Index: /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs
===================================================================
--- /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs	(revision 16599)
+++ /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs	(revision 16600)
@@ -250,5 +250,5 @@
       var targetValues = new double[rows.Length, targetVars.Length];
 
-      
+
       // collect values of all target variables                                            
       var colIdx = 0;
@@ -275,5 +275,5 @@
       var theta = new double[paramNodes.Count + latentVariables.Length * episodes.Count()];
       for (int i = 0; i < theta.Length; i++)
-        theta[i] = random.NextDouble() * 2.0e-2 - 1.0e-2;
+        theta[i] = random.NextDouble() * 2.0e-1 - 1.0e-1;
 
       optTheta = new double[0];
@@ -284,7 +284,6 @@
         alglib.minlmsetcond(state, 0.0, 0.0, 0.0, maxParameterOptIterations);
         alglib.minlmsetgradientcheck(state, 1.0e-3);
-        //TODO: create a type
-        var myState = new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver };
-        alglib.minlmoptimize(state, EvaluateObjectiveVector, EvaluateObjectiveVectorAndJacobian, null, myState); 
+        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);
@@ -330,5 +329,5 @@
         nmse = targetValues.Length;
       }
-      
+
     }
 
@@ -342,31 +341,15 @@
     alglib.ndimensional_jac jac;
 
-    private static void EvaluateObjectiveVector(double[] x, double[] fi, object obj) {
-      EvaluateObjectiveVectorAndJacobian(x, fi, null, obj);
-    }
-
-    private static void EvaluateObjectiveVectorAndJacobian(double[] x, double[] fi, double[,] jac, object obj) {
-      var trees = (ISymbolicExpressionTree[])((object[])obj)[0];
-      var targetVariables = (string[])((object[])obj)[1];
-      var problemData = (IRegressionProblemData)((object[])obj)[2];
-      var targetValues = (double[,])((object[])obj)[3];
-      var episodes = (IntRange[])((object[])obj)[4];
-      var numericIntegrationSteps = (int)((object[])obj)[5];
-      var latentVariables = (string[])((object[])obj)[6];
-      var odeSolver = (string)((object[])obj)[7];
+    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);
+    }
+
+    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
-      // TODO: stop integration early for diverging solutions
-      predicted = Integrate(
-          trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
-          problemData.Dataset,
-          problemData.AllowedInputVariables.ToArray(),
-          targetVariables,
-          latentVariables,
-          episodes,
-          x,
-          odeSolver,
-          numericIntegrationSteps).ToArray();
+      predicted = Integrate(optimizationData, x).ToArray();
 
       // clear all result data structures
@@ -376,5 +359,5 @@
       }
 
-      if (predicted.Length != targetValues.GetLength(0)) {
+      if (predicted.Length != optimizationData.targetValues.GetLength(0)) {
         return;
       }
@@ -398,8 +381,8 @@
       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 < targetVariables.Length; c++) {
+        for (int c = 0; c < optimizationData.targetVariables.Length; c++) {
 
           var y_pred_f = y_pred[c].Item1;
-          var y = targetValues[r, c];
+          var y = optimizationData.targetValues[r, c];
 
           fi[i] = (y - y_pred_f);
@@ -451,14 +434,6 @@
           var episodes = new[] { episode };
           var optTheta = ((DoubleArray)bestIndividualAndQuality.Item1["OptTheta_" + eIdx]).ToArray(); // see evaluate
-          var trainingPrediction = Integrate(
-                                   trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
-                                   problemData.Dataset,
-                                   problemData.AllowedInputVariables.ToArray(),
-                                   targetVars,
-                                   latentVariables,
-                                   episodes,
-                                   optTheta,
-                                   OdeSolver,
-                                   NumericIntegrationSteps).ToArray();
+          var optimizationData = new OptimizationData(trees, targetVars, problemData, null, episodes, NumericIntegrationSteps, latentVariables, OdeSolver);
+          var trainingPrediction = Integrate(optimizationData, optTheta).ToArray();
           trainingPredictions.Add(trainingPrediction);
           eIdx++;
@@ -492,14 +467,6 @@
       } else {
         var optTheta = ((DoubleArray)bestIndividualAndQuality.Item1["OptTheta"]).ToArray(); // see evaluate
-        var trainingPrediction = Integrate(
-                                   trees,  // we assume trees contain expressions for the change of each target variable over time dy/dt
-                                   problemData.Dataset,
-                                   problemData.AllowedInputVariables.ToArray(),
-                                   targetVars,
-                                   latentVariables,
-                                   TrainingEpisodes,
-                                   optTheta,
-                                   OdeSolver,
-                                   NumericIntegrationSteps).ToArray();
+        var optimizationData = new OptimizationData(trees, targetVars, problemData, null, TrainingEpisodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver);
+        var trainingPrediction = Integrate(optimizationData, optTheta).ToArray();
         // for target values and latent variables
         var trainingRows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start));
@@ -564,14 +531,6 @@
         var testList = new ItemList<DataTable>();
         var testRows = ProblemData.TestIndices.ToArray();
-        var testPrediction = Integrate(
-         trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
-         problemData.Dataset,
-         problemData.AllowedInputVariables.ToArray(),
-         targetVars,
-         latentVariables,
-         new IntRange[] { ProblemData.TestPartition },
-         optTheta,
-         OdeSolver,
-         NumericIntegrationSteps).ToArray();
+        var testOptimizationData = new OptimizationData(trees, targetVars, problemData, null, new IntRange[] { ProblemData.TestPartition }, NumericIntegrationSteps, latentVariables, OdeSolver);
+        var testPrediction = Integrate(testOptimizationData, optTheta).ToArray();
 
         for (int colIdx = 0; colIdx < trees.Length; colIdx++) {
@@ -686,13 +645,19 @@
     // for a solver with the necessary features see: https://computation.llnl.gov/projects/sundials/cvodes
 
-    public static IEnumerable<Tuple<double, Vector>[]> Integrate(
-      ISymbolicExpressionTree[] trees, IDataset dataset,
-      string[] inputVariables, string[] targetVariables, string[] latentVariables, IEnumerable<IntRange> episodes,
-      double[] parameterValues,
-      string odeSolver, int numericIntegrationSteps = 100) {
+    public static IEnumerable<Tuple<double, Vector>[]> Integrate(OptimizationData optimizationData, double[] parameterValues) {
+
+      var trees = optimizationData.trees;
+      var dataset = optimizationData.problemData.Dataset;
+      var inputVariables = optimizationData.problemData.AllowedInputVariables.ToArray();
+      var targetVariables = optimizationData.targetVariables;
+      var latentVariables = optimizationData.latentVariables;
+      var episodes = optimizationData.episodes;
+      var odeSolver = optimizationData.odeSolver;
+      var numericIntegrationSteps = optimizationData.numericIntegrationSteps;
+
 
       // TODO: numericIntegrationSteps is only relevant for the HeuristicLab solver
       var episodeIdx = 0;
-      foreach (var episode in episodes) {
+      foreach (var episode in optimizationData.episodes) {
         var rows = Enumerable.Range(episode.Start, episode.End - episode.Start);
 
@@ -1129,22 +1094,22 @@
       ISymbolicExpressionTreeNode node,
       Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>> nodeValues,      // contains value and gradient vector for a node (variables and constants only)
-      out double f,
-      out Vector g
+      out double z,
+      out Vector dz
       ) {
-      double fl, fr;
-      Vector gl, gr;
+      double f, g;
+      Vector df, dg;
       switch (node.Symbol.Name) {
         case "+": {
-            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);
+            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);
             break;
           }
         case "*": {
-            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'
+            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'
             break;
           }
@@ -1152,52 +1117,52 @@
         case "-": {
             if (node.SubtreeCount == 1) {
-              InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
-              f = -fl;
-              g = gl.Scale(-1.0);
+              InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
+              z = -f;
+              dz = -df;
             } else {
-              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);
+              InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
+              InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
+
+              z = f - g;
+              dz = df - dg;
             }
             break;
           }
         case "%": {
-            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
-            InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
+            InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
+            InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
 
             // protected division
-            if (fr.IsAlmost(0.0)) {
-              f = 0;
-              g = Vector.Zero;
+            if (g.IsAlmost(0.0)) {
+              z = 0;
+              dz = Vector.Zero;
             } else {
-              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'
+              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'
             }
             break;
           }
         case "sin": {
-            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
-            f = Math.Sin(fl);
-            g = gl.Scale(Math.Cos(fl));
+            InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
+            z = Math.Sin(f);
+            dz = Math.Cos(f) * df;
             break;
           }
         case "cos": {
-            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
-            f = Math.Cos(fl);
-            g = gl.Scale(-Math.Sin(fl));
+            InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
+            z = Math.Cos(f);
+            dz = -Math.Sin(f) * df;
             break;
           }
         case "sqr": {
-            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
-            f = fl * fl;
-            g = gl.Scale(2.0 * fl);
+            InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
+            z = f * f;
+            dz = 2.0 * f * df;
             break;
           }
         default: {
             var t = nodeValues[node];
-            f = t.Item1;
-            g = Vector.CreateNew(t.Item2);
+            z = t.Item1;
+            dz = Vector.CreateNew(t.Item2);
             break;
           }
@@ -1455,4 +1420,26 @@
       return ProblemData;
     }
+
+    public class OptimizationData {
+      public readonly ISymbolicExpressionTree[] trees;
+      public readonly string[] targetVariables;
+      public readonly IRegressionProblemData problemData;
+      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) {
+        this.trees = trees;
+        this.targetVariables = targetVars;
+        this.problemData = problemData;
+        this.targetValues = targetValues;
+        this.episodes = episodes;
+        this.numericIntegrationSteps = numericIntegrationSteps;
+        this.latentVariables = latentVariables;
+        this.odeSolver = odeSolver;
+      }
+    }
     #endregion
 
Index: /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Solution.cs
===================================================================
--- /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Solution.cs	(revision 16599)
+++ /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Solution.cs	(revision 16600)
@@ -91,14 +91,6 @@
       var random = new FastRandom(12345);
       Problem.OptimizeForEpisodes(trees, problemData, targetVars, latentVariables, random, new[] { forecastEpisode }, 100, numericIntegrationSteps, odeSolver, out optL0, out snmse);
-      var predictions = Problem.Integrate(
-         trees,
-         problemData.Dataset,
-         problemData.AllowedInputVariables.ToArray(),
-         targetVars,
-         latentVariables,
-         new[] { forecastEpisode },
-         optL0,
-         odeSolver,
-         numericIntegrationSteps).ToArray();
+      var optimizationData = new Problem.OptimizationData(trees, targetVars, problemData, null, new[] { forecastEpisode }, numericIntegrationSteps, latentVariables, odeSolver);
+      var predictions = Problem.Integrate(optimizationData, optL0).ToArray();
       return predictions.Select(p => p.Select(pi => pi.Item1).ToArray()).ToArray();
     }
Index: /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Vector.cs
===================================================================
--- /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Vector.cs	(revision 16599)
+++ /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Vector.cs	(revision 16600)
@@ -16,17 +16,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 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;
+    //   }
+    // }
 
 
@@ -38,42 +38,30 @@
       }
       Trace.Assert(a.arr.Length == b.arr.Length);
+      var res = new double[a.arr.Length];
       for (int i = 0; i < a.arr.Length; i++)
-        a.arr[i] -= b.arr[i];
-      return a;
+        res[i] = a.arr[i] - b.arr[i];
+      return new Vector(res);
     }
 
     public static Vector operator -(Vector a, Vector b) {
-      if (b == Zero) return a;
-      if (a == Zero) return -b;
-      Trace.Assert(a.arr.Length == b.arr.Length);
-      var res = new double[a.arr.Length];
-      for (int i = 0; i < res.Length; i++)
-        res[i] = a.arr[i] - b.arr[i];
-      return new Vector(res);
+      return Subtract(a, b);
     }
     public static Vector operator -(Vector v) {
-      if (v == Zero) return Zero;
-      for (int i = 0; i < v.arr.Length; i++)
-        v.arr[i] = -v.arr[i];
-      return v;
+      return v.Scale(-1.0);
     }
 
     public static Vector operator *(double s, Vector v) {
-      if (v == Zero) return Zero;
-      if (s == 0.0) return Zero;
-      var res = new double[v.arr.Length];
-      for (int i = 0; i < res.Length; i++)
-        res[i] = s * v.arr[i];
-      return new Vector(res);
+      return v.Scale(s);
     }
 
     public Vector Scale(double s) {
-      if (this != Zero) {
-        for (int i = 0; i < arr.Length; i++) {
-          arr[i] *= 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;
       }
+      return new Vector(res);
 
-      return this;
     }
 
@@ -81,4 +69,5 @@
       return s * v;
     }
+
     public static Vector operator *(Vector u, Vector v) {
       if (v == Zero) return Zero;
@@ -89,4 +78,5 @@
       return new Vector(res);
     }
+
     public static Vector operator /(double s, Vector v) {
       if (s == 0.0) return Zero;
@@ -94,9 +84,10 @@
       var res = new double[v.arr.Length];
       for (int i = 0; i < res.Length; i++)
-        res[i] = 1.0 / v.arr[i];
+        res[i] = s / v.arr[i];
       return new Vector(res);
     }
+
     public static Vector operator /(Vector v, double s) {
-      return v * 1.0 / s;
+      return v.Scale(1.0 / s);
     }
 
