from __future__ import print_function import gpflow import tensorflow as tf import os import numpy as np import cProfile import csv tol=1e-11 nRepeats = 50 predict_limits = [-4., 4.] inducing_points_limits = [-1., 9] hold_out_limits = [0.20, 0.60] optimization_limits = [18., 25.] def readCsvFile(fileName): reader = csv.reader(open(fileName,'r')) dataList = [] for row in reader: dataList.append([float(elem) for elem in row]) return np.array(dataList) def getTrainingTestData(): overallX = readCsvFile('data/snelson_train_inputs') overallY = readCsvFile('data/snelson_train_outputs') trainIndeces = [] testIndeces = [] nPoints = overallX.shape[0] for index in range(nPoints): if index % 4 == 0: trainIndeces.append(index) else: testIndeces.append(index) xtrain = overallX[trainIndeces,:] xtest = overallX[testIndeces, :] ytrain = overallY[trainIndeces, :] ytest = overallY[testIndeces, :] return xtrain,ytrain,xtest,ytest def getLogPredictiveDensities(targetValues, means, variances): assert(targetValues.flatten().shape == targetValues.shape) assert(means.flatten().shape == means.shape) assert(variances.flatten().shape == variances.shape) assert(len(targetValues) == len(means)) assert(len(variances) == len(means)) deltas = targetValues - means mahalanobisTerms = -0.5*deltas**2/variances normalizationTerms = -0.5 * np.log(variances) - 0.5 * np.log(2.*np.pi) return mahalanobisTerms + normalizationTerms def getKernel(): return gpflow.kernels.RBF(1) def getRegressionModel(X,Y): m = gpflow.models.GPR(X, Y, kern=getKernel()) m.likelihood.variance = 1. m.kern.lengthscales = 1. m.kern.variance = 1. return m def getSparseModel(X,Y,isFITC=False): if not(isFITC): m = gpflow.models.SGPR(X, Y, kern=getKernel(), Z=X.copy()) else: m = gpflow.models.GPRFITC(X, Y, kern=getKernel(), Z=X.copy()) return m def printModelParameters(model): print("Likelihood variance ", model.likelihood.variance, "\n") print("Kernel variance ", model.kern.variance, "\n") print("Kernel lengthscale ", model.kern.lengthscales, "\n") def plotPredictions(ax, model, color, label): xtest = np.sort(readCsvFile('data/snelson_test_inputs')) predMean, predVar = model.predict_y(xtest) ax.plot(xtest, predMean, color, label=label) ax.plot(xtest, predMean + 2.*np.sqrt(predVar),color) ax.plot(xtest, predMean - 2.*np.sqrt(predVar), color) def trainSparseModel(xtrain,ytrain,exact_model,isFITC, xtest, ytest): sparse_model = getSparseModel(xtrain,ytrain,isFITC) sparse_model.likelihood.variance = exact_model.likelihood.variance.read_value().copy() sparse_model.kern.lengthscales = exact_model.kern.lengthscales.read_value().copy() sparse_model.kern.variance = exact_model.kern.variance.read_value().copy() callback = cb(sparse_model, xtest, ytest) for repeatIndex in range(nRepeats): print("repeatIndex ",repeatIndex) sparse_model.optimize(disp=False, maxiter= 2000, tol=tol, callback=callback) return sparse_model, callback def plotComparisonFigure(xtrain, sparse_model,exact_model, ax_predictions, ax_inducing_points, ax_optimization, iterations, log_likelihoods,hold_out_likelihood, title): plotPredictions(ax_predictions, exact_model, 'g', label='Exact') plotPredictions(ax_predictions, sparse_model, 'b', label='Approximate') ax_predictions.legend(loc=9) ax_predictions.plot( sparse_model.feature.Z.value , -1.*np.ones( xtrain.shape ), 'ko' ) ax_predictions.set_ylim( predict_limits ) ax_inducing_points.plot( xtrain, sparse_model.feature.Z.value, 'bo' ) xs= np.linspace( ax_inducing_points.get_xlim()[0], ax_inducing_points.get_xlim()[1], 200 ) ax_inducing_points.plot( xs, xs, 'g' ) ax_inducing_points.set_xlabel('Optimal inducing point position') ax_inducing_points.set_ylabel('Learnt inducing point position') ax_inducing_points.set_ylim(inducing_points_limits) ax_optimization.plot(iterations, -1.*np.array(log_likelihoods), 'g-') ax_optimization.set_ylim(optimization_limits) ax2 = ax_optimization.twinx() ax2.plot(iterations, -1.*np.array(hold_out_likelihood), 'b-') ax_optimization.set_xlabel('Minimization iterations') ax_optimization.set_ylabel('Minimization objective', color='g') ax2.set_ylim(hold_out_limits) ax2.set_ylabel('Hold out negative log likelihood', color='b') class cb(): def __init__(self, model, xtest, ytest, holdout_inverval=100): self.model = model self.holdout_inverval = holdout_inverval self.xtest = xtest self.ytest = ytest self.log_likelihoods = [] self.hold_out_likelihood = [] self.n_iters = [] self.counter = 0 def __call__(self, info): if (self.counter%self.holdout_inverval) == 0 or (self.counter <= 10): predictive_mean, predictive_variance = self.model.predict_y(self.xtest) self.n_iters.append(self.counter) self.log_likelihoods.append(self.model.compute_log_likelihood()) self.hold_out_likelihood.append(getLogPredictiveDensities(self.ytest.flatten() , predictive_mean.flatten(), predictive_variance.flatten()).mean()) self.counter+=1 def stretch(lenNIters, initialValues): stretched = np.ones(lenNIters) * initialValues[-1] stretched[0:len(initialValues)] = initialValues return stretched def snelsonDemo(): from matplotlib import pyplot as plt from IPython import embed xtrain,ytrain,xtest,ytest = getTrainingTestData() #run exact inference on training data. exact_model = getRegressionModel(xtrain,ytrain) exact_model.optimize(maxiter= 2000000, tol=tol) figA, axes = plt.subplots(1,1) inds = np.argsort(xtrain.flatten()) axes.plot(xtrain[inds,:], ytrain[inds,:], 'ro') plotPredictions(axes, exact_model, 'g', None) figB, axes = plt.subplots(3,2) #run sparse model on training data intialized from exact optimal solution. VFEmodel, VFEcb = trainSparseModel(xtrain,ytrain,exact_model,False,xtest,ytest) FITCmodel, FITCcb = trainSparseModel(xtrain,ytrain,exact_model,True,xtest,ytest) print("Exact model parameters \n") printModelParameters(exact_model) print("Sparse model parameters for VFE optimization \n") printModelParameters(VFEmodel) print("Sparse model parameters for FITC optimization \n") printModelParameters(FITCmodel) VFEiters = FITCcb.n_iters VFElog_likelihoods = stretch(len(VFEiters), VFEcb.log_likelihoods) VFEhold_out_likelihood = stretch(len(VFEiters), VFEcb.hold_out_likelihood) plotComparisonFigure(xtrain, VFEmodel, exact_model, axes[0,0], axes[1,0], axes[2,0], VFEiters, VFElog_likelihoods.tolist(), VFEhold_out_likelihood.tolist(), "VFE") plotComparisonFigure(xtrain, FITCmodel, exact_model, axes[0,1], axes[1,1], axes[2,1],FITCcb.n_iters, FITCcb.log_likelihoods, FITCcb.hold_out_likelihood , "FITC") axes[0,0].set_title('VFE', loc='center', fontdict = {'fontsize': 22}) axes[0,1].set_title('FITC', loc='center', fontdict = {'fontsize': 22}) embed() if __name__ == '__main__': snelsonDemo()