pr.py @ db4ed82

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

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

SasView

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

Download in other formats: