Python fitting.SpecfitFuns类代码示例

    def calculateShiftsFFT(self, stack, reference, offsets=None, widths=None, crop=False):
        if DEBUG:
            print("Offsets = ", offsets)
            print("Widths = ", widths)
        data = stack.data
        if offsets is None:
            offsets = [0.0, 0.0]
        if widths is None:
            widths = [reference.shape[0], reference.shape[1]]
        fft2Function = numpy.fft.fft2
        if 1:
            DTYPE = numpy.float32
        else:
            DTYPE = numpy.float64
        image2 = numpy.zeros((widths[0], widths[1]), dtype=DTYPE)
        shape = image2.shape

        USE_APODIZATION_WINDOW = False
        apo = [10, 10]
        if USE_APODIZATION_WINDOW:
            # use apodization window
            window = numpy.outer(SpecfitFuns.slit([0.5, shape[0] * 0.5, shape[0] - 4 * apo[0], apo[0]],
                                          numpy.arange(float(shape[0]))),
                                 SpecfitFuns.slit([0.5, shape[1] * 0.5, shape[1] - 4 * apo[1], apo[1]],
                                          numpy.arange(float(shape[1])))).astype(DTYPE)
        else:
            window = numpy.zeros((shape[0], shape[1]), dtype=DTYPE)
            window[apo[0]:shape[0] - apo[0], apo[1]:shape[1] - apo[1]] = 1
        image2[:,:] = window * reference[offsets[0]:offsets[0]+widths[0],
                                         offsets[1]:offsets[1]+widths[1]]
        image2fft2 = fft2Function(image2)
        mcaIndex = stack.info.get('McaIndex')
        shifts = numpy.zeros((data.shape[mcaIndex], 2), numpy.float)
        image1 = numpy.zeros(image2.shape, dtype=DTYPE)
        total = float(data.shape[mcaIndex])
        if mcaIndex == 0:
            for i in range(data.shape[mcaIndex]):
                image1[:,:] = window * data[i][offsets[0]:offsets[0]+widths[0],
                                               offsets[1]:offsets[1]+widths[1]]

                image1fft2 = fft2Function(image1)
                shifts[i] = ImageRegistration.measure_offset_from_ffts(image2fft2,
                                                                       image1fft2)
                if DEBUG:
                    print("Index = %d shift = %.4f, %.4f" % (i, shifts[i][0], shifts[i][1]))
                self._progress = (100 * i) / total
        elif mcaIndex in [2, -1]:
            for i in range(data.shape[mcaIndex]):
                image1[:,:] = window * data[:,:,i][offsets[0]:offsets[0]+widths[0],
                                               offsets[1]:offsets[1]+widths[1]]

                image1fft2 = fft2Function(image1)
                shifts[i] = ImageRegistration.measure_offset_from_ffts(image2fft2,
                                                                       image1fft2)
                if DEBUG:
                    print("Index = %d shift = %.4f, %.4f" % (i, shifts[i][0], shifts[i][1]))
                self._progress = (100 * i) / total
        else:
            raise IndexError("Only stacks of images or spectra supported. 1D index should be 0 or 2")
        return shifts

def estimateXANESEdge(spectrum, energy=None, full=False):
    if energy is None:
        x = numpy.arange(len(spectrum)).astype(numpy.float)
    else:
        x = numpy.array(energy, dtype=numpy.float, copy=True)
    y = numpy.array(spectrum, dtype=numpy.float, copy=True)

    # make sure data are sorted
    idx = energy.argsort(kind='mergesort')
    x = numpy.take(energy, idx)
    y = numpy.take(spectrum, idx)

    # make sure data are strictly increasing
    delta = x[1:] - x[:-1]
    dmin = delta.min()
    dmax = delta.max()
    if delta.min() <= 1.0e-10:
        # force data are strictly increasing
        # although we do not consider last point
        idx = numpy.nonzero(delta>0)[0]
        x = numpy.take(x, idx)
        y = numpy.take(y, idx)
        delta = None

    sortedX = x
    sortedY = y

    # use a regularly spaced spectrum
    if dmax != dmin:
        # choose the number of points or deduce it from
        # the input data length?
        nchannels = 2 * len(spectrum)
        xi = numpy.linspace(x[1], x[-2], nchannels).reshape(-1, 1)
        x.shape = -1
        y.shape = -1
        y = SpecfitFuns.interpol([x], y, xi, y.min())
        x = xi
        x.shape = -1
        y.shape = -1

    # take the first derivative
    npoints = 7
    xPrime = x[npoints:-npoints]
    yPrime = SGModule.getSavitzkyGolay(y, npoints=npoints, degree=2, order=1)

    # get the index at maximum value
    iMax = numpy.argmax(yPrime)

    # get the center of mass
    w = 2 * npoints
    selection = yPrime[iMax-w:iMax+w+1]
    edge = (selection * xPrime[iMax-w:iMax+w+1]).sum(dtype=numpy.float)/\
           selection.sum(dtype=numpy.float)

    if full:
        # return intermediate information
        return edge, sortedX, sortedY, xPrime, yPrime
    else:
        # return the corresponding x value
        return edge

    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

 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])

def estimateXANESEdge(spectrum, energy=None, npoints=5, full=False,
                      sanitize=True):
    if sanitize:
        if energy is None:
            x = numpy.arange(len(spectrum)).astype(numpy.float)
        else:
            x = numpy.array(energy, dtype=numpy.float, copy=False)
        y = numpy.array(spectrum, dtype=numpy.float, copy=False)
        # make sure data are sorted
        idx = energy.argsort(kind='mergesort')
        x = numpy.take(energy, idx)
        y = numpy.take(spectrum, idx)

        # make sure data are strictly increasing
        delta = x[1:] - x[:-1]
        dmin = delta.min()
        dmax = delta.max()
        if delta.min() <= 1.0e-10:
            # force data are strictly increasing
            # although we do not consider last point
            idx = numpy.nonzero(delta>0)[0]
            x = numpy.take(x, idx)
            y = numpy.take(y, idx)
            delta = None

        # use a regularly spaced spectrum
        if dmax != dmin:
            # choose the number of points or deduce it from
            # the input data length?
            nchannels = 10 * x.size
            xi = numpy.linspace(x[1], x[-2], nchannels).reshape(-1, 1)
            x.shape = -1
            y.shape = -1
            y = SpecfitFuns.interpol([x], y, xi, y.min())
            x = xi
    else:
        # take views
        x = energy[:]
        y = spectrum[:]

    x.shape = -1
    y.shape = -1

    # Sorted and regularly spaced values
    sortedX = x
    sortedY = y    

    ddict = getE0SavitzkyGolay(sortedX,
                               sortedY,
                               points=npoints,
                               full=full)
    if full:
        # return intermediate information
        return ddict["edge"], sortedX, sortedY, ddict["xPrime"], ddict["yPrime"]
    else:
        # return the corresponding x value
        return ddict

    def interpolate(self, factor=1.):
        """
        Input
        -----
        factor : float
            factor used to oversample existing data, use
            with caution.

        Interpolates all existing curves to an equidistant
        x-range using the either the active or the first
        curve do determine the number of data points.
        Use this method instead of self.getAllCurves() when
        performin FFT related tasks.

        Returns
        -------
        interpCurves : ndarray
            Array containing the interpolated curves shown
            in the plot window.
            Format: [(x0, y0, legend0, info0), ...]
        """
        curves = self.getAllCurves()
        if len(curves) < 1:
            if DEBUG == 1:
                print('interpolate -- no curves present')
            raise ValueError("At least 1 curve needed")
            return

        activeCurve = self.getActiveCurve()
        if not activeCurve:
            activeCurve = curves[0]
        else:
            activeLegend = activeCurve[2]
            idx = list.index([curve[2] for curve in curves],
                             activeLegend)
            activeCurve = curves[idx]
        activeX, activeY, activeLegend, activeInfo = activeCurve[0:4]

        # Determine average spaceing between Datapoints
        step = numpy.average(numpy.diff(activeX))
        xmin, xmax = self.getXLimits([x for (x,y,leg,info) in curves],
                                     overlap=False)
        num  = factor * numpy.ceil((xmax-xmin)/step)

        # Create equidistant x-range, exclude first and last point
        xeq = numpy.linspace(xmin, xmax, num, endpoint=False)[:-1]

        # Interpolate on sections of xeq
        interpCurves = []
        for (x,y,legend,info) in curves:
            idx = numpy.nonzero((x.min()<xeq) & (xeq<x.max()))[0]
            xi = numpy.take(xeq, idx)
            yi = SpecfitFuns.interpol([x], y, xi.reshape(-1,1), y.min())
            yi.shape = -1
            interpCurves += [(xi, yi, legend, info)]
        return interpCurves

 def periodic_gauss(self, pars, x):
     """
     Fit function periodic_gauss(pars, x)
     pars = [npeaks, delta, height, position, fwhm]
     """
     newpars = numpy.zeros((pars[0], 3), numpy.float)
     for i in range(int(pars[0])):
         newpars[i, 0] = pars[2]
         newpars[i, 1] = pars[3] + i * pars[1]
         newpars[:, 2] = pars[4]
     return SpecfitFuns.gauss(newpars,x)

 def _fitBkgSubtract(spectra, config=None, anchorslist=None, fitmodel=None):
     """Subtract brackground from data and add it to fit model
     """
     for k in range(spectra.shape[1]):
         # obtain the smoothed spectrum
         background = SpecfitFuns.SavitskyGolay(spectra[:, k],
                                 config['fit']['stripfilterwidth'])
         lastAnchor = 0
         for anchor in anchorslist:
             if (anchor > lastAnchor) and (anchor < background.size):
                 background[lastAnchor:anchor] =\
                         SpecfitFuns.snip1d(background[lastAnchor:anchor],
                                            config['fit']['snipwidth'],
                                            0)
                 lastAnchor = anchor
         if lastAnchor < background.size:
             background[lastAnchor:] =\
                     SpecfitFuns.snip1d(background[lastAnchor:],
                                        config['fit']['snipwidth'],
                                        0)
         spectra[:, k] -= background
         if fitmodel is not None:
             fitmodel[:, k] = background

    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

