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

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

Typos fixed

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