# Path: rbf/gproc.py
# def gpiso(phi, eps, var, dim=None):
#     '''
#     Creates an isotropic Gaussian process which has a covariance function that
#     is described by a radial basis function.
# 
#     Parameters
#     ----------
#     phi : str or RBF instance
#         Radial basis function describing the covariance function.
# 
#     eps : float
#         Shape parameter.
# 
#     var : float
#         Variance.
# 
#     dim : int, optional
#         Fixes the spatial dimensions of the domain.
# 
#     Returns
#     -------
#     GaussianProcess
# 
#     Notes
#     -----
#     Not all radial basis functions are positive definite, which means that it
#     is possible to instantiate an invalid `GaussianProcess`. The method
#     `is_positive_definite` provides a necessary but not sufficient test for
#     positive definiteness. Examples of predefined `RBF` instances which are
#     positive definite include: `rbf.basis.se`, `rbf.basis.ga`, `rbf.basis.exp`,
#     `rbf.basis.iq`, `rbf.basis.imq`.
# 
#     '''
#     phi = rbf.basis.get_rbf(phi)
# 
#     def isotropic_covariance(x1, x2, diff1, diff2):
#         diff = diff1 + diff2
#         coeff = var*(-1)**sum(diff2)
#         out = coeff*phi(x1, x2, eps=eps, diff=diff)
#         return out
# 
#     def isotropic_variance(x, diff):
#         coeff = var*(-1)**sum(diff)
#         value = coeff*phi.center_value(eps=eps, diff=2*diff)
#         out = np.full(x.shape[0], value)
#         return out
# 
#     out = GaussianProcess(
#         covariance=isotropic_covariance,
#         variance=isotropic_variance,
#         dim=dim,
#         differentiable=True
#         )
#     return out


import logging
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from rbf.gproc import gpiso
''' 
This script demonstrates how to perform Gaussian process regression on a noisy
data set. It also demonstrates drawing samples of the prior and posterior to
provide the user with an intuitive understanding of their distributions.
'''


logging.basicConfig(level=logging.DEBUG)
np.random.seed(1)

# observation points
y = np.linspace(-7.5, 7.5, 25) 
# interpolation points
x = np.linspace(-7.5, 7.5, 1000) 
# true signal
u_true = np.exp(-0.3*np.abs(x))*np.sin(x)
# observation noise covariance
dsigma = np.full(25, 0.1)
dcov = np.diag(dsigma**2)
# noisy observations of the signal
d = np.exp(-0.3*np.abs(y))*np.sin(y) + np.random.normal(0.0, dsigma)

def negative_log_likelihood(params):
    log_eps, log_var = params
    gp = gpiso('se', eps=10**log_eps, var=10**log_var) 
    out = -gp.log_likelihood(y[:, None], d, dcov=dcov)
    return out

log_eps, log_var = minimize(negative_log_likelihood, [0.0, 0.0]).x
# create a prior GaussianProcess using the most likely variance and lengthscale

# Path: rbf/gproc.py
# def gpiso(phi, eps, var, dim=None):
#     '''
#     Creates an isotropic Gaussian process which has a covariance function that
#     is described by a radial basis function.
# 
#     Parameters
#     ----------
#     phi : str or RBF instance
#         Radial basis function describing the covariance function.
# 
#     eps : float
#         Shape parameter.
# 
#     var : float
#         Variance.
# 
#     dim : int, optional
#         Fixes the spatial dimensions of the domain.
# 
#     Returns
#     -------
#     GaussianProcess
# 
#     Notes
#     -----
#     Not all radial basis functions are positive definite, which means that it
#     is possible to instantiate an invalid `GaussianProcess`. The method
#     `is_positive_definite` provides a necessary but not sufficient test for
#     positive definiteness. Examples of predefined `RBF` instances which are
#     positive definite include: `rbf.basis.se`, `rbf.basis.ga`, `rbf.basis.exp`,
#     `rbf.basis.iq`, `rbf.basis.imq`.
# 
#     '''
#     phi = rbf.basis.get_rbf(phi)
# 
#     def isotropic_covariance(x1, x2, diff1, diff2):
#         diff = diff1 + diff2
#         coeff = var*(-1)**sum(diff2)
#         out = coeff*phi(x1, x2, eps=eps, diff=diff)
#         return out
# 
#     def isotropic_variance(x, diff):
#         coeff = var*(-1)**sum(diff)
#         value = coeff*phi.center_value(eps=eps, diff=2*diff)
#         out = np.full(x.shape[0], value)
#         return out
# 
#     out = GaussianProcess(
#         covariance=isotropic_covariance,
#         variance=isotropic_variance,
#         dim=dim,
#         differentiable=True
#         )
#     return out


# Path: docs/scripts/gproc.a.py
import logging
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from rbf.gproc import gpiso

''' 
This script demonstrates how to perform Gaussian process regression on a noisy
data set. It also demonstrates drawing samples of the prior and posterior to
provide the user with an intuitive understanding of their distributions.
'''


logging.basicConfig(level=logging.DEBUG)
np.random.seed(1)

# observation points
y = np.linspace(-7.5, 7.5, 25) 
# interpolation points
x = np.linspace(-7.5, 7.5, 1000) 
# true signal
u_true = np.exp(-0.3*np.abs(x))*np.sin(x)
# observation noise covariance
dsigma = np.full(25, 0.1)
dcov = np.diag(dsigma**2)
# noisy observations of the signal
d = np.exp(-0.3*np.abs(y))*np.sin(y) + np.random.normal(0.0, dsigma)

def negative_log_likelihood(params):
    log_eps, log_var = params
    gp = gpiso('se', eps=10**log_eps, var=10**log_var) 
    out = -gp.log_likelihood(y[:, None], d, dcov=dcov)
    return out

log_eps, log_var = minimize(negative_log_likelihood, [0.0, 0.0]).x
# create a prior GaussianProcess using the most likely variance and lengthscale

gp_prior = gpiso('se', eps=10**log_eps, var=10**log_var)