source: sasview/src/sas/sascalc/pr/invertor.py @ 13da5f5

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

Added error check to invert fucntion

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