Changeset 16250 for branches/2925_AutoDiffForDynamicalModels

Timestamp:

10/22/18 18:35:27 (6 years ago)

Author:

gkronber

Message:

#2925 preparations to switch between HL ODE solver and CVODES

File:

• branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs (modified) (15 diffs)

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

@@ -40 +40 @@
 namespace HeuristicLab.Problems.DynamicalSystemsModelling {
 
-  
+
 
 
@@ -50 +50 @@
   [StorableClass]
   public sealed class Problem : SingleObjectiveBasicProblem<MultiEncoding>, IRegressionProblem, IProblemInstanceConsumer<IRegressionProblemData>, IProblemInstanceExporter<IRegressionProblemData> {
-
-
-
-
     #region parameter names
     private const string ProblemDataParameterName = "Data";
@@ -64 +60 @@
     private const string TrainingEpisodesParameterName = "Training episodes";
     private const string OptimizeParametersForEpisodesParameterName = "Optimize parameters for episodes";
+    private const string OdeSolverParameterName = "ODE Solver";
     #endregion
 
@@ -95 +92 @@
     public IFixedValueParameter<BoolValue> OptimizeParametersForEpisodesParameter {
       get { return (IFixedValueParameter<BoolValue>)Parameters[OptimizeParametersForEpisodesParameterName]; }
+    }
+    public IConstrainedValueParameter<StringValue> OdeSolverParameter {
+      get { return (IConstrainedValueParameter<StringValue>)Parameters[OdeSolverParameterName]; }
     }
     #endregion
@@ -130 +130 @@
     public bool OptimizeParametersForEpisodes {
       get { return OptimizeParametersForEpisodesParameter.Value.Value; }
+    }
+
+    public string OdeSolver {
+      get { return OdeSolverParameter.Value.Value; }
+      set {
+        var matchingValue = OdeSolverParameter.ValidValues.FirstOrDefault(v => v.Value == value);
+        if (matchingValue == null) throw new ArgumentOutOfRangeException();
+        else OdeSolverParameter.Value = matchingValue;
+      }
     }
 
@@ -173 +182 @@
       Parameters.Add(new ValueParameter<ItemList<IntRange>>(TrainingEpisodesParameterName, "A list of ranges that should be used for training, each range represents an independent episode. This overrides the TrainingSet parameter in ProblemData.", new ItemList<IntRange>()));
       Parameters.Add(new FixedValueParameter<BoolValue>(OptimizeParametersForEpisodesParameterName, "Flag to select if parameters should be optimized globally or for each episode individually.", new BoolValue(false)));
+
+      var solversStr = new string[] { "HeuristicLab", "COVDES" };
+      var solvers = new ItemSet<StringValue>(
+        solversStr.Select(s => new StringValue(s).AsReadOnly())
+        );
+      Parameters.Add(new ConstrainedValueParameter<StringValue>(OdeSolverParameterName, "The solver to use for solving the initial value ODE problems", solvers));
+
       RegisterEventHandlers();
       InitAllParameters();
@@ -208 +224 @@
     }
 
-    private void OptimizeForEpisodes(ISymbolicExpressionTree[] trees, IRandom random, IEnumerable<IntRange> episodes, out double[] optTheta, out double nmse) {
+    private void OptimizeForEpisodes(
+      ISymbolicExpressionTree[] trees,
+      IRandom random,
+      IEnumerable<IntRange> episodes,
+      out double[] optTheta,
+      out double nmse) {
       var rows = episodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start)).ToArray();
       var problemData = ProblemData;
@@ -243 +264 @@
         alglib.minlbfgssetcond(state, 0.0, 0.0, 0.0, MaximumParameterOptimizationIterations);
         alglib.minlbfgsoptimize(state, EvaluateObjectiveAndGradient, null,
-          new object[] { trees, targetVars, problemData, nodeIdx, targetValues, episodes.ToArray(), NumericIntegrationSteps, latentVariables }); //TODO: create a type
+          new object[] { trees, targetVars, problemData, nodeIdx, targetValues, episodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver }); //TODO: create a type
+
         alglib.minlbfgsresults(state, out optTheta, out report);
 
@@ -292 +314 @@
       var numericIntegrationSteps = (int)((object[])obj)[6];
       var latentVariables = (string[])((object[])obj)[7];
-
-      var predicted = Integrate(
-        trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
-        problemData.Dataset,
-        problemData.AllowedInputVariables.ToArray(),
-        targetVariables,
-        latentVariables,
-        episodes,
-        nodeIdx,
-        x, numericIntegrationSteps).ToArray();
+      var odeSolver = (string)((object[])obj)[8];
+
+
+      Tuple<double, Vector>[][] predicted = null; // one array of predictions for each episode
+      predicted = Integrate(
+          trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
+          problemData.Dataset,
+          problemData.AllowedInputVariables.ToArray(),
+          targetVariables,
+          latentVariables,
+          episodes,
+          nodeIdx,
+          x,
+          odeSolver,
+          numericIntegrationSteps).ToArray();
 
 
@@ -379 +406 @@
                                    nodeIdx,
                                    optTheta,
+                                   OdeSolver,
                                    NumericIntegrationSteps).ToArray();
           trainingPredictions.Add(trainingPrediction);
@@ -420 +448 @@
                                    nodeIdx,
                                    optTheta,
+                                   OdeSolver,
                                    NumericIntegrationSteps).ToArray();
         // only for actual target values
@@ -444 +473 @@
          nodeIdx,
          optTheta,
+         OdeSolver,
          NumericIntegrationSteps).ToArray();
 
@@ -577 +607 @@
     private static IEnumerable<Tuple<double, Vector>[]> Integrate(
       ISymbolicExpressionTree[] trees, IDataset dataset, string[] inputVariables, string[] targetVariables, string[] latentVariables, IEnumerable<IntRange> episodes,
-      Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx, double[] parameterValues, int numericIntegrationSteps = 100) {
-
-      int NUM_STEPS = numericIntegrationSteps;
-      double h = 1.0 / NUM_STEPS;
+      Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx, double[] parameterValues,
+      string odeSolver, int numericIntegrationSteps = 100) {
+
+      // TODO: numericIntegrationSteps is only relevant for the HeuristicLab solver
 
       foreach (var episode in episodes) {
@@ -603 +633 @@
           variableValues.Add(latentVar, Tuple.Create(0.0, Vector.Zero)); // we don't have observations for latent variables -> assume zero as starting value
         }
-        var calculatedVariables = targetVariables.Concat(latentVariables); // TODO: must conincide with the order of trees in the encoding
+        var calculatedVariables = targetVariables.Concat(latentVariables).ToArray(); // TODO: must conincide with the order of trees in the encoding
 
         foreach (var t in rows.Skip(1)) {
-          for (int step = 0; step < NUM_STEPS; step++) {
-            var deltaValues = new Dictionary<string, Tuple<double, Vector>>();
-            foreach (var tup in trees.Zip(calculatedVariables, Tuple.Create)) {
-              var tree = tup.Item1;
-              var targetVarName = tup.Item2;
-              // skip programRoot and startSymbol
-              var res = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), variableValues, nodeIdx, parameterValues);
-              deltaValues.Add(targetVarName, res);
-            }
-
-            // update variableValues for next step
-            foreach (var kvp in deltaValues) {
-              var oldVal = variableValues[kvp.Key];
-              variableValues[kvp.Key] = Tuple.Create(
-                oldVal.Item1 + h * kvp.Value.Item1,
-                oldVal.Item2 + h * kvp.Value.Item2
-              );
-            }
-          }
+          if (odeSolver == "HeuristicLab")
+            IntegrateHL(trees, nodeIdx, calculatedVariables, variableValues, parameterValues, numericIntegrationSteps);
+          else if (odeSolver == "CVODES")
+            IntegrateCVODES(trees, nodeIdx, calculatedVariables, variableValues, parameterValues);
+          else throw new InvalidOperationException("Unknown ODE solver " + odeSolver);
 
           // only return the target variables for calculation of errors
@@ -635 +651 @@
             variableValues[varName] = Tuple.Create(dataset.GetDoubleValue(varName, t), Vector.Zero);
           }
+        }
+      }
+    }
+
+    private static void IntegrateCVODES(
+      ISymbolicExpressionTree[] trees, // f(u,p) in tree representation
+      Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx,
+      string[] calculatedVariables, // names of elements of u
+      Dictionary<string, Tuple<double, Vector>> variableValues,  //  u(t): output of the system
+      double[] parameterValues // p
+      ) {
+      // the RHS of the ODE
+      CVODES.CVRhsFunc f = CreateOdeRhs(trees, nodeIdx, calculatedVariables, parameterValues);
+      // XXX Jacobian function
+    }
+
+    private static CVODES.CVRhsFunc CreateOdeRhs(
+      ISymbolicExpressionTree[] trees,
+      Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx,
+      string[] calculatedVariables,
+      double[] parameterValues) {
+      var variableValues = new Dictionary<string, Tuple<double, Vector>>();
+      // init variableValues with zeros
+      for (int i = 0; i < calculatedVariables.Length; i++) {
+        var backingArr = new double[parameterValues.Length];
+        variableValues.Add(calculatedVariables[i], Tuple.Create(0.0, new Vector(backingArr)));
+      }
+      return (double t,
+              IntPtr y, // N_Vector, current value of y (input)
+              IntPtr ydot, // N_Vector (calculated value of y' (output)
+              IntPtr user_data // optional user data, (unused here)
+              ) => {
+                // copy y array to variableValues dictionary
+                for(int i=0;i<calculatedVariables.Length;i++) {
+                  var yi = CVODES.NV_Get_Ith_S(y, (long)i);
+                  variableValues[calculatedVariables[i]] = Tuple.Create(yi, Vector.Zero);
+                }
+                for (int i = 0; i < trees.Length; i++) {
+                  var tree = trees[i];
+                  var res_i = InterpretRec(tree.Root, variableValues, nodeIdx, parameterValues); // TODO: perf - we don't need AutoDiff-based gradients here
+                  CVODES.NV_Set_Ith_S(ydot, i, res_i.Item1);
+                }
+                return 0;
+              };
+    }
+
+    private static void IntegrateHL(ISymbolicExpressionTree[] trees, Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx,
+      string[] calculatedVariables,
+      Dictionary<string, Tuple<double, Vector>> variableValues,
+      double[] parameterValues,
+      int numericIntegrationSteps) {
+      double h = 1.0 / numericIntegrationSteps;
+      for (int step = 0; step < numericIntegrationSteps; step++) {
+        var deltaValues = new Dictionary<string, Tuple<double, Vector>>();
+        foreach (var tup in trees.Zip(calculatedVariables, Tuple.Create)) {
+          var tree = tup.Item1;
+          var targetVarName = tup.Item2;
+          // skip programRoot and startSymbol
+          var res = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), variableValues, nodeIdx, parameterValues);
+          deltaValues.Add(targetVarName, res);
+        }
+
+        // update variableValues for next step, trapecoid integration
+        foreach (var kvp in deltaValues) {
+          var oldVal = variableValues[kvp.Key];
+          variableValues[kvp.Key] = Tuple.Create(
+            oldVal.Item1 + h * kvp.Value.Item1,
+            oldVal.Item2 + h * kvp.Value.Item2
+          );
         }
       }