source: sasview/src/sas/pr/invertor.py @ c76bdc3

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 c76bdc3 was c76bdc3, checked in by Doucet, Mathieu <doucetm@…>, 10 years ago

pylint fixes

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