Changeset 16682

Timestamp:

03/13/19 23:55:38 (6 years ago)

Author:

gkronber

Message:

#2994: worked on auto diff for intervals and vectors

Location:

branches/2994-AutoDiffForIntervals

Files:

• HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Crossovers/SymbolicDataAnalysisExpressionProbabilisticFunctionalCrossover.cs (modified) (1 diff)
• HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs (modified) (24 diffs)
• Tests/AutoDiffTest.cs (modified) (4 diffs)

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

@@ -116 +116 @@
       // remove branches with an infinite or NaN behavioral distance
       for (int i = weights.Count - 1; i >= 0; --i) {
-        if (Double.IsNaN(weights[i]) || Double.IsInfinity(weights[i])) {
+        if (double.IsNaN(weights[i]) || double.IsInfinity(weights[i])) {
           weights.RemoveAt(i);
           allowedBranches.RemoveAt(i);

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

@@ -192 +192 @@
   }
 
-  public class VectorAutoDiffEvaluator : Interpreter<VectorDual<DoubleVector>> {
+  public class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<DoubleVector>> {
     private const int BATCHSIZE = 128;
     [ThreadStatic]
@@ -241 +241 @@
 
         // TRANSPOSE into JAC
-        // here we assume that we usually calculate the gradient over _all_ parameters (otherwise it would be better to propagate only relevant vectors in evaluation)
-        for (int j = 0; j < nParams; ++j) 
-          code[0].value.Gradient[j].CopyColumnTo(jac, j, rowIndex, BATCHSIZE); XXX CONTINUE HERE
+        var g = code[0].value.Gradient;
+        for (int j = 0; j < nParams; ++j) {
+          g.Elements[j].CopyColumnTo(jac, j, rowIndex, BATCHSIZE);
+        }
       }
 
@@ -249 +250 @@
         Evaluate(code);
         code[0].value.Value.CopyTo(fi, roundedTotal, remainingRows);
+
+        var g = code[0].value.Gradient;
         for (int j = 0; j < nParams; ++j)
-          code[0].value.Gradient[j].CopyColumnTo(jac, j, roundedTotal, remainingRows);
+          g.Elements[j].CopyColumnTo(jac, j, roundedTotal, remainingRows);
       }
     }
@@ -256 +259 @@
     protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
       var zero = new DoubleVector(BATCHSIZE);
-      //       var g = new Vector<DoubleVector>(node2paramIdx.Count);
-      //       for (int i = 0; i < BATCHSIZE; ++i) g[i] = new DoubleVector(BATCHSIZE);
-      instruction.value = new VectorDual<DoubleVector>(zero, null);
+      instruction.value = new MultivariateDual<DoubleVector>(zero);
     }
 
     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 DoubleVector(g_arr);
+        instruction.value = new MultivariateDual<DoubleVector>(new DoubleVector(BATCHSIZE), paramIdx, g); // only a single column for the gradient
+      } else {
+        instruction.value = new MultivariateDual<DoubleVector>(new DoubleVector(BATCHSIZE));
+      }
+
       instruction.dblVal = constant.Value;
-
-      var g_arr = new double[BATCHSIZE];
-      if (node2paramIdx.ContainsKey(constant)) {
-        for (int i = 0; i < BATCHSIZE; i++) g_arr[i] = 1.0;
-      }
-      var g = new DoubleVector(g_arr);
-
-      instruction.value = new VectorDual<DoubleVector>(new DoubleVector(BATCHSIZE), new Vector<DoubleVector>(new[] { g })); // only a single column for the gradient
       instruction.value.Value.AssignConstant(instruction.dblVal);
     }
@@ -286 +288 @@
       if (node2paramIdx.ContainsKey(variable)) {
         paramIdx = node2paramIdx[variable];
+        var f = new DoubleVector(BATCHSIZE);
+        var g = new DoubleVector(BATCHSIZE);
+        instruction.value = new MultivariateDual<DoubleVector>(f, paramIdx, g);
+      } else {
+        var f = new DoubleVector(BATCHSIZE);
+        instruction.value = new MultivariateDual<DoubleVector>(f);
       }
 
       instruction.dblVal = variable.Weight;
       instruction.data = new object[] { data, paramIdx };
-
-      var f = new DoubleVector(BATCHSIZE);
-      var g = new DoubleVector(BATCHSIZE);
-      instruction.value = new VectorDual<DoubleVector>(f, new Vector<DoubleVector>(new[] { g }));
     }
 
@@ -305 +309 @@
       if (paramIdx >= 0) {
         // update gradient with variable values
-        var g = a.value.Gradient[0];
+        var g = a.value.Gradient.Elements[paramIdx];
         for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) {
           g[i] = data[rows[i]];
@@ -316 +320 @@
   public class IntervalEvaluator : Interpreter<AlgebraicInterval> {
     [ThreadStatic]
-    private Dictionary<string, AlgebraicInterval> intervals;
-
-    public AlgebraicInterval Evaluate(ISymbolicExpressionTree tree, Dictionary<string, AlgebraicInterval> intervals) {
+    private Dictionary<string, Interval> intervals;
+
+    public Interval Evaluate(ISymbolicExpressionTree tree, Dictionary<string, Interval> intervals) {
       this.intervals = intervals;
       var code = Compile(tree);
       Evaluate(code);
-      return code[0].value;
-    }
-
+      return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);
+    }
+
+    public Interval Evaluate(ISymbolicExpressionTree tree, Dictionary<string, Interval> intervals, ISymbolicExpressionTreeNode[] paramNodes, out double[] lowerGradient, out double[] upperGradient) {
+      this.intervals = intervals;
+      var code = Compile(tree);
+      Evaluate(code);
+      lowerGradient = new double[paramNodes.Length];
+      upperGradient = new double[paramNodes.Length];
+      var l = code[0].value.LowerBound;
+      var u = code[0].value.UpperBound;
+      for (int i = 0; i < paramNodes.Length; ++i) {
+        lowerGradient[i] = l.Gradient.Elements[paramNodes[i]];
+        upperGradient[i] = u.Gradient.Elements[paramNodes[i]];
+      }
+      return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);
+    }
 
     protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
@@ -333 +351 @@
     protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
       instruction.dblVal = constant.Value;
-      instruction.value = new AlgebraicInterval(constant.Value, constant.Value);
+      instruction.value = new AlgebraicInterval(
+        new MultivariateDual<Double>(constant.Value, constant, 1.0),
+        new MultivariateDual<Double>(constant.Value, constant, 1.0) // use node as key
+        );
     }
 
     protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
       instruction.dblVal = variable.Weight;
-      instruction.value = intervals[variable.VariableName];
+      instruction.value = new AlgebraicInterval(
+        low: new MultivariateDual<Double>(intervals[variable.VariableName].LowerBound, variable, intervals[variable.VariableName].LowerBound),  // bounds change by variable value d/dc (c I(var)) = I(var)
+        high: new MultivariateDual<Double>(intervals[variable.VariableName].UpperBound, variable, intervals[variable.VariableName].UpperBound)
+        );
     }
 
@@ -347 +371 @@
 
   public interface IAlgebraicType<T> {
+    T AssignAbs(T a); // set this to assign abs(a)
     T Assign(T a); // assign this to same value as a (copy!)
     T AssignNeg(T a); // set this to negative(a)
@@ -361 +386 @@
     T AssignIntPower(T a, int p);
     T AssignIntRoot(T a, int r);
+    T AssignSgn(T a); // set this to sign(a)
     T Clone();
   }
 
+  public static class Algebraic {
+    public static T Abs<T>(this T a) where T : IAlgebraicType<T> { a.AssignAbs(a); return a; }
+    public static T Neg<T>(this T a) where T : IAlgebraicType<T> { a.AssignNeg(a); return a; }
+    public static T Inv<T>(this T a) where T : IAlgebraicType<T> { a.AssignInv(a); return a; }
+    public static T Log<T>(this T a) where T : IAlgebraicType<T> { a.AssignLog(a); return a; }
+    public static T Exp<T>(this T a) where T : IAlgebraicType<T> { a.AssignExp(a); return a; }
+    public static T Sin<T>(this T a) where T : IAlgebraicType<T> { a.AssignSin(a); return a; }
+    public static T Cos<T>(this T a) where T : IAlgebraicType<T> { a.AssignCos(a); return a; }
+    public static T Sgn<T>(this T a) where T : IAlgebraicType<T> { a.AssignSgn(a); return a; }
+    public static T IntPower<T>(this T a, int p) where T : IAlgebraicType<T> { a.AssignIntPower(a, p); return a; }
+    public static T IntRoot<T>(this T a, int r) where T : IAlgebraicType<T> { a.AssignIntRoot(a, r); return a; }
+
+    public static T Max<T>(T a, T b) where T : IAlgebraicType<T> {
+      // ((a + b) + abs(b - a)) / 2
+      return a.Clone().Add(b).Add(b.Clone().Sub(a).Abs()).Scale(1.0 / 2.0);
+    }
+    public static T Min<T>(T a, T b) where T : IAlgebraicType<T> {
+      // ((a + b) - abs(a - b)) / 2
+      return a.Clone().Add(b).Sub(a.Clone().Sub(b).Abs()).Scale(1.0 / 2.0);
+    }
+  }
+
+
+  // algebraic type wrapper for a double value
+  public class Double : IAlgebraicType<Double> {
+    public static implicit operator Double(double value) { return new Double(value); }
+    public static implicit operator double(Double value) { return value.Value; }
+    public double Value;
+    public Double() { }
+    public Double(double value) { this.Value = value; }
+    public Double Add(Double a) { Value += a.Value; return this; }
+    public Double Assign(Double a) { Value = a.Value; return this; }
+
+    public Double AssignAbs(Double a) { Value = Math.Abs(a.Value); return this; }
+    public Double AssignCos(Double a) { Value = Math.Cos(a.Value); return this; }
+    public Double AssignExp(Double a) { Value = Math.Exp(a.Value); return this; }
+    public Double AssignIntPower(Double a, int p) { Value = Math.Pow(a.Value, p); return this; }
+    public Double AssignIntRoot(Double a, int r) { Value = Math.Pow(a.Value, 1.0 / r); return this; }
+    public Double AssignInv(Double a) { Value = 1.0 / a.Value; return this; }
+    public Double AssignLog(Double a) { Value = Math.Log(a.Value); return this; }
+    public Double AssignNeg(Double a) { Value = -a.Value; return this; }
+    public Double AssignSin(Double a) { Value = Math.Sin(a.Value); return this; }
+    public Double AssignSgn(Double a) { Value = Math.Sign(a.Value); return this; }
+    public Double Clone() { return new Double(Value); }
+    public Double Div(Double a) { Value /= a.Value; return this; }
+    public Double Mul(Double a) { Value *= a.Value; return this; }
+    public Double Scale(double s) { Value *= s; return this; }
+    public Double Sub(Double a) { Value -= a.Value; return this; }
+
+  }
 
   // a simple vector as an algebraic type
@@ -374 +450 @@
       arr = new double[length];
     }
+
+    public DoubleVector() { }
 
     /// <summary>
@@ -398 +476 @@
     }
 
-
-    public DoubleVector Add(DoubleVector a) {
-      for (int i = 0; i < arr.Length; ++i) {
-        arr[i] += a.arr[i];
-      }
-      return this;
-    }
-
     public void AssignConstant(double constantValue) {
       for (int i = 0; i < arr.Length; ++i) {
@@ -412 +482 @@
     }
 
-    public DoubleVector Assign(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
-    public DoubleVector AssignCos(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }
-    public DoubleVector Div(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }
-    public DoubleVector AssignExp(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }
+    public DoubleVector Assign(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
+    public DoubleVector AssignCos(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }
+    public DoubleVector Div(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }
+    public DoubleVector AssignExp(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }
     public DoubleVector AssignIntPower(DoubleVector a, int p) { throw new NotImplementedException(); }
     public DoubleVector AssignIntRoot(DoubleVector a, int r) { throw new NotImplementedException(); }
-    public DoubleVector AssignInv(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }
-    public DoubleVector AssignLog(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }
+    public DoubleVector AssignInv(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }
+    public DoubleVector AssignLog(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }
+    public DoubleVector Add(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; }
     public DoubleVector Mul(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; }
-    public DoubleVector AssignNeg(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = -a.arr[i]; } return this; }
+    public DoubleVector AssignNeg(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = -a.arr[i]; } return this; }
     public DoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; }
-    public DoubleVector AssignSin(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }
+    public DoubleVector AssignSin(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }
     public DoubleVector Sub(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] -= a.arr[i]; } return this; }
-
+    public DoubleVector AssignAbs(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Abs(a.arr[i]); } return this; }
+    public DoubleVector AssignSgn(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sign(a.arr[i]); } return this; }
     public DoubleVector Clone() {
       var v = new DoubleVector(this.arr.Length);
@@ -437 +509 @@
 
   }
-
 
   // vectors of algebraic types
@@ -506 +577 @@
     public Vector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }
     public Vector<T> Sub(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }
+    public Vector<T> AssignAbs(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignAbs(a.elems[i]); } return this; }
+    public Vector<T> AssignSgn(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSgn(a.elems[i]); } return this; }
+  }
+
+
+  /// <summary>
+  /// A sparse vector of algebraic types. Elements are accessed via a key of type K
+  /// </summary>
+  /// <typeparam name="K">Key type</typeparam>
+  /// <typeparam name="T">Element type</typeparam>
+  public class SparseVector<K, T> : IAlgebraicType<SparseVector<K, T>> where T : IAlgebraicType<T> {
+
+    private Dictionary<K, T> elems;
+
+    public IReadOnlyDictionary<K, T> Elements => elems;
+
+    public SparseVector(SparseVector<K, T> original) {
+      elems = original.elems.ToDictionary(kvp => kvp.Key, kvp => kvp.Value.Clone());
+    }
+
+    public SparseVector(K[] keys, T[] values) {
+      if (keys.Length != values.Length) throw new ArgumentException("lengths of keys and values doesn't match in SparseVector");
+      elems = new Dictionary<K, T>(keys.Length);
+      for (int i = 0; i < keys.Length; ++i) {
+        elems.Add(keys[i], values[i]);
+      }
+    }
+
+    public SparseVector() {
+      this.elems = new Dictionary<K, T>();
+    }
+
+    public SparseVector<K, T> Add(SparseVector<K, T> a) {
+      foreach (var kvp in a.elems) {
+        if (elems.TryGetValue(kvp.Key, out var value)) {
+          value.Add(kvp.Value);
+        } else {
+          elems.Add(kvp.Key, kvp.Value.Clone());
+        }
+      }
+      return this;
+    }
+
+    public SparseVector<K, T> Scale(T s) {
+      foreach (var kvp in elems) {
+        kvp.Value.Mul(s);
+      }
+      return this;
+    }
+
+    public SparseVector<K, T> Scale(double s) {
+      foreach (var kvp in elems) {
+        kvp.Value.Scale(s);
+      }
+      return this;
+    }
+
+
+    public SparseVector<K, T> Assign(SparseVector<K, T> a) {
+      elems.Clear();
+      elems = a.elems.ToDictionary(kvp => kvp.Key, kvp => kvp.Value.Clone());
+      return this;
+    }
+
+    public SparseVector<K, T> AssignCos(SparseVector<K, T> a) {
+      throw new NotImplementedException();
+    }
+
+    public SparseVector<K, T> AssignExp(SparseVector<K, T> a) {
+      throw new NotImplementedException();
+    }
+
+    public SparseVector<K, T> AssignIntPower(SparseVector<K, T> a, int p) {
+      throw new NotImplementedException();
+    }
+
+    public SparseVector<K, T> AssignIntRoot(SparseVector<K, T> a, int r) {
+      throw new NotImplementedException();
+    }
+
+    public SparseVector<K, T> AssignInv(SparseVector<K, T> a) {
+      throw new NotImplementedException();
+    }
+
+    public SparseVector<K, T> AssignLog(SparseVector<K, T> a) {
+      throw new NotImplementedException();
+    }
+
+    public SparseVector<K, T> AssignNeg(SparseVector<K, T> a) {
+      throw new NotImplementedException();
+    }
+
+    public SparseVector<K, T> AssignSin(SparseVector<K, T> a) {
+      throw new NotImplementedException();
+    }
+
+    public SparseVector<K, T> Clone() {
+      return new SparseVector<K, T>(this);
+    }
+
+    public SparseVector<K, T> Div(SparseVector<K, T> a) {
+      throw new NotImplementedException();
+    }
+
+    public SparseVector<K, T> Mul(SparseVector<K, T> a) {
+      throw new NotImplementedException();
+    }
+
+    public SparseVector<K, T> Sub(SparseVector<K, T> a) {
+      foreach (var kvp in a.elems) {
+        if (elems.TryGetValue(kvp.Key, out var value)) {
+          value.Sub(kvp.Value);
+        } else {
+          elems.Add(kvp.Key, kvp.Value.Clone().Neg());
+        }
+      }
+      return this;
+    }
+
+    public SparseVector<K, T> AssignAbs(SparseVector<K, T> a) {
+      elems.Clear();
+      foreach (var kvp in a.elems) {
+        elems.Add(kvp.Key, kvp.Value.Clone().Abs());
+      }
+      return this;
+    }
+
+    public SparseVector<K, T> AssignSgn(SparseVector<K, T> a) {
+      elems.Clear();
+      foreach (var kvp in a.elems) {
+        elems.Add(kvp.Key, kvp.Value.Clone().Sgn());
+      }
+      return this;
+    }
   }
 
 
   public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {
-    private double low;
-    private double high;
-
-    public double LowerBound => low;
-    public double UpperBound => high;
+    private MultivariateDual<Double> low;
+    private MultivariateDual<Double> high;
+
+    public MultivariateDual<Double> LowerBound => low.Clone();
+    public MultivariateDual<Double> UpperBound => high.Clone();
 
     public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }
 
+    public AlgebraicInterval(MultivariateDual<Double> low, MultivariateDual<Double> high) {
+      this.low = low.Clone();
+      this.high = high.Clone();
+    }
+
     public AlgebraicInterval(double low, double high) {
-      this.low = low;
-      this.high = high;
+      this.low = new MultivariateDual<Double>(new Double(low));
+      this.high = new MultivariateDual<Double>(new Double(high));
     }
 
     public AlgebraicInterval Add(AlgebraicInterval a) {
-      low += a.low;
-      high += a.high;
+      low.Add(a.low);
+      high.Add(a.high);
       return this;
     }
@@ -540 +750 @@
 
     public AlgebraicInterval Div(AlgebraicInterval a) {
-      if (Contains(0.0)) {
-        if (a.low.IsAlmost(0.0))
-          Mul(new AlgebraicInterval(1.0 / a.high, double.PositiveInfinity));
-        else if (a.high.IsAlmost(0.0))
-          Mul(new AlgebraicInterval(double.NegativeInfinity, 1.0 / a.low));
-        else {
-          low = double.NegativeInfinity;
-          high = double.PositiveInfinity;
-        }
+      if (a.Contains(0.0)) {
+        if (a.low.Value.Value.IsAlmost(0.0) && a.high.Value.Value.IsAlmost(0.0)) {
+          low = new MultivariateDual<Double>(double.NegativeInfinity);
+          high = new MultivariateDual<Double>(double.PositiveInfinity);
+        } else if (a.low.Value.Value.IsAlmost(0.0))
+          Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<Double>(double.PositiveInfinity)));
+        else
+          Mul(new AlgebraicInterval(new MultivariateDual<Double>(double.NegativeInfinity), a.low.Clone().Inv()));
       } else {
-        Mul(new AlgebraicInterval(1.0 / a.high, 1.0 / a.low));
+        Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
       }
       return this;
@@ -556 +765 @@
 
     public AlgebraicInterval AssignExp(AlgebraicInterval a) {
-      low = Math.Exp(a.low);
-      high = Math.Exp(a.high);
+      low.AssignExp(a.low);
+      high.AssignExp(a.high);
       return this;
     }
