Changeset 15403 for branches/MCTS-SymbReg-2796/HeuristicLab.Algorithms.DataAnalysis

Timestamp:

10/04/17 22:00:52 (7 years ago)

Author:

gkronber

Message:

#2796 worked on MCTS symb reg

Location:

branches/MCTS-SymbReg-2796

Files:

• . (modified) (1 prop)
• HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/Automaton.cs (modified) (3 diffs)
• HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/ConstraintHandler.cs (modified) (1 diff)
• HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/ExpressionEvaluator.cs (modified) (1 diff)
• HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/MctsSymbolicRegressionAlgorithm.cs (modified) (13 diffs)
• HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/MctsSymbolicRegressionStatic.cs (modified) (14 diffs)

branches/MCTS-SymbReg-2796/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/Automaton.cs

@@ -96 +96 @@
     }
 
-    // reverse notation ops
-    // Expr -> c 0 Term { '+' Term } '+' '*' c '+' 'exit'
+    // reverse notation 
+    // Expr -> 0 Term { '+' Term } '+' 'exit'
     // Term -> c Fact { '*' Fact } '*'
     // Fact -> VarFact | ExpFact | LogFact | InvFact
@@ -120 +120 @@
       actionStrings = new List<string>[nStates, nStates];
 
-      // Expr -> c 0 Term { '+' Term } '+' '*' c '+' 'exit'
+      // Expr -> 0 Term { '+' Term } '+' 'exit'
       AddTransition(StateExpr, StateTermStart, () => {
         codeGenerator.Reset();
-        codeGenerator.Emit1(OpCodes.LoadParamN);
         codeGenerator.Emit1(OpCodes.LoadConst0);
         constraintHandler.Reset();
-      }, "c 0, Reset");
+      }, "0, Reset");
       AddTransition(StateTermEnd, StateExprEnd, () => {
         codeGenerator.Emit1(OpCodes.Add);
-        codeGenerator.Emit1(OpCodes.Mul);
-        codeGenerator.Emit1(OpCodes.LoadParamN);
-        codeGenerator.Emit1(OpCodes.Add);
         codeGenerator.Emit1(OpCodes.Exit);
-      }, "+*c+ exit");
+      }, "+ exit");
       if (allowMultipleTerms)
         AddTransition(StateTermEnd, StateTermStart, () => {
@@ -353 +349 @@
     private readonly int[] followStatesBuf = new int[1000];
     public void FollowStates(int state, out int[] buf, out int nElements) {
-      // for loop instead of where iterator
       var fs = followStates[state];
       int j = 0;

branches/MCTS-SymbReg-2796/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/ConstraintHandler.cs

@@ -29 +29 @@
   // This class restricts the set of allowed transitions of the automaton to prevent exploration of duplicate expressions.
   // It would be possible to implement this class in such a way that the search never visits a duplicate expression. However,
-  // it seems very intricate to detect this robustly and in all cases while generating an expression because 
-  // some for of lookahead is necessary. 
-  // Instead the constraint handler only catches the obvious duplicates directly, but does not guarantee that the search always produces a valid expression.
+  // it seems very intricate to detect this robustly while generating an expression because 
+  // some form of lookahead is necessary. 
+  // Instead, the constraint handler only catches the obvious duplicates directly, but does not guarantee that the search always produces a valid expression.
   // The ratio of the number of unsuccessful searches, that need backtracking should be tracked in the MCTS alg (MctsSymbolicRegressionStatic)
 

branches/MCTS-SymbReg-2796/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/ExpressionEvaluator.cs

@@ -161 +161 @@
                 var minFx = fxc.Min() - consts[curParamIdx]; // stack[topOfStack] is f(x) + c
 
-                var delta = 1.0 - minFx - consts[curParamIdx];
-                // adjust c so that minFx + c = 1 ... log(minFx + c) = 0
+                // adjust c so that minFx + c = e and therefore log(minFx + c) = log(e) = 1
+                // this initialization works in combination with the gradient check (instead of initializing such that log(minFx + c) = 0
+                var delta = Math.E - minFx - consts[curParamIdx];
                 consts[curParamIdx] += delta;
 

branches/MCTS-SymbReg-2796/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/MctsSymbolicRegressionAlgorithm.cs

@@ -22 +22 @@
 using System;
 using System.Linq;
-using System.Runtime.CompilerServices;
 using System.Threading;
 using HeuristicLab.Algorithms.DataAnalysis.MctsSymbolicRegression.Policies;
@@ -36 +35 @@
 
 namespace HeuristicLab.Algorithms.DataAnalysis.MctsSymbolicRegression {
-  [Item("MCTS Symbolic Regression", "Monte carlo tree search for symbolic regression. Useful mainly as a base learner in gradient boosting.")]
+  // TODO: support pause (persisting/cloning the state)
+  [Item("MCTS Symbolic Regression", "Monte carlo tree search for symbolic regression.")]
   [StorableClass]
   [Creatable(CreatableAttribute.Categories.DataAnalysisRegression, Priority = 250)]
@@ -53 +53 @@
     private const string CreateSolutionParameterName = "CreateSolution";
     private const string PunishmentFactorParameterName = "PunishmentFactor";
-
-    private const string VariableProductFactorName = "product(xi)";
-    private const string ExpFactorName = "exp(c * product(xi))";
-    private const string LogFactorName = "log(c + sum(c*product(xi))";
-    private const string InvFactorName = "1 / (1 + sum(c*product(xi))";
-    private const string FactorSumsName = "sum of multiple terms";
+    private const string CollectParetoOptimalSolutionsParameterName = "CollectParetoOptimalSolutions";
+    private const string LambdaParameterName = "Lambda";
+
+    private const string VariableProductFactorName = "x * y * ...";
+    private const string ExpFactorName = "exp(c * x * y ...)";
+    private const string LogFactorName = "log(c + c1 x + c2 x + ...)";
+    private const string InvFactorName = "1 / (1 + c1 x + c2 x + ...)";
+    private const string FactorSumsName = "t1(x) + t2(x) + ... ";
     #endregion
 
@@ -94 +96 @@
     public IFixedValueParameter<BoolValue> CreateSolutionParameter {
       get { return (IFixedValueParameter<BoolValue>)Parameters[CreateSolutionParameterName]; }
+    }
+    public IFixedValueParameter<BoolValue> CollectParetoOptimalSolutionsParameter {
+      get { return (IFixedValueParameter<BoolValue>)Parameters[CollectParetoOptimalSolutionsParameterName]; }
+    }
+    public IFixedValueParameter<DoubleValue> LambdaParameter {
+      get { return (IFixedValueParameter<DoubleValue>)Parameters[LambdaParameterName]; }
     }
     #endregion
@@ -136 +144 @@
       get { return CreateSolutionParameter.Value.Value; }
       set { CreateSolutionParameter.Value.Value = value; }
+    }
+    public bool CollectParetoOptimalSolutions {
+      get { return CollectParetoOptimalSolutionsParameter.Value.Value; }
+      set { CollectParetoOptimalSolutionsParameter.Value.Value = value; }
+    }
+    public double Lambda {
+      get { return LambdaParameter.Value.Value; }
+      set { LambdaParameter.Value.Value = value; }
     }
     #endregion
@@ -177 +193 @@
       Parameters.Add(new FixedValueParameter<IntValue>(ConstantOptimizationIterationsParameterName,
         "Number of iterations for constant optimization. A small number of iterations should be sufficient for most models. " +
-        "Set to 0 to disable constants optimization.", new IntValue(10)));
+        "Set to 0 to let the algorithm stop automatically when it converges. Set to -1 to disable constants optimization.", new IntValue(10)));
       Parameters.Add(new FixedValueParameter<BoolValue>(ScaleVariablesParameterName,
-        "Set to true to scale all input variables to the range [0..1]", new BoolValue(false)));
+        "Set to true to all input variables to the range [0..1]", new BoolValue(true)));
       Parameters[ScaleVariablesParameterName].Hidden = true;
       Parameters.Add(new FixedValueParameter<DoubleValue>(PunishmentFactorParameterName, "Estimations of models can be bounded. The estimation limits are calculated in the following way (lb = mean(y) - punishmentFactor*range(y), ub = mean(y) + punishmentFactor*range(y))", new DoubleValue(10)));
@@ -187 +203 @@
       Parameters[UpdateIntervalParameterName].Hidden = true;
       Parameters.Add(new FixedValueParameter<BoolValue>(CreateSolutionParameterName,
-        "Flag that indicates if a solution should be produced at the end of the run", new BoolValue(true)));
+        "Optionally produce a solution at the end of the run", new BoolValue(true)));
       Parameters[CreateSolutionParameterName].Hidden = true;
+
+      Parameters.Add(new FixedValueParameter<BoolValue>(CollectParetoOptimalSolutionsParameterName,
+        "Optionally collect a set of Pareto-optimal solutions minimizing error and complexity.", new BoolValue(false)));
+      Parameters[CollectParetoOptimalSolutionsParameterName].Hidden = true;
+
+      Parameters.Add(new FixedValueParameter<DoubleValue>(LambdaParameterName,
+        "Lambda is the factor for the regularization term in the objective function (Obj = (y - f(x,p))² + lambda * |p|²)", new DoubleValue(0.0)));
     }
 
@@ -195 +218 @@
     }
 
+    // TODO: support pause and restart
     protected override void Run(CancellationToken cancellationToken) {
       // Set up the algorithm
       if (SetSeedRandomly) Seed = new System.Random().Next();
+      var collectPareto = CollectParetoOptimalSolutions;
 
       // Set up the results display
@@ -229 +254 @@
       var gradEvals = new IntValue();
       Results.Add(new Result("Gradient evaluations", gradEvals));
-      var paretoBestModelsResult = new Result("ParetoBestModels", typeof(ItemList<ISymbolicRegressionSolution>));
-      Results.Add(paretoBestModelsResult);
-
+
+      Result paretoBestModelsResult = new Result("ParetoBestModels", typeof(ItemList<ISymbolicRegressionSolution>));
+      if (collectPareto) {
+        Results.Add(paretoBestModelsResult);
+      }
 
       // same as in SymbolicRegressionSingleObjectiveProblem
@@ -246 +273 @@
       var problemData = (IRegressionProblemData)Problem.ProblemData.Clone();
       if (!AllowedFactors.CheckedItems.Any()) throw new ArgumentException("At least on type of factor must be allowed");
-      var state = MctsSymbolicRegressionStatic.CreateState(problemData, (uint)Seed, MaxVariableReferences, ScaleVariables, ConstantOptimizationIterations,
-        Policy,
+      var state = MctsSymbolicRegressionStatic.CreateState(problemData, (uint)Seed, MaxVariableReferences, ScaleVariables, 
+        ConstantOptimizationIterations, Lambda,
+        Policy, collectPareto,
         lowerLimit, upperLimit,
         allowProdOfVars: AllowedFactors.CheckedItems.Any(s => s.Value.Value == VariableProductFactorName),
@@ -261 +289 @@
       double curBestQ = 0.0;
       int n = 0;
+
+      // cancelled before we acutally started
+      cancellationToken.ThrowIfCancellationRequested();
+
       // Loop until iteration limit reached or canceled.
-      for (int i = 0; i < Iterations && !state.Done; i++) {
-        cancellationToken.ThrowIfCancellationRequested();
-
+      for (int i = 0; i < Iterations && !state.Done && !cancellationToken.IsCancellationRequested; i++) {
         var q = MctsSymbolicRegressionStatic.MakeStep(state);
         sumQ += q; // sum of qs in the last updateinterval iterations
@@ -286 +316 @@
           totalRollouts.Value = state.TotalRollouts;
 
-          paretoBestModelsResult.Value = new ItemList<ISymbolicRegressionSolution>(state.ParetoBestModels);
+          if (collectPareto) {
+            paretoBestModelsResult.Value = new ItemList<ISymbolicRegressionSolution>(state.ParetoBestModels);
+          }
 
           table.Rows["Best quality"].Values.Add(bestQuality.Value);
@@ -296 +328 @@
       }
 
-      // final results
+      // final results (assumes that at least one iteration was calculated)
       if (n > 0) {
         if (bestQ > bestQuality.Value) {

branches/MCTS-SymbReg-2796/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/MctsSymbolicRegressionStatic.cs

@@ -35 +35 @@
 namespace HeuristicLab.Algorithms.DataAnalysis.MctsSymbolicRegression {
   public static class MctsSymbolicRegressionStatic {
-    // TODO: SGD with adagrad instead of lbfgs?
-    // TODO: check Taylor expansion capabilities (ln(x), sqrt(x), exp(x)) in combination with GBT
-    // TODO: optimize for 3 targets concurrently (y, 1/y, exp(y), and log(y))? Would simplify the number of possible expressions again
+    // OBJECTIVES:
+    // 1) solve toy problems without numeric constants (to show that structure search is effective / efficient)
+    //    - e.g. Keijzer, Nguyen ... where no numeric constants are involved
+    //    - assumptions: 
+    //      - we don't know the necessary operations or functions -> all available functions could be necessary
+    //      - but we do not need to tune numeric constants -> no scaling of input variables x!
+    // 2) Solve toy problems with numeric constants to make the algorithm invariant concerning variable scale.
+    //    This is important for real world applications.
+    //    - e.g. Korns or Vladislavleva problems where numeric constants are involved
+    //    - assumptions:
+    //      - any numeric constant is possible (a-priori we might assume that small abs. constants are more likely)
+    //      - standardization of variables is possible (or might be necessary) as we adjust numeric parameters of the expression anyway
+    //      - to simplify the problem we can restrict the set of functions e.g. we assume which functions are necessary for the problem instance
+    //        -> several steps: (a) polyinomials, (b) rational polynomials, (c) exponential or logarithmic functions, rational functions with exponential and logarithmic parts
+    // 3) efficiency and effectiveness for real-world problems
+    //    - e.g. Tower problem
+    //    - (1) and (2) combined, structure search must be effective in combination with numeric optimization of constants
+    //    
+
+    // TODO NEXT: check if transformation of y is correct and works
+    // TODO: transform target y to zero-mean and remove linear scaling parameters
+    // TODO: include offset for variables as parameter
+    // TODO: why does LM optimization converge so slowly with exp(x), log(x), and 1/x allowed?
+    // TODO: support e(-x) and possibly (1/-x)
+    // TODO: is it OK to initialize all constants to 1?
     #region static API
 
@@ -76 +98 @@
       private readonly double[] scalingFactor;
       private readonly double[] scalingOffset;
+      private readonly double yStdDev; // for scaling parameters (e.g. stopping condition for LM)
       private readonly int constOptIterations;
+      private readonly double lambda; // weight of penalty term for regularization
       private readonly double lowerEstimationLimit, upperEstimationLimit;
+      private readonly bool collectParetoOptimalModels;
       private readonly List<ISymbolicRegressionSolution> paretoBestModels = new List<ISymbolicRegressionSolution>();
       private readonly List<double[]> paretoFront = new List<double[]>(); // matching the models
@@ -99 +124 @@
       private readonly double[][] gradBuf;
 
-      public State(IRegressionProblemData problemData, uint randSeed, int maxVariables, bool scaleVariables, int constOptIterations,
+      public State(IRegressionProblemData problemData, uint randSeed, int maxVariables, bool scaleVariables,
+        int constOptIterations, double lambda,
         IPolicy treePolicy = null,
+        bool collectParetoOptimalModels = false,
         double lowerEstimationLimit = double.MinValue, double upperEstimationLimit = double.MaxValue,
         bool allowProdOfVars = true,
@@ -108 +135 @@
         bool allowMultipleTerms = false) {
 
+        if (lambda < 0) throw new ArgumentException("Lambda must be larger or equal zero", "lambda");
+
         this.problemData = problemData;
         this.constOptIterations = constOptIterations;
+        this.lambda = lambda;
         this.evalFun = this.Eval;
         this.lowerEstimationLimit = lowerEstimationLimit;
         this.upperEstimationLimit = upperEstimationLimit;
+        this.collectParetoOptimalModels = collectParetoOptimalModels;
 
         random = new MersenneTwister(randSeed);
@@ -128 +159 @@
         this.x = x;
         this.y = y;
+        this.yStdDev = HeuristicLab.Common.EnumerableStatisticExtensions.StandardDeviation(y);
         this.testX = testX;
         this.testY = testY;
@@ -181 +213 @@
           var t = new SymbolicExpressionTree(treeGen.Exec(bestCode, bestConsts, bestNParams, scalingFactor, scalingOffset));
           var model = new SymbolicRegressionModel(problemData.TargetVariable, t, interpreter, lowerEstimationLimit, upperEstimationLimit);
-
-          // model has already been scaled linearly in Eval
+          model.Scale(problemData); // apply linear scaling
           return model;
         }
@@ -212 +243 @@
           Array.Copy(optConsts, bestConsts, bestNParams);
         }
-        // multi-objective best
-        var complexity = // SymbolicDataAnalysisModelComplexityCalculator.CalculateComplexity() TODO
-          Array.FindIndex(code, (opc)=>opc==(byte)OpCodes.Exit);
-        UpdateParetoFront(q, complexity, code, optConsts, nParams, scalingFactor, scalingOffset);
-
+        if (collectParetoOptimalModels) {
+          // multi-objective best
+          var complexity = // SymbolicDataAnalysisModelComplexityCalculator.CalculateComplexity() TODO: implement Kommenda's tree complexity directly in the evaluator
+            Array.FindIndex(code, (opc) => opc == (byte)OpCodes.Exit);  // use length of expression as surrogate for complexity
+          UpdateParetoFront(q, complexity, code, optConsts, nParams, scalingFactor, scalingOffset);
+        }
         return q;
       }
+
+      private void Eval(byte[] code, int nParams, out double rsq, out double[] optConsts) {
+        // we make a first pass to determine a valid starting configuration for all constants
+        // constant c in log(c + f(x)) is adjusted to guarantee that x is positive (see expression evaluator)
+        // scale and offset are set to optimal starting configuration
+        // assumes scale is the first param and offset is the last param
+
+        // reset constants
+        Array.Copy(ones, constsBuf, nParams);
+        evaluator.Exec(code, x, constsBuf, predBuf, adjustOffsetForLogAndExp: true);
+        funcEvaluations++;
+
+        if (nParams == 0 || constOptIterations < 0) {
+          // if we don't need to optimize parameters then we are done
+          // changing scale and offset does not influence r²
+          rsq = RSq(y, predBuf);
+          optConsts = constsBuf;
+        } else {
+          // optimize constants using the starting point calculated above
+          OptimizeConstsLm(code, constsBuf, nParams, 0.0, nIters: constOptIterations);
+
+          evaluator.Exec(code, x, constsBuf, predBuf);
+          funcEvaluations++;
+
+          rsq = RSq(y, predBuf);
+          optConsts = constsBuf;
+        }
+      }
+
+
+
+      #region helpers
+      private static double RSq(IEnumerable<double> x, IEnumerable<double> y) {
+        OnlineCalculatorError error;
+        double r = OnlinePearsonsRCalculator.Calculate(x, y, out error);
+        return error == OnlineCalculatorError.None ? r * r : 0.0;
+      }
+
+
+      private void OptimizeConstsLm(byte[] code, double[] consts, int nParams, double epsF = 0.0, int nIters = 100) {
+        double[] optConsts = new double[nParams]; // allocate a smaller buffer for constants opt (TODO perf?)
+        Array.Copy(consts, optConsts, nParams);
+
+        // direct usage of LM is recommended in alglib manual for better performance than the lsfit interface (which uses lm internally).
+        alglib.minlmstate state;
+        alglib.minlmreport rep = null;
+        alglib.minlmcreatevj(y.Length + 1, optConsts, out state); // +1 for penalty term
+        // Using the change of the gradient as stopping criterion is recommended in alglib manual.
+        // However, the most recent version of alglib (as of Oct 2017) only supports epsX as stopping criterion
+        alglib.minlmsetcond(state, epsg: 1E-6 * yStdDev, epsf: epsF, epsx: 0.0, maxits: nIters);
+        // alglib.minlmsetgradientcheck(state, 1E-5);
+        alglib.minlmoptimize(state, Func, FuncAndJacobian, null, code);
+        alglib.minlmresults(state, out optConsts, out rep);
+        funcEvaluations += rep.nfunc;
+        gradEvaluations += rep.njac * nParams;
+
+        if (rep.terminationtype < 0) throw new ArgumentException("lm failed: termination type = " + rep.terminationtype);
+
+        // only use optimized constants if successful
+        if (rep.terminationtype >= 0) {
+          Array.Copy(optConsts, consts, optConsts.Length);
+        }
+      }
+
+      private void Func(double[] arg, double[] fi, object obj) {
+        var code = (byte[])obj;
+        int n = predBuf.Length;
+        evaluator.Exec(code, x, arg, predBuf); // gradients are nParams x vLen
+        for (int r = 0; r < n; r++) {
+          var res = predBuf[r] - y[r];
+          fi[r] = res;
+        }
+
+        var penaltyIdx = fi.Length - 1;
+        fi[penaltyIdx] = 0.0;
+        // calc length of parameter vector for regularization
+        var aa = 0.0;
+        for (int i = 0; i < arg.Length; i++) {
+          aa += arg[i] * arg[i];
+        }
+        if (lambda > 0 && aa > 0) {
+          // scale lambda using stdDev(y) to make the parameter independent of the scale of y
+          // scale lambda using n to make parameter independent of the number of training points
+          // take the root because LM squares the result
+          fi[penaltyIdx] = Math.Sqrt(n * lambda / yStdDev * aa);
+        }
+      }
+
+      private void FuncAndJacobian(double[] arg, double[] fi, double[,] jac, object obj) {
+        int n = predBuf.Length;
+        int nParams = arg.Length;
+        var code = (byte[])obj;
+        evaluator.ExecGradient(code, x, arg, predBuf, gradBuf); // gradients are nParams x vLen
+        for (int r = 0; r < n; r++) {
+          var res = predBuf[r] - y[r];
+          fi[r] = res;
+
+          for (int k = 0; k < nParams; k++) {
+            jac[r, k] = gradBuf[k][r];
+          }
+        }
+        // calc length of parameter vector for regularization
+        double aa = 0.0;
+        for (int i = 0; i < arg.Length; i++) {
+          aa += arg[i] * arg[i];
+        }
+
+        var penaltyIdx = fi.Length - 1;
+        if (lambda > 0 && aa > 0) {
+          fi[penaltyIdx] = 0.0;
+          // scale lambda using stdDev(y) to make the parameter independent of the scale of y
+          // scale lambda using n to make parameter independent of the number of training points
+          // take the root because alglib LM squares the result
+          fi[penaltyIdx] = Math.Sqrt(n * lambda / yStdDev * aa);
+
+          for (int i = 0; i < arg.Length; i++) {
+            jac[penaltyIdx, i] = 0.5 / fi[penaltyIdx] * 2 * n * lambda / yStdDev * arg[i];
+          }
+        } else {
+          fi[penaltyIdx] = 0.0;
+          for (int i = 0; i < arg.Length; i++) {
+            jac[penaltyIdx, i] = 0.0;
+          }
+        }
+      }
+
 
       private void UpdateParetoFront(double q, int complexity, byte[] code, double[] param, int nParam,
@@ -233 +391 @@
           }
         }
-        if(isNonDominated) {
+        if (isNonDominated) {
           paretoFront.Add(cur);
 
@@ -242 +400 @@
           var t = new SymbolicExpressionTree(treeGen.Exec(code, param, nParam, scalingFactor, scalingOffset));
           var model = new SymbolicRegressionModel(problemData.TargetVariable, t, interpreter, lowerEstimationLimit, upperEstimationLimit);
+          model.Scale(problemData); // apply linear scaling
 
           var sol = model.CreateRegressionSolution(this.problemData);
           sol.Name = string.Format("{0:N5} {1}", q, complexity);
-          
+
           paretoBestModels.Add(sol);
         }
-        for(int i=paretoFront.Count-2;i>=0;i--) {
+        for (int i = paretoFront.Count - 2; i >= 0; i--) {
           var @ref = paretoFront[i];
           var domRes = DominationCalculator<int>.Dominates(cur, @ref, max, true);
-          if(domRes == DominationResult.Dominates) {
+          if (domRes == DominationResult.Dominates) {
             paretoFront.RemoveAt(i);
             paretoBestModels.RemoveAt(i);
@@ -257 +416 @@
         }
       }
-
-      private void Eval(byte[] code, int nParams, out double rsq, out double[] optConsts) {
-        // we make a first pass to determine a valid starting configuration for all constants
-        // constant c in log(c + f(x)) is adjusted to guarantee that x is positive (see expression evaluator)
-        // scale and offset are set to optimal starting configuration
-        // assumes scale is the first param and offset is the last param
-        double alpha;
-        double beta;
-
-        // reset constants
-        Array.Copy(ones, constsBuf, nParams);
-        evaluator.Exec(code, x, constsBuf, predBuf, adjustOffsetForLogAndExp: true);
-        funcEvaluations++;
-
-        // calc opt scaling (alpha*f(x) + beta)
-        OnlineCalculatorError error;
-        OnlineLinearScalingParameterCalculator.Calculate(predBuf, y, out alpha, out beta, out error);
-        if (error == OnlineCalculatorError.None) {
-          constsBuf[0] *= beta;
-          constsBuf[nParams - 1] = constsBuf[nParams - 1] * beta + alpha;
-        }
-        if (nParams <= 2 || constOptIterations <= 0) {
-          // if we don't need to optimize parameters then we are done
-          // changing scale and offset does not influence r²
-          rsq = RSq(y, predBuf);
-          optConsts = constsBuf;
-        } else {
-          // optimize constants using the starting point calculated above
-          OptimizeConstsLm(code, constsBuf, nParams, 0.0, nIters: constOptIterations);
-
-          evaluator.Exec(code, x, constsBuf, predBuf);
-          funcEvaluations++;
-
-          rsq = RSq(y, predBuf);
-          optConsts = constsBuf;
-        }
-      }
-
-
-
-      #region helpers
-      private static double RSq(IEnumerable<double> x, IEnumerable<double> y) {
-        OnlineCalculatorError error;
-        double r = OnlinePearsonsRCalculator.Calculate(x, y, out error);
-        return error == OnlineCalculatorError.None ? r * r : 0.0;
-      }
-
-
-      private void OptimizeConstsLm(byte[] code, double[] consts, int nParams, double epsF = 0.0, int nIters = 100) {
-        double[] optConsts = new double[nParams]; // allocate a smaller buffer for constants opt (TODO perf?)
-        Array.Copy(consts, optConsts, nParams);
-
-        alglib.minlmstate state;
-        alglib.minlmreport rep = null;
-        alglib.minlmcreatevj(y.Length, optConsts, out state);
-        alglib.minlmsetcond(state, 0.0, epsF, 0.0, nIters);
-        //alglib.minlmsetgradientcheck(state, 0.000001);
-        alglib.minlmoptimize(state, Func, FuncAndJacobian, null, code);
-        alglib.minlmresults(state, out optConsts, out rep);
-        funcEvaluations += rep.nfunc;
-        gradEvaluations += rep.njac * nParams;
-
-        if (rep.terminationtype < 0) throw new ArgumentException("lm failed: termination type = " + rep.terminationtype);
-
-        // only use optimized constants if successful
-        if (rep.terminationtype >= 0) {
-          Array.Copy(optConsts, consts, optConsts.Length);
-        }
-      }
-
-      private void Func(double[] arg, double[] fi, object obj) {
-        var code = (byte[])obj;
-        evaluator.Exec(code, x, arg, predBuf); // gradients are nParams x vLen
-        for (int r = 0; r < predBuf.Length; r++) {
-          var res = predBuf[r] - y[r];
-          fi[r] = res;
-        }
-      }
-      private void FuncAndJacobian(double[] arg, double[] fi, double[,] jac, object obj) {
-        int nParams = arg.Length;
-        var code = (byte[])obj;
-        evaluator.ExecGradient(code, x, arg, predBuf, gradBuf); // gradients are nParams x vLen
-        for (int r = 0; r < predBuf.Length; r++) {
-          var res = predBuf[r] - y[r];
-          fi[r] = res;
-
-          for (int k = 0; k < nParams; k++) {
-            jac[r, k] = gradBuf[k][r];
-          }
-        }
-      }
       #endregion
     }
 
+
+    /// <summary>
+    /// Static method to initialize a state for the algorithm
+    /// </summary>
+    /// <param name="problemData">The problem data</param>
+    /// <param name="randSeed">Random seed.</param>
+    /// <param name="maxVariables">Maximum number of variable references that are allowed in the expression.</param>
+    /// <param name="scaleVariables">Optionally scale input variables to the interval [0..1] (recommended)</param>
+    /// <param name="constOptIterations">Maximum number of iterations for constants optimization (Levenberg-Marquardt)</param>
+    /// <param name="lambda">Penalty factor for regularization (0..inf.), small penalty disabled regularization.</param>
+    /// <param name="policy">Tree search policy (random, ucb, eps-greedy, ...)</param>
+    /// <param name="collectParameterOptimalModels">Optionally collect all Pareto-optimal solutions having minimal length and error.</param>
+    /// <param name="lowerEstimationLimit">Optionally limit the result of the expression to this lower value.</param>
+    /// <param name="upperEstimationLimit">Optionally limit the result of the expression to this upper value.</param>
+    /// <param name="allowProdOfVars">Allow products of expressions.</param>
+    /// <param name="allowExp">Allow expressions with exponentials.</param>
+    /// <param name="allowLog">Allow expressions with logarithms</param>
+    /// <param name="allowInv">Allow expressions with 1/x</param>
+    /// <param name="allowMultipleTerms">Allow expressions which are sums of multiple terms.</param>
+    /// <returns></returns>
+
     public static IState CreateState(IRegressionProblemData problemData, uint randSeed, int maxVariables = 3,
-      bool scaleVariables = true, int constOptIterations = 0,
+      bool scaleVariables = true, int constOptIterations = -1, double lambda = 0.0,
       IPolicy policy = null,
+      bool collectParameterOptimalModels = false,
       double lowerEstimationLimit = double.MinValue, double upperEstimationLimit = double.MaxValue,
       bool allowProdOfVars = true,
@@ -361 +451 @@
       bool allowMultipleTerms = false
       ) {
-      return new State(problemData, randSeed, maxVariables, scaleVariables, constOptIterations,
-        policy,
+      return new State(problemData, randSeed, maxVariables, scaleVariables, constOptIterations, lambda,
+        policy, collectParameterOptimalModels,
         lowerEstimationLimit, upperEstimationLimit,
         allowProdOfVars, allowExp, allowLog, allowInv, allowMultipleTerms);
@@ -481 +571 @@
       var i = 0;
       if (scaleVariables) {
-        scalingFactor = new double[xs.Length];
-        scalingOffset = new double[xs.Length];
+        scalingFactor = new double[xs.Length + 1];
+        scalingOffset = new double[xs.Length + 1];
       } else {
         scalingFactor = null;
@@ -502 +592 @@
       }
 
+      if (scaleVariables) {
+        // transform target variable to zero-mean
+        scalingFactor[i] = 1.0;
+        scalingOffset[i] = -problemData.Dataset.GetDoubleValues(problemData.TargetVariable, rows).Average();
+      }
+
       GenerateData(problemData, rows, scalingFactor, scalingOffset, out xs, out y);
     }
@@ -518 +614 @@
       }
 
-      y = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, rows).ToArray();
+      {
+        var sf = scalingFactor == null ? 1.0 : scalingFactor[i];
+        var offset = scalingFactor == null ? 0.0 : scalingOffset[i];
+        y = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, rows).Select(yi => yi * sf + offset).ToArray();
+      }
     }
   }