def test():
    import numpy
    from PyMca5.PyMcaMath.fitting import SpecfitFuns
    x = numpy.arange(1000.)
    data = numpy.zeros((50, 1000), numpy.float)

    #the peaks to be fitted
    p0 = [100., 300., 50.,
          200., 500., 30.,
          300., 800., 65]

    #generate the data to be fitted
    for i in range(data.shape[0]):
        nPeaks = 3 - i % 3
        data[i,:] = SpecfitFuns.gauss(p0[:3*nPeaks],x)

    oldShape = data.shape
    data.shape = 1,oldShape[0], oldShape[1]

    instance = StackSimpleFit()
    instance.setData(x, data)
    # TODO: Generate this file "on-the-fly" to be able to test everywhere
    instance.setConfigurationFile(r"C:\StackSimpleFit.cfg")
    instance.processStack()

 def apvoigt(self,pars,x):
    """
    Fit function.
    """
    return SpecfitFuns.apvoigt(pars,x)

 def slit(self,pars,x):
     """
     A fit function.
     """
     return 0.5*SpecfitFuns.slit(pars,x)

 def stepup(self,pars,x):
     """
     A fit function.
     """
     return SpecfitFuns.upstep(pars,x)

 def stepdown(self,pars,x):
     """
     A fit function.
     """
     return SpecfitFuns.downstep(pars,x)

 def splitpvoigt(self,pars,x):
    """
    Asymmetric pseudovoigt.
    """
    return SpecfitFuns.splitpvoigt(pars,x)