def __init__(self, knots_x, knots_y, n_points):
        self.knots_x = knots_x
        self.knots_y = knots_y

        tck = splrep(knots_x, knots_y)
        self.spline_x = np.linspace(knots_x[0], knots_x[-1], n_points)
        self.spline_y = splev(self.spline_x, tck)

def _obj_beam_fit(params, beams, x_nodes):
    """
    params = [2.953e-03, 1.156e+00, 1.297e+01, 9.747e-01 ,  9.970e-01 ,  8.509e-01, 1.076e+00 ,  1.487e+00 ,  7.864e-01,   1.072e+00 ,  1.015e+00]
    """
    import time
    from scipy.interpolate.fitpack import splev, splrep
    import pysynphot as S
    
    ### Spline continuum + gaussian line
    l0 = 6563.*(1+params[1])
    if (l0 < 1.12e4) | (l0 > 1.63e4):
        return -np.inf
        
    line = S.GaussianSource(params[2], l0, 10)
    tck = splrep(x_nodes, params[3:], k=3, s=0)
    xcon = np.arange(0.9e4,1.8e4,0.01e4)
    ycon = splev(xcon, tck, der=0, ext=0)
    spec = S.ArraySpectrum(xcon, ycon, fluxunits='flam', keepneg=True)+line
    
    lnprob = 0
    for key in beams.keys():
        beam = beams[key]
        modelf = beam.compute_model(beam.clip_thumb, xspec=spec.wave, yspec=spec.flux, in_place=False)
        lnprob += -0.5*np.sum(((beam.cutout_scif-params[0]-modelf)**2/beam.cutout_varf)[beam.cutout_maskf])
    
    if ~np.isfinite(lnprob):
        lnprob = -np.inf
        
    print params, lnprob
    #time.sleep(0.2)
    
    return lnprob

def _loss(params, beams, x_nodes):
    import time
    from scipy.interpolate.fitpack import splev, splrep
    import pysynphot as S
    
    ### Spline continuum + gaussian line
    # l0 = 6563.*(1+params[1])
    # if (l0 < 1.12e4) | (l0 > 1.63e4):
    #     l0 = 1.e4
    #     
    # line = S.GaussianSource(params[2], l0, 10)
    tck = splrep(x_nodes, params, k=3, s=0)
    xcon = np.arange(0.9e4,1.8e4,0.01e4)
    ycon = splev(xcon, tck, der=0, ext=0)
    spec = S.ArraySpectrum(xcon, ycon, fluxunits='flam', keepneg=True)#+line
    
    lnprob = 0
    dof = 0
    for key in beams.keys():
        beam = beams[key]
        modelf = beam.compute_model(beam.clip_thumb, xspec=spec.wave, yspec=spec.flux, in_place=False)
        lnprob += np.sum(((beam.cutout_scif-modelf)**2/beam.cutout_varf)[beam.cutout_maskf])
        dof += beam.cutout_maskf.sum()
        
    print params, lnprob, lnprob/(dof-len(params))
    #time.sleep(0.2)
    
    return lnprob

 def check_4(self,f=f1,per=0,s=0,a=0,b=2*pi,N=20,xb=None,xe=None,
           ia=0,ib=2*pi,dx=0.2*pi):
     if xb is None: xb=a
     if xe is None: xe=b
     x=a+(b-a)*arange(N+1,dtype=float)/float(N)    # nodes
     x1=a+(b-a)*arange(1,N,dtype=float)/float(N-1) # middle points of the nodes
     v,v1=f(x),f(x1)
     nk=[]
     put(" u = %s   N = %d"%(repr(round(dx,3)),N))
     put("  k  :  [x(u), %s(x(u))]  Error of splprep  Error of splrep "%(f(0,None)))
     for k in range(1,6):
         tckp,u=splprep([x,v],s=s,per=per,k=k,nest=-1)
         tck=splrep(x,v,s=s,per=per,k=k)
         uv=splev(dx,tckp)
         err1 = abs(uv[1]-f(uv[0]))
         err2 = abs(splev(uv[0],tck)-f(uv[0]))
         assert_(err1 < 1e-2)
         assert_(err2 < 1e-2)
         put("  %d  :  %s    %.1e           %.1e"%\
               (k,repr([round(z,3) for z in uv]),
                err1,
                err2))
     put("Derivatives of parametric cubic spline at u (first function):")
     k=3
     tckp,u=splprep([x,v],s=s,per=per,k=k,nest=-1)
     for d in range(1,k+1):
         uv=splev(dx,tckp,d)
         put(" %s "%(repr(uv[0])))

 def make_diff_func(self) :
     ''' everytime parameters are changed, the diffusion 
         function must be REMADE !! 
     '''
     if self.difftype == 'const':
         self.diff  = lambda y: zeros(len(y), float) + self.param[0]
         self.diff_u= lambda y: zeros(len(y), float) + 0.             
     elif self.difftype == 'power':
         self.diff  = lambda y: y**self.param[0] \
                 * (self.param[1]+self.param[2]*y+self.param[3]*y**2) 
         self.diff_u= lambda y: self.param[0]*y**(self.param[0]-1.) \
                 * (self.param[1]+self.param[2]*y+self.param[3]*y**2) + \
                 y**self.param[0] * (self.param[2]+2.*self.param[3]*y)
     elif self.difftype == 'bspline':
         #create an interp object
         from scipy.interpolate.fitpack import splrep,splev
         #paramuval should be list of u values
         #param should be list of D(u) values
         #create the data of the cubic bspline:
         # no smoother so s = 0
         self.splinedata = splrep(self.paramuval, self.param,
                 xb = None, xe = None, s=0,
                 k = 3, full_output = 0, quiet = 1)
         self.diff   = lambda y: splev(y, self.splinedata, der = 0)
         self.diff_u = lambda y: splev(y, self.splinedata, der = 1)
     elif self.difftype == 'bspline1storder':
         #create an interp object
         from scipy.interpolate.fitpack import splrep,splev
         #paramuval should be list of u values
         #param should be list of D(u) values
         #create the data of the linear bspline:
         # no smoother so s = 0
         self.splinedata = splrep(self.paramuval, self.param,
                 xb = None, xe = None, s=0,
                 k = 1, full_output = 0, quiet = 1)
         self.diff   = lambda y: splev(y, self.splinedata, der = 0)
         self.diff_u = lambda y: splev(y, self.splinedata, der = 1)
     elif self.difftype == '1D2pint':
         #interpolation with 2 points (=piecewise linear)
         self.diff   = GridUtils.GridFunc1D([self.paramuval],self.param)
         self.diff_u = None
     else:
         print ('Unknown diffusion type given', self.difftype)
         sys.exit()
     #value of diffusion is again in agreement with par
     self.modified = False

    def __init__(self):
        # non-uniform grid, just to make it sure
        x = np.linspace(0, 1, 100)**3
        y = np.sin(20 * x)
        self.spl = splrep(x, y)

        # double check that knots are non-uniform
        assert_(np.diff(self.spl[0]).ptp() > 0)

 def test_1d_shape(self):
     x = [1,2,3,4,5]
     y = [4,5,6,7,8]
     tck = splrep(x, y)
     z = splev([1], tck)
     assert_equal(z.shape, (1,))
     z = splev(1, tck)
     assert_equal(z.shape, ())

 def test_2d_shape(self):
     x = [1, 2, 3, 4, 5]
     y = [4, 5, 6, 7, 8]
     tck = splrep(x, y)
     t = np.array([[1.0, 1.5, 2.0, 2.5], [3.0, 3.5, 4.0, 4.5]])
     z = splev(t, tck)
     z0 = splev(t[0], tck)
     z1 = splev(t[1], tck)
     assert_equal(z, np.row_stack((z0, z1)))

