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 


79  [Storable]


80  private double[,] l;


81 


82  [Storable]


83  private double[,] x;


84  [Storable]


85  private Scaling inputScaling;


86 


87 


88  [StorableConstructor]


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


90  private GaussianProcessModel(GaussianProcessModel original, Cloner cloner)


91  : base(original, cloner) {


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


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


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


95  this.negativeLogLikelihood = original.negativeLogLikelihood;


96  this.targetVariable = original.targetVariable;


97  this.sqrSigmaNoise = original.sqrSigmaNoise;


98 


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


100  this.allowedInputVariables = original.allowedInputVariables;


101  this.alpha = original.alpha;


102  this.l = original.l;


103  this.x = original.x;


104  }


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


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


107  : base() {


108  this.name = ItemName;


109  this.description = ItemDescription;


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


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


112  this.targetVariable = targetVariable;


113  this.allowedInputVariables = allowedInputVariables.ToArray();


114 


115 


116  int nVariables = this.allowedInputVariables.Length;


117  this.meanFunction.SetParameter(hyp


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


119  .ToArray());


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


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


122  .ToArray());


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


124 


125  CalculateModel(ds, rows);


126  }


127 


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


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


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


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


132 


133  int n = x.GetLength(0);


134  l = new double[n, n];


135 


136  meanFunction.SetData(x);


137 


138  // calculate means and covariances


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


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


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


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


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


144  }


145  }


146 


147  // cholesky decomposition


148  int info;


149  alglib.densesolverreport denseSolveRep;


150 


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


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


153 


154  // calculate sum of diagonal elements for likelihood


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


156 


157  // solve for alpha


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


159 


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


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


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


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


164 


165  // derivatives


166  int nAllowedVariables = x.GetLength(1);


167 


168  alglib.matinvreport matInvRep;


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


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


171 


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


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


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


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


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


177  }


178 


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


180 


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


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


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


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


185  }


186 


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


188  if (covGradients.Length > 0) {


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


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


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


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


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


194  }


195  }


196 


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


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


199  // diag


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


201  }


202  }


203  }


204 


205  hyperparameterGradients =


206  meanGradients


207  .Concat(covGradients)


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


209 


210  }


211 


212 


213  public override IDeepCloneable Clone(Cloner cloner) {


214  return new GaussianProcessModel(this, cloner);


215  }


216 


217  #region IRegressionModel Members


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


219  return GetEstimatedValuesHelper(dataset, rows);


220  }


221  public GaussianProcessRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {


222  return new GaussianProcessRegressionSolution(this, problemData);


223  }


224  IRegressionSolution IRegressionModel.CreateRegressionSolution(IRegressionProblemData problemData) {


225  return CreateRegressionSolution(problemData);


226  }


227  #endregion


228 


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


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


231  int newN = newX.GetLength(0);


232  int n = x.GetLength(0);


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


234  meanFunction.SetData(newX);


235  var ms = meanFunction.GetMean(newX);


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


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


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


239  }


240  }


241 


242  return Enumerable.Range(0, newN)


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


244  }


245 


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


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


248  int newN = newX.GetLength(0);


249  int n = x.GetLength(0);


250 


251  var kss = new double[newN];


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


253 


254  // for stddev


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


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


257 


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


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


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


261  }


262  }


263 


264  // for stddev


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


266 


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


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


269  kss[i] = sumV;


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


271  }


272  return kss;


273  }


274  }


275  }

