Changeset 16610 for branches/2925_AutoDiffForDynamicalModels

Timestamp:

02/16/19 06:53:44 (6 years ago)

Author:

gkronber

Message:

#2925: added support for multiple training episodes (and log, exp functions)

Location:

branches/2925_AutoDiffForDynamicalModels

Files:

• AutoDiffForDynamicalModelsTest/AutoDiffForDynamicalModelsTest.csproj (modified) (1 diff)
• AutoDiffForDynamicalModelsTest/TestOdeIdentification.cs (modified) (1 diff)
• HeuristicLab.Problems.DynamicalSystemsModelling/3.3/HeuristicLab.Problems.DynamicalSystemsModelling-3.3.csproj (modified) (1 diff)
• HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs (modified) (22 diffs)
• HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Solution.cs (modified) (1 diff)

branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest/AutoDiffForDynamicalModelsTest.csproj

@@ -76 +76 @@
   </ItemGroup>
   <ItemGroup>
-    <Compile Include="TestCvodes.cs" />
     <Compile Include="Properties\AssemblyInfo.cs" />
     <Compile Include="TestOdeIdentification.cs" />

branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest/TestOdeIdentification.cs

@@ -42 +42 @@
       Assert.AreEqual(6.8350792038369173E-20, ((DoubleValue)alg.Results["SNMSE"].Value).Value, 1E-8);
     }
+    [TestMethod]
+    public void RunOdeIdentificationMultipleEpisodes() {
+      var alg = new OdeParameterIdentification();
+      var dynProb = new Problem();
+      var parser = new HeuristicLab.Problems.Instances.DataAnalysis.TableFileParser();
+      // var fileName = @"C:\reps\HEAL\EuroCAST - Kronberger\DataGeneration\test.csv";
+      var fileName = @"C:\reps\HEAL\EuroCAST - Kronberger\DataGeneration\s-system.csv";
+      parser.Parse(fileName, true);
+      var prov = new HeuristicLab.Problems.Instances.DataAnalysis.RegressionCSVInstanceProvider();
+      var regressionProblemData = prov.ImportData(fileName);
+      regressionProblemData.TrainingPartition.Start = 0;
+      regressionProblemData.TrainingPartition.End = 61;
+      regressionProblemData.TestPartition.Start = 61;
+      var allowedInputs = new List<string>(); // empty
+      var targets = new List<string>(new[] { "y1", "y2", "y3", "y4", "y5" });
+      foreach (var checkedItem in regressionProblemData.InputVariables) {
+        regressionProblemData.InputVariables.SetItemCheckedState(checkedItem, allowedInputs.Contains(checkedItem.Value));
+      }
+      dynProb.Load(regressionProblemData);
+
+      foreach (var checkedItem in dynProb.TargetVariables) {
+        dynProb.TargetVariables.SetItemCheckedState(checkedItem, targets.Contains(checkedItem.Value));
+      }
+
+      dynProb.TrainingEpisodesParameter.Value.Add(new IntRange(0, 30));
+      dynProb.TrainingEpisodesParameter.Value.Add(new IntRange(30, 61));
+      dynProb.NumericIntegrationStepsParameter.Value.Value = 1;
+
+      var rand = new FastRandom(1234);
+      var expressions = new string[] {
+"1.0*(y3*exp(log(y5)*0.1)) + 1.0*(y1*y1)",
+"1.0*(y1*y1) + 1.0*(y2*y2)",
+"1.0*exp(log(y2)*0.1) + 1.0*exp(log(y2)*0.1)*(y3*y3)",
+"1.0*(y1*y1)*exp(log(y5)*0.1) + 1.0*(y4*y4)",
+"1.0*(y4*y4)+1.0*(y5*y5)"
+      };
+
+      alg.CreateSolution(dynProb, expressions, 200, rand);
+      Assert.AreEqual(0.0, ((DoubleValue)alg.Results["SNMSE"].Value).Value, 1E-8);
+    }
   }
 }

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

@@ -113 +113 @@
   </ItemGroup>
   <ItemGroup>
-    <Compile Include="CVODES.cs" />
     <Compile Include="OdeParameterIdentification.cs" />
     <Compile Include="Plugin.cs" />

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

@@ -278 +278 @@
         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 targetValuesDiff = new List<double>();
+        foreach (var ep in episodes) {
+          var episodeRows = Enumerable.Range(ep.Start, ep.Size);
+          var targetValues = problemData.Dataset.GetDoubleValues(targetVars[treeIdx], episodeRows).ToArray();
+          targetValuesDiff.AddRange(targetValues.Skip(1).Zip(targetValues, (t1, t0) => t1 - t0));// TODO: smoothing or multi-pole);
+        }
+        var adjustedEpisodes = episodes.Select(ep => new IntRange(ep.Start, ep.End - 1)); // because we lose the last row in the differencing step
+        var myState = new OptimizationData(new[] { t },
+          targetVars,
+          problemData.AllowedInputVariables.Concat(targetVars).ToArray(),
+          problemData, new[] { targetValuesDiff.ToArray() }, adjustedEpisodes.ToArray(), -99, latentVariables, string.Empty); // TODO
         var paramCount = myState.nodeValueLookup.ParameterCount;
 
@@ -295 +300 @@
             var upperBounds = Enumerable.Repeat(100.0, p.Length).ToArray();
             Array.Copy(initialTheta[treeIdx], p, p.Length);
-            alglib.minlmcreatevj(targetValuesDiff.Length, p, out state);
+            alglib.minlmcreatevj(targetValuesDiff.Count, p, out state);
             alglib.minlmsetcond(state, 0.0, 0.0, 0.0, maxParameterOptIterations);
             alglib.minlmsetbc(state, lowerBounds, upperBounds);
@@ -325 +330 @@
     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 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 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);
+      var myState = new OptimizationData(trees, targetVars, problemData.AllowedInputVariables.ToArray(), problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver);
       optTheta = initialTheta;
 
@@ -378 +382 @@
       var nodeValueLookup = optimizationData.nodeValueLookup;
       var ds = problemData.Dataset;
+      var variables = optimizationData.variables;
+
+      nodeValueLookup.UpdateParamValues(x);
+
       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]));
+        foreach (var variable in variables) {
+          nodeValueLookup.SetVariableValue(variable, ds.GetDoubleValue(variable, rows[trainIdx])); // TODO: perf
         }
         // interpret all trees
@@ -404 +409 @@
       var ds = problemData.Dataset;
       var rows = optimizationData.rows;
+      var variables = optimizationData.variables;
 
       var nodeValueLookup = optimizationData.nodeValueLookup;
@@ -412 +418 @@
       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]));
+        foreach (var variable in variables) {
+          nodeValueLookup.SetVariableValue(variable, ds.GetDoubleValue(variable, rows[trainIdx])); // TODO: perf
         }
 
@@ -511 +517 @@
         foreach (var episode in TrainingEpisodes) {
           var episodes = new[] { episode };
-          var optimizationData = new OptimizationData(trees, targetVars, problemData, null, episodes, NumericIntegrationSteps, latentVariables, OdeSolver);
+          var optimizationData = new OptimizationData(trees, targetVars, problemData.AllowedInputVariables.ToArray(), problemData, null, episodes, NumericIntegrationSteps, latentVariables, OdeSolver);
           var trainingPrediction = Integrate(optimizationData).ToArray();
           trainingPredictions.Add(trainingPrediction);
@@ -543 +549 @@
         results["Models"].Value = models;
       } else {
-        var optimizationData = new OptimizationData(trees, targetVars, problemData, null, TrainingEpisodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver);
+        var optimizationData = new OptimizationData(trees, targetVars, problemData.AllowedInputVariables.ToArray(), problemData, null, TrainingEpisodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver);
         var trainingPrediction = Integrate(optimizationData).ToArray();
 
@@ -549 +555 @@
         var numParams = optimizationData.nodeValueLookup.ParameterCount;
         // for target values and latent variables
-        var trainingRows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start));
+        var trainingRows = optimizationData.rows;
         for (int colIdx = 0; colIdx < trees.Length; colIdx++) {
           // is target variable
@@ -587 +593 @@
         var targetValues = targetVars.Select(v => problemData.Dataset.GetDoubleValues(v, trainingRows).ToArray()).ToArray();
         int r = 0;
-        double invN = 1.0 / trainingRows.Count();
+
         foreach (var y_pred in trainingPrediction) {
           // calculate objective function gradient
@@ -598 +604 @@
             var res = (y - y_pred_f);
             var ressq = res * res;
-            f_i += ressq * invN;
-            g_i.Add(y_pred[colIdx].Item2.Scale(-2.0 * res * invN));
+            f_i += ressq * optimizationData.inverseStandardDeviation[colIdx];
+            g_i.Add(y_pred[colIdx].Item2.Scale(-2.0 * res * optimizationData.inverseStandardDeviation[colIdx]));
           }
           seRow.Values.Add(f_i);
@@ -610 +616 @@
         var testList = new ItemList<DataTable>();
         var testRows = ProblemData.TestIndices.ToArray();
-        var testOptimizationData = new OptimizationData(trees, targetVars, problemData, null, new IntRange[] { ProblemData.TestPartition }, NumericIntegrationSteps, latentVariables, OdeSolver);
+        var testOptimizationData = new OptimizationData(trees, targetVars, problemData.AllowedInputVariables.ToArray(), problemData, null, new IntRange[] { ProblemData.TestPartition }, NumericIntegrationSteps, latentVariables, OdeSolver);
         var testPrediction = Integrate(testOptimizationData).ToArray();
 
@@ -701 +707 @@
       double[,] jac = new double[n, d];
       Integrate(optimizationData, fi, jac);
-      for (int i = 0;i < optimizationData.rows.Length;i++) {
+      for (int i = 0; i < optimizationData.rows.Length; i++) {
         var res = new Tuple<double, Vector>[nTargets];
-        for (int j = 0;j<nTargets;j++) {
+        for (int j = 0; j < nTargets; j++) {
           res[j] = Tuple.Create(fi[i * nTargets + j], Vector.CreateFromMatrixRow(jac, i * nTargets + j));
         }
@@ -721 +727 @@
       var calculatedVariables = targetVariables.Concat(latentVariables).ToArray(); // TODO: must conincide with the order of trees in the encoding
 
+
+
       var nodeValues = optimizationData.nodeValueLookup;
 
       // TODO: numericIntegrationSteps is only relevant for the HeuristicLab solver
+      var outputRowIdx = 0;
       var episodeIdx = 0;
       foreach (var episode in optimizationData.episodes) {
@@ -729 +738 @@
 
         var t0 = rows.First();
-        var outputRowIdx = 0;
 
         // initialize values for inputs and targets from dataset
@@ -736 +744 @@
           nodeValues.SetVariableValue(varName, y0, Vector.Zero);
         }
-        foreach(var varName in targetVariables) {
+        foreach (var varName in targetVariables) {
           var y0 = dataset.GetDoubleValue(varName, t0);
           nodeValues.SetVariableValue(varName, y0, Vector.Zero);
@@ -1141 +1149 @@
         return nodeValues.NodeValue(node);
       } else if (node.Symbol is Addition) {
+        Debug.Assert(node.SubtreeCount == 2);
         var f = InterpretRec(node.GetSubtree(0), nodeValues);
         var g = InterpretRec(node.GetSubtree(1), nodeValues);
         return f + g;
       } else if (node.Symbol is Multiplication) {
+        Debug.Assert(node.SubtreeCount == 2);
         var f = InterpretRec(node.GetSubtree(0), nodeValues);
         var g = InterpretRec(node.GetSubtree(1), nodeValues);
         return f * g;
       } else if (node.Symbol is Subtraction) {
+        Debug.Assert(node.SubtreeCount <= 2);
         if (node.SubtreeCount == 1) {
           var f = InterpretRec(node.GetSubtree(0), nodeValues);
@@ -1159 +1170 @@
         }
       } else if (node.Symbol is Division) {
-        var f = InterpretRec(node.GetSubtree(0), nodeValues);
-        var g = InterpretRec(node.GetSubtree(1), nodeValues);
-
-        // protected division
-        if (g.IsAlmost(0.0)) {
-          return 0;
+        Debug.Assert(node.SubtreeCount <= 2);
+
+        if (node.SubtreeCount == 1) {
+          var f = InterpretRec(node.GetSubtree(0), nodeValues);
+          // protected division
+          if (f.IsAlmost(0.0)) {
+            return 0;
+          } else {
+            return 1.0 / f;
+          }
         } else {
-          return f / g;
+          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;
+          }
         }
       } else if (node.Symbol is Sine) {
+        Debug.Assert(node.SubtreeCount == 1);
+
         var f = InterpretRec(node.GetSubtree(0), nodeValues);
         return Math.Sin(f);
       } else if (node.Symbol is Cosine) {
+        Debug.Assert(node.SubtreeCount == 1);
+
         var f = InterpretRec(node.GetSubtree(0), nodeValues);
         return Math.Cos(f);
       } else if (node.Symbol is Square) {
+        Debug.Assert(node.SubtreeCount == 1);
+
         var f = InterpretRec(node.GetSubtree(0), nodeValues);
         return f * f;
+      } else if (node.Symbol is Exponential) {
+        Debug.Assert(node.SubtreeCount == 1);
+
+        var f = InterpretRec(node.GetSubtree(0), nodeValues);
+        return Math.Exp(f);
+      } else if (node.Symbol is Logarithm) {
+        Debug.Assert(node.SubtreeCount == 1);
+
+        var f = InterpretRec(node.GetSubtree(0), nodeValues);
+        return Math.Log(f);
       } else throw new NotSupportedException("unsupported symbol");
     }
@@ -1215 +1254 @@
 
       } else if (node.Symbol is Division) {
-        InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
-        InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
-
-        // protected division
-        if (g.IsAlmost(0.0)) {
-          z = 0;
-          dz = Vector.Zero;
+        if (node.SubtreeCount == 1) {
+          InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
+          // protected division
+          if (f.IsAlmost(0.0)) {
+            z = 0;
+            dz = Vector.Zero;
+          } else {
+            z = 1.0 / f;
+            dz = df.Scale(z * z);
+          }
         } else {
-          var inv_g = 1.0 / g;
-          z = f * inv_g;
-
-          dz = dg.Scale(-f * inv_g * inv_g).Add(df.Scale(inv_g));
+          InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
+          InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
+
+          // protected division
+          if (g.IsAlmost(0.0)) {
+            z = 0;
+            dz = Vector.Zero;
+          } else {
+            var inv_g = 1.0 / g;
+            z = f * inv_g;
+
+            dz = dg.Scale(-f * inv_g * inv_g).Add(df.Scale(inv_g));
+          }
         }
 
@@ -1242 +1293 @@
         z = f * f;
         dz = df.Scale(2.0 * f);
+      } else if (node.Symbol is Exponential) {
+        InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
+        z = Math.Exp(f);
+        dz = df.Scale(Math.Exp(f));
+      } else if (node.Symbol is Logarithm) {
+        InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
+        z = Math.Log(f);
+        dz = df.Scale(1.0 / f);
       } else {
         throw new NotSupportedException("unsupported symbol");
@@ -1518 +1577 @@
       public readonly NodeValueLookup nodeValueLookup;
       public readonly int[] rows;
-
-      public OptimizationData(ISymbolicExpressionTree[] trees, string[] targetVars, IRegressionProblemData problemData,
+      internal readonly string[] variables;
+
+      public OptimizationData(ISymbolicExpressionTree[] trees, string[] targetVars, string[] inputVariables,
+        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;
         this.problemData = problemData;
         this.targetValues = targetValues;
-        if (targetValues != null)
-          this.inverseStandardDeviation = targetValues.Select(vec => 1.0 / vec.StandardDeviation()).ToArray();
-        else
-          this.inverseStandardDeviation = Enumerable.Repeat(1.0, trees.Length).ToArray();
+        this.variables = inputVariables;
+        // if (targetValues != null) {
+        //   this.inverseStandardDeviation = new double[targetValues.Length];
+        //   for (int i = 0; i < targetValues.Length; i++) {
+        //     // calculate variance for each episode separately and calc the average 
+        //     var epStartIdx = 0;
+        //     var stdevs = new List<double>();
+        //     foreach (var ep in episodes) {
+        //       var epValues = targetValues[i].Skip(epStartIdx).Take(ep.Size);
+        //       stdevs.Add(epValues.StandardDeviation());
+        //       epStartIdx += ep.Size;
+        //     }
+        //     inverseStandardDeviation[i] = 1.0 / stdevs.Average();
+        //   }
+        // } else
+        this.inverseStandardDeviation = Enumerable.Repeat(1.0, trees.Length).ToArray();
         this.episodes = episodes;
         this.numericIntegrationSteps = numericIntegrationSteps;
@@ -1589 +1661 @@
         } else {
           var fakeNode = new VariableTreeNode(new Variable());
+          fakeNode.Weight = 1.0;
+          fakeNode.VariableName = variableName;
           var newNodeList = new List<ISymbolicExpressionTreeNode>();
           newNodeList.Add(fakeNode);

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

@@ -91 +91 @@
       // var random = new FastRandom(12345);
       // snmse = Problem.OptimizeForEpisodes(trees, problemData, targetVars, latentVariables, random, new[] { forecastEpisode }, 100, numericIntegrationSteps, odeSolver);
-      var optimizationData = new Problem.OptimizationData(trees, targetVars, problemData, null, new[] { forecastEpisode }, numericIntegrationSteps, latentVariables, odeSolver);
+      var optimizationData = new Problem.OptimizationData(trees, targetVars, problemData.AllowedInputVariables.ToArray(), problemData, null, new[] { forecastEpisode }, numericIntegrationSteps, latentVariables, odeSolver);
       // 
       //