Changeset 34f3ad0 in sasview for pr_inversion/src/sans/pr


Ignore:
Timestamp:
Apr 27, 2012 9:07:00 AM (13 years ago)
Author:
Mathieu Doucet <doucetm@…>
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:
8a621ac
Parents:
1f9f3c8a
Message:

make pylint happier

Location:
pr_inversion/src/sans/pr
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • pr_inversion/src/sans/pr/distance_explorer.py

    rd9621b17 r34f3ad0  
    1616over that range. 
    1717""" 
     18import sys 
    1819 
    19 import numpy 
    20 import math 
    21 import sys 
    2220 
    2321class Results: 
    … …  
    4038        self.bck = [] 
    4139        self.d_max = [] 
    42         ## List of errors found during the last exploration         
     40        ## List of errors found during the last exploration 
    4341        self.errors = [] 
     42         
    4443         
    4544class DistExplorer(object): 
    … …  
    5554         
    5655        """ 
    57         self.pr_state  = pr_state 
     56        self.pr_state = pr_state 
    5857        self._default_min = 0.8 * self.pr_state.d_max 
    5958        self._default_max = 1.2 * self.pr_state.d_max 
    60  
    6159         
    6260    def __call__(self, dmin=None, dmax=None, npts=10): 
    … …  
    8078         
    8179        # Loop over d_max values 
    82         for i in range(npts):     
    83             d = dmin + i * (dmax - dmin)/(npts-1.0) 
     80        for i in range(npts): 
     81            d = dmin + i * (dmax - dmin) / (npts - 1.0) 
    8482            self.pr_state.d_max = d 
    8583            try: 
    86                 out, cov = self.pr_state.invert(self.pr_state.nfunc)     
     84                out, cov = self.pr_state.invert(self.pr_state.nfunc) 
    8785             
    8886                # Store results 
    … …  
    10098                results.pos.append(pos) 
    10199                results.pos_err.append(pos_err) 
    102                 results.osc.append(osc)            
     100                results.osc.append(osc) 
    103101            except: 
    104102                # This inversion failed, skip this D_max value 
    … …  
    108106             
    109107        return results 
    110          

  • pr_inversion/src/sans/pr/invertor.py

    rea59943 r34f3ad0  
    2020    Future work: extend this function to allow topic selection 
    2121    """ 
    22     info_txt  = "The inversion approach is based on Moore, J. Appl. Cryst. "  
     22    info_txt  = "The inversion approach is based on Moore, J. Appl. Cryst. " 
    2323    info_txt += "(1980) 13, 168-175.\n\n" 
    2424    info_txt += "P(r) is set to be equal to an expansion of base functions " 
    … …  
    6262    In the following i refers to the ith base function coefficient. 
    6363    The matrix has its entries j in its first Npts rows set to 
    64         A[j][i] = (Fourier transformed base function for point j)  
     64        A[j][i] = (Fourier transformed base function for point j) 
    6565         
    6666    We them choose a number of r-points, n_r, to evaluate the second 
    6767    derivative of P(r) at. This is used as our regularization term. 
    6868    For a vector r of length n_r, the following n_r rows are set to 
    69         A[j+Npts][i] = (2nd derivative of P(r), d**2(P(r))/d(r)**2,  
     69        A[j+Npts][i] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, 
    7070        evaluated at r[j]) 
    7171         
    … …  
    101101    info = {} 
    102102     
    103      
    104103    def __init__(self): 
    105104        Cinvertor.__init__(self) 
    … …  
    109108        restore the state of invertor for pickle 
    110109        """ 
    111         (self.__dict__,self.alpha, self.d_max, 
    112          self.q_min, self.q_max,  
     110        (self.__dict__, self.alpha, self.d_max, 
     111         self.q_min, self.q_max, 
    113112         self.x, self.y, 
    114113         self.err, self.has_bck, 
    … …  
    120119        """ 
    121120 
    122         state = (self.__dict__,  
     121        state = (self.__dict__, 
    123122                 self.alpha, self.d_max, 
    124123                 self.q_min, self.q_max, 
    125124                 self.x, self.y, 
    126                  self.err, self.has_bck,  
     125                 self.err, self.has_bck, 
    127126                 self.slit_height, self.slit_width, 
    128127                 ) 
    129         return (Invertor,tuple(), state, None, None) 
     128        return (Invertor, tuple(), state, None, None) 
    130129     
    131130    def __setattr__(self, name, value): 
    … …  
    224223         
    225224        invertor = Invertor() 
    226         invertor.chi2    = self.chi2  
    227         invertor.elapsed = self.elapsed  
    228         invertor.nfunc   = self.nfunc  
     225        invertor.chi2    = self.chi2 
     226        invertor.elapsed = self.elapsed 
     227        invertor.nfunc   = self.nfunc 
    229228        invertor.alpha   = self.alpha 
    230229        invertor.d_max   = self.d_max 
    … …  
    237236        invertor.has_bck = self.has_bck 
    238237        invertor.slit_height = self.slit_height 
    239         invertor.slit_width  = self.slit_width 
     238        invertor.slit_width = self.slit_width 
    240239         
    241240        invertor.info = copy.deepcopy(self.info) 
    … …  
    254253        In the following i refers to the ith base function coefficient. 
    255254        The matrix has its entries j in its first Npts rows set to 
    256             A[i][j] = (Fourier transformed base function for point j)  
     255            A[i][j] = (Fourier transformed base function for point j) 
    257256             
    258257        We them choose a number of r-points, n_r, to evaluate the second 
    … …  
    271270        :param nfunc: number of base functions to use. 
    272271        :param nr: number of r points to evaluate the 2nd derivative at for the reg. term. 
    273         :return: c_out, c_cov - the coefficients with covariance matrix  
     272        :return: c_out, c_cov - the coefficients with covariance matrix 
    274273         
    275274        """ 
    … …  
    295294        The minimization function is set to 
    296295        sum_i[ (I_obs(q_i) - I_theo(q_i))/err**2 ] + alpha * reg_term, 
    297         where the reg_term is given by Svergun: it is the integral of  
     296        where the reg_term is given by Svergun: it is the integral of 
    298297        the square of the first derivative 
    299298        of P(r), d(P(r))/dr, integrated over the full range of r. 
    300299         
    301300        :param nfunc: number of base functions to use. 
    302         :param nr: number of r points to evaluate the 2nd derivative at  
     301        :param nr: number of r points to evaluate the 2nd derivative at 
    303302            for the reg. term. 
    304303         
    305         :return: c_out, c_cov - the coefficients with covariance matrix  
    306          
    307         """ 
    308          
    309         #from scipy import optimize 
    310         #import time 
    311          
     304        :return: c_out, c_cov - the coefficients with covariance matrix 
     305         
     306        """ 
    312307        self.nfunc = nfunc 
    313308        # First, check that the current data is valid 
    … …  
    333328         
    334329        if cov_x is None: 
    335             cov_x = numpy.ones([nfunc,nfunc]) 
     330            cov_x = numpy.ones([nfunc, nfunc]) 
    336331            cov_x *= math.fabs(chisqr) 
    337332        return out, cov_x 
    … …  
    343338        with a set of base functions. 
    344339         
    345         This method is provided as a test.  
    346         """ 
    347         #from scipy import optimize 
    348          
     340        This method is provided as a test. 
     341        """ 
    349342        # First, check that the current data is valid 
    350343        if self.is_valid() <= 0: 
    … …  
    354347        p = numpy.ones(nfunc) 
    355348        t_0 = time.time() 
    356         out, cov_x, info, mesg, success = optimize.leastsq(self.pr_residuals, p, 
     349        out, cov_x, _, _, _ = optimize.leastsq(self.pr_residuals, p, 
    357350                                                            full_output=1) 
    358351         
    … …  
    371364     
    372365    def pr_err(self, c, c_cov, r): 
    373         """     
     366        """ 
    374367        Returns the value of P(r) for a given r, and base function 
    375368        coefficients, with error. 
    … …  
    403396        In the following i refers to the ith base function coefficient. 
    404397        The matrix has its entries j in its first Npts rows set to 
    405             A[i][j] = (Fourier transformed base function for point j)  
     398            A[i][j] = (Fourier transformed base function for point j) 
    406399             
    407400        We them choose a number of r-points, n_r, to evaluate the second 
    … …  
    449442            nfunc += 1 
    450443 
    451         a = numpy.zeros([npts+nq, nfunc]) 
    452         b = numpy.zeros(npts+nq) 
     444        a = numpy.zeros([npts + nq, nfunc]) 
     445        b = numpy.zeros(npts + nq) 
    453446        err = numpy.zeros([nfunc, nfunc]) 
    454447         
    … …  
    461454              
    462455        # Perform the inversion (least square fit) 
    463         c, chi2, rank, n = lstsq(a, b) 
     456        c, chi2, _, _ = lstsq(a, b) 
    464457        # Sanity check 
    465458        try: 
    … …  
    469462        self.chi2 = chi2 
    470463                 
    471         inv_cov = numpy.zeros([nfunc,nfunc]) 
     464        inv_cov = numpy.zeros([nfunc, nfunc]) 
    472465        # Get the covariance matrix, defined as inv_cov = a_transposed * a 
    473466        self._get_invcov_matrix(nfunc, nr, a, inv_cov) 
    … …  
    477470                     
    478471        if math.fabs(self.alpha) > 0: 
    479             new_alpha = sum_sig/(sum_reg/self.alpha) 
     472            new_alpha = sum_sig / (sum_reg / self.alpha) 
    480473        else: 
    481474            new_alpha = 0.0 
    … …  
    484477        try: 
    485478            cov = numpy.linalg.pinv(inv_cov) 
    486             err = math.fabs(chi2/float(npts-nfunc)) * cov 
     479            err = math.fabs(chi2 / float(npts - nfunc)) * cov 
    487480        except: 
    488481            # We were not able to estimate the errors 
    … …  
    509502            self.cov = err_0 
    510503             
     504        # Store computation time 
     505        self.elapsed = time.time() - t_0 
     506         
    511507        return self.out, self.cov 
    512508         
    … …  
    516512        number of terms 
    517513         
    518         :param isquit_func: reference to thread function to call to  
     514        :param isquit_func: reference to thread function to call to 
    519515                            check whether the computation needs to 
    520516                            be stopped. 
    … …  
    529525        except: 
    530526            # If we fail, estimate alpha and return the default 
    531             # number of terms  
    532             best_alpha, message, elapsed =self.estimate_alpha(self.nfunc) 
     527            # number of terms 
     528            best_alpha, _, _ = self.estimate_alpha(self.nfunc) 
    533529            return self.nfunc, best_alpha, "Could not estimate number of terms" 
    534530                     
    … …  
    547543        """ 
    548544        #import time 
    549         try:             
     545        try: 
    550546            pr = self.clone() 
    551547             
    … …  
    560556                  
    561557            # Perform inversion to find the largest alpha 
    562             out, cov = pr.invert(nfunc) 
     558            out, _ = pr.invert(nfunc) 
    563559            elapsed = time.time() - starttime 
    564560            initial_alpha = pr.alpha 
    … …  
    567563            # Try the inversion with the estimated alpha 
    568564            pr.alpha = pr.suggested_alpha 
    569             out, cov = pr.invert(nfunc) 
     565            out, _ = pr.invert(nfunc) 
    570566     
    571567            npeaks = pr.get_peaks(out) 
    … …  
    573569            # just return the estimate 
    574570            if npeaks > 1: 
    575                 #message = "Your P(r) is not smooth,  
     571                #message = "Your P(r) is not smooth, 
    576572                #please check your inversion parameters" 
    577573                message = None 
    … …  
    587583                for i in range(10): 
    588584                    pr.alpha = (0.33)**(i+1) * alpha 
    589                     out, cov = pr.invert(nfunc) 
     585                    out, _ = pr.invert(nfunc) 
    590586                     
    591587                    peaks = pr.get_peaks(out) 
    … …  
    599595                # just return that 
    600596                if not found and initial_peaks == 1 and \ 
    601                     initial_alpha<best_alpha: 
     597                    initial_alpha < best_alpha: 
    602598                    best_alpha = initial_alpha 
    603599                     
    … …  
    608604                    message = None 
    609605                elif best_alpha >= 0.5 * pr.suggested_alpha: 
    610                     # best alpha is too big, return a  
     606                    # best alpha is too big, return a 
    611607                    # reasonable value 
    612608                    message  = "The estimated alpha for your system is too " 
    613                     messsage += "large. " 
     609                    message += "large. " 
    614610                    message += "Try increasing your maximum distance." 
    615611                 
    … …  
    620616            return 0, message, elapsed 
    621617     
    622          
    623618    def to_file(self, path, npts=100): 
    624619        """ 
    … …  
    660655        file.close() 
    661656      
    662          
    663657    def from_file(self, path): 
    664658        """ 
    … …  
    676670                fd = open(path, 'r') 
    677671                 
    678                 buff    = fd.read() 
    679                 lines   = buff.split('\n') 
     672                buff = fd.read() 
     673                lines = buff.split('\n') 
    680674                for line in lines: 
    681675                    if line.startswith('#d_max='): 
    … …  
    736730                        self.out[i] = float(toks2[0]) 
    737731                         
    738                         self.cov[i][i] = float(toks2[1])                         
     732                        self.cov[i][i] = float(toks2[1]) 
    739733             
    740734            except: 
    … …  
    742736                raise RuntimeError, msg 
    743737        else: 
    744             msg = "Invertor.from_file: '%s' is not a file" % str(path)  
     738            msg = "Invertor.from_file: '%s' is not a file" % str(path) 
    745739            raise RuntimeError, msg 
    746          
    747           
    748 if __name__ == "__main__": 
    749     o = Invertor() 
    750  
    751      
    752      
    753      
    754      

  • pr_inversion/src/sans/pr/num_term.py

    r3fd8390 r34f3ad0  
    1  
    2 #import unittest 
    31import math 
    42import numpy 
    5 #import sys 
    6 #import string 
    73import copy 
    84from sans.pr.invertor import Invertor 
    95 
    10 class Num_terms():     
     6 
     7class Num_terms(): 
    118    """ 
    129    """ 
    … …  
    6259    def get0_out(self): 
    6360        """ 
    64         """  
     61        """ 
    6562        inver = self.invertor 
    6663        self.osc_list = [] 
    … …  
    7067            if self.isquit_func != None: 
    7168                self.isquit_func() 
    72             best_alpha, message, elapsed = inver.estimate_alpha(k) 
     69            best_alpha, message, _ = inver.estimate_alpha(k) 
    7370            inver.alpha = best_alpha 
    7471            inver.out, inver.cov = inver.lstsq(k) 
    … …  
    8279            self.mess_list.append(message) 
    8380          
    84         #print "osc", self.osc_list 
    85         #print "err", self.err_list 
    86         #print "alpha", self.alpha_list 
    87         #new_ls = [] 
    8881        new_osc1 = [] 
    8982        new_osc2 = [] 
    … …  
    10194            if self.err_list[i] < 0.8 and self.err_list[i] >= 0.7: 
    10295                new_osc3.append(self.osc_list[i]) 
    103                 falg7 = True 
     96                flag7 = True 
    10497                  
    10598        if flag9 == True: 
    … …  
    110103            self.dataset = new_osc3 
    111104          
    112         #print "dataset", self.dataset 
    113105        return self.dataset 
    114106         
    … …  
    148140        try: 
    149141            self.isquit_func = isquit_func 
    150             #self.nterm_max = len(self.invertor.x) 
    151             #self.nterm_max = 32 
    152142            nts = self.compare_err() 
    153             #print "nts", nts 
    154143            div = len(nts) 
    155144            tem = float(div)/2.0 
    … …  
    163152            return self.nterm_min, self.alpha_list[10], self.mess_list[10] 
    164153 
     154 
    165155#For testing 
    166156def load(path): 
    167     #import numpy, math, sys 
    168157    # Read the data from the data file 
    169158    data_x   = numpy.zeros(0) 
    … …  
    187176                        scale = 0.05 * math.sqrt(test_y) 
    188177                        #scale = 0.05/math.sqrt(y) 
    189                         min_err = 0.01*y 
     178                        min_err = 0.01 * y 
    190179                    err = scale * math.sqrt(test_y) + min_err 
    191180                    #err = 0 
    … …  
    204193    x, y, erro = load("test/Cyl_A_D102.txt") 
    205194    i.d_max = 102.0 
    206     i.nfunc = 10     
     195    i.nfunc = 10 
    207196    #i.q_max = 0.4 
    208197    #i.q_min = 0.07 
    209     i.x   = x 
    210     i.y   = y 
    211     i.err = erro    
     198    i.x = x 
     199    i.y = y 
     200    i.err = erro 
    212201    #i.out, i.cov = i.lstsq(10) 
    213202    # Testing estimator 
    214203    est = Num_terms(i) 
    215204    print est.num_terms() 
    216      
    217          
    218          
    219          
    220          
    221          
Note: See TracChangeset for help on using the changeset viewer.