Changeset 13651 for trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression

Timestamp:

03/05/16 08:25:08 (9 years ago)

Author:

gkronber

Message:

#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

Location:

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression

Files:

• Automaton.cs (modified) (7 diffs)
• ConstraintHandler.cs (modified) (4 diffs)
• Disassembler.cs (modified) (2 diffs)
• ExpressionEvaluator.cs (modified) (7 diffs)
• MctsSymbolicRegressionAlgorithm.cs (modified) (4 diffs)
• MctsSymbolicRegressionStatic.cs (modified) (12 diffs)

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/Automaton.cs

@@ -26 +26 @@
 namespace HeuristicLab.Algorithms.DataAnalysis.MctsSymbolicRegression {
   // this is the core class for generating expressions. 
-  // the automaton determines which expressions are allowed
+  // it represents a finite state automaton, each state transition can be associated with an action (e.g. to produce code).
+  // the automaton determines the possible structures for expressions.
+  // 
+  // to understand this code it is worthwile to generate a graphical visualization of the automaton (see PrintAutomaton).
+  // If the code is compiled in debug mode the automaton produces a Graphviz file into the folder of the application
+  // whenever an instance of the automaton is constructed.
+  //
+  // This class relies on two other classes: 
+  // - CodeGenerator to produce code for a stack-based evaluator and
+  // - ConstraintHandler to restrict the allowed set of expressions. 
+  //
+  // The ConstraintHandler extends the automaton and adds semantic restrictions for expressions produced by the automaton.
+  // 
+  // 
   internal class Automaton {
     public const int StateExpr = 1;
@@ -52 +65 @@
     public const int StateInvTFStart = 23;
     public const int StateInvTFEnd = 24;
-    private const int FirstDynamicState = 25;
+    public const int FirstDynamicState = 25;
+    // more states for individual variables are created dynamically
 
     private const int StartState = StateExpr;
@@ -221 +235 @@
         () => {
           codeGenerator.Emit1(OpCodes.LoadConst0);
-        },
-        "0");
+          constraintHandler.StartNewTermInPoly();
+        },
+        "0, StartTermInPoly");
       AddTransition(StateLogTEnd, StateLogFactorEnd,
         () => {
@@ -271 +286 @@
         () => {
           codeGenerator.Emit1(OpCodes.LoadConst1);
-        },
-        "c");
+          constraintHandler.StartNewTermInPoly();
+        },
+        "c, StartTermInPoly");
       AddTransition(StateInvTEnd, StateInvFactorEnd,
         () => {
@@ -337 +353 @@
     private readonly int[] followStatesBuf = new int[1000];
     public void FollowStates(int state, out int[] buf, out int nElements) {
-      // return followStates[state]
-      //   .Where(s => s < FirstDynamicState || s >= minVarIdx) // for variables we only allow non-decreasing state sequences
-      //   // the following states imply an additional variable being added to the expression
-      //   // F, Sum, Prod
-      //   .Where(s => (s != StateF && s != StateSum && s != StateProd) || variablesRemaining > 0);
-
       // for loop instead of where iterator
       var fs = followStates[state];
       int j = 0;
-      //Console.Write(stateNames[CurrentState] + " allowed: ");
       for (int i = 0; i < fs.Count; i++) {
         var s = fs[i];
         if (constraintHandler.IsAllowedFollowState(state, s)) {
-          //Console.Write(s + " ");
           followStatesBuf[j++] = s;
         }
       }
-      //Console.WriteLine();
       buf = followStatesBuf;
       nElements = j;
@@ -361 +368 @@
 
     public void Goto(int targetState) {
-      //Console.WriteLine("->{0}", stateNames[targetState]);
-      // Contract.Assert(FollowStates(CurrentState).Contains(targetState));
-
       if (actions[CurrentState, targetState] != null)
         actions[CurrentState, targetState].ForEach(a => a()); // execute all actions
@@ -370 +374 @@
 
     public bool IsFinalState(int s) {
-      return s == StateExprEnd;
+      return s == StateExprEnd && !constraintHandler.IsInvalidExpression;
     }
 

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/ConstraintHandler.cs

@@ -2 +2 @@
 /* HeuristicLab
  * Copyright (C) 2002-2015 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
- * and the BEACON Center for the Study of Evolution in Action.
  * 
  * This file is part of HeuristicLab.
@@ -21 +20 @@
 #endregion
 
-
+using System;
+using System.Collections.Generic;
 using System.Diagnostics.Contracts;
+using System.Linq;
 
 namespace HeuristicLab.Algorithms.DataAnalysis.MctsSymbolicRegression {
 
-  // more states for individual variables are created dynamically
+  // 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.
+  // The ratio of the number of unsuccessful searches, that need backtracking should be tracked in the MCTS alg (MctsSymbolicRegressionStatic)
+
+  // All changes to this class should be tested through unit tests. It is important that the ConstraintHandler is not too restrictive.
+
+  // the constraints are derived from a canonical form for expressions.
+  // overall we can enforce a limited number of variable references
+  // 
+  // an expression is a sum of terms t_1 ... t_n where terms are ordered according to a relation t_i (<=)_term t_j for each pair t_i, t_j and i <= j
+  // a term is a product of factors where factors are ordered according to relation f_i (<=)_factor f_j for each pair f_i,f_j and i <= j
+
+  // we want to enforce lower-order terms before higher-order terms in expressions (based on number of variable references)
+  // factors can have different types (variable, exp, log, inverse)
+
+  // (<=)_term  [IsSmallerOrEqualTerm(t_i, t_j)]
+  //   1.  NumberOfVarRefs(t_i) < NumberOfVarRefs(t_j)  --> true           enforce terms with non-decreasing number of var refs
+  //   2.  NumberOfVarRefs(t_i) > NumberOfVarRefs(t_j)  --> false
+  //   3.  NumFactors(t_i) > NumFactors(t_j)            --> true           enforce terms with non-increasing number of factors
+  //   4.  NumFactors(t_i) < NumFactors(t_j)            --> false
+  //   5.  for all k factors: Factor(k, t_i) (<=)_factor  Factor(k, t_j) --> true // factors must be non-decreasing
+  //   6.  all factors are (=)_factor                   --> true
+  //   7.  else false
+
+  // (<=)_factor  [IsSmallerOrEqualFactor(f_i, f_j)]
+  //   1.  FactorType(t_i) < FactorType(t_j)  --> true           enforce terms with non-decreasing factor type (var < exp < log < inv)
+  //   2.  FactorType(t_i) > FactorType(t_j)  --> false
+  //   3.  Compare the two factors specifically
+  //     - variables: varIdx_i <= varIdx_j (only one var reference)
+  //     - exp: number of variable references and then varIdx_i <= varIdx_j for each position
+  //     - log: number of variable references and ...
+  //     - inv: number of variable references and ...
+  //
+
+  // for log and inverse factors we allow all polynomials as argument
+  // a polynomial is a sum of terms t_1 ... t_n where terms are ordered according to a relation t_i (<=)_poly t_j for each pair t_i, t_j and i <= j
+
+  // (<=)_poly  [IsSmallerOrEqualPoly(t_i, t_j)]
+  //  1. NumberOfVarRefs(t_i) < NumberOfVarRefs(t_j)         --> true // enforce non-decreasing number of var refs
+  //  2. NumberOfVarRefs(t_i) > NumberOfVarRefs(t_j)         --> false // enforce non-decreasing number of var refs
+  //  3. for all k variables: VarIdx(k,t_i) > VarIdx(k, t_j) --> false // enforce non-decreasing variable idx
+
+
+  // we store the following to make comparsions:
+  // - prevTerm (complete & containing all factors)
+  // - curTerm  (incomplete & containing all completed factors)
+  // - curFactor (incomplete)
   internal class ConstraintHandler {
     private int nVars;
     private readonly int maxVariables;
-
-    public int prevTermFirstVariableState;
-    public int curTermFirstVariableState;
-    public int prevTermFirstFactorType;
-    public int curTermFirstFactorType;
-    public int prevFactorType;
-    public int curFactorType;
-    public int prevFactorFirstVariableState;
-    public int curFactorFirstVariableState;
-    public int prevVariableRef;
+    private bool invalidExpression;
+
+    public bool IsInvalidExpression {
+      get { return invalidExpression; }
+    }
+
+
+    private TermInformation prevTerm;
+    private TermInformation curTerm;
+    private FactorInformation curFactor;
+
+
+    private class TermInformation {
+      public int numVarReferences { get { return factors.Sum(f => f.numVarReferences); } }
+      public List<FactorInformation> factors = new List<FactorInformation>();
+    }
+
+    private class FactorInformation {
+      public int numVarReferences = 0;
+      public int factorType; // use the state number to represent types
+
+      // for variable factors
+      public int variableState = -1;
+
+      // for exp  factors
+      public List<int> expVariableStates = new List<int>();
+
+      // for log and inv factors
+      public List<List<int>> polyVariableStates = new List<List<int>>();
+    }
 
 
@@ -46 +116 @@
     }
 
-    // 1) an expression is a sum of terms t_1 ... t_n
-    //    FirstFactorType(t_i) <= FirstFactorType(t_j) for each pair t_i, t_j where i < j 
-    //    FirstVarReference(t_i) <= FirstVarReference(t_j) for each pair t_i, t_j where i < j and FirstFactorType(t_i) = FirstFactorType(t_j)
-    // 2) a term is a product of factors, each factor is either a variable factor, an exp factor, a log factor or an inverse factor
-    //    FactorType(f_i) <= FactorType(f_j) for each pair of factors f_i, f_j and i < j
-    //    FirstVarReference(f_i) <= FirstVarReference(f_j) for each pair of factors f_i, f_j and i < j and FactorType(f_i) = FactorType(f_j)
-    // 3) a variable factor is a product of variable references v1...vn
-    //    VarIdx(v_i) <= VarIdx(v_j) for each pair of variable references v_i, v_j and i < j
-    //    (IMPLICIT) FirstVarReference(t) <= VarIdx(v_i) for each variable reference v_i in term t
-    // 4) an exponential factor is the exponential of a product of variables v1...vn
-    //    VarIdx(v_i) <= VarIdx(v_j) for each pair of variable references v_i, v_j and i < j
-    //    (IMPLICIT) FirstVarReference(t) <= VarIdx(v_i) for each variable reference v_i in term t
-    // 5) a log factor is a sum of terms t_i where each term is a product of variables
-    //    FirstVarReference(t_i) <= FirstVarReference(t_j) for each pair of terms t_i, t_j and i < j 
-    //    for each term t: VarIdx(v_i) <= VarIdx(v_j) for each pair of variable references v_i, v_j and i < j in t
+    // the order relations for terms and factors
+
+    private static int CompareTerms(TermInformation a, TermInformation b) {
+      if (a.numVarReferences < b.numVarReferences) return -1;
+      if (a.numVarReferences > b.numVarReferences) return 1;
+
+      if (a.factors.Count > b.factors.Count) return -1;  // terms with more factors should be ordered first
+      if (a.factors.Count < b.factors.Count) return +1;
+
+      var aFactors = a.factors.GetEnumerator();
+      var bFactors = b.factors.GetEnumerator();
+      while (aFactors.MoveNext() & bFactors.MoveNext()) {
+        var c = CompareFactors(aFactors.Current, bFactors.Current);
+        if (c < 0) return -1;
+        if (c > 0) return 1;
+      }
+      // all factors are the same => terms are the same
+      return 0;
+    }
+
+    private static int CompareFactors(FactorInformation a, FactorInformation b) {
+      if (a.factorType < b.factorType) return -1;
+      if (a.factorType > b.factorType) return +1;
+      // same factor types
+      if (a.factorType == Automaton.StateVariableFactorStart) {
+        return a.variableState.CompareTo(b.variableState);
+      } else if (a.factorType == Automaton.StateExpFactorStart) {
+        return CompareStateLists(a.expVariableStates, b.expVariableStates);
+      } else {
+        if (a.numVarReferences < b.numVarReferences) return -1;
+        if (a.numVarReferences > b.numVarReferences) return +1;
+        if (a.polyVariableStates.Count > b.polyVariableStates.Count) return -1; // more terms in the poly should be ordered first
+        if (a.polyVariableStates.Count < b.polyVariableStates.Count) return +1;
+        // log and inv
+        var aTerms = a.polyVariableStates.GetEnumerator();
+        var bTerms = b.polyVariableStates.GetEnumerator();
+        while (aTerms.MoveNext() & bTerms.MoveNext()) {
+          var c = CompareStateLists(aTerms.Current, bTerms.Current);
+          if (c != 0) return c;
+        }
+        return 0; // all terms in the polynomial are the same
+      }
+    }
+
+    private static int CompareStateLists(List<int> a, List<int> b) {
+      if (a.Count < b.Count) return -1;
+      if (a.Count > b.Count) return +1;
+      for (int i = 0; i < a.Count; i++) {
+        if (a[i] < b[i]) return -1;
+        if (a[i] > b[i]) return +1;
+      }
+      return 0; // all states are the same
+    }
+
+
+    private bool IsNewTermAllowed() {
+      // next term must have at least as many variable references as the previous term
+      return prevTerm == null || nVars + prevTerm.numVarReferences <= maxVariables;
+    }
+
+    private bool IsNewFactorAllowed() {
+      // next factor must have a larger or equal type compared to the previous factor.
+      // if the types are the same it must have at least as many variable references.
+      // so if the prevFactor is any other than invFactor (last possible type) then we only need to be able to add one variable
+      // otherwise we need to be able to add at least as many variables as the previous factor
+      return !curTerm.factors.Any() ||
+             (nVars + curTerm.factors.Last().numVarReferences <= maxVariables);
+    }
+
+    private bool IsAllowedAsNextFactorType(int followState) {
+      // IsNewTermAllowed already ensures that we can add a term with enough variable references
+
+      // enforce constraints within terms (compare to prev factor)
+      if (curTerm.factors.Any()) {
+        // enforce non-decreasing factor types
+        if (curTerm.factors.Last().factorType > followState) return false;
+        // when the factor type is the same, starting a new factor is only allowed if we can add at least the number of variables of the prev factor
+        if (curTerm.factors.Last().factorType == followState && nVars + curTerm.factors.Last().numVarReferences > maxVariables) return false;
+      }
+
+      // enforce constraints on terms (compare to prev term)
+      // meaning that we must ensure non-decreasing terms
+      if (prevTerm != null) {
+        // a factor type is only allowed if we can then produce a term that is larger or equal to the prev term
+        // (1) if we the number of variable references still remaining is larger than the number of variable references in the prev term 
+        //     then it is always possible to build a larger term
+        // (2) otherwise we try to build the largest possible term starting from current factors in the term.
+        //     
+
+        var numVarRefsRemaining = maxVariables - nVars;
+        Contract.Assert(!curTerm.factors.Any() || curTerm.factors.Last().numVarReferences <= numVarRefsRemaining);
+
+        if (prevTerm.numVarReferences < numVarRefsRemaining) return true;
+
+        // variable factors must be handled differently because they can only contain one variable reference
+        if (followState == Automaton.StateVariableFactorStart) {
+          // append the variable factor and the maximum possible state from the previous factor to create a larger factor
+          var varF = CreateLargestPossibleFactor(Automaton.StateVariableFactorStart, 1);
+          var maxF = CreateLargestPossibleFactor(prevTerm.factors.Max(f => f.factorType), numVarRefsRemaining - 1);
+          var origFactorCount = curTerm.factors.Count;
+          // add this factor to the current term
+          curTerm.factors.Add(varF);
+          curTerm.factors.Add(maxF);
+          var c = CompareTerms(prevTerm, curTerm);
+          // restore term
+          curTerm.factors.RemoveRange(origFactorCount, 2);
+          // if the prev term is still larger then this followstate is not allowed
+          if (c > 0) {
+            return false;
+          }
+        } else {
+          var newF = CreateLargestPossibleFactor(followState, numVarRefsRemaining);
+
+          var origFactorCount = curTerm.factors.Count;
+          // add this factor to the current term
+          curTerm.factors.Add(newF);
+          var c = CompareTerms(prevTerm, curTerm);
+          // restore term
+          curTerm.factors.RemoveAt(origFactorCount);
+          // if the prev term is still larger then this followstate is not allowed
+          if (c > 0) {
+            return false;
+          }
+        }
+      }
+      return true;
+    }
+
+    // largest possible factor of the given kind
+    private FactorInformation CreateLargestPossibleFactor(int factorType, int numVarRefs) {
+      var newF = new FactorInformation();
+      newF.factorType = factorType;
+      if (factorType == Automaton.StateVariableFactorStart) {
+        newF.variableState = int.MaxValue;
+        newF.numVarReferences = 1;
+      } else if (factorType == Automaton.StateExpFactorStart) {
+        for (int i = 0; i < numVarRefs; i++)
+          newF.expVariableStates.Add(int.MaxValue);
+        newF.numVarReferences = numVarRefs;
+      } else if (factorType == Automaton.StateInvFactorStart || factorType == Automaton.StateLogFactorStart) {
+        for (int i = 0; i < numVarRefs; i++) {
+          newF.polyVariableStates.Add(new List<int>());
+          newF.polyVariableStates[i].Add(int.MaxValue);
+        }
+        newF.numVarReferences = numVarRefs;
+      }
+      return newF;
+    }
+
+    private bool IsAllowedAsNextVariableFactor(int variableState) {
+      Contract.Assert(variableState >= Automaton.FirstDynamicState);
+      return !curTerm.factors.Any() || curTerm.factors.Last().variableState <= variableState;
+    }
+
+    private bool IsAllowedAsNextInExp(int variableState) {
+      Contract.Assert(variableState >= Automaton.FirstDynamicState);
+      if (curFactor.expVariableStates.Any() && curFactor.expVariableStates.Last() > variableState) return false;
+      if (curTerm.factors.Any()) {
+        // try and compare with prev factor      
+        curFactor.numVarReferences++;
+        curFactor.expVariableStates.Add(variableState);
+        var c = CompareFactors(curTerm.factors.Last(), curFactor);
+        curFactor.numVarReferences--;
+        curFactor.expVariableStates.RemoveAt(curFactor.expVariableStates.Count - 1);
+        return c <= 0;
+      }
+      return true;
+    }
+
+    private bool IsNewTermAllowedInPoly() {
+      return nVars + curFactor.polyVariableStates.Last().Count() <= maxVariables;
+    }
+
+    private bool IsAllowedAsNextInPoly(int variableState) {
+      Contract.Assert(variableState >= Automaton.FirstDynamicState);
+      return !curFactor.polyVariableStates.Any() ||
+             !curFactor.polyVariableStates.Last().Any() ||
+              curFactor.polyVariableStates.Last().Last() <= variableState;
+    }
+    private bool IsTermCompleteInPoly() {
+      var nTerms = curFactor.polyVariableStates.Count;
+      return nTerms == 1 ||
+             curFactor.polyVariableStates[nTerms - 2].Count <= curFactor.polyVariableStates[nTerms - 1].Count;
+
+    }
+    private bool IsCompleteExp() {
+      return !curTerm.factors.Any() || CompareFactors(curTerm.factors.Last(), curFactor) <= 0;
+    }
+
     public bool IsAllowedFollowState(int currentState, int followState) {
-      // the following states are always allowed
+      // an invalid action was taken earlier on => nothing can be done anymore
+      if (invalidExpression) return false;
+      // states that have no alternative are always allowed
+      // some ending states are only allowed if enough variables have been used in the term
       if (
-        followState == Automaton.StateVariableFactorEnd ||
-        followState == Automaton.StateExpFEnd ||
-        followState == Automaton.StateExpFactorEnd ||
-        followState == Automaton.StateLogTFEnd ||
-        followState == Automaton.StateLogTEnd ||
-        followState == Automaton.StateLogFactorEnd ||
-        followState == Automaton.StateInvTFEnd ||
-        followState == Automaton.StateInvTEnd ||
-        followState == Automaton.StateInvFactorEnd ||
-        followState == Automaton.StateFactorEnd ||
-        followState == Automaton.StateTermEnd ||
-        followState == Automaton.StateExprEnd
+        currentState == Automaton.StateTermStart ||           // no alternative
+        currentState == Automaton.StateExpFactorStart ||
+        currentState == Automaton.StateLogFactorStart ||
+        currentState == Automaton.StateInvFactorStart ||
+        followState == Automaton.StateVariableFactorEnd ||    // no alternative
+        followState == Automaton.StateExpFEnd ||              // no alternative
+        followState == Automaton.StateLogTFEnd ||             // no alternative
+        followState == Automaton.StateInvTFEnd ||             // no alternative
+        followState == Automaton.StateFactorEnd ||            // always allowed because no alternative
+        followState == Automaton.StateExprEnd                 // we could also constrain the minimum number of terms here
       ) return true;
 
 
-      // all other states are only allowed if we can add more variables
-      if (nVars >= maxVariables) return false;
-
-      // the following states are always allowed when we can add more variables
+      // starting a new term is only allowed if we can add a term with at least the number of variables of the prev term
+      if (followState == Automaton.StateTermStart && !IsNewTermAllowed()) return false;
+      if (followState == Automaton.StateFactorStart && !IsNewFactorAllowed()) return false;
+      if (currentState == Automaton.StateFactorStart && !IsAllowedAsNextFactorType(followState)) return false;
+      if (followState == Automaton.StateTermEnd && prevTerm != null && CompareTerms(prevTerm, curTerm) > 0) return false;
+
+      // all of these states add at least one variable
       if (
-        followState == Automaton.StateTermStart ||
-        followState == Automaton.StateFactorStart ||
-        followState == Automaton.StateExpFStart ||
-        followState == Automaton.StateLogTStart ||
-        followState == Automaton.StateLogTFStart ||
-        followState == Automaton.StateInvTStart ||
-        followState == Automaton.StateInvTFStart
-        ) return true;
-
-      // enforce non-decreasing factor types 
-      if (currentState == Automaton.StateFactorStart) {
-        if (curFactorType < 0) {
-          //    FirstFactorType(t_i) <= FirstFactorType(t_j) for each pair t_i, t_j where i < j
-          return prevTermFirstFactorType <= followState;
-        } else {
-          // FactorType(f_i) <= FactorType(f_j) for each pair of factors f_i, f_j and i < j
-          return curFactorType <= followState;
-        }
-      }
-      // enforce non-decreasing variables references in variable and exp factors
-      if (currentState == Automaton.StateVariableFactorStart || currentState == Automaton.StateExpFStart || currentState == Automaton.StateLogTFStart || currentState == Automaton.StateInvTFStart) {
-        if (prevVariableRef > followState) return false; // never allow decreasing variables
-        if (prevFactorType < 0) {
-          // FirstVarReference(t_i) <= FirstVarReference(t_j) for each pair t_i, t_j where i < j
-          return prevTermFirstVariableState <= followState;
-        } else if (prevFactorType == curFactorType) {
-          // (FirstVarReference(f_i) <= FirstVarReference(f_j) for each pair of factors f_i, f_j and i < j and FactorType(f_i) = FactorType(f_j)
-          return prevFactorFirstVariableState <= followState;
-        }
-      }
-
-
-      return true;
+          followState == Automaton.StateVariableFactorStart ||
+          followState == Automaton.StateExpFactorStart || followState == Automaton.StateExpFStart ||
+          followState == Automaton.StateLogFactorStart || followState == Automaton.StateLogTStart ||
+          followState == Automaton.StateLogTFStart ||
+          followState == Automaton.StateInvFactorStart || followState == Automaton.StateInvTStart ||
+          followState == Automaton.StateInvTFStart) {
+        if (nVars + 1 > maxVariables) return false;
+      }
+
+      if (currentState == Automaton.StateVariableFactorStart && !IsAllowedAsNextVariableFactor(followState)) return false;
+      else if (currentState == Automaton.StateExpFStart && !IsAllowedAsNextInExp(followState)) return false;
+      else if (followState == Automaton.StateLogTStart && !IsNewTermAllowedInPoly()) return false;
+      else if (currentState == Automaton.StateLogTFStart && !IsAllowedAsNextInPoly(followState)) return false;
+      else if (followState == Automaton.StateInvTStart && !IsNewTermAllowedInPoly()) return false;
+      else if (currentState == Automaton.StateInvTFStart && !IsAllowedAsNextInPoly(followState)) return false;
+      // finishing an exponential factor is only allowed when the number of variable references is large enough
+      else if (followState == Automaton.StateExpFactorEnd && !IsCompleteExp()) return false;
+      // finishing a polynomial (in log or inv) is only allowed when the number of variable references is large enough
+      else if (followState == Automaton.StateInvTEnd && !IsTermCompleteInPoly()) return false;
+      else if (followState == Automaton.StateLogTEnd && !IsTermCompleteInPoly()) return false;
+
+      else if (nVars > maxVariables) return false;
+      else return true;
     }
 
@@ -122 +360 @@
     public void Reset() {
       nVars = 0;
-
-
-      prevTermFirstVariableState = -1;
-      curTermFirstVariableState = -1;
-      prevTermFirstFactorType = -1;
-      curTermFirstFactorType = -1;
-      prevVariableRef = -1;
-      prevFactorType = -1;
-      curFactorType = -1;
-      curFactorFirstVariableState = -1;
-      prevFactorFirstVariableState = -1;
+      prevTerm = null;
+      curTerm = null;
+      curFactor = null;
+      invalidExpression = false;
     }
 
     public void StartTerm() {
-      // reset factor type. in each term we can start with each type of factor
-      prevTermFirstVariableState = curTermFirstVariableState;
-      curTermFirstVariableState = -1;
-
-      prevTermFirstFactorType = curTermFirstFactorType;
-      curTermFirstFactorType = -1;
-
-
-      prevFactorType = -1;
-      curFactorType = -1;
-
-      curFactorFirstVariableState = -1;
-      prevFactorFirstVariableState = -1;
+      curTerm = new TermInformation();
     }
 
     public void StartFactor(int state) {
-      prevFactorType = curFactorType;
-      curFactorType = -1;
-
-      prevFactorFirstVariableState = curFactorFirstVariableState;
-      curFactorFirstVariableState = -1;
-
-
-      // store the first factor type
-      if (curTermFirstFactorType < 0) {
-        curTermFirstFactorType = state;
-      }
-      curFactorType = state;
-
-      // reset variable references. in each factor we can start with each variable reference
-      prevVariableRef = -1;
+      curFactor = new FactorInformation();
+      curFactor.factorType = state;
     }
 
 
     public void AddVarToCurrentFactor(int state) {
-
-      Contract.Assert(prevVariableRef <= state);
-
-      // store the first variable reference for each factor 
-      if (curFactorFirstVariableState < 0) {
-        curFactorFirstVariableState = state;
-
-        // store the first variable reference for each term
-        if (curTermFirstVariableState < 0) {
-          curTermFirstVariableState = state;
-        }
-      }
-      prevVariableRef = state;
+      Contract.Assert(Automaton.FirstDynamicState <= state);
+      Contract.Assert(curTerm != null);
+      Contract.Assert(curFactor != null);
 
       nVars++;
+      curFactor.numVarReferences++;
+
+      if (curFactor.factorType == Automaton.StateVariableFactorStart) {
+        Contract.Assert(curFactor.variableState < 0); // not set before
+        curFactor.variableState = state;
+      } else if (curFactor.factorType == Automaton.StateExpFactorStart) {
+        curFactor.expVariableStates.Add(state);
+      } else if (curFactor.factorType == Automaton.StateLogFactorStart ||
+                 curFactor.factorType == Automaton.StateInvFactorStart) {
+        curFactor.polyVariableStates.Last().Add(state);
+      } else throw new InvalidProgramException();
+    }
+
+    public void StartNewTermInPoly() {
+      curFactor.polyVariableStates.Add(new List<int>());
     }
 
     public void EndFactor() {
-      Contract.Assert(prevFactorFirstVariableState <= curFactorFirstVariableState);
-      Contract.Assert(prevFactorType <= curFactorType);
+      // enforce non-decreasing factors 
+      if (curTerm.factors.Any() && CompareFactors(curTerm.factors.Last(), curFactor) > 0)
+        invalidExpression = true;
+      curTerm.factors.Add(curFactor);
+      curFactor = null;
     }
 
     public void EndTerm() {
-
-      Contract.Assert(prevFactorType <= curFactorType);
-      Contract.Assert(prevTermFirstVariableState <= curTermFirstVariableState);
+      // enforce non-decreasing terms (TODO: equal terms should not be allowed)
+      if (prevTerm != null && CompareTerms(prevTerm, curTerm) > 0)
+        invalidExpression = true;
+      prevTerm = curTerm;
+      curTerm = null;
     }
   }

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/Disassembler.cs

@@ -23 +23 @@
 
 namespace HeuristicLab.Algorithms.DataAnalysis.MctsSymbolicRegression {
-#if DEBUG
   internal class Disassembler {
     public static string CodeToString(byte[] code, double[] consts) {
@@ -51 +50 @@
     }
   }
-#endif
 }

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/ExpressionEvaluator.cs

@@ -28 +28 @@
   internal class ExpressionEvaluator {
     // manages it's own vector buffers
-    private readonly List<double[]> vectorBuffers = new List<double[]>();
-    private readonly List<double[]> scalarBuffers = new List<double[]>(); // scalars are vectors of length 1 (to allow mixing scalars and vectors on the same stack)
+    private readonly double[][] vectorBuffers;
+    private readonly double[][] scalarBuffers; // scalars are vectors of length 1 (to allow mixing scalars and vectors on the same stack)
+    private int lastVecBufIdx;
+    private int lastScalarBufIdx;
 
 
     private double[] GetVectorBuffer() {
-      var v = vectorBuffers[vectorBuffers.Count - 1];
-      vectorBuffers.RemoveAt(vectorBuffers.Count - 1);
-      return v;
+      return vectorBuffers[--lastVecBufIdx];
     }
     private double[] GetScalarBuffer() {
-      var v = scalarBuffers[scalarBuffers.Count - 1];
-      scalarBuffers.RemoveAt(scalarBuffers.Count - 1);
-      return v;
+      return scalarBuffers[--lastScalarBufIdx];
     }
 
     private void ReleaseBuffer(double[] buf) {
-      (buf.Length == 1 ? scalarBuffers : vectorBuffers).Add(buf);
+      if (buf.Length == 1) {
+        scalarBuffers[lastScalarBufIdx++] = buf;
+      } else {
+        vectorBuffers[lastVecBufIdx++] = buf;
+      }
     }
 
@@ -73 +75 @@
 
       // preallocate buffers 
+      vectorBuffers = new double[MaxStackSize * (1 + MaxParams)][];
+      scalarBuffers = new double[MaxStackSize * (1 + MaxParams)][];
       for (int i = 0; i < MaxStackSize; i++) {
         ReleaseBuffer(new double[vLen]);
@@ -94 +98 @@
       short arg;
       // checked at the end to make sure we do not leak buffers
-      int initialScalarCount = scalarBuffers.Count;
-      int initialVectorCount = vectorBuffers.Count;
+      int initialScalarCount = lastScalarBufIdx;
+      int initialVectorCount = lastVecBufIdx;
 
       while (true) {
@@ -179 +183 @@
 
                 var f = 1.0 / (maxFx * consts[curParamIdx]);
-                // adjust c so that maxFx*c = 1 
+                // adjust c so that maxFx*c = 1 TODO: this is not ideal as enforce positive argument to exp()
                 consts[curParamIdx] *= f;
 
@@ -211 +215 @@
             }
             ReleaseBuffer(r);
-            Contract.Assert(vectorBuffers.Count == initialVectorCount);
-            Contract.Assert(scalarBuffers.Count == initialScalarCount);
+            Contract.Assert(lastVecBufIdx == initialVectorCount);
+            Contract.Assert(lastScalarBufIdx == initialScalarCount);
             return;
         }
@@ -232 +236 @@
 
       // checked at the end to make sure we do not leak buffers
-      int initialScalarCount = scalarBuffers.Count;
-      int initialVectorCount = vectorBuffers.Count;
+      int initialScalarCount = lastScalarBufIdx;
+      int initialVectorCount = lastVecBufIdx;
 
       while (true) {
@@ -400 +404 @@
             }
 
-            Contract.Assert(vectorBuffers.Count == initialVectorCount);
-            Contract.Assert(scalarBuffers.Count == initialScalarCount);
+            Contract.Assert(lastVecBufIdx == initialVectorCount);
+            Contract.Assert(lastScalarBufIdx == initialScalarCount);
             return; // break loop
         }

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/MctsSymbolicRegressionAlgorithm.cs

@@ -221 +221 @@
       Results.Add(new Result("Average quality", avgQuality));
 
+      var totalRollouts = new IntValue();
+      Results.Add(new Result("Total rollouts", totalRollouts));
+      var effRollouts = new IntValue();
+      Results.Add(new Result("Effective rollouts", effRollouts));
+      var funcEvals = new IntValue();
+      Results.Add(new Result("Function evaluations", funcEvals));
+      var gradEvals = new IntValue();
+      Results.Add(new Result("Gradient evaluations", gradEvals));
+
+
       // same as in SymbolicRegressionSingleObjectiveProblem
       var y = Problem.ProblemData.Dataset.GetDoubleValues(Problem.ProblemData.TargetVariable,
@@ -266 +276 @@
           curBestQ = 0.0;
 
+          funcEvals.Value = state.FuncEvaluations;
+          gradEvals.Value = state.GradEvaluations;
+          effRollouts.Value = state.EffectiveRollouts;
+          totalRollouts.Value = state.TotalRollouts;
+
           table.Rows["Best quality"].Values.Add(bestQuality.Value);
           table.Rows["Current best quality"].Values.Add(curQuality.Value);
@@ -280 +295 @@
         avgQuality.Value = sumQ / n;
 
+        funcEvals.Value = state.FuncEvaluations;
+        gradEvals.Value = state.GradEvaluations;
+        effRollouts.Value = state.EffectiveRollouts;
+        totalRollouts.Value = state.TotalRollouts;
+
         table.Rows["Best quality"].Values.Add(bestQuality.Value);
         table.Rows["Current best quality"].Values.Add(curQuality.Value);
         table.Rows["Average quality"].Values.Add(avgQuality.Value);
         iterations.Value = iterations.Value + n;
+
       }
 
@@ -289 +310 @@
       Results.Add(new Result("Best solution quality (train)", new DoubleValue(state.BestSolutionTrainingQuality)));
       Results.Add(new Result("Best solution quality (test)", new DoubleValue(state.BestSolutionTestQuality)));
+
 
       // produce solution 

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/MctsSymbolicRegressionStatic.cs

@@ -44 +44 @@
       double BestSolutionTrainingQuality { get; }
       double BestSolutionTestQuality { get; }
+      int TotalRollouts { get; }
+      int EffectiveRollouts { get; }
+      int FuncEvaluations { get; }
+      int GradEvaluations { get; } // number of gradient evaluations (* num parameters) to get a value representative of the effort comparable to the number of function evaluations
+      // TODO other stats on LM optimizer might be interesting here
     }
 
@@ -57 +62 @@
       internal readonly List<Tree> bestChildrenBuf;
       internal readonly Func<byte[], int, double> evalFun;
+      // MCTS might get stuck. Track statistics on the number of effective rollouts
+      internal int totalRollouts;
+      internal int effectiveRollouts;
 
 
@@ -77 +85 @@
       private int bestNParams;
       private double[] bestConsts;
+
+      // stats
+      private int funcEvaluations;
+      private int gradEvaluations;
 
       // buffers
@@ -173 +185 @@
         }
       }
+
+      public int TotalRollouts { get { return totalRollouts; } }
+      public int EffectiveRollouts { get { return effectiveRollouts; } }
+      public int FuncEvaluations { get { return funcEvaluations; } }
+      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
+
       #endregion
 
@@ -204 +222 @@
         Array.Copy(ones, constsBuf, nParams);
         evaluator.Exec(code, x, constsBuf, predBuf, adjustOffsetForLogAndExp: true);
+        funcEvaluations++;
 
         // calc opt scaling (alpha*f(x) + beta)
@@ -220 +239 @@
           // 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;
@@ -237 +259 @@
 
       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
+        double[] optConsts = new double[nParams]; // allocate a smaller buffer for constants opt (TODO perf?)
         Array.Copy(consts, optConsts, nParams);
 
@@ -247 +269 @@
         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);
@@ -311 +335 @@
       var rand = mctsState.random;
       double c = mctsState.c;
-
-      automaton.Reset();
-      return TreeSearchRec(rand, tree, c, automaton, eval, bestChildrenBuf);
-    }
-
-    private static double TreeSearchRec(IRandom rand, Tree tree, double c, Automaton automaton, Func<byte[], int, double> eval, List<Tree> bestChildrenBuf) {
+      double q = 0;
+      bool success = false;
+      do {
+        automaton.Reset();
+        success = TryTreeSearchRec(rand, tree, c, automaton, eval, bestChildrenBuf, out q);
+        mctsState.totalRollouts++;
+      } while (!success && !tree.done);
+      mctsState.effectiveRollouts++;
+      return q;
+    }
+
+    // tree search might fail because of constraints for expressions
+    // in this case we get stuck we just restart
+    // see ConstraintHandler.cs for more info
+    private static bool TryTreeSearchRec(IRandom rand, Tree tree, double c, Automaton automaton, Func<byte[], int, double> eval, List<Tree> bestChildrenBuf,
+      out double q) {
       Tree selectedChild = null;
-      double q;
       Contract.Assert(tree.state == automaton.CurrentState);
       Contract.Assert(!tree.done);
@@ -332 +365 @@
           tree.visits++;
           tree.sumQuality += q;
-          return q;
+          return true; // we reached a final state
         } else {
           // EXPAND 
@@ -338 +371 @@
           int nFs;
           automaton.FollowStates(automaton.CurrentState, out possibleFollowStates, out nFs);
-
+          if (nFs == 0) {
+            // stuck in a dead end (no final state and no allowed follow states)
+            q = 0;
+            tree.done = true;
+            tree.children = null;
+            return false;
+          }
           tree.children = new Tree[nFs];
-          for (int i = 0; i < tree.children.Length; i++) tree.children[i] = new Tree() { children = null, done = false, state = possibleFollowStates[i], visits = 0 };
+          for (int i = 0; i < tree.children.Length; i++)
+            tree.children[i] = new Tree() { children = null, done = false, state = possibleFollowStates[i], visits = 0 };
 
           selectedChild = SelectFinalOrRandom(automaton, tree, rand);
@@ -351 +391 @@
       // make selected step and recurse
       automaton.Goto(selectedChild.state);
-      q = TreeSearchRec(rand, selectedChild, c, automaton, eval, bestChildrenBuf);
-
-      tree.sumQuality += q;
-      tree.visits++;
+      var success = TryTreeSearchRec(rand, selectedChild, c, automaton, eval, bestChildrenBuf, out q);
+      if (success) {
+        // only update if successful
+        tree.sumQuality += q;
+        tree.visits++;
+      }
 
       // tree.done = tree.children.All(ch => ch.done);
       tree.done = true; for (int i = 0; i < tree.children.Length && tree.done; i++) tree.done = tree.children[i].done;
       if (tree.done) {
-        tree.children = null; // cut of the sub-branch if it has been fully explored
+        tree.children = null; // cut off the sub-branch if it has been fully explored
         // TODO: update all qualities and visits to remove the information gained from this whole branch
       }
-      return q;
+      return success;
     }