pr.py @ 52725d6

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

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

SasView

source: sasview/prview/perspectives/pr/pr.py @ 52725d6

Download in other formats: