pr.py @ 9a11937

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

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

SasView

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

Download in other formats: