Changeset 16660 for branches/2925_AutoDiffForDynamicalModels

Timestamp:

03/07/19 17:00:39 (6 years ago)

Author:

gkronber

Message:

#2925: re-introduced partial support for latent variables

Location:

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

Files:

• Problem.cs (modified) (20 diffs)
• Solution.cs (modified) (3 diffs)
• SolutionView.cs (modified) (1 diff)

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

@@ -187 +187 @@
       InitAllParameters();
 
-      // 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
-
     }
 
@@ -235 +231 @@
       string odeSolver) {
 
-      // extract constants from tree
-      var constantNodes = trees.Select(t => t.IterateNodesPrefix().OfType<ConstantTreeNode>().ToArray()).ToArray();
+      // extract constants from trees (without trees for latent variables)
+      var targetVariableTrees = trees.Take(targetVars.Length).ToArray();
+      var latentVariableTrees = trees.Skip(targetVars.Length).ToArray();
+      var constantNodes = targetVariableTrees.Select(t => t.IterateNodesPrefix().OfType<ConstantTreeNode>().ToArray()).ToArray();
       var initialTheta = constantNodes.Select(nodes => nodes.Select(n => n.Value).ToArray()).ToArray();
 
@@ -242 +240 @@
       double nmse = PreTuneParameters(trees, problemData, targetVars, latentVariables, random, episodes, maxParameterOptIterations,
         initialTheta, out double[] pretunedParameters);
+
+      // extend parameter vector to include parameters for latent variable trees
+      pretunedParameters = pretunedParameters
+        .Concat(latentVariableTrees
+        .SelectMany(t => t.IterateNodesPrefix().OfType<ConstantTreeNode>().Select(n => n.Value)))
+        .ToArray();
 
       // optimize parameters using integration of f(x,y) to calculate y(t)
@@ -276 +280 @@
       var maxTreeNmse = 100 * episodes.Sum(ep => ep.Size);
 
+      var targetTrees = trees.Take(targetVars.Length).ToArray();
+      var latentTrees = trees.Take(latentVariables.Length).ToArray();
+
+      {
+        // first calculate values of latent variables by integration
+        var inputVariables = targetVars.Concat(latentTrees.SelectMany(t => t.IterateNodesPrefix().OfType<VariableTreeNode>().Select(n => n.VariableName))).Except(latentVariables).Distinct();
+        var myState = new OptimizationData(latentTrees, targetVars, inputVariables.ToArray(), problemData, null, episodes.ToArray(), 10, latentVariables, "HeuristicLab");
+
+        var fi = new double[myState.rows.Length * targetVars.Length];
+        var jac = new double[myState.rows.Length * targetVars.Length, myState.nodeValueLookup.ParameterCount];
+        var latentValues = new double[myState.rows.Length, latentVariables.Length];
+        Integrate(myState, fi, jac, latentValues);
+
+        // add integrated latent variables to dataset
+        var modifiedDataset = ((Dataset)problemData.Dataset).ToModifiable();
+        foreach (var variable in latentVariables) {
+          modifiedDataset.AddVariable(variable, Enumerable.Repeat(0.0, modifiedDataset.Rows).ToList()); // empty column
+        }
+        int predIdx = 0;
+        foreach (var ep in episodes) {
+          for (int r = ep.Start; r < ep.End; r++) {
+            for (int latVarIdx = 0; latVarIdx < latentVariables.Length; latVarIdx++) {
+              modifiedDataset.SetVariableValue(latentValues[predIdx, latVarIdx], latentVariables[latVarIdx], r);
+            }
+            predIdx++;
+          }
+        }
+
+        problemData = new RegressionProblemData(modifiedDataset, problemData.AllowedInputVariables, problemData.TargetVariable);
+      }
       // NOTE: the order of values in parameter matches prefix order of constant nodes in trees
-      for (int treeIdx = 0; treeIdx < trees.Length; treeIdx++) {
-        var t = trees[treeIdx];
+      for (int treeIdx = 0; treeIdx < targetTrees.Length; treeIdx++) {
+        var t = targetTrees[treeIdx];
 
         var targetValuesDiff = new List<double>();
@@ -290 +324 @@
         // data for input variables is assumed to be known
         // input variables in pretuning are all target variables and all variable names that occur in the tree
-        var inputVariables = targetVars.Concat(t.IterateNodesPrefix().OfType<VariableTreeNode>().Select(n => n.VariableName));
+        var inputVariables = targetVars.Concat(t.IterateNodesPrefix().OfType<VariableTreeNode>().Select(n => n.VariableName)).Distinct();
 
         var myState = new OptimizationData(new[] { t },
@@ -345 +379 @@
       IEnumerable<IntRange> episodes, int maxParameterOptIterations, double[] initialTheta, int numericIntegrationSteps, string odeSolver, out double[] optTheta) {
       var rowsForDataExtraction = episodes.SelectMany(e => Enumerable.Range(e.Start, e.Size)).ToArray();
-      var targetValues = new double[trees.Length][];
-      for (int treeIdx = 0; treeIdx < trees.Length; treeIdx++) {
+      var targetValues = new double[targetVars.Length][];
+      for (int treeIdx = 0; treeIdx < targetVars.Length; treeIdx++) {
         var t = trees[treeIdx];
 
@@ -354 +388 @@
       // data for input variables is assumed to be known
       // input variables are all variable names that occur in the trees except for target variables (we assume that trees have been generated correctly)
-      var inputVariables = trees.SelectMany(t => t.IterateNodesPrefix().OfType<VariableTreeNode>().Select(n => n.VariableName)).Except(targetVars);
+      var inputVariables = trees.SelectMany(t => t.IterateNodesPrefix().OfType<VariableTreeNode>().Select(n => n.VariableName))
+        .Except(targetVars)
+        .Except(latentVariables)
+        .Distinct();
 
       var myState = new OptimizationData(trees, targetVars, inputVariables.ToArray(), problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver);
@@ -495 +532 @@
       nodeValueLookup.UpdateParamValues(x);
 
-      Integrate(optimizationData, fi, jac);
+      Integrate(optimizationData, fi, jac, null);
       var trees = optimizationData.trees;
 
       // update result with error
       for (int trainIdx = 0; trainIdx < rows.Length; trainIdx++) {
-        for (int i = 0; i < trees.Length; i++) {
+        for (int i = 0; i < optimizationData.targetVariables.Length; i++) {
           var tree = trees[i];
           var y = optimizationData.targetValues[i][trainIdx];
@@ -583 +620 @@
         // data for input variables is assumed to be known
         // input variables are all variable names that occur in the trees except for target variables (we assume that trees have been generated correctly)
-        var inputVariables = trees.SelectMany(t => t.IterateNodesPrefix().OfType<VariableTreeNode>().Select(n => n.VariableName)).Except(targetVars);
+        var inputVariables = trees
+          .SelectMany(t => t.IterateNodesPrefix().OfType<VariableTreeNode>().Select(n => n.VariableName))
+          .Except(targetVars)
+          .Except(latentVariables)
+          .Distinct();
 
         var optimizationData = new OptimizationData(trees, targetVars, inputVariables.ToArray(), problemData, null, TrainingEpisodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver);
-        var trainingPrediction = Integrate(optimizationData).ToArray();
-
-
         var numParams = optimizationData.nodeValueLookup.ParameterCount;
+
+        var fi = new double[optimizationData.rows.Length * targetVars.Length];
+        var jac = new double[optimizationData.rows.Length * targetVars.Length, numParams];
+        var latentValues = new double[optimizationData.rows.Length, latentVariables.Length];
+        Integrate(optimizationData, fi, jac, latentValues);
+
+
         // for target values and latent variables
         var trainingRows = optimizationData.rows;
@@ -598 +643 @@
             var trainingDataTable = new DataTable(targetVar + " prediction (training)");
             var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, trainingRows));
-            var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, trainingPrediction.Select(arr => arr[colIdx].Item1).ToArray());
+            var idx = Enumerable.Range(0, trainingRows.Length).Select(i => i * targetVars.Length + colIdx);
+            var pred = idx.Select(i => fi[i]);
+            var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, pred.ToArray());
             trainingDataTable.Rows.Add(actualValuesRow);
             trainingDataTable.Rows.Add(predictedValuesRow);
 
             for (int paramIdx = 0; paramIdx < numParams; paramIdx++) {
-              var paramSensitivityRow = new DataRow($"∂{targetVar}/∂θ{paramIdx}", $"Sensitivities of parameter {paramIdx}", trainingPrediction.Select(arr => arr[colIdx].Item2[paramIdx]).ToArray());
+              var paramSensitivityRow = new DataRow($"∂{targetVar}/∂θ{paramIdx}", $"Sensitivities of parameter {paramIdx}", idx.Select(i => jac[i, paramIdx]).ToArray());
               paramSensitivityRow.VisualProperties.SecondYAxis = true;
               trainingDataTable.Rows.Add(paramSensitivityRow);
@@ -611 +658 @@
             var latentVar = latentVariables[colIdx - targetVars.Length];
             var trainingDataTable = new DataTable(latentVar + " prediction (training)");
-            var predictedValuesRow = new DataRow(latentVar + " pred.", "Predicted values for " + latentVar, trainingPrediction.Select(arr => arr[colIdx].Item1).ToArray());
+            var idx = Enumerable.Range(0, trainingRows.Length);
+            var pred = idx.Select(i => latentValues[i, colIdx - targetVars.Length]);
+            var predictedValuesRow = new DataRow(latentVar + " pred.", "Predicted values for " + latentVar, pred.ToArray());
             var emptyRow = new DataRow(latentVar);
             trainingDataTable.Rows.Add(emptyRow);
@@ -630 +679 @@
         int r = 0;
 
-        foreach (var y_pred in trainingPrediction) {
-          // calculate objective function gradient
-          double f_i = 0.0;
-          Vector g_i = Vector.CreateNew(new double[numParams]);
-          for (int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
-            var y_pred_f = y_pred[colIdx].Item1;
-            var y = targetValues[colIdx][r];
-
-            var res = (y - y_pred_f) * optimizationData.inverseStandardDeviation[colIdx];
-            var ressq = res * res;
-            f_i += ressq;
-            g_i.Add(y_pred[colIdx].Item2.Scale(-2.0 * res));
-          }
-          seRow.Values.Add(f_i);
-          for (int j = 0; j < g_i.Length; j++) gradientRows[j].Values.Add(g_i[j]);
-          r++;
-        }
-        results["Squared error and gradient"].Value = errorTable;
+        // foreach (var y_pred in trainingPrediction) {
+        //   // calculate objective function gradient
+        //   double f_i = 0.0;
+        //   Vector g_i = Vector.CreateNew(new double[numParams]);
+        //   for (int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
+        //     var y_pred_f = y_pred[colIdx].Item1;
+        //     var y = targetValues[colIdx][r];
+        // 
+        //     var res = (y - y_pred_f) * optimizationData.inverseStandardDeviation[colIdx];
+        //     var ressq = res * res;
+        //     f_i += ressq;
+        //     g_i.Add(y_pred[colIdx].Item2.Scale(-2.0 * res));
+        //   }
+        //   seRow.Values.Add(f_i);
+        //   for (int j = 0; j < g_i.Length; j++) gradientRows[j].Values.Add(g_i[j]);
+        //   r++;
+        // }
+        // results["Squared error and gradient"].Value = errorTable;
 
         // TODO: DRY for training and test
@@ -667 +716 @@
 
           } else {
-            var latentVar = latentVariables[colIdx - targetVars.Length];
-            var testDataTable = new DataTable(latentVar + " prediction (test)");
-            var predictedValuesRow = new DataRow(latentVar + " pred.", "Predicted values for " + latentVar, testPrediction.Select(arr => arr[colIdx].Item1).ToArray());
-            var emptyRow = new DataRow(latentVar);
-            testDataTable.Rows.Add(emptyRow);
-            testDataTable.Rows.Add(predictedValuesRow);
-            testList.Add(testDataTable);
+            // var latentVar = latentVariables[colIdx - targetVars.Length];
+            // var testDataTable = new DataTable(latentVar + " prediction (test)");
+            // var predictedValuesRow = new DataRow(latentVar + " pred.", "Predicted values for " + latentVar, testPrediction.Select(arr => arr[colIdx].Item1).ToArray());
+            // var emptyRow = new DataRow(latentVar);
+            // testDataTable.Rows.Add(emptyRow);
+            // testDataTable.Rows.Add(predictedValuesRow);
+            // testList.Add(testDataTable);
           }
         }
@@ -742 +791 @@
       double[] fi = new double[n];
       double[,] jac = new double[n, d];
-      Integrate(optimizationData, fi, jac);
+      Integrate(optimizationData, fi, jac, null);
       for (int i = 0; i < optimizationData.rows.Length; i++) {
         var res = new Tuple<double, Vector>[nTargets];
@@ -752 +801 @@
     }
 
-    public static void Integrate(OptimizationData optimizationData, double[] fi, double[,] jac) {
+    public static void Integrate(OptimizationData optimizationData, double[] fi, double[,] jac, double[,] latentValues) {
       var trees = optimizationData.trees;
       var dataset = optimizationData.problemData.Dataset;
@@ -795 +844 @@
           outputRowIdx++;
         }
+
+        var latentValueRowIdx = 0;
+        var latentValueColIdx = 0;
+        foreach (var varName in latentVariables) {
+          var y0 = 0.0; // assume we start at zero
+          nodeValues.SetVariableValue(varName, y0, Vector.Zero);
+
+          if (latentValues != null) {
+            latentValues[latentValueRowIdx, latentValueColIdx++] = y0;
+          }
+        }
+        latentValueColIdx = 0; latentValueRowIdx++;
 
         { // CODE BELOW DOESN'T WORK ANYMORE
@@ -822 +883 @@
           prevT = t;
 
-          for (int i = 0; i < calculatedVariables.Length; i++) {
-            var targetVar = calculatedVariables[i];
+          // update output for target variables (TODO: if we want to visualize the latent variables then we need to provide a separate output)
+          for (int i = 0; i < targetVariables.Length; i++) {
+            var targetVar = targetVariables[i];
             var yt = nodeValues.GetVariableValue(targetVar);
 
@@ -829 +891 @@
             if (double.IsNaN(yt.Item1) || double.IsInfinity(yt.Item1)) {
               for (; outputRowIdx < fi.Length; outputRowIdx++) {
-                var prevIdx = outputRowIdx - calculatedVariables.Length;
+                var prevIdx = outputRowIdx - targetVariables.Length;
                 fi[outputRowIdx] = fi[prevIdx]; // current <- prev
                 if (jac != null) for (int j = 0; j < jac.GetLength(1); j++) jac[outputRowIdx, j] = jac[prevIdx, j];
@@ -840 +902 @@
             g.CopyTo(jac, outputRowIdx);
             outputRowIdx++;
+          }
+          if (latentValues != null) {
+            foreach (var latentVariable in latentVariables) {
+              var lt = nodeValues.GetVariableValue(latentVariable).Item1;
+              latentValues[latentValueRowIdx, latentValueColIdx++] = lt;
+            }
+            latentValueRowIdx++; latentValueColIdx = 0;
           }
 
@@ -1451 +1520 @@
 
     #region  helper
+
+    private static IEnumerable<T> EveryNth<T>(IEnumerable<T> xs, int step) {
+      var e = xs.GetEnumerator();
+      while (e.MoveNext()) {
+        for (int i = 0; i < step; i++) {
+          if (!e.MoveNext()) yield break;
+        }
+        yield return e.Current;
+      }
+    }
 
     private void InitAllParameters() {

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

@@ -86 +86 @@
     }
 
-    public IEnumerable<double[]> Predict(IntRange episode, int forecastHorizon, out double snmse) {
+    public IEnumerable<double[]> Predict(IntRange episode, int forecastHorizon) {
       var forecastEpisode = new IntRange(episode.Start, episode.End + forecastHorizon);
       // 
@@ -92 +92 @@
       // snmse = Problem.OptimizeForEpisodes(trees, problemData, targetVars, latentVariables, random, new[] { forecastEpisode }, 100, numericIntegrationSteps, odeSolver);
 
-      var inputVariables = trees.SelectMany(t => t.IterateNodesPrefix().OfType<VariableTreeNode>().Select(n => n.VariableName)).Except(targetVars);
+      var inputVariables = trees.SelectMany(t => t.IterateNodesPrefix().OfType<VariableTreeNode>().Select(n => n.VariableName))
+        .Except(targetVars)
+        .Except(latentVariables)
+        .Distinct();
 
       var optimizationData = new Problem.OptimizationData(trees, targetVars, inputVariables.ToArray(), problemData, null, new[] { forecastEpisode }, numericIntegrationSteps, latentVariables, odeSolver);
@@ -99 +102 @@
       // var theta = Problem.ExtractParametersFromTrees(trees);
 
-      snmse = 0.0; // TODO
-      return Problem.Integrate(optimizationData).Select(p => p.Select(pi => pi.Item1).ToArray()).ToArray();
+
+      var fi = new double[forecastEpisode.Size * targetVars.Length];
+      var jac = new double[forecastEpisode.Size * targetVars.Length, optimizationData.nodeValueLookup.ParameterCount];
+      var latentValues = new double[forecastEpisode.Size, LatentVariables.Length];
+      Problem.Integrate(optimizationData, fi, jac, latentValues);
+      for (int i = 0; i < forecastEpisode.Size; i++) {
+        var res = new double[targetVars.Length + latentVariables.Length];
+        for (int j = 0; j < targetVars.Length; j++) {
+          res[j] = fi[i * targetVars.Length + j];
+        }
+        for (int j = 0; j < latentVariables.Length; j++) {
+          res[targetVars.Length + j] = latentValues[i, j];
+        }
+        yield return res;
+      }
     }
   }

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

@@ -73 +73 @@
 
       double snmse;
-      var predictions = Content.Predict(predictionEpisode, forecastHorizon, out snmse).ToArray();
-      errorLabel.Text = $"SNMSE: {snmse:e4}";
+      var predictions = Content.Predict(predictionEpisode, forecastHorizon).ToArray();
+      errorLabel.Text = $"SNMSE: {0.0:e4}";
       var trainingPredictions = predictions.Take(predictionEpisode.Size).ToArray();
       var forecastPredictions = predictions.Skip(predictionEpisode.Size).ToArray();