pr.py @ 9b85e93

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

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

SasView

source: sasview/prview/perspectives/pr/pr.py @ 9b85e93

Download in other formats: