pr.py @ 9a5097c

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.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since 9a5097c was 9a5097c, checked in by andyfaff, 7 years ago
MAINT: import numpy as np
Property mode set to `100644`
File size: 47.0 KB

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

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

SasView

source: sasview/src/sas/sasgui/perspectives/pr/pr.py @ 9a5097c

Download in other formats: