pr.py @ 75df58b

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

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

SasView

source: sasview/prview/perspectives/pr/pr.py @ 75df58b

Download in other formats: