Python SpecfitFuns.subac方法代码示例

# 需要导入模块: from PyMca5.PyMcaMath.fitting import SpecfitFuns [as 别名]
# 或者: from PyMca5.PyMcaMath.fitting.SpecfitFuns import subac [as 别名]
    def guess_fwhm(self, **kw):
        if 'y' in kw:
            y=kw['y']
        else:
            return self.config['FwhmPoints']
        if 'x' in kw:
            x=kw['x']
        else:
            x=numpy.arange(len(y))*1.0

        #set at least a default value for the fwhm
        fwhm=4

        zz=SpecfitFuns.subac(y,1.000,1000)
        yfit=y-zz

        #now I should do some sort of peak search ...
        maximum = max(yfit)
        idx = numpy.nonzero(yfit == maximum)[0]
        pos = numpy.take(x,idx)[-1]
        posindex = idx[-1]
        height = yfit[posindex]
        imin = posindex
        while ((yfit[imin] > 0.5*height) & (imin >0)):
            imin=imin - 1
        imax=posindex
        while ((yfit[imax] > 0.5*height) & (imax <(len(yfit)-1))):
            imax=imax + 1
        fwhm=max(imax-imin-1,fwhm)

        return fwhm

# 需要导入模块: from PyMca5.PyMcaMath.fitting import SpecfitFuns [as 别名]
# 或者: from PyMca5.PyMcaMath.fitting.SpecfitFuns import subac [as 别名]
 def bkg_internal(self,pars,x):
    """
    Internal Background
    """
    #fast
    #return self.zz
    #slow: recalculate the background as function of the parameters
    #yy=SpecfitFuns.subac(self.ydata*self.fitconfig['Yscaling'],
    #                     pars[0],pars[1])
    yy=SpecfitFuns.subac(self.ydata*1.0,
                         pars[0],pars[1])
    nrx=shape(x)[0]
    nry=shape(yy)[0]
    if nrx == nry:
         return SpecfitFuns.subac(yy,pars[0],pars[1])
    else:
         return SpecfitFuns.subac(numpy.take(yy,numpy.arange(0,nry,2)),pars[0],pars[1])

# 需要导入模块: from PyMca5.PyMcaMath.fitting import SpecfitFuns [as 别名]
# 或者: from PyMca5.PyMcaMath.fitting.SpecfitFuns import subac [as 别名]
    def estimate_linear(self, xx, yy, zzz, xscaling=1.0, yscaling=None):
        # compute strip bg and use it to estimate the linear bg parameters
        zz = SpecfitFuns.subac(yy, 1.000, 10000)
        n = float(len(zz))
        Sy = numpy.sum(zz)
        Sx = float(numpy.sum(xx))
        Sxx = float(numpy.sum(xx * xx))
        Sxy = float(numpy.sum(xx * zz))

        deno = n * Sxx - (Sx * Sx)
        if deno != 0:
            bg = (Sxx * Sy - Sx * Sxy) / deno
            slope = (n * Sxy - Sx * Sy) / deno
        else:
            bg = 0.0
            slope = 0.0
        estimated_par = [bg, slope]
        # code = 0: FREE
        constraints = [[0, 0], [0, 0], [0, 0]]
        return estimated_par, constraints

# 需要导入模块: from PyMca5.PyMcaMath.fitting import SpecfitFuns [as 别名]
# 或者: from PyMca5.PyMcaMath.fitting.SpecfitFuns import subac [as 别名]
    def estimate_gauss(self,xx,yy,zzz,xscaling=1.0,yscaling=None):
       if yscaling == None:
            try:
                yscaling=self.config['Yscaling']
            except:
                yscaling=1.0
       if yscaling == 0:
            yscaling=1.0
       fittedpar=[]
       zz=SpecfitFuns.subac(yy,1.000,10000)

       npoints = len(zz)
       if self.config['AutoFwhm']:
            search_fwhm=self.guess_fwhm(x=xx,y=yy)
       else:
            search_fwhm=int(float(self.config['FwhmPoints']))
       search_sens=float(self.config['Sensitivity'])
       search_mca=int(float(self.config['McaMode']))

       if search_fwhm < 3:
            search_fwhm = 3
            self.config['FwhmPoints']=3

       if search_sens < 1:
            search_sens = 1
            self.config['Sensitivity']=1

       if npoints > 1.5*search_fwhm:
            peaks=self.seek(yy,fwhm=search_fwhm,
                               sensitivity=search_sens,
                               yscaling=yscaling,
                               mca=search_mca)
            #print "estimate peaks = ",peaks
            #peaks=self.seek(yy-zz,fwhm=search_fwhm,
            #                   sensitivity=search_sens,
            #                   yscaling=yscaling,
            #                   mca=search_mca)
       else:
            peaks = []
       if not len(peaks):
            mca = int(float(self.config.get('McaMode', 0)))
            forcePeak =  int(float(self.config.get('ForcePeakPresence', 0)))
            if (not mca) and forcePeak:
                delta = yy - zz
                peaks  = [int(numpy.nonzero(delta == delta.max())[0])]

       #print "peaks = ",peaks
       #print "peaks subac = ",self.seek(yy-zz,fwhm=search_fwhm,
       #                        sensitivity=search_sens,
       #                        yscaling=yscaling,
       #                        mca=search_mca)
       largest_index = 0
       if len(peaks) > 0:
            j = 0
            for i in peaks:
                if j == 0:
                    sig=5*abs(xx[npoints-1]-xx[0])/npoints
                    peakpos = xx[int(i)]
                    if abs(peakpos) < 1.0e-16:
                        peakpos = 0.0
                    param = numpy.array([yy[int(i)] - zz[int(i)], peakpos ,sig])
                    largest = param
                else:
                    param2 = numpy.array([yy[int(i)] - zz [int(i)], xx[int(i)] ,sig])
                    param = numpy.concatenate((param,param2))
                    if (param2[0] > largest[0]):
                        largest = param2
                        largest_index = j
                j = j + 1
            xw = numpy.resize(xx,(npoints,1))
            if (0):
                sy = numpy.sqrt(abs(yy))
                yw = numpy.resize(yy-zz,(npoints,1))
                sy = numpy.resize(sy,(npoints,1))
                datawork = numpy.concatenate((xw,yw,sy),1)
            else:
                yw = numpy.resize(yy-zz,(npoints,1))
                datawork = numpy.concatenate((xw,yw),1)
            cons = numpy.zeros((3,len(param)),numpy.float)
            cons [0] [0:len(param):3] = CPOSITIVE
            #force peaks to stay around their position
            if (1):
                cons [0] [1:len(param):3] = CQUOTED
                #This does not work!!!!
                #FWHM should be given in terms of X not of points!
                #cons [1] [1:len(param):3] = param[1:len(param):3]-0.5*search_fwhm
                #cons [2] [1:len(param):3] = param[1:len(param):3]+0.5*search_fwhm
                if len(xw) > search_fwhm:
                    fwhmx = numpy.fabs(xw[int(search_fwhm)]-xw[0])
                    cons [1] [1:len(param):3] = param[1:len(param):3]-0.5*fwhmx
                    cons [2] [1:len(param):3] = param[1:len(param):3]+0.5*fwhmx
                else:
                    cons [1] [1:len(param):3] = numpy.ones(shape(param[1:len(param):3]),numpy.float)*min(xw)
                    cons [2] [1:len(param):3] = numpy.ones(shape(param[1:len(param):3]),numpy.float)*max(xw)

            if 0:
                cons [0] [2:len(param):3] = CFACTOR
                cons [1] [2:len(param):3] = 2
                cons [2] [2:len(param):3] = 1.0
                cons [0] [2] = CPOSITIVE
#.........这里部分代码省略.........

