pr.py @ 380ff18

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 380ff18 was 380ff18, checked in by Jae Cho <jhjcho@…>, 13 years ago
need to call the complete()
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	import DataLoader
32
33	from pr_widgets import load_error
34	from sans.guiframe.plugin_base import PluginBase
35
36
37	PR_FIT_LABEL = r"$P_{fit}(r)$"
38	PR_LOADED_LABEL = r"$P_{loaded}(r)$"
39	IQ_DATA_LABEL = r"$I_{obs}(q)$"
40	IQ_FIT_LABEL = r"$I_{fit}(q)$"
41	IQ_SMEARED_LABEL = r"$I_{smeared}(q)$"
42	GROUP_ID_IQ_DATA = r"$I_{obs}(q)$"
43	GROUP_ID_PR_FIT = r"$P_{fit}(r)$"
44
45
46
47	class Plugin(PluginBase):
48	"""
49	"""
50	DEFAULT_ALPHA = 0.0001
51	DEFAULT_NFUNC = 10
52	DEFAULT_DMAX = 140.0
53
54	def __init__(self, standalone=True):
55	PluginBase.__init__(self, name="Pr Inversion", standalone=standalone)
56	## Simulation window manager
57	self.simview = None
58
59	## State data
60	self.alpha = self.DEFAULT_ALPHA
61	self.nfunc = self.DEFAULT_NFUNC
62	self.max_length = self.DEFAULT_DMAX
63	self.q_min = None
64	self.q_max = None
65	self.has_bck = False
66	self.slit_height = 0
67	self.slit_width = 0
68	## Remember last plottable processed
69	self.last_data = "sphere_60_q0_2.txt"
70	self._current_file_data = None
71	## Time elapsed for last computation [sec]
72	# Start with a good default
73	self.elapsed = 0.022
74	self.iq_data_shown = False
75
76	## Current invertor
77	self.invertor = None
78	self.pr = None
79	self.data_id = IQ_DATA_LABEL
80	# Copy of the last result in case we need to display it.
81	self._last_pr = None
82	self._last_out = None
83	self._last_cov = None
84	## Calculation thread
85	self.calc_thread = None
86	## Estimation thread
87	self.estimation_thread = None
88	## Result panel
89	self.control_panel = None
90	## Currently views plottable
91	self.current_plottable = None
92	## Number of P(r) points to display on the output plot
93	self._pr_npts = 51
94	## Flag to let the plug-in know that it is running standalone
95	self.standalone = standalone
96	self._normalize_output = False
97	self._scale_output_unity = False
98
99	## List of added P(r) plots
100	self._added_plots = {}
101	self._default_Iq = {}
102	self.list_plot_id = []
103
104	# Associate the inversion state reader with .prv files
105	from inversion_state import Reader
106
107	# Create a CanSAS/Pr reader
108	self.state_reader = Reader(self.set_state)
109	self._extensions = '.prv'
110	l = Loader()
111	l.associate_file_reader('.prv', self.state_reader)
112	#l.associate_file_reader(".svs", self.state_reader)
113
114	# Log startup
115	logging.info("Pr(r) plug-in started")
116
117	def delete_data(self, data_id):
118	"""
119	delete the data association with prview
120	"""
121	if self.calc_thread is not None and self.calc_thread.isrunning():
122	msg = "Data in use in Prview\n"
123	if self.estimation_thread is not None and \
124	self.estimation_thread.isrunning():
125	msg = "Data in use in Prview\n"
126	self.control_panel.clear_panel()
127
128	def get_data(self):
129	"""
130	"""
131	return self.current_plottable
132
133	def set_state(self, state=None, datainfo=None):
134	"""
135	Call-back method for the inversion state reader.
136	This method is called when a .prv file is loaded.
137
138	:param state: InversionState object
139	:param datainfo: Data1D object [optional]
140
141	"""
142	try:
143	if datainfo.__class__.__name__ == 'list':
144	if len(datainfo) >= 1:
145	data = datainfo[0]
146	else:
147	data = None
148	else:
149	data = datainfo
150	if data is None:
151	raise RuntimeError, "Pr.set_state: datainfo parameter cannot be None in standalone mode"
152
153	# Ensuring that plots are coordinated correctly
154	t = time.localtime(data.meta_data['prstate'].timestamp)
155	time_str = time.strftime("%b %d %H:%M", t)
156
157	# Check that no time stamp is already appended
158	max_char = data.meta_data['prstate'].file.find("[")
159	if max_char < 0:
160	max_char = len(data.meta_data['prstate'].file)
161
162	datainfo.meta_data['prstate'].file = data.meta_data['prstate'].file[0:max_char] +' [' + time_str + ']'
163	data.filename = data.meta_data['prstate'].file
164	# TODO:
165	#remove this call when state save all information about the gui data
166	# such as ID , Group_ID, etc...
167	#make self.current_plottable = datainfo directly
168	self.current_plottable = self.parent.create_gui_data(data,None)
169	self.current_plottable.group_id = data.meta_data['prstate'].file
170
171	# Make sure the user sees the P(r) panel after loading
172	#self.parent.set_perspective(self.perspective)
173	self.on_perspective(event=None)
174	# Load the P(r) results
175	#state = self.state_reader.get_state()
176	data_dict = {self.current_plottable.id:self.current_plottable}
177	self.parent.add_data(data_list=data_dict)
178	wx.PostEvent(self.parent, NewPlotEvent(plot=self.current_plottable,
179	title=self.current_plottable.title))
180	self.control_panel.set_state(state)
181	except:
182	logging.error("prview.set_state: %s" % sys.exc_value)
183
184
185	def help(self, evt):
186	"""
187	Show a general help dialog.
188
189	:TODO: replace the text with a nice image
190
191	"""
192	from inversion_panel import HelpDialog
193	dialog = HelpDialog(None, -1)
194	if dialog.ShowModal() == wx.ID_OK:
195	dialog.Destroy()
196	else:
197	dialog.Destroy()
198
199	def _fit_pr(self, evt):
200	"""
201	"""
202	from sans.pr.invertor import Invertor
203	# Generate P(r) for sphere
204	radius = 60.0
205	d_max = 2*radius
206
207	r = pylab.arange(0.01, d_max, d_max/51.0)
208	M = len(r)
209	y = numpy.zeros(M)
210	pr_err = numpy.zeros(M)
211
212	sum = 0.0
213	for j in range(M):
214	value = self.pr_theory(r[j], radius)
215	sum += value
216	y[j] = value
217	pr_err[j] = math.sqrt(y[j])
218
219	y = y/sum*d_max/len(r)
220
221	# Perform fit
222	pr = Invertor()
223	pr.d_max = d_max
224	pr.alpha = 0
225	pr.x = r
226	pr.y = y
227	pr.err = pr_err
228	out, cov = pr.pr_fit()
229	for i in range(len(out)):
230	print "%g +- %g" % (out[i], math.sqrt(cov[i][i]))
231
232	# Show input P(r)
233	title = "Pr"
234	new_plot = Data1D(pr.x, pr.y, dy=pr.err)
235	new_plot.name = "P_{obs}(r)"
236	new_plot.xaxis("\\rm{r}", 'A')
237	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
238	new_plot.group_id = "P_{obs}(r)"
239	new_plot.id = "P_{obs}(r)"
240	new_plot.title = title
241	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
242	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
243
244	# Show P(r) fit
245	self.show_pr(out, pr)
246
247	# Show I(q) fit
248	q = pylab.arange(0.001, 0.1, 0.01/51.0)
249	self.show_iq(out, pr, q)
250
251	def show_shpere(self, x, radius=70.0, x_range=70.0):
252	"""
253	"""
254	# Show P(r)
255	y_true = numpy.zeros(len(x))
256
257	sum_true = 0.0
258	for i in range(len(x)):
259	y_true[i] = self.pr_theory(x[i], radius)
260	sum_true += y_true[i]
261
262	y_true = y_true/sum_true*x_range/len(x)
263
264	# Show the theory P(r)
265	new_plot = Data1D(x, y_true)
266	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
267	new_plot.name = "P_{true}(r)"
268	new_plot.xaxis("\\rm{r}", 'A')
269	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
270	new_plot.id = "P_{true}(r)"
271	new_plot.group_id = "P_{true}(r)"
272	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
273	#Put this call in plottables/guitools
274	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot,
275	title="Sphere P(r)"))
276
277
278	def get_npts(self):
279	"""
280	Returns the number of points in the I(q) data
281	"""
282	try:
283	return len(self.pr.x)
284	except:
285	return 0
286
287	def show_iq(self, out, pr, q=None):
288	"""
289	"""
290	qtemp = pr.x
291	if not q==None:
292	qtemp = q
293
294	# Make a plot
295	maxq = -1
296	for q_i in qtemp:
297	if q_i>maxq:
298	maxq=q_i
299
300	minq = 0.001
301
302	# Check for user min/max
303	if not pr.q_min==None:
304	minq = pr.q_min
305	if not pr.q_max==None:
306	maxq = pr.q_max
307
308	x = pylab.arange(minq, maxq, maxq/301.0)
309	y = numpy.zeros(len(x))
310	err = numpy.zeros(len(x))
311	for i in range(len(x)):
312	value = pr.iq(out, x[i])
313	y[i] = value
314	try:
315	err[i] = math.sqrt(math.fabs(value))
316	except:
317	err[i] = 1.0
318	print "Error getting error", value, x[i]
319
320	new_plot = Data1D(x, y)
321	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
322	new_plot.name = IQ_FIT_LABEL
323	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
324	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
325	title = "I(q)"
326	new_plot.title = title
327
328	# If we have a group ID, use it
329	if pr.info.has_key("plot_group_id"):
330	new_plot.group_id = pr.info["plot_group_id"]
331	new_plot.id = IQ_FIT_LABEL
332	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
333
334	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
335
336	# If we have used slit smearing, plot the smeared I(q) too
337	if pr.slit_width>0 or pr.slit_height>0:
338	x = pylab.arange(minq, maxq, maxq/301.0)
339	y = numpy.zeros(len(x))
340	err = numpy.zeros(len(x))
341	for i in range(len(x)):
342	value = pr.iq_smeared(out, x[i])
343	y[i] = value
344	try:
345	err[i] = math.sqrt(math.fabs(value))
346	except:
347	err[i] = 1.0
348	print "Error getting error", value, x[i]
349
350	new_plot = Data1D(x, y)
351	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
352	new_plot.name = IQ_SMEARED_LABEL
353	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
354	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
355	# If we have a group ID, use it
356	if pr.info.has_key("plot_group_id"):
357	new_plot.group_id = pr.info["plot_group_id"]
358	new_plot.id = IQ_SMEARED_LABEL
359	new_plot.title = title
360	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
361	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
362
363	def _on_pr_npts(self, evt):
364	"""
365	Redisplay P(r) with a different number of points
366	"""
367	from inversion_panel import PrDistDialog
368	dialog = PrDistDialog(None, -1)
369	dialog.set_content(self._pr_npts)
370	if dialog.ShowModal() == wx.ID_OK:
371	self._pr_npts= dialog.get_content()
372	dialog.Destroy()
373	self.show_pr(self._last_out, self._last_pr, self._last_cov)
374	else:
375	dialog.Destroy()
376
377
378	def show_pr(self, out, pr, cov=None):
379	"""
380	"""
381	# Show P(r)
382	x = pylab.arange(0.0, pr.d_max, pr.d_max/self._pr_npts)
383
384	y = numpy.zeros(len(x))
385	dy = numpy.zeros(len(x))
386	y_true = numpy.zeros(len(x))
387
388	sum = 0.0
389	pmax = 0.0
390	cov2 = numpy.ascontiguousarray(cov)
391
392	for i in range(len(x)):
393	if cov2==None:
394	value = pr.pr(out, x[i])
395	else:
396	(value, dy[i]) = pr.pr_err(out, cov2, x[i])
397	sum += value*pr.d_max/len(x)
398
399	# keep track of the maximum P(r) value
400	if value>pmax:
401	pmax = value
402
403	y[i] = value
404
405	if self._normalize_output==True:
406	y = y/sum
407	dy = dy/sum
408	elif self._scale_output_unity==True:
409	y = y/pmax
410	dy = dy/pmax
411
412	if cov2==None:
413	new_plot = Data1D(x, y)
414	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
415	else:
416	new_plot = Data1D(x, y, dy=dy)
417	new_plot.name = PR_FIT_LABEL
418	new_plot.xaxis("\\rm{r}", 'A')
419	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
420	new_plot.title = "P(r) fit"
421	new_plot.id = PR_FIT_LABEL
422	# Make sure that the plot is linear
423	new_plot.xtransform = "x"
424	new_plot.ytransform = "y"
425	new_plot.group_id = GROUP_ID_PR_FIT
426	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
427	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="P(r) fit"))
428	return x, pr.d_max
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	new_plot.group_id = self._added_plots[plot].group_id
677	new_plot.id = self._added_plots[plot].id
678	new_plot.title = self._added_plots[plot].title
679	new_plot.name = self._added_plots[plot].name
680	new_plot.xaxis("\\rm{r}", 'A')
681	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
682	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
683	wx.PostEvent(self.parent,
684	NewPlotEvent(plot=new_plot, update=True,
685	title=self._added_plots[plot].name))
686
687	def _on_scale_unity(self, evt):
688	"""
689	Scale the maximum P(r) value on each displayed plot to 1.
690
691	:param evt: Menu event
692
693	"""
694	self._scale_output_unity = True
695	self._normalize_output = False
696
697	self.show_pr(self._last_out, self._last_pr, self._last_cov)
698
699	# Now scale the added plots too
700	for plot in self._added_plots:
701	_max = 0
702	for y in self._added_plots[plot].y:
703	if y>_max:
704	_max = y
705	y = self._added_plots[plot].y/_max
706
707	new_plot = Data1D(self._added_plots[plot].x, y)
708	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
709	new_plot.name = self._added_plots[plot].name
710	new_plot.xaxis("\\rm{r}", 'A')
711	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
712	self.parent.update_theory(data_id=self.data_id,theory=new_plot)
713	wx.PostEvent(self.parent,
714	NewPlotEvent(plot=new_plot, update=True,
715	title=self._added_plots[plot].name))
716
717
718	def _on_add_data(self, evt):
719	"""
720	Add a data curve to the plot
721
722	:WARNING: this will be removed once guiframe.plotting has
723	its full functionality
724	"""
725	path = self.choose_file()
726	if path==None:
727	return
728
729	#x, y, err = self.parent.load_ascii_1D(path)
730	# Use data loader to load file
731	try:
732	dataread = Loader().load(path)
733	x = None
734	y = None
735	err = None
736	if dataread.__class__.__name__ == 'Data1D':
737	x = dataread.x
738	y = dataread.y
739	err = dataread.dy
740	else:
741	if isinstance(dataread, list) and len(dataread)>0:
742	x = dataread[0].x
743	y = dataread[0].y
744	err = dataread[0].dy
745	msg = "PrView only allows a single data set at a time. "
746	msg += "Only the first data set was loaded."
747	wx.PostEvent(self.parent, StatusEvent(status=msg))
748	else:
749	msg = "This tool can only read 1D data"
750	wx.PostEvent(self.parent, StatusEvent(status=msg))
751	return
752
753	except:
754	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
755	return
756
757	filename = os.path.basename(path)
758	new_plot = Data1D(x, y)
759	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
760	new_plot.name = filename
761	new_plot.xaxis("\\rm{r}", 'A')
762	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
763	# Store a ref to the plottable for later use
764	self._added_plots[filename] = new_plot
765	self._default_Iq[filename] = numpy.copy(y)
766	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=filename))
767
768	def start_thread(self):
769	"""
770	"""
771	from pr_thread import CalcPr
772	from copy import deepcopy
773
774	# If a thread is already started, stop it
775	if self.calc_thread != None and self.calc_thread.isrunning():
776	self.calc_thread.stop()
777
778	pr = self.pr.clone()
779	self.calc_thread = CalcPr(pr, self.nfunc,
780	error_func=self._thread_error,
781	completefn=self._completed, updatefn=None)
782	self.calc_thread.queue()
783	self.calc_thread.ready(2.5)
784
785	def _thread_error(self, error):
786	"""
787	"""
788	wx.PostEvent(self.parent, StatusEvent(status=error))
789
790	def _estimate_completed(self, alpha, message, elapsed):
791	"""
792	Parameter estimation completed,
793	display the results to the user
794
795	:param alpha: estimated best alpha
796	:param elapsed: computation time
797
798	"""
799	# Save useful info
800	self.elapsed = elapsed
801	self.control_panel.alpha_estimate = alpha
802	if not message==None:
803	wx.PostEvent(self.parent, StatusEvent(status=str(message)))
804	self.perform_estimateNT()
805
806	def _estimateNT_completed(self, nterms, alpha, message, elapsed):
807	"""
808	Parameter estimation completed,
809	display the results to the user
810
811	:param alpha: estimated best alpha
812	:param nterms: estimated number of terms
813	:param elapsed: computation time
814
815	"""
816	# Save useful info
817	self.elapsed = elapsed
818	self.control_panel.nterms_estimate = nterms
819	self.control_panel.alpha_estimate = alpha
820	if not message==None:
821	wx.PostEvent(self.parent, StatusEvent(status=str(message)))
822
823	def _completed(self, out, cov, pr, elapsed):
824	"""
825	wxCallAfter Method called with the results when the inversion
826	is done
827
828	:param out: output coefficient for the base functions
829	:param cov: covariance matrix
830	:param pr: Invertor instance
831	:param elapsed: time spent computing
832	"""
833	wx.CallAfter(self,_completed_call, out, cov, pr, elapsed)
834
835	def _completed_call(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 = pr.info["plot_group_id"]
891	new_plot.id = self.data_id
892	self.parent.update_theory(data_id=self.data_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 = GROUP_ID_IQ_DATA
940	new_plot.id = self.data_id
941	new_plot.title = "I(q)"
942	wx.PostEvent(self.parent,
943	NewPlotEvent(plot=new_plot, title="I(q)", reset=reset))
944
945	self.current_plottable = new_plot
946	# Get Q range
947	self.control_panel.q_min = self.pr.x.min()
948	self.control_panel.q_max = self.pr.x.max()
949
950	def save_data(self, filepath, prstate=None):
951	"""
952	Save data in provided state object.
953
954	:TODO: move the state code away from inversion_panel and move it here.
955	Then remove the "prstate" input and make this method private.
956
957	:param filepath: path of file to write to
958	:param prstate: P(r) inversion state
959
960	"""
961	#TODO: do we need this or can we use DataLoader.loader.save directly?
962
963	# Add output data and coefficients to state
964	prstate.coefficients = self._last_out
965	prstate.covariance = self._last_cov
966
967	# Write the output to file
968	# First, check that the data is of the right type
969	if issubclass(self.current_plottable.__class__,
970	DataLoader.data_info.Data1D):
971	self.state_reader.write(filepath, self.current_plottable, prstate)
972	else:
973	msg = "pr.save_data: the data being saved is not a"
974	msg += " DataLoader.data_info.Data1D object"
975	raise RuntimeError, msg
976
977
978	def setup_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None,
979	bck=False, height=0, width=0):
980	"""
981	"""
982	self.alpha = alpha
983	self.nfunc = nfunc
984	self.max_length = d_max
985	self.q_min = q_min
986	self.q_max = q_max
987	self.has_bck = bck
988	self.slit_height = height
989	self.slit_width = width
990
991	try:
992	pr = self._create_plot_pr()
993	if not pr==None:
994	self.pr = pr
995	self.perform_inversion()
996	except:
997	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
998
999	def estimate_plot_inversion(self, alpha, nfunc, d_max,
1000	q_min=None, q_max=None,
1001	bck=False, height=0, width=0):
1002	"""
1003	"""
1004	self.alpha = alpha
1005	self.nfunc = nfunc
1006	self.max_length = d_max
1007	self.q_min = q_min
1008	self.q_max = q_max
1009	self.has_bck = bck
1010	self.slit_height = height
1011	self.slit_width = width
1012
1013	try:
1014	pr = self._create_plot_pr()
1015	if not pr==None:
1016	self.pr = pr
1017	self.perform_estimate()
1018	except:
1019	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1020
1021	def _create_plot_pr(self, estimate=False):
1022	"""
1023	Create and prepare invertor instance from
1024	a plottable data set.
1025
1026	:param path: path of the file to read in
1027
1028	"""
1029	# Sanity check
1030	if self.current_plottable is None:
1031	msg = "Please load a valid data set before proceeding."
1032	wx.PostEvent(self.parent, StatusEvent(status=msg))
1033	return None
1034
1035	# Get the data from the chosen data set and perform inversion
1036	pr = Invertor()
1037	pr.d_max = self.max_length
1038	pr.alpha = self.alpha
1039	pr.q_min = self.q_min
1040	pr.q_max = self.q_max
1041	pr.x = self.current_plottable.x
1042	pr.y = self.current_plottable.y
1043	pr.has_bck = self.has_bck
1044	pr.slit_height = self.slit_height
1045	pr.slit_width = self.slit_width
1046
1047	# Keep track of the plot window title to ensure that
1048	# we can overlay the plots
1049	pr.info["plot_group_id"] = self.current_plottable.group_id
1050
1051	# Fill in errors if none were provided
1052	err = self.current_plottable.dy
1053	all_zeros = True
1054	if err == None:
1055	err = numpy.zeros(len(pr.y))
1056	else:
1057	for i in range(len(err)):
1058	if err[i]>0:
1059	all_zeros = False
1060
1061	if all_zeros:
1062	scale = None
1063	min_err = 0.0
1064	for i in range(len(pr.y)):
1065	# Scale the error so that we can fit over several decades of Q
1066	if scale==None:
1067	scale = 0.05*math.sqrt(pr.y[i])
1068	min_err = 0.01*pr.y[i]
1069	err[i] = scale*math.sqrt( math.fabs(pr.y[i]) ) + min_err
1070	message = "The loaded file had no error bars, "
1071	message += "statistical errors are assumed."
1072	wx.PostEvent(self.parent, StatusEvent(status=message))
1073
1074	pr.err = err
1075
1076	return pr
1077
1078
1079	def setup_file_inversion(self, alpha, nfunc, d_max, data,
1080	path=None, q_min=None, q_max=None,
1081	bck=False, height=0, width=0):
1082	"""
1083	"""
1084	self.alpha = alpha
1085	self.nfunc = nfunc
1086	self.max_length = d_max
1087	self.q_min = q_min
1088	self.q_max = q_max
1089	self.has_bck = bck
1090	self.slit_height = height
1091	self.slit_width = width
1092
1093	try:
1094	#pr = self._create_file_pr(path)
1095	pr = self._create_file_pr(data)
1096	if not pr==None:
1097	self.pr = pr
1098	self.perform_inversion()
1099	except:
1100	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1101
1102	def estimate_file_inversion(self, alpha, nfunc, d_max, data,
1103	path=None, q_min=None, q_max=None,
1104	bck=False, height=0, width=0):
1105	"""
1106	"""
1107	self.alpha = alpha
1108	self.nfunc = nfunc
1109	self.max_length = d_max
1110	self.q_min = q_min
1111	self.q_max = q_max
1112	self.has_bck = bck
1113	self.slit_height = height
1114	self.slit_width = width
1115
1116	try:
1117	pr = self._create_file_pr(data)
1118	#pr = self._create_file_pr(path)
1119	if not pr is None:
1120	self.pr = pr
1121	self.perform_estimate()
1122	except:
1123	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1124
1125
1126	def _create_file_pr(self, data):
1127	"""
1128	Create and prepare invertor instance from
1129	a file data set.
1130
1131	:param path: path of the file to read in
1132
1133	"""
1134	# Load data
1135	#if os.path.isfile(path):
1136	"""
1137	if self._current_file_data is not None \
1138	and self._current_file_data.path==path:
1139	# Protect against corrupted data from
1140	# previous failed load attempt
1141	if self._current_file_data.x is None:
1142	return None
1143	x = self._current_file_data.x
1144	y = self._current_file_data.y
1145	err = self._current_file_data.err
1146
1147	message = "The data from this file has already been loaded."
1148	wx.PostEvent(self.parent, StatusEvent(status=message))
1149	else:
1150	"""
1151	# Reset the status bar so that we don't get mixed up
1152	# with old messages.
1153	#TODO: refactor this into a proper status handling
1154	wx.PostEvent(self.parent, StatusEvent(status=''))
1155	try:
1156	class FileData:
1157	x = None
1158	y = None
1159	err = None
1160	path = None
1161	def __init__(self, path):
1162	self.path = path
1163
1164	self._current_file_data = FileData(data.path)
1165	self._current_file_data.x = data.x
1166	self._current_file_data.y = data.y
1167	self._current_file_data.err = data.dy
1168	x, y, err = data.x, data.y, data.dy
1169	except:
1170	load_error(sys.exc_value)
1171	return None
1172
1173	# If the file contains no data, just return
1174	if x is None or len(x) == 0:
1175	load_error("The loaded file contains no data")
1176	return None
1177
1178	# If we have not errors, add statistical errors
1179	if y is not None:
1180	if err == None or numpy.all(err) == 0:
1181	err = numpy.zeros(len(y))
1182	scale = None
1183	min_err = 0.0
1184	for i in range(len(y)):
1185	# Scale the error so that we can fit over several decades of Q
1186	if scale == None:
1187	scale = 0.05 * math.sqrt(y[i])
1188	min_err = 0.01 * y[i]
1189	err[i] = scale * math.sqrt(math.fabs(y[i])) + min_err
1190	message = "The loaded file had no error bars, "
1191	message += "statistical errors are assumed."
1192	wx.PostEvent(self.parent, StatusEvent(status=message))
1193
1194	try:
1195	# Get the data from the chosen data set and perform inversion
1196	pr = Invertor()
1197	pr.d_max = self.max_length
1198	pr.alpha = self.alpha
1199	pr.q_min = self.q_min
1200	pr.q_max = self.q_max
1201	pr.x = x
1202	pr.y = y
1203	pr.err = err
1204	pr.has_bck = self.has_bck
1205	pr.slit_height = self.slit_height
1206	pr.slit_width = self.slit_width
1207	return pr
1208	except:
1209	load_error(sys.exc_value)
1210	return None
1211
1212	def perform_estimate(self):
1213	"""
1214	"""
1215	from pr_thread import EstimatePr
1216	from copy import deepcopy
1217
1218	# If a thread is already started, stop it
1219	if self.estimation_thread != None and \
1220	self.estimation_thread.isrunning():
1221	self.estimation_thread.stop()
1222
1223	pr = self.pr.clone()
1224	self.estimation_thread = EstimatePr(pr, self.nfunc,
1225	error_func=self._thread_error,
1226	completefn = self._estimate_completed,
1227	updatefn = None)
1228	self.estimation_thread.queue()
1229	self.estimation_thread.ready(2.5)
1230
1231	def perform_estimateNT(self):
1232	"""
1233	"""
1234	from pr_thread import EstimateNT
1235	from copy import deepcopy
1236
1237	# If a thread is already started, stop it
1238	if self.estimation_thread != None and self.estimation_thread.isrunning():
1239	self.estimation_thread.stop()
1240
1241	pr = self.pr.clone()
1242	# Skip the slit settings for the estimation
1243	# It slows down the application and it doesn't change the estimates
1244	pr.slit_height = 0.0
1245	pr.slit_width = 0.0
1246	self.estimation_thread = EstimateNT(pr, self.nfunc,
1247	error_func=self._thread_error,
1248	completefn = self._estimateNT_completed,
1249	updatefn = None)
1250	self.estimation_thread.queue()
1251	self.estimation_thread.ready(2.5)
1252
1253	def perform_inversion(self):
1254	"""
1255	"""
1256	# Time estimate
1257	#estimated = self.elapsedself.nfunc*2
1258	#message = "Computation time may take up to %g seconds" % self.elapsed
1259	#wx.PostEvent(self.parent, StatusEvent(status=message))
1260
1261	# Start inversion thread
1262	self.start_thread()
1263	return
1264
1265	out, cov = self.pr.lstsq(self.nfunc)
1266
1267	# Save useful info
1268	self.elapsed = self.pr.elapsed
1269
1270	for i in range(len(out)):
1271	try:
1272	print "%d: %g +- %g" % (i, out[i],
1273	math.sqrt(math.fabs(cov[i][i])))
1274	except:
1275	print "%d: %g +- ?" % (i, out[i])
1276
1277	# Make a plot of I(q) data
1278	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
1279	new_plot.name = "I_{obs}(q)"
1280	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
1281	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
1282	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
1283	# Show I(q) fit
1284	self.show_iq(out, self.pr)
1285	# Show P(r) fit
1286	x_values, x_range = self.show_pr(out, self.pr, cov=cov)
1287
1288	def _on_context_inversion(self, event):
1289	"""
1290	"""
1291	panel = event.GetEventObject()
1292	Plugin.on_perspective(self, event=event)
1293
1294	# If we have more than one displayed plot, make the user choose
1295	if len(panel.plots) >= 1 and \
1296	panel.graph.selected_plottable in panel.plots:
1297	dataset = panel.plots[panel.graph.selected_plottable].name
1298	else:
1299	logging.info("Prview Error: No data is available")
1300	return
1301
1302	# Store a reference to the current plottable
1303	# If we have a suggested value, use it.
1304	try:
1305	estimate = float(self.control_panel.alpha_estimate)
1306	self.control_panel.alpha = estimate
1307	except:
1308	self.control_panel.alpha = self.alpha
1309	logging.info("Prview :Alpha Not estimate yet")
1310	pass
1311	try:
1312	estimate = int(self.control_panel.nterms_estimate)
1313	self.control_panel.nfunc = estimate
1314	except:
1315	self.control_panel.nfunc = self.nfunc
1316	logging.info("Prview : ntemrs Not estimate yet")
1317	pass
1318
1319	self.current_plottable = panel.plots[panel.graph.selected_plottable]
1320	self.set_data([self.current_plottable])
1321	self.control_panel.plotname = dataset
1322	#self.control_panel.nfunc = self.nfunc
1323	self.control_panel.d_max = self.max_length
1324	#self.parent.set_perspective(self.perspective)
1325	self.control_panel._on_invert(None)
1326
1327	def get_panels(self, parent):
1328	"""
1329	Create and return a list of panel objects
1330	"""
1331	from inversion_panel import InversionControl
1332
1333	self.parent = parent
1334	self.control_panel = InversionControl(self.parent, -1,
1335	style=wx.RAISED_BORDER,
1336	standalone=self.standalone)
1337	self.control_panel.set_manager(self)
1338	self.control_panel.nfunc = self.nfunc
1339	self.control_panel.d_max = self.max_length
1340	self.control_panel.alpha = self.alpha
1341	self.perspective = []
1342	self.perspective.append(self.control_panel.window_name)
1343
1344	return [self.control_panel]
1345
1346	def set_data(self, data_list=None):
1347	"""
1348	receive a list of data to compute pr
1349	"""
1350	if data_list is None:
1351	data_list = []
1352	if len(data_list) >= 1:
1353	if len(data_list) == 1:
1354	data = data_list[0]
1355	else:
1356	msg = "Pr panel does not allow multiple Data.\n"
1357	msg += "Please select one!\n"
1358	from pr_widgets import DataDialog
1359	dlg = DataDialog(data_list=data_list, text=msg)
1360	if dlg.ShowModal() == wx.ID_OK:
1361	data = dlg.get_data()
1362	if data is None:
1363	return
1364	if issubclass(data.__class__, Data1D):
1365	try:
1366	wx.PostEvent(self.parent, NewPlotEvent(action='remove',
1367	group_id=GROUP_ID_IQ_DATA,
1368	id=self.data_id))
1369
1370	self.data_id = data.id
1371	self.control_panel._change_file(evt=None, data=data)
1372	except:
1373	msg = "Prview Set_data: " + str(sys.exc_value)
1374	wx.PostEvent(self.parent, StatusEvent(status=msg,
1375	info="error"))
1376	else:
1377	msg = "Pr cannot be computed for data of "
1378	msg += "type %s" % (data_list[0].__class__.__name__)
1379	wx.PostEvent(self.parent,
1380	StatusEvent(status=msg, info='error'))
1381	else:
1382	msg = "Pr contain no data"
1383	wx.PostEvent(self.parent, StatusEvent(status=msg, info='warning'))
1384
1385	def post_init(self):
1386	"""
1387	Post initialization call back to close the loose ends
1388	[Somehow openGL needs this call]
1389	"""
1390	if self.standalone:
1391	self.parent.set_perspective(self.perspective)
1392
1393	if __name__ == "__main__":
1394	i = Plugin()
1395	print i.perform_estimateNT()
1396
1397
1398
1399

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

SasView

source: sasview/prview/perspectives/pr/pr.py @ 380ff18

Download in other formats: