Changeset 17176 for branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/ConstrainedConstantOptimizationEvaluator.cs

Timestamp:

07/28/19 19:48:12 (5 years ago)

Author:

gkronber

Message:

#2994 linear scaling for const opt with constraints

File:

• branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/ConstrainedConstantOptimizationEvaluator.cs (modified) (7 diffs)

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

@@ -223 +223 @@
       if (!updateVariableWeights) throw new NotSupportedException("not updating variable weights is not supported");
       if (!updateConstantsInTree) throw new NotSupportedException("not updating tree parameters is not supported");
-      if (applyLinearScaling) throw new NotSupportedException("linear scaling is not supported");
+      if (!applyLinearScaling) throw new NotSupportedException("application without linear scaling is not supported");
 
       // we always update constants, so we don't need to calculate initial quality
@@ -234 +234 @@
       var dataIntervals = problemData.VariableRanges.GetIntervals();
 
+      // buffers 
+      var target = problemData.TargetVariableTrainingValues.ToArray();
+      var targetStDev = target.StandardDeviationPop();
+      var targetVariance = targetStDev * targetStDev;
+      var targetMean = target.Average();
+      var pred = interpreter.GetSymbolicExpressionTreeValues(tree, problemData.Dataset, problemData.TrainingIndices).ToArray();
+      var predStDev = pred.StandardDeviationPop();
+      var predMean = pred.Average();
+
+      var scalingFactor = targetStDev / predStDev;
+      var offset = targetMean - predMean * scalingFactor;
+
+      ISymbolicExpressionTree scaledTree = null;
+      if (applyLinearScaling) scaledTree = CopyAndScaleTree(tree, scalingFactor, offset);
+
       // convert constants to variables named theta...
-      var treeForDerivation = ReplaceConstWithVar(tree, out List<string> thetaNames, out List<double> thetaValues); // copies the tree
+      var treeForDerivation = ReplaceConstWithVar(scaledTree, out List<string> thetaNames, out List<double> thetaValues); // copies the tree
 
       // create trees for relevant derivatives
@@ -288 +303 @@
       }
 
-      // buffers for calculate_jacobian
-      var target = problemData.TargetVariableTrainingValues.ToArray();
-      var targetVariance = target.VariancePop();
       var fi_eval = new double[target.Length];
       var jac_eval = new double[target.Length, thetaValues.Count];
@@ -337 +349 @@
           alglib.minnscreate(thetaValues.Count, thetaValues.ToArray(), out state);
           alglib.minnssetbc(state, thetaValues.Select(_ => -10000.0).ToArray(), thetaValues.Select(_ => +10000.0).ToArray());
-          alglib.minnssetcond(state, 1E-7, maxIterations);
+          alglib.minnssetcond(state, 0, maxIterations);
           var s = Enumerable.Repeat(1d, thetaValues.Count).ToArray();  // scale is set to unit scale
           alglib.minnssetscale(state, s);
@@ -352 +364 @@
 
           if (rep.terminationtype > 0) {
+            // update parameters in tree
+            var pIdx = 0;
+            // here we lose the two last parameters (for linear scaling)
+            foreach (var node in tree.IterateNodesPostfix()) {
+              if (node is ConstantTreeNode constTreeNode) {
+                constTreeNode.Value = xOpt[pIdx++];
+              } else if (node is VariableTreeNode varTreeNode) {
+                varTreeNode.Weight = xOpt[pIdx++];
+              }
+            }
+            // note: we keep the optimized constants even when the tree is worse.
+            // assert that we lose the last two parameters 
+            if (pIdx != xOpt.Length - 2) throw new InvalidProgramException();
+          }
+          if (Math.Abs(rep.nlcerr) > 0.01) return targetVariance; // constraints are violated
+        } catch (ArithmeticException) {
+          return targetVariance;
+        } catch (alglib.alglibexception) {
+          // eval MSE of original tree
+          return targetVariance;
+        }
+      } else if (solver.Contains("minnlc")) {
+        alglib.minnlcstate state;
+        alglib.minnlcreport rep;
+        alglib.optguardreport optGuardRep;
+        try {
+          alglib.minnlccreate(thetaValues.Count, thetaValues.ToArray(), out state);
+          alglib.minnlcsetalgoslp(state);        // SLP is more robust but slower
+          alglib.minnlcsetbc(state, thetaValues.Select(_ => -10000.0).ToArray(), thetaValues.Select(_ => +10000.0).ToArray());
+          alglib.minnlcsetcond(state, 0, maxIterations);
+          var s = Enumerable.Repeat(1d, thetaValues.Count).ToArray();  // scale is set to unit scale
+          alglib.minnlcsetscale(state, s);
+
+          // set non-linear constraints: 0 equality constraints, constraintTrees inequality constraints
+          alglib.minnlcsetnlc(state, 0, constraintTrees.Count);
+          alglib.minnlcoptguardsmoothness(state, 1);
+
+          alglib.minnlcoptimize(state, calculate_jacobian, null, null);
+          alglib.minnlcresults(state, out double[] xOpt, out rep);
+          alglib.minnlcoptguardresults(state, out optGuardRep);
+          if (optGuardRep.nonc0suspected) throw new InvalidProgramException("optGuardRep.nonc0suspected");
+          if (optGuardRep.nonc1suspected) {
+            alglib.minnlcoptguardnonc1test1results(state, out alglib.optguardnonc1test1report strrep, out alglib.optguardnonc1test1report lngrep);
+            throw new InvalidProgramException("optGuardRep.nonc1suspected");
+          }
+
+          // counter.FunctionEvaluations += rep.nfev; TODO
+          counter.GradientEvaluations += rep.nfev;
+
+          if (rep.terminationtype != -8) {
             // update parameters in tree
             var pIdx = 0;
@@ -362 +424 @@
             }
             // note: we keep the optimized constants even when the tree is worse.
-          }
-          if (Math.Abs(rep.nlcerr) > 0.01) return targetVariance; // constraints are violated
-        } catch (ArithmeticException) {
-          return targetVariance;
-        } catch (alglib.alglibexception) {
-          // eval MSE of original tree
-          return targetVariance;
-        }
-      } else if (solver.Contains("minnlc")) {
-        alglib.minnlcstate state;
-        alglib.minnlcreport rep;
-        alglib.optguardreport optGuardRep;
-        try {
-          alglib.minnlccreate(thetaValues.Count, thetaValues.ToArray(), out state);
-          alglib.minnlcsetalgoslp(state);        // SLP is more robust but slower
-          alglib.minnlcsetbc(state, thetaValues.Select(_ => -10000.0).ToArray(), thetaValues.Select(_ => +10000.0).ToArray());
-          alglib.minnlcsetcond(state, 1E-7, maxIterations);
-          var s = Enumerable.Repeat(1d, thetaValues.Count).ToArray();  // scale is set to unit scale
-          alglib.minnlcsetscale(state, s);
-
-          // set non-linear constraints: 0 equality constraints, constraintTrees inequality constraints
-          alglib.minnlcsetnlc(state, 0, constraintTrees.Count);
-          alglib.minnlcoptguardsmoothness(state, 1);
-
-          alglib.minnlcoptimize(state, calculate_jacobian, null, null);
-          alglib.minnlcresults(state, out double[] xOpt, out rep);
-          alglib.minnlcoptguardresults(state, out optGuardRep);
-          if (optGuardRep.nonc0suspected) throw new InvalidProgramException("optGuardRep.nonc0suspected");
-          if (optGuardRep.nonc1suspected) throw new InvalidProgramException("optGuardRep.nonc1suspected");
-
-          // counter.FunctionEvaluations += rep.nfev; TODO
-          counter.GradientEvaluations += rep.nfev;
-
-          if (rep.terminationtype != -8) {
-            // update parameters in tree
-            var pIdx = 0;
-            foreach (var node in tree.IterateNodesPostfix()) {
-              if (node is ConstantTreeNode constTreeNode) {
-                constTreeNode.Value = xOpt[pIdx++];
-              } else if (node is VariableTreeNode varTreeNode) {
-                varTreeNode.Weight = xOpt[pIdx++];
-              }
-            }
-
-            // note: we keep the optimized constants even when the tree is worse.
+            // assert that we lose the last two parameters 
+            if (pIdx != xOpt.Length - 2) throw new InvalidProgramException();
+
           }
           if (Math.Abs(rep.nlcerr) > 0.01) return targetVariance; // constraints are violated
@@ -421 +441 @@
 
       // evaluate tree with updated constants
-      var residualVariance = SymbolicRegressionSingleObjectiveMeanSquaredErrorEvaluator.Calculate(interpreter, tree, lowerEstimationLimit, upperEstimationLimit, problemData, rows, applyLinearScaling: false);
+      var residualVariance = SymbolicRegressionSingleObjectiveMeanSquaredErrorEvaluator.Calculate(interpreter, scaledTree, lowerEstimationLimit, upperEstimationLimit, problemData, rows, applyLinearScaling: false);
       return Math.Min(residualVariance, targetVariance);
+    }
+
+    private static ISymbolicExpressionTree CopyAndScaleTree(ISymbolicExpressionTree tree, double scalingFactor, double offset) {
+      var m = (ISymbolicExpressionTree)tree.Clone();
+
+      var add = MakeNode<Addition>(MakeNode<Multiplication>(m.Root.GetSubtree(0).GetSubtree(0), CreateConstant(scalingFactor)), CreateConstant(offset));
+      m.Root.GetSubtree(0).RemoveSubtree(0);
+      m.Root.GetSubtree(0).AddSubtree(add);
+      return m;
     }