pr.py @ 6675651

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.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since 6675651 was 5777106, checked in by Mathieu Doucet <doucetm@…>, 11 years ago
Moving things around. Will definitely not build.
Property mode set to `100644`
File size: 48.8 KB

Line
1
2	################################################################################
3	#This software was developed by the University of Tennessee as part of the
4	#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
5	#project funded by the US National Science Foundation.
6	#
7	#See the license text in license.txt
8	#
9	#copyright 2009, University of Tennessee
10	################################################################################
11
12
13	# Make sure the option of saving each curve is available
14	# Use the I(q) curve as input and compare the output to P(r)
15
16	import os
17	import sys
18	import wx
19	import logging
20	import time
21	import math
22	import numpy
23	import pylab
24	from sans.guiframe.gui_manager import MDIFrame
25	from sans.guiframe.dataFitting import Data1D
26	from sans.guiframe.events import NewPlotEvent
27	from sans.guiframe.events import StatusEvent
28	from sans.guiframe.gui_style import GUIFRAME_ID
29	from sans.pr.invertor import Invertor
30	from sans.dataloader.loader import Loader
31	import sans.dataloader
32
33	from pr_widgets import load_error
34	from sans.guiframe.plugin_base import PluginBase
35
36
37	PR_FIT_LABEL = r"$P_{fit}(r)$"
38	PR_LOADED_LABEL = r"$P_{loaded}(r)$"
39	IQ_DATA_LABEL = r"$I_{obs}(q)$"
40	IQ_FIT_LABEL = r"$I_{fit}(q)$"
41	IQ_SMEARED_LABEL = r"$I_{smeared}(q)$"
42	GROUP_ID_IQ_DATA = r"$I_{obs}(q)$"
43	GROUP_ID_PR_FIT = r"$P_{fit}(r)$"
44
45
46
47	class Plugin(PluginBase):
48	"""
49	P(r) inversion perspective
50	"""
51	DEFAULT_ALPHA = 0.0001
52	DEFAULT_NFUNC = 10
53	DEFAULT_DMAX = 140.0
54
55	def __init__(self, standalone=True):
56	PluginBase.__init__(self, name="Pr Inversion", standalone=standalone)
57	## Simulation window manager
58	self.simview = None
59
60	## State data
61	self.alpha = self.DEFAULT_ALPHA
62	self.nfunc = self.DEFAULT_NFUNC
63	self.max_length = self.DEFAULT_DMAX
64	self.q_min = None
65	self.q_max = None
66	self.has_bck = False
67	self.slit_height = 0
68	self.slit_width = 0
69	## Remember last plottable processed
70	self.last_data = "sphere_60_q0_2.txt"
71	self._current_file_data = None
72	## Time elapsed for last computation [sec]
73	# Start with a good default
74	self.elapsed = 0.022
75	self.iq_data_shown = False
76
77	## Current invertor
78	self.invertor = None
79	self.pr = None
80	self.data_id = IQ_DATA_LABEL
81	# Copy of the last result in case we need to display it.
82	self._last_pr = None
83	self._last_out = None
84	self._last_cov = None
85	## Calculation thread
86	self.calc_thread = None
87	## Estimation thread
88	self.estimation_thread = None
89	## Result panel
90	self.control_panel = None
91	## Currently views plottable
92	self.current_plottable = None
93	## Number of P(r) points to display on the output plot
94	self._pr_npts = 51
95	## Flag to let the plug-in know that it is running standalone
96	self.standalone = standalone
97	self._normalize_output = False
98	self._scale_output_unity = False
99
100	## List of added P(r) plots
101	self._added_plots = {}
102	self._default_Iq = {}
103	self.list_plot_id = []
104
105	# Associate the inversion state reader with .prv files
106	from inversion_state import Reader
107
108	# Create a CanSAS/Pr reader
109	self.state_reader = Reader(self.set_state)
110	self._extensions = '.prv'
111	l = Loader()
112	l.associate_file_reader('.prv', self.state_reader)
113	#l.associate_file_reader(".svs", self.state_reader)
114
115	# Log startup
116	logging.info("Pr(r) plug-in started")
117
118	def delete_data(self, data_id):
119	"""
120	delete the data association with prview
121	"""
122	self.control_panel.clear_panel()
123
124	def get_data(self):
125	"""
126	Returns the current data
127	"""
128	return self.current_plottable
129
130	def set_state(self, state=None, datainfo=None):
131	"""
132	Call-back method for the inversion state reader.
133	This method is called when a .prv file is loaded.
134
135	:param state: InversionState object
136	:param datainfo: Data1D object [optional]
137
138	"""
139	try:
140	if datainfo.__class__.__name__ == 'list':
141	if len(datainfo) >= 1:
142	data = datainfo[0]
143	else:
144	data = None
145	else:
146	data = datainfo
147	if data is None:
148	msg = "Pr.set_state: datainfo parameter cannot "
149	msg += "be None in standalone mode"
150	raise RuntimeError, msg
151
152	# Ensuring that plots are coordinated correctly
153	t = time.localtime(data.meta_data['prstate'].timestamp)
154	time_str = time.strftime("%b %d %H:%M", t)
155
156	# Check that no time stamp is already appended
157	max_char = data.meta_data['prstate'].file.find("[")
158	if max_char < 0:
159	max_char = len(data.meta_data['prstate'].file)
160
161	datainfo.meta_data['prstate'].file =\
162	data.meta_data['prstate'].file[0:max_char]\
163	+ ' [' + time_str + ']'
164
165	data.filename = data.meta_data['prstate'].file
166	# TODO:
167	#remove this call when state save all information about the gui data
168	# such as ID , Group_ID, etc...
169	#make self.current_plottable = datainfo directly
170	self.current_plottable = self.parent.create_gui_data(data, None)
171	self.current_plottable.group_id = data.meta_data['prstate'].file
172
173	# Make sure the user sees the P(r) panel after loading
174	#self.parent.set_perspective(self.perspective)
175	self.on_perspective(event=None)
176	# Load the P(r) results
177	#state = self.state_reader.get_state()
178	data_dict = {self.current_plottable.id:self.current_plottable}
179	self.parent.add_data(data_list=data_dict)
180	wx.PostEvent(self.parent, NewPlotEvent(plot=self.current_plottable,
181	title=self.current_plottable.title))
182	self.control_panel.set_state(state)
183	except:
184	logging.error("prview.set_state: %s" % sys.exc_value)
185
186
187	def help(self, evt):
188	"""
189	Show a general help dialog.
190
191	:TODO: replace the text with a nice image
192
193	"""
194	from inversion_panel import HelpDialog
195	dialog = HelpDialog(None, -1)
196	if dialog.ShowModal() == wx.ID_OK:
197	dialog.Destroy()
198	else:
199	dialog.Destroy()
200
201	def _fit_pr(self, evt):
202	"""
203	"""
204	# Generate P(r) for sphere
205	radius = 60.0
206	d_max = 2*radius
207
208	r = pylab.arange(0.01, d_max, d_max/51.0)
209	M = len(r)
210	y = numpy.zeros(M)
211	pr_err = numpy.zeros(M)
212
213	sum = 0.0
214	for j in range(M):
215	value = self.pr_theory(r[j], radius)
216	sum += value
217	y[j] = value
218	pr_err[j] = math.sqrt(y[j])
219
220	y = y/sum*d_max/len(r)
221
222	# Perform fit
223	pr = Invertor()
224	pr.d_max = d_max
225	pr.alpha = 0
226	pr.x = r
227	pr.y = y
228	pr.err = pr_err
229	out, cov = pr.pr_fit()
230	for i in range(len(out)):
231	print "%g +- %g" % (out[i], math.sqrt(cov[i][i]))
232
233	# Show input P(r)
234	title = "Pr"
235	new_plot = Data1D(pr.x, pr.y, dy=pr.err)
236	new_plot.name = "P_{obs}(r)"
237	new_plot.xaxis("\\rm{r}", 'A')
238	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
239	new_plot.group_id = "P_{obs}(r)"
240	new_plot.id = "P_{obs}(r)"
241	new_plot.title = title
242	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
243	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
244
245	# Show P(r) fit
246	self.show_pr(out, pr)
247
248	# Show I(q) fit
249	q = pylab.arange(0.001, 0.1, 0.01/51.0)
250	self.show_iq(out, pr, q)
251
252	def show_shpere(self, x, radius=70.0, x_range=70.0):
253	"""
254	"""
255	# Show P(r)
256	y_true = numpy.zeros(len(x))
257
258	sum_true = 0.0
259	for i in range(len(x)):
260	y_true[i] = self.pr_theory(x[i], radius)
261	sum_true += y_true[i]
262
263	y_true = y_true/sum_true*x_range/len(x)
264
265	# Show the theory P(r)
266	new_plot = Data1D(x, y_true)
267	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
268	new_plot.name = "P_{true}(r)"
269	new_plot.xaxis("\\rm{r}", 'A')
270	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
271	new_plot.id = "P_{true}(r)"
272	new_plot.group_id = "P_{true}(r)"
273	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
274	#Put this call in plottables/guitools
275	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot,
276	title="Sphere P(r)"))
277
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	"""
291	qtemp = pr.x
292	if not q == None:
293	qtemp = q
294
295	# Make a plot
296	maxq = -1
297	for q_i in qtemp:
298	if q_i>maxq:
299	maxq=q_i
300
301	minq = 0.001
302
303	# Check for user min/max
304	if not pr.q_min == None:
305	minq = pr.q_min
306	if not pr.q_max == None:
307	maxq = pr.q_max
308
309	x = pylab.arange(minq, maxq, maxq/301.0)
310	y = numpy.zeros(len(x))
311	err = numpy.zeros(len(x))
312	for i in range(len(x)):
313	value = pr.iq(out, x[i])
314	y[i] = value
315	try:
316	err[i] = math.sqrt(math.fabs(value))
317	except:
318	err[i] = 1.0
319	print "Error getting error", value, x[i]
320
321	new_plot = Data1D(x, y)
322	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
323	new_plot.name = IQ_FIT_LABEL
324	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
325	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
326	title = "I(q)"
327	new_plot.title = title
328
329	# If we have a group ID, use it
330	if pr.info.has_key("plot_group_id"):
331	new_plot.group_id = pr.info["plot_group_id"]
332	new_plot.id = IQ_FIT_LABEL
333	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
334
335	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
336
337	# If we have used slit smearing, plot the smeared I(q) too
338	if pr.slit_width > 0 or pr.slit_height > 0:
339	x = pylab.arange(minq, maxq, maxq/301.0)
340	y = numpy.zeros(len(x))
341	err = numpy.zeros(len(x))
342	for i in range(len(x)):
343	value = pr.iq_smeared(out, x[i])
344	y[i] = value
345	try:
346	err[i] = math.sqrt(math.fabs(value))
347	except:
348	err[i] = 1.0
349	print "Error getting error", value, x[i]
350
351	new_plot = Data1D(x, y)
352	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
353	new_plot.name = IQ_SMEARED_LABEL
354	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
355	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
356	# If we have a group ID, use it
357	if pr.info.has_key("plot_group_id"):
358	new_plot.group_id = pr.info["plot_group_id"]
359	new_plot.id = IQ_SMEARED_LABEL
360	new_plot.title = title
361	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
362	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
363
364	def _on_pr_npts(self, evt):
365	"""
366	Redisplay P(r) with a different number of points
367	"""
368	from inversion_panel import PrDistDialog
369	dialog = PrDistDialog(None, -1)
370	dialog.set_content(self._pr_npts)
371	if dialog.ShowModal() == wx.ID_OK:
372	self._pr_npts= dialog.get_content()
373	dialog.Destroy()
374	self.show_pr(self._last_out, self._last_pr, self._last_cov)
375	else:
376	dialog.Destroy()
377
378
379	def show_pr(self, out, pr, cov=None):
380	"""
381	"""
382	# Show P(r)
383	x = pylab.arange(0.0, pr.d_max, pr.d_max/self._pr_npts)
384
385	y = numpy.zeros(len(x))
386	dy = numpy.zeros(len(x))
387	y_true = numpy.zeros(len(x))
388
389	sum = 0.0
390	pmax = 0.0
391	cov2 = numpy.ascontiguousarray(cov)
392
393	for i in range(len(x)):
394	if cov2==None:
395	value = pr.pr(out, x[i])
396	else:
397	(value, dy[i]) = pr.pr_err(out, cov2, x[i])
398	sum += value*pr.d_max/len(x)
399
400	# keep track of the maximum P(r) value
401	if value>pmax:
402	pmax = value
403
404	y[i] = value
405
406	if self._normalize_output==True:
407	y = y/sum
408	dy = dy/sum
409	elif self._scale_output_unity==True:
410	y = y/pmax
411	dy = dy/pmax
412
413	if cov2==None:
414	new_plot = Data1D(x, y)
415	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
416	else:
417	new_plot = Data1D(x, y, dy=dy)
418	new_plot.name = PR_FIT_LABEL
419	new_plot.xaxis("\\rm{r}", 'A')
420	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
421	new_plot.title = "P(r) fit"
422	new_plot.id = PR_FIT_LABEL
423	# Make sure that the plot is linear
424	new_plot.xtransform = "x"
425	new_plot.ytransform = "y"
426	new_plot.group_id = GROUP_ID_PR_FIT
427	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
428	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="P(r) fit"))
429	return x, pr.d_max
430
431	def load(self, data):
432	"""
433	Load data. This will eventually be replaced
434	by our standard DataLoader class.
435	"""
436	class FileData:
437	x = None
438	y = None
439	err = None
440	path = None
441
442	def __init__(self, path):
443	self.path = path
444
445	self._current_file_data = FileData(data.path)
446
447	# Use data loader to load file
448	#dataread = Loader().load(path)
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 = numpy.zeros(0)
485	data_y = numpy.zeros(0)
486	data_err = numpy.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 = numpy.append(data_x, x)
509	data_y = numpy.append(data_y, y)
510	data_err = numpy.append(data_err, err)
511	except:
512	pass
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 = numpy.zeros(0)
533	data_y = numpy.zeros(0)
534	data_err = numpy.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 = numpy.append(data_x, x)
560	data_y = numpy.append(data_y, y)
561	data_err = numpy.append(data_err, err)
562	except:
563	pass
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.5r/R)2) ((1.0-0.5r/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,
623	IQ_SMEARED_LABEL]:
624	return []
625	elif item.id == graph.selected_plottable:
626	if not self.standalone and issubclass(item.__class__, Data1D):
627	return [["Compute P(r)",
628	"Compute P(r) from distribution",
629	self._on_context_inversion]]
630
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 = numpy.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	sum = numpy.sum(self._added_plots[plot].y)
671	npts = len(self._added_plots[plot].x)
672	sum *= self._added_plots[plot].x[npts-1]/npts
673	y = self._added_plots[plot].y/sum
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	"""
721	from pr_thread import CalcPr
722
723	# If a thread is already started, stop it
724	if self.calc_thread != None and self.calc_thread.isrunning():
725	self.calc_thread.stop()
726
727	pr = self.pr.clone()
728	self.calc_thread = CalcPr(pr, self.nfunc,
729	error_func=self._thread_error,
730	completefn=self._completed, updatefn=None)
731	self.calc_thread.queue()
732	self.calc_thread.ready(2.5)
733
734	def _thread_error(self, error):
735	"""
736	"""
737	wx.PostEvent(self.parent, StatusEvent(status=error))
738
739	def _estimate_completed(self, alpha, message, elapsed):
740	"""
741	Parameter estimation completed,
742	display the results to the user
743
744	:param alpha: estimated best alpha
745	:param elapsed: computation time
746
747	"""
748	# Save useful info
749	self.elapsed = elapsed
750	self.control_panel.alpha_estimate = alpha
751	if not message==None:
752	wx.PostEvent(self.parent, StatusEvent(status=str(message)))
753	self.perform_estimateNT()
754
755	def _estimateNT_completed(self, nterms, alpha, message, elapsed):
756	"""
757	Parameter estimation completed,
758	display the results to the user
759
760	:param alpha: estimated best alpha
761	:param nterms: estimated number of terms
762	:param elapsed: computation time
763
764	"""
765	# Save useful info
766	self.elapsed = elapsed
767	self.control_panel.nterms_estimate = nterms
768	self.control_panel.alpha_estimate = alpha
769	if not message==None:
770	wx.PostEvent(self.parent, StatusEvent(status=str(message)))
771
772	def _completed(self, out, cov, pr, elapsed):
773	"""
774	wxCallAfter Method called with the results when the inversion
775	is done
776
777	:param out: output coefficient for the base functions
778	:param cov: covariance matrix
779	:param pr: Invertor instance
780	:param elapsed: time spent computing
781	"""
782	# Ensure hat you have all inputs are ready at the time call happens:
783	# Without CallAfter, it will freeze with wx >= 2.9.
784	wx.CallAfter(self._completed_call, out, cov, pr, elapsed)
785
786	def _completed_call(self, out, cov, pr, elapsed):
787	"""
788	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	"""
797	# Save useful info
798	self.elapsed = elapsed
799	# Keep a copy of the last result
800	self._last_pr = pr.clone()
801	self._last_out = out
802	self._last_cov = cov
803
804	# Save Pr invertor
805	self.pr = pr
806
807	#message = "Computation completed in"
808	#message += %g seconds [chi2=%g]" % (elapsed, pr.chi2)
809	#wx.PostEvent(self.parent, StatusEvent(status=message))
810
811	cov = numpy.ascontiguousarray(cov)
812
813	# Show result on control panel
814	self.control_panel.chi2 = pr.chi2
815	self.control_panel.elapsed = elapsed
816	self.control_panel.oscillation = pr.oscillations(out)
817	self.control_panel.positive = pr.get_positive(out)
818	self.control_panel.pos_err = pr.get_pos_err(out, cov)
819	self.control_panel.rg = pr.rg(out)
820	self.control_panel.iq0 = pr.iq0(out)
821	self.control_panel.bck = pr.background
822
823	# Show I(q) fit
824	self.show_iq(out, self.pr)
825
826	# Show P(r) fit
827	x_values, x_range = self.show_pr(out, self.pr, cov)
828
829	def show_data(self, path=None, data=None, reset=False):
830	"""
831	Show data read from a file
832
833	:param path: file path
834	:param reset: if True all other plottables will be cleared
835
836	"""
837	#if path is not None:
838	if data is not None:
839	try:
840	pr = self._create_file_pr(data)
841	except:
842	status = "Problem reading data: %s" % sys.exc_value
843	wx.PostEvent(self.parent, StatusEvent(status=status))
844	raise RuntimeError, status
845
846	# If the file contains nothing, just return
847	if pr is None:
848	raise RuntimeError, "Loaded data is invalid"
849
850	self.pr = pr
851
852	# Make a plot of I(q) data
853	if self.pr.err == None:
854	new_plot = Data1D(self.pr.x, self.pr.y)
855	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
856	else:
857	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
858	new_plot.name = IQ_DATA_LABEL
859	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
860	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
861	new_plot.interactive = True
862	new_plot.group_id = GROUP_ID_IQ_DATA
863	new_plot.id = self.data_id
864	new_plot.title = "I(q)"
865	wx.PostEvent(self.parent,
866	NewPlotEvent(plot=new_plot, title="I(q)", reset=reset))
867
868	self.current_plottable = new_plot
869	# Get Q range
870	self.control_panel.q_min = min(self.pr.x)
871	self.control_panel.q_max = max(self.pr.x)
872
873	def save_data(self, filepath, prstate=None):
874	"""
875	Save data in provided state object.
876
877	:TODO: move the state code away from inversion_panel and move it here.
878	Then remove the "prstate" input and make this method private.
879
880	:param filepath: path of file to write to
881	:param prstate: P(r) inversion state
882
883	"""
884	#TODO: do we need this or can we use DataLoader.loader.save directly?
885
886	# Add output data and coefficients to state
887	prstate.coefficients = self._last_out
888	prstate.covariance = self._last_cov
889
890	# Write the output to file
891	# First, check that the data is of the right type
892	if issubclass(self.current_plottable.__class__,
893	sans.dataloader.data_info.Data1D):
894	self.state_reader.write(filepath, self.current_plottable, prstate)
895	else:
896	msg = "pr.save_data: the data being saved is not a"
897	msg += " sans.dataloader.data_info.Data1D object"
898	raise RuntimeError, msg
899
900
901	def setup_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None,
902	bck=False, height=0, width=0):
903	"""
904	"""
905	self.alpha = alpha
906	self.nfunc = nfunc
907	self.max_length = d_max
908	self.q_min = q_min
909	self.q_max = q_max
910	self.has_bck = bck
911	self.slit_height = height
912	self.slit_width = width
913
914	try:
915	pr = self._create_plot_pr()
916	if not pr==None:
917	self.pr = pr
918	self.perform_inversion()
919	except:
920	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
921
922	def estimate_plot_inversion(self, alpha, nfunc, d_max,
923	q_min=None, q_max=None,
924	bck=False, height=0, width=0):
925	"""
926	"""
927	self.alpha = alpha
928	self.nfunc = nfunc
929	self.max_length = d_max
930	self.q_min = q_min
931	self.q_max = q_max
932	self.has_bck = bck
933	self.slit_height = height
934	self.slit_width = width
935
936	try:
937	pr = self._create_plot_pr()
938	if not pr==None:
939	self.pr = pr
940	self.perform_estimate()
941	except:
942	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
943
944	def _create_plot_pr(self, estimate=False):
945	"""
946	Create and prepare invertor instance from
947	a plottable data set.
948
949	:param path: path of the file to read in
950
951	"""
952	# Sanity check
953	if self.current_plottable is None:
954	msg = "Please load a valid data set before proceeding."
955	wx.PostEvent(self.parent, StatusEvent(status=msg))
956	return None
957
958	# Get the data from the chosen data set and perform inversion
959	pr = Invertor()
960	pr.d_max = self.max_length
961	pr.alpha = self.alpha
962	pr.q_min = self.q_min
963	pr.q_max = self.q_max
964	pr.x = self.current_plottable.x
965	pr.y = self.current_plottable.y
966	pr.has_bck = self.has_bck
967	pr.slit_height = self.slit_height
968	pr.slit_width = self.slit_width
969
970	# Keep track of the plot window title to ensure that
971	# we can overlay the plots
972	pr.info["plot_group_id"] = self.current_plottable.group_id
973
974	# Fill in errors if none were provided
975	err = self.current_plottable.dy
976	all_zeros = True
977	if err == None:
978	err = numpy.zeros(len(pr.y))
979	else:
980	for i in range(len(err)):
981	if err[i]>0:
982	all_zeros = False
983
984	if all_zeros:
985	scale = None
986	min_err = 0.0
987	for i in range(len(pr.y)):
988	# Scale the error so that we can fit over several decades of Q
989	if scale==None:
990	scale = 0.05*math.sqrt(pr.y[i])
991	min_err = 0.01*pr.y[i]
992	err[i] = scale*math.sqrt( math.fabs(pr.y[i]) ) + min_err
993	message = "The loaded file had no error bars, "
994	message += "statistical errors are assumed."
995	wx.PostEvent(self.parent, StatusEvent(status=message))
996
997	pr.err = err
998
999	return pr
1000
1001
1002	def setup_file_inversion(self, alpha, nfunc, d_max, data,
1003	path=None, q_min=None, q_max=None,
1004	bck=False, height=0, width=0):
1005	"""
1006	"""
1007	self.alpha = alpha
1008	self.nfunc = nfunc
1009	self.max_length = d_max
1010	self.q_min = q_min
1011	self.q_max = q_max
1012	self.has_bck = bck
1013	self.slit_height = height
1014	self.slit_width = width
1015
1016	try:
1017	#pr = self._create_file_pr(path)
1018	pr = self._create_file_pr(data)
1019	if not pr==None:
1020	self.pr = pr
1021	self.perform_inversion()
1022	except:
1023	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1024
1025	def estimate_file_inversion(self, alpha, nfunc, d_max, data,
1026	path=None, q_min=None, q_max=None,
1027	bck=False, height=0, width=0):
1028	"""
1029	"""
1030	self.alpha = alpha
1031	self.nfunc = nfunc
1032	self.max_length = d_max
1033	self.q_min = q_min
1034	self.q_max = q_max
1035	self.has_bck = bck
1036	self.slit_height = height
1037	self.slit_width = width
1038
1039	try:
1040	pr = self._create_file_pr(data)
1041	#pr = self._create_file_pr(path)
1042	if not pr is None:
1043	self.pr = pr
1044	self.perform_estimate()
1045	except:
1046	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1047
1048
1049	def _create_file_pr(self, data):
1050	"""
1051	Create and prepare invertor instance from
1052	a file data set.
1053
1054	:param path: path of the file to read in
1055
1056	"""
1057	# Load data
1058	#if os.path.isfile(path):
1059	"""
1060	if self._current_file_data is not None \
1061	and self._current_file_data.path==path:
1062	# Protect against corrupted data from
1063	# previous failed load attempt
1064	if self._current_file_data.x is None:
1065	return None
1066	x = self._current_file_data.x
1067	y = self._current_file_data.y
1068	err = self._current_file_data.err
1069
1070	message = "The data from this file has already been loaded."
1071	wx.PostEvent(self.parent, StatusEvent(status=message))
1072	else:
1073	"""
1074	# Reset the status bar so that we don't get mixed up
1075	# with old messages.
1076	#TODO: refactor this into a proper status handling
1077	wx.PostEvent(self.parent, StatusEvent(status=''))
1078	try:
1079	class FileData:
1080	x = None
1081	y = None
1082	err = None
1083	path = None
1084	def __init__(self, path):
1085	self.path = path
1086
1087	self._current_file_data = FileData(data.path)
1088	self._current_file_data.x = data.x
1089	self._current_file_data.y = data.y
1090	self._current_file_data.err = data.dy
1091	x, y, err = data.x, data.y, data.dy
1092	except:
1093	load_error(sys.exc_value)
1094	return None
1095
1096	# If the file contains no data, just return
1097	if x is None or len(x) == 0:
1098	load_error("The loaded file contains no data")
1099	return None
1100
1101	# If we have not errors, add statistical errors
1102	if y is not None:
1103	if err == None or numpy.all(err) == 0:
1104	err = numpy.zeros(len(y))
1105	scale = None
1106	min_err = 0.0
1107	for i in range(len(y)):
1108	# Scale the error so that we can fit over several decades of Q
1109	if scale == None:
1110	scale = 0.05 * math.sqrt(y[i])
1111	min_err = 0.01 * y[i]
1112	err[i] = scale * math.sqrt(math.fabs(y[i])) + min_err
1113	message = "The loaded file had no error bars, "
1114	message += "statistical errors are assumed."
1115	wx.PostEvent(self.parent, StatusEvent(status=message))
1116
1117	try:
1118	# Get the data from the chosen data set and perform inversion
1119	pr = Invertor()
1120	pr.d_max = self.max_length
1121	pr.alpha = self.alpha
1122	pr.q_min = self.q_min
1123	pr.q_max = self.q_max
1124	pr.x = x
1125	pr.y = y
1126	pr.err = err
1127	pr.has_bck = self.has_bck
1128	pr.slit_height = self.slit_height
1129	pr.slit_width = self.slit_width
1130	return pr
1131	except:
1132	load_error(sys.exc_value)
1133	return None
1134
1135	def perform_estimate(self):
1136	"""
1137	"""
1138	from pr_thread import EstimatePr
1139	from copy import deepcopy
1140
1141	# If a thread is already started, stop it
1142	if self.estimation_thread != None and \
1143	self.estimation_thread.isrunning():
1144	self.estimation_thread.stop()
1145
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	"""
1157	from pr_thread import EstimateNT
1158	from copy import deepcopy
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
1164	pr = self.pr.clone()
1165	# Skip the slit settings for the estimation
1166	# It slows down the application and it doesn't change the estimates
1167	pr.slit_height = 0.0
1168	pr.slit_width = 0.0
1169	self.estimation_thread = EstimateNT(pr, self.nfunc,
1170	error_func=self._thread_error,
1171	completefn = self._estimateNT_completed,
1172	updatefn = None)
1173	self.estimation_thread.queue()
1174	self.estimation_thread.ready(2.5)
1175
1176	def perform_inversion(self):
1177	"""
1178	"""
1179	# Time estimate
1180	#estimated = self.elapsedself.nfunc*2
1181	#message = "Computation time may take up to %g seconds" % self.elapsed
1182	#wx.PostEvent(self.parent, StatusEvent(status=message))
1183
1184	# Start inversion thread
1185	self.start_thread()
1186	return
1187
1188	out, cov = self.pr.lstsq(self.nfunc)
1189
1190	# Save useful info
1191	self.elapsed = self.pr.elapsed
1192
1193	for i in range(len(out)):
1194	try:
1195	print "%d: %g +- %g" % (i, out[i],
1196	math.sqrt(math.fabs(cov[i][i])))
1197	except:
1198	print "%d: %g +- ?" % (i, out[i])
1199
1200	# Make a plot of I(q) data
1201	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
1202	new_plot.name = "I_{obs}(q)"
1203	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
1204	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
1205	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
1206	# Show I(q) fit
1207	self.show_iq(out, self.pr)
1208	# Show P(r) fit
1209	x_values, x_range = self.show_pr(out, self.pr, cov=cov)
1210
1211	def _on_context_inversion(self, event):
1212	"""
1213	"""
1214	panel = event.GetEventObject()
1215	Plugin.on_perspective(self, event=event)
1216
1217	# If we have more than one displayed plot, make the user choose
1218	if len(panel.plots) >= 1 and \
1219	panel.graph.selected_plottable in panel.plots:
1220	dataset = panel.plots[panel.graph.selected_plottable].name
1221	else:
1222	logging.info("Prview Error: No data is available")
1223	return
1224
1225	# Store a reference to the current plottable
1226	# If we have a suggested value, use it.
1227	try:
1228	estimate = float(self.control_panel.alpha_estimate)
1229	self.control_panel.alpha = estimate
1230	except:
1231	self.control_panel.alpha = self.alpha
1232	logging.info("Prview :Alpha Not estimate yet")
1233	pass
1234	try:
1235	estimate = int(self.control_panel.nterms_estimate)
1236	self.control_panel.nfunc = estimate
1237	except:
1238	self.control_panel.nfunc = self.nfunc
1239	logging.info("Prview : ntemrs Not estimate yet")
1240	pass
1241
1242	self.current_plottable = panel.plots[panel.graph.selected_plottable]
1243	self.set_data([self.current_plottable])
1244	self.control_panel.plotname = dataset
1245	#self.control_panel.nfunc = self.nfunc
1246	self.control_panel.d_max = self.max_length
1247	#self.parent.set_perspective(self.perspective)
1248	self.control_panel._on_invert(None)
1249
1250	def get_panels(self, parent):
1251	"""
1252	Create and return a list of panel objects
1253	"""
1254	from inversion_panel import InversionControl
1255
1256	self.parent = parent
1257	self.frame = MDIFrame(self.parent, None, 'None', (100, 200))
1258	self.control_panel = InversionControl(self.frame, -1,
1259	style=wx.RAISED_BORDER,
1260	standalone=self.standalone)
1261	self.frame.set_panel(self.control_panel)
1262	self._frame_set_helper()
1263	self.control_panel.set_manager(self)
1264	self.control_panel.nfunc = self.nfunc
1265	self.control_panel.d_max = self.max_length
1266	self.control_panel.alpha = self.alpha
1267	self.perspective = []
1268	self.perspective.append(self.control_panel.window_name)
1269
1270	return [self.control_panel]
1271
1272	def set_data(self, data_list=None):
1273	"""
1274	receive a list of data to compute pr
1275	"""
1276	if data_list is None:
1277	data_list = []
1278	if len(data_list) >= 1:
1279	if len(data_list) == 1:
1280	data = data_list[0]
1281	else:
1282	data_1d_list = []
1283	data_2d_list = []
1284	error_msg = ""
1285	# separate data into data1d and data2d list
1286	for data in data_list:
1287	if data is not None:
1288	if issubclass(data.__class__, Data1D):
1289	data_1d_list.append(data)
1290	else:
1291	error_msg += " %s type %s \n" % (str(data.name),
1292	str(data.__class__.__name__))
1293	data_2d_list.append(data)
1294	if len(data_2d_list) > 0:
1295	msg = "PrView does not support the following data types:\n"
1296	msg += error_msg
1297	if len(data_1d_list) == 0:
1298	wx.PostEvent(self.parent,
1299	StatusEvent(status=msg, info='error'))
1300	return
1301	msg += "Prview does not allow multiple data!\n"
1302	msg += "Please select one.\n"
1303	if len(data_list) > 1:
1304	from pr_widgets import DataDialog
1305	dlg = DataDialog(data_list=data_1d_list, text=msg)
1306	if dlg.ShowModal() == wx.ID_OK:
1307	data = dlg.get_data()
1308	else:
1309	data = None
1310	dlg.Destroy()
1311	if data is None:
1312	msg += "PrView receives no data. \n"
1313	wx.PostEvent(self.parent,
1314	StatusEvent(status=msg, info='error'))
1315	return
1316	if issubclass(data.__class__, Data1D):
1317	try:
1318	wx.PostEvent(self.parent, NewPlotEvent(action='remove',
1319	group_id=GROUP_ID_IQ_DATA,
1320	id=self.data_id))
1321
1322	self.data_id = data.id
1323	self.control_panel._change_file(evt=None, data=data)
1324	except:
1325	msg = "Prview Set_data: " + str(sys.exc_value)
1326	wx.PostEvent(self.parent, StatusEvent(status=msg,
1327	info="error"))
1328	else:
1329	msg = "Pr cannot be computed for data of "
1330	msg += "type %s" % (data_list[0].__class__.__name__)
1331	wx.PostEvent(self.parent,
1332	StatusEvent(status=msg, info='error'))
1333	else:
1334	msg = "Pr contain no data"
1335	wx.PostEvent(self.parent, StatusEvent(status=msg, info='warning'))
1336
1337	def post_init(self):
1338	"""
1339	Post initialization call back to close the loose ends
1340	[Somehow openGL needs this call]
1341	"""
1342	if self.standalone:
1343	self.parent.set_perspective(self.perspective)
1344
1345	if __name__ == "__main__":
1346	i = Plugin()
1347	print i.perform_estimateNT()
1348
1349
1350
1351

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

SasView

source: sasview/src/sans/perspectives/pr/pr.py @ 6675651

Download in other formats: