Index: branches/3040_VectorBasedGP/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/SymbolicDataAnalysisExpressionTreeVectorInterpreter.cs
===================================================================
--- branches/3040_VectorBasedGP/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/SymbolicDataAnalysisExpressionTreeInterpreter.cs	(revision 17448)
+++ branches/3040_VectorBasedGP/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/SymbolicDataAnalysisExpressionTreeVectorInterpreter.cs	(revision 17455)
@@ -20,8 +20,7 @@
 #endregion
 
-using DoubleVector = MathNet.Numerics.LinearAlgebra.Vector<double>;
-
 using System;
 using System.Collections.Generic;
+using HeuristicLab.Analysis;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -30,13 +29,14 @@
 using HeuristicLab.Parameters;
 using HEAL.Attic;
+using MathNet.Numerics.LinearAlgebra;
 using MathNet.Numerics.Statistics;
 
+using DoubleVector = MathNet.Numerics.LinearAlgebra.Vector<double>;
+
 namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
-  [StorableType("FB94F333-B32A-44FB-A561-CBDE76693D20")]
-  [Item("SymbolicDataAnalysisExpressionTreeInterpreter", "Interpreter for symbolic expression trees including automatically defined functions.")]
-  public class SymbolicDataAnalysisExpressionTreeInterpreter : ParameterizedNamedItem,
-    ISymbolicDataAnalysisExpressionTreeInterpreter {
-    private const string CheckExpressionsWithIntervalArithmeticParameterName = "CheckExpressionsWithIntervalArithmetic";
-    private const string CheckExpressionsWithIntervalArithmeticParameterDescription = "Switch that determines if the interpreter checks the validity of expressions with interval arithmetic before evaluating the expression.";
+  [StorableType("DE68A1D9-5AFC-4DDD-AB62-29F3B8FC28E0")]
+  [Item("SymbolicDataAnalysisExpressionTreeVectorInterpreter", "Interpreter for symbolic expression trees including vector arithmetic.")]
+  public class SymbolicDataAnalysisExpressionTreeVectorInterpreter : ParameterizedNamedItem, ISymbolicDataAnalysisExpressionTreeInterpreter {
+
     private const string EvaluatedSolutionsParameterName = "EvaluatedSolutions";
 
@@ -50,8 +50,4 @@
 
     #region parameter properties
-    public IFixedValueParameter<BoolValue> CheckExpressionsWithIntervalArithmeticParameter {
-      get { return (IFixedValueParameter<BoolValue>)Parameters[CheckExpressionsWithIntervalArithmeticParameterName]; }
-    }
-
     public IFixedValueParameter<IntValue> EvaluatedSolutionsParameter {
       get { return (IFixedValueParameter<IntValue>)Parameters[EvaluatedSolutionsParameterName]; }
@@ -60,9 +56,4 @@
 
     #region properties
-    public bool CheckExpressionsWithIntervalArithmetic {
-      get { return CheckExpressionsWithIntervalArithmeticParameter.Value.Value; }
-      set { CheckExpressionsWithIntervalArithmeticParameter.Value.Value = value; }
-    }
-
     public int EvaluatedSolutions {
       get { return EvaluatedSolutionsParameter.Value.Value; }
@@ -72,23 +63,20 @@
 
     [StorableConstructor]
-    protected SymbolicDataAnalysisExpressionTreeInterpreter(StorableConstructorFlag _) : base(_) { }
-
-    protected SymbolicDataAnalysisExpressionTreeInterpreter(SymbolicDataAnalysisExpressionTreeInterpreter original,
-      Cloner cloner)
+    protected SymbolicDataAnalysisExpressionTreeVectorInterpreter(StorableConstructorFlag _) : base(_) { }
+
+    protected SymbolicDataAnalysisExpressionTreeVectorInterpreter(SymbolicDataAnalysisExpressionTreeVectorInterpreter original, Cloner cloner)
       : base(original, cloner) { }
 
     public override IDeepCloneable Clone(Cloner cloner) {
-      return new SymbolicDataAnalysisExpressionTreeInterpreter(this, cloner);
-    }
-
-    public SymbolicDataAnalysisExpressionTreeInterpreter()
-      : base("SymbolicDataAnalysisExpressionTreeInterpreter", "Interpreter for symbolic expression trees including automatically defined functions.") {
-      Parameters.Add(new FixedValueParameter<BoolValue>(CheckExpressionsWithIntervalArithmeticParameterName, "Switch that determines if the interpreter checks the validity of expressions with interval arithmetic before evaluating the expression.", new BoolValue(false)));
+      return new SymbolicDataAnalysisExpressionTreeVectorInterpreter(this, cloner);
+    }
+
+    public SymbolicDataAnalysisExpressionTreeVectorInterpreter()
+      : base("SymbolicDataAnalysisExpressionTreeVectorInterpreter", "Interpreter for symbolic expression trees including vector arithmetic.") {
       Parameters.Add(new FixedValueParameter<IntValue>(EvaluatedSolutionsParameterName, "A counter for the total number of solutions the interpreter has evaluated", new IntValue(0)));
     }
 
-    protected SymbolicDataAnalysisExpressionTreeInterpreter(string name, string description)
+    protected SymbolicDataAnalysisExpressionTreeVectorInterpreter(string name, string description)
       : base(name, description) {
-      Parameters.Add(new FixedValueParameter<BoolValue>(CheckExpressionsWithIntervalArithmeticParameterName, "Switch that determines if the interpreter checks the validity of expressions with interval arithmetic before evaluating the expression.", new BoolValue(false)));
       Parameters.Add(new FixedValueParameter<IntValue>(EvaluatedSolutionsParameterName, "A counter for the total number of solutions the interpreter has evaluated", new IntValue(0)));
     }
@@ -96,18 +84,5 @@
     [StorableHook(HookType.AfterDeserialization)]
     private void AfterDeserialization() {
-      var evaluatedSolutions = new IntValue(0);
-      var checkExpressionsWithIntervalArithmetic = new BoolValue(false);
-      if (Parameters.ContainsKey(EvaluatedSolutionsParameterName)) {
-        var evaluatedSolutionsParameter = (IValueParameter<IntValue>)Parameters[EvaluatedSolutionsParameterName];
-        evaluatedSolutions = evaluatedSolutionsParameter.Value;
-        Parameters.Remove(EvaluatedSolutionsParameterName);
-      }
-      Parameters.Add(new FixedValueParameter<IntValue>(EvaluatedSolutionsParameterName, "A counter for the total number of solutions the interpreter has evaluated", evaluatedSolutions));
-      if (Parameters.ContainsKey(CheckExpressionsWithIntervalArithmeticParameterName)) {
-        var checkExpressionsWithIntervalArithmeticParameter = (IValueParameter<BoolValue>)Parameters[CheckExpressionsWithIntervalArithmeticParameterName];
-        Parameters.Remove(CheckExpressionsWithIntervalArithmeticParameterName);
-        checkExpressionsWithIntervalArithmetic = checkExpressionsWithIntervalArithmeticParameter.Value;
-      }
-      Parameters.Add(new FixedValueParameter<BoolValue>(CheckExpressionsWithIntervalArithmeticParameterName, CheckExpressionsWithIntervalArithmeticParameterDescription, checkExpressionsWithIntervalArithmetic));
+
     }
 
@@ -121,10 +96,5 @@
 
     private readonly object syncRoot = new object();
-    public IEnumerable<double> GetSymbolicExpressionTreeValues(ISymbolicExpressionTree tree, IDataset dataset,
-      IEnumerable<int> rows) {
-      if (CheckExpressionsWithIntervalArithmetic) {
-        throw new NotSupportedException("Interval arithmetic is not yet supported in the symbolic data analysis interpreter.");
-      }
-
+    public IEnumerable<double> GetSymbolicExpressionTreeValues(ISymbolicExpressionTree tree, IDataset dataset, IEnumerable<int> rows) {
       lock (syncRoot) {
         EvaluatedSolutions++; // increment the evaluated solutions counter
@@ -134,5 +104,8 @@
       foreach (var rowEnum in rows) {
         int row = rowEnum;
-        yield return Evaluate(dataset, ref row, state);
+        var result = Evaluate(dataset, ref row, state);
+        if (!result.IsScalar)
+          throw new InvalidOperationException("Result of the tree is not a scalar.");
+        yield return result.Scalar;
         state.Reset();
       }
@@ -152,7 +125,4 @@
           var factorTreeNode = instr.dynamicNode as BinaryFactorVariableTreeNode;
           instr.data = dataset.GetReadOnlyStringValues(factorTreeNode.VariableName);
-        } else if (instr.opCode == OpCodes.VectorVariable) {
-          var vectorVariableTreeNode = (VectorVariableTreeNode)instr.dynamicNode;
-          instr.data = dataset.GetReadOnlyDoubleVectorValues(vectorVariableTreeNode.VariableName);
         } else if (instr.opCode == OpCodes.LagVariable) {
           var laggedVariableTreeNode = (LaggedVariableTreeNode)instr.dynamicNode;
@@ -168,424 +138,204 @@
     }
 
-    public virtual double Evaluate(IDataset dataset, ref int row, InterpreterState state) {
+
+    public struct EvaluationResult {
+      public double Scalar { get; }
+      public bool IsScalar => !double.IsNaN(Scalar);
+
+      public DoubleVector Vector { get; }
+      public bool IsVector => !(Vector.Count == 1 && double.IsNaN(Vector[0]));
+
+      public bool IsNaN => !IsScalar && !IsVector;
+
+      public EvaluationResult(double scalar) {
+        Scalar = scalar;
+        Vector = NaNVector;
+      }
+      public EvaluationResult(DoubleVector vector) {
+        Vector = vector;
+        Scalar = double.NaN;
+      }
+
+      public override string ToString() {
+        if (IsScalar) return Scalar.ToString();
+        if (IsVector) return Vector.ToVectorString();
+        return "NaN";
+      }
+
+      public static readonly EvaluationResult NaN = new EvaluationResult(double.NaN);
+      private static readonly DoubleVector NaNVector = DoubleVector.Build.Dense(1, double.NaN);
+    }
+
+    private static EvaluationResult ArithmeticApply(EvaluationResult lhs, EvaluationResult rhs,
+      Func<double, double, double> ssFunc = null,
+      Func<double, DoubleVector, DoubleVector> svFunc = null,
+      Func<DoubleVector, double, DoubleVector> vsFunc = null,
+      Func<DoubleVector, DoubleVector, DoubleVector> vvFunc = null) {
+      if (lhs.IsScalar && rhs.IsScalar && ssFunc != null) return new EvaluationResult(ssFunc(lhs.Scalar, rhs.Scalar));
+      if (lhs.IsScalar && rhs.IsVector && svFunc != null) return new EvaluationResult(svFunc(lhs.Scalar, rhs.Vector));
+      if (lhs.IsVector && rhs.IsScalar && vsFunc != null) return new EvaluationResult(vsFunc(lhs.Vector, rhs.Scalar));
+      if (lhs.IsVector && rhs.IsVector && vvFunc != null) return new EvaluationResult(vvFunc(lhs.Vector, rhs.Vector));
+      throw new NotSupportedException($"Unsupported combination of argument types: ({lhs}) / ({rhs})");
+    }
+
+    private static EvaluationResult FunctionApply(EvaluationResult val,
+      Func<double, double> sFunc = null,
+      Func<DoubleVector, DoubleVector> vFunc = null) {
+      if (val.IsScalar && sFunc != null) return new EvaluationResult(sFunc(val.Scalar));
+      if (val.IsVector && vFunc != null) return new EvaluationResult(vFunc(val.Vector));
+      throw new NotSupportedException($"Unsupported argument type ({val})");
+    }
+
+    public virtual EvaluationResult Evaluate(IDataset dataset, ref int row, InterpreterState state) {
       Instruction currentInstr = state.NextInstruction();
       switch (currentInstr.opCode) {
         case OpCodes.Add: {
-            double s = Evaluate(dataset, ref row, state);
+            var cur = Evaluate(dataset, ref row, state);
             for (int i = 1; i < currentInstr.nArguments; i++) {
-              s += Evaluate(dataset, ref row, state);
+              var op = Evaluate(dataset, ref row, state);
+              cur = ArithmeticApply(cur, op,
+                (s1, s2) => s1 + s2,
+                (s1, v2) => s1 + v2,
+                (v1, s2) => v1 + s2,
+                (v1, v2) => v1 + v2);
             }
-            return s;
+            return cur;
           }
         case OpCodes.Sub: {
-            double s = Evaluate(dataset, ref row, state);
+            var cur = Evaluate(dataset, ref row, state);
             for (int i = 1; i < currentInstr.nArguments; i++) {
-              s -= Evaluate(dataset, ref row, state);
+              var op = Evaluate(dataset, ref row, state);
+              cur = ArithmeticApply(cur, op,
+                (s1, s2) => s1 - s2,
+                (s1, v2) => s1 - v2,
+                (v1, s2) => v1 - s2,
+                (v1, v2) => v1 - v2);
             }
-            if (currentInstr.nArguments == 1) { s = -s; }
-            return s;
+            return cur;
           }
         case OpCodes.Mul: {
-            double p = Evaluate(dataset, ref row, state);
+            var cur = Evaluate(dataset, ref row, state);
             for (int i = 1; i < currentInstr.nArguments; i++) {
-              p *= Evaluate(dataset, ref row, state);
+              var op = Evaluate(dataset, ref row, state);
+              cur = ArithmeticApply(cur, op,
+                (s1, s2) => s1 * s2,
+                (s1, v2) => s1 * v2,
+                (v1, s2) => v1 * s2,
+                (v1, v2) => v1.PointwiseMultiply(v2));
             }
-            return p;
+            return cur;
           }
         case OpCodes.Div: {
-            double p = Evaluate(dataset, ref row, state);
+            var cur = Evaluate(dataset, ref row, state);
             for (int i = 1; i < currentInstr.nArguments; i++) {
-              p /= Evaluate(dataset, ref row, state);
+              var op = Evaluate(dataset, ref row, state);
+              cur = ArithmeticApply(cur, op,
+                (s1, s2) => s1 / s2,
+                (s1, v2) => s1 / v2,
+                (v1, s2) => v1 / s2,
+                (v1, v2) => v1 / v2);
             }
-            if (currentInstr.nArguments == 1) { p = 1.0 / p; }
-            return p;
-          }
-        case OpCodes.Average: {
-            double sum = Evaluate(dataset, ref row, state);
-            for (int i = 1; i < currentInstr.nArguments; i++) {
-              sum += Evaluate(dataset, ref row, state);
-            }
-            return sum / currentInstr.nArguments;
+            return cur;
           }
         case OpCodes.Absolute: {
-            return Math.Abs(Evaluate(dataset, ref row, state));
+            var cur = Evaluate(dataset, ref row, state);
+            return FunctionApply(cur, Math.Abs, DoubleVector.Abs);
           }
         case OpCodes.Tanh: {
-            return Math.Tanh(Evaluate(dataset, ref row, state));
+            var cur = Evaluate(dataset, ref row, state);
+            return FunctionApply(cur, Math.Tanh, DoubleVector.Tanh);
           }
         case OpCodes.Cos: {
-            return Math.Cos(Evaluate(dataset, ref row, state));
+            var cur = Evaluate(dataset, ref row, state);
+            return FunctionApply(cur, Math.Cos, DoubleVector.Cos);
           }
         case OpCodes.Sin: {
-            return Math.Sin(Evaluate(dataset, ref row, state));
+            var cur = Evaluate(dataset, ref row, state);
+            return FunctionApply(cur, Math.Sin, DoubleVector.Sin);
           }
         case OpCodes.Tan: {
-            return Math.Tan(Evaluate(dataset, ref row, state));
+            var cur = Evaluate(dataset, ref row, state);
+            return FunctionApply(cur, Math.Tan, DoubleVector.Tan);
           }
         case OpCodes.Square: {
-            return Math.Pow(Evaluate(dataset, ref row, state), 2);
+            var cur = Evaluate(dataset, ref row, state);
+            return FunctionApply(cur,
+              s => Math.Pow(s, 2),
+              v => v.PointwisePower(2));
           }
         case OpCodes.Cube: {
-            return Math.Pow(Evaluate(dataset, ref row, state), 3);
+            var cur = Evaluate(dataset, ref row, state);
+            return FunctionApply(cur,
+              s => Math.Pow(s, 3),
+              v => v.PointwisePower(3));
           }
         case OpCodes.Power: {
-            double x = Evaluate(dataset, ref row, state);
-            double y = Math.Round(Evaluate(dataset, ref row, state));
-            return Math.Pow(x, y);
+            var x = Evaluate(dataset, ref row, state);
+            var y = Evaluate(dataset, ref row, state);
+            return ArithmeticApply(x, y,
+              (s1, s2) => Math.Pow(s1, Math.Round(s2)),
+              (s1, v2) => DoubleVector.Build.Dense(v2.Count, s1).PointwisePower(DoubleVector.Round(v2)),
+              (v1, s2) => v1.PointwisePower(Math.Round(s2)),
+              (v1, v2) => v1.PointwisePower(DoubleVector.Round(v2)));
           }
         case OpCodes.SquareRoot: {
-            return Math.Sqrt(Evaluate(dataset, ref row, state));
+            var cur = Evaluate(dataset, ref row, state);
+            return FunctionApply(cur,
+              s => Math.Sqrt(s),
+              v => DoubleVector.Sqrt(v));
           }
         case OpCodes.CubeRoot: {
-            var arg = Evaluate(dataset, ref row, state);
-            return arg < 0 ? -Math.Pow(-arg, 1.0 / 3.0) : Math.Pow(arg, 1.0 / 3.0);
+            var cur = Evaluate(dataset, ref row, state);
+            return FunctionApply(cur,
+              s => s < 0 ? -Math.Pow(-s, 1.0 / 3.0) : Math.Pow(s, 1.0 / 3.0),
+              v => v.Map(s => s < 0 ? -Math.Pow(-s, 1.0 / 3.0) : Math.Pow(s, 1.0 / 3.0)));
           }
         case OpCodes.Root: {
-            double x = Evaluate(dataset, ref row, state);
-            double y = Math.Round(Evaluate(dataset, ref row, state));
-            return Math.Pow(x, 1 / y);
+            var x = Evaluate(dataset, ref row, state);
+            var y = Evaluate(dataset, ref row, state);
+            return ArithmeticApply(x, y,
+              (s1, s2) => Math.Pow(s1, 1.0 / Math.Round(s2)),
+              (s1, v2) => DoubleVector.Build.Dense(v2.Count, s1).PointwisePower(1.0 / DoubleVector.Round(v2)),
+              (v1, s2) => v1.PointwisePower(1.0 / Math.Round(s2)),
+              (v1, v2) => v1.PointwisePower(1.0 / DoubleVector.Round(v2)));
           }
         case OpCodes.Exp: {
-            return Math.Exp(Evaluate(dataset, ref row, state));
+            var cur = Evaluate(dataset, ref row, state);
+            return FunctionApply(cur,
+              s => Math.Exp(s),
+              v => DoubleVector.Exp(v));
           }
         case OpCodes.Log: {
-            return Math.Log(Evaluate(dataset, ref row, state));
-          }
-        case OpCodes.Gamma: {
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) { return double.NaN; } else { return alglib.gammafunction(x); }
-          }
-        case OpCodes.Psi: {
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) return double.NaN;
-            else if (x <= 0 && (Math.Floor(x) - x).IsAlmost(0)) return double.NaN;
-            return alglib.psi(x);
-          }
-        case OpCodes.Dawson: {
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) { return double.NaN; }
-            return alglib.dawsonintegral(x);
-          }
-        case OpCodes.ExponentialIntegralEi: {
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) { return double.NaN; }
-            return alglib.exponentialintegralei(x);
-          }
-        case OpCodes.SineIntegral: {
-            double si, ci;
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) return double.NaN;
-            else {
-              alglib.sinecosineintegrals(x, out si, out ci);
-              return si;
-            }
-          }
-        case OpCodes.CosineIntegral: {
-            double si, ci;
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) return double.NaN;
-            else {
-              alglib.sinecosineintegrals(x, out si, out ci);
-              return ci;
-            }
-          }
-        case OpCodes.HyperbolicSineIntegral: {
-            double shi, chi;
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) return double.NaN;
-            else {
-              alglib.hyperbolicsinecosineintegrals(x, out shi, out chi);
-              return shi;
-            }
-          }
-        case OpCodes.HyperbolicCosineIntegral: {
-            double shi, chi;
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) return double.NaN;
-            else {
-              alglib.hyperbolicsinecosineintegrals(x, out shi, out chi);
-              return chi;
-            }
-          }
-        case OpCodes.FresnelCosineIntegral: {
-            double c = 0, s = 0;
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) return double.NaN;
-            else {
-              alglib.fresnelintegral(x, ref c, ref s);
-              return c;
-            }
-          }
-        case OpCodes.FresnelSineIntegral: {
-            double c = 0, s = 0;
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) return double.NaN;
-            else {
-              alglib.fresnelintegral(x, ref c, ref s);
-              return s;
-            }
-          }
-        case OpCodes.AiryA: {
-            double ai, aip, bi, bip;
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) return double.NaN;
-            else {
-              alglib.airy(x, out ai, out aip, out bi, out bip);
-              return ai;
-            }
-          }
-        case OpCodes.AiryB: {
-            double ai, aip, bi, bip;
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) return double.NaN;
-            else {
-              alglib.airy(x, out ai, out aip, out bi, out bip);
-              return bi;
-            }
-          }
-        case OpCodes.Norm: {
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) return double.NaN;
-            else return alglib.normaldistribution(x);
-          }
-        case OpCodes.Erf: {
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) return double.NaN;
-            else return alglib.errorfunction(x);
-          }
-        case OpCodes.Bessel: {
-            var x = Evaluate(dataset, ref row, state);
-            if (double.IsNaN(x)) return double.NaN;
-            else return alglib.besseli0(x);
-          }
-
-        case OpCodes.AnalyticQuotient: {
-            var x1 = Evaluate(dataset, ref row, state);
-            var x2 = Evaluate(dataset, ref row, state);
-            return x1 / Math.Pow(1 + x2 * x2, 0.5);
-          }
-        case OpCodes.IfThenElse: {
-            double condition = Evaluate(dataset, ref row, state);
-            double result;
-            if (condition > 0.0) {
-              result = Evaluate(dataset, ref row, state); state.SkipInstructions();
-            } else {
-              state.SkipInstructions(); result = Evaluate(dataset, ref row, state);
-            }
-            return result;
-          }
-        case OpCodes.AND: {
-            double result = Evaluate(dataset, ref row, state);
-            for (int i = 1; i < currentInstr.nArguments; i++) {
-              if (result > 0.0) result = Evaluate(dataset, ref row, state);
-              else {
-                state.SkipInstructions();
-              }
-            }
-            return result > 0.0 ? 1.0 : -1.0;
-          }
-        case OpCodes.OR: {
-            double result = Evaluate(dataset, ref row, state);
-            for (int i = 1; i < currentInstr.nArguments; i++) {
-              if (result <= 0.0) result = Evaluate(dataset, ref row, state);
-              else {
-                state.SkipInstructions();
-              }
-            }
-            return result > 0.0 ? 1.0 : -1.0;
-          }
-        case OpCodes.NOT: {
-            return Evaluate(dataset, ref row, state) > 0.0 ? -1.0 : 1.0;
-          }
-        case OpCodes.XOR: {
-            //mkommend: XOR on multiple inputs is defined as true if the number of positive signals is odd
-            // this is equal to a consecutive execution of binary XOR operations.
-            int positiveSignals = 0;
-            for (int i = 0; i < currentInstr.nArguments; i++) {
-              if (Evaluate(dataset, ref row, state) > 0.0) { positiveSignals++; }
-            }
-            return positiveSignals % 2 != 0 ? 1.0 : -1.0;
-          }
-        case OpCodes.GT: {
-            double x = Evaluate(dataset, ref row, state);
-            double y = Evaluate(dataset, ref row, state);
-            if (x > y) { return 1.0; } else { return -1.0; }
-          }
-        case OpCodes.LT: {
-            double x = Evaluate(dataset, ref row, state);
-            double y = Evaluate(dataset, ref row, state);
-            if (x < y) { return 1.0; } else { return -1.0; }
-          }
-        case OpCodes.TimeLag: {
-            var timeLagTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
-            row += timeLagTreeNode.Lag;
-            double result = Evaluate(dataset, ref row, state);
-            row -= timeLagTreeNode.Lag;
-            return result;
-          }
-        case OpCodes.Integral: {
-            int savedPc = state.ProgramCounter;
-            var timeLagTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
-            double sum = 0.0;
-            for (int i = 0; i < Math.Abs(timeLagTreeNode.Lag); i++) {
-              row += Math.Sign(timeLagTreeNode.Lag);
-              sum += Evaluate(dataset, ref row, state);
-              state.ProgramCounter = savedPc;
-            }
-            row -= timeLagTreeNode.Lag;
-            sum += Evaluate(dataset, ref row, state);
-            return sum;
-          }
-
-        //mkommend: derivate calculation taken from: 
-        //http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
-        //one sided smooth differentiatior, N = 4
-        // y' = 1/8h (f_i + 2f_i-1, -2 f_i-3 - f_i-4)
-        case OpCodes.Derivative: {
-            int savedPc = state.ProgramCounter;
-            double f_0 = Evaluate(dataset, ref row, state); row--;
-            state.ProgramCounter = savedPc;
-            double f_1 = Evaluate(dataset, ref row, state); row -= 2;
-            state.ProgramCounter = savedPc;
-            double f_3 = Evaluate(dataset, ref row, state); row--;
-            state.ProgramCounter = savedPc;
-            double f_4 = Evaluate(dataset, ref row, state);
-            row += 4;
-
-            return (f_0 + 2 * f_1 - 2 * f_3 - f_4) / 8; // h = 1
-          }
-        case OpCodes.Call: {
-            // evaluate sub-trees
-            double[] argValues = new double[currentInstr.nArguments];
-            for (int i = 0; i < currentInstr.nArguments; i++) {
-              argValues[i] = Evaluate(dataset, ref row, state);
-            }
-            // push on argument values on stack 
-            state.CreateStackFrame(argValues);
-
-            // save the pc
-            int savedPc = state.ProgramCounter;
-            // set pc to start of function  
-            state.ProgramCounter = (ushort)currentInstr.data;
-            // evaluate the function
-            double v = Evaluate(dataset, ref row, state);
-
-            // delete the stack frame
-            state.RemoveStackFrame();
-
-            // restore the pc => evaluation will continue at point after my subtrees  
-            state.ProgramCounter = savedPc;
-            return v;
-          }
-        case OpCodes.Arg: {
-            return state.GetStackFrameValue((ushort)currentInstr.data);
+            var cur = Evaluate(dataset, ref row, state);
+            return FunctionApply(cur,
+              s => Math.Log(s),
+              v => DoubleVector.Log(v));
           }
         case OpCodes.Variable: {
-            if (row < 0 || row >= dataset.Rows) return double.NaN;
+            if (row < 0 || row >= dataset.Rows) return EvaluationResult.NaN;
             var variableTreeNode = (VariableTreeNode)currentInstr.dynamicNode;
-            return ((IList<double>)currentInstr.data)[row] * variableTreeNode.Weight;
+            if (currentInstr.data is IList<double> doubleList)
+              return new EvaluationResult(doubleList[row] * variableTreeNode.Weight);
+            if (currentInstr.data is IList<DoubleVector> doubleVectorList)
+              return new EvaluationResult(doubleVectorList[row] * variableTreeNode.Weight);
+            throw new NotSupportedException($"Unsupported type of variable: {currentInstr.data.GetType().GetPrettyName()}");
           }
         case OpCodes.BinaryFactorVariable: {
-            if (row < 0 || row >= dataset.Rows) return double.NaN;
+            if (row < 0 || row >= dataset.Rows) return EvaluationResult.NaN;
             var factorVarTreeNode = currentInstr.dynamicNode as BinaryFactorVariableTreeNode;
-            return ((IList<string>)currentInstr.data)[row] == factorVarTreeNode.VariableValue ? factorVarTreeNode.Weight : 0;
+            return new EvaluationResult(((IList<string>)currentInstr.data)[row] == factorVarTreeNode.VariableValue ? factorVarTreeNode.Weight : 0);
           }
         case OpCodes.FactorVariable: {
-            if (row < 0 || row >= dataset.Rows) return double.NaN;
+            if (row < 0 || row >= dataset.Rows) return EvaluationResult.NaN;
             var factorVarTreeNode = currentInstr.dynamicNode as FactorVariableTreeNode;
-            return factorVarTreeNode.GetValue(((IList<string>)currentInstr.data)[row]);
-          }
-        case OpCodes.LagVariable: {
-            var laggedVariableTreeNode = (LaggedVariableTreeNode)currentInstr.dynamicNode;
-            int actualRow = row + laggedVariableTreeNode.Lag;
-            if (actualRow < 0 || actualRow >= dataset.Rows) { return double.NaN; }
-            return ((IList<double>)currentInstr.data)[actualRow] * laggedVariableTreeNode.Weight;
+            return new EvaluationResult(factorVarTreeNode.GetValue(((IList<string>)currentInstr.data)[row]));
           }
         case OpCodes.Constant: {
             var constTreeNode = (ConstantTreeNode)currentInstr.dynamicNode;
-            return constTreeNode.Value;
-          }
-
-        //mkommend: this symbol uses the logistic function f(x) = 1 / (1 + e^(-alpha * x) ) 
-        //to determine the relative amounts of the true and false branch see http://en.wikipedia.org/wiki/Logistic_function
-        case OpCodes.VariableCondition: {
-            if (row < 0 || row >= dataset.Rows) return double.NaN;
-            var variableConditionTreeNode = (VariableConditionTreeNode)currentInstr.dynamicNode;
-            if (!variableConditionTreeNode.Symbol.IgnoreSlope) {
-              double variableValue = ((IList<double>)currentInstr.data)[row];
-              double x = variableValue - variableConditionTreeNode.Threshold;
-              double p = 1 / (1 + Math.Exp(-variableConditionTreeNode.Slope * x));
-
-              double trueBranch = Evaluate(dataset, ref row, state);
-              double falseBranch = Evaluate(dataset, ref row, state);
-
-              return trueBranch * p + falseBranch * (1 - p);
-            } else {
-              // strict threshold
-              double variableValue = ((IList<double>)currentInstr.data)[row];
-              if (variableValue <= variableConditionTreeNode.Threshold) {
-                var left = Evaluate(dataset, ref row, state);
-                state.SkipInstructions();
-                return left;
-              } else {
-                state.SkipInstructions();
-                return Evaluate(dataset, ref row, state);
-              }
-            }
-          }
-
-        case OpCodes.VectorSum: {
-            DoubleVector v = VectorEvaluate(dataset, ref row, state);
-            return v.Sum();
-          }
-        case OpCodes.VectorMean: {
-            DoubleVector v = VectorEvaluate(dataset, ref row, state);
-            return v.Mean();
+            return new EvaluationResult(constTreeNode.Value);
           }
 
         default:
