Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
10/22/18 18:35:27 (6 years ago)
Author:
gkronber
Message:

#2925 preparations to switch between HL ODE solver and CVODES

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs

    r16249 r16250  
    4040namespace HeuristicLab.Problems.DynamicalSystemsModelling {
    4141
    42  
     42
    4343
    4444
     
    5050  [StorableClass]
    5151  public sealed class Problem : SingleObjectiveBasicProblem<MultiEncoding>, IRegressionProblem, IProblemInstanceConsumer<IRegressionProblemData>, IProblemInstanceExporter<IRegressionProblemData> {
    52 
    53 
    54 
    55 
    5652    #region parameter names
    5753    private const string ProblemDataParameterName = "Data";
     
    6460    private const string TrainingEpisodesParameterName = "Training episodes";
    6561    private const string OptimizeParametersForEpisodesParameterName = "Optimize parameters for episodes";
     62    private const string OdeSolverParameterName = "ODE Solver";
    6663    #endregion
    6764
     
    9592    public IFixedValueParameter<BoolValue> OptimizeParametersForEpisodesParameter {
    9693      get { return (IFixedValueParameter<BoolValue>)Parameters[OptimizeParametersForEpisodesParameterName]; }
     94    }
     95    public IConstrainedValueParameter<StringValue> OdeSolverParameter {
     96      get { return (IConstrainedValueParameter<StringValue>)Parameters[OdeSolverParameterName]; }
    9797    }
    9898    #endregion
     
    130130    public bool OptimizeParametersForEpisodes {
    131131      get { return OptimizeParametersForEpisodesParameter.Value.Value; }
     132    }
     133
     134    public string OdeSolver {
     135      get { return OdeSolverParameter.Value.Value; }
     136      set {
     137        var matchingValue = OdeSolverParameter.ValidValues.FirstOrDefault(v => v.Value == value);
     138        if (matchingValue == null) throw new ArgumentOutOfRangeException();
     139        else OdeSolverParameter.Value = matchingValue;
     140      }
    132141    }
    133142
     
    173182      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>()));
    174183      Parameters.Add(new FixedValueParameter<BoolValue>(OptimizeParametersForEpisodesParameterName, "Flag to select if parameters should be optimized globally or for each episode individually.", new BoolValue(false)));
     184
     185      var solversStr = new string[] { "HeuristicLab", "COVDES" };
     186      var solvers = new ItemSet<StringValue>(
     187        solversStr.Select(s => new StringValue(s).AsReadOnly())
     188        );
     189      Parameters.Add(new ConstrainedValueParameter<StringValue>(OdeSolverParameterName, "The solver to use for solving the initial value ODE problems", solvers));
     190
    175191      RegisterEventHandlers();
    176192      InitAllParameters();
     
    208224    }
    209225
    210     private void OptimizeForEpisodes(ISymbolicExpressionTree[] trees, IRandom random, IEnumerable<IntRange> episodes, out double[] optTheta, out double nmse) {
     226    private void OptimizeForEpisodes(
     227      ISymbolicExpressionTree[] trees,
     228      IRandom random,
     229      IEnumerable<IntRange> episodes,
     230      out double[] optTheta,
     231      out double nmse) {
    211232      var rows = episodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start)).ToArray();
    212233      var problemData = ProblemData;
     
    243264        alglib.minlbfgssetcond(state, 0.0, 0.0, 0.0, MaximumParameterOptimizationIterations);
    244265        alglib.minlbfgsoptimize(state, EvaluateObjectiveAndGradient, null,
    245           new object[] { trees, targetVars, problemData, nodeIdx, targetValues, episodes.ToArray(), NumericIntegrationSteps, latentVariables }); //TODO: create a type
     266          new object[] { trees, targetVars, problemData, nodeIdx, targetValues, episodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver }); //TODO: create a type
     267
    246268        alglib.minlbfgsresults(state, out optTheta, out report);
    247269
     
    292314      var numericIntegrationSteps = (int)((object[])obj)[6];
    293315      var latentVariables = (string[])((object[])obj)[7];
    294 
    295       var predicted = Integrate(
    296         trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
    297         problemData.Dataset,
    298         problemData.AllowedInputVariables.ToArray(),
    299         targetVariables,
    300         latentVariables,
    301         episodes,
    302         nodeIdx,
    303         x, numericIntegrationSteps).ToArray();
     316      var odeSolver = (string)((object[])obj)[8];
     317
     318
     319      Tuple<double, Vector>[][] predicted = null; // one array of predictions for each episode
     320      predicted = Integrate(
     321          trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
     322          problemData.Dataset,
     323          problemData.AllowedInputVariables.ToArray(),
     324          targetVariables,
     325          latentVariables,
     326          episodes,
     327          nodeIdx,
     328          x,
     329          odeSolver,
     330          numericIntegrationSteps).ToArray();
    304331
    305332
     
    379406                                   nodeIdx,
    380407                                   optTheta,
     408                                   OdeSolver,
    381409                                   NumericIntegrationSteps).ToArray();
    382410          trainingPredictions.Add(trainingPrediction);
     
    420448                                   nodeIdx,
    421449                                   optTheta,
     450                                   OdeSolver,
    422451                                   NumericIntegrationSteps).ToArray();
    423452        // only for actual target values
     
    444473         nodeIdx,
    445474         optTheta,
     475         OdeSolver,
    446476         NumericIntegrationSteps).ToArray();
    447477
     
    577607    private static IEnumerable<Tuple<double, Vector>[]> Integrate(
    578608      ISymbolicExpressionTree[] trees, IDataset dataset, string[] inputVariables, string[] targetVariables, string[] latentVariables, IEnumerable<IntRange> episodes,
    579       Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx, double[] parameterValues, int numericIntegrationSteps = 100) {
    580 
    581       int NUM_STEPS = numericIntegrationSteps;
    582       double h = 1.0 / NUM_STEPS;
     609      Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx, double[] parameterValues,
     610      string odeSolver, int numericIntegrationSteps = 100) {
     611
     612      // TODO: numericIntegrationSteps is only relevant for the HeuristicLab solver
    583613
    584614      foreach (var episode in episodes) {
     
    603633          variableValues.Add(latentVar, Tuple.Create(0.0, Vector.Zero)); // we don't have observations for latent variables -> assume zero as starting value
    604634        }
    605         var calculatedVariables = targetVariables.Concat(latentVariables); // TODO: must conincide with the order of trees in the encoding
     635        var calculatedVariables = targetVariables.Concat(latentVariables).ToArray(); // TODO: must conincide with the order of trees in the encoding
    606636
    607637        foreach (var t in rows.Skip(1)) {
    608           for (int step = 0; step < NUM_STEPS; step++) {
    609             var deltaValues = new Dictionary<string, Tuple<double, Vector>>();
    610             foreach (var tup in trees.Zip(calculatedVariables, Tuple.Create)) {
    611               var tree = tup.Item1;
    612               var targetVarName = tup.Item2;
    613               // skip programRoot and startSymbol
    614               var res = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), variableValues, nodeIdx, parameterValues);
    615               deltaValues.Add(targetVarName, res);
    616             }
    617 
    618             // update variableValues for next step
    619             foreach (var kvp in deltaValues) {
    620               var oldVal = variableValues[kvp.Key];
    621               variableValues[kvp.Key] = Tuple.Create(
    622                 oldVal.Item1 + h * kvp.Value.Item1,
    623                 oldVal.Item2 + h * kvp.Value.Item2
    624               );
    625             }
    626           }
     638          if (odeSolver == "HeuristicLab")
     639            IntegrateHL(trees, nodeIdx, calculatedVariables, variableValues, parameterValues, numericIntegrationSteps);
     640          else if (odeSolver == "CVODES")
     641            IntegrateCVODES(trees, nodeIdx, calculatedVariables, variableValues, parameterValues);
     642          else throw new InvalidOperationException("Unknown ODE solver " + odeSolver);
    627643
    628644          // only return the target variables for calculation of errors
     
    635651            variableValues[varName] = Tuple.Create(dataset.GetDoubleValue(varName, t), Vector.Zero);
    636652          }
     653        }
     654      }
     655    }
     656
     657    private static void IntegrateCVODES(
     658      ISymbolicExpressionTree[] trees, // f(u,p) in tree representation
     659      Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx,
     660      string[] calculatedVariables, // names of elements of u
     661      Dictionary<string, Tuple<double, Vector>> variableValues,  //  u(t): output of the system
     662      double[] parameterValues // p
     663      ) {
     664      // the RHS of the ODE
     665      CVODES.CVRhsFunc f = CreateOdeRhs(trees, nodeIdx, calculatedVariables, parameterValues);
     666      // XXX Jacobian function
     667    }
     668
     669    private static CVODES.CVRhsFunc CreateOdeRhs(
     670      ISymbolicExpressionTree[] trees,
     671      Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx,
     672      string[] calculatedVariables,
     673      double[] parameterValues) {
     674      var variableValues = new Dictionary<string, Tuple<double, Vector>>();
     675      // init variableValues with zeros
     676      for (int i = 0; i < calculatedVariables.Length; i++) {
     677        var backingArr = new double[parameterValues.Length];
     678        variableValues.Add(calculatedVariables[i], Tuple.Create(0.0, new Vector(backingArr)));
     679      }
     680      return (double t,
     681              IntPtr y, // N_Vector, current value of y (input)
     682              IntPtr ydot, // N_Vector (calculated value of y' (output)
     683              IntPtr user_data // optional user data, (unused here)
     684              ) => {
     685                // copy y array to variableValues dictionary
     686                for(int i=0;i<calculatedVariables.Length;i++) {
     687                  var yi = CVODES.NV_Get_Ith_S(y, (long)i);
     688                  variableValues[calculatedVariables[i]] = Tuple.Create(yi, Vector.Zero);
     689                }
     690                for (int i = 0; i < trees.Length; i++) {
     691                  var tree = trees[i];
     692                  var res_i = InterpretRec(tree.Root, variableValues, nodeIdx, parameterValues); // TODO: perf - we don't need AutoDiff-based gradients here
     693                  CVODES.NV_Set_Ith_S(ydot, i, res_i.Item1);
     694                }
     695                return 0;
     696              };
     697    }
     698
     699    private static void IntegrateHL(ISymbolicExpressionTree[] trees, Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx,
     700      string[] calculatedVariables,
     701      Dictionary<string, Tuple<double, Vector>> variableValues,
     702      double[] parameterValues,
     703      int numericIntegrationSteps) {
     704      double h = 1.0 / numericIntegrationSteps;
     705      for (int step = 0; step < numericIntegrationSteps; step++) {
     706        var deltaValues = new Dictionary<string, Tuple<double, Vector>>();
     707        foreach (var tup in trees.Zip(calculatedVariables, Tuple.Create)) {
     708          var tree = tup.Item1;
     709          var targetVarName = tup.Item2;
     710          // skip programRoot and startSymbol
     711          var res = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), variableValues, nodeIdx, parameterValues);
     712          deltaValues.Add(targetVarName, res);
     713        }
     714
     715        // update variableValues for next step, trapecoid integration
     716        foreach (var kvp in deltaValues) {
     717          var oldVal = variableValues[kvp.Key];
     718          variableValues[kvp.Key] = Tuple.Create(
     719            oldVal.Item1 + h * kvp.Value.Item1,
     720            oldVal.Item2 + h * kvp.Value.Item2
     721          );
    637722        }
    638723      }
Note: See TracChangeset for help on using the changeset viewer.