pr.py @ b659551

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 b659551 was b659551, checked in by Mathieu Doucet <doucetm@…>, 16 years ago
Added reset button, added option to change number of points on the P(r) plot, fixed minor problems with displaying error bars.
Property mode set to `100644`
File size: 22.6 KB

Line
1	#TODO: Use simview to generate P(r) and I(q) pairs in sansview.
2	# Make sure the option of saving each curve is available
3	# Use the I(q) curve as input and compare the output to P(r)
4
5	import os
6	import sys
7	import wx
8	from sans.guitools.plottables import Data1D, Theory1D
9	from sans.guicomm.events import NewPlotEvent, StatusEvent
10	import math, numpy
11	from sans.pr.invertor import Invertor
12
13	class Plugin:
14
15	DEFAULT_ALPHA = 0.0001
16	DEFAULT_NFUNC = 10
17	DEFAULT_DMAX = 140.0
18
19	def __init__(self):
20	## Plug-in name
21	self.sub_menu = "Pr inversion"
22
23	## Reference to the parent window
24	self.parent = None
25
26	## Simulation window manager
27	self.simview = None
28
29	## List of panels for the simulation perspective (names)
30	self.perspective = []
31
32	## State data
33	self.alpha = self.DEFAULT_ALPHA
34	self.nfunc = self.DEFAULT_NFUNC
35	self.max_length = self.DEFAULT_DMAX
36	self.q_min = None
37	self.q_max = None
38	## Remember last plottable processed
39	self.last_data = "sphere_60_q0_2.txt"
40	## Time elapsed for last computation [sec]
41	# Start with a good default
42	self.elapsed = 0.022
43
44	## Current invertor
45	self.invertor = None
46	## Calculation thread
47	self.calc_thread = None
48	## Estimation thread
49	self.estimation_thread = None
50	## Result panel
51	self.control_panel = None
52	## Currently views plottable
53	self.current_plottable = None
54	## Number of P(r) points to display on the output plot
55	self._pr_npts = 51
56
57	def populate_menu(self, id, owner):
58	"""
59	Create a menu for the plug-in
60	"""
61	return []
62	import wx
63	shapes = wx.Menu()
64
65	id = wx.NewId()
66	shapes.Append(id, '&Sphere test')
67	wx.EVT_MENU(owner, id, self._fit_pr)
68
69	return [(id, shapes, "P(r)")]
70
71	def help(self, evt):
72	"""
73	Show a general help dialog.
74	TODO: replace the text with a nice image
75	"""
76	from inversion_panel import HelpDialog
77	dialog = HelpDialog(None, -1)
78	if dialog.ShowModal() == wx.ID_OK:
79	dialog.Destroy()
80	else:
81	dialog.Destroy()
82
83	def _fit_pr(self, evt):
84	from sans.pr.invertor import Invertor
85	import numpy
86	import pylab
87	import math
88	from sans.guicomm.events import NewPlotEvent
89	from sans.guitools.plottables import Data1D, Theory1D
90
91	# Generate P(r) for sphere
92	radius = 60.0
93	d_max = 2*radius
94
95
96	r = pylab.arange(0.01, d_max, d_max/51.0)
97	M = len(r)
98	y = numpy.zeros(M)
99	pr_err = numpy.zeros(M)
100
101	sum = 0.0
102	for j in range(M):
103	value = self.pr_theory(r[j], radius)
104	sum += value
105	y[j] = value
106	pr_err[j] = math.sqrt(y[j])
107
108
109	y = y/sum*d_max/len(r)
110
111
112
113	# Perform fit
114	pr = Invertor()
115	pr.d_max = d_max
116	pr.alpha = 0
117	pr.x = r
118	pr.y = y
119	pr.err = pr_err
120	out, cov = pr.pr_fit()
121	for i in range(len(out)):
122	print "%g +- %g" % (out[i], math.sqrt(cov[i][i]))
123
124
125	# Show input P(r)
126	new_plot = Data1D(pr.x, pr.y, dy=pr.err)
127	new_plot.name = "P_{obs}(r)"
128	new_plot.xaxis("\\rm{r}", 'A')
129	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
130	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Pr"))
131
132	# Show P(r) fit
133	self.show_pr(out, pr)
134
135	# Show I(q) fit
136	q = pylab.arange(0.001, 0.1, 0.01/51.0)
137	self.show_iq(out, pr, q)
138
139
140	def show_shpere(self, x, radius=70.0, x_range=70.0):
141	import numpy
142	import pylab
143	import math
144	from sans.guicomm.events import NewPlotEvent
145	from sans.guitools.plottables import Data1D, Theory1D
146	# Show P(r)
147	y_true = numpy.zeros(len(x))
148
149	sum_true = 0.0
150	for i in range(len(x)):
151	y_true[i] = self.pr_theory(x[i], radius)
152	sum_true += y_true[i]
153
154	y_true = y_true/sum_true*x_range/len(x)
155
156	# Show the theory P(r)
157	new_plot = Theory1D(x, y_true)
158	new_plot.name = "P_{true}(r)"
159	new_plot.xaxis("\\rm{r}", 'A')
160	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
161
162
163	#Put this call in plottables/guitools
164	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Sphere P(r)"))
165
166	def get_npts(self):
167	"""
168	Returns the number of points in the I(q) data
169	"""
170	try:
171	return len(self.pr.x)
172	except:
173	return 0
174
175	def show_iq(self, out, pr, q=None):
176	import numpy
177	import pylab
178	import math
179	from sans.guicomm.events import NewPlotEvent
180	from sans.guitools.plottables import Data1D, Theory1D
181
182	qtemp = pr.x
183	if not q==None:
184	qtemp = q
185
186	# Make a plot
187	maxq = -1
188	for q_i in qtemp:
189	if q_i>maxq:
190	maxq=q_i
191
192	minq = 0.001
193
194	# Check for user min/max
195	if not pr.q_min==None:
196	minq = pr.q_min
197	if not pr.q_max==None:
198	maxq = pr.q_max
199
200	x = pylab.arange(minq, maxq, maxq/301.0)
201	y = numpy.zeros(len(x))
202	err = numpy.zeros(len(x))
203	for i in range(len(x)):
204	value = pr.iq(out, x[i])
205	y[i] = value
206	try:
207	err[i] = math.sqrt(math.fabs(value))
208	except:
209	err[i] = 1.0
210	print "Error getting error", value, x[i]
211
212	new_plot = Theory1D(x, y)
213	new_plot.name = "I_{fit}(q)"
214	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
215	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
216	#new_plot.group_id = "test group"
217	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
218
219
220	def _on_pr_npts(self, evt):
221	"""
222	Redisplay P(r) with a different number of points
223	"""
224	from inversion_panel import PrDistDialog
225	dialog = PrDistDialog(None, -1)
226	dialog.set_content(self._pr_npts)
227	if dialog.ShowModal() == wx.ID_OK:
228	self._pr_npts= dialog.get_content()
229	dialog.Destroy()
230	self.show_pr(self.pr.out, self.pr, self.pr.cov)
231	else:
232	dialog.Destroy()
233
234
235	def show_pr(self, out, pr, cov=None):
236	import numpy
237	import pylab
238	import math
239	from sans.guicomm.events import NewPlotEvent
240	from sans.guitools.plottables import Data1D, Theory1D
241
242	# Show P(r)
243	x = pylab.arange(0.0, pr.d_max, pr.d_max/self._pr_npts)
244
245	y = numpy.zeros(len(x))
246	dy = numpy.zeros(len(x))
247	y_true = numpy.zeros(len(x))
248
249	sum = 0.0
250	cov2 = numpy.ascontiguousarray(cov)
251
252	for i in range(len(x)):
253	if cov2==None:
254	value = pr.pr(out, x[i])
255	else:
256	(value, dy[i]) = pr.pr_err(out, cov2, x[i])
257	sum += value*pr.d_max/len(x)
258	y[i] = value
259
260	y = y/sum
261	dy = dy/sum
262
263	if cov2==None:
264	new_plot = Theory1D(x, y)
265	else:
266	new_plot = Data1D(x, y, dy=dy)
267	new_plot.name = "P_{fit}(r)"
268	new_plot.xaxis("\\rm{r}", 'A')
269	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
270
271	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="P(r) fit"))
272
273	return x, pr.d_max
274
275
276	def choose_file(self):
277	"""
278
279	"""
280	#TODO: this should be in a common module
281	return self.parent.choose_file()
282
283	def load(self, path = "sphere_60_q0_2.txt"):
284	import numpy, math, sys
285	# Read the data from the data file
286	data_x = numpy.zeros(0)
287	data_y = numpy.zeros(0)
288	data_err = numpy.zeros(0)
289	scale = None
290	min_err = 0.0
291	if not path == None:
292	input_f = open(path,'r')
293	buff = input_f.read()
294	lines = buff.split('\n')
295	for line in lines:
296	try:
297	toks = line.split()
298	x = float(toks[0])
299	y = float(toks[1])
300	if len(toks)>2:
301	err = float(toks[2])
302	else:
303	if scale==None:
304	scale = 0.05*math.sqrt(y)
305	#scale = 0.05/math.sqrt(y)
306	min_err = 0.01*y
307	err = scale*math.sqrt(y)+min_err
308	#err = 0
309
310	data_x = numpy.append(data_x, x)
311	data_y = numpy.append(data_y, y)
312	data_err = numpy.append(data_err, err)
313	except:
314	pass
315
316	return data_x, data_y, data_err
317
318	def pr_theory(self, r, R):
319	"""
320
321	"""
322	if r<=2*R:
323	return 12.0* ((0.5r/R)2) ((1.0-0.5r/R)2) ( 2.0 + 0.5*r/R )
324	else:
325	return 0.0
326
327	def get_context_menu(self, graph=None):
328	"""
329	Get the context menu items available for P(r)
330	@param graph: the Graph object to which we attach the context menu
331	@return: a list of menu items with call-back function
332	"""
333	# Look whether this Graph contains P(r) data
334	for item in graph.plottables:
335	if item.name=="P_{fit}(r)":
336
337	return [["Compute P(r)", "Compute P(r) from distribution", self._on_context_inversion],
338	["Add P(r) data", "Load a data file and display it on this plot", self._on_add_data],
339	["Change number of P(r) points", "Change the number of points on the P(r) output", self._on_pr_npts]]
340
341	return [["Compute P(r)", "Compute P(r) from distribution", self._on_context_inversion]]
342
343	def _on_add_data(self, evt):
344	"""
345	Add a data curve to the plot
346	WARNING: this will be removed once guiframe.plotting has its full functionality
347	"""
348	path = self.choose_file()
349	if path==None:
350	return
351
352	x, y, err = self.parent.load_ascii_1D(path)
353
354	#new_plot = Data1D(x, y, dy=err)
355	new_plot = Theory1D(x, y)
356	new_plot.name = "P_{loaded}(r)"
357	new_plot.xaxis("\\rm{r}", 'A')
358	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
359
360	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="P(r) fit"))
361
362
363
364	def start_thread(self):
365	from pr_thread import CalcPr
366	from copy import deepcopy
367
368	# If a thread is already started, stop it
369	if self.calc_thread != None and self.calc_thread.isrunning():
370	self.calc_thread.stop()
371
372	pr = self.pr.clone()
373	self.calc_thread = CalcPr(pr, self.nfunc, error_func=self._thread_error, completefn=self._completed, updatefn=None)
374	self.calc_thread.queue()
375	self.calc_thread.ready(2.5)
376
377	def _thread_error(self, error):
378	wx.PostEvent(self.parent, StatusEvent(status=error))
379
380	def _estimate_completed(self, alpha, message, elapsed):
381	"""
382	Parameter estimation completed,
383	display the results to the user
384	@param alpha: estimated best alpha
385	@param elapsed: computation time
386	"""
387	# Save useful info
388	self.elapsed = elapsed
389	self.control_panel.alpha_estimate = alpha
390	if not message==None:
391	wx.PostEvent(self.parent, StatusEvent(status=str(message)))
392
393	def _completed(self, out, cov, pr, elapsed):
394	"""
395	Method called with the results when the inversion
396	is done
397
398	@param out: output coefficient for the base functions
399	@param cov: covariance matrix
400	@param pr: Invertor instance
401	@param elapsed: time spent computing
402	"""
403	from copy import deepcopy
404	# Save useful info
405	self.elapsed = elapsed
406	# Save Pr invertor
407	self.pr = pr
408
409	message = "Computation completed in %g seconds [chi2=%g]" % (elapsed, pr.chi2)
410	wx.PostEvent(self.parent, StatusEvent(status=message))
411
412	cov = numpy.ascontiguousarray(cov)
413
414	# Show result on control panel
415	self.control_panel.chi2 = pr.chi2
416	self.control_panel.elapsed = elapsed
417	self.control_panel.oscillation = pr.oscillations(out)
418	#print "OSCILL", pr.oscillations(out)
419	print "PEAKS:", pr.get_peaks(out)
420	self.control_panel.positive = pr.get_positive(out)
421	self.control_panel.pos_err = pr.get_pos_err(out, cov)
422
423	for i in range(len(out)):
424	try:
425	print "%d: %g +- %g" % (i, out[i], math.sqrt(math.fabs(cov[i][i])))
426	except:
427	print sys.exc_value
428	print "%d: %g +- ?" % (i, out[i])
429
430	# Make a plot of I(q) data
431	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
432	new_plot.name = "I_{obs}(q)"
433	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
434	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
435	#new_plot.group_id = "test group"
436	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
437
438	# Show I(q) fit
439	self.show_iq(out, self.pr)
440
441	# Show P(r) fit
442	x_values, x_range = self.show_pr(out, self.pr, cov)
443
444	# Popup result panel
445	#result_panel = InversionResults(self.parent, -1, style=wx.RAISED_BORDER)
446
447	def show_data(self, path=None):
448	if not path==None:
449	self._create_file_pr(path)
450
451	# Make a plot of I(q) data
452	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
453	new_plot.name = "I_{obs}(q)"
454	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
455	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
456	#new_plot.group_id = "test group"
457	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
458
459
460	def setup_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None):
461	self.alpha = alpha
462	self.nfunc = nfunc
463	self.max_length = d_max
464	self.q_min = q_min
465	self.q_max = q_max
466
467	try:
468	self._create_plot_pr()
469	self.perform_inversion()
470	except:
471	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
472
473	def estimate_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None):
474	self.alpha = alpha
475	self.nfunc = nfunc
476	self.max_length = d_max
477	self.q_min = q_min
478	self.q_max = q_max
479
480	try:
481	self._create_plot_pr()
482	self.perform_estimate()
483	except:
484	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
485
486	def _create_plot_pr(self):
487	"""
488	Create and prepare invertor instance from
489	a plottable data set.
490	@param path: path of the file to read in
491	"""
492	# Get the data from the chosen data set and perform inversion
493	pr = Invertor()
494	pr.d_max = self.max_length
495	pr.alpha = self.alpha
496	pr.q_min = self.q_min
497	pr.q_max = self.q_max
498	pr.x = self.current_plottable.x
499	pr.y = self.current_plottable.y
500
501	# Fill in errors if none were provided
502	if self.current_plottable.dy == None:
503	print "no error", self.current_plottable.name
504	y = numpy.zeros(len(pr.y))
505	for i in range(len(pr.y)):
506	y[i] = math.sqrt(pr.y[i])
507	pr.err = y
508	else:
509	pr.err = self.current_plottable.dy
510
511	self.pr = pr
512
513
514	def setup_file_inversion(self, alpha, nfunc, d_max, path, q_min=None, q_max=None):
515	self.alpha = alpha
516	self.nfunc = nfunc
517	self.max_length = d_max
518	self.q_min = q_min
519	self.q_max = q_max
520
521	try:
522	if self._create_file_pr(path):
523	self.perform_inversion()
524	except:
525	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
526
527	def estimate_file_inversion(self, alpha, nfunc, d_max, path, q_min=None, q_max=None):
528	self.alpha = alpha
529	self.nfunc = nfunc
530	self.max_length = d_max
531	self.q_min = q_min
532	self.q_max = q_max
533
534	try:
535	if self._create_file_pr(path):
536	self.perform_estimate()
537	except:
538	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
539
540
541	def _create_file_pr(self, path):
542	"""
543	Create and prepare invertor instance from
544	a file data set.
545	@param path: path of the file to read in
546	"""
547	# Load data
548	if os.path.isfile(path):
549	x, y, err = self.load(path)
550
551	# Get the data from the chosen data set and perform inversion
552	pr = Invertor()
553	pr.d_max = self.max_length
554	pr.alpha = self.alpha
555	pr.q_min = self.q_min
556	pr.q_max = self.q_max
557	pr.x = x
558	pr.y = y
559	pr.err = err
560
561	self.pr = pr
562	return True
563	return False
564
565	def perform_estimate(self):
566	from pr_thread import EstimatePr
567	from copy import deepcopy
568
569	wx.PostEvent(self.parent, StatusEvent(status=''))
570	# If a thread is already started, stop it
571	if self.estimation_thread != None and self.estimation_thread.isrunning():
572	self.estimation_thread.stop()
573
574	pr = self.pr.clone()
575	self.estimation_thread = EstimatePr(pr, self.nfunc, error_func=self._thread_error,
576	completefn = self._estimate_completed,
577	updatefn = None)
578	self.estimation_thread.queue()
579	self.estimation_thread.ready(2.5)
580
581	def perform_inversion(self):
582
583	# Time estimate
584	#estimated = self.elapsedself.nfunc*2
585	message = "Computation time may take up to %g seconds" % self.elapsed
586	wx.PostEvent(self.parent, StatusEvent(status=message))
587
588	# Start inversion thread
589	self.start_thread()
590	return
591
592	out, cov = self.pr.lstsq(self.nfunc)
593
594	# Save useful info
595	self.elapsed = self.pr.elapsed
596
597	for i in range(len(out)):
598	try:
599	print "%d: %g +- %g" % (i, out[i], math.sqrt(math.fabs(cov[i][i])))
600	except:
601	print "%d: %g +- ?" % (i, out[i])
602
603
604
605	# Make a plot of I(q) data
606	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
607	new_plot.name = "I_{obs}(q)"
608	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
609	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
610	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
611
612	# Show I(q) fit
613	self.show_iq(out, self.pr)
614
615	# Show P(r) fit
616	x_values, x_range = self.show_pr(out, self.pr, cov=cov)
617
618
619
620	def _on_context_inversion(self, event):
621	panel = event.GetEventObject()
622
623	from inversion_panel import InversionDlg
624
625	# If we have more than one displayed plot, make the user choose
626	if len(panel.plots)>1:
627	dialog = InversionDlg(None, -1, "P(r) Inversion", panel.plots, pars=False)
628	dialog.set_content(self.last_data, self.nfunc, self.alpha, self.max_length)
629	if dialog.ShowModal() == wx.ID_OK:
630	dataset = dialog.get_content()
631	dialog.Destroy()
632	else:
633	dialog.Destroy()
634	return
635	elif len(panel.plots)==1:
636	dataset = panel.plots.keys()[0]
637	else:
638	print "Error: No data is available"
639	return
640
641	# Store a reference to the current plottable
642	self.current_plottable = panel.plots[dataset]
643	self.control_panel.plotname = dataset
644	self.control_panel.nfunc = self.nfunc
645	self.control_panel.d_max = self.max_length
646	self.control_panel.alpha = self.alpha
647	self.parent.set_perspective(self.perspective)
648
649	def get_panels(self, parent):
650	"""
651	Create and return a list of panel objects
652	"""
653	from inversion_panel import InversionControl
654
655	self.parent = parent
656	self.control_panel = InversionControl(self.parent, -1, style=wx.RAISED_BORDER)
657	self.control_panel.set_manager(self)
658	self.control_panel.nfunc = self.nfunc
659	self.control_panel.d_max = self.max_length
660	self.control_panel.alpha = self.alpha
661
662	self.perspective = []
663	self.perspective.append(self.control_panel.window_name)
664	return [self.control_panel]
665
666	def get_perspective(self):
667	"""
668	Get the list of panel names for this perspective
669	"""
670	return self.perspective
671
672	def on_perspective(self, event):
673	"""
674	Call back function for the perspective menu item.
675	We notify the parent window that the perspective
676	has changed.
677	"""
678	self.parent.set_perspective(self.perspective)
679
680	def post_init(self):
681	"""
682	Post initialization call back to close the loose ends
683	[Somehow openGL needs this call]
684	"""
685	self.parent.set_perspective(self.perspective)
686

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

SasView

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

Download in other formats: