pr.py @ 91d910d

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

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

SasView

source: sasview/prview/perspectives/pr/pr.py @ 91d910d

Download in other formats: