Changeset 17303 for branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic

Timestamp:

10/03/19 12:30:19 (5 years ago)

Author:

gkronber

Message:

#2994 continued refactoring and extended unit tests. Interval calculation still fails for some edge cases (mainly for undefined behaviour). VectorEvaluator and VectorAutoDiffEvaluator produce the same results as the LinearInterpreter. TODO: check gradient calculation

Location:

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

Files:

• HeuristicLab.Problems.DataAnalysis.Symbolic-3.4.csproj (modified) (1 diff)
• Interpreter/Algebraic.cs (modified) (1 diff)
• Interpreter/AlgebraicDouble.cs (modified) (3 diffs)
• Interpreter/AlgebraicDoubleVector.cs (modified) (2 diffs)
• Interpreter/AlgebraicInterval.cs (modified) (6 diffs)
• Interpreter/AlgebraicSparseVector.cs (modified) (1 diff)
• Interpreter/Dual.cs (modified) (2 diffs)
• Interpreter/IAlgebraicType.cs (modified) (2 diffs)
• Interpreter/MultivariateDual.cs (modified) (2 diffs)
• Interpreter/VectorAutoDiffEvaluator.cs (modified) (5 diffs)
• Interpreter/VectorOfAlgebraic.cs (added)

branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/HeuristicLab.Problems.DataAnalysis.Symbolic-3.4.csproj

@@ -182 +182 @@
     <Compile Include="Interpreter\VectorAutoDiffEvaluator.cs" />
     <Compile Include="Interpreter\VectorEvaluator.cs" />
+    <Compile Include="Interpreter\VectorOfAlgebraic.cs" />
     <Compile Include="Selectors\DiversitySelector.cs" />
     <Compile Include="SymbolicDataAnalysisExpressionTreeAverageSimilarityCalculator.cs" />

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

@@ -12 +12 @@
     public static T IntRoot<T>(this T a, int r) where T : IAlgebraicType<T> { a.AssignIntRoot(a.Clone(), r); return a; }
 
-    internal static T Min<T>(T a, T b) where T : IAlgebraicType<T> { return a.Clone().AssignMin(b); }
-    internal static T Max<T>(T a, T b) where T : IAlgebraicType<T> { return a.Clone().AssignMax(b); }
-
     // public static T Max<T>(T a, T b) where T : IAlgebraicType<T> {
     //   // ((a + b) + abs(b - a)) / 2

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

@@ -5 +5 @@
   // algebraic type wrapper for a double value
   [DebuggerDisplay("{Value}")]
-  public sealed class AlgebraicDouble : IAlgebraicType<AlgebraicDouble> {
+  public sealed class AlgebraicDouble : IComparableAlgebraicType<AlgebraicDouble> {
     public static implicit operator AlgebraicDouble(double value) { return new AlgebraicDouble(value); }
     public static implicit operator double(AlgebraicDouble value) { return value.Value; }
@@ -35 +35 @@
     public AlgebraicDouble AssignIntPower(AlgebraicDouble a, int p) { Value = Math.Pow(a.Value, p); return this; }
     public AlgebraicDouble AssignIntRoot(AlgebraicDouble a, int r) { Value = IntRoot(a.Value, r); return this; }
-    public AlgebraicDouble AssignMin(AlgebraicDouble other) { Value = Math.Min(Value, other.Value); return this; }
-    public AlgebraicDouble AssignMax(AlgebraicDouble other) { Value = Math.Max(Value, other.Value); return this; }
+    //public AlgebraicDouble AssignMin(AlgebraicDouble other) { Value = Math.Min(Value, other.Value); return this; }
+    //public AlgebraicDouble AssignMax(AlgebraicDouble other) { Value = Math.Max(Value, other.Value); return this; }
     public AlgebraicDouble AssignAbs(AlgebraicDouble a) { Value = Math.Abs(a.Value); return this; }
     public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = double.IsNaN(a.Value) ? double.NaN : Math.Sign(a.Value); return this; }
@@ -51 +51 @@
       return Value.ToString();
     }
+
+    public int CompareTo(AlgebraicDouble other) {
+      return Value.CompareTo(other.Value);
+    }
   }
 }

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

@@ -6 +6 @@
   // a simple vector as an algebraic type
   [DebuggerDisplay("DoubleVector(len={Length}): {string.}")]
+  [Obsolete("Use VectorOfAlgebraic instead")]
   public class AlgebraicDoubleVector : IAlgebraicType<AlgebraicDoubleVector> {
     private double[] arr;
@@ -40 +41 @@
     public AlgebraicDoubleVector AssignIntPower(AlgebraicDoubleVector a, int p) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Pow(a.arr[i], p); } return this; }
     public AlgebraicDoubleVector AssignIntRoot(AlgebraicDoubleVector a, int r) { for (int i = 0; i < arr.Length; ++i) { arr[i] = IntRoot(a.arr[i], r); } return this; }
-    public AlgebraicDoubleVector AssignMin(AlgebraicDoubleVector other) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Min(arr[i], other.arr[i]); } return this; }
-    public AlgebraicDoubleVector AssignMax(AlgebraicDoubleVector other) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Max(arr[i], other.arr[i]); } return this; }
+    //public AlgebraicDoubleVector AssignMin(AlgebraicDoubleVector other) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Min(arr[i], other.arr[i]); } return this; }
+    //public AlgebraicDoubleVector AssignMax(AlgebraicDoubleVector other) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Max(arr[i], other.arr[i]); } return this; }
     public AlgebraicDoubleVector AssignAbs(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Abs(a.arr[i]); } return this; }
     public AlgebraicDoubleVector AssignSgn(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sign(a.arr[i]); } return this; }

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

@@ -47 +47 @@
       var v4 = high.Clone().Mul(a.high);
 
-      low = Min(Min(v1, v2), Min(v3, v4));
-      high = Max(Max(v1, v2), Max(v3, v4));
-
-      return this;
-    }
-
-
-    private static MultivariateDual<AlgebraicDouble> Min(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
-      return a.Value < b.Value ? a : b;
-    }
-    private static MultivariateDual<AlgebraicDouble> Max(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
-      return a.Value > b.Value ? a : b;
+      AssignLowAndHigh(v1, v2, v3, v4);
+
+      return this;
     }
 
@@ -133 +124 @@
         if (a.Contains(0.0)) {
           low = new MultivariateDual<AlgebraicDouble>(0.0);
-          high = a.low.IntPower(p).AssignMax(a.high.IntPower(p));
-        } else {
-          var lowPower = a.low.IntPower(p);
-          var highPower = a.high.IntPower(p);
-          low = lowPower.AssignMin(highPower);
-          high = lowPower.AssignMax(highPower);
+          AssignMax(high, a.low.IntPower(p), a.high.IntPower(p));
+        } else {
+          AssignLowAndHigh(a.low.IntPower(p), a.high.IntPower(p));
         }
       } else {
@@ -148 +136 @@
           var lowPower = a.low.IntPower(p);
           var highPower = a.high.IntPower(p);
-          low = lowPower.AssignMin(highPower);
-          high = lowPower.AssignMax(highPower);
+          AssignMin(low, lowPower, highPower);
+          AssignMax(high, lowPower, highPower);
         }
       }
@@ -217 +205 @@
 
       // assume low and high are in the same quadrant
-      low = Algebraic.Min(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin());
-      high = Algebraic.Max(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin());
+      AssignLowAndHigh(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin()); // AssignLowAndHigh determines the lower and higher value 
 
       // override min and max if necessary
@@ -260 +247 @@
     public AlgebraicInterval AssignAbs(AlgebraicInterval a) {
       if (a.Contains(0.0)) {
+        low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound
+        AssignMax(high, a.low.Clone().Abs(), a.high.Clone().Abs());
+      } else {
         var abslow = a.low.Clone().Abs();
         var abshigh = a.high.Clone().Abs();
-        a.high.Assign(Algebraic.Max(abslow, abshigh));
-        a.low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound
-      } else {
-        var abslow = a.low.Clone().Abs();
-        var abshigh = a.high.Clone().Abs();
-        a.low.Assign(Algebraic.Min(abslow, abshigh));
-        a.high.Assign(Algebraic.Max(abslow, abshigh));
+        AssignLowAndHigh(abslow, abshigh);
       }
       return this;
@@ -279 +263 @@
     }
 
