def generalized_gauss(collocation_points, verbose=False):
    p, result1, result2 = args.degree, [], []
    interp = BarycentricInterpolator(collocation_points)
    for n in range(p//2):
        filename = 'weights{}.txt'.format(n+1)
        with open(filename) as f:
            header = f.readline()
            x = float(header.split()[-1])
            # check if quadrature is for the given collacation point
            assert math.isclose(x, collocation_points[n], rel_tol=fixed.exact_tol)
        y, w = np.loadtxt(filename, unpack=True)
        check_generalized(lambda f: np.sum(w*f(y)), x, verbose)
        I, b1, b2 = np.eye(collocation_points.size), [], []
        for i in range(collocation_points.size):
            interp.set_yi(I[i])
            b1.append(interp(y))
            b2.append(interp(1-y[::-1]))
        result1.append(( y, w, np.array(b1) ))
        result2.append(( 1-y[::-1], w[::-1], np.array(b2) ))
    return result1 + result2[::-1]

 def test_append(self):
     P = BarycentricInterpolator(self.xs[:3],self.ys[:3])
     P.add_xi(self.xs[3:],self.ys[3:])
     assert_almost_equal(self.true_poly(self.test_xs),P(self.test_xs))

 def test_delayed(self):
     P = BarycentricInterpolator(self.xs)
     P.set_yi(self.ys)
     assert_almost_equal(self.true_poly(self.test_xs),P(self.test_xs))

npts = 7
noiseAmp = 0.03
xdata = np.linspace(0.2, 10., npts) + 0.3*np.random.randn(npts)
xdata = np.sort(xdata)
ydata = f(xdata) * (1.0 + noiseAmp * np.random.randn(npts))
np.savetxt('logCurveData9.txt', zip(xdata, ydata), fmt='%10.2f')

# Create x array spanning data set plus 5%
xmin, xmax = frangeAdd(xdata, 0.05)
x = np.linspace(xmin, xmax, 200)

# Plot data and "fit"
plt.plot(xdata, ydata, 'or', label='data')
plt.plot(x, f(x), 'k-', label='fitting function')

# Create y array from cubic spline
f_cubic = interp1d(xdata, ydata, kind='cubic', bounds_error=False)
plt.plot(x, f_cubic(x), 'k-', label='cubic spline')

# Create y array from univariate spline
f_univar = UnivariateSpline(xdata, ydata, w=None, bbox=[xmin, xmax], k=3)
plt.plot(x, f_univar(x), label='univariate spline')

# Create y array from barycentric interpolation
f_bary = BarycentricInterpolator(xdata, ydata)
plt.plot(x, f_bary.__call__(x), label='barycentric interp')

plt.legend(loc='best')

plt.show()

 def interpolator(self, x, values):
     """
     Returns a polynomial with vector coefficients which interpolates the values at the Chebyshev points x
     """
     # hacking the barycentric interpolator by computing the weights in advance
     p = Bary([0.])
     N = len(values)
     weights = np.ones(N)
     weights[0] = .5
     weights[1::2] = -1
     weights[-1] *= .5
     p.wi = weights
     p.xi = x
     p.set_yi(values)
     return p

t0_a= 10.8
sigma_t_a= 0.2
sigma= 0.04
active_space_a= np.array(map(lambda x: 1./(1+np.exp((x-L)/sigma)),x))*1./(1+np.exp(t0_a/sigma_t_a))
active_const_a= np.array(map(lambda x: 1.,x))*1./(1+np.exp(t0_a/sigma_t_a))
t0_i= 12.2
sigma_t_i= .8
sigma= 0.04
active_space_i= np.array(map(lambda x: 1./(1+np.exp((x-L)/sigma)),x))*1./(1+np.exp(t0_i/sigma_t_i))
active_const_i= np.array(map(lambda x: 1.,x))*1./(1+np.exp(t0_i/sigma_t_i))

rng= np.linspace(-1,1,1000)
x_step= 2./1000
interpolator= BI(x)
avg_h_blade= []
avg_Q_minus_blade= []
avg_T_minus_blade= []
avg_shear_blade= []

zeta= zeta_h*active_space_a + zeta_b*active_const_a
zetabar= zetabar_h*active_space_i + zetabar_b*active_const_i
lambda2= lambda2_h*active_space_a + lambda2_b*active_const_a
#xi = xi_h*active_space_a + xi_b*active_const_a

n_images= 150
image_step= 200

for i in np.arange(n_images*image_step):
    if i%image_step == 0: