1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022012 Heuristic and Evolutionary Algorithms Laboratory (HEAL)


4  *


5  * This file is part of HeuristicLab.


6  *


7  * HeuristicLab is free software: you can redistribute it and/or modify


8  * it under the terms of the GNU General Public License as published by


9  * the Free Software Foundation, either version 3 of the License, or


10  * (at your option) any later version.


11  *


12  * HeuristicLab is distributed in the hope that it will be useful,


13  * but WITHOUT ANY WARRANTY; without even the implied warranty of


14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


15  * GNU General Public License for more details.


16  *


17  * You should have received a copy of the GNU General Public License


18  * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.


19  */


20  #endregion


21 


22  using System;


23  using System.Collections.Generic;


24  using System.Linq;


25  using HeuristicLab.Common;


26  using HeuristicLab.Core;


27  using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;


28  using HeuristicLab.Problems.DataAnalysis;


29 


30  namespace HeuristicLab.Algorithms.DataAnalysis {


31  /// <summary>


32  /// Represents a Gaussian process model.


33  /// </summary>


34  [StorableClass]


35  [Item("GaussianProcessModel", "Represents a Gaussian process posterior.")]


36  public sealed class GaussianProcessModel : NamedItem, IGaussianProcessModel {


37  [Storable]


38  private double negativeLogLikelihood;


39  public double NegativeLogLikelihood {


40  get { return negativeLogLikelihood; }


41  }


42 


43  [Storable]


44  private double[] hyperparameterGradients;


45  public double[] HyperparameterGradients {


46  get {


47  var copy = new double[hyperparameterGradients.Length];


48  Array.Copy(hyperparameterGradients, copy, copy.Length);


49  return copy;


50  }


51  }


52 


53  [Storable]


54  private ICovarianceFunction covarianceFunction;


55  public ICovarianceFunction CovarianceFunction {


56  get { return covarianceFunction; }


57  }


58  [Storable]


59  private IMeanFunction meanFunction;


60  public IMeanFunction MeanFunction {


61  get { return meanFunction; }


62  }


63  [Storable]


64  private string targetVariable;


65  public string TargetVariable {


66  get { return targetVariable; }


67  }


68  [Storable]


69  private string[] allowedInputVariables;


70  public string[] AllowedInputVariables {


71  get { return allowedInputVariables; }


72  }


73 


74  [Storable]


75  private double[] alpha;


76  [Storable]


77  private double sqrSigmaNoise;


78  public double SigmaNoise {


79  get { return Math.Sqrt(sqrSigmaNoise); }


80  }


81 


82  [Storable]


83  private double[,] l;


84 


85  [Storable]


86  private double[,] x;


87  [Storable]


88  private Scaling inputScaling;


89 


90 


91  [StorableConstructor]


92  private GaussianProcessModel(bool deserializing) : base(deserializing) { }


93  private GaussianProcessModel(GaussianProcessModel original, Cloner cloner)


94  : base(original, cloner) {


95  this.meanFunction = cloner.Clone(original.meanFunction);


96  this.covarianceFunction = cloner.Clone(original.covarianceFunction);


97  this.inputScaling = cloner.Clone(original.inputScaling);


98  this.negativeLogLikelihood = original.negativeLogLikelihood;


99  this.targetVariable = original.targetVariable;


100  this.sqrSigmaNoise = original.sqrSigmaNoise;


101 


102  // shallow copies of arrays because they cannot be modified


103  this.allowedInputVariables = original.allowedInputVariables;


104  this.alpha = original.alpha;


105  this.l = original.l;


106  this.x = original.x;


107  }


108  public GaussianProcessModel(Dataset ds, string targetVariable, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows,


109  IEnumerable<double> hyp, IMeanFunction meanFunction, ICovarianceFunction covarianceFunction)


110  : base() {


111  this.name = ItemName;


112  this.description = ItemDescription;


113  this.meanFunction = (IMeanFunction)meanFunction.Clone();


114  this.covarianceFunction = (ICovarianceFunction)covarianceFunction.Clone();


115  this.targetVariable = targetVariable;


116  this.allowedInputVariables = allowedInputVariables.ToArray();


117 


118 


119  int nVariables = this.allowedInputVariables.Length;


120  this.meanFunction.SetParameter(hyp


121  .Take(this.meanFunction.GetNumberOfParameters(nVariables))


122  .ToArray());


123  this.covarianceFunction.SetParameter(hyp.Skip(this.meanFunction.GetNumberOfParameters(nVariables))


124  .Take(this.covarianceFunction.GetNumberOfParameters(nVariables))


125  .ToArray());


126  sqrSigmaNoise = Math.Exp(2.0 * hyp.Last());


127 


128  CalculateModel(ds, rows);


129  }


130 


131  private void CalculateModel(Dataset ds, IEnumerable<int> rows) {


132  inputScaling = new Scaling(ds, allowedInputVariables, rows);


133  x = AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputVariables, rows, inputScaling);


134  var y = ds.GetDoubleValues(targetVariable, rows);


135 


136  int n = x.GetLength(0);


137  l = new double[n, n];


138 


139  // calculate means and covariances


140  double[] m = meanFunction.GetMean(x);


141  for (int i = 0; i < n; i++) {


142  for (int j = i; j < n; j++) {


143  l[j, i] = covarianceFunction.GetCovariance(x, i, j) / sqrSigmaNoise;


144  if (j == i) l[j, i] += 1.0;


145  }


146  }


147 


148  // cholesky decomposition


149  int info;


150  alglib.densesolverreport denseSolveRep;


151 


152  var res = alglib.trfac.spdmatrixcholesky(ref l, n, false);


153  if (!res) throw new ArgumentException("Matrix is not positive semidefinite");


154 


155  // calculate sum of diagonal elements for likelihood


156  double diagSum = Enumerable.Range(0, n).Select(i => Math.Log(l[i, i])).Sum();


157 


158  // solve for alpha


159  double[] ym = y.Zip(m, (a, b) => a  b).ToArray();


160 


161  alglib.spdmatrixcholeskysolve(l, n, false, ym, out info, out denseSolveRep, out alpha);


162  for (int i = 0; i < alpha.Length; i++)


163  alpha[i] = alpha[i] / sqrSigmaNoise;


164  negativeLogLikelihood = 0.5 * Util.ScalarProd(ym, alpha) + diagSum + (n / 2.0) * Math.Log(2.0 * Math.PI * sqrSigmaNoise);


165 


166  // derivatives


167  int nAllowedVariables = x.GetLength(1);


168 


169  alglib.matinvreport matInvRep;


170  double[,] lCopy = new double[l.GetLength(0), l.GetLength(1)];


171  Array.Copy(l, lCopy, lCopy.Length);


172 


173  alglib.spdmatrixcholeskyinverse(ref lCopy, n, false, out info, out matInvRep);


174  if (info != 1) throw new ArgumentException("Can't invert matrix to calculate gradients.");


175  for (int i = 0; i < n; i++) {


176  for (int j = 0; j <= i; j++)


177  lCopy[i, j] = lCopy[i, j] / sqrSigmaNoise  alpha[i] * alpha[j];


178  }


179 


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


181 


182  double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];


183  for (int i = 0; i < meanGradients.Length; i++) {


184  var meanGrad = meanFunction.GetGradients(i, x);


185  meanGradients[i] = Util.ScalarProd(meanGrad, alpha);


186  }


187 


188  double[] covGradients = new double[covarianceFunction.GetNumberOfParameters(nAllowedVariables)];


189  if (covGradients.Length > 0) {


190  for (int i = 0; i < n; i++) {


191  for (int j = 0; j < i; j++) {


192  var g = covarianceFunction.GetGradient(x, i, j).ToArray();


193  for (int k = 0; k < covGradients.Length; k++) {


194  covGradients[k] += lCopy[i, j] * g[k];


195  }


196  }


197 


198  var gDiag = covarianceFunction.GetGradient(x, i, i).ToArray();


199  for (int k = 0; k < covGradients.Length; k++) {


200  // diag


201  covGradients[k] += 0.5 * lCopy[i, i] * gDiag[k];


202  }


203  }


204  }


205 


206  hyperparameterGradients =


207  meanGradients


208  .Concat(covGradients)


209  .Concat(new double[] { noiseGradient }).ToArray();


210 


211  }


212 


213 


214  public override IDeepCloneable Clone(Cloner cloner) {


215  return new GaussianProcessModel(this, cloner);


216  }


217 


218  #region IRegressionModel Members


219  public IEnumerable<double> GetEstimatedValues(Dataset dataset, IEnumerable<int> rows) {


220  return GetEstimatedValuesHelper(dataset, rows);


221  }


222  public GaussianProcessRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {


223  return new GaussianProcessRegressionSolution(this, new RegressionProblemData(problemData));


224  }


225  IRegressionSolution IRegressionModel.CreateRegressionSolution(IRegressionProblemData problemData) {


226  return CreateRegressionSolution(problemData);


227  }


228  #endregion


229 


230 


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


232  var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);


233  int newN = newX.GetLength(0);


234  int n = x.GetLength(0);


235  var Ks = new double[newN, n];


236  var ms = meanFunction.GetMean(newX);


237  for (int i = 0; i < newN; i++) {


238  for (int j = 0; j < n; j++) {


239  Ks[i, j] = covarianceFunction.GetCrossCovariance(x, newX, j, i);


240  }


241  }


242 


243  return Enumerable.Range(0, newN)


244  .Select(i => ms[i] + Util.ScalarProd(Util.GetRow(Ks, i), alpha));


245  }


246 


247  public IEnumerable<double> GetEstimatedVariance(Dataset dataset, IEnumerable<int> rows) {


248  var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);


249  int newN = newX.GetLength(0);


250  int n = x.GetLength(0);


251 


252  var kss = new double[newN];


253  double[,] sWKs = new double[n, newN];


254 


255  // for stddev


256  for (int i = 0; i < newN; i++)


257  kss[i] = covarianceFunction.GetCovariance(newX, i, i);


258 


259  for (int i = 0; i < newN; i++) {


260  for (int j = 0; j < n; j++) {


261  sWKs[j, i] = covarianceFunction.GetCrossCovariance(x, newX, j, i) / Math.Sqrt(sqrSigmaNoise);


262  }


263  }


264 


265  // for stddev


266  alglib.ablas.rmatrixlefttrsm(n, newN, l, 0, 0, false, false, 0, ref sWKs, 0, 0);


267 


268  for (int i = 0; i < newN; i++) {


269  var sumV = Util.ScalarProd(Util.GetCol(sWKs, i), Util.GetCol(sWKs, i));


270  kss[i] = sumV;


271  if (kss[i] < 0) kss[i] = 0;


272  }


273  return kss;


274  }


275  }


276  }

