pr.py @ df1a6ed

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalcmagnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since df1a6ed was 861f1880, checked in by GitHub <noreply@…>, 8 years ago

Merge pull request #64 from andyfaff/has_key

MAINT: remove has_key occurrences

Property mode set to 100644

File size: 47.0 KB

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

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

SasView

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

Download in other formats: