invertor.py @ 4a0536a

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 4a0536a was f71287f4, checked in by Mathieu Doucet <doucetm@…>, 16 years ago
Allow user to set q_min/q_max
Property mode set to `100644`
File size: 16.1 KB

Line
1	from sans.pr.core.pr_inversion import Cinvertor
2	import numpy
3	import sys
4
5	class Invertor(Cinvertor):
6
7	## Chisqr of the last computation
8	chi2 = 0
9	## Time elapsed for last computation
10	elapsed = 0
11	## Alpha to get the reg term the same size as the signal
12	suggested_alpha = 0
13	## Last number of base functions used
14	nfunc = 0
15	## Last output values
16	out = None
17	## Last errors on output values
18	cov = None
19
20	def __init__(self):
21	Cinvertor.__init__(self)
22
23	def __setattr__(self, name, value):
24	"""
25	Set the value of an attribute.
26	Access the parent class methods for
27	x, y, err and d_max.
28	"""
29	if name=='x':
30	if 0.0 in value:
31	raise ValueError, "Invertor: one of your q-values is zero. Delete that entry before proceeding"
32	return self.set_x(value)
33	elif name=='y':
34	return self.set_y(value)
35	elif name=='err':
36	return self.set_err(value)
37	elif name=='d_max':
38	return self.set_dmax(value)
39	elif name=='q_min':
40	if value==None:
41	return self.set_qmin(-1.0)
42	return self.set_qmin(value)
43	elif name=='q_max':
44	if value==None:
45	return self.set_qmax(-1.0)
46	return self.set_qmax(value)
47	elif name=='alpha':
48	return self.set_alpha(value)
49
50	return Cinvertor.__setattr__(self, name, value)
51
52	def __getattr__(self, name):
53	"""
54	Return the value of an attribute
55	For the moment x, y, err and d_max are write-only
56	TODO: change that!
57	"""
58	import numpy
59	if name=='x':
60	out = numpy.ones(self.get_nx())
61	self.get_x(out)
62	return out
63	elif name=='y':
64	out = numpy.ones(self.get_ny())
65	self.get_y(out)
66	return out
67	elif name=='err':
68	out = numpy.ones(self.get_nerr())
69	self.get_err(out)
70	return out
71	elif name=='d_max':
72	return self.get_dmax()
73	elif name=='q_min':
74	qmin = self.get_qmin()
75	if qmin<0:
76	return None
77	return qmin
78	elif name=='q_max':
79	qmax = self.get_qmax()
80	if qmax<0:
81	return None
82	return qmax
83	elif name=='alpha':
84	return self.get_alpha()
85	elif name in self.__dict__:
86	return self.__dict__[name]
87	return None
88
89	def clone(self):
90	"""
91	Return a clone of this instance
92	"""
93	invertor = Invertor()
94	invertor.chi2 = self.chi2
95	invertor.elapsed = self.elapsed
96	invertor.alpha = self.alpha
97	invertor.d_max = self.d_max
98	invertor.q_min = self.q_min
99	invertor.q_max = self.q_max
100
101	invertor.x = self.x
102	invertor.y = self.y
103	invertor.err = self.err
104
105	return invertor
106
107	def invert(self, nfunc=5):
108	"""
109	Perform inversion to P(r)
110	"""
111	from scipy import optimize
112	import time
113
114	self.nfunc = nfunc
115	# First, check that the current data is valid
116	if self.is_valid()<=0:
117	raise RuntimeError, "Invertor.invert: Data array are of different length"
118
119	p = numpy.ones(nfunc)
120	t_0 = time.time()
121	out, cov_x, info, mesg, success = optimize.leastsq(self.residuals, p, full_output=1, warning=True)
122
123	# Compute chi^2
124	res = self.residuals(out)
125	chisqr = 0
126	for i in range(len(res)):
127	chisqr += res[i]
128
129	self.chi2 = chisqr
130
131	# Store computation time
132	self.elapsed = time.time() - t_0
133
134	return out, cov_x
135
136	def pr_fit(self, nfunc=5):
137	"""
138	Perform inversion to P(r)
139	"""
140	from scipy import optimize
141
142	# First, check that the current data is valid
143	if self.is_valid()<=0:
144	raise RuntimeError, "Invertor.invert: Data arrays are of different length"
145
146	p = numpy.ones(nfunc)
147	t_0 = time.time()
148	out, cov_x, info, mesg, success = optimize.leastsq(self.pr_residuals, p, full_output=1, warning=True)
149
150	# Compute chi^2
151	res = self.pr_residuals(out)
152	chisqr = 0
153	for i in range(len(res)):
154	chisqr += res[i]
155
156	self.chisqr = chisqr
157
158	# Store computation time
159	self.elapsed = time.time() - t_0
160
161	return out, cov_x
162
163	def pr_err(self, c, c_cov, r):
164	import math
165	c_err = numpy.zeros(len(c))
166	for i in range(len(c)):
167	try:
168	c_err[i] = math.sqrt(math.fabs(c_cov[i][i]))
169	except:
170	import sys
171	print sys.exc_value
172	print "oups", c_cov[i][i]
173	c_err[i] = c[i]
174
175	return self.get_pr_err(c, c_err, r)
176
177	def _accept_q(self, q):
178	"""
179	Check q-value against user-defined range
180	"""
181	if not self.q_min==None and q<self.q_min:
182	return False
183	if not self.q_max==None and q>self.q_max:
184	return False
185	return True
186
187	def lstsq(self, nfunc=5):
188	#TODO: do this on the C side
189	#
190	# To make sure an array is contiguous:
191	# blah = numpy.ascontiguousarray(blah_original)
192	# ... before passing it to C
193
194	import math
195	from scipy.linalg.basic import lstsq
196
197	self.nfunc = nfunc
198	# a -- An M x N matrix.
199	# b -- An M x nrhs matrix or M vector.
200	npts = len(self.x)
201	nq = 20
202	sqrt_alpha = math.sqrt(self.alpha)
203
204	a = numpy.zeros([npts+nq, nfunc])
205	b = numpy.zeros(npts+nq)
206	err = numpy.zeros(nfunc)
207
208	for j in range(nfunc):
209	for i in range(npts):
210	if self._accept_q(self.x[i]):
211	a[i][j] = self.basefunc_ft(self.d_max, j+1, self.x[i])/self.err[i]
212	for i_q in range(nq):
213	r = self.d_max/nq*i_q
214	#a[i_q+npts][j] = sqrt_alpha * 1.0/nqself.d_max2.0math.fabs(math.sin(math.pi(j+1)r/self.d_max) + math.pi(j+1)r/self.d_max math.cos(math.pi(j+1)r/self.d_max))
215	a[i_q+npts][j] = sqrt_alpha * 1.0/nqself.d_max2.0(2.0math.pi(j+1)/self.d_maxmath.cos(math.pi(j+1)r/self.d_max) + math.pi*2(j+1)*2r/self.d_max*2 math.sin(math.pi(j+1)r/self.d_max))
216
217	for i in range(npts):
218	if self._accept_q(self.x[i]):
219	b[i] = self.y[i]/self.err[i]
220
221	c, chi2, rank, n = lstsq(a, b)
222	self.chi2 = chi2
223
224	at = numpy.transpose(a)
225	inv_cov = numpy.zeros([nfunc,nfunc])
226	for i in range(nfunc):
227	for j in range(nfunc):
228	inv_cov[i][j] = 0.0
229	for k in range(npts):
230	if self._accept_q(self.x[i]):
231	inv_cov[i][j] = at[i][k]*a[k][j]
232
233	# Compute the reg term size for the output
234	sum_sig = 0.0
235	sum_reg = 0.0
236	for j in range(nfunc):
237	for i in range(npts):
238	if self._accept_q(self.x[i]):
239	sum_sig += (a[i][j])**2
240	for i in range(nq):
241	sum_reg += (a[i+npts][j])**2
242
243	if math.fabs(self.alpha)>0:
244	new_alpha = sum_sig/(sum_reg/self.alpha)
245	else:
246	new_alpha = 0.0
247	self.suggested_alpha = new_alpha
248
249	try:
250	err = math.fabs(chi2/(npts-nfunc))* inv_cov
251	except:
252	print "Error estimating uncertainties"
253
254	# Keep a copy of the last output
255	self.out = c
256	self.cov = err
257
258	return c, err
259
260	def svd(self, nfunc=5):
261	import math, time
262	# Ac - b = 0
263
264	A = numpy.zeros([nfunc, nfunc])
265	y = numpy.zeros(nfunc)
266
267	t_0 = time.time()
268	for i in range(nfunc):
269	# A
270	for j in range(nfunc):
271	A[i][j] = 0.0
272	for k in range(len(self.x)):
273	err = self.err[k]
274	A[i][j] += 1.0/err/err*self.basefunc_ft(self.d_max, j+1, self.x[k]) \
275	*self.basefunc_ft(self.d_max, i+1, self.x[k]);
276	#print A[i][j]
277	#A[i][j] -= self.alpha(math.cos(math.pi(i+j)) - math.cos(math.pi*(i-j)));
278	if i==j:
279	A[i][j] += -1.0*self.alpha
280	elif i-j==1 or i-j==-1:
281	A[i][j] += 1.0*self.alpha
282	#print " ",A[i][j]
283	# y
284	y[i] = 0.0
285	for k in range(len(self.x)):
286	y[i] = self.y[k]/self.err[k]/self.err[k]*self.basefunc_ft(self.d_max, i+1, self.x[k])
287
288	print time.time()-t_0, 'secs'
289
290	# use numpy.pinv(A)
291	#inv_A = numpy.linalg.inv(A)
292	#c = y*inv_A
293	print y
294	c = numpy.linalg.solve(A, y)
295
296
297	print c
298
299	err = numpy.zeros(len(c))
300	return c, err
301
302	def estimate_alpha(self, nfunc):
303	"""
304	Returns a reasonable guess for the
305	regularization constant alpha
306
307	@return: alpha, message, elapsed
308
309	where alpha is the estimate for alpha,
310	message is a message for the user,
311	elapsed is the computation time
312	"""
313	import time
314	try:
315	pr = self.clone()
316
317	# T_0 for computation time
318	starttime = time.time()
319
320	# If the current alpha is zero, try
321	# another value
322	if pr.alpha<=0:
323	pr.alpha = 0.0001
324
325	# Perform inversion to find the largest alpha
326	out, cov = pr.lstsq(nfunc)
327	elapsed = time.time()-starttime
328	initial_alpha = pr.alpha
329	initial_peaks = pr.get_peaks(out)
330
331	# Try the inversion with the estimated alpha
332	pr.alpha = pr.suggested_alpha
333	out, cov = pr.lstsq(nfunc)
334
335	npeaks = pr.get_peaks(out)
336	# if more than one peak to start with
337	# just return the estimate
338	if npeaks>1:
339	message = "Your P(r) is not smooth, please check your inversion parameters"
340	return pr.suggested_alpha, message, elapsed
341	else:
342
343	# Look at smaller values
344	# We assume that for the suggested alpha, we have 1 peak
345	# if not, send a message to change parameters
346	alpha = pr.suggested_alpha
347	best_alpha = pr.suggested_alpha
348	found = False
349	for i in range(10):
350	pr.alpha = (0.33)*(i+1)alpha
351	out, cov = pr.lstsq(nfunc)
352
353	peaks = pr.get_peaks(out)
354	print pr.alpha, peaks
355	if peaks>1:
356	found = True
357	break
358	best_alpha = pr.alpha
359
360	# If we didn't find a turning point for alpha and
361	# the initial alpha already had only one peak,
362	# just return that
363	if not found and initial_peaks==1 and initial_alpha<best_alpha:
364	best_alpha = initial_alpha
365
366	# Check whether the size makes sense
367	message=''
368
369	if not found:
370	message = "None"
371	elif best_alpha>=0.5*pr.suggested_alpha:
372	# best alpha is too big, return a
373	# reasonable value
374	message = "The estimated alpha for your system is too large. "
375	message += "Try increasing your maximum distance."
376
377	return best_alpha, message, elapsed
378
379	except:
380	message = "Invertor.estimate_alpha: %s" % sys.exc_value
381	return 0, message, elapsed
382
383
384	def to_file(self, path, npts=100):
385	"""
386	Save the state to a file that will be readable
387	by SliceView.
388	@param path: path of the file to write
389	@param npts: number of P(r) points to be written
390	"""
391	import pylab
392
393	file = open(path, 'w')
394	file.write("#d_max=%g\n" % self.d_max)
395	file.write("#nfunc=%g\n" % self.nfunc)
396	file.write("#alpha=%g\n" % self.alpha)
397	file.write("#chi2=%g\n" % self.chi2)
398	file.write("#elapsed=%g\n" % self.elapsed)
399	file.write("#alpha_estimate=%g\n" % self.suggested_alpha)
400	if not self.out==None:
401	if len(self.out)==len(self.cov):
402	for i in range(len(self.out)):
403	file.write("#C_%i=%s+-%s\n" % (i, str(self.out[i]), str(self.cov[i][i])))
404	file.write("<r> <Pr> <dPr>\n")
405	r = pylab.arange(0.0, self.d_max, self.d_max/npts)
406
407	for r_i in r:
408	(value, err) = self.pr_err(self.out, self.cov, r_i)
409	file.write("%g %g %g\n" % (r_i, value, err))
410
411	file.close()
412
413	def from_file(self, path):
414	"""
415	Load the state of the Invertor from a file,
416	to be able to generate P(r) from a set of
417	parameters.
418	@param path: path of the file to load
419	"""
420	import os
421	import re
422	if os.path.isfile(path):
423	try:
424	fd = open(path, 'r')
425
426	buff = fd.read()
427	lines = buff.split('\n')
428	for line in lines:
429	if line.startswith('#d_max='):
430	toks = line.split('=')
431	self.d_max = float(toks[1])
432	elif line.startswith('#nfunc='):
433	toks = line.split('=')
434	self.nfunc = int(toks[1])
435	self.out = numpy.zeros(self.nfunc)
436	self.cov = numpy.zeros([self.nfunc, self.nfunc])
437	elif line.startswith('#alpha='):
438	toks = line.split('=')
439	self.alpha = float(toks[1])
440	elif line.startswith('#chi2='):
441	toks = line.split('=')
442	self.chi2 = float(toks[1])
443	elif line.startswith('#elapsed='):
444	toks = line.split('=')
445	self.elapsed = float(toks[1])
446	elif line.startswith('#alpha_estimate='):
447	toks = line.split('=')
448	self.suggested_alpha = float(toks[1])
449
450	# Now read in the parameters
451	elif line.startswith('#C_'):
452	toks = line.split('=')
453	p = re.compile('#C_([0-9]+)')
454	m = p.search(toks[0])
455	toks2 = toks[1].split('+-')
456	i = int(m.group(1))
457	self.out[i] = float(toks2[0])
458
459	self.cov[i][i] = float(toks2[1])
460
461	except:
462	raise RuntimeError, "Invertor.from_file: corrupted file\n%s" % sys.exc_value
463	else:
464	raise RuntimeError, "Invertor.from_file: '%s' is not a file" % str(path)
465
466
467
468
469	if __name__ == "__main__":
470	o = Invertor()
471
472
473
474
475

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

SasView

source: sasview/pr_inversion/invertor.py @ 4a0536a

Download in other formats: