Changeset 16604 for branches/2925_AutoDiffForDynamicalModels

Timestamp:

02/13/19 17:29:49 (6 years ago)

Author:

gkronber

Message:

#2925: efficiency improvements

Location:

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

Files:

• Problem.cs (modified) (18 diffs)
• Vector.cs (modified) (5 diffs)

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

@@ -428 +428 @@
           var y = optimizationData.targetValues[i][trainIdx];
           fi[termIdx] = (y - f) * optimizationData.inverseStandardDeviation[i]; // scale of NMSE
-          if (jac != null && g != Vector.Zero) for (int j = 0; j < g.Length; j++) jac[termIdx, j] = -g[j] * optimizationData.inverseStandardDeviation[i]; 
+          if (jac != null && g != Vector.Zero) for (int j = 0; j < g.Length; j++) jac[termIdx, j] = -g[j] * optimizationData.inverseStandardDeviation[i];
 
           termIdx++;
@@ -457 +457 @@
       nodeValueLookup.UpdateParamValues(x);
 
-      var predictionEnumerator = Integrate(optimizationData).GetEnumerator();
+      Integrate(optimizationData, fi, jac);
       var trees = optimizationData.trees;
 
-      predictionEnumerator.MoveNext();
-      var currentPrediction = predictionEnumerator.Current;
+      // update result with error
       for (int trainIdx = 0; trainIdx < rows.Length; trainIdx++) {
-        var pred_i = currentPrediction;
-
         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) * optimizationData.inverseStandardDeviation[i];  // scale for normalized squared error
-          if (jac != null && g != Vector.Zero) for (int j = 0; j < g.Length; j++) jac[outputIdx, j] = -g[j] * optimizationData.inverseStandardDeviation[i];
-
+          fi[outputIdx] = (y - fi[outputIdx]) * optimizationData.inverseStandardDeviation[i];  // scale for normalized squared error
+          if (jac != null) for (int j = 0; j < x.Length; j++) jac[outputIdx, j] = -jac[outputIdx, j] * optimizationData.inverseStandardDeviation[i];
           outputIdx++;
         }
-        if (predictionEnumerator.MoveNext()) currentPrediction = predictionEnumerator.Current; // if we stopped early => just use the last prediction until the end
       }
     }
@@ -606 +599 @@
             var ressq = res * res;
             f_i += ressq * invN;
-            g_i = g_i - 2.0 * res * y_pred[colIdx].Item2 * invN;
+            g_i.Add(y_pred[colIdx].Item2.Scale(-2.0 * res * invN));
           }
           seRow.Values.Add(f_i);
@@ -652 +645 @@
         var clonedTrees = new List<ISymbolicExpressionTree>();
         for (int idx = 0; idx < trees.Length; idx++) {
-          clonedTrees.Add((ISymbolicExpressionTree)trees[idx]);
+          clonedTrees.Add((ISymbolicExpressionTree)trees[idx].Clone());
         }
         var ds = problemData.Dataset;
@@ -689 +682 @@
     }
 
-
-    public static double[] ExtractParametersFromTrees(ISymbolicExpressionTree[] trees) {
-      return trees
-        .SelectMany(t => t.IterateNodesPrefix().OfType<ConstantTreeNode>().Select(n => n.Value))
-        .ToArray();
-    }
-
-
-
     #region interpretation
 
@@ -711 +695 @@
 
     public static IEnumerable<Tuple<double, Vector>[]> Integrate(OptimizationData optimizationData) {
-
+      var nTargets = optimizationData.targetVariables.Length;
+      var n = optimizationData.rows.Length * optimizationData.targetVariables.Length;
+      var d = optimizationData.nodeValueLookup.ParameterCount;
+      double[] fi = new double[n];
+      double[,] jac = new double[n, d];
+      Integrate(optimizationData, fi, jac);
+      for (int i = 0;i < optimizationData.rows.Length;i++) {
+        var res = new Tuple<double, Vector>[nTargets];
+        for (int j = 0;j<nTargets;j++) {
+          res[j] = Tuple.Create(fi[i * nTargets + j], Vector.CreateFromMatrixRow(jac, i * nTargets + j));
+        }
+        yield return res;
+      }
+    }
+
+    public static void Integrate(OptimizationData optimizationData, double[] fi, double[,] jac) {
       var trees = optimizationData.trees;
       var dataset = optimizationData.problemData.Dataset;
@@ -730 +729 @@
 
         var t0 = rows.First();
+        var outputRowIdx = 0;
 
         // initialize values for inputs and targets from dataset
-        foreach (var varName in inputVariables.Concat(targetVariables)) {
-          nodeValues.SetVariableValue(varName, dataset.GetDoubleValue(varName, t0), Vector.Zero);
+        foreach (var varName in inputVariables) {
+          var y0 = dataset.GetDoubleValue(varName, t0);
+          nodeValues.SetVariableValue(varName, y0, Vector.Zero);
+        }
+        foreach(var varName in targetVariables) {
+          var y0 = dataset.GetDoubleValue(varName, t0);
+          nodeValues.SetVariableValue(varName, y0, Vector.Zero);
+
+          // output starting value
+          fi[outputRowIdx] = y0;
+          Vector.Zero.CopyTo(jac, outputRowIdx);
+
+          outputRowIdx++;
         }
 
@@ -752 +763 @@
         }
 
-        // return first value as stored in the dataset
-        yield return calculatedVariables
-          .Select(calcVarName => nodeValues.GetVariableValue(calcVarName))
-          .ToArray();
-
         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)) {
@@ -767 +773 @@
           prevT = t;
 
-
-          var res = calculatedVariables
-            .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;
+          for (int i = 0; i < calculatedVariables.Length; i++) {
+            var targetVar = calculatedVariables[i];
+            var yt = nodeValues.GetVariableValue(targetVar);
+
+            // fill up remaining rows with last valid value if there are invalid values
+            if (double.IsNaN(yt.Item1) || double.IsInfinity(yt.Item1)) {
+              for (; outputRowIdx < fi.Length; outputRowIdx++) {
+                var prevIdx = outputRowIdx - calculatedVariables.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];
+              }
+              return;
+            };
+
+            fi[outputRowIdx] = yt.Item1;
+            var g = yt.Item2;
+            g.CopyTo(jac, outputRowIdx);
+            outputRowIdx++;
+          }
 
           // update for next time step (only the inputs)
@@ -1087 +1104 @@
       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)
       NodeValueLookup nodeValues,
-      // double[] parameterValues,
       int numericIntegrationSteps) {
 
@@ -1114 +1129 @@
           var varName = calculatedVariables[i];
           var oldVal = nodeValues.GetVariableValue(varName);
-          nodeValues.SetVariableValue(varName, oldVal.Item1 + h * deltaF[i], oldVal.Item2 + deltaG[i].Scale(h));
+          nodeValues.SetVariableValue(varName, oldVal.Item1 + h * deltaF[i], oldVal.Item2.Add(deltaG[i].Scale(h)));
         }
       }
@@ -1123 +1138 @@
       if (node is ConstantTreeNode) {
         return ((ConstantTreeNode)node).Value;
-      } else if(node is VariableTreeNode) {
+      } else if (node is VariableTreeNode) {
         return nodeValues.NodeValue(node);
       } else if (node.Symbol is Addition) {
@@ -1175 +1190 @@
       if (node.Symbol is Constant || node.Symbol is Variable) {
         z = nodeValues.NodeValue(node);
-        dz = Vector.CreateNew(nodeValues.NodeGradient(node));
+        dz = Vector.CreateNew(nodeValues.NodeGradient(node)); // original gradient vectors are never changed by evaluation
       } else if (node.Symbol is Addition) {
         InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
         InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
         z = f + g;
-        dz = df + dg; // Vector.AddTo(gl, gr);
-
+        dz = df.Add(dg);
       } else if (node.Symbol is Multiplication) {
         InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
         InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
         z = f * g;
-        dz = df * g + f * dg;  // Vector.AddTo(gl.Scale(fr), gr.Scale(fl)); // f'*g + f*g'
+        dz = df.Scale(g).Add(dg.Scale(f));  // f'*g + f*g'
 
       } else if (node.Symbol is Subtraction) {
@@ -1192 +1206 @@
           InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
           z = -f;
-          dz = -df;
+          dz = df.Scale(-1.0);
         } else {
           InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
           InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
-
           z = f - g;
-          dz = df - dg;
+          dz = df.Subtract(dg);
         }
 
@@ -1210 +1223 @@
           dz = Vector.Zero;
         } else {
-          z = f / g;
-          dz = -f / (g * g) * dg + df / g; // Vector.AddTo(dg.Scale(f * -1.0 / (g * g)), df.Scale(1.0 / g)); // (f/g)' = f * (1/g)' + 1/g * f' = f * -1/g² * g' + 1/g * f'
+          var inv_g = 1.0 / g;
+          z = f * inv_g;
+
+          dz = dg.Scale(-f * inv_g * inv_g).Add(df.Scale(inv_g));
         }
 
@@ -1217 +1232 @@
         InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
         z = Math.Sin(f);
-        dz = Math.Cos(f) * df;
+        dz = df.Scale(Math.Cos(f));
 
       } else if (node.Symbol is Cosine) {
         InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
         z = Math.Cos(f);
-        dz = -Math.Sin(f) * df;
+        dz = df.Scale(-Math.Sin(f));
       } else if (node.Symbol is Square) {
         InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
         z = f * f;
-        dz = 2.0 * f * df;
+        dz = df.Scale(2.0 * f);
       } else {
         throw new NotSupportedException("unsupported symbol");
@@ -1522 +1537 @@
         this.odeSolver = odeSolver;
         this.nodeValueLookup = new NodeValueLookup(trees);
-        this.rows =  episodes.SelectMany(ep => Enumerable.Range(ep.Start, ep.Size)).ToArray();
+        this.rows = episodes.SelectMany(ep => Enumerable.Range(ep.Start, ep.Size)).ToArray();
       }
     }
@@ -1541 +1556 @@
         this.constantNodes = trees.SelectMany(t => t.IterateNodesPrefix().OfType<ConstantTreeNode>()).ToArray();
         constantGradientVectors = new Vector[constantNodes.Length];
-        for(int paramIdx = 0; paramIdx < constantNodes.Length; paramIdx++) {
+        for (int paramIdx = 0; paramIdx < constantNodes.Length; paramIdx++) {
           constantGradientVectors[paramIdx] = Vector.CreateIndicator(length: constantNodes.Length, idx: paramIdx);
 

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

@@ -17 +17 @@
     }
 
-    // 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 Vector Add(Vector b) {
+      var a = this;
+      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) {
+    public Vector Subtract(Vector b) {
+      var a = this;
       if (b == Zero) return a;
       if (a == Zero) {
@@ -39 +41 @@
       }
       Trace.Assert(a.arr.Length == b.arr.Length);
-      var res = new double[a.arr.Length];
       for (int i = 0; i < a.arr.Length; i++)
-        res[i] = a.arr[i] - b.arr[i];
-      return new Vector(res);
+        a.arr[i] -= b.arr[i];
+      return a;
     }
 
-    public static Vector operator -(Vector a, Vector b) {
-      return Subtract(a, b);
-    }
-    public static Vector operator -(Vector v) {
-      return v.Scale(-1.0);
-    }
+    // public static Vector operator -(Vector a, Vector b) {
+    //   return Subtract(a, b);
+    // }
+    // public static Vector operator -(Vector v) {
+    //   return v.Scale(-1.0);
+    // }
 
-    public static Vector operator *(double s, Vector v) {
-      return v.Scale(s);
-    }
+    //public static Vector operator *(double s, Vector v) {
+    //  return v.Scale(s);
+    //}
 
     public Vector Scale(double s) {
       if (this == Zero) return Zero;
 
-      var res = new double[arr.Length];
       for (int i = 0; i < arr.Length; i++) {
-        res[i] = arr[i] * s;
+        arr[i] *= s;
       }
-      return new Vector(res);
-
+      return this;
     }
 
-    public static Vector operator *(Vector v, double s) {
-      return s * v;
-    }
+    // public static Vector operator *(Vector v, double s) {
+    //   return s * v;
+    // }
 
     public static Vector operator *(Vector u, Vector v) {
@@ -89 +88 @@
     }
 
-    public static Vector operator /(Vector v, double s) {
-      return v.Scale(1.0 / s);
+    // public static Vector operator /(Vector v, double s) {
+    //   return v.Scale(1.0 / s);
+    // }
+
+    public Vector Sin() {
+      for (int i = 0; i < arr.Length; i++) arr[i] = Math.Sin(arr[i]);
+      return this;
     }
-
-    public static Vector Sin(Vector s) {
-      var res = new double[s.arr.Length];
-      for (int i = 0; i < res.Length; i++) res[i] = Math.Sin(s.arr[i]);
-      return new Vector(res);
-    }
-    public static Vector Cos(Vector s) {
-      var res = new double[s.arr.Length];
-      for (int i = 0; i < res.Length; i++) res[i] = Math.Cos(s.arr[i]);
-      return new Vector(res);
+    public Vector Cos() {
+      for (int i = 0; i < arr.Length; i++) arr[i] = Math.Cos(arr[i]);
+      return this;
     }
 
@@ -137 +134 @@
     }
 
+    public void CopyTo(double[,] target, int rowIdx) {
+      if (target == null) return;
+      if (this == Zero) {
+        for (int j = 0; j < target.GetLength(1); j++) target[rowIdx, j] = 0.0;
+      } else {
+        for (int j = 0; j < target.GetLength(1); j++) target[rowIdx, j] = arr[j];
+      }
+    }
     /// <summary>
     /// Creates a new vector
@@ -165 +170 @@
       return new Vector(gArr);
     }
+
+    internal static Vector CreateFromMatrixRow(double[,] jac, int rowIdx) {
+      var arr = new double[jac.GetLength(1)];
+      for (int i = 0; i < arr.Length; i++) arr[i] = jac[rowIdx, i];
+      return new Vector(arr);
+    }
   }
 }