Changeset d84a90c in sasview


Ignore:
Timestamp:
Jun 3, 2010 4:47:46 PM (15 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:
7116b6e0
Parents:
aa36f96
Message:

working on documentation

Location:
pr_inversion
Files:
53 added
3 edited
1 moved

Legend:

Unmodified
Added
Removed

  • pr_inversion/distance_explorer.py

    re83c75a rd84a90c  
     1 
     2################################################################################ 
     3#This software was developed by the University of Tennessee as part of the 
     4#Distributed Data Analysis of Neutron Scattering Experiments (DANSE) 
     5#project funded by the US National Science Foundation.  
     6# 
     7#See the license text in license.txt 
     8# 
     9#copyright 2009, University of Tennessee 
     10################################################################################ 
     11 
    112""" 
    2     Module to explore the P(r) inversion results for a range 
    3     of D_max value. User picks a number of points and a range of 
    4     distances, then get a series of outputs as a function of D_max 
    5     over that range. 
    6      
    7     This software was developed by the University of Tennessee as part of the 
    8     Distributed Data Analysis of Neutron Scattering Experiments (DANSE) 
    9     project funded by the US National Science Foundation.  
     13Module to explore the P(r) inversion results for a range 
     14of D_max value. User picks a number of points and a range of 
     15distances, then get a series of outputs as a function of D_max 
     16over that range. 
     17""" 
    1018 
    11     See the license text in license.txt 
    12  
    13     copyright 2009, University of Tennessee 
    14 """ 
    1519import numpy 
    1620import math 
    … …  
    1923class Results: 
    2024    """ 
    21         Class to hold the inversion output parameters 
    22         as a function of D_max 
     25    Class to hold the inversion output parameters 
     26    as a function of D_max 
    2327    """ 
    2428    def __init__(self): 
    2529        """ 
    26             Initialization. Create empty arrays 
    27             and dictionary of labels. 
     30        Initialization. Create empty arrays 
     31        and dictionary of labels. 
    2832        """ 
    2933        # Array of output for each inversion 
    … …  
    4145class DistExplorer(object): 
    4246    """ 
    43         The explorer class 
     47    The explorer class 
    4448    """ 
    4549     
    4650    def __init__(self, pr_state): 
    4751        """ 
    48             Initialization. 
    49              
    50             @param pr_state: sans.pr.invertor.Invertor object 
     52        Initialization. 
     53         
     54        :param pr_state: sans.pr.invertor.Invertor object 
     55         
    5156        """ 
    5257        self.pr_state  = pr_state 
    … …  
    5762    def __call__(self, dmin=None, dmax=None, npts=10): 
    5863        """ 
    59             Compute the outputs as a function of D_max. 
    60              
    61             @param dmin: minimum value for D_max 
    62             @param dmax: maximum value for D_max 
    63             @param npts: number of points for D_max 
     64        Compute the outputs as a function of D_max. 
     65         
     66        :param dmin: minimum value for D_max 
     67        :param dmax: maximum value for D_max 
     68        :param npts: number of points for D_max 
     69         
    6470        """ 
    6571        # Take care of the defaults if needed 

  • pr_inversion/invertor.py

    r97d69d9 rd84a90c  
    11""" 
    2     Module to perform P(r) inversion. 
    3     The module contains the Invertor class. 
     2Module to perform P(r) inversion. 
     3The module contains the Invertor class. 
    44""" 
    55from sans.pr.core.pr_inversion import Cinvertor 
    … …  
    1111def help(): 
    1212    """ 
    13         Provide general online help text 
    14         Future work: extend this function to allow topic selection 
     13    Provide general online help text 
     14    Future work: extend this function to allow topic selection 
    1515    """ 
    1616    info_txt  = "The inversion approach is based on Moore, J. Appl. Cryst. (1980) 13, 168-175.\n\n" 
    … …  
    3535class Invertor(Cinvertor): 
    3636    """ 
    37         Invertor class to perform P(r) inversion 
    38          
    39         The problem is solved by posing the problem as  Ax = b, 
    40         where x is the set of coefficients we are looking for. 
    41          
    42         Npts is the number of points. 
    43          
    44         In the following i refers to the ith base function coefficient. 
    45         The matrix has its entries j in its first Npts rows set to 
    46             A[j][i] = (Fourier transformed base function for point j)  
    47              
    48         We them choose a number of r-points, n_r, to evaluate the second 
    49         derivative of P(r) at. This is used as our regularization term. 
    50         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]) 
    52              
    53         The vector b has its first Npts entries set to 
    54             b[j] = (I(q) observed for point j) 
    55              
    56         The following n_r entries are set to zero. 
    57          
    58         The result is found by using scipy.linalg.basic.lstsq to invert 
    59         the matrix and find the coefficients x. 
    60          
    61         Methods inherited from Cinvertor: 
    62         - get_peaks(pars): returns the number of P(r) peaks 
    63         - oscillations(pars): returns the oscillation parameters for the output P(r) 
    64         - get_positive(pars): returns the fraction of P(r) that is above zero 
    65         - get_pos_err(pars): returns the fraction of P(r) that is 1-sigma above zero 
     37    Invertor class to perform P(r) inversion 
     38     
     39    The problem is solved by posing the problem as  Ax = b, 
     40    where x is the set of coefficients we are looking for. 
     41     
     42    Npts is the number of points. 
     43     
     44    In the following i refers to the ith base function coefficient. 
     45    The matrix has its entries j in its first Npts rows set to 
     46        A[j][i] = (Fourier transformed base function for point j)  
     47         
     48    We them choose a number of r-points, n_r, to evaluate the second 
     49    derivative of P(r) at. This is used as our regularization term. 
     50    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]) 
     52         
     53    The vector b has its first Npts entries set to 
     54        b[j] = (I(q) observed for point j) 
     55         
     56    The following n_r entries are set to zero. 
     57     
     58    The result is found by using scipy.linalg.basic.lstsq to invert 
     59    the matrix and find the coefficients x. 
     60     
     61    Methods inherited from Cinvertor: 
     62    - get_peaks(pars): returns the number of P(r) peaks 
     63    - oscillations(pars): returns the oscillation parameters for the output P(r) 
     64    - get_positive(pars): returns the fraction of P(r) that is above zero 
     65    - get_pos_err(pars): returns the fraction of P(r) that is 1-sigma above zero 
    6666    """ 
    6767    ## Chisqr of the last computation 
    … …  
    8888    def __setattr__(self, name, value): 
    8989        """ 
    90             Set the value of an attribute. 
    91             Access the parent class methods for 
    92             x, y, err, d_max, q_min, q_max and alpha 
     90        Set the value of an attribute. 
     91        Access the parent class methods for 
     92        x, y, err, d_max, q_min, q_max and alpha 
    9393        """ 
    9494        if   name=='x': 
    … …  
    129129    def __getattr__(self, name): 
    130130        """ 
    131            Return the value of an attribute 
     131        Return the value of an attribute 
    132132        """ 
    133133        import numpy 
    … …  
    174174    def clone(self): 
    175175        """ 
    176             Return a clone of this instance 
     176        Return a clone of this instance 
    177177        """ 
    178178        import copy 
    … …  
    200200    def invert(self, nfunc=10, nr=20): 
    201201        """ 
    202             Perform inversion to P(r) 
    203              
    204             The problem is solved by posing the problem as  Ax = b, 
    205             where x is the set of coefficients we are looking for. 
    206              
    207             Npts is the number of points. 
    208              
    209             In the following i refers to the ith base function coefficient. 
    210             The matrix has its entries j in its first Npts rows set to 
    211                 A[i][j] = (Fourier transformed base function for point j)  
    212                  
    213             We them choose a number of r-points, n_r, to evaluate the second 
    214             derivative of P(r) at. This is used as our regularization term. 
    215             For a vector r of length n_r, the following n_r rows are set to 
    216                 A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, evaluated at r[j]) 
    217                  
    218             The vector b has its first Npts entries set to 
    219                 b[j] = (I(q) observed for point j) 
    220                  
    221             The following n_r entries are set to zero. 
    222              
    223             The result is found by using scipy.linalg.basic.lstsq to invert 
    224             the matrix and find the coefficients x. 
    225              
    226             @param nfunc: number of base functions to use. 
    227             @param nr: number of r points to evaluate the 2nd derivative at for the reg. term. 
    228             @return: c_out, c_cov - the coefficients with covariance matrix  
     202        Perform inversion to P(r) 
     203         
     204        The problem is solved by posing the problem as  Ax = b, 
     205        where x is the set of coefficients we are looking for. 
     206         
     207        Npts is the number of points. 
     208         
     209        In the following i refers to the ith base function coefficient. 
     210        The matrix has its entries j in its first Npts rows set to 
     211            A[i][j] = (Fourier transformed base function for point j)  
     212             
     213        We them choose a number of r-points, n_r, to evaluate the second 
     214        derivative of P(r) at. This is used as our regularization term. 
     215        For a vector r of length n_r, the following n_r rows are set to 
     216            A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, evaluated at r[j]) 
     217             
     218        The vector b has its first Npts entries set to 
     219            b[j] = (I(q) observed for point j) 
     220             
     221        The following n_r entries are set to zero. 
     222         
     223        The result is found by using scipy.linalg.basic.lstsq to invert 
     224        the matrix and find the coefficients x. 
     225         
     226        :param nfunc: number of base functions to use. 
     227        :param nr: number of r points to evaluate the 2nd derivative at for the reg. term. 
     228        :return: c_out, c_cov - the coefficients with covariance matrix  
     229         
    229230        """ 
    230231        # Reset the background value before proceeding 
    … …  
    234235    def iq(self, out, q): 
    235236        """ 
    236             Function to call to evaluate the scattering intensity 
    237             @param args: c-parameters, and q 
    238             @return: I(q) 
     237        Function to call to evaluate the scattering intensity 
     238         
     239        :param args: c-parameters, and q 
     240        :return: I(q) 
     241         
    239242        """ 
    240243        return Cinvertor.iq(self, out, q)+self.background 
    … …  
    242245    def invert_optimize(self, nfunc=10, nr=20): 
    243246        """ 
    244             Slower version of the P(r) inversion that uses scipy.optimize.leastsq. 
    245              
    246             This probably produce more reliable results, but is much slower. 
    247             The minimization function is set to sum_i[ (I_obs(q_i) - I_theo(q_i))/err**2 ] + alpha * reg_term, 
    248             where the reg_term is given by Svergun: it is the integral of the square of the first derivative 
    249             of P(r), d(P(r))/dr, integrated over the full range of r. 
    250              
    251             @param nfunc: number of base functions to use. 
    252             @param nr: number of r points to evaluate the 2nd derivative at for the reg. term. 
    253             @return: c_out, c_cov - the coefficients with covariance matrix  
     247        Slower version of the P(r) inversion that uses scipy.optimize.leastsq. 
     248         
     249        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 
     252        of P(r), d(P(r))/dr, integrated over the full range of r. 
     253         
     254        :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. 
     256         
     257        :return: c_out, c_cov - the coefficients with covariance matrix  
     258         
    254259        """ 
    255260         
    … …  
    281286    def pr_fit(self, nfunc=5): 
    282287        """ 
    283             This is a direct fit to a given P(r). It assumes that the y data 
    284             is set to some P(r) distribution that we are trying to reproduce 
    285             with a set of base functions. 
    286              
    287             This method is provided as a test.  
     288        This is a direct fit to a given P(r). It assumes that the y data 
     289        is set to some P(r) distribution that we are trying to reproduce 
     290        with a set of base functions. 
     291         
     292        This method is provided as a test.  
    288293        """ 
    289294        from scipy import optimize 
    … …  
    312317    def pr_err(self, c, c_cov, r): 
    313318        """     
    314             Returns the value of P(r) for a given r, and base function 
    315             coefficients, with error. 
    316              
    317             @param c: base function coefficients 
    318             @param c_cov: covariance matrice of the base function coefficients 
    319             @param r: r-value to evaluate P(r) at 
    320             @return: P(r) 
     319        Returns the value of P(r) for a given r, and base function 
     320        coefficients, with error. 
     321         
     322        :param c: base function coefficients 
     323        :param c_cov: covariance matrice of the base function coefficients 
     324        :param r: r-value to evaluate P(r) at 
     325         
     326        :return: P(r) 
     327         
    321328        """ 
    322329        return self.get_pr_err(c, c_cov, r) 
    … …  
    324331    def _accept_q(self, q): 
    325332        """ 
    326             Check q-value against user-defined range 
     333        Check q-value against user-defined range 
    327334        """ 
    328335        if not self.q_min==None and q<self.q_min: 
    … …  
    334341    def lstsq(self, nfunc=5, nr=20): 
    335342        """ 
    336             The problem is solved by posing the problem as  Ax = b, 
    337             where x is the set of coefficients we are looking for. 
    338              
    339             Npts is the number of points. 
    340              
    341             In the following i refers to the ith base function coefficient. 
    342             The matrix has its entries j in its first Npts rows set to 
    343                 A[i][j] = (Fourier transformed base function for point j)  
    344                  
    345             We them choose a number of r-points, n_r, to evaluate the second 
    346             derivative of P(r) at. This is used as our regularization term. 
    347             For a vector r of length n_r, the following n_r rows are set to 
    348                 A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, evaluated at r[j]) 
    349                  
    350             The vector b has its first Npts entries set to 
    351                 b[j] = (I(q) observed for point j) 
    352                  
    353             The following n_r entries are set to zero. 
    354              
    355             The result is found by using scipy.linalg.basic.lstsq to invert 
    356             the matrix and find the coefficients x. 
    357              
    358             @param nfunc: number of base functions to use. 
    359             @param nr: number of r points to evaluate the 2nd derivative at for the reg. term. 
    360  
    361             If the result does not allow us to compute the covariance matrix, 
    362             a matrix filled with zeros will be returned. 
     343        The problem is solved by posing the problem as  Ax = b, 
     344        where x is the set of coefficients we are looking for. 
     345         
     346        Npts is the number of points. 
     347         
     348        In the following i refers to the ith base function coefficient. 
     349        The matrix has its entries j in its first Npts rows set to 
     350            A[i][j] = (Fourier transformed base function for point j)  
     351             
     352        We them choose a number of r-points, n_r, to evaluate the second 
     353        derivative of P(r) at. This is used as our regularization term. 
     354        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]) 
     356             
     357        The vector b has its first Npts entries set to 
     358            b[j] = (I(q) observed for point j) 
     359             
     360        The following n_r entries are set to zero. 
     361         
     362        The result is found by using scipy.linalg.basic.lstsq to invert 
     363        the matrix and find the coefficients x. 
     364         
     365        :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. 
     367 
     368        If the result does not allow us to compute the covariance matrix, 
     369        a matrix filled with zeros will be returned. 
    363370 
    364371        """ 
    … …  
    445452    def estimate_numterms(self, isquit_func=None): 
    446453        """ 
    447             Returns a reasonable guess for the 
    448             number of terms 
    449             @param isquit_func: reference to thread function to call to  
    450                                 check whether the computation needs to 
    451                                 be stopped. 
    452              
    453             @return: number of terms, alpha, message 
     454        Returns a reasonable guess for the 
     455        number of terms 
     456         
     457        :param isquit_func: reference to thread function to call to  
     458                            check whether the computation needs to 
     459                            be stopped. 
     460         
     461        :return: number of terms, alpha, message 
     462         
    454463        """ 
    455464        from num_term import Num_terms 
    … …  
    465474    def estimate_alpha(self, nfunc): 
    466475        """ 
    467             Returns a reasonable guess for the 
    468             regularization constant alpha 
    469              
    470             @param nfunc: number of terms to use in the expansion. 
    471             @return: alpha, message, elapsed 
    472              
    473             where alpha is the estimate for alpha, 
    474             message is a message for the user, 
    475             elapsed is the computation time 
     476        Returns a reasonable guess for the 
     477        regularization constant alpha 
     478         
     479        :param nfunc: number of terms to use in the expansion. 
     480         
     481        :return: alpha, message, elapsed 
     482         
     483        where alpha is the estimate for alpha, 
     484        message is a message for the user, 
     485        elapsed is the computation time 
    476486        """ 
    477487        import time 
    … …  
    549559    def to_file(self, path, npts=100): 
    550560        """ 
    551             Save the state to a file that will be readable 
    552             by SliceView. 
    553             @param path: path of the file to write 
    554             @param npts: number of P(r) points to be written 
     561        Save the state to a file that will be readable 
     562        by SliceView. 
     563         
     564        :param path: path of the file to write 
     565        :param npts: number of P(r) points to be written 
     566         
    555567        """ 
    556568        file = open(path, 'w') 
    … …  
    586598    def from_file(self, path): 
    587599        """ 
    588             Load the state of the Invertor from a file, 
    589             to be able to generate P(r) from a set of 
    590             parameters. 
    591             @param path: path of the file to load 
     600        Load the state of the Invertor from a file, 
     601        to be able to generate P(r) from a set of 
     602        parameters. 
     603         
     604        :param path: path of the file to load 
     605         
    592606        """ 
    593607        import os 

  • pr_inversion/num_term.py

    r97d69d9 rd84a90c  
    33 
    44class Num_terms():     
    5      
    6      def __init__(self, invertor): 
    7          self.invertor = invertor 
    8          self.nterm_min = 10 
    9          self.nterm_max = len(self.invertor.x) 
    10          if self.nterm_max>50: 
    11              self.nterm_max=50 
    12          self.isquit_func = None 
    13           
    14          self.osc_list = [] 
    15          self.err_list = [] 
    16          self.alpha_list = [] 
    17          self.mess_list = [] 
    18           
    19          self.dataset = [] 
     5    """ 
     6    """ 
     7    def __init__(self, invertor): 
     8        """ 
     9        """ 
     10        self.invertor = invertor 
     11        self.nterm_min = 10 
     12        self.nterm_max = len(self.invertor.x) 
     13        if self.nterm_max>50: 
     14            self.nterm_max=50 
     15        self.isquit_func = None 
     16          
     17        self.osc_list = [] 
     18        self.err_list = [] 
     19        self.alpha_list = [] 
     20        self.mess_list = [] 
     21          
     22        self.dataset = [] 
    2023      
    21      def is_odd(self, n): 
    22          return bool(n%2) 
    23  
    24      def sort_osc(self): 
    25          import copy 
    26          osc = copy.deepcopy(self.dataset) 
    27          lis = [] 
    28          for i in range(len(osc)): 
    29              osc.sort() 
    30              re = osc.pop(0) 
    31              lis.append(re) 
    32          return lis 
     24    def is_odd(self, n): 
     25        """ 
     26        """ 
     27        return bool(n%2) 
     28 
     29    def sort_osc(self): 
     30        """ 
     31        """ 
     32        import copy 
     33        osc = copy.deepcopy(self.dataset) 
     34        lis = [] 
     35        for i in range(len(osc)): 
     36            osc.sort() 
     37            re = osc.pop(0) 
     38            lis.append(re) 
     39        return lis 
    3340            
    34      def median_osc(self): 
    35          osc = self.sort_osc() 
    36          dv = len(osc) 
    37          med = float(dv) / 2.0 
    38          odd = self.is_odd(dv) 
    39          medi = 0 
    40          for i in range(dv): 
    41              if odd == True: 
    42                  medi = osc[int(med)] 
    43              else: 
    44                  medi = osc[int(med) - 1] 
    45          return medi 
    46  
    47      def get0_out(self): 
     41    def median_osc(self): 
     42        """ 
     43        """ 
     44        osc = self.sort_osc() 
     45        dv = len(osc) 
     46        med = float(dv) / 2.0 
     47        odd = self.is_odd(dv) 
     48        medi = 0 
     49        for i in range(dv): 
     50            if odd == True: 
     51                medi = osc[int(med)] 
     52            else: 
     53                medi = osc[int(med) - 1] 
     54        return medi 
     55 
     56    def get0_out(self): 
     57        """ 
     58        """  
    4859        inver = self.invertor 
    4960        self.osc_list = [] 
    … …  
    96107        return self.dataset 
    97108         
    98      def ls_osc(self): 
     109    def ls_osc(self): 
     110        """ 
     111        """ 
    99112        # Generate data 
    100113        ls_osc = self.get0_out() 
    … …  
    109122        return ls 
    110123 
    111      def compare_err(self): 
     124    def compare_err(self): 
     125        """ 
     126        """ 
    112127        ls = self.ls_osc() 
    113128        #print "ls", ls 
    … …  
    122137        return nt_ls 
    123138 
    124      def num_terms(self, isquit_func=None): 
    125          try: 
    126              self.isquit_func = isquit_func 
    127              #self.nterm_max = len(self.invertor.x) 
    128              #self.nterm_max = 32 
    129              nts = self.compare_err() 
    130              #print "nts", nts 
    131              div = len(nts) 
    132              tem = float(div)/2.0 
    133              odd = self.is_odd(div) 
    134              if odd == True: 
    135                  nt = nts[int(tem)] 
    136              else: 
    137                  nt = nts[int(tem) - 1] 
    138              return nt, self.alpha_list[nt - 10], self.mess_list[nt-10] 
    139          except: 
    140              return self.nterm_min, self.alpha_list[10], self.mess_list[10] 
     139    def num_terms(self, isquit_func=None): 
     140        """ 
     141        """ 
     142        try: 
     143            self.isquit_func = isquit_func 
     144            #self.nterm_max = len(self.invertor.x) 
     145            #self.nterm_max = 32 
     146            nts = self.compare_err() 
     147            #print "nts", nts 
     148            div = len(nts) 
     149            tem = float(div)/2.0 
     150            odd = self.is_odd(div) 
     151            if odd == True: 
     152                nt = nts[int(tem)] 
     153            else: 
     154                nt = nts[int(tem) - 1] 
     155            return nt, self.alpha_list[nt - 10], self.mess_list[nt-10] 
     156        except: 
     157            return self.nterm_min, self.alpha_list[10], self.mess_list[10] 
    141158 
    142159#For testing 
    143160def load(path): 
    144         import numpy, math, sys 
    145         # Read the data from the data file 
    146         data_x   = numpy.zeros(0) 
    147         data_y   = numpy.zeros(0) 
    148         data_err = numpy.zeros(0) 
    149         scale    = None 
    150         min_err  = 0.0 
    151         if not path == None: 
    152             input_f = open(path,'r') 
    153             buff    = input_f.read() 
    154             lines   = buff.split('\n') 
    155             for line in lines: 
    156                 try: 
    157                     toks = line.split() 
    158                     x = float(toks[0]) 
    159                     y = float(toks[1]) 
    160                     if len(toks)>2: 
    161                         err = float(toks[2]) 
    162                     else: 
    163                         if scale==None: 
    164                             scale = 0.05*math.sqrt(y) 
    165                             #scale = 0.05/math.sqrt(y) 
    166                             min_err = 0.01*y 
    167                         err = scale*math.sqrt(y)+min_err 
    168                         #err = 0 
    169                          
    170                     data_x = numpy.append(data_x, x) 
    171                     data_y = numpy.append(data_y, y) 
    172                     data_err = numpy.append(data_err, err) 
    173                 except: 
    174                     pass 
    175                     
    176         return data_x, data_y, data_err 
    177      
     161    import numpy, math, sys 
     162    # Read the data from the data file 
     163    data_x   = numpy.zeros(0) 
     164    data_y   = numpy.zeros(0) 
     165    data_err = numpy.zeros(0) 
     166    scale    = None 
     167    min_err  = 0.0 
     168    if not path == None: 
     169        input_f = open(path,'r') 
     170        buff    = input_f.read() 
     171        lines   = buff.split('\n') 
     172        for line in lines: 
     173            try: 
     174                toks = line.split() 
     175                x = float(toks[0]) 
     176                y = float(toks[1]) 
     177                if len(toks)>2: 
     178                    err = float(toks[2]) 
     179                else: 
     180                    if scale==None: 
     181                        scale = 0.05*math.sqrt(y) 
     182                        #scale = 0.05/math.sqrt(y) 
     183                        min_err = 0.01*y 
     184                    err = scale*math.sqrt(y)+min_err 
     185                    #err = 0 
     186                     
     187                data_x = numpy.append(data_x, x) 
     188                data_y = numpy.append(data_y, y) 
     189                data_err = numpy.append(data_err, err) 
     190            except: 
     191                pass 
     192                
     193    return data_x, data_y, data_err 
     194 
    178195 
    179196if __name__ == "__main__": 
Note: See TracChangeset for help on using the changeset viewer.