Index: /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs
===================================================================
--- /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs	(revision 16154)
+++ /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs	(revision 16155)
@@ -117,4 +117,5 @@
     private const string NumericIntegrationStepsParameterName = "Steps for numeric integration";
     private const string TrainingEpisodesParameterName = "Training episodes";
+    private const string OptimizeParametersForEpisodesParameterName = "Optimize parameters for episodes";
     #endregion
 
@@ -145,4 +146,7 @@
     public IValueParameter<ItemList<IntRange>> TrainingEpisodesParameter {
       get { return (IValueParameter<ItemList<IntRange>>)Parameters[TrainingEpisodesParameterName]; }
+    }
+    public IFixedValueParameter<BoolValue> OptimizeParametersForEpisodesParameter {
+      get { return (IFixedValueParameter<BoolValue>)Parameters[OptimizeParametersForEpisodesParameterName]; }
     }
     #endregion
@@ -177,4 +181,7 @@
     public IEnumerable<IntRange> TrainingEpisodes {
       get { return TrainingEpisodesParameter.Value; }
+    }
+    public bool OptimizeParametersForEpisodes {
+      get { return OptimizeParametersForEpisodesParameter.Value.Value; }
     }
 
@@ -193,4 +200,7 @@
     [StorableHook(HookType.AfterDeserialization)]
     private void AfterDeserialization() {
+      if(!Parameters.ContainsKey(OptimizeParametersForEpisodesParameterName)) {
+        Parameters.Add(new FixedValueParameter<BoolValue>(OptimizeParametersForEpisodesParameterName, "Flag to select if parameters should be optimized globally or for each episode individually.", new BoolValue(false)));
+      }
       RegisterEventHandlers();
     }
@@ -216,5 +226,5 @@
       Parameters.Add(new FixedValueParameter<IntValue>(NumericIntegrationStepsParameterName, "Number of steps in the numeric integration that are taken from one row to the next (set to 1 to 100). More steps makes the algorithm slower, less steps worsens the accuracy of the numeric integration scheme.", new IntValue(10)));
       Parameters.Add(new ValueParameter<ItemList<IntRange>>(TrainingEpisodesParameterName, "A list of ranges that should be used for training, each range represents an independent episode. This overrides the TrainingSet parameter in ProblemData.", new ItemList<IntRange>()));
-
+      Parameters.Add(new FixedValueParameter<BoolValue>(OptimizeParametersForEpisodesParameterName, "Flag to select if parameters should be optimized globally or for each episode individually.", new BoolValue(false)));
       RegisterEventHandlers();
       InitAllParameters();
@@ -225,10 +235,33 @@
     }
 
-
     public override double Evaluate(Individual individual, IRandom random) {
       var trees = individual.Values.Select(v => v.Value).OfType<ISymbolicExpressionTree>().ToArray(); // extract all trees from individual
 
+      if(OptimizeParametersForEpisodes) {
+        int eIdx = 0;
+        double totalNMSE = 0.0;
+        int totalSize = 0;
+        foreach(var episode in TrainingEpisodes) {
+          double[] optTheta;
+          double nmse;
+          OptimizeForEpisodes(trees, random, new[] { episode }, out optTheta, out nmse);
+          individual["OptTheta_" + eIdx] = new DoubleArray(optTheta); // write back optimized parameters so that we can use them in the Analysis method
+          eIdx++;
+          totalNMSE += nmse * episode.Size;
+          totalSize += episode.Size;
+        }
+        return totalNMSE / totalSize;
+      } else {
+        double[] optTheta;
+        double nmse;
+        OptimizeForEpisodes(trees, random, TrainingEpisodes, out optTheta, out nmse);
+        individual["OptTheta"] = new DoubleArray(optTheta); // write back optimized parameters so that we can use them in the Analysis method
+        return nmse;
+      }
+    }
+
+    private void OptimizeForEpisodes(ISymbolicExpressionTree[] trees, IRandom random, IEnumerable<IntRange> episodes, out double[] optTheta, out double nmse) {
+      var rows = episodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start)).ToArray();
       var problemData = ProblemData;
-      var rows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start)).ToArray();
       var targetVars = TargetVariables.CheckedItems.Select(i => i.Value).ToArray();
       var latentVariables = Enumerable.Range(1, NumberOfLatentVariables).Select(i => "λ" + i).ToArray(); // TODO: must coincide with the variables which are actually defined in the grammar and also for which we actually have trees
@@ -256,5 +289,5 @@
       var theta = nodeIdx.Select(_ => random.NextDouble() * 2.0 - 1.0).ToArray(); // init params randomly from Unif(-1,1)
 
-      double[] optTheta = new double[0];
+      optTheta = new double[0];
       if(theta.Length > 0) {
         alglib.minlbfgsstate state;
@@ -263,5 +296,5 @@
         alglib.minlbfgssetcond(state, 0.0, 0.0, 0.0, MaximumParameterOptimizationIterations);
         alglib.minlbfgsoptimize(state, EvaluateObjectiveAndGradient, null,
-          new object[] { trees, targetVars, problemData, nodeIdx, targetValues, TrainingEpisodes.ToArray(), NumericIntegrationSteps, latentVariables }); //TODO: create a type
+          new object[] { trees, targetVars, problemData, nodeIdx, targetValues, episodes.ToArray(), NumericIntegrationSteps, latentVariables }); //TODO: create a type
         alglib.minlbfgsresults(state, out optTheta, out report);
 
@@ -292,16 +325,13 @@
                           * NFEV countains number of function calculations
          */
-        if(report.terminationtype < 0) return double.MaxValue;
+        if(report.terminationtype < 0) { nmse = 10E6; return; }
       }
 
       // perform evaluation for optimal theta to get quality value
       double[] grad = new double[optTheta.Length];
-      double optQuality = double.NaN;
-      EvaluateObjectiveAndGradient(optTheta, ref optQuality, grad,
-        new object[] { trees, targetVars, problemData, nodeIdx, targetValues, TrainingEpisodes.ToArray(), NumericIntegrationSteps, latentVariables });
-      if(double.IsNaN(optQuality) || double.IsInfinity(optQuality)) return 10E6; // return a large value (TODO: be consistent by using NMSE)
-
-      individual["OptTheta"] = new DoubleArray(optTheta); // write back optimized parameters so that we can use them in the Analysis method
-      return optQuality;
+      nmse = double.NaN;
+      EvaluateObjectiveAndGradient(optTheta, ref nmse, grad,
+        new object[] { trees, targetVars, problemData, nodeIdx, targetValues, episodes.ToArray(), NumericIntegrationSteps, latentVariables });
+      if(double.IsNaN(nmse) || double.IsInfinity(nmse)) { nmse = 10E6; return; } // return a large value (TODO: be consistent by using NMSE)
     }
 
@@ -371,11 +401,9 @@
       }
 
+      var bestIndividualAndQuality = this.GetBestIndividual(individuals, qualities);
+      var trees = bestIndividualAndQuality.Item1.Values.Select(v => v.Value).OfType<ISymbolicExpressionTree>().ToArray(); // extract all trees from individual
+
       // TODO extract common functionality from Evaluate and Analyze
-      var bestIndividualAndQuality = this.GetBestIndividual(individuals, qualities);
-      var optTheta = ((DoubleArray)bestIndividualAndQuality.Item1["OptTheta"]).ToArray(); // see evaluate
-      var trees = bestIndividualAndQuality.Item1.Values.Select(v => v.Value).OfType<ISymbolicExpressionTree>().ToArray(); // extract all trees from individual
       var nodeIdx = new Dictionary<ISymbolicExpressionTreeNode, int>();
