source: sasview/src/sas/sascalc/pr/invertor.py @ 4b7d322

ESS_GUIESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalc
Last change on this file since 4b7d322 was eeea6a3, checked in by wojciech, 6 years ago

Simulated error handling added to more logical place

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