Python numpy.polyder函数代码示例

    def calculateCentroidMeasurements(self):
        self.X[self.badFrames, :] = ma.masked
        if not self.useSmoothingFilterDerivatives:
            self.v[1:-1] = (self.X[2:, :] - self.X[0:-2])/(2.0/self.frameRate)
        else:
            # use a cubic polynomial filter to estimate the velocity
            self.v = ma.zeros(self.X.shape)
            halfWindow = int(np.round(self.filterWindow/2.*self.frameRate))
            for i in xrange(halfWindow, self.v.shape[0]-halfWindow):
                start = i-halfWindow
                mid = i
                finish = i+halfWindow+1
                if not np.any(self.X.mask[start:finish,:]):
                    px = np.polyder(np.polyfit(self.t[start:finish]-self.t[mid],
                                               self.X[start:finish, 0], 3))
                    py = np.polyder(np.polyfit(self.t[start:finish]-self.t[mid],
                                               self.X[start:finish, 1], 3))
                    self.v[i,:] = [np.polyval(px, 0), np.polyval(py, 0)]
                else:
                    self.v[i,:] = ma.masked

        self.s = ma.sqrt(ma.sum(ma.power(self.v, 2), axis=1))
        self.phi = ma.arctan2(self.v[:, 1], self.v[:, 0])
        self.t[self.badFrames] = ma.masked
        self.X[self.badFrames, :] = ma.masked
        self.v[self.badFrames, :] = ma.masked
        self.s[self.badFrames] = ma.masked
        self.phi[self.badFrames] = ma.masked

def sin_poly(P):
  poly_size = len(P)+1

  temp = np.poly1d([0])
  coeffCos = np.poly1d([0])
  coeffSin = np.poly1d([0])
  tmpCoeffCos = np.poly1d([0])
  tmpCoeffSin = np.poly1d([0])

  cos0 = math.cos(P[0])
  sin0 = math.sin(P[0])

  coeffSin[0] = 1.0
  temp[0] = sin0

  # compute the derivative of P
  dP = np.polyder(P)

  facti = 1
  for i in xrange(1,poly_size):
    facti *= i

    tmpCoeffCos = np.polyder(coeffCos) + coeffSin * dP
    tmpCoeffSin = np.polyder(coeffSin) - coeffCos * dP

    coeffCos = tmpCoeffCos
    coeffSin = tmpCoeffSin

    temp[i] = (coeffCos[0] * cos0 + coeffSin[0] * sin0)/float(facti)

  return temp

def deltafactor_polyfit(volumes, energies):
    """
    This is the routine used to compute V0, B0, B1 in the deltafactor code.

    Taken from deltafactor/eosfit.py
    """
    fitdata = np.polyfit(volumes**(-2./3.), energies, 3, full=True)
    ssr = fitdata[1]
    sst = np.sum((energies - np.average(energies))**2.)
    residuals0 = ssr/sst
    deriv0 = np.poly1d(fitdata[0])
    deriv1 = np.polyder(deriv0, 1)
    deriv2 = np.polyder(deriv1, 1)
    deriv3 = np.polyder(deriv2, 1)

    v0 = 0
    x = 0
    for x in np.roots(deriv1):
        if x > 0 and deriv2(x) > 0:
            v0 = x**(-3./2.)
            break
    else:
        raise EOSError("No minimum could be found")

    derivV2 = 4./9. * x**5. * deriv2(x)
    derivV3 = (-20./9. * x**(13./2.) * deriv2(x) - 8./27. * x**(15./2.) * deriv3(x))
    b0 = derivV2 / x**(3./2.)
    b1 = -1 - x**(-3./2.) * derivV3 / derivV2

    #print('deltafactor polyfit:')
    #print('e0, b0, b1, v0')
    #print(fitdata[0], b0, b1, v0)

    n = collections.namedtuple("DeltaFitResults", "v0 b0 b1 poly1d")
    return n(v0, b0, b1, fitdata[0])

    def fit_sjeos(self):
        """Calculate volume, energy, and bulk modulus.

        Returns the optimal volume, the minimum energy, and the bulk
        modulus.  Notice that the ASE units for the bulk modulus is
        eV/Angstrom^3 - to get the value in GPa, do this::

          v0, e0, B = eos.fit()
          print(B / kJ * 1.0e24, 'GPa')

        """

        fit0 = np.poly1d(np.polyfit(self.v**-(1 / 3), self.e, 3))
        fit1 = np.polyder(fit0, 1)
        fit2 = np.polyder(fit1, 1)

        self.v0 = None
        for t in np.roots(fit1):
            if isinstance(t, float) and t > 0 and fit2(t) > 0:
                self.v0 = t**-3
                break

        if self.v0 is None:
            raise ValueError('No minimum!')

        self.e0 = fit0(t)
        self.B = t**5 * fit2(t) / 9
        self.fit0 = fit0

        return self.v0, self.e0, self.B

