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


28  using HeuristicLab.Parameters;


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


30 


31  namespace HeuristicLab.Algorithms.DataAnalysis {


32  [StorableClass]


33  [Item(Name = "CovarianceRationalQuadraticArd",


34  Description = "Rational quadratic covariance function with automatic relevance determination for Gaussian processes.")]


35  public sealed class CovarianceRationalQuadraticArd : ParameterizedNamedItem, ICovarianceFunction {


36  public IValueParameter<DoubleValue> ScaleParameter {


37  get { return (IValueParameter<DoubleValue>)Parameters["Scale"]; }


38  }


39 


40  public IValueParameter<DoubleArray> InverseLengthParameter {


41  get { return (IValueParameter<DoubleArray>)Parameters["InverseLength"]; }


42  }


43 


44  public IValueParameter<DoubleValue> ShapeParameter {


45  get { return (IValueParameter<DoubleValue>)Parameters["Shape"]; }


46  }


47 


48  [StorableConstructor]


49  private CovarianceRationalQuadraticArd(bool deserializing)


50  : base(deserializing) {


51  }


52 


53  private CovarianceRationalQuadraticArd(CovarianceRationalQuadraticArd original, Cloner cloner)


54  : base(original, cloner) {


55  }


56 


57  public CovarianceRationalQuadraticArd()


58  : base() {


59  Name = ItemName;


60  Description = ItemDescription;


61 


62  Parameters.Add(new OptionalValueParameter<DoubleValue>("Scale", "The scale parameter of the rational quadratic covariance function with ARD."));


63  Parameters.Add(new OptionalValueParameter<DoubleArray>("InverseLength", "The inverse length parameter for automatic relevance determination."));


64  Parameters.Add(new OptionalValueParameter<DoubleValue>("Shape", "The shape parameter (alpha) of the rational quadratic covariance function with ARD."));


65  }


66 


67  public override IDeepCloneable Clone(Cloner cloner) {


68  return new CovarianceRationalQuadraticArd(this, cloner);


69  }


70 


71  public int GetNumberOfParameters(int numberOfVariables) {


72  return


73  (ScaleParameter.Value != null ? 0 : 1) +


74  (ShapeParameter.Value != null ? 0 : 1) +


75  (InverseLengthParameter.Value != null ? 0 : numberOfVariables);


76  }


77 


78  public void SetParameter(double[] p) {


79  double scale, shape;


80  double[] inverseLength;


81  GetParameterValues(p, out scale, out shape, out inverseLength);


82  ScaleParameter.Value = new DoubleValue(scale);


83  ShapeParameter.Value = new DoubleValue(shape);


84  InverseLengthParameter.Value = new DoubleArray(inverseLength);


85  }


86 


87  private void GetParameterValues(double[] p, out double scale, out double shape, out double[] inverseLength) {


88  int c = 0;


89  // gather parameter values


90  if (ScaleParameter.Value != null) {


91  scale = ScaleParameter.Value.Value;


92  } else {


93  scale = Math.Exp(2 * p[c]);


94  c++;


95  }


96  if (ShapeParameter.Value != null) {


97  shape = ShapeParameter.Value.Value;


98  } else {


99  shape = Math.Exp(p[c]);


100  c++;


101  }


102  if (InverseLengthParameter.Value != null) {


103  inverseLength = InverseLengthParameter.Value.ToArray();


104  } else {


105  inverseLength = p.Skip(2).Select(e => 1.0 / Math.Exp(e)).ToArray();


106  c += inverseLength.Length;


107  }


108  if (p.Length != c) throw new ArgumentException("The length of the parameter vector does not match the number of free parameters for CovarianceRationalQuadraticArd", "p");


109  }


110 


111  public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {


112  double scale, shape;


113  double[] inverseLength;


114  GetParameterValues(p, out scale, out shape, out inverseLength);


115  // create functions


116  var cov = new ParameterizedCovarianceFunction();


117  cov.Covariance = (x, i, j) => {


118  double d = i == j


119  ? 0.0


120  : Util.SqrDist(x, i, j, inverseLength, columnIndices);


121  return scale * Math.Pow(1 + 0.5 * d / shape, shape);


122  };


123  cov.CrossCovariance = (x, xt, i, j) => {


124  double d = Util.SqrDist(x, i, xt, j, inverseLength, columnIndices);


125  return scale * Math.Pow(1 + 0.5 * d / shape, shape);


126  };


127  cov.CovarianceGradient = (x, i, j) => GetGradient(x, i, j, columnIndices, scale, shape, inverseLength);


128  return cov;


129  }


130 


131  private static IEnumerable<double> GetGradient(double[,] x, int i, int j, IEnumerable<int> columnIndices, double scale, double shape, double[] inverseLength) {


132  if (columnIndices == null) columnIndices = Enumerable.Range(0, x.GetLength(1));


133  double d = i == j


134  ? 0.0


135  : Util.SqrDist(x, i, j, inverseLength, columnIndices);


136  double b = 1 + 0.5 * d / shape;


137  int k = 0;


138  foreach (var columnIndex in columnIndices) {


139  yield return scale * Math.Pow(b, shape  1) * Util.SqrDist(x[i, columnIndex] * inverseLength[k], x[j, columnIndex] * inverseLength[k]);


140  k++;


141  }


142  yield return 2 * scale * Math.Pow(b, shape);


143  yield return scale * Math.Pow(b, shape) * (0.5 * d / b  shape * Math.Log(b));


144  }


145  }


146  }