def compute_colors(N):
    xref = np.linspace(0, 1, CMRref.shape[0])
    x = np.linspace(0, 1, N)
    cmap = np.zeros((N, 3))

    for i in range(3):
        tck = splrep(xref, CMRref[:, i], s=0)  # cubic spline (default) without smoothing
        cmap[:, i] = splev(x, tck)

    # Limit to range [0,1]
    cmap -= np.min(cmap)
    cmap /= np.max(cmap)

    return cmap

    def test_extrapolation_modes(self):
        # test extrapolation modes
        #    * if ext=0, return the extrapolated value.
        #    * if ext=1, return 0
        #    * if ext=2, raise a ValueError
        #    * if ext=3, return the boundary value.
        x = [1,2,3]
        y = [0,2,4]
        tck = splrep(x, y, k=1)

        rstl = [[-2, 6], [0, 0], None, [0, 4]]
        for ext in (0, 1, 3):
            assert_array_almost_equal(splev([0, 4], tck, ext=ext), rstl[ext])

        assert_raises(ValueError, splev, [0, 4], tck, ext=2)

    def check_1(self, f=f1, per=0, s=0, a=0, b=2 * pi, N=20, at=0, xb=None, xe=None):
        if xb is None:
            xb = a
        if xe is None:
            xe = b
        x = a + (b - a) * arange(N + 1, dtype=float) / float(N)  # nodes
        x1 = a + (b - a) * arange(1, N, dtype=float) / float(N - 1)  # middle points of the nodes
        v, v1 = f(x), f(x1)
        nk = []

        def err_est(k, d):
            # Assume f has all derivatives < 1
            h = 1.0 / float(N)
            tol = 5 * h ** (0.75 * (k - d))
            if s > 0:
                tol += 1e5 * s
            return tol

        for k in range(1, 6):
            tck = splrep(x, v, s=s, per=per, k=k, xe=xe)
            if at:
                t = tck[0][k:-k]
            else:
                t = x1
            nd = []
            for d in range(k + 1):
                tol = err_est(k, d)
                err = norm2(f(t, d) - splev(t, tck, d)) / norm2(f(t, d))
                assert_(err < tol, (k, d, err, tol))
                nd.append((err, tol))
            nk.append(nd)
        put(
            "\nf = %s  s=S_k(x;t,c)  x in [%s, %s] > [%s, %s]"
            % (f(None), repr(round(xb, 3)), repr(round(xe, 3)), repr(round(a, 3)), repr(round(b, 3)))
        )
        if at:
            str = "at knots"
        else:
            str = "at the middle of nodes"
        put(" per=%d s=%s Evaluation %s" % (per, repr(s), str))
        put(" k :  |f-s|^2  |f'-s'| |f''-.. |f'''-. |f''''- |f'''''")
        k = 1
        for l in nk:
            put(" %d : " % k)
            for r in l:
                put(" %.1e  %.1e" % r)
            put("\n")
            k = k + 1

 def check_3(self,f=f1,per=0,s=0,a=0,b=2*pi,N=20,xb=None,xe=None,
           ia=0,ib=2*pi,dx=0.2*pi):
     if xb is None: xb=a
     if xe is None: xe=b
     x=a+(b-a)*arange(N+1,dtype=float)/float(N)    # nodes
     v=f(x)
     nk=[]
     put("  k  :     Roots of s(x) approx %s  x in [%s,%s]:"%\
           (f(None),repr(round(a,3)),repr(round(b,3))))
     for k in range(1,6):
         tck=splrep(x,v,s=s,per=per,k=k,xe=xe)
         roots = sproot(tck)
         if k == 3:
             assert_allclose(roots, pi*array([1, 2, 3, 4]),
                             rtol=1e-3)
         put('  %d  : %s'%(k,repr(roots.tolist())))

def interp_masked1d(marr, kind='linear'):
    """

    Interpolates masked values in an array according to the given method.

    Parameters
    ----------
    marr : MaskedArray
        Array to fill
    kind : {'constant', 'linear', 'cubic', quintic'}, optional
        Type of interpolation

    """
    if np.ndim(marr) > 1:
        raise ValueError("array must be 1 dimensional!")
    #
    marr = marray(marr, copy=True)
    if getmask(marr) is nomask:
        return marr
    #
    unmaskedIndices = (~marr._mask).nonzero()[0]
    if unmaskedIndices.size < 2:
        return marr
    #
    kind = kind.lower()
    if kind == 'constant':
        return forward_fill(marr)
    try:
        k = {'linear' : 1,
             'cubic' : 3,
             'quintic' : 5}[kind.lower()]
    except KeyError:
        raise ValueError("Unsupported interpolation type.")

    first_unmasked, last_unmasked = flatnotmasked_edges(marr)

    vals = marr.data[unmaskedIndices]

    from scipy.interpolate import fitpack
    tck = fitpack.splrep(unmaskedIndices, vals, k=k)

    maskedIndices = marr._mask.nonzero()[0]
    interpIndices = maskedIndices[(maskedIndices > first_unmasked) & \
                                  (maskedIndices < last_unmasked)]
    marr[interpIndices] = fitpack.splev(interpIndices, tck).astype(marr.dtype)
    return marr

    def check_2(self, f=f1, per=0, s=0, a=0, b=2 * pi, N=20, xb=None, xe=None, ia=0, ib=2 * pi, dx=0.2 * pi):
        if xb is None:
            xb = a
        if xe is None:
            xe = b
        x = a + (b - a) * arange(N + 1, dtype=float) / float(N)  # nodes
        v = f(x)

        def err_est(k, d):
            # Assume f has all derivatives < 1
            h = 1.0 / float(N)
            tol = 5 * h ** (0.75 * (k - d))
            if s > 0:
                tol += 1e5 * s
            return tol

        nk = []
        for k in range(1, 6):
            tck = splrep(x, v, s=s, per=per, k=k, xe=xe)
            nk.append([splint(ia, ib, tck), spalde(dx, tck)])
        put(
            "\nf = %s  s=S_k(x;t,c)  x in [%s, %s] > [%s, %s]"
            % (f(None), repr(round(xb, 3)), repr(round(xe, 3)), repr(round(a, 3)), repr(round(b, 3)))
        )
        put(
            " per=%d s=%s N=%d [a, b] = [%s, %s]  dx=%s"
            % (per, repr(s), N, repr(round(ia, 3)), repr(round(ib, 3)), repr(round(dx, 3)))
        )
        put(" k :  int(s,[a,b]) Int.Error   Rel. error of s^(d)(dx) d = 0, .., k")
        k = 1
        for r in nk:
            if r[0] < 0:
                sr = "-"
            else:
                sr = " "
            put(" %d   %s%.8f   %.1e " % (k, sr, abs(r[0]), abs(r[0] - (f(ib, -1) - f(ia, -1)))))
            d = 0
            for dr in r[1]:
                err = abs(1 - dr / f(dx, d))
                tol = err_est(k, d)
                assert_(err < tol, (k, d))
                put(" %.1e %.1e" % (err, tol))
                d = d + 1
            put("\n")
            k = k + 1

def interp_masked1d(marr, kind='linear'):
    """interp_masked1d(marr, king='linear')

Interpolates masked values in marr according to method kind.
kind must be one of 'constant', 'linear', 'cubic', quintic'
"""
    if numeric.ndim(marr) > 1: 
        raise ValueError("array must be 1 dimensional!")
    #
    marr = marray(marr, copy=True)
    if getmask(marr) is nomask: 
        return marr
    #
    unmaskedIndices = (~marr._mask).nonzero()[0]
    if unmaskedIndices.size < 2: 
        return marr
    #    
    kind = kind.lower()
    if kind == 'constant': 
        return forward_fill(marr)
    try:
        k = {'linear' : 1,
             'cubic' : 3,
             'quintic' : 5}[kind.lower()]
    except KeyError:
        raise ValueError("Unsupported interpolation type.")
    
    first_unmasked, last_unmasked = flatnotmasked_edges(marr)
    
    vals = marr.data[unmaskedIndices]
    
    tck = fitpack.splrep(unmaskedIndices, vals, k=k)
    
    maskedIndices = marr._mask.nonzero()[0]
    interpIndices = maskedIndices[(maskedIndices > first_unmasked) & \
                                  (maskedIndices < last_unmasked)]
    marr[interpIndices] = fitpack.splev(interpIndices, tck).astype(marr.dtype)
    return marr

def test_splev_der_k():
    # regression test for gh-2188: splev(x, tck, der=k) gives garbage or crashes
    # for x outside of knot range

    # test case from gh-2188
    tck = (np.array([0., 0., 2.5, 2.5]),
           np.array([-1.56679978, 2.43995873, 0., 0.]),
           1)
    t, c, k = tck
    x = np.array([-3, 0, 2.5, 3])

    # an explicit form of the linear spline
    assert_allclose(splev(x, tck), c[0] + (c[1] - c[0]) * x/t[2])
    assert_allclose(splev(x, tck, 1), (c[1]-c[0]) / t[2])

    # now check a random spline vs splder
    np.random.seed(1234)
    x = np.sort(np.random.random(30))
    y = np.random.random(30)
    t, c, k = splrep(x, y)

    x = [t[0] - 1., t[-1] + 1.]
    tck2 = splder((t, c, k), k)
    assert_allclose(splev(x, (t, c, k), k), splev(x, tck2))

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate.fitpack import splrep, splev

