pr.py @ 66ff250

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

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

SasView

source: sasview/prview/perspectives/pr/pr.py @ 66ff250

Download in other formats: