pr.py @ f784b54

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 f784b54 was b1a9e8b, checked in by Gervaise Alina <gervyh@…>, 14 years ago
remove .svs extensions for default reader
Property mode set to `100644`
File size: 49.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.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 DataLoader.loader import Loader
31	import 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	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	"""
834	from copy import deepcopy
835	# Save useful info
836	self.elapsed = elapsed
837	# Keep a copy of the last result
838	self._last_pr = pr.clone()
839	self._last_out = out
840	self._last_cov = cov
841
842	# Save Pr invertor
843	self.pr = pr
844
845	#message = "Computation completed in"
846	#message += %g seconds [chi2=%g]" % (elapsed, pr.chi2)
847	#wx.PostEvent(self.parent, StatusEvent(status=message))
848
849	cov = numpy.ascontiguousarray(cov)
850
851	# Show result on control panel
852	self.control_panel.chi2 = pr.chi2
853	self.control_panel.elapsed = elapsed
854	self.control_panel.oscillation = pr.oscillations(out)
855	#print "OSCILL", pr.oscillations(out)
856	#print "PEAKS:", pr.get_peaks(out)
857	self.control_panel.positive = pr.get_positive(out)
858	self.control_panel.pos_err = pr.get_pos_err(out, cov)
859	self.control_panel.rg = pr.rg(out)
860	self.control_panel.iq0 = pr.iq0(out)
861	self.control_panel.bck = pr.background
862
863	if False:
864	for i in range(len(out)):
865	try:
866	print "%d: %g +- %g" % (i, out[i],
867	math.sqrt(math.fabs(cov[i][i])))
868	except:
869	print sys.exc_value
870	print "%d: %g +- ?" % (i, out[i])
871
872	# Make a plot of I(q) data
873	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
874	new_plot.name = IQ_DATA_LABEL
875	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
876	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
877	if pr.info.has_key("plot_group_id"):
878	new_plot.group_id = pr.info["plot_group_id"]
879	new_plot.id = self.data_id
880	self.parent.update_theory(data_id=self.data_id,
881	theory=new_plot)
882	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
883
884	# Show I(q) fit
885	self.show_iq(out, self.pr)
886
887	# Show P(r) fit
888	x_values, x_range = self.show_pr(out, self.pr, cov)
889
890	# Popup result panel
891	#result_panel = InversionResults(self.parent,
892	#-1, style=wx.RAISED_BORDER)
893
894	def show_data(self, path=None, data=None, reset=False):
895	"""
896	Show data read from a file
897
898	:param path: file path
899	:param reset: if True all other plottables will be cleared
900
901	"""
902	#if path is not None:
903	if data is not None:
904	try:
905	pr = self._create_file_pr(data)
906	except:
907	status = "Problem reading data: %s" % sys.exc_value
908	wx.PostEvent(self.parent, StatusEvent(status=status))
909	raise RuntimeError, status
910
911	# If the file contains nothing, just return
912	if pr is None:
913	raise RuntimeError, "Loaded data is invalid"
914
915	self.pr = pr
916
917	# Make a plot of I(q) data
918	if self.pr.err == None:
919	new_plot = Data1D(self.pr.x, self.pr.y)
920	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
921	else:
922	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
923	new_plot.name = IQ_DATA_LABEL
924	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
925	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
926	new_plot.interactive = True
927	new_plot.group_id = GROUP_ID_IQ_DATA
928	new_plot.id = self.data_id
929	new_plot.title = "I(q)"
930	wx.PostEvent(self.parent,
931	NewPlotEvent(plot=new_plot, title="I(q)", reset=reset))
932
933	self.current_plottable = new_plot
934	# Get Q range
935	self.control_panel.q_min = self.pr.x.min()
936	self.control_panel.q_max = self.pr.x.max()
937
938	def save_data(self, filepath, prstate=None):
939	"""
940	Save data in provided state object.
941
942	:TODO: move the state code away from inversion_panel and move it here.
943	Then remove the "prstate" input and make this method private.
944
945	:param filepath: path of file to write to
946	:param prstate: P(r) inversion state
947
948	"""
949	#TODO: do we need this or can we use DataLoader.loader.save directly?
950
951	# Add output data and coefficients to state
952	prstate.coefficients = self._last_out
953	prstate.covariance = self._last_cov
954
955	# Write the output to file
956	# First, check that the data is of the right type
957	if issubclass(self.current_plottable.__class__,
958	DataLoader.data_info.Data1D):
959	self.state_reader.write(filepath, self.current_plottable, prstate)
960	else:
961	msg = "pr.save_data: the data being saved is not a"
962	msg += " DataLoader.data_info.Data1D object"
963	raise RuntimeError, msg
964
965
966	def setup_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None,
967	bck=False, height=0, width=0):
968	"""
969	"""
970	self.alpha = alpha
971	self.nfunc = nfunc
972	self.max_length = d_max
973	self.q_min = q_min
974	self.q_max = q_max
975	self.has_bck = bck
976	self.slit_height = height
977	self.slit_width = width
978
979	try:
980	pr = self._create_plot_pr()
981	if not pr==None:
982	self.pr = pr
983	self.perform_inversion()
984	except:
985	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
986
987	def estimate_plot_inversion(self, alpha, nfunc, d_max,
988	q_min=None, q_max=None,
989	bck=False, height=0, width=0):
990	"""
991	"""
992	self.alpha = alpha
993	self.nfunc = nfunc
994	self.max_length = d_max
995	self.q_min = q_min
996	self.q_max = q_max
997	self.has_bck = bck
998	self.slit_height = height
999	self.slit_width = width
1000
1001	try:
1002	pr = self._create_plot_pr()
1003	if not pr==None:
1004	self.pr = pr
1005	self.perform_estimate()
1006	except:
1007	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1008
1009	def _create_plot_pr(self, estimate=False):
1010	"""
1011	Create and prepare invertor instance from
1012	a plottable data set.
1013
1014	:param path: path of the file to read in
1015
1016	"""
1017	# Sanity check
1018	if self.current_plottable is None:
1019	msg = "Please load a valid data set before proceeding."
1020	wx.PostEvent(self.parent, StatusEvent(status=msg))
1021	return None
1022
1023	# Get the data from the chosen data set and perform inversion
1024	pr = Invertor()
1025	pr.d_max = self.max_length
1026	pr.alpha = self.alpha
1027	pr.q_min = self.q_min
1028	pr.q_max = self.q_max
1029	pr.x = self.current_plottable.x
1030	pr.y = self.current_plottable.y
1031	pr.has_bck = self.has_bck
1032	pr.slit_height = self.slit_height
1033	pr.slit_width = self.slit_width
1034
1035	# Keep track of the plot window title to ensure that
1036	# we can overlay the plots
1037	pr.info["plot_group_id"] = self.current_plottable.group_id
1038
1039	# Fill in errors if none were provided
1040	err = self.current_plottable.dy
1041	all_zeros = True
1042	if err == None:
1043	err = numpy.zeros(len(pr.y))
1044	else:
1045	for i in range(len(err)):
1046	if err[i]>0:
1047	all_zeros = False
1048
1049	if all_zeros:
1050	scale = None
1051	min_err = 0.0
1052	for i in range(len(pr.y)):
1053	# Scale the error so that we can fit over several decades of Q
1054	if scale==None:
1055	scale = 0.05*math.sqrt(pr.y[i])
1056	min_err = 0.01*pr.y[i]
1057	err[i] = scale*math.sqrt( math.fabs(pr.y[i]) ) + min_err
1058	message = "The loaded file had no error bars, "
1059	message += "statistical errors are assumed."
1060	wx.PostEvent(self.parent, StatusEvent(status=message))
1061
1062	pr.err = err
1063
1064	return pr
1065
1066
1067	def setup_file_inversion(self, alpha, nfunc, d_max, data,
1068	path=None, q_min=None, q_max=None,
1069	bck=False, height=0, width=0):
1070	"""
1071	"""
1072	self.alpha = alpha
1073	self.nfunc = nfunc
1074	self.max_length = d_max
1075	self.q_min = q_min
1076	self.q_max = q_max
1077	self.has_bck = bck
1078	self.slit_height = height
1079	self.slit_width = width
1080
1081	try:
1082	#pr = self._create_file_pr(path)
1083	pr = self._create_file_pr(data)
1084	if not pr==None:
1085	self.pr = pr
1086	self.perform_inversion()
1087	except:
1088	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1089
1090	def estimate_file_inversion(self, alpha, nfunc, d_max, data,
1091	path=None, q_min=None, q_max=None,
1092	bck=False, height=0, width=0):
1093	"""
1094	"""
1095	self.alpha = alpha
1096	self.nfunc = nfunc
1097	self.max_length = d_max
1098	self.q_min = q_min
1099	self.q_max = q_max
1100	self.has_bck = bck
1101	self.slit_height = height
1102	self.slit_width = width
1103
1104	try:
1105	pr = self._create_file_pr(data)
1106	#pr = self._create_file_pr(path)
1107	if not pr is None:
1108	self.pr = pr
1109	self.perform_estimate()
1110	except:
1111	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1112
1113
1114	def _create_file_pr(self, data):
1115	"""
1116	Create and prepare invertor instance from
1117	a file data set.
1118
1119	:param path: path of the file to read in
1120
1121	"""
1122	# Load data
1123	#if os.path.isfile(path):
1124	"""
1125	if self._current_file_data is not None \
1126	and self._current_file_data.path==path:
1127	# Protect against corrupted data from
1128	# previous failed load attempt
1129	if self._current_file_data.x is None:
1130	return None
1131	x = self._current_file_data.x
1132	y = self._current_file_data.y
1133	err = self._current_file_data.err
1134
1135	message = "The data from this file has already been loaded."
1136	wx.PostEvent(self.parent, StatusEvent(status=message))
1137	else:
1138	"""
1139	# Reset the status bar so that we don't get mixed up
1140	# with old messages.
1141	#TODO: refactor this into a proper status handling
1142	wx.PostEvent(self.parent, StatusEvent(status=''))
1143	try:
1144	class FileData:
1145	x = None
1146	y = None
1147	err = None
1148	path = None
1149	def __init__(self, path):
1150	self.path = path
1151
1152	self._current_file_data = FileData(data.path)
1153	self._current_file_data.x = data.x
1154	self._current_file_data.y = data.y
1155	self._current_file_data.err = data.dy
1156	x, y, err = data.x, data.y, data.dy
1157	except:
1158	load_error(sys.exc_value)
1159	return None
1160
1161	# If the file contains no data, just return
1162	if x is None or len(x) == 0:
1163	load_error("The loaded file contains no data")
1164	return None
1165
1166	# If we have not errors, add statistical errors
1167	if y is not None:
1168	if err == None or numpy.all(err) == 0:
1169	err = numpy.zeros(len(y))
1170	scale = None
1171	min_err = 0.0
1172	for i in range(len(y)):
1173	# Scale the error so that we can fit over several decades of Q
1174	if scale == None:
1175	scale = 0.05 * math.sqrt(y[i])
1176	min_err = 0.01 * y[i]
1177	err[i] = scale * math.sqrt(math.fabs(y[i])) + min_err
1178	message = "The loaded file had no error bars, "
1179	message += "statistical errors are assumed."
1180	wx.PostEvent(self.parent, StatusEvent(status=message))
1181
1182	try:
1183	# Get the data from the chosen data set and perform inversion
1184	pr = Invertor()
1185	pr.d_max = self.max_length
1186	pr.alpha = self.alpha
1187	pr.q_min = self.q_min
1188	pr.q_max = self.q_max
1189	pr.x = x
1190	pr.y = y
1191	pr.err = err
1192	pr.has_bck = self.has_bck
1193	pr.slit_height = self.slit_height
1194	pr.slit_width = self.slit_width
1195	return pr
1196	except:
1197	load_error(sys.exc_value)
1198	return None
1199
1200	def perform_estimate(self):
1201	"""
1202	"""
1203	from pr_thread import EstimatePr
1204	from copy import deepcopy
1205
1206	# If a thread is already started, stop it
1207	if self.estimation_thread != None and \
1208	self.estimation_thread.isrunning():
1209	self.estimation_thread.stop()
1210
1211	pr = self.pr.clone()
1212	self.estimation_thread = EstimatePr(pr, self.nfunc,
1213	error_func=self._thread_error,
1214	completefn = self._estimate_completed,
1215	updatefn = None)
1216	self.estimation_thread.queue()
1217	self.estimation_thread.ready(2.5)
1218
1219	def perform_estimateNT(self):
1220	"""
1221	"""
1222	from pr_thread import EstimateNT
1223	from copy import deepcopy
1224
1225	# If a thread is already started, stop it
1226	if self.estimation_thread != None and self.estimation_thread.isrunning():
1227	self.estimation_thread.stop()
1228
1229	pr = self.pr.clone()
1230	# Skip the slit settings for the estimation
1231	# It slows down the application and it doesn't change the estimates
1232	pr.slit_height = 0.0
1233	pr.slit_width = 0.0
1234	self.estimation_thread = EstimateNT(pr, self.nfunc,
1235	error_func=self._thread_error,
1236	completefn = self._estimateNT_completed,
1237	updatefn = None)
1238	self.estimation_thread.queue()
1239	self.estimation_thread.ready(2.5)
1240
1241	def perform_inversion(self):
1242	"""
1243	"""
1244	# Time estimate
1245	#estimated = self.elapsedself.nfunc*2
1246	#message = "Computation time may take up to %g seconds" % self.elapsed
1247	#wx.PostEvent(self.parent, StatusEvent(status=message))
1248
1249	# Start inversion thread
1250	self.start_thread()
1251	return
1252
1253	out, cov = self.pr.lstsq(self.nfunc)
1254
1255	# Save useful info
1256	self.elapsed = self.pr.elapsed
1257
1258	for i in range(len(out)):
1259	try:
1260	print "%d: %g +- %g" % (i, out[i],
1261	math.sqrt(math.fabs(cov[i][i])))
1262	except:
1263	print "%d: %g +- ?" % (i, out[i])
1264
1265	# Make a plot of I(q) data
1266	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
1267	new_plot.name = "I_{obs}(q)"
1268	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
1269	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
1270	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
1271	# Show I(q) fit
1272	self.show_iq(out, self.pr)
1273	# Show P(r) fit
1274	x_values, x_range = self.show_pr(out, self.pr, cov=cov)
1275
1276	def _on_context_inversion(self, event):
1277	"""
1278	"""
1279	panel = event.GetEventObject()
1280	Plugin.on_perspective(self, event=event)
1281
1282	# If we have more than one displayed plot, make the user choose
1283	if len(panel.plots) >= 1 and \
1284	panel.graph.selected_plottable in panel.plots:
1285	dataset = panel.plots[panel.graph.selected_plottable].name
1286	else:
1287	logging.info("Prview Error: No data is available")
1288	return
1289
1290	# Store a reference to the current plottable
1291	# If we have a suggested value, use it.
1292	try:
1293	estimate = float(self.control_panel.alpha_estimate)
1294	self.control_panel.alpha = estimate
1295	except:
1296	self.control_panel.alpha = self.alpha
1297	logging.info("Prview :Alpha Not estimate yet")
1298	pass
1299	try:
1300	estimate = int(self.control_panel.nterms_estimate)
1301	self.control_panel.nfunc = estimate
1302	except:
1303	self.control_panel.nfunc = self.nfunc
1304	logging.info("Prview : ntemrs Not estimate yet")
1305	pass
1306
1307	self.current_plottable = panel.plots[panel.graph.selected_plottable]
1308	self.set_data([self.current_plottable])
1309	self.control_panel.plotname = dataset
1310	#self.control_panel.nfunc = self.nfunc
1311	self.control_panel.d_max = self.max_length
1312	#self.parent.set_perspective(self.perspective)
1313	self.control_panel._on_invert(None)
1314
1315	def get_panels(self, parent):
1316	"""
1317	Create and return a list of panel objects
1318	"""
1319	from inversion_panel import InversionControl
1320
1321	self.parent = parent
1322	self.control_panel = InversionControl(self.parent, -1,
1323	style=wx.RAISED_BORDER,
1324	standalone=self.standalone)
1325	self.control_panel.set_manager(self)
1326	self.control_panel.nfunc = self.nfunc
1327	self.control_panel.d_max = self.max_length
1328	self.control_panel.alpha = self.alpha
1329	self.perspective = []
1330	self.perspective.append(self.control_panel.window_name)
1331
1332	return [self.control_panel]
1333
1334	def set_data(self, data_list=None):
1335	"""
1336	receive a list of data to compute pr
1337	"""
1338	if data_list is None:
1339	data_list = []
1340	if len(data_list) >= 1:
1341	if len(data_list) == 1:
1342	data = data_list[0]
1343	else:
1344	msg = "Pr panel does not allow multiple Data.\n"
1345	msg += "Please select one!\n"
1346	from pr_widgets import DataDialog
1347	dlg = DataDialog(data_list=data_list, text=msg)
1348	if dlg.ShowModal() == wx.ID_OK:
1349	data = dlg.get_data()
1350	if data is None:
1351	return
1352	if issubclass(data.__class__, Data1D):
1353	try:
1354	wx.PostEvent(self.parent, NewPlotEvent(action='remove',
1355	group_id=GROUP_ID_IQ_DATA,
1356	id=self.data_id))
1357
1358	self.data_id = data.id
1359	self.control_panel._change_file(evt=None, data=data)
1360	except:
1361	msg = "Prview Set_data: " + str(sys.exc_value)
1362	wx.PostEvent(self.parent, StatusEvent(status=msg,
1363	info="error"))
1364	else:
1365	msg = "Pr cannot be computed for data of "
1366	msg += "type %s" % (data_list[0].__class__.__name__)
1367	wx.PostEvent(self.parent,
1368	StatusEvent(status=msg, info='error'))
1369	else:
1370	msg = "Pr contain no data"
1371	wx.PostEvent(self.parent, StatusEvent(status=msg, info='warning'))
1372
1373	def post_init(self):
1374	"""
1375	Post initialization call back to close the loose ends
1376	[Somehow openGL needs this call]
1377	"""
1378	if self.standalone:
1379	self.parent.set_perspective(self.perspective)
1380
1381	if __name__ == "__main__":
1382	i = Plugin()
1383	print i.perform_estimateNT()
1384
1385
1386
1387

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

SasView

source: sasview/prview/perspectives/pr/pr.py @ f784b54

Download in other formats: