Changeset 16597

Timestamp:

02/11/19 14:15:47 (6 years ago)

Author:

gkronber

Message:

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

Location:

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

Files:

• OdeParameterIdentification.cs (modified) (1 diff)
• Problem.cs (modified) (17 diffs)
• SolutionView.Designer.cs (modified) (1 diff)
• Vector.cs (modified) (5 diffs)

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

@@ -140 +140 @@
     #region nonlinear regression
     protected override void Run(CancellationToken cancellationToken) {
-      IRegressionSolution bestSolution = null;
       if (SetSeedRandomly) Seed = (new System.Random()).Next();
       var rand = new MersenneTwister((uint)Seed);

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

@@ -274 +274 @@
         alglib.minlbfgscreate(Math.Min(theta.Length, 5), theta, out state);
         alglib.minlbfgssetcond(state, 0.0, 0.0, 0.0, maxParameterOptIterations);
-        //alglib.minlbfgssetgradientcheck(state, 1e-6);
+        // alglib.minlbfgssetgradientcheck(state, 1e-4);
         alglib.minlbfgsoptimize(state, EvaluateObjectiveAndGradient, null,
           new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver }); //TODO: create a type
@@ -306 +306 @@
                           * NFEV countains number of function calculations
          */
-        if (report.terminationtype < 0) { nmse = 10E6; return; }
+        if (report.terminationtype < 0) { nmse = 10.0; return; }
       }
 
@@ -314 +314 @@
       EvaluateObjectiveAndGradient(optTheta, ref nmse, grad,
         new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver });
-      if (double.IsNaN(nmse) || double.IsInfinity(nmse)) { nmse = 10E6; return; } // return a large value (TODO: be consistent by using NMSE)
+      if (double.IsNaN(nmse) || double.IsInfinity(nmse)) { nmse = 10.0; return; } // return a large value (TODO: be consistent by using NMSE)
     }
 
@@ -342 +342 @@
 
       if (predicted.Length != targetValues.GetLength(0)) {
-        f = double.MaxValue;
+        f = 10.0; // TODO
         Array.Clear(grad, 0, grad.Length);
         return;
@@ -349 +349 @@
       // 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.Variance())
-        .Select(v => 1.0 / v)
-        .ToArray();
+      // 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
@@ -370 +373 @@
           var res = (y - y_pred_f);
           var ressq = res * res;
-          f += ressq * invN * invVar[c];
-          g += -2.0 * res * y_pred[c].Item2 * invN * invVar[c];
+          f += ressq * invN * invVar[c] /* * Math.Exp(-0.2 * r) */ ;
+          g += -2.0 * res * y_pred[c].Item2 * invN * invVar[c] /* * Math.Exp(-0.2 * r) */;
         }
         r++;
@@ -396 +399 @@
       if (!results.ContainsKey("Solution")) {
         results.Add(new Result("Solution", typeof(Solution)));
+      }
+      if (!results.ContainsKey("Squared error and gradient")) {
+        results.Add(new Result("Squared error and gradient", typeof(DataTable)));
       }
 
@@ -477 +483 @@
             trainingDataTable.Rows.Add(actualValuesRow);
             trainingDataTable.Rows.Add(predictedValuesRow);
+
+            for (int paramIdx = 0; paramIdx < optTheta.Length; paramIdx++) {
+              var paramSensitivityRow = new DataRow($"∂{targetVar}/∂θ{paramIdx}", $"Sensitivities of parameter {paramIdx}", trainingPrediction.Select(arr => arr[colIdx].Item2[paramIdx]).ToArray());
+              paramSensitivityRow.VisualProperties.SecondYAxis = true;
+              trainingDataTable.Rows.Add(paramSensitivityRow);
+            }
             trainingList.Add(trainingDataTable);
           } else {
@@ -488 +500 @@
           }
         }
+
+        var errorTable = new DataTable("Squared error and gradient");
+        var seRow = new DataRow("Squared error");
+        var gradientRows = Enumerable.Range(0, optTheta.Length).Select(i => new DataRow($"∂SE/∂θ{i}")).ToArray();
+        errorTable.Rows.Add(seRow);
+        foreach (var gRow in gradientRows) {
+          gRow.VisualProperties.SecondYAxis = true;
+          errorTable.Rows.Add(gRow);
+        }
+        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
+          double f_i = 0.0;
+          Vector g_i = Vector.CreateNew(new double[optTheta.Length]);
+          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);
+            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)*/;
+          }
+          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
         var testList = new ItemList<DataTable>();
@@ -659 +702 @@
             IntegrateHL(trees, calculatedVariables, variableValues, parameterValues, numericIntegrationSteps);
           else if (odeSolver == "CVODES")
-            IntegrateCVODES(trees, calculatedVariables, variableValues, parameterValues, t - prevT);
+            throw new NotImplementedException();
+          // IntegrateCVODES(trees, calculatedVariables, variableValues, parameterValues, t - prevT);
           else throw new InvalidOperationException("Unknown ODE solver " + odeSolver);
           prevT = t;
@@ -687 +731 @@
     #region CVODES
 
-
+    /* 
     /// <summary>
     ///  Here we use CVODES to solve the ODE. Forward sensitivities are used to calculate the gradient for parameter optimization
@@ -981 +1025 @@
         };
     }
+    */
     #endregion
 
@@ -1012 +1057 @@
       }
 
+      double[] deltaF = new double[calculatedVariables.Length];
+      Vector[] deltaG = new Vector[calculatedVariables.Length];
 
       double h = 1.0 / numericIntegrationSteps;
       for (int step = 0; step < numericIntegrationSteps; step++) {
-        var deltaValues = new Dictionary<string, Tuple<double, Vector>>();
+        //var deltaValues = new Dictionary<string, Tuple<double, Vector>>();
         for (int i = 0; i < trees.Length; i++) {
           var tree = trees[i];
@@ -1021 +1068 @@
 
           // Root.GetSubtree(0).GetSubtree(0) skips programRoot and startSymbol
-          var res = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues);
-          deltaValues.Add(targetVarName, res);
+          double f; Vector g;
+          InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues, out f, out g);
+          deltaF[i] = f;
+          deltaG[i] = g;
         }
 
         // update variableValues for next step, trapezoid integration
-        foreach (var kvp in deltaValues) {
-          var oldVal = variableValues[kvp.Key];
+        for (int i = 0; i < trees.Length; i++) {
+          var varName = calculatedVariables[i];
+          var oldVal = variableValues[varName];
           var newVal = Tuple.Create(
-            oldVal.Item1 + h * kvp.Value.Item1,
-            oldVal.Item2 + h * kvp.Value.Item2
+            oldVal.Item1 + h * deltaF[i],
+            oldVal.Item2 + deltaG[i].Scale(h)
           );
-          variableValues[kvp.Key] = newVal;
-        }
-
-
+          variableValues[varName] = newVal;
+        }
+
+        // TODO perf
         foreach (var node in nodeValues.Keys.ToArray()) {
           if (node.SubtreeCount == 0 && !IsConstantNode(node)) {
@@ -1046 +1096 @@
     }
 
-    private static Tuple<double, Vector> InterpretRec(
+    private static void InterpretRec(
       ISymbolicExpressionTreeNode node,
-      Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>> nodeValues      // contains value and gradient vector for a node (variables and constants only)
-        ) {
-
+      Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>> nodeValues,      // contains value and gradient vector for a node (variables and constants only)
+      out double f,
+      out Vector g
+      ) {
+      double fl, fr;
+      Vector gl, gr;
       switch (node.Symbol.Name) {
         case "+": {
-            var l = InterpretRec(node.GetSubtree(0), nodeValues);
-            var r = InterpretRec(node.GetSubtree(1), nodeValues);
-
-            return Tuple.Create(l.Item1 + r.Item1, l.Item2 + r.Item2);
+            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
+            InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
+            f = fl + fr;
+            g = Vector.AddTo(gl, gr);
+            break;
           }
         case "*": {
-            var l = InterpretRec(node.GetSubtree(0), nodeValues);
-            var r = InterpretRec(node.GetSubtree(1), nodeValues);
-
-            return Tuple.Create(l.Item1 * r.Item1, l.Item2 * r.Item1 + l.Item1 * r.Item2);
+            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
+            InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
+            f = fl * fr;
+            g = Vector.AddTo(gl.Scale(fr), gr.Scale(fl)); // f'*g + f*g'
+            break;
           }
 
         case "-": {
-            var l = InterpretRec(node.GetSubtree(0), nodeValues);
-            var r = InterpretRec(node.GetSubtree(1), nodeValues);
-
-            return Tuple.Create(l.Item1 - r.Item1, l.Item2 - r.Item2);
+            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
+            InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
+            f = fl - fr;
+            g = Vector.Subtract(gl, gr);
+            break;
           }
         case "%": {
-            var l = InterpretRec(node.GetSubtree(0), nodeValues);
-            var r = InterpretRec(node.GetSubtree(1), nodeValues);
+            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
+            InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
 
             // protected division
-            if (r.Item1.IsAlmost(0.0)) {
-              return Tuple.Create(0.0, Vector.Zero);
+            if (fr.IsAlmost(0.0)) {
+              f = 0;
+              g = Vector.Zero;
             } else {
-              return Tuple.Create(
-                l.Item1 / r.Item1,
-                l.Item1 * -1.0 / (r.Item1 * r.Item1) * r.Item2 + 1.0 / r.Item1 * l.Item2 // (f/g)' = f * (1/g)' + 1/g * f' = f * -1/g² * g' + 1/g * f'
-                );
+              f = fl / fr;
+              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'
             }
+            break;
           }
         case "sin": {
-            var x = InterpretRec(node.GetSubtree(0), nodeValues);
-            return Tuple.Create(
-              Math.Sin(x.Item1),
-              Vector.Cos(x.Item2) * x.Item2
-            );
+            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
+            f = Math.Sin(fl);
+            g = gl.Scale(Math.Cos(fl));
+            break;
           }
         case "cos": {
-            var x = InterpretRec(node.GetSubtree(0), nodeValues);
-            return Tuple.Create(
-              Math.Cos(x.Item1),
-              -Vector.Sin(x.Item2) * x.Item2
-            );
+            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
+            f = Math.Cos(fl);
+            g = gl.Scale(-Math.Sin(fl));
+            break;
           }
         case "sqr": {
-            var x = InterpretRec(node.GetSubtree(0), nodeValues);
-            return Tuple.Create(
-              x.Item1 * x.Item1,
-              2.0 * x.Item1 * x.Item2
-            );
+            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
+            f = fl * fl;
+            g = gl.Scale(2.0 * fl);
+            break;
           }
         default: {
-            return nodeValues[node];  // value and gradient for constants and variables must be set by the caller
+            var t = nodeValues[node];
+            f = t.Item1;
+            g = Vector.CreateNew(t.Item2);
+            break;
           }
       }
@@ -1223 +1279 @@
       var newVariablesList = new CheckedItemList<StringValue>(ProblemData.Dataset.VariableNames.Select(str => new StringValue(str).AsReadOnly()).ToArray()).AsReadOnly();
       var matchingItems = newVariablesList.Where(item => currentlySelectedVariables.Contains(item.Value)).ToArray();
-      foreach (var matchingItem in matchingItems) {
-        newVariablesList.SetItemCheckedState(matchingItem, true);
+      foreach (var item in newVariablesList) {
+        if (currentlySelectedVariables.Contains(item.Value)) {
+          newVariablesList.SetItemCheckedState(item, true);
+        } else {
+          newVariablesList.SetItemCheckedState(item, false);
+        }
       }
       TargetVariablesParameter.Value = newVariablesList;
@@ -1244 +1304 @@
       // whenever ProblemData is changed we create a new grammar with the necessary symbols
       var g = new SimpleSymbolicExpressionGrammar();
-      g.AddSymbols(FunctionSet.CheckedItems.OrderBy(i => i.Index).Select(i => i.Value.Value).ToArray(), 2, 2);
-
-      // TODO
-      //g.AddSymbols(new[] {
-      //  "exp",
-      //  "log", // log( <expr> ) // TODO: init a theta to ensure the value is always positive
-      //  "exp_minus" // exp((-1) * <expr>
-      //}, 1, 1);
+      var unaryFunc = new string[] { "sin", "cos", "sqr" };
+      var binaryFunc = new string[] { "+", "-", "*", "%" };
+      foreach (var func in unaryFunc) {
+        if (FunctionSet.CheckedItems.Any(ci => ci.Value.Value == func)) g.AddSymbol(func, 1, 1);
+      }
+      foreach (var func in binaryFunc) {
+        if (FunctionSet.CheckedItems.Any(ci => ci.Value.Value == func)) g.AddSymbol(func, 2, 2);
+      }
 
       foreach (var variableName in ProblemData.AllowedInputVariables.Union(TargetVariables.CheckedItems.Select(i => i.Value.Value)))

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

@@ -90 +90 @@
       this.forecastTextbox.Size = new System.Drawing.Size(100, 20);
       this.forecastTextbox.TabIndex = 2;
-      this.forecastTextbox.Text = "180";
+      this.forecastTextbox.Text = "1";
       // 
       // forecastLabel

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

@@ -9 +9 @@
       if (a == Zero) return b;
       if (b == Zero) return a;
-      Debug.Assert(a.arr.Length == b.arr.Length);
+      Trace.Assert(a.arr.Length == b.arr.Length);
       var res = new double[a.arr.Length];
       for (int i = 0; i < res.Length; i++)
@@ -16 +16 @@
     }
 
+    public static Vector AddTo(Vector a, Vector b) {
+      if (b == Zero) return a;
+      if (a == Zero) {
+        var vArr = new double[b.Length];
+        b.CopyTo(vArr);
+        return new Vector(vArr);
+      } else {
+        Trace.Assert(a.arr.Length == b.arr.Length);
+        for (int i = 0; i < a.arr.Length; i++)
+          a.arr[i] += b.arr[i];
+        return a;
+      }
+    }
+
+
+    public static Vector Subtract(Vector a, Vector b) {
+      if (b == Zero) return a;
+      if (a == Zero) {
+        var vArr = new double[b.Length];
+        a = new Vector(vArr);
+      }
+      Trace.Assert(a.arr.Length == b.arr.Length);
+      for (int i = 0; i < a.arr.Length; i++)
+        a.arr[i] -= b.arr[i];
+      return a;
+    }
+
     public static Vector operator -(Vector a, Vector b) {
       if (b == Zero) return a;
       if (a == Zero) return -b;
-      Debug.Assert(a.arr.Length == b.arr.Length);
+      Trace.Assert(a.arr.Length == b.arr.Length);
       var res = new double[a.arr.Length];
       for (int i = 0; i < res.Length; i++)
@@ -39 +66 @@
         res[i] = s * v.arr[i];
       return new Vector(res);
+    }
+
+    public Vector Scale(double s) {
+      if (this != Zero) {
+        for (int i = 0; i < arr.Length; i++) {
+          arr[i] *= s;
+        }
+      }
+
+      return this;
     }
 
@@ -95 +132 @@
     private readonly double[] arr; // backing array;
 
+    /// <summary>
+    /// 
+    /// </summary>
+    /// <param name="v">Is not cloned!</param>
     public Vector(double[] v) {
       this.arr = v;
@@ -100 +141 @@
 
     public void CopyTo(double[] target) {
-      Debug.Assert(arr.Length <= target.Length);
+      Trace.Assert(arr.Length <= target.Length);
       Array.Copy(arr, target, arr.Length);
+    }
+
+
+    /// <summary>
+    /// Creates a new vector
+    /// </summary>
+    /// <param name="a"> is cloned</param>
+    /// <returns></returns>
+    public static Vector CreateNew(double[] a) {
+      var arr = new double[a.Length];
+      Array.Copy(a, arr, arr.Length);
+      return new Vector(arr);
+    }
+
+    /// <summary>
+    /// Creates a new vector
+    /// </summary>
+    /// <param name="v"> is cloned</param>
+    /// <returns></returns>
+    public static Vector CreateNew(Vector v) {
+      if (v == Zero) return Zero;
+      var arr = new double[v.arr.Length];
+      Array.Copy(v.arr, arr, arr.Length);
+      return new Vector(arr);
     }
   }