-    public AlgebraicInterval AssignMin(AlgebraicInterval other) {
-      low.AssignMin(other.low);
-      high.AssignMin(other.high);
-      return this;
-    }
-
-    public AlgebraicInterval AssignMax(AlgebraicInterval other) {
-      low.AssignMax(other.low);
-      high.AssignMax(other.high);
-      return this;
-    }
+
+    #region helper
+    private void AssignMin(MultivariateDual<AlgebraicDouble> dest, MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
+      if (a.Value.CompareTo(b.Value) <= 0) {
+        dest.Assign(a);
+      } else {
+        dest.Assign(b);
+      }
+    }
+    private void AssignMax(MultivariateDual<AlgebraicDouble> dest, MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
+      if (a.Value.CompareTo(b.Value) <= 0) {
+        dest.Assign(b);
+      } else {
+        dest.Assign(a);
+      }
+    }
+
+    // determines the smaller and larger value and sets low and high accordingly
+    private void AssignLowAndHigh(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
+      // we must make sure that low and high are different objects when a == b
+      if (a.Value.CompareTo(b.Value) == 0) {
+        low.Assign(a);
+        high.Assign(b);
+      } else {
+        AssignMin(low, a, b);
+        AssignMax(high, a, b);
+      }
+    }
+    private void AssignLowAndHigh(params MultivariateDual<AlgebraicDouble>[] xs) {
+      if (xs.Length <= 2) throw new ArgumentException("need at least 3 arguments");
+      AssignLowAndHigh(xs[0], xs[1]);
+      for(int i=2;i<xs.Length;i++) {
+        // we must make sure that low and high are different objects when a == b
+        if (low.Value.CompareTo(xs[i].Value) > 0) low.Assign(xs[i]);
+        if (high.Value.CompareTo(xs[i].Value) < 0) high.Assign(xs[i]);
+      }
+    }
+    #endregion
+
   }
 }

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

@@ -109 +109 @@
     }
 
-    public AlgebraicSparseVector<K, T> AssignMin(AlgebraicSparseVector<K, T> other) {
-      // assumes that keys without a matching key in other are zero and vice versa
-      foreach (var kvp in elems) if (!other.elems.ContainsKey(kvp.Key)) kvp.Value.AssignMin(kvp.Value.Zero); // min(v, 0)
-      foreach (var kvp in other.elems) {
-        if (elems.TryGetValue(kvp.Key, out T value))
-          value.AssignMin(kvp.Value);
-        else
-          elems.Add(kvp.Key, kvp.Value.Zero.AssignMin(kvp.Value));
-      }
-      return this;
-    }
-
-    public AlgebraicSparseVector<K, T> AssignMax(AlgebraicSparseVector<K, T> other) {
-      // assumes that keys without a matching key in other are zero and vice versa
-      foreach (var kvp in elems) if (!other.elems.ContainsKey(kvp.Key)) kvp.Value.AssignMax(kvp.Value.Zero); // max(v, 0)
-      foreach (var kvp in other.elems) {
-        if (elems.TryGetValue(kvp.Key, out T value))
-          value.AssignMax(kvp.Value);
-        else
-          elems.Add(kvp.Key, kvp.Value.Zero.AssignMax(kvp.Value));
-      }
-      return this;
-    }
-
-
     public AlgebraicSparseVector<K, T> Clone() {
       return new AlgebraicSparseVector<K, T>(this);

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

@@ -4 +4 @@
 namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
   public class Dual<V> : IAlgebraicType<Dual<V>>
-    where V : IAlgebraicType<V> {
+    where V : IComparableAlgebraicType<V> {
     [DebuggerBrowsable(DebuggerBrowsableState.Never)]
     private V v;
@@ -54 +54 @@
 
     public Dual<V> AssignMin(Dual<V> other) {
-      throw new NotImplementedException();
+      if(v.CompareTo(other.v) > 0) {
+        v.Assign(other.v);
+        dv.Assign(other.dv);
+      }
+      return this;
     }
 
     public Dual<V> AssignMax(Dual<V> other) {
-      throw new NotImplementedException();
+      if (v.CompareTo(other.v) <= 0) {
+        v.Assign(other.v);
+        dv.Assign(other.dv);
+      }
+      return this;
     }
   }

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

@@ -1 @@
-namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
+using System;
+
+namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
   public interface IAlgebraicType<T> {
     T Zero { get; } // Zero and One must create new objects
@@ -21 +23 @@
     T AssignIntRoot(T a, int r);
     T AssignSgn(T a); // set this to sign(a)
-    T AssignMin(T other); // set this min(this, other)
-    T AssignMax(T other); // set this max(this, other)
+    //T AssignMin(T other); // set this min(this, other)
+    //T AssignMax(T other); // set this max(this, other)
     T Clone();
   }
+
+  public interface IComparableAlgebraicType<T> : IAlgebraicType<T>, IComparable<T> {
+
+  }
 }

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

@@ -17 +17 @@
     public AlgebraicSparseVector<object, V> Gradient => dv; // <key,value> partial derivative identified via the key
 
+
     private MultivariateDual(MultivariateDual<V> orig) { this.v = orig.v.Clone(); this.dv = orig.dv.Clone(); }
+
+    /// <summary>
+    /// Constructor which sets value and derivative to zero
+    /// </summary>
+    public MultivariateDual() {
+      v = new V(); // assumed to be zero
+      dv = new AlgebraicSparseVector<object, V>();
+    }
 
     /// <summary>
@@ -83 +92 @@
     public MultivariateDual<V> AssignAbs(MultivariateDual<V> a) { v.AssignAbs(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Sgn()); return this; }      // abs(f(x))' = f(x)*f'(x) / |f(x)|  doesn't work for intervals      
     public MultivariateDual<V> AssignSgn(MultivariateDual<V> a) { v.AssignSgn(a.v); dv = a.dv.Zero; return this; } // sign(f(x))' = 0;     
-
-    public MultivariateDual<V> AssignMin(MultivariateDual<V> other) {
-      throw new NotImplementedException();
-    }
-
-    public MultivariateDual<V> AssignMax(MultivariateDual<V> other) {
-      throw new NotImplementedException();
-    }
   }
 }

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

