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

Timestamp:

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

Author:

gkronber

Message:

#2789: worked on nonlinear regression with constraints

File:

• branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/TreeToDiffSharpConverter.cs (moved) (moved from branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/TreeToAutoDiffTermConverter.cs) (7 diffs)

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();