pr.py @ aa4b8379

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

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

SasView

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

Download in other formats: