Changeset 16601 for branches

Timestamp:

02/13/19 10:14:03 (6 years ago)

Author:

gkronber

Message:

#2925: refactored dynamical modelling code

Location:

branches/2925_AutoDiffForDynamicalModels

Files:

• AutoDiffForDynamicalModelsTest/AutoDiffForDynamicalModelsTest.csproj (modified) (3 diffs)
• AutoDiffForDynamicalModelsTest/Genetic Algorithm (GA).hl (deleted)
• AutoDiffForDynamicalModelsTest/TestGA.cs (deleted)
• AutoDiffForDynamicalModelsTest/TestOdeIdentification.cs (added)
• HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs (modified) (23 diffs)
• HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Vector.cs (modified) (2 diffs)

branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest/AutoDiffForDynamicalModelsTest.csproj

@@ -78 +78 @@
     <Compile Include="TestCvodes.cs" />
     <Compile Include="Properties\AssemblyInfo.cs" />
-    <Compile Include="TestGA.cs" />
+    <Compile Include="TestOdeIdentification.cs" />
   </ItemGroup>
   <ItemGroup>
-    <None Include="Genetic Algorithm %28GA%29.hl">
-      <CopyToOutputDirectory>PreserveNewest</CopyToOutputDirectory>
-    </None>
     <None Include="packages.config" />
   </ItemGroup>
   <ItemGroup>
+    <ProjectReference Include="..\HeuristicLab.Algorithms.DataAnalysis\3.4\HeuristicLab.Algorithms.DataAnalysis-3.4.csproj">
+      <Project>{2E782078-FA81-4B70-B56F-74CE38DAC6C8}</Project>
+      <Name>HeuristicLab.Algorithms.DataAnalysis-3.4</Name>
+    </ProjectReference>
     <ProjectReference Include="..\HeuristicLab.Algorithms.GeneticAlgorithm\3.3\HeuristicLab.Algorithms.GeneticAlgorithm-3.3.csproj">
       <Project>{A51DA44F-CB35-4F6F-99F5-2A2E904AB93B}</Project>
       <Name>HeuristicLab.Algorithms.GeneticAlgorithm-3.3</Name>
+    </ProjectReference>
+    <ProjectReference Include="..\HeuristicLab.Collections\3.3\HeuristicLab.Collections-3.3.csproj">
+      <Project>{958B43BC-CC5C-4FA2-8628-2B3B01D890B6}</Project>
+      <Name>HeuristicLab.Collections-3.3</Name>
     </ProjectReference>
     <ProjectReference Include="..\HeuristicLab.Common\3.3\HeuristicLab.Common-3.3.csproj">
@@ -99 +104 @@
       <Name>HeuristicLab.Core-3.3</Name>
     </ProjectReference>
+    <ProjectReference Include="..\HeuristicLab.Data\3.3\HeuristicLab.Data-3.3.csproj">
+      <Project>{BBAB9DF5-5EF3-4BA8-ADE9-B36E82114937}</Project>
+      <Name>HeuristicLab.Data-3.3</Name>
+    </ProjectReference>
     <ProjectReference Include="..\HeuristicLab.Optimization\3.3\HeuristicLab.Optimization-3.3.csproj">
       <Project>{14AB8D24-25BC-400C-A846-4627AA945192}</Project>
@@ -106 +115 @@
       <Project>{edcc8202-4463-4122-b01f-21b664343dfb}</Project>
       <Name>HeuristicLab.Problems.DynamicalSystemsModelling-3.3</Name>
+    </ProjectReference>
+    <ProjectReference Include="..\HeuristicLab.Problems.Instances.DataAnalysis\3.3\HeuristicLab.Problems.Instances.DataAnalysis-3.3.csproj">
+      <Project>{94c7714e-29d4-4d6d-b213-2c18d627ab75}</Project>
+      <Name>HeuristicLab.Problems.Instances.DataAnalysis-3.3</Name>
     </ProjectReference>
   </ItemGroup>

branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs

@@ -39 +39 @@
 
 namespace HeuristicLab.Problems.DynamicalSystemsModelling {
-  // Eine weitere Möglichkeit ist spline-smoothing der Daten (über Zeit) um damit für jede Zielvariable
-  // einen bereinigten (Rauschen) Wert und die Ableitung dy/dt für alle Beobachtungen zu bekommen
-  // danach kann man die Regression direkt für dy/dt anwenden (mit bereinigten y als inputs)
   [Item("Dynamical Systems Modelling Problem", "TODO")]
   [Creatable(CreatableAttribute.Categories.GeneticProgrammingProblems, Priority = 900)]
@@ -179 +176 @@
       Parameters.Add(new FixedValueParameter<BoolValue>(OptimizeParametersForEpisodesParameterName, "Flag to select if parameters should be optimized globally or for each episode individually.", new BoolValue(false)));
 
-      var solversStr = new string[] { "HeuristicLab", "CVODES" };
+      var solversStr = new string[] { "HeuristicLab" /* , "CVODES" */};
       var solvers = new ItemSet<StringValue>(
         solversStr.Select(s => new StringValue(s).AsReadOnly())
@@ -188 +185 @@
       InitAllParameters();
 
-      // TODO: do not clear selection of target variables when the input variables are changed (keep selected target variables)
       // TODO: UI hangs when selecting / deselecting input variables because the encoding is updated on each item
       // TODO: use training range as default training episode
       // TODO: write back optimized parameters to solution?
       // TODO: optimization of starting values for latent variables in CVODES solver
+      // TODO: allow to specify the name for the time variable in the dataset and allow variable step-sizes
+      // TODO: check grammars (input variables) after cloning
 
     }
@@ -201 +199 @@
                                                                                                       // retreive optimized parameters from nodes?
 
-      var problemData = Standardize(ProblemData);
+      var problemData = ProblemData;
       var targetVars = TargetVariables.CheckedItems.OrderBy(i => i.Index).Select(i => i.Value.Value).ToArray();
       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
       if (OptimizeParametersForEpisodes) {
+        throw new NotImplementedException();
         int eIdx = 0;
         double totalNMSE = 0.0;
@@ -225 +224 @@
         return nmse;
       }
-    }
-
-    private IRegressionProblemData Standardize(IRegressionProblemData problemData) {
-      // var standardizedDataset = new Dataset(
-      //   problemData.Dataset.DoubleVariables,
-      //   problemData.Dataset.DoubleVariables.Select(v => Standardize(problemData.Dataset.GetReadOnlyDoubleValues(v)).ToList()));
-      // return new RegressionProblemData(standardizedDataset, problemData.AllowedInputVariables, problemData.TargetVariable);
-      return new RegressionProblemData(problemData.Dataset, problemData.AllowedInputVariables, problemData.TargetVariable);
     }
 
@@ -247 +238 @@
       out double[] optTheta,
       out double nmse) {
-      var rows = episodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start)).ToArray();
-      var targetValues = new double[rows.Length, targetVars.Length];
-
-
-      // collect values of all target variables                                            
-      var colIdx = 0;
-      foreach (var targetVar in targetVars) {
-        int rowIdx = 0;
-        foreach (var value in problemData.Dataset.GetDoubleValues(targetVar, rows)) {
-          targetValues[rowIdx, colIdx] = value;
-          rowIdx++;
-        }
-        colIdx++;
-      }
-
+
+
+      // optimize parameters by fitting f(x,y) to calculated differences dy/dt(t)
+      nmse = PreTuneParameters(trees, problemData, targetVars, latentVariables, random, episodes, maxParameterOptIterations, out optTheta);
+
+      // optimize parameters using integration of f(x,y) to calculate y(t)
+      nmse = OptimizeParameters(trees, problemData, targetVars, latentVariables, episodes, maxParameterOptIterations, optTheta, numericIntegrationSteps, odeSolver, out optTheta);
+
+      if (double.IsNaN(nmse) || double.IsInfinity(nmse)) nmse = 100 * trees.Length * episodes.Sum(ep => ep.Size);
+    }
+
+    private static double PreTuneParameters(
+      ISymbolicExpressionTree[] trees,
+      IRegressionProblemData problemData,
+      string[] targetVars,
+      string[] latentVariables,
+      IRandom random,
+      IEnumerable<IntRange> episodes,
+      int maxParameterOptIterations,
+      out double[] optTheta) {
+      var thetas = new List<double>();
+      double nmse = 0.0;
       // NOTE: the order of values in parameter matches prefix order of constant nodes in trees
-      var paramNodes = new List<ISymbolicExpressionTreeNode>();
-      foreach (var t in trees) {
-        foreach (var n in t.IterateNodesPrefix()) {
-          if (IsConstantNode(n))
-            paramNodes.Add(n);
-        }
-      }
-      // init params randomly
-      // theta contains parameter values for trees and then the initial values for latent variables (a separate vector for each episode)
-      // inital values for latent variables are also optimized
-      var theta = new double[paramNodes.Count + latentVariables.Length * episodes.Count()];
-      for (int i = 0; i < theta.Length; i++)
-        theta[i] = random.NextDouble() * 2.0e-1 - 1.0e-1;
-
-      optTheta = new double[0];
-      if (theta.Length > 0) {
-        alglib.minlmstate state;
-        alglib.minlmreport report;
-        alglib.minlmcreatevj(targetValues.Length, theta, out state);
-        alglib.minlmsetcond(state, 0.0, 0.0, 0.0, maxParameterOptIterations);
-        alglib.minlmsetgradientcheck(state, 1.0e-3);
-        var myState = new OptimizationData(trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver);
-        alglib.minlmoptimize(state, EvaluateObjectiveVector, EvaluateObjectiveVectorAndJacobian, null, myState);
-
-        alglib.minlmresults(state, out optTheta, out report);
-
-        /*************************************************************************
-         Levenberg-Marquardt algorithm results
-         
-         INPUT PARAMETERS:
-             State   -   algorithm state
-         
-         OUTPUT PARAMETERS:
-             X       -   array[0..N-1], solution
-             Rep     -   optimization  report;  includes  termination   codes   and
-                         additional information. Termination codes are listed below,
-                         see comments for this structure for more info.
-                         Termination code is stored in rep.terminationtype field:
-                         * -8    optimizer detected NAN/INF values either in the
-                                 function itself, or in its Jacobian
-                         * -7    derivative correctness check failed;
-                                 see rep.funcidx, rep.varidx for
-                                 more information.
-                         * -3    constraints are inconsistent
-                         *  2    relative step is no more than EpsX.
-                         *  5    MaxIts steps was taken
-                         *  7    stopping conditions are too stringent,
-                                 further improvement is impossible
-                         *  8    terminated by user who called minlmrequesttermination().
-                                 X contains point which was "current accepted" when
-                                 termination request was submitted.
-         
-           -- ALGLIB --
-              Copyright 10.03.2009 by Bochkanov Sergey
-         *************************************************************************/
-        if (report.terminationtype < 0) { nmse = 10.0; return; }
-
-        nmse = state.f; //TODO check
-
-        // var myState = new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver };
-        // EvaluateObjectiveVector(optTheta, ref nmse, grad,myState);
-        if (double.IsNaN(nmse) || double.IsInfinity(nmse)) { nmse = 10.0; return; } // return a large value (TODO: be consistent by using NMSE)
-      } else {
-        // no parameters
-        nmse = targetValues.Length;
-      }
-
-    }
-
-    // private static IEnumerable<double> Standardize(IEnumerable<double> xs) {
-    //   var m = xs.Average();
-    //   var s = xs.StandardDeviationPop();
-    //   return xs.Select(xi => (xi - m) / s);
-    // }
-
-    alglib.ndimensional_fvec fvec;
-    alglib.ndimensional_jac jac;
-
+      for (int treeIdx = 0; treeIdx < trees.Length; treeIdx++) {
+        var t = trees[treeIdx];
+
+        // TODO: calculation of differences for multiple episodes
+        var rowsForDataExtraction = episodes.SelectMany(e => Enumerable.Range(e.Start, e.End + 1 - e.Start)).ToArray();
+        var targetValues = problemData.Dataset.GetDoubleValues(targetVars[treeIdx], rowsForDataExtraction).ToArray();
+        var targetValuesDiff = targetValues.Skip(1).Zip(targetValues, (t1, t0) => t1 - t0).ToArray();// TODO: smoothing or multi-pole
+
+        var myState = new OptimizationData(new[] { t }, targetVars, problemData, new[] { targetValuesDiff }, episodes.ToArray(), -99, latentVariables, string.Empty); // TODO
+        var paramCount = myState.nodeValueLookup.ParameterCount;
+
+        // init params randomly
+        // theta contains parameter values for trees and then the initial values for latent variables (a separate vector for each episode)
+        // inital values for latent variables are also optimized
+        var theta = new double[paramCount + latentVariables.Length * episodes.Count()];
+        for (int i = 0; i < theta.Length; i++)
+          theta[i] = random.NextDouble() * 2.0e-1 - 1.0e-1;
+
+        optTheta = new double[0];
+        if (theta.Length > 0) {
+          alglib.minlmstate state;
+          alglib.minlmreport report;
+          alglib.minlmcreatevj(targetValuesDiff.Length, theta, out state);
+          alglib.minlmsetcond(state, 0.0, 0.0, 0.0, maxParameterOptIterations);
+          // alglib.minlmsetgradientcheck(state, 1.0e-3);
+          alglib.minlmoptimize(state, EvaluateObjectiveVector, EvaluateObjectiveVectorAndJacobian, null, myState);
+
+          alglib.minlmresults(state, out optTheta, out report);
+
+
+          if (report.terminationtype < 0) { throw new InvalidOperationException("there was a problem in the optimizer"); }
+
+          thetas.AddRange(optTheta);
+        }
+        nmse += EvaluateMSE(optTheta, myState);
+      } // foreach tree
+      optTheta = thetas.ToArray();
+
+      var maxNmse = 100 * trees.Length * episodes.Sum(ep => ep.Size);
+      if (double.IsNaN(nmse) || double.IsInfinity(nmse) || nmse > maxNmse) nmse = maxNmse;
+      return nmse;
+    }
+
+
+    // similar to above but this time we integrate and optimize all parameters for all targets concurrently
+    private static double OptimizeParameters(ISymbolicExpressionTree[] trees, IRegressionProblemData problemData, string[] targetVars, string[] latentVariables,
+      IEnumerable<IntRange> episodes, int maxParameterOptIterations, double[] initialTheta, int numericIntegrationSteps, string odeSolver, out double[] optTheta) {
+      var rowsForDataExtraction = episodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start)).ToArray();
+      var targetValues = new double[trees.Length][];
+      for (int treeIdx = 0; treeIdx < trees.Length; treeIdx++) {
+        var t = trees[treeIdx];
+
+        // TODO: calculation of differences for multiple episodes
+        targetValues[treeIdx] = problemData.Dataset.GetDoubleValues(targetVars[treeIdx], rowsForDataExtraction).ToArray();
+      }
+
+      var myState = new OptimizationData(trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver);
+      optTheta = initialTheta;
+
+      if (initialTheta.Length > 0) {
+
+        try {
+          alglib.minlmstate state;
+          alglib.minlmreport report;
+          alglib.minlmcreatevj(rowsForDataExtraction.Length * trees.Length, initialTheta, out state);
+          alglib.minlmsetcond(state, 0.0, 0.0, 0.0, maxParameterOptIterations);
+          // alglib.minlmsetgradientcheck(state, 1.0e-3);
+          alglib.minlmoptimize(state, IntegrateAndEvaluateObjectiveVector, IntegrateAndEvaluateObjectiveVectorAndJacobian, null, myState);
+
+          alglib.minlmresults(state, out optTheta, out report);
+
+          if (report.terminationtype < 0) {
+            // there was a problem: reset theta and evaluate for inital values
+            optTheta = initialTheta;
+          }
+        } catch (alglib.alglibexception) {
+          optTheta = initialTheta;
+        }
+      }
+      var nmse = EvaluateIntegratedMSE(optTheta, myState);
+      var maxNmse = 100 * targetValues.Length * rowsForDataExtraction.Length;
+      if (double.IsNaN(nmse) || double.IsInfinity(nmse) || nmse > maxNmse) nmse = maxNmse;
+      return nmse;
+    }
+
+
+    // helper
+    public static double EvaluateMSE(double[] x, OptimizationData optimizationData) {
+      var fi = new double[optimizationData.rows.Count()];
+      EvaluateObjectiveVector(x, fi, optimizationData);
+      return fi.Sum(fii => fii * fii) / fi.Length;
+    }
     public static void EvaluateObjectiveVector(double[] x, double[] fi, object optimizationData) { EvaluateObjectiveVector(x, fi, (OptimizationData)optimizationData); } // for alglib
     public static void EvaluateObjectiveVector(double[] x, double[] fi, OptimizationData optimizationData) {
-      EvaluateObjectiveVectorAndJacobian(x, fi, null, optimizationData);
+      var rows = optimizationData.rows.ToArray();
+      var problemData = optimizationData.problemData;
+      var nodeValueLookup = optimizationData.nodeValueLookup;
+      var ds = problemData.Dataset;
+      int outputIdx = 0;
+
+      nodeValueLookup.UpdateParamValues(x);
+
+      for (int trainIdx = 0; trainIdx < rows.Length; trainIdx++) {
+        // update variable values
+        foreach (var variable in problemData.AllowedInputVariables.Concat(optimizationData.targetVariables)) {
+          nodeValueLookup.SetVariableValue(variable, ds.GetDoubleValue(variable, rows[trainIdx]));
+        }
+        // interpret all trees
+        for (int treeIdx = 0; treeIdx < optimizationData.trees.Length; treeIdx++) {
+          var tree = optimizationData.trees[treeIdx];
+          var pred = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValueLookup);
+          var y = optimizationData.targetValues[treeIdx][trainIdx];
+          fi[outputIdx++] = y - pred;
+        }
+      }
     }
 
     public static void EvaluateObjectiveVectorAndJacobian(double[] x, double[] fi, double[,] jac, object optimizationData) { EvaluateObjectiveVectorAndJacobian(x, fi, jac, (OptimizationData)optimizationData); } // for alglib
     public static void EvaluateObjectiveVectorAndJacobian(double[] x, double[] fi, double[,] jac, OptimizationData optimizationData) {
-
-
-      Tuple<double, Vector>[][] predicted = null; // one array of predictions for each episode
-      predicted = Integrate(optimizationData, x).ToArray();
-
-      // clear all result data structures
-      for (int j = 0; j < fi.Length; j++) {
-        fi[j] = 10.0;
-        if (jac != null) Array.Clear(jac, 0, jac.Length);
-      }
-
-      if (predicted.Length != optimizationData.targetValues.GetLength(0)) {
-        return;
-      }
-
-      // for normalized MSE = 1/variance(t) * MSE(t, pred)
-      // TODO: Perf. (by standardization of target variables before evaluation of all trees)      
-      // var invVar = Enumerable.Range(0, targetVariables.Length)
-      //   .Select(c => Enumerable.Range(0, targetValues.GetLength(0)).Select(row => targetValues[row, c])) // column vectors
-      //   .Select(vec => vec.StandardDeviation()) // TODO: variance of stddev
-      //   .Select(v => 1.0 / v)
-      //   .ToArray();
-
-      // double[] invVar = Enumerable.Repeat(1.0, targetVariables.Length).ToArray();
-
-
-      // objective function is NMSE
-      int n = predicted.Length;
-      double invN = 1.0 / n;
-      int i = 0;
-      int r = 0;
-      foreach (var y_pred in predicted) {
-        // y_pred contains the predicted values for target variables first and then predicted values for latent variables
-        for (int c = 0; c < optimizationData.targetVariables.Length; c++) {
-
-          var y_pred_f = y_pred[c].Item1;
-          var y = optimizationData.targetValues[r, c];
-
-          fi[i] = (y - y_pred_f);
-          var g = y_pred[c].Item2;
-          if (jac != null && g != Vector.Zero) for (int j = 0; j < g.Length; j++) jac[i, j] = -g[j];
-          i++; // we put the errors for each target variable after each other 
-        }
-        r++;
+      // extract variable values from dataset
+      var variableValues = new Dictionary<string, Tuple<double, Vector>>();
+      var problemData = optimizationData.problemData;
+      var ds = problemData.Dataset;
+      var rows = optimizationData.episodes.SelectMany(ep => Enumerable.Range(ep.Start, ep.Size)).ToArray();
+
+      var nodeValueLookup = optimizationData.nodeValueLookup;
+      nodeValueLookup.UpdateParamValues(x);
+
+      // TODO add integration of latent variables
+      int termIdx = 0;
+
+      for (int trainIdx = 0; trainIdx < rows.Length; trainIdx++) {
+        // update variable values
+        foreach (var variable in problemData.AllowedInputVariables.Concat(optimizationData.targetVariables)) {
+          nodeValueLookup.SetVariableValue(variable, ds.GetDoubleValue(variable, rows[trainIdx]));
+        }
+
+        var calculatedVariables = optimizationData.targetVariables;
+
+        var trees = optimizationData.trees;
+        for (int i = 0; i < trees.Length; i++) {
+          var tree = trees[i];
+          var targetVarName = calculatedVariables[i];
+
+          double f; Vector g;
+          InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValueLookup, out f, out g);
+
+          var y = optimizationData.targetValues[i][trainIdx];
+          fi[termIdx] = y - f;
+          if (jac != null && g != Vector.Zero) for (int j = 0; j < g.Length; j++) jac[termIdx, j] = -g[j];
+
+          termIdx++;
+        }
+      }
+
+    }
+
+    // helper
+    public static double EvaluateIntegratedMSE(double[] x, OptimizationData optimizationData) {
+      var fi = new double[optimizationData.rows.Count() * optimizationData.targetVariables.Length];
+      IntegrateAndEvaluateObjectiveVector(x, fi, optimizationData);
+      return fi.Sum(fii => fii * fii) / fi.Length;
+    }
+    public static void IntegrateAndEvaluateObjectiveVector(double[] x, double[] fi, object optimizationData) { IntegrateAndEvaluateObjectiveVector(x, fi, (OptimizationData)optimizationData); } // for alglib
+    public static void IntegrateAndEvaluateObjectiveVector(double[] x, double[] fi, OptimizationData optimizationData) {
+      IntegrateAndEvaluateObjectiveVectorAndJacobian(x, fi, null, optimizationData);
+    }
+
+    public static void IntegrateAndEvaluateObjectiveVectorAndJacobian(double[] x, double[] fi, double[,] jac, object optimizationData) { IntegrateAndEvaluateObjectiveVectorAndJacobian(x, fi, jac, (OptimizationData)optimizationData); } // for alglib
+    public static void IntegrateAndEvaluateObjectiveVectorAndJacobian(double[] x, double[] fi, double[,] jac, OptimizationData optimizationData) {
+      var rows = optimizationData.rows.ToArray();
+      var problemData = optimizationData.problemData;
+      var nodeValueLookup = optimizationData.nodeValueLookup;
+      var ds = problemData.Dataset;
+      int outputIdx = 0;
+
+      var prediction = Integrate(optimizationData, x).ToArray();
+      var trees = optimizationData.trees;
+
+      for (int trainIdx = 0; trainIdx < rows.Length; trainIdx++) {
+        var pred_i = prediction[Math.Min(trainIdx, prediction.Length - 1)]; // if we stopped earlier in the integration then just use the last generated value
+
+        for (int i = 0; i < trees.Length; i++) {
+          var tree = trees[i];
+          var y = optimizationData.targetValues[i][trainIdx];
+          var f = pred_i[i].Item1;
+          var g = pred_i[i].Item2;
+          fi[outputIdx] = y - f;
+          if (jac != null && g != Vector.Zero) for (int j = 0; j < g.Length; j++) jac[outputIdx, j] = -g[j];
+
+          outputIdx++;
+        }
       }
     }
@@ -422 +483 @@
       results["SNMSE"].Value = new DoubleValue(bestIndividualAndQuality.Item2);
 
-      var problemData = Standardize(ProblemData);
+      var problemData = ProblemData;
       var targetVars = TargetVariables.CheckedItems.OrderBy(i => i.Index).Select(i => i.Value.Value).ToArray();
       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
@@ -519 +580 @@
             var res = (y - y_pred_f);
             var ressq = res * res;
-            f_i += ressq * invN /* * Math.Exp(-0.2 * r) */;
-            g_i = g_i - 2.0 * res * y_pred[colIdx].Item2 * invN /* * Math.Exp(-0.2 * r)*/;
+            f_i += ressq * invN;
+            g_i = g_i - 2.0 * res * y_pred[colIdx].Item2 * invN;
           }
           seRow.Values.Add(f_i);
@@ -603 +664 @@
             ref nextParIdx));
 
-          // var shownTree = (SymbolicExpressionTree)tree.Clone();
-          // var constantsNodeOrig = tree.IterateNodesPrefix().Where(IsConstantNode);
-          // var constantsNodeShown = shownTree.IterateNodesPrefix().Where(IsConstantNode);
-          // 
-          // foreach (var n in constantsNodeOrig.Zip(constantsNodeShown, (original, shown) => new { original, shown })) {
-          //   double constantsVal = optTheta[nodeIdx[n.original]];
-          // 
-          //   ConstantTreeNode replacementNode = new ConstantTreeNode(new Constant()) { Value = constantsVal };
-          // 
-          //   var parentNode = n.shown.Parent;
-          //   int replacementIndex = parentNode.IndexOfSubtree(n.shown);
-          //   parentNode.RemoveSubtree(replacementIndex);
-          //   parentNode.InsertSubtree(replacementIndex, replacementNode);
-          // }
 
           var origTreeVar = new HeuristicLab.Core.Variable(varName + "(original)");
