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

Last change on this file since 16398 was 16398, checked in by gkronber, 8 months ago

#2925: small changes

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