source: sasview/src/sas/sascalc/pr/invertor.py @ 2469df7

magnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 2469df7 was 2469df7, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

lint: update 'if x==True/False?' to 'if x/not x:'

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