source: branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Algorithms.DataAnalysis.Experimental/ToeplitzGaussianProcessModel.cs @ 11598

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2012 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
 *
 * This file is part of HeuristicLab.
 *
 * HeuristicLab is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * HeuristicLab is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
 */
#endregion

using System;
using System.Collections.Generic;
using System.Linq;
using System.Runtime.InteropServices;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Problems.DataAnalysis;
using HeuristicLabEigen;
using ILNumerics;

namespace HeuristicLab.Algorithms.DataAnalysis {
  /// <summary>
  /// Represents a Gaussian process model.
  /// </summary>
  [StorableClass]
  [Item("ToeplitzGaussianProcessModel", "Gaussian process model implemented using ILNumerics.")]
  public sealed class ToeplitzGaussianProcessModel : NamedItem, IGaussianProcessModel {

    [DllImport("toeplitz.dll", CallingConvention = CallingConvention.Cdecl)]
    public static extern void r4to_sl_(ref int n, float[] a, float[] x, float[] b, ref int job);



    [Storable]
    private double negativeLogLikelihood;
    public double NegativeLogLikelihood {
      get { return negativeLogLikelihood; }
    }

    [Storable]
    private double[] hyperparameterGradients;
    public double[] HyperparameterGradients {
      get {
        var copy = new double[hyperparameterGradients.Length];
        Array.Copy(hyperparameterGradients, copy, copy.Length);
        return copy;
      }
    }

    [Storable]
    private ICovarianceFunction covarianceFunction;
    public ICovarianceFunction CovarianceFunction {
      get { return covarianceFunction; }
    }
    [Storable]
    private IMeanFunction meanFunction;
    public IMeanFunction MeanFunction {
      get { return meanFunction; }
    }
    [Storable]
    private string targetVariable;
    public string TargetVariable {
      get { return targetVariable; }
    }
    [Storable]
    private string[] allowedInputVariables;
    public string[] AllowedInputVariables {
      get { return allowedInputVariables; }
    }


    [Storable]
    private double sqrSigmaNoise;
    public double SigmaNoise {
      get { return Math.Sqrt(sqrSigmaNoise); }
    }

    [Storable]
    private double[] meanParameter;
    [Storable]
    private double[] covarianceParameter;
    [Storable]
    private float[] alpha;


    [Storable]
    private double[,] x;
    [Storable]
    private Scaling inputScaling;

    [StorableConstructor]
    private ToeplitzGaussianProcessModel(bool deserializing) : base(deserializing) { }
    private ToeplitzGaussianProcessModel(ToeplitzGaussianProcessModel original, Cloner cloner)
      : base(original, cloner) {
      this.meanFunction = cloner.Clone(original.meanFunction);
      this.covarianceFunction = cloner.Clone(original.covarianceFunction);
      this.inputScaling = cloner.Clone(original.inputScaling);
      this.negativeLogLikelihood = original.negativeLogLikelihood;
      this.targetVariable = original.targetVariable;
      this.sqrSigmaNoise = original.sqrSigmaNoise;
      if (original.meanParameter != null) {
        this.meanParameter = (double[])original.meanParameter.Clone();
      }
      if (original.covarianceParameter != null) {
        this.covarianceParameter = (double[])original.covarianceParameter.Clone();
      }

      // shallow copies of arrays because they cannot be modified
      this.allowedInputVariables = original.allowedInputVariables; 
      this.alpha = original.alpha;
      this.x = original.x;
    }
    public ToeplitzGaussianProcessModel(Dataset ds, string targetVariable, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows,
      IEnumerable<double> hyp, IMeanFunction meanFunction, ICovarianceFunction covarianceFunction)
      : base() {
      this.name = ItemName;
      this.description = ItemDescription;
      this.meanFunction = (IMeanFunction)meanFunction.Clone();
      this.covarianceFunction = (ICovarianceFunction)covarianceFunction.Clone();
      this.targetVariable = targetVariable;
      this.allowedInputVariables = allowedInputVariables.ToArray();


      int nVariables = this.allowedInputVariables.Length;
      meanParameter = hyp
        .Take(this.meanFunction.GetNumberOfParameters(nVariables))
        .ToArray();

      covarianceParameter = hyp.Skip(this.meanFunction.GetNumberOfParameters(nVariables))
                                             .Take(this.covarianceFunction.GetNumberOfParameters(nVariables))
                                             .ToArray();
      sqrSigmaNoise = Math.Exp(2.0 * hyp.Last());

      CalculateModel(ds, rows);
    }

    private void CalculateModel(Dataset ds, IEnumerable<int> rows) {
      inputScaling = new Scaling(ds, allowedInputVariables, rows);
      x = AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputVariables, rows, inputScaling);
      var y = ds.GetDoubleValues(targetVariable, rows);

      int n = x.GetLength(0);
      var l = new float[2 * n - 1];

      // calculate means and covariances
      var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, Enumerable.Range(0, x.GetLength(1)));
      var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
      var l0 = (float)(cov.Covariance(x, 0, 0) / sqrSigmaNoise + 1.0);
      l[0] = 1.0f;
      for (int i = 1; i < n; i++) {
        l[i] = (float)((cov.Covariance(x, 0, i) / sqrSigmaNoise)/ l0);
      }
      for (int i = n; i < 2 * n - 1; i++) {
        l[i] = (float)((cov.Covariance(x, i - n + 1, 0) / sqrSigmaNoise) / l0);
      }


      int info = 0;
      alpha = new float[n];

      // solve for alpha
      float[] ym = y.Zip(Enumerable.Range(0, x.GetLength(0)).Select(r => mean.Mean(x, r)), (a, b) => a - b)
        .Select(e => (float)e)
        .ToArray();
      double nll;

      //unsafe {
      //  fixed (double* ap = &alpha[0])
      //  fixed (double* ymp = &ym[0])
      //  fixed (double* invlP = &invL[0])
      //  fixed (double* lp = &l[0]) {
      //    myEigen.Solve(lp, ymp, ap, invlP, sqrSigmaNoise, n, &nll, &info);
      //  }
      //}
      int job = 0;
      r4to_sl_(ref n, l, ym, alpha, ref job);
      alpha = alpha.Select(a => (float)(a / sqrSigmaNoise / l0)).ToArray();


      var yDotAlpha = ym.Zip(alpha, (f, f1) => f * f1).Sum();

      var lambda = toeplitz_cholesky_diag(n, l);
      var logDetK = 2 * lambda.Sum(f => Math.Log(f * Math.Sqrt(l0)));


      this.negativeLogLikelihood = 0.5 * yDotAlpha + 0.5 * logDetK + n / 2.0 * Math.Log(2.0 * Math.PI * sqrSigmaNoise);

      //double noiseGradient = sqrSigmaNoise * Enumerable.Range(0, n).Select(i => invL[i + i * n]).Sum();

      // derivatives
      int nAllowedVariables = x.GetLength(1);
      //
      //double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];
      //for (int k = 0; k < meanGradients.Length; k++) {
      //  var meanGrad = Enumerable.Range(0, alpha.Length)
      //  .Select(r => mean.Gradient(x, r, k));
      //  meanGradients[k] = -meanGrad.Zip(alpha, (d1, f) => d1 * f).Sum();
      //}
      //
      //double[] covGradients = new double[covarianceFunction.GetNumberOfParameters(nAllowedVariables)];
      //if (covGradients.Length > 0) {
      //  for (int i = 0; i < n; i++) {
      //    for (int j = 0; j < i; j++) {
      //      var g = cov.CovarianceGradient(x, i, j).ToArray();
      //      for (int k = 0; k < covGradients.Length; k++) {
      //        covGradients[k] += invL[j + i * n] * g[k];
      //      }
      //    }
      //
      //    var gDiag = cov.CovarianceGradient(x, i, i).ToArray();
      //    for (int k = 0; k < covGradients.Length; k++) {
      //      // diag
      //      covGradients[k] += 0.5 * invL[i + i * n] * gDiag[k];
      //    }
      //  }
      //}
      //
      hyperparameterGradients =
        Enumerable.Repeat(0.0,
                          meanFunction.GetNumberOfParameters(nAllowedVariables) +
                          covarianceFunction.GetNumberOfParameters(nAllowedVariables) + 1).ToArray();
      //hyperparameterGradients =
      //  meanGradients
      //  .Concat(covGradients)
      //  .Concat(new double[] { noiseGradient }).ToArray();

    }

    public static double[] toeplitz_cholesky_diag(int n, float[] a) {
      double div;
      double[] g;
      double g1j;
      double g2j;
      int i;
      int j;
      double[] l;
      double rho;

      l = new double[n];


      g = new double[2 * n];

      for (j = 0; j < n; j++) {
        g[0 + j * 2] = a[j];
      }
      g[1 + 0 * 2] = 0.0;
      for (j = 1; j < n; j++) {
        g[1 + j * 2] = a[j + n - 1];
      }

      l[0] = g[0];
      //for (i = 0; i < n; i++) {
      //  l[i, 0] = g[0 + i * 2];
      //}
      for (j = n - 1; 1 <= j; j--) {
        g[0 + j * 2] = g[0 + (j - 1) * 2];
      }
      g[0 + 0 * 2] = 0.0;

      for (i = 1; i < n; i++) {
        rho = -g[1 + i * 2] / g[0 + i * 2];
        div = Math.Sqrt((1.0 - rho) * (1.0 + rho));

        for (j = i; j < n; j++) {
          g1j = g[0 + j * 2];
          g2j = g[1 + j * 2];
          g[0 + j * 2] = (g1j + rho * g2j) / div;
          g[1 + j * 2] = (rho * g1j + g2j) / div;
        }

        l[i] = g[i * 2];
        //for (j = i; j < n; j++) {
        //  l[j, i] = g[0 + j * 2];
        //}
        for (j = n - 1; i < j; j--) {
          g[0 + j * 2] = g[0 + (j - 1) * 2];
        }
        g[0 + i * 2] = 0.0;
      }


      return l;
    }

    public override IDeepCloneable Clone(Cloner cloner) {
      return new ToeplitzGaussianProcessModel(this, cloner);
    }

    // is called by the solution creator to set all parameter values of the covariance and mean function
    // to the optimized values (necessary to make the values visible in the GUI)
    public void FixParameters() {
      covarianceFunction.SetParameter(covarianceParameter);
      meanFunction.SetParameter(meanParameter);
      covarianceParameter = new double[0];
      meanParameter = new double[0];
    }

    #region IRegressionModel Members
    public IEnumerable<double> GetEstimatedValues(Dataset dataset, IEnumerable<int> rows) {
      return GetEstimatedValuesHelper(dataset, rows);
    }
    public GaussianProcessRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
      return new GaussianProcessRegressionSolution(this, new RegressionProblemData(problemData));
    }
    IRegressionSolution IRegressionModel.CreateRegressionSolution(IRegressionProblemData problemData) {
      return CreateRegressionSolution(problemData);
    }
    #endregion


    private IEnumerable<double> GetEstimatedValuesHelper(Dataset dataset, IEnumerable<int> rows) {

      var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);
      int newN = newX.GetLength(0);
      int n = x.GetLength(0);
      var Ks = new double[newN, n];
      var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, Enumerable.Range(0, newX.GetLength(1)));
      var ms = Enumerable.Range(0, newX.GetLength(0))
      .Select(r => mean.Mean(newX, r))
      .ToArray();
      var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
      for (int i = 0; i < newN; i++) {
        for (int j = 0; j < n; j++) {
          Ks[i, j] = cov.CrossCovariance(x, newX, j, i);
        }
      }

      return Enumerable.Range(0, newN)
        .Select(i => ms[i] + Util.ScalarProd(Util.GetRow(Ks, i), alpha.Select(a=>(double)a)));
    }

    public IEnumerable<double> GetEstimatedVariance(Dataset dataset, IEnumerable<int> rows) {
      return rows.Select(r => sqrSigmaNoise);
    }
  }
}