# 需要导入模块: from PyMca5.PyMcaMath.fitting import SpecfitFuns [as 别名]
# 或者: from PyMca5.PyMcaMath.fitting.SpecfitFuns import subac [as 别名]
    def estimate_periodic_gauss(self,xx,yy,zzz,xscaling=1.0,yscaling=None):
        if yscaling == None:
            try:
                yscaling=self.config['Yscaling']
            except:
                yscaling=1.0
        if yscaling == 0:
            yscaling=1.0
        fittedpar=[]
        zz=SpecfitFuns.subac(yy,1.000,10000)

        npoints = len(zz)
        if self.config['AutoFwhm']:
            search_fwhm=self.guess_fwhm(x=xx,y=yy)
        else:
            search_fwhm=int(float(self.config['FwhmPoints']))
        search_sens=float(self.config['Sensitivity'])
        search_mca=int(float(self.config['McaMode']))

        if search_fwhm < 3:
            search_fwhm = 3
            self.config['FwhmPoints']=3

        if search_sens < 1:
            search_sens = 1
            self.config['Sensitivity']=1

        if npoints > 1.5*search_fwhm:
            peaks=self.seek(yy,fwhm=search_fwhm,
                               sensitivity=search_sens,
                               yscaling=yscaling,
                               mca=search_mca)
        else:
             peaks = []
        npeaks = len(peaks)
        if not npeaks:
            fittedpar = []
            cons = numpy.zeros((3,len(fittedpar)), numpy.float)
            return fittedpar, cons

        fittedpar = [0.0, 0.0, 0.0, 0.0, 0.0]

        #The number of peaks
        fittedpar[0] = npeaks

        #The separation between peaks in x units
        delta  = 0.0
        height = 0.0
        for i in range(npeaks):
            height += yy[int(peaks[i])] - zz [int(peaks[i])]
            if i != ((npeaks)-1):
                delta += (xx[int(peaks[i+1])] - xx[int(peaks[i])])

        #delta between peaks
        if npeaks > 1:
            fittedpar[1] = delta/(npeaks-1)

        #starting height
        fittedpar[2] = height/npeaks

        #position of the first peak
        fittedpar[3] = xx[int(peaks[0])]

        #Estimate the fwhm
        fittedpar[4] = search_fwhm

        #setup constraints
        cons = numpy.zeros((3, 5), numpy.float)
        cons [0] [0] = CFIXED     #the number of gaussians
        if npeaks == 1:
            cons[0][1] = CFIXED  #the delta between peaks
        else:
            cons[0][1] = CFREE   #the delta between peaks
        j = 2
        #Setup height area constrains
        if self.config['NoConstrainsFlag'] == 0:
            if self.config['HeightAreaFlag']:
                #POSITIVE = 1
                cons[0] [j] = 1
                cons[1] [j] = 0
                cons[2] [j] = 0
        j=j+1

        #Setup position constrains
        if self.config['NoConstrainsFlag'] == 0:
            if self.config['PositionFlag']:
                        #QUOTED = 2
                        cons[0][j]=2
                        cons[1][j]=min(xx)
                        cons[2][j]=max(xx)
        j=j+1

        #Setup positive FWHM constrains
        if self.config['NoConstrainsFlag'] == 0:
            if self.config['PosFwhmFlag']:
                        #POSITIVE=1
                        cons[0][j]=1
                        cons[1][j]=0
                        cons[2][j]=0
        j=j+1
#.........这里部分代码省略.........

# 需要导入模块: from PyMca5.PyMcaMath.fitting import SpecfitFuns [as 别名]
# 或者: from PyMca5.PyMcaMath.fitting.SpecfitFuns import subac [as 别名]
    def update(self):
        if self._y is None:
            return

        pars = self.getParameters()

        #smoothed data
        y = numpy.ravel(numpy.array(self._y)).astype(numpy.float)
        ysmooth = SpecfitFuns.SavitskyGolay(y, pars['stripfilterwidth'])        
        f=[0.25,0.5,0.25]
        ysmooth[1:-1] = numpy.convolve(ysmooth,f,mode=0)
        ysmooth[0] = 0.5 *(ysmooth[0] + ysmooth[1])
        ysmooth[-1] = 0.5 * (ysmooth[-1] + ysmooth[-2])

        #loop for anchors
        x = self._x
        niter = pars['stripiterations']
        anchorslist = []
        if pars['stripanchorsflag']:
            if pars['stripanchorslist'] is not None:
                ravelled = x
                for channel in pars['stripanchorslist']:
                    if channel <= ravelled[0]:continue
                    index = numpy.nonzero(ravelled >= channel)[0]
                    if len(index):
                        index = min(index)
                        if index > 0:
                            anchorslist.append(index)
        if niter > 1000:
            stripBackground = SpecfitFuns.subac(ysmooth,
                                            pars['stripconstant'],
                                            niter,
                                            pars['stripwidth'],
                                            anchorslist)
            #final smoothing
            stripBackground = SpecfitFuns.subac(stripBackground,
                                  pars['stripconstant'],
                                  500,1,
                                  anchorslist)
        elif niter > 0:
            stripBackground = SpecfitFuns.subac(ysmooth,
                                            pars['stripconstant'],
                                            niter,
                                            pars['stripwidth'],
                                            anchorslist)
        else:
            stripBackground = 0.0 * ysmooth + ysmooth.min()

        if len(anchorslist) == 0:
            anchorslist = [0, len(ysmooth)-1]
        anchorslist.sort()
        snipBackground = 0.0 * ysmooth
        lastAnchor = 0
        width = pars['snipwidth']
        for anchor in anchorslist:
            if (anchor > lastAnchor) and (anchor < len(ysmooth)):
                snipBackground[lastAnchor:anchor] =\
                            SpecfitFuns.snip1d(ysmooth[lastAnchor:anchor], width, 0)
                lastAnchor = anchor
        if lastAnchor < len(ysmooth):
            snipBackground[lastAnchor:] =\
                            SpecfitFuns.snip1d(ysmooth[lastAnchor:], width, 0)

        self.graphWidget.addCurve(x, y, \
                                  legend='Input Data',\
                                  replace=True,
                                  replot=False)
        self.graphWidget.addCurve(x, stripBackground,\
                                  legend='Strip Background',\
                                  replot=False)
        self.graphWidget.addCurve(x, snipBackground,\
                                  legend='SNIP Background',
                                  replot=True)