pr.py @ ea5551f

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 ea5551f was ea5551f, checked in by Mathieu Doucet <doucetm@…>, 16 years ago
Fixed problem with silent failure when loading file with Q=0 points.
Property mode set to `100644`
File size: 35.5 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	IQ_FIT_LABEL = "I_{fit}(q)"
18	IQ_SMEARED_LABEL = "I_{smeared}(q)"
19
20	import wx.lib
21	(NewPrFileEvent, EVT_PR_FILE) = wx.lib.newevent.NewEvent()
22
23
24	class Plugin:
25
26	DEFAULT_ALPHA = 0.0001
27	DEFAULT_NFUNC = 10
28	DEFAULT_DMAX = 140.0
29
30	def __init__(self):
31	## Plug-in name
32	self.sub_menu = "Pr inversion"
33
34	## Reference to the parent window
35	self.parent = None
36
37	## Simulation window manager
38	self.simview = None
39
40	## List of panels for the simulation perspective (names)
41	self.perspective = []
42
43	## State data
44	self.alpha = self.DEFAULT_ALPHA
45	self.nfunc = self.DEFAULT_NFUNC
46	self.max_length = self.DEFAULT_DMAX
47	self.q_min = None
48	self.q_max = None
49	self.has_bck = False
50	self.slit_height = 0
51	self.slit_width = 0
52	## Remember last plottable processed
53	self.last_data = "sphere_60_q0_2.txt"
54	self._current_file_data = None
55	## Time elapsed for last computation [sec]
56	# Start with a good default
57	self.elapsed = 0.022
58	self.iq_data_shown = False
59
60	## Current invertor
61	self.invertor = None
62	self.pr = None
63	# Copy of the last result in case we need to display it.
64	self._last_pr = None
65	self._last_out = None
66	self._last_cov = None
67	## Calculation thread
68	self.calc_thread = None
69	## Estimation thread
70	self.estimation_thread = None
71	## Result panel
72	self.control_panel = None
73	## Currently views plottable
74	self.current_plottable = None
75	## Number of P(r) points to display on the output plot
76	self._pr_npts = 51
77	## Flag to let the plug-in know that it is running standalone
78	self.standalone = True
79	self._normalize_output = False
80	self._scale_output_unity = False
81
82	# Log startup
83	logging.info("Pr(r) plug-in started")
84
85
86
87	def populate_menu(self, id, owner):
88	"""
89	Create a menu for the plug-in
90	"""
91	return []
92
93	def help(self, evt):
94	"""
95	Show a general help dialog.
96	TODO: replace the text with a nice image
97	"""
98	from inversion_panel import HelpDialog
99	dialog = HelpDialog(None, -1)
100	if dialog.ShowModal() == wx.ID_OK:
101	dialog.Destroy()
102	else:
103	dialog.Destroy()
104
105	def _fit_pr(self, evt):
106	from sans.pr.invertor import Invertor
107	import numpy
108	import pylab
109	import math
110	from sans.guicomm.events import NewPlotEvent
111	from sans.guitools.plottables import Data1D, Theory1D
112
113	# Generate P(r) for sphere
114	radius = 60.0
115	d_max = 2*radius
116
117
118	r = pylab.arange(0.01, d_max, d_max/51.0)
119	M = len(r)
120	y = numpy.zeros(M)
121	pr_err = numpy.zeros(M)
122
123	sum = 0.0
124	for j in range(M):
125	value = self.pr_theory(r[j], radius)
126	sum += value
127	y[j] = value
128	pr_err[j] = math.sqrt(y[j])
129
130
131	y = y/sum*d_max/len(r)
132
133
134
135	# Perform fit
136	pr = Invertor()
137	pr.d_max = d_max
138	pr.alpha = 0
139	pr.x = r
140	pr.y = y
141	pr.err = pr_err
142	out, cov = pr.pr_fit()
143	for i in range(len(out)):
144	print "%g +- %g" % (out[i], math.sqrt(cov[i][i]))
145
146
147	# Show input P(r)
148	new_plot = Data1D(pr.x, pr.y, dy=pr.err)
149	new_plot.name = "P_{obs}(r)"
150	new_plot.xaxis("\\rm{r}", 'A')
151	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
152	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Pr"))
153
154	# Show P(r) fit
155	self.show_pr(out, pr)
156
157	# Show I(q) fit
158	q = pylab.arange(0.001, 0.1, 0.01/51.0)
159	self.show_iq(out, pr, q)
160
161
162	def show_shpere(self, x, radius=70.0, x_range=70.0):
163	import numpy
164	import pylab
165	import math
166	from sans.guicomm.events import NewPlotEvent
167	from sans.guitools.plottables import Data1D, Theory1D
168	# Show P(r)
169	y_true = numpy.zeros(len(x))
170
171	sum_true = 0.0
172	for i in range(len(x)):
173	y_true[i] = self.pr_theory(x[i], radius)
174	sum_true += y_true[i]
175
176	y_true = y_true/sum_true*x_range/len(x)
177
178	# Show the theory P(r)
179	new_plot = Theory1D(x, y_true)
180	new_plot.name = "P_{true}(r)"
181	new_plot.xaxis("\\rm{r}", 'A')
182	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
183
184
185	#Put this call in plottables/guitools
186	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Sphere P(r)"))
187
188	def get_npts(self):
189	"""
190	Returns the number of points in the I(q) data
191	"""
192	try:
193	return len(self.pr.x)
194	except:
195	return 0
196
197	def show_iq(self, out, pr, q=None):
198	import numpy
199	import pylab
200	import math
201	from sans.guicomm.events import NewPlotEvent
202	from sans.guitools.plottables import Data1D, Theory1D
203
204	qtemp = pr.x
205	if not q==None:
206	qtemp = q
207
208	# Make a plot
209	maxq = -1
210	for q_i in qtemp:
211	if q_i>maxq:
212	maxq=q_i
213
214	minq = 0.001
215
216	# Check for user min/max
217	if not pr.q_min==None:
218	minq = pr.q_min
219	if not pr.q_max==None:
220	maxq = pr.q_max
221
222	x = pylab.arange(minq, maxq, maxq/301.0)
223	y = numpy.zeros(len(x))
224	err = numpy.zeros(len(x))
225	for i in range(len(x)):
226	value = pr.iq(out, x[i])
227	y[i] = value
228	try:
229	err[i] = math.sqrt(math.fabs(value))
230	except:
231	err[i] = 1.0
232	print "Error getting error", value, x[i]
233
234	new_plot = Theory1D(x, y)
235	new_plot.name = IQ_FIT_LABEL
236	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
237	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
238	#new_plot.group_id = "test group"
239	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="I(q)"))
240
241	# If we have used slit smearing, plot the smeared I(q) too
242	if pr.slit_width>0 or pr.slit_height>0:
243	x = pylab.arange(minq, maxq, maxq/301.0)
244	y = numpy.zeros(len(x))
245	err = numpy.zeros(len(x))
246	for i in range(len(x)):
247	value = pr.iq_smeared(out, x[i])
248	y[i] = value
249	try:
250	err[i] = math.sqrt(math.fabs(value))
251	except:
252	err[i] = 1.0
253	print "Error getting error", value, x[i]
254
255	new_plot = Theory1D(x, y)
256	new_plot.name = IQ_SMEARED_LABEL
257	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
258	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
259	#new_plot.group_id = "test group"
260	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="I(q)"))
261
262
263	def _on_pr_npts(self, evt):
264	"""
265	Redisplay P(r) with a different number of points
266	"""
267	from inversion_panel import PrDistDialog
268	dialog = PrDistDialog(None, -1)
269	dialog.set_content(self._pr_npts)
270	if dialog.ShowModal() == wx.ID_OK:
271	self._pr_npts= dialog.get_content()
272	dialog.Destroy()
273	self.show_pr(self._last_out, self._last_pr, self._last_cov)
274	else:
275	dialog.Destroy()
276
277
278	def show_pr(self, out, pr, cov=None):
279	import numpy
280	import pylab
281	import math
282	from sans.guicomm.events import NewPlotEvent
283	from sans.guitools.plottables import Data1D, Theory1D
284
285	# Show P(r)
286	x = pylab.arange(0.0, pr.d_max, pr.d_max/self._pr_npts)
287
288	y = numpy.zeros(len(x))
289	dy = numpy.zeros(len(x))
290	y_true = numpy.zeros(len(x))
291
292	sum = 0.0
293	pmax = 0.0
294	cov2 = numpy.ascontiguousarray(cov)
295
296	for i in range(len(x)):
297	if cov2==None:
298	value = pr.pr(out, x[i])
299	else:
300	(value, dy[i]) = pr.pr_err(out, cov2, x[i])
301	sum += value*pr.d_max/len(x)
302
303	# keep track of the maximum P(r) value
304	if value>pmax:
305	pmax = value
306
307	y[i] = value
308
309	if self._normalize_output==True:
310	y = y/sum
311	dy = dy/sum
312	elif self._scale_output_unity==True:
313	y = y/pmax
314	dy = dy/pmax
315
316	if cov2==None:
317	new_plot = Theory1D(x, y)
318	else:
319	new_plot = Data1D(x, y, dy=dy)
320	new_plot.name = "P_{fit}(r)"
321	new_plot.xaxis("\\rm{r}", 'A')
322	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
323
324	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="P(r) fit"))
325
326	return x, pr.d_max
327
328
329	def choose_file(self):
330	"""
331
332	"""
333	#TODO: this should be in a common module
334	return self.parent.choose_file()
335
336
337	def load(self, path):
338	"""
339	Load data. This will eventually be replaced
340	by our standard DataLoader class.
341	"""
342	#TODO: use DataLoader
343	class FileData:
344	x = None
345	y = None
346	err = None
347	path = None
348
349	def __init__(self, path):
350	self.path = path
351
352	self._current_file_data = FileData(path)
353	basename = os.path.basename(path)
354	root, ext = os.path.splitext(basename)
355	if ext.lower()=='.abs':
356	x, y, err = self.load_abs(path)
357
358	else:
359	x, y, err = self.load_columns(path)
360
361	self._current_file_data.x = x
362	self._current_file_data.y = y
363	self._current_file_data.err = err
364	return x, y, err
365
366	def load_columns(self, path = "sphere_60_q0_2.txt"):
367	"""
368	Load 2- or 3- column ascii
369	"""
370	import numpy, math, sys
371	# Read the data from the data file
372	data_x = numpy.zeros(0)
373	data_y = numpy.zeros(0)
374	data_err = numpy.zeros(0)
375	scale = None
376	min_err = 0.0
377	if not path == None:
378	input_f = open(path,'r')
379	buff = input_f.read()
380	lines = buff.split('\n')
381	for line in lines:
382	try:
383	toks = line.split()
384	x = float(toks[0])
385	y = float(toks[1])
386	if len(toks)>2:
387	err = float(toks[2])
388	else:
389	if scale==None:
390	scale = 0.05*math.sqrt(y)
391	#scale = 0.05/math.sqrt(y)
392	min_err = 0.01*y
393	err = scale*math.sqrt(y)+min_err
394	#err = 0
395
396	data_x = numpy.append(data_x, x)
397	data_y = numpy.append(data_y, y)
398	data_err = numpy.append(data_err, err)
399	except:
400	pass
401
402	if not scale==None:
403	message = "The loaded file had no error bars, statistical errors are assumed."
404	wx.PostEvent(self.parent, StatusEvent(status=message))
405	else:
406	wx.PostEvent(self.parent, StatusEvent(status=''))
407
408	return data_x, data_y, data_err
409
410	def load_abs(self, path):
411	"""
412	Load an IGOR .ABS reduced file
413	@param path: file path
414	@return: x, y, err vectors
415	"""
416	import numpy, math, sys
417	# Read the data from the data file
418	data_x = numpy.zeros(0)
419	data_y = numpy.zeros(0)
420	data_err = numpy.zeros(0)
421	scale = None
422	min_err = 0.0
423
424	data_started = False
425	if not path == None:
426	input_f = open(path,'r')
427	buff = input_f.read()
428	lines = buff.split('\n')
429	for line in lines:
430	if data_started==True:
431	try:
432	toks = line.split()
433	x = float(toks[0])
434	y = float(toks[1])
435	if len(toks)>2:
436	err = float(toks[2])
437	else:
438	if scale==None:
439	scale = 0.05*math.sqrt(y)
440	#scale = 0.05/math.sqrt(y)
441	min_err = 0.01*y
442	err = scale*math.sqrt(y)+min_err
443	#err = 0
444
445	data_x = numpy.append(data_x, x)
446	data_y = numpy.append(data_y, y)
447	data_err = numpy.append(data_err, err)
448	except:
449	pass
450	elif line.find("The 6 columns")>=0:
451	data_started = True
452
453	if not scale==None:
454	message = "The loaded file had no error bars, statistical errors are assumed."
455	wx.PostEvent(self.parent, StatusEvent(status=message))
456	else:
457	wx.PostEvent(self.parent, StatusEvent(status=''))
458
459	return data_x, data_y, data_err
460
461
462
463	def pr_theory(self, r, R):
464	"""
465
466	"""
467	if r<=2*R:
468	return 12.0* ((0.5r/R)2) ((1.0-0.5r/R)2) ( 2.0 + 0.5*r/R )
469	else:
470	return 0.0
471
472	def get_context_menu(self, graph=None):
473	"""
474	Get the context menu items available for P(r)
475	@param graph: the Graph object to which we attach the context menu
476	@return: a list of menu items with call-back function
477	"""
478	# Look whether this Graph contains P(r) data
479	#if graph.selected_plottable==IQ_DATA_LABEL:
480	for item in graph.plottables:
481	if item.name==PR_FIT_LABEL:
482	m_list = [["Add P(r) data", "Load a data file and display it on this plot", self._on_add_data],
483	["Change number of P(r) points", "Change the number of points on the P(r) output", self._on_pr_npts]]
484
485	if self._scale_output_unity==True or self._normalize_output==True:
486	m_list.append(["Disable P(r) scaling",
487	"Let the output P(r) keep the scale of the data",
488	self._on_disable_scaling])
489
490	if self._scale_output_unity==False:
491	m_list.append(["Scale P_max(r) to unity",
492	"Scale P(r) so that its maximum is 1",
493	self._on_scale_unity])
494
495	if self._normalize_output==False:
496	m_list.append(["Normalize P(r) to unity",
497	"Normalize the integral of P(r) to 1",
498	self._on_normalize])
499
500	return m_list
501	#return [["Add P(r) data", "Load a data file and display it on this plot", self._on_add_data],
502	# ["Change number of P(r) points", "Change the number of points on the P(r) output", self._on_pr_npts]]
503
504	elif item.name==graph.selected_plottable:
505	return [["Compute P(r)", "Compute P(r) from distribution", self._on_context_inversion]]
506
507	return []
508
509	def _on_disable_scaling(self, evt):
510	"""
511	Disable P(r) scaling
512	@param evt: Menu event
513	"""
514	self._normalize_output = False
515	self._scale_output_unity = False
516	self.show_pr(self._last_out, self._last_pr, self._last_cov)
517
518	def _on_normalize(self, evt):
519	"""
520	Switch normalization ON/OFF
521	@param evt: Menu event
522	"""
523	self._normalize_output = True
524	self._scale_output_unity = False
525
526	self.show_pr(self._last_out, self._last_pr, self._last_cov)
527
528	def _on_scale_unity(self, evt):
529	"""
530	Switch normalization ON/OFF
531	@param evt: Menu event
532	"""
533	self._scale_output_unity = True
534	self._normalize_output = False
535
536	self.show_pr(self._last_out, self._last_pr, self._last_cov)
537
538	def _on_add_data(self, evt):
539	"""
540	Add a data curve to the plot
541	WARNING: this will be removed once guiframe.plotting has its full functionality
542	"""
543	path = self.choose_file()
544	if path==None:
545	return
546
547	x, y, err = self.parent.load_ascii_1D(path)
548
549	#new_plot = Data1D(x, y, dy=err)
550	new_plot = Theory1D(x, y)
551	new_plot.name = "P_{loaded}(r)"
552	new_plot.xaxis("\\rm{r}", 'A')
553	new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
554
555	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="P(r) fit"))
556
557
558
559	def start_thread(self):
560	from pr_thread import CalcPr
561	from copy import deepcopy
562
563	# If a thread is already started, stop it
564	if self.calc_thread != None and self.calc_thread.isrunning():
565	self.calc_thread.stop()
566
567	pr = self.pr.clone()
568	self.calc_thread = CalcPr(pr, self.nfunc, error_func=self._thread_error, completefn=self._completed, updatefn=None)
569	self.calc_thread.queue()
570	self.calc_thread.ready(2.5)
571
572	def _thread_error(self, error):
573	wx.PostEvent(self.parent, StatusEvent(status=error))
574
575	def _estimate_completed(self, alpha, message, elapsed):
576	"""
577	Parameter estimation completed,
578	display the results to the user
579	@param alpha: estimated best alpha
580	@param elapsed: computation time
581	"""
582	# Save useful info
583	self.elapsed = elapsed
584	self.control_panel.alpha_estimate = alpha
585	if not message==None:
586	wx.PostEvent(self.parent, StatusEvent(status=str(message)))
587
588	self.perform_estimateNT()
589
590
591
592	def _estimateNT_completed(self, nterms, alpha, message, elapsed):
593	"""
594	Parameter estimation completed,
595	display the results to the user
596	@param alpha: estimated best alpha
597	@param nterms: estimated number of terms
598	@param elapsed: computation time
599	"""
600	# Save useful info
601	self.elapsed = elapsed
602	self.control_panel.nterms_estimate = nterms
603	self.control_panel.alpha_estimate = alpha
604	if not message==None:
605	wx.PostEvent(self.parent, StatusEvent(status=str(message)))
606
607	def _completed(self, out, cov, pr, elapsed):
608	"""
609	Method called with the results when the inversion
610	is done
611
612	@param out: output coefficient for the base functions
613	@param cov: covariance matrix
614	@param pr: Invertor instance
615	@param elapsed: time spent computing
616	"""
617	from copy import deepcopy
618	# Save useful info
619	self.elapsed = elapsed
620	# Keep a copy of the last result
621	self._last_pr = pr.clone()
622	self._last_out = out
623	self._last_cov = cov
624
625	# Save Pr invertor
626	self.pr = pr
627
628	#message = "Computation completed in %g seconds [chi2=%g]" % (elapsed, pr.chi2)
629	#wx.PostEvent(self.parent, StatusEvent(status=message))
630
631	cov = numpy.ascontiguousarray(cov)
632
633	# Show result on control panel
634	self.control_panel.chi2 = pr.chi2
635	self.control_panel.elapsed = elapsed
636	self.control_panel.oscillation = pr.oscillations(out)
637	#print "OSCILL", pr.oscillations(out)
638	#print "PEAKS:", pr.get_peaks(out)
639	self.control_panel.positive = pr.get_positive(out)
640	self.control_panel.pos_err = pr.get_pos_err(out, cov)
641	self.control_panel.rg = pr.rg(out)
642	self.control_panel.iq0 = pr.iq0(out)
643	self.control_panel.bck = pr.background
644
645	if False:
646	for i in range(len(out)):
647	try:
648	print "%d: %g +- %g" % (i, out[i], math.sqrt(math.fabs(cov[i][i])))
649	except:
650	print sys.exc_value
651	print "%d: %g +- ?" % (i, out[i])
652
653	# Make a plot of I(q) data
654	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
655	new_plot.name = "I_{obs}(q)"
656	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
657	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
658	#new_plot.group_id = "test group"
659	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
660
661	# Show I(q) fit
662	self.show_iq(out, self.pr)
663
664	# Show P(r) fit
665	x_values, x_range = self.show_pr(out, self.pr, cov)
666
667	# Popup result panel
668	#result_panel = InversionResults(self.parent, -1, style=wx.RAISED_BORDER)
669
670	def show_data(self, path=None, reset=False):
671	"""
672	Show data read from a file
673	@param path: file path
674	@param reset: if True all other plottables will be cleared
675	"""
676
677
678	if path is not None:
679	try:
680	pr = self._create_file_pr(path)
681	self.pr = pr
682	except:
683	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
684	return
685
686	# Make a plot of I(q) data
687	if self.pr.err==None:
688	new_plot = Theory1D(self.pr.x, self.pr.y)
689	else:
690	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
691	new_plot.name = "I_{obs}(q)"
692	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
693	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
694	new_plot.interactive = True
695	#new_plot.group_id = "test group"
696	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="I(q)", reset=reset))
697
698	# Get Q range
699	self.control_panel.q_min = self.pr.x.min()
700	self.control_panel.q_max = self.pr.x.max()
701
702
703
704	def setup_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None,
705	bck=False, height=0, width=0):
706	self.alpha = alpha
707	self.nfunc = nfunc
708	self.max_length = d_max
709	self.q_min = q_min
710	self.q_max = q_max
711	self.has_bck = bck
712	self.slit_height = height
713	self.slit_width = width
714
715	try:
716	pr = self._create_plot_pr()
717	if not pr==None:
718	self.pr = pr
719	self.perform_inversion()
720	except:
721	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
722
723	def estimate_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None,
724	bck=False, height=0, width=0):
725	self.alpha = alpha
726	self.nfunc = nfunc
727	self.max_length = d_max
728	self.q_min = q_min
729	self.q_max = q_max
730	self.has_bck = bck
731	self.slit_height = height
732	self.slit_width = width
733
734	try:
735	pr = self._create_plot_pr()
736	if not pr==None:
737	self.pr = pr
738	self.perform_estimate()
739	except:
740	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
741
742	def _create_plot_pr(self, estimate=False):
743	"""
744	Create and prepare invertor instance from
745	a plottable data set.
746	@param path: path of the file to read in
747	"""
748	# Get the data from the chosen data set and perform inversion
749	pr = Invertor()
750	pr.d_max = self.max_length
751	pr.alpha = self.alpha
752	pr.q_min = self.q_min
753	pr.q_max = self.q_max
754	pr.x = self.current_plottable.x
755	pr.y = self.current_plottable.y
756	pr.has_bck = self.has_bck
757	pr.slit_height = self.slit_height
758	pr.slit_width = self.slit_width
759
760	# Fill in errors if none were provided
761	if self.current_plottable.dy == None:
762	print "no error", self.current_plottable.name
763	y = numpy.zeros(len(pr.y))
764	for i in range(len(pr.y)):
765	y[i] = math.sqrt(pr.y[i])
766	pr.err = y
767	else:
768	pr.err = self.current_plottable.dy
769
770	#self.pr = pr
771	return pr
772
773
774	def setup_file_inversion(self, alpha, nfunc, d_max, path, q_min=None, q_max=None,
775	bck=False, height=0, width=0):
776	self.alpha = alpha
777	self.nfunc = nfunc
778	self.max_length = d_max
779	self.q_min = q_min
780	self.q_max = q_max
781	self.has_bck = bck
782	self.slit_height = height
783	self.slit_width = width
784
785	try:
786	pr = self._create_file_pr(path)
787	if not pr==None:
788	self.pr = pr
789	self.perform_inversion()
790	except:
791	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
792
793	def estimate_file_inversion(self, alpha, nfunc, d_max, path, q_min=None, q_max=None,
794	bck=False, height=0, width=0):
795	self.alpha = alpha
796	self.nfunc = nfunc
797	self.max_length = d_max
798	self.q_min = q_min
799	self.q_max = q_max
800	self.has_bck = bck
801	self.slit_height = height
802	self.slit_width = width
803
804	try:
805	pr = self._create_file_pr(path)
806	if not pr==None:
807	self.pr = pr
808	self.perform_estimate()
809	except:
810	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
811
812
813	def _create_file_pr(self, path):
814	"""
815	Create and prepare invertor instance from
816	a file data set.
817	@param path: path of the file to read in
818	"""
819	from sans.guiframe.data_loader import load_ascii_1D
820	import numpy
821
822	# Load data
823	if os.path.isfile(path):
824
825	if self._current_file_data is not None \
826	and self._current_file_data.path==path:
827	x = self._current_file_data.x
828	y = self._current_file_data.y
829	err = self._current_file_data.err
830	else:
831	x, y, err = self.load(path)
832	#x, y, err = load_ascii_1D(path)
833
834	# If we have not errors, add statistical errors
835	if err==None:
836	err = numpy.zeros(len(y))
837	for i in range(len(y)):
838	err[i] = math.sqrt( math.fabs(y[i]) )
839	message = "The loaded file had no error bars, statistical errors are assumed."
840	wx.PostEvent(self.parent, StatusEvent(status=message))
841	else:
842	wx.PostEvent(self.parent, StatusEvent(status=''))
843
844	# Get the data from the chosen data set and perform inversion
845	pr = Invertor()
846	pr.d_max = self.max_length
847	pr.alpha = self.alpha
848	pr.q_min = self.q_min
849	pr.q_max = self.q_max
850	pr.x = x
851	pr.y = y
852	pr.err = err
853	pr.has_bck = self.has_bck
854	pr.slit_height = self.slit_height
855	pr.slit_width = self.slit_width
856	return pr
857	return None
858
859	def perform_estimate(self):
860	from pr_thread import EstimatePr
861	from copy import deepcopy
862
863	# If a thread is already started, stop it
864	if self.estimation_thread != None and self.estimation_thread.isrunning():
865	self.estimation_thread.stop()
866
867	pr = self.pr.clone()
868	self.estimation_thread = EstimatePr(pr, self.nfunc, error_func=self._thread_error,
869	completefn = self._estimate_completed,
870	updatefn = None)
871	self.estimation_thread.queue()
872	self.estimation_thread.ready(2.5)
873
874	def perform_estimateNT(self):
875	from pr_thread import EstimateNT
876	from copy import deepcopy
877
878	# If a thread is already started, stop it
879	if self.estimation_thread != None and self.estimation_thread.isrunning():
880	self.estimation_thread.stop()
881
882	pr = self.pr.clone()
883	# Skip the slit settings for the estimation
884	# It slows down the application and it doesn't change the estimates
885	pr.slit_height = 0.0
886	pr.slit_width = 0.0
887	self.estimation_thread = EstimateNT(pr, self.nfunc, error_func=self._thread_error,
888	completefn = self._estimateNT_completed,
889	updatefn = None)
890	self.estimation_thread.queue()
891	self.estimation_thread.ready(2.5)
892
893	def perform_inversion(self):
894
895	# Time estimate
896	#estimated = self.elapsedself.nfunc*2
897	#message = "Computation time may take up to %g seconds" % self.elapsed
898	#wx.PostEvent(self.parent, StatusEvent(status=message))
899
900	# Start inversion thread
901	self.start_thread()
902	return
903
904	out, cov = self.pr.lstsq(self.nfunc)
905
906	# Save useful info
907	self.elapsed = self.pr.elapsed
908
909	for i in range(len(out)):
910	try:
911	print "%d: %g +- %g" % (i, out[i], math.sqrt(math.fabs(cov[i][i])))
912	except:
913	print "%d: %g +- ?" % (i, out[i])
914
915
916
917	# Make a plot of I(q) data
918	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
919	new_plot.name = "I_{obs}(q)"
920	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
921	new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
922	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
923
924	# Show I(q) fit
925	self.show_iq(out, self.pr)
926
927	# Show P(r) fit
928	x_values, x_range = self.show_pr(out, self.pr, cov=cov)
929
930
931
932	def _on_context_inversion(self, event):
933	panel = event.GetEventObject()
934
935	from inversion_panel import InversionDlg
936
937	# If we have more than one displayed plot, make the user choose
938	if len(panel.plots)>1 and panel.graph.selected_plottable in panel.plots:
939	dataset = panel.graph.selected_plottable
940	if False:
941	dialog = InversionDlg(None, -1, "P(r) Inversion", panel.plots, pars=False)
942	dialog.set_content(self.last_data, self.nfunc, self.alpha, self.max_length)
943	if dialog.ShowModal() == wx.ID_OK:
944	dataset = dialog.get_content()
945	dialog.Destroy()
946	else:
947	dialog.Destroy()
948	return
949	elif len(panel.plots)==1:
950	dataset = panel.plots.keys()[0]
951	else:
952	print "Error: No data is available"
953	return
954
955	# Store a reference to the current plottable
956	# If we have a suggested value, use it.
957	try:
958	estimate = float(self.control_panel.alpha_estimate)
959	self.control_panel.alpha = estimate
960	except:
961	self.control_panel.alpha = self.alpha
962	print "No estimate yet"
963	pass
964	try:
965	estimate = int(self.control_panel.nterms_estimate)
966	self.control_panel.nfunc = estimate
967	except:
968	self.control_panel.nfunc = self.nfunc
969	print "No estimate yet"
970	pass
971
972	self.current_plottable = panel.plots[dataset]
973	self.control_panel.plotname = dataset
974	#self.control_panel.nfunc = self.nfunc
975	self.control_panel.d_max = self.max_length
976	self.parent.set_perspective(self.perspective)
977	self.control_panel._on_invert(None)
978
979	def get_panels(self, parent):
980	"""
981	Create and return a list of panel objects
982	"""
983	from inversion_panel import InversionControl
984
985	self.parent = parent
986	self.control_panel = InversionControl(self.parent, -1,
987	style=wx.RAISED_BORDER,
988	standalone=self.standalone)
989	self.control_panel.set_manager(self)
990	self.control_panel.nfunc = self.nfunc
991	self.control_panel.d_max = self.max_length
992	self.control_panel.alpha = self.alpha
993
994	self.perspective = []
995	self.perspective.append(self.control_panel.window_name)
996
997	self.parent.Bind(EVT_PR_FILE, self._on_new_file)
998
999	return [self.control_panel]
1000
1001	def _on_new_file(self, evt):
1002	"""
1003	Called when the application manager posted an
1004	EVT_PR_FILE event. Just prompt the control
1005	panel to load a new data file.
1006	"""
1007	self.control_panel._change_file(None)
1008
1009	def get_perspective(self):
1010	"""
1011	Get the list of panel names for this perspective
1012	"""
1013	return self.perspective
1014
1015	def on_perspective(self, event):
1016	"""
1017	Call back function for the perspective menu item.
1018	We notify the parent window that the perspective
1019	has changed.
1020	"""
1021	self.parent.set_perspective(self.perspective)
1022
1023	def post_init(self):
1024	"""
1025	Post initialization call back to close the loose ends
1026	[Somehow openGL needs this call]
1027	"""
1028	self.parent.set_perspective(self.perspective)
1029
1030	if __name__ == "__main__":
1031	i = Plugin()
1032	print i.perform_estimateNT()
1033
1034
1035
1036

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

SasView

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

Download in other formats: