Python numpy 模块，polyder() 实例源码

def calc_approx(time, altitude, div = 5, dt = 1.0, deg = 100, lim = 1):
    time_list, altitude_list, velocity_y = [],[],[]

    for i in range(div):
        f = trendline.approx_func(
            time[max(int((i - lim) * len(time) / div), 0):min(len(time), int((i + lim) * len(time) / div))],
            altitude[max(int((i - lim) * len(time) / div), 0):min(len(time), int((i + lim) * len(time) / div))], deg)
        der = np.polyder(f)

        #print(time[int(i * len(time) / div)], min(time[-1], time[int((i + 1) * len(time) / div) - 1]))

        for t in np.arange(time[int(i * len(time) / div)], min(time[-1], time[int((i + 1) * len(time) / div) - 1]), dt):
        #for t in time:
        #    if t >= time[int(i * len(time) / div)] and t < min(time[-1], time[int((i + 1) * len(time) / div) - 1]):
            time_list.append(t)
            altitude_list.append(f(t))
            velocity_y.append(der(t))

    return time_list, velocity_y, altitude_list

def calc_adot(cp,order=1):
    a=[]
    for i in range(len(cp)):
        ptmp = []
        tmp = 0
        for j in range(len(cp)):
            if j != i:
                row = []
                row.insert(0,1/(cp[i]-cp[j]))
                row.insert(1,-cp[j]/(cp[i]-cp[j]))
                ptmp.insert(tmp,row)
                tmp += 1
        p=[1]
        for j in range(len(cp)-1):
            p = conv(p,ptmp[j])
        pder = numpy.polyder(p,order)
        arow = []
        for j in range(len(cp)):
            arow.append(numpy.polyval(pder,cp[j]))
        a.append(arow)
    return a

def test_polyder_return_type(self):
        # Ticket #1249
        assert_(isinstance(np.polyder(np.poly1d([1]), 0), np.poly1d))
        assert_(isinstance(np.polyder([1], 0), np.ndarray))
        assert_(isinstance(np.polyder(np.poly1d([1]), 1), np.poly1d))
        assert_(isinstance(np.polyder([1], 1), np.ndarray))

def compute_grun_along_one_direction(nq,modes,ngeo,cgeo,celldmsx,freqgeo,rangegeo,xindex=0):
    """
    Compute the Gruneisen parameters along one direction.
    This function uses a 1-dimensional polynomial of fourth degree to fit the 
    frequencies along a certain direction (along a and c axis in hexagonal systems
    for example). 
    """
    # set a numpy array of volumes for the fit (n=5)
    xtemp=[]
    for igeo in rangegeo:
        xtemp.append(celldmsx[igeo,xindex])
    x=np.array(xtemp)

    grun=[]
    for iq in range(0,nq):
        grunq=[]
        for ifreq in range(0,modes):
            ytemp=[]
            for igeo in rangegeo:
                ytemp.append(freqgeo[igeo,iq,ifreq])
            y=np.array(ytemp)    
            z=np.polyfit(x, y, 4)
            p=np.poly1d(z)
            pderiv=np.polyder(p)
            if freqgeo[cgeo[xindex],iq,ifreq]<1E-3:
                grunq.append(0.0)
            else:
                grunq.append(pderiv(celldmsx[cgeo[xindex],xindex])/freqgeo[cgeo[xindex],iq,ifreq]) #*celldmsx[cgeo[xindex],xindex])
        grun.append(grunq)

    return np.array(grun)    


################################################################################

def test_polyder_return_type(self):
        # Ticket #1249
        assert_(isinstance(np.polyder(np.poly1d([1]), 0), np.poly1d))
        assert_(isinstance(np.polyder([1], 0), np.ndarray))
        assert_(isinstance(np.polyder(np.poly1d([1]), 1), np.poly1d))
        assert_(isinstance(np.polyder([1], 1), np.ndarray))

def basis_polys():
    """ Returns the basis polynomials and their derivatives and antiderivatives
    """
    nodes, _, _ = quad()
    ? = [lagrange(nodes,eye(N+1)[i]) for i in range(N+1)]
    ?Der = [[polyder(?_p, m=a) for ?_p in ?] for a in range(N+1)]
    ?Int = [polyint(?_p) for ?_p in ?]
    return ?, ?Der, ?Int

def test_polyder_return_type(self):
        # Ticket #1249
        assert_(isinstance(np.polyder(np.poly1d([1]), 0), np.poly1d))
        assert_(isinstance(np.polyder([1], 0), np.ndarray))
        assert_(isinstance(np.polyder(np.poly1d([1]), 1), np.poly1d))
        assert_(isinstance(np.polyder([1], 1), np.ndarray))

def test_polyder_return_type(self):
        # Ticket #1249
        assert_(isinstance(np.polyder(np.poly1d([1]), 0), np.poly1d))
        assert_(isinstance(np.polyder([1], 0), np.ndarray))
        assert_(isinstance(np.polyder(np.poly1d([1]), 1), np.poly1d))
        assert_(isinstance(np.polyder([1], 1), np.ndarray))

def test_polyder_return_type(self):
        # Ticket #1249
        assert_(isinstance(np.polyder(np.poly1d([1]), 0), np.poly1d))
        assert_(isinstance(np.polyder([1], 0), np.ndarray))
        assert_(isinstance(np.polyder(np.poly1d([1]), 1), np.poly1d))
        assert_(isinstance(np.polyder([1], 1), np.ndarray))

def polyval(fit, points, der=0, avg=False):
    """Evaluate polynomial generated by ``polyfit()`` on `points`.

    Parameters
    ----------
    fit, points : see polyfit()
    der : int, optional
        Derivative order. Only for 1D, uses np.polyder().
    avg : bool, optional
        Internal hack, only used by ``avgpolyval()``.

    Notes
    -----
    For 1D we provide "analytic" derivatives using np.polyder(). For ND, we
    didn't implement an equivalent machinery. For 2D, you might get away with
    fitting a bispline (see Interpol2D) and use it's derivs. For ND, try rbf.py's RBF
    interpolator which has at least 1st derivatives for arbitrary dimensions.

    See Also
    --------
    :class:`PolyFit`, :class:`PolyFit1D`, :func:`polyfit`
    """
    pscale, pmin = fit['pscale'], fit['pmin']
    vscale, vmin = fit['vscale'], fit['vmin']
    if der > 0:
        assert points.shape[1] == 1, "deriv only for 1d poly (ndim=1)"
        # ::-1 b/c numpy stores poly coeffs in reversed order
        dcoeffs = np.polyder(fit['coeffs'][::-1], m=der)
        return np.polyval(dcoeffs, (points[:,0] - pmin[0,0]) / pscale[0,0]) / \
            pscale[0,0]**der * vscale
    else:
        vand = vander((points - pmin) / pscale, fit['deg'])
        if avg:
            return np.dot(vand, fit['coeffs']) * vscale
        else:
            return np.dot(vand, fit['coeffs']) * vscale + vmin

def calc_approx_func(time, altitude, div = 5, der = 0 ,deg = 100, lim = 1):
    dict = {}

    for i in range(div):
        f = trendline.approx_func(
            time[max(int((i - lim) * len(time) / div), 0):min(len(time), int((i + lim) * len(time) / div))],
            altitude[max(int((i - lim) * len(time) / div), 0):min(len(time), int((i + lim) * len(time) / div))], deg)

        for j in range(der):
            f = np.polyder(f)

        dict[min(time[-1], time[int((i + 1) * len(time) / div) - 1])] = f

    return dict

def test_polyder_return_type(self):
        # Ticket #1249
        assert_(isinstance(np.polyder(np.poly1d([1]), 0), np.poly1d))
        assert_(isinstance(np.polyder([1], 0), np.ndarray))
        assert_(isinstance(np.polyder(np.poly1d([1]), 1), np.poly1d))
        assert_(isinstance(np.polyder([1], 1), np.ndarray))

def birch_murnaghan_fit(energies, volumes):
    """
    least squares fit of a Birch-Murnaghan equation of state curve. From delta project
    containing in its columns the volumes in A^3/atom and energies in eV/atom
    # The following code is based on the source code of eos.py from the Atomic
    # Simulation Environment (ASE) <https://wiki.fysik.dtu.dk/ase/>.
    :params energies: list (numpy arrays!) of total energies eV/atom
    :params volumes: list (numpy arrays!) of volumes in A^3/atom

    #volume, bulk_modulus, bulk_deriv, residuals = Birch_Murnaghan_fit(data)
    """
    fitdata = np.polyfit(volumes[:]**(-2./3.), energies[:], 3, full=True)
    ssr = fitdata[1]
    sst = np.sum((energies[:] - np.average(energies[:]))**2.)
    #print(fitdata)
    #print(ssr)
    #print(sst)
    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./2.)
            break

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

    derivV2 = 4./9. * x**5. * deriv2(x)
    derivV3 = (-20./9. * x**(13./2.) * deriv2(x) -
               8./27. * x**(15./2.) * deriv3(x))
    bulk_modulus0 = derivV2 / x**(3./2.)
    #print('bulk modulus 0: {} '.format(bulk_modulus0))
    bulk_deriv0 = -1 - x**(-3./2.) * derivV3 / derivV2

    return volume0, bulk_modulus0, bulk_deriv0, residuals0