pr.py @ 8a7d922

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

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

SasView

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

Download in other formats: