pr.py @ b7ce5ad

magnetic_scatt

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

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

SasView

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

Download in other formats: