#!/usr/bin/python

import numpy as N

from scipy.interpolate import splrep,splev
from scipy.optimize import fmin_cobyla,leastsq,fmin
from scipy import io

def func(a):
    sp = splrep(x,a)
    return splev(xdata, sp)
        
def minimize(a):
    res = N.abs(ydata-func(a))
    return res.sum()

class leftturn(object):
    """constraint which forces a positive second derivative"""
    def __init__(self, n):
        self.n = n

    def __call__(self, y, x, xdata, ydata):
        sp = splrep(x,y)
        dy = splev([xdata[self.n]],sp,1)
        return -dy

class smaller(object):
    """constraint which forces the model to be smaller than the data"""
    def __init__(self, n):
        self.n = n

    def __call__(self, y, x, xdata, ydata):
        sp = splrep(x,y)
        ny = splev([xdata[self.n]],sp)
        return ydata[self.n]-ny
    
n = 100
ord_spline = 5

xdata = N.linspace(0,10,n)
ydata = N.exp(-xdata/2.0)+N.random.normal(0,.03,xdata.shape)+0.1*N.sin(xdata*2.0)

x = N.linspace(xdata.min(),xdata.max(),ord_spline)

res = fmin_cobyla(minimize, N.zeros(x.shape), cons=[smaller(i) for i in range(n)], consargs=(x,xdata,ydata))
y = func(res)

import pylab as P

P.plot(xdata,ydata)
P.plot(xdata,y)
P.plot(xdata,ydata-y)
P.plot(x,res,'ro')
P.show()