Index: /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/OdeParameterIdentification.cs
===================================================================
--- /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/OdeParameterIdentification.cs	(revision 16602)
+++ /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/OdeParameterIdentification.cs	(revision 16603)
@@ -189,69 +189,4 @@
       problem.Analyze(new[] { ind }, new[] { quality }, Results, rand);
     }
-
-    // private ISymbolicExpressionTree Convert(ISymbolicExpressionTree tree) {
-    //   return new SymbolicExpressionTree(Convert(tree.Root));
-    // }
-
-
-    // for translation from symbolic expressions to simple symbols
-    // private static Dictionary<Type, string> sym2str = new Dictionary<Type, string>() {
-    //   {typeof(Addition), "+" },
-    //   {typeof(Subtraction), "-" },
-    //   {typeof(Multiplication), "*" },
-    //   {typeof(Sine), "sin" },
-    //   {typeof(Cosine), "cos" },
-    //   {typeof(Square), "sqr" },
-    // };
-
-    // private ISymbolicExpressionTreeNode Convert(ISymbolicExpressionTreeNode node) {
-    //   if (sym2str.ContainsKey(node.Symbol.GetType())) {
-    //     var children = node.Subtrees.Select(st => Convert(st)).ToArray();
-    //     return Make(sym2str[node.Symbol.GetType()], children);
-    //   } else if (node.Symbol is ProgramRootSymbol) {
-    //     var child = Convert(node.GetSubtree(0));
-    //     node.RemoveSubtree(0);
-    //     node.AddSubtree(child);
-    //     return node;
-    //   } else if (node.Symbol is StartSymbol) {
-    //     var child = Convert(node.GetSubtree(0));
-    //     node.RemoveSubtree(0);
-    //     node.AddSubtree(child);
-    //     return node;
-    //   } else if (node.Symbol is Division) {
-    //     var children = node.Subtrees.Select(st => Convert(st)).ToArray();
-    //     if (children.Length == 1) {
-    //       return Make("%", new[] { new SimpleSymbol("θ", 0).CreateTreeNode(), children[0] });
-    //     } else if (children.Length != 2) throw new ArgumentException("Division is not supported for multiple arguments");
-    //     else return Make("%", children);
-    //   } else if (node.Symbol is Constant) {
-    //     return new SimpleSymbol("θ", 0).CreateTreeNode();
-    //   } else if (node.Symbol is DataAnalysis.Symbolic.Variable) {
-    //     var varNode = node as VariableTreeNode;
-    //     if (!varNode.Weight.IsAlmost(1.0)) throw new ArgumentException("Variable weights are not supported");
-    //     return new SimpleSymbol(varNode.VariableName, 0).CreateTreeNode();
-    //   } else throw new ArgumentException("Unsupported symbol: " + node.Symbol.Name);
-    // }
-
-    // private ISymbolicExpressionTreeNode Make(string op, ISymbolicExpressionTreeNode[] children) {
-    //   if (children.Length == 1) {
-    //     var s = new SimpleSymbol(op, 1).CreateTreeNode();
-    //     s.AddSubtree(children.First());
-    //     return s;
-    //   } else {
-    //     var s = new SimpleSymbol(op, 2).CreateTreeNode();
-    //     var c0 = children[0];
-    //     var c1 = children[1];
-    //     s.AddSubtree(c0);
-    //     s.AddSubtree(c1);
-    //     for (int i = 2; i < children.Length; i++) {
-    //       var sn = new SimpleSymbol(op, 2).CreateTreeNode();
-    //       sn.AddSubtree(s);
-    //       sn.AddSubtree(children[i]);
-    //       s = sn;
-    //     }
-    //     return s;
-    //   }
-    // }
     #endregion
   }
Index: /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs
===================================================================
--- /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs	(revision 16602)
+++ /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs	(revision 16603)
@@ -234,4 +234,5 @@
       string odeSolver) {
 
+      // extract constants from tree
       var constantNodes = trees.Select(t => t.IterateNodesPrefix().OfType<ConstantTreeNode>().ToArray()).ToArray();
       var initialTheta = constantNodes.Select(nodes => nodes.Select(n => n.Value).ToArray()).ToArray();
@@ -373,5 +374,5 @@
     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) {
-      var rows = optimizationData.rows.ToArray();
+      var rows = optimizationData.rows;
       var problemData = optimizationData.problemData;
       var nodeValueLookup = optimizationData.nodeValueLookup;
@@ -402,10 +403,9 @@
       var problemData = optimizationData.problemData;
       var ds = problemData.Dataset;
-      var rows = optimizationData.episodes.SelectMany(ep => Enumerable.Range(ep.Start, ep.Size)).ToArray();
+      var rows = optimizationData.rows;
 
       var nodeValueLookup = optimizationData.nodeValueLookup;
       nodeValueLookup.UpdateParamValues(x);
 
-      // TODO add integration of latent variables
       int termIdx = 0;
 
@@ -427,6 +427,6 @@
 
           var y = optimizationData.targetValues[i][trainIdx];
-          fi[termIdx] = y - f;
-          if (jac != null && g != Vector.Zero) for (int j = 0; j < g.Length; j++) jac[termIdx, j] = -g[j];
+          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]; 
 
           termIdx++;
@@ -455,9 +455,13 @@
       int outputIdx = 0;
 
-      var prediction = Integrate(optimizationData, x).ToArray();
+      nodeValueLookup.UpdateParamValues(x);
+
+      var predictionEnumerator = Integrate(optimizationData).GetEnumerator();
       var trees = optimizationData.trees;
 
+      predictionEnumerator.MoveNext();
+      var currentPrediction = predictionEnumerator.Current;
       for (int trainIdx = 0; trainIdx < rows.Length; trainIdx++) {
-        var pred_i = prediction[Math.Min(trainIdx, prediction.Length - 1)]; // if we stopped earlier in the integration then just use the last generated value
+        var pred_i = currentPrediction;
 
         for (int i = 0; i < trees.Length; i++) {
@@ -466,9 +470,10 @@
           var f = pred_i[i].Item1;
           var g = pred_i[i].Item2;
-          fi[outputIdx] = y - f;
-          if (jac != null && g != Vector.Zero) for (int j = 0; j < g.Length; j++) jac[outputIdx, j] = -g[j];
+          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];
 
           outputIdx++;
         }
+        if (predictionEnumerator.MoveNext()) currentPrediction = predictionEnumerator.Current; // if we stopped early => just use the last prediction until the end
       }
     }
@@ -513,7 +518,6 @@
         foreach (var episode in TrainingEpisodes) {
           var episodes = new[] { episode };
-          var optTheta = ((DoubleArray)bestIndividualAndQuality.Item1["OptTheta_" + eIdx]).ToArray(); // see evaluate
           var optimizationData = new OptimizationData(trees, targetVars, problemData, null, episodes, NumericIntegrationSteps, latentVariables, OdeSolver);
-          var trainingPrediction = Integrate(optimizationData, optTheta).ToArray();
+          var trainingPrediction = Integrate(optimizationData).ToArray();
           trainingPredictions.Add(trainingPrediction);
           eIdx++;
@@ -546,7 +550,9 @@
         results["Models"].Value = models;
       } else {
-        var optTheta = Problem.ExtractParametersFromTrees(trees);
         var optimizationData = new OptimizationData(trees, targetVars, problemData, null, TrainingEpisodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver);
-        var trainingPrediction = Integrate(optimizationData, optTheta).ToArray();
+        var trainingPrediction = Integrate(optimizationData).ToArray();
+
+
+        var numParams = optimizationData.nodeValueLookup.ParameterCount;
         // for target values and latent variables
         var trainingRows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start));
@@ -561,5 +567,5 @@
             trainingDataTable.Rows.Add(predictedValuesRow);
 
-            for (int paramIdx = 0; paramIdx < optTheta.Length; paramIdx++) {
+            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());
               paramSensitivityRow.VisualProperties.SecondYAxis = true;
@@ -580,5 +586,5 @@
         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();
+        var gradientRows = Enumerable.Range(0, numParams).Select(i => new DataRow($"∂SE/∂θ{i}")).ToArray();
         errorTable.Rows.Add(seRow);
         foreach (var gRow in gradientRows) {
@@ -592,5 +598,5 @@
           // calculate objective function gradient
           double f_i = 0.0;
-          Vector g_i = Vector.CreateNew(new double[optTheta.Length]);
+          Vector g_i = Vector.CreateNew(new double[numParams]);
           for (int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
             var y_pred_f = y_pred[colIdx].Item1;
@@ -612,5 +618,5 @@
         var testRows = ProblemData.TestIndices.ToArray();
         var testOptimizationData = new OptimizationData(trees, targetVars, problemData, null, new IntRange[] { ProblemData.TestPartition }, NumericIntegrationSteps, latentVariables, OdeSolver);
-        var testPrediction = Integrate(testOptimizationData, optTheta).ToArray();
+        var testPrediction = Integrate(testOptimizationData).ToArray();
 
         for (int colIdx = 0; colIdx < trees.Length; colIdx++) {
@@ -644,21 +650,11 @@
         var models = new VariableCollection();    // to store target var names and original version of tree
 
-        var optimizedTrees = new List<ISymbolicExpressionTree>();
-        int nextParIdx = 0;
+        var clonedTrees = new List<ISymbolicExpressionTree>();
         for (int idx = 0; idx < trees.Length; idx++) {
-          var tree = trees[idx];
-          // optimizedTrees.Add(new SymbolicExpressionTree(FixParameters(tree.Root, optTheta.ToArray(), ref nextParIdx)));
-          optimizedTrees.Add(tree);
+          clonedTrees.Add((ISymbolicExpressionTree)trees[idx]);
         }
         var ds = problemData.Dataset;
-        var newVarNames = Enumerable.Range(0, nextParIdx).Select(i => "c_" + i).ToArray();
-        var allVarNames = ds.DoubleVariables.Concat(newVarNames);
-        var newVarValues = Enumerable.Range(0, nextParIdx).Select(i => "c_" + i).ToArray();
-        var allVarValues = ds.DoubleVariables.Select(varName => ds.GetDoubleValues(varName).ToList())
-          .Concat(Enumerable.Range(0, nextParIdx).Select(i => Enumerable.Repeat(optTheta[i], ds.Rows).ToList()))
-          .ToList();
-        var newDs = new Dataset(allVarNames, allVarValues);
-        var newProblemData = new RegressionProblemData(newDs, problemData.AllowedInputVariables.Concat(newVarValues).ToArray(), problemData.TargetVariable);
-        results["Solution"].Value = new Solution(optimizedTrees.ToArray(),
+        var newProblemData = new RegressionProblemData((IDataset)ds.Clone(), problemData.AllowedInputVariables, problemData.TargetVariable);
+        results["Solution"].Value = new Solution(clonedTrees.ToArray(),
                    // optTheta,
                    newProblemData,
@@ -670,5 +666,4 @@
 
 
-        nextParIdx = 0;
         for (int idx = 0; idx < trees.Length; idx++) {
           var varName = string.Empty;
@@ -715,5 +710,5 @@
     // for a solver with the necessary features see: https://computation.llnl.gov/projects/sundials/cvodes
 
-    public static IEnumerable<Tuple<double, Vector>[]> Integrate(OptimizationData optimizationData, double[] parameterValues) {
+    public static IEnumerable<Tuple<double, Vector>[]> Integrate(OptimizationData optimizationData) {
 
       var trees = optimizationData.trees;
@@ -727,8 +722,5 @@
       var calculatedVariables = targetVariables.Concat(latentVariables).ToArray(); // TODO: must conincide with the order of trees in the encoding
 
-
       var nodeValues = optimizationData.nodeValueLookup;
-      nodeValues.UpdateParamValues(parameterValues);
-
 
       // TODO: numericIntegrationSteps is only relevant for the HeuristicLab solver
@@ -744,18 +736,18 @@
         }
 
-        { // CODE BELOW MUST BE TESTED
-          if (latentVariables.Length > 0) throw new NotImplementedException();
-
-          // add value entries for latent variables which are also integrated
-          // initial values are at the end of the parameter vector
-          // separate initial values for each episode
-          var initialValueIdx = parameterValues.Length - episodes.Count() * latentVariables.Length + episodeIdx * latentVariables.Length;
-          foreach (var latentVar in latentVariables) {
-            var arr = new double[parameterValues.Length]; // backing array
-            arr[initialValueIdx] = 1.0;
-            var g = new Vector(arr);
-            nodeValues.SetVariableValue(latentVar, parameterValues[initialValueIdx], g); // we don't have observations for latent variables therefore we optimize the initial value for each episode 
-            initialValueIdx++;
-          }
+        { // CODE BELOW DOESN'T WORK ANYMORE
+          // if (latentVariables.Length > 0) throw new NotImplementedException();
+          // 
+          // // add value entries for latent variables which are also integrated
+          // // initial values are at the end of the parameter vector
+          // // separate initial values for each episode
+          // var initialValueIdx = parameterValues.Length - episodes.Count() * latentVariables.Length + episodeIdx * latentVariables.Length;
+          // foreach (var latentVar in latentVariables) {
+          //   var arr = new double[parameterValues.Length]; // backing array
+          //   arr[initialValueIdx] = 1.0;
+          //   var g = new Vector(arr);
+          //   nodeValues.SetVariableValue(latentVar, parameterValues[initialValueIdx], g); // we don't have observations for latent variables therefore we optimize the initial value for each episode 
+          //   initialValueIdx++;
+          // }
         }
 
@@ -1129,5 +1121,7 @@
     // TODO: use an existing interpreter implementation instead
     private static double InterpretRec(ISymbolicExpressionTreeNode node, NodeValueLookup nodeValues) {
-      if (node.Symbol is Constant || node.Symbol is Variable) {
+      if (node is ConstantTreeNode) {
+        return ((ConstantTreeNode)node).Value;
+      } else if(node is VariableTreeNode) {
         return nodeValues.NodeValue(node);
       } else if (node.Symbol is Addition) {
@@ -1370,8 +1364,18 @@
       var g = CreateGrammar();
       foreach (var targetVar in TargetVariables.CheckedItems) {
-        encoding = encoding.Add(new SymbolicExpressionTreeEncoding(targetVar + "_tree", g, MaximumLength, MaximumLength)); // only limit by length
+        var e = new SymbolicExpressionTreeEncoding(targetVar + "_tree", g, MaximumLength, MaximumLength);
+        var multiManipulator = e.Operators.Where(op => op is MultiSymbolicExpressionTreeManipulator).First();
+        var filteredOperators = e.Operators.Where(op => !(op is IManipulator)).ToArray();
+        // make sure our multi-manipulator is the only manipulator
+        e.Operators = new IOperator[] { multiManipulator }.Concat(filteredOperators);
+        encoding = encoding.Add(e); // only limit by length
       }
       for (int i = 1; i <= NumberOfLatentVariables; i++) {
-        encoding = encoding.Add(new SymbolicExpressionTreeEncoding("λ" + i + "_tree", g, MaximumLength, MaximumLength));
+        var e = new SymbolicExpressionTreeEncoding("λ" + i + "_tree", g, MaximumLength, MaximumLength);
+        var multiManipulator = e.Operators.Where(op => op is MultiSymbolicExpressionTreeManipulator).First();
+        var filteredOperators = e.Operators.Where(op => !(op is IManipulator)).ToArray();
+        // make sure our multi-manipulator is the only manipulator
+        e.Operators = new IOperator[] { multiManipulator }.Concat(filteredOperators);
+        encoding = encoding.Add(e);
       }
       Encoding = encoding;
@@ -1492,4 +1496,5 @@
       public readonly IRegressionProblemData problemData;
       public readonly double[][] targetValues;
+      public readonly double[] inverseStandardDeviation;
       public readonly IntRange[] episodes;
       public readonly int numericIntegrationSteps;
@@ -1497,5 +1502,5 @@
       public readonly string odeSolver;
       public readonly NodeValueLookup nodeValueLookup;
-      public IEnumerable<int> rows => episodes.SelectMany(ep => Enumerable.Range(ep.Start, ep.Size));
+      public readonly int[] rows;
 
       public OptimizationData(ISymbolicExpressionTree[] trees, string[] targetVars, IRegressionProblemData problemData,
@@ -1508,4 +1513,8 @@
         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.episodes = episodes;
         this.numericIntegrationSteps = numericIntegrationSteps;
@@ -1513,4 +1522,5 @@
         this.odeSolver = odeSolver;
         this.nodeValueLookup = new NodeValueLookup(trees);
+        this.rows =  episodes.SelectMany(ep => Enumerable.Range(ep.Start, ep.Size)).ToArray();
       }
     }
@@ -1519,5 +1529,8 @@
       private readonly Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>> node2val = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
       private readonly Dictionary<string, List<ISymbolicExpressionTreeNode>> name2nodes = new Dictionary<string, List<ISymbolicExpressionTreeNode>>();
-      private readonly Dictionary<int, ISymbolicExpressionTreeNode> paramIdx2node = new Dictionary<int, ISymbolicExpressionTreeNode>();
+      private readonly ConstantTreeNode[] constantNodes;
+      private readonly Vector[] constantGradientVectors;
+
+      // private readonly Dictionary<int, ISymbolicExpressionTreeNode> paramIdx2node = new Dictionary<int, ISymbolicExpressionTreeNode>();
 
       public double NodeValue(ISymbolicExpressionTreeNode node) => node2val[node].Item1;
@@ -1525,11 +1538,12 @@
 
       public NodeValueLookup(ISymbolicExpressionTree[] trees) {
-        int paramIdx = 0;
-
-        var constantNodes = trees.SelectMany(t => t.IterateNodesPrefix().Where(IsConstantNode)).ToArray();
-        foreach (var n in constantNodes) {
-          paramIdx2node[paramIdx] = n;
-          SetParamValue(paramIdx, GetConstantValue(n), Vector.CreateIndicator(length: constantNodes.Length, idx: paramIdx));
-          paramIdx++;
+
+        this.constantNodes = trees.SelectMany(t => t.IterateNodesPrefix().OfType<ConstantTreeNode>()).ToArray();
+        constantGradientVectors = new Vector[constantNodes.Length];
+        for(int paramIdx = 0; paramIdx < constantNodes.Length; paramIdx++) {
+          constantGradientVectors[paramIdx] = Vector.CreateIndicator(length: constantNodes.Length, idx: paramIdx);
+
+          var node = constantNodes[paramIdx];
+          node2val[node] = Tuple.Create(node.Value, constantGradientVectors[paramIdx]);
         }
 
@@ -1547,10 +1561,5 @@
       }
 
-      public int ParameterCount => paramIdx2node.Count;
-
-
-      public void SetParamValue(int paramIdx, double val, Vector dVal) {
-        node2val[paramIdx2node[paramIdx]] = Tuple.Create(val, dVal);
-      }
+      public int ParameterCount => constantNodes.Length;
 
       public void SetVariableValue(string variableName, double val) {
@@ -1573,8 +1582,7 @@
 
       internal void UpdateParamValues(double[] x) {
-        Trace.Assert(x.Length == paramIdx2node.Count);
-        for (int paramIdx = 0; paramIdx < x.Length; paramIdx++) {
-          var n = paramIdx2node[paramIdx];
-          node2val[n] = Tuple.Create(x[paramIdx], node2val[n].Item2); // prevent allocation of new Vector
+        for (int i = 0; i < x.Length; i++) {
+          constantNodes[i].Value = x[i];
+          node2val[constantNodes[i]] = Tuple.Create(x[i], constantGradientVectors[i]);
         }
       }
Index: /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Solution.cs
===================================================================
--- /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Solution.cs	(revision 16602)
+++ /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Solution.cs	(revision 16603)
@@ -88,14 +88,14 @@
     public IEnumerable<double[]> Predict(IntRange episode, int forecastHorizon, out double snmse) {
       var forecastEpisode = new IntRange(episode.Start, episode.End + forecastHorizon);
+      // 
+      // 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 theta = Problem.ExtractParametersFromTrees(trees);
 
-      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 theta = Problem.ExtractParametersFromTrees(trees);
-
-      var predictions = Problem.Integrate(optimizationData, theta).ToArray();
-      return predictions.Select(p => p.Select(pi => pi.Item1).ToArray()).ToArray();
+      snmse = 0.0; // TODO
+      return Problem.Integrate(optimizationData).Select(p => p.Select(pi => pi.Item1).ToArray()).ToArray();
     }
   }