from util.point import Point
from walking.helpers import bezier_curves
from walking.helpers.bezier_curves import LineSegment


x = (0, 2.0, 6.0, 8.0)
y = (-47.7, -44.7, -44.7, -47.7)

tck = splrep(x, y)

x2 = np.linspace(0, 8)
y2 = splev(x2, tck)

points = [Point((x_, 0, y_)) for (x_, y_) in zip(x, y)]
l1 = LineSegment(points[0], points[1])
l2 = LineSegment(points[1], points[2])
l3 = LineSegment(points[2], points[3])

yy2 = [bezier_curves.get_interpolated_position(x2_/8, l1, l2, l3)[2] for x2_ in x2]

plt.plot(x, y, 'o', x, y, '--', x2, y2, x2, yy2)

plt.show()

diodes=[]
for rad in ['0','1']:
	for diode in ['0','1']:
		diodes.append(str(horn)+rad+diode)

for diode in diodes:
	#first do 91-952, read response file then go thru ODs
	fl=glob(adc_dir+'/'+'ADCDX12_%s*_952_*.pic' % diode)
	if len(fl)>0:
		respfile=fl[0]
		if os.path.exists(respfile):
			resp = read_LFI_response(respfile)
			horn = resp.keys['horn']
			rad  = resp.keys['radiometer']
			det  = resp.keys['detector'] 
			sky_spline = fit.splrep(resp.sky_volt_out,resp.sky_volt_in,s=0.0)
			ref_spline = fit.splrep(resp.load_volt_out,resp.load_volt_in,s=0.0)
			for od in range(91,953):
				rawfile=rawdata_dir+'od%s_rca%s.fits' % (od,diode)
				corrfile=adc_corr_data_dir+'od%s_rca%s.fits' % (od,diode)
				if not(os.path.exists(corrfile)):
					if os.path.exists(rawfile):
						hdulist = fits.open(rawfile)
						data = hdulist[1].data
						sky_volt  = data.field('SKY')
						load_volt = data.field('REF')						
						sky_cor = fit.splev(sky_volt,sky_spline)
						load_cor = fit.splev(load_volt,ref_spline)
						data.field('SKY')[:] = sky_cor 
						data.field('REF')[:] = load_cor
						hdulist[1].header.update('hierarch instrument','LFI_ADC_corr')

