Index: /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Algebraic.cs
===================================================================
--- /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Algebraic.cs	(revision 17295)
+++ /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Algebraic.cs	(revision 17295)
@@ -0,0 +1,26 @@
+﻿namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
+  public static class Algebraic {
+    public static T Abs<T>(this T a) where T : IAlgebraicType<T> { a.AssignAbs(a.Clone()); return a; }
+    public static T Neg<T>(this T a) where T : IAlgebraicType<T> { a.AssignNeg(a.Clone()); return a; }
+    public static T Inv<T>(this T a) where T : IAlgebraicType<T> { a.AssignInv(a.Clone()); return a; }
+    public static T Log<T>(this T a) where T : IAlgebraicType<T> { a.AssignLog(a.Clone()); return a; }
+    public static T Exp<T>(this T a) where T : IAlgebraicType<T> { a.AssignExp(a.Clone()); return a; }
+    public static T Sin<T>(this T a) where T : IAlgebraicType<T> { a.AssignSin(a.Clone()); return a; }
+    public static T Cos<T>(this T a) where T : IAlgebraicType<T> { a.AssignCos(a.Clone()); return a; }
+    public static T Sgn<T>(this T a) where T : IAlgebraicType<T> { a.AssignSgn(a.Clone()); return a; }
+    public static T IntPower<T>(this T a, int p) where T : IAlgebraicType<T> { a.AssignIntPower(a.Clone(), p); return a; }
+    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
+    //   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);
+    // }
+  }
+}
Index: /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/AlgebraicDouble.cs
===================================================================
--- /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/AlgebraicDouble.cs	(revision 17295)
+++ /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/AlgebraicDouble.cs	(revision 17295)
@@ -0,0 +1,54 @@
+﻿using System;
+using System.Diagnostics;
+
+namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
+  // algebraic type wrapper for a double value
+  [DebuggerDisplay("{Value}")]
+  public sealed class AlgebraicDouble : IAlgebraicType<AlgebraicDouble> {
+    public static implicit operator AlgebraicDouble(double value) { return new AlgebraicDouble(value); }
+    public static implicit operator double(AlgebraicDouble value) { return value.Value; }
+    public double Value;
+
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    public AlgebraicDouble Zero => new AlgebraicDouble(0.0);
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    public AlgebraicDouble One => new AlgebraicDouble(1.0);
+
+    public bool IsInfinity => IsNegativeInfinity || IsPositiveInfinity;
+    public bool IsNegativeInfinity => double.IsNegativeInfinity(Value);
+    public bool IsPositiveInfinity => double.IsPositiveInfinity(Value);
+    public AlgebraicDouble() { }
+    public AlgebraicDouble(double value) { this.Value = value; }
+    public AlgebraicDouble Assign(AlgebraicDouble a) { Value = a.Value; return this; }
+    public AlgebraicDouble Add(AlgebraicDouble a) { Value += a.Value; return this; }
+    public AlgebraicDouble Sub(AlgebraicDouble a) { Value -= a.Value; return this; }
+    public AlgebraicDouble Mul(AlgebraicDouble a) { Value *= a.Value; return this; }
+    public AlgebraicDouble Div(AlgebraicDouble a) { Value /= a.Value; return this; }
+    public AlgebraicDouble Scale(double s) { Value *= s; return this; }
+    public AlgebraicDouble AssignInv(AlgebraicDouble a) { Value = 1.0 / a.Value; return this; }
+    public AlgebraicDouble AssignNeg(AlgebraicDouble a) { Value = -a.Value; return this; }
+    public AlgebraicDouble AssignSin(AlgebraicDouble a) { Value = Math.Sin(a.Value); return this; }
+    public AlgebraicDouble AssignCos(AlgebraicDouble a) { Value = Math.Cos(a.Value); return this; }
+    public AlgebraicDouble AssignTanh(AlgebraicDouble a) { Value = Math.Tanh(a.Value); return this; }
+    public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = Math.Log(a.Value); return this; }
+    public AlgebraicDouble AssignExp(AlgebraicDouble a) { Value = Math.Exp(a.Value); return this; }
+    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; }
+
+    // helper
+    private static double IntRoot(double value, int r) {
+      if (r % 2 == 0) return Math.Pow(value, 1.0 / r);
+      else return value < 0 ? -Math.Pow(-value, 1.0 / r) : Math.Pow(value, 1.0 / r);
+    }
+
+    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; }
+    public AlgebraicDouble Clone() { return new AlgebraicDouble(Value); }
+
+    public override string ToString() {
+      return Value.ToString();
+    }
+  }
+}
Index: /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/AlgebraicDoubleVector.cs
===================================================================
--- /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/AlgebraicDoubleVector.cs	(revision 17295)
+++ /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/AlgebraicDoubleVector.cs	(revision 17295)
@@ -0,0 +1,85 @@
+﻿using System;
+using System.Diagnostics;
+
+namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
+  // a simple vector as an algebraic type
+  [DebuggerDisplay("DoubleVector(len={Length}): {string.}")]
+  public class AlgebraicDoubleVector : IAlgebraicType<AlgebraicDoubleVector> {
+    private double[] arr;
+    public double this[int idx] { get { return arr[idx]; } set { arr[idx] = value; } }
+    public int Length => arr.Length;
+
+    public AlgebraicDoubleVector(int length) { arr = new double[length]; }
+
+    public AlgebraicDoubleVector() { }
+
+    /// <summary>
+    /// 
+    /// </summary>
+    /// <param name="arr">array is not copied</param>
+    public AlgebraicDoubleVector(double[] arr) { this.arr = arr; }
+
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    public AlgebraicDoubleVector Zero => new AlgebraicDoubleVector(new double[this.Length]); // must return vector of same length as this (therefore Zero is not static)
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    public AlgebraicDoubleVector One => new AlgebraicDoubleVector(new double[this.Length]).AssignConstant(1.0); // must return vector of same length as this (therefore Zero is not static)
+    public AlgebraicDoubleVector Assign(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
+    public AlgebraicDoubleVector Add(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; }
+    public AlgebraicDoubleVector Sub(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] -= a.arr[i]; } return this; }
+    public AlgebraicDoubleVector Mul(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; }
+    public AlgebraicDoubleVector Div(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }
+    public AlgebraicDoubleVector AssignNeg(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = -a.arr[i]; } return this; }
+    public AlgebraicDoubleVector AssignInv(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }
+    public AlgebraicDoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; }
+    public AlgebraicDoubleVector AssignLog(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }
+    public AlgebraicDoubleVector AssignSin(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }
+    public AlgebraicDoubleVector AssignExp(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }
+    public AlgebraicDoubleVector AssignCos(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }
+    public AlgebraicDoubleVector AssignTanh(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Tanh(a.arr[i]); } return this; }
+    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; }
+
+    // helper
+    private double IntRoot(double v, int r) {
+      if (r % 2 == 0) return Math.Pow(v, 1.0 / r);
+      else return v < 0 ? -Math.Pow(-v, 1.0 / r) : Math.Pow(v, 1.0 / r);
+    }
+
+    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; }
+
+    public AlgebraicDoubleVector Clone() {
+      var v = new AlgebraicDoubleVector(this.arr.Length);
+      Array.Copy(arr, v.arr, v.arr.Length);
+      return v;
+    }
+
+    public AlgebraicDoubleVector AssignConstant(double constantValue) {
+      for (int i = 0; i < arr.Length; ++i) {
+        arr[i] = constantValue;
+      }
+      return this;
+    }
+
+    public void CopyTo(double[] dest, int idx, int length) {
+      Array.Copy(arr, 0, dest, idx, length);
+    }
+
+    public void CopyFrom(double[] data, int rowIndex) {
+      Array.Copy(data, rowIndex, arr, 0, Math.Min(arr.Length, data.Length - rowIndex));
+    }
+    public void CopyRowTo(double[,] dest, int row) {
+      for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j];
+    }
+
+    internal void CopyColumnTo(double[,] dest, int column, int row, int len) {
+      for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j];
+    }
+
+    public override string ToString() {
+      return "{" + string.Join(", ", arr.Take(Math.Max(5, arr.Length))) + (arr.Length > 5 ? "..." : string.Empty) + "}";
+    }
+  }
+}
Index: /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/AlgebraicInterval.cs
===================================================================
--- /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/AlgebraicInterval.cs	(revision 17295)
+++ /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/AlgebraicInterval.cs	(revision 17295)
@@ -0,0 +1,293 @@
+﻿using System;
+using System.Diagnostics;
+using System.Linq;
+
+namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
+  // this is our own implementation of interval arithmetic
+  // for a well worked out definition of interval operations for IEEE reals see:
+  // Stahl: Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations, Dissertation, JKU, 1995
+  [DebuggerDisplay("[{low.Value}..{high.Value}]")]
+  public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    private MultivariateDual<AlgebraicDouble> low;
+    public MultivariateDual<AlgebraicDouble> LowerBound => low.Clone();
+
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    private MultivariateDual<AlgebraicDouble> high;
+    public MultivariateDual<AlgebraicDouble> UpperBound => high.Clone();
+
+
+    public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }
+
+    public AlgebraicInterval(MultivariateDual<AlgebraicDouble> low, MultivariateDual<AlgebraicDouble> high) {
+      this.low = low.Clone();
+      this.high = high.Clone();
+    }
+
+    public AlgebraicInterval(double low, double high) {
+      this.low = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(low));
+      this.high = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(high));
+    }
+
+    [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);
+      high.Add(a.high);
+      return this;
+    }
+
+    public AlgebraicInterval Mul(AlgebraicInterval a) {
+      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 = 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;
+    }
+
+    public AlgebraicInterval Assign(AlgebraicInterval a) {
+      low = a.low;
+      high = a.high;
+      return this;
+    }
+
+    public AlgebraicInterval AssignCos(AlgebraicInterval a) {
+      return AssignSin(a.Clone().Add(new AlgebraicInterval(Math.PI / 2, Math.PI / 2)));
+    }
+
+    public AlgebraicInterval Div(AlgebraicInterval a) {
+      if (a.Contains(0.0)) {
+        if (a.low.Value.Value == 0 && a.high.Value.Value == 0) {
+          if (this.low.Value >= 0) {
+            // pos / 0
+          } else if (this.high.Value <= 0) {
+            // neg / 0
+          } else {
+            low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
+            high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
+          }
+        } else if (a.low.Value.Value >= 0) {
+          // a is positive 
+          Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));
+        } else if (a.high.Value <= 0) {
+          // a is negative
+          Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv()));
+        } else {
+          // a is interval over zero
+          low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
+          high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
+        }
+      } else {
+        Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
+      }
+      return this;
+    }
+
+    public AlgebraicInterval AssignExp(AlgebraicInterval a) {
+      low.AssignExp(a.low);
+      high.AssignExp(a.high);
+      return this;
+    }
+
+    // tanh is a bijective function
+    public AlgebraicInterval AssignTanh(AlgebraicInterval a) {
+      low.AssignTanh(a.low);
+      high.AssignTanh(a.high);
+      return this;
+    }
+
+    public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {
+      if (p < 0) {  // x^-3 == 1/(x^3)
+        AssignIntPower(a, -p);
+        return AssignInv(this);
+      } else if (p == 0) {
+        if (a.Contains(0.0)) {
+          // => 0^0 = 0 ; might not be relevant
+          low = new MultivariateDual<AlgebraicDouble>(0.0);
+          high = new MultivariateDual<AlgebraicDouble>(1.0);
+          return this;
+        } else {
+          // => 1
+          low = new MultivariateDual<AlgebraicDouble>(1.0);
+          high = new MultivariateDual<AlgebraicDouble>(1.0);
+          return this;
+        }
+      } else if (p == 1) return this;
+      else if (p % 2 == 0) {
+        // p is even => interval must be positive
+        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);
+        }
+      } else {
+        // p is uneven
+        if (a.Contains(0.0)) {
+          low.AssignIntPower(a.low, p);
+          high.AssignIntPower(a.high, p);
+        } else {
+          var lowPower = a.low.IntPower(p);
+          var highPower = a.high.IntPower(p);
+          low = lowPower.AssignMin(highPower);
+          high = lowPower.AssignMax(highPower);
+        }
+      }
+      return this;
+    }
+
+    public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {
+      if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; }
+      if (r == 1) return this;
+      if (r < 0) {
+        // x^ (-1/2) = 1 / (x^(1/2))
+        AssignIntRoot(a, -r);
+        return AssignInv(this);
+      } else {
+        // root only defined for positive arguments for even roots
+        if (r % 2 == 0 && a.LowerBound.Value.Value < 0) {
+          low = new MultivariateDual<AlgebraicDouble>(double.NaN);
+          high = new MultivariateDual<AlgebraicDouble>(double.NaN);
+          return this;
+        } else {
+          low.AssignIntRoot(a.low, r);
+          high.AssignIntRoot(a.high, r);
+          return this;
+        }
+      }
+    }
+
+    public AlgebraicInterval AssignInv(AlgebraicInterval a) {
+      low = new MultivariateDual<AlgebraicDouble>(1.0);
+      high = new MultivariateDual<AlgebraicDouble>(1.0);
+      return Div(a);
+    }
+
+    public AlgebraicInterval AssignLog(AlgebraicInterval a) {
+      low.AssignLog(a.low);
+      high.AssignLog(a.high);
+      return this;
+    }
+
+    public AlgebraicInterval AssignNeg(AlgebraicInterval a) {
+      low.AssignNeg(a.high);
+      high.AssignNeg(a.low);
+      return this;
+    }
+
+    public AlgebraicInterval Scale(double s) {
+      low.Scale(s);
+      high.Scale(s);
+      if (s < 0) {
+        var t = low;
+        low = high;
+        high = t;
+      }
+      return this;
+    }
+
+    public AlgebraicInterval AssignSin(AlgebraicInterval a) {
+      var lower = a.LowerBound.Value.Value;
+      var size = a.UpperBound.Value.Value - lower;
+      if (size < 0) throw new InvalidProgramException(); // ASSERT interval >= 0;
+
+      if (size >= Math.PI * 2) {
+        low = new MultivariateDual<AlgebraicDouble>(-1.0); // zero gradient
+        high = new MultivariateDual<AlgebraicDouble>(1.0);
+        return this;
+      }
+
+      // 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());
+
+      // override min and max if necessary
+
+      // shift interval 'a' into the range [-2pi .. 2pi] without changing the size of the interval to simplify the checks
+      lower = lower % (2 * Math.PI); // lower in [-2pi .. 2pi]      
+
+      // handle min = -1 and max = 1 cases explicitly
+      var pi_2 = Math.PI / 2.0;
+      var maxima = new double[] { -3 * pi_2, pi_2 };
+      var minima = new double[] { -pi_2, 3 * pi_2 };
+
+      // override min and max if necessary
+      if (maxima.Any(m => lower < m && lower + size > m)) {
+        // max = 1
+        high = new MultivariateDual<AlgebraicDouble>(1.0); // zero gradient
+      }
+
+      if (minima.Any(m => lower < m && lower + size > m)) {
+        // min = -1;
+        low = new MultivariateDual<AlgebraicDouble>(-1.0); // zero gradient
+      }
+      return this;
+    }
+
+    public AlgebraicInterval Sub(AlgebraicInterval a) {
+      // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1]
+      low.Sub(a.high);
+      high.Sub(a.low);
+      return this;
+    }
+
+    public AlgebraicInterval Clone() {
+      return new AlgebraicInterval(low, high);
+    }
+
+    public bool Contains(double val) {
+      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<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));
+      }
+      return this;
+    }
+
+    public AlgebraicInterval AssignSgn(AlgebraicInterval a) {
+      low.AssignSgn(a.low);
+      high.AssignSgn(a.high);
+      return this;
+    }
+
+    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;
+    }
+  }
+}
Index: /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/AlgebraicSparseVector.cs
===================================================================
--- /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/AlgebraicSparseVector.cs	(revision 17295)
+++ /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/AlgebraicSparseVector.cs	(revision 17295)
@@ -0,0 +1,211 @@
+﻿using System;
+using System.Collections.Generic;
+using System.Diagnostics;
+using System.Linq;
+
+namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
+  /*
+  // vectors of algebraic types
+  public sealed class AlgebraicVector<T> : IAlgebraicType<AlgebraicVector<T>> where T : IAlgebraicType<T>, new() {
+    private T[] elems;
+
+    public T this[int idx] { get { return elems[idx]; } set { elems[idx] = value; } }
+
+    public int Length => elems.Length;
+
+    private AlgebraicVector() { }
+
+    public AlgebraicVector(int len) { elems = new T[len]; }
+
+    /// <summary>
+    /// 
+    /// </summary>
+    /// <param name="elems">The array is copied (element-wise clone)</param>
+    public AlgebraicVector(T[] elems) {
+      this.elems = new T[elems.Length];
+      for (int i = 0; i < elems.Length; ++i) { this.elems[i] = elems[i].Clone(); }
+    }
+
+    /// <summary>
+    /// 
+    /// </summary>
+    /// <param name="elems">Array is not copied!</param>
+    /// <returns></returns>
+    public AlgebraicVector<T> FromArray(T[] elems) {
+      var v = new AlgebraicVector<T>();
+      v.elems = elems;
+      return v;
+    }
+
+    public void CopyTo(T[] dest) {
+      if (dest.Length != elems.Length) throw new InvalidOperationException("arr lengths do not match in Vector<T>.Copy");
+      Array.Copy(elems, dest, dest.Length);
+    }
+
+    public AlgebraicVector<T> Clone() { return new AlgebraicVector<T>(elems); }
+
+
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    public AlgebraicVector<T> Zero => new AlgebraicVector<T>(Length);
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    public AlgebraicVector<T> One { get { var v = new AlgebraicVector<T>(Length); for (int i = 0; i < elems.Length; ++i) elems[i] = new T().One; return v; } }
+    public AlgebraicVector<T> Assign(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }
+    public AlgebraicVector<T> Add(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; }
+    public AlgebraicVector<T> Sub(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }
+    public AlgebraicVector<T> Mul(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; }
+    public AlgebraicVector<T> Div(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; }
+    public AlgebraicVector<T> AssignNeg(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; }
+    public AlgebraicVector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; }
+    public AlgebraicVector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }
+    public AlgebraicVector<T> AssignInv(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; }
+    public AlgebraicVector<T> AssignLog(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; }
+    public AlgebraicVector<T> AssignExp(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; }
+    public AlgebraicVector<T> AssignSin(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; }
+    public AlgebraicVector<T> AssignCos(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignCos(a.elems[i]); } return this; }
+    public AlgebraicVector<T> AssignIntPower(AlgebraicVector<T> a, int p) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignIntPower(a.elems[i], p); } return this; }
+    public AlgebraicVector<T> AssignIntRoot(AlgebraicVector<T> a, int r) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignIntRoot(a.elems[i], r); } return this; }
+    public AlgebraicVector<T> AssignAbs(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignAbs(a.elems[i]); } return this; }
+    public AlgebraicVector<T> AssignSgn(AlgebraicVector<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>
+  [DebuggerDisplay("SparseVector: {ToString()}")]
+  public sealed class AlgebraicSparseVector<K, T> : IAlgebraicType<AlgebraicSparseVector<K, T>> where T : IAlgebraicType<T> {
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    private Dictionary<K, T> elems;
+    public IReadOnlyDictionary<K, T> Elements => elems;
+
+
+    public AlgebraicSparseVector(AlgebraicSparseVector<K, T> original) {
+      elems = original.elems.ToDictionary(kvp => kvp.Key, kvp => kvp.Value.Clone());
+    }
+
+    /// <summary>
+    /// 
+    /// </summary>
+    /// <param name="keys"></param>
+    /// <param name="values">values are cloned</param>
+    public AlgebraicSparseVector(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].Clone());
+      }
+    }
+
+    public AlgebraicSparseVector() {
+      this.elems = new Dictionary<K, T>();
+    }
+
+    // keep only elements from a
+    private void AssignFromSource(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {
+      // remove elems from this which don't occur in a
+      List<K> keysToRemove = new List<K>();
+      foreach (var kvp in elems) {
+        if (!a.elems.ContainsKey(kvp.Key)) keysToRemove.Add(kvp.Key);
+      }
+      foreach (var o in keysToRemove) elems.Remove(o); // -> zero
+
+      foreach (var kvp in a.elems) {
+        if (elems.TryGetValue(kvp.Key, out T value))
+          mapAssign(kvp.Value, value);
+        else
+          elems.Add(kvp.Key, mapAssign(kvp.Value, kvp.Value.Zero));
+      }
+    }
+
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    public AlgebraicSparseVector<K, T> Zero => new AlgebraicSparseVector<K, T>();
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    public AlgebraicSparseVector<K, T> One => throw new NotSupportedException();
+
+    public AlgebraicSparseVector<K, T> Scale(T s) { foreach (var kvp in elems) { kvp.Value.Mul(s); } return this; }
+    public AlgebraicSparseVector<K, T> Scale(double s) { foreach (var kvp in elems) { kvp.Value.Scale(s); } return this; }
+
+    public AlgebraicSparseVector<K, T> Assign(AlgebraicSparseVector<K, T> a) { elems.Clear(); AssignFromSource(a, (src, dest) => dest.Assign(src)); return this; }
+    public AlgebraicSparseVector<K, T> AssignInv(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignInv(src)); return this; }
+    public AlgebraicSparseVector<K, T> AssignNeg(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignNeg(src)); return this; }
+    public AlgebraicSparseVector<K, T> AssignLog(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignLog(src)); return this; }
+    public AlgebraicSparseVector<K, T> AssignExp(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignExp(src)); return this; }
+    public AlgebraicSparseVector<K, T> AssignIntPower(AlgebraicSparseVector<K, T> a, int p) { AssignFromSource(a, (src, dest) => dest.AssignIntPower(src, p)); return this; }
+    public AlgebraicSparseVector<K, T> AssignIntRoot(AlgebraicSparseVector<K, T> a, int r) { AssignFromSource(a, (src, dest) => dest.AssignIntRoot(src, r)); return this; }
+    public AlgebraicSparseVector<K, T> AssignSin(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSin(src)); return this; }
+    public AlgebraicSparseVector<K, T> AssignCos(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignCos(src)); return this; }
+    public AlgebraicSparseVector<K, T> AssignTanh(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignTanh(src)); return this; }
+    public AlgebraicSparseVector<K, T> AssignAbs(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignAbs(src)); return this; }
+    public AlgebraicSparseVector<K, T> AssignSgn(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSgn(src)); return this; }
+    public AlgebraicSparseVector<K, T> Add(AlgebraicSparseVector<K, T> a) {
+      foreach (var kvp in a.elems) {
+        if (elems.TryGetValue(kvp.Key, out T value))
+          value.Add(kvp.Value);
+        else
+          elems.Add(kvp.Key, kvp.Value.Clone()); // 0 + a
+      }
+      return this;
+    }
+
+    public AlgebraicSparseVector<K, T> Sub(AlgebraicSparseVector<K, T> a) {
+      foreach (var kvp in a.elems) {
+        if (elems.TryGetValue(kvp.Key, out T value))
+          value.Sub(kvp.Value);
+        else
+          elems.Add(kvp.Key, kvp.Value.Zero.Sub(kvp.Value)); // 0 - a
+      }
+      return this;
+    }
+
+    public AlgebraicSparseVector<K, T> Mul(AlgebraicSparseVector<K, T> a) {
+      var keys = elems.Keys.ToArray();
+      foreach (var k in keys) if (!a.elems.ContainsKey(k)) elems.Remove(k); // 0 * a => 0
+      foreach (var kvp in a.elems) {
+        if (elems.TryGetValue(kvp.Key, out T value))
+          value.Mul(kvp.Value); // this * a
+      }
+      return this;
+    }
+
+    public AlgebraicSparseVector<K, T> Div(AlgebraicSparseVector<K, T> a) {
+      return Mul(a.Clone().Inv());
+    }
+
+    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);
+    }
+
+    public override string ToString() {
+      return "[" + string.Join(" ", elems.Select(kvp => kvp.Key + ": " + kvp.Value)) + "]";
+    }
+  }
+}
Index: /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Dual.cs
===================================================================
--- /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Dual.cs	(revision 17295)
+++ /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Dual.cs	(revision 17295)
@@ -0,0 +1,63 @@
+﻿using System;
+using System.Diagnostics;
+
+namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
+  public class Dual<V> : IAlgebraicType<Dual<V>>
+    where V : IAlgebraicType<V> {
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    private V v;
+    public V Value => v;
+
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    private V dv;
+    public V Derivative => dv;
+
+    public Dual(V v, V dv) { this.v = v; this.dv = dv; }
+
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    public Dual<V> Zero => new Dual<V>(Value.Zero, Derivative.Zero);
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    public Dual<V> One => new Dual<V>(Value.One, Derivative.Zero);
+
+    public Dual<V> Assign(Dual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }
+    public Dual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }
+    public Dual<V> Add(Dual<V> a) { v.Add(a.v); dv.Add(a.dv); return this; }
+    public Dual<V> Sub(Dual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }
+    public Dual<V> AssignNeg(Dual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }
+    public Dual<V> AssignInv(Dual<V> a) { v.AssignInv(a.v); dv.AssignNeg(a.dv).Mul(v).Mul(v); return this; } // (1/f(x))' = - f(x)' / f(x)^2
+
+    // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);
+    public Dual<V> Mul(Dual<V> a) {
+      var t1 = a.dv.Clone().Mul(v);
+      var t2 = dv.Clone().Mul(a.v);
+      dv.Assign(t1).Add(t2);
+
+      v.Mul(a.v);
+      return this;
+    }
+    public Dual<V> Div(Dual<V> a) { Mul(a.Inv()); return this; }
+
+    public Dual<V> AssignExp(Dual<V> a) { v.AssignExp(a.v); dv.Assign(a.dv).Mul(v); return this; } // exp(f(x)) = exp(f(x))*f(x)'
+    public Dual<V> AssignLog(Dual<V> a) { v.AssignLog(a.v); dv.Assign(a.dv).Div(a.v); return this; }     // log(x)' = 1/f(x) * f(x)' 
+
+    public Dual<V> AssignIntPower(Dual<V> a, int p) { v.AssignIntPower(a.v, p); dv.Assign(a.dv).Scale(p).Mul(a.v.Clone().IntPower(p - 1)); return this; }
+    public Dual<V> AssignIntRoot(Dual<V> a, int r) { v.AssignIntRoot(a.v, r); dv.Assign(a.dv).Scale(1.0 / r).Mul(a.v.IntRoot(r - 1)); return this; }
+
+    public Dual<V> AssignSin(Dual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Mul(a.v.Clone().Cos()); return this; }
+    public Dual<V> AssignCos(Dual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Mul(a.v.Clone().Sin()); return this; }
+    public Dual<V> AssignTanh(Dual<V> a) { v.AssignTanh(a.v); dv.Assign(a.dv.Mul(v.Clone().IntPower(2).Neg().Add(Value.One))); return this; }
+
+    public Dual<V> AssignAbs(Dual<V> a) { v.AssignAbs(a.v); dv.Assign(a.dv).Mul(a.v.Clone().Sgn()); return this; }       // abs(f(x))' = f(x)*f'(x) / |f(x)|      
+    public Dual<V> AssignSgn(Dual<V> a) { v.AssignSgn(a.v); dv.Assign(a.dv.Zero); return this; }
+
+    public Dual<V> Clone() { return new Dual<V>(v.Clone(), dv.Clone()); }
+
+    public Dual<V> AssignMin(Dual<V> other) {
+      throw new NotImplementedException();
+    }
+
+    public Dual<V> AssignMax(Dual<V> other) {
+      throw new NotImplementedException();
+    }
+  }
+}
Index: /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/IAlgebraicType.cs
===================================================================
--- /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/IAlgebraicType.cs	(revision 17295)
+++ /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/IAlgebraicType.cs	(revision 17295)
@@ -0,0 +1,27 @@
+﻿namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
+  public interface IAlgebraicType<T> {
+    T Zero { get; } // Zero and One must create new objects
+    T One { get; }
+
+    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)
+    T AssignInv(T a); // set this to inv(a);
+    T Scale(double s); // scale this with s
+    T Add(T a); // add a to this
+    T Sub(T a); // subtract a from this
+    T Mul(T a); // multiply this with a
+    T Div(T a); // divide this by a
+    T AssignLog(T a); // set this to log a
+    T AssignExp(T a); // set this to exp(a)
+    T AssignSin(T a); // set this to sin(a)
+    T AssignCos(T a); // set this to cos(a)
+    T AssignTanh(T a); // set this to tanh(a)
+    T AssignIntPower(T a, int p);
+    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 Clone();
+  }
+}
Index: /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs
===================================================================
--- /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs	(revision 17294)
+++ /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs	(revision 17295)
@@ -1,6 +1,3 @@
 ﻿using System;
-using System.Collections.Generic;
-using System.Diagnostics;
-using System.Linq;
 using HeuristicLab.Common;
 using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
@@ -156,1104 +153,3 @@
 
   }
-
-
-  public sealed class VectorEvaluator : Interpreter<AlgebraicDoubleVector> {
-    private const int BATCHSIZE = 128;
-    [ThreadStatic]
-    private Dictionary<string, double[]> cachedData;
-
-    [ThreadStatic]
-    private IDataset dataset;
-
-    [ThreadStatic]
-    private int rowIndex;
-
-    [ThreadStatic]
-    private int[] rows;
-
-    private void InitCache(IDataset dataset) {
-      this.dataset = dataset;
-      cachedData = new Dictionary<string, double[]>();
-      foreach (var v in dataset.DoubleVariables) {
-        cachedData[v] = dataset.GetReadOnlyDoubleValues(v).ToArray();
-      }
-    }
-
-    public double[] Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows) {
-      if (cachedData == null || this.dataset != dataset) {
-        InitCache(dataset);
-      }
-
-      this.rows = rows;
-      var code = Compile(tree);
-      var remainingRows = rows.Length % BATCHSIZE;
-      var roundedTotal = rows.Length - remainingRows;
-
-      var result = new double[rows.Length];
-
-      for (rowIndex = 0; rowIndex < roundedTotal; rowIndex += BATCHSIZE) {
-        Evaluate(code);
-        code[0].value.CopyTo(result, rowIndex, BATCHSIZE);
-      }
-
-      if (remainingRows > 0) {
-        Evaluate(code);
-        code[0].value.CopyTo(result, roundedTotal, remainingRows);
-      }
-
-      return result;
-    }
-
-    protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
-      instruction.dblVal = constant.Value;
-      instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
-      instruction.value.AssignConstant(instruction.dblVal);
-    }
-
-    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
-      instruction.dblVal = variable.Weight;
-      instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
-      if (cachedData.ContainsKey(variable.VariableName)) {
-        instruction.data = cachedData[variable.VariableName];
-      } else {
-        instruction.data = dataset.GetDoubleValues(variable.VariableName).ToArray();
-        cachedData[variable.VariableName] = (double[])instruction.data;
-      }
-    }
-
-    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
-      instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
-    }
-
-    protected override void LoadVariable(Instruction a) {
-      var data = (double[])a.data;
-      for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) a.value[i - rowIndex] = data[rows[i]];
-      a.value.Scale(a.dblVal);
-    }
-  }
-
-  public sealed class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<AlgebraicDoubleVector>> {
-    private const int BATCHSIZE = 128;
-    [ThreadStatic]
-    private Dictionary<string, double[]> cachedData;
-
-    [ThreadStatic]
-    private IDataset dataset;
-
-    [ThreadStatic]
-    private int rowIndex;
-
-    [ThreadStatic]
-    private int[] rows;
-
-    [ThreadStatic]
-    private Dictionary<ISymbolicExpressionTreeNode, int> node2paramIdx;
-
-    private void InitCache(IDataset dataset) {
-      this.dataset = dataset;
-      cachedData = new Dictionary<string, double[]>();
-      foreach (var v in dataset.DoubleVariables) {
-        cachedData[v] = dataset.GetDoubleValues(v).ToArray();
-      }
-    }
-
-    /// <summary>
-    /// 
-    /// </summary>
-    /// <param name="tree"></param>
-    /// <param name="dataset"></param>
-    /// <param name="rows"></param>
-    /// <param name="parameterNodes"></param>
-    /// <param name="fi">Function output. Must be allocated by the caller.</param>
-    /// <param name="jac">Jacobian matrix. Must be allocated by the caller.</param>
-    public void Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows, ISymbolicExpressionTreeNode[] parameterNodes, double[] fi, double[,] jac) {
-      if (cachedData == null || this.dataset != dataset) {
-        InitCache(dataset);
-      }
-
-      int nParams = parameterNodes.Length;
-      node2paramIdx = new Dictionary<ISymbolicExpressionTreeNode, int>();
-      for (int i = 0; i < parameterNodes.Length; i++) node2paramIdx.Add(parameterNodes[i], i);
-
-      var code = Compile(tree);
-
-      var remainingRows = rows.Length % BATCHSIZE;
-      var roundedTotal = rows.Length - remainingRows;
-
-      this.rows = rows;
-
-      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;
-          }
-        }
-      }
-
-      if (remainingRows > 0) {
-        Evaluate(code);
-        code[0].value.Value.CopyTo(fi, roundedTotal, remainingRows);
-
-        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;
-          }
-      }
-    }
-
-    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
-      var zero = new AlgebraicDoubleVector(BATCHSIZE);
-      instruction.value = new MultivariateDual<AlgebraicDoubleVector>(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 AlgebraicDoubleVector(g_arr);
-        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE), paramIdx, g); // only a single column for the gradient
-      } else {
-        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE));
-      }
-
-      instruction.dblVal = constant.Value;
-      instruction.value.Value.AssignConstant(instruction.dblVal);
-    }
-
-    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
-      double[] data;
-      if (cachedData.ContainsKey(variable.VariableName)) {
-        data = cachedData[variable.VariableName];
-      } else {
-        data = dataset.GetReadOnlyDoubleValues(variable.VariableName).ToArray();
-        cachedData[variable.VariableName] = (double[])instruction.data;
-      }
-
-      var paramIdx = -1;
-      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);
-      } else {
-        var f = new AlgebraicDoubleVector(BATCHSIZE);
-        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(f);
-      }
-
-      instruction.dblVal = variable.Weight;
-      instruction.data = new object[] { data, paramIdx };
-    }
-
-    protected override void LoadVariable(Instruction a) {
-      var paramIdx = (int)((object[])a.data)[1];
-      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);
-
-      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]];
-        }
-      }
-    }
-  }
-
-
-  public sealed class IntervalEvaluator : Interpreter<AlgebraicInterval> {
-    [ThreadStatic]
-    private IDictionary<string, Interval> intervals;
-
-    public Interval Evaluate(ISymbolicExpressionTree tree, IDictionary<string, Interval> intervals) {
-      this.intervals = intervals;
-      var code = Compile(tree);
-      Evaluate(code);
-      if (code[0].value.LowerBound.Value.Value > code[0].value.UpperBound.Value.Value) throw new InvalidProgramException($"lower: {code[0].value.LowerBound.Value.Value} > upper: {code[0].value.UpperBound.Value.Value}");
-      return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);
-    }
-
-    public Interval Evaluate(ISymbolicExpressionTree tree, IDictionary<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) {
-        if (paramNodes[i] == null) continue;
-        if (l.Gradient.Elements.TryGetValue(paramNodes[i], out AlgebraicDouble value)) lowerGradient[i] = value;
-        if (u.Gradient.Elements.TryGetValue(paramNodes[i], out value)) upperGradient[i] = value;
-      }
-      return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);
-    }
-
-    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
-      instruction.value = new AlgebraicInterval(0, 0);
-    }
-
-
-    protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
-      instruction.dblVal = constant.Value;
-      instruction.value = new AlgebraicInterval(
-        new MultivariateDual<AlgebraicDouble>(instruction.dblVal, constant, 1.0),
-        new MultivariateDual<AlgebraicDouble>(instruction.dblVal, constant, 1.0) // use node as key
-        );
-    }
-
-    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
-      instruction.dblVal = variable.Weight;
-      var v1 = instruction.dblVal * intervals[variable.VariableName].LowerBound;
-      var v2 = instruction.dblVal * intervals[variable.VariableName].UpperBound;
-      var lower = Math.Min(v1, v2);
-      var upper = Math.Max(v1, v2);
-      // we assume that the for variable nodes ( v(x,w) = w * x ) the gradient is returned for parameter w
-      instruction.value = new AlgebraicInterval(
-        low: new MultivariateDual<AlgebraicDouble>(v: lower, key: variable, dv: intervals[variable.VariableName].LowerBound),
-        high: new MultivariateDual<AlgebraicDouble>(v: upper, key: variable, dv: intervals[variable.VariableName].UpperBound)
-        );
-    }
-
-    protected override void LoadVariable(Instruction a) {
-      // nothing to do
-    }
-  }
-
-  public interface IAlgebraicType<T> {
-    T Zero { get; } // Zero and One must create new objects
-    T One { get; }
-
-    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)
-    T AssignInv(T a); // set this to inv(a);
-    T Scale(double s); // scale this with s
-    T Add(T a); // add a to this
-    T Sub(T a); // subtract a from this
-    T Mul(T a); // multiply this with a
-    T Div(T a); // divide this by a
-    T AssignLog(T a); // set this to log a
-    T AssignExp(T a); // set this to exp(a)
-    T AssignSin(T a); // set this to sin(a)
-    T AssignCos(T a); // set this to cos(a)
-    T AssignTanh(T a); // set this to tanh(a)
-    T AssignIntPower(T a, int p);
-    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 Clone();
-  }
-
-  public static class Algebraic {
-    public static T Abs<T>(this T a) where T : IAlgebraicType<T> { a.AssignAbs(a.Clone()); return a; }
-    public static T Neg<T>(this T a) where T : IAlgebraicType<T> { a.AssignNeg(a.Clone()); return a; }
-    public static T Inv<T>(this T a) where T : IAlgebraicType<T> { a.AssignInv(a.Clone()); return a; }
-    public static T Log<T>(this T a) where T : IAlgebraicType<T> { a.AssignLog(a.Clone()); return a; }
-    public static T Exp<T>(this T a) where T : IAlgebraicType<T> { a.AssignExp(a.Clone()); return a; }
-    public static T Sin<T>(this T a) where T : IAlgebraicType<T> { a.AssignSin(a.Clone()); return a; }
-    public static T Cos<T>(this T a) where T : IAlgebraicType<T> { a.AssignCos(a.Clone()); return a; }
-    public static T Sgn<T>(this T a) where T : IAlgebraicType<T> { a.AssignSgn(a.Clone()); return a; }
-    public static T IntPower<T>(this T a, int p) where T : IAlgebraicType<T> { a.AssignIntPower(a.Clone(), p); return a; }
-    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
-    //   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
-  [DebuggerDisplay("{Value}")]
-  public sealed class AlgebraicDouble : IAlgebraicType<AlgebraicDouble> {
-    public static implicit operator AlgebraicDouble(double value) { return new AlgebraicDouble(value); }
-    public static implicit operator double(AlgebraicDouble value) { return value.Value; }
-    public double Value;
-
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    public AlgebraicDouble Zero => new AlgebraicDouble(0.0);
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    public AlgebraicDouble One => new AlgebraicDouble(1.0);
-
-    public bool IsInfinity => IsNegativeInfinity || IsPositiveInfinity;
-    public bool IsNegativeInfinity => double.IsNegativeInfinity(Value);
-    public bool IsPositiveInfinity => double.IsPositiveInfinity(Value);
-    public AlgebraicDouble() { }
-    public AlgebraicDouble(double value) { this.Value = value; }
-    public AlgebraicDouble Assign(AlgebraicDouble a) { Value = a.Value; return this; }
-    public AlgebraicDouble Add(AlgebraicDouble a) { Value += a.Value; return this; }
-    public AlgebraicDouble Sub(AlgebraicDouble a) { Value -= a.Value; return this; }
-    public AlgebraicDouble Mul(AlgebraicDouble a) { Value *= a.Value; return this; }
-    public AlgebraicDouble Div(AlgebraicDouble a) { Value /= a.Value; return this; }
-    public AlgebraicDouble Scale(double s) { Value *= s; return this; }
-    public AlgebraicDouble AssignInv(AlgebraicDouble a) { Value = 1.0 / a.Value; return this; }
-    public AlgebraicDouble AssignNeg(AlgebraicDouble a) { Value = -a.Value; return this; }
-    public AlgebraicDouble AssignSin(AlgebraicDouble a) { Value = Math.Sin(a.Value); return this; }
-    public AlgebraicDouble AssignCos(AlgebraicDouble a) { Value = Math.Cos(a.Value); return this; }
-    public AlgebraicDouble AssignTanh(AlgebraicDouble a) { Value = Math.Tanh(a.Value); return this; }
-    public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = Math.Log(a.Value); return this; }
-    public AlgebraicDouble AssignExp(AlgebraicDouble a) { Value = Math.Exp(a.Value); return this; }
-    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; }
-
-    // helper
-    private static double IntRoot(double value, int r) {
-      if (r % 2 == 0) return Math.Pow(value, 1.0 / r);
-      else return value < 0 ? -Math.Pow(-value, 1.0 / r) : Math.Pow(value, 1.0 / r);
-    }
-
-    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; }
-    public AlgebraicDouble Clone() { return new AlgebraicDouble(Value); }
-
-    public override string ToString() {
-      return Value.ToString();
-    }
-  }
-
-  // a simple vector as an algebraic type
-  [DebuggerDisplay("DoubleVector(len={Length}): {string.}")]
-  public class AlgebraicDoubleVector : IAlgebraicType<AlgebraicDoubleVector> {
-    private double[] arr;
-    public double this[int idx] { get { return arr[idx]; } set { arr[idx] = value; } }
-    public int Length => arr.Length;
-
-    public AlgebraicDoubleVector(int length) { arr = new double[length]; }
-
-    public AlgebraicDoubleVector() { }
-
-    /// <summary>
-    /// 
-    /// </summary>
-    /// <param name="arr">array is not copied</param>
-    public AlgebraicDoubleVector(double[] arr) { this.arr = arr; }
-
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    public AlgebraicDoubleVector Zero => new AlgebraicDoubleVector(new double[this.Length]); // must return vector of same length as this (therefore Zero is not static)
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    public AlgebraicDoubleVector One => new AlgebraicDoubleVector(new double[this.Length]).AssignConstant(1.0); // must return vector of same length as this (therefore Zero is not static)
-    public AlgebraicDoubleVector Assign(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
-    public AlgebraicDoubleVector Add(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; }
-    public AlgebraicDoubleVector Sub(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] -= a.arr[i]; } return this; }
-    public AlgebraicDoubleVector Mul(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; }
-    public AlgebraicDoubleVector Div(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }
-    public AlgebraicDoubleVector AssignNeg(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = -a.arr[i]; } return this; }
-    public AlgebraicDoubleVector AssignInv(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }
-    public AlgebraicDoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; }
-    public AlgebraicDoubleVector AssignLog(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }
-    public AlgebraicDoubleVector AssignSin(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }
-    public AlgebraicDoubleVector AssignExp(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }
-    public AlgebraicDoubleVector AssignCos(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }
-    public AlgebraicDoubleVector AssignTanh(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Tanh(a.arr[i]); } return this; }
-    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; }
-
-    // helper
-    private double IntRoot(double v, int r) {
-      if (r % 2 == 0) return Math.Pow(v, 1.0 / r);
-      else return v < 0 ? -Math.Pow(-v, 1.0 / r) : Math.Pow(v, 1.0 / r);
-    }
-
-    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; }
-
-    public AlgebraicDoubleVector Clone() {
-      var v = new AlgebraicDoubleVector(this.arr.Length);
-      Array.Copy(arr, v.arr, v.arr.Length);
-      return v;
-    }
-
-    public AlgebraicDoubleVector AssignConstant(double constantValue) {
-      for (int i = 0; i < arr.Length; ++i) {
-        arr[i] = constantValue;
-      }
-      return this;
-    }
-
-    public void CopyTo(double[] dest, int idx, int length) {
-      Array.Copy(arr, 0, dest, idx, length);
-    }
-
-    public void CopyFrom(double[] data, int rowIndex) {
-      Array.Copy(data, rowIndex, arr, 0, Math.Min(arr.Length, data.Length - rowIndex));
-    }
-    public void CopyRowTo(double[,] dest, int row) {
-      for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j];
-    }
-
-    internal void CopyColumnTo(double[,] dest, int column, int row, int len) {
-      for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j];
-    }
-
-    public override string ToString() {
-      return "{" + string.Join(", ", arr.Take(Math.Max(5, arr.Length))) + (arr.Length > 5 ? "..." : string.Empty) + "}";
-    }
-  }
-
-
-  /*
-  // vectors of algebraic types
-  public sealed class AlgebraicVector<T> : IAlgebraicType<AlgebraicVector<T>> where T : IAlgebraicType<T>, new() {
-    private T[] elems;
-
-    public T this[int idx] { get { return elems[idx]; } set { elems[idx] = value; } }
-
-    public int Length => elems.Length;
-
-    private AlgebraicVector() { }
-
-    public AlgebraicVector(int len) { elems = new T[len]; }
-
-    /// <summary>
-    /// 
-    /// </summary>
-    /// <param name="elems">The array is copied (element-wise clone)</param>
-    public AlgebraicVector(T[] elems) {
-      this.elems = new T[elems.Length];
-      for (int i = 0; i < elems.Length; ++i) { this.elems[i] = elems[i].Clone(); }
-    }
-
-    /// <summary>
-    /// 
-    /// </summary>
-    /// <param name="elems">Array is not copied!</param>
-    /// <returns></returns>
-    public AlgebraicVector<T> FromArray(T[] elems) {
-      var v = new AlgebraicVector<T>();
-      v.elems = elems;
-      return v;
-    }
-
-    public void CopyTo(T[] dest) {
-      if (dest.Length != elems.Length) throw new InvalidOperationException("arr lengths do not match in Vector<T>.Copy");
-      Array.Copy(elems, dest, dest.Length);
-    }
-
-    public AlgebraicVector<T> Clone() { return new AlgebraicVector<T>(elems); }
-
-
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    public AlgebraicVector<T> Zero => new AlgebraicVector<T>(Length);
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    public AlgebraicVector<T> One { get { var v = new AlgebraicVector<T>(Length); for (int i = 0; i < elems.Length; ++i) elems[i] = new T().One; return v; } }
-    public AlgebraicVector<T> Assign(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }
-    public AlgebraicVector<T> Add(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; }
-    public AlgebraicVector<T> Sub(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }
-    public AlgebraicVector<T> Mul(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; }
-    public AlgebraicVector<T> Div(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; }
-    public AlgebraicVector<T> AssignNeg(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; }
-    public AlgebraicVector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; }
-    public AlgebraicVector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }
-    public AlgebraicVector<T> AssignInv(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; }
-    public AlgebraicVector<T> AssignLog(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; }
-    public AlgebraicVector<T> AssignExp(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; }
-    public AlgebraicVector<T> AssignSin(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; }
-    public AlgebraicVector<T> AssignCos(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignCos(a.elems[i]); } return this; }
-    public AlgebraicVector<T> AssignIntPower(AlgebraicVector<T> a, int p) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignIntPower(a.elems[i], p); } return this; }
-    public AlgebraicVector<T> AssignIntRoot(AlgebraicVector<T> a, int r) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignIntRoot(a.elems[i], r); } return this; }
-    public AlgebraicVector<T> AssignAbs(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignAbs(a.elems[i]); } return this; }
-    public AlgebraicVector<T> AssignSgn(AlgebraicVector<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>
-  [DebuggerDisplay("SparseVector: {ToString()}")]
-  public sealed class AlgebraicSparseVector<K, T> : IAlgebraicType<AlgebraicSparseVector<K, T>> where T : IAlgebraicType<T> {
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    private Dictionary<K, T> elems;
-    public IReadOnlyDictionary<K, T> Elements => elems;
-
-
-    public AlgebraicSparseVector(AlgebraicSparseVector<K, T> original) {
-      elems = original.elems.ToDictionary(kvp => kvp.Key, kvp => kvp.Value.Clone());
-    }
-
-    /// <summary>
-    /// 
-    /// </summary>
-    /// <param name="keys"></param>
-    /// <param name="values">values are cloned</param>
-    public AlgebraicSparseVector(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].Clone());
-      }
-    }
-
-    public AlgebraicSparseVector() {
-      this.elems = new Dictionary<K, T>();
-    }
-
-    // keep only elements from a
-    private void AssignFromSource(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {
-      // remove elems from this which don't occur in a
-      List<K> keysToRemove = new List<K>();
-      foreach (var kvp in elems) {
-        if (!a.elems.ContainsKey(kvp.Key)) keysToRemove.Add(kvp.Key);
-      }
-      foreach (var o in keysToRemove) elems.Remove(o); // -> zero
-
-      foreach (var kvp in a.elems) {
-        if (elems.TryGetValue(kvp.Key, out T value))
-          mapAssign(kvp.Value, value);
-        else
-          elems.Add(kvp.Key, mapAssign(kvp.Value, kvp.Value.Zero));
-      }
-    }
-
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    public AlgebraicSparseVector<K, T> Zero => new AlgebraicSparseVector<K, T>();
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    public AlgebraicSparseVector<K, T> One => throw new NotSupportedException();
-
-    public AlgebraicSparseVector<K, T> Scale(T s) { foreach (var kvp in elems) { kvp.Value.Mul(s); } return this; }
-    public AlgebraicSparseVector<K, T> Scale(double s) { foreach (var kvp in elems) { kvp.Value.Scale(s); } return this; }
-
-    public AlgebraicSparseVector<K, T> Assign(AlgebraicSparseVector<K, T> a) { elems.Clear(); AssignFromSource(a, (src, dest) => dest.Assign(src)); return this; }
-    public AlgebraicSparseVector<K, T> AssignInv(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignInv(src)); return this; }
-    public AlgebraicSparseVector<K, T> AssignNeg(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignNeg(src)); return this; }
-    public AlgebraicSparseVector<K, T> AssignLog(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignLog(src)); return this; }
-    public AlgebraicSparseVector<K, T> AssignExp(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignExp(src)); return this; }
-    public AlgebraicSparseVector<K, T> AssignIntPower(AlgebraicSparseVector<K, T> a, int p) { AssignFromSource(a, (src, dest) => dest.AssignIntPower(src, p)); return this; }
-    public AlgebraicSparseVector<K, T> AssignIntRoot(AlgebraicSparseVector<K, T> a, int r) { AssignFromSource(a, (src, dest) => dest.AssignIntRoot(src, r)); return this; }
-    public AlgebraicSparseVector<K, T> AssignSin(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSin(src)); return this; }
-    public AlgebraicSparseVector<K, T> AssignCos(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignCos(src)); return this; }
-    public AlgebraicSparseVector<K, T> AssignTanh(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignTanh(src)); return this; }
-    public AlgebraicSparseVector<K, T> AssignAbs(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignAbs(src)); return this; }
-    public AlgebraicSparseVector<K, T> AssignSgn(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSgn(src)); return this; }
-    public AlgebraicSparseVector<K, T> Add(AlgebraicSparseVector<K, T> a) {
-      foreach (var kvp in a.elems) {
-        if (elems.TryGetValue(kvp.Key, out T value))
-          value.Add(kvp.Value);
-        else
-          elems.Add(kvp.Key, kvp.Value.Clone()); // 0 + a
-      }
-      return this;
-    }
-
-    public AlgebraicSparseVector<K, T> Sub(AlgebraicSparseVector<K, T> a) {
-      foreach (var kvp in a.elems) {
-        if (elems.TryGetValue(kvp.Key, out T value))
-          value.Sub(kvp.Value);
-        else
-          elems.Add(kvp.Key, kvp.Value.Zero.Sub(kvp.Value)); // 0 - a
-      }
-      return this;
-    }
-
-    public AlgebraicSparseVector<K, T> Mul(AlgebraicSparseVector<K, T> a) {
-      var keys = elems.Keys.ToArray();
-      foreach (var k in keys) if (!a.elems.ContainsKey(k)) elems.Remove(k); // 0 * a => 0
-      foreach (var kvp in a.elems) {
-        if (elems.TryGetValue(kvp.Key, out T value))
-          value.Mul(kvp.Value); // this * a
-      }
-      return this;
-    }
-
-    public AlgebraicSparseVector<K, T> Div(AlgebraicSparseVector<K, T> a) {
-      return Mul(a.Clone().Inv());
-    }
-
-    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);
-    }
-
-    public override string ToString() {
-      return "[" + string.Join(" ", elems.Select(kvp => kvp.Key + ": " + kvp.Value)) + "]";
-    }
-  }
-
-  // this is our own implementation of interval arithmetic
-  // for a well worked out definition of interval operations for IEEE reals see:
-  // Stahl: Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations, Dissertation, JKU, 1995
-  [DebuggerDisplay("[{low.Value}..{high.Value}]")]
-  public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    private MultivariateDual<AlgebraicDouble> low;
-    public MultivariateDual<AlgebraicDouble> LowerBound => low.Clone();
-
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    private MultivariateDual<AlgebraicDouble> high;
-    public MultivariateDual<AlgebraicDouble> UpperBound => high.Clone();
-
-
-    public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }
-
-    public AlgebraicInterval(MultivariateDual<AlgebraicDouble> low, MultivariateDual<AlgebraicDouble> high) {
-      this.low = low.Clone();
-      this.high = high.Clone();
-    }
-
-    public AlgebraicInterval(double low, double high) {
-      this.low = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(low));
-      this.high = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(high));
-    }
-
-    [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);
-      high.Add(a.high);
-      return this;
-    }
-
-    public AlgebraicInterval Mul(AlgebraicInterval a) {
-      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 = 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;
-    }
-
-    public AlgebraicInterval Assign(AlgebraicInterval a) {
-      low = a.low;
-      high = a.high;
-      return this;
-    }
-
-    public AlgebraicInterval AssignCos(AlgebraicInterval a) {
-      return AssignSin(a.Clone().Add(new AlgebraicInterval(Math.PI / 2, Math.PI / 2)));
-    }
-
-    public AlgebraicInterval Div(AlgebraicInterval a) {
-      if (a.Contains(0.0)) {
-        if (a.low.Value.Value == 0 && a.high.Value.Value == 0) {
-          if (this.low.Value >= 0) {
-            // pos / 0
-          } else if (this.high.Value <= 0) {
-            // neg / 0
-          } else {
-            low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
-            high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
-          }
-        } else if (a.low.Value.Value >= 0) {
-          // a is positive 
-          Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));
-        } else if (a.high.Value <= 0) {
-          // a is negative
-          Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv()));
-        } else {
-          // a is interval over zero
-          low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
-          high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
-        }
-      } else {
-        Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
-      }
-      return this;
-    }
-
-    public AlgebraicInterval AssignExp(AlgebraicInterval a) {
-      low.AssignExp(a.low);
-      high.AssignExp(a.high);
-      return this;
-    }
-
-    // tanh is a bijective function
-    public AlgebraicInterval AssignTanh(AlgebraicInterval a) {
-      low.AssignTanh(a.low);
-      high.AssignTanh(a.high);
-      return this;
-    }
-
-    public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {
-      if (p < 0) {  // x^-3 == 1/(x^3)
-        AssignIntPower(a, -p);
-        return AssignInv(this);
-      } else if (p == 0) {
-        if (a.Contains(0.0)) {
-          // => 0^0 = 0 ; might not be relevant
-          low = new MultivariateDual<AlgebraicDouble>(0.0);
-          high = new MultivariateDual<AlgebraicDouble>(1.0);
-          return this;
-        } else {
-          // => 1
-          low = new MultivariateDual<AlgebraicDouble>(1.0);
-          high = new MultivariateDual<AlgebraicDouble>(1.0);
-          return this;
-        }
-      } else if (p == 1) return this;
-      else if (p % 2 == 0) {
-        // p is even => interval must be positive
-        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);
-        }
-      } else {
-        // p is uneven
-        if (a.Contains(0.0)) {
-          low.AssignIntPower(a.low, p);
-          high.AssignIntPower(a.high, p);
-        } else {
-          var lowPower = a.low.IntPower(p);
-          var highPower = a.high.IntPower(p);
-          low = lowPower.AssignMin(highPower);
-          high = lowPower.AssignMax(highPower);
-        }
-      }
-      return this;
-    }
-
-    public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {
-      if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; }
-      if (r == 1) return this;
-      if (r < 0) {
-        // x^ (-1/2) = 1 / (x^(1/2))
-        AssignIntRoot(a, -r);
-        return AssignInv(this);
-      } else {
-        // root only defined for positive arguments for even roots
-        if (r % 2 == 0 && a.LowerBound.Value.Value < 0) {
-          low = new MultivariateDual<AlgebraicDouble>(double.NaN);
-          high = new MultivariateDual<AlgebraicDouble>(double.NaN);
-          return this;
-        } else {
-          low.AssignIntRoot(a.low, r);
-          high.AssignIntRoot(a.high, r);
-          return this;
-        }
-      }
-    }
-
-    public AlgebraicInterval AssignInv(AlgebraicInterval a) {
-      low = new MultivariateDual<AlgebraicDouble>(1.0);
-      high = new MultivariateDual<AlgebraicDouble>(1.0);
-      return Div(a);
-    }
-
-    public AlgebraicInterval AssignLog(AlgebraicInterval a) {
-      low.AssignLog(a.low);
-      high.AssignLog(a.high);
-      return this;
-    }
-
-    public AlgebraicInterval AssignNeg(AlgebraicInterval a) {
-      low.AssignNeg(a.high);
-      high.AssignNeg(a.low);
-      return this;
-    }
-
-    public AlgebraicInterval Scale(double s) {
-      low.Scale(s);
-      high.Scale(s);
-      if (s < 0) {
-        var t = low;
-        low = high;
-        high = t;
-      }
-      return this;
-    }
-
-    public AlgebraicInterval AssignSin(AlgebraicInterval a) {
-      var lower = a.LowerBound.Value.Value;
-      var size = a.UpperBound.Value.Value - lower;
-      if (size < 0) throw new InvalidProgramException(); // ASSERT interval >= 0;
-
-      if (size >= Math.PI * 2) {
-        low = new MultivariateDual<AlgebraicDouble>(-1.0); // zero gradient
-        high = new MultivariateDual<AlgebraicDouble>(1.0);
-        return this;
-      }
-
-      // 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());
-
-      // override min and max if necessary
-
-      // shift interval 'a' into the range [-2pi .. 2pi] without changing the size of the interval to simplify the checks
-      lower = lower % (2 * Math.PI); // lower in [-2pi .. 2pi]      
-
-      // handle min = -1 and max = 1 cases explicitly
-      var pi_2 = Math.PI / 2.0;
-      var maxima = new double[] { -3 * pi_2, pi_2 };
-      var minima = new double[] { -pi_2, 3 * pi_2 };
-
-      // override min and max if necessary
-      if (maxima.Any(m => lower < m && lower + size > m)) {
-        // max = 1
-        high = new MultivariateDual<AlgebraicDouble>(1.0); // zero gradient
-      }
-
-      if (minima.Any(m => lower < m && lower + size > m)) {
-        // min = -1;
-        low = new MultivariateDual<AlgebraicDouble>(-1.0); // zero gradient
-      }
-      return this;
-    }
-
-    public AlgebraicInterval Sub(AlgebraicInterval a) {
-      // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1]
-      low.Sub(a.high);
-      high.Sub(a.low);
-      return this;
-    }
-
-    public AlgebraicInterval Clone() {
-      return new AlgebraicInterval(low, high);
-    }
-
-    public bool Contains(double val) {
-      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<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));
-      }
-      return this;
-    }
-
-    public AlgebraicInterval AssignSgn(AlgebraicInterval a) {
-      low.AssignSgn(a.low);
-      high.AssignSgn(a.high);
-      return this;
-    }
-
-    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;
-    }
-  }
-
-  public class Dual<V> : IAlgebraicType<Dual<V>>
-    where V : IAlgebraicType<V> {
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    private V v;
-    public V Value => v;
-
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    private V dv;
-    public V Derivative => dv;
-
-    public Dual(V v, V dv) { this.v = v; this.dv = dv; }
-
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    public Dual<V> Zero => new Dual<V>(Value.Zero, Derivative.Zero);
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    public Dual<V> One => new Dual<V>(Value.One, Derivative.Zero);
-
-    public Dual<V> Assign(Dual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }
-    public Dual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }
-    public Dual<V> Add(Dual<V> a) { v.Add(a.v); dv.Add(a.dv); return this; }
-    public Dual<V> Sub(Dual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }
-    public Dual<V> AssignNeg(Dual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }
-    public Dual<V> AssignInv(Dual<V> a) { v.AssignInv(a.v); dv.AssignNeg(a.dv).Mul(v).Mul(v); return this; } // (1/f(x))' = - f(x)' / f(x)^2
-
-    // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);
-    public Dual<V> Mul(Dual<V> a) {
-      var t1 = a.dv.Clone().Mul(v);
-      var t2 = dv.Clone().Mul(a.v);
-      dv.Assign(t1).Add(t2);
-
-      v.Mul(a.v);
-      return this;
-    }
-    public Dual<V> Div(Dual<V> a) { Mul(a.Inv()); return this; }
-
-    public Dual<V> AssignExp(Dual<V> a) { v.AssignExp(a.v); dv.Assign(a.dv).Mul(v); return this; } // exp(f(x)) = exp(f(x))*f(x)'
-    public Dual<V> AssignLog(Dual<V> a) { v.AssignLog(a.v); dv.Assign(a.dv).Div(a.v); return this; }     // log(x)' = 1/f(x) * f(x)' 
-
-    public Dual<V> AssignIntPower(Dual<V> a, int p) { v.AssignIntPower(a.v, p); dv.Assign(a.dv).Scale(p).Mul(a.v.Clone().IntPower(p - 1)); return this; }
-    public Dual<V> AssignIntRoot(Dual<V> a, int r) { v.AssignIntRoot(a.v, r); dv.Assign(a.dv).Scale(1.0 / r).Mul(a.v.IntRoot(r - 1)); return this; }
-
-    public Dual<V> AssignSin(Dual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Mul(a.v.Clone().Cos()); return this; }
-    public Dual<V> AssignCos(Dual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Mul(a.v.Clone().Sin()); return this; }
-    public Dual<V> AssignTanh(Dual<V> a) { v.AssignTanh(a.v); dv.Assign(a.dv.Mul(v.Clone().IntPower(2).Neg().Add(Value.One))); return this; }
-
-    public Dual<V> AssignAbs(Dual<V> a) { v.AssignAbs(a.v); dv.Assign(a.dv).Mul(a.v.Clone().Sgn()); return this; }       // abs(f(x))' = f(x)*f'(x) / |f(x)|      
-    public Dual<V> AssignSgn(Dual<V> a) { v.AssignSgn(a.v); dv.Assign(a.dv.Zero); return this; }
-
-    public Dual<V> Clone() { return new Dual<V>(v.Clone(), dv.Clone()); }
-
-    public Dual<V> AssignMin(Dual<V> other) {
-      throw new NotImplementedException();
-    }
-
-    public Dual<V> AssignMax(Dual<V> other) {
-      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>
-  [DebuggerDisplay("v={Value}; dv={dv}")]
-  public class MultivariateDual<V> : IAlgebraicType<MultivariateDual<V>> where V : IAlgebraicType<V>, new() {
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    private V v;
-    public V Value => v;
-
-    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
-    private AlgebraicSparseVector<object, V> dv;
-    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 without partial derivative
-    /// </summary>
-    /// <param name="v"></param>
-    public MultivariateDual(V v) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<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 AlgebraicSparseVector<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 AlgebraicSparseVector<object, V>(new[] { key }, new[] { dv }); }
-
-    /// <summary>
-    /// Constructor with a given value and gradient. For internal use.
-    /// </summary>
-    /// <param name="v">The value (not cloned).</param>
-    /// <param name="gradient">The gradient (not cloned).</param>
-    internal MultivariateDual(V v, AlgebraicSparseVector<object, V> gradient) { this.v = v; this.dv = gradient; }
-
-    public MultivariateDual<V> Clone() { return new MultivariateDual<V>(this); }
-
-    public MultivariateDual<V> Zero => new MultivariateDual<V>(Value.Zero, Gradient.Zero);
-    public MultivariateDual<V> One => new MultivariateDual<V>(Value.One, Gradient.Zero);
-
-    public MultivariateDual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }
-
-    public MultivariateDual<V> Add(MultivariateDual<V> a) { v.Add(a.v); dv.Add(a.dv); return this; }
-    public MultivariateDual<V> Sub(MultivariateDual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }
-    public MultivariateDual<V> Assign(MultivariateDual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }
-    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).Add(t2);
-
-      v.Mul(a.v);
-      return this;
-    }
-    public MultivariateDual<V> Div(MultivariateDual<V> a) { v.Div(a.v); dv.Mul(a.dv.Inv()); return this; }
-    public MultivariateDual<V> AssignNeg(MultivariateDual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }
-    public MultivariateDual<V> AssignInv(MultivariateDual<V> a) { v.AssignInv(a.v); dv.AssignNeg(a.dv).Scale(v).Scale(v); return this; }   // (1/f(x))' = - f(x)' / f(x)^2
-
-    public MultivariateDual<V> AssignSin(MultivariateDual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Cos()); return this; }
-    public MultivariateDual<V> AssignCos(MultivariateDual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Scale(a.v.Clone().Sin()); return this; }
-    public MultivariateDual<V> AssignTanh(MultivariateDual<V> a) { v.AssignTanh(a.v); dv.Assign(a.dv.Scale(v.Clone().IntPower(2).Neg().Add(Value.One))); return this; }     // tanh(f(x))' = f(x)'sech²(f(x)) = f(x)'(1 - tanh²(f(x)))
-
-    public MultivariateDual<V> AssignIntPower(MultivariateDual<V> a, int p) { v.AssignIntPower(a.v, p); dv.Assign(a.dv).Scale(p).Scale(a.v.Clone().IntPower(p - 1)); return this; }
-    public MultivariateDual<V> AssignIntRoot(MultivariateDual<V> a, int r) { v.AssignIntRoot(a.v, r); dv.Assign(a.dv).Scale(1.0 / r).Scale(a.v.IntRoot(r - 1)); return this; }
-
-    public MultivariateDual<V> AssignExp(MultivariateDual<V> a) { v.AssignExp(a.v); dv.Assign(a.dv).Scale(v); return this; } // exp(f(x)) = exp(f(x))*f(x)'     
-    public MultivariateDual<V> AssignLog(MultivariateDual<V> a) { v.AssignLog(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Inv()); return this; }     // log(x)' = 1/f(x) * f(x)' 
-
-    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) {
-      XXX
-    }
-
-    public MultivariateDual<V> AssignMax(MultivariateDual<V> other) {
-      XXX
-    }
-  }
 }
Index: /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/IntervalEvaluator.cs
===================================================================
--- /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/IntervalEvaluator.cs	(revision 17295)
+++ /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/IntervalEvaluator.cs	(revision 17295)
@@ -0,0 +1,64 @@
+﻿using System;
+using System.Collections.Generic;
+using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
+
+namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
+  public sealed class IntervalEvaluator : Interpreter<AlgebraicInterval> {
+    [ThreadStatic]
+    private IDictionary<string, Interval> intervals;
+
+    public Interval Evaluate(ISymbolicExpressionTree tree, IDictionary<string, Interval> intervals) {
+      this.intervals = intervals;
+      var code = Compile(tree);
+      Evaluate(code);
+      if (code[0].value.LowerBound.Value.Value > code[0].value.UpperBound.Value.Value) throw new InvalidProgramException($"lower: {code[0].value.LowerBound.Value.Value} > upper: {code[0].value.UpperBound.Value.Value}");
+      return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);
+    }
+
+    public Interval Evaluate(ISymbolicExpressionTree tree, IDictionary<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) {
+        if (paramNodes[i] == null) continue;
+        if (l.Gradient.Elements.TryGetValue(paramNodes[i], out AlgebraicDouble value)) lowerGradient[i] = value;
+        if (u.Gradient.Elements.TryGetValue(paramNodes[i], out value)) upperGradient[i] = value;
+      }
+      return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);
+    }
+
+    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
+      instruction.value = new AlgebraicInterval(0, 0);
+    }
+
+
+    protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
+      instruction.dblVal = constant.Value;
+      instruction.value = new AlgebraicInterval(
+        new MultivariateDual<AlgebraicDouble>(instruction.dblVal, constant, 1.0),
+        new MultivariateDual<AlgebraicDouble>(instruction.dblVal, constant, 1.0) // use node as key
+        );
+    }
+
+    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
+      instruction.dblVal = variable.Weight;
+      var v1 = instruction.dblVal * intervals[variable.VariableName].LowerBound;
+      var v2 = instruction.dblVal * intervals[variable.VariableName].UpperBound;
+      var lower = Math.Min(v1, v2);
+      var upper = Math.Max(v1, v2);
+      // we assume that the for variable nodes ( v(x,w) = w * x ) the gradient is returned for parameter w
+      instruction.value = new AlgebraicInterval(
+        low: new MultivariateDual<AlgebraicDouble>(v: lower, key: variable, dv: intervals[variable.VariableName].LowerBound),
+        high: new MultivariateDual<AlgebraicDouble>(v: upper, key: variable, dv: intervals[variable.VariableName].UpperBound)
+        );
+    }
+
+    protected override void LoadVariable(Instruction a) {
+      // nothing to do
+    }
+  }
+}
Index: /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/MultivariateDual.cs
===================================================================
--- /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/MultivariateDual.cs	(revision 17295)
+++ /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/MultivariateDual.cs	(revision 17295)
@@ -0,0 +1,93 @@
+﻿using System.Diagnostics;
+
+namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
+  /// <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>
+  [DebuggerDisplay("v={Value}; dv={dv}")]
+  public class MultivariateDual<V> : IAlgebraicType<MultivariateDual<V>> where V : IAlgebraicType<V>, new() {
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    private V v;
+    public V Value => v;
+
+    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
+    private AlgebraicSparseVector<object, V> dv;
+    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 without partial derivative
+    /// </summary>
+    /// <param name="v"></param>
+    public MultivariateDual(V v) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<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 AlgebraicSparseVector<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 AlgebraicSparseVector<object, V>(new[] { key }, new[] { dv }); }
+
+    /// <summary>
+    /// Constructor with a given value and gradient. For internal use.
+    /// </summary>
+    /// <param name="v">The value (not cloned).</param>
+    /// <param name="gradient">The gradient (not cloned).</param>
+    internal MultivariateDual(V v, AlgebraicSparseVector<object, V> gradient) { this.v = v; this.dv = gradient; }
+
+    public MultivariateDual<V> Clone() { return new MultivariateDual<V>(this); }
+
+    public MultivariateDual<V> Zero => new MultivariateDual<V>(Value.Zero, Gradient.Zero);
+    public MultivariateDual<V> One => new MultivariateDual<V>(Value.One, Gradient.Zero);
+
+    public MultivariateDual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }
+
+    public MultivariateDual<V> Add(MultivariateDual<V> a) { v.Add(a.v); dv.Add(a.dv); return this; }
+    public MultivariateDual<V> Sub(MultivariateDual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }
+    public MultivariateDual<V> Assign(MultivariateDual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }
+    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).Add(t2);
+
+      v.Mul(a.v);
+      return this;
+    }
+    public MultivariateDual<V> Div(MultivariateDual<V> a) { v.Div(a.v); dv.Mul(a.dv.Inv()); return this; }
+    public MultivariateDual<V> AssignNeg(MultivariateDual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }
+    public MultivariateDual<V> AssignInv(MultivariateDual<V> a) { v.AssignInv(a.v); dv.AssignNeg(a.dv).Scale(v).Scale(v); return this; }   // (1/f(x))' = - f(x)' / f(x)^2
+
+    public MultivariateDual<V> AssignSin(MultivariateDual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Cos()); return this; }
+    public MultivariateDual<V> AssignCos(MultivariateDual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Scale(a.v.Clone().Sin()); return this; }
+    public MultivariateDual<V> AssignTanh(MultivariateDual<V> a) { v.AssignTanh(a.v); dv.Assign(a.dv.Scale(v.Clone().IntPower(2).Neg().Add(Value.One))); return this; }     // tanh(f(x))' = f(x)'sech²(f(x)) = f(x)'(1 - tanh²(f(x)))
+
+    public MultivariateDual<V> AssignIntPower(MultivariateDual<V> a, int p) { v.AssignIntPower(a.v, p); dv.Assign(a.dv).Scale(p).Scale(a.v.Clone().IntPower(p - 1)); return this; }
+    public MultivariateDual<V> AssignIntRoot(MultivariateDual<V> a, int r) { v.AssignIntRoot(a.v, r); dv.Assign(a.dv).Scale(1.0 / r).Scale(a.v.IntRoot(r - 1)); return this; }
+
+    public MultivariateDual<V> AssignExp(MultivariateDual<V> a) { v.AssignExp(a.v); dv.Assign(a.dv).Scale(v); return this; } // exp(f(x)) = exp(f(x))*f(x)'     
+    public MultivariateDual<V> AssignLog(MultivariateDual<V> a) { v.AssignLog(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Inv()); return this; }     // log(x)' = 1/f(x) * f(x)' 
+
+    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) {
+      XXX
+    }
+
+    public MultivariateDual<V> AssignMax(MultivariateDual<V> other) {
+      XXX
+    }
+  }
+}
Index: /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/VectorAutoDiffEvaluator.cs
===================================================================
--- /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/VectorAutoDiffEvaluator.cs	(revision 17295)
+++ /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/VectorAutoDiffEvaluator.cs	(revision 17295)
@@ -0,0 +1,145 @@
+﻿using System;
+using System.Collections.Generic;
+using System.Linq;
+using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
+
+namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
+  public sealed class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<AlgebraicDoubleVector>> {
+    private const int BATCHSIZE = 128;
+    [ThreadStatic]
+    private Dictionary<string, double[]> cachedData;
+
+    [ThreadStatic]
+    private IDataset dataset;
+
+    [ThreadStatic]
+    private int rowIndex;
+
+    [ThreadStatic]
+    private int[] rows;
+
+    [ThreadStatic]
+    private Dictionary<ISymbolicExpressionTreeNode, int> node2paramIdx;
+
+    private void InitCache(IDataset dataset) {
+      this.dataset = dataset;
+      cachedData = new Dictionary<string, double[]>();
+      foreach (var v in dataset.DoubleVariables) {
+        cachedData[v] = dataset.GetDoubleValues(v).ToArray();
+      }
+    }
+
+    /// <summary>
+    /// 
+    /// </summary>
+    /// <param name="tree"></param>
+    /// <param name="dataset"></param>
+    /// <param name="rows"></param>
+    /// <param name="parameterNodes"></param>
+    /// <param name="fi">Function output. Must be allocated by the caller.</param>
+    /// <param name="jac">Jacobian matrix. Must be allocated by the caller.</param>
+    public void Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows, ISymbolicExpressionTreeNode[] parameterNodes, double[] fi, double[,] jac) {
+      if (cachedData == null || this.dataset != dataset) {
+        InitCache(dataset);
+      }
+
+      int nParams = parameterNodes.Length;
+      node2paramIdx = new Dictionary<ISymbolicExpressionTreeNode, int>();
+      for (int i = 0; i < parameterNodes.Length; i++) node2paramIdx.Add(parameterNodes[i], i);
+
+      var code = Compile(tree);
+
+      var remainingRows = rows.Length % BATCHSIZE;
+      var roundedTotal = rows.Length - remainingRows;
+
+      this.rows = rows;
+
+      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;
+          }
+        }
+      }
+
+      if (remainingRows > 0) {
+        Evaluate(code);
+        code[0].value.Value.CopyTo(fi, roundedTotal, remainingRows);
+
+        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;
+          }
+      }
+    }
+
+    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
+      var zero = new AlgebraicDoubleVector(BATCHSIZE);
+      instruction.value = new MultivariateDual<AlgebraicDoubleVector>(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 AlgebraicDoubleVector(g_arr);
+        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE), paramIdx, g); // only a single column for the gradient
+      } else {
+        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE));
+      }
+
+      instruction.dblVal = constant.Value;
+      instruction.value.Value.AssignConstant(instruction.dblVal);
+    }
+
+    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
+      double[] data;
+      if (cachedData.ContainsKey(variable.VariableName)) {
+        data = cachedData[variable.VariableName];
+      } else {
+        data = dataset.GetReadOnlyDoubleValues(variable.VariableName).ToArray();
+        cachedData[variable.VariableName] = (double[])instruction.data;
+      }
+
+      var paramIdx = -1;
+      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);
+      } else {
+        var f = new AlgebraicDoubleVector(BATCHSIZE);
+        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(f);
+      }
+
+      instruction.dblVal = variable.Weight;
+      instruction.data = new object[] { data, paramIdx };
+    }
+
+    protected override void LoadVariable(Instruction a) {
+      var paramIdx = (int)((object[])a.data)[1];
+      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);
+
+      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]];
+        }
+      }
+    }
+  }
+}
Index: /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/VectorEvaluator.cs
===================================================================
--- /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/VectorEvaluator.cs	(revision 17295)
+++ /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/VectorEvaluator.cs	(revision 17295)
@@ -0,0 +1,81 @@
+﻿using System;
+using System.Collections.Generic;
+using System.Linq;
+using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
+
+namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
+  public sealed class VectorEvaluator : Interpreter<AlgebraicDoubleVector> {
+    private const int BATCHSIZE = 128;
+    [ThreadStatic]
+    private Dictionary<string, double[]> cachedData;
+
+    [ThreadStatic]
+    private IDataset dataset;
+
+    [ThreadStatic]
+    private int rowIndex;
+
+    [ThreadStatic]
+    private int[] rows;
+
+    private void InitCache(IDataset dataset) {
+      this.dataset = dataset;
+      cachedData = new Dictionary<string, double[]>();
+      foreach (var v in dataset.DoubleVariables) {
+        cachedData[v] = dataset.GetReadOnlyDoubleValues(v).ToArray();
+      }
+    }
+
+    public double[] Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows) {
+      if (cachedData == null || this.dataset != dataset) {
+        InitCache(dataset);
+      }
+
+      this.rows = rows;
+      var code = Compile(tree);
+      var remainingRows = rows.Length % BATCHSIZE;
+      var roundedTotal = rows.Length - remainingRows;
+
+      var result = new double[rows.Length];
+
+      for (rowIndex = 0; rowIndex < roundedTotal; rowIndex += BATCHSIZE) {
+        Evaluate(code);
+        code[0].value.CopyTo(result, rowIndex, BATCHSIZE);
+      }
+
+      if (remainingRows > 0) {
+        Evaluate(code);
+        code[0].value.CopyTo(result, roundedTotal, remainingRows);
+      }
+
+      return result;
+    }
+
+    protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
+      instruction.dblVal = constant.Value;
+      instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
+      instruction.value.AssignConstant(instruction.dblVal);
+    }
+
+    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
+      instruction.dblVal = variable.Weight;
+      instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
+      if (cachedData.ContainsKey(variable.VariableName)) {
+        instruction.data = cachedData[variable.VariableName];
+      } else {
+        instruction.data = dataset.GetDoubleValues(variable.VariableName).ToArray();
+        cachedData[variable.VariableName] = (double[])instruction.data;
+      }
+    }
+
+    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
+      instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
+    }
+
+    protected override void LoadVariable(Instruction a) {
+      var data = (double[])a.data;
+      for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) a.value[i - rowIndex] = data[rows[i]];
+      a.value.Scale(a.dblVal);
+    }
+  }
+}