def BM(energies):

    fitdata = np.polyfit(energies[:, 0] ** (-2.0 / 3.0), energies[:, 1], 3, full=True)
    ssr = fitdata[1]
    sst = np.sum((energies[:, 1] - np.average(energies[:, 1])) ** 2.0)
    residuals0 = ssr / sst
    deriv0 = np.poly1d(fitdata[0])
    deriv1 = np.polyder(deriv0, 1)
    deriv2 = np.polyder(deriv1, 1)
    deriv3 = np.polyder(deriv2, 1)

    volume0 = 0
    x = 0
    for x in np.roots(deriv1):
        if x > 0 and deriv2(x) > 0:
            volume0 = x ** (-3.0 / 2.0)
            break

    if volume0 == 0:
        print("Error: No minimum could be found")
        exit()

    derivV2 = 4.0 / 9.0 * x ** 5.0 * deriv2(x)
    derivV3 = -20.0 / 9.0 * x ** (13.0 / 2.0) * deriv2(x) - 8.0 / 27.0 * x ** (15.0 / 2.0) * deriv3(x)
    bulk_modulus0 = derivV2 / x ** (3.0 / 2.0)
    bulk_deriv0 = -1 - x ** (-3.0 / 2.0) * derivV3 / derivV2

    return volume0, bulk_modulus0, bulk_deriv0, residuals0

def findStepWithPoly(x,data,kind="stepUp",excludePoints=100,order=20,fitrange=100):
  """ Look for a step in the data
      Data is 1D array
      The 'kind' keywords should be either 'stepUp' or 'stepDown'
      the 'excludePoints' keyword is used to limit the search 'excludePoints' 
      away from the extremes
      The data are fit with a polynomial of order 'order', then the maximum
      (or minimum) of derivative is located.
      After this first attempt, the position is refined by a second polyfit
      in a range [-fitrange,+fitrange] around the first guess
  """
  if (kind == "stepUp"):
    use = "max"
  else:
    use = "min"
  poly    = np.polyfit(x,data,order)
  polyder = np.polyder(poly)
  x_poly1 = findDerPeak(x,np.polyval(polyder,x),use=use,excludePoints=excludePoints)
  # find closest
  idx = np.abs(x - x_poly1).argmin()
  idx = slice(idx-fitrange,idx+fitrange)
  poly    = np.polyfit(x[idx],data[idx],order)
  polyder = np.polyder(poly)
  x_poly2 = findDerPeak(x[idx],np.polyval(polyder,x[idx]),use=use,excludePoints=10)
  return x_poly2

def calc_max_curvature_point(x_arr, y_arr):
    rot_x = []
    rot_y = []

    # Rotate all of the x and y co-ordinates around the origin by 90deg
    for current_x, current_y in zip(x_arr, y_arr):
        out_x, out_y = rotate(current_x, current_y, 90)
        rot_x.append(out_x)
        rot_y.append(out_y)

    pol = get_best_poly_fit(rot_x, rot_y, 7)

    # Differentiate the polynomial
    deriv_one = numpy.polyder(pol)
    deriv_two = numpy.polyder(deriv_one)

    # Roots of the 1st derivative - that is, X-values for where 1st deriv = 0
    roots = numpy.roots(deriv_one)

    if roots.size == 1:
        result = rotate(roots, numpy.polyval(pol, roots), -90)
    else:
        prev_val = 0
        selected_root = -1

        for root in roots:
            val = numpy.polyval(deriv_two, root)
            if val > prev_val:
                selected_root = root

        result = rotate(selected_root, numpy.polyval(pol, selected_root), -90)

    return [float(result[0]), float(result[1])]

def solve_for_nearest( px,py,rx,ry ):
  dpx = polyder(px)
  dpy = polyder(py)
  cp = polymul( dpx, px ) + polymul( dpy, py )
  cp = polyadd( cp, -rx*dpx )      
  cp = polyadd( cp, -ry*dpy )
  t = roots(cp)
  t = real(t[isreal(t)])
  t = t[ (t>=0) * (t<=1) ]

  ##tt = linspace(0,1,100)
  ##from pylab import plot
  ##plot( polyval(px,tt), polyval(py,tt), 'k', hold = 0 )
  ##plot( [rx],[ry], 'r.' )
  ##plot( polyval(px,t[isreal(t)*(real(t)>=0)*(real(t)<=1)]),
  ##      polyval(py,t[isreal(t)*(real(t)>=0)*(real(t)<=1)]), 'o' )
  ##pdb.set_trace()

  if len(t):
    if len(t) == 1:
      return t[0]
    else:
      ux = polyval( px, t )
      uy = polyval( py, t )
      d = hypot( ux - rx, uy - ry )
      return t[ d==d.min() ][0]
  else:
    t = array([0.0,1.0])
    ux = polyval( px, t )
    uy = polyval( py, t )
    d = hypot( ux - rx, uy - ry )
    if d[0] < d[1]:
      return 0.0
    else:
      return 1.0

def solve_pathlength_func_bendrad(p, edge_y, edge_dydx, D, L0, BendRad, n_grid=100):
    """A helper function for solve_pathlength where a fixed minimum bend radius 
    is desired.
    
    Parameters
    ----------
    BendRad: float
        Enforced minimum bend radius.
    """
    pathlength_poly = np.poly1d(p)
    p_d = np.polyder(pathlength_poly)
    p_d_d = np.polyder(p_d)
    L = integrate.quad(polynomial_pathlength, 0, D, args=(p_d))

    # Find the minimum radius of curvature along the curve.
    # Given the the curve is non-differentiable, try brute-force with n_grid points
    # along the curve.
    x_vect = np.meshgrid(0, D, n_grid)

    a = bend_radius(x_vect, p_d, p_d_d)

    # Find the minimum radius of curvature.
    SmallCurve = np.min(a)

    retpar = [
        pathlength_poly(0) - edge_y[0],
        pathlength_poly(D) - edge_y[1],
        p_d(0) - edge_dydx[0],
        p_d(D) - edge_dydx[1],
        L[0] - L0,
        SmallCurve - BendRad,
    ]
    # print(retpar); pdb.set_trace()
    return retpar

    def splitBimodal(self, x, y, largepoly=30):
        p = np.polyfit(x, y, largepoly) # polynomial coefficients for fit

        extrema = np.roots(np.polyder(p))
        extrema = extrema[np.isreal(extrema)]
        extrema = extrema[(extrema - x[1]) * (x[-2] - extrema) > 0] # exclude the endpoints due false maxima during fitting
        try:
            root_vals = [sum([p[::-1][i]*(root**i) for i in range(len(p))]) for root in extrema]
            peaks = extrema[np.argpartition(root_vals, -2)][-2:] # find two peaks of bimodal distribution

            mid, = np.where((x - peaks[0])* (peaks[1] - x) > 0)
             # want data points between the peaks
        except:
            warnings.warn("Peak finding failed!")
            return None

        try:
            p_mid = np.polyfit(x[mid], y[mid], 2) # fit middle section to a parabola
            midpoint = np.roots(np.polyder(p_mid))[0]
        except:
            warnings.warn("Polynomial fit between peaks of distribution poorly conditioned. Falling back on using the minimum! May result in inaccurate split determination.")
            if len(mid) == 0:
                return None

            midx = np.argmin(y[mid])
            midpoint = x[mid][midx]

        return midpoint

def poly_sigma(time, y):
    """
    Calculates velocity and acceleration by fitting a polynomial over the prior N elements
    defined by poly_window. The sigma value defines an averaging window inside of poly_window
    in which all points inside of sigma window are averaged and treated as a single point for the
    polynomial fitting process.
    """
    y_d = np.zeros(time.shape)
    y_dd = np.zeros(time.shape)
    window = poly_window

    for i in range(window * sigma, time.shape[0]):
        y_history = y[i - window * sigma + 1:i + 1]
        y_hist_avg = np.mean(y_history.reshape(-1, sigma), axis=1)
        params = np.polyfit(
            x=time[i - window * sigma + 1:i + 1:sigma], y=y_hist_avg, deg=degree)

        p = np.poly1d(params)
        p_d = np.polyder(p)
        p_dd = np.polyder(p_d)

        y_d[i] = p_d(time[i])
        y_dd[i] = p_dd(time[i])

    y_d = y_d / encoder_resolution
    y_dd = y_dd / encoder_resolution

    return y_d, y_dd

    def _set_params(self):
        """
        Overriden to account for the fact the fit with volume**(2/3) instead
        of volume.
        """
        deriv0 = np.poly1d(self.eos_params)
        deriv1 = np.polyder(deriv0, 1)
        deriv2 = np.polyder(deriv1, 1)
        deriv3 = np.polyder(deriv2, 1)

        for x in np.roots(deriv1):
            if x > 0 and deriv2(x) > 0:
                v0 = x**(-3./2.)
                break
        else:
            raise EOSError("No minimum could be found")

        derivV2 = 4./9. * x**5. * deriv2(x)
        derivV3 = (-20./9. * x**(13./2.) * deriv2(x) - 8./27. *
                   x**(15./2.) * deriv3(x))
        b0 = derivV2 / x**(3./2.)
        b1 = -1 - x**(-3./2.) * derivV3 / derivV2

        # e0, b0, b1, v0
        self._params = [deriv0(v0**(-2./3.)), b0, b1, v0]

def fit_edge_hist(bins, counts, fwhm_guess=10.0):
    if len(bins) == len(counts)+1: bins = bins[:-1]+0.5*(bins[1]-bins[0]) # convert bin edge to bin centers if neccesary
    pfit = np.polyfit(bins, counts, 3)
    edgeGuess = np.roots(np.polyder(pfit,2))
    try:
        preGuessX, postGuessX = np.sort(np.roots(np.polyder(pfit,1)))
    except:
        raise ValueError("failed to generate guesses")
    use = bins>(edgeGuess+2*fwhm_guess)
    if np.sum(use)>4:
        pfit2 = np.polyfit(bins[use], counts[use],1)
        slope_guess = pfit2[0]
    else:
        slope_guess=1
    pGuess = np.array([edgeGuess, np.polyval(pfit,preGuessX), np.polyval(pfit,postGuessX),fwhm_guess,slope_guess],dtype='float64')

    try:
        pOut = curve_fit(edge_model, bins, counts, pGuess)
    except:
        return (0,0,0,0,0,0)
    (edgeCenter, preHeight, postHeight, fwhm, bgSlope) = pOut[0]
    model_counts = edge_model(bins, edgeCenter, preHeight, postHeight, fwhm, bgSlope)
    num_degree_of_freedom = float(len(bins)-1-5) # num points - 1 - number of fitted parameters
    chi2 = np.sum(((counts - model_counts)**2)/model_counts)/num_degree_of_freedom
    return (edgeCenter, preHeight, postHeight, fwhm, bgSlope, chi2)

def second_deriv_from_profile(xinterp,F):
    """Calculate second derivative at extrema"""
    dFdx = np.polyder(F,m=1)
    d2Fdx2 = np.polyder(F,m=2)
    minidx, maxidx = extrema_from_profile(xinterp,F)
    omegamin = d2Fdx2(xinterp[minidx])
    omegamax = d2Fdx2(xinterp[maxidx])
    return omegamin, omegamax

def polykruemm(p,x):
    y = np.polyval(p,x)
    dy = np.polyval(np.polyder(p),x)
    ddy = np.polyval(np.polyder(p,2),x)
    kappa = np.abs(ddy/(1+dy**2)**(3/2))
    mx = x - dy*(1+dy**2)/ddy
    my = y + (1+dy**2)/ddy
    return kappa,mx,my