Changeset 17325 for branches/2994-AutoDiffForIntervals

Timestamp:

10/09/19 11:13:11 (5 years ago)

Author:

gkronber

Message:

#2994: worked on ConstrainedNLS

Location:

branches/2994-AutoDiffForIntervals

Files:

• HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/ConstrainedNLS.cs (modified) (6 diffs)
• HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/ConstrainedNLSInternal.cs (modified) (15 diffs)
• HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/NLOptEvaluator.cs (modified) (2 diffs)
• HeuristicLab.Tests/HeuristicLab.Problems.DataAnalysis.Symbolic-3.4/IntervalEvaluatorAutoDiffTest.cs (modified) (4 diffs)

branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/ConstrainedNLS.cs

@@ -10 +10 @@
 using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
 using HeuristicLab.Analysis;
+using System.Collections.Generic;
 
 namespace HeuristicLab.Problems.DataAnalysis.Symbolic.Regression {
@@ -157 +158 @@
       Results.AddOrUpdateResult("Qualities", qualitiesTable);
 
-      var curConstraintValue = new DoubleValue(0);
-      Results.AddOrUpdateResult("Current Constraint Value", curConstraintValue);
-      var curConstraintIdx = new IntValue(0);
-      Results.AddOrUpdateResult("Current Constraint Index", curConstraintIdx);
-
-      var curConstraintRow = new DataRow("Constraint Value");
-      var constraintsTable = new DataTable("Constraints");
-
-      constraintsTable.Rows.Add(curConstraintRow);
+      var constraintRows = new List<IndexedDataRow<int>>(); // for access via index
+      var constraintsTable = new IndexedDataTable<int>("Constraints");
       Results.AddOrUpdateResult("Constraints", constraintsTable);
+      foreach (var constraint in problem.ProblemData.IntervalConstraints.Constraints.Where(c => c.Enabled)) {
+        if (constraint.Interval.LowerBound > double.NegativeInfinity) {
+          var constraintRow = new IndexedDataRow<int>("-" + constraint.Expression + " < " + (-constraint.Interval.LowerBound));
+          constraintRows.Add(constraintRow);
+          constraintsTable.Rows.Add(constraintRow);
+        }
+        if (constraint.Interval.UpperBound < double.PositiveInfinity) {
+          var constraintRow = new IndexedDataRow<int>(constraint.Expression + " < " + (constraint.Interval.UpperBound));
+          constraintRows.Add(constraintRow);
+          constraintsTable.Rows.Add(constraintRow);
+        }
+      }
+
+      var parametersTable = new IndexedDataTable<int>("Parameters");
 
       #endregion
@@ -175 +183 @@
       var formatter = new InfixExpressionFormatter();
       var constraintDescriptions = state.ConstraintDescriptions.ToArray();
-      foreach(var constraintTree in state.constraintTrees) {
+      foreach (var constraintTree in state.constraintTrees) {
         // HACK to remove parameter nodes which occurr multiple times
         var reparsedTree = parser.Parse(formatter.Format(constraintTree));
@@ -186 +194 @@
       state.ConstraintEvaluated += State_ConstraintEvaluated;
 
-      state.Optimize();
+      state.Optimize(ConstrainedNLSInternal.OptimizationMode.UpdateParametersAndKeepLinearScaling);
       bestError.Value = state.BestError;
       curQualityRow.Values.Add(state.CurError);
@@ -205 +213 @@
         curQualityRow.Values.Add(state.CurError);
         bestQualityRow.Values.Add(bestError.Value);
+
+        // on the first call create the data rows
+        if(!parametersTable.Rows.Any()) {
+          for(int i=0;i<state.BestSolution.Length;i++) {
+            parametersTable.Rows.Add(new IndexedDataRow<int>("p" + i));
+          }
+        }
+        for (int i = 0; i < state.BestSolution.Length; i++) {
+          parametersTable.Rows["p" + i].Values.Add(Tuple.Create(functionEvaluations.Value, state.BestSolution[i])); // TODO: remove access via string
+        }
       }
 
       // local function
       void State_ConstraintEvaluated(int constraintIdx, double value) {
-        curConstraintIdx.Value = constraintIdx;
-        curConstraintValue.Value = value;
-        curConstraintRow.Values.Add(value);
+        constraintRows[constraintIdx].Values.Add(Tuple.Create(functionEvaluations.Value, value));
       }
     }
@@ -217 +233 @@
     private static ISymbolicRegressionSolution CreateSolution(ISymbolicExpressionTree tree, IRegressionProblemData problemData) {
       var model = new SymbolicRegressionModel(problemData.TargetVariable, tree, new SymbolicDataAnalysisExpressionTreeLinearInterpreter());
-      // model.Scale(problemData);
+      // model.CreateRegressionSolution produces a new ProblemData and recalculates intervals ==> use SymbolicRegressionSolution.ctor instead
+      var sol = new SymbolicRegressionSolution(model, (IRegressionProblemData)problemData.Clone());
+      // NOTE: the solution has slightly different derivative values because simplification of derivatives can be done differently when parameter values are fixed.
+
       // var sol = model.CreateRegressionSolution((IRegressionProblemData)problemData.Clone());
-      // model.CreateRegressionSolution produces a new ProblemData and recalculates intervals
-
-      var sol = new SymbolicRegressionSolution(model, (IRegressionProblemData)problemData.Clone());
-
-      // ==> set variable ranges to same range as in original problemData
-      // foreach(var interval in problemData.VariableRanges.GetIntervals()) {
-      //   sol.ProblemData.VariableRanges.SetInterval
-      // }
+
       return sol;
     }

branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/ConstrainedNLSInternal.cs

@@ -6 +6 @@
 using HeuristicLab.Common;
 using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
+using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression.Extensions;
 
 namespace HeuristicLab.Problems.DataAnalysis.Symbolic.Regression {
@@ -40 +41 @@
     private double[] bestConstraintValues;
     public double[] BestConstraintValues => bestConstraintValues;
+
+    private bool disposed = false;
 
 
@@ -120 +123 @@
       }
 
+      // all trees are linearly scaled (to improve GP performance)
       #region linear scaling
       var predStDev = pred.StandardDeviationPop();
@@ -134 +138 @@
 
       // convert constants to variables named theta...
-      var treeForDerivation = ReplaceConstWithVar(scaledTree, out List<string> thetaNames, out thetaValues); // copies the tree
+      var treeForDerivation = ReplaceAndExtractParameters(scaledTree, out List<string> thetaNames, out thetaValues); // copies the tree
 
       // create trees for relevant derivatives
@@ -220 +224 @@
 
     ~ConstrainedNLSInternal() {
-      Dispose();
-    }
-
-
-    internal void Optimize() {
+      Dispose(false);
+    }
+
+
+    public enum OptimizationMode { ReadOnly, UpdateParameters, UpdateParametersAndKeepLinearScaling };
+
+    internal void Optimize(OptimizationMode mode) {
       if (invalidProblem) return;
       var x = thetaValues.ToArray();  /* initial guess */
@@ -233 +239 @@
         // throw new InvalidOperationException($"NLOpt failed {res} {NLOpt.nlopt_get_errmsg(nlopt)}");
         return;
-      } else if (minf <= bestError) {
+      } else /*if ( minf <= bestError ) */{
         bestSolution = x;
         bestError = minf;
@@ -245 +251 @@
 
         // update parameters in tree
-        var pIdx = 0;
-        // here we lose the two last parameters (for linear scaling)
-        foreach (var node in scaledTree.IterateNodesPostfix()) {
-          if (node is ConstantTreeNode constTreeNode) {
-            constTreeNode.Value = x[pIdx++];
-          } else if (node is VariableTreeNode varTreeNode) {
-            varTreeNode.Weight = x[pIdx++];
-          }
-        }
-        if (pIdx != x.Length) throw new InvalidProgramException();
-      }
-      bestTree = scaledTree;
-    }
+        UpdateParametersInTree(scaledTree, x);
+
+        if (mode == OptimizationMode.UpdateParameters) {
+          // update original expression (when called from evaluator we want to write back optimized parameters)
+          expr.Root.GetSubtree(0).RemoveSubtree(0); // delete old tree
+          expr.Root.GetSubtree(0).InsertSubtree(0,
+            scaledTree.Root.GetSubtree(0).GetSubtree(0).GetSubtree(0).GetSubtree(0) // insert the optimized sub-tree (without scaling nodes)
+            );
+        } else if (mode == OptimizationMode.UpdateParametersAndKeepLinearScaling) {
+          expr.Root.GetSubtree(0).RemoveSubtree(0); // delete old tree
+          expr.Root.GetSubtree(0).InsertSubtree(0, scaledTree.Root.GetSubtree(0).GetSubtree(0)); // insert the optimized sub-tree (including scaling nodes)
+        }
+      }
+      bestTree = expr;
+    }
+
 
     double CalculateObjective(uint dim, double[] curX, double[] grad, IntPtr data) {
@@ -314 +323 @@
       }
 
-      // UpdateBestSolution(sse / target.Length, curX);
+      UpdateBestSolution(sse / target.Length, curX);
       RaiseFunctionEvaluated();
 
@@ -427 +436 @@
       UpdateConstraintViolations(constraintData.Idx, interval.UpperBound);
       if (double.IsNaN(interval.UpperBound)) {
-        if(grad!=null)Array.Clear(grad, 0, grad.Length);
+        if (grad != null) Array.Clear(grad, 0, grad.Length);
         return double.MaxValue;
       } else return interval.UpperBound;
@@ -463 +472 @@
     }
 
-    private static void UpdateConstants(ISymbolicExpressionTreeNode[] nodes, double[] constants) {
-      if (nodes.Length != constants.Length) throw new InvalidOperationException();
-      for (int i = 0; i < nodes.Length; i++) {
-        if (nodes[i] is VariableTreeNode varNode) varNode.Weight = constants[i];
-        else if (nodes[i] is ConstantTreeNode constNode) constNode.Value = constants[i];
-      }
-    }
 
     private NLOpt.nlopt_algorithm GetSolver(string solver) {
@@ -514 +516 @@
     }
 
+    // determines the nodes over which we can calculate the partial derivative
+    // this is different from the vector of all parameters because not every tree contains all parameters
     private static ISymbolicExpressionTreeNode[] GetParameterNodes(ISymbolicExpressionTree tree, List<ConstantTreeNode>[] allNodes) {
       // TODO better solution necessary
@@ -553 +557 @@
     }
 
-    private static ISymbolicExpressionTree ReplaceConstWithVar(ISymbolicExpressionTree tree, out List<string> thetaNames, out List<double> thetaValues) {
+
+
+
+    private void UpdateParametersInTree(ISymbolicExpressionTree scaledTree, double[] x) {
+      var pIdx = 0;
+      // here we lose the two last parameters (for linear scaling)
+      foreach (var node in scaledTree.IterateNodesPostfix()) {
+        if (node is ConstantTreeNode constTreeNode) {
+          constTreeNode.Value = x[pIdx++];
+        } else if (node is VariableTreeNode varTreeNode) {
+          if (varTreeNode.Weight != 1.0) // see ReplaceAndExtractParameters
+            varTreeNode.Weight = x[pIdx++];
+        }
+      }
+      if (pIdx != x.Length) throw new InvalidProgramException();
+    }
+
+    private static ISymbolicExpressionTree ReplaceAndExtractParameters(ISymbolicExpressionTree tree, out List<string> thetaNames, out List<double> thetaValues) {
       thetaNames = new List<string>();
       thetaValues = new List<double>();
@@ -578 +599 @@
         }
         if (node is VariableTreeNode varTreeNode) {
+          if (varTreeNode.Weight == 1) continue; // NOTE: here we assume that we do not tune variable weights when they are originally exactly 1 because we assume that the tree has been parsed and the tree explicitly has the structure w * var
+
           var thetaVar = (VariableTreeNode)new Problems.DataAnalysis.Symbolic.Variable().CreateTreeNode();
           thetaVar.Weight = 1;
@@ -626 +649 @@
 
     public void Dispose() {
+      Dispose(true);
+      GC.SuppressFinalize(this);
+    }
+
+    protected virtual void Dispose(bool disposing) {
+      if (disposed)
+        return;
+
+      if (disposing) {
+        // Free any other managed objects here.
+      }
+
+      // Free any unmanaged objects here.
       if (nlopt != IntPtr.Zero) {
         NLOpt.nlopt_destroy(nlopt);
@@ -637 +673 @@
           }
       }
+
+      disposed = true;
     }
     #endregion

branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/NLOptEvaluator.cs

@@ -199 +199 @@
     }
 
