source: branches/HeuristicLab.RegressionSolutionGradientView/HeuristicLab.Algorithms.DataAnalysis/3.4/GBM/GradientBoostingRegressionAlgorithm.cs @ 13948

Last change on this file since 13948 was 13948, checked in by pfleck, 5 years ago

#2597

  • Merged recent trunk changes.
  • Adapted VariablesUsedForPrediction property for RegressionSolutionTargetResponseGradientView.
  • Fixed a reference (.dll to project ref).
File size: 23.0 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 StoreRunsParameterName = "StoreRuns";
67    private const string RegressionAlgorithmSolutionResultParameterName = "RegressionAlgorithmResult";
68
69    #endregion
70
71    #region ParameterProperties
72
73    public IFixedValueParameter<IntValue> IterationsParameter {
74      get { return (IFixedValueParameter<IntValue>)Parameters[IterationsParameterName]; }
75    }
76
77    public IFixedValueParameter<DoubleValue> NuParameter {
78      get { return (IFixedValueParameter<DoubleValue>)Parameters[NuParameterName]; }
79    }
80
81    public IFixedValueParameter<DoubleValue> RParameter {
82      get { return (IFixedValueParameter<DoubleValue>)Parameters[RParameterName]; }
83    }
84
85    public IFixedValueParameter<DoubleValue> MParameter {
86      get { return (IFixedValueParameter<DoubleValue>)Parameters[MParameterName]; }
87    }
88
89    // regression algorithms are currently: DataAnalysisAlgorithms, BasicAlgorithms and engine algorithms with no common interface
90    public IConstrainedValueParameter<IAlgorithm> RegressionAlgorithmParameter {
91      get { return (IConstrainedValueParameter<IAlgorithm>)Parameters[RegressionAlgorithmParameterName]; }
92    }
93
94    public IFixedValueParameter<StringValue> RegressionAlgorithmSolutionResultParameter {
95      get { return (IFixedValueParameter<StringValue>)Parameters[RegressionAlgorithmSolutionResultParameterName]; }
96    }
97
98    public IFixedValueParameter<IntValue> SeedParameter {
99      get { return (IFixedValueParameter<IntValue>)Parameters[SeedParameterName]; }
100    }
101
102    public FixedValueParameter<BoolValue> SetSeedRandomlyParameter {
103      get { return (FixedValueParameter<BoolValue>)Parameters[SetSeedRandomlyParameterName]; }
104    }
105
106    public IFixedValueParameter<BoolValue> CreateSolutionParameter {
107      get { return (IFixedValueParameter<BoolValue>)Parameters[CreateSolutionParameterName]; }
108    }
109    public IFixedValueParameter<BoolValue> StoreRunsParameter {
110      get { return (IFixedValueParameter<BoolValue>)Parameters[StoreRunsParameterName]; }
111    }
112
113    #endregion
114
115    #region Properties
116
117    public int Iterations {
118      get { return IterationsParameter.Value.Value; }
119      set { IterationsParameter.Value.Value = value; }
120    }
121
122    public int Seed {
123      get { return SeedParameter.Value.Value; }
124      set { SeedParameter.Value.Value = value; }
125    }
126
127    public bool SetSeedRandomly {
128      get { return SetSeedRandomlyParameter.Value.Value; }
129      set { SetSeedRandomlyParameter.Value.Value = value; }
130    }
131
132    public double Nu {
133      get { return NuParameter.Value.Value; }
134      set { NuParameter.Value.Value = value; }
135    }
136
137    public double R {
138      get { return RParameter.Value.Value; }
139      set { RParameter.Value.Value = value; }
140    }
141
142    public double M {
143      get { return MParameter.Value.Value; }
144      set { MParameter.Value.Value = value; }
145    }
146
147    public bool CreateSolution {
148      get { return CreateSolutionParameter.Value.Value; }
149      set { CreateSolutionParameter.Value.Value = value; }
150    }
151
152    public bool StoreRuns {
153      get { return StoreRunsParameter.Value.Value; }
154      set { StoreRunsParameter.Value.Value = value; }
155    }
156
157    public IAlgorithm RegressionAlgorithm {
158      get { return RegressionAlgorithmParameter.Value; }
159    }
160
161    public string RegressionAlgorithmResult {
162      get { return RegressionAlgorithmSolutionResultParameter.Value.Value; }
163      set { RegressionAlgorithmSolutionResultParameter.Value.Value = value; }
164    }
165
166    #endregion
167
168    [StorableConstructor]
169    protected GradientBoostingRegressionAlgorithm(bool deserializing)
170      : base(deserializing) {
171    }
172
173    protected GradientBoostingRegressionAlgorithm(GradientBoostingRegressionAlgorithm original, Cloner cloner)
174      : base(original, cloner) {
175    }
176
177    public override IDeepCloneable Clone(Cloner cloner) {
178      return new GradientBoostingRegressionAlgorithm(this, cloner);
179    }
180
181    public GradientBoostingRegressionAlgorithm() {
182      Problem = new RegressionProblem(); // default problem
183      var mctsSymbReg = new MctsSymbolicRegressionAlgorithm();
184      mctsSymbReg.Iterations = 10000;
185      mctsSymbReg.StoreAlgorithmInEachRun = false;
186      var sgp = CreateOSGP();
187      var regressionAlgs = new ItemSet<IAlgorithm>(new IAlgorithm[] {
188        new RandomForestRegression(),
189        sgp,
190        mctsSymbReg
191      });
192      foreach (var alg in regressionAlgs) alg.Prepare();
193
194
195      Parameters.Add(new FixedValueParameter<IntValue>(IterationsParameterName,
196        "Number of iterations", new IntValue(100)));
197      Parameters.Add(new FixedValueParameter<IntValue>(SeedParameterName,
198        "The random seed used to initialize the new pseudo random number generator.", new IntValue(0)));
199      Parameters.Add(new FixedValueParameter<BoolValue>(SetSeedRandomlyParameterName,
200        "True if the random seed should be set to a random value, otherwise false.", new BoolValue(true)));
201      Parameters.Add(new FixedValueParameter<DoubleValue>(NuParameterName,
202        "The learning rate nu when updating predictions in GBM (0 < nu <= 1)", new DoubleValue(0.5)));
203      Parameters.Add(new FixedValueParameter<DoubleValue>(RParameterName,
204        "The fraction of rows that are sampled randomly for the base learner in each iteration (0 < r <= 1)",
205        new DoubleValue(1)));
206      Parameters.Add(new FixedValueParameter<DoubleValue>(MParameterName,
207        "The fraction of variables that are sampled randomly for the base learner in each iteration (0 < m <= 1)",
208        new DoubleValue(0.5)));
209      Parameters.Add(new ConstrainedValueParameter<IAlgorithm>(RegressionAlgorithmParameterName,
210        "The regression algorithm to use as a base learner", regressionAlgs, mctsSymbReg));
211      Parameters.Add(new FixedValueParameter<StringValue>(RegressionAlgorithmSolutionResultParameterName,
212        "The name of the solution produced by the regression algorithm", new StringValue("Solution")));
213      Parameters[RegressionAlgorithmSolutionResultParameterName].Hidden = true;
214      Parameters.Add(new FixedValueParameter<BoolValue>(CreateSolutionParameterName,
215        "Flag that indicates if a solution should be produced at the end of the run", new BoolValue(true)));
216      Parameters[CreateSolutionParameterName].Hidden = true;
217      Parameters.Add(new FixedValueParameter<BoolValue>(StoreRunsParameterName,
218        "Flag that indicates if the results of the individual runs should be stored for detailed analysis", new BoolValue(false)));
219      Parameters[StoreRunsParameterName].Hidden = true;
220    }
221
222    protected override void Run(CancellationToken cancellationToken) {
223      // Set up the algorithm
224      if (SetSeedRandomly) Seed = new System.Random().Next();
225      var rand = new MersenneTwister((uint)Seed);
226
227      // Set up the results display
228      var iterations = new IntValue(0);
229      Results.Add(new Result("Iterations", iterations));
230
231      var table = new DataTable("Qualities");
232      table.Rows.Add(new DataRow("R² (train)"));
233      table.Rows.Add(new DataRow("R² (test)"));
234      Results.Add(new Result("Qualities", table));
235      var curLoss = new DoubleValue();
236      var curTestLoss = new DoubleValue();
237      Results.Add(new Result("R² (train)", curLoss));
238      Results.Add(new Result("R² (test)", curTestLoss));
239      var runCollection = new RunCollection();
240      if (StoreRuns)
241        Results.Add(new Result("Runs", runCollection));
242
243      // init
244      var problemData = Problem.ProblemData;
245      var targetVarName = problemData.TargetVariable;
246      var activeVariables = problemData.AllowedInputVariables.Concat(new string[] { problemData.TargetVariable });
247      var modifiableDataset = new ModifiableDataset(
248        activeVariables,
249        activeVariables.Select(v => problemData.Dataset.GetDoubleValues(v).ToList()));
250
251      var trainingRows = problemData.TrainingIndices;
252      var testRows = problemData.TestIndices;
253      var yPred = new double[trainingRows.Count()];
254      var yPredTest = new double[testRows.Count()];
255      var y = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TrainingIndices).ToArray();
256      var curY = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TrainingIndices).ToArray();
257
258      var yTest = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TestIndices).ToArray();
259      var curYTest = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TestIndices).ToArray();
260      var nu = Nu;
261      var mVars = (int)Math.Ceiling(M * problemData.AllowedInputVariables.Count());
262      var rRows = (int)Math.Ceiling(R * problemData.TrainingIndices.Count());
263      var alg = RegressionAlgorithm;
264      List<IRegressionModel> models = new List<IRegressionModel>();
265      try {
266
267        // Loop until iteration limit reached or canceled.
268        for (int i = 0; i < Iterations; i++) {
269          cancellationToken.ThrowIfCancellationRequested();
270
271          modifiableDataset.RemoveVariable(targetVarName);
272          modifiableDataset.AddVariable(targetVarName, curY.Concat(curYTest));
273
274          SampleTrainingData(rand, modifiableDataset, rRows, problemData.Dataset, curY, problemData.TargetVariable, problemData.TrainingIndices); // all training indices from the original problem data are allowed
275          var modifiableProblemData = new RegressionProblemData(modifiableDataset,
276            problemData.AllowedInputVariables.SampleRandomWithoutRepetition(rand, mVars),
277            problemData.TargetVariable);
278          modifiableProblemData.TrainingPartition.Start = 0;
279          modifiableProblemData.TrainingPartition.End = rRows;
280          modifiableProblemData.TestPartition.Start = problemData.TestPartition.Start;
281          modifiableProblemData.TestPartition.End = problemData.TestPartition.End;
282
283          if (!TrySetProblemData(alg, modifiableProblemData))
284            throw new NotSupportedException("The algorithm cannot be used with GBM.");
285
286          IRegressionModel model;
287          IRun run;
288
289          // try to find a model. The algorithm might fail to produce a model. In this case we just retry until the iterations are exhausted
290          if (TryExecute(alg, rand.Next(), RegressionAlgorithmResult, out model, out run)) {
291            int row = 0;
292            // update predictions for training and test
293            // update new targets (in the case of squared error loss we simply use negative residuals)
294            foreach (var pred in model.GetEstimatedValues(problemData.Dataset, trainingRows)) {
295              yPred[row] = yPred[row] + nu * pred;
296              curY[row] = y[row] - yPred[row];
297              row++;
298            }
299            row = 0;
300            foreach (var pred in model.GetEstimatedValues(problemData.Dataset, testRows)) {
301              yPredTest[row] = yPredTest[row] + nu * pred;
302              curYTest[row] = yTest[row] - yPredTest[row];
303              row++;
304            }
305            // determine quality
306            OnlineCalculatorError error;
307            var trainR = OnlinePearsonsRCalculator.Calculate(yPred, y, out error);
308            var testR = OnlinePearsonsRCalculator.Calculate(yPredTest, yTest, out error);
309
310            // iteration results
311            curLoss.Value = error == OnlineCalculatorError.None ? trainR * trainR : 0.0;
312            curTestLoss.Value = error == OnlineCalculatorError.None ? testR * testR : 0.0;
313
314            models.Add(model);
315
316
317          }
318
319          if (StoreRuns)
320            runCollection.Add(run);
321          table.Rows["R² (train)"].Values.Add(curLoss.Value);
322          table.Rows["R² (test)"].Values.Add(curTestLoss.Value);
323          iterations.Value = i + 1;
324        }
325
326        // produce solution
327        if (CreateSolution) {
328          // when all our models are symbolic models we can easily combine them to a single model
329          if (models.All(m => m is ISymbolicRegressionModel)) {
330            Results.Add(new Result("Solution", CreateSymbolicSolution(models, Nu, (IRegressionProblemData)problemData.Clone())));
331          }
332          // just produce an ensemble solution for now (TODO: correct scaling or linear regression for ensemble model weights)
333
334          var ensembleSolution = CreateEnsembleSolution(models, (IRegressionProblemData)problemData.Clone());
335          Results.Add(new Result("EnsembleSolution", ensembleSolution));
336        }
337      }
338      finally {
339        // reset everything
340        alg.Prepare(true);
341      }
342    }
343
344    private static IRegressionEnsembleSolution CreateEnsembleSolution(List<IRegressionModel> models,
345      IRegressionProblemData problemData) {
346      var rows = problemData.TrainingPartition.Size;
347      var features = models.Count;
348      double[,] inputMatrix = new double[rows, features + 1];
349      //add model estimates
350      for (int m = 0; m < models.Count; m++) {
351        var model = models[m];
352        var estimates = model.GetEstimatedValues(problemData.Dataset, problemData.TrainingIndices);
353        int estimatesCounter = 0;
354        foreach (var estimate in estimates) {
355          inputMatrix[estimatesCounter, m] = estimate;
356          estimatesCounter++;
357        }
358      }
359
360      //add target
361      var targets = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TrainingIndices);
362      int targetCounter = 0;
363      foreach (var target in targets) {
364        inputMatrix[targetCounter, models.Count] = target;
365        targetCounter++;
366      }
367
368      alglib.linearmodel lm = new alglib.linearmodel();
369      alglib.lrreport ar = new alglib.lrreport();
370      double[] coefficients;
371      int retVal = 1;
372      alglib.lrbuildz(inputMatrix, rows, features, out retVal, out lm, out ar);
373      if (retVal != 1) throw new ArgumentException("Error in calculation of linear regression solution");
374
375      alglib.lrunpack(lm, out coefficients, out features);
376
377      var ensembleModel = new RegressionEnsembleModel(models, coefficients.Take(models.Count)) { AverageModelEstimates = false };
378      var ensembleSolution = (IRegressionEnsembleSolution)ensembleModel.CreateRegressionSolution(problemData);
379      return ensembleSolution;
380    }
381
382
383    private IAlgorithm CreateOSGP() {
384      // configure strict osgp
385      var alg = new OffspringSelectionGeneticAlgorithm.OffspringSelectionGeneticAlgorithm();
386      var prob = new SymbolicRegressionSingleObjectiveProblem();
387      prob.MaximumSymbolicExpressionTreeDepth.Value = 7;
388      prob.MaximumSymbolicExpressionTreeLength.Value = 15;
389      alg.Problem = prob;
390      alg.SuccessRatio.Value = 1.0;
391      alg.ComparisonFactorLowerBound.Value = 1.0;
392      alg.ComparisonFactorUpperBound.Value = 1.0;
393      alg.MutationProbability.Value = 0.15;
394      alg.PopulationSize.Value = 200;
395      alg.MaximumSelectionPressure.Value = 100;
396      alg.MaximumEvaluatedSolutions.Value = 20000;
397      alg.SelectorParameter.Value = alg.SelectorParameter.ValidValues.OfType<GenderSpecificSelector>().First();
398      alg.MutatorParameter.Value = alg.MutatorParameter.ValidValues.OfType<MultiSymbolicExpressionTreeManipulator>().First();
399      alg.StoreAlgorithmInEachRun = false;
400      return alg;
401    }
402
403    private void SampleTrainingData(MersenneTwister rand, ModifiableDataset ds, int rRows,
404      IDataset sourceDs, double[] curTarget, string targetVarName, IEnumerable<int> trainingIndices) {
405      var selectedRows = trainingIndices.SampleRandomWithoutRepetition(rand, rRows).ToArray();
406      int t = 0;
407      object[] srcRow = new object[ds.Columns];
408      var varNames = ds.DoubleVariables.ToArray();
409      foreach (var r in selectedRows) {
410        // take all values from the original dataset
411        for (int c = 0; c < srcRow.Length; c++) {
412          var col = sourceDs.GetReadOnlyDoubleValues(varNames[c]);
413          srcRow[c] = col[r];
414        }
415        ds.ReplaceRow(t, srcRow);
416        // but use the updated target values
417        ds.SetVariableValue(curTarget[r], targetVarName, t);
418        t++;
419      }
420    }
421
422    private static ISymbolicRegressionSolution CreateSymbolicSolution(List<IRegressionModel> models, double nu, IRegressionProblemData problemData) {
423      var symbModels = models.OfType<ISymbolicRegressionModel>();
424      var lowerLimit = symbModels.Min(m => m.LowerEstimationLimit);
425      var upperLimit = symbModels.Max(m => m.UpperEstimationLimit);
426      var interpreter = new SymbolicDataAnalysisExpressionTreeLinearInterpreter();
427      var progRootNode = new ProgramRootSymbol().CreateTreeNode();
428      var startNode = new StartSymbol().CreateTreeNode();
429
430      var addNode = new Addition().CreateTreeNode();
431      var mulNode = new Multiplication().CreateTreeNode();
432      var scaleNode = (ConstantTreeNode)new Constant().CreateTreeNode(); // all models are scaled using the same nu
433      scaleNode.Value = nu;
434
435      foreach (var m in symbModels) {
436        var relevantPart = m.SymbolicExpressionTree.Root.GetSubtree(0).GetSubtree(0); // skip root and start
437        addNode.AddSubtree((ISymbolicExpressionTreeNode)relevantPart.Clone());
438      }
439
440      mulNode.AddSubtree(addNode);
441      mulNode.AddSubtree(scaleNode);
442      startNode.AddSubtree(mulNode);
443      progRootNode.AddSubtree(startNode);
444      var t = new SymbolicExpressionTree(progRootNode);
445      var combinedModel = new SymbolicRegressionModel(problemData.TargetVariable, t, interpreter, lowerLimit, upperLimit);
446      var sol = new SymbolicRegressionSolution(combinedModel, problemData);
447      return sol;
448    }
449
450    private static bool TrySetProblemData(IAlgorithm alg, IRegressionProblemData problemData) {
451      var prob = alg.Problem as IRegressionProblem;
452      // there is already a problem and it is compatible -> just set problem data
453      if (prob != null) {
454        prob.ProblemDataParameter.Value = problemData;
455        return true;
456      } else return false;
457    }
458
459    private static bool TryExecute(IAlgorithm alg, int seed, string regressionAlgorithmResultName, out IRegressionModel model, out IRun run) {
460      model = null;
461      SetSeed(alg, seed);
462      using (var wh = new AutoResetEvent(false)) {
463        Exception ex = null;
464        EventHandler<EventArgs<Exception>> handler = (sender, args) => {
465          ex = args.Value;
466          wh.Set();
467        };
468        EventHandler handler2 = (sender, args) => wh.Set();
469        alg.ExceptionOccurred += handler;
470        alg.Stopped += handler2;
471        try {
472          alg.Prepare();
473          alg.Start();
474          wh.WaitOne();
475
476          if (ex != null) throw new AggregateException(ex);
477          run = alg.Runs.Last();
478          alg.Runs.Clear();
479          var sols = alg.Results.Select(r => r.Value).OfType<IRegressionSolution>();
480          if (!sols.Any()) return false;
481          var sol = sols.First();
482          if (sols.Skip(1).Any()) {
483            // more than one solution => use regressionAlgorithmResult
484            if (alg.Results.ContainsKey(regressionAlgorithmResultName)) {
485              sol = (IRegressionSolution)alg.Results[regressionAlgorithmResultName].Value;
486            }
487          }
488          var symbRegSol = sol as SymbolicRegressionSolution;
489          // only accept symb reg solutions that do not hit the estimation limits
490          // NaN evaluations would not be critical but are problematic if we want to combine all symbolic models into a single symbolic model
491          if (symbRegSol == null ||
492            (symbRegSol.TrainingLowerEstimationLimitHits == 0 && symbRegSol.TrainingUpperEstimationLimitHits == 0 &&
493             symbRegSol.TestLowerEstimationLimitHits == 0 && symbRegSol.TestUpperEstimationLimitHits == 0) &&
494            symbRegSol.TrainingNaNEvaluations == 0 && symbRegSol.TestNaNEvaluations == 0) {
495            model = sol.Model;
496          }
497        }
498        finally {
499          alg.ExceptionOccurred -= handler;
500          alg.Stopped -= handler2;
501        }
502      }
503      return model != null;
504    }
505
506    private static void SetSeed(IAlgorithm alg, int seed) {
507      // no common interface for algs that use a PRNG -> use naming convention to set seed
508      var paramItem = alg as IParameterizedItem;
509
510      if (paramItem.Parameters.ContainsKey("SetSeedRandomly")) {
511        ((BoolValue)paramItem.Parameters["SetSeedRandomly"].ActualValue).Value = false;
512        ((IntValue)paramItem.Parameters["Seed"].ActualValue).Value = seed;
513      } else {
514        throw new ArgumentException("Base learner does not have a seed parameter (algorithm {0})", alg.Name);
515      }
516
517    }
518  }
519}
Note: See TracBrowser for help on using the repository browser.