Changeset 15313

Timestamp:

08/08/17 19:51:31 (8 years ago)

Author:

gkronber

Message:

#2789: worked on nonlinear regression with constraints

Location:

branches/MathNetNumerics-Exploration-2789

Files:

• HeuristicLab.Algorithms.DataAnalysis.Experimental/DiffSharp.dll (added)
• HeuristicLab.Algorithms.DataAnalysis.Experimental/FSharp.Quotations.Evaluator.dll (added)
• HeuristicLab.Algorithms.DataAnalysis.Experimental/HeuristicLab.Algorithms.DataAnalysis.Experimental.csproj (modified) (5 diffs)
• HeuristicLab.Algorithms.DataAnalysis.Experimental/SymbolicRegressionConstantOptimizationEvaluator.cs (modified) (10 diffs)
• HeuristicLab.Algorithms.DataAnalysis.Experimental/TreeToDiffSharpConverter.cs (moved) (moved from branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/TreeToAutoDiffTermConverter.cs) (7 diffs)
• HeuristicLab.ExtLibs.MathNetNumerics/Plugin.cs.frame (modified) (1 diff)
• MathNetNumerics-Exploration.sln (modified) (2 diffs)
• Test (added)
• Test/Properties (added)
• Test/Properties/AssemblyInfo.cs (added)
• Test/Test.csproj (added)
• Test/UnitTest1.cs (added)

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/HeuristicLab.Algorithms.DataAnalysis.Experimental.csproj

@@ -10 +10 @@
     <RootNamespace>HeuristicLab.Algorithms.DataAnalysis.Experimental</RootNamespace>
     <AssemblyName>HeuristicLab.Algorithms.DataAnalysis.Experimental</AssemblyName>
-    <TargetFrameworkVersion>v4.5.1</TargetFrameworkVersion>
+    <TargetFrameworkVersion>v4.6</TargetFrameworkVersion>
     <FileAlignment>512</FileAlignment>
     <TargetFrameworkProfile />
+    <NuGetPackageImportStamp>
+    </NuGetPackageImportStamp>
   </PropertyGroup>
   <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Debug|AnyCPU' ">
@@ -39 +41 @@
       <SpecificVersion>False</SpecificVersion>
       <HintPath>..\..\..\trunk\sources\bin\AutoDiff-1.0.dll</HintPath>
+    </Reference>
+    <Reference Include="DiffSharp, Version=1.0.0.0, Culture=neutral, processorArchitecture=MSIL">
+      <SpecificVersion>False</SpecificVersion>
+      <HintPath>.\DiffSharp.dll</HintPath>
     </Reference>
     <Reference Include="HeuristicLab.Algorithms.DataAnalysis-3.4">
@@ -101 +107 @@
     </Reference>
     <Reference Include="HeuristicLab.Problems.Instances-3.3, Version=3.3.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec" />
+    <Reference Include="HeuristicLab.Random-3.3, Version=3.3.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
+      <SpecificVersion>False</SpecificVersion>
+      <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Random-3.3.dll</HintPath>
+    </Reference>
     <Reference Include="System" />
     <Reference Include="System.Core" />
@@ -117 +127 @@
     <Compile Include="Properties\AssemblyInfo.cs" />
     <Compile Include="SymbolicRegressionConstantOptimizationEvaluator.cs" />
-    <Compile Include="TreeToAutoDiffTermConverter.cs">
-      <SubType>Code</SubType>
-    </Compile>
+    <Compile Include="TreeToDiffSharpConverter.cs" />
   </ItemGroup>
   <ItemGroup>
@@ -126 +134 @@
   <ItemGroup>
     <Content Include="alglibnet2.dll">
+      <CopyToOutputDirectory>PreserveNewest</CopyToOutputDirectory>
+    </Content>
+    <Content Include="DiffSharp.dll">
+      <CopyToOutputDirectory>PreserveNewest</CopyToOutputDirectory>
+    </Content>
+    <Content Include="FSharp.Quotations.Evaluator.dll">
       <CopyToOutputDirectory>PreserveNewest</CopyToOutputDirectory>
     </Content>

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/SymbolicRegressionConstantOptimizationEvaluator.cs

@@ -23 +23 @@
 using System.Collections.Generic;
 using System.Linq;
-using System.Runtime.Remoting.Channels;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -33 +32 @@
 using HeuristicLab.Problems.DataAnalysis.Symbolic;
 using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;
+using DiffSharp.Interop.Float64;
 
 namespace HeuristicLab.Algorithms.DataAnalysis.Experimental {
@@ -168 +168 @@
       // A dictionary is used to find parameters
       double[] initialConstants;
-      var parameters = new List<TreeToAutoDiffTermConverter.DataForVariable>();
-
-      TreeToAutoDiffTermConverter.ParametricFunction func;
-      TreeToAutoDiffTermConverter.ParametricFunctionGradient func_grad;
-      TreeToAutoDiffTermConverter.ParametricFunctionGradient func_grad_for_vars;
-      if (!TreeToAutoDiffTermConverter.TryConvertToAutoDiff(tree, updateVariableWeights, out parameters, out initialConstants, 
-        out func, out func_grad, out func_grad_for_vars))
+      var parameters = new List<TreeToDiffSharpConverter.DataForVariable>();
+
+      Func<DV, D> func;
+      if (!TreeToDiffSharpConverter.TryConvertToDiffSharp(tree, updateVariableWeights, out parameters, out initialConstants,
+        out func))
         throw new NotSupportedException("Could not optimize constants of symbolic expression tree due to not supported symbols used in the tree.");
       if (parameters.Count == 0) return 0.0; // gkronber: constant expressions always have a R² of 0.0 
@@ -182 +180 @@
       // extract inital constants
       double[] c = new double[initialConstants.Length + 2];
+      double[] s = new double[c.Length];
       {
-        c[0] = 0.0;
-        c[1] = 1.0;
+        c[0] = 1.0;
+        c[1] = 0.0;
         Array.Copy(initialConstants, 0, c, 2, initialConstants.Length);
-      }
+
+        // s[0] = 1.0;
+        // s[1] = 1.0;
+        // Array.Copy(initialConstants.Select(ci=>Math.Abs(ci)).ToArray()
+        //   , 0, s, 2, initialConstants.Length);
+      }
+      s = Enumerable.Repeat(1.0, c.Length).ToArray();
+
       double[] originalConstants = (double[])c.Clone();
       double originalQuality = SymbolicRegressionSingleObjectivePearsonRSquaredEvaluator.Calculate(interpreter, tree, lowerEstimationLimit, upperEstimationLimit, problemData, rows, applyLinearScaling);
@@ -195 +201 @@
       IDataset ds = problemData.Dataset;
       double[,] x = new double[rows.Count(), parameters.Count];
-      double[,] constraints = new double[rows.Count(), parameters.Count + 1]; // +1 for constraint for f(x)
-      string[,] comp = string[rows.Count(), parameters.Count + 1];
+      int col = 0;
       int row = 0;
-      foreach (var r in rows) {
-        int col = 0;
-        foreach (var info in parameterEntries) {
+      foreach (var info in parameterEntries) {
+        row = 0;
+        foreach (var r in rows) {
           if (ds.VariableHasType<double>(info.variableName)) {
             x[row, col] = ds.GetDoubleValue(info.variableName, r + info.lag);
@@ -207 +212 @@
           } else throw new InvalidProgramException("found a variable of unknown type");
 
-          // find the matching df/dx column
-          var colName = string.Format("df/d({0})", info.variableName);
-          constraints[row, col] = ds.GetDoubleValue(colName, r);
-
-          var compColName = string.Format("df/d({0}) constraint-type")
-          comp[row, col] = ds.GetStringValue(compColName, r);
-          col++;
-        }
+          row++;
+        }
+        col++;
+      }
+
+      var target = problemData.TargetVariable;
+      var constraintRows = Enumerable.Range(0, problemData.Dataset.Rows).Where(rIdx => double.IsNaN(ds.GetDoubleValue(target, rIdx)));
+      double[,] constraints = new double[constraintRows.Count(), parameters.Count + 1]; // +1 for constraint for f(x)
+      string[,] comp = new string[constraintRows.Count(), parameters.Count + 1];
+      int eqConstraints = 0;
+      int ieqConstraints = 0;
+      col = 0;
+      foreach (var info in parameterEntries) {
+        row = 0;
+        // find the matching df/dx column
+        var colName = string.Format("df/d({0})", info.variableName);
+        var compColName = string.Format("df/d({0}) constraint-type", info.variableName);
+
+        if (ds.VariableNames.Contains(colName)) {
+          foreach (var r in constraintRows) {
+            constraints[row, col] = ds.GetDoubleValue(colName, r);
+            comp[row, col] = ds.GetStringValue(compColName, r);
+
+            if (comp[row, col] == "EQ") eqConstraints++;
+            else if (comp[row, col] == "LEQ" || comp[row, col] == "GEQ") ieqConstraints++;
+
+            row++;
+          }
+        }
+        col++;
+      }
+      // f(x) constraint
+      row = 0;
+      col = constraints.GetLength(1) - 1;
+      foreach (var r in constraintRows) {
         constraints[row, col] = ds.GetDoubleValue("f(x)", r);
         comp[row, col] = ds.GetStringValue("f(x) constraint-type", r);
+        if (comp[row, col] == "EQ") eqConstraints++;
+        else if (comp[row, col] == "LEQ" || comp[row, col] == "GEQ") ieqConstraints++;
         row++;
       }
+
+
       double[] y = ds.GetDoubleValues(problemData.TargetVariable, rows).ToArray();
       int n = x.GetLength(0);
@@ -224 +260 @@
       int k = c.Length;
 
-      // alglib.ndimensional_pfunc function_cx_1_func = CreatePFunc(func);
-      // alglib.ndimensional_pgrad function_cx_1_grad = CreatePGrad(func_grad);
-      alglib.ndimensional_jac jac = CreateJac(x, y, constraints, comp, func, func_grad, func_grad_for_vars);
-      double[] s = c.Select(ci=>Math.Max(Math.Abs(ci), 1E-6)).ToArray(); // use absolute value of variables as scale
+      alglib.ndimensional_jac jac = CreateJac(x, y, constraints, comp, func, AD.Grad(func));
       double rho = 1000;
       int outeriters = 3;
@@ -243 +276 @@
         alglib.minnlcsetprecexactlowrank(state, updateFreq);
         // TODO set constraints;
-        alglib.minnlcsetnlc(state, 0, 2);
+        alglib.minnlcsetnlc(state, eqConstraints, ieqConstraints);
         alglib.minnlcoptimize(state, jac, null, null);
         alglib.minnlcresults(state, out c, out rep);
@@ -286 +319 @@
       double[,] x, // x longer than y 
       double[] y, // only targets
+      double[,] constraints, // df/d(xi), same order as for x
       string[,] comparison, // {LEQ, GEQ, EQ }
-      double[,] constraints, // df/d(xi), same order as for x
-      TreeToAutoDiffTermConverter.ParametricFunction func,
-      TreeToAutoDiffTermConverter.ParametricFunctionGradient func_grad,
-      TreeToAutoDiffTermConverter.ParametricFunctionGradient func_grad_for_vars) {
+      Func<DV, D> func,
+      Func<DV, DV> func_grad) {
       return (double[] c, double[] fi, double[,] jac, object o) => {
         // objective function is sum of squared errors
@@ -298 +330 @@
         Array.Clear(fi, 0, fi.Length);
         Array.Clear(jac, 0, jac.Length);
-        double[] xi = new double[nParams];
+        var p = new double[nParams + c.Length];
+        Array.Copy(c, 0, p, nParams, c.Length); // copy c to the end of the function parameters vector
         for (int rowIdx = 0; rowIdx < nRows; rowIdx++) {
-          // copy row 
+          // copy x_i to the beginning of the function parameters vector
           for (int cIdx = 0; cIdx < nParams; cIdx++)
-            xi[cIdx] = x[rowIdx, cIdx];
-          var fg = func_grad(c, xi);
-          double f = fg.Item2;
-          double[] g = fg.Item1;
+            p[cIdx] = x[rowIdx, cIdx];
+
+          double f = (double)func(p);
+          double[] g = (double[])func_grad(p);
           var e = y[rowIdx] - f;
           fi[0] += e * e;
           // update gradient
           for (int colIdx = 0; colIdx < c.Length; colIdx++) {
-            jac[0, colIdx] += -2 * e * g[colIdx];
+            jac[0, colIdx] += -2 * e * g[nParams + colIdx]; // skip the elements for the variable values
           }
         }
 
-        // constraints
-        var nConstraintPoints = constraints.GetLength(0);
-
-        // TODO: AutoDiff for gradients d/d(c) d/d(xi) f(xi,c) for all xi
-        // converter should produce the differentials for all variables as functions which can be differentiated wrt the parameters c
-
+        int fidx = 1;
+        int constraintRows = constraints.GetLength(0);
+
+        // eq constraints
+        for (int rowIdx = 0; rowIdx < constraintRows; rowIdx++) {
+          for (int colIdx = 0; colIdx < constraints.GetLength(1); colIdx++) {
+            if (comparison[rowIdx, colIdx] == "EQ") {
+              throw new NotSupportedException();
+            }
+          }
+        }
+        // ineq constraints
+        for (int rowIdx = 0; rowIdx < constraintRows; rowIdx++) {
+          for (int colIdx = 0; colIdx < constraints.GetLength(1); colIdx++) {
+            // there is a constraint value
+            if (!double.IsNaN(constraints[rowIdx, colIdx]) && !string.IsNullOrEmpty(comparison[rowIdx, colIdx])) {
+              var factor = (comparison[rowIdx, colIdx] == "LEQ") ? 1.0
+              : comparison[rowIdx, colIdx] == "GEQ" ? -1.0 : 0.0;
+
+              // f(x) constraint
+              if (colIdx == constraints.GetLength(1) - 1) {
+                // copy x_i to the beginning of the function parameters vector
+                for (int cIdx = 0; cIdx < nParams; cIdx++)
+                  p[cIdx] = x[rowIdx, cIdx];
+
+                fi[fidx] = factor * ((double)(func(p)) - constraints[rowIdx, colIdx]);
+                var g = (double[])func_grad(p);
+                for (int jacIdx = 0; jacIdx < c.Length; jacIdx++) {
+                  jac[fidx, jacIdx] = factor * g[nParams + jacIdx]; // skip the elements for the variable values
+                }
+                fidx++;
+              } else {
+                // df / dxi constraint
+                var g = (double[])func_grad(p);
+                fi[fidx] = factor * g[colIdx];
+                var hess = AD.Hessian(func, p);
+                for (int jacIdx = 0; jacIdx < c.Length; jacIdx++) {
+                  jac[fidx, jacIdx] = factor * (double)hess[nParams + colIdx, nParams + jacIdx]; // skip the elements for the variable values
+                }
+                fidx++;
+              }
+            }
+          }
+        }
       };
     }
 
-    // private static alglib.ndimensional_pfunc CreatePFunc(TreeToAutoDiffTermConverter.ParametricFunction func) {
-    //   return (double[] c, double[] x, ref double fx, object o) => {
-    //     fx = func(c, x);
-    //   };
-    // }
-    // 
-    // private static alglib.ndimensional_pgrad CreatePGrad(TreeToAutoDiffTermConverter.ParametricFunctionGradient func_grad) {
-    //   return (double[] c, double[] x, ref double fx, double[] grad, object o) => {
-    //     var tupel = func_grad(c, x);
-    //     fx = tupel.Item2;
-    //     Array.Copy(tupel.Item1, grad, grad.Length);
-    //   };
-    // }
     public static bool CanOptimizeConstants(ISymbolicExpressionTree tree) {
       return TreeToAutoDiffTermConverter.IsCompatible(tree);

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/TreeToDiffSharpConverter.cs

@@ -24 +24 @@
 using System.Linq;
 using System.Runtime.Serialization;
-using AutoDiff;
 using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
+using HeuristicLab.Problems.DataAnalysis.Symbolic;
+using DiffSharp.Interop.Float64;
+using System.Linq.Expressions;
+using System.Reflection;
 
 namespace HeuristicLab.Algorithms.DataAnalysis.Experimental {
-  public class TreeToAutoDiffTermConverter {
-    public delegate double ParametricFunction(double[] vars, double[] @params);
-
-    public delegate Tuple<double[], double> ParametricFunctionGradient(double[] vars, double[] @params);
+  public class TreeToDiffSharpConverter {
+    public delegate double ParametricFunction(double[] vars);
+
+    public delegate Tuple<double[], double> ParametricFunctionGradient(double[] vars);
 
     #region helper class
@@ -59 +62 @@
     #endregion
 
-    #region derivations of functions
-    // create function factory for arctangent
-    private static readonly Func<Term, UnaryFunc> arctan = UnaryFunc.Factory(
-      eval: Math.Atan,
-      diff: x => 1 / (1 + x * x));
-
-    private static readonly Func<Term, UnaryFunc> sin = UnaryFunc.Factory(
-      eval: Math.Sin,
-      diff: Math.Cos);
-
-    private static readonly Func<Term, UnaryFunc> cos = UnaryFunc.Factory(
-      eval: Math.Cos,
-      diff: x => -Math.Sin(x));
-
-    private static readonly Func<Term, UnaryFunc> tan = UnaryFunc.Factory(
-      eval: Math.Tan,
-      diff: x => 1 + Math.Tan(x) * Math.Tan(x));
-
-    private static readonly Func<Term, UnaryFunc> erf = UnaryFunc.Factory(
-      eval: alglib.errorfunction,
-      diff: x => 2.0 * Math.Exp(-(x * x)) / Math.Sqrt(Math.PI));
-
-    private static readonly Func<Term, UnaryFunc> norm = UnaryFunc.Factory(
-      eval: alglib.normaldistribution,
-      diff: x => -(Math.Exp(-(x * x)) * Math.Sqrt(Math.Exp(x * x)) * x) / Math.Sqrt(2 * Math.PI));
-
-    #endregion
-
-    public static bool TryConvertToAutoDiff(ISymbolicExpressionTree tree, bool makeVariableWeightsVariable,
+
+    public static bool TryConvertToDiffSharp(ISymbolicExpressionTree tree, bool makeVariableWeightsVariable,
       out List<DataForVariable> parameters, out double[] initialConstants,
-      out ParametricFunction func,
-      out ParametricFunctionGradient func_grad,
-      out ParametricFunctionGradient func_grad_for_vars) {
+      out Func<DV, D> func) {
 
       // use a transformator object which holds the state (variable list, parameter list, ...) for recursive transformation of the tree
-      var transformator = new TreeToAutoDiffTermConverter(makeVariableWeightsVariable);
-      AutoDiff.Term term;
+      var transformator = new TreeToDiffSharpConverter(makeVariableWeightsVariable);
       try {
-        term = transformator.ConvertToAutoDiff(tree.Root.GetSubtree(0));
+
+        // the list of variable names represents the names for dv[0] ... dv[d-1] where d is the number of input variables
+        // the remaining entries of d represent the parameter values
+        transformator.ExtractParameters(tree.Root.GetSubtree(0));
+
+        var lambda = transformator.CreateDelegate(tree, transformator.parameters);
+        func = lambda.Compile();
+
         var parameterEntries = transformator.parameters.ToArray(); // guarantee same order for keys and values
-        var compiledTerm = term.Compile(
-          transformator.variables.ToArray(),
-          parameterEntries.Select(kvp => kvp.Value).ToArray());
-        
         parameters = new List<DataForVariable>(parameterEntries.Select(kvp => kvp.Key));
         initialConstants = transformator.initialConstants.ToArray();
-        func = (vars, @params) => compiledTerm.Evaluate(vars, @params);
-        func_grad = (vars, @params) => compiledTerm.Differentiate(vars, @params);
-        func_grad_for_vars = (vars, @params) => compiledTerm.Differentiate(@params,vars);
         return true;
       } catch (ConversionException) {
         func = null;
-        func_grad = null;
-        func_grad_for_vars = null;
         parameters = null;
         initialConstants = null;
@@ -119 +90 @@
     }
 
+    public Expression<Func<DV, D>> CreateDelegate(ISymbolicExpressionTree tree, Dictionary<DataForVariable, int> parameters) {
+      paramIdx = parameters.Count; // first non-variable parameter
+      var dv = Expression.Parameter(typeof(DV));
+      var expr = MakeExpr(tree.Root.GetSubtree(0), parameters, dv);
+      var lambda = Expression.Lambda<Func<DV, D>>(expr, dv);
+      return lambda;
+    }
+
     // state for recursive transformation of trees 
-    private readonly
-    List<double> initialConstants;
-    private readonly Dictionary<DataForVariable, AutoDiff.Variable> parameters;
-    private readonly List<AutoDiff.Variable> variables;
+    private readonly List<double> initialConstants;
+    private readonly Dictionary<DataForVariable, int> parameters;
     private readonly bool makeVariableWeightsVariable;
-
-    private TreeToAutoDiffTermConverter(bool makeVariableWeightsVariable) {
+    private int paramIdx;
+
+    private TreeToDiffSharpConverter(bool makeVariableWeightsVariable) {
       this.makeVariableWeightsVariable = makeVariableWeightsVariable;
       this.initialConstants = new List<double>();
-      this.parameters = new Dictionary<DataForVariable, AutoDiff.Variable>();
-      this.variables = new List<AutoDiff.Variable>();
-    }
-
-    private AutoDiff.Term ConvertToAutoDiff(ISymbolicExpressionTreeNode node) {
-      if (node.Symbol is Constant) {
+      this.parameters = new Dictionary<DataForVariable, int>();
+    }
+
+    private void ExtractParameters(ISymbolicExpressionTreeNode node) {
+      if (node.Symbol is HeuristicLab.Problems.DataAnalysis.Symbolic.Constant) {
         initialConstants.Add(((ConstantTreeNode)node).Value);
-        var var = new AutoDiff.Variable();
-        variables.Add(var);
-        return var;
-      }
-      if (node.Symbol is Variable || node.Symbol is BinaryFactorVariable) {
+      } else if (node.Symbol is HeuristicLab.Problems.DataAnalysis.Symbolic.Variable || node.Symbol is BinaryFactorVariable) {
+        var varNode = node as VariableTreeNodeBase;
+        var factorVarNode = node as BinaryFactorVariableTreeNode;
+        // factor variable values are only 0 or 1 and set in x accordingly
+        var varValue = factorVarNode != null ? factorVarNode.VariableValue : string.Empty;
+        FindOrCreateParameter(parameters, varNode.VariableName, varValue);
+
+        if (makeVariableWeightsVariable) {
+          initialConstants.Add(varNode.Weight);
+        }
+      } else if (node.Symbol is FactorVariable) {
+        var factorVarNode = node as FactorVariableTreeNode;
+        var products = new List<D>();
+        foreach (var variableValue in factorVarNode.Symbol.GetVariableValues(factorVarNode.VariableName)) {
+          FindOrCreateParameter(parameters, factorVarNode.VariableName, variableValue);
+
+          initialConstants.Add(factorVarNode.GetValue(variableValue));
+        }
+      } else if (node.Symbol is LaggedVariable) {
+        var varNode = node as LaggedVariableTreeNode;
+        FindOrCreateParameter(parameters, varNode.VariableName, string.Empty, varNode.Lag);
+
+        if (makeVariableWeightsVariable) {
+          initialConstants.Add(varNode.Weight);
+        }
+      } else if (node.Symbol is Addition) {
+        foreach (var subTree in node.Subtrees) {
+          ExtractParameters(subTree);
+        }
+      } else if (node.Symbol is Subtraction) {
+        for (int i = 0; i < node.SubtreeCount; i++) {
+          ExtractParameters(node.GetSubtree(i));
+        }
+      } else if (node.Symbol is Multiplication) {
+        foreach (var subTree in node.Subtrees) {
+          ExtractParameters(subTree);
+        }
+      } else if (node.Symbol is Division) {
+        foreach (var subTree in node.Subtrees) {
+          ExtractParameters(subTree);
+        }
+      } else if (node.Symbol is Logarithm) {
+        ExtractParameters(node.GetSubtree(0));
+      } else if (node.Symbol is Exponential) {
+        ExtractParameters(node.GetSubtree(0));
+      } else if (node.Symbol is Square) {
+        ExtractParameters(node.GetSubtree(0));
+      } else if (node.Symbol is SquareRoot) {
+        ExtractParameters(node.GetSubtree(0));
+      } else if (node.Symbol is Sine) {
+        ExtractParameters(node.GetSubtree(0));
+      } else if (node.Symbol is Cosine) {
+        ExtractParameters(node.GetSubtree(0));
+      } else if (node.Symbol is Tangent) {
+        ExtractParameters(node.GetSubtree(0));
+      } else if (node.Symbol is StartSymbol) {
+        ExtractParameters(node.GetSubtree(0));
+      } else throw new ConversionException();
+    }
+
+    private Func<DV, D> CreateDiffSharpFunc(ISymbolicExpressionTreeNode node, Dictionary<DataForVariable, int> parameters) {
+      this.paramIdx = parameters.Count; // first idx of non-variable parameter      
+      var f = CreateDiffSharpFunc(node, parameters);
+      return (DV paramValues) => f(paramValues);
+    }
+
+    private static readonly MethodInfo DvIndexer = typeof(DV).GetMethod("get_Item", new[] { typeof(int) });
+    private static readonly MethodInfo d_Add_d = typeof(D).GetMethod("op_Addition", new[] { typeof(D), typeof(D) });
+    private static readonly MethodInfo d_Neg = typeof(D).GetMethod("Neg", new[] { typeof(D) });
+    private static readonly MethodInfo d_Mul_d = typeof(D).GetMethod("op_Multiply", new[] { typeof(D), typeof(D) });
+    private static readonly MethodInfo d_Mul_f = typeof(D).GetMethod("op_Multiply", new[] { typeof(D), typeof(double) });
+    private static readonly MethodInfo d_Div_d = typeof(D).GetMethod("op_Division", new[] { typeof(D), typeof(D) });
+    private static readonly MethodInfo f_Div_d = typeof(D).GetMethod("op_Division", new[] { typeof(double), typeof(D) });
+    private static readonly MethodInfo d_Sub_d = typeof(D).GetMethod("op_Subtraction", new[] { typeof(D), typeof(D) });
+    private static readonly MethodInfo d_Pow_f = typeof(D).GetMethod("Pow", new[] { typeof(D), typeof(double) });
+    private static readonly MethodInfo d_Log = typeof(D).GetMethod("Log", new[] { typeof(D) });
+    private static readonly MethodInfo d_Exp = typeof(D).GetMethod("Exp", new[] { typeof(D) });
+
+
+
+    private Expression MakeExpr(ISymbolicExpressionTreeNode node, Dictionary<DataForVariable, int> parameters, ParameterExpression dv) {
+      if (node.Symbol is HeuristicLab.Problems.DataAnalysis.Symbolic.Constant) {
+        return Expression.Call(dv, DvIndexer, Expression.Constant(paramIdx++));
+      }
+      if (node.Symbol is HeuristicLab.Problems.DataAnalysis.Symbolic.Variable || node.Symbol is BinaryFactorVariable) {
         var varNode = node as VariableTreeNodeBase;
         var factorVarNode = node as BinaryFactorVariableTreeNode;
@@ -148 +205 @@
 
         if (makeVariableWeightsVariable) {
-          initialConstants.Add(varNode.Weight);
-          var w = new AutoDiff.Variable();
-          variables.Add(w);
-          return AutoDiff.TermBuilder.Product(w, par);
+          var w = Expression.Call(dv, DvIndexer, Expression.Constant(paramIdx++));
+          var v = Expression.Call(dv, DvIndexer, Expression.Constant(par));
+          return Expression.Call(d_Mul_d, w, v);
         } else {
-          return varNode.Weight * par;
+          var w = Expression.Constant(varNode.Weight);
+          var v = Expression.Call(dv, DvIndexer, Expression.Constant(par));
+          return Expression.Call(d_Mul_f, v, w);
         }
       }
       if (node.Symbol is FactorVariable) {
         var factorVarNode = node as FactorVariableTreeNode;
-        var products = new List<Term>();
-        foreach (var variableValue in factorVarNode.Symbol.GetVariableValues(factorVarNode.VariableName)) {
+        var products = new List<D>();
+        var firstValue = factorVarNode.Symbol.GetVariableValues(factorVarNode.VariableName).First();
+        var parForFirstValue = FindOrCreateParameter(parameters, factorVarNode.VariableName, firstValue);
+        var weightForFirstValue = Expression.Call(dv, DvIndexer, Expression.Constant(paramIdx++));
+        var valForFirstValue = Expression.Call(dv, DvIndexer, Expression.Constant(parForFirstValue));
+        var res = Expression.Call(d_Mul_d, weightForFirstValue, valForFirstValue);
+
+        foreach (var variableValue in factorVarNode.Symbol.GetVariableValues(factorVarNode.VariableName).Skip(1)) {
           var par = FindOrCreateParameter(parameters, factorVarNode.VariableName, variableValue);
-
-          initialConstants.Add(factorVarNode.GetValue(variableValue));
-          var wVar = new AutoDiff.Variable();
-          variables.Add(wVar);
-
-          products.Add(AutoDiff.TermBuilder.Product(wVar, par));
-        }
-        return AutoDiff.TermBuilder.Sum(products);
-      }
-      if (node.Symbol is LaggedVariable) {
-        var varNode = node as LaggedVariableTreeNode;
-        var par = FindOrCreateParameter(parameters, varNode.VariableName, string.Empty, varNode.Lag);
-
-        if (makeVariableWeightsVariable) {
-          initialConstants.Add(varNode.Weight);
-          var w = new AutoDiff.Variable();
-          variables.Add(w);
-          return AutoDiff.TermBuilder.Product(w, par);
+      
+          var weight = Expression.Call(dv, DvIndexer, Expression.Constant(paramIdx++));
+          var v = Expression.Call(dv, DvIndexer, Expression.Constant(par));
+
+          res = Expression.Call(d_Add_d, res, Expression.Call(d_Mul_d, weight, v));
+        }
+        return res;
+      }
+      // if (node.Symbol is LaggedVariable) {
+      //   var varNode = node as LaggedVariableTreeNode;
+      //   var par = FindOrCreateParameter(parameters, varNode.VariableName, string.Empty, varNode.Lag);
+      // 
+      //   if (makeVariableWeightsVariable) {
+      //     initialConstants.Add(varNode.Weight);
+      //     var w = paramValues[paramIdx++];
+      //     return w * paramValues[par];
+      //   } else {
+      //     return varNode.Weight * paramValues[par];
+      //   }
+      // }
+      if (node.Symbol is Addition) {
+        var f = MakeExpr(node.Subtrees.First(), parameters, dv);
+
+        foreach (var subTree in node.Subtrees.Skip(1)) {
+          f = Expression.Call(d_Add_d, f, MakeExpr(subTree, parameters, dv));
+        }
+        return f;
+      }
+      if (node.Symbol is Subtraction) {
+        if (node.SubtreeCount == 1) {
+          return Expression.Call(d_Neg, MakeExpr(node.Subtrees.First(), parameters, dv));
         } else {
-          return varNode.Weight * par;
-        }
-      }
-      if (node.Symbol is Addition) {
-        List<AutoDiff.Term> terms = new List<Term>();
-        foreach (var subTree in node.Subtrees) {
-          terms.Add(ConvertToAutoDiff(subTree));
-        }
-        return AutoDiff.TermBuilder.Sum(terms);
-      }
-      if (node.Symbol is Subtraction) {
-        List<AutoDiff.Term> terms = new List<Term>();
-        for (int i = 0; i < node.SubtreeCount; i++) {
-          AutoDiff.Term t = ConvertToAutoDiff(node.GetSubtree(i));
-          if (i > 0) t = -t;
-          terms.Add(t);
-        }
-        if (terms.Count == 1) return -terms[0];
-        else return AutoDiff.TermBuilder.Sum(terms);
+          var f = MakeExpr(node.Subtrees.First(), parameters, dv);
+
+          foreach (var subTree in node.Subtrees.Skip(1)) {
+            f = Expression.Call(d_Sub_d, f, MakeExpr(subTree, parameters, dv));
+          }
+          return f;
+        }
       }
       if (node.Symbol is Multiplication) {
-        List<AutoDiff.Term> terms = new List<Term>();
-        foreach (var subTree in node.Subtrees) {
-          terms.Add(ConvertToAutoDiff(subTree));
-        }
-        if (terms.Count == 1) return terms[0];
-        else return terms.Aggregate((a, b) => new AutoDiff.Product(a, b));
+        var f = MakeExpr(node.Subtrees.First(), parameters, dv);
+        foreach (var subTree in node.Subtrees.Skip(1)) {
+          f = Expression.Call(d_Mul_d, f, MakeExpr(subTree, parameters, dv));
+        }
+        return f;
       }
       if (node.Symbol is Division) {
-        List<AutoDiff.Term> terms = new List<Term>();
-        foreach (var subTree in node.Subtrees) {
-          terms.Add(ConvertToAutoDiff(subTree));
-        }
-        if (terms.Count == 1) return 1.0 / terms[0];
-        else return terms.Aggregate((a, b) => new AutoDiff.Product(a, 1.0 / b));
+        if (node.SubtreeCount == 1) {
+          return Expression.Call(f_Div_d, Expression.Constant(1.0), MakeExpr(node.Subtrees.First(), parameters, dv));
+        } else {
+          var f = MakeExpr(node.Subtrees.First(), parameters, dv);
+
+          foreach (var subTree in node.Subtrees.Skip(1)) {
+            f = Expression.Call(d_Div_d, f, MakeExpr(subTree, parameters, dv));
+          }
+          return f;
+        }
       }
       if (node.Symbol is Logarithm) {
-        return AutoDiff.TermBuilder.Log(
-          ConvertToAutoDiff(node.GetSubtree(0)));
+        return Expression.Call(d_Log, MakeExpr(node.GetSubtree(0), parameters, dv));
       }
       if (node.Symbol is Exponential) {
-        return AutoDiff.TermBuilder.Exp(
-          ConvertToAutoDiff(node.GetSubtree(0)));
+        return Expression.Call(d_Exp, MakeExpr(node.GetSubtree(0), parameters, dv));
       }
       if (node.Symbol is Square) {
-        return AutoDiff.TermBuilder.Power(
-          ConvertToAutoDiff(node.GetSubtree(0)), 2.0);
+        return Expression.Call(d_Pow_f, MakeExpr(node.GetSubtree(0), parameters, dv), Expression.Constant(2.0));
       }
       if (node.Symbol is SquareRoot) {
-        return AutoDiff.TermBuilder.Power(
-          ConvertToAutoDiff(node.GetSubtree(0)), 0.5);
-      }
-      if (node.Symbol is Sine) {
-        return sin(
-          ConvertToAutoDiff(node.GetSubtree(0)));
-      }
-      if (node.Symbol is Cosine) {
-        return cos(
-          ConvertToAutoDiff(node.GetSubtree(0)));
-      }
-      if (node.Symbol is Tangent) {
-        return tan(
-          ConvertToAutoDiff(node.GetSubtree(0)));
-      }
-      if (node.Symbol is Erf) {
-        return erf(
-          ConvertToAutoDiff(node.GetSubtree(0)));
-      }
-      if (node.Symbol is Norm) {
-        return norm(
-          ConvertToAutoDiff(node.GetSubtree(0)));
-      }
+        return Expression.Call(d_Pow_f, MakeExpr(node.GetSubtree(0), parameters, dv), Expression.Constant(0.5));
+      }
+      // if (node.Symbol is Sine) {
+      //   return AD.Sin(CreateDiffSharpFunc(node.GetSubtree(0), parameters, paramValues));
+      // }
+      // if (node.Symbol is Cosine) {
+      //   return AD.Cos(CreateDiffSharpFunc(node.GetSubtree(0), parameters, paramValues));
+      // }
+      // if (node.Symbol is Tangent) {
+      //   return AD.Tan(CreateDiffSharpFunc(node.GetSubtree(0), parameters, paramValues));
+      // }
       if (node.Symbol is StartSymbol) {
-        var alpha = new AutoDiff.Variable();
-        var beta = new AutoDiff.Variable();
-        variables.Add(beta);
-        variables.Add(alpha);
-        return ConvertToAutoDiff(node.GetSubtree(0)) * alpha + beta;
+        var alpha = Expression.Call(dv, DvIndexer, Expression.Constant(paramIdx++));
+        var beta = Expression.Call(dv, DvIndexer, Expression.Constant(paramIdx++));
+
+        return Expression.Call(d_Add_d, beta,
+          Expression.Call(d_Mul_d, alpha, MakeExpr(node.GetSubtree(0), parameters, dv)));
       }
       throw new ConversionException();
@@ -265 +318 @@
     // for each factor variable value we need a parameter which represents a binary indicator for that variable & value combination
     // each binary indicator is only necessary once. So we only create a parameter if this combination is not yet available.
-    private static Term FindOrCreateParameter(Dictionary<DataForVariable, AutoDiff.Variable> parameters,
+    private static int FindOrCreateParameter(Dictionary<DataForVariable, int> parameters,
       string varName, string varValue = "", int lag = 0) {
       var data = new DataForVariable(varName, varValue, lag);
-
-      AutoDiff.Variable par = null;
-      if (!parameters.TryGetValue(data, out par)) {
-        // not found -> create new parameter and entries in names and values lists
-        par = new AutoDiff.Variable();
-        parameters.Add(data, par);
-      }
-      return par;
+      int idx = -1;
+      if (parameters.TryGetValue(data, out idx)) return idx;
+      else parameters[data] = parameters.Count;
+      return idx;
     }
 
@@ -282 +331 @@
         from n in tree.Root.GetSubtree(0).IterateNodesPrefix()
         where
-          !(n.Symbol is Variable) &&
+          !(n.Symbol is HeuristicLab.Problems.DataAnalysis.Symbolic.Variable) &&
           !(n.Symbol is BinaryFactorVariable) &&
           !(n.Symbol is FactorVariable) &&
           !(n.Symbol is LaggedVariable) &&
-          !(n.Symbol is Constant) &&
+          !(n.Symbol is HeuristicLab.Problems.DataAnalysis.Symbolic.Constant) &&
           !(n.Symbol is Addition) &&
           !(n.Symbol is Subtraction) &&
@@ -298 +347 @@
           !(n.Symbol is Cosine) &&
           !(n.Symbol is Tangent) &&
-          !(n.Symbol is Erf) &&
-          !(n.Symbol is Norm) &&
           !(n.Symbol is StartSymbol)
         select n).Any();

branches/MathNetNumerics-Exploration-2789/HeuristicLab.ExtLibs.MathNetNumerics/Plugin.cs.frame

@@ -26 +26 @@
   [PluginFile("HeuristicLab.ExtLibs.MathNetNumerics.dll", PluginFileType.Assembly)]
   [PluginFile("MathNet.Numerics.dll", PluginFileType.Assembly)]
-  [PluginFile("x64/libiomp5md.dll", PluginFileType.NativeDll)]
-  [PluginFile("x64/MathNet.Numerics.MKL.dll", PluginFileType.NativeDll)]
-  [PluginFile("x86/libiomp5md.dll", PluginFileType.NativeDll)]
-  [PluginFile("x86/MathNet.Numerics.MKL.dll", PluginFileType.NativeDll)]
+  //[PluginFile("x64/libiomp5md.dll", PluginFileType.NativeDll)]
+  //[PluginFile("x64/MathNet.Numerics.MKL.dll", PluginFileType.NativeDll)]
+  //[PluginFile("x86/libiomp5md.dll", PluginFileType.NativeDll)]
+  //[PluginFile("x86/MathNet.Numerics.MKL.dll", PluginFileType.NativeDll)]
   public class Plugin : PluginBase {
   }

branches/MathNetNumerics-Exploration-2789/MathNetNumerics-Exploration.sln

@@ -7 +7 @@
 EndProject
 Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "HeuristicLab.Algorithms.DataAnalysis.Experimental", "HeuristicLab.Algorithms.DataAnalysis.Experimental\HeuristicLab.Algorithms.DataAnalysis.Experimental.csproj", "{D97F91B5-B71B-412D-90A5-177EED948A3A}"
+EndProject
+Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Test", "Test\Test.csproj", "{F5F7CF4D-F560-4B26-9552-21F2DE441F1F}"
 EndProject
 Global
@@ -22 +24 @@
     {D97F91B5-B71B-412D-90A5-177EED948A3A}.Release|Any CPU.ActiveCfg = Release|Any CPU
     {D97F91B5-B71B-412D-90A5-177EED948A3A}.Release|Any CPU.Build.0 = Release|Any CPU
+    {F5F7CF4D-F560-4B26-9552-21F2DE441F1F}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+    {F5F7CF4D-F560-4B26-9552-21F2DE441F1F}.Debug|Any CPU.Build.0 = Debug|Any CPU
+    {F5F7CF4D-F560-4B26-9552-21F2DE441F1F}.Release|Any CPU.ActiveCfg = Release|Any CPU
+    {F5F7CF4D-F560-4B26-9552-21F2DE441F1F}.Release|Any CPU.Build.0 = Release|Any CPU
   EndGlobalSection
   GlobalSection(SolutionProperties) = preSolution