@@ -570 +779 @@
 
     public AlgebraicInterval AssignInv(AlgebraicInterval a) {
-      low = 1.0;
-      high = 1.0;
+      low = new MultivariateDual<Double>(1.0);
+      high = new MultivariateDual<Double>(1.0);
       return Div(a);
     }
 
     public AlgebraicInterval AssignLog(AlgebraicInterval a) {
-      low = Math.Log(a.low);
-      high = Math.Log(a.high);
+      low.AssignLog(a.low);
+      high.AssignLog(a.high);
       return this;
     }
 
     public AlgebraicInterval Mul(AlgebraicInterval a) {
-      double v1 = low * a.low;
-      double v2 = low * a.high;
-      double v3 = high * a.low;
-      double v4 = high * a.high;
-
-      low = Math.Min(Math.Min(v1, v2), Math.Min(v3, v4));
-      high = Math.Max(Math.Max(v1, v2), Math.Max(v3, v4));
+      var v1 = low.Clone().Mul(a.low);
+      var v2 = low.Clone().Mul(a.high);
+      var v3 = high.Clone().Mul(a.low);
+      var v4 = high.Clone().Mul(a.high);
+
+      low = Algebraic.Min(Algebraic.Min(v1, v2), Algebraic.Min(v3, v4));
+      high = Algebraic.Max(Algebraic.Max(v1, v2), Algebraic.Max(v3, v4));
       return this;
     }
 
     public AlgebraicInterval AssignNeg(AlgebraicInterval a) {
-      low = 0;
-      high = 0;
-      return Sub(a);
+      throw new NotImplementedException();
     }
 
     public AlgebraicInterval Scale(double s) {
+      low.Scale(s);
+      high.Scale(s);
       if (s < 0) {
-        low = s * high;
-        high = s * low;
-      } else {
-        low = s * low;
-        high = s * high;
+        var t = low;
+        low = high;
+        high = t;
       }
       return this;
@@ -614 +821 @@
 
     public AlgebraicInterval Sub(AlgebraicInterval a) {
-      low -= a.high;
-      high -= a.low;
+      // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1]
+      low.Sub(a.high);
+      high.Sub(a.low);
       return this;
     }
@@ -624 +832 @@
 
     public bool Contains(double val) {
-      return low <= val && val <= high;
-    }
-  }
+      return LowerBound.Value.Value <= val && val <= UpperBound.Value.Value;
+    }
+
+    public AlgebraicInterval AssignAbs(AlgebraicInterval a) {
+      if (a.Contains(0.0)) {
+        var abslow = a.low.Clone().Abs();
+        var abshigh = a.high.Clone().Abs();
+        a.high.Assign(Algebraic.Max(abslow, abshigh));
+        a.low.Assign(new MultivariateDual<Double>(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));
+      }
+      return this;
+    }
+
+    public AlgebraicInterval AssignSgn(AlgebraicInterval a) {
+      throw new NotImplementedException();
+    }
+  }
+
 
   public class Dual<V> : IAlgebraicType<Dual<V>>
@@ -633 +861 @@
     private V dv;
 
+    public V Value => v;
+    public V Derivative => dv;
+
     public Dual(V v, V dv) {
       this.v = v;
@@ -720 +951 @@
       throw new NotImplementedException();
     }