def setup_fit():
    
    snlim = 0
    snlim = 1./np.sqrt(len(beams))
    
    Rmax = 5
    
    for key in FLT.keys():
        beam = beams[key]
        beam.cutout_sci = beam.get_cutout(FLT[key].im['SCI'].data)*1
        beam.cutout_dq = beam.get_cutout(FLT[key].im['DQ'].data)*1
        beam.cutout_err = beam.get_cutout(FLT[key].im['ERR'].data)*1
        
        beam.cutout_mask = (beam.cutout_dq-(beam.cutout_dq & 512) == 0) & (beam.cutout_err > 0) & (beam.cutout_err < 10)
        
        sh = beam.cutout_sci.shape
        yp, xp = np.indices(beam.thumb.shape)
        r = np.sqrt((xp-sh[0]/2)**2+(yp-sh[0]/2)**2)
        beam.norm = beam.thumb[r <= Rmax].sum()/1.e-17
        beam.clip_thumb = beam.thumb*(r <= Rmax)
        beam.compute_model(beam.clip_thumb)
        
        sn = beam.model/beam.cutout_err
        #beam.mask &= (sn > 0.5) & (beam.err > 0)
        beam.cutout_mask &= (sn > snlim) & (beam.cutout_err > 0)
        
        beam.cutout_sci[~beam.cutout_mask] = 0
        beam.cutout_err[~beam.cutout_mask] = 0
        beam.cutout_var = beam.cutout_err**2
        
        beam.cutout_scif = beam.cutout_sci.flatten()
        beam.cutout_varf = beam.cutout_var.flatten()
        beam.cutout_maskf = beam.cutout_mask.flatten()
                
        ds9.view((beam.cutout_scif-beam.cutout_modelf).reshape(beam.cutout_sci.shape)*beam.cutout_mask)
        
    
    ##
    x_nodes = np.arange(1.0e4,1.71e4,0.1e4)
    x_nodes = np.arange(1.0e4,1.71e4,0.03e4)
    x_nodes = np.arange(1.05e4,1.71e4,0.075e4)
    
    init = np.append([0, 1.148, 500], np.ones(len(x_nodes)))
    init = np.append([0, 1.2078, 500], np.ones(len(x_nodes)))
    step_sig = np.append([0.01,0.1,50], np.ones(len(x_nodes))*0.1)    

    init = np.append([0, 1.0, 500], np.ones(len(x_nodes)))
    step_sig = np.append([0.0,0.2,150], np.ones(len(x_nodes))*0.1)    

    ### don't fit redshift
    #init = np.append([0, 1.0, 0], np.ones(len(x_nodes)))
    #step_sig = np.append([0.0,0.,0], np.ones(len(x_nodes))*0.2)    
    
    obj_fun = _obj_beam_fit
    obj_args = [beams, x_nodes]
        
    ndim, nwalkers = len(init), len(init)*2
    p0 = [(init+np.random.normal(size=ndim)*step_sig) 
          for i in xrange(nwalkers)]
    
    NTHREADS, NSTEP = 2, 5
    sampler = emcee.EnsembleSampler(nwalkers, ndim, obj_fun, args = obj_args, 
                                    threads=NTHREADS)
    #
    t0 = time.time()
    result = sampler.run_mcmc(p0, NSTEP)
    t1 = time.time()
    print 'Sampler: %.1f s' %(t1-t0)
    
    param_names = ['bg', 'z', 'Ha']
    param_names.extend(['spl%d' %i for i in range(len(init)-3)])
    
    import unicorn.interlace_fit
    chain = unicorn.interlace_fit.emceeChain(chain=sampler.chain, param_names=param_names)
    
    #obj_fun(chain.map, beams)
    #obj_fun(init, beams, x_nodes)
    params = chain.map*1.
    #params = init
    
    ### Show spectra
    from scipy.interpolate.fitpack import splev, splrep

    # NDRAW=100
    # draw = chain.draw_random(NDRAW)
    # for i in range(NDRAW):
    #     line = S.GaussianSource(draw[i,2], 6563.*(1+draw[i,1]), 10)
    #     tck = splrep(x_nodes, draw[i,3:], k=3, s=0)
    #     xcon = np.arange(0.9e4,1.8e4,0.01e4)
    #     ycon = splev(xcon, tck, der=0, ext=0)
    #     spec = S.ArraySpectrum(xcon, ycon, fluxunits='flam', keepneg=True)+line
    #     plt.plot(spec.wave, spec.flux, alpha=0.1, color='red')
    #     
    line = S.GaussianSource(params[2], 6563.*(1+params[1]), 10)
    tck = splrep(x_nodes, params[3:], k=3, s=0)
    xcon = np.arange(0.9e4,1.8e4,0.01e4)
    ycon = splev(xcon, tck, der=0, ext=0)
    pspec = S.ArraySpectrum(xcon, ycon, fluxunits='flam', keepneg=True)+line

    for key in beams.keys():
#.........这里部分代码省略.........

    yi = np.zeros(np.shape(xi))
    yi_hd = np.zeros(np.shape(xi_hd))
    yi[:] = f(xi[:])
    yi_hd = f(xi_hd[:])

    #xi = np.linspace(-2, 2, 101)
    #yi = 10 * xi / (1 + 100 * np.power(xi, 2))
    data = np.c_[xi, yi]

    p = poly.NewtonPolynomial(points=data)
    d4p = deriv.dr_dxr(p, r=4)
    norm = norm.Norm(d4p)
    T, eps = cut.cutab(norm, xi, eps, r)

    tck = fitpack.splrep(xi, yi, t=T[1:-1])
    fit = fitpack.splev(xi, tck)
    error = np.abs(fit - yi)
    print '*' * 10 + " FIRST ERROR RATIO " + '*' * 10 + '\n\n'
    print "error / tol = %f" % (error.max() / tol)

    while error.max() > tol:
        T, eps = cut.cutab(norm, xi, eps, r)
        print "\n\nRun %d" % run_count
        print "eps = %e in outer loop" % eps
        eps = eps / 2
        tck = fitpack.splrep(xi, yi, t=T[1:-1])
        fit = fitpack.splev(xi, tck)
        fit_hd = fitpack.splev(xi_hd, tck)
        fT = np.zeros(np.shape(T))
        fT[:] = f(T[:])