pr.py @ 9e0dd49

magnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249unittest-saveload

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

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

SasView

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

Download in other formats: