'''This program computes the derivatives of noisy functions numerically in three different ways.
'''

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage.filters import gaussian_filter1d
from scipy.interpolate import InterpolatedUnivariateSpline

'''data'''
n = 801
x = np.linspace(-2,2,n)
y = -2*x**3 + 4*x

'''random numbers'''
np.random.seed(0)
dy = np.random.normal(0,1,n)/200

'''noise added to data'''
y += dy

''' left figure '''
plt.subplot(131)
plt.plot(x,y,'b-')

'''numerical derivative of noisy function'''
# We compute (y[i+1] - y[i]) / (x[i+1] - x[i]) for i in range(n-1)
v = (y[1:] - y[:-1]) / (x[1:] - x[:-1])
plt.plot(x[:-1], v ,'r-', linewidth=0.5)
plt.axis([-2.1,2.1,-10,10])

''' central figure '''
plt.subplot(132)
plt.plot(x,y,'b')

'''filtering the data'''
u = gaussian_filter1d(y,7, mode='nearest')

'''numerical derivative after filtering the data'''
v = (u[1:] - u[:-1]) / (x[1:] - x[:-1])
plt.plot(x[:-1], v,'r-')
plt.axis([-2.1,2.1,-10,10])

'''right figure'''
plt.subplot(133)
plt.plot(x,y,'b-')

'''smoothing of data by means of splines'''
spline = InterpolatedUnivariateSpline(x, y)
spline.set_smoothing_factor(0.5)
u = spline(x)

'''numerical derivative of splines'''
v = (u[1:] - u[:-1]) / (x[1:] - x[:-1])
plt.plot(x[:-1], v,'r-')
plt.axis([-2.1,2.1,-10,10])

plt.show()