pr.py @ 62c6dc0

ESS_GUIESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalc

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

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

SasView

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

Download in other formats: