Changeset 16250
- Timestamp:
- 10/22/18 18:35:27 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs
r16249 r16250 40 40 namespace HeuristicLab.Problems.DynamicalSystemsModelling { 41 41 42 42 43 43 44 44 … … 50 50 [StorableClass] 51 51 public sealed class Problem : SingleObjectiveBasicProblem<MultiEncoding>, IRegressionProblem, IProblemInstanceConsumer<IRegressionProblemData>, IProblemInstanceExporter<IRegressionProblemData> { 52 53 54 55 56 52 #region parameter names 57 53 private const string ProblemDataParameterName = "Data"; … … 64 60 private const string TrainingEpisodesParameterName = "Training episodes"; 65 61 private const string OptimizeParametersForEpisodesParameterName = "Optimize parameters for episodes"; 62 private const string OdeSolverParameterName = "ODE Solver"; 66 63 #endregion 67 64 … … 95 92 public IFixedValueParameter<BoolValue> OptimizeParametersForEpisodesParameter { 96 93 get { return (IFixedValueParameter<BoolValue>)Parameters[OptimizeParametersForEpisodesParameterName]; } 94 } 95 public IConstrainedValueParameter<StringValue> OdeSolverParameter { 96 get { return (IConstrainedValueParameter<StringValue>)Parameters[OdeSolverParameterName]; } 97 97 } 98 98 #endregion … … 130 130 public bool OptimizeParametersForEpisodes { 131 131 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 } 132 141 } 133 142 … … 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>())); 174 183 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 175 191 RegisterEventHandlers(); 176 192 InitAllParameters(); … … 208 224 } 209 225 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) { 211 232 var rows = episodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start)).ToArray(); 212 233 var problemData = ProblemData; … … 243 264 alglib.minlbfgssetcond(state, 0.0, 0.0, 0.0, MaximumParameterOptimizationIterations); 244 265 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 246 268 alglib.minlbfgsresults(state, out optTheta, out report); 247 269 … … 292 314 var numericIntegrationSteps = (int)((object[])obj)[6]; 293 315 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(); 304 331 305 332 … … 379 406 nodeIdx, 380 407 optTheta, 408 OdeSolver, 381 409 NumericIntegrationSteps).ToArray(); 382 410 trainingPredictions.Add(trainingPrediction); … … 420 448 nodeIdx, 421 449 optTheta, 450 OdeSolver, 422 451 NumericIntegrationSteps).ToArray(); 423 452 // only for actual target values … … 444 473 nodeIdx, 445 474 optTheta, 475 OdeSolver, 446 476 NumericIntegrationSteps).ToArray(); 447 477 … … 577 607 private static IEnumerable<Tuple<double, Vector>[]> Integrate( 578 608 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 583 613 584 614 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 604 634 } 605 var calculatedVariables = targetVariables.Concat(latentVariables) ; // TODO: must conincide with the order of trees in the encoding635 var calculatedVariables = targetVariables.Concat(latentVariables).ToArray(); // TODO: must conincide with the order of trees in the encoding 606 636 607 637 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); 627 643 628 644 // only return the target variables for calculation of errors … … 635 651 variableValues[varName] = Tuple.Create(dataset.GetDoubleValue(varName, t), Vector.Zero); 636 652 } 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 ); 637 722 } 638 723 }
Note: See TracChangeset
for help on using the changeset viewer.