[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 | |
---|
[fba204bf] | 229 | def add_errors(self, yvalues): |
---|
| 230 | """ |
---|
| 231 | Adds errors to data set is they are not avaialble |
---|
| 232 | :return: |
---|
| 233 | """ |
---|
| 234 | stats_errors = np.zeros(len(yvalues)) |
---|
| 235 | for i in range(len(yvalues)): |
---|
| 236 | # Scale the error so that we can fit over several decades of Q |
---|
| 237 | scale = 0.05 * np.sqrt(yvalues[i]) |
---|
| 238 | min_err = 0.01 * yvalues[i] |
---|
| 239 | stats_errors[i] = scale * np.sqrt(np.fabs(yvalues[i])) + min_err |
---|
[6701a0b] | 240 | logger.warning("Simulated errors have been added to the data set\n") |
---|
[fba204bf] | 241 | return stats_errors |
---|
| 242 | |
---|
[51f14603] | 243 | def clone(self): |
---|
| 244 | """ |
---|
| 245 | Return a clone of this instance |
---|
| 246 | """ |
---|
| 247 | #import copy |
---|
[3350ad6] | 248 | |
---|
[51f14603] | 249 | invertor = Invertor() |
---|
[3350ad6] | 250 | invertor.chi2 = self.chi2 |
---|
[51f14603] | 251 | invertor.elapsed = self.elapsed |
---|
[3350ad6] | 252 | invertor.nfunc = self.nfunc |
---|
| 253 | invertor.alpha = self.alpha |
---|
| 254 | invertor.d_max = self.d_max |
---|
| 255 | invertor.q_min = self.q_min |
---|
| 256 | invertor.q_max = self.q_max |
---|
| 257 | |
---|
[51f14603] | 258 | invertor.x = self.x |
---|
| 259 | invertor.y = self.y |
---|
[eeea6a3] | 260 | invertor.err = self.err |
---|
[cb62bd5] | 261 | invertor.est_bck = self.est_bck |
---|
| 262 | invertor.background = self.background |
---|
[51f14603] | 263 | invertor.slit_height = self.slit_height |
---|
| 264 | invertor.slit_width = self.slit_width |
---|
[3350ad6] | 265 | |
---|
[51f14603] | 266 | invertor.info = copy.deepcopy(self.info) |
---|
[3350ad6] | 267 | |
---|
[51f14603] | 268 | return invertor |
---|
[3350ad6] | 269 | |
---|
[51f14603] | 270 | def invert(self, nfunc=10, nr=20): |
---|
| 271 | """ |
---|
| 272 | Perform inversion to P(r) |
---|
[3350ad6] | 273 | |
---|
[51f14603] | 274 | The problem is solved by posing the problem as Ax = b, |
---|
| 275 | where x is the set of coefficients we are looking for. |
---|
[3350ad6] | 276 | |
---|
[51f14603] | 277 | Npts is the number of points. |
---|
[3350ad6] | 278 | |
---|
[51f14603] | 279 | In the following i refers to the ith base function coefficient. |
---|
| 280 | The matrix has its entries j in its first Npts rows set to :: |
---|
| 281 | |
---|
| 282 | A[i][j] = (Fourier transformed base function for point j) |
---|
[3350ad6] | 283 | |
---|
[6701a0b] | 284 | We then choose a number of r-points, n_r, to evaluate the second |
---|
[51f14603] | 285 | derivative of P(r) at. This is used as our regularization term. |
---|
| 286 | For a vector r of length n_r, the following n_r rows are set to :: |
---|
| 287 | |
---|
| 288 | A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, evaluated at r[j]) |
---|
[3350ad6] | 289 | |
---|
[51f14603] | 290 | The vector b has its first Npts entries set to :: |
---|
| 291 | |
---|
| 292 | b[j] = (I(q) observed for point j) |
---|
[3350ad6] | 293 | |
---|
[51f14603] | 294 | The following n_r entries are set to zero. |
---|
[3350ad6] | 295 | |
---|
[51f14603] | 296 | The result is found by using scipy.linalg.basic.lstsq to invert |
---|
| 297 | the matrix and find the coefficients x. |
---|
[3350ad6] | 298 | |
---|
[51f14603] | 299 | :param nfunc: number of base functions to use. |
---|
| 300 | :param nr: number of r points to evaluate the 2nd derivative at for the reg. term. |
---|
| 301 | :return: c_out, c_cov - the coefficients with covariance matrix |
---|
| 302 | """ |
---|
| 303 | # Reset the background value before proceeding |
---|
[cb62bd5] | 304 | # self.background = 0.0 |
---|
| 305 | if not self.est_bck: |
---|
| 306 | self.y -= self.background |
---|
| 307 | out, cov = self.lstsq(nfunc, nr=nr) |
---|
| 308 | if not self.est_bck: |
---|
| 309 | self.y += self.background |
---|
| 310 | return out, cov |
---|
[3350ad6] | 311 | |
---|
[51f14603] | 312 | def iq(self, out, q): |
---|
| 313 | """ |
---|
| 314 | Function to call to evaluate the scattering intensity |
---|
[3350ad6] | 315 | |
---|
[51f14603] | 316 | :param args: c-parameters, and q |
---|
| 317 | :return: I(q) |
---|
[3350ad6] | 318 | |
---|
[51f14603] | 319 | """ |
---|
| 320 | return Cinvertor.iq(self, out, q) + self.background |
---|
[3350ad6] | 321 | |
---|
[51f14603] | 322 | def invert_optimize(self, nfunc=10, nr=20): |
---|
| 323 | """ |
---|
| 324 | Slower version of the P(r) inversion that uses scipy.optimize.leastsq. |
---|
[3350ad6] | 325 | |
---|
[51f14603] | 326 | This probably produce more reliable results, but is much slower. |
---|
| 327 | The minimization function is set to |
---|
| 328 | sum_i[ (I_obs(q_i) - I_theo(q_i))/err**2 ] + alpha * reg_term, |
---|
| 329 | where the reg_term is given by Svergun: it is the integral of |
---|
| 330 | the square of the first derivative |
---|
| 331 | of P(r), d(P(r))/dr, integrated over the full range of r. |
---|
[3350ad6] | 332 | |
---|
[51f14603] | 333 | :param nfunc: number of base functions to use. |
---|
| 334 | :param nr: number of r points to evaluate the 2nd derivative at |
---|
| 335 | for the reg. term. |
---|
[3350ad6] | 336 | |
---|
[51f14603] | 337 | :return: c_out, c_cov - the coefficients with covariance matrix |
---|
[3350ad6] | 338 | |
---|
[51f14603] | 339 | """ |
---|
| 340 | self.nfunc = nfunc |
---|
| 341 | # First, check that the current data is valid |
---|
| 342 | if self.is_valid() <= 0: |
---|
| 343 | msg = "Invertor.invert: Data array are of different length" |
---|
[574adc7] | 344 | raise RuntimeError(msg) |
---|
[3350ad6] | 345 | |
---|
[9a5097c] | 346 | p = np.ones(nfunc) |
---|
[51f14603] | 347 | t_0 = time.time() |
---|
[3350ad6] | 348 | out, cov_x, _, _, _ = optimize.leastsq(self.residuals, p, full_output=1) |
---|
| 349 | |
---|
[51f14603] | 350 | # Compute chi^2 |
---|
| 351 | res = self.residuals(out) |
---|
| 352 | chisqr = 0 |
---|
| 353 | for i in range(len(res)): |
---|
| 354 | chisqr += res[i] |
---|
[3350ad6] | 355 | |
---|
[51f14603] | 356 | self.chi2 = chisqr |
---|
| 357 | |
---|
| 358 | # Store computation time |
---|
| 359 | self.elapsed = time.time() - t_0 |
---|
[3350ad6] | 360 | |
---|
[51f14603] | 361 | if cov_x is None: |
---|
[9a5097c] | 362 | cov_x = np.ones([nfunc, nfunc]) |
---|
[51f14603] | 363 | cov_x *= math.fabs(chisqr) |
---|
| 364 | return out, cov_x |
---|
[3350ad6] | 365 | |
---|
[51f14603] | 366 | def pr_fit(self, nfunc=5): |
---|
| 367 | """ |
---|
| 368 | This is a direct fit to a given P(r). It assumes that the y data |
---|
| 369 | is set to some P(r) distribution that we are trying to reproduce |
---|
| 370 | with a set of base functions. |
---|
[3350ad6] | 371 | |
---|
[51f14603] | 372 | This method is provided as a test. |
---|
| 373 | """ |
---|
| 374 | # First, check that the current data is valid |
---|
| 375 | if self.is_valid() <= 0: |
---|
| 376 | msg = "Invertor.invert: Data arrays are of different length" |
---|
[574adc7] | 377 | raise RuntimeError(msg) |
---|
[3350ad6] | 378 | |
---|
[9a5097c] | 379 | p = np.ones(nfunc) |
---|
[51f14603] | 380 | t_0 = time.time() |
---|
[3350ad6] | 381 | out, cov_x, _, _, _ = optimize.leastsq(self.pr_residuals, p, full_output=1) |
---|
| 382 | |
---|
[51f14603] | 383 | # Compute chi^2 |
---|
| 384 | res = self.pr_residuals(out) |
---|
| 385 | chisqr = 0 |
---|
| 386 | for i in range(len(res)): |
---|
| 387 | chisqr += res[i] |
---|
[3350ad6] | 388 | |
---|
[51f14603] | 389 | self.chisqr = chisqr |
---|
[3350ad6] | 390 | |
---|
[51f14603] | 391 | # Store computation time |
---|
| 392 | self.elapsed = time.time() - t_0 |
---|
| 393 | |
---|
| 394 | return out, cov_x |
---|
[3350ad6] | 395 | |
---|
[51f14603] | 396 | def pr_err(self, c, c_cov, r): |
---|
| 397 | """ |
---|
| 398 | Returns the value of P(r) for a given r, and base function |
---|
| 399 | coefficients, with error. |
---|
[3350ad6] | 400 | |
---|
[51f14603] | 401 | :param c: base function coefficients |
---|
| 402 | :param c_cov: covariance matrice of the base function coefficients |
---|
| 403 | :param r: r-value to evaluate P(r) at |
---|
[3350ad6] | 404 | |
---|
[51f14603] | 405 | :return: P(r) |
---|
[3350ad6] | 406 | |
---|
[51f14603] | 407 | """ |
---|
| 408 | return self.get_pr_err(c, c_cov, r) |
---|
[3350ad6] | 409 | |
---|
[51f14603] | 410 | def _accept_q(self, q): |
---|
| 411 | """ |
---|
| 412 | Check q-value against user-defined range |
---|
| 413 | """ |
---|
[ac07a3a] | 414 | if self.q_min is not None and q < self.q_min: |
---|
[51f14603] | 415 | return False |
---|
[ac07a3a] | 416 | if self.q_max is not None and q > self.q_max: |
---|
[51f14603] | 417 | return False |
---|
| 418 | return True |
---|
[3350ad6] | 419 | |
---|
[51f14603] | 420 | def lstsq(self, nfunc=5, nr=20): |
---|
| 421 | """ |
---|
| 422 | The problem is solved by posing the problem as Ax = b, |
---|
| 423 | where x is the set of coefficients we are looking for. |
---|
[3350ad6] | 424 | |
---|
[51f14603] | 425 | Npts is the number of points. |
---|
[3350ad6] | 426 | |
---|
[51f14603] | 427 | In the following i refers to the ith base function coefficient. |
---|
| 428 | The matrix has its entries j in its first Npts rows set to :: |
---|
| 429 | |
---|
| 430 | A[i][j] = (Fourier transformed base function for point j) |
---|
[3350ad6] | 431 | |
---|
[51f14603] | 432 | We them choose a number of r-points, n_r, to evaluate the second |
---|
| 433 | derivative of P(r) at. This is used as our regularization term. |
---|
| 434 | For a vector r of length n_r, the following n_r rows are set to :: |
---|
| 435 | |
---|
| 436 | A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, |
---|
| 437 | evaluated at r[j]) |
---|
[3350ad6] | 438 | |
---|
[51f14603] | 439 | The vector b has its first Npts entries set to :: |
---|
| 440 | |
---|
| 441 | b[j] = (I(q) observed for point j) |
---|
[3350ad6] | 442 | |
---|
[51f14603] | 443 | The following n_r entries are set to zero. |
---|
[3350ad6] | 444 | |
---|
[51f14603] | 445 | The result is found by using scipy.linalg.basic.lstsq to invert |
---|
| 446 | the matrix and find the coefficients x. |
---|
[3350ad6] | 447 | |
---|
[51f14603] | 448 | :param nfunc: number of base functions to use. |
---|
| 449 | :param nr: number of r points to evaluate the 2nd derivative at for the reg. term. |
---|
| 450 | |
---|
| 451 | If the result does not allow us to compute the covariance matrix, |
---|
| 452 | a matrix filled with zeros will be returned. |
---|
| 453 | |
---|
| 454 | """ |
---|
| 455 | # Note: To make sure an array is contiguous: |
---|
[9a5097c] | 456 | # blah = np.ascontiguousarray(blah_original) |
---|
[51f14603] | 457 | # ... before passing it to C |
---|
[3350ad6] | 458 | |
---|
[51f14603] | 459 | if self.is_valid() < 0: |
---|
| 460 | msg = "Invertor: invalid data; incompatible data lengths." |
---|
[574adc7] | 461 | raise RuntimeError(msg) |
---|
[3350ad6] | 462 | |
---|
[51f14603] | 463 | self.nfunc = nfunc |
---|
| 464 | # a -- An M x N matrix. |
---|
| 465 | # b -- An M x nrhs matrix or M vector. |
---|
| 466 | npts = len(self.x) |
---|
[3350ad6] | 467 | nq = nr |
---|
[51f14603] | 468 | sqrt_alpha = math.sqrt(math.fabs(self.alpha)) |
---|
| 469 | if sqrt_alpha < 0.0: |
---|
| 470 | nq = 0 |
---|
| 471 | |
---|
| 472 | # If we need to fit the background, add a term |
---|
[b8080e1] | 473 | if self.est_bck: |
---|
[51f14603] | 474 | nfunc_0 = nfunc |
---|
| 475 | nfunc += 1 |
---|
| 476 | |
---|
[9a5097c] | 477 | a = np.zeros([npts + nq, nfunc]) |
---|
| 478 | b = np.zeros(npts + nq) |
---|
| 479 | err = np.zeros([nfunc, nfunc]) |
---|
[3350ad6] | 480 | |
---|
[51f14603] | 481 | # Construct the a matrix and b vector that represent the problem |
---|
| 482 | t_0 = time.time() |
---|
| 483 | try: |
---|
| 484 | self._get_matrix(nfunc, nq, a, b) |
---|
[d04ac05] | 485 | except Exception as exc: |
---|
| 486 | raise RuntimeError("Invertor: could not invert I(Q)\n %s" % str(exc)) |
---|
[3350ad6] | 487 | |
---|
[51f14603] | 488 | # Perform the inversion (least square fit) |
---|
| 489 | c, chi2, _, _ = lstsq(a, b) |
---|
| 490 | # Sanity check |
---|
| 491 | try: |
---|
| 492 | float(chi2) |
---|
| 493 | except: |
---|
[6da860a] | 494 | chi2 = -1.0 |
---|
[51f14603] | 495 | self.chi2 = chi2 |
---|
[3350ad6] | 496 | |
---|
[9a5097c] | 497 | inv_cov = np.zeros([nfunc, nfunc]) |
---|
[51f14603] | 498 | # Get the covariance matrix, defined as inv_cov = a_transposed * a |
---|
| 499 | self._get_invcov_matrix(nfunc, nr, a, inv_cov) |
---|
[3350ad6] | 500 | |
---|
[51f14603] | 501 | # Compute the reg term size for the output |
---|
| 502 | sum_sig, sum_reg = self._get_reg_size(nfunc, nr, a) |
---|
[3350ad6] | 503 | |
---|
[51f14603] | 504 | if math.fabs(self.alpha) > 0: |
---|
| 505 | new_alpha = sum_sig / (sum_reg / self.alpha) |
---|
| 506 | else: |
---|
| 507 | new_alpha = 0.0 |
---|
| 508 | self.suggested_alpha = new_alpha |
---|
[3350ad6] | 509 | |
---|
[51f14603] | 510 | try: |
---|
[9a5097c] | 511 | cov = np.linalg.pinv(inv_cov) |
---|
[51f14603] | 512 | err = math.fabs(chi2 / float(npts - nfunc)) * cov |
---|
[b8080e1] | 513 | except: |
---|
[51f14603] | 514 | # We were not able to estimate the errors |
---|
| 515 | # Return an empty error matrix |
---|
[b8080e1] | 516 | logger.error(sys.exc_value) |
---|
[3350ad6] | 517 | |
---|
[51f14603] | 518 | # Keep a copy of the last output |
---|
[b8080e1] | 519 | if not self.est_bck: |
---|
[51f14603] | 520 | self.out = c |
---|
| 521 | self.cov = err |
---|
| 522 | else: |
---|
| 523 | self.background = c[0] |
---|
[3350ad6] | 524 | |
---|
[9a5097c] | 525 | err_0 = np.zeros([nfunc, nfunc]) |
---|
| 526 | c_0 = np.zeros(nfunc) |
---|
[3350ad6] | 527 | |
---|
[51f14603] | 528 | for i in range(nfunc_0): |
---|
[3350ad6] | 529 | c_0[i] = c[i + 1] |
---|
[51f14603] | 530 | for j in range(nfunc_0): |
---|
[3350ad6] | 531 | err_0[i][j] = err[i + 1][j + 1] |
---|
| 532 | |
---|
[51f14603] | 533 | self.out = c_0 |
---|
| 534 | self.cov = err_0 |
---|
[3350ad6] | 535 | |
---|
[51f14603] | 536 | # Store computation time |
---|
| 537 | self.elapsed = time.time() - t_0 |
---|
[3350ad6] | 538 | |
---|
[51f14603] | 539 | return self.out, self.cov |
---|
[3350ad6] | 540 | |
---|
[51f14603] | 541 | def estimate_numterms(self, isquit_func=None): |
---|
| 542 | """ |
---|
| 543 | Returns a reasonable guess for the |
---|
| 544 | number of terms |
---|
[3350ad6] | 545 | |
---|
[51f14603] | 546 | :param isquit_func: |
---|
| 547 | reference to thread function to call to check whether the computation needs to |
---|
| 548 | be stopped. |
---|
[3350ad6] | 549 | |
---|
[51f14603] | 550 | :return: number of terms, alpha, message |
---|
[3350ad6] | 551 | |
---|
[51f14603] | 552 | """ |
---|
[b8080e1] | 553 | from .num_term import NTermEstimator |
---|
[5f8fc78] | 554 | estimator = NTermEstimator(self.clone()) |
---|
[51f14603] | 555 | try: |
---|
| 556 | return estimator.num_terms(isquit_func) |
---|
[b8080e1] | 557 | except: |
---|
[51f14603] | 558 | # If we fail, estimate alpha and return the default |
---|
| 559 | # number of terms |
---|
| 560 | best_alpha, _, _ = self.estimate_alpha(self.nfunc) |
---|
[b8080e1] | 561 | logger.warning("Invertor.estimate_numterms: %s" % sys.exc_value) |
---|
[51f14603] | 562 | return self.nfunc, best_alpha, "Could not estimate number of terms" |
---|
[3350ad6] | 563 | |
---|
[51f14603] | 564 | def estimate_alpha(self, nfunc): |
---|
| 565 | """ |
---|
| 566 | Returns a reasonable guess for the |
---|
| 567 | regularization constant alpha |
---|
[3350ad6] | 568 | |
---|
[51f14603] | 569 | :param nfunc: number of terms to use in the expansion. |
---|
[3350ad6] | 570 | |
---|
[51f14603] | 571 | :return: alpha, message, elapsed |
---|
[3350ad6] | 572 | |
---|
[51f14603] | 573 | where alpha is the estimate for alpha, |
---|
| 574 | message is a message for the user, |
---|
| 575 | elapsed is the computation time |
---|
| 576 | """ |
---|
| 577 | #import time |
---|
| 578 | try: |
---|
| 579 | pr = self.clone() |
---|
[3350ad6] | 580 | |
---|
[51f14603] | 581 | # T_0 for computation time |
---|
| 582 | starttime = time.time() |
---|
| 583 | elapsed = 0 |
---|
[3350ad6] | 584 | |
---|
[51f14603] | 585 | # If the current alpha is zero, try |
---|
| 586 | # another value |
---|
| 587 | if pr.alpha <= 0: |
---|
| 588 | pr.alpha = 0.0001 |
---|
[3350ad6] | 589 | |
---|
[51f14603] | 590 | # Perform inversion to find the largest alpha |
---|
| 591 | out, _ = pr.invert(nfunc) |
---|
| 592 | elapsed = time.time() - starttime |
---|
| 593 | initial_alpha = pr.alpha |
---|
| 594 | initial_peaks = pr.get_peaks(out) |
---|
[3350ad6] | 595 | |
---|
[51f14603] | 596 | # Try the inversion with the estimated alpha |
---|
| 597 | pr.alpha = pr.suggested_alpha |
---|
| 598 | out, _ = pr.invert(nfunc) |
---|
[3350ad6] | 599 | |
---|
[51f14603] | 600 | npeaks = pr.get_peaks(out) |
---|
| 601 | # if more than one peak to start with |
---|
| 602 | # just return the estimate |
---|
| 603 | if npeaks > 1: |
---|
| 604 | #message = "Your P(r) is not smooth, |
---|
| 605 | #please check your inversion parameters" |
---|
| 606 | message = None |
---|
| 607 | return pr.suggested_alpha, message, elapsed |
---|
| 608 | else: |
---|
[3350ad6] | 609 | |
---|
[51f14603] | 610 | # Look at smaller values |
---|
| 611 | # We assume that for the suggested alpha, we have 1 peak |
---|
| 612 | # if not, send a message to change parameters |
---|
| 613 | alpha = pr.suggested_alpha |
---|
| 614 | best_alpha = pr.suggested_alpha |
---|
| 615 | found = False |
---|
| 616 | for i in range(10): |
---|
[3350ad6] | 617 | pr.alpha = (0.33) ** (i + 1) * alpha |
---|
[51f14603] | 618 | out, _ = pr.invert(nfunc) |
---|
[3350ad6] | 619 | |
---|
[51f14603] | 620 | peaks = pr.get_peaks(out) |
---|
| 621 | if peaks > 1: |
---|
| 622 | found = True |
---|
| 623 | break |
---|
| 624 | best_alpha = pr.alpha |
---|
[3350ad6] | 625 | |
---|
[51f14603] | 626 | # If we didn't find a turning point for alpha and |
---|
| 627 | # the initial alpha already had only one peak, |
---|
| 628 | # just return that |
---|
| 629 | if not found and initial_peaks == 1 and \ |
---|
| 630 | initial_alpha < best_alpha: |
---|
| 631 | best_alpha = initial_alpha |
---|
[3350ad6] | 632 | |
---|
[51f14603] | 633 | # Check whether the size makes sense |
---|
| 634 | message = '' |
---|
[3350ad6] | 635 | |
---|
[51f14603] | 636 | if not found: |
---|
| 637 | message = None |
---|
| 638 | elif best_alpha >= 0.5 * pr.suggested_alpha: |
---|
| 639 | # best alpha is too big, return a |
---|
| 640 | # reasonable value |
---|
[3350ad6] | 641 | message = "The estimated alpha for your system is too " |
---|
[51f14603] | 642 | message += "large. " |
---|
| 643 | message += "Try increasing your maximum distance." |
---|
[3350ad6] | 644 | |
---|
[51f14603] | 645 | return best_alpha, message, elapsed |
---|
[3350ad6] | 646 | |
---|
[b8080e1] | 647 | except: |
---|
| 648 | message = "Invertor.estimate_alpha: %s" % sys.exc_value |
---|
[51f14603] | 649 | return 0, message, elapsed |
---|
[3350ad6] | 650 | |
---|
[51f14603] | 651 | def to_file(self, path, npts=100): |
---|
| 652 | """ |
---|
| 653 | Save the state to a file that will be readable |
---|
| 654 | by SliceView. |
---|
[3350ad6] | 655 | |
---|
[51f14603] | 656 | :param path: path of the file to write |
---|
| 657 | :param npts: number of P(r) points to be written |
---|
[3350ad6] | 658 | |
---|
[51f14603] | 659 | """ |
---|
| 660 | file = open(path, 'w') |
---|
| 661 | file.write("#d_max=%g\n" % self.d_max) |
---|
| 662 | file.write("#nfunc=%g\n" % self.nfunc) |
---|
| 663 | file.write("#alpha=%g\n" % self.alpha) |
---|
| 664 | file.write("#chi2=%g\n" % self.chi2) |
---|
| 665 | file.write("#elapsed=%g\n" % self.elapsed) |
---|
| 666 | file.write("#qmin=%s\n" % str(self.q_min)) |
---|
| 667 | file.write("#qmax=%s\n" % str(self.q_max)) |
---|
| 668 | file.write("#slit_height=%g\n" % self.slit_height) |
---|
| 669 | file.write("#slit_width=%g\n" % self.slit_width) |
---|
| 670 | file.write("#background=%g\n" % self.background) |
---|
[b8080e1] | 671 | if self.est_bck: |
---|
[51f14603] | 672 | file.write("#has_bck=1\n") |
---|
| 673 | else: |
---|
| 674 | file.write("#has_bck=0\n") |
---|
| 675 | file.write("#alpha_estimate=%g\n" % self.suggested_alpha) |
---|
[45dffa69] | 676 | if self.out is not None: |
---|
[51f14603] | 677 | if len(self.out) == len(self.cov): |
---|
| 678 | for i in range(len(self.out)): |
---|
| 679 | file.write("#C_%i=%s+-%s\n" % (i, str(self.out[i]), |
---|
[3350ad6] | 680 | str(self.cov[i][i]))) |
---|
[51f14603] | 681 | file.write("<r> <Pr> <dPr>\n") |
---|
[9a5097c] | 682 | r = np.arange(0.0, self.d_max, self.d_max / npts) |
---|
[3350ad6] | 683 | |
---|
[51f14603] | 684 | for r_i in r: |
---|
| 685 | (value, err) = self.pr_err(self.out, self.cov, r_i) |
---|
| 686 | file.write("%g %g %g\n" % (r_i, value, err)) |
---|
[3350ad6] | 687 | |
---|
[51f14603] | 688 | file.close() |
---|
[3350ad6] | 689 | |
---|
[51f14603] | 690 | def from_file(self, path): |
---|
| 691 | """ |
---|
| 692 | Load the state of the Invertor from a file, |
---|
| 693 | to be able to generate P(r) from a set of |
---|
| 694 | parameters. |
---|
[3350ad6] | 695 | |
---|
[51f14603] | 696 | :param path: path of the file to load |
---|
[3350ad6] | 697 | |
---|
[51f14603] | 698 | """ |
---|
| 699 | #import os |
---|
| 700 | #import re |
---|
| 701 | if os.path.isfile(path): |
---|
| 702 | try: |
---|
| 703 | fd = open(path, 'r') |
---|
[3350ad6] | 704 | |
---|
[51f14603] | 705 | buff = fd.read() |
---|
| 706 | lines = buff.split('\n') |
---|
| 707 | for line in lines: |
---|
| 708 | if line.startswith('#d_max='): |
---|
| 709 | toks = line.split('=') |
---|
| 710 | self.d_max = float(toks[1]) |
---|
| 711 | elif line.startswith('#nfunc='): |
---|
| 712 | toks = line.split('=') |
---|
| 713 | self.nfunc = int(toks[1]) |
---|
[9a5097c] | 714 | self.out = np.zeros(self.nfunc) |
---|
| 715 | self.cov = np.zeros([self.nfunc, self.nfunc]) |
---|
[51f14603] | 716 | elif line.startswith('#alpha='): |
---|
| 717 | toks = line.split('=') |
---|
| 718 | self.alpha = float(toks[1]) |
---|
| 719 | elif line.startswith('#chi2='): |
---|
| 720 | toks = line.split('=') |
---|
| 721 | self.chi2 = float(toks[1]) |
---|
| 722 | elif line.startswith('#elapsed='): |
---|
| 723 | toks = line.split('=') |
---|
| 724 | self.elapsed = float(toks[1]) |
---|
| 725 | elif line.startswith('#alpha_estimate='): |
---|
| 726 | toks = line.split('=') |
---|
| 727 | self.suggested_alpha = float(toks[1]) |
---|
| 728 | elif line.startswith('#qmin='): |
---|
| 729 | toks = line.split('=') |
---|
| 730 | try: |
---|
| 731 | self.q_min = float(toks[1]) |
---|
| 732 | except: |
---|
| 733 | self.q_min = None |
---|
| 734 | elif line.startswith('#qmax='): |
---|
| 735 | toks = line.split('=') |
---|
| 736 | try: |
---|
| 737 | self.q_max = float(toks[1]) |
---|
| 738 | except: |
---|
| 739 | self.q_max = None |
---|
| 740 | elif line.startswith('#slit_height='): |
---|
| 741 | toks = line.split('=') |
---|
| 742 | self.slit_height = float(toks[1]) |
---|
| 743 | elif line.startswith('#slit_width='): |
---|
| 744 | toks = line.split('=') |
---|
| 745 | self.slit_width = float(toks[1]) |
---|
| 746 | elif line.startswith('#background='): |
---|
| 747 | toks = line.split('=') |
---|
| 748 | self.background = float(toks[1]) |
---|
| 749 | elif line.startswith('#has_bck='): |
---|
| 750 | toks = line.split('=') |
---|
[b8080e1] | 751 | self.est_bck = int(toks[1]) == 1 |
---|
[3350ad6] | 752 | |
---|
[51f14603] | 753 | # Now read in the parameters |
---|
| 754 | elif line.startswith('#C_'): |
---|
| 755 | toks = line.split('=') |
---|
| 756 | p = re.compile('#C_([0-9]+)') |
---|
| 757 | m = p.search(toks[0]) |
---|
| 758 | toks2 = toks[1].split('+-') |
---|
| 759 | i = int(m.group(1)) |
---|
| 760 | self.out[i] = float(toks2[0]) |
---|
[3350ad6] | 761 | |
---|
[51f14603] | 762 | self.cov[i][i] = float(toks2[1]) |
---|
[3350ad6] | 763 | |
---|
[b8080e1] | 764 | except: |
---|
| 765 | msg = "Invertor.from_file: corrupted file\n%s" % sys.exc_value |
---|
[574adc7] | 766 | raise RuntimeError(msg) |
---|
[51f14603] | 767 | else: |
---|
| 768 | msg = "Invertor.from_file: '%s' is not a file" % str(path) |
---|
[574adc7] | 769 | raise RuntimeError(msg) |
---|