source: sasview/src/sas/sascalc/pr/invertor.py @ 57a91fc

ESS_GUIESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalc
Last change on this file since 57a91fc was 57a91fc, checked in by Piotr Rozyczko <piotr.rozyczko@…>, 5 years ago

cherrypicked changes from master to sascalc (except PR#190)

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