Index: /trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessCovarianceOptimizationProblem.cs
===================================================================
--- /trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessCovarianceOptimizationProblem.cs (revision 12946)
+++ /trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessCovarianceOptimizationProblem.cs (revision 12946)
@@ -0,0 +1,346 @@
+#region License Information
+/* HeuristicLab
+ * Copyright (C) 2002-2015 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 .
+ */
+#endregion
+
+using System;
+using System.Collections.Generic;
+using System.Diagnostics.Contracts;
+using System.Drawing.Text;
+using System.Linq;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Data;
+using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
+using HeuristicLab.Optimization;
+using HeuristicLab.Parameters;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+using HeuristicLab.Problems.DataAnalysis;
+using HeuristicLab.Problems.Instances;
+
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+ [Item("Gaussian Process Covariance Optimization Problem", "")]
+ [Creatable(CreatableAttribute.Categories.GeneticProgrammingProblems, Priority = 300)]
+ [StorableClass]
+ public sealed class GaussianProcessCovarianceOptimizationProblem : SymbolicExpressionTreeProblem, IRegressionProblem, IProblemInstanceConsumer, IProblemInstanceExporter {
+ #region static variables and ctor
+ private static readonly CovarianceMaternIso maternIso1;
+ private static readonly CovarianceMaternIso maternIso3;
+ private static readonly CovarianceMaternIso maternIso5;
+ private static readonly CovariancePiecewisePolynomial piecewisePoly0;
+ private static readonly CovariancePiecewisePolynomial piecewisePoly1;
+ private static readonly CovariancePiecewisePolynomial piecewisePoly2;
+ private static readonly CovariancePiecewisePolynomial piecewisePoly3;
+ private static readonly CovariancePolynomial poly2;
+ private static readonly CovariancePolynomial poly3;
+ private static readonly CovarianceSpectralMixture spectralMixture1;
+ private static readonly CovarianceSpectralMixture spectralMixture3;
+ private static readonly CovarianceSpectralMixture spectralMixture5;
+ private static readonly CovarianceLinear linear;
+ private static readonly CovarianceLinearArd linearArd;
+ private static readonly CovarianceNeuralNetwork neuralNetwork;
+ private static readonly CovariancePeriodic periodic;
+ private static readonly CovarianceRationalQuadraticIso ratQuadraticIso;
+ private static readonly CovarianceRationalQuadraticArd ratQuadraticArd;
+ private static readonly CovarianceSquaredExponentialArd sqrExpArd;
+ private static readonly CovarianceSquaredExponentialIso sqrExpIso;
+
+ static GaussianProcessCovarianceOptimizationProblem() {
+ // cumbersome initialization because of ConstrainedValueParameters
+ maternIso1 = new CovarianceMaternIso(); SetConstrainedValueParameter(maternIso1.DParameter, 1);
+ maternIso3 = new CovarianceMaternIso(); SetConstrainedValueParameter(maternIso3.DParameter, 3);
+ maternIso5 = new CovarianceMaternIso(); SetConstrainedValueParameter(maternIso5.DParameter, 5);
+
+ piecewisePoly0 = new CovariancePiecewisePolynomial(); SetConstrainedValueParameter(piecewisePoly0.VParameter, 0);
+ piecewisePoly1 = new CovariancePiecewisePolynomial(); SetConstrainedValueParameter(piecewisePoly1.VParameter, 1);
+ piecewisePoly2 = new CovariancePiecewisePolynomial(); SetConstrainedValueParameter(piecewisePoly2.VParameter, 2);
+ piecewisePoly3 = new CovariancePiecewisePolynomial(); SetConstrainedValueParameter(piecewisePoly3.VParameter, 3);
+
+ poly2 = new CovariancePolynomial(); poly2.DegreeParameter.Value.Value = 2;
+ poly3 = new CovariancePolynomial(); poly3.DegreeParameter.Value.Value = 3;
+
+ spectralMixture1 = new CovarianceSpectralMixture(); spectralMixture1.QParameter.Value.Value = 1;
+ spectralMixture3 = new CovarianceSpectralMixture(); spectralMixture3.QParameter.Value.Value = 3;
+ spectralMixture5 = new CovarianceSpectralMixture(); spectralMixture5.QParameter.Value.Value = 5;
+
+ linear = new CovarianceLinear();
+ linearArd = new CovarianceLinearArd();
+ neuralNetwork = new CovarianceNeuralNetwork();
+ periodic = new CovariancePeriodic();
+ ratQuadraticArd = new CovarianceRationalQuadraticArd();
+ ratQuadraticIso = new CovarianceRationalQuadraticIso();
+ sqrExpArd = new CovarianceSquaredExponentialArd();
+ sqrExpIso = new CovarianceSquaredExponentialIso();
+ }
+
+ private static void SetConstrainedValueParameter(IConstrainedValueParameter param, int val) {
+ param.Value = param.ValidValues.Single(v => v.Value == val);
+ }
+
+ #endregion
+
+
+ #region parameter names
+
+ private const string ProblemDataParameterName = "ProblemData";
+ private const string ConstantOptIterationsParameterName = "Constant optimization steps";
+ private const string RestartsParameterName = "Restarts";
+
+ #endregion
+
+ #region Parameter Properties
+ IParameter IDataAnalysisProblem.ProblemDataParameter { get { return ProblemDataParameter; } }
+
+ public IValueParameter ProblemDataParameter {
+ get { return (IValueParameter)Parameters[ProblemDataParameterName]; }
+ }
+ public IFixedValueParameter ConstantOptIterationsParameter {
+ get { return (IFixedValueParameter)Parameters[ConstantOptIterationsParameterName]; }
+ }
+ public IFixedValueParameter RestartsParameter {
+ get { return (IFixedValueParameter)Parameters[RestartsParameterName]; }
+ }
+ #endregion
+
+ #region Properties
+
+ public IRegressionProblemData ProblemData {
+ get { return ProblemDataParameter.Value; }
+ set { ProblemDataParameter.Value = value; }
+ }
+ IDataAnalysisProblemData IDataAnalysisProblem.ProblemData { get { return ProblemData; } }
+
+ public int ConstantOptIterations {
+ get { return ConstantOptIterationsParameter.Value.Value; }
+ set { ConstantOptIterationsParameter.Value.Value = value; }
+ }
+
+ public int Restarts {
+ get { return RestartsParameter.Value.Value; }
+ set { RestartsParameter.Value.Value = value; }
+ }
+ #endregion
+
+ public override bool Maximization {
+ get { return true; } // return log likelihood (instead of negative log likelihood as in GPR
+ }
+
+ public GaussianProcessCovarianceOptimizationProblem()
+ : base() {
+ Parameters.Add(new ValueParameter(ProblemDataParameterName, "The data for the regression problem", new RegressionProblemData()));
+ Parameters.Add(new FixedValueParameter(ConstantOptIterationsParameterName, "Number of optimization steps for hyperparameter values", new IntValue(50)));
+ Parameters.Add(new FixedValueParameter(RestartsParameterName, "The number of random restarts for constant optimization.", new IntValue(10)));
+ Parameters["Restarts"].Hidden = true;
+ var g = new SimpleSymbolicExpressionGrammar();
+ g.AddSymbols(new string[] { "Sum", "Product" }, 2, 2);
+ g.AddTerminalSymbols(new string[]
+ {
+ "Linear",
+ "LinearArd",
+ "MaternIso1",
+ "MaternIso3",
+ "MaternIso5",
+ "NeuralNetwork",
+ "Periodic",
+ "PiecewisePolynomial0",
+ "PiecewisePolynomial1",
+ "PiecewisePolynomial2",
+ "PiecewisePolynomial3",
+ "Polynomial2",
+ "Polynomial3",
+ "RationalQuadraticArd",
+ "RationalQuadraticIso",
+ "SpectralMixture1",
+ "SpectralMixture3",
+ "SpectralMixture5",
+ "SquaredExponentialArd",
+ "SquaredExponentialIso"
+ });
+ base.Encoding = new SymbolicExpressionTreeEncoding(g, 10, 5);
+ }
+
+
+ public override double Evaluate(ISymbolicExpressionTree tree, IRandom random) {
+
+ var meanFunction = new MeanConst();
+ var problemData = ProblemData;
+ var ds = problemData.Dataset;
+ var targetVariable = problemData.TargetVariable;
+ var allowedInputVariables = problemData.AllowedInputVariables.ToArray();
+ var nVars = allowedInputVariables.Length;
+ var trainingRows = problemData.TrainingIndices.ToArray();
+
+ // use the same covariance function for each restart
+ var covarianceFunction = TreeToCovarianceFunction(tree);
+
+ // allocate hyperparameters
+ var hyperParameters = new double[meanFunction.GetNumberOfParameters(nVars) + covarianceFunction.GetNumberOfParameters(nVars) + 1]; // mean + cov + noise
+ double[] bestHyperParameters = new double[hyperParameters.Length];
+ var bestObjValue = new double[1] { double.MinValue };
+
+ // data that is necessary for the objective function
+ var data = Tuple.Create(ds, targetVariable, allowedInputVariables, trainingRows, (IMeanFunction)meanFunction, covarianceFunction, bestObjValue);
+
+ for (int t = 0; t < Restarts; t++) {
+ var prevBest = bestObjValue[0];
+ var prevBestHyperParameters = new double[hyperParameters.Length];
+ Array.Copy(bestHyperParameters, prevBestHyperParameters, bestHyperParameters.Length);
+
+ // initialize hyperparameters
+ hyperParameters[0] = ds.GetDoubleValues(targetVariable).Average(); // mean const
+
+ for (int i = 0; i < covarianceFunction.GetNumberOfParameters(nVars); i++) {
+ hyperParameters[1 + i] = random.NextDouble() * 2.0 - 1.0;
+ }
+ hyperParameters[hyperParameters.Length - 1] = 1.0; // s² = exp(2), TODO: other inits better?
+
+ // use alglib.bfgs for hyper-parameter optimization ...
+ double epsg = 0;
+ double epsf = 0.00001;
+ double epsx = 0;
+ double stpmax = 1;
+ int maxits = ConstantOptIterations;
+ alglib.mincgstate state;
+ alglib.mincgreport rep;
+
+ alglib.mincgcreate(hyperParameters, out state);
+ alglib.mincgsetcond(state, epsg, epsf, epsx, maxits);
+ alglib.mincgsetstpmax(state, stpmax);
+ alglib.mincgoptimize(state, ObjectiveFunction, null, data);
+
+ alglib.mincgresults(state, out bestHyperParameters, out rep);
+
+ if (rep.terminationtype < 0) {
+ // error -> restore previous best quality
+ bestObjValue[0] = prevBest;
+ Array.Copy(prevBestHyperParameters, bestHyperParameters, prevBestHyperParameters.Length);
+ }
+ }
+
+ return bestObjValue[0];
+ }
+
+ public void ObjectiveFunction(double[] x, ref double func, double[] grad, object obj) {
+ // we want to optimize the model likelihood by changing the hyperparameters and also return the gradient for each hyperparameter
+ var data = (Tuple)obj;
+ var ds = data.Item1;
+ var targetVariable = data.Item2;
+ var allowedInputVariables = data.Item3;
+ var trainingRows = data.Item4;
+ var meanFunction = data.Item5;
+ var covarianceFunction = data.Item6;
+ var bestObjValue = data.Item7;
+ var hyperParameters = x; // the decision variable vector
+
+ try {
+ var model = new GaussianProcessModel(ds, targetVariable, allowedInputVariables, trainingRows, hyperParameters, meanFunction, covarianceFunction);
+
+ func = model.NegativeLogLikelihood; // mincgoptimize, so we return negative likelihood
+ bestObjValue[0] = Math.Max(bestObjValue[0], -func); // problem itself is a maximization problem
+ var gradients = model.HyperparameterGradients;
+ Array.Copy(gradients, grad, gradients.Length);
+ } catch (Exception) {
+ // building the GaussianProcessModel might fail, in this case we return the worst possible objective value
+ func = 1.0E+300;
+ Array.Clear(grad, 0, grad.Length);
+ }
+ }
+
+ private ICovarianceFunction TreeToCovarianceFunction(ISymbolicExpressionTree tree) {
+ return TreeToCovarianceFunction(tree.Root.GetSubtree(0).GetSubtree(0)); // skip programroot and startsymbol
+ }
+
+ private ICovarianceFunction TreeToCovarianceFunction(ISymbolicExpressionTreeNode node) {
+ switch (node.Symbol.Name) {
+ case "Sum": {
+ var sum = new CovarianceSum();
+ sum.Terms.Add(TreeToCovarianceFunction(node.GetSubtree(0)));
+ sum.Terms.Add(TreeToCovarianceFunction(node.GetSubtree(1)));
+ return sum;
+ }
+ case "Product": {
+ var prod = new CovarianceProduct();
+ prod.Factors.Add(TreeToCovarianceFunction(node.GetSubtree(0)));
+ prod.Factors.Add(TreeToCovarianceFunction(node.GetSubtree(1)));
+ return prod;
+ }
+ // covFunction is cloned by the model so we can reuse instances of terminal covariance functions
+ case "Linear": return linear;
+ case "LinearArd": return linearArd;
+ case "MaternIso1": return maternIso1;
+ case "MaternIso3": return maternIso3;
+ case "MaternIso5": return maternIso5;
+ case "NeuralNetwork": return neuralNetwork;
+ case "Periodic": return periodic;
+ case "PiecewisePolynomial0": return piecewisePoly0;
+ case "PiecewisePolynomial1": return piecewisePoly1;
+ case "PiecewisePolynomial2": return piecewisePoly2;
+ case "PiecewisePolynomial3": return piecewisePoly3;
+ case "Polynomial2": return poly2;
+ case "Polynomial3": return poly3;
+ case "RationalQuadraticArd": return ratQuadraticArd;
+ case "RationalQuadraticIso": return ratQuadraticIso;
+ case "SpectralMixture1": return spectralMixture1;
+ case "SpectralMixture3": return spectralMixture3;
+ case "SpectralMixture5": return spectralMixture5;
+ case "SquaredExponentialArd": return sqrExpArd;
+ case "SquaredExponentialIso": return sqrExpIso;
+ default: throw new InvalidProgramException(string.Format("Found invalid symbol {0}", node.Symbol.Name));
+ }
+ }
+
+
+ // persistence
+ [StorableConstructor]
+ private GaussianProcessCovarianceOptimizationProblem(bool deserializing) : base(deserializing) { }
+ [StorableHook(HookType.AfterDeserialization)]
+ private void AfterDeserialization() {
+ }
+
+ // cloning
+ private GaussianProcessCovarianceOptimizationProblem(GaussianProcessCovarianceOptimizationProblem original, Cloner cloner)
+ : base(original, cloner) {
+ }
+ public override IDeepCloneable Clone(Cloner cloner) {
+ return new GaussianProcessCovarianceOptimizationProblem(this, cloner);
+ }
+
+ public void Load(IRegressionProblemData data) {
+ this.ProblemData = data;
+ OnProblemDataChanged();
+ }
+
+ public IRegressionProblemData Export() {
+ return ProblemData;
+ }
+
+ #region events
+ public event EventHandler ProblemDataChanged;
+
+
+ private void OnProblemDataChanged() {
+ var handler = ProblemDataChanged;
+ if (handler != null)
+ handler(this, EventArgs.Empty);
+ }
+ #endregion
+
+ }
+}
Index: /trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/HeuristicLab.Algorithms.DataAnalysis-3.4.csproj
===================================================================
--- /trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/HeuristicLab.Algorithms.DataAnalysis-3.4.csproj (revision 12945)
+++ /trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/HeuristicLab.Algorithms.DataAnalysis-3.4.csproj (revision 12946)
@@ -189,4 +189,5 @@
Code
+