pr.py @ 14cd91b1

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

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

SasView

source: sasview/prview/perspectives/pr/pr.py @ 14cd91b1

Download in other formats: