1 | from sans.pr.core.pr_inversion import Cinvertor |
---|
2 | import numpy |
---|
3 | |
---|
4 | class Invertor(Cinvertor): |
---|
5 | |
---|
6 | ## Chisqr of the last computation |
---|
7 | chi2 = 0 |
---|
8 | ## Time elapsed for last computation |
---|
9 | elapsed = 0 |
---|
10 | ## Alpha to get the reg term the same size as the signal |
---|
11 | suggested_alpha = 0 |
---|
12 | |
---|
13 | def __init__(self): |
---|
14 | Cinvertor.__init__(self) |
---|
15 | |
---|
16 | def __setattr__(self, name, value): |
---|
17 | """ |
---|
18 | Set the value of an attribute. |
---|
19 | Access the parent class methods for |
---|
20 | x, y, err and d_max. |
---|
21 | """ |
---|
22 | if name=='x': |
---|
23 | if 0.0 in value: |
---|
24 | raise ValueError, "Invertor: one of your q-values is zero. Delete that entry before proceeding" |
---|
25 | return self.set_x(value) |
---|
26 | elif name=='y': |
---|
27 | return self.set_y(value) |
---|
28 | elif name=='err': |
---|
29 | return self.set_err(value) |
---|
30 | elif name=='d_max': |
---|
31 | return self.set_dmax(value) |
---|
32 | elif name=='alpha': |
---|
33 | return self.set_alpha(value) |
---|
34 | |
---|
35 | return Cinvertor.__setattr__(self, name, value) |
---|
36 | |
---|
37 | def __getattr__(self, name): |
---|
38 | """ |
---|
39 | Return the value of an attribute |
---|
40 | For the moment x, y, err and d_max are write-only |
---|
41 | TODO: change that! |
---|
42 | """ |
---|
43 | import numpy |
---|
44 | if name=='x': |
---|
45 | out = numpy.ones(self.get_nx()) |
---|
46 | self.get_x(out) |
---|
47 | return out |
---|
48 | elif name=='y': |
---|
49 | out = numpy.ones(self.get_ny()) |
---|
50 | self.get_y(out) |
---|
51 | return out |
---|
52 | elif name=='err': |
---|
53 | out = numpy.ones(self.get_nerr()) |
---|
54 | self.get_err(out) |
---|
55 | return out |
---|
56 | elif name=='d_max': |
---|
57 | return self.get_dmax() |
---|
58 | elif name=='alpha': |
---|
59 | return self.get_alpha() |
---|
60 | elif name in self.__dict__: |
---|
61 | return self.__dict__[name] |
---|
62 | return None |
---|
63 | |
---|
64 | def clone(self): |
---|
65 | """ |
---|
66 | Return a clone of this instance |
---|
67 | """ |
---|
68 | invertor = Invertor() |
---|
69 | invertor.chi2 = self.chi2 |
---|
70 | invertor.elapsed = self.elapsed |
---|
71 | invertor.alpha = self.alpha |
---|
72 | invertor.d_max = self.d_max |
---|
73 | |
---|
74 | invertor.x = self.x |
---|
75 | invertor.y = self.y |
---|
76 | invertor.err = self.err |
---|
77 | |
---|
78 | return invertor |
---|
79 | |
---|
80 | def invert(self, nfunc=5): |
---|
81 | """ |
---|
82 | Perform inversion to P(r) |
---|
83 | """ |
---|
84 | from scipy import optimize |
---|
85 | import time |
---|
86 | |
---|
87 | # First, check that the current data is valid |
---|
88 | if self.is_valid()<=0: |
---|
89 | raise RuntimeError, "Invertor.invert: Data array are of different length" |
---|
90 | |
---|
91 | p = numpy.ones(nfunc) |
---|
92 | t_0 = time.time() |
---|
93 | out, cov_x, info, mesg, success = optimize.leastsq(self.residuals, p, full_output=1, warning=True) |
---|
94 | |
---|
95 | # Compute chi^2 |
---|
96 | res = self.residuals(out) |
---|
97 | chisqr = 0 |
---|
98 | for i in range(len(res)): |
---|
99 | chisqr += res[i] |
---|
100 | |
---|
101 | self.chi2 = chisqr |
---|
102 | |
---|
103 | # Store computation time |
---|
104 | self.elapsed = time.time() - t_0 |
---|
105 | |
---|
106 | return out, cov_x |
---|
107 | |
---|
108 | def pr_fit(self, nfunc=5): |
---|
109 | """ |
---|
110 | Perform inversion to P(r) |
---|
111 | """ |
---|
112 | from scipy import optimize |
---|
113 | |
---|
114 | # First, check that the current data is valid |
---|
115 | if self.is_valid()<=0: |
---|
116 | raise RuntimeError, "Invertor.invert: Data arrays are of different length" |
---|
117 | |
---|
118 | p = numpy.ones(nfunc) |
---|
119 | t_0 = time.time() |
---|
120 | out, cov_x, info, mesg, success = optimize.leastsq(self.pr_residuals, p, full_output=1, warning=True) |
---|
121 | |
---|
122 | # Compute chi^2 |
---|
123 | res = self.pr_residuals(out) |
---|
124 | chisqr = 0 |
---|
125 | for i in range(len(res)): |
---|
126 | chisqr += res[i] |
---|
127 | |
---|
128 | self.chisqr = chisqr |
---|
129 | |
---|
130 | # Store computation time |
---|
131 | self.elapsed = time.time() - t_0 |
---|
132 | |
---|
133 | return out, cov_x |
---|
134 | |
---|
135 | def pr_err(self, c, c_cov, r): |
---|
136 | import math |
---|
137 | c_err = numpy.zeros(len(c)) |
---|
138 | for i in range(len(c)): |
---|
139 | try: |
---|
140 | c_err[i] = math.sqrt(math.fabs(c_cov[i][i])) |
---|
141 | except: |
---|
142 | import sys |
---|
143 | print sys.exc_value |
---|
144 | print "oups", c_cov[i][i] |
---|
145 | c_err[i] = c[i] |
---|
146 | |
---|
147 | return self.get_pr_err(c, c_err, r) |
---|
148 | |
---|
149 | def lstsq(self, nfunc=5): |
---|
150 | import math |
---|
151 | from scipy.linalg.basic import lstsq |
---|
152 | |
---|
153 | # a -- An M x N matrix. |
---|
154 | # b -- An M x nrhs matrix or M vector. |
---|
155 | npts = len(self.x) |
---|
156 | nq = 20 |
---|
157 | sqrt_alpha = math.sqrt(self.alpha) |
---|
158 | |
---|
159 | a = numpy.zeros([npts+nq, nfunc]) |
---|
160 | b = numpy.zeros(npts+nq) |
---|
161 | err = numpy.zeros(nfunc) |
---|
162 | |
---|
163 | for j in range(nfunc): |
---|
164 | for i in range(npts): |
---|
165 | a[i][j] = self.basefunc_ft(self.d_max, j+1, self.x[i])/self.err[i] |
---|
166 | for i_q in range(nq): |
---|
167 | r = self.d_max/nq*i_q |
---|
168 | #a[i_q+npts][j] = sqrt_alpha * 1.0/nq*self.d_max*2.0*math.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)) |
---|
169 | a[i_q+npts][j] = sqrt_alpha * 1.0/nq*self.d_max*2.0*(2.0*math.pi*(j+1)/self.d_max*math.cos(math.pi*(j+1)*r/self.d_max) + math.pi**2*(j+1)**2*r/self.d_max**2 * math.sin(math.pi*(j+1)*r/self.d_max)) |
---|
170 | |
---|
171 | for i in range(npts): |
---|
172 | b[i] = self.y[i]/self.err[i] |
---|
173 | |
---|
174 | c, chi2, rank, n = lstsq(a, b) |
---|
175 | self.chi2 = chi2 |
---|
176 | |
---|
177 | at = numpy.transpose(a) |
---|
178 | inv_cov = numpy.zeros([nfunc,nfunc]) |
---|
179 | for i in range(nfunc): |
---|
180 | for j in range(nfunc): |
---|
181 | inv_cov[i][j] = 0.0 |
---|
182 | for k in range(npts): |
---|
183 | inv_cov[i][j] = at[i][k]*a[k][j] |
---|
184 | |
---|
185 | # Compute the reg term size for the output |
---|
186 | sum_sig = 0.0 |
---|
187 | sum_reg = 0.0 |
---|
188 | for j in range(nfunc): |
---|
189 | for i in range(npts): |
---|
190 | sum_sig += (a[i][j])**2 |
---|
191 | for i in range(nq): |
---|
192 | sum_reg += (a[i_q+npts][j])**2 |
---|
193 | |
---|
194 | if math.fabs(self.alpha)>0: |
---|
195 | new_alpha = sum_sig/(sum_reg/self.alpha) |
---|
196 | else: |
---|
197 | new_alpha = 0.0 |
---|
198 | self.suggested_alpha = new_alpha |
---|
199 | |
---|
200 | try: |
---|
201 | err = math.fabs(chi2/(npts-nfunc))* inv_cov |
---|
202 | except: |
---|
203 | print "Error estimating uncertainties" |
---|
204 | |
---|
205 | |
---|
206 | return c, err |
---|
207 | |
---|
208 | def svd(self, nfunc=5): |
---|
209 | import math, time |
---|
210 | # Ac - b = 0 |
---|
211 | |
---|
212 | A = numpy.zeros([nfunc, nfunc]) |
---|
213 | y = numpy.zeros(nfunc) |
---|
214 | |
---|
215 | t_0 = time.time() |
---|
216 | for i in range(nfunc): |
---|
217 | # A |
---|
218 | for j in range(nfunc): |
---|
219 | A[i][j] = 0.0 |
---|
220 | for k in range(len(self.x)): |
---|
221 | err = self.err[k] |
---|
222 | A[i][j] += 1.0/err/err*self.basefunc_ft(self.d_max, j+1, self.x[k]) \ |
---|
223 | *self.basefunc_ft(self.d_max, i+1, self.x[k]); |
---|
224 | #print A[i][j] |
---|
225 | #A[i][j] -= self.alpha*(math.cos(math.pi*(i+j)) - math.cos(math.pi*(i-j))); |
---|
226 | if i==j: |
---|
227 | A[i][j] += -1.0*self.alpha |
---|
228 | elif i-j==1 or i-j==-1: |
---|
229 | A[i][j] += 1.0*self.alpha |
---|
230 | #print " ",A[i][j] |
---|
231 | # y |
---|
232 | y[i] = 0.0 |
---|
233 | for k in range(len(self.x)): |
---|
234 | y[i] = self.y[k]/self.err[k]/self.err[k]*self.basefunc_ft(self.d_max, i+1, self.x[k]) |
---|
235 | |
---|
236 | print time.time()-t_0, 'secs' |
---|
237 | |
---|
238 | # use numpy.pinv(A) |
---|
239 | #inv_A = numpy.linalg.inv(A) |
---|
240 | #c = y*inv_A |
---|
241 | print y |
---|
242 | c = numpy.linalg.solve(A, y) |
---|
243 | |
---|
244 | |
---|
245 | print c |
---|
246 | |
---|
247 | err = numpy.zeros(len(c)) |
---|
248 | return c, err |
---|
249 | |
---|
250 | |
---|
251 | |
---|
252 | |
---|
253 | |
---|
254 | if __name__ == "__main__": |
---|
255 | o = Invertor() |
---|
256 | |
---|
257 | |
---|
258 | |
---|
259 | |
---|
260 | |
---|