!date
import matplotlib.pyplot as plt, seaborn as sns
%matplotlib inline

import pymc.gp, sklearn

import numpy as np
from sklearn import gaussian_process

def f(x):
    return x * np.sin(x)
X = np.atleast_2d([1., 3., 5., 6., 7., 8.]).T
y = f(X).ravel()

gp = gaussian_process.GaussianProcess(theta0=1e-2, thetaL=1e-4, thetaU=1e-1, nugget=.00001)
gp.fit(X, y)  

x = np.atleast_2d(np.linspace(0, 10, 1000)).T
y_pred, sigma2_pred = gp.predict(x, eval_MSE=True)
sigma_pred = np.sqrt(sigma2_pred)

plt.plot(X.ravel(), y, 's')
plt.plot(x, y_pred)

gp = gaussian_process.GaussianProcess(corr=sklearn.gaussian_process.correlation_models.generalized_exponential,
    theta0=[1e-2,10], thetaL=[1e-4,10], thetaU=[1e-1,10], nugget=1)
gp.fit(X, y)

x = np.atleast_2d(np.linspace(0, 10, 1000)).T
y_pred, sigma2_pred = gp.predict(x, eval_MSE=True)
sigma_pred = np.sqrt(sigma2_pred)

plt.plot(X.ravel(), y, 's')
plt.plot(x, y_pred)

pymc.gp.cov_funs.matern.euclidean(X, np.array([0]), diff_degree=1.4, amp=0.4, scale=10)

def iso_matern(theta, d):
    """
    Matern correlation model.

        theta, d --> r(diff_degree=theta[0], amp=theta[1], scale=theta[2]; d)

    Parameters
    ----------
    theta : array_like
        An array with shape 3 (isotropic) giving the
        autocorrelation parameters (diff_degree, amp, scale)
    d : array_like
        An array with shape (n_eval, n_features) giving the componentwise
        distances between locations x and x' at which the correlation model
        should be evaluated.
    Returns
    -------
    r : array_like
        An array with shape (n_eval, ) with the values of the autocorrelation
        model.
    """
    print theta
    #return sklearn.gaussian_process.correlation_models.generalized_exponential(theta[:2], d)
    theta = np.asarray(theta, dtype=np.float).reshape(3)
    #theta[0] = 2
    #theta[1] = .4
    d = np.asarray(d, dtype=np.float)

    assert theta.shape == (3,), 'expected theta.shape == (3,),  with diff_degree=theta[0], amp=theta[1], scale=theta[2]'

    r = pymc.gp.cov_funs.matern.euclidean(
        d, np.array([0]),
        diff_degree=theta[0], amp=theta[1], scale=theta[2])
    return np.array(r).reshape(d.shape[0],)
iso_matern([1.4,.4,10], X).shape

gp = gaussian_process.GaussianProcess(corr=iso_matern,
    theta0=[1.5,1,10], nugget=.001)
gp.fit(X, y)

x = np.atleast_2d(np.linspace(0, 10, 1000)).T
y_pred, sigma2_pred = gp.predict(x, eval_MSE=True)
sigma_pred = np.sqrt(sigma2_pred)

plt.plot(X.ravel(), y, 's')

plt.plot(x, y_pred)

gp = gaussian_process.GaussianProcess(corr=iso_matern,
    theta0=[2.5,.5,1], thetaL=[1., .1, 1], thetaU=[10., 1., 1], nugget=y/10.)
gp.fit(X, y)

x = np.atleast_2d(np.linspace(0, 10, 1000)).T
y_pred, sigma2_pred = gp.predict(x, eval_MSE=True)
sigma_pred = np.sqrt(sigma2_pred)

plt.plot(X.ravel(), y, 's')

plt.plot(x, y_pred)