1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022015 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.Diagnostics.Contracts;


25  using System.Linq;


26  using HeuristicLab.Common;


27  using HeuristicLab.Core;


28  using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;


29  using HeuristicLab.Problems.DataAnalysis;


30  using HeuristicLab.Problems.DataAnalysis.Symbolic;


31  using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;


32  using HeuristicLab.Random;


33 


34  namespace 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 


178  var t = new SymbolicExpressionTree(treeGen.Exec(bestCode, bestConsts, bestNParams, scalingFactor, scalingOffset));


179  var model = new SymbolicRegressionModel(t, interpreter, lowerEstimationLimit, upperEstimationLimit);


180 


181  // model has already been scaled linearly in Eval


182  return model;


183  }


184  }


185 


186  public int TotalRollouts { get { return totalRollouts; } }


187  public int EffectiveRollouts { get { return effectiveRollouts; } }


188  public int FuncEvaluations { get { return funcEvaluations; } }


189  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


190 


191  #endregion


192 


193  private double Eval(byte[] code, int nParams) {


194  double[] optConsts;


195  double q;


196  Eval(code, nParams, out q, out optConsts);


197 


198  if (q > bestRSq) {


199  bestRSq = q;


200  bestNParams = nParams;


201  this.bestCode = new byte[code.Length];


202  this.bestConsts = new double[bestNParams];


203 


204  Array.Copy(code, bestCode, code.Length);


205  Array.Copy(optConsts, bestConsts, bestNParams);


206  }


207 


208  return q;


209  }


210 


211  private void Eval(byte[] code, int nParams, out double rsq, out double[] optConsts) {


212  // we make a first pass to determine a valid starting configuration for all constants


213  // constant c in log(c + f(x)) is adjusted to guarantee that x is positive (see expression evaluator)


214  // scale and offset are set to optimal starting configuration


215  // assumes scale is the first param and offset is the last param


216  double alpha;


217  double beta;


218 


219  // reset constants


220  Array.Copy(ones, constsBuf, nParams);


221  evaluator.Exec(code, x, constsBuf, predBuf, adjustOffsetForLogAndExp: true);


222  funcEvaluations++;


223 


224  // calc opt scaling (alpha*f(x) + beta)


225  OnlineCalculatorError error;


226  OnlineLinearScalingParameterCalculator.Calculate(predBuf, y, out alpha, out beta, out error);


227  if (error == OnlineCalculatorError.None) {


228  constsBuf[0] *= beta;


229  constsBuf[nParams  1] = constsBuf[nParams  1] * beta + alpha;


230  }


231  if (nParams <= 2  constOptIterations <= 0) {


232  // if we don't need to optimize parameters then we are done


233  // changing scale and offset does not influence r²


234  rsq = RSq(y, predBuf);


235  optConsts = constsBuf;


236  } else {


237  // optimize constants using the starting point calculated above


238  OptimizeConstsLm(code, constsBuf, nParams, 0.0, nIters: constOptIterations);


239 


240  evaluator.Exec(code, x, constsBuf, predBuf);


241  funcEvaluations++;


242 


243  rsq = RSq(y, predBuf);


244  optConsts = constsBuf;


245  }


246  }


247 


248 


249 


250  #region helpers


251  private static double RSq(IEnumerable<double> x, IEnumerable<double> y) {


252  OnlineCalculatorError error;


253  double r = OnlinePearsonsRCalculator.Calculate(x, y, out error);


254  return error == OnlineCalculatorError.None ? r * r : 0.0;


255  }


256 


257 


258  private void OptimizeConstsLm(byte[] code, double[] consts, int nParams, double epsF = 0.0, int nIters = 100) {


259  double[] optConsts = new double[nParams]; // allocate a smaller buffer for constants opt (TODO perf?)


260  Array.Copy(consts, optConsts, nParams);


261 


262  alglib.minlmstate state;


263  alglib.minlmreport rep = null;


264  alglib.minlmcreatevj(y.Length, optConsts, out state);


265  alglib.minlmsetcond(state, 0.0, epsF, 0.0, nIters);


266  //alglib.minlmsetgradientcheck(state, 0.000001);


267  alglib.minlmoptimize(state, Func, FuncAndJacobian, null, code);


268  alglib.minlmresults(state, out optConsts, out rep);


269  funcEvaluations += rep.nfunc;


270  gradEvaluations += rep.njac * nParams;


271 


272  if (rep.terminationtype < 0) throw new ArgumentException("lm failed: termination type = " + rep.terminationtype);


273 


274  // only use optimized constants if successful


275  if (rep.terminationtype >= 0) {


276  Array.Copy(optConsts, consts, optConsts.Length);


277  }


278  }


