Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/MctsSymbolicRegressionStatic.cs @ 13651

Last change on this file since 13651 was 13651, checked in by gkronber, 8 years ago

#2581:

  • added unit tests for the number of different expressions
  • fixed problems in Automaton and constraintHandler that lead to duplicate expressions
  • added possibility for MCTS to handle dead-ends in the search tree (when it is not possible to construct a valid new expression)
  • added statistics on function and gradient evaluations
File size: 20.8 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2015 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
22using System;
23using System.Collections.Generic;
24using System.Diagnostics.Contracts;
25using System.Linq;
26using HeuristicLab.Common;
27using HeuristicLab.Core;
28using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
29using HeuristicLab.Problems.DataAnalysis;
30using HeuristicLab.Problems.DataAnalysis.Symbolic;
31using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;
32using HeuristicLab.Random;
33
34namespace HeuristicLab.Algorithms.DataAnalysis.MctsSymbolicRegression {
35  public static class MctsSymbolicRegressionStatic {
36    // TODO: SGD with adagrad instead of lbfgs?
37    // TODO: check Taylor expansion capabilities (ln(x), sqrt(x), exp(x)) in combination with GBT
38    // TODO: optimize for 3 targets concurrently (y, 1/y, exp(y), and log(y))? Would simplify the number of possible expressions again
39    #region static API
40
41    public interface IState {
42      bool Done { get; }
43      ISymbolicRegressionModel BestModel { get; }
44      double BestSolutionTrainingQuality { get; }
45      double BestSolutionTestQuality { get; }
46      int TotalRollouts { get; }
47      int EffectiveRollouts { get; }
48      int FuncEvaluations { get; }
49      int GradEvaluations { get; } // number of gradient evaluations (* num parameters) to get a value representative of the effort comparable to the number of function evaluations
50      // TODO other stats on LM optimizer might be interesting here
51    }
52
53    // created through factory method
54    private class State : IState {
55      private const int MaxParams = 100;
56
57      // state variables used by MCTS
58      internal readonly Automaton automaton;
59      internal IRandom random { get; private set; }
60      internal readonly double c;
61      internal readonly Tree tree;
62      internal readonly List<Tree> bestChildrenBuf;
63      internal readonly Func<byte[], int, double> evalFun;
64      // MCTS might get stuck. Track statistics on the number of effective rollouts
65      internal int totalRollouts;
66      internal int effectiveRollouts;
67
68
69      // state variables used only internally (for eval function)
70      private readonly IRegressionProblemData problemData;
71      private readonly double[][] x;
72      private readonly double[] y;
73      private readonly double[][] testX;
74      private readonly double[] testY;
75      private readonly double[] scalingFactor;
76      private readonly double[] scalingOffset;
77      private readonly int constOptIterations;
78      private readonly double lowerEstimationLimit, upperEstimationLimit;
79
80      private readonly ExpressionEvaluator evaluator, testEvaluator;
81
82      // values for best solution
83      private double bestRSq;
84      private byte[] bestCode;
85      private int bestNParams;
86      private double[] bestConsts;
87
88      // stats
89      private int funcEvaluations;
90      private int gradEvaluations;
91
92      // buffers
93      private readonly double[] ones; // vector of ones (as default params)
94      private readonly double[] constsBuf;
95      private readonly double[] predBuf, testPredBuf;
96      private readonly double[][] gradBuf;
97
98      public State(IRegressionProblemData problemData, uint randSeed, int maxVariables, double c, bool scaleVariables, int constOptIterations,
99        double lowerEstimationLimit = double.MinValue, double upperEstimationLimit = double.MaxValue,
100        bool allowProdOfVars = true,
101        bool allowExp = true,
102        bool allowLog = true,
103        bool allowInv = true,
104        bool allowMultipleTerms = false) {
105
106        this.problemData = problemData;
107        this.c = c;
108        this.constOptIterations = constOptIterations;
109        this.evalFun = this.Eval;
110        this.lowerEstimationLimit = lowerEstimationLimit;
111        this.upperEstimationLimit = upperEstimationLimit;
112
113        random = new MersenneTwister(randSeed);
114
115        // prepare data for evaluation
116        double[][] x;
117        double[] y;
118        double[][] testX;
119        double[] testY;
120        double[] scalingFactor;
121        double[] scalingOffset;
122        // get training and test datasets (scale linearly based on training set if required)
123        GenerateData(problemData, scaleVariables, problemData.TrainingIndices, out x, out y, out scalingFactor, out scalingOffset);
124        GenerateData(problemData, problemData.TestIndices, scalingFactor, scalingOffset, out testX, out testY);
125        this.x = x;
126        this.y = y;
127        this.testX = testX;
128        this.testY = testY;
129        this.scalingFactor = scalingFactor;
130        this.scalingOffset = scalingOffset;
131        this.evaluator = new ExpressionEvaluator(y.Length, lowerEstimationLimit, upperEstimationLimit);
132        // we need a separate evaluator because the vector length for the test dataset might differ
133        this.testEvaluator = new ExpressionEvaluator(testY.Length, lowerEstimationLimit, upperEstimationLimit);
134
135        this.automaton = new Automaton(x, maxVariables, allowProdOfVars, allowExp, allowLog, allowInv, allowMultipleTerms);
136        this.tree = new Tree() { state = automaton.CurrentState };
137
138        // reset best solution
139        this.bestRSq = 0;
140        // code for default solution (constant model)
141        this.bestCode = new byte[] { (byte)OpCodes.LoadConst0, (byte)OpCodes.Exit };
142        this.bestNParams = 0;
143        this.bestConsts = null;
144
145        // init buffers
146        this.ones = Enumerable.Repeat(1.0, MaxParams).ToArray();
147        constsBuf = new double[MaxParams];
148        this.bestChildrenBuf = new List<Tree>(2 * x.Length); // the number of follow states in the automaton is O(number of variables) 2 * number of variables should be sufficient (capacity is increased if necessary anyway)
149        this.predBuf = new double[y.Length];
150        this.testPredBuf = new double[testY.Length];
151
152        this.gradBuf = Enumerable.Range(0, MaxParams).Select(_ => new double[y.Length]).ToArray();
153      }
154
155      #region IState inferface
156      public bool Done { get { return tree != null && tree.done; } }
157
158      public double BestSolutionTrainingQuality {
159        get {
160          evaluator.Exec(bestCode, x, bestConsts, predBuf);
161          return RSq(y, predBuf);
162        }
163      }
164
165      public double BestSolutionTestQuality {
166        get {
167          testEvaluator.Exec(bestCode, testX, bestConsts, testPredBuf);
168          return RSq(testY, testPredBuf);
169        }
170      }
171
172      // takes the code of the best solution and creates and equivalent symbolic regression model
173      public ISymbolicRegressionModel BestModel {
174        get {
175          var treeGen = new SymbolicExpressionTreeGenerator(problemData.AllowedInputVariables.ToArray());
176          var interpreter = new SymbolicDataAnalysisExpressionTreeLinearInterpreter();
177          var simplifier = new SymbolicDataAnalysisExpressionTreeSimplifier();
178
179          var t = new SymbolicExpressionTree(treeGen.Exec(bestCode, bestConsts, bestNParams, scalingFactor, scalingOffset));
180          var simpleT = simplifier.Simplify(t);
181          var model = new SymbolicRegressionModel(simpleT, interpreter, lowerEstimationLimit, upperEstimationLimit);
182
183          // model has already been scaled linearly in Eval
184          return model;
185        }
186      }
187
188      public int TotalRollouts { get { return totalRollouts; } }
189      public int EffectiveRollouts { get { return effectiveRollouts; } }
190      public int FuncEvaluations { get { return funcEvaluations; } }
191      public int GradEvaluations { get { return gradEvaluations; } } // number of gradient evaluations (* num parameters) to get a value representative of the effort comparable to the number of function evaluations
192
193      #endregion
194
195      private double Eval(byte[] code, int nParams) {
196        double[] optConsts;
197        double q;
198        Eval(code, nParams, out q, out optConsts);
199
200        if (q > bestRSq) {
201          bestRSq = q;
202          bestNParams = nParams;
203          this.bestCode = new byte[code.Length];
204          this.bestConsts = new double[bestNParams];
205
206          Array.Copy(code, bestCode, code.Length);
207          Array.Copy(optConsts, bestConsts, bestNParams);
208        }
209
210        return q;
211      }
212
213      private void Eval(byte[] code, int nParams, out double rsq, out double[] optConsts) {
214        // we make a first pass to determine a valid starting configuration for all constants
215        // constant c in log(c + f(x)) is adjusted to guarantee that x is positive (see expression evaluator)
216        // scale and offset are set to optimal starting configuration
217        // assumes scale is the first param and offset is the last param
218        double alpha;
219        double beta;
220
221        // reset constants
222        Array.Copy(ones, constsBuf, nParams);
223        evaluator.Exec(code, x, constsBuf, predBuf, adjustOffsetForLogAndExp: true);
224        funcEvaluations++;
225
226        // calc opt scaling (alpha*f(x) + beta)
227        OnlineCalculatorError error;
228        OnlineLinearScalingParameterCalculator.Calculate(predBuf, y, out alpha, out beta, out error);
229        if (error == OnlineCalculatorError.None) {
230          constsBuf[0] *= beta;
231          constsBuf[nParams - 1] = constsBuf[nParams - 1] * beta + alpha;
232        }
233        if (nParams <= 2 || constOptIterations <= 0) {
234          // if we don't need to optimize parameters then we are done
235          // changing scale and offset does not influence r²
236          rsq = RSq(y, predBuf);
237          optConsts = constsBuf;
238        } else {
239          // optimize constants using the starting point calculated above
240          OptimizeConstsLm(code, constsBuf, nParams, 0.0, nIters: constOptIterations);
241
242          evaluator.Exec(code, x, constsBuf, predBuf);
243          funcEvaluations++;
244
245          rsq = RSq(y, predBuf);
246          optConsts = constsBuf;
247        }
248      }
249
250
251
252      #region helpers
253      private static double RSq(IEnumerable<double> x, IEnumerable<double> y) {
254        OnlineCalculatorError error;
255        double r = OnlinePearsonsRCalculator.Calculate(x, y, out error);
256        return error == OnlineCalculatorError.None ? r * r : 0.0;
257      }
258
259
260      private void OptimizeConstsLm(byte[] code, double[] consts, int nParams, double epsF = 0.0, int nIters = 100) {
261        double[] optConsts = new double[nParams]; // allocate a smaller buffer for constants opt (TODO perf?)
262        Array.Copy(consts, optConsts, nParams);
263
264        alglib.minlmstate state;
265        alglib.minlmreport rep = null;
266        alglib.minlmcreatevj(y.Length, optConsts, out state);
267        alglib.minlmsetcond(state, 0.0, epsF, 0.0, nIters);
268        //alglib.minlmsetgradientcheck(state, 0.000001);
269        alglib.minlmoptimize(state, Func, FuncAndJacobian, null, code);
270        alglib.minlmresults(state, out optConsts, out rep);
271        funcEvaluations += rep.nfunc;
272        gradEvaluations += rep.njac * nParams;
273
274        if (rep.terminationtype < 0) throw new ArgumentException("lm failed: termination type = " + rep.terminationtype);
275
276        // only use optimized constants if successful
277        if (rep.terminationtype >= 0) {
278          Array.Copy(optConsts, consts, optConsts.Length);
279        }
280      }
281
282      private void Func(double[] arg, double[] fi, object obj) {
283        // 0.5 * MSE and gradient
284        var code = (byte[])obj;
285        evaluator.Exec(code, x, arg, predBuf); // gradients are nParams x vLen
286        for (int r = 0; r < predBuf.Length; r++) {
287          var res = predBuf[r] - y[r];
288          fi[r] = res;
289        }
290      }
291      private void FuncAndJacobian(double[] arg, double[] fi, double[,] jac, object obj) {
292        int nParams = arg.Length;
293        var code = (byte[])obj;
294        evaluator.ExecGradient(code, x, arg, predBuf, gradBuf); // gradients are nParams x vLen
295        for (int r = 0; r < predBuf.Length; r++) {
296          var res = predBuf[r] - y[r];
297          fi[r] = res;
298
299          for (int k = 0; k < nParams; k++) {
300            jac[r, k] = gradBuf[k][r];
301          }
302        }
303      }
304      #endregion
305    }
306
307    public static IState CreateState(IRegressionProblemData problemData, uint randSeed, int maxVariables = 3, double c = 1.0,
308      bool scaleVariables = true, int constOptIterations = 0, double lowerEstimationLimit = double.MinValue, double upperEstimationLimit = double.MaxValue,
309      bool allowProdOfVars = true,
310      bool allowExp = true,
311      bool allowLog = true,
312      bool allowInv = true,
313      bool allowMultipleTerms = false
314      ) {
315      return new State(problemData, randSeed, maxVariables, c, scaleVariables, constOptIterations,
316        lowerEstimationLimit, upperEstimationLimit,
317        allowProdOfVars, allowExp, allowLog, allowInv, allowMultipleTerms);
318    }
319
320    // returns the quality of the evaluated solution
321    public static double MakeStep(IState state) {
322      var mctsState = state as State;
323      if (mctsState == null) throw new ArgumentException("state");
324      if (mctsState.Done) throw new NotSupportedException("The tree search has enumerated all possible solutions.");
325
326      return TreeSearch(mctsState);
327    }
328    #endregion
329
330    private static double TreeSearch(State mctsState) {
331      var automaton = mctsState.automaton;
332      var tree = mctsState.tree;
333      var eval = mctsState.evalFun;
334      var bestChildrenBuf = mctsState.bestChildrenBuf;
335      var rand = mctsState.random;
336      double c = mctsState.c;
337      double q = 0;
338      bool success = false;
339      do {
340        automaton.Reset();
341        success = TryTreeSearchRec(rand, tree, c, automaton, eval, bestChildrenBuf, out q);
342        mctsState.totalRollouts++;
343      } while (!success && !tree.done);
344      mctsState.effectiveRollouts++;
345      return q;
346    }
347
348    // tree search might fail because of constraints for expressions
349    // in this case we get stuck we just restart
350    // see ConstraintHandler.cs for more info
351    private static bool TryTreeSearchRec(IRandom rand, Tree tree, double c, Automaton automaton, Func<byte[], int, double> eval, List<Tree> bestChildrenBuf,
352      out double q) {
353      Tree selectedChild = null;
354      Contract.Assert(tree.state == automaton.CurrentState);
355      Contract.Assert(!tree.done);
356      if (tree.children == null) {
357        if (automaton.IsFinalState(tree.state)) {
358          // final state
359          tree.done = true;
360
361          // EVALUATE
362          byte[] code; int nParams;
363          automaton.GetCode(out code, out nParams);
364          q = eval(code, nParams);
365          tree.visits++;
366          tree.sumQuality += q;
367          return true; // we reached a final state
368        } else {
369          // EXPAND
370          int[] possibleFollowStates;
371          int nFs;
372          automaton.FollowStates(automaton.CurrentState, out possibleFollowStates, out nFs);
373          if (nFs == 0) {
374            // stuck in a dead end (no final state and no allowed follow states)
375            q = 0;
376            tree.done = true;
377            tree.children = null;
378            return false;
379          }
380          tree.children = new Tree[nFs];
381          for (int i = 0; i < tree.children.Length; i++)
382            tree.children[i] = new Tree() { children = null, done = false, state = possibleFollowStates[i], visits = 0 };
383
384          selectedChild = SelectFinalOrRandom(automaton, tree, rand);
385        }
386      } else {
387        // tree.children != null
388        // UCT selection within tree
389        selectedChild = SelectUct(tree, rand, c, bestChildrenBuf);
390      }
391      // make selected step and recurse
392      automaton.Goto(selectedChild.state);
393      var success = TryTreeSearchRec(rand, selectedChild, c, automaton, eval, bestChildrenBuf, out q);
394      if (success) {
395        // only update if successful
396        tree.sumQuality += q;
397        tree.visits++;
398      }
399
400      // tree.done = tree.children.All(ch => ch.done);
401      tree.done = true; for (int i = 0; i < tree.children.Length && tree.done; i++) tree.done = tree.children[i].done;
402      if (tree.done) {
403        tree.children = null; // cut off the sub-branch if it has been fully explored
404        // TODO: update all qualities and visits to remove the information gained from this whole branch
405      }
406      return success;
407    }
408
409    private static Tree SelectUct(Tree tree, IRandom rand, double c, List<Tree> bestChildrenBuf) {
410      // determine total tries of still active children
411      int totalTries = 0;
412      bestChildrenBuf.Clear();
413      for (int i = 0; i < tree.children.Length; i++) {
414        var ch = tree.children[i];
415        if (ch.done) continue;
416        if (ch.visits == 0) bestChildrenBuf.Add(ch);
417        else totalTries += tree.children[i].visits;
418      }
419      // if there are unvisited children select a random child
420      if (bestChildrenBuf.Any()) {
421        return bestChildrenBuf[rand.Next(bestChildrenBuf.Count)];
422      }
423      Contract.Assert(totalTries > 0); // the tree is not done yet so there is at least on child that is not done
424      double logTotalTries = Math.Log(totalTries);
425      var bestQ = double.NegativeInfinity;
426      for (int i = 0; i < tree.children.Length; i++) {
427        var ch = tree.children[i];
428        if (ch.done) continue;
429        var childQ = ch.AverageQuality + c * Math.Sqrt(logTotalTries / ch.visits);
430        if (childQ > bestQ) {
431          bestChildrenBuf.Clear();
432          bestChildrenBuf.Add(ch);
433          bestQ = childQ;
434        } else if (childQ >= bestQ) {
435          bestChildrenBuf.Add(ch);
436        }
437      }
438      return bestChildrenBuf.Count > 0 ? bestChildrenBuf[rand.Next(bestChildrenBuf.Count)] : bestChildrenBuf[0];
439    }
440
441    private static Tree SelectFinalOrRandom(Automaton automaton, Tree tree, IRandom rand) {
442      // if one of the new children leads to a final state then go there
443      // otherwise choose a random child
444      int selectedChildIdx = -1;
445      // find first final state if there is one
446      for (int i = 0; i < tree.children.Length; i++) {
447        if (automaton.IsFinalState(tree.children[i].state)) {
448          selectedChildIdx = i;
449          break;
450        }
451      }
452      // no final state -> select a random child
453      if (selectedChildIdx == -1) {
454        selectedChildIdx = rand.Next(tree.children.Length);
455      }
456      return tree.children[selectedChildIdx];
457    }
458
459    // scales data and extracts values from dataset into arrays
460    private static void GenerateData(IRegressionProblemData problemData, bool scaleVariables, IEnumerable<int> rows,
461      out double[][] xs, out double[] y, out double[] scalingFactor, out double[] scalingOffset) {
462      xs = new double[problemData.AllowedInputVariables.Count()][];
463
464      var i = 0;
465      if (scaleVariables) {
466        scalingFactor = new double[xs.Length];
467        scalingOffset = new double[xs.Length];
468      } else {
469        scalingFactor = null;
470        scalingOffset = null;
471      }
472      foreach (var var in problemData.AllowedInputVariables) {
473        if (scaleVariables) {
474          var minX = problemData.Dataset.GetDoubleValues(var, rows).Min();
475          var maxX = problemData.Dataset.GetDoubleValues(var, rows).Max();
476          var range = maxX - minX;
477
478          // scaledX = (x - min) / range
479          var sf = 1.0 / range;
480          var offset = -minX / range;
481          scalingFactor[i] = sf;
482          scalingOffset[i] = offset;
483          i++;
484        }
485      }
486
487      GenerateData(problemData, rows, scalingFactor, scalingOffset, out xs, out y);
488    }
489
490    // extract values from dataset into arrays
491    private static void GenerateData(IRegressionProblemData problemData, IEnumerable<int> rows, double[] scalingFactor, double[] scalingOffset,
492     out double[][] xs, out double[] y) {
493      xs = new double[problemData.AllowedInputVariables.Count()][];
494
495      int i = 0;
496      foreach (var var in problemData.AllowedInputVariables) {
497        var sf = scalingFactor == null ? 1.0 : scalingFactor[i];
498        var offset = scalingFactor == null ? 0.0 : scalingOffset[i];
499        xs[i++] =
500          problemData.Dataset.GetDoubleValues(var, rows).Select(xi => xi * sf + offset).ToArray();
501      }
502
503      y = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, rows).ToArray();
504    }
505  }
506}
Note: See TracBrowser for help on using the repository browser.