pr.py @ 3cd5806

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

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

SasView

source: sasview/prview/perspectives/pr/pr.py @ 3cd5806

Download in other formats: