source: trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GBM/GradientBoostingRegressionAlgorithm.cs @ 13721

Last change on this file since 13721 was 13721, checked in by mkommend, 5 years ago

#2591: Changed all GP covariance and mean functions to use int[] for column indices instead of IEnumerable<int>. Changed GP utils, GPModel and StudentTProcessModell as well to use fewer iterators and adapted unit tests to new interface.

File size: 21.5 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2015 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 * and the BEACON Center for the Study of Evolution in Action.
5 *
6 * This file is part of HeuristicLab.
7 *
8 * HeuristicLab is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * HeuristicLab is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
20 */
21#endregion
22
23using System;
24using System.Collections.Generic;
25using System.Linq;
26using System.Threading;
27using HeuristicLab.Analysis;
28using HeuristicLab.Common;
29using HeuristicLab.Core;
30using HeuristicLab.Data;
31using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
32using HeuristicLab.Optimization;
33using HeuristicLab.Parameters;
34using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
35using HeuristicLab.Problems.DataAnalysis;
36using HeuristicLab.Problems.DataAnalysis.Symbolic;
37using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;
38using HeuristicLab.Random;
39using HeuristicLab.Selection;
40
41namespace HeuristicLab.Algorithms.DataAnalysis.MctsSymbolicRegression {
42  [Item("Gradient Boosting Machine Regression (GBM)",
43    "Gradient boosting for any regression base learner (e.g. MCTS symbolic regression)")]
44  [StorableClass]
45  [Creatable(CreatableAttribute.Categories.DataAnalysisRegression, Priority = 350)]
46  public class GradientBoostingRegressionAlgorithm : BasicAlgorithm {
47    public override Type ProblemType {
48      get { return typeof(IRegressionProblem); }
49    }
50
51    public new IRegressionProblem Problem {
52      get { return (IRegressionProblem)base.Problem; }
53      set { base.Problem = value; }
54    }
55
56    #region ParameterNames
57
58    private const string IterationsParameterName = "Iterations";
59    private const string NuParameterName = "Nu";
60    private const string MParameterName = "M";
61    private const string RParameterName = "R";
62    private const string RegressionAlgorithmParameterName = "RegressionAlgorithm";
63    private const string SeedParameterName = "Seed";
64    private const string SetSeedRandomlyParameterName = "SetSeedRandomly";
65    private const string CreateSolutionParameterName = "CreateSolution";
66    private const string RegressionAlgorithmSolutionResultParameterName = "RegressionAlgorithmResult";
67
68    #endregion
69
70    #region ParameterProperties
71
72    public IFixedValueParameter<IntValue> IterationsParameter {
73      get { return (IFixedValueParameter<IntValue>)Parameters[IterationsParameterName]; }
74    }
75
76    public IFixedValueParameter<DoubleValue> NuParameter {
77      get { return (IFixedValueParameter<DoubleValue>)Parameters[NuParameterName]; }
78    }
79
80    public IFixedValueParameter<DoubleValue> RParameter {
81      get { return (IFixedValueParameter<DoubleValue>)Parameters[RParameterName]; }
82    }
83
84    public IFixedValueParameter<DoubleValue> MParameter {
85      get { return (IFixedValueParameter<DoubleValue>)Parameters[MParameterName]; }
86    }
87
88    // regression algorithms are currently: DataAnalysisAlgorithms, BasicAlgorithms and engine algorithms with no common interface
89    public IConstrainedValueParameter<IAlgorithm> RegressionAlgorithmParameter {
90      get { return (IConstrainedValueParameter<IAlgorithm>)Parameters[RegressionAlgorithmParameterName]; }
91    }
92
93    public IFixedValueParameter<StringValue> RegressionAlgorithmSolutionResultParameter {
94      get { return (IFixedValueParameter<StringValue>)Parameters[RegressionAlgorithmSolutionResultParameterName]; }
95    }
96
97    public IFixedValueParameter<IntValue> SeedParameter {
98      get { return (IFixedValueParameter<IntValue>)Parameters[SeedParameterName]; }
99    }
100
101    public FixedValueParameter<BoolValue> SetSeedRandomlyParameter {
102      get { return (FixedValueParameter<BoolValue>)Parameters[SetSeedRandomlyParameterName]; }
103    }
104
105    public IFixedValueParameter<BoolValue> CreateSolutionParameter {
106      get { return (IFixedValueParameter<BoolValue>)Parameters[CreateSolutionParameterName]; }
107    }
108
109    #endregion
110
111    #region Properties
112
113    public int Iterations {
114      get { return IterationsParameter.Value.Value; }
115      set { IterationsParameter.Value.Value = value; }
116    }
117
118    public int Seed {
119      get { return SeedParameter.Value.Value; }
120      set { SeedParameter.Value.Value = value; }
121    }
122
123    public bool SetSeedRandomly {
124      get { return SetSeedRandomlyParameter.Value.Value; }
125      set { SetSeedRandomlyParameter.Value.Value = value; }
126    }
127
128    public double Nu {
129      get { return NuParameter.Value.Value; }
130      set { NuParameter.Value.Value = value; }
131    }
132
133    public double R {
134      get { return RParameter.Value.Value; }
135      set { RParameter.Value.Value = value; }
136    }
137
138    public double M {
139      get { return MParameter.Value.Value; }
140      set { MParameter.Value.Value = value; }
141    }
142
143    public bool CreateSolution {
144      get { return CreateSolutionParameter.Value.Value; }
145      set { CreateSolutionParameter.Value.Value = value; }
146    }
147
148    public IAlgorithm RegressionAlgorithm {
149      get { return RegressionAlgorithmParameter.Value; }
150    }
151
152    public string RegressionAlgorithmResult {
153      get { return RegressionAlgorithmSolutionResultParameter.Value.Value; }
154      set { RegressionAlgorithmSolutionResultParameter.Value.Value = value; }
155    }
156
157    #endregion
158
159    [StorableConstructor]
160    protected GradientBoostingRegressionAlgorithm(bool deserializing)
161      : base(deserializing) {
162    }
163
164    protected GradientBoostingRegressionAlgorithm(GradientBoostingRegressionAlgorithm original, Cloner cloner)
165      : base(original, cloner) {
166    }
167
168    public override IDeepCloneable Clone(Cloner cloner) {
169      return new GradientBoostingRegressionAlgorithm(this, cloner);
170    }
171
172    public GradientBoostingRegressionAlgorithm() {
173      Problem = new RegressionProblem(); // default problem
174      var mctsSymbReg = new MctsSymbolicRegressionAlgorithm();
175      mctsSymbReg.Iterations = 10000;
176      mctsSymbReg.StoreAlgorithmInEachRun = false;
177      var sgp = CreateOSGP();
178      var regressionAlgs = new ItemSet<IAlgorithm>(new IAlgorithm[] {
179        new RandomForestRegression(),
180        sgp,
181        mctsSymbReg
182      });
183      foreach (var alg in regressionAlgs) alg.Prepare();
184
185
186      Parameters.Add(new FixedValueParameter<IntValue>(IterationsParameterName,
187        "Number of iterations", new IntValue(100)));
188      Parameters.Add(new FixedValueParameter<IntValue>(SeedParameterName,
189        "The random seed used to initialize the new pseudo random number generator.", new IntValue(0)));
190      Parameters.Add(new FixedValueParameter<BoolValue>(SetSeedRandomlyParameterName,
191        "True if the random seed should be set to a random value, otherwise false.", new BoolValue(true)));
192      Parameters.Add(new FixedValueParameter<DoubleValue>(NuParameterName,
193        "The learning rate nu when updating predictions in GBM (0 < nu <= 1)", new DoubleValue(0.5)));
194      Parameters.Add(new FixedValueParameter<DoubleValue>(RParameterName,
195        "The fraction of rows that are sampled randomly for the base learner in each iteration (0 < r <= 1)",
196        new DoubleValue(1)));
197      Parameters.Add(new FixedValueParameter<DoubleValue>(MParameterName,
198        "The fraction of variables that are sampled randomly for the base learner in each iteration (0 < m <= 1)",
199        new DoubleValue(0.5)));
200      Parameters.Add(new ConstrainedValueParameter<IAlgorithm>(RegressionAlgorithmParameterName,
201        "The regression algorithm to use as a base learner", regressionAlgs, mctsSymbReg));
202      Parameters.Add(new FixedValueParameter<StringValue>(RegressionAlgorithmSolutionResultParameterName,
203        "The name of the solution produced by the regression algorithm", new StringValue("Solution")));
204      Parameters[RegressionAlgorithmSolutionResultParameterName].Hidden = true;
205      Parameters.Add(new FixedValueParameter<BoolValue>(CreateSolutionParameterName,
206        "Flag that indicates if a solution should be produced at the end of the run", new BoolValue(true)));
207      Parameters[CreateSolutionParameterName].Hidden = true;
208    }
209
210    protected override void Run(CancellationToken cancellationToken) {
211      // Set up the algorithm
212      if (SetSeedRandomly) Seed = new System.Random().Next();
213      var rand = new MersenneTwister((uint)Seed);
214
215      // Set up the results display
216      var iterations = new IntValue(0);
217      Results.Add(new Result("Iterations", iterations));
218
219      var table = new DataTable("Qualities");
220      table.Rows.Add(new DataRow("Loss (train)"));
221      table.Rows.Add(new DataRow("Loss (test)"));
222      Results.Add(new Result("Qualities", table));
223      var curLoss = new DoubleValue();
224      var curTestLoss = new DoubleValue();
225      Results.Add(new Result("Loss (train)", curLoss));
226      Results.Add(new Result("Loss (test)", curTestLoss));
227      var runCollection = new RunCollection();
228      Results.Add(new Result("Runs", runCollection));
229
230      // init
231      var problemData = Problem.ProblemData;
232      var targetVarName = Problem.ProblemData.TargetVariable;
233      var activeVariables = problemData.AllowedInputVariables.Concat(new string[] { problemData.TargetVariable });
234      var modifiableDataset = new ModifiableDataset(
235        activeVariables,
236        activeVariables.Select(v => problemData.Dataset.GetDoubleValues(v).ToList()));
237
238      var trainingRows = problemData.TrainingIndices;
239      var testRows = problemData.TestIndices;
240      var yPred = new double[trainingRows.Count()];
241      var yPredTest = new double[testRows.Count()];
242      var y = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TrainingIndices).ToArray();
243      var curY = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TrainingIndices).ToArray();
244
245      var yTest = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TestIndices).ToArray();
246      var curYTest = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TestIndices).ToArray();
247      var nu = Nu;
248      var mVars = (int)Math.Ceiling(M * problemData.AllowedInputVariables.Count());
249      var rRows = (int)Math.Ceiling(R * problemData.TrainingIndices.Count());
250      var alg = RegressionAlgorithm;
251      List<IRegressionModel> models = new List<IRegressionModel>();
252      try {
253
254        // Loop until iteration limit reached or canceled.
255        for (int i = 0; i < Iterations; i++) {
256          cancellationToken.ThrowIfCancellationRequested();
257
258          modifiableDataset.RemoveVariable(targetVarName);
259          modifiableDataset.AddVariable(targetVarName, curY.Concat(curYTest));
260
261          SampleTrainingData(rand, modifiableDataset, rRows, problemData.Dataset, curY, problemData.TargetVariable, problemData.TrainingIndices); // all training indices from the original problem data are allowed
262          var modifiableProblemData = new RegressionProblemData(modifiableDataset,
263            problemData.AllowedInputVariables.SampleRandomWithoutRepetition(rand, mVars),
264            problemData.TargetVariable);
265          modifiableProblemData.TrainingPartition.Start = 0;
266          modifiableProblemData.TrainingPartition.End = rRows;
267          modifiableProblemData.TestPartition.Start = problemData.TestPartition.Start;
268          modifiableProblemData.TestPartition.End = problemData.TestPartition.End;
269
270          if (!TrySetProblemData(alg, modifiableProblemData))
271            throw new NotSupportedException("The algorithm cannot be used with GBM.");
272
273          IRegressionModel model;
274          IRun run;
275          // try to find a model. The algorithm might fail to produce a model. In this case we just retry until the iterations are exhausted
276          if (TryExecute(alg, RegressionAlgorithmResult, out model, out run)) {
277            int row = 0;
278            // update predictions for training and test
279            // update new targets (in the case of squared error loss we simply use negative residuals)
280            foreach (var pred in model.GetEstimatedValues(problemData.Dataset, trainingRows)) {
281              yPred[row] = yPred[row] + nu * pred;
282              curY[row] = y[row] - yPred[row];
283              row++;
284            }
285            row = 0;
286            foreach (var pred in model.GetEstimatedValues(problemData.Dataset, testRows)) {
287              yPredTest[row] = yPredTest[row] + nu * pred;
288              curYTest[row] = yTest[row] - yPredTest[row];
289              row++;
290            }
291            // determine quality
292            OnlineCalculatorError error;
293            var trainR = OnlinePearsonsRCalculator.Calculate(yPred, y, out error);
294            var testR = OnlinePearsonsRCalculator.Calculate(yPredTest, yTest, out error);
295
296            // iteration results
297            curLoss.Value = error == OnlineCalculatorError.None ? trainR * trainR : 0.0;
298            curTestLoss.Value = error == OnlineCalculatorError.None ? testR * testR : 0.0;
299
300            models.Add(model);
301
302
303          }
304
305          runCollection.Add(run);
306          table.Rows["Loss (train)"].Values.Add(curLoss.Value);
307          table.Rows["Loss (test)"].Values.Add(curTestLoss.Value);
308          iterations.Value = i + 1;
309        }
310
311        // produce solution
312        if (CreateSolution) {
313          // when all our models are symbolic models we can easily combine them to a single model
314          if (models.All(m => m is ISymbolicRegressionModel)) {
315            Results.Add(new Result("Solution", CreateSymbolicSolution(models, Nu, (IRegressionProblemData)problemData.Clone())));
316          }
317          // just produce an ensemble solution for now (TODO: correct scaling or linear regression for ensemble model weights)
318
319          var ensembleSolution = CreateEnsembleSolution(models, (IRegressionProblemData)problemData.Clone());
320          Results.Add(new Result("EnsembleSolution", ensembleSolution));
321        }
322      }
323      finally {
324        // reset everything
325        alg.Prepare(true);
326      }
327    }
328
329    private static IRegressionEnsembleSolution CreateEnsembleSolution(List<IRegressionModel> models,
330      IRegressionProblemData problemData) {
331      var rows = problemData.TrainingPartition.Size;
332      var features = models.Count;
333      double[,] inputMatrix = new double[rows, features + 1];
334
335      //add model estimates
336      for (int m = 0; m < models.Count; m++) {
337        var model = models[m];
338        var estimates = model.GetEstimatedValues(problemData.Dataset, problemData.TrainingIndices);
339        int estimatesCounter = 0;
340        foreach (var estimate in estimates) {
341          inputMatrix[estimatesCounter, m] = estimate;
342          estimatesCounter++;
343        }
344      }
345
346      //add target
347      var targets = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TrainingIndices);
348      int targetCounter = 0;
349      foreach (var target in targets) {
350        inputMatrix[targetCounter, models.Count] = target;
351        targetCounter++;
352      }
353
354      alglib.linearmodel lm = new alglib.linearmodel();
355      alglib.lrreport ar = new alglib.lrreport();
356      double[] coefficients;
357      int retVal = 1;
358      alglib.lrbuildz(inputMatrix, rows, features, out retVal, out lm, out ar);
359      if (retVal != 1) throw new ArgumentException("Error in calculation of linear regression solution");
360
361      alglib.lrunpack(lm, out coefficients, out features);
362
363      var ensembleModel = new RegressionEnsembleModel(models, coefficients.Take(models.Count)) { AverageModelEstimates = false };
364      var ensembleSolution = ensembleModel.CreateRegressionSolution((IRegressionProblemData)problemData.Clone());
365      return ensembleSolution;
366    }
367
368
369    private IAlgorithm CreateOSGP() {
370      // configure strict osgp
371      var alg = new OffspringSelectionGeneticAlgorithm.OffspringSelectionGeneticAlgorithm();
372      var prob = new SymbolicRegressionSingleObjectiveProblem();
373      prob.MaximumSymbolicExpressionTreeDepth.Value = 7;
374      prob.MaximumSymbolicExpressionTreeLength.Value = 15;
375      alg.Problem = prob;
376      alg.SuccessRatio.Value = 1.0;
377      alg.ComparisonFactorLowerBound.Value = 1.0;
378      alg.ComparisonFactorUpperBound.Value = 1.0;
379      alg.MutationProbability.Value = 0.15;
380      alg.PopulationSize.Value = 200;
381      alg.MaximumSelectionPressure.Value = 100;
382      alg.MaximumEvaluatedSolutions.Value = 20000;
383      alg.SelectorParameter.Value = alg.SelectorParameter.ValidValues.OfType<GenderSpecificSelector>().First();
384      alg.MutatorParameter.Value = alg.MutatorParameter.ValidValues.OfType<MultiSymbolicExpressionTreeManipulator>().First();
385      alg.StoreAlgorithmInEachRun = false;
386      return alg;
387    }
388
389    private void SampleTrainingData(MersenneTwister rand, ModifiableDataset ds, int rRows,
390      IDataset sourceDs, double[] curTarget, string targetVarName, IEnumerable<int> trainingIndices) {
391      var selectedRows = trainingIndices.SampleRandomWithoutRepetition(rand, rRows).ToArray();
392      int t = 0;
393      object[] srcRow = new object[ds.Columns];
394      var varNames = ds.DoubleVariables.ToArray();
395      foreach (var r in selectedRows) {
396        // take all values from the original dataset
397        for (int c = 0; c < srcRow.Length; c++) {
398          var col = sourceDs.GetReadOnlyDoubleValues(varNames[c]);
399          srcRow[c] = col[r];
400        }
401        ds.ReplaceRow(t, srcRow);
402        // but use the updated target values
403        ds.SetVariableValue(curTarget[r], targetVarName, t);
404        t++;
405      }
406    }
407
408    private static ISymbolicRegressionSolution CreateSymbolicSolution(List<IRegressionModel> models, double nu, IRegressionProblemData problemData) {
409      var symbModels = models.OfType<ISymbolicRegressionModel>();
410      var lowerLimit = symbModels.Min(m => m.LowerEstimationLimit);
411      var upperLimit = symbModels.Max(m => m.UpperEstimationLimit);
412      var interpreter = new SymbolicDataAnalysisExpressionTreeLinearInterpreter();
413      var progRootNode = new ProgramRootSymbol().CreateTreeNode();
414      var startNode = new StartSymbol().CreateTreeNode();
415
416      var addNode = new Addition().CreateTreeNode();
417      var mulNode = new Multiplication().CreateTreeNode();
418      var scaleNode = (ConstantTreeNode)new Constant().CreateTreeNode(); // all models are scaled using the same nu
419      scaleNode.Value = nu;
420
421      foreach (var m in symbModels) {
422        var relevantPart = m.SymbolicExpressionTree.Root.GetSubtree(0).GetSubtree(0); // skip root and start
423        addNode.AddSubtree((ISymbolicExpressionTreeNode)relevantPart.Clone());
424      }
425
426      mulNode.AddSubtree(addNode);
427      mulNode.AddSubtree(scaleNode);
428      startNode.AddSubtree(mulNode);
429      progRootNode.AddSubtree(startNode);
430      var t = new SymbolicExpressionTree(progRootNode);
431      var combinedModel = new SymbolicRegressionModel(t, interpreter, lowerLimit, upperLimit);
432      var sol = new SymbolicRegressionSolution(combinedModel, problemData);
433      return sol;
434    }
435
436    private static bool TrySetProblemData(IAlgorithm alg, IRegressionProblemData problemData) {
437      var prob = alg.Problem as IRegressionProblem;
438      // there is already a problem and it is compatible -> just set problem data
439      if (prob != null) {
440        prob.ProblemDataParameter.Value = problemData;
441        return true;
442      } else return false;
443    }
444
445    private static bool TryExecute(IAlgorithm alg, string regressionAlgorithmResultName, out IRegressionModel model, out IRun run) {
446      model = null;
447      using (var wh = new AutoResetEvent(false)) {
448        EventHandler<EventArgs<Exception>> handler = (sender, args) => wh.Set();
449        EventHandler handler2 = (sender, args) => wh.Set();
450        alg.ExceptionOccurred += handler;
451        alg.Stopped += handler2;
452        try {
453          alg.Prepare();
454          alg.Start();
455          wh.WaitOne();
456
457          run = alg.Runs.Last();
458          var sols = alg.Results.Select(r => r.Value).OfType<IRegressionSolution>();
459          if (!sols.Any()) return false;
460          var sol = sols.First();
461          if (sols.Skip(1).Any()) {
462            // more than one solution => use regressionAlgorithmResult
463            if (alg.Results.ContainsKey(regressionAlgorithmResultName)) {
464              sol = (IRegressionSolution)alg.Results[regressionAlgorithmResultName].Value;
465            }
466          }
467          var symbRegSol = sol as SymbolicRegressionSolution;
468          // only accept symb reg solutions that do not hit the estimation limits
469          // NaN evaluations would not be critical but are problematic if we want to combine all symbolic models into a single symbolic model
470          if (symbRegSol == null ||
471            (symbRegSol.TrainingLowerEstimationLimitHits == 0 && symbRegSol.TrainingUpperEstimationLimitHits == 0 &&
472             symbRegSol.TestLowerEstimationLimitHits == 0 && symbRegSol.TestUpperEstimationLimitHits == 0) &&
473            symbRegSol.TrainingNaNEvaluations == 0 && symbRegSol.TestNaNEvaluations == 0) {
474            model = sol.Model;
475          }
476        }
477        finally {
478          alg.ExceptionOccurred -= handler;
479          alg.Stopped -= handler2;
480        }
481      }
482      return model != null;
483    }
484  }
485}
Note: See TracBrowser for help on using the repository browser.