Python fitpack.splrep函数代码示例

    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