@@ -5 +5 @@
 
 namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
-  public sealed class VectorAutoDiffEvaluator : InterpreterBase<MultivariateDual<AlgebraicDoubleVector>> {
+  public sealed class VectorAutoDiffEvaluator : InterpreterBase<VectorOfAlgebraic<MultivariateDual<AlgebraicDouble>>> {
     private const int BATCHSIZE = 128;
     [ThreadStatic]
@@ -57 +57 @@
       for (rowIndex = 0; rowIndex < roundedTotal; rowIndex += BATCHSIZE) {
         Evaluate(code);
-        code[0].value.Value.CopyTo(fi, rowIndex, BATCHSIZE);
 
-        // TRANSPOSE into JAC
-        var g = code[0].value.Gradient;
-        for (int j = 0; j < nParams; ++j) {
-          if (g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {
-            v.CopyColumnTo(jac, j, rowIndex, BATCHSIZE);
-          } else {
-            for (int r = 0; r < BATCHSIZE; r++) jac[rowIndex + r, j] = 0.0;
+        // code[0].value.Value.CopyTo(fi, rowIndex, BATCHSIZE);
+        var v = code[0].value;
+        for (int k = 0; k < BATCHSIZE; k++) {
+          fi[rowIndex + k] = v[k].Value.Value;
+
+          // copy gradient to Jacobian
+          var g = v[k].Gradient;
+          for (int j = 0; j < nParams; ++j) {
+            if (g.Elements.TryGetValue(j, out AlgebraicDouble gj)) {
+              jac[rowIndex + k, j] = gj.Value;
+            } else {
+              jac[rowIndex + k, j] = 0.0;
+            }
           }
         }
@@ -72 +77 @@
       if (remainingRows > 0) {
         Evaluate(code);
-        code[0].value.Value.CopyTo(fi, roundedTotal, remainingRows);
+        // code[0].value.Value.CopyTo(fi, roundedTotal, remainingRows);
+        var v = code[0].value;
+        for (int k = 0; k < remainingRows; k++) {
+          fi[roundedTotal + k] = v[k].Value.Value;
 
-        var g = code[0].value.Gradient;
-        for (int j = 0; j < nParams; ++j)
-          if (g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {
-            v.CopyColumnTo(jac, j, roundedTotal, remainingRows);
-          } else {
-            for (int r = 0; r < remainingRows; r++) jac[roundedTotal + r, j] = 0.0;
+          var g = v[k].Gradient;
+          for (int j = 0; j < nParams; ++j) {
+            if (g.Elements.TryGetValue(j, out AlgebraicDouble gj)) {
+              jac[roundedTotal + k, j] = gj.Value;
+            } else {
+              jac[roundedTotal + k, j] = 0.0;
+            }
           }
+        }
       }
     }
 
     protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
-      var zero = new AlgebraicDoubleVector(BATCHSIZE);
-      instruction.value = new MultivariateDual<AlgebraicDoubleVector>(zero);
+      instruction.value = new VectorOfAlgebraic<MultivariateDual<AlgebraicDouble>>(BATCHSIZE).Zero; // XXX zero needed?
     }
 
     protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
-      var g_arr = new double[BATCHSIZE];
       if (node2paramIdx.TryGetValue(constant, out var paramIdx)) {
-        for (int i = 0; i < BATCHSIZE; i++) g_arr[i] = 1.0;
-        var g = new AlgebraicDoubleVector(g_arr);
-        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE), paramIdx, g); // only a single column for the gradient
+        instruction.value = new VectorOfAlgebraic<MultivariateDual<AlgebraicDouble>>(BATCHSIZE);
+        for (int k = 0; k < BATCHSIZE; k++) {
+          instruction.value[k] = new MultivariateDual<AlgebraicDouble>(constant.Value, paramIdx, 1.0); // gradient is 1.0 for all elements
+        }
       } else {
-        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE));
+        instruction.value = new VectorOfAlgebraic<MultivariateDual<AlgebraicDouble>>(BATCHSIZE);
+        for (int k = 0; k < BATCHSIZE; k++) {
+          instruction.value[k] = new MultivariateDual<AlgebraicDouble>(constant.Value); // zero gradient
+        }
       }
 
-      instruction.dblVal = constant.Value;
-      instruction.value.Value.AssignConstant(instruction.dblVal);
+      instruction.dblVal = constant.Value; // also store the parameter value in the instruction (not absolutely necessary, will not be used)
     }
 
@@ -115 +126 @@
       if (node2paramIdx.ContainsKey(variable)) {
         paramIdx = node2paramIdx[variable];
-        var f = new AlgebraicDoubleVector(BATCHSIZE);
-        var g = new AlgebraicDoubleVector(BATCHSIZE);
-        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(f, paramIdx, g);
+        instruction.value = new VectorOfAlgebraic<MultivariateDual<AlgebraicDouble>>(BATCHSIZE);
+        for(int k=0;k<BATCHSIZE;k++) {
+          instruction.value[k] = new MultivariateDual<AlgebraicDouble>(0.0, paramIdx, 0.0); // values are set in LoadVariable()
+        }
       } else {
         var f = new AlgebraicDoubleVector(BATCHSIZE);
-        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(f);
+        instruction.value = new VectorOfAlgebraic<MultivariateDual<AlgebraicDouble>>(BATCHSIZE);
+        for (int k = 0; k < BATCHSIZE; k++) {
+          instruction.value[k] = new MultivariateDual<AlgebraicDouble>(0.0); // values are set in LoadVariable()
+        }
       }
 
-      instruction.dblVal = variable.Weight;
+      instruction.dblVal = variable.Weight; 
       instruction.data = new object[] { data, paramIdx };
     }
@@ -131 +146 @@
       var data = (double[])((object[])a.data)[0];
 
-      for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) a.value.Value[i - rowIndex] = data[rows[i]];
-      a.value.Scale(a.dblVal);
+      for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) {
+        a.value[i - rowIndex].Value.Assign(a.dblVal * data[rows[i]]);
+      }
 
       if (paramIdx >= 0) {
         // update gradient with variable values
-        var g = a.value.Gradient.Elements[paramIdx];
         for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) {
-          g[i - rowIndex] = data[rows[i]];
+          a.value[i - rowIndex].Gradient.Elements[paramIdx].Assign(data[rows[i]]);
         }
       }