pr.py @ 8cd029b

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

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

SasView

source: sasview/prview/perspectives/pr/pr.py @ 8cd029b

Download in other formats: