# 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 numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fmin
from rbf.gproc import gpiso
''' 
This script demonstrate how to optimize the hyperparameters for a 
Gaussian process based on the marginal likelihood. Optimization is 
performed in two ways, first with a grid search method and then with a 
downhill simplex method.
'''
np.random.seed(3)

# True signal which we want to recover. This is a exponential function
# with mean=0.0, variance=2.0, and time-scale=0.1. For graphical
# purposes, we will only estimate the variance and time-scale.
eps = 0.1
var = 2.0
gp = gpiso('exp', eps=eps, var=var)

n = 500 # number of observations
time = np.linspace(-5.0, 5.0, n)[:,None] # observation points
data = gp.sample(time) # signal which we want to describe

# find the optimal hyperparameter with a brute force grid search
eps_search = 10**np.linspace(-2, 1, 30)
var_search = 10**np.linspace(-2, 2, 30)
log_likelihoods = np.zeros((30, 30))
for i, eps_test in enumerate(eps_search): 
  for j, var_test in enumerate(var_search): 

# 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.f.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fmin
from rbf.gproc import gpiso

''' 
This script demonstrate how to optimize the hyperparameters for a 
Gaussian process based on the marginal likelihood. Optimization is 
performed in two ways, first with a grid search method and then with a 
downhill simplex method.
'''
np.random.seed(3)

# True signal which we want to recover. This is a exponential function
# with mean=0.0, variance=2.0, and time-scale=0.1. For graphical
# purposes, we will only estimate the variance and time-scale.
eps = 0.1
var = 2.0
gp = gpiso('exp', eps=eps, var=var)

n = 500 # number of observations
time = np.linspace(-5.0, 5.0, n)[:,None] # observation points
data = gp.sample(time) # signal which we want to describe

# find the optimal hyperparameter with a brute force grid search
eps_search = 10**np.linspace(-2, 1, 30)
var_search = 10**np.linspace(-2, 2, 30)
log_likelihoods = np.zeros((30, 30))
for i, eps_test in enumerate(eps_search): 
  for j, var_test in enumerate(var_search): 

    gp = gpiso('exp', eps=eps_test, var=var_test)