@@ -655 +702 @@
       var odeSolver = optimizationData.odeSolver;
       var numericIntegrationSteps = optimizationData.numericIntegrationSteps;
+      var calculatedVariables = targetVariables.Concat(latentVariables).ToArray(); // TODO: must conincide with the order of trees in the encoding
+
+
+      var nodeValues = optimizationData.nodeValueLookup;
+      nodeValues.UpdateParamValues(parameterValues);
 
 
@@ -660 +712 @@
       var episodeIdx = 0;
       foreach (var episode in optimizationData.episodes) {
-        var rows = Enumerable.Range(episode.Start, episode.End - episode.Start);
-
-        // integrate forward starting with known values for the target in t0
-
-        var variableValues = new Dictionary<string, Tuple<double, Vector>>();
+        var rows = Enumerable.Range(episode.Start, episode.End - episode.Start).ToArray();
+
         var t0 = rows.First();
-        foreach (var varName in inputVariables) {
-          variableValues.Add(varName, Tuple.Create(dataset.GetDoubleValue(varName, t0), Vector.Zero));
-        }
-        foreach (var varName in targetVariables) {
-          variableValues.Add(varName, Tuple.Create(dataset.GetDoubleValue(varName, t0), Vector.Zero));
-        }
-        // add value entries for latent variables which are also integrated
-        // initial values are at the end of the parameter vector
-        // separete initial values for each episode
-        var initialValueIdx = parameterValues.Length - episodes.Count() * latentVariables.Length + episodeIdx * latentVariables.Length;
-        foreach (var latentVar in latentVariables) {
-          var arr = new double[parameterValues.Length]; // backing array
-          arr[initialValueIdx] = 1.0;
-          var g = new Vector(arr);
-          variableValues.Add(latentVar,
-            Tuple.Create(parameterValues[initialValueIdx], g)); // we don't have observations for latent variables therefore we optimize the initial value for each episode 
-          initialValueIdx++;
-        }
-
-        var calculatedVariables = targetVariables.Concat(latentVariables).ToArray(); // TODO: must conincide with the order of trees in the encoding
+
+        // initialize values for inputs and targets from dataset
+        foreach (var varName in inputVariables.Concat(targetVariables)) {
+          nodeValues.SetVariableValue(varName, dataset.GetDoubleValue(varName, t0), Vector.Zero);
+        }
+
+        { // CODE BELOW MUST BE TESTED
+          if (latentVariables.Length > 0) throw new NotImplementedException();
+
+          // add value entries for latent variables which are also integrated
+          // initial values are at the end of the parameter vector
+          // separate initial values for each episode
+          var initialValueIdx = parameterValues.Length - episodes.Count() * latentVariables.Length + episodeIdx * latentVariables.Length;
+          foreach (var latentVar in latentVariables) {
+            var arr = new double[parameterValues.Length]; // backing array
+            arr[initialValueIdx] = 1.0;
+            var g = new Vector(arr);
+            nodeValues.SetVariableValue(latentVar, parameterValues[initialValueIdx], g); // we don't have observations for latent variables therefore we optimize the initial value for each episode 
+            initialValueIdx++;
+          }
+        }
 
         // return first value as stored in the dataset
         yield return calculatedVariables
-          .Select(calcVarName => variableValues[calcVarName])
+          .Select(calcVarName => nodeValues.GetVariableValue(calcVarName))
           .ToArray();
 
-        var prevT = rows.First(); // TODO: here we should use a variable for t if it is available. Right now we assume equidistant measurements.
+        var prevT = t0; // TODO: here we should use a variable for t if it is available. Right now we assume equidistant measurements.
         foreach (var t in rows.Skip(1)) {
           if (odeSolver == "HeuristicLab")
-            IntegrateHL(trees, calculatedVariables, variableValues, parameterValues, numericIntegrationSteps);
+            IntegrateHL(trees, calculatedVariables, nodeValues, numericIntegrationSteps); // integrator updates nodeValues
           else if (odeSolver == "CVODES")
             throw new NotImplementedException();
@@ -702 +752 @@
           prevT = t;
 
-          // This check doesn't work with the HeuristicLab integrator if there are input variables 
-          //if (variableValues.Count == targetVariables.Length) {
-          // only return the target variables for calculation of errors
+
           var res = calculatedVariables
-            .Select(targetVar => variableValues[targetVar])
+            .Select(targetVar => nodeValues.GetVariableValue(targetVar))
             .ToArray();
+
+          // stop early if there are undefined values
           if (res.Any(ri => double.IsNaN(ri.Item1) || double.IsInfinity(ri.Item1))) yield break;
           yield return res;
-          //} else {
-          //  yield break; // stop early on problems
-          //}
-
-
-          // update for next time step
+
+          // update for next time step (only the inputs)
           foreach (var varName in inputVariables) {
-            variableValues[varName] = Tuple.Create(dataset.GetDoubleValue(varName, t), Vector.Zero);
+            nodeValues.SetVariableValue(varName, dataset.GetDoubleValue(varName, t), Vector.Zero);
           }
         }
@@ -1026 +1072 @@
       ISymbolicExpressionTree[] trees,
       string[] calculatedVariables, // names of integrated variables
-      Dictionary<string, Tuple<double, Vector>> variableValues, //y (intput and output) input: y(t0), output: y(t0+1)
-      double[] parameterValues,
+                                    // Dictionary<string, Tuple<double, Vector>> variableValues, //y (intput and output) input: y(t0), output: y(t0+1)
+      NodeValueLookup nodeValues,
+      // double[] parameterValues,
       int numericIntegrationSteps) {
 
-      var nodeValues = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
-
-      // the nodeValues for parameters are constant over time
-      // TODO: this needs to be done only once for each iteration in gradient descent (whenever parameter values change)
-      // NOTE: the order must be the same as above (prefix order for nodes)
-      int pIdx = 0;
-      foreach (var tree in trees) {
-        foreach (var node in tree.Root.IterateNodesPrefix()) {
-          if (IsConstantNode(node)) {
-            var gArr = new double[parameterValues.Length]; // backing array
-            gArr[pIdx] = 1.0;
-            var g = new Vector(gArr);
-            nodeValues.Add(node, new Tuple<double, Vector>(parameterValues[pIdx], g));
-            pIdx++;
-          } else if (node.SubtreeCount == 0) {
-            // for (latent) variables get the values from variableValues
-            var varName = node.Symbol.Name;
-            nodeValues.Add(node, variableValues[varName]);
-          }
-        }
-      }
 
       double[] deltaF = new double[calculatedVariables.Length];
@@ -1057 +1083 @@
       double h = 1.0 / numericIntegrationSteps;
       for (int step = 0; step < numericIntegrationSteps; step++) {
-        //var deltaValues = new Dictionary<string, Tuple<double, Vector>>();
+
+        // evaluate all trees
         for (int i = 0; i < trees.Length; i++) {
           var tree = trees[i];
-          var targetVarName = calculatedVariables[i];
 
           // Root.GetSubtree(0).GetSubtree(0) skips programRoot and startSymbol
@@ -1072 +1098 @@
         for (int i = 0; i < trees.Length; i++) {
           var varName = calculatedVariables[i];
-          var oldVal = variableValues[varName];
-          var newVal = Tuple.Create(
-            oldVal.Item1 + h * deltaF[i],
-            oldVal.Item2 + deltaG[i].Scale(h)
-          );
-          variableValues[varName] = newVal;
-        }
-
-        // TODO perf
-        foreach (var node in nodeValues.Keys.ToArray()) {
-          if (node.SubtreeCount == 0 && !IsConstantNode(node)) {
-            // update values for (latent) variables
-            var varName = node.Symbol.Name;
-            nodeValues[node] = variableValues[varName];
-          }
-        }
+          var oldVal = nodeValues.GetVariableValue(varName);
+          nodeValues.SetVariableValue(varName, oldVal.Item1 + h * deltaF[i], oldVal.Item2 + deltaG[i].Scale(h));
+        }
+      }
+    }
+
+    private static double InterpretRec(ISymbolicExpressionTreeNode node, NodeValueLookup nodeValues) {
+      switch (node.Symbol.Name) {
+        case "+": {
+            var f = InterpretRec(node.GetSubtree(0), nodeValues);
+            var g = InterpretRec(node.GetSubtree(1), nodeValues);
+            return f + g;
+          }
+        case "*": {
+            var f = InterpretRec(node.GetSubtree(0), nodeValues);
+            var g = InterpretRec(node.GetSubtree(1), nodeValues);
+            return f * g;
+          }
+
+        case "-": {
+            if (node.SubtreeCount == 1) {
+              var f = InterpretRec(node.GetSubtree(0), nodeValues);
+              return -f;
+            } else {
+              var f = InterpretRec(node.GetSubtree(0), nodeValues);
+              var g = InterpretRec(node.GetSubtree(1), nodeValues);
+
+              return f - g;
+            }
+          }
+        case "%": {
+            var f = InterpretRec(node.GetSubtree(0), nodeValues);
+            var g = InterpretRec(node.GetSubtree(1), nodeValues);
+
+            // protected division
+            if (g.IsAlmost(0.0)) {
+              return 0;
+            } else {
+              return f / g;
+            }
+          }
+        case "sin": {
+            var f = InterpretRec(node.GetSubtree(0), nodeValues);
+            return Math.Sin(f);
+          }
+        case "cos": {
+            var f = InterpretRec(node.GetSubtree(0), nodeValues);
+            return Math.Cos(f);
+          }
+        case "sqr": {
+            var f = InterpretRec(node.GetSubtree(0), nodeValues);
+            return f * f;
+          }
+        default: {
+            return nodeValues.NodeValue(node);
+          }
       }
     }
@@ -1093 +1159 @@
     private static void InterpretRec(
       ISymbolicExpressionTreeNode node,
-      Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>> nodeValues,      // contains value and gradient vector for a node (variables and constants only)
+       NodeValueLookup nodeValues,      // contains value and gradient vector for a node (variables and constants only)
       out double z,
       out Vector dz
@@ -1162 +1228 @@
           }
         default: {
-            var t = nodeValues[node];
-            z = t.Item1;
-            dz = Vector.CreateNew(t.Item2);
+            z = nodeValues.NodeValue(node);
+            dz = Vector.CreateNew(nodeValues.NodeGradient(node));
             break;
           }
@@ -1265 +1330 @@
       return n.Symbol.Name[0] == 'θ';
     }
+    private static double GetConstantValue(ISymbolicExpressionTreeNode n) {
+      return 0.0; // TODO: needs to be updated when we write back values to the tree
+    }
     private static bool IsLatentVariableNode(ISymbolicExpressionTreeNode n) {
       return n.Symbol.Name[0] == 'λ';
@@ -1270 +1338 @@
     private static bool IsVariableNode(ISymbolicExpressionTreeNode n) {
       return (n.SubtreeCount == 0) && !IsConstantNode(n) && !IsLatentVariableNode(n);
+    }
+    private static string GetVariableName(ISymbolicExpressionTreeNode n) {
+      return n.Symbol.Name;
     }
 
@@ -1410 +1481 @@
     #endregion
 
+
     #region Import & Export
     public void Load(IRegressionProblemData data) {
@@ -1420 +1492 @@
       return ProblemData;
     }
+    #endregion
+
+
+    // TODO: for integration we only need a part of the data that we need for optimization
 
     public class OptimizationData {
@@ -1425 +1501 @@
       public readonly string[] targetVariables;
       public readonly IRegressionProblemData problemData;
-      public readonly double[,] targetValues;
+      public readonly double[][] targetValues;
       public readonly IntRange[] episodes;
       public readonly int numericIntegrationSteps;
       public readonly string[] latentVariables;
       public readonly string odeSolver;
-
-      public OptimizationData(ISymbolicExpressionTree[] trees, string[] targetVars, IRegressionProblemData problemData, double[,] targetValues, IntRange[] episodes, int numericIntegrationSteps, string[] latentVariables, string odeSolver) {
+      public readonly NodeValueLookup nodeValueLookup;
+      public IEnumerable<int> rows => episodes.SelectMany(ep => Enumerable.Range(ep.Start, ep.Size));
+
+      public OptimizationData(ISymbolicExpressionTree[] trees, string[] targetVars, IRegressionProblemData problemData,
+        double[][] targetValues,
+        IntRange[] episodes,
+        int numericIntegrationSteps, string[] latentVariables, string odeSolver) {
+        if (episodes.Length > 1) throw new NotSupportedException("multiple episodes are not yet implemented");
         this.trees = trees;
         this.targetVariables = targetVars;
@@ -1440 +1522 @@
         this.latentVariables = latentVariables;
         this.odeSolver = odeSolver;
-      }
-    }
-    #endregion
-
+        this.nodeValueLookup = new NodeValueLookup(trees);
+      }
+    }
+
+    public class NodeValueLookup {
+      private readonly Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>> node2val = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
+      private readonly Dictionary<string, List<ISymbolicExpressionTreeNode>> name2nodes = new Dictionary<string, List<ISymbolicExpressionTreeNode>>();
+      private readonly Dictionary<int, ISymbolicExpressionTreeNode> paramIdx2node = new Dictionary<int, ISymbolicExpressionTreeNode>();
+
+      public double NodeValue(ISymbolicExpressionTreeNode node) => node2val[node].Item1;
+      public Vector NodeGradient(ISymbolicExpressionTreeNode node) => node2val[node].Item2;
+
+      public NodeValueLookup(ISymbolicExpressionTree[] trees) {
+        int paramIdx = 0;
+
+        var constantNodes = trees.SelectMany(t => t.IterateNodesPrefix().Where(IsConstantNode)).ToArray();
+        foreach (var n in constantNodes) {
+          paramIdx2node[paramIdx] = n;
+          SetParamValue(paramIdx, GetConstantValue(n), Vector.CreateIndicator(length: constantNodes.Length, idx: paramIdx));
+          paramIdx++;
+        }
+
+        foreach (var tree in trees) {
+          foreach (var node in tree.IterateNodesPrefix().Where(IsVariableNode)) {
+            var varName = GetVariableName(node);
+            if (!name2nodes.TryGetValue(varName, out List<ISymbolicExpressionTreeNode> nodes)) {
+              nodes = new List<ISymbolicExpressionTreeNode>();
+              name2nodes.Add(varName, nodes);
+            }
+            nodes.Add(node);
+            SetVariableValue(varName, 0.0);
+          }
+        }
+      }
+
+      public int ParameterCount => paramIdx2node.Count;
+
+
+      public void SetParamValue(int paramIdx, double val, Vector dVal) {
+        node2val[paramIdx2node[paramIdx]] = Tuple.Create(val, dVal);
+      }
+
+      public void SetVariableValue(string variableName, double val) {
+        SetVariableValue(variableName, val, Vector.Zero);
+      }
+      public Tuple<double, Vector> GetVariableValue(string variableName) {
+        return node2val[name2nodes[variableName].First()];
+      }
+      public void SetVariableValue(string variableName, double val, Vector dVal) {
+        if (name2nodes.TryGetValue(variableName, out List<ISymbolicExpressionTreeNode> nodes)) {
+          nodes.ForEach(n => node2val[n] = Tuple.Create(val, dVal));
+        } else {
+          var fakeNode = new SimpleSymbol(variableName, 0).CreateTreeNode();
+          var newNodeList = new List<ISymbolicExpressionTreeNode>();
+          newNodeList.Add(fakeNode);
+          name2nodes.Add(variableName, newNodeList);
+          node2val[fakeNode] = Tuple.Create(val, dVal);
+        }
+      }
+
+      internal void UpdateParamValues(double[] x) {
+        Trace.Assert(x.Length == paramIdx2node.Count);
+        for (int paramIdx = 0; paramIdx < x.Length; paramIdx++) {
+          var n = paramIdx2node[paramIdx];
+          node2val[n] = Tuple.Create(x[paramIdx], node2val[n].Item2); // prevent allocation of new Vector
+        }
+      }
+    }
   }
 }

branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Vector.cs

@@ -1 +1 @@
 using System;
 using System.Diagnostics;
+using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
 
 namespace HeuristicLab.Problems.DynamicalSystemsModelling {
@@ -158 +159 @@
       return new Vector(arr);
     }
+
+    internal static Vector CreateIndicator(int length, int idx) {
+      var gArr = new double[length]; // backing array
+      gArr[idx] = 1.0;
+      return new Vector(gArr);
+    }
   }
 }