source: sasview/src/sas/sascalc/pr/invertor.py @ 6701a0b

ESS_GUIESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalc
Last change on this file since 6701a0b was 6701a0b, checked in by wojciech, 6 years ago

Typos fixed

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