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.Linq;


24  using HeuristicLab.Common;


25  using HeuristicLab.Core;


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


27 


28  namespace HeuristicLab.Algorithms.DataAnalysis {


29  [StorableClass]


30  [Item(Name = "CovarianceRQiso",


31  Description = "Isotropic rational quadratic covariance function for Gaussian processes.")]


32  public class CovarianceRQiso : Item, ICovarianceFunction {


33  [Storable]


34  private double[,] x;


35  [Storable]


36  private double[,] xt;


37  [Storable]


38  private double sf2;


39  public double Scale { get { return sf2; } }


40  [Storable]


41  private double l;


42  public double Length { get { return l; } }


43  [Storable]


44  private double alpha;


45  public double Shape { get { return alpha; } }


46  [Storable]


47  private bool symmetric;


48  private double[,] d2;


49 


50  [StorableConstructor]


51  protected CovarianceRQiso(bool deserializing)


52  : base(deserializing) {


53  }


54 


55  protected CovarianceRQiso(CovarianceRQiso original, Cloner cloner)


56  : base(original, cloner) {


57  if (original.x != null) {


58  this.x = new double[original.x.GetLength(0), original.x.GetLength(1)];


59  Array.Copy(original.x, this.x, x.Length);


60 


61  this.xt = new double[original.xt.GetLength(0), original.xt.GetLength(1)];


62  Array.Copy(original.xt, this.xt, xt.Length);


63 


64  this.d2 = new double[original.d2.GetLength(0), original.d2.GetLength(1)];


65  Array.Copy(original.d2, this.d2, d2.Length);


66  this.sf2 = original.sf2;


67  }


68  this.sf2 = original.sf2;


69  this.l = original.l;


70  this.alpha = original.alpha;


71  this.symmetric = original.symmetric;


72  }


73 


74  public CovarianceRQiso()


75  : base() {


76  }


77 


78  public override IDeepCloneable Clone(Cloner cloner) {


79  return new CovarianceRQiso(this, cloner);


80  }


81 


82  public int GetNumberOfParameters(int numberOfVariables) {


83  return 3;


84  }


85 


86  public void SetParameter(double[] hyp) {


87  this.l = Math.Exp(hyp[0]);


88  this.sf2 = Math.Exp(2 * hyp[1]);


89  this.alpha = Math.Exp(hyp[2]);


90  d2 = null;


91  }


92  public void SetData(double[,] x) {


93  SetData(x, x);


94  this.symmetric = true;


95  }


96 


97 


98  public void SetData(double[,] x, double[,] xt) {


99  this.symmetric = false;


100  this.x = x;


101  this.xt = xt;


102  d2 = null;


103  }


104 


105  public double GetCovariance(int i, int j) {


106  if (d2 == null) CalculateSquaredDistances();


107  return sf2 * Math.Pow(1 + 0.5 * d2[i, j] / alpha, alpha);


108  }


109 


110  public double GetGradient(int i, int j, int k) {


111  switch (k) {


112  case 0: return sf2 * Math.Pow(1 + 0.5 * d2[i, j] / alpha, alpha  1) * d2[i, j];


113  case 1: return 2 * sf2 * Math.Pow((1 + 0.5 * d2[i, j] / alpha), (alpha));


114  case 2: {


115  double g = (1 + 0.5 * d2[i, j] / alpha);


116  g = sf2 * Math.Pow(g, alpha) * (0.5 * d2[i, j] / g  alpha * Math.Log(g));


117  return g;


118  }


119  default: throw new ArgumentException("CovarianceRQiso has three hyperparameters", "k");


120  }


121  }


122 


123  private void CalculateSquaredDistances() {


124  if (x.GetLength(1) != xt.GetLength(1)) throw new InvalidOperationException();


125  int rows = x.GetLength(0);


126  int cols = xt.GetLength(0);


127  d2 = new double[rows, cols];


128  double lInv = 1.0 / l;


129  if (symmetric) {


130  for (int i = 0; i < rows; i++) {


131  for (int j = i; j < rows; j++) {


132  d2[i, j] = Util.SqrDist(Util.GetRow(x, i).Select(e => e * lInv), Util.GetRow(xt, j).Select(e => e * lInv));


133  d2[j, i] = d2[i, j];


134  }


135  }


136  } else {


137  for (int i = 0; i < rows; i++) {


138  for (int j = 0; j < cols; j++) {


139  d2[i, j] = Util.SqrDist(Util.GetRow(x, i).Select(e => e * lInv), Util.GetRow(xt, j).Select(e => e * lInv));


140  }


141  }


142  }


143  }


144  }


145  }

