[SciPy-User] Crossing of Splines

[denis]

On Jun 8, 8:22 pm, Marco <gae... at gmail.com> wrote:
> Hi all!
>
> I have 2 different datasets which I fit using interpolate.splrep().
>
> I am interested in finding the point where the splines cross: as of
> now I use interpolate.splev() to evaluate each spline and then look
> for a zero in the difference of the evaluated splines.


Following Anne's idea, if you extend class UnivariateSpline to use
given knots,
you could then subtract the whole splines / all their coefficients
and run roots() on the difference, along the lines

    a = myUniSpline( x, ya, s=s, knots= )
    b = myUniSpline( x, yb, s=s, knots= )
    a.minus( b.getcoeffs() )  # spline a - b
    roots = a.roots()

Marco, are you interpolating, s=0, or smoothing, s > 0 ?
If smoothing, knots wobble around with s and y (try the little test
below) --
good, adaptive, for interpolating 1 function, not so good for your
problem..

cheers
  -- denis

    """ scipy UnivariateSpline sensitivity to s """

    from __future__ import division
    import sys
    import numpy as np
    from scipy.interpolate import UnivariateSpline  # $scipy/
interpolate/fitpack2.py
    import pylab as pl

    np.set_printoptions( 2, threshold=10, edgeitems=3, suppress=True )

    N = 100
    H = 5
    cycle = 4
    plot = 0
    exec "\n".join( sys.argv[1:] )  # N= ...

    x = np.arange(N+1)
    xup = np.arange( 0, N + 1e-6, 1/H )
    if cycle == 0:
        y = np.zeros(N+1)
        y[N//2] = 1
    else:
        y = np.sin( 2*np.pi * np.arange(N+1) / cycle )

 
#...............................................................................
    title = "UnivariateSpline N=%d H=%d cycle=%.2g" % (N, H, cycle)
    print title
    for s in (0, .1, .5, 10, None ):
        uspline = UnivariateSpline( x, y, s=s )  # s=0 interpolates
        yup = uspline( xup )
        res = uspline.get_residual()  # == s == |y - yup[::H]|**2
        label = "s: %s  res: %.2g" % (s, res)
        print label
        knots = uspline.get_knots()
        print "%d knots: %s" % (len(knots), knots )
        roots = uspline.roots()
        print "%d roots: %s" % (len(roots), roots )
        print ""
        if plot:
            pl.plot( xup, yup, label=label )
            if cycle == 0:
                pl.xlim( N//2 - 10, N//2 + 10 )

    if plot:
        pl.title(title)
        pl.legend()
        pl.show()