pr.py @ baf1ef3

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 baf1ef3 was d80f7e1, checked in by Jae Cho <jhjcho@…>, 13 years ago
updated dataloader calls
Property mode set to `100644`
File size: 50.5 KB

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

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

SasView

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

Download in other formats: