Changeset 1db4a53 in sasview for pr_inversion/invertor.py


Ignore:
Timestamp:
Nov 15, 2010 5:31:33 PM (14 years ago)
Author:
Gervaise Alina <gervyh@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
519c693
Parents:
6996942
Message:

working on pylint

File:
1 edited

Legend:

Unmodified
Added
Removed

  • pr_inversion/invertor.py

    rd84a90c r1db4a53  
    33The module contains the Invertor class. 
    44""" 
    5 from sans.pr.core.pr_inversion import Cinvertor 
     5 
    66import numpy 
    77import sys 
    8 import math, time 
     8import math 
     9import time 
     10import copy 
     11import os 
     12import re 
    913from numpy.linalg import lstsq 
     14from scipy import optimize 
     15from sans.pr.core.pr_inversion import Cinvertor 
    1016 
    1117def help(): 
    … …  
    1420    Future work: extend this function to allow topic selection 
    1521    """ 
    16     info_txt  = "The inversion approach is based on Moore, J. Appl. Cryst. (1980) 13, 168-175.\n\n" 
    17     info_txt += "P(r) is set to be equal to an expansion of base functions of the type " 
    18     info_txt += "phi_n(r) = 2*r*sin(pi*n*r/D_max). The coefficient of each base functions " 
    19     info_txt += "in the expansion is found by performing a least square fit with the " 
     22    info_txt  = "The inversion approach is based on Moore, J. Appl. Cryst. "  
     23    info_txt += "(1980) 13, 168-175.\n\n" 
     24    info_txt += "P(r) is set to be equal to an expansion of base functions " 
     25    info_txt += "of the type " 
     26    info_txt += "phi_n(r) = 2*r*sin(pi*n*r/D_max). The coefficient of each " 
     27    info_txt += "base functions " 
     28    info_txt += "in the expansion is found by performing a least square fit " 
     29    info_txt += "with the " 
    2030    info_txt += "following fit function:\n\n" 
    21     info_txt += "chi**2 = sum_i[ I_meas(q_i) - I_th(q_i) ]**2/error**2 + Reg_term\n\n" 
    22     info_txt += "where I_meas(q) is the measured scattering intensity and I_th(q) is " 
    23     info_txt += "the prediction from the Fourier transform of the P(r) expansion. " 
    24     info_txt += "The Reg_term term is a regularization term set to the second derivative " 
    25     info_txt += "d**2P(r)/dr**2 integrated over r. It is used to produce a smooth P(r) output.\n\n" 
     31    info_txt += "chi**2 = sum_i[ I_meas(q_i) - I_th(q_i) ]**2/error**2 +" 
     32    info_txt += "Reg_term\n\n" 
     33    info_txt += "where I_meas(q) is the measured scattering intensity and " 
     34    info_txt += "I_th(q) is " 
     35    info_txt += "the prediction from the Fourier transform of the P(r) " 
     36    info_txt += "expansion. " 
     37    info_txt += "The Reg_term term is a regularization term set to the second" 
     38    info_txt += " derivative " 
     39    info_txt += "d**2P(r)/dr**2 integrated over r. It is used to produce " 
     40    info_txt += "a smooth P(r) output.\n\n" 
    2641    info_txt += "The following are user inputs:\n\n" 
    27     info_txt += "   - Number of terms: the number of base functions in the P(r) expansion.\n\n" 
    28     info_txt += "   - Regularization constant: a multiplicative constant to set the size of " 
     42    info_txt += "   - Number of terms: the number of base functions in the P(r)" 
     43    info_txt += " expansion.\n\n" 
     44    info_txt += "   - Regularization constant: a multiplicative constant " 
     45    info_txt += "to set the size of " 
    2946    info_txt += "the regularization term.\n\n" 
    30     info_txt += "   - Maximum distance: the maximum distance between any two points in the system.\n" 
     47    info_txt += "   - Maximum distance: the maximum distance between any " 
     48    info_txt += "two points in the system.\n" 
    3149      
    3250    return info_txt 
    … …  
    4967    derivative of P(r) at. This is used as our regularization term. 
    5068    For a vector r of length n_r, the following n_r rows are set to 
    51         A[j+Npts][i] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, evaluated at r[j]) 
     69        A[j+Npts][i] = (2nd derivative of P(r), d**2(P(r))/d(r)**2,  
     70        evaluated at r[j]) 
    5271         
    5372    The vector b has its first Npts entries set to 
    … …  
    92111        x, y, err, d_max, q_min, q_max and alpha 
    93112        """ 
    94         if   name=='x': 
     113        if   name == 'x': 
    95114            if 0.0 in value: 
    96                 raise ValueError, "Invertor: one of your q-values is zero. Delete that entry before proceeding" 
     115                msg = "Invertor: one of your q-values is zero. " 
     116                msg += "Delete that entry before proceeding" 
     117                raise ValueError, msg 
    97118            return self.set_x(value) 
    98         elif name=='y': 
     119        elif name == 'y': 
    99120            return self.set_y(value) 
    100         elif name=='err': 
     121        elif name == 'err': 
    101122            value2 = abs(value) 
    102123            return self.set_err(value2) 
    103         elif name=='d_max': 
     124        elif name == 'd_max': 
    104125            return self.set_dmax(value) 
    105         elif name=='q_min': 
    106             if value==None: 
     126        elif name == 'q_min': 
     127            if value == None: 
    107128                return self.set_qmin(-1.0) 
    108129            return self.set_qmin(value) 
    109         elif name=='q_max': 
    110             if value==None: 
     130        elif name == 'q_max': 
     131            if value == None: 
    111132                return self.set_qmax(-1.0) 
    112133            return self.set_qmax(value) 
    113         elif name=='alpha': 
     134        elif name == 'alpha': 
    114135            return self.set_alpha(value) 
    115         elif name=='slit_height': 
     136        elif name == 'slit_height': 
    116137            return self.set_slit_height(value) 
    117         elif name=='slit_width': 
     138        elif name == 'slit_width': 
    118139            return self.set_slit_width(value) 
    119         elif name=='has_bck': 
    120             if value==True: 
     140        elif name == 'has_bck': 
     141            if value == True: 
    121142                return self.set_has_bck(1) 
    122             elif value==False: 
     143            elif value == False: 
    123144                return self.set_has_bck(0) 
    124145            else: 
    … …  
    131152        Return the value of an attribute 
    132153        """ 
    133         import numpy 
    134         if   name=='x': 
     154        #import numpy 
     155        if name == 'x': 
    135156            out = numpy.ones(self.get_nx()) 
    136157            self.get_x(out) 
    137158            return out 
    138         elif name=='y': 
     159        elif name == 'y': 
    139160            out = numpy.ones(self.get_ny()) 
    140161            self.get_y(out) 
    141162            return out 
    142         elif name=='err': 
     163        elif name == 'err': 
    143164            out = numpy.ones(self.get_nerr()) 
    144165            self.get_err(out) 
    145166            return out 
    146         elif name=='d_max': 
     167        elif name == 'd_max': 
    147168            return self.get_dmax() 
    148         elif name=='q_min': 
     169        elif name == 'q_min': 
    149170            qmin = self.get_qmin() 
    150             if qmin<0: 
     171            if qmin < 0: 
    151172                return None 
    152173            return qmin 
    153         elif name=='q_max': 
     174        elif name == 'q_max': 
    154175            qmax = self.get_qmax() 
    155             if qmax<0: 
     176            if qmax < 0: 
    156177                return None 
    157178            return qmax 
    158         elif name=='alpha': 
     179        elif name == 'alpha': 
    159180            return self.get_alpha() 
    160         elif name=='slit_height': 
     181        elif name == 'slit_height': 
    161182            return self.get_slit_height() 
    162         elif name=='slit_width': 
     183        elif name == 'slit_width': 
    163184            return self.get_slit_width() 
    164         elif name=='has_bck': 
     185        elif name == 'has_bck': 
    165186            value = self.get_has_bck() 
    166             if value==1: 
     187            if value == 1: 
    167188                return True 
    168189            else: 
    … …  
    176197        Return a clone of this instance 
    177198        """ 
    178         import copy 
     199        #import copy 
    179200         
    180201        invertor = Invertor() 
    … …  
    241262         
    242263        """ 
    243         return Cinvertor.iq(self, out, q)+self.background 
     264        return Cinvertor.iq(self, out, q) + self.background 
    244265     
    245266    def invert_optimize(self, nfunc=10, nr=20): 
    … …  
    248269         
    249270        This probably produce more reliable results, but is much slower. 
    250         The minimization function is set to sum_i[ (I_obs(q_i) - I_theo(q_i))/err**2 ] + alpha * reg_term, 
    251         where the reg_term is given by Svergun: it is the integral of the square of the first derivative 
     271        The minimization function is set to 
     272        sum_i[ (I_obs(q_i) - I_theo(q_i))/err**2 ] + alpha * reg_term, 
     273        where the reg_term is given by Svergun: it is the integral of  
     274        the square of the first derivative 
    252275        of P(r), d(P(r))/dr, integrated over the full range of r. 
    253276         
    254277        :param nfunc: number of base functions to use. 
    255         :param nr: number of r points to evaluate the 2nd derivative at for the reg. term. 
     278        :param nr: number of r points to evaluate the 2nd derivative at  
     279            for the reg. term. 
    256280         
    257281        :return: c_out, c_cov - the coefficients with covariance matrix  
    … …  
    259283        """ 
    260284         
    261         from scipy import optimize 
    262         import time 
     285        #from scipy import optimize 
     286        #import time 
    263287         
    264288        self.nfunc = nfunc 
    265289        # First, check that the current data is valid 
    266         if self.is_valid()<=0: 
    267             raise RuntimeError, "Invertor.invert: Data array are of different length" 
     290        if self.is_valid() <= 0: 
     291            msg = "Invertor.invert: Data array are of different length" 
     292            raise RuntimeError, msg 
    268293         
    269294        p = numpy.ones(nfunc) 
    270295        t_0 = time.time() 
    271         out, cov_x, info, mesg, success = optimize.leastsq(self.residuals, p, full_output=1, warning=True) 
     296        out, cov_x, _, _, _ = optimize.leastsq(self.residuals, 
     297                                                            p, full_output=1, 
     298                                                             warning=True) 
    272299         
    273300        # Compute chi^2 
    … …  
    292319        This method is provided as a test.  
    293320        """ 
    294         from scipy import optimize 
     321        #from scipy import optimize 
    295322         
    296323        # First, check that the current data is valid 
    297         if self.is_valid()<=0: 
    298             raise RuntimeError, "Invertor.invert: Data arrays are of different length" 
     324        if self.is_valid() <= 0: 
     325            msg = "Invertor.invert: Data arrays are of different length" 
     326            raise RuntimeError, msg 
    299327         
    300328        p = numpy.ones(nfunc) 
    301329        t_0 = time.time() 
    302         out, cov_x, info, mesg, success = optimize.leastsq(self.pr_residuals, p, full_output=1, warning=True) 
     330        out, cov_x, info, mesg, success = optimize.leastsq(self.pr_residuals, p, 
     331                                                            full_output=1, 
     332                                                            warning=True) 
    303333         
    304334        # Compute chi^2 
    … …  
    333363        Check q-value against user-defined range 
    334364        """ 
    335         if not self.q_min==None and q<self.q_min: 
     365        if not self.q_min == None and q < self.q_min: 
    336366            return False 
    337         if not self.q_max==None and q>self.q_max: 
     367        if not self.q_max == None and q > self.q_max: 
    338368            return False 
    339369        return True 
    … …  
    353383        derivative of P(r) at. This is used as our regularization term. 
    354384        For a vector r of length n_r, the following n_r rows are set to 
    355             A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, evaluated at r[j]) 
     385            A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, 
     386            evaluated at r[j]) 
    356387             
    357388        The vector b has its first Npts entries set to 
    … …  
    364395         
    365396        :param nfunc: number of base functions to use. 
    366         :param nr: number of r points to evaluate the 2nd derivative at for the reg. term. 
     397        :param nr: number of r points to evaluate the 2nd derivative at 
     398            for the reg. term. 
    367399 
    368400        If the result does not allow us to compute the covariance matrix, 
    … …  
    374406        # ... before passing it to C 
    375407         
    376         if self.is_valid()<0: 
    377             raise RuntimeError, "Invertor: invalid data; incompatible data lengths." 
     408        if self.is_valid() < 0: 
     409            msg = "Invertor: invalid data; incompatible data lengths." 
     410            raise RuntimeError, msg 
    378411         
    379412        self.nfunc = nfunc 
    … …  
    383416        nq   = nr 
    384417        sqrt_alpha = math.sqrt(math.fabs(self.alpha)) 
    385         if sqrt_alpha<0.0: 
     418        if sqrt_alpha < 0.0: 
    386419            nq = 0 
    387420 
    388421        # If we need to fit the background, add a term 
    389         if self.has_bck==True: 
     422        if self.has_bck == True: 
    390423            nfunc_0 = nfunc 
    391424            nfunc += 1 
    … …  
    415448        sum_sig, sum_reg = self._get_reg_size(nfunc, nr, a) 
    416449                     
    417         if math.fabs(self.alpha)>0: 
     450        if math.fabs(self.alpha) > 0: 
    418451            new_alpha = sum_sig/(sum_reg/self.alpha) 
    419452        else: 
    … …  
    430463             
    431464        # Keep a copy of the last output 
    432         if self.has_bck==False: 
     465        if self.has_bck == False: 
    433466            self.background = 0 
    434467            self.out = c 
    … …  
    485518        elapsed is the computation time 
    486519        """ 
    487         import time 
     520        #import time 
    488521        try:             
    489522            pr = self.clone() 
    … …  
    495528            # If the current alpha is zero, try 
    496529            # another value 
    497             if pr.alpha<=0: 
     530            if pr.alpha <= 0: 
    498531                pr.alpha = 0.0001 
    499532                  
    500533            # Perform inversion to find the largest alpha 
    501534            out, cov = pr.invert(nfunc) 
    502             elapsed = time.time()-starttime 
     535            elapsed = time.time() - starttime 
    503536            initial_alpha = pr.alpha 
    504537            initial_peaks = pr.get_peaks(out) 
    … …  
    511544            # if more than one peak to start with 
    512545            # just return the estimate 
    513             if npeaks>1: 
    514                 #message = "Your P(r) is not smooth, please check your inversion parameters" 
     546            if npeaks > 1: 
     547                #message = "Your P(r) is not smooth,  
     548                #please check your inversion parameters" 
    515549                message = None 
    516550                return pr.suggested_alpha, message, elapsed 
    … …  
    524558                found = False 
    525559                for i in range(10): 
    526                     pr.alpha = (0.33)**(i+1)*alpha 
     560                    pr.alpha = (0.33)**(i+1) * alpha 
    527561                    out, cov = pr.invert(nfunc) 
    528562                     
    529563                    peaks = pr.get_peaks(out) 
    530                     if peaks>1: 
     564                    if peaks > 1: 
    531565                        found = True 
    532566                        break 
    … …  
    536570                # the initial alpha already had only one peak, 
    537571                # just return that 
    538                 if not found and initial_peaks==1 and initial_alpha<best_alpha: 
     572                if not found and initial_peaks == 1 and \ 
     573                    initial_alpha<best_alpha: 
    539574                    best_alpha = initial_alpha 
    540575                     
    541576                # Check whether the size makes sense 
    542                 message='' 
     577                message = '' 
    543578                 
    544579                if not found: 
    545580                    message = None 
    546                 elif best_alpha>=0.5*pr.suggested_alpha: 
     581                elif best_alpha >= 0.5 * pr.suggested_alpha: 
    547582                    # best alpha is too big, return a  
    548583                    # reasonable value 
    549                     message  = "The estimated alpha for your system is too large. " 
     584                    message  = "The estimated alpha for your system is too " 
     585                    messsage += "large. " 
    550586                    message += "Try increasing your maximum distance." 
    551587                 
    … …  
    577613        file.write("#slit_width=%g\n" % self.slit_width) 
    578614        file.write("#background=%g\n" % self.background) 
    579         if self.has_bck==True: 
     615        if self.has_bck == True: 
    580616            file.write("#has_bck=1\n") 
    581617        else: 
    582618            file.write("#has_bck=0\n") 
    583619        file.write("#alpha_estimate=%g\n" % self.suggested_alpha) 
    584         if not self.out==None: 
    585             if len(self.out)==len(self.cov): 
     620        if not self.out == None: 
     621            if len(self.out) == len(self.cov): 
    586622                for i in range(len(self.out)): 
    587                     file.write("#C_%i=%s+-%s\n" % (i, str(self.out[i]), str(self.cov[i][i]))) 
     623                    file.write("#C_%i=%s+-%s\n" % (i, str(self.out[i]), 
     624                                                    str(self.cov[i][i]))) 
    588625        file.write("<r>  <Pr>  <dPr>\n") 
    589626        r = numpy.arange(0.0, self.d_max, self.d_max/npts) 
    … …  
    605642         
    606643        """ 
    607         import os 
    608         import re 
     644        #import os 
     645        #import re 
    609646        if os.path.isfile(path): 
    610647            try: 
    … …  
    657694                    elif line.startswith('#has_bck='): 
    658695                        toks = line.split('=') 
    659                         if int(toks[1])==1: 
    660                             self.has_bck=True 
     696                        if int(toks[1]) == 1: 
     697                            self.has_bck = True 
    661698                        else: 
    662                             self.has_bck=False 
     699                            self.has_bck = False 
    663700             
    664701                    # Now read in the parameters 
    … …  
    674711             
    675712            except: 
    676                 raise RuntimeError, "Invertor.from_file: corrupted file\n%s" % sys.exc_value 
     713                msg = "Invertor.from_file: corrupted file\n%s" % sys.exc_value 
     714                raise RuntimeError, msg 
    677715        else: 
    678             raise RuntimeError, "Invertor.from_file: '%s' is not a file" % str(path)  
    679          
    680          
    681      
    682      
     716            msg = "Invertor.from_file: '%s' is not a file" % str(path)  
     717            raise RuntimeError, msg 
     718         
     719          
    683720if __name__ == "__main__": 
    684721    o = Invertor() 
Note: See TracChangeset for help on using the changeset viewer.