-          throw new NotSupportedException();
-      }
-    }
-
-    public virtual DoubleVector VectorEvaluate(IDataset dataset, ref int row, InterpreterState state) {
-      Instruction currentInstr = state.NextInstruction();
-      switch (currentInstr.opCode) {
-        case OpCodes.VectorAdd: {
-            DoubleVector s = VectorEvaluate(dataset, ref row, state);
-            for (int i = 1; i < currentInstr.nArguments; i++) {
-              s += VectorEvaluate(dataset, ref row, state);
-            }
-            return s;
-          }
-        case OpCodes.VectorSub: {
-            DoubleVector s = VectorEvaluate(dataset, ref row, state);
-            for (int i = 1; i < currentInstr.nArguments; i++) {
-              s -= VectorEvaluate(dataset, ref row, state);
-            }
-            return s;
-          }
-        case OpCodes.VectorMul: {
-            DoubleVector s = VectorEvaluate(dataset, ref row, state);
-            for (int i = 1; i < currentInstr.nArguments; i++) {
-              s = s.PointwiseMultiply(VectorEvaluate(dataset, ref row, state));
-            }
-            return s;
-          }
-        case OpCodes.VectorDiv: {
-            DoubleVector s = VectorEvaluate(dataset, ref row, state);
-            for (int i = 1; i < currentInstr.nArguments; i++) {
-              s /= VectorEvaluate(dataset, ref row, state);
-            }
-            return s;
-          }
-
-        case OpCodes.VectorVariable: {
-            if (row < 0 || row >= dataset.Rows) return DoubleVector.Build.Dense(new[] { double.NaN });
-            var vectorVarTreeNode = currentInstr.dynamicNode as VectorVariableTreeNode;
-            return ((IList<DoubleVector>)currentInstr.data)[row] * vectorVarTreeNode.Weight;
-          }
-        default:
-          throw new NotSupportedException();
+          throw new NotSupportedException($"Unsupported OpCode: {currentInstr.opCode}");
       }
     }
