[3350ad6] | 1 | # pylint: disable=invalid-name |
---|
[51f14603] | 2 | """ |
---|
| 3 | Module to perform P(r) inversion. |
---|
| 4 | The module contains the Invertor class. |
---|
[bc3e38c] | 5 | |
---|
| 6 | FIXME: The way the Invertor interacts with its C component should be cleaned up |
---|
[51f14603] | 7 | """ |
---|
| 8 | |
---|
[9a5097c] | 9 | import numpy as np |
---|
[51f14603] | 10 | import sys |
---|
| 11 | import math |
---|
| 12 | import time |
---|
| 13 | import copy |
---|
| 14 | import os |
---|
| 15 | import re |
---|
[3350ad6] | 16 | import logging |
---|
[51f14603] | 17 | from numpy.linalg import lstsq |
---|
| 18 | from scipy import optimize |
---|
[b699768] | 19 | from sas.sascalc.pr.core.pr_inversion import Cinvertor |
---|
[51f14603] | 20 | |
---|
[463e7ffc] | 21 | logger = logging.getLogger(__name__) |
---|
[c155a16] | 22 | |
---|
[51f14603] | 23 | def help(): |
---|
| 24 | """ |
---|
| 25 | Provide general online help text |
---|
| 26 | Future work: extend this function to allow topic selection |
---|
| 27 | """ |
---|
[3350ad6] | 28 | info_txt = "The inversion approach is based on Moore, J. Appl. Cryst. " |
---|
[51f14603] | 29 | info_txt += "(1980) 13, 168-175.\n\n" |
---|
| 30 | info_txt += "P(r) is set to be equal to an expansion of base functions " |
---|
| 31 | info_txt += "of the type " |
---|
| 32 | info_txt += "phi_n(r) = 2*r*sin(pi*n*r/D_max). The coefficient of each " |
---|
| 33 | info_txt += "base functions " |
---|
| 34 | info_txt += "in the expansion is found by performing a least square fit " |
---|
| 35 | info_txt += "with the " |
---|
| 36 | info_txt += "following fit function:\n\n" |
---|
| 37 | info_txt += "chi**2 = sum_i[ I_meas(q_i) - I_th(q_i) ]**2/error**2 +" |
---|
| 38 | info_txt += "Reg_term\n\n" |
---|
| 39 | info_txt += "where I_meas(q) is the measured scattering intensity and " |
---|
| 40 | info_txt += "I_th(q) is " |
---|
| 41 | info_txt += "the prediction from the Fourier transform of the P(r) " |
---|
| 42 | info_txt += "expansion. " |
---|
| 43 | info_txt += "The Reg_term term is a regularization term set to the second" |
---|
| 44 | info_txt += " derivative " |
---|
| 45 | info_txt += "d**2P(r)/dr**2 integrated over r. It is used to produce " |
---|
| 46 | info_txt += "a smooth P(r) output.\n\n" |
---|
| 47 | info_txt += "The following are user inputs:\n\n" |
---|
| 48 | info_txt += " - Number of terms: the number of base functions in the P(r)" |
---|
| 49 | info_txt += " expansion.\n\n" |
---|
| 50 | info_txt += " - Regularization constant: a multiplicative constant " |
---|
| 51 | info_txt += "to set the size of " |
---|
| 52 | info_txt += "the regularization term.\n\n" |
---|
| 53 | info_txt += " - Maximum distance: the maximum distance between any " |
---|
| 54 | info_txt += "two points in the system.\n" |
---|
[3350ad6] | 55 | |
---|
[51f14603] | 56 | return info_txt |
---|
[3350ad6] | 57 | |
---|
[51f14603] | 58 | |
---|
| 59 | class Invertor(Cinvertor): |
---|
| 60 | """ |
---|
| 61 | Invertor class to perform P(r) inversion |
---|
[3350ad6] | 62 | |
---|
[51f14603] | 63 | The problem is solved by posing the problem as Ax = b, |
---|
| 64 | where x is the set of coefficients we are looking for. |
---|
[3350ad6] | 65 | |
---|
[51f14603] | 66 | Npts is the number of points. |
---|
[3350ad6] | 67 | |
---|
[51f14603] | 68 | In the following i refers to the ith base function coefficient. |
---|
| 69 | The matrix has its entries j in its first Npts rows set to :: |
---|
| 70 | |
---|
| 71 | A[j][i] = (Fourier transformed base function for point j) |
---|
[3350ad6] | 72 | |
---|
[6701a0b] | 73 | We then choose a number of r-points, n_r, to evaluate the second |
---|
[51f14603] | 74 | derivative of P(r) at. This is used as our regularization term. |
---|
| 75 | For a vector r of length n_r, the following n_r rows are set to :: |
---|
| 76 | |
---|
| 77 | A[j+Npts][i] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, |
---|
| 78 | evaluated at r[j]) |
---|
[3350ad6] | 79 | |
---|
[51f14603] | 80 | The vector b has its first Npts entries set to :: |
---|
| 81 | |
---|
| 82 | b[j] = (I(q) observed for point j) |
---|
[3350ad6] | 83 | |
---|
[51f14603] | 84 | The following n_r entries are set to zero. |
---|
[3350ad6] | 85 | |
---|
[51f14603] | 86 | The result is found by using scipy.linalg.basic.lstsq to invert |
---|
| 87 | the matrix and find the coefficients x. |
---|
[3350ad6] | 88 | |
---|
[51f14603] | 89 | Methods inherited from Cinvertor: |
---|
| 90 | |
---|
| 91 | * ``get_peaks(pars)``: returns the number of P(r) peaks |
---|
| 92 | * ``oscillations(pars)``: returns the oscillation parameters for the output P(r) |
---|
| 93 | * ``get_positive(pars)``: returns the fraction of P(r) that is above zero |
---|
| 94 | * ``get_pos_err(pars)``: returns the fraction of P(r) that is 1-sigma above zero |
---|
| 95 | """ |
---|
| 96 | ## Chisqr of the last computation |
---|
[3350ad6] | 97 | chi2 = 0 |
---|
[51f14603] | 98 | ## Time elapsed for last computation |
---|
| 99 | elapsed = 0 |
---|
| 100 | ## Alpha to get the reg term the same size as the signal |
---|
| 101 | suggested_alpha = 0 |
---|
| 102 | ## Last number of base functions used |
---|
| 103 | nfunc = 10 |
---|
| 104 | ## Last output values |
---|
| 105 | out = None |
---|
| 106 | ## Last errors on output values |
---|
| 107 | cov = None |
---|
| 108 | ## Background value |
---|
| 109 | background = 0 |
---|
| 110 | ## Information dictionary for application use |
---|
| 111 | info = {} |
---|
[3350ad6] | 112 | |
---|
[51f14603] | 113 | def __init__(self): |
---|
| 114 | Cinvertor.__init__(self) |
---|
[3350ad6] | 115 | |
---|
[51f14603] | 116 | def __setstate__(self, state): |
---|
| 117 | """ |
---|
| 118 | restore the state of invertor for pickle |
---|
| 119 | """ |
---|
| 120 | (self.__dict__, self.alpha, self.d_max, |
---|
| 121 | self.q_min, self.q_max, |
---|
| 122 | self.x, self.y, |
---|
[cb62bd5] | 123 | self.err, self.est_bck, |
---|
[51f14603] | 124 | self.slit_height, self.slit_width) = state |
---|
[3350ad6] | 125 | |
---|
[51f14603] | 126 | def __reduce_ex__(self, proto): |
---|
| 127 | """ |
---|
| 128 | Overwrite the __reduce_ex__ |
---|
| 129 | """ |
---|
| 130 | |
---|
| 131 | state = (self.__dict__, |
---|
| 132 | self.alpha, self.d_max, |
---|
| 133 | self.q_min, self.q_max, |
---|
| 134 | self.x, self.y, |
---|
[cb62bd5] | 135 | self.err, self.est_bck, |
---|
[51f14603] | 136 | self.slit_height, self.slit_width, |
---|
[3350ad6] | 137 | ) |
---|
[51f14603] | 138 | return (Invertor, tuple(), state, None, None) |
---|
[3350ad6] | 139 | |
---|
[51f14603] | 140 | def __setattr__(self, name, value): |
---|
| 141 | """ |
---|
| 142 | Set the value of an attribute. |
---|
| 143 | Access the parent class methods for |
---|
| 144 | x, y, err, d_max, q_min, q_max and alpha |
---|
| 145 | """ |
---|
[fba204bf] | 146 | if name == 'x': |
---|
[51f14603] | 147 | if 0.0 in value: |
---|
| 148 | msg = "Invertor: one of your q-values is zero. " |
---|
| 149 | msg += "Delete that entry before proceeding" |
---|
[574adc7] | 150 | raise ValueError(msg) |
---|
[51f14603] | 151 | return self.set_x(value) |
---|
| 152 | elif name == 'y': |
---|
| 153 | return self.set_y(value) |
---|
| 154 | elif name == 'err': |
---|
| 155 | value2 = abs(value) |
---|
| 156 | return self.set_err(value2) |
---|
| 157 | elif name == 'd_max': |
---|
[2c60f304] | 158 | if value <= 0.0: |
---|
| 159 | msg = "Invertor: d_max must be greater than zero." |
---|
| 160 | msg += "Correct that entry before proceeding" |
---|
[574adc7] | 161 | raise ValueError(msg) |
---|
[51f14603] | 162 | return self.set_dmax(value) |
---|
| 163 | elif name == 'q_min': |
---|
[235f514] | 164 | if value is None: |
---|
[51f14603] | 165 | return self.set_qmin(-1.0) |
---|
| 166 | return self.set_qmin(value) |
---|
| 167 | elif name == 'q_max': |
---|
[235f514] | 168 | if value is None: |
---|
[51f14603] | 169 | return self.set_qmax(-1.0) |
---|
| 170 | return self.set_qmax(value) |
---|
| 171 | elif name == 'alpha': |
---|
| 172 | return self.set_alpha(value) |
---|
| 173 | elif name == 'slit_height': |
---|
| 174 | return self.set_slit_height(value) |
---|
| 175 | elif name == 'slit_width': |
---|
| 176 | return self.set_slit_width(value) |
---|
[cb62bd5] | 177 | elif name == 'est_bck': |
---|
[51f14603] | 178 | if value == True: |
---|
[cb62bd5] | 179 | return self.set_est_bck(1) |
---|
[51f14603] | 180 | elif value == False: |
---|
[cb62bd5] | 181 | return self.set_est_bck(0) |
---|
[51f14603] | 182 | else: |
---|
[574adc7] | 183 | raise ValueError("Invertor: est_bck can only be True or False") |
---|
[3350ad6] | 184 | |
---|
[51f14603] | 185 | return Cinvertor.__setattr__(self, name, value) |
---|
[3350ad6] | 186 | |
---|
[51f14603] | 187 | def __getattr__(self, name): |
---|
| 188 | """ |
---|
| 189 | Return the value of an attribute |
---|
| 190 | """ |
---|
| 191 | #import numpy |
---|
| 192 | if name == 'x': |
---|
[9a5097c] | 193 | out = np.ones(self.get_nx()) |
---|
[51f14603] | 194 | self.get_x(out) |
---|
| 195 | return out |
---|
| 196 | elif name == 'y': |
---|
[9a5097c] | 197 | out = np.ones(self.get_ny()) |
---|
[51f14603] | 198 | self.get_y(out) |
---|
| 199 | return out |
---|
| 200 | elif name == 'err': |
---|
[9a5097c] | 201 | out = np.ones(self.get_nerr()) |
---|
[51f14603] | 202 | self.get_err(out) |
---|
| 203 | return out |
---|
| 204 | elif name == 'd_max': |
---|
| 205 | return self.get_dmax() |
---|
| 206 | elif name == 'q_min': |
---|
| 207 | qmin = self.get_qmin() |
---|
| 208 | if qmin < 0: |
---|
| 209 | return None |
---|
| 210 | return qmin |
---|
| 211 | elif name == 'q_max': |
---|
| 212 | qmax = self.get_qmax() |
---|
| 213 | if qmax < 0: |
---|
| 214 | return None |
---|
| 215 | return qmax |
---|
| 216 | elif name == 'alpha': |
---|
| 217 | return self.get_alpha() |
---|
| 218 | elif name == 'slit_height': |
---|
| 219 | return self.get_slit_height() |
---|
| 220 | elif name == 'slit_width': |
---|
| 221 | return self.get_slit_width() |
---|
[cb62bd5] | 222 | elif name == 'est_bck': |
---|
| 223 | value = self.get_est_bck() |
---|
[b8080e1] | 224 | return value == 1 |
---|
[51f14603] | 225 | elif name in self.__dict__: |
---|
| 226 | return self.__dict__[name] |
---|
| 227 | return None |
---|
[3350ad6] | 228 | |
---|
[aed159f] | 229 | def add_errors(self, sigma=0.05): |
---|
[fba204bf] | 230 | """ |
---|
[aed159f] | 231 | Adds errors to data set is they are not available. |
---|
| 232 | Uses $\Delta y = \sigma | y |$. |
---|
[fba204bf] | 233 | """ |
---|
[aed159f] | 234 | self.dy = sigma * np.fabs(self.y) |
---|
[fba204bf] | 235 | |
---|
[51f14603] | 236 | def clone(self): |
---|
| 237 | """ |
---|
| 238 | Return a clone of this instance |
---|
| 239 | """ |
---|
| 240 | #import copy |
---|
[3350ad6] | 241 | |
---|
[51f14603] | 242 | invertor = Invertor() |
---|
[3350ad6] | 243 | invertor.chi2 = self.chi2 |
---|
[51f14603] | 244 | invertor.elapsed = self.elapsed |
---|
[3350ad6] | 245 | invertor.nfunc = self.nfunc |
---|
| 246 | invertor.alpha = self.alpha |
---|
| 247 | invertor.d_max = self.d_max |
---|
| 248 | invertor.q_min = self.q_min |
---|
| 249 | invertor.q_max = self.q_max |
---|
| 250 | |
---|
[51f14603] | 251 | invertor.x = self.x |
---|
| 252 | invertor.y = self.y |
---|
[eeea6a3] | 253 | invertor.err = self.err |
---|
[cb62bd5] | 254 | invertor.est_bck = self.est_bck |
---|
| 255 | invertor.background = self.background |
---|
[51f14603] | 256 | invertor.slit_height = self.slit_height |
---|
| 257 | invertor.slit_width = self.slit_width |
---|
[3350ad6] | 258 | |
---|
[51f14603] | 259 | invertor.info = copy.deepcopy(self.info) |
---|
[3350ad6] | 260 | |
---|
[51f14603] | 261 | return invertor |
---|
[3350ad6] | 262 | |
---|
[51f14603] | 263 | def invert(self, nfunc=10, nr=20): |
---|
| 264 | """ |
---|
| 265 | Perform inversion to P(r) |
---|
[3350ad6] | 266 | |
---|
[51f14603] | 267 | The problem is solved by posing the problem as Ax = b, |
---|
| 268 | where x is the set of coefficients we are looking for. |
---|
[3350ad6] | 269 | |
---|
[51f14603] | 270 | Npts is the number of points. |
---|
[3350ad6] | 271 | |
---|
[51f14603] | 272 | In the following i refers to the ith base function coefficient. |
---|
| 273 | The matrix has its entries j in its first Npts rows set to :: |
---|
| 274 | |
---|
| 275 | A[i][j] = (Fourier transformed base function for point j) |
---|
[3350ad6] | 276 | |
---|
[6701a0b] | 277 | We then choose a number of r-points, n_r, to evaluate the second |
---|
[51f14603] | 278 | derivative of P(r) at. This is used as our regularization term. |
---|
| 279 | For a vector r of length n_r, the following n_r rows are set to :: |
---|
| 280 | |
---|
| 281 | A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, evaluated at r[j]) |
---|
[3350ad6] | 282 | |
---|
[51f14603] | 283 | The vector b has its first Npts entries set to :: |
---|
| 284 | |
---|
| 285 | b[j] = (I(q) observed for point j) |
---|
[3350ad6] | 286 | |
---|
[51f14603] | 287 | The following n_r entries are set to zero. |
---|
[3350ad6] | 288 | |
---|
[51f14603] | 289 | The result is found by using scipy.linalg.basic.lstsq to invert |
---|
| 290 | the matrix and find the coefficients x. |
---|
[3350ad6] | 291 | |
---|
[51f14603] | 292 | :param nfunc: number of base functions to use. |
---|
| 293 | :param nr: number of r points to evaluate the 2nd derivative at for the reg. term. |
---|
| 294 | :return: c_out, c_cov - the coefficients with covariance matrix |
---|
| 295 | """ |
---|
| 296 | # Reset the background value before proceeding |
---|
[cb62bd5] | 297 | # self.background = 0.0 |
---|
| 298 | if not self.est_bck: |
---|
| 299 | self.y -= self.background |
---|
| 300 | out, cov = self.lstsq(nfunc, nr=nr) |
---|
| 301 | if not self.est_bck: |
---|
| 302 | self.y += self.background |
---|
| 303 | return out, cov |
---|
[3350ad6] | 304 | |
---|
[51f14603] | 305 | def iq(self, out, q): |
---|
| 306 | """ |
---|
| 307 | Function to call to evaluate the scattering intensity |
---|
[3350ad6] | 308 | |
---|
[51f14603] | 309 | :param args: c-parameters, and q |
---|
| 310 | :return: I(q) |
---|
[3350ad6] | 311 | |
---|
[51f14603] | 312 | """ |
---|
| 313 | return Cinvertor.iq(self, out, q) + self.background |
---|
[3350ad6] | 314 | |
---|
[51f14603] | 315 | def invert_optimize(self, nfunc=10, nr=20): |
---|
| 316 | """ |
---|
| 317 | Slower version of the P(r) inversion that uses scipy.optimize.leastsq. |
---|
[3350ad6] | 318 | |
---|
[51f14603] | 319 | This probably produce more reliable results, but is much slower. |
---|
| 320 | The minimization function is set to |
---|
| 321 | sum_i[ (I_obs(q_i) - I_theo(q_i))/err**2 ] + alpha * reg_term, |
---|
| 322 | where the reg_term is given by Svergun: it is the integral of |
---|
| 323 | the square of the first derivative |
---|
| 324 | of P(r), d(P(r))/dr, integrated over the full range of r. |
---|
[3350ad6] | 325 | |
---|
[51f14603] | 326 | :param nfunc: number of base functions to use. |
---|
| 327 | :param nr: number of r points to evaluate the 2nd derivative at |
---|
| 328 | for the reg. term. |
---|
[3350ad6] | 329 | |
---|
[51f14603] | 330 | :return: c_out, c_cov - the coefficients with covariance matrix |
---|
[3350ad6] | 331 | |
---|
[51f14603] | 332 | """ |
---|
| 333 | self.nfunc = nfunc |
---|
| 334 | # First, check that the current data is valid |
---|
| 335 | if self.is_valid() <= 0: |
---|
| 336 | msg = "Invertor.invert: Data array are of different length" |
---|
[574adc7] | 337 | raise RuntimeError(msg) |
---|
[3350ad6] | 338 | |
---|
[9a5097c] | 339 | p = np.ones(nfunc) |
---|
[51f14603] | 340 | t_0 = time.time() |
---|
[3350ad6] | 341 | out, cov_x, _, _, _ = optimize.leastsq(self.residuals, p, full_output=1) |
---|
| 342 | |
---|
[51f14603] | 343 | # Compute chi^2 |
---|
| 344 | res = self.residuals(out) |
---|
| 345 | chisqr = 0 |
---|
| 346 | for i in range(len(res)): |
---|
| 347 | chisqr += res[i] |
---|
[3350ad6] | 348 | |
---|
[51f14603] | 349 | self.chi2 = chisqr |
---|
| 350 | |
---|
| 351 | # Store computation time |
---|
| 352 | self.elapsed = time.time() - t_0 |
---|
[3350ad6] | 353 | |
---|
[51f14603] | 354 | if cov_x is None: |
---|
[9a5097c] | 355 | cov_x = np.ones([nfunc, nfunc]) |
---|
[51f14603] | 356 | cov_x *= math.fabs(chisqr) |
---|
| 357 | return out, cov_x |
---|
[3350ad6] | 358 | |
---|
[51f14603] | 359 | def pr_fit(self, nfunc=5): |
---|
| 360 | """ |
---|
| 361 | This is a direct fit to a given P(r). It assumes that the y data |
---|
| 362 | is set to some P(r) distribution that we are trying to reproduce |
---|
| 363 | with a set of base functions. |
---|
[3350ad6] | 364 | |
---|
[51f14603] | 365 | This method is provided as a test. |
---|
| 366 | """ |
---|
| 367 | # First, check that the current data is valid |
---|
| 368 | if self.is_valid() <= 0: |
---|
| 369 | msg = "Invertor.invert: Data arrays are of different length" |
---|
[574adc7] | 370 | raise RuntimeError(msg) |
---|
[3350ad6] | 371 | |
---|
[9a5097c] | 372 | p = np.ones(nfunc) |
---|
[51f14603] | 373 | t_0 = time.time() |
---|
[3350ad6] | 374 | out, cov_x, _, _, _ = optimize.leastsq(self.pr_residuals, p, full_output=1) |
---|
| 375 | |
---|
[51f14603] | 376 | # Compute chi^2 |
---|
| 377 | res = self.pr_residuals(out) |
---|
| 378 | chisqr = 0 |
---|
| 379 | for i in range(len(res)): |
---|
| 380 | chisqr += res[i] |
---|
[3350ad6] | 381 | |
---|
[51f14603] | 382 | self.chisqr = chisqr |
---|
[3350ad6] | 383 | |
---|
[51f14603] | 384 | # Store computation time |
---|
| 385 | self.elapsed = time.time() - t_0 |
---|
| 386 | |
---|
| 387 | return out, cov_x |
---|
[3350ad6] | 388 | |
---|
[51f14603] | 389 | def pr_err(self, c, c_cov, r): |
---|
| 390 | """ |
---|
| 391 | Returns the value of P(r) for a given r, and base function |
---|
| 392 | coefficients, with error. |
---|
[3350ad6] | 393 | |
---|
[51f14603] | 394 | :param c: base function coefficients |
---|
| 395 | :param c_cov: covariance matrice of the base function coefficients |
---|
| 396 | :param r: r-value to evaluate P(r) at |
---|
[3350ad6] | 397 | |
---|
[51f14603] | 398 | :return: P(r) |
---|
[3350ad6] | 399 | |
---|
[51f14603] | 400 | """ |
---|
| 401 | return self.get_pr_err(c, c_cov, r) |
---|
[3350ad6] | 402 | |
---|
[51f14603] | 403 | def _accept_q(self, q): |
---|
| 404 | """ |
---|
| 405 | Check q-value against user-defined range |
---|
| 406 | """ |
---|
[ac07a3a] | 407 | if self.q_min is not None and q < self.q_min: |
---|
[51f14603] | 408 | return False |
---|
[ac07a3a] | 409 | if self.q_max is not None and q > self.q_max: |
---|
[51f14603] | 410 | return False |
---|
| 411 | return True |
---|
[3350ad6] | 412 | |
---|
[51f14603] | 413 | def lstsq(self, nfunc=5, nr=20): |
---|
| 414 | """ |
---|
| 415 | The problem is solved by posing the problem as Ax = b, |
---|
| 416 | where x is the set of coefficients we are looking for. |
---|
[3350ad6] | 417 | |
---|
[51f14603] | 418 | Npts is the number of points. |
---|
[3350ad6] | 419 | |
---|
[51f14603] | 420 | In the following i refers to the ith base function coefficient. |
---|
| 421 | The matrix has its entries j in its first Npts rows set to :: |
---|
| 422 | |
---|
| 423 | A[i][j] = (Fourier transformed base function for point j) |
---|
[3350ad6] | 424 | |
---|
[51f14603] | 425 | We them choose a number of r-points, n_r, to evaluate the second |
---|
| 426 | derivative of P(r) at. This is used as our regularization term. |
---|
| 427 | For a vector r of length n_r, the following n_r rows are set to :: |
---|
| 428 | |
---|
| 429 | A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, |
---|
| 430 | evaluated at r[j]) |
---|
[3350ad6] | 431 | |
---|
[51f14603] | 432 | The vector b has its first Npts entries set to :: |
---|
| 433 | |
---|
| 434 | b[j] = (I(q) observed for point j) |
---|
[3350ad6] | 435 | |
---|
[51f14603] | 436 | The following n_r entries are set to zero. |
---|
[3350ad6] | 437 | |
---|
[51f14603] | 438 | The result is found by using scipy.linalg.basic.lstsq to invert |
---|
| 439 | the matrix and find the coefficients x. |
---|
[3350ad6] | 440 | |
---|
[51f14603] | 441 | :param nfunc: number of base functions to use. |
---|
| 442 | :param nr: number of r points to evaluate the 2nd derivative at for the reg. term. |
---|
| 443 | |
---|
| 444 | If the result does not allow us to compute the covariance matrix, |
---|
| 445 | a matrix filled with zeros will be returned. |
---|
| 446 | |
---|
| 447 | """ |
---|
| 448 | # Note: To make sure an array is contiguous: |
---|
[9a5097c] | 449 | # blah = np.ascontiguousarray(blah_original) |
---|
[51f14603] | 450 | # ... before passing it to C |
---|
[3350ad6] | 451 | |
---|
[51f14603] | 452 | if self.is_valid() < 0: |
---|
| 453 | msg = "Invertor: invalid data; incompatible data lengths." |
---|
[574adc7] | 454 | raise RuntimeError(msg) |
---|
[3350ad6] | 455 | |
---|
[51f14603] | 456 | self.nfunc = nfunc |
---|
| 457 | # a -- An M x N matrix. |
---|
| 458 | # b -- An M x nrhs matrix or M vector. |
---|
| 459 | npts = len(self.x) |
---|
[3350ad6] | 460 | nq = nr |
---|
[51f14603] | 461 | sqrt_alpha = math.sqrt(math.fabs(self.alpha)) |
---|
| 462 | if sqrt_alpha < 0.0: |
---|
| 463 | nq = 0 |
---|
| 464 | |
---|
| 465 | # If we need to fit the background, add a term |
---|
[b8080e1] | 466 | if self.est_bck: |
---|
[51f14603] | 467 | nfunc_0 = nfunc |
---|
| 468 | nfunc += 1 |
---|
| 469 | |
---|
[9a5097c] | 470 | a = np.zeros([npts + nq, nfunc]) |
---|
| 471 | b = np.zeros(npts + nq) |
---|
| 472 | err = np.zeros([nfunc, nfunc]) |
---|
[3350ad6] | 473 | |
---|
[51f14603] | 474 | # Construct the a matrix and b vector that represent the problem |
---|
| 475 | t_0 = time.time() |
---|
| 476 | try: |
---|
| 477 | self._get_matrix(nfunc, nq, a, b) |
---|
[d04ac05] | 478 | except Exception as exc: |
---|
| 479 | raise RuntimeError("Invertor: could not invert I(Q)\n %s" % str(exc)) |
---|
[3350ad6] | 480 | |
---|
[51f14603] | 481 | # Perform the inversion (least square fit) |
---|
| 482 | c, chi2, _, _ = lstsq(a, b) |
---|
| 483 | # Sanity check |
---|
| 484 | try: |
---|
| 485 | float(chi2) |
---|
| 486 | except: |
---|
[6da860a] | 487 | chi2 = -1.0 |
---|
[51f14603] | 488 | self.chi2 = chi2 |
---|
[3350ad6] | 489 | |
---|
[9a5097c] | 490 | inv_cov = np.zeros([nfunc, nfunc]) |
---|
[51f14603] | 491 | # Get the covariance matrix, defined as inv_cov = a_transposed * a |
---|
| 492 | self._get_invcov_matrix(nfunc, nr, a, inv_cov) |
---|
[3350ad6] | 493 | |
---|
[51f14603] | 494 | # Compute the reg term size for the output |
---|
| 495 | sum_sig, sum_reg = self._get_reg_size(nfunc, nr, a) |
---|
[3350ad6] | 496 | |
---|
[51f14603] | 497 | if math.fabs(self.alpha) > 0: |
---|
| 498 | new_alpha = sum_sig / (sum_reg / self.alpha) |
---|
| 499 | else: |
---|
| 500 | new_alpha = 0.0 |
---|
| 501 | self.suggested_alpha = new_alpha |
---|
[3350ad6] | 502 | |
---|
[51f14603] | 503 | try: |
---|
[9a5097c] | 504 | cov = np.linalg.pinv(inv_cov) |
---|
[51f14603] | 505 | err = math.fabs(chi2 / float(npts - nfunc)) * cov |
---|
[b8080e1] | 506 | except: |
---|
[51f14603] | 507 | # We were not able to estimate the errors |
---|
| 508 | # Return an empty error matrix |
---|
[aed159f] | 509 | logger.error(sys.exc_info()[1]) |
---|
[3350ad6] | 510 | |
---|
[51f14603] | 511 | # Keep a copy of the last output |
---|
[b8080e1] | 512 | if not self.est_bck: |
---|
[51f14603] | 513 | self.out = c |
---|
| 514 | self.cov = err |
---|
| 515 | else: |
---|
| 516 | self.background = c[0] |
---|
[3350ad6] | 517 | |
---|
[9a5097c] | 518 | err_0 = np.zeros([nfunc, nfunc]) |
---|
| 519 | c_0 = np.zeros(nfunc) |
---|
[3350ad6] | 520 | |
---|
[51f14603] | 521 | for i in range(nfunc_0): |
---|
[3350ad6] | 522 | c_0[i] = c[i + 1] |
---|
[51f14603] | 523 | for j in range(nfunc_0): |
---|
[3350ad6] | 524 | err_0[i][j] = err[i + 1][j + 1] |
---|
| 525 | |
---|
[51f14603] | 526 | self.out = c_0 |
---|
| 527 | self.cov = err_0 |
---|
[3350ad6] | 528 | |
---|
[51f14603] | 529 | # Store computation time |
---|
| 530 | self.elapsed = time.time() - t_0 |
---|
[3350ad6] | 531 | |
---|
[51f14603] | 532 | return self.out, self.cov |
---|
[3350ad6] | 533 | |
---|
[51f14603] | 534 | def estimate_numterms(self, isquit_func=None): |
---|
| 535 | """ |
---|
| 536 | Returns a reasonable guess for the |
---|
| 537 | number of terms |
---|
[3350ad6] | 538 | |
---|
[51f14603] | 539 | :param isquit_func: |
---|
| 540 | reference to thread function to call to check whether the computation needs to |
---|
| 541 | be stopped. |
---|
[3350ad6] | 542 | |
---|
[51f14603] | 543 | :return: number of terms, alpha, message |
---|
[3350ad6] | 544 | |
---|
[51f14603] | 545 | """ |
---|
[b8080e1] | 546 | from .num_term import NTermEstimator |
---|
[5f8fc78] | 547 | estimator = NTermEstimator(self.clone()) |
---|
[51f14603] | 548 | try: |
---|
| 549 | return estimator.num_terms(isquit_func) |
---|
[b8080e1] | 550 | except: |
---|
[51f14603] | 551 | # If we fail, estimate alpha and return the default |
---|
| 552 | # number of terms |
---|
| 553 | best_alpha, _, _ = self.estimate_alpha(self.nfunc) |
---|
[aed159f] | 554 | logger.warning("Invertor.estimate_numterms: %s" % sys.exc_info()[1]) |
---|
[51f14603] | 555 | return self.nfunc, best_alpha, "Could not estimate number of terms" |
---|
[3350ad6] | 556 | |
---|
[51f14603] | 557 | def estimate_alpha(self, nfunc): |
---|
| 558 | """ |
---|
| 559 | Returns a reasonable guess for the |
---|
| 560 | regularization constant alpha |
---|
[3350ad6] | 561 | |
---|
[51f14603] | 562 | :param nfunc: number of terms to use in the expansion. |
---|
[3350ad6] | 563 | |
---|
[51f14603] | 564 | :return: alpha, message, elapsed |
---|
[3350ad6] | 565 | |
---|
[51f14603] | 566 | where alpha is the estimate for alpha, |
---|
| 567 | message is a message for the user, |
---|
| 568 | elapsed is the computation time |
---|
| 569 | """ |
---|
| 570 | #import time |
---|
| 571 | try: |
---|
| 572 | pr = self.clone() |
---|
[3350ad6] | 573 | |
---|
[51f14603] | 574 | # T_0 for computation time |
---|
| 575 | starttime = time.time() |
---|
| 576 | elapsed = 0 |
---|
[3350ad6] | 577 | |
---|
[51f14603] | 578 | # If the current alpha is zero, try |
---|
| 579 | # another value |
---|
| 580 | if pr.alpha <= 0: |
---|
| 581 | pr.alpha = 0.0001 |
---|
[3350ad6] | 582 | |
---|
[51f14603] | 583 | # Perform inversion to find the largest alpha |
---|
| 584 | out, _ = pr.invert(nfunc) |
---|
| 585 | elapsed = time.time() - starttime |
---|
| 586 | initial_alpha = pr.alpha |
---|
| 587 | initial_peaks = pr.get_peaks(out) |
---|
[3350ad6] | 588 | |
---|
[51f14603] | 589 | # Try the inversion with the estimated alpha |
---|
| 590 | pr.alpha = pr.suggested_alpha |
---|
| 591 | out, _ = pr.invert(nfunc) |
---|
[3350ad6] | 592 | |
---|
[51f14603] | 593 | npeaks = pr.get_peaks(out) |
---|
| 594 | # if more than one peak to start with |
---|
| 595 | # just return the estimate |
---|
| 596 | if npeaks > 1: |
---|
| 597 | #message = "Your P(r) is not smooth, |
---|
| 598 | #please check your inversion parameters" |
---|
| 599 | message = None |
---|
| 600 | return pr.suggested_alpha, message, elapsed |
---|
| 601 | else: |
---|
[3350ad6] | 602 | |
---|
[51f14603] | 603 | # Look at smaller values |
---|
| 604 | # We assume that for the suggested alpha, we have 1 peak |
---|
| 605 | # if not, send a message to change parameters |
---|
| 606 | alpha = pr.suggested_alpha |
---|
| 607 | best_alpha = pr.suggested_alpha |
---|
| 608 | found = False |
---|
| 609 | for i in range(10): |
---|
[3350ad6] | 610 | pr.alpha = (0.33) ** (i + 1) * alpha |
---|
[51f14603] | 611 | out, _ = pr.invert(nfunc) |
---|
[3350ad6] | 612 | |
---|
[51f14603] | 613 | peaks = pr.get_peaks(out) |
---|
| 614 | if peaks > 1: |
---|
| 615 | found = True |
---|
| 616 | break |
---|
| 617 | best_alpha = pr.alpha |
---|
[3350ad6] | 618 | |
---|
[51f14603] | 619 | # If we didn't find a turning point for alpha and |
---|
| 620 | # the initial alpha already had only one peak, |
---|
| 621 | # just return that |
---|
| 622 | if not found and initial_peaks == 1 and \ |
---|
| 623 | initial_alpha < best_alpha: |
---|
| 624 | best_alpha = initial_alpha |
---|
[3350ad6] | 625 | |
---|
[51f14603] | 626 | # Check whether the size makes sense |
---|
| 627 | message = '' |
---|
[3350ad6] | 628 | |
---|
[51f14603] | 629 | if not found: |
---|
| 630 | message = None |
---|
| 631 | elif best_alpha >= 0.5 * pr.suggested_alpha: |
---|
| 632 | # best alpha is too big, return a |
---|
| 633 | # reasonable value |
---|
[3350ad6] | 634 | message = "The estimated alpha for your system is too " |
---|
[51f14603] | 635 | message += "large. " |
---|
| 636 | message += "Try increasing your maximum distance." |
---|
[3350ad6] | 637 | |
---|
[51f14603] | 638 | return best_alpha, message, elapsed |
---|
[3350ad6] | 639 | |
---|
[b8080e1] | 640 | except: |
---|
[aed159f] | 641 | message = "Invertor.estimate_alpha: %s" % sys.exc_info()[1] |
---|
[51f14603] | 642 | return 0, message, elapsed |
---|
[3350ad6] | 643 | |
---|
[51f14603] | 644 | def to_file(self, path, npts=100): |
---|
| 645 | """ |
---|
| 646 | Save the state to a file that will be readable |
---|
| 647 | by SliceView. |
---|
[3350ad6] | 648 | |
---|
[51f14603] | 649 | :param path: path of the file to write |
---|
| 650 | :param npts: number of P(r) points to be written |
---|
[3350ad6] | 651 | |
---|
[51f14603] | 652 | """ |
---|
| 653 | file = open(path, 'w') |
---|
| 654 | file.write("#d_max=%g\n" % self.d_max) |
---|
| 655 | file.write("#nfunc=%g\n" % self.nfunc) |
---|
| 656 | file.write("#alpha=%g\n" % self.alpha) |
---|
| 657 | file.write("#chi2=%g\n" % self.chi2) |
---|
| 658 | file.write("#elapsed=%g\n" % self.elapsed) |
---|
| 659 | file.write("#qmin=%s\n" % str(self.q_min)) |
---|
| 660 | file.write("#qmax=%s\n" % str(self.q_max)) |
---|
| 661 | file.write("#slit_height=%g\n" % self.slit_height) |
---|
| 662 | file.write("#slit_width=%g\n" % self.slit_width) |
---|
| 663 | file.write("#background=%g\n" % self.background) |
---|
[b8080e1] | 664 | if self.est_bck: |
---|
[51f14603] | 665 | file.write("#has_bck=1\n") |
---|
| 666 | else: |
---|
| 667 | file.write("#has_bck=0\n") |
---|
| 668 | file.write("#alpha_estimate=%g\n" % self.suggested_alpha) |
---|
[45dffa69] | 669 | if self.out is not None: |
---|
[51f14603] | 670 | if len(self.out) == len(self.cov): |
---|
| 671 | for i in range(len(self.out)): |
---|
| 672 | file.write("#C_%i=%s+-%s\n" % (i, str(self.out[i]), |
---|
[3350ad6] | 673 | str(self.cov[i][i]))) |
---|
[51f14603] | 674 | file.write("<r> <Pr> <dPr>\n") |
---|
[9a5097c] | 675 | r = np.arange(0.0, self.d_max, self.d_max / npts) |
---|
[3350ad6] | 676 | |
---|
[51f14603] | 677 | for r_i in r: |
---|
| 678 | (value, err) = self.pr_err(self.out, self.cov, r_i) |
---|
| 679 | file.write("%g %g %g\n" % (r_i, value, err)) |
---|
[3350ad6] | 680 | |
---|
[51f14603] | 681 | file.close() |
---|
[3350ad6] | 682 | |
---|
[51f14603] | 683 | def from_file(self, path): |
---|
| 684 | """ |
---|
| 685 | Load the state of the Invertor from a file, |
---|
| 686 | to be able to generate P(r) from a set of |
---|
| 687 | parameters. |
---|
[3350ad6] | 688 | |
---|
[51f14603] | 689 | :param path: path of the file to load |
---|
[3350ad6] | 690 | |
---|
[51f14603] | 691 | """ |
---|
| 692 | #import os |
---|
| 693 | #import re |
---|
| 694 | if os.path.isfile(path): |
---|
| 695 | try: |
---|
| 696 | fd = open(path, 'r') |
---|
[3350ad6] | 697 | |
---|
[51f14603] | 698 | buff = fd.read() |
---|
| 699 | lines = buff.split('\n') |
---|
| 700 | for line in lines: |
---|
| 701 | if line.startswith('#d_max='): |
---|
| 702 | toks = line.split('=') |
---|
| 703 | self.d_max = float(toks[1]) |
---|
| 704 | elif line.startswith('#nfunc='): |
---|
| 705 | toks = line.split('=') |
---|
| 706 | self.nfunc = int(toks[1]) |
---|
[9a5097c] | 707 | self.out = np.zeros(self.nfunc) |
---|
| 708 | self.cov = np.zeros([self.nfunc, self.nfunc]) |
---|
[51f14603] | 709 | elif line.startswith('#alpha='): |
---|
| 710 | toks = line.split('=') |
---|
| 711 | self.alpha = float(toks[1]) |
---|
| 712 | elif line.startswith('#chi2='): |
---|
| 713 | toks = line.split('=') |
---|
| 714 | self.chi2 = float(toks[1]) |
---|
| 715 | elif line.startswith('#elapsed='): |
---|
| 716 | toks = line.split('=') |
---|
| 717 | self.elapsed = float(toks[1]) |
---|
| 718 | elif line.startswith('#alpha_estimate='): |
---|
| 719 | toks = line.split('=') |
---|
| 720 | self.suggested_alpha = float(toks[1]) |
---|
| 721 | elif line.startswith('#qmin='): |
---|
| 722 | toks = line.split('=') |
---|
| 723 | try: |
---|
| 724 | self.q_min = float(toks[1]) |
---|
| 725 | except: |
---|
| 726 | self.q_min = None |
---|
| 727 | elif line.startswith('#qmax='): |
---|
| 728 | toks = line.split('=') |
---|
| 729 | try: |
---|
| 730 | self.q_max = float(toks[1]) |
---|
| 731 | except: |
---|
| 732 | self.q_max = None |
---|
| 733 | elif line.startswith('#slit_height='): |
---|
| 734 | toks = line.split('=') |
---|
| 735 | self.slit_height = float(toks[1]) |
---|
| 736 | elif line.startswith('#slit_width='): |
---|
| 737 | toks = line.split('=') |
---|
| 738 | self.slit_width = float(toks[1]) |
---|
| 739 | elif line.startswith('#background='): |
---|
| 740 | toks = line.split('=') |
---|
| 741 | self.background = float(toks[1]) |
---|
| 742 | elif line.startswith('#has_bck='): |
---|
| 743 | toks = line.split('=') |
---|
[b8080e1] | 744 | self.est_bck = int(toks[1]) == 1 |
---|
[3350ad6] | 745 | |
---|
[51f14603] | 746 | # Now read in the parameters |
---|
| 747 | elif line.startswith('#C_'): |
---|
| 748 | toks = line.split('=') |
---|
| 749 | p = re.compile('#C_([0-9]+)') |
---|
| 750 | m = p.search(toks[0]) |
---|
| 751 | toks2 = toks[1].split('+-') |
---|
| 752 | i = int(m.group(1)) |
---|
| 753 | self.out[i] = float(toks2[0]) |
---|
[3350ad6] | 754 | |
---|
[51f14603] | 755 | self.cov[i][i] = float(toks2[1]) |
---|
[3350ad6] | 756 | |
---|
[b8080e1] | 757 | except: |
---|
[aed159f] | 758 | msg = "Invertor.from_file: corrupted file\n%s" % sys.exc_info()[1] |
---|
[574adc7] | 759 | raise RuntimeError(msg) |
---|
[51f14603] | 760 | else: |
---|
| 761 | msg = "Invertor.from_file: '%s' is not a file" % str(path) |
---|
[574adc7] | 762 | raise RuntimeError(msg) |
---|