Changeset 1db4a53 in sasview for pr_inversion


Ignore:
Timestamp:
Nov 15, 2010 3: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

Location:
pr_inversion
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • pr_inversion/distance_explorer.py

    rd84a90c r1db4a53  
    5656        """ 
    5757        self.pr_state  = pr_state 
    58         self._default_min = 0.8*self.pr_state.d_max 
    59         self._default_max = 1.2*self.pr_state.d_max 
     58        self._default_min = 0.8 * self.pr_state.d_max 
     59        self._default_max = 1.2 * self.pr_state.d_max 
    6060 
    6161         
    … …  
    103103            except: 
    104104                # This inversion failed, skip this D_max value 
    105                 results.errors.append("ExploreDialog: inversion failed for D_max=%s\n %s" % (str(d), sys.exc_value)) 
     105                msg = "ExploreDialog: inversion failed for " 
     106                msg += "D_max=%s\n %s" % (str(d), sys.exc_value) 
     107                results.errors.append(msg) 
    106108             
    107109        return results 

  • 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() 

  • pr_inversion/num_term.py

    rd84a90c r1db4a53  
    1 import unittest, math, numpy, sys, string 
     1 
     2#import unittest 
     3import math 
     4import numpy 
     5#import sys 
     6#import string 
     7import copy 
    28from sans.pr.invertor import Invertor 
    39 
    … …  
    1117        self.nterm_min = 10 
    1218        self.nterm_max = len(self.invertor.x) 
    13         if self.nterm_max>50: 
    14             self.nterm_max=50 
     19        if self.nterm_max > 50: 
     20            self.nterm_max = 50 
    1521        self.isquit_func = None 
    1622          
    … …  
    2531        """ 
    2632        """ 
    27         return bool(n%2) 
     33        return bool(n % 2) 
    2834 
    2935    def sort_osc(self): 
    3036        """ 
    3137        """ 
    32         import copy 
     38        #import copy 
    3339        osc = copy.deepcopy(self.dataset) 
    3440        lis = [] 
    … …  
    6975            osc = inver.oscillations(inver.out) 
    7076            err = inver.get_pos_err(inver.out, inver.cov) 
    71             if osc>10.0: 
     77            if osc > 10.0: 
    7278                break 
    7379            self.osc_list.append(osc) 
    … …  
    7985        #print "err", self.err_list 
    8086        #print "alpha", self.alpha_list 
    81         new_ls = [] 
     87        #new_ls = [] 
    8288        new_osc1 = [] 
    83         new_osc2= [] 
     89        new_osc2 = [] 
    8490        new_osc3 = [] 
    85         flag9=False 
    86         flag8=False 
    87         flag7=False 
     91        flag9 = False 
     92        flag8 = False 
     93        flag7 = False 
    8894        for i in range(len(self.err_list)): 
    89             if self.err_list[i] <= 1.0 and self.err_list[i] >=0.9: 
     95            if self.err_list[i] <= 1.0 and self.err_list[i] >= 0.9: 
    9096                new_osc1.append(self.osc_list[i]) 
    91                 flag9=True 
    92             if self.err_list[i] < 0.9 and self.err_list[i] >=0.8: 
     97                flag9 = True 
     98            if self.err_list[i] < 0.9 and self.err_list[i] >= 0.8: 
    9399                new_osc2.append(self.osc_list[i]) 
    94                 flag8=True 
    95             if self.err_list[i] <0.8 and self.err_list[i] >= 0.7: 
     100                flag8 = True 
     101            if self.err_list[i] < 0.8 and self.err_list[i] >= 0.7: 
    96102                new_osc3.append(self.osc_list[i]) 
    97                 falg7=True 
     103                falg7 = True 
    98104                  
    99         if flag9==True: 
     105        if flag9 == True: 
    100106            self.dataset = new_osc1 
    101         elif flag8==True: 
     107        elif flag8 == True: 
    102108            self.dataset = new_osc2 
    103109        else: 
    … …  
    159165#For testing 
    160166def load(path): 
    161     import numpy, math, sys 
     167    #import numpy, math, sys 
    162168    # Read the data from the data file 
    163169    data_x   = numpy.zeros(0) 
    … …  
    173179            try: 
    174180                toks = line.split() 
    175                 x = float(toks[0]) 
    176                 y = float(toks[1]) 
    177                 if len(toks)>2: 
     181                test_x = float(toks[0]) 
     182                test_y = float(toks[1]) 
     183                if len(toks) > 2: 
    178184                    err = float(toks[2]) 
    179185                else: 
    180                     if scale==None: 
    181                         scale = 0.05*math.sqrt(y) 
     186                    if scale == None: 
     187                        scale = 0.05 * math.sqrt(test_y) 
    182188                        #scale = 0.05/math.sqrt(y) 
    183189                        min_err = 0.01*y 
    184                     err = scale*math.sqrt(y)+min_err 
     190                    err = scale * math.sqrt(test_y) + min_err 
    185191                    #err = 0 
    186192                     
    187                 data_x = numpy.append(data_x, x) 
    188                 data_y = numpy.append(data_y, y) 
     193                data_x = numpy.append(data_x, test_x) 
     194                data_y = numpy.append(data_y, test_y) 
    189195                data_err = numpy.append(data_err, err) 
    190196            except: 
Note: See TracChangeset for help on using the changeset viewer.