source: sasview/pr_inversion/src/sans/pr/invertor.py @ aa9ea2e

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since aa9ea2e was 007f9b3, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Re #5 Fix p(r) unit tests

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