Free cookie consent management tool by TermsFeed Policy Generator

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

Last change on this file since 16597 was 16597, checked in by gkronber, 5 years ago

#2925: made some adaptations while debugging parameter identification for dynamical models

File size: 67.5 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      var problemData = ProblemData;
204      var targetVars = TargetVariables.CheckedItems.OrderBy(i => i.Index).Select(i => i.Value.Value).ToArray();
205      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
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, problemData, targetVars, latentVariables, random, new[] { episode }, MaximumParameterOptimizationIterations, NumericIntegrationSteps, OdeSolver, 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, problemData, targetVars, latentVariables, random, TrainingEpisodes, MaximumParameterOptimizationIterations, NumericIntegrationSteps, OdeSolver, 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    public static void OptimizeForEpisodes(
230      ISymbolicExpressionTree[] trees,
231      IRegressionProblemData problemData,
232      string[] targetVars,
233      string[] latentVariables,
234      IRandom random,
235      IEnumerable<IntRange> episodes,
236      int maxParameterOptIterations,
237      int numericIntegrationSteps,
238      string odeSolver,
239      out double[] optTheta,
240      out double nmse) {
241      var rows = episodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start)).ToArray();
242      var targetValues = new double[rows.Length, targetVars.Length];
243
244      // collect values of all target variables
245      var colIdx = 0;
246      foreach (var targetVar in targetVars) {
247        int rowIdx = 0;
248        foreach (var value in problemData.Dataset.GetDoubleValues(targetVar, rows)) {
249          targetValues[rowIdx, colIdx] = value;
250          rowIdx++;
251        }
252        colIdx++;
253      }
254
255      // NOTE: the order of values in parameter matches prefix order of constant nodes in trees
256      var paramNodes = new List<ISymbolicExpressionTreeNode>();
257      foreach (var t in trees) {
258        foreach (var n in t.IterateNodesPrefix()) {
259          if (IsConstantNode(n))
260            paramNodes.Add(n);
261        }
262      }
263      // init params randomly
264      // theta contains parameter values for trees and then the initial values for latent variables (a separate vector for each episode)
265      // inital values for latent variables are also optimized
266      var theta = new double[paramNodes.Count + latentVariables.Length * episodes.Count()];
267      for (int i = 0; i < theta.Length; i++)
268        theta[i] = random.NextDouble() * 2.0 - 1.0;
269
270      optTheta = new double[0];
271      if (theta.Length > 0) {
272        alglib.minlbfgsstate state;
273        alglib.minlbfgsreport report;
274        alglib.minlbfgscreate(Math.Min(theta.Length, 5), theta, out state);
275        alglib.minlbfgssetcond(state, 0.0, 0.0, 0.0, maxParameterOptIterations);
276        // alglib.minlbfgssetgradientcheck(state, 1e-4);
277        alglib.minlbfgsoptimize(state, EvaluateObjectiveAndGradient, null,
278          new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver }); //TODO: create a type
279
280        alglib.minlbfgsresults(state, out optTheta, out report);
281
282        /*
283         *
284         *         L-BFGS algorithm results
285
286          INPUT PARAMETERS:
287              State   -   algorithm state
288
289          OUTPUT PARAMETERS:
290              X       -   array[0..N-1], solution
291              Rep     -   optimization report:
292                          * Rep.TerminationType completetion code:
293                              * -7    gradient verification failed.
294                                      See MinLBFGSSetGradientCheck() for more information.
295                              * -2    rounding errors prevent further improvement.
296                                      X contains best point found.
297                              * -1    incorrect parameters were specified
298                              *  1    relative function improvement is no more than
299                                      EpsF.
300                              *  2    relative step is no more than EpsX.
301                              *  4    gradient norm is no more than EpsG
302                              *  5    MaxIts steps was taken
303                              *  7    stopping conditions are too stringent,
304                                      further improvement is impossible
305                          * Rep.IterationsCount contains iterations count
306                          * NFEV countains number of function calculations
307         */
308        if (report.terminationtype < 0) { nmse = 10.0; return; }
309      }
310
311      // perform evaluation for optimal theta to get quality value
312      double[] grad = new double[optTheta.Length];
313      nmse = double.NaN;
314      EvaluateObjectiveAndGradient(optTheta, ref nmse, grad,
315        new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver });
316      if (double.IsNaN(nmse) || double.IsInfinity(nmse)) { nmse = 10.0; return; } // return a large value (TODO: be consistent by using NMSE)
317    }
318
319    private static void EvaluateObjectiveAndGradient(double[] x, ref double f, double[] grad, object obj) {
320      var trees = (ISymbolicExpressionTree[])((object[])obj)[0];
321      var targetVariables = (string[])((object[])obj)[1];
322      var problemData = (IRegressionProblemData)((object[])obj)[2];
323      var targetValues = (double[,])((object[])obj)[3];
324      var episodes = (IntRange[])((object[])obj)[4];
325      var numericIntegrationSteps = (int)((object[])obj)[5];
326      var latentVariables = (string[])((object[])obj)[6];
327      var odeSolver = (string)((object[])obj)[7];
328
329
330      Tuple<double, Vector>[][] predicted = null; // one array of predictions for each episode
331      // TODO: stop integration early for diverging solutions
332      predicted = Integrate(
333          trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
334          problemData.Dataset,
335          problemData.AllowedInputVariables.ToArray(),
336          targetVariables,
337          latentVariables,
338          episodes,
339          x,
340          odeSolver,
341          numericIntegrationSteps).ToArray();
342
343      if (predicted.Length != targetValues.GetLength(0)) {
344        f = 10.0; // TODO
345        Array.Clear(grad, 0, grad.Length);
346        return;
347      }
348
349      // for normalized MSE = 1/variance(t) * MSE(t, pred)
350      // TODO: Perf. (by standardization of target variables before evaluation of all trees)     
351      // var invVar = Enumerable.Range(0, targetVariables.Length)
352      //   .Select(c => Enumerable.Range(0, targetValues.GetLength(0)).Select(row => targetValues[row, c])) // column vectors
353      //   .Select(vec => vec.StandardDeviation()) // TODO: variance of stddev
354      //   .Select(v => 1.0 / v)
355      //   .ToArray();
356
357      double[] invVar = Enumerable.Repeat(1.0, targetVariables.Length).ToArray();
358
359
360      // objective function is NMSE
361      f = 0.0;
362      int n = predicted.Length;
363      double invN = 1.0 / n;
364      var g = Vector.Zero;
365      int r = 0;
366      foreach (var y_pred in predicted) {
367        // y_pred contains the predicted values for target variables first and then predicted values for latent variables
368        for (int c = 0; c < targetVariables.Length; c++) {
369
370          var y_pred_f = y_pred[c].Item1;
371          var y = targetValues[r, c];
372
373          var res = (y - y_pred_f);
374          var ressq = res * res;
375          f += ressq * invN * invVar[c] /* * Math.Exp(-0.2 * r) */ ;
376          g += -2.0 * res * y_pred[c].Item2 * invN * invVar[c] /* * Math.Exp(-0.2 * r) */;
377        }
378        r++;
379      }
380
381      g.CopyTo(grad);
382    }
383
384    public override void Analyze(Individual[] individuals, double[] qualities, ResultCollection results, IRandom random) {
385      base.Analyze(individuals, qualities, results, random);
386
387      if (!results.ContainsKey("Prediction (training)")) {
388        results.Add(new Result("Prediction (training)", typeof(ReadOnlyItemList<DataTable>)));
389      }
390      if (!results.ContainsKey("Prediction (test)")) {
391        results.Add(new Result("Prediction (test)", typeof(ReadOnlyItemList<DataTable>)));
392      }
393      if (!results.ContainsKey("Models")) {
394        results.Add(new Result("Models", typeof(VariableCollection)));
395      }
396      if (!results.ContainsKey("SNMSE")) {
397        results.Add(new Result("SNMSE", typeof(DoubleValue)));
398      }
399      if (!results.ContainsKey("Solution")) {
400        results.Add(new Result("Solution", typeof(Solution)));
401      }
402      if (!results.ContainsKey("Squared error and gradient")) {
403        results.Add(new Result("Squared error and gradient", typeof(DataTable)));
404      }
405
406      var bestIndividualAndQuality = this.GetBestIndividual(individuals, qualities);
407      var trees = bestIndividualAndQuality.Item1.Values.Select(v => v.Value).OfType<ISymbolicExpressionTree>().ToArray(); // extract all trees from individual
408
409      results["SNMSE"].Value = new DoubleValue(bestIndividualAndQuality.Item2);
410
411      var problemData = ProblemData;
412      var targetVars = TargetVariables.CheckedItems.OrderBy(i => i.Index).Select(i => i.Value.Value).ToArray();
413      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
414
415      var trainingList = new ItemList<DataTable>();
416
417      if (OptimizeParametersForEpisodes) {
418        var eIdx = 0;
419        var trainingPredictions = new List<Tuple<double, Vector>[][]>();
420        foreach (var episode in TrainingEpisodes) {
421          var episodes = new[] { episode };
422          var optTheta = ((DoubleArray)bestIndividualAndQuality.Item1["OptTheta_" + eIdx]).ToArray(); // see evaluate
423          var trainingPrediction = Integrate(
424                                   trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
425                                   problemData.Dataset,
426                                   problemData.AllowedInputVariables.ToArray(),
427                                   targetVars,
428                                   latentVariables,
429                                   episodes,
430                                   optTheta,
431                                   OdeSolver,
432                                   NumericIntegrationSteps).ToArray();
433          trainingPredictions.Add(trainingPrediction);
434          eIdx++;
435        }
436
437        // only for target values
438        var trainingRows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start));
439        for (int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
440          var targetVar = targetVars[colIdx];
441          var trainingDataTable = new DataTable(targetVar + " prediction (training)");
442          var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, trainingRows));
443          var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, trainingPredictions.SelectMany(arr => arr.Select(row => row[colIdx].Item1)).ToArray());
444          trainingDataTable.Rows.Add(actualValuesRow);
445          trainingDataTable.Rows.Add(predictedValuesRow);
446          trainingList.Add(trainingDataTable);
447        }
448        results["Prediction (training)"].Value = trainingList.AsReadOnly();
449
450
451        var models = new VariableCollection();
452
453        foreach (var tup in targetVars.Zip(trees, Tuple.Create)) {
454          var targetVarName = tup.Item1;
455          var tree = tup.Item2;
456
457          var origTreeVar = new HeuristicLab.Core.Variable(targetVarName + "(original)");
458          origTreeVar.Value = (ISymbolicExpressionTree)tree.Clone();
459          models.Add(origTreeVar);
460        }
461        results["Models"].Value = models;
462      } else {
463        var optTheta = ((DoubleArray)bestIndividualAndQuality.Item1["OptTheta"]).ToArray(); // see evaluate
464        var trainingPrediction = Integrate(
465                                   trees,  // we assume trees contain expressions for the change of each target variable over time dy/dt
466                                   problemData.Dataset,
467                                   problemData.AllowedInputVariables.ToArray(),
468                                   targetVars,
469                                   latentVariables,
470                                   TrainingEpisodes,
471                                   optTheta,
472                                   OdeSolver,
473                                   NumericIntegrationSteps).ToArray();
474        // for target values and latent variables
475        var trainingRows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start));
476        for (int colIdx = 0; colIdx < trees.Length; colIdx++) {
477          // is target variable
478          if (colIdx < targetVars.Length) {
479            var targetVar = targetVars[colIdx];
480            var trainingDataTable = new DataTable(targetVar + " prediction (training)");
481            var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, trainingRows));
482            var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, trainingPrediction.Select(arr => arr[colIdx].Item1).ToArray());
483            trainingDataTable.Rows.Add(actualValuesRow);
484            trainingDataTable.Rows.Add(predictedValuesRow);
485
486            for (int paramIdx = 0; paramIdx < optTheta.Length; paramIdx++) {
487              var paramSensitivityRow = new DataRow($"∂{targetVar}/∂θ{paramIdx}", $"Sensitivities of parameter {paramIdx}", trainingPrediction.Select(arr => arr[colIdx].Item2[paramIdx]).ToArray());
488              paramSensitivityRow.VisualProperties.SecondYAxis = true;
489              trainingDataTable.Rows.Add(paramSensitivityRow);
490            }
491            trainingList.Add(trainingDataTable);
492          } else {
493            var latentVar = latentVariables[colIdx - targetVars.Length];
494            var trainingDataTable = new DataTable(latentVar + " prediction (training)");
495            var predictedValuesRow = new DataRow(latentVar + " pred.", "Predicted values for " + latentVar, trainingPrediction.Select(arr => arr[colIdx].Item1).ToArray());
496            var emptyRow = new DataRow(latentVar);
497            trainingDataTable.Rows.Add(emptyRow);
498            trainingDataTable.Rows.Add(predictedValuesRow);
499            trainingList.Add(trainingDataTable);
500          }
501        }
502
503        var errorTable = new DataTable("Squared error and gradient");
504        var seRow = new DataRow("Squared error");
505        var gradientRows = Enumerable.Range(0, optTheta.Length).Select(i => new DataRow($"∂SE/∂θ{i}")).ToArray();
506        errorTable.Rows.Add(seRow);
507        foreach (var gRow in gradientRows) {
508          gRow.VisualProperties.SecondYAxis = true;
509          errorTable.Rows.Add(gRow);
510        }
511        var targetValues = targetVars.Select(v => problemData.Dataset.GetDoubleValues(v, trainingRows).ToArray()).ToArray();
512        int r = 0;
513        double invN = 1.0 / trainingRows.Count();
514        foreach (var y_pred in trainingPrediction) {
515          // calculate objective function gradient
516          double f_i = 0.0;
517          Vector g_i = Vector.CreateNew(new double[optTheta.Length]);
518          for (int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
519            var y_pred_f = y_pred[colIdx].Item1;
520            var y = targetValues[colIdx][r];
521
522            var res = (y - y_pred_f);
523            var ressq = res * res;
524            f_i += ressq * invN /* * Math.Exp(-0.2 * r) */;
525            g_i = g_i - 2.0 * res * y_pred[colIdx].Item2 * invN /* * Math.Exp(-0.2 * r)*/;
526          }
527          seRow.Values.Add(f_i);
528          for (int j = 0; j < g_i.Length; j++) gradientRows[j].Values.Add(g_i[j]);
529          r++;
530        }
531        results["Squared error and gradient"].Value = errorTable;
532
533        // TODO: DRY for training and test
534        var testList = new ItemList<DataTable>();
535        var testRows = ProblemData.TestIndices.ToArray();
536        var testPrediction = Integrate(
537         trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
538         problemData.Dataset,
539         problemData.AllowedInputVariables.ToArray(),
540         targetVars,
541         latentVariables,
542         new IntRange[] { ProblemData.TestPartition },
543         optTheta,
544         OdeSolver,
545         NumericIntegrationSteps).ToArray();
546
547        for (int colIdx = 0; colIdx < trees.Length; colIdx++) {
548          // is target variable
549          if (colIdx < targetVars.Length) {
550            var targetVar = targetVars[colIdx];
551            var testDataTable = new DataTable(targetVar + " prediction (test)");
552            var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, testRows));
553            var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, testPrediction.Select(arr => arr[colIdx].Item1).ToArray());
554            testDataTable.Rows.Add(actualValuesRow);
555            testDataTable.Rows.Add(predictedValuesRow);
556            testList.Add(testDataTable);
557
558          } else {
559            var latentVar = latentVariables[colIdx - targetVars.Length];
560            var testDataTable = new DataTable(latentVar + " prediction (test)");
561            var predictedValuesRow = new DataRow(latentVar + " pred.", "Predicted values for " + latentVar, testPrediction.Select(arr => arr[colIdx].Item1).ToArray());
562            var emptyRow = new DataRow(latentVar);
563            testDataTable.Rows.Add(emptyRow);
564            testDataTable.Rows.Add(predictedValuesRow);
565            testList.Add(testDataTable);
566          }
567        }
568
569        results["Prediction (training)"].Value = trainingList.AsReadOnly();
570        results["Prediction (test)"].Value = testList.AsReadOnly();
571
572
573        #region simplification of models
574        // TODO the dependency of HeuristicLab.Problems.DataAnalysis.Symbolic is not ideal
575        var models = new VariableCollection();    // to store target var names and original version of tree
576
577        var optimizedTrees = new List<ISymbolicExpressionTree>();
578        int nextParIdx = 0;
579        for (int idx = 0; idx < trees.Length; idx++) {
580          var tree = trees[idx];
581          optimizedTrees.Add(new SymbolicExpressionTree(FixParameters(tree.Root, optTheta.ToArray(), ref nextParIdx)));
582        }
583        var ds = problemData.Dataset;
584        var newVarNames = Enumerable.Range(0, nextParIdx).Select(i => "c_" + i).ToArray();
585        var allVarNames = ds.DoubleVariables.Concat(newVarNames);
586        var newVarValues = Enumerable.Range(0, nextParIdx).Select(i => "c_" + i).ToArray();
587        var allVarValues = ds.DoubleVariables.Select(varName => ds.GetDoubleValues(varName).ToList())
588          .Concat(Enumerable.Range(0, nextParIdx).Select(i => Enumerable.Repeat(optTheta[i], ds.Rows).ToList()))
589          .ToList();
590        var newDs = new Dataset(allVarNames, allVarValues);
591        var newProblemData = new RegressionProblemData(newDs, problemData.AllowedInputVariables.Concat(newVarValues).ToArray(), problemData.TargetVariable);
592        results["Solution"].Value = new Solution(optimizedTrees.ToArray(),
593                   // optTheta,
594                   newProblemData,
595                   targetVars,
596                   latentVariables,
597                   TrainingEpisodes,
598                   OdeSolver,
599                   NumericIntegrationSteps);
600
601
602        nextParIdx = 0;
603        for (int idx = 0; idx < trees.Length; idx++) {
604          var varName = string.Empty;
605          if (idx < targetVars.Length) {
606            varName = targetVars[idx];
607          } else {
608            varName = latentVariables[idx - targetVars.Length];
609          }
610          var tree = trees[idx];
611
612          // when we reference HeuristicLab.Problems.DataAnalysis.Symbolic we can translate symbols
613          var shownTree = new SymbolicExpressionTree(TranslateTreeNode(tree.Root, optTheta.ToArray(),
614            ref nextParIdx));
615
616          // var shownTree = (SymbolicExpressionTree)tree.Clone();
617          // var constantsNodeOrig = tree.IterateNodesPrefix().Where(IsConstantNode);
618          // var constantsNodeShown = shownTree.IterateNodesPrefix().Where(IsConstantNode);
619          //
620          // foreach (var n in constantsNodeOrig.Zip(constantsNodeShown, (original, shown) => new { original, shown })) {
621          //   double constantsVal = optTheta[nodeIdx[n.original]];
622          //
623          //   ConstantTreeNode replacementNode = new ConstantTreeNode(new Constant()) { Value = constantsVal };
624          //
625          //   var parentNode = n.shown.Parent;
626          //   int replacementIndex = parentNode.IndexOfSubtree(n.shown);
627          //   parentNode.RemoveSubtree(replacementIndex);
628          //   parentNode.InsertSubtree(replacementIndex, replacementNode);
629          // }
630
631          var origTreeVar = new HeuristicLab.Core.Variable(varName + "(original)");
632          origTreeVar.Value = (ISymbolicExpressionTree)tree.Clone();
633          models.Add(origTreeVar);
634          var simplifiedTreeVar = new HeuristicLab.Core.Variable(varName + "(simplified)");
635          simplifiedTreeVar.Value = TreeSimplifier.Simplify(shownTree);
636          models.Add(simplifiedTreeVar);
637
638        }
639
640        results["Models"].Value = models;
641        #endregion
642      }
643    }
644
645
646    #region interpretation
647
648    // the following uses auto-diff to calculate the gradient w.r.t. the parameters forward in time.
649    // 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.
650
651    // a comparison of three potential calculation methods for the gradient is given in:
652    // 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
653    // "Our comparison establishes that the adjoint method is computationally more efficient for numerical estimation of parametric gradients
654    // for state-space models — both linear and non-linear, as in the case of a dynamical causal model (DCM)"
655
656    // for a solver with the necessary features see: https://computation.llnl.gov/projects/sundials/cvodes
657
658    public static IEnumerable<Tuple<double, Vector>[]> Integrate(
659      ISymbolicExpressionTree[] trees, IDataset dataset,
660      string[] inputVariables, string[] targetVariables, string[] latentVariables, IEnumerable<IntRange> episodes,
661      double[] parameterValues,
662      string odeSolver, int numericIntegrationSteps = 100) {
663
664      // TODO: numericIntegrationSteps is only relevant for the HeuristicLab solver
665      var episodeIdx = 0;
666      foreach (var episode in episodes) {
667        var rows = Enumerable.Range(episode.Start, episode.End - episode.Start);
668
669        // integrate forward starting with known values for the target in t0
670
671        var variableValues = new Dictionary<string, Tuple<double, Vector>>();
672        var t0 = rows.First();
673        foreach (var varName in inputVariables) {
674          variableValues.Add(varName, Tuple.Create(dataset.GetDoubleValue(varName, t0), Vector.Zero));
675        }
676        foreach (var varName in targetVariables) {
677          variableValues.Add(varName, Tuple.Create(dataset.GetDoubleValue(varName, t0), Vector.Zero));
678        }
679        // add value entries for latent variables which are also integrated
680        // initial values are at the end of the parameter vector
681        // separete initial values for each episode
682        var initialValueIdx = parameterValues.Length - episodes.Count() * latentVariables.Length + episodeIdx * latentVariables.Length;
683        foreach (var latentVar in latentVariables) {
684          var arr = new double[parameterValues.Length]; // backing array
685          arr[initialValueIdx] = 1.0;
686          var g = new Vector(arr);
687          variableValues.Add(latentVar,
688            Tuple.Create(parameterValues[initialValueIdx], g)); // we don't have observations for latent variables therefore we optimize the initial value for each episode
689          initialValueIdx++;
690        }
691
692        var calculatedVariables = targetVariables.Concat(latentVariables).ToArray(); // TODO: must conincide with the order of trees in the encoding
693
694        // return first value as stored in the dataset
695        yield return calculatedVariables
696          .Select(calcVarName => variableValues[calcVarName])
697          .ToArray();
698
699        var prevT = rows.First(); // TODO: here we should use a variable for t if it is available. Right now we assume equidistant measurements.
700        foreach (var t in rows.Skip(1)) {
701          if (odeSolver == "HeuristicLab")
702            IntegrateHL(trees, calculatedVariables, variableValues, parameterValues, numericIntegrationSteps);
703          else if (odeSolver == "CVODES")
704            throw new NotImplementedException();
705          // IntegrateCVODES(trees, calculatedVariables, variableValues, parameterValues, t - prevT);
706          else throw new InvalidOperationException("Unknown ODE solver " + odeSolver);
707          prevT = t;
708
709          // This check doesn't work with the HeuristicLab integrator if there are input variables
710          //if (variableValues.Count == targetVariables.Length) {
711          // only return the target variables for calculation of errors
712          var res = calculatedVariables
713            .Select(targetVar => variableValues[targetVar])
714            .ToArray();
715          if (res.Any(ri => double.IsNaN(ri.Item1) || double.IsInfinity(ri.Item1))) yield break;
716          yield return res;
717          //} else {
718          //  yield break; // stop early on problems
719          //}
720
721
722          // update for next time step
723          foreach (var varName in inputVariables) {
724            variableValues[varName] = Tuple.Create(dataset.GetDoubleValue(varName, t), Vector.Zero);
725          }
726        }
727        episodeIdx++;
728      }
729    }
730
731    #region CVODES
732
733    /*
734    /// <summary>
735    ///  Here we use CVODES to solve the ODE. Forward sensitivities are used to calculate the gradient for parameter optimization
736    /// </summary>
737    /// <param name="trees">Each equation in the ODE represented as a tree</param>
738    /// <param name="calculatedVariables">The names of the calculated variables</param>
739    /// <param name="variableValues">The start values of the calculated variables as well as their sensitivites over parameters</param>
740    /// <param name="parameterValues">The current parameter values</param>
741    /// <param name="t">The time t up to which we need to integrate.</param>
742    private static void IntegrateCVODES(
743      ISymbolicExpressionTree[] trees, // f(y,p) in tree representation
744      string[] calculatedVariables, // names of elements of y
745      Dictionary<string, Tuple<double, Vector>> variableValues,  //  y (intput and output) input: y(t0), output: y(t0+t)
746      double[] parameterValues, // p
747      double t // duration t for which we want to integrate
748      ) {
749
750      // the RHS of the ODE
751      // dy/dt = f(y_t,x_t,p)
752      CVODES.CVRhsFunc f = CreateOdeRhs(trees, calculatedVariables, parameterValues);
753      // the Jacobian ∂f/∂y
754      CVODES.CVDlsJacFunc jac = CreateJac(trees, calculatedVariables, parameterValues);
755
756      // the RHS for the forward sensitivities (∂f/∂y)s_i(t) + ∂f/∂p_i
757      CVODES.CVSensRhsFn sensF = CreateSensitivityRhs(trees, calculatedVariables, parameterValues);
758
759      // setup solver
760      int numberOfEquations = trees.Length;
761      IntPtr y = IntPtr.Zero;
762      IntPtr cvode_mem = IntPtr.Zero;
763      IntPtr A = IntPtr.Zero;
764      IntPtr yS0 = IntPtr.Zero;
765      IntPtr linearSolver = IntPtr.Zero;
766      var ns = parameterValues.Length; // number of parameters
767
768      try {
769        y = CVODES.N_VNew_Serial(numberOfEquations);
770        // init y to current values of variables
771        // y must be initialized before calling CVodeInit
772        for (int i = 0; i < calculatedVariables.Length; i++) {
773          CVODES.NV_Set_Ith_S(y, i, variableValues[calculatedVariables[i]].Item1);
774        }
775
776        cvode_mem = CVODES.CVodeCreate(CVODES.MultistepMethod.CV_ADAMS, CVODES.NonlinearSolverIteration.CV_FUNCTIONAL);
777
778        var flag = CVODES.CVodeInit(cvode_mem, f, 0.0, y);
779        Debug.Assert(CVODES.CV_SUCCESS == flag);
780
781        double relTol = 1.0e-2;
782        double absTol = 1.0;
783        flag = CVODES.CVodeSStolerances(cvode_mem, relTol, absTol);  // TODO: probably need to adjust absTol per variable
784        Debug.Assert(CVODES.CV_SUCCESS == flag);
785
786        A = CVODES.SUNDenseMatrix(numberOfEquations, numberOfEquations);
787        Debug.Assert(A != IntPtr.Zero);
788
789        linearSolver = CVODES.SUNDenseLinearSolver(y, A);
790        Debug.Assert(linearSolver != IntPtr.Zero);
791
792        flag = CVODES.CVDlsSetLinearSolver(cvode_mem, linearSolver, A);
793        Debug.Assert(CVODES.CV_SUCCESS == flag);
794
795        flag = CVODES.CVDlsSetJacFn(cvode_mem, jac);
796        Debug.Assert(CVODES.CV_SUCCESS == flag);
797
798        yS0 = CVODES.N_VCloneVectorArray_Serial(ns, y); // clone the output vector for each parameter
799        unsafe {
800          // set to initial sensitivities supplied by caller
801          for (int pIdx = 0; pIdx < ns; pIdx++) {
802            var yS0_i = *((IntPtr*)yS0.ToPointer() + pIdx);
803            for (var varIdx = 0; varIdx < calculatedVariables.Length; varIdx++) {
804              CVODES.NV_Set_Ith_S(yS0_i, varIdx, variableValues[calculatedVariables[varIdx]].Item2[pIdx]);
805            }
806          }
807        }
808
809        flag = CVODES.CVodeSensInit(cvode_mem, ns, CVODES.CV_SIMULTANEOUS, sensF, yS0);
810        Debug.Assert(CVODES.CV_SUCCESS == flag);
811
812        flag = CVODES.CVodeSensEEtolerances(cvode_mem);
813        Debug.Assert(CVODES.CV_SUCCESS == flag);
814
815        // make one forward integration step
816        double tout = 0.0; // first output time
817        flag = CVODES.CVode(cvode_mem, t, y, ref tout, CVODES.CV_NORMAL);
818        if (flag == CVODES.CV_SUCCESS) {
819          Debug.Assert(t == tout);
820
821          // get sensitivities
822          flag = CVODES.CVodeGetSens(cvode_mem, ref tout, yS0);
823          Debug.Assert(CVODES.CV_SUCCESS == flag);
824
825          // update variableValues based on integration results
826          for (int varIdx = 0; varIdx < calculatedVariables.Length; varIdx++) {
827            var yi = CVODES.NV_Get_Ith_S(y, varIdx);
828            var gArr = new double[parameterValues.Length];
829            for (var pIdx = 0; pIdx < parameterValues.Length; pIdx++) {
830              unsafe {
831                var yS0_pi = *((IntPtr*)yS0.ToPointer() + pIdx);
832                gArr[pIdx] = CVODES.NV_Get_Ith_S(yS0_pi, varIdx);
833              }
834            }
835            variableValues[calculatedVariables[varIdx]] = Tuple.Create(yi, new Vector(gArr));
836          }
837        } else {
838          variableValues.Clear();   // indicate problems by not returning new values
839        }
840
841        // cleanup all allocated objects
842      } finally {
843        if (y != IntPtr.Zero) CVODES.N_VDestroy_Serial(y);
844        if (cvode_mem != IntPtr.Zero) CVODES.CVodeFree(ref cvode_mem);
845        if (linearSolver != IntPtr.Zero) CVODES.SUNLinSolFree(linearSolver);
846        if (A != IntPtr.Zero) CVODES.SUNMatDestroy(A);
847        if (yS0 != IntPtr.Zero) CVODES.N_VDestroyVectorArray_Serial(yS0, ns);
848      }
849    }
850
851
852    private static CVODES.CVRhsFunc CreateOdeRhs(
853      ISymbolicExpressionTree[] trees,
854      string[] calculatedVariables,
855      double[] parameterValues) {
856      // we don't need to calculate a gradient here
857      return (double t,
858              IntPtr y, // N_Vector, current value of y (input)
859              IntPtr ydot, // N_Vector (calculated value of y' (output)
860              IntPtr user_data // optional user data, (unused here)
861              ) => {
862                // TODO: perf
863                var nodeValues = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
864
865                int pIdx = 0;
866                foreach (var tree in trees) {
867                  foreach (var n in tree.IterateNodesPrefix()) {
868                    if (IsConstantNode(n)) {
869                      nodeValues.Add(n, Tuple.Create(parameterValues[pIdx], Vector.Zero)); // here we do not need a gradient
870                      pIdx++;
871                    } else if (n.SubtreeCount == 0) {
872                      // for variables and latent variables get the value from variableValues
873                      var varName = n.Symbol.Name;
874                      var varIdx = Array.IndexOf(calculatedVariables, varName); // TODO: perf!
875                      if (varIdx < 0) throw new InvalidProgramException();
876                      var y_i = CVODES.NV_Get_Ith_S(y, (long)varIdx);
877                      nodeValues.Add(n, Tuple.Create(y_i, Vector.Zero)); // no gradient needed
878                    }
879                  }
880                }
881                for (int i = 0; i < trees.Length; i++) {
882                  var tree = trees[i];
883                  var res_i = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues);
884                  CVODES.NV_Set_Ith_S(ydot, i, res_i.Item1);
885                }
886                return 0;
887              };
888    }
889
890    private static CVODES.CVDlsJacFunc CreateJac(
891      ISymbolicExpressionTree[] trees,
892      string[] calculatedVariables,
893      double[] parameterValues) {
894
895      return (
896        double t, // current time (input)
897        IntPtr y, // N_Vector, current value of y (input)
898        IntPtr fy, // N_Vector, current value of f (input)
899        IntPtr Jac, // SUNMatrix ∂f/∂y (output, rows i contains are ∂f_i/∂y vector)
900        IntPtr user_data, // optional (unused here)
901        IntPtr tmp1, // N_Vector, optional (unused here)
902        IntPtr tmp2, // N_Vector, optional (unused here)
903        IntPtr tmp3 // N_Vector, optional (unused here)
904      ) => {
905        // here we need to calculate partial derivatives for the calculated variables y
906        var nodeValues = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
907        int pIdx = 0;
908        foreach (var tree in trees) {
909          foreach (var n in tree.IterateNodesPrefix()) {
910            if (IsConstantNode(n)) {
911              nodeValues.Add(n, Tuple.Create(parameterValues[pIdx], Vector.Zero)); // here we need a gradient over y which is zero for parameters
912              pIdx++;
913            } else if (n.SubtreeCount == 0) {
914              // for variables and latent variables we use values supplied in y and init gradient vectors accordingly
915              var varName = n.Symbol.Name;
916              var varIdx = Array.IndexOf(calculatedVariables, varName); // TODO: perf!
917              if (varIdx < 0) throw new InvalidProgramException();
918
919              var y_i = CVODES.NV_Get_Ith_S(y, (long)varIdx);
920              var gArr = new double[CVODES.NV_LENGTH_S(y)]; // backing array
921              gArr[varIdx] = 1.0;
922              var g = new Vector(gArr);
923              nodeValues.Add(n, Tuple.Create(y_i, g));
924            }
925          }
926        }
927
928        for (int i = 0; i < trees.Length; i++) {
929          var tree = trees[i];
930          var res = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues);
931          var g = res.Item2;
932          for (int j = 0; j < calculatedVariables.Length; j++) {
933            CVODES.SUNDenseMatrix_Set(Jac, i, j, g[j]);
934          }
935        }
936        return 0; // on success
937      };
938    }
939
940
941    // to calculate sensitivities RHS for all equations at once
942    // must compute (∂f/∂y)s_i(t) + ∂f/∂p_i and store in ySdot.
943    // Index i refers to parameters, dimensionality of matrix and vectors is number of equations
944    private static CVODES.CVSensRhsFn CreateSensitivityRhs(ISymbolicExpressionTree[] trees, string[] calculatedVariables, double[] parameterValues) {
945      return (
946              int Ns, // number of parameters
947              double t, // current time
948              IntPtr y, // N_Vector y(t) (input)
949              IntPtr ydot, // N_Vector dy/dt(t) (input)
950              IntPtr yS, // N_Vector*, one vector for each parameter (input)
951              IntPtr ySdot, // N_Vector*, one vector for each parameter (output)
952              IntPtr user_data, // optional (unused here)
953              IntPtr tmp1, // N_Vector, optional (unused here)
954              IntPtr tmp2 // N_Vector, optional (unused here)
955        ) => {
956          // here we need to calculate partial derivatives for the calculated variables y as well as for the parameters
957          var nodeValues = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
958          var d = calculatedVariables.Length + parameterValues.Length; // dimensionality of gradient
959          // first collect variable values
960          foreach (var tree in trees) {
961            foreach (var n in tree.IterateNodesPrefix()) {
962              if (IsVariableNode(n)) {
963                // for variables and latent variables we use values supplied in y and init gradient vectors accordingly
964                var varName = n.Symbol.Name;
965                var varIdx = Array.IndexOf(calculatedVariables, varName); // TODO: perf!
966                if (varIdx < 0) throw new InvalidProgramException();
967
968                var y_i = CVODES.NV_Get_Ith_S(y, (long)varIdx);
969                var gArr = new double[d]; // backing array
970                gArr[varIdx] = 1.0;
971                var g = new Vector(gArr);
972                nodeValues.Add(n, Tuple.Create(y_i, g));
973              }
974            }
975          }
976          // then collect constants
977          int pIdx = 0;
978          foreach (var tree in trees) {
979            foreach (var n in tree.IterateNodesPrefix()) {
980              if (IsConstantNode(n)) {
981                var gArr = new double[d];
982                gArr[calculatedVariables.Length + pIdx] = 1.0;
983                var g = new Vector(gArr);
984                nodeValues.Add(n, Tuple.Create(parameterValues[pIdx], g));
985                pIdx++;
986              }
987            }
988          }
989          // gradient vector is [∂f/∂y_1, ∂f/∂y_2, ... ∂f/∂yN, ∂f/∂p_1 ... ∂f/∂p_K]
990
991
992          for (pIdx = 0; pIdx < Ns; pIdx++) {
993            unsafe {
994              var sDot_pi = *((IntPtr*)ySdot.ToPointer() + pIdx);
995              CVODES.N_VConst_Serial(0.0, sDot_pi);
996            }
997          }
998
999          for (int i = 0; i < trees.Length; i++) {
1000            var tree = trees[i];
1001            var res = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues);
1002            var g = res.Item2;
1003
1004
1005            // update ySdot = (∂f/∂y)s_i(t) + ∂f/∂p_i
1006
1007            for (pIdx = 0; pIdx < Ns; pIdx++) {
1008              unsafe {
1009                var sDot_pi = *((IntPtr*)ySdot.ToPointer() + pIdx);
1010                var s_pi = *((IntPtr*)yS.ToPointer() + pIdx);
1011
1012                var v = CVODES.NV_Get_Ith_S(sDot_pi, i);
1013                // (∂f/∂y)s_i(t)
1014                var p = 0.0;
1015                for (int yIdx = 0; yIdx < calculatedVariables.Length; yIdx++) {
1016                  p += g[yIdx] * CVODES.NV_Get_Ith_S(s_pi, yIdx);
1017                }
1018                // + ∂f/∂p_i
1019                CVODES.NV_Set_Ith_S(sDot_pi, i, v + p + g[calculatedVariables.Length + pIdx]);
1020              }
1021            }
1022
1023          }
1024          return 0; // on success
1025        };
1026    }
1027    */
1028    #endregion
1029
1030    private static void IntegrateHL(
1031      ISymbolicExpressionTree[] trees,
1032      string[] calculatedVariables, // names of integrated variables
1033      Dictionary<string, Tuple<double, Vector>> variableValues, //y (intput and output) input: y(t0), output: y(t0+1)
1034      double[] parameterValues,
1035      int numericIntegrationSteps) {
1036
1037      var nodeValues = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
1038
1039      // the nodeValues for parameters are constant over time
1040      // TODO: this needs to be done only once for each iteration in gradient descent (whenever parameter values change)
1041      // NOTE: the order must be the same as above (prefix order for nodes)
1042      int pIdx = 0;
1043      foreach (var tree in trees) {
1044        foreach (var node in tree.Root.IterateNodesPrefix()) {
1045          if (IsConstantNode(node)) {
1046            var gArr = new double[parameterValues.Length]; // backing array
1047            gArr[pIdx] = 1.0;
1048            var g = new Vector(gArr);
1049            nodeValues.Add(node, new Tuple<double, Vector>(parameterValues[pIdx], g));
1050            pIdx++;
1051          } else if (node.SubtreeCount == 0) {
1052            // for (latent) variables get the values from variableValues
1053            var varName = node.Symbol.Name;
1054            nodeValues.Add(node, variableValues[varName]);
1055          }
1056        }
1057      }
1058
1059      double[] deltaF = new double[calculatedVariables.Length];
1060      Vector[] deltaG = new Vector[calculatedVariables.Length];
1061
1062      double h = 1.0 / numericIntegrationSteps;
1063      for (int step = 0; step < numericIntegrationSteps; step++) {
1064        //var deltaValues = new Dictionary<string, Tuple<double, Vector>>();
1065        for (int i = 0; i < trees.Length; i++) {
1066          var tree = trees[i];
1067          var targetVarName = calculatedVariables[i];
1068
1069          // Root.GetSubtree(0).GetSubtree(0) skips programRoot and startSymbol
1070          double f; Vector g;
1071          InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues, out f, out g);
1072          deltaF[i] = f;
1073          deltaG[i] = g;
1074        }
1075
1076        // update variableValues for next step, trapezoid integration
1077        for (int i = 0; i < trees.Length; i++) {
1078          var varName = calculatedVariables[i];
1079          var oldVal = variableValues[varName];
1080          var newVal = Tuple.Create(
1081            oldVal.Item1 + h * deltaF[i],
1082            oldVal.Item2 + deltaG[i].Scale(h)
1083          );
1084          variableValues[varName] = newVal;
1085        }
1086
1087        // TODO perf
1088        foreach (var node in nodeValues.Keys.ToArray()) {
1089          if (node.SubtreeCount == 0 && !IsConstantNode(node)) {
1090            // update values for (latent) variables
1091            var varName = node.Symbol.Name;
1092            nodeValues[node] = variableValues[varName];
1093          }
1094        }
1095      }
1096    }
1097
1098    private static void InterpretRec(
1099      ISymbolicExpressionTreeNode node,
1100      Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>> nodeValues,      // contains value and gradient vector for a node (variables and constants only)
1101      out double f,
1102      out Vector g
1103      ) {
1104      double fl, fr;
1105      Vector gl, gr;
1106      switch (node.Symbol.Name) {
1107        case "+": {
1108            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
1109            InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
1110            f = fl + fr;
1111            g = Vector.AddTo(gl, gr);
1112            break;
1113          }
1114        case "*": {
1115            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
1116            InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
1117            f = fl * fr;
1118            g = Vector.AddTo(gl.Scale(fr), gr.Scale(fl)); // f'*g + f*g'
1119            break;
1120          }
1121
1122        case "-": {
1123            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
1124            InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
1125            f = fl - fr;
1126            g = Vector.Subtract(gl, gr);
1127            break;
1128          }
1129        case "%": {
1130            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
1131            InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
1132
1133            // protected division
1134            if (fr.IsAlmost(0.0)) {
1135              f = 0;
1136              g = Vector.Zero;
1137            } else {
1138              f = fl / fr;
1139              g = Vector.AddTo(gr.Scale(fl * -1.0 / (fr * fr)), gl.Scale(1.0 / fr)); // (f/g)' = f * (1/g)' + 1/g * f' = f * -1/g² * g' + 1/g * f'
1140            }
1141            break;
1142          }
1143        case "sin": {
1144            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
1145            f = Math.Sin(fl);
1146            g = gl.Scale(Math.Cos(fl));
1147            break;
1148          }
1149        case "cos": {
1150            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
1151            f = Math.Cos(fl);
1152            g = gl.Scale(-Math.Sin(fl));
1153            break;
1154          }
1155        case "sqr": {
1156            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
1157            f = fl * fl;
1158            g = gl.Scale(2.0 * fl);
1159            break;
1160          }
1161        default: {
1162            var t = nodeValues[node];
1163            f = t.Item1;
1164            g = Vector.CreateNew(t.Item2);
1165            break;
1166          }
1167      }
1168    }
1169    #endregion
1170
1171    #region events
1172    /*
1173     * Dependencies between parameters:
1174     *
1175     * ProblemData
1176     *    |
1177     *    V
1178     * TargetVariables   FunctionSet    MaximumLength    NumberOfLatentVariables
1179     *               |   |                 |                   |
1180     *               V   V                 |                   |
1181     *             Grammar <---------------+-------------------
1182     *                |
1183     *                V
1184     *            Encoding
1185     */
1186    private void RegisterEventHandlers() {
1187      ProblemDataParameter.ValueChanged += ProblemDataParameter_ValueChanged;
1188      if (ProblemDataParameter.Value != null) ProblemDataParameter.Value.Changed += ProblemData_Changed;
1189
1190      TargetVariablesParameter.ValueChanged += TargetVariablesParameter_ValueChanged;
1191      if (TargetVariablesParameter.Value != null) TargetVariablesParameter.Value.CheckedItemsChanged += CheckedTargetVariablesChanged;
1192
1193      FunctionSetParameter.ValueChanged += FunctionSetParameter_ValueChanged;
1194      if (FunctionSetParameter.Value != null) FunctionSetParameter.Value.CheckedItemsChanged += CheckedFunctionsChanged;
1195
1196      MaximumLengthParameter.Value.ValueChanged += MaximumLengthChanged;
1197
1198      NumberOfLatentVariablesParameter.Value.ValueChanged += NumLatentVariablesChanged;
1199    }
1200
1201    private void NumLatentVariablesChanged(object sender, EventArgs e) {
1202      UpdateGrammarAndEncoding();
1203    }
1204
1205    private void MaximumLengthChanged(object sender, EventArgs e) {
1206      UpdateGrammarAndEncoding();
1207    }
1208
1209    private void FunctionSetParameter_ValueChanged(object sender, EventArgs e) {
1210      FunctionSetParameter.Value.CheckedItemsChanged += CheckedFunctionsChanged;
1211    }
1212
1213    private void CheckedFunctionsChanged(object sender, CollectionItemsChangedEventArgs<IndexedItem<StringValue>> e) {
1214      UpdateGrammarAndEncoding();
1215    }
1216
1217    private void TargetVariablesParameter_ValueChanged(object sender, EventArgs e) {
1218      TargetVariablesParameter.Value.CheckedItemsChanged += CheckedTargetVariablesChanged;
1219    }
1220
1221    private void CheckedTargetVariablesChanged(object sender, CollectionItemsChangedEventArgs<IndexedItem<StringValue>> e) {
1222      UpdateGrammarAndEncoding();
1223    }
1224
1225    private void ProblemDataParameter_ValueChanged(object sender, EventArgs e) {
1226      ProblemDataParameter.Value.Changed += ProblemData_Changed;
1227      OnProblemDataChanged();
1228      OnReset();
1229    }
1230
1231    private void ProblemData_Changed(object sender, EventArgs e) {
1232      OnProblemDataChanged();
1233      OnReset();
1234    }
1235
1236    private void OnProblemDataChanged() {
1237      UpdateTargetVariables();        // implicitly updates other dependent parameters
1238      var handler = ProblemDataChanged;
1239      if (handler != null) handler(this, EventArgs.Empty);
1240    }
1241
1242    #endregion
1243
1244    #region  helper
1245
1246    private void InitAllParameters() {
1247      UpdateTargetVariables(); // implicitly updates the grammar and the encoding     
1248    }
1249
1250    private ReadOnlyCheckedItemList<StringValue> CreateFunctionSet() {
1251      var l = new CheckedItemList<StringValue>();
1252      l.Add(new StringValue("+").AsReadOnly());
1253      l.Add(new StringValue("*").AsReadOnly());
1254      l.Add(new StringValue("%").AsReadOnly());
1255      l.Add(new StringValue("-").AsReadOnly());
1256      l.Add(new StringValue("sin").AsReadOnly());
1257      l.Add(new StringValue("cos").AsReadOnly());
1258      l.Add(new StringValue("sqr").AsReadOnly());
1259      return l.AsReadOnly();
1260    }
1261
1262    private static bool IsConstantNode(ISymbolicExpressionTreeNode n) {
1263      return n.Symbol.Name[0] == 'θ';
1264    }
1265    private static bool IsLatentVariableNode(ISymbolicExpressionTreeNode n) {
1266      return n.Symbol.Name[0] == 'λ';
1267    }
1268    private static bool IsVariableNode(ISymbolicExpressionTreeNode n) {
1269      return (n.SubtreeCount == 0) && !IsConstantNode(n) && !IsLatentVariableNode(n);
1270    }
1271
1272
1273    private void UpdateTargetVariables() {
1274      var currentlySelectedVariables = TargetVariables.CheckedItems
1275        .OrderBy(i => i.Index)
1276        .Select(i => i.Value.Value)
1277        .ToArray();
1278
1279      var newVariablesList = new CheckedItemList<StringValue>(ProblemData.Dataset.VariableNames.Select(str => new StringValue(str).AsReadOnly()).ToArray()).AsReadOnly();
1280      var matchingItems = newVariablesList.Where(item => currentlySelectedVariables.Contains(item.Value)).ToArray();
1281      foreach (var item in newVariablesList) {
1282        if (currentlySelectedVariables.Contains(item.Value)) {
1283          newVariablesList.SetItemCheckedState(item, true);
1284        } else {
1285          newVariablesList.SetItemCheckedState(item, false);
1286        }
1287      }
1288      TargetVariablesParameter.Value = newVariablesList;
1289    }
1290
1291    private void UpdateGrammarAndEncoding() {
1292      var encoding = new MultiEncoding();
1293      var g = CreateGrammar();
1294      foreach (var targetVar in TargetVariables.CheckedItems) {
1295        encoding = encoding.Add(new SymbolicExpressionTreeEncoding(targetVar + "_tree", g, MaximumLength, MaximumLength)); // only limit by length
1296      }
1297      for (int i = 1; i <= NumberOfLatentVariables; i++) {
1298        encoding = encoding.Add(new SymbolicExpressionTreeEncoding("λ" + i + "_tree", g, MaximumLength, MaximumLength));
1299      }
1300      Encoding = encoding;
1301    }
1302
1303    private ISymbolicExpressionGrammar CreateGrammar() {
1304      // whenever ProblemData is changed we create a new grammar with the necessary symbols
1305      var g = new SimpleSymbolicExpressionGrammar();
1306      var unaryFunc = new string[] { "sin", "cos", "sqr" };
1307      var binaryFunc = new string[] { "+", "-", "*", "%" };
1308      foreach (var func in unaryFunc) {
1309        if (FunctionSet.CheckedItems.Any(ci => ci.Value.Value == func)) g.AddSymbol(func, 1, 1);
1310      }
1311      foreach (var func in binaryFunc) {
1312        if (FunctionSet.CheckedItems.Any(ci => ci.Value.Value == func)) g.AddSymbol(func, 2, 2);
1313      }
1314
1315      foreach (var variableName in ProblemData.AllowedInputVariables.Union(TargetVariables.CheckedItems.Select(i => i.Value.Value)))
1316        g.AddTerminalSymbol(variableName);
1317
1318      // generate symbols for numeric parameters for which the value is optimized using AutoDiff
1319      // we generate multiple symbols to balance the probability for selecting a numeric parameter in the generation of random trees
1320      var numericConstantsFactor = 2.0;
1321      for (int i = 0; i < numericConstantsFactor * (ProblemData.AllowedInputVariables.Count() + TargetVariables.CheckedItems.Count()); i++) {
1322        g.AddTerminalSymbol("θ" + i); // numeric parameter for which the value is optimized using AutoDiff
1323      }
1324
1325      // generate symbols for latent variables
1326      for (int i = 1; i <= NumberOfLatentVariables; i++) {
1327        g.AddTerminalSymbol("λ" + i); // numeric parameter for which the value is optimized using AutoDiff
1328      }
1329
1330      return g;
1331    }
1332
1333
1334
1335
1336
1337    private ISymbolicExpressionTreeNode FixParameters(ISymbolicExpressionTreeNode n, double[] parameterValues, ref int nextParIdx) {
1338      ISymbolicExpressionTreeNode translatedNode = null;
1339      if (n.Symbol is StartSymbol) {
1340        translatedNode = new StartSymbol().CreateTreeNode();
1341      } else if (n.Symbol is ProgramRootSymbol) {
1342        translatedNode = new ProgramRootSymbol().CreateTreeNode();
1343      } else if (n.Symbol.Name == "+") {
1344        translatedNode = new SimpleSymbol("+", 2).CreateTreeNode();
1345      } else if (n.Symbol.Name == "-") {
1346        translatedNode = new SimpleSymbol("-", 2).CreateTreeNode();
1347      } else if (n.Symbol.Name == "*") {
1348        translatedNode = new SimpleSymbol("*", 2).CreateTreeNode();
1349      } else if (n.Symbol.Name == "%") {
1350        translatedNode = new SimpleSymbol("%", 2).CreateTreeNode();
1351      } else if (n.Symbol.Name == "sin") {
1352        translatedNode = new SimpleSymbol("sin", 1).CreateTreeNode();
1353      } else if (n.Symbol.Name == "cos") {
1354        translatedNode = new SimpleSymbol("cos", 1).CreateTreeNode();
1355      } else if (n.Symbol.Name == "sqr") {
1356        translatedNode = new SimpleSymbol("sqr", 1).CreateTreeNode();
1357      } else if (IsConstantNode(n)) {
1358        translatedNode = new SimpleSymbol("c_" + nextParIdx, 0).CreateTreeNode();
1359        nextParIdx++;
1360      } else {
1361        translatedNode = new SimpleSymbol(n.Symbol.Name, n.SubtreeCount).CreateTreeNode();
1362      }
1363      foreach (var child in n.Subtrees) {
1364        translatedNode.AddSubtree(FixParameters(child, parameterValues, ref nextParIdx));
1365      }
1366      return translatedNode;
1367    }
1368
1369
1370    private ISymbolicExpressionTreeNode TranslateTreeNode(ISymbolicExpressionTreeNode n, double[] parameterValues, ref int nextParIdx) {
1371      ISymbolicExpressionTreeNode translatedNode = null;
1372      if (n.Symbol is StartSymbol) {
1373        translatedNode = new StartSymbol().CreateTreeNode();
1374      } else if (n.Symbol is ProgramRootSymbol) {
1375        translatedNode = new ProgramRootSymbol().CreateTreeNode();
1376      } else if (n.Symbol.Name == "+") {
1377        translatedNode = new Addition().CreateTreeNode();
1378      } else if (n.Symbol.Name == "-") {
1379        translatedNode = new Subtraction().CreateTreeNode();
1380      } else if (n.Symbol.Name == "*") {
1381        translatedNode = new Multiplication().CreateTreeNode();
1382      } else if (n.Symbol.Name == "%") {
1383        translatedNode = new Division().CreateTreeNode();
1384      } else if (n.Symbol.Name == "sin") {
1385        translatedNode = new Sine().CreateTreeNode();
1386      } else if (n.Symbol.Name == "cos") {
1387        translatedNode = new Cosine().CreateTreeNode();
1388      } else if (n.Symbol.Name == "sqr") {
1389        translatedNode = new Square().CreateTreeNode();
1390      } else if (IsConstantNode(n)) {
1391        var constNode = (ConstantTreeNode)new Constant().CreateTreeNode();
1392        constNode.Value = parameterValues[nextParIdx];
1393        nextParIdx++;
1394        translatedNode = constNode;
1395      } else {
1396        // assume a variable name
1397        var varName = n.Symbol.Name;
1398        var varNode = (VariableTreeNode)new Variable().CreateTreeNode();
1399        varNode.Weight = 1.0;
1400        varNode.VariableName = varName;
1401        translatedNode = varNode;
1402      }
1403      foreach (var child in n.Subtrees) {
1404        translatedNode.AddSubtree(TranslateTreeNode(child, parameterValues, ref nextParIdx));
1405      }
1406      return translatedNode;
1407    }
1408    #endregion
1409
1410    #region Import & Export
1411    public void Load(IRegressionProblemData data) {
1412      Name = data.Name;
1413      Description = data.Description;
1414      ProblemData = data;
1415    }
1416
1417    public IRegressionProblemData Export() {
1418      return ProblemData;
1419    }
1420    #endregion
1421
1422  }
1423}
Note: See TracBrowser for help on using the repository browser.