pr.py @ 59d542c

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since 59d542c was 59d542c, checked in by Jae Cho <jhjcho@…>, 13 years ago
fixed error bar changes
Property mode set to `100644`
File size: 49.9 KB

Line
1
2	################################################################################
3	#This software was developed by the University of Tennessee as part of the
4	#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
5	#project funded by the US National Science Foundation.
6	#
7	#See the license text in license.txt
8	#
9	#copyright 2009, University of Tennessee
10	################################################################################
11
12
13	# Make sure the option of saving each curve is available
14	# Use the I(q) curve as input and compare the output to P(r)
15
16	import os
17	import sys
18	import wx
19	import logging
20	import time
21	import copy
22	import math
23	import numpy
24	import pylab
25	from sans.guiframe.dataFitting import Data1D
26	from sans.guiframe.events import NewPlotEvent
27	from sans.guiframe.events import StatusEvent
28	from sans.guiframe.gui_style import GUIFRAME_ID
29	from sans.pr.invertor import Invertor
30	from DataLoader.loader import Loader
31
32	from pr_widgets import load_error
33	from sans.guiframe.plugin_base import PluginBase
34
35
36	PR_FIT_LABEL = r"$P_{fit}(r)$"
37	PR_LOADED_LABEL = r"$P_{loaded}(r)$"
38	IQ_DATA_LABEL = r"$I_{obs}(q)$"
39	IQ_FIT_LABEL = r"$I_{fit}(q)$"
40	IQ_SMEARED_LABEL = r"$I_{smeared}(q)$"
41	GROUP_ID_IQ_DATA = r"$I_{obs}(q)$"
42	GROUP_ID_PR_FIT = r"$P_{fit}(r)$"
43
44
45
46	class Plugin(PluginBase):
47	"""
48	"""
49	DEFAULT_ALPHA = 0.0001
50	DEFAULT_NFUNC = 10
51	DEFAULT_DMAX = 140.0
52
53	def __init__(self, standalone=True):
54	PluginBase.__init__(self, name="Pr Inversion", standalone=standalone)
55	## Simulation window manager
56	self.simview = None
57
58	## State data
59	self.alpha = self.DEFAULT_ALPHA
60	self.nfunc = self.DEFAULT_NFUNC
61	self.max_length = self.DEFAULT_DMAX
62	self.q_min = None
63	self.q_max = None
64	self.has_bck = False
65	self.slit_height = 0
66	self.slit_width = 0
67	## Remember last plottable processed
68	self.last_data = "sphere_60_q0_2.txt"
69	self._current_file_data = None
70	## Time elapsed for last computation [sec]
71	# Start with a good default
72	self.elapsed = 0.022
73	self.iq_data_shown = False
74
75	## Current invertor
76	self.invertor = None
77	self.pr = None
78	self.data_id = IQ_DATA_LABEL
79	# Copy of the last result in case we need to display it.
80	self._last_pr = None
81	self._last_out = None
82	self._last_cov = None
83	## Calculation thread
84	self.calc_thread = None
85	## Estimation thread
86	self.estimation_thread = None
87	## Result panel
88	self.control_panel = None
89	## Currently views plottable
90	self.current_plottable = None
91	## Number of P(r) points to display on the output plot
92	self._pr_npts = 51
93	## Flag to let the plug-in know that it is running standalone
94	self.standalone = standalone
95	self._normalize_output = False
96	self._scale_output_unity = False
97
98	## List of added P(r) plots
99	self._added_plots = {}
100	self._default_Iq = {}
101	self.list_plot_id = []
102
103	# Associate the inversion state reader with .prv files
104	from inversion_state import Reader
105
106	# Create a CanSAS/Pr reader
107	self.state_reader = Reader(self.set_state)
108	self._extensions = '.prv'
109	l = Loader()
110	l.associate_file_reader('.prv', self.state_reader)
111	l.associate_file_reader(".svs", self.state_reader)
112
113	# Log startup
114	logging.info("Pr(r) plug-in started")
115
116	def delete_data(self, data_id):
117	"""
118	delete the data association with prview
119	"""
120	if self.calc_thread is not None and self.calc_thread.isrunning():
121	msg = "Data in use in Prview\n"
122	if self.estimation_thread is not None and \
123	self.estimation_thread.isrunning():
124	msg = "Data in use in Prview\n"
125	self.control_panel.clear_panel()
126
127	def get_data(self):
128	"""
129	"""
130	return self.current_plottable
131
132	def set_state(self, state=None, datainfo=None):
133	"""
134	Call-back method for the inversion state reader.
135	This method is called when a .prv file is loaded.
136
137	:param state: InversionState object
138	:param datainfo: Data1D object [optional]
139
140	"""
141	try:
142	if datainfo.__class__.__name__ == 'list':
143	if len(datainfo) >= 1:
144	data = datainfo[0]
145	else:
146	data = None
147	else:
148	data = datainfo
149	if data is None:
150	raise RuntimeError, "Pr.set_state: datainfo parameter cannot be None in standalone mode"
151
152	# Ensuring that plots are coordinated correctly
153	t = time.localtime(data.meta_data['prstate'].timestamp)
154	time_str = time.strftime("%b %d %H:%M", t)
155
156	# Check that no time stamp is already appended
157	max_char = data.meta_data['prstate'].file.find("[")
158	if max_char < 0:
159	max_char = len(data.meta_data['prstate'].file)
160
161	datainfo.meta_data['prstate'].file = data.meta_data['prstate'].file[0:max_char] +' [' + time_str + ']'
162	data.filename = data.meta_data['prstate'].file
163	# TODO:
164	#remove this call when state save all information about the gui data
165	# such as ID , Group_ID, etc...
166	#make self.current_plottable = datainfo directly
167	self.current_plottable = self.parent.create_gui_data(data,None)
168	self.current_plottable.group_id = data.meta_data['prstate'].file
169
170	# Make sure the user sees the P(r) panel after loading
171	#self.parent.set_perspective(self.perspective)
172	self.on_perspective(event=None)
173	# Load the P(r) results
174	#state = self.state_reader.get_state()
175	data_dict = {self.current_plottable.id:self.current_plottable}
176	self.parent.add_data(data_list=data_dict)
177	wx.PostEvent(self.parent, NewPlotEvent(plot=self.current_plottable,
178	title=self.current_plottable.title))
179	self.control_panel.set_state(state)
180	except:
181	logging.error("prview.set_state: %s" % sys.exc_value)
182
183
184	def help(self, evt):
185	"""
186	Show a general help dialog.
187
188	:TODO: replace the text with a nice image
189
190	"""
191	from inversion_panel import HelpDialog
192	dialog = HelpDialog(None, -1)
193	if dialog.ShowModal() == wx.ID_OK:
194	dialog.Destroy()
195	else:
196	dialog.Destroy()
197
198	def _fit_pr(self, evt):
199	"""
200	"""
201	from sans.pr.invertor import Invertor
202	# Generate P(r) for sphere
203	radius = 60.0
204	d_max = 2*radius
205
206	r = pylab.arange(0.01, d_max, d_max/51.0)
207	M = len(r)
208	y = numpy.zeros(M)
209	pr_err = numpy.zeros(M)
210
211	sum = 0.0
212	for j in range(M):
213	value = self.pr_theory(r[j], radius)
214	sum += value
215	y[j] = value
216	pr_err[j] = math.sqrt(y[j])
217
218	y = y/sum*d_max/len(r)
219
220	# Perform fit
221	pr = Invertor()
222	pr.d_max = d_max
223	pr.alpha = 0
224	pr.x = r
225	pr.y = y
226	pr.err = pr_err
227	out, cov = pr.pr_fit()
228	for i in range(len(out)):
229	print "%g +- %g" % (out[i], math.sqrt(cov[i][i]))
230
231	# Show input P(r)
232	title = "Pr"
233	new_plot = Data1D(pr.x, pr.y, dy=pr.err)
234	new_plot.name = "P_{obs}(r)"
235	new_plot.xaxis("\\rm{r}", 'A')
236	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
237	new_plot.group_id = "P_{obs}(r)"
238	new_plot.id = "P_{obs}(r)"
239	new_plot.title = title
240	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
241	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
242
243	# Show P(r) fit
244	self.show_pr(out, pr)
245
246	# Show I(q) fit
247	q = pylab.arange(0.001, 0.1, 0.01/51.0)
248	self.show_iq(out, pr, q)
249
250	def show_shpere(self, x, radius=70.0, x_range=70.0):
251	"""
252	"""
253	# Show P(r)
254	y_true = numpy.zeros(len(x))
255
256	sum_true = 0.0
257	for i in range(len(x)):
258	y_true[i] = self.pr_theory(x[i], radius)
259	sum_true += y_true[i]
260
261	y_true = y_true/sum_true*x_range/len(x)
262
263	# Show the theory P(r)
264	new_plot = Data1D(x, y_true)
265	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
266	new_plot.name = "P_{true}(r)"
267	new_plot.xaxis("\\rm{r}", 'A')
268	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
269	new_plot.id = "P_{true}(r)"
270	new_plot.group_id = "P_{true}(r)"
271	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
272	#Put this call in plottables/guitools
273	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot,
274	title="Sphere P(r)"))
275
276
277	def get_npts(self):
278	"""
279	Returns the number of points in the I(q) data
280	"""
281	try:
282	return len(self.pr.x)
283	except:
284	return 0
285
286	def show_iq(self, out, pr, q=None):
287	"""
288	"""
289	qtemp = pr.x
290	if not q==None:
291	qtemp = q
292
293	# Make a plot
294	maxq = -1
295	for q_i in qtemp:
296	if q_i>maxq:
297	maxq=q_i
298
299	minq = 0.001
300
301	# Check for user min/max
302	if not pr.q_min==None:
303	minq = pr.q_min
304	if not pr.q_max==None:
305	maxq = pr.q_max
306
307	x = pylab.arange(minq, maxq, maxq/301.0)
308	y = numpy.zeros(len(x))
309	err = numpy.zeros(len(x))
310	for i in range(len(x)):
311	value = pr.iq(out, x[i])
312	y[i] = value
313	try:
314	err[i] = math.sqrt(math.fabs(value))
315	except:
316	err[i] = 1.0
317	print "Error getting error", value, x[i]
318
319	new_plot = Data1D(x, y)
320	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
321	new_plot.name = IQ_FIT_LABEL
322	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
323	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
324	title = "I(q)"
325	new_plot.title = title
326
327	# If we have a group ID, use it
328	if pr.info.has_key("plot_group_id"):
329	new_plot.group_id = pr.info["plot_group_id"]
330	new_plot.id = IQ_FIT_LABEL
331	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
332
333	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
334
335	# If we have used slit smearing, plot the smeared I(q) too
336	if pr.slit_width>0 or pr.slit_height>0:
337	x = pylab.arange(minq, maxq, maxq/301.0)
338	y = numpy.zeros(len(x))
339	err = numpy.zeros(len(x))
340	for i in range(len(x)):
341	value = pr.iq_smeared(out, x[i])
342	y[i] = value
343	try:
344	err[i] = math.sqrt(math.fabs(value))
345	except:
346	err[i] = 1.0
347	print "Error getting error", value, x[i]
348
349	new_plot = Data1D(x, y)
350	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
351	new_plot.name = IQ_SMEARED_LABEL
352	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
353	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
354	# If we have a group ID, use it
355	if pr.info.has_key("plot_group_id"):
356	new_plot.group_id = pr.info["plot_group_id"]
357	new_plot.id = IQ_SMEARED_LABEL
358	new_plot.title = title
359	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
360	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
361
362	def _on_pr_npts(self, evt):
363	"""
364	Redisplay P(r) with a different number of points
365	"""
366	from inversion_panel import PrDistDialog
367	dialog = PrDistDialog(None, -1)
368	dialog.set_content(self._pr_npts)
369	if dialog.ShowModal() == wx.ID_OK:
370	self._pr_npts= dialog.get_content()
371	dialog.Destroy()
372	self.show_pr(self._last_out, self._last_pr, self._last_cov)
373	else:
374	dialog.Destroy()
375
376
377	def show_pr(self, out, pr, cov=None):
378	"""
379	"""
380	# Show P(r)
381	x = pylab.arange(0.0, pr.d_max, pr.d_max/self._pr_npts)
382
383	y = numpy.zeros(len(x))
384	dy = numpy.zeros(len(x))
385	y_true = numpy.zeros(len(x))
386
387	sum = 0.0
388	pmax = 0.0
389	cov2 = numpy.ascontiguousarray(cov)
390
391	for i in range(len(x)):
392	if cov2==None:
393	value = pr.pr(out, x[i])
394	else:
395	(value, dy[i]) = pr.pr_err(out, cov2, x[i])
396	sum += value*pr.d_max/len(x)
397
398	# keep track of the maximum P(r) value
399	if value>pmax:
400	pmax = value
401
402	y[i] = value
403
404	if self._normalize_output==True:
405	y = y/sum
406	dy = dy/sum
407	elif self._scale_output_unity==True:
408	y = y/pmax
409	dy = dy/pmax
410
411	if cov2==None:
412	new_plot = Data1D(x, y)
413	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
414	else:
415	new_plot = Data1D(x, y, dy=dy)
416	new_plot.name = PR_FIT_LABEL
417	new_plot.xaxis("\\rm{r}", 'A')
418	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
419	new_plot.title = "P(r) fit"
420	new_plot.id = PR_FIT_LABEL
421	# Make sure that the plot is linear
422	new_plot.xtransform = "x"
423	new_plot.ytransform = "y"
424	new_plot.group_id = GROUP_ID_PR_FIT
425	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
426	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="P(r) fit"))
427	return x, pr.d_max
428
429	def load(self, data):
430	"""
431	Load data. This will eventually be replaced
432	by our standard DataLoader class.
433	"""
434	class FileData:
435	x = None
436	y = None
437	err = None
438	path = None
439
440	def __init__(self, path):
441	self.path = path
442
443	self._current_file_data = FileData(data.path)
444
445	# Use data loader to load file
446	#dataread = Loader().load(path)
447	dataread = data
448	# Notify the user if we could not read the file
449	if dataread is None:
450	raise RuntimeError, "Invalid data"
451
452	x = None
453	y = None
454	err = None
455	if dataread.__class__.__name__ == 'Data1D':
456	x = dataread.x
457	y = dataread.y
458	err = dataread.dy
459	else:
460	if isinstance(dataread, list) and len(dataread)>0:
461	x = dataread[0].x
462	y = dataread[0].y
463	err = dataread[0].dy
464	msg = "PrView only allows a single data set at a time. "
465	msg += "Only the first data set was loaded."
466	wx.PostEvent(self.parent, StatusEvent(status=msg))
467	else:
468	if dataread is None:
469	return x, y, err
470	raise RuntimeError, "This tool can only read 1D data"
471
472	self._current_file_data.x = x
473	self._current_file_data.y = y
474	self._current_file_data.err = err
475	return x, y, err
476
477	def load_columns(self, path = "sphere_60_q0_2.txt"):
478	"""
479	Load 2- or 3- column ascii
480	"""
481	import numpy, math, sys
482	# Read the data from the data file
483	data_x = numpy.zeros(0)
484	data_y = numpy.zeros(0)
485	data_err = numpy.zeros(0)
486	scale = None
487	min_err = 0.0
488	if not path == None:
489	input_f = open(path,'r')
490	buff = input_f.read()
491	lines = buff.split('\n')
492	for line in lines:
493	try:
494	toks = line.split()
495	x = float(toks[0])
496	y = float(toks[1])
497	if len(toks)>2:
498	err = float(toks[2])
499	else:
500	if scale==None:
501	scale = 0.05*math.sqrt(y)
502	#scale = 0.05/math.sqrt(y)
503	min_err = 0.01*y
504	err = scale*math.sqrt(y)+min_err
505	#err = 0
506
507	data_x = numpy.append(data_x, x)
508	data_y = numpy.append(data_y, y)
509	data_err = numpy.append(data_err, err)
510	except:
511	pass
512
513	if not scale==None:
514	message = "The loaded file had no error bars, statistical errors are assumed."
515	wx.PostEvent(self.parent, StatusEvent(status=message))
516	else:
517	wx.PostEvent(self.parent, StatusEvent(status=''))
518
519	return data_x, data_y, data_err
520
521	def load_abs(self, path):
522	"""
523	Load an IGOR .ABS reduced file
524
525	:param path: file path
526
527	:return: x, y, err vectors
528
529	"""
530	import numpy, math, sys
531	# Read the data from the data file
532	data_x = numpy.zeros(0)
533	data_y = numpy.zeros(0)
534	data_err = numpy.zeros(0)
535	scale = None
536	min_err = 0.0
537
538	data_started = False
539	if not path == None:
540	input_f = open(path,'r')
541	buff = input_f.read()
542	lines = buff.split('\n')
543	for line in lines:
544	if data_started==True:
545	try:
546	toks = line.split()
547	x = float(toks[0])
548	y = float(toks[1])
549	if len(toks)>2:
550	err = float(toks[2])
551	else:
552	if scale==None:
553	scale = 0.05*math.sqrt(y)
554	#scale = 0.05/math.sqrt(y)
555	min_err = 0.01*y
556	err = scale*math.sqrt(y)+min_err
557	#err = 0
558
559	data_x = numpy.append(data_x, x)
560	data_y = numpy.append(data_y, y)
561	data_err = numpy.append(data_err, err)
562	except:
563	pass
564	elif line.find("The 6 columns")>=0:
565	data_started = True
566
567	if not scale==None:
568	message = "The loaded file had no error bars, statistical errors are assumed."
569	wx.PostEvent(self.parent, StatusEvent(status=message))
570	else:
571	wx.PostEvent(self.parent, StatusEvent(status=''))
572
573	return data_x, data_y, data_err
574
575	def pr_theory(self, r, R):
576	"""
577	"""
578	if r<=2*R:
579	return 12.0* ((0.5r/R)2) ((1.0-0.5r/R)2) ( 2.0 + 0.5*r/R )
580	else:
581	return 0.0
582
583	def get_context_menu(self, plotpanel=None):
584	"""
585	Get the context menu items available for P(r)
586
587	:param graph: the Graph object to which we attach the context menu
588
589	:return: a list of menu items with call-back function
590
591	"""
592	graph = plotpanel.graph
593	# Look whether this Graph contains P(r) data
594	if graph.selected_plottable not in plotpanel.plots:
595	return []
596	item = plotpanel.plots[graph.selected_plottable]
597	if item.id == PR_FIT_LABEL:
598	#add_data_hint = "Load a data file and display it on this plot"
599	#["Add P(r) data",add_data_hint , self._on_add_data],
600	change_n_hint = "Change the number of"
601	change_n_hint += " points on the P(r) output"
602	change_n_label = "Change number of P(r) points"
603	m_list = [[change_n_label, change_n_hint , self._on_pr_npts]]
604
605	if self._scale_output_unity or self._normalize_output:
606	hint = "Let the output P(r) keep the scale of the data"
607	m_list.append(["Disable P(r) scaling", hint,
608	self._on_disable_scaling])
609	if not self._scale_output_unity:
610	m_list.append(["Scale P_max(r) to unity",
611	"Scale P(r) so that its maximum is 1",
612	self._on_scale_unity])
613	if not self._normalize_output:
614	m_list.append(["Normalize P(r) to unity",
615	"Normalize the integral of P(r) to 1",
616	self._on_normalize])
617
618	return m_list
619
620	elif item.id in [PR_LOADED_LABEL, IQ_DATA_LABEL, IQ_FIT_LABEL,
621	IQ_SMEARED_LABEL]:
622	return []
623	elif item.id == graph.selected_plottable:
624	if not self.standalone and issubclass(item.__class__, Data1D):
625	return [["Compute P(r)",
626	"Compute P(r) from distribution",
627	self._on_context_inversion]]
628
629	return []
630
631	def _on_disable_scaling(self, evt):
632	"""
633	Disable P(r) scaling
634
635	:param evt: Menu event
636
637	"""
638	self._normalize_output = False
639	self._scale_output_unity = False
640	self.show_pr(self._last_out, self._last_pr, self._last_cov)
641
642	# Now replot the original added data
643	for plot in self._added_plots:
644	self._added_plots[plot].y = numpy.copy(self._default_Iq[plot])
645	wx.PostEvent(self.parent,
646	NewPlotEvent(plot=self._added_plots[plot],
647	title=self._added_plots[plot].name,
648	update=True))
649
650	# Need the update flag in the NewPlotEvent to protect against
651	# the plot no longer being there...
652
653	def _on_normalize(self, evt):
654	"""
655	Normalize the area under the P(r) curve to 1.
656	This operation is done for all displayed plots.
657
658	:param evt: Menu event
659
660	"""
661	self._normalize_output = True
662	self._scale_output_unity = False
663
664	self.show_pr(self._last_out, self._last_pr, self._last_cov)
665
666	# Now scale the added plots too
667	for plot in self._added_plots:
668	sum = numpy.sum(self._added_plots[plot].y)
669	npts = len(self._added_plots[plot].x)
670	sum *= self._added_plots[plot].x[npts-1]/npts
671	y = self._added_plots[plot].y/sum
672
673	new_plot = Data1D(self._added_plots[plot].x, y)
674	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
675	new_plot.group_id = self._added_plots[plot].group_id
676	new_plot.id = self._added_plots[plot].id
677	new_plot.title = self._added_plots[plot].title
678	new_plot.name = self._added_plots[plot].name
679	new_plot.xaxis("\\rm{r}", 'A')
680	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
681	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
682	wx.PostEvent(self.parent,
683	NewPlotEvent(plot=new_plot, update=True,
684	title=self._added_plots[plot].name))
685
686	def _on_scale_unity(self, evt):
687	"""
688	Scale the maximum P(r) value on each displayed plot to 1.
689
690	:param evt: Menu event
691
692	"""
693	self._scale_output_unity = True
694	self._normalize_output = False
695
696	self.show_pr(self._last_out, self._last_pr, self._last_cov)
697
698	# Now scale the added plots too
699	for plot in self._added_plots:
700	_max = 0
701	for y in self._added_plots[plot].y:
702	if y>_max:
703	_max = y
704	y = self._added_plots[plot].y/_max
705
706	new_plot = Data1D(self._added_plots[plot].x, y)
707	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
708	new_plot.name = self._added_plots[plot].name
709	new_plot.xaxis("\\rm{r}", 'A')
710	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
711	self.parent.update_theory(data_id=self.data_id,theory=new_plot)
712	wx.PostEvent(self.parent,
713	NewPlotEvent(plot=new_plot, update=True,
714	title=self._added_plots[plot].name))
715
716
717	def _on_add_data(self, evt):
718	"""
719	Add a data curve to the plot
720
721	:WARNING: this will be removed once guiframe.plotting has
722	its full functionality
723	"""
724	path = self.choose_file()
725	if path==None:
726	return
727
728	#x, y, err = self.parent.load_ascii_1D(path)
729	# Use data loader to load file
730	try:
731	dataread = Loader().load(path)
732	x = None
733	y = None
734	err = None
735	if dataread.__class__.__name__ == 'Data1D':
736	x = dataread.x
737	y = dataread.y
738	err = dataread.dy
739	else:
740	if isinstance(dataread, list) and len(dataread)>0:
741	x = dataread[0].x
742	y = dataread[0].y
743	err = dataread[0].dy
744	msg = "PrView only allows a single data set at a time. "
745	msg += "Only the first data set was loaded."
746	wx.PostEvent(self.parent, StatusEvent(status=msg))
747	else:
748	msg = "This tool can only read 1D data"
749	wx.PostEvent(self.parent, StatusEvent(status=msg))
750	return
751
752	except:
753	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
754	return
755
756	filename = os.path.basename(path)
757	new_plot = Data1D(x, y)
758	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
759	new_plot.name = filename
760	new_plot.xaxis("\\rm{r}", 'A')
761	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
762	# Store a ref to the plottable for later use
763	self._added_plots[filename] = new_plot
764	self._default_Iq[filename] = numpy.copy(y)
765	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=filename))
766
767	def start_thread(self):
768	"""
769	"""
770	from pr_thread import CalcPr
771	from copy import deepcopy
772
773	# If a thread is already started, stop it
774	if self.calc_thread != None and self.calc_thread.isrunning():
775	self.calc_thread.stop()
776
777	pr = self.pr.clone()
778	self.calc_thread = CalcPr(pr, self.nfunc,
779	error_func=self._thread_error,
780	completefn=self._completed, updatefn=None)
781	self.calc_thread.queue()
782	self.calc_thread.ready(2.5)
783
784	def _thread_error(self, error):
785	"""
786	"""
787	wx.PostEvent(self.parent, StatusEvent(status=error))
788
789	def _estimate_completed(self, alpha, message, elapsed):
790	"""
791	Parameter estimation completed,
792	display the results to the user
793
794	:param alpha: estimated best alpha
795	:param elapsed: computation time
796
797	"""
798	# Save useful info
799	self.elapsed = elapsed
800	self.control_panel.alpha_estimate = alpha
801	if not message==None:
802	wx.PostEvent(self.parent, StatusEvent(status=str(message)))
803	self.perform_estimateNT()
804
805	def _estimateNT_completed(self, nterms, alpha, message, elapsed):
806	"""
807	Parameter estimation completed,
808	display the results to the user
809
810	:param alpha: estimated best alpha
811	:param nterms: estimated number of terms
812	:param elapsed: computation time
813
814	"""
815	# Save useful info
816	self.elapsed = elapsed
817	self.control_panel.nterms_estimate = nterms
818	self.control_panel.alpha_estimate = alpha
819	if not message==None:
820	wx.PostEvent(self.parent, StatusEvent(status=str(message)))
821
822	def _completed(self, out, cov, pr, elapsed):
823	"""
824	Method called with the results when the inversion
825	is done
826
827	:param out: output coefficient for the base functions
828	:param cov: covariance matrix
829	:param pr: Invertor instance
830	:param elapsed: time spent computing
831
832	"""
833	from copy import deepcopy
834	# Save useful info
835	self.elapsed = elapsed
836	# Keep a copy of the last result
837	self._last_pr = pr.clone()
838	self._last_out = out
839	self._last_cov = cov
840
841	# Save Pr invertor
842	self.pr = pr
843
844	#message = "Computation completed in"
845	#message += %g seconds [chi2=%g]" % (elapsed, pr.chi2)
846	#wx.PostEvent(self.parent, StatusEvent(status=message))
847
848	cov = numpy.ascontiguousarray(cov)
849
850	# Show result on control panel
851	self.control_panel.chi2 = pr.chi2
852	self.control_panel.elapsed = elapsed
853	self.control_panel.oscillation = pr.oscillations(out)
854	#print "OSCILL", pr.oscillations(out)
855	#print "PEAKS:", pr.get_peaks(out)
856	self.control_panel.positive = pr.get_positive(out)
857	self.control_panel.pos_err = pr.get_pos_err(out, cov)
858	self.control_panel.rg = pr.rg(out)
859	self.control_panel.iq0 = pr.iq0(out)
860	self.control_panel.bck = pr.background
861
862	if False:
863	for i in range(len(out)):
864	try:
865	print "%d: %g +- %g" % (i, out[i],
866	math.sqrt(math.fabs(cov[i][i])))
867	except:
868	print sys.exc_value
869	print "%d: %g +- ?" % (i, out[i])
870
871	# Make a plot of I(q) data
872	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
873	new_plot.name = IQ_DATA_LABEL
874	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
875	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
876	if pr.info.has_key("plot_group_id"):
877	new_plot.group_id = pr.info["plot_group_id"]
878	new_plot.id = self.data_id
879	self.parent.update_theory(data_id=self.data_id,
880	theory=new_plot)
881	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
882
883	# Show I(q) fit
884	self.show_iq(out, self.pr)
885
886	# Show P(r) fit
887	x_values, x_range = self.show_pr(out, self.pr, cov)
888
889	# Popup result panel
890	#result_panel = InversionResults(self.parent,
891	#-1, style=wx.RAISED_BORDER)
892
893	def show_data(self, path=None, data=None, reset=False):
894	"""
895	Show data read from a file
896
897	:param path: file path
898	:param reset: if True all other plottables will be cleared
899
900	"""
901	#if path is not None:
902	if data is not None:
903	try:
904	pr = self._create_file_pr(data)
905	except:
906	status = "Problem reading data: %s" % sys.exc_value
907	wx.PostEvent(self.parent, StatusEvent(status=status))
908	raise RuntimeError, status
909
910	# If the file contains nothing, just return
911	if pr is None:
912	raise RuntimeError, "Loaded data is invalid"
913
914	self.pr = pr
915
916	# Make a plot of I(q) data
917	if self.pr.err == None:
918	new_plot = Data1D(self.pr.x, self.pr.y)
919	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
920	else:
921	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
922	new_plot.name = IQ_DATA_LABEL
923	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
924	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
925	new_plot.interactive = True
926	new_plot.group_id = GROUP_ID_IQ_DATA
927	new_plot.id = self.data_id
928	new_plot.title = "I(q)"
929	wx.PostEvent(self.parent,
930	NewPlotEvent(plot=new_plot, title="I(q)", reset=reset))
931
932	self.current_plottable = new_plot
933	# Get Q range
934	self.control_panel.q_min = self.pr.x.min()
935	self.control_panel.q_max = self.pr.x.max()
936
937	def save_data(self, filepath, prstate=None):
938	"""
939	Save data in provided state object.
940
941	:TODO: move the state code away from inversion_panel and move it here.
942	Then remove the "prstate" input and make this method private.
943
944	:param filepath: path of file to write to
945	:param prstate: P(r) inversion state
946
947	"""
948	#TODO: do we need this or can we use DataLoader.loader.save directly?
949
950	# Add output data and coefficients to state
951	prstate.coefficients = self._last_out
952	prstate.covariance = self._last_cov
953
954	# Write the output to file
955	# First, check that the data is of the right type
956	if issubclass(self.current_plottable.__class__,
957	DataLoader.data_info.Data1D):
958	self.state_reader.write(filepath, self.current_plottable, prstate)
959	else:
960	msg = "pr.save_data: the data being saved is not a"
961	msg += " DataLoader.data_info.Data1D object"
962	raise RuntimeError, msg
963
964
965	def setup_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None,
966	bck=False, height=0, width=0):
967	"""
968	"""
969	self.alpha = alpha
970	self.nfunc = nfunc
971	self.max_length = d_max
972	self.q_min = q_min
973	self.q_max = q_max
974	self.has_bck = bck
975	self.slit_height = height
976	self.slit_width = width
977
978	try:
979	pr = self._create_plot_pr()
980	if not pr==None:
981	self.pr = pr
982	self.perform_inversion()
983	except:
984	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
985
986	def estimate_plot_inversion(self, alpha, nfunc, d_max,
987	q_min=None, q_max=None,
988	bck=False, height=0, width=0):
989	"""
990	"""
991	self.alpha = alpha
992	self.nfunc = nfunc
993	self.max_length = d_max
994	self.q_min = q_min
995	self.q_max = q_max
996	self.has_bck = bck
997	self.slit_height = height
998	self.slit_width = width
999
1000	try:
1001	pr = self._create_plot_pr()
1002	if not pr==None:
1003	self.pr = pr
1004	self.perform_estimate()
1005	except:
1006	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1007
1008	def _create_plot_pr(self, estimate=False):
1009	"""
1010	Create and prepare invertor instance from
1011	a plottable data set.
1012
1013	:param path: path of the file to read in
1014
1015	"""
1016	# Sanity check
1017	if self.current_plottable is None:
1018	msg = "Please load a valid data set before proceeding."
1019	wx.PostEvent(self.parent, StatusEvent(status=msg))
1020	return None
1021
1022	# Get the data from the chosen data set and perform inversion
1023	pr = Invertor()
1024	pr.d_max = self.max_length
1025	pr.alpha = self.alpha
1026	pr.q_min = self.q_min
1027	pr.q_max = self.q_max
1028	pr.x = self.current_plottable.x
1029	pr.y = self.current_plottable.y
1030	pr.has_bck = self.has_bck
1031	pr.slit_height = self.slit_height
1032	pr.slit_width = self.slit_width
1033
1034	# Keep track of the plot window title to ensure that
1035	# we can overlay the plots
1036	pr.info["plot_group_id"] = self.current_plottable.group_id
1037
1038	# Fill in errors if none were provided
1039	err = self.current_plottable.dy
1040	all_zeros = True
1041	if err == None:
1042	err = numpy.zeros(len(pr.y))
1043	else:
1044	for i in range(len(err)):
1045	if err[i]>0:
1046	all_zeros = False
1047
1048	if all_zeros:
1049	scale = None
1050	min_err = 0.0
1051	for i in range(len(pr.y)):
1052	# Scale the error so that we can fit over several decades of Q
1053	if scale==None:
1054	scale = 0.05*math.sqrt(pr.y[i])
1055	min_err = 0.01*pr.y[i]
1056	err[i] = scale*math.sqrt( math.fabs(pr.y[i]) ) + min_err
1057	message = "The loaded file had no error bars, "
1058	message += "statistical errors are assumed."
1059	wx.PostEvent(self.parent, StatusEvent(status=message))
1060
1061	pr.err = err
1062
1063	return pr
1064
1065
1066	def setup_file_inversion(self, alpha, nfunc, d_max, data,
1067	path=None, q_min=None, q_max=None,
1068	bck=False, height=0, width=0):
1069	"""
1070	"""
1071	self.alpha = alpha
1072	self.nfunc = nfunc
1073	self.max_length = d_max
1074	self.q_min = q_min
1075	self.q_max = q_max
1076	self.has_bck = bck
1077	self.slit_height = height
1078	self.slit_width = width
1079
1080	try:
1081	#pr = self._create_file_pr(path)
1082	pr = self._create_file_pr(data)
1083	if not pr==None:
1084	self.pr = pr
1085	self.perform_inversion()
1086	except:
1087	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1088
1089	def estimate_file_inversion(self, alpha, nfunc, d_max, data,
1090	path=None, q_min=None, q_max=None,
1091	bck=False, height=0, width=0):
1092	"""
1093	"""
1094	self.alpha = alpha
1095	self.nfunc = nfunc
1096	self.max_length = d_max
1097	self.q_min = q_min
1098	self.q_max = q_max
1099	self.has_bck = bck
1100	self.slit_height = height
1101	self.slit_width = width
1102
1103	try:
1104	pr = self._create_file_pr(data)
1105	#pr = self._create_file_pr(path)
1106	if not pr is None:
1107	self.pr = pr
1108	self.perform_estimate()
1109	except:
1110	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1111
1112
1113	def _create_file_pr(self, data):
1114	"""
1115	Create and prepare invertor instance from
1116	a file data set.
1117
1118	:param path: path of the file to read in
1119
1120	"""
1121	# Load data
1122	#if os.path.isfile(path):
1123	"""
1124	if self._current_file_data is not None \
1125	and self._current_file_data.path==path:
1126	# Protect against corrupted data from
1127	# previous failed load attempt
1128	if self._current_file_data.x is None:
1129	return None
1130	x = self._current_file_data.x
1131	y = self._current_file_data.y
1132	err = self._current_file_data.err
1133
1134	message = "The data from this file has already been loaded."
1135	wx.PostEvent(self.parent, StatusEvent(status=message))
1136	else:
1137	"""
1138	# Reset the status bar so that we don't get mixed up
1139	# with old messages.
1140	#TODO: refactor this into a proper status handling
1141	wx.PostEvent(self.parent, StatusEvent(status=''))
1142	try:
1143	class FileData:
1144	x = None
1145	y = None
1146	err = None
1147	path = None
1148	def __init__(self, path):
1149	self.path = path
1150
1151	self._current_file_data = FileData(data.path)
1152	self._current_file_data.x = data.x
1153	self._current_file_data.y = data.y
1154	self._current_file_data.err = data.dy
1155	x, y, err = data.x, data.y, data.dy
1156	except:
1157	load_error(sys.exc_value)
1158	return None
1159
1160	# If the file contains no data, just return
1161	if x is None or len(x) == 0:
1162	load_error("The loaded file contains no data")
1163	return None
1164
1165	# If we have not errors, add statistical errors
1166	if y is not None:
1167	if err == None or numpy.all(err) == 0:
1168	err = numpy.zeros(len(y))
1169	scale = None
1170	min_err = 0.0
1171	for i in range(len(y)):
1172	# Scale the error so that we can fit over several decades of Q
1173	if scale == None:
1174	scale = 0.05 * math.sqrt(y[i])
1175	min_err = 0.01 * y[i]
1176	err[i] = scale * math.sqrt(math.fabs(y[i])) + min_err
1177	message = "The loaded file had no error bars, "
1178	message += "statistical errors are assumed."
1179	wx.PostEvent(self.parent, StatusEvent(status=message))
1180
1181	try:
1182	# Get the data from the chosen data set and perform inversion
1183	pr = Invertor()
1184	pr.d_max = self.max_length
1185	pr.alpha = self.alpha
1186	pr.q_min = self.q_min
1187	pr.q_max = self.q_max
1188	pr.x = x
1189	pr.y = y
1190	pr.err = err
1191	pr.has_bck = self.has_bck
1192	pr.slit_height = self.slit_height
1193	pr.slit_width = self.slit_width
1194	return pr
1195	except:
1196	load_error(sys.exc_value)
1197	return None
1198
1199	def perform_estimate(self):
1200	"""
1201	"""
1202	from pr_thread import EstimatePr
1203	from copy import deepcopy
1204
1205	# If a thread is already started, stop it
1206	if self.estimation_thread != None and \
1207	self.estimation_thread.isrunning():
1208	self.estimation_thread.stop()
1209
1210	pr = self.pr.clone()
1211	self.estimation_thread = EstimatePr(pr, self.nfunc,
1212	error_func=self._thread_error,
1213	completefn = self._estimate_completed,
1214	updatefn = None)
1215	self.estimation_thread.queue()
1216	self.estimation_thread.ready(2.5)
1217
1218	def perform_estimateNT(self):
1219	"""
1220	"""
1221	from pr_thread import EstimateNT
1222	from copy import deepcopy
1223
1224	# If a thread is already started, stop it
1225	if self.estimation_thread != None and self.estimation_thread.isrunning():
1226	self.estimation_thread.stop()
1227
1228	pr = self.pr.clone()
1229	# Skip the slit settings for the estimation
1230	# It slows down the application and it doesn't change the estimates
1231	pr.slit_height = 0.0
1232	pr.slit_width = 0.0
1233	self.estimation_thread = EstimateNT(pr, self.nfunc,
1234	error_func=self._thread_error,
1235	completefn = self._estimateNT_completed,
1236	updatefn = None)
1237	self.estimation_thread.queue()
1238	self.estimation_thread.ready(2.5)
1239
1240	def perform_inversion(self):
1241	"""
1242	"""
1243	# Time estimate
1244	#estimated = self.elapsedself.nfunc*2
1245	#message = "Computation time may take up to %g seconds" % self.elapsed
1246	#wx.PostEvent(self.parent, StatusEvent(status=message))
1247
1248	# Start inversion thread
1249	self.start_thread()
1250	return
1251
1252	out, cov = self.pr.lstsq(self.nfunc)
1253
1254	# Save useful info
1255	self.elapsed = self.pr.elapsed
1256
1257	for i in range(len(out)):
1258	try:
1259	print "%d: %g +- %g" % (i, out[i],
1260	math.sqrt(math.fabs(cov[i][i])))
1261	except:
1262	print "%d: %g +- ?" % (i, out[i])
1263
1264	# Make a plot of I(q) data
1265	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
1266	new_plot.name = "I_{obs}(q)"
1267	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
1268	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
1269	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
1270	# Show I(q) fit
1271	self.show_iq(out, self.pr)
1272	# Show P(r) fit
1273	x_values, x_range = self.show_pr(out, self.pr, cov=cov)
1274
1275	def _on_context_inversion(self, event):
1276	"""
1277	"""
1278	panel = event.GetEventObject()
1279	Plugin.on_perspective(self, event=event)
1280
1281	# If we have more than one displayed plot, make the user choose
1282	if len(panel.plots) >= 1 and \
1283	panel.graph.selected_plottable in panel.plots:
1284	dataset = panel.plots[panel.graph.selected_plottable].name
1285	else:
1286	logging.info("Prview Error: No data is available")
1287	return
1288
1289	# Store a reference to the current plottable
1290	# If we have a suggested value, use it.
1291	try:
1292	estimate = float(self.control_panel.alpha_estimate)
1293	self.control_panel.alpha = estimate
1294	except:
1295	self.control_panel.alpha = self.alpha
1296	logging.info("Prview :Alpha Not estimate yet")
1297	pass
1298	try:
1299	estimate = int(self.control_panel.nterms_estimate)
1300	self.control_panel.nfunc = estimate
1301	except:
1302	self.control_panel.nfunc = self.nfunc
1303	logging.info("Prview : ntemrs Not estimate yet")
1304	pass
1305
1306	self.current_plottable = panel.plots[panel.graph.selected_plottable]
1307	self.set_data([self.current_plottable])
1308	self.control_panel.plotname = dataset
1309	#self.control_panel.nfunc = self.nfunc
1310	self.control_panel.d_max = self.max_length
1311	#self.parent.set_perspective(self.perspective)
1312	self.control_panel._on_invert(None)
1313
1314	def get_panels(self, parent):
1315	"""
1316	Create and return a list of panel objects
1317	"""
1318	from inversion_panel import InversionControl
1319
1320	self.parent = parent
1321	self.control_panel = InversionControl(self.parent, -1,
1322	style=wx.RAISED_BORDER,
1323	standalone=self.standalone)
1324	self.control_panel.set_manager(self)
1325	self.control_panel.nfunc = self.nfunc
1326	self.control_panel.d_max = self.max_length
1327	self.control_panel.alpha = self.alpha
1328	self.perspective = []
1329	self.perspective.append(self.control_panel.window_name)
1330
1331	return [self.control_panel]
1332
1333	def set_data(self, data_list=None):
1334	"""
1335	receive a list of data to compute pr
1336	"""
1337	if data_list is None:
1338	data_list = []
1339	if len(data_list) >= 1:
1340	if len(data_list) == 1:
1341	data = data_list[0]
1342	else:
1343	msg = "Pr panel does not allow multiple Data.\n"
1344	msg += "Please select one!\n"
1345	from pr_widgets import DataDialog
1346	dlg = DataDialog(data_list=data_list, text=msg)
1347	if dlg.ShowModal() == wx.ID_OK:
1348	data = dlg.get_data()
1349	if data is None:
1350	return
1351	if issubclass(data.__class__, Data1D):
1352	try:
1353	wx.PostEvent(self.parent, NewPlotEvent(action='remove',
1354	group_id=GROUP_ID_IQ_DATA,
1355	id=self.data_id))
1356
1357	self.data_id = data.id
1358	self.control_panel._change_file(evt=None, data=data)
1359	except:
1360	msg = "Prview Set_data: " + str(sys.exc_value)
1361	wx.PostEvent(self.parent, StatusEvent(status=msg,
1362	info="error"))
1363	else:
1364	msg = "Pr cannot be computed for data of "
1365	msg += "type %s" % (data_list[0].__class__.__name__)
1366	wx.PostEvent(self.parent,
1367	StatusEvent(status=msg, info='error'))
1368	else:
1369	msg = "Pr contain no data"
1370	wx.PostEvent(self.parent, StatusEvent(status=msg, info='warning'))
1371
1372	def post_init(self):
1373	"""
1374	Post initialization call back to close the loose ends
1375	[Somehow openGL needs this call]
1376	"""
1377	if self.standalone:
1378	self.parent.set_perspective(self.perspective)
1379
1380	if __name__ == "__main__":
1381	i = Plugin()
1382	print i.perform_estimateNT()
1383
1384
1385
1386

Note: See TracBrowser for help on using the repository browser.

SasView

source: sasview/prview/perspectives/pr/pr.py @ 59d542c

Download in other formats: