pr.py @ 552b3d5

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 552b3d5 was e88ebfd, checked in by Gervaise Alina <gervyh@…>, 14 years ago
working on guiframe
Property mode set to `100644`
File size: 49.5 KB

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

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

SasView

source: sasview/prview/perspectives/pr/pr.py @ 552b3d5

Download in other formats: