package func;

import dist.AbstractConditionalDistribution;
import dist.Distribution;
import dist.UnivariateGaussian;
import shared.DataSet;
import shared.Instance;
import util.linalg.CholeskyFactorization;
import util.linalg.DenseVector;
import util.linalg.Matrix;
import util.linalg.RectangularMatrix;
import util.linalg.Vector;
import func.svm.Kernel;
import func.svm.LinearKernel;

/**
 * A gaussian process regression
 * @author Andrew Guillory gtg008g@mail.gatech.edu
 * @version 1.0
 */
public class GaussianProcessRegression extends AbstractConditionalDistribution implements FunctionApproximater {
    /**
     * The kernel
     */
    private Kernel kernel;
    /**
     * The noise sigma
     */
    private double noiseSigma;
    /**
     * The kernel matrix
     */
    private Matrix c;
    /**
     * The cholesky factorization of the kernel matrix
     */
    private CholeskyFactorization cf;
    /**
     * The a values
     */
    private Vector a;
    
    /**
     * Make a new gaussian process regression
     * @param k the kernel to use
     * @param noise the noise sigma value
     */
    public GaussianProcessRegression(Kernel k, double noise) {
        this.kernel = k;
        this.noiseSigma = noise;
    }
    /**
     * Make a new default gaussian process regression
     */
    public GaussianProcessRegression() {
        this(new LinearKernel(), 1);
    }
    /**
     * @see func.FunctionApproximater#estimate(shared.DataSet)
     */
    public void estimate(DataSet set) {
        // make the kernel matrix
        c = new RectangularMatrix(set.size(), set.size());
        kernel.setExamples(set);
        for (int i = 0; i < c.m(); i++) {
            for (int j = 0; j < c.n(); j++) {
                c.set(i,j, kernel.value(i,j));
            }
        }
        // add in the noise
        c.plusEquals(RectangularMatrix.eye(set.size()).times(noiseSigma * noiseSigma));
        // make the target vector
        Vector t = new DenseVector(set.size());
        for (int i = 0; i < t.size(); i++) {
            t.set(i, set.get(i).getLabel().getContinuous());
        }
        // solve for the weights
        cf = new CholeskyFactorization(c);
        a = cf.solve(t);
    }
    
    /**
     * @see func.FunctionApproximater#value(shared.Instance)
     */
    public Instance value(Instance i) {
        return distributionFor(i).mode();
    }
    
    /**
     * @see dist.ConditionalDistribution#distributionFor(shared.Instance)
     */
    public Distribution distributionFor(Instance instance) {
        Vector k = new DenseVector(c.m());
        for (int i = 0; i < k.size(); i++) {
            k.set(i, kernel.value(i, instance));
        }        
        double mean = a.dotProduct(k);
        double sigma = Math.sqrt(
            kernel.value(instance, instance) - k.dotProduct(cf.solve(k)));
        return new UnivariateGaussian(mean, sigma);
    }
    

}