Changeset 16999 for branches

Timestamp:

05/31/19 13:55:51 (5 years ago)

Author:

gkronber

Message:

#2925: Added optimization of weights for variables and added an integration method which uses CVODES to integrate over the whole episode (without input variables)

Location:

branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3

Files:

• HeuristicLab.Problems.DynamicalSystemsModelling-3.3.csproj (modified) (1 diff)
• Problem.cs (modified) (21 diffs)
• ProblemInstanceProvider.cs (modified) (3 diffs)

branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/HeuristicLab.Problems.DynamicalSystemsModelling-3.3.csproj

@@ -58 +58 @@
     <CodeAnalysisRuleSet>AllRules.ruleset</CodeAnalysisRuleSet>
     <Prefer32Bit>false</Prefer32Bit>
+    <AllowUnsafeBlocks>true</AllowUnsafeBlocks>
   </PropertyGroup>
   <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Debug|x64' ">

branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs

@@ -229 +229 @@
       Parameters.Add(new FixedValueParameter<DoubleValue>("Numeric differences smoothing", "Determines the amount of smoothing for the numeric differences which are calculated for pre-tuning. Values from -8 to 8 are reasonable. Use very low value if the data contains no noise. Default: 2.", new DoubleValue(2.0)));
 
-      var solversStr = new string[] { "HeuristicLab", "CVODES" };
+      var solversStr = new string[] { "HeuristicLab", "CVODES", "CVODES (full)" };
       var solvers = new ItemSet<StringValue>(
         solversStr.Select(s => new StringValue(s).AsReadOnly())
@@ -238 +238 @@
       InitAllParameters();
 
-      // TODO: use training range as default training episode
       // TODO: optimization of starting values for latent variables in CVODES solver
       // TODO: allow to specify the name for the time variable in the dataset and allow variable step-sizes
@@ -249 +248 @@
       var targetVars = TargetVariables.CheckedItems.OrderBy(i => i.Index).Select(i => i.Value.Value).ToArray();
       var latentVariables = Enumerable.Range(1, NumberOfLatentVariables).Select(i => "λ" + i).ToArray(); // TODO: must coincide with the variables which are actually defined in the grammar and also for which we actually have trees
+      if (latentVariables.Any()) throw new NotSupportedException("latent variables are not supported"); // not sure if everything still works in combination with latent variables
       if (OptimizeParametersForEpisodes) {
         throw new NotImplementedException();
@@ -272 +272 @@
         double nmse = OptimizeForEpisodes(trees, problemData, targetVars, latentVariables, random, trainingEpisodes, MaximumPretuningParameterOptimizationIterations, NumericIntegrationSteps, OdeSolver, MaximumOdeParameterOptimizationIterations,
           PretuningErrorWeight.Value.Value, OdeErrorWeight.Value.Value, NumericDifferencesSmoothing);
-        // individual["OptTheta"] = new DoubleArray(optTheta); // write back optimized parameters so that we can use them in the Analysis method
         return nmse;
       }
@@ -298 +297 @@
       var targetVariableTrees = trees.Take(targetVars.Length).ToArray();
       var latentVariableTrees = trees.Skip(targetVars.Length).ToArray();
-      var constantNodes = targetVariableTrees.Select(t => t.IterateNodesPrefix().OfType<ConstantTreeNode>().ToArray()).ToArray();
-      var initialTheta = constantNodes.Select(nodes => nodes.Select(n => n.Value).ToArray()).ToArray();
+      // var constantNodes = targetVariableTrees.Select(t => t.IterateNodesPrefix().OfType<ConstantTreeNode>().ToArray()).ToArray();
+      // var initialTheta = constantNodes.Select(nodes => nodes.Select(n => n.Value).ToArray()).ToArray();
+      var constantNodes = targetVariableTrees.Select(
+        t => t.IterateNodesPrefix()
+        .Where(n => n.SubtreeCount == 0) // select leaves
+        .ToArray()).ToArray();
+      var initialTheta = constantNodes.Select(
+        a => a.Select(
+          n => {
+            if (n is VariableTreeNode varTreeNode) {
+              return varTreeNode.Weight;
+            } else if (n is ConstantTreeNode constTreeNode) {
+              return constTreeNode.Value;
+            } else throw new InvalidProgramException();
+          }).ToArray()).ToArray();
 
       // optimize parameters by fitting f(x,y) to calculated differences dy/dt(t)
@@ -313 +325 @@
       pretunedParameters = pretunedParameters
         .Concat(latentVariableTrees
-        .SelectMany(t => t.IterateNodesPrefix().OfType<ConstantTreeNode>().Select(n => n.Value)))
+        .SelectMany(t => t.IterateNodesPrefix()
+        .Where(n => n.SubtreeCount == 0)
+        .Select(n => {
+          if (n is VariableTreeNode varTreeNode) {
+            return varTreeNode.Weight;
+          } else if (n is ConstantTreeNode constTreeNode) {
+            return constTreeNode.Value;
+          } else throw new InvalidProgramException();
+        })))
         .ToArray();
 
@@ -332 +352 @@
       var paramIdx = 0;
       for (var treeIdx = 0; treeIdx < constantNodes.Length; treeIdx++) {
-        for (int i = 0; i < constantNodes[treeIdx].Length; i++)
-          constantNodes[treeIdx][i].Value = optTheta[paramIdx++];
+        for (int i = 0; i < constantNodes[treeIdx].Length; i++) {
+          if (constantNodes[treeIdx][i] is VariableTreeNode varTreeNode) {
+            varTreeNode.Weight = optTheta[paramIdx++];
+          } else if (constantNodes[treeIdx][i] is ConstantTreeNode constTreeNode) {
+            constTreeNode.Value = optTheta[paramIdx++];
+          }
+        }
       }
       return nmse;
@@ -359 +384 @@
       if (latentVariables.Length > 0) {
         var inputVariables = targetVars.Concat(latentTrees.SelectMany(t => t.IterateNodesPrefix().OfType<VariableTreeNode>().Select(n => n.VariableName))).Except(latentVariables).Distinct();
-        var myState = new OptimizationData(latentTrees, targetVars, inputVariables.ToArray(), problemData, null, episodes.ToArray(), 10, latentVariables, "HeuristicLab");
+        var myState = new OptimizationData(latentTrees, targetVars, inputVariables.ToArray(), problemData, null, episodes.ToArray(), 10, latentVariables, "NONE");
 
         var fi = new double[myState.rows.Length * targetVars.Length];
@@ -527 +552 @@
         foreach (var variable in variables) {
           // in this problem we also allow fixed numeric parameters (represented as variables with the value as name)
-          if (double.TryParse(variable, NumberStyles.Float, CultureInfo.InvariantCulture, out double value)) {
-            nodeValueLookup.SetVariableValue(variable, value); // TODO: Perf we don't need to set this for each index
-          } else {
-            nodeValueLookup.SetVariableValue(variable, ds.GetDoubleValue(variable, rows[trainIdx])); // TODO: perf
-          }
+          // if (double.TryParse(variable, NumberStyles.Float, CultureInfo.InvariantCulture, out double value)) {
+          //   nodeValueLookup.SetVariableValue(variable, value); // TODO: Perf we don't need to set this for each index
+          // } else {
+          nodeValueLookup.SetVariableValue(variable, ds.GetDoubleValue(variable, rows[trainIdx])); // TODO: perf
+          // }
         }
         // interpret all trees
@@ -561 +586 @@
         foreach (var variable in variables) {
           // in this problem we also allow fixed numeric parameters (represented as variables with the value as name)
-          if (double.TryParse(variable, NumberStyles.Float, CultureInfo.InvariantCulture, out double value)) {
-            nodeValueLookup.SetVariableValue(variable, value); // TODO: Perf we don't need to set this for each index
-          } else {
-            nodeValueLookup.SetVariableValue(variable, ds.GetDoubleValue(variable, rows[trainIdx])); // TODO: perf
-          }
+          // if (double.TryParse(variable, NumberStyles.Float, CultureInfo.InvariantCulture, out double value)) {
+          //   nodeValueLookup.SetVariableValue(variable, value); // TODO: Perf we don't need to set this for each index
+          // } else {
+          nodeValueLookup.SetVariableValue(variable, ds.GetDoubleValue(variable, rows[trainIdx])); // TODO: perf
+          // }
         }
 
@@ -658 +683 @@
 
       var bestIndividualAndQuality = this.GetBestIndividual(individuals, qualities);
-      var trees = bestIndividualAndQuality.Item1.Values.Select(v => v.Value).OfType<ISymbolicExpressionTree>().ToArray(); // extract all trees from individual
+      var trees = bestIndividualAndQuality.Item1.Values.Select(v => v.Value)
+        .OfType<ISymbolicExpressionTree>().ToArray(); // extract all trees from individual
 
       results["SNMSE"].Value = new DoubleValue(bestIndividualAndQuality.Item2);
@@ -966 +992 @@
         foreach (var varName in inputVariables) {
           // in this problem we also allow fixed numeric parameters (represented as variables with the value as name)
-          if (double.TryParse(varName, NumberStyles.Float, CultureInfo.InvariantCulture, out double value)) {
-            nodeValues.SetVariableValue(varName, value, Vector.Zero);
-          } else {
-            var y0 = dataset.GetDoubleValue(varName, t0);
-            nodeValues.SetVariableValue(varName, y0, Vector.Zero);
-          }
+          // if (double.TryParse(varName, NumberStyles.Float, CultureInfo.InvariantCulture, out double value)) {
+          //   nodeValues.SetVariableValue(varName, value, Vector.Zero);
+          // } else {
+          var y0 = dataset.GetDoubleValue(varName, t0);
+          nodeValues.SetVariableValue(varName, y0, Vector.Zero);
+          //}
         }
         foreach (var varName in targetVariables) {
@@ -1013 +1039 @@
 
         var prevT = t0; // TODO: here we should use a variable for t if it is available. Right now we assume equidistant measurements.
-        foreach (var t in rows.Skip(1)) {
-          if (odeSolver == "HeuristicLab")
-            IntegrateHL(trees, calculatedVariables, nodeValues, numericIntegrationSteps); // integrator updates nodeValues
-          else if (odeSolver == "CVODES")
-            IntegrateCVODES(trees, calculatedVariables, nodeValues);
-          else throw new InvalidOperationException("Unknown ODE solver " + odeSolver);
-          prevT = t;
-
-          // update output for target variables (TODO: if we want to visualize the latent variables then we need to provide a separate output)
-          for (int i = 0; i < targetVariables.Length; i++) {
-            var targetVar = targetVariables[i];
-            var yt = nodeValues.GetVariableValue(targetVar);
-
-            // fill up remaining rows with last valid value if there are invalid values
-            if (double.IsNaN(yt.Item1) || double.IsInfinity(yt.Item1)) {
-              for (; outputRowIdx < fi.Length; outputRowIdx++) {
-                var prevIdx = outputRowIdx - targetVariables.Length;
-                fi[outputRowIdx] = fi[prevIdx]; // current <- prev
-                if (jac != null) for (int j = 0; j < jac.GetLength(1); j++) jac[outputRowIdx, j] = jac[prevIdx, j];
+        if (odeSolver == "CVODES (full)") {
+          IntegrateCVODES(trees, calculatedVariables, nodeValues, rows, fi, jac);
+        } else {
+
+          foreach (var t in rows.Skip(1)) {
+            if (odeSolver == "HeuristicLab")
+              IntegrateHL(trees, calculatedVariables, nodeValues, numericIntegrationSteps); // integrator updates nodeValues
+            else if (odeSolver == "CVODES")
+              IntegrateCVODES(trees, calculatedVariables, nodeValues);
+            else throw new InvalidOperationException("Unknown ODE solver " + odeSolver);
+            prevT = t;
+
+            // update output for target variables (TODO: if we want to visualize the latent variables then we need to provide a separate output)
+            for (int i = 0; i < targetVariables.Length; i++) {
+              var targetVar = targetVariables[i];
+              var yt = nodeValues.GetVariableValue(targetVar);
+
+              // fill up remaining rows with last valid value if there are invalid values
+              if (double.IsNaN(yt.Item1) || double.IsInfinity(yt.Item1)) {
+                for (; outputRowIdx < fi.Length; outputRowIdx++) {
+                  var prevIdx = outputRowIdx - targetVariables.Length;
+                  fi[outputRowIdx] = fi[prevIdx]; // current <- prev
+                  if (jac != null) for (int j = 0; j < jac.GetLength(1); j++) jac[outputRowIdx, j] = jac[prevIdx, j];
+                }
+                return;
+              };
+
+              fi[outputRowIdx] = yt.Item1;
+              var g = yt.Item2;
+              g.CopyTo(jac, outputRowIdx);
+              outputRowIdx++;
+            }
+            if (latentValues != null) {
+              foreach (var latentVariable in latentVariables) {
+                var lt = nodeValues.GetVariableValue(latentVariable).Item1;
+                latentValues[latentValueRowIdx, latentValueColIdx++] = lt;
               }
-              return;
-            };
-
-            fi[outputRowIdx] = yt.Item1;
-            var g = yt.Item2;
-            g.CopyTo(jac, outputRowIdx);
-            outputRowIdx++;
-          }
-          if (latentValues != null) {
-            foreach (var latentVariable in latentVariables) {
-              var lt = nodeValues.GetVariableValue(latentVariable).Item1;
-              latentValues[latentValueRowIdx, latentValueColIdx++] = lt;
+              latentValueRowIdx++; latentValueColIdx = 0;
             }
-            latentValueRowIdx++; latentValueColIdx = 0;
-          }
-
-          // update for next time step (only the inputs)
-          foreach (var varName in inputVariables) {
-            // in this problem we also allow fixed numeric parameters (represented as variables with the value as name)
-            if (double.TryParse(varName, NumberStyles.Float, CultureInfo.InvariantCulture, out double value)) {
-              // value is unchanged
-            } else {
+
+            // update for next time step (only the inputs)
+            foreach (var varName in inputVariables) {
+              // in this problem we also allow fixed numeric parameters (represented as variables with the value as name)
+              // if (double.TryParse(varName, NumberStyles.Float, CultureInfo.InvariantCulture, out double value)) {
+              //   // value is unchanged
+              // } else {
               nodeValues.SetVariableValue(varName, dataset.GetDoubleValue(varName, t), Vector.Zero);
+              // }
             }
           }
@@ -1183 +1214 @@
     }
 
+    /// <summary>
+    ///  Here we use CVODES to solve the ODE. Forward sensitivities are used to calculate the gradient for parameter optimization
+    /// </summary>
+    /// <param name="trees">Each equation in the ODE represented as a tree</param>
+    /// <param name="calculatedVariables">The names of the calculated variables</param>
+    /// <param name="t">The time t up to which we need to integrate.</param>
+    private static void IntegrateCVODES(
+      ISymbolicExpressionTree[] trees, // f(y,p) in tree representation
+      string[] calculatedVariables, // names of elements of y
+      NodeValueLookup nodeValues,
+      IEnumerable<int> rows,
+      double[] fi,
+      double[,] jac
+      ) {
+
+      // the RHS of the ODE
+      // dy/dt = f(y_t,x_t,p)
+      CVODES.CVRhsFunc f = CreateOdeRhs(trees, calculatedVariables, nodeValues);
+
+      var calcSens = jac != null;
+      // the Jacobian ∂f/∂y
+      CVODES.CVDlsJacFunc jacF = CreateJac(trees, calculatedVariables, nodeValues);
+
+      // the RHS for the forward sensitivities (∂f/∂y)s_i(t) + ∂f/∂p_i
+      CVODES.CVSensRhsFn sensF = CreateSensitivityRhs(trees, calculatedVariables, nodeValues);
+
+      // setup solver
+      int numberOfEquations = trees.Length;
+      IntPtr y = IntPtr.Zero;
+      IntPtr cvode_mem = IntPtr.Zero;
+      IntPtr A = IntPtr.Zero;
+      IntPtr yS0 = IntPtr.Zero;
+      IntPtr linearSolver = IntPtr.Zero;
+      var ns = nodeValues.ParameterCount; // number of parameters
+
+      try {
+        y = CVODES.N_VNew_Serial(numberOfEquations);
+        // init y to current values of variables
+        // y must be initialized before calling CVodeInit
+        for (int i = 0; i < calculatedVariables.Length; i++) {
+          CVODES.NV_Set_Ith_S(y, i, nodeValues.GetVariableValue(calculatedVariables[i]).Item1);
+        }
+
+        cvode_mem = CVODES.CVodeCreate(CVODES.MultistepMethod.CV_ADAMS, CVODES.NonlinearSolverIteration.CV_FUNCTIONAL);
+
+        var flag = CVODES.CVodeInit(cvode_mem, f, rows.First(), y);
+        Assert(CVODES.CV_SUCCESS == flag);
+
+        flag = CVODES.CVodeSetErrHandlerFn(cvode_mem, errorFunction, IntPtr.Zero);
+        Assert(CVODES.CV_SUCCESS == flag);
+        double relTol = 1.0e-2;
+        double absTol = 1.0;
+        flag = CVODES.CVodeSStolerances(cvode_mem, relTol, absTol);  // TODO: probably need to adjust absTol per variable
+        Assert(CVODES.CV_SUCCESS == flag);
+
+        A = CVODES.SUNDenseMatrix(numberOfEquations, numberOfEquations);
+        Assert(A != IntPtr.Zero);
+
+        linearSolver = CVODES.SUNDenseLinearSolver(y, A);
+        Assert(linearSolver != IntPtr.Zero);
+
+        flag = CVODES.CVDlsSetLinearSolver(cvode_mem, linearSolver, A);
+        Assert(CVODES.CV_SUCCESS == flag);
+
+        flag = CVODES.CVDlsSetJacFn(cvode_mem, jacF);
+        Assert(CVODES.CV_SUCCESS == flag);
+
+        if (calcSens) {
+
+          yS0 = CVODES.N_VCloneVectorArray_Serial(ns, y); // clone the output vector for each parameter
+          unsafe {
+            // set to initial sensitivities supplied by caller
+            for (int pIdx = 0; pIdx < ns; pIdx++) {
+              var yS0_i = *((IntPtr*)yS0.ToPointer() + pIdx);
+              for (var varIdx = 0; varIdx < calculatedVariables.Length; varIdx++) {
+                CVODES.NV_Set_Ith_S(yS0_i, varIdx, nodeValues.GetVariableValue(calculatedVariables[varIdx]).Item2[pIdx]); // TODO: perf
+              }
+            }
+          }
+
+          flag = CVODES.CVodeSensInit(cvode_mem, ns, CVODES.CV_SIMULTANEOUS, sensF, yS0);
+          Assert(CVODES.CV_SUCCESS == flag);
+
+          flag = CVODES.CVodeSensEEtolerances(cvode_mem);
+          Assert(CVODES.CV_SUCCESS == flag);
+        }
+        // integrate
+        int outputIdx = calculatedVariables.Length; // values at t0 do not need to be set.
+        foreach (var tout in rows.Skip(1)) {
+          double tret = 0;
+          flag = CVODES.CVode(cvode_mem, tout, y, ref tret, CVODES.CV_NORMAL);
+          if (flag == CVODES.CV_SUCCESS) {
+            // Assert(1.0 == tout);
+            if (calcSens) {
+              // get sensitivities
+              flag = CVODES.CVodeGetSens(cvode_mem, ref tret, yS0);
+              Assert(CVODES.CV_SUCCESS == flag);
+            }
+            // update variableValues based on integration results
+            for (int varIdx = 0; varIdx < calculatedVariables.Length; varIdx++) {
+              var yi = CVODES.NV_Get_Ith_S(y, varIdx);
+              fi[outputIdx] = yi;
+              if (calcSens) {
+                // var gArr = new double[ns];
+                for (var pIdx = 0; pIdx < ns; pIdx++) {
+                  unsafe {
+                    var yS0_pi = *((IntPtr*)yS0.ToPointer() + pIdx);
+                    jac[outputIdx, pIdx] = CVODES.NV_Get_Ith_S(yS0_pi, varIdx);
+                  }
+                }
+              }
+              outputIdx++;
+            }
+
+          } else {
+            // fill up remaining values
+            while (outputIdx < fi.Length) {
+              fi[outputIdx] = fi[outputIdx - calculatedVariables.Length];
+              if (calcSens) {
+                for (var pIdx = 0; pIdx < ns; pIdx++) {
+                  jac[outputIdx, pIdx] = jac[outputIdx - calculatedVariables.Length, pIdx];
+                }
+              }
+              outputIdx++;
+            }
+            return;
+          }
+        }
+
+        // cleanup all allocated objects
+      } finally {
+        if (y != IntPtr.Zero) CVODES.N_VDestroy_Serial(y);
+        if (cvode_mem != IntPtr.Zero) CVODES.CVodeFree(ref cvode_mem);
+        if (linearSolver != IntPtr.Zero) CVODES.SUNLinSolFree(linearSolver);
+        if (A != IntPtr.Zero) CVODES.SUNMatDestroy(A);
+        if (yS0 != IntPtr.Zero) CVODES.N_VDestroyVectorArray_Serial(yS0, ns);
+      }
+    }
+
     private static void errorFunction(int errorCode, IntPtr module, IntPtr function, IntPtr msg, IntPtr ehdata) {
       var moduleStr = Marshal.PtrToStringAnsi(module);
@@ -1188 +1358 @@
       var msgStr = Marshal.PtrToStringAnsi(msg);
       string type = errorCode == 0 ? "Warning" : "Error";
-      throw new InvalidProgramException($"{type}: {msgStr} Module: {moduleStr} Function: {functionStr}");
+      // throw new InvalidProgramException($"{type}: {msgStr} Module: {moduleStr} Function: {functionStr}");
     }
 
@@ -1257 +1427 @@
           InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues, out double z, out Vector dz);
           for (int j = 0; j < calculatedVariables.Length; j++) {
-            CVODES.SUNDenseMatrix_Set(Jac, i, j, dz[j]);
+            CVODES.SUNDenseMatrix_Set(Jac, i, j, dz[j]);  //TODO: must set as in SensitivityRhs!
           }
         }
@@ -1375 +1545 @@
     // TODO: use an existing interpreter implementation instead
     private static double InterpretRec(ISymbolicExpressionTreeNode node, NodeValueLookup nodeValues) {
-      if (node is ConstantTreeNode) {
-        return ((ConstantTreeNode)node).Value;
-      } else if (node is VariableTreeNode) {
-        return nodeValues.NodeValue(node);
+      if (node is ConstantTreeNode constTreeNode) {
+        return nodeValues.ConstantNodeValue(constTreeNode);
+      } else if (node is VariableTreeNode varTreeNode) {
+        return nodeValues.VariableNodeValue(varTreeNode);
       } else if (node.Symbol is Addition) {
         var f = InterpretRec(node.GetSubtree(0), nodeValues);
@@ -1476 +1646 @@
       double f, g;
       Vector df, dg;
-      if (node.Symbol is Constant || node.Symbol is Variable) {
-        z = nodeValues.NodeValue(node);
-        dz = Vector.CreateNew(nodeValues.NodeGradient(node)); // original gradient vectors are never changed by evaluation
+      if (node is ConstantTreeNode constTreeNode) {
+        var val = nodeValues.ConstantNodeValueAndGradient(constTreeNode);
+        z = val.Item1;
+        dz = val.Item2;
+      } else if (node is VariableTreeNode varTreeNode) {
+        var val = nodeValues.VariableNodeValueAndGradient(varTreeNode);
+        z = val.Item1;
+        dz = val.Item2;
       } else if (node.Symbol is Addition) {
         InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
@@ -1775 +1950 @@
     }
 
-    private static bool IsConstantNode(ISymbolicExpressionTreeNode n) {
-      return n is ConstantTreeNode;
-    }
-    private static double GetConstantValue(ISymbolicExpressionTreeNode n) {
-      return ((ConstantTreeNode)n).Value;
-    }
-    private static bool IsLatentVariableNode(ISymbolicExpressionTreeNode n) {
-      return n.Symbol.Name[0] == 'λ';
-    }
-    private static bool IsVariableNode(ISymbolicExpressionTreeNode n) {
-      return (n.SubtreeCount == 0) && !IsConstantNode(n) && !IsLatentVariableNode(n);
-    }
-    private static string GetVariableName(ISymbolicExpressionTreeNode n) {
-      return ((VariableTreeNode)n).VariableName;
-    }
+    // private static bool IsLatentVariableNode(ISymbolicExpressionTreeNode n) {
+    //   return n.Symbol.Name[0] == 'λ';
+    // }
 
     private void UpdateTargetVariables() {
@@ -1893 +2056 @@
       // configure initialization of variables
       var varSy = (Variable)grammar.GetSymbol("Variable");
-      // fix variable weights to 1.0
-      varSy.WeightMu = 1.0;
-      varSy.WeightSigma = 0.0;
+      // init variables to a small value and allow manipulation
+      varSy.WeightMu = 0.0;
+      varSy.WeightSigma = 1e-1;
       varSy.WeightManipulatorMu = 0.0;
-      varSy.WeightManipulatorSigma = 0.0;
-      varSy.MultiplicativeWeightManipulatorSigma = 0.0;
+      varSy.WeightManipulatorSigma = 1.0;
+      varSy.MultiplicativeWeightManipulatorSigma = 1.0;
 
       foreach (var f in FunctionSet) {
@@ -1975 +2138 @@
     public class NodeValueLookup {
       private readonly Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>> node2val = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
-      private readonly Dictionary<string, List<ISymbolicExpressionTreeNode>> name2nodes = new Dictionary<string, List<ISymbolicExpressionTreeNode>>();
-      private readonly ConstantTreeNode[] constantNodes;
+      private readonly ISymbolicExpressionTreeNode[] leafNodes;
       private readonly Vector[] constantGradientVectors;
-
-
-      public double NodeValue(ISymbolicExpressionTreeNode node) => node2val[node].Item1;
-      public Vector NodeGradient(ISymbolicExpressionTreeNode node) => node2val[node].Item2;
+      private readonly Dictionary<string, Tuple<double, Vector>> variableValues = new Dictionary<string, Tuple<double, Vector>>();
+
+      // accessors for current values of constant and variable nodes. For variable nodes we also need to account for the variable weight
+      public double ConstantNodeValue(ConstantTreeNode node) => node2val[node].Item1;
+      public Tuple<double, Vector> ConstantNodeValueAndGradient(ConstantTreeNode node) { var v = node2val[node]; return Tuple.Create(v.Item1, Vector.CreateNew(v.Item2)); }
+      public double VariableNodeValue(VariableTreeNode node) => variableValues[node.VariableName].Item1 * node2val[node].Item1;
+      public Tuple<double, Vector> VariableNodeValueAndGradient(VariableTreeNode node) {
+        // (f*g)' = (f'*g)+(g'*f)        
+        var g = node2val[node];
+        var f = variableValues[node.VariableName];
+
+        return Tuple.Create(
+          g.Item1 * f.Item1,
+          Vector.CreateNew(f.Item2).Scale(g.Item1).Add(Vector.CreateNew(g.Item2).Scale(f.Item1)));
+      }
 
       public NodeValueLookup(ISymbolicExpressionTree[] trees, bool variableGradient = false) {
-        this.constantNodes = trees.SelectMany(t => t.IterateNodesPrefix().OfType<ConstantTreeNode>()).ToArray();
+        this.leafNodes = trees.SelectMany(t => t.IterateNodesPrefix().Where(n => n.SubtreeCount==0)).ToArray();
         if (!variableGradient) {
-          constantGradientVectors = new Vector[constantNodes.Length];
-          for (int paramIdx = 0; paramIdx < constantNodes.Length; paramIdx++) {
-            constantGradientVectors[paramIdx] = Vector.CreateIndicator(length: constantNodes.Length, idx: paramIdx);
-
-            var node = constantNodes[paramIdx];
-            node2val[node] = Tuple.Create(node.Value, constantGradientVectors[paramIdx]);
-          }
-
-          foreach (var tree in trees) {
-            foreach (var node in tree.IterateNodesPrefix().Where(IsVariableNode)) {
-              var varName = GetVariableName(node);
-              if (!name2nodes.TryGetValue(varName, out List<ISymbolicExpressionTreeNode> nodes)) {
-                nodes = new List<ISymbolicExpressionTreeNode>();
-                name2nodes.Add(varName, nodes);
-              }
-              nodes.Add(node);
-              SetVariableValue(varName, 0.0);  // this value is updated in the prediction loop
-            }
-          }
-        }
-        else {
+          constantGradientVectors = new Vector[leafNodes.Length];
+          for (int paramIdx = 0; paramIdx < leafNodes.Length; paramIdx++) {
+            constantGradientVectors[paramIdx] = Vector.CreateIndicator(length: leafNodes.Length, idx: paramIdx);
+
+            var node = leafNodes[paramIdx];
+            if (node is ConstantTreeNode constTreeNode) {
+              node2val[node] = Tuple.Create(constTreeNode.Value, constantGradientVectors[paramIdx]);
+            } else if (node is VariableTreeNode varTreeNode) {
+              node2val[node] = Tuple.Create(varTreeNode.Weight, constantGradientVectors[paramIdx]);
+            } else throw new InvalidProgramException();
+          }
+        } else {
           // variable gradient means we want to calculate the gradient over the target variables instead of parameters
-          for (int paramIdx = 0; paramIdx < constantNodes.Length; paramIdx++) {
-            var node = constantNodes[paramIdx];
-            node2val[node] = Tuple.Create(node.Value, Vector.Zero);
-          }
-
-          foreach (var tree in trees) {
-            foreach (var node in tree.IterateNodesPrefix().Where(IsVariableNode)) {
-              var varName = GetVariableName(node);
-              if (!name2nodes.TryGetValue(varName, out List<ISymbolicExpressionTreeNode> nodes)) {
-                nodes = new List<ISymbolicExpressionTreeNode>();
-                name2nodes.Add(varName, nodes);
-              }
-              nodes.Add(node);
-              SetVariableValue(varName, 0.0);  // this value is updated in the prediction loop
-            }
-          }
-        }
-      }
-
-      public int ParameterCount => constantNodes.Length;
+          for (int paramIdx = 0; paramIdx < leafNodes.Length; paramIdx++) {
+            var node = leafNodes[paramIdx];
+            if (node is ConstantTreeNode constTreeNode) {
+              node2val[node] = Tuple.Create(constTreeNode.Value, Vector.Zero);
+            } else if (node is VariableTreeNode varTreeNode) {
+              node2val[node] = Tuple.Create(varTreeNode.Weight, Vector.Zero);
+            } else throw new InvalidProgramException();
+          }
+        }
+      }
+
+      public int ParameterCount => leafNodes.Length;
 
       public void SetVariableValue(string variableName, double val) {
         SetVariableValue(variableName, val, Vector.Zero);
       }
+      /// <summary>
+      /// returns the current value for variable variableName
+      /// </summary>
+      /// <param name="variableName"></param>
+      /// <returns></returns>
       public Tuple<double, Vector> GetVariableValue(string variableName) {
-        return node2val[name2nodes[variableName].First()];
-      }
+        return variableValues[variableName];
+      }
+
+      /// <summary>
+      /// sets the current value for variable variableName
+      /// </summary>
+      /// <param name="variableName"></param>
+      /// <param name="val"></param>
+      /// <param name="dVal"></param>
       public void SetVariableValue(string variableName, double val, Vector dVal) {
-        if (name2nodes.TryGetValue(variableName, out List<ISymbolicExpressionTreeNode> nodes)) {
-          nodes.ForEach(n => node2val[n] = Tuple.Create(val, dVal));
-        } else {
-          var fakeNode = new VariableTreeNode(new Variable());
-          fakeNode.Weight = 1.0;
-          fakeNode.VariableName = variableName;
-          var newNodeList = new List<ISymbolicExpressionTreeNode>();
-          newNodeList.Add(fakeNode);
-          name2nodes.Add(variableName, newNodeList);
-          node2val[fakeNode] = Tuple.Create(val, dVal);
-        }
+        variableValues[variableName] = Tuple.Create(val, dVal);
+        // if (name2nodes.TryGetValue(variableName, out List<ISymbolicExpressionTreeNode> nodes)) {
+        //   nodes.ForEach(n => node2val[n] = Tuple.Create(val, dVal));
+        // } else {
+        //   var fakeNode = new VariableTreeNode(new Variable());
+        //   fakeNode.Weight = 1.0;
+        //   fakeNode.VariableName = variableName;
+        //   var newNodeList = new List<ISymbolicExpressionTreeNode>();
+        //   newNodeList.Add(fakeNode);
+        //   name2nodes.Add(variableName, newNodeList);
+        //   node2val[fakeNode] = Tuple.Create(val, dVal);
+        // }
       }
 
       internal void UpdateParamValues(double[] x) {
         for (int i = 0; i < x.Length; i++) {
-          constantNodes[i].Value = x[i];
-          node2val[constantNodes[i]] = Tuple.Create(x[i], constantGradientVectors[i]);
+          node2val[leafNodes[i]] = Tuple.Create(x[i], constantGradientVectors[i]);
         }
       }

branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/ProblemInstanceProvider.cs

@@ -293 +293 @@
         TargetVariables = new[] { "x1", "x2", "v1", "v2" },
         InputVariables = new string[] { },
-        TrainingEpisodes = new IntRange[] { new IntRange(0, 820) },
+        TrainingEpisodes = new IntRange[] { new IntRange(0, 200) },
         TestEpisodes = new IntRange[] { },
         FileName = "double_linear_h_1_equidistant.txt",
@@ -306 +306 @@
         TargetVariables = new[] { "x1", "x2", "v1", "v2" },
         InputVariables = new string[] { },
-        TrainingEpisodes = new IntRange[] { new IntRange(0, 500) },
+        TrainingEpisodes = new IntRange[] { new IntRange(0, 150) },
         TestEpisodes = new IntRange[] { },
         FileName = "real_double_linear_h_1_equidistant.txt",
@@ -332 +332 @@
         TargetVariables = new[] { "theta1", "theta2", "omega1", "omega2" },
         InputVariables = new string[] { },
-        TrainingEpisodes = new IntRange[] { new IntRange(0, 886) },
+        TrainingEpisodes = new IntRange[] { new IntRange(0, 200) },
         TestEpisodes = new IntRange[] {new IntRange(886, 1731) },
         FileName = "real_double_pend_h_1_equidistant.txt",