pr.py @ 868846b

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 868846b was cb62bd5, checked in by lewis, 7 years ago
Use manually inputted background level in Pr calculation
Property mode set to `100644`
File size: 46.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	from __future__ import print_function
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_value)
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_value)
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_value)
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_value
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_value))
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_value))
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_value))
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_value))
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_value)
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_value)
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_value)
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 @ 868846b

Download in other formats: