pr.py @ d16339c

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

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

SasView

source: sasview/src/sas/perspectives/pr/pr.py @ d16339c

Download in other formats: