Ignore:
Timestamp:
07/29/10 19:13:46 (14 years ago)
Author:
Kana Sugimoto
Message:

New Development: Yes

JIRA Issue: No (test merging alma branch)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s):

Description:


Location:
branches/mergetest
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/mergetest

  • branches/mergetest/python/asapfitter.py

    r1739 r1779  
    11import _asap
    22from asap import rcParams
    3 from asap import print_log_dec
     3from asap import print_log, print_log_dec
    44from asap import _n_bools
    55from asap import mask_and
     6from asap import asaplog
    67
    78class fitter:
    … …  
    5960            msg = "Please give a correct scan"
    6061            if rcParams['verbose']:
    61                 print msg
     62                #print msg
     63                asaplog.push(msg)
     64                print_log('ERROR')
    6265                return
    6366            else:
    … …  
    7982            lpoly:   use polynomial of the order given with linear least squares fit
    8083            gauss:   fit the number of gaussian specified
     84            lorentz: fit the number of lorentzian specified
    8185        Example:
    82             fitter.set_function(gauss=2) # will fit two gaussians
    8386            fitter.set_function(poly=3)  # will fit a 3rd order polynomial via nonlinear method
    8487            fitter.set_function(lpoly=3)  # will fit a 3rd order polynomial via linear method
     88            fitter.set_function(gauss=2) # will fit two gaussians
     89            fitter.set_function(lorentz=2) # will fit two lorentzians
    8590        """
    8691        #default poly order 0
    … …  
    102107            self.components = [ 3 for i in range(n) ]
    103108            self.uselinear = False
     109        elif kwargs.has_key('lorentz'):
     110            n = kwargs.get('lorentz')
     111            self.fitfunc = 'lorentz'
     112            self.fitfuncs = [ 'lorentz' for i in range(n) ]
     113            self.components = [ 3 for i in range(n) ]
     114            self.uselinear = False
    104115        else:
    105116            msg = "Invalid function type."
    106117            if rcParams['verbose']:
    107                 print msg
     118                #print msg
     119                asaplog.push(msg)
     120                print_log('ERROR')
    108121                return
    109122            else:
    … …  
    114127        return
    115128
    116     @print_log_dec
     129    #@print_log_dec
    117130    def fit(self, row=0, estimate=False):
    118131        """
    … …  
    135148            msg = "Fitter not yet initialised. Please set data & fit function"
    136149            if rcParams['verbose']:
    137                 print msg
     150                #print msg
     151                asaplog.push(msg)
     152                print_log('ERROR')
    138153                return
    139154            else:
    … …  
    145160                self.y = self.data._getspectrum(row)
    146161                self.mask = mask_and(self.mask, self.data._getmask(row))
    147                 from asap import asaplog
    148162                asaplog.push("Fitting:")
    149163                i = row
    … …  
    155169                asaplog.push(out,False)
    156170        self.fitter.setdata(self.x, self.y, self.mask)
    157         if self.fitfunc == 'gauss':
     171        if self.fitfunc == 'gauss' or self.fitfunc == 'lorentz':
    158172            ps = self.fitter.getparameters()
    159173            if len(ps) == 0 or estimate:
    … …  
    171185        except RuntimeError, msg:
    172186            if rcParams['verbose']:
    173                 print msg
     187                #print msg
     188                print_log()
     189                asaplog.push(str(msg))
     190                print_log('ERROR')
    174191            else:
    175192                raise
    176193        self._fittedrow = row
    177194        self.fitted = True
     195        print_log()
    178196        return
    179197
    … …  
    204222                self.data._addfit(fit,self._fittedrow)
    205223
    206     @print_log_dec
     224    #@print_log_dec
    207225    def set_parameters(self,*args,**kwargs):
    208226        """
    … …  
    228246            msg = "Please specify a fitting function first."
    229247            if rcParams['verbose']:
    230                 print msg
     248                #print msg
     249                asaplog.push(msg)
     250                print_log('ERROR')
    231251                return
    232252            else:
    233253                raise RuntimeError(msg)
    234         if self.fitfunc == "gauss" and component is not None:
     254        if (self.fitfunc == "gauss" or self.fitfunc == 'lorentz') and component is not None:
    235255            if not self.fitted and sum(self.fitter.getparameters()) == 0:
    236256                pars = _n_bools(len(self.components)*3, False)
    … …  
    247267        if fixed is not None:
    248268            self.fitter.setfixedparameters(fixed)
     269        print_log()
    249270        return
    250271
    … …  
    269290            msg = "Function only operates on Gaussian components."
    270291            if rcParams['verbose']:
    271                 print msg
     292                #print msg
     293                asaplog.push(msg)
     294                print_log('ERROR')
    272295                return
    273296            else:
    … …  
    280303            msg = "Please select a valid  component."
    281304            if rcParams['verbose']:
    282                 print msg
     305                #print msg
     306                asaplog.push(msg)
     307                print_log('ERROR')
    283308                return
    284309            else:
    285310                raise ValueError(msg)
    286311
     312    def set_lorentz_parameters(self, peak, centre, fwhm,
     313                             peakfixed=0, centrefixed=0,
     314                             fwhmfixed=0,
     315                             component=0):
     316        """
     317        Set the Parameters of a 'Lorentzian' component, set with set_function.
     318        Parameters:
     319            peak, centre, fwhm:  The gaussian parameters
     320            peakfixed,
     321            centrefixed,
     322            fwhmfixed:           Optional parameters to indicate if
     323                                 the paramters should be held fixed during
     324                                 the fitting process. The default is to keep
     325                                 all parameters flexible.
     326            component:           The number of the component (Default is the
     327                                 component 0)
     328        """
     329        if self.fitfunc != "lorentz":
     330            msg = "Function only operates on Lorentzian components."
     331            if rcParams['verbose']:
     332                #print msg
     333                asaplog.push(msg)
     334                print_log('ERROR')
     335                return
     336            else:
     337                raise ValueError(msg)
     338        if 0 <= component < len(self.components):
     339            d = {'params':[peak, centre, fwhm],
     340                 'fixed':[peakfixed, centrefixed, fwhmfixed]}
     341            self.set_parameters(d, component)
     342        else:
     343            msg = "Please select a valid  component."
     344            if rcParams['verbose']:
     345                #print msg
     346                asaplog.push(msg)
     347                print_log('ERROR')
     348                return
     349            else:
     350                raise ValueError(msg)
     351
    287352    def get_area(self, component=None):
    288353        """
    289         Return the area under the fitted gaussian component.
    290         Parameters:
    291               component:   the gaussian component selection,
     354        Return the area under the fitted gaussian/lorentzian component.
     355        Parameters:
     356              component:   the gaussian/lorentzian component selection,
    292357                           default (None) is the sum of all components
    293358        Note:
    294               This will only work for gaussian fits.
     359              This will only work for gaussian/lorentzian fits.
    295360        """
    296361        if not self.fitted: return
    297         if self.fitfunc == "gauss":
     362        if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
    298363            pars = list(self.fitter.getparameters())
    299364            from math import log,pi,sqrt
    300             fac = sqrt(pi/log(16.0))
     365            if self.fitfunc == "gauss":
     366                fac = sqrt(pi/log(16.0))
     367            elif self.fitfunc == "lorentz":
     368                fac = pi/2.0
    301369            areas = []
    302370            for i in range(len(self.components)):
    … …  
    321389            msg = "Not yet fitted."
    322390            if rcParams['verbose']:
    323                 print msg
     391                #print msg
     392                asaplog.push(msg)
     393                print_log('ERROR')
    324394                return
    325395            else:
    … …  
    328398        cerrs = errs
    329399        if component is not None:
    330             if self.fitfunc == "gauss":
     400            if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
    331401                i = 3*component
    332402                if i < len(errs):
    … …  
    344414            msg = "Not yet fitted."
    345415            if rcParams['verbose']:
    346                 print msg
     416                #print msg
     417                asaplog.push(msg)
     418                print_log('ERROR')
    347419                return
    348420            else:
    … …  
    353425        area = []
    354426        if component is not None:
    355             if self.fitfunc == "gauss":
     427            if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
    356428                i = 3*component
    357429                cpars = pars[i:i+3]
    … …  
    368440            cfixed = fixed
    369441            cerrs = errs
    370             if self.fitfunc == "gauss":
     442            if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
    371443                for c in range(len(self.components)):
    372444                  a = self.get_area(c)
    … …  
    374446        fpars = self._format_pars(cpars, cfixed, errors and cerrs, area)
    375447        if rcParams['verbose']:
    376             print fpars
     448            #print fpars
     449            asaplog.push(fpars)
     450            print_log()
    377451        return {'params':cpars, 'fixed':cfixed, 'formatted': fpars,
    378452                'errors':cerrs}
    … …  
    391465                c+=1
    392466            out = out[:-1]  # remove trailing ','
    393         elif self.fitfunc == 'gauss':
     467        elif self.fitfunc == 'gauss' or self.fitfunc == 'lorentz':
    394468            i = 0
    395469            c = 0
    … …  
    415489        fixed = self.fitter.getfixedparameters()
    416490        if rcParams['verbose']:
    417             print self._format_pars(pars,fixed,None)
     491            #print self._format_pars(pars,fixed,None)
     492            asaplog.push(self._format_pars(pars,fixed,None))
     493            print_log()
    418494        return pars
    419495
    … …  
    425501            msg = "Not yet fitted."
    426502            if rcParams['verbose']:
    427                 print msg
     503                #print msg
     504                asaplog.push(msg)
     505                print_log('ERROR')
    428506                return
    429507            else:
    … …  
    438516            msg = "Not yet fitted."
    439517            if rcParams['verbose']:
    440                 print msg
     518                #print msg
     519                asaplog.push(msg)
     520                print_log('ERROR')
    441521                return
    442522            else:
    … …  
    444524        ch2 = self.fitter.getchi2()
    445525        if rcParams['verbose']:
    446             print 'Chi^2 = %3.3f' % (ch2)
     526            #print 'Chi^2 = %3.3f' % (ch2)
     527            asaplog.push( 'Chi^2 = %3.3f' % (ch2) )
     528            print_log()
    447529        return ch2
    448530
    … …  
    454536            msg = "Not yet fitted."
    455537            if rcParams['verbose']:
    456                 print msg
     538                #print msg
     539                asaplog.push(msg)
     540                print_log('ERROR')
    457541                return
    458542            else:
    … …  
    460544        return self.fitter.getfit()
    461545
    462     @print_log_dec
     546    #@print_log_dec
    463547    def commit(self):
    464548        """
    … …  
    468552            msg = "Not yet fitted."
    469553            if rcParams['verbose']:
    470                 print msg
     554                #print msg
     555                asaplog.push(msg)
     556                print_log('ERROR')
    471557                return
    472558            else:
    … …  
    476562            msg = "Not a scantable"
    477563            if rcParams['verbose']:
    478                 print msg
     564                #print msg
     565                asaplog.push(msg)
     566                print_log('ERROR')
    479567                return
    480568            else:
    … …  
    482570        scan = self.data.copy()
    483571        scan._setspectrum(self.fitter.getresidual())
     572        print_log()
    484573        return scan
    485574
    486     @print_log_dec
     575    #@print_log_dec
    487576    def plot(self, residual=False, components=None, plotparms=False,
    488577             filename=None):
    … …  
    525614
    526615        colours = ["#777777","#dddddd","red","orange","purple","green","magenta", "cyan"]
     616        nomask=True
     617        for i in range(len(m)):
     618            nomask = nomask and m[i]
     619        label0='Masked Region'
     620        label1='Spectrum'
     621        if ( nomask ):
     622            label0=label1
     623        else:
     624            y = ma.masked_array( self.y, mask = m )
     625            self._p.palette(1,colours)
     626            self._p.set_line( label = label1 )
     627            self._p.plot( self.x, y )
    527628        self._p.palette(0,colours)
    528         self._p.set_line(label='Spectrum')
     629        self._p.set_line(label=label0)
    529630        y = ma.masked_array(self.y,mask=logical_not(m))
    530631        self._p.plot(self.x, y)
    531632        if residual:
    532             self._p.palette(1)
     633            self._p.palette(7)
    533634            self._p.set_line(label='Residual')
    534635            y = ma.masked_array(self.get_residual(),
    … …  
    571672        if (not rcParams['plotter.gui']):
    572673            self._p.save(filename)
    573 
    574     @print_log_dec
     674        print_log()
     675
     676    #@print_log_dec
    575677    def auto_fit(self, insitu=None, plot=False):
    576678        """
    … …  
    583685            msg = "Data is not a scantable"
    584686            if rcParams['verbose']:
    585                 print msg
     687                #print msg
     688                asaplog.push(msg)
     689                print_log('ERROR')
    586690                return
    587691            else:
    … …  
    593697            scan = self.data
    594698        rows = xrange(scan.nrow())
    595         from asap import asaplog
     699        # Save parameters of baseline fits as a class attribute.
     700        # NOTICE: This does not reflect changes in scantable!
     701        if len(rows) > 0: self.blpars=[]
    596702        asaplog.push("Fitting:")
    597703        for r in rows:
    … …  
    608714            self.fit()
    609715            x = self.get_parameters()
     716            fpar = self.get_parameters()
    610717            if plot:
    611718                self.plot(residual=True)
    612719                x = raw_input("Accept fit ([y]/n): ")
    613720                if x.upper() == 'N':
     721                    self.blpars.append(None)
    614722                    continue
    615723            scan._setspectrum(self.fitter.getresidual(), r)
     724            self.blpars.append(fpar)
    616725        if plot:
    617726            self._p.unmap()
    618727            self._p = None
     728        print_log()
    619729        return scan
Note: See TracChangeset for help on using the changeset viewer.