279 


280  private void Func(double[] arg, double[] fi, object obj) {


281  var code = (byte[])obj;


282  evaluator.Exec(code, x, arg, predBuf); // gradients are nParams x vLen


283  for (int r = 0; r < predBuf.Length; r++) {


284  var res = predBuf[r]  y[r];


285  fi[r] = res;


286  }


287  }


288  private void FuncAndJacobian(double[] arg, double[] fi, double[,] jac, object obj) {


289  int nParams = arg.Length;


290  var code = (byte[])obj;


291  evaluator.ExecGradient(code, x, arg, predBuf, gradBuf); // gradients are nParams x vLen


292  for (int r = 0; r < predBuf.Length; r++) {


293  var res = predBuf[r]  y[r];


294  fi[r] = res;


295 


296  for (int k = 0; k < nParams; k++) {


297  jac[r, k] = gradBuf[k][r];


298  }


299  }


300  }


301  #endregion


302  }


303 


304  public static IState CreateState(IRegressionProblemData problemData, uint randSeed, int maxVariables = 3, double c = 1.0,


305  bool scaleVariables = true, int constOptIterations = 0, double lowerEstimationLimit = double.MinValue, double upperEstimationLimit = double.MaxValue,


306  bool allowProdOfVars = true,


307  bool allowExp = true,


308  bool allowLog = true,


309  bool allowInv = true,


310  bool allowMultipleTerms = false


311  ) {


312  return new State(problemData, randSeed, maxVariables, c, scaleVariables, constOptIterations,


313  lowerEstimationLimit, upperEstimationLimit,


314  allowProdOfVars, allowExp, allowLog, allowInv, allowMultipleTerms);


315  }


316 


317  // returns the quality of the evaluated solution


318  public static double MakeStep(IState state) {


319  var mctsState = state as State;


320  if (mctsState == null) throw new ArgumentException("state");


321  if (mctsState.Done) throw new NotSupportedException("The tree search has enumerated all possible solutions.");


322 


323  return TreeSearch(mctsState);


324  }


325  #endregion


326 


327  private static double TreeSearch(State mctsState) {


328  var automaton = mctsState.automaton;


329  var tree = mctsState.tree;


330  var eval = mctsState.evalFun;


331  var bestChildrenBuf = mctsState.bestChildrenBuf;


332  var rand = mctsState.random;


333  double c = mctsState.c;


334  double q = 0;


335  bool success = false;


336  do {


337  automaton.Reset();


338  success = TryTreeSearchRec(rand, tree, c, automaton, eval, bestChildrenBuf, out q);


339  mctsState.totalRollouts++;


340  } while (!success && !tree.done);


341  mctsState.effectiveRollouts++;


342  return q;


343  }


344 


345  // tree search might fail because of constraints for expressions


346  // in this case we get stuck we just restart


347  // see ConstraintHandler.cs for more info


348  private static bool TryTreeSearchRec(IRandom rand, Tree tree, double c, Automaton automaton, Func<byte[], int, double> eval, List<Tree> bestChildrenBuf,


349  out double q) {


350  Tree selectedChild = null;


351  Contract.Assert(tree.state == automaton.CurrentState);


352  Contract.Assert(!tree.done);


353  if (tree.children == null) {


354  if (automaton.IsFinalState(tree.state)) {


355  // final state


356  tree.done = true;


357 


358  // EVALUATE


359  byte[] code; int nParams;


360  automaton.GetCode(out code, out nParams);


361  q = eval(code, nParams);


362  tree.visits++;


363  tree.sumQuality += q;


364  return true; // we reached a final state


365  } else {


366  // EXPAND


367  int[] possibleFollowStates;


368  int nFs;


369  automaton.FollowStates(automaton.CurrentState, out possibleFollowStates, out nFs);


370  if (nFs == 0) {


371  // stuck in a dead end (no final state and no allowed follow states)


372  q = 0;


373  tree.done = true;


374  tree.children = null;


375  tree.visits = 1;


376  return false;


377  }


378  tree.children = new Tree[nFs];


379  for (int i = 0; i < tree.children.Length; i++)


380  tree.children[i] = new Tree() { children = null, done = false, state = possibleFollowStates[i], visits = 0 };


381 


382  selectedChild = SelectFinalOrRandom(automaton, tree, rand);


383  }


384  } else {


385  // tree.children != null


386  // UCT selection within tree


387  selectedChild = SelectUct(tree, rand, c, bestChildrenBuf);


388  }


389  // make selected step and recurse


390  automaton.Goto(selectedChild.state);


391  var success = TryTreeSearchRec(rand, selectedChild, c, automaton, eval, bestChildrenBuf, out q);


392  if (success) {


393  // only update if successful


394  tree.sumQuality += q;


395  tree.visits++;


396  }


397 


398  // tree.done = tree.children.All(ch => ch.done);


399  tree.done = true; for (int i = 0; i < tree.children.Length && tree.done; i++) tree.done = tree.children[i].done;


400  if (tree.done) {


401  tree.children = null; // cut off the subbranch if it has been fully explored


402  // TODO: update all qualities and visits to remove the information gained from this whole branch


403  }


404  return success;


405  }


406 


407  private static Tree SelectUct(Tree tree, IRandom rand, double c, List<Tree> bestChildrenBuf) {


408  // determine total tries of still active children


409  int totalTries = 0;


410  bestChildrenBuf.Clear();


411  for (int i = 0; i < tree.children.Length; i++) {


412  var ch = tree.children[i];


413  if (ch.done) continue;


414  if (ch.visits == 0) bestChildrenBuf.Add(ch);


415  else totalTries += tree.children[i].visits;


416  }


417  // if there are unvisited children select a random child


418  if (bestChildrenBuf.Any()) {


419  return bestChildrenBuf[rand.Next(bestChildrenBuf.Count)];


420  }


421  Contract.Assert(totalTries > 0); // the tree is not done yet so there is at least on child that is not done


422  double logTotalTries = Math.Log(totalTries);


423  var bestQ = double.NegativeInfinity;


424  for (int i = 0; i < tree.children.Length; i++) {


425  var ch = tree.children[i];


426  if (ch.done) continue;


427  var childQ = ch.AverageQuality + c * Math.Sqrt(logTotalTries / ch.visits);


428  if (childQ > bestQ) {


429  bestChildrenBuf.Clear();


430  bestChildrenBuf.Add(ch);


431  bestQ = childQ;


432  } else if (childQ >= bestQ) {


433  bestChildrenBuf.Add(ch);


434  }


435  }


436  return bestChildrenBuf[rand.Next(bestChildrenBuf.Count)];


437  }


438 


439  private static Tree SelectFinalOrRandom(Automaton automaton, Tree tree, IRandom rand) {


440  // if one of the new children leads to a final state then go there


441  // otherwise choose a random child


442  int selectedChildIdx = 1;


443  // find first final state if there is one


444  for (int i = 0; i < tree.children.Length; i++) {


445  if (automaton.IsFinalState(tree.children[i].state)) {


446  selectedChildIdx = i;


447  break;


448  }


449  }


450  // no final state > select a random child


451  if (selectedChildIdx == 1) {


452  selectedChildIdx = rand.Next(tree.children.Length);


453  }


454  return tree.children[selectedChildIdx];


455  }


456 


457  // scales data and extracts values from dataset into arrays


458  private static void GenerateData(IRegressionProblemData problemData, bool scaleVariables, IEnumerable<int> rows,


459  out double[][] xs, out double[] y, out double[] scalingFactor, out double[] scalingOffset) {


460  xs = new double[problemData.AllowedInputVariables.Count()][];


461 


462  var i = 0;


463  if (scaleVariables) {


464  scalingFactor = new double[xs.Length];


465  scalingOffset = new double[xs.Length];


466  } else {


467  scalingFactor = null;


468  scalingOffset = null;


469  }


470  foreach (var var in problemData.AllowedInputVariables) {


471  if (scaleVariables) {


472  var minX = problemData.Dataset.GetDoubleValues(var, rows).Min();


473  var maxX = problemData.Dataset.GetDoubleValues(var, rows).Max();


474  var range = maxX  minX;


475 


476  // scaledX = (x  min) / range


477  var sf = 1.0 / range;


478  var offset = minX / range;


479  scalingFactor[i] = sf;


480  scalingOffset[i] = offset;


481  i++;


482  }


483  }


484 


485  GenerateData(problemData, rows, scalingFactor, scalingOffset, out xs, out y);


486  }


487 


488  // extract values from dataset into arrays


489  private static void GenerateData(IRegressionProblemData problemData, IEnumerable<int> rows, double[] scalingFactor, double[] scalingOffset,


490  out double[][] xs, out double[] y) {


491  xs = new double[problemData.AllowedInputVariables.Count()][];


492 


493  int i = 0;


494  foreach (var var in problemData.AllowedInputVariables) {


495  var sf = scalingFactor == null ? 1.0 : scalingFactor[i];


496  var offset = scalingFactor == null ? 0.0 : scalingOffset[i];


497  xs[i++] =


498  problemData.Dataset.GetDoubleValues(var, rows).Select(xi => xi * sf + offset).ToArray();


499  }


500 


501  y = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, rows).ToArray();


502  }


503  }


504  }