-
-
       foreach(var tree in trees) {
         foreach(var node in tree.Root.IterateNodesPrefix().Where(n => IsConstantNode(n))) {
@@ -388,91 +416,137 @@
 
       var trainingList = new ItemList<DataTable>();
-      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,
-       TrainingEpisodes,
-       nodeIdx,
-       optTheta,
-       NumericIntegrationSteps).ToArray();
-
-      // only for actual target values
-      var trainingRows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start));
-      for(int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
-        var targetVar = targetVars[colIdx];
-        var trainingDataTable = new DataTable(targetVar + " prediction (training)");
-        var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, trainingRows));
-        var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, trainingPrediction.Select(arr => arr[colIdx].Item1).ToArray());
-        trainingDataTable.Rows.Add(actualValuesRow);
-        trainingDataTable.Rows.Add(predictedValuesRow);
-        trainingList.Add(trainingDataTable);
-      }
-
-      // TODO: DRY for training and test
-      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 },
-       nodeIdx,
-       optTheta,
-       NumericIntegrationSteps).ToArray();
-
-      for(int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
-        var targetVar = targetVars[colIdx];
-        var testDataTable = new DataTable(targetVar + " prediction (test)");
-        var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, testRows));
-        var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, testPrediction.Select(arr => arr[colIdx].Item1).ToArray());
-        testDataTable.Rows.Add(actualValuesRow);
-        testDataTable.Rows.Add(predictedValuesRow);
-        testList.Add(testDataTable);
-      }
-
-      results["Prediction (training)"].Value = trainingList.AsReadOnly();
-      results["Prediction (test)"].Value = testList.AsReadOnly();
-
-      #region simplification of models
-      // TODO the dependency of HeuristicLab.Problems.DataAnalysis.Symbolic is not ideal
-      var models = new VariableCollection();    // to store target var names and original version of tree
-
-      foreach(var tup in targetVars.Zip(trees, Tuple.Create)) {
-        var targetVarName = tup.Item1;
-        var tree = tup.Item2;
-
-        // when we reference HeuristicLab.Problems.DataAnalysis.Symbolic we can translate symbols
-        int nextParIdx = 0;
-        var shownTree = new SymbolicExpressionTree(TranslateTreeNode(tree.Root, optTheta, ref nextParIdx));
-
-        // var shownTree = (SymbolicExpressionTree)tree.Clone();
-        // var constantsNodeOrig = tree.IterateNodesPrefix().Where(IsConstantNode);
-        // var constantsNodeShown = shownTree.IterateNodesPrefix().Where(IsConstantNode);
-        // 
-        // foreach (var n in constantsNodeOrig.Zip(constantsNodeShown, (original, shown) => new { original, shown })) {
-        //   double constantsVal = optTheta[nodeIdx[n.original]];
-        // 
-        //   ConstantTreeNode replacementNode = new ConstantTreeNode(new Constant()) { Value = constantsVal };
-        // 
-        //   var parentNode = n.shown.Parent;
-        //   int replacementIndex = parentNode.IndexOfSubtree(n.shown);
-        //   parentNode.RemoveSubtree(replacementIndex);
-        //   parentNode.InsertSubtree(replacementIndex, replacementNode);
-        // }
-
-        var origTreeVar = new HeuristicLab.Core.Variable(targetVarName + "(original)");
-        origTreeVar.Value = (ISymbolicExpressionTree)tree.Clone();
-        models.Add(origTreeVar);
-        var simplifiedTreeVar = new HeuristicLab.Core.Variable(targetVarName + "(simplified)");
-        simplifiedTreeVar.Value = TreeSimplifier.Simplify(shownTree);
-        models.Add(simplifiedTreeVar);
-
-      }
-      results["Models"].Value = models;
-      #endregion
+
+      if(OptimizeParametersForEpisodes) {
+        var eIdx = 0;
+        var trainingPredictions = new List<Tuple<double, Vector>[][]>();
+        foreach(var episode in TrainingEpisodes) {
+          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,
+                                   nodeIdx,
+                                   optTheta,
+                                   NumericIntegrationSteps).ToArray();
+          trainingPredictions.Add(trainingPrediction);
+          eIdx++;
+        }
+
+        // only for actual target values
+        var trainingRows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start));
+        for(int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
+          var targetVar = targetVars[colIdx];
+          var trainingDataTable = new DataTable(targetVar + " prediction (training)");
+          var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, trainingRows));
+          var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, trainingPredictions.SelectMany(arr => arr.Select(row => row[colIdx].Item1)).ToArray());
+          trainingDataTable.Rows.Add(actualValuesRow);
+          trainingDataTable.Rows.Add(predictedValuesRow);
+          trainingList.Add(trainingDataTable);
+        }
+        results["Prediction (training)"].Value = trainingList.AsReadOnly();
+
+
+        var models = new VariableCollection();
+
+        foreach(var tup in targetVars.Zip(trees, Tuple.Create)) {
+          var targetVarName = tup.Item1;
+          var tree = tup.Item2;
+
+          var origTreeVar = new HeuristicLab.Core.Variable(targetVarName + "(original)");
+          origTreeVar.Value = (ISymbolicExpressionTree)tree.Clone();
+          models.Add(origTreeVar);
+        }
+        results["Models"].Value = models;
+      } 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 y'(t)
+                                   problemData.Dataset,
+                                   problemData.AllowedInputVariables.ToArray(),
+                                   targetVars,
+                                   latentVariables,
+                                   TrainingEpisodes,
+                                   nodeIdx,
+                                   optTheta,
+                                   NumericIntegrationSteps).ToArray();
+        // only for actual target values
+        var trainingRows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start));
+        for(int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
+          var targetVar = targetVars[colIdx];
+          var trainingDataTable = new DataTable(targetVar + " prediction (training)");
+          var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, trainingRows));
+          var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, trainingPrediction.Select(arr => arr[colIdx].Item1).ToArray());
+          trainingDataTable.Rows.Add(actualValuesRow);
+          trainingDataTable.Rows.Add(predictedValuesRow);
+          trainingList.Add(trainingDataTable);
+        }
+        // TODO: DRY for training and test
+        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 },
+         nodeIdx,
+         optTheta,
+         NumericIntegrationSteps).ToArray();
+
+        for(int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
+          var targetVar = targetVars[colIdx];
+          var testDataTable = new DataTable(targetVar + " prediction (test)");
+          var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, testRows));
+          var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, testPrediction.Select(arr => arr[colIdx].Item1).ToArray());
+          testDataTable.Rows.Add(actualValuesRow);
+          testDataTable.Rows.Add(predictedValuesRow);
+          testList.Add(testDataTable);
+        }
+
+        results["Prediction (training)"].Value = trainingList.AsReadOnly();
+        results["Prediction (test)"].Value = testList.AsReadOnly();
+        #region simplification of models
+        // TODO the dependency of HeuristicLab.Problems.DataAnalysis.Symbolic is not ideal
+        var models = new VariableCollection();    // to store target var names and original version of tree
+
+        foreach(var tup in targetVars.Zip(trees, Tuple.Create)) {
+          var targetVarName = tup.Item1;
+          var tree = tup.Item2;
+
+          // when we reference HeuristicLab.Problems.DataAnalysis.Symbolic we can translate symbols
+          int nextParIdx = 0;
+          var shownTree = new SymbolicExpressionTree(TranslateTreeNode(tree.Root, optTheta, ref nextParIdx));
+
+          // var shownTree = (SymbolicExpressionTree)tree.Clone();
+          // var constantsNodeOrig = tree.IterateNodesPrefix().Where(IsConstantNode);
+          // var constantsNodeShown = shownTree.IterateNodesPrefix().Where(IsConstantNode);
+          // 
+          // foreach (var n in constantsNodeOrig.Zip(constantsNodeShown, (original, shown) => new { original, shown })) {
+          //   double constantsVal = optTheta[nodeIdx[n.original]];
+          // 
+          //   ConstantTreeNode replacementNode = new ConstantTreeNode(new Constant()) { Value = constantsVal };
+          // 
+          //   var parentNode = n.shown.Parent;
+          //   int replacementIndex = parentNode.IndexOfSubtree(n.shown);
+          //   parentNode.RemoveSubtree(replacementIndex);
+          //   parentNode.InsertSubtree(replacementIndex, replacementNode);
+          // }
+
+          var origTreeVar = new HeuristicLab.Core.Variable(targetVarName + "(original)");
+          origTreeVar.Value = (ISymbolicExpressionTree)tree.Clone();
+          models.Add(origTreeVar);
+          var simplifiedTreeVar = new HeuristicLab.Core.Variable(targetVarName + "(simplified)");
+          simplifiedTreeVar.Value = TreeSimplifier.Simplify(shownTree);
+          models.Add(simplifiedTreeVar);
+
+        }
+        results["Models"].Value = models;
+        #endregion
+      }
     }
 
