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

def matrixN(tau, rows=-1, last_value=1.0):
    n = tau.shape[0]
    if rows == -1:
        rows = n
    N = np.zeros((rows, n))
    # construct the lagrange polynomials
    circulating_one = np.asarray([1.0]+[0.0]*(n-1))
    lag_pol = []
    for i in range(n):
        lag_pol.append(intpl.lagrange(tau, np.roll(circulating_one, i)))
        N[:, i] = -np.ones(rows)*lag_pol[-1](last_value)
    return N

 def updateApproximation(self):
   self.logger.debug('')
   self.logger.debug('entering updateApproximation')
   self.logger.debug('self.x = \n%s'  % self.x)
   self.logger.debug('self.fx = \n%s' % self.fx)
   if len(self.x[0]) == 1: # we are dealing with the one-dimensional case
     self.logger.debug('self.x[0] == 1')
     xIn = self.sorted_x.copy()
     fIn = self.sorted_fx.copy()
     withinDistance = ( self.sortedNormedDistances[1] / self.sortedNormedDistances ) < 10
     withinDistance[0] = True # need at least two values
     withinDistance[1] = True # 
     if self.protected_ix:
       self.logger.debug('protected exists')
       try: 
         lenProtected_ix = len(self.protected_ix)
       except:
         lenProtected_ix = 1
       actual_nIntpoints_max = self.nIntpoints_max - lenProtected_ix
       withinDistance[self.protected_ix] = True # protected_ix must always be carried over
     else:
       self.logger.debug('protected doesnt exist')
       actual_nIntpoints_max = self.nIntpoints_max
     actual_nIntpoints = np.min([actual_nIntpoints_max, len(xIn)])
     self.logger.debug('actual_nIntpoints = %s' % actual_nIntpoints)
     self.logger.debug('len(xIn) = %s' % len(xIn))
     withinMaxIntpoints = np.array([ True ] * actual_nIntpoints + [ False ] * ( len(xIn) - actual_nIntpoints ) )
     withinMaxIntpoints[self.protected_ix] = True
     self.logger.debug('withinMaxIntpoints = %s' % withinMaxIntpoints)
     self.logger.debug('withinDistance = %s' % withinDistance)
     self.intPoints_ix = withinMaxIntpoints * withinDistance
     xIn = xIn[self.intPoints_ix]
     fIn = fIn[self.intPoints_ix]
     self.logger.debug('after nIntpoints  xIn = %s \n' % xIn)
     self.logger.debug('after nIntpoints  fIn = %s \n' % fIn)
     self.logger.debug('transforming fIn:')
     fIn = np.log((1. - self.target_fx) + fIn)
     self.logger.debug('sending xIn = %s to lagrange\n' % xIn)
     self.logger.debug('sending fIn = %s to lagrange\n' % fIn)
     ### create approximation ###
     self.gx = si.lagrange(xIn[:, 0], fIn)
     ### -------------------- ###
   else:
     self.logger.debug('multi-dimensional case not implemented... exiting')
     os._exit(1)
   if self.makePlots:
     self.plot()
   self.logger.debug('done updating approximation')
   self.logger.debug('')

    def _interpolate_solution(self, ts):

        h = self.tf / self.nk
        stage_starts = pd.Series(h * np.arange(self.nk))
        stages = stage_starts.searchsorted(ts, side='right') - 1

        out = np.empty((len(ts), self.nx))
    
        for ki in range(self.nk):
            for ni in range(self.nx):
                interp = lagrange(self.col_vars['tau_root'], 
                                  self.var.x_op[ki, :, ni])

                out[stages == ki, ni] = interp(
                    (ts[stages == ki] - stage_starts[ki])/h)

        return out

def interpolate_to_t_end(nodes_on_unit, values):
    """
    Assume a GaussLegendre nodes, we are interested in the value at the end of
    the interval, but we now only the values in the interior of the interval.
    We compute the value by legendre interpolation.
    :param nodes_on_unit: nodes transformed to the unit interval
    :param values: values on those nodes
    :return: interpolation to the end of the interval
    """
    n = nodes_on_unit.shape[0]
    circulating_one = np.asarray([1.0]+[0.0]*(n-1))
    lag_pol = []
    result = np.zeros(values[0].shape)
    for i in range(n):
        lag_pol.append(intpl.lagrange(nodes_on_unit, np.roll(circulating_one, i)))
        result += values[i]*lag_pol[-1](1.0)
    return result

def integrate1d_romberg(f,a,b,N):
    if a == b:
        return 0
    k = np.log2(N)
    if not k % 1 == 0:
        raise ValueError("N has to be a power of 2")
    k = int(k)
    h = b - a
    tf = h*(f(a)+f(b))/2.
    h_list = [h]
    tf_list = [tf]
    for i in range(0,k):
        tf = 0.5*tf + 0.5 * h * sum(f(a+(np.arange(0,2**i) + 0.5)*h))
        h *= 0.5
        tf_list.append(tf)
        h_list.append(h)
    
    tf = np.array(tf_list)
    h = np.array(h_list)
    return lagrange(h**2, tf)(0)