
import numpy      as np
import GPflow     as gpflow
import tensorflow as tf
import os

def getData( ):
    rng = np.random.RandomState(1)
    N   = 30
    X   = rng.rand(N,1)
    Y   = np.sin(12*X) + 0.66*np.cos(25*X) + rng.randn(N,1)*0.1 + 3
    return X, Y

def getRegressionModel( X , Y ):
    # build the GPR object
    k                     = gpflow.kernels.Matern52(1)
    meanf                 = gpflow.mean_functions.Linear(1,0)
    m                     = gpflow.models.GPR(X, Y, k, meanf)
    m.likelihood.variance = 0.01
    print("Here are the parameters before optimization")
    m
    return m

def optimizeModel( m ):
    opt = gpflow.train.ScipyOptimizer()
    opt.minimize(m)
    print("Here are the parameters after optimization")
    m

def setModelPriors( m ):
    # we'll choose rather arbitrary priors.
    m.kern.lengthscales.prior   = gpflow.priors.Gamma(1., 1.)
    m.kern.variance.prior       = gpflow.priors.Gamma(1., 1.)
    m.likelihood.variance.prior = gpflow.priors.Gamma(1., 1.)
    m.mean_function.A.prior     = gpflow.priors.Gaussian(0., 10.)
    m.mean_function.b.prior     = gpflow.priors.Gaussian(0., 10.)
    print("model with priors ", m)

def getSamples( m ):
    hmc     = gpflow.train.HMC()
    samples = hmc.sample(m, num_samples=100, epsilon = 0.1)
    return samples

def runExperiments( sampling=True, outputGraphs=False ):
    X,Y = getData()
    m   = getRegressionModel(X,Y)
    optimizeModel(m)
    if sampling:
        setModelPriors(m)
        samples = getSamples(m)

if __name__ == '__main__':
    runExperiments()
    # cProfile.run('runExperiments(plotting=False)')