def fitEdge(og,e,VERBOSE=False):
    path = og[e[0]][e[1]]['wpath']
    #pstart=og.node[e[0]]['wcrd']
    #pend = og.node[e[1]]['wcrd']
    #path.extend([pend])
    #p = [pstart]
    #p.extend(path)
    ae = np.array(path)
    
    if( ae.shape[0] > 3 ): # there are not enough points to fit with
        if( VERBOSE ): print("refitting edge with %d points"%len(path))

        s = ae.shape[0]/2.0

        fit2 = splprep(ae.transpose(),task=0,full_output =1, s=s)[0]
        u = np.array(list(range(ae.shape[0]+1))).\
                astype(np.float64)/(ae.shape[0])
        # location of spline points
        og[e[0]][e[1]]['d0'] = splev(u,fit2[0],der=0)
        # first derivative (tangent) of spline
        og[e[0]][e[1]]['d1'] =np.array( splev(u,fit2[0],der=1))
        # second derivative (curvature) of spline
        og[e[0]][e[1]]['d2'] = np.array(splev(u,fit2[0],der=2))
    else:
        if( VERBOSE ): print("no or insufficient points to fit")
        og[e[0]][e[1]]['d0'] = None
        og[e[0]][e[1]]['d1'] = None
        og[e[0]][e[1]]['d2'] = None

 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 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 test_splantider_vs_splint(self):
        # Check antiderivative vs. FITPACK
        spl2 = splantider(self.spl)

        # no extrapolation, splint assumes function is zero outside
        # range
        xx = np.linspace(0, 1, 20)

        for x1 in xx:
            for x2 in xx:
                y1 = splint(x1, x2, self.spl)
                y2 = splev(x2, spl2) - splev(x1, spl2)
                assert_allclose(y1, y2)

    def test_splder_vs_splev(self):
        # Check derivative vs. FITPACK

        for n in range(3+1):
            # Also extrapolation!
            xx = np.linspace(-1, 2, 2000)
            if n == 3:
                # ... except that FITPACK extrapolates strangely for
                # order 0, so let's not check that.
                xx = xx[(xx >= 0) & (xx <= 1)]

            dy = splev(xx, self.spl, n)
            spl2 = splder(self.spl, n)
            dy2 = splev(xx, spl2)
            assert_allclose(dy, dy2)

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 _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 __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 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 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 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))

 def fitPath(self,maxiter=40, s=20000):
     """fits a least squares spline through the path
     find descriptions of maxiter and s from scipy documentation"""
     try:
         pth = path.getCrds()
         pp = [pth[:,0],pth[:,1],pth[:,2]]
         if( len(pp[0]) <= 3 ): # there are not enough points to fit with
             self.splinePoints = pp
             self.splineTangents = None
             self.splineCurvature = None
             return 
         else:
             s = len(pp[0])/2.0
             cont = 1
             j=0
             while( j < maxiter ):
                 fit2 = splprep(pp,task=0,full_output =1, s=s)[0]
                 u = na.array(range(pp[0].shape[0]+1)).\
                         astype(na.float64)/(pp[0].shape[0])
                 # location of spline points
                 self.splinePoints = splev(u,fit2[0],der=0)
                 # first derivative (tangent) of spline
                 valsd = splev(u,fit2[0],der=1)
                 self.splineTangents = na.array(valsd).astype(na.float64)
                 # second derivative (curvature) of spline
                 valsc = splev(u,fit2[0],der=2)
                 self.splineCurvature = na.array(valsc).astype(na.float64)
                 success = 1 # this doesn't make much sense
                 if( success ): # how should I be determining success?
                     break
                 else:
                     s = s*0.9
                 j += 1
             return True
     except Exception, error:
         print "failed in fitPath()",error

    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 spli_basex(p, x ,knots=None , deg = 3 , order = 0 ):
    """Vandermonde type matrix for splines.

    Returns a matrix whose columns are the values of the b-splines of deg
    `deg` associated with the knot sequence `knots` evaluated at the points
    `x`.

    Parameters
    ----------
    p : array_like
        Parameter array containing:
         - the number of knots
         - the lower bound of the approximation
         - the upper bound of the approximation
    x : array_like
        Points at which to evaluate the b-splines.
    deg : int
        Degree of the splines.
        Default: cubic splines
    knots : array_like
        List of knots. The convention here is that the interior knots have
        been extended at both ends by ``deg + 1`` extra knots - see augbreaks.
        If not given the default is equidistant grid
    order : int
        Evaluate the derivative of the spline
    Returns
    -------
    vander : ndarray
        Vandermonde like matrix of shape (m,n), where ``m = len(x)`` and
        ``n = len(augbreaks) - deg - 1``

    Notes
    -----
    The knots exending the interior points are usually taken to be the same
    as the endpoints of the interval on which the spline will be evaluated.

    """
    n , a , b = p[0] , p[1] , p[2]
    if knots is None:
        knots = np.linspace(a , b , n + 1 - deg)
    augbreaks = np.concatenate(( a * np.ones((deg)),knots, b * np.ones((deg))))
    m = len(augbreaks) - deg - 1
    v = np.empty((m, len(x)))
    d = np.eye(m, len(augbreaks))
    for i in range(m):
        v[i] = spl.splev(x, (augbreaks, d[i], deg),order)
    return v.T

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_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)
     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)
         if k == 3:
             roots = sproot(tck)
             assert_allclose(splev(roots, tck), 0, atol=1e-10, rtol=1e-10)
             assert_allclose(roots, pi * array([1, 2, 3, 4]), rtol=1e-3)
             put("  %d  : %s" % (k, repr(roots.tolist())))
         else:
             assert_raises(ValueError, sproot, tck)

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

			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')
						hdulist[1].header.update('hierarch ADC_correction','True')
						print 'ADC correcting %s -> %s' % (rawfile,corrfile)
						hdulist.writeto(corrfile)
	if len(fl)==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 os.path.exists(rawfile):
				if not(os.path.exists(corrfile)):
					shutil.copy(rawfile,corrfile)