Changeset 17200

Timestamp:

08/12/19 10:50:02 (5 years ago)

Author:

gkronber

Message:

#2994 implemented checking of autodiff with numeric differentials and fixed a bug in IntervalEvaluator.

Location:

branches/2994-AutoDiffForIntervals

Files:

• HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/ConstrainedNLSInternal.cs (modified) (5 diffs)
• HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs (modified) (4 diffs)

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

@@ -1 +1 @@
 using System;
 using System.Collections.Generic;
+using System.Diagnostics;
 using System.Linq;
 using System.Runtime.InteropServices;
@@ -193 +194 @@
       bestSolution = x;
       bestError = minf;
-      
+
       if (res < 0) {
         throw new InvalidOperationException($"NLOpt failed {res} {NLOpt.nlopt_get_errmsg(nlopt)}");
@@ -232 +233 @@
         // average
         for (int j = 0; j < grad.Length; j++) { grad[j] /= target.Length; }
+
+        #region check gradient
+        var eval = new VectorEvaluator();
+        for (int i = 0; i < dim; i++) {
+          // make two additional evaluations
+          var xForNumDiff = (double[])curX.Clone();
+          double delta = Math.Abs(xForNumDiff[i] * 1e-5);
+          xForNumDiff[i] += delta;
+          UpdateThetaValues(xForNumDiff);
+          var evalHigh = eval.Evaluate(preparedTree, problemData.Dataset, trainingRows);
+          var mseHigh = MSE(target, evalHigh);
+          xForNumDiff[i] = curX[i] - delta;
+          UpdateThetaValues(xForNumDiff);
+          var evalLow = eval.Evaluate(preparedTree, problemData.Dataset, trainingRows);
+          var mseLow = MSE(target, evalLow);
+
+          var numericDiff = (mseHigh - mseLow) / (2 * delta);
+          var autoDiff = grad[i];
+          if ((Math.Abs(autoDiff) < 1e-10 && Math.Abs(numericDiff) > 1e-2)
+            || (Math.Abs(autoDiff) >= 1e-10 && Math.Abs((numericDiff - autoDiff) / numericDiff) > 1e-2))
+            throw new InvalidProgramException();
+        }
+        #endregion
       } else {
         var eval = new VectorEvaluator();
         var prediction = eval.Evaluate(preparedTree, problemData.Dataset, trainingRows);
 
-        // calc sum of squared errors and gradient
+        // calc sum of squared errors
         sse = 0.0;
         for (int i = 0; i < target.Length; i++) {
@@ -249 +273 @@
       if (double.IsNaN(sse)) return double.MaxValue;
       return sse / target.Length;
+    }
+
+    private double MSE(double[] a, double[] b) {
+      Trace.Assert(a.Length == b.Length);
+      var sse = 0.0;
+      for (int i = 0; i < a.Length; i++) sse += (a[i] - b[i]) * (a[i] - b[i]);
+      return sse / a.Length;
     }
 
@@ -280 +311 @@
       // we transformed this to a constraint c(x) <= 0, so only the upper bound is relevant for us
       if (grad != null) for (int j = 0; j < grad.Length; j++) { grad[j] = upperGradient[j]; }
+
+      #region check gradient
+      for (int i = 0; i < dim; i++) {
+        // make two additional evaluations
+        var xForNumDiff = (double[])curX.Clone();
+        double delta = Math.Abs(xForNumDiff[i] * 1e-5);
+        xForNumDiff[i] += delta;
+        UpdateThetaValues(xForNumDiff);
+        var evalHigh = intervalEvaluator.Evaluate(constraintData.Tree, dataIntervals, constraintData.ParameterNodes,
+          out double[] unusedLowerGradientHigh, out double[] unusedUpperGradientHigh);
+
+        xForNumDiff[i] = curX[i] - delta;
+        UpdateThetaValues(xForNumDiff);
+        var evalLow = intervalEvaluator.Evaluate(constraintData.Tree, dataIntervals, constraintData.ParameterNodes,
+          out double[] unusedLowerGradientLow, out double[] unusedUpperGradientLow);
+        
+        var numericDiff = (evalHigh.UpperBound - evalLow.UpperBound) / (2 * delta);
+        var autoDiff = grad[i];
+        // XXX GRADIENT FOR UPPER BOUND FOR FUNCTION IS ZERO FOR THE FIRST SET OF VARIABLES? GRADIENT FOR ADDITION INCORRECT?!
+        if ((Math.Abs(autoDiff) < 1e-10 && Math.Abs(numericDiff) > 1e-2)
+          || (Math.Abs(autoDiff) >= 1e-10 && Math.Abs((numericDiff - autoDiff) / numericDiff) > 1e-2))
+          throw new InvalidProgramException();
+      }
+      #endregion
+
+
       UpdateConstraintViolations(constraintData.Idx, interval.UpperBound);
       if (double.IsNaN(interval.UpperBound)) return double.MaxValue;

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

@@ -788 +788 @@
     [DebuggerBrowsable(DebuggerBrowsableState.Never)]
     public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0);
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
     public AlgebraicInterval One => new AlgebraicInterval(1.0, 1.0);
+
     public AlgebraicInterval Add(AlgebraicInterval a) {
       low.Add(a.low);
@@ -818 +820 @@
     public AlgebraicInterval Div(AlgebraicInterval a) {
       if (a.Contains(0.0)) {
-        if (a.low.Value.Value.IsAlmost(0.0) && a.high.Value.Value.IsAlmost(0.0)) {
+        if (a.low.Value.Value == 0 && a.high.Value.Value == 0) {
           low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
           high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
-        } else if (a.low.Value.Value.IsAlmost(0.0))
+        } else if (a.low.Value.Value == 0)
           Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));
         else
@@ -945 +947 @@
       //divide the interval by PI/2 so that the optima lie at x element of N (0,1,2,3,4,...)
       double Pihalf = Math.PI / 2;
-      var scaled = this.Clone().Scale(1.0 / Pihalf);
+      var scaled = a.Clone().Scale(1.0 / Pihalf);
       //move to positive scale
       if (scaled.LowerBound.Value.Value < 0) {
@@ -965 +967 @@
       }
 
-      low = new MultivariateDual<AlgebraicDouble>(sinValues.Min());
+      low = new MultivariateDual<AlgebraicDouble>(sinValues.Min()); // TODO, XXX FIX CALCULATION TO SUPPORT GRADIENTS!
       high = new MultivariateDual<AlgebraicDouble>(sinValues.Max());
       return this;