-  }
-
-  public class VectorDual<V> : IAlgebraicType<VectorDual<V>> where V : IAlgebraicType<V> {
+
+    public Dual<V> AssignAbs(Dual<V> a) {
+      v.AssignAbs(a.v);
+      // abs(f(x))' = f(x)*f'(x) / |f(x)|      
+      dv.Assign(a.dv).Mul(a.v.Clone().Sgn());
+      return this;
+    }
+
+    public Dual<V> AssignSgn(Dual<V> a) {
+      throw new NotImplementedException();
+    }
+  }
+
+  /// <summary>
+  /// An algebraic type which has a value as well as the partial derivatives of the value over multiple variables.
+  /// </summary>
+  /// <typeparam name="V"></typeparam>
+  public class MultivariateDual<V> : IAlgebraicType<MultivariateDual<V>> where V : IAlgebraicType<V>, new() {
     private V v;
     public V Value => v;
 
-    private Vector<V> dv;
-    public Vector<V> Gradient => dv; // column-oriented (dv[0] is the vector for the first parameter)
-
-
-    public VectorDual(V v, Vector<V> dv) {
-      this.v = v;
-      this.dv = dv;
-    }
-
-    public VectorDual<V> Clone() {
-      return new VectorDual<V>(v.Clone(), dv.Clone());
-    }
-
-    public VectorDual<V> Add(VectorDual<V> a) {
+    private SparseVector<object, V> dv;
+    public SparseVector<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 without partial derivative
+    /// </summary>
+    /// <param name="v"></param>
+    public MultivariateDual(V v) {
+      this.v = v.Clone();
+      this.dv = new SparseVector<object, V>();
+    }
+
+    /// <summary>
+    /// Constructor for multiple partial derivatives
+    /// </summary>
+    /// <param name="v"></param>
+    /// <param name="keys"></param>
+    /// <param name="dv"></param>
+    public MultivariateDual(V v, object[] keys, V[] dv) {
+      this.v = v.Clone();
+      this.dv = new SparseVector<object, V>(keys, dv);
+    }
+
+    /// <summary>
+    /// Constructor for a single partial derivative 
+    /// </summary>
+    /// <param name="v"></param>
+    /// <param name="key"></param>
+    /// <param name="dv"></param>
+    public MultivariateDual(V v, object key, V dv) {
+      this.v = v.Clone();
+      this.dv = new SparseVector<object, V>(new[] { key }, new[] { dv });
+    }
+
+    public MultivariateDual<V> Clone() {
+      return new MultivariateDual<V>(this);
+    }
+
+    public MultivariateDual<V> Add(MultivariateDual<V> a) {
       v.Add(a.v);
-      dv.Concat(a.dv);
-      return this;
-    }
-
-    public VectorDual<V> Assign(VectorDual<V> a) {
+      dv.Add(a.dv);
+      return this;
+    }
+
+    public MultivariateDual<V> Assign(MultivariateDual<V> a) {
       v.Assign(a.v);
-      dv = a.dv.Clone();
-      return this;
-    }
-
-    public VectorDual<V> AssignCos(VectorDual<V> a) {
-      v.AssignCos(a.v);
-      V t = default(V);
-      t.AssignNeg(t.AssignSin(a.v));
-      dv = a.dv.Clone().Scale(t);
-      return this;
-    }
-
-    public VectorDual<V> AssignExp(VectorDual<V> a) {
-      v.AssignExp(a.v);
-      dv = a.dv.Clone().Scale(v);
-      return this;
-    }
-
-    public VectorDual<V> AssignIntPower(VectorDual<V> a, int p) {
-      throw new NotImplementedException();
-    }
-
-    public VectorDual<V> AssignIntRoot(VectorDual<V> a, int r) {
-      throw new NotImplementedException();
-    }
-
-    public VectorDual<V> AssignInv(VectorDual<V> a) {
-      throw new NotImplementedException();
-    }
-
-    public VectorDual<V> AssignLog(VectorDual<V> a) {
-      throw new NotImplementedException();
-    }
-
-    public VectorDual<V> AssignNeg(VectorDual<V> a) {
-      throw new NotImplementedException();
-    }
-
-    public VectorDual<V> AssignSin(VectorDual<V> a) {
-      throw new NotImplementedException();
-    }
-
-    public VectorDual<V> Div(VectorDual<V> a) {
-      throw new NotImplementedException();
-    }
-
-    public VectorDual<V> Mul(VectorDual<V> a) {
+      dv.Assign(a.dv);
+      return this;
+    }
+
+    public MultivariateDual<V> AssignCos(MultivariateDual<V> a) {
+      throw new NotImplementedException();
+    }
+
+    public MultivariateDual<V> AssignExp(MultivariateDual<V> a) {
+      throw new NotImplementedException();
+    }
+
+    public MultivariateDual<V> AssignIntPower(MultivariateDual<V> a, int p) {
+      throw new NotImplementedException();
+    }
+
+    public MultivariateDual<V> AssignIntRoot(MultivariateDual<V> a, int r) {
+      throw new NotImplementedException();
+    }
+
+    public MultivariateDual<V> AssignInv(MultivariateDual<V> a) {
+      throw new NotImplementedException();
+    }
+
+    public MultivariateDual<V> AssignLog(MultivariateDual<V> a) {
+      throw new NotImplementedException();
+    }
+
+    public MultivariateDual<V> AssignNeg(MultivariateDual<V> a) {
+      throw new NotImplementedException();
+    }
+
+    public MultivariateDual<V> AssignSin(MultivariateDual<V> a) {
+      throw new NotImplementedException();
+    }
+
+    public MultivariateDual<V> Div(MultivariateDual<V> a) {
+      throw new NotImplementedException();
+    }
+
+    public MultivariateDual<V> Mul(MultivariateDual<V> a) {
       // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);
 
       var t1 = a.dv.Clone().Scale(v);
       var t2 = dv.Clone().Scale(a.v);
-      dv.Assign(t1).Concat(t2); // TODO: CHECK ORDER
+      dv.Assign(t1).Add(t2);
 
       v.Mul(a.v);
@@ -804 +1074 @@
     }
 
-    public VectorDual<V> Scale(double s) {
+    public MultivariateDual<V> Scale(double s) {
       v.Scale(s);
       dv.Scale(s);
@@ -810 +1080 @@
     }
 
-    public VectorDual<V> Sub(VectorDual<V> a) {
+    public MultivariateDual<V> Sub(MultivariateDual<V> a) {
       v.Sub(a.v);
-      dv.Concat(a.dv.AssignNeg(a.dv));
+      dv.Sub(a.dv);
+      return this;
+    }
+
+    public MultivariateDual<V> AssignAbs(MultivariateDual<V> a) {
+      v.AssignAbs(a.v);
+      // abs(f(x))' = f(x)*f'(x) / |f(x)|  doesn't work for intervals      
+
+      dv.Assign(a.dv).Scale(a.v.Clone().Sgn());
+      return this;
+    }
+
+    public MultivariateDual<V> AssignSgn(MultivariateDual<V> a) {
+      v.AssignSgn(a.v);
+      // sign(f(x)) = 0;
+      dv = new SparseVector<object, V>();
       return this;
     }

branches/2994-AutoDiffForIntervals/Tests/AutoDiffTest.cs

@@ -11 +11 @@
   public class AutoDiffTest {
     [TestMethod]
-    public void TestMethod1() {
+    public void Test() {
       {
+        // eval
         var parser = new InfixExpressionParser();
         var t = parser.Parse("2.0*x+y");
 
+        // interval eval
         var evaluator = new IntervalEvaluator();
-        var intervals = new Dictionary<string, AlgebraicInterval>();
-        intervals.Add("x", new AlgebraicInterval(-1.0, 1.0));
-        intervals.Add("y", new AlgebraicInterval(2.0, 10.0));
+        var intervals = new Dictionary<string, Interval>();
+        intervals.Add("x", new Interval(-1.0, 1.0));
+        intervals.Add("y", new Interval(2.0, 10.0));
         var resultInterval = evaluator.Evaluate(t, intervals);
         Assert.AreEqual(0, resultInterval.LowerBound);
@@ -26 +28 @@
 
       {
+        // vector eval
         var parser = new InfixExpressionParser();
         var t = parser.Parse("2.0*x+y");
@@ -54 +57 @@
 
       {
+        // vector eval and auto-diff
         var parser = new InfixExpressionParser();
         var t = parser.Parse("2.0*x+y");
@@ -98 +102 @@
       }
 
+      {
+        // interval eval and auto-diff
+        var parser = new InfixExpressionParser();
+        var t = parser.Parse("2.0*x+y");
+        var p0 = t.IterateNodesPostfix().First(n => n is ConstantTreeNode);
+        var p1 = t.IterateNodesPostfix().First(n => (n is VariableTreeNode var) && var.VariableName == "y");
+        var paramNodes = new ISymbolicExpressionTreeNode[] { p0, p1 };
+
+        var evaluator = new IntervalEvaluator();
+        var intervals = new Dictionary<string, Interval>();
+        intervals.Add("x", new Interval(-1.0, 1.0));
+        intervals.Add("y", new Interval(2.0, 10.0));
+        var resultInterval = evaluator.Evaluate(t, intervals, paramNodes, out double[] lowerGradient, out double[] upperGradient);
+        Assert.AreEqual(0, resultInterval.LowerBound);
+        Assert.AreEqual(12, resultInterval.UpperBound);
+
+        Assert.AreEqual(-1, lowerGradient[0]);
+        Assert.AreEqual(2, lowerGradient[1]);
+        Assert.AreEqual(1, upperGradient[0]);
+        Assert.AreEqual(10, upperGradient[1]);
+      }
+
     }
   }