source: branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs @ 16329

Last change on this file since 16329 was 16329, checked in by gkronber, 9 months ago

#2925: made several extensions in relation to blood glucose prediction

  • added sqr function,
  • added outputs for latent variables (linechart and model),
  • added optimization of initial values for latent variables (for each episode separately)
  • TODO: test with CVODES (so far only our own integration scheme has been tested)
File size: 61.0 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2018 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22using System;
23using System.Collections.Generic;
24using System.Diagnostics;
25using System.Linq;
26using HeuristicLab.Analysis;
27using HeuristicLab.Collections;
28using HeuristicLab.Common;
29using HeuristicLab.Core;
30using HeuristicLab.Data;
31using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
32using HeuristicLab.Optimization;
33using HeuristicLab.Parameters;
34using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
35using HeuristicLab.Problems.DataAnalysis;
36using HeuristicLab.Problems.DataAnalysis.Symbolic;
37using HeuristicLab.Problems.Instances;
38using Variable = HeuristicLab.Problems.DataAnalysis.Symbolic.Variable;
39
40namespace HeuristicLab.Problems.DynamicalSystemsModelling {
41
42
43
44
45  // Eine weitere Möglichkeit ist spline-smoothing der Daten (über Zeit) um damit für jede Zielvariable
46  // einen bereinigten (Rauschen) Wert und die Ableitung dy/dt für alle Beobachtungen zu bekommen
47  // danach kann man die Regression direkt für dy/dt anwenden (mit bereinigten y als inputs)
48  [Item("Dynamical Systems Modelling Problem", "TODO")]
49  [Creatable(CreatableAttribute.Categories.GeneticProgrammingProblems, Priority = 900)]
50  [StorableClass]
51  public sealed class Problem : SingleObjectiveBasicProblem<MultiEncoding>, IRegressionProblem, IProblemInstanceConsumer<IRegressionProblemData>, IProblemInstanceExporter<IRegressionProblemData> {
52    #region parameter names
53    private const string ProblemDataParameterName = "Data";
54    private const string TargetVariablesParameterName = "Target variables";
55    private const string FunctionSetParameterName = "Function set";
56    private const string MaximumLengthParameterName = "Size limit";
57    private const string MaximumParameterOptimizationIterationsParameterName = "Max. parameter optimization iterations";
58    private const string NumberOfLatentVariablesParameterName = "Number of latent variables";
59    private const string NumericIntegrationStepsParameterName = "Steps for numeric integration";
60    private const string TrainingEpisodesParameterName = "Training episodes";
61    private const string OptimizeParametersForEpisodesParameterName = "Optimize parameters for episodes";
62    private const string OdeSolverParameterName = "ODE Solver";
63    #endregion
64
65    #region Parameter Properties
66    IParameter IDataAnalysisProblem.ProblemDataParameter { get { return ProblemDataParameter; } }
67
68    public IValueParameter<IRegressionProblemData> ProblemDataParameter {
69      get { return (IValueParameter<IRegressionProblemData>)Parameters[ProblemDataParameterName]; }
70    }
71    public IValueParameter<ReadOnlyCheckedItemList<StringValue>> TargetVariablesParameter {
72      get { return (IValueParameter<ReadOnlyCheckedItemList<StringValue>>)Parameters[TargetVariablesParameterName]; }
73    }
74    public IValueParameter<ReadOnlyCheckedItemList<StringValue>> FunctionSetParameter {
75      get { return (IValueParameter<ReadOnlyCheckedItemList<StringValue>>)Parameters[FunctionSetParameterName]; }
76    }
77    public IFixedValueParameter<IntValue> MaximumLengthParameter {
78      get { return (IFixedValueParameter<IntValue>)Parameters[MaximumLengthParameterName]; }
79    }
80    public IFixedValueParameter<IntValue> MaximumParameterOptimizationIterationsParameter {
81      get { return (IFixedValueParameter<IntValue>)Parameters[MaximumParameterOptimizationIterationsParameterName]; }
82    }
83    public IFixedValueParameter<IntValue> NumberOfLatentVariablesParameter {
84      get { return (IFixedValueParameter<IntValue>)Parameters[NumberOfLatentVariablesParameterName]; }
85    }
86    public IFixedValueParameter<IntValue> NumericIntegrationStepsParameter {
87      get { return (IFixedValueParameter<IntValue>)Parameters[NumericIntegrationStepsParameterName]; }
88    }
89    public IValueParameter<ItemList<IntRange>> TrainingEpisodesParameter {
90      get { return (IValueParameter<ItemList<IntRange>>)Parameters[TrainingEpisodesParameterName]; }
91    }
92    public IFixedValueParameter<BoolValue> OptimizeParametersForEpisodesParameter {
93      get { return (IFixedValueParameter<BoolValue>)Parameters[OptimizeParametersForEpisodesParameterName]; }
94    }
95    public IConstrainedValueParameter<StringValue> OdeSolverParameter {
96      get { return (IConstrainedValueParameter<StringValue>)Parameters[OdeSolverParameterName]; }
97    }
98    #endregion
99
100    #region Properties
101    public IRegressionProblemData ProblemData {
102      get { return ProblemDataParameter.Value; }
103      set { ProblemDataParameter.Value = value; }
104    }
105    IDataAnalysisProblemData IDataAnalysisProblem.ProblemData { get { return ProblemData; } }
106
107    public ReadOnlyCheckedItemList<StringValue> TargetVariables {
108      get { return TargetVariablesParameter.Value; }
109    }
110
111    public ReadOnlyCheckedItemList<StringValue> FunctionSet {
112      get { return FunctionSetParameter.Value; }
113    }
114
115    public int MaximumLength {
116      get { return MaximumLengthParameter.Value.Value; }
117    }
118    public int MaximumParameterOptimizationIterations {
119      get { return MaximumParameterOptimizationIterationsParameter.Value.Value; }
120    }
121    public int NumberOfLatentVariables {
122      get { return NumberOfLatentVariablesParameter.Value.Value; }
123    }
124    public int NumericIntegrationSteps {
125      get { return NumericIntegrationStepsParameter.Value.Value; }
126    }
127    public IEnumerable<IntRange> TrainingEpisodes {
128      get { return TrainingEpisodesParameter.Value; }
129    }
130    public bool OptimizeParametersForEpisodes {
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      }
141    }
142
143    #endregion
144
145    public event EventHandler ProblemDataChanged;
146
147    public override bool Maximization {
148      get { return false; } // we minimize NMSE
149    }
150
151    #region item cloning and persistence
152    // persistence
153    [StorableConstructor]
154    private Problem(bool deserializing) : base(deserializing) { }
155    [StorableHook(HookType.AfterDeserialization)]
156    private void AfterDeserialization() {
157      if (!Parameters.ContainsKey(OptimizeParametersForEpisodesParameterName)) {
158        Parameters.Add(new FixedValueParameter<BoolValue>(OptimizeParametersForEpisodesParameterName, "Flag to select if parameters should be optimized globally or for each episode individually.", new BoolValue(false)));
159      }
160      RegisterEventHandlers();
161    }
162
163    // cloning
164    private Problem(Problem original, Cloner cloner)
165      : base(original, cloner) {
166      RegisterEventHandlers();
167    }
168    public override IDeepCloneable Clone(Cloner cloner) { return new Problem(this, cloner); }
169    #endregion
170
171    public Problem()
172      : base() {
173      var targetVariables = new CheckedItemList<StringValue>().AsReadOnly(); // HACK: it would be better to provide a new class derived from IDataAnalysisProblem
174      var functions = CreateFunctionSet();
175      Parameters.Add(new ValueParameter<IRegressionProblemData>(ProblemDataParameterName, "The data captured from the dynamical system. Use CSV import functionality to import data.", new RegressionProblemData()));
176      Parameters.Add(new ValueParameter<ReadOnlyCheckedItemList<StringValue>>(TargetVariablesParameterName, "Target variables (overrides setting in ProblemData)", targetVariables));
177      Parameters.Add(new ValueParameter<ReadOnlyCheckedItemList<StringValue>>(FunctionSetParameterName, "The list of allowed functions", functions));
178      Parameters.Add(new FixedValueParameter<IntValue>(MaximumLengthParameterName, "The maximally allowed length of each expression. Set to a small value (5 - 25). Default = 10", new IntValue(10)));
179      Parameters.Add(new FixedValueParameter<IntValue>(MaximumParameterOptimizationIterationsParameterName, "The maximum number of iterations for optimization of parameters (using L-BFGS). More iterations makes the algorithm slower, fewer iterations might prevent convergence in the optimization scheme. Default = 100", new IntValue(100)));
180      Parameters.Add(new FixedValueParameter<IntValue>(NumberOfLatentVariablesParameterName, "Latent variables (unobserved variables) allow us to produce expressions which are integrated up and can be used in other expressions. They are handled similarly to target variables in forward simulation / integration. The difference to target variables is that there are no data to which the calculated values of latent variables are compared. Set to a small value (0 .. 5) as necessary (default = 0)", new IntValue(0)));
181      Parameters.Add(new FixedValueParameter<IntValue>(NumericIntegrationStepsParameterName, "Number of steps in the numeric integration that are taken from one row to the next (set to 1 to 100). More steps makes the algorithm slower, less steps worsens the accuracy of the numeric integration scheme.", new IntValue(10)));
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>()));
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", "CVODES" };
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, solvers.First()));
190
191      RegisterEventHandlers();
192      InitAllParameters();
193
194      // TODO: do not clear selection of target variables when the input variables are changed (keep selected target variables)
195      // TODO: UI hangs when selecting / deselecting input variables because the encoding is updated on each item
196      // TODO: use training range as default training episode
197      // TODO: write back optimized parameters to solution?
198
199    }
200
201    public override double Evaluate(Individual individual, IRandom random) {
202      var trees = individual.Values.Select(v => v.Value).OfType<ISymbolicExpressionTree>().ToArray(); // extract all trees from individual
203      // write back optimized parameters to tree nodes instead of the separate OptTheta variable
204      // retreive optimized parameters from nodes?
205
206      if (OptimizeParametersForEpisodes) {
207        int eIdx = 0;
208        double totalNMSE = 0.0;
209        int totalSize = 0;
210        foreach (var episode in TrainingEpisodes) {
211          double[] optTheta;
212          double nmse;
213          OptimizeForEpisodes(trees, random, new[] { episode }, out optTheta, out nmse);
214          individual["OptTheta_" + eIdx] = new DoubleArray(optTheta); // write back optimized parameters so that we can use them in the Analysis method
215          eIdx++;
216          totalNMSE += nmse * episode.Size;
217          totalSize += episode.Size;
218        }
219        return totalNMSE / totalSize;
220      } else {
221        double[] optTheta;
222        double nmse;
223        OptimizeForEpisodes(trees, random, TrainingEpisodes, out optTheta, out nmse);
224        individual["OptTheta"] = new DoubleArray(optTheta); // write back optimized parameters so that we can use them in the Analysis method
225        return nmse;
226      }
227    }
228
229    private void OptimizeForEpisodes(
230      ISymbolicExpressionTree[] trees,
231      IRandom random,
232      IEnumerable<IntRange> episodes,
233      out double[] optTheta,
234      out double nmse) {
235      var rows = episodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start)).ToArray();
236      var problemData = ProblemData;
237      var targetVars = TargetVariables.CheckedItems.OrderBy(i => i.Index).Select(i => i.Value.Value).ToArray();
238      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
239      var targetValues = new double[rows.Length, targetVars.Length];
240
241      // collect values of all target variables
242      var colIdx = 0;
243      foreach (var targetVar in targetVars) {
244        int rowIdx = 0;
245        foreach (var value in problemData.Dataset.GetDoubleValues(targetVar, rows)) {
246          targetValues[rowIdx, colIdx] = value;
247          rowIdx++;
248        }
249        colIdx++;
250      }
251
252      // NOTE: the order of values in parameter matches prefix order of constant nodes in trees
253      var paramNodes = new List<ISymbolicExpressionTreeNode>();
254      foreach (var t in trees) {
255        foreach (var n in t.IterateNodesPrefix()) {
256          if (IsConstantNode(n))
257            paramNodes.Add(n);
258        }
259      }
260      // init params randomly from Unif(-1e-5, 1e-5)
261      // theta contains parameter values for trees and then the initial values for latent variables (a separate vector for each episode)
262      // inital values for latent variables are also optimized
263      var theta = new double[paramNodes.Count + latentVariables.Length * episodes.Count()];
264      for (int i = 0; i < theta.Length; i++)
265        theta[i] = random.NextDouble() * 2.0e-2 - 1.0e-2;
266
267      optTheta = new double[0];
268      if (theta.Length > 0) {
269        alglib.minlbfgsstate state;
270        alglib.minlbfgsreport report;
271        alglib.minlbfgscreate(Math.Min(theta.Length, 5), theta, out state);
272        alglib.minlbfgssetcond(state, 0.0, 0.0, 0.0, MaximumParameterOptimizationIterations);
273        //alglib.minlbfgssetgradientcheck(state, 1e-6);
274        alglib.minlbfgsoptimize(state, EvaluateObjectiveAndGradient, null,
275          new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver }); //TODO: create a type
276
277        alglib.minlbfgsresults(state, out optTheta, out report);
278
279        /*
280         *
281         *         L-BFGS algorithm results
282
283          INPUT PARAMETERS:
284              State   -   algorithm state
285
286          OUTPUT PARAMETERS:
287              X       -   array[0..N-1], solution
288              Rep     -   optimization report:
289                          * Rep.TerminationType completetion code:
290                              * -7    gradient verification failed.
291                                      See MinLBFGSSetGradientCheck() for more information.
292                              * -2    rounding errors prevent further improvement.
293                                      X contains best point found.
294                              * -1    incorrect parameters were specified
295                              *  1    relative function improvement is no more than
296                                      EpsF.
297                              *  2    relative step is no more than EpsX.
298                              *  4    gradient norm is no more than EpsG
299                              *  5    MaxIts steps was taken
300                              *  7    stopping conditions are too stringent,
301                                      further improvement is impossible
302                          * Rep.IterationsCount contains iterations count
303                          * NFEV countains number of function calculations
304         */
305        if (report.terminationtype < 0) { nmse = 10E6; return; }
306      }
307
308      // perform evaluation for optimal theta to get quality value
309      double[] grad = new double[optTheta.Length];
310      nmse = double.NaN;
311      EvaluateObjectiveAndGradient(optTheta, ref nmse, grad,
312        new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver });
313      if (double.IsNaN(nmse) || double.IsInfinity(nmse)) { nmse = 10E6; return; } // return a large value (TODO: be consistent by using NMSE)
314    }
315
316    private static void EvaluateObjectiveAndGradient(double[] x, ref double f, double[] grad, object obj) {
317      var trees = (ISymbolicExpressionTree[])((object[])obj)[0];
318      var targetVariables = (string[])((object[])obj)[1];
319      var problemData = (IRegressionProblemData)((object[])obj)[2];
320      var targetValues = (double[,])((object[])obj)[3];
321      var episodes = (IntRange[])((object[])obj)[4];
322      var numericIntegrationSteps = (int)((object[])obj)[5];
323      var latentVariables = (string[])((object[])obj)[6];
324      var odeSolver = (string)((object[])obj)[7];
325
326
327      Tuple<double, Vector>[][] predicted = null; // one array of predictions for each episode
328      // TODO: stop integration early for diverging solutions
329      predicted = Integrate(
330          trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
331          problemData.Dataset,
332          problemData.AllowedInputVariables.ToArray(),
333          targetVariables,
334          latentVariables,
335          episodes,
336          x,
337          odeSolver,
338          numericIntegrationSteps).ToArray();
339
340      if (predicted.Length != targetValues.GetLength(0)) {
341        f = double.MaxValue;
342        Array.Clear(grad, 0, grad.Length);
343        return;
344      }
345
346      // for normalized MSE = 1/variance(t) * MSE(t, pred)
347      // TODO: Perf. (by standardization of target variables before evaluation of all trees)     
348      var invVar = Enumerable.Range(0, targetVariables.Length)
349        .Select(c => Enumerable.Range(0, targetValues.GetLength(0)).Select(row => targetValues[row, c])) // column vectors
350        .Select(vec => vec.Variance())
351        .Select(v => 1.0 / v)
352        .ToArray();
353
354      // objective function is NMSE
355      f = 0.0;
356      int n = predicted.Length;
357      double invN = 1.0 / n;
358      var g = Vector.Zero;
359      int r = 0;
360      foreach (var y_pred in predicted) {
361        // y_pred contains the predicted values for target variables first and then predicted values for latent variables
362        for (int c = 0; c < targetVariables.Length; c++) {
363
364          var y_pred_f = y_pred[c].Item1;
365          var y = targetValues[r, c];
366
367          var res = (y - y_pred_f);
368          var ressq = res * res;
369          f += ressq * invN * invVar[c];
370          g += -2.0 * res * y_pred[c].Item2 * invN * invVar[c];
371        }
372        r++;
373      }
374
375      g.CopyTo(grad);
376    }
377
378    public override void Analyze(Individual[] individuals, double[] qualities, ResultCollection results, IRandom random) {
379      base.Analyze(individuals, qualities, results, random);
380
381      if (!results.ContainsKey("Prediction (training)")) {
382        results.Add(new Result("Prediction (training)", typeof(ReadOnlyItemList<DataTable>)));
383      }
384      if (!results.ContainsKey("Prediction (test)")) {
385        results.Add(new Result("Prediction (test)", typeof(ReadOnlyItemList<DataTable>)));
386      }
387      if (!results.ContainsKey("Models")) {
388        results.Add(new Result("Models", typeof(VariableCollection)));
389      }
390
391      var bestIndividualAndQuality = this.GetBestIndividual(individuals, qualities);
392      var trees = bestIndividualAndQuality.Item1.Values.Select(v => v.Value).OfType<ISymbolicExpressionTree>().ToArray(); // extract all trees from individual
393
394      var problemData = ProblemData;
395      var targetVars = TargetVariables.CheckedItems.OrderBy(i => i.Index).Select(i => i.Value.Value).ToArray();
396      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
397
398      var trainingList = new ItemList<DataTable>();
399
400      if (OptimizeParametersForEpisodes) {
401        var eIdx = 0;
402        var trainingPredictions = new List<Tuple<double, Vector>[][]>();
403        foreach (var episode in TrainingEpisodes) {
404          var episodes = new[] { episode };
405          var optTheta = ((DoubleArray)bestIndividualAndQuality.Item1["OptTheta_" + eIdx]).ToArray(); // see evaluate
406          var trainingPrediction = Integrate(
407                                   trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
408                                   problemData.Dataset,
409                                   problemData.AllowedInputVariables.ToArray(),
410                                   targetVars,
411                                   latentVariables,
412                                   episodes,
413                                   optTheta,
414                                   OdeSolver,
415                                   NumericIntegrationSteps).ToArray();
416          trainingPredictions.Add(trainingPrediction);
417          eIdx++;
418        }
419
420        // only for target values
421        var trainingRows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start));
422        for (int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
423          var targetVar = targetVars[colIdx];
424          var trainingDataTable = new DataTable(targetVar + " prediction (training)");
425          var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, trainingRows));
426          var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, trainingPredictions.SelectMany(arr => arr.Select(row => row[colIdx].Item1)).ToArray());
427          trainingDataTable.Rows.Add(actualValuesRow);
428          trainingDataTable.Rows.Add(predictedValuesRow);
429          trainingList.Add(trainingDataTable);
430        }
431        results["Prediction (training)"].Value = trainingList.AsReadOnly();
432
433
434        var models = new VariableCollection();
435
436        foreach (var tup in targetVars.Zip(trees, Tuple.Create)) {
437          var targetVarName = tup.Item1;
438          var tree = tup.Item2;
439
440          var origTreeVar = new HeuristicLab.Core.Variable(targetVarName + "(original)");
441          origTreeVar.Value = (ISymbolicExpressionTree)tree.Clone();
442          models.Add(origTreeVar);
443        }
444        results["Models"].Value = models;
445      } else {
446        var optTheta = ((DoubleArray)bestIndividualAndQuality.Item1["OptTheta"]).ToArray(); // see evaluate
447        var trainingPrediction = Integrate(
448                                   trees,  // we assume trees contain expressions for the change of each target variable over time dy/dt
449                                   problemData.Dataset,
450                                   problemData.AllowedInputVariables.ToArray(),
451                                   targetVars,
452                                   latentVariables,
453                                   TrainingEpisodes,
454                                   optTheta,
455                                   OdeSolver,
456                                   NumericIntegrationSteps).ToArray();
457        // for target values and latent variables
458        var trainingRows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start));
459        for (int colIdx = 0; colIdx < trees.Length; colIdx++) {
460          // is target variable
461          if (colIdx < targetVars.Length) {
462            var targetVar = targetVars[colIdx];
463            var trainingDataTable = new DataTable(targetVar + " prediction (training)");
464            var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, trainingRows));
465            var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, trainingPrediction.Select(arr => arr[colIdx].Item1).ToArray());
466            trainingDataTable.Rows.Add(actualValuesRow);
467            trainingDataTable.Rows.Add(predictedValuesRow);
468            trainingList.Add(trainingDataTable);
469          } else {
470            var latentVar = latentVariables[colIdx - targetVars.Length];
471            var trainingDataTable = new DataTable(latentVar + " prediction (training)");
472            var predictedValuesRow = new DataRow(latentVar + " pred.", "Predicted values for " + latentVar, trainingPrediction.Select(arr => arr[colIdx].Item1).ToArray());
473            var emptyRow = new DataRow(latentVar);
474            trainingDataTable.Rows.Add(emptyRow);
475            trainingDataTable.Rows.Add(predictedValuesRow);
476            trainingList.Add(trainingDataTable);
477          }
478        }
479        // TODO: DRY for training and test
480        var testList = new ItemList<DataTable>();
481        var testRows = ProblemData.TestIndices.ToArray();
482        var testPrediction = Integrate(
483         trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
484         problemData.Dataset,
485         problemData.AllowedInputVariables.ToArray(),
486         targetVars,
487         latentVariables,
488         new IntRange[] { ProblemData.TestPartition },
489         optTheta,
490         OdeSolver,
491         NumericIntegrationSteps).ToArray();
492
493        for (int colIdx = 0; colIdx < trees.Length; colIdx++) {
494          // is target variable
495          if (colIdx < targetVars.Length) {
496            var targetVar = targetVars[colIdx];
497            var testDataTable = new DataTable(targetVar + " prediction (test)");
498            var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, testRows));
499            var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, testPrediction.Select(arr => arr[colIdx].Item1).ToArray());
500            testDataTable.Rows.Add(actualValuesRow);
501            testDataTable.Rows.Add(predictedValuesRow);
502            testList.Add(testDataTable);
503
504          } else {
505            var latentVar = latentVariables[colIdx - targetVars.Length];
506            var testDataTable = new DataTable(latentVar + " prediction (test)");
507            var predictedValuesRow = new DataRow(latentVar + " pred.", "Predicted values for " + latentVar, testPrediction.Select(arr => arr[colIdx].Item1).ToArray());
508            var emptyRow = new DataRow(latentVar);
509            testDataTable.Rows.Add(emptyRow);
510            testDataTable.Rows.Add(predictedValuesRow);
511            testList.Add(testDataTable);
512          }
513        }
514
515        results["Prediction (training)"].Value = trainingList.AsReadOnly();
516        results["Prediction (test)"].Value = testList.AsReadOnly();
517        #region simplification of models
518        // TODO the dependency of HeuristicLab.Problems.DataAnalysis.Symbolic is not ideal
519        var models = new VariableCollection();    // to store target var names and original version of tree
520
521        int nextParIdx = 0;
522        for (int idx = 0; idx < trees.Length; idx++) {
523          var varName = string.Empty;
524          if (idx < targetVars.Length) {
525            varName = targetVars[idx];
526          } else {
527            varName = latentVariables[idx - targetVars.Length];
528          }
529          var tree = trees[idx];
530
531          // when we reference HeuristicLab.Problems.DataAnalysis.Symbolic we can translate symbols
532          var shownTree = new SymbolicExpressionTree(TranslateTreeNode(tree.Root, optTheta.ToArray(),
533            ref nextParIdx));
534
535          // var shownTree = (SymbolicExpressionTree)tree.Clone();
536          // var constantsNodeOrig = tree.IterateNodesPrefix().Where(IsConstantNode);
537          // var constantsNodeShown = shownTree.IterateNodesPrefix().Where(IsConstantNode);
538          //
539          // foreach (var n in constantsNodeOrig.Zip(constantsNodeShown, (original, shown) => new { original, shown })) {
540          //   double constantsVal = optTheta[nodeIdx[n.original]];
541          //
542          //   ConstantTreeNode replacementNode = new ConstantTreeNode(new Constant()) { Value = constantsVal };
543          //
544          //   var parentNode = n.shown.Parent;
545          //   int replacementIndex = parentNode.IndexOfSubtree(n.shown);
546          //   parentNode.RemoveSubtree(replacementIndex);
547          //   parentNode.InsertSubtree(replacementIndex, replacementNode);
548          // }
549
550          var origTreeVar = new HeuristicLab.Core.Variable(varName + "(original)");
551          origTreeVar.Value = (ISymbolicExpressionTree)tree.Clone();
552          models.Add(origTreeVar);
553          var simplifiedTreeVar = new HeuristicLab.Core.Variable(varName + "(simplified)");
554          simplifiedTreeVar.Value = TreeSimplifier.Simplify(shownTree);
555          models.Add(simplifiedTreeVar);
556
557        }
558        results["Models"].Value = models;
559        #endregion
560      }
561    }
562
563
564    #region interpretation
565
566    // the following uses auto-diff to calculate the gradient w.r.t. the parameters forward in time.
567    // this is basically the method described in Gronwall T. Note on the derivatives with respect to a parameter of the solutions of a system of differential equations. Ann. Math. 1919;20:292–296.
568
569    // a comparison of three potential calculation methods for the gradient is given in:
570    // Sengupta, B., Friston, K. J., & Penny, W. D. (2014). Efficient gradient computation for dynamical models. Neuroimage, 98(100), 521–527. http://doi.org/10.1016/j.neuroimage.2014.04.040
571    // "Our comparison establishes that the adjoint method is computationally more efficient for numerical estimation of parametric gradients
572    // for state-space models — both linear and non-linear, as in the case of a dynamical causal model (DCM)"
573
574    // for a solver with the necessary features see: https://computation.llnl.gov/projects/sundials/cvodes
575
576    private static IEnumerable<Tuple<double, Vector>[]> Integrate(
577      ISymbolicExpressionTree[] trees, IDataset dataset,
578      string[] inputVariables, string[] targetVariables, string[] latentVariables, IEnumerable<IntRange> episodes,
579      double[] parameterValues,
580      string odeSolver, int numericIntegrationSteps = 100) {
581
582      // TODO: numericIntegrationSteps is only relevant for the HeuristicLab solver
583      var episodeIdx = 0;
584      foreach (var episode in episodes) {
585        var rows = Enumerable.Range(episode.Start, episode.End - episode.Start);
586
587        // integrate forward starting with known values for the target in t0
588
589        var variableValues = new Dictionary<string, Tuple<double, Vector>>();
590        var t0 = rows.First();
591        foreach (var varName in inputVariables) {
592          variableValues.Add(varName, Tuple.Create(dataset.GetDoubleValue(varName, t0), Vector.Zero));
593        }
594        foreach (var varName in targetVariables) {
595          variableValues.Add(varName, Tuple.Create(dataset.GetDoubleValue(varName, t0), Vector.Zero));
596        }
597        // add value entries for latent variables which are also integrated
598        // initial values are at the end of the parameter vector
599        // separete initial values for each episode
600        var initialValueIdx = parameterValues.Length - episodes.Count() * latentVariables.Length + episodeIdx * latentVariables.Length;
601        foreach (var latentVar in latentVariables) {
602          var arr = new double[parameterValues.Length]; // backing array
603          arr[initialValueIdx] = 1.0;
604          var g = new Vector(arr);
605          variableValues.Add(latentVar,
606            Tuple.Create(parameterValues[initialValueIdx], g)); // we don't have observations for latent variables therefore we optimize the initial value for each episode
607          initialValueIdx++;
608        }
609
610        var calculatedVariables = targetVariables.Concat(latentVariables).ToArray(); // TODO: must conincide with the order of trees in the encoding
611
612        // return first value as stored in the dataset
613        yield return calculatedVariables
614          .Select(calcVarName => variableValues[calcVarName])
615          .ToArray();
616
617        var prevT = rows.First(); // TODO: here we should use a variable for t if it is available. Right now we assume equidistant measurements.
618        foreach (var t in rows.Skip(1)) {
619          if (odeSolver == "HeuristicLab")
620            IntegrateHL(trees, calculatedVariables, variableValues, parameterValues, numericIntegrationSteps);
621          else if (odeSolver == "CVODES")
622            IntegrateCVODES(trees, calculatedVariables, variableValues, parameterValues, t - prevT);
623          else throw new InvalidOperationException("Unknown ODE solver " + odeSolver);
624          prevT = t;
625
626          // This check doesn't work with the HeuristicLab integrator if there are input variables
627          //if (variableValues.Count == targetVariables.Length) {
628          // only return the target variables for calculation of errors
629          var res = calculatedVariables
630            .Select(targetVar => variableValues[targetVar])
631            .ToArray();
632          if (res.Any(ri => double.IsNaN(ri.Item1) || double.IsInfinity(ri.Item1))) yield break;
633          yield return res;
634          //} else {
635          //  yield break; // stop early on problems
636          //}
637
638
639          // update for next time step
640          foreach (var varName in inputVariables) {
641            variableValues[varName] = Tuple.Create(dataset.GetDoubleValue(varName, t), Vector.Zero);
642          }
643        }
644        episodeIdx++;
645      }
646    }
647
648
649    /// <summary>
650    ///  Here we use CVODES to solve the ODE. Forward sensitivities are used to calculate the gradient for parameter optimization
651    /// </summary>
652    /// <param name="trees">Each equation in the ODE represented as a tree</param>
653    /// <param name="calculatedVariables">The names of the calculated variables</param>
654    /// <param name="variableValues">The start values of the calculated variables as well as their sensitivites over parameters</param>
655    /// <param name="parameterValues">The current parameter values</param>
656    /// <param name="t">The time t up to which we need to integrate.</param>
657    private static void IntegrateCVODES(
658      ISymbolicExpressionTree[] trees, // f(y,p) in tree representation
659      string[] calculatedVariables, // names of elements of y
660      Dictionary<string, Tuple<double, Vector>> variableValues,  //  y (intput and output) input: y(t0), output: y(t0+t)
661      double[] parameterValues, // p
662      double t // duration t for which we want to integrate
663      ) {
664
665      // the RHS of the ODE
666      // dy/dt = f(y_t,x_t,p)
667      CVODES.CVRhsFunc f = CreateOdeRhs(trees, calculatedVariables, parameterValues);
668      // the Jacobian ∂f/∂y
669      CVODES.CVDlsJacFunc jac = CreateJac(trees, calculatedVariables, parameterValues);
670
671      // the RHS for the forward sensitivities (∂f/∂y)s_i(t) + ∂f/∂p_i
672      CVODES.CVSensRhsFn sensF = CreateSensitivityRhs(trees, calculatedVariables, parameterValues);
673
674      // setup solver
675      int numberOfEquations = trees.Length;
676      IntPtr y = IntPtr.Zero;
677      IntPtr cvode_mem = IntPtr.Zero;
678      IntPtr A = IntPtr.Zero;
679      IntPtr yS0 = IntPtr.Zero;
680      IntPtr linearSolver = IntPtr.Zero;
681      var ns = parameterValues.Length; // number of parameters
682
683      try {
684        y = CVODES.N_VNew_Serial(numberOfEquations);
685        // init y to current values of variables
686        // y must be initialized before calling CVodeInit
687        for (int i = 0; i < calculatedVariables.Length; i++) {
688          CVODES.NV_Set_Ith_S(y, i, variableValues[calculatedVariables[i]].Item1);
689        }
690
691        cvode_mem = CVODES.CVodeCreate(CVODES.MultistepMethod.CV_ADAMS, CVODES.NonlinearSolverIteration.CV_FUNCTIONAL);
692
693        var flag = CVODES.CVodeInit(cvode_mem, f, 0.0, y);
694        Debug.Assert(CVODES.CV_SUCCESS == flag);
695
696        double relTol = 1.0e-2;
697        double absTol = 1.0;
698        flag = CVODES.CVodeSStolerances(cvode_mem, relTol, absTol);  // TODO: probably need to adjust absTol per variable
699        Debug.Assert(CVODES.CV_SUCCESS == flag);
700
701        A = CVODES.SUNDenseMatrix(numberOfEquations, numberOfEquations);
702        Debug.Assert(A != IntPtr.Zero);
703
704        linearSolver = CVODES.SUNDenseLinearSolver(y, A);
705        Debug.Assert(linearSolver != IntPtr.Zero);
706
707        flag = CVODES.CVDlsSetLinearSolver(cvode_mem, linearSolver, A);
708        Debug.Assert(CVODES.CV_SUCCESS == flag);
709
710        flag = CVODES.CVDlsSetJacFn(cvode_mem, jac);
711        Debug.Assert(CVODES.CV_SUCCESS == flag);
712
713        yS0 = CVODES.N_VCloneVectorArray_Serial(ns, y); // clone the output vector for each parameter
714        unsafe {
715          // set to initial sensitivities supplied by caller
716          for (int pIdx = 0; pIdx < ns; pIdx++) {
717            var yS0_i = *((IntPtr*)yS0.ToPointer() + pIdx);
718            for (var varIdx = 0; varIdx < calculatedVariables.Length; varIdx++) {
719              CVODES.NV_Set_Ith_S(yS0_i, varIdx, variableValues[calculatedVariables[varIdx]].Item2[pIdx]);
720            }
721          }
722        }
723
724        flag = CVODES.CVodeSensInit(cvode_mem, ns, CVODES.CV_SIMULTANEOUS, sensF, yS0);
725        Debug.Assert(CVODES.CV_SUCCESS == flag);
726
727        flag = CVODES.CVodeSensEEtolerances(cvode_mem);
728        Debug.Assert(CVODES.CV_SUCCESS == flag);
729
730        // make one forward integration step
731        double tout = 0.0; // first output time
732        flag = CVODES.CVode(cvode_mem, t, y, ref tout, CVODES.CV_NORMAL);
733        if (flag == CVODES.CV_SUCCESS) {
734          Debug.Assert(t == tout);
735
736          // get sensitivities
737          flag = CVODES.CVodeGetSens(cvode_mem, ref tout, yS0);
738          Debug.Assert(CVODES.CV_SUCCESS == flag);
739
740          // update variableValues based on integration results
741          for (int varIdx = 0; varIdx < calculatedVariables.Length; varIdx++) {
742            var yi = CVODES.NV_Get_Ith_S(y, varIdx);
743            var gArr = new double[parameterValues.Length];
744            for (var pIdx = 0; pIdx < parameterValues.Length; pIdx++) {
745              unsafe {
746                var yS0_pi = *((IntPtr*)yS0.ToPointer() + pIdx);
747                gArr[pIdx] = CVODES.NV_Get_Ith_S(yS0_pi, varIdx);
748              }
749            }
750            variableValues[calculatedVariables[varIdx]] = Tuple.Create(yi, new Vector(gArr));
751          }
752        } else {
753          variableValues.Clear();   // indicate problems by not returning new values
754        }
755
756        // cleanup all allocated objects
757      } finally {
758        if (y != IntPtr.Zero) CVODES.N_VDestroy_Serial(y);
759        if (cvode_mem != IntPtr.Zero) CVODES.CVodeFree(ref cvode_mem);
760        if (linearSolver != IntPtr.Zero) CVODES.SUNLinSolFree(linearSolver);
761        if (A != IntPtr.Zero) CVODES.SUNMatDestroy(A);
762        if (yS0 != IntPtr.Zero) CVODES.N_VDestroyVectorArray_Serial(yS0, ns);
763      }
764    }
765
766
767    private static CVODES.CVRhsFunc CreateOdeRhs(
768      ISymbolicExpressionTree[] trees,
769      string[] calculatedVariables,
770      double[] parameterValues) {
771      // we don't need to calculate a gradient here -> no nodes are selected for
772      // --> no nodes are selected to be included in the gradient calculation
773      var nodeIdx = new Dictionary<ISymbolicExpressionTreeNode, int>();
774      return (double t,
775              IntPtr y, // N_Vector, current value of y (input)
776              IntPtr ydot, // N_Vector (calculated value of y' (output)
777              IntPtr user_data // optional user data, (unused here)
778              ) => {
779                // TODO: perf
780                var nodeValues = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
781
782                int pIdx = 0;
783                foreach (var tree in trees) {
784                  foreach (var n in tree.IterateNodesPrefix()) {
785                    if (IsConstantNode(n)) {
786                      nodeValues.Add(n, Tuple.Create(parameterValues[pIdx], Vector.Zero)); // here we do not need a gradient
787                      pIdx++;
788                    } else if (n.SubtreeCount == 0) {
789                      // for variables and latent variables get the value from variableValues
790                      var varName = n.Symbol.Name;
791                      var varIdx = Array.IndexOf(calculatedVariables, varName); // TODO: perf!
792                      if (varIdx < 0) throw new InvalidProgramException();
793                      var y_i = CVODES.NV_Get_Ith_S(y, (long)varIdx);
794                      nodeValues.Add(n, Tuple.Create(y_i, Vector.Zero)); // no gradient needed
795                    }
796                  }
797                }
798                for (int i = 0; i < trees.Length; i++) {
799                  var tree = trees[i];
800                  var res_i = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues);
801                  CVODES.NV_Set_Ith_S(ydot, i, res_i.Item1);
802                }
803                return 0;
804              };
805    }
806
807    private static CVODES.CVDlsJacFunc CreateJac(
808      ISymbolicExpressionTree[] trees,
809      string[] calculatedVariables,
810      double[] parameterValues) {
811
812      return (
813        double t, // current time (input)
814        IntPtr y, // N_Vector, current value of y (input)
815        IntPtr fy, // N_Vector, current value of f (input)
816        IntPtr Jac, // SUNMatrix ∂f/∂y (output, rows i contains are ∂f_i/∂y vector)
817        IntPtr user_data, // optional (unused here)
818        IntPtr tmp1, // N_Vector, optional (unused here)
819        IntPtr tmp2, // N_Vector, optional (unused here)
820        IntPtr tmp3 // N_Vector, optional (unused here)
821      ) => {
822        // here we need to calculate partial derivatives for the calculated variables y
823        var nodeValues = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
824        int pIdx = 0;
825        foreach (var tree in trees) {
826          foreach (var n in tree.IterateNodesPrefix()) {
827            if (IsConstantNode(n)) {
828              nodeValues.Add(n, Tuple.Create(parameterValues[pIdx], Vector.Zero)); // here we need a gradient over y which is zero for parameters
829              pIdx++;
830            } else if (n.SubtreeCount == 0) {
831              // for variables and latent variables we use values supplied in y and init gradient vectors accordingly
832              var varName = n.Symbol.Name;
833              var varIdx = Array.IndexOf(calculatedVariables, varName); // TODO: perf!
834              if (varIdx < 0) throw new InvalidProgramException();
835
836              var y_i = CVODES.NV_Get_Ith_S(y, (long)varIdx);
837              var gArr = new double[CVODES.NV_LENGTH_S(y)]; // backing array
838              gArr[varIdx] = 1.0;
839              var g = new Vector(gArr);
840              nodeValues.Add(n, Tuple.Create(y_i, g));
841            }
842          }
843        }
844
845        for (int i = 0; i < trees.Length; i++) {
846          var tree = trees[i];
847          var res = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues);
848          var g = res.Item2;
849          for (int j = 0; j < calculatedVariables.Length; j++) {
850            CVODES.SUNDenseMatrix_Set(Jac, i, j, g[j]);
851          }
852        }
853        return 0; // on success
854      };
855    }
856
857
858    // to calculate sensitivities RHS for all equations at once
859    // must compute (∂f/∂y)s_i(t) + ∂f/∂p_i and store in ySdot.
860    // Index i refers to parameters, dimensionality of matrix and vectors is number of equations
861    private static CVODES.CVSensRhsFn CreateSensitivityRhs(ISymbolicExpressionTree[] trees, string[] calculatedVariables, double[] parameterValues) {
862      return (
863              int Ns, // number of parameters
864              double t, // current time
865              IntPtr y, // N_Vector y(t) (input)
866              IntPtr ydot, // N_Vector dy/dt(t) (input)
867              IntPtr yS, // N_Vector*, one vector for each parameter (input)
868              IntPtr ySdot, // N_Vector*, one vector for each parameter (output)
869              IntPtr user_data, // optional (unused here)
870              IntPtr tmp1, // N_Vector, optional (unused here)
871              IntPtr tmp2 // N_Vector, optional (unused here)
872        ) => {
873          // here we need to calculate partial derivatives for the calculated variables y as well as for the parameters
874          var nodeValues = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
875          var d = calculatedVariables.Length + parameterValues.Length; // dimensionality of gradient
876          // first collect variable values
877          foreach (var tree in trees) {
878            foreach (var n in tree.IterateNodesPrefix()) {
879              if (IsVariableNode(n)) {
880                // for variables and latent variables we use values supplied in y and init gradient vectors accordingly
881                var varName = n.Symbol.Name;
882                var varIdx = Array.IndexOf(calculatedVariables, varName); // TODO: perf!
883                if (varIdx < 0) throw new InvalidProgramException();
884
885                var y_i = CVODES.NV_Get_Ith_S(y, (long)varIdx);
886                var gArr = new double[d]; // backing array
887                gArr[varIdx] = 1.0;
888                var g = new Vector(gArr);
889                nodeValues.Add(n, Tuple.Create(y_i, g));
890              }
891            }
892          }
893          // then collect constants
894          int pIdx = 0;
895          foreach (var tree in trees) {
896            foreach (var n in tree.IterateNodesPrefix()) {
897              if (IsConstantNode(n)) {
898                var gArr = new double[d];
899                gArr[calculatedVariables.Length + pIdx] = 1.0;
900                var g = new Vector(gArr);
901                nodeValues.Add(n, Tuple.Create(parameterValues[pIdx], g));
902                pIdx++;
903              }
904            }
905          }
906          // gradient vector is [∂f/∂y_1, ∂f/∂y_2, ... ∂f/∂yN, ∂f/∂p_1 ... ∂f/∂p_K]
907
908
909          for (pIdx = 0; pIdx < Ns; pIdx++) {
910            unsafe {
911              var sDot_pi = *((IntPtr*)ySdot.ToPointer() + pIdx);
912              CVODES.N_VConst_Serial(0.0, sDot_pi);
913            }
914          }
915
916          for (int i = 0; i < trees.Length; i++) {
917            var tree = trees[i];
918            var res = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues);
919            var g = res.Item2;
920
921
922            // update ySdot = (∂f/∂y)s_i(t) + ∂f/∂p_i
923
924            for (pIdx = 0; pIdx < Ns; pIdx++) {
925              unsafe {
926                var sDot_pi = *((IntPtr*)ySdot.ToPointer() + pIdx);
927                var s_pi = *((IntPtr*)yS.ToPointer() + pIdx);
928
929                var v = CVODES.NV_Get_Ith_S(sDot_pi, i);
930                // (∂f/∂y)s_i(t)
931                var p = 0.0;
932                for (int yIdx = 0; yIdx < calculatedVariables.Length; yIdx++) {
933                  p += g[yIdx] * CVODES.NV_Get_Ith_S(s_pi, yIdx);
934                }
935                // + ∂f/∂p_i
936                CVODES.NV_Set_Ith_S(sDot_pi, i, v + p + g[calculatedVariables.Length + pIdx]);
937              }
938            }
939
940          }
941          return 0; // on success
942        };
943    }
944
945    private static void IntegrateHL(
946      ISymbolicExpressionTree[] trees,
947      string[] calculatedVariables, // names of integrated variables
948      Dictionary<string, Tuple<double, Vector>> variableValues, //y (intput and output) input: y(t0), output: y(t0+1)
949      double[] parameterValues,
950      int numericIntegrationSteps) {
951
952      var nodeValues = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
953
954      // the nodeValues for parameters are constant over time
955      // TODO: this needs to be done only once for each iteration in gradient descent (whenever parameter values change)
956      // NOTE: the order must be the same as above (prefix order for nodes)
957      int pIdx = 0;
958      foreach (var tree in trees) {
959        foreach (var node in tree.Root.IterateNodesPrefix()) {
960          if (IsConstantNode(node)) {
961            var gArr = new double[parameterValues.Length]; // backing array
962            gArr[pIdx] = 1.0;
963            var g = new Vector(gArr);
964            nodeValues.Add(node, new Tuple<double, Vector>(parameterValues[pIdx], g));
965            pIdx++;
966          } else if (node.SubtreeCount == 0) {
967            // for (latent) variables get the values from variableValues
968            var varName = node.Symbol.Name;
969            nodeValues.Add(node, variableValues[varName]);
970          }
971        }
972      }
973
974
975      double h = 1.0 / numericIntegrationSteps;
976      for (int step = 0; step < numericIntegrationSteps; step++) {
977        var deltaValues = new Dictionary<string, Tuple<double, Vector>>();
978        for (int i = 0; i < trees.Length; i++) {
979          var tree = trees[i];
980          var targetVarName = calculatedVariables[i];
981
982          // Root.GetSubtree(0).GetSubtree(0) skips programRoot and startSymbol
983          var res = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues);
984          deltaValues.Add(targetVarName, res);
985        }
986
987        // update variableValues for next step, trapezoid integration
988        foreach (var kvp in deltaValues) {
989          var oldVal = variableValues[kvp.Key];
990          var newVal = Tuple.Create(
991            oldVal.Item1 + h * kvp.Value.Item1,
992            oldVal.Item2 + h * kvp.Value.Item2
993          );
994          variableValues[kvp.Key] = newVal;
995        }
996        // update nodeValues from variableValues
997        // TODO: perf using dictionary with list of nodes for each variable
998        foreach (var tree in trees) {
999          foreach (var node in tree.Root.IterateNodesPrefix().Where(n => IsVariableNode(n))) {
1000            var varName = node.Symbol.Name;
1001            nodeValues[node] = variableValues[varName];
1002          }
1003        }
1004      }
1005    }
1006
1007    private static Tuple<double, Vector> InterpretRec(
1008      ISymbolicExpressionTreeNode node,
1009      Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>> nodeValues      // contains value and gradient vector for a node (variables and constants only)
1010        ) {
1011
1012      switch (node.Symbol.Name) {
1013        case "+": {
1014            var l = InterpretRec(node.GetSubtree(0), nodeValues);
1015            var r = InterpretRec(node.GetSubtree(1), nodeValues);
1016
1017            return Tuple.Create(l.Item1 + r.Item1, l.Item2 + r.Item2);
1018          }
1019        case "*": {
1020            var l = InterpretRec(node.GetSubtree(0), nodeValues);
1021            var r = InterpretRec(node.GetSubtree(1), nodeValues);
1022
1023            return Tuple.Create(l.Item1 * r.Item1, l.Item2 * r.Item1 + l.Item1 * r.Item2);
1024          }
1025
1026        case "-": {
1027            var l = InterpretRec(node.GetSubtree(0), nodeValues);
1028            var r = InterpretRec(node.GetSubtree(1), nodeValues);
1029
1030            return Tuple.Create(l.Item1 - r.Item1, l.Item2 - r.Item2);
1031          }
1032        case "%": {
1033            var l = InterpretRec(node.GetSubtree(0), nodeValues);
1034            var r = InterpretRec(node.GetSubtree(1), nodeValues);
1035
1036            // protected division
1037            if (r.Item1.IsAlmost(0.0)) {
1038              return Tuple.Create(0.0, Vector.Zero);
1039            } else {
1040              return Tuple.Create(
1041                l.Item1 / r.Item1,
1042                l.Item1 * -1.0 / (r.Item1 * r.Item1) * r.Item2 + 1.0 / r.Item1 * l.Item2 // (f/g)' = f * (1/g)' + 1/g * f' = f * -1/g² * g' + 1/g * f'
1043                );
1044            }
1045          }
1046        case "sin": {
1047            var x = InterpretRec(node.GetSubtree(0), nodeValues);
1048            return Tuple.Create(
1049              Math.Sin(x.Item1),
1050              Vector.Cos(x.Item2) * x.Item2
1051            );
1052          }
1053        case "cos": {
1054            var x = InterpretRec(node.GetSubtree(0), nodeValues);
1055            return Tuple.Create(
1056              Math.Cos(x.Item1),
1057              -Vector.Sin(x.Item2) * x.Item2
1058            );
1059          }
1060        case "sqr": {
1061            var x = InterpretRec(node.GetSubtree(0), nodeValues);
1062            return Tuple.Create(
1063              x.Item1 * x.Item1,
1064              2.0 * x.Item1 * x.Item2
1065            );
1066          }
1067        default: {
1068            return nodeValues[node];  // value and gradient for constants and variables must be set by the caller
1069          }
1070      }
1071    }
1072    #endregion
1073
1074    #region events
1075    /*
1076     * Dependencies between parameters:
1077     *
1078     * ProblemData
1079     *    |
1080     *    V
1081     * TargetVariables   FunctionSet    MaximumLength    NumberOfLatentVariables
1082     *               |   |                 |                   |
1083     *               V   V                 |                   |
1084     *             Grammar <---------------+-------------------
1085     *                |
1086     *                V
1087     *            Encoding
1088     */
1089    private void RegisterEventHandlers() {
1090      ProblemDataParameter.ValueChanged += ProblemDataParameter_ValueChanged;
1091      if (ProblemDataParameter.Value != null) ProblemDataParameter.Value.Changed += ProblemData_Changed;
1092
1093      TargetVariablesParameter.ValueChanged += TargetVariablesParameter_ValueChanged;
1094      if (TargetVariablesParameter.Value != null) TargetVariablesParameter.Value.CheckedItemsChanged += CheckedTargetVariablesChanged;
1095
1096      FunctionSetParameter.ValueChanged += FunctionSetParameter_ValueChanged;
1097      if (FunctionSetParameter.Value != null) FunctionSetParameter.Value.CheckedItemsChanged += CheckedFunctionsChanged;
1098
1099      MaximumLengthParameter.Value.ValueChanged += MaximumLengthChanged;
1100
1101      NumberOfLatentVariablesParameter.Value.ValueChanged += NumLatentVariablesChanged;
1102    }
1103
1104    private void NumLatentVariablesChanged(object sender, EventArgs e) {
1105      UpdateGrammarAndEncoding();
1106    }
1107
1108    private void MaximumLengthChanged(object sender, EventArgs e) {
1109      UpdateGrammarAndEncoding();
1110    }
1111
1112    private void FunctionSetParameter_ValueChanged(object sender, EventArgs e) {
1113      FunctionSetParameter.Value.CheckedItemsChanged += CheckedFunctionsChanged;
1114    }
1115
1116    private void CheckedFunctionsChanged(object sender, CollectionItemsChangedEventArgs<IndexedItem<StringValue>> e) {
1117      UpdateGrammarAndEncoding();
1118    }
1119
1120    private void TargetVariablesParameter_ValueChanged(object sender, EventArgs e) {
1121      TargetVariablesParameter.Value.CheckedItemsChanged += CheckedTargetVariablesChanged;
1122    }
1123
1124    private void CheckedTargetVariablesChanged(object sender, CollectionItemsChangedEventArgs<IndexedItem<StringValue>> e) {
1125      UpdateGrammarAndEncoding();
1126    }
1127
1128    private void ProblemDataParameter_ValueChanged(object sender, EventArgs e) {
1129      ProblemDataParameter.Value.Changed += ProblemData_Changed;
1130      OnProblemDataChanged();
1131      OnReset();
1132    }
1133
1134    private void ProblemData_Changed(object sender, EventArgs e) {
1135      OnProblemDataChanged();
1136      OnReset();
1137    }
1138
1139    private void OnProblemDataChanged() {
1140      UpdateTargetVariables();        // implicitly updates other dependent parameters
1141      var handler = ProblemDataChanged;
1142      if (handler != null) handler(this, EventArgs.Empty);
1143    }
1144
1145    #endregion
1146
1147    #region  helper
1148
1149    private void InitAllParameters() {
1150      UpdateTargetVariables(); // implicitly updates the grammar and the encoding     
1151    }
1152
1153    private ReadOnlyCheckedItemList<StringValue> CreateFunctionSet() {
1154      var l = new CheckedItemList<StringValue>();
1155      l.Add(new StringValue("+").AsReadOnly());
1156      l.Add(new StringValue("*").AsReadOnly());
1157      l.Add(new StringValue("%").AsReadOnly());
1158      l.Add(new StringValue("-").AsReadOnly());
1159      l.Add(new StringValue("sin").AsReadOnly());
1160      l.Add(new StringValue("cos").AsReadOnly());
1161      l.Add(new StringValue("sqr").AsReadOnly());
1162      return l.AsReadOnly();
1163    }
1164
1165    private static bool IsConstantNode(ISymbolicExpressionTreeNode n) {
1166      return n.Symbol.Name.StartsWith("θ");
1167    }
1168    private static bool IsLatentVariableNode(ISymbolicExpressionTreeNode n) {
1169      return n.Symbol.Name.StartsWith("λ");
1170    }
1171    private static bool IsVariableNode(ISymbolicExpressionTreeNode n) {
1172      return (n.SubtreeCount == 0) && !IsConstantNode(n) && !IsLatentVariableNode(n);
1173    }
1174
1175
1176    private void UpdateTargetVariables() {
1177      var currentlySelectedVariables = TargetVariables.CheckedItems
1178        .OrderBy(i => i.Index)
1179        .Select(i => i.Value.Value)
1180        .ToArray();
1181
1182      var newVariablesList = new CheckedItemList<StringValue>(ProblemData.Dataset.VariableNames.Select(str => new StringValue(str).AsReadOnly()).ToArray()).AsReadOnly();
1183      var matchingItems = newVariablesList.Where(item => currentlySelectedVariables.Contains(item.Value)).ToArray();
1184      foreach (var matchingItem in matchingItems) {
1185        newVariablesList.SetItemCheckedState(matchingItem, true);
1186      }
1187      TargetVariablesParameter.Value = newVariablesList;
1188    }
1189
1190    private void UpdateGrammarAndEncoding() {
1191      var encoding = new MultiEncoding();
1192      var g = CreateGrammar();
1193      foreach (var targetVar in TargetVariables.CheckedItems) {
1194        encoding = encoding.Add(new SymbolicExpressionTreeEncoding(targetVar + "_tree", g, MaximumLength, MaximumLength)); // only limit by length
1195      }
1196      for (int i = 1; i <= NumberOfLatentVariables; i++) {
1197        encoding = encoding.Add(new SymbolicExpressionTreeEncoding("λ" + i + "_tree", g, MaximumLength, MaximumLength));
1198      }
1199      Encoding = encoding;
1200    }
1201
1202    private ISymbolicExpressionGrammar CreateGrammar() {
1203      // whenever ProblemData is changed we create a new grammar with the necessary symbols
1204      var g = new SimpleSymbolicExpressionGrammar();
1205      g.AddSymbols(FunctionSet.CheckedItems.OrderBy(i => i.Index).Select(i => i.Value.Value).ToArray(), 2, 2);
1206
1207      // TODO
1208      //g.AddSymbols(new[] {
1209      //  "exp",
1210      //  "log", // log( <expr> ) // TODO: init a theta to ensure the value is always positive
1211      //  "exp_minus" // exp((-1) * <expr>
1212      //}, 1, 1);
1213
1214      foreach (var variableName in ProblemData.AllowedInputVariables.Union(TargetVariables.CheckedItems.Select(i => i.Value.Value)))
1215        g.AddTerminalSymbol(variableName);
1216
1217      // generate symbols for numeric parameters for which the value is optimized using AutoDiff
1218      // we generate multiple symbols to balance the probability for selecting a numeric parameter in the generation of random trees
1219      var numericConstantsFactor = 2.0;
1220      for (int i = 0; i < numericConstantsFactor * (ProblemData.AllowedInputVariables.Count() + TargetVariables.CheckedItems.Count()); i++) {
1221        g.AddTerminalSymbol("θ" + i); // numeric parameter for which the value is optimized using AutoDiff
1222      }
1223
1224      // generate symbols for latent variables
1225      for (int i = 1; i <= NumberOfLatentVariables; i++) {
1226        g.AddTerminalSymbol("λ" + i); // numeric parameter for which the value is optimized using AutoDiff
1227      }
1228
1229      return g;
1230    }
1231
1232
1233    private ISymbolicExpressionTreeNode TranslateTreeNode(ISymbolicExpressionTreeNode n, double[] parameterValues, ref int nextParIdx) {
1234      ISymbolicExpressionTreeNode translatedNode = null;
1235      if (n.Symbol is StartSymbol) {
1236        translatedNode = new StartSymbol().CreateTreeNode();
1237      } else if (n.Symbol is ProgramRootSymbol) {
1238        translatedNode = new ProgramRootSymbol().CreateTreeNode();
1239      } else if (n.Symbol.Name == "+") {
1240        translatedNode = new Addition().CreateTreeNode();
1241      } else if (n.Symbol.Name == "-") {
1242        translatedNode = new Subtraction().CreateTreeNode();
1243      } else if (n.Symbol.Name == "*") {
1244        translatedNode = new Multiplication().CreateTreeNode();
1245      } else if (n.Symbol.Name == "%") {
1246        translatedNode = new Division().CreateTreeNode();
1247      } else if (n.Symbol.Name == "sin") {
1248        translatedNode = new Sine().CreateTreeNode();
1249      } else if (n.Symbol.Name == "cos") {
1250        translatedNode = new Cosine().CreateTreeNode();
1251      } else if (n.Symbol.Name == "sqr") {
1252        translatedNode = new Square().CreateTreeNode();
1253      } else if (IsConstantNode(n)) {
1254        var constNode = (ConstantTreeNode)new Constant().CreateTreeNode();
1255        constNode.Value = parameterValues[nextParIdx];
1256        nextParIdx++;
1257        translatedNode = constNode;
1258      } else {
1259        // assume a variable name
1260        var varName = n.Symbol.Name;
1261        var varNode = (VariableTreeNode)new Variable().CreateTreeNode();
1262        varNode.Weight = 1.0;
1263        varNode.VariableName = varName;
1264        translatedNode = varNode;
1265      }
1266      foreach (var child in n.Subtrees) {
1267        translatedNode.AddSubtree(TranslateTreeNode(child, parameterValues, ref nextParIdx));
1268      }
1269      return translatedNode;
1270    }
1271    #endregion
1272
1273    #region Import & Export
1274    public void Load(IRegressionProblemData data) {
1275      Name = data.Name;
1276      Description = data.Description;
1277      ProblemData = data;
1278    }
1279
1280    public IRegressionProblemData Export() {
1281      return ProblemData;
1282    }
1283    #endregion
1284
1285  }
1286}
Note: See TracBrowser for help on using the repository browser.