Python interpolate.lagrange函数代码示例

def interpolation_matrix_1d(fine_grid, coarse_grid, k=2, return_type="csc", periodic=False, T=1.0):
    """
    We construct the interpolation matrix between two 1d grids, using lagrange interpolation.

    :param fine_grid: a one dimensional 1d array containing the nodes of the fine grid
    :param coarse_grid: a one dimensional 1d array containing the nodes of the coarse grid
    :param k: order of the restriction
    :return: a interpolation matrix
    """
    M = np.zeros((fine_grid.size, coarse_grid.size))
    n_f = fine_grid.size

    for i, p in zip(range(n_f), fine_grid):
        if periodic:
            nn,cont_arr = next_neighbors_periodic(p, coarse_grid, k, T)
            circulating_one = np.asarray([1.0]+[0.0]*(k-1))
            lag_pol = []
            for l in range(k):
                lag_pol.append(intpl.lagrange(cont_arr, np.roll(circulating_one, l)))
            M[i, nn] = np.asarray(map(lambda x: x(p), lag_pol))
        else:
            nn = next_neighbors(p, coarse_grid, k)
        # construct the lagrange polynomials for the k neighbors
            circulating_one = np.asarray([1.0]+[0.0]*(k-1))
            lag_pol = []
            for l in range(k):
                lag_pol.append(intpl.lagrange(coarse_grid[nn], np.roll(circulating_one, l)))
            M[i, nn] = np.asarray(map(lambda x: x(p), lag_pol))
    return to_sparse(M, return_type)

def print_table_at_segment_for_g_and_fs(x_segment, xs):
    fs = map(lambda x: mega_f(x), xs[:])
    gs = map(lambda x: mega_g(x), xs[:])
    #using the x_i = x_i to make python to close over the value of x_i and not, the name that is "lazy watched" once
    omega = [(lambda x, x_i=x_i : x - x_i) for x_i in xs]
    f_lagrange = construct_Lagrange_polynomial(xs, fs, omega)
    g_lagrange = construct_Lagrange_polynomial(xs, gs, omega)
    p = sci.lagrange(xs, fs)
    q = sci.lagrange(xs, gs)
    #checking the results
    for i in range(len(xs)):
        assert abs(eval_Lagrange(g_lagrange, xs[i]) - gs[i]) < float_difference
        assert p(xs[i]) - eval_Lagrange(f_lagrange, xs[i]) < float_difference
        assert q(xs[i]) - eval_Lagrange(g_lagrange, xs[i]) < float_difference  
    # our interpolation degree 
    n = len(omega)
    # prepare lambda functions for absolute of (n+1 derivatives) of f and g
    abs_nth_der_f = lambda x, n=n, der1=nth_dif_of_cos(n + 1, float(1)/3), der2=nth_dif_of_sin(n + 1, float(1)/2):    abs(der1(x) - der2(x))
    abs_nth_der_g = lambda x, n=n, g1=nth_dif_of_cos(n + 1, float(5)) :     abs(g1(x))
    #calculate their max values
    max_f_der_value = stupid_max(abs_nth_der_f, x_segment)
    assert max_f_der_value >= 0
    max_g_der_value = stupid_max(abs_nth_der_g, x_segment)
    assert max_g_der_value
    print abs_nth_der_f(0) >= 0
    print "Max values: |f`", n,"|", max_f_der_value, " ; |g`",n,"|", max_g_der_value
    #construct As
    A_f = lambda x, omega=omega : (abs(numpy.prod(map(lambda factor : factor(x), omega))) * max_f_der_value) / float(factorial(n + 1))
    A_g = lambda x : (abs(numpy.prod(map(lambda factor : factor(x), omega))) * max_g_der_value) / float(factorial(n + 1))
    # prepare yourself for iterating
    h = float(x_segment[-1] - x_segment[0]) / float(part_amount)
    print h
    #Print the target tables
    header = "-" * 11 + "Table for F at" + str(x_segment) + "-" *11
    print header
    print "-" * len(header)
    print_table(x_segment, h, mega_f, f_lagrange, A_f)
    print "-" * len(header)
    print
    header = "-" * 11 + "Table for G at" + str(x_segment) + "-" * 11
    print header
    print "-" * len(header)
    print_table(x_segment, h, mega_g, g_lagrange, A_g)
    print "-" * len(header)
    #plotting
    f_eval = lambda x : mega_f(x)
    g_eval = lambda x : mega_g(x)
    f_lagr_eval = lambda x: eval_Lagrange(f_lagrange, x)
    g_lagr_eval = lambda x: eval_Lagrange(g_lagrange, x)
    #f_lagr_eval = lambda x, f=f_lagrange: numpy.sum(map(lambda factor: factor(x), f_lagrange))
    mpmath.plot([f_eval, f_lagr_eval], x_segment)
    mpmath.plot([g_eval, g_lagr_eval], x_segment)

def test_interp_telles():
    # The problem I solve:
    # exact is from mathematica
    sing_pt = 1.005;
    denom = lambda x: (sing_pt - x) ** 2;
    numer = lambda x: x ** 3;
    f = lambda x: numer(x) / denom(x);
    exact_f = lambda s: (s * (-4 + 6 * s ** 2 - 3 * s * (-1 + s ** 2) * \
                                (log((-1 - s) / (1 - s))))) / (-1 + s ** 2)
    exact = exact_f(sing_pt)

    # Solved with standard Telles quadrature
    x_nearest = 1.0;
    D = sing_pt - 1.0;
    N = 14;
    [tx, tw] = telles_quasi_singular(N, x_nearest, D);
    est_telles = np.sum(f(tx) * tw)

    # Solved with gauss quadrature
    [gx, gw] = gaussxw(N)
    est_gauss = np.sum(f(gx) * gw)

    # Solved with interpolation and Telles quadrature
    # X = Interpolation Points
    X = gx;
    # Y = Value of function at the interpolation points
    Y = f(gx) * denom(gx);
    # WARNING, WARNING, WARNING: This implementation of lagrange interpolation
    # is super unstable. I just downloaded it from somewhere online. A
    # reimplementation using barycentric Lagrange interpolation is necessary to
    # go above N = appx 20.
    P = spi.lagrange(X, Y)

    est_interp_telles = sum(P(tx) / denom(tx) * tw)
    np.testing.assert_almost_equal(est_telles, est_interp_telles)

def return_polynomial_coefficients(curve_list):
	xdata = [x[0] for x in curve_list]
	ydata = [x[1] for x in curve_list]
	np.set_printoptions(precision=6)
	np.set_printoptions(suppress=True)
	p = interpolate.lagrange(xdata, ydata)
	return p

def ployinterp_column(s,n,k=5):
    # 取数
    y = s[list(range(n-k,n))+list(range(n+1,n+1+k))]
    # 删除空值
    y = y[y.notnull()]
    # 插值并返回插值结果
    return lagrange(y.index,list(y))(n)

    def _pd(self):
        m2, m3, nd, ld, type = self._m2, self._m3, self._nd,  self._ld, self._type
        if type == 1:
            dx = ld/(nd-1)
            p02 = Point(0, 1.0)
            dy1 = m2*dx
            p12 = p02 + Point(dx, dy1, 0.0, 0.0)
            pn2 = Point( ld, self._Aout)
            dy2 = m3*dx
            pn_12 = pn2 - Point(dx, dy2, 0.0, 0.0)
            pint =  [p02, p12, pn_12, pn2]
            n = len(pint)
            x = np.ones(n)
            y = np.ones(n)
            for i in xrange(n):
                x[i], y[i], nul1, nul2 = pint[i].split()
            xp = np.ones(nd)
            for i in xrange(nd):
                xp[i] =  i*dx     
        
        
        from scipy import interpolate
        f = interpolate.lagrange(x, y)
        yp = f(xp)
#         print x, y
#         print xp, yp
            
#         f = interpolate.interp1d(x, y)
        pct = []
        for i in xrange(nd):
            pct.append(Point(xp[i], yp[i]))
#         pct.append(pn1)  
        return pct

    def _interpolate_boundary_constraints(self, ts):

        stage_starts = [0.]
        for i in xrange(self.nk-1):
            stage_starts += [self.var.h_op[self._get_stage_index(i)] +
                             stage_starts[-1]]
        stage_starts = pd.Series(stage_starts)
        stages = stage_starts.searchsorted(ts, side='right') - 1

        for ki in range(self.nk):
            for ji in xrange(1, self.d+1):

                x = {met : var_op for met, var_op in 
                     zip(self.boundary_species, self.var.x_op[k,j])}

                interp = lagrange(self.col_vars['tau_root'], 
                                  (self.col_vars['C'].T.dot(
                                      self.var.x_op[ki, :, ni]) /
                                   self.var.h_op[self._get_stage_index(ki)]))

                out[stages == ki, ni] = interp(
                    (ts[stages == ki] - stage_starts[ki]) /
                    self.var.h_op[self._get_stage_index(ki)])

        return out

def intgl_simp38(f,a,b,steps=-1,h=1):
    if steps>0:
        xis = np.linspace(a,b,steps+1)
        h  = xis[1]-xis[0]
    fxis = f(xis)
    wis = np.zeros(steps+1)
    pcs = []; fpcs = []
    for i in xrange(0,steps-2,3):
        wis[i:i+4] += [1,3,3,1]
        pcs.append(xis[i:i+4])
        fpcs.append(fxis[i:i+4])
    wis *= 3*h/8
    if steps%3==2:
        wis[-3:] += [h/3,4*h/3,h/3]
        pcs.append(xis[-3:])
        fpcs.append(fxis[-3:])
    elif steps%3==1:
        wis[-2:] += [h/2,h/2]
        pcs.append(xis[-2:])
        fpcs.append(fxis[-2:])
    fapprox = lambda x: np.piecewise(x,
                                     [np.logical_and(p[0]<=x,x<=p[-1]) for p in pcs],
                                     [lagrange(pcs[i],fpcs[i]) for i in xrange(len(pcs))])# np.interp(x,xis,fxis)
    # fapprox = lambda x: np.interp(x,xis,fxis)
    return (sum(fxis*wis),xis,fxis,wis,fapprox) # h/2 * sum(np.array([f(x) for x in xs]) * np.array([1]+[2]*(len(xs)-2)+[1]))

 def _find_coefficients(self):
   polynomials = []
   for curve in self.curves:
     xdata = [x[0] for x in curve]
     ydata = [x[1] for x in curve]
     p = interpolate.lagrange(xdata, ydata)
     polynomials.append(p)
   return polynomials

def intgl_glquad(f,a,b,n):
    ans = spintegrate.fixed_quad(f,a,b,n=n)
    xis,wis = np.polynomial.legendre.leggauss(n)
    # print(xis,wis)
    xis = (b+a)/2 + (b-a)/2*xis
    fxis = f(xis)
    wis = (b-a)/2*wis
    return (ans[0],xis,fxis,wis,lagrange(xis,fxis))

    def lagrangeBasis(self, x):
        "Returns lagrange interpolant which has value 1 in x"

        xGrid = self.getList()
        yGrid = [0]*len(self.wnodes)
        i = xGrid.index(x)
        yGrid[i] = 1
        return lagrange(xGrid, yGrid)

def card_poly(k):
	x = []
	y = []
	for i in range(k + 1):
		n = k + 2 + i
		x.append(n)
		y.append(Ank.card_ank(n, k))
	p = si.lagrange(x, y)
	return p

 def _interpol_init_guess(self):
     Pig = self._Pig
     n = len(Pig)
     x = np.ones(n)
     y = np.ones(n)
     for i in xrange(n):
         x[i], y[i], nul1, nul2 = Pig[i].split()
     from scipy import interpolate
     f = interpolate.lagrange(x, y)
     self._yp = f(self._xp)     

 def interpolate(self, mode=InterpolateMode.LAGRANGE):
     xy = self.asarray()
     if mode == InterpolateMode.LAGRANGE or mode is None:
         delegate = interpolate.lagrange(xy.T[0],
                                         xy.T[1])
     else:
         kind = InterpolateMode(mode).to_string()
         delegate = interpolate.interp1d(xy.T[0],
                                         xy.T[1], kind=kind)
     self.delegate = delegate
     return self

def lazy():
    from scipy import interpolate

    x = list(range(1, 11))
    y = [u(n) for n in range(1, 11)]

    FITs = 0
    for d in range(1, 11):
        z = interpolate.lagrange(x[:d], y[:d])
        FITs += round(z(d+1))

    return FITs