from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

def mspline(x, k, i, T):
    """Evaluate the kth order m-spline at x.                                  
                                                                              
    Arguments:                                                                
    ----------                                                                
    x : the point at which the spline is evaluated                            
    k : the order of the spline (k>=1)                                        
    i : the interval on which the spline is defined (i<len(T)-k)              
    T : the knots                                                             
                                                                              
    Returns:                                                                  
    --------                                                                  
    The value of the spline at x                                              
                                                                              
    Note: Use mspline_knots(...) to create T.                                 
    """

    if x<T[i]: return 0.
    if x>=T[i+k]: return 0.

    if k==1: return 1/(T[i+1]-T[i])

    d1 = x-T[i]
    d2 = T[i+k]-x

    v = d1*mspline(x, k-1, i, T) + d2*mspline(x, k-1, i+1, T)
    v *= k
    v /= (k-1)*(T[i+k]-T[i])

    return v

def ispline(x, k, i, T):
    """Evaluate the kth order i-spline at x.                                  
                                                                              
    Arguments:                                                                
    ----------                                                                
    x : the point at which the spline is evaluated                            
    k : the order of the spline                                               
    i : the interval on which the spline is defined (i<len(T)-k-1)            
    T : the knots                                                             
                                                                              
    Returns:                                                                  
    --------                                                                  
    The value of the spline at x                                              
                                                                              
    Note: Use ispline_knots(...) to create T.                                 
    """

    j = np.flatnonzero(T<=x)[-1]
    if j<i: return 0.
    if j-k+1>i: return 1.

    v = 0.
    for m in range(i,j+1):
        v += mspline(x, k+1, m, T)*(T[m+k+1]-T[m])/(k+1)

    return v

def mspline_knots(xmin, xmax, k, m):
    """Create an order list of points in the inteval (xmin,xmax)              
    suitable for use as knots in m-splines and i-splines.                     
                                                                              
    Arguments:                                                                
    ----------                                                                
    xmin : the lower bound of the interval                                    
    xmax : the upper bount of the interval                                    
    k : the order of the m-spline for which the knots will be used            
    m : the number of intervals into which (xmin,xmax) will be partitioned    
    """

    T = [xmin]*(k-1)
    h = (xmax-xmin)/m
    T.extend(np.arange(xmin, xmax, h))
    T.extend([xmax]*k)

    return T

def ispline_knots(xmin, xmax, k, m):
    """Create an order list of points in the inteval (xmin,xmax)              
    suitable for use as knots in m-splines and i-splines.                     
                                                                              
    Arguments:                                                                
    ----------                                                                
    xmin : the lower bound of the interval                                    
    xmax : the upper bount of the interval                                    
    k : the order of the m-spline for which the knots will be used            
    m : the number of intervals into which (xmin,xmax) will be partitioned    
    """

    return mspline_knots(xmin, xmax, k+1, m)

def test_ispline():

    plt.clf()
    plt.axes([0.1,0.1,0.7,0.8])
    k = 3
    T = ispline_knots(0, 1, k, 5)

    A = []
    for i in range(1,len(T)-k-1):
        F = []
        for x in np.arange(0,1.01,0.01):
            F.append([x,ispline(x, k, i, T)])
        F = np.array(F)
        a = plt.plot(F[:,0], F[:,1], '-')
        A.append(a)
        if i==0: print F
        plt.hold(True)
    b = plt.figlegend(A, ("i=1", "i=2", "i=3", "i=4", "i=5", "i=6", "i=7"),\
                   'center right')
    b.draw_frame(False)
    plt.xlabel(r"$x$")
    plt.ylabel(r"$M_i(x