-    private static void GetParameterNodes(ISymbolicExpressionTree tree, out List<ISymbolicExpressionTreeNode> thetaNodes, out List<double> thetaValues) {
-      thetaNodes = new List<ISymbolicExpressionTreeNode>();
-      thetaValues = new List<double>();
-
-      var nodes = tree.IterateNodesPrefix().ToArray();
-      for (int i = 0; i < nodes.Length; ++i) {
-        var node = nodes[i];
-        if (node is VariableTreeNode variableTreeNode) {
-          thetaValues.Add(variableTreeNode.Weight);
-          thetaNodes.Add(node);
-        } else if (node is ConstantTreeNode constantTreeNode) {
-          thetaNodes.Add(node);
-          thetaValues.Add(constantTreeNode.Value);
-        }
-      }
-    }
-
 
     public static double OptimizeConstants(ISymbolicDataAnalysisExpressionTreeInterpreter interpreter,
@@ -230 +213 @@
 
       using (var state = new ConstrainedNLSInternal(solver, tree, maxIterations, problemData, 0, 0, 0)) {
-        state.Optimize();
+        state.Optimize(ConstrainedNLSInternal.OptimizationMode.UpdateParameters);
         return state.BestError;
       }

branches/2994-AutoDiffForIntervals/HeuristicLab.Tests/HeuristicLab.Problems.DataAnalysis.Symbolic-3.4/IntervalEvaluatorAutoDiffTest.cs

@@ -61 +61 @@
       var eval = new IntervalEvaluator();
       var parser = new InfixExpressionParser();
+      var intervals = new Dictionary<string, Interval>() {
+        { "x", new Interval(1, 2) },
+        { "unit", new Interval(0, 1) },
+        { "neg", new Interval(-1, 0) },
+      };
       var t = parser.Parse("sqr(x)");
       var paramNodes = t.IterateNodesPostfix().Where(n => n.SubtreeCount == 0).ToArray();
-      var intervals = new Dictionary<string, Interval>() {
-        { "x", new Interval(1, 2) },
-        { "y", new Interval(0, 1) }
-      };
-      var r = eval.Evaluate(t, intervals, paramNodes, out double[] lg, out double[] ug);
-      // TODO
-      // Assert.AreEqual(XXX, r.LowerBound);
-      // Assert.AreEqual(XXX, r.UpperBound);
-      // 
-      // Assert.AreEqual(XXX, lg[0]); // x
-      // Assert.AreEqual(XXX, ug[0]);
-      // 
-      // for  { "x", new Interval(1, 2) },
-      //   { "y", new Interval(0, 1) },
-      // 
-      // 0 <> -2,50012500572888E-05 for y in SQR(LOG('y'))
-      // 0 <> 2, 49987500573946E-05 for x in SQR(LOG('x'))
+      var r = eval.Evaluate(t, intervals, paramNodes, out double[] lg, out double[] ug);
+      Assert.AreEqual(1, r.LowerBound);
+      Assert.AreEqual(4, r.UpperBound);
+
+      Assert.AreEqual(2.0, lg[0]); // x
+      Assert.AreEqual(8.0, ug[0]);
+
+      t = parser.Parse("sqr(log(unit))");
+      paramNodes = t.IterateNodesPostfix().Where(n => n.SubtreeCount == 0).ToArray();
+      r = eval.Evaluate(t, intervals, paramNodes, out lg, out ug);
+      Assert.AreEqual(0.0, r.LowerBound);
+      Assert.AreEqual(double.PositiveInfinity, r.UpperBound);
+
+      Assert.AreEqual(0.0, lg[0]); // x
+      Assert.AreEqual(double.NaN, ug[0]);
+
     }
 
@@ -126 +130 @@
       var eval = new IntervalEvaluator();
       var parser = new InfixExpressionParser();
+      var intervals = new Dictionary<string, Interval>() {
+        { "x", new Interval(3, 4) },
+        { "z", new Interval(1, 2) }
+      };
       var t = parser.Parse("cos(x)");
       var paramNodes = t.IterateNodesPostfix().Where(n => n.SubtreeCount == 0).ToArray();
-      var intervals = new Dictionary<string, Interval>() {
-        { "x", new Interval(3, 4) },
-      };
       var r = eval.Evaluate(t, intervals, paramNodes, out double[] lg, out double[] ug);
       Assert.AreEqual(-1, r.LowerBound); //  3..4 crosses pi and cos(pi) == -1
@@ -137 +142 @@
       Assert.AreEqual(0, lg[0]); // x
       Assert.AreEqual(-4 * Math.Sin(4), ug[0]);
+
+      t = parser.Parse("LOG(COS('z'))");
+      paramNodes = t.IterateNodesPostfix().Where(n => n.SubtreeCount == 0).ToArray();
+      r = eval.Evaluate(t, intervals, paramNodes, out  lg, out  ug);
+      Assert.AreEqual(double.NaN, r.LowerBound); 
+      Assert.AreEqual(Math.Log(Math.Cos(1)), r.UpperBound); 
+
+      Assert.AreEqual(-2 * Math.Sin(2) / Math.Cos(2), lg[0], 1e-5); // x
+      Assert.AreEqual(-1 * Math.Sin(1) / Math.Cos(1), ug[0], 1e-5);
+      
     }
 
@@ -182 +197 @@
       r = eval.Evaluate(t, intervals, paramNodes, out lg, out ug);
 
-      Assert.AreEqual(0.5 * Math.Sqrt(1e-10), lg[0], 1e-6); // z          --> lim x -> 0 (sqrt(x)) = 0
+      Assert.AreEqual(0.5 * Math.Sqrt(1e-10), lg[0], 1e-5); // --> lim x -> 0 (sqrt(x)) = 0
       Assert.AreEqual(0.5, ug[0], 1e-5);
+
+      t = parser.Parse("sqrt(y - z)"); // 1..2 - 0..1
+      paramNodes = t.IterateNodesPostfix().Where(n => n.SubtreeCount == 0).ToArray();
+      r = eval.Evaluate(t, intervals, paramNodes, out lg, out ug);
+      Assert.AreEqual(0, r.LowerBound);
+      Assert.AreEqual(Math.Sqrt(2), r.UpperBound);
+
+      Assert.AreEqual(double.PositiveInfinity, lg[0], 1e-5); // y
+      Assert.AreEqual(1/ Math.Sqrt(2)  , ug[0], 1e-5);
+      Assert.AreEqual(double.NegativeInfinity, lg[1], 1e-5); // z
+      Assert.AreEqual(0.0   , ug[1], 1e-5);
     }