pr.py @ 4b73c3e

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

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

SasView

source: sasview/inversionview/src/sans/perspectives/pr/pr.py @ 4b73c3e

Download in other formats: