Index: /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs
===================================================================
--- /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs	(revision 15969)
+++ /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs	(revision 15970)
@@ -112,4 +112,6 @@
     private const string MaximumLengthParameterName = "Size limit";
     private const string MaximumParameterOptimizationIterationsParameterName = "Max. parameter optimization iterations";
+    private const string NumberOfLatentVariablesParameterName = "Number of latent variables";
+    private const string NumericIntegrationStepsParameterName = "Steps for numeric integration";
     #endregion
 
@@ -131,4 +133,10 @@
     public IFixedValueParameter<IntValue> MaximumParameterOptimizationIterationsParameter {
       get { return (IFixedValueParameter<IntValue>)Parameters[MaximumParameterOptimizationIterationsParameterName]; }
+    }
+    public IFixedValueParameter<IntValue> NumberOfLatentVariablesParameter {
+      get { return (IFixedValueParameter<IntValue>)Parameters[NumberOfLatentVariablesParameterName]; }
+    }
+    public IFixedValueParameter<IntValue> NumericIntegrationStepsParameter {
+      get { return (IFixedValueParameter<IntValue>)Parameters[NumericIntegrationStepsParameterName]; }
     }
     #endregion
@@ -155,4 +163,11 @@
       get { return MaximumParameterOptimizationIterationsParameter.Value.Value; }
     }
+    public int NumberOfLatentVariables {
+      get { return NumberOfLatentVariablesParameter.Value.Value; }
+    }
+    public int NumericIntegrationSteps {
+      get { return NumericIntegrationStepsParameter.Value.Value; }
+    }
+
     #endregion                                                                                      
 
@@ -184,9 +199,11 @@
       var targetVariables = new CheckedItemCollection<StringValue>().AsReadOnly(); // HACK: it would be better to provide a new class derived from IDataAnalysisProblem
       var functions = CreateFunctionSet();
-      Parameters.Add(new ValueParameter<IRegressionProblemData>(ProblemDataParameterName, "The data captured from the dynamical system", new RegressionProblemData()));
+      Parameters.Add(new ValueParameter<IRegressionProblemData>(ProblemDataParameterName, "The data captured from the dynamical system. Use CSV import functionality to import data.", new RegressionProblemData()));
       Parameters.Add(new ValueParameter<ReadOnlyCheckedItemCollection<StringValue>>(TargetVariablesParameterName, "Target variables (overrides setting in ProblemData)", targetVariables));
       Parameters.Add(new ValueParameter<ReadOnlyCheckedItemCollection<StringValue>>(FunctionSetParameterName, "The list of allowed functions", functions));
-      Parameters.Add(new FixedValueParameter<IntValue>(MaximumLengthParameterName, "The maximally allowed length of each expression", new IntValue(20)));
-      Parameters.Add(new FixedValueParameter<IntValue>(MaximumParameterOptimizationIterationsParameterName, "The maximum number of iterations for optimization of parameters (using L-BFGS)", new IntValue(100)));
+      Parameters.Add(new FixedValueParameter<IntValue>(MaximumLengthParameterName, "The maximally allowed length of each expression. Set to a small value (5 - 25). Default = 10", new IntValue(10)));
+      Parameters.Add(new FixedValueParameter<IntValue>(MaximumParameterOptimizationIterationsParameterName, "The maximum number of iterations for optimization of parameters (using L-BFGS). More iterations makes the algorithm slower, fewer iterations might prevent convergence in the optimization scheme. Default = 100", new IntValue(100)));
+      Parameters.Add(new FixedValueParameter<IntValue>(NumberOfLatentVariablesParameterName, "Latent variables (unobserved variables) allow us to produce expressions which are integrated up and can be used in other expressions. They are handled similarly to target variables in forward simulation / integration. The difference to target variables is that there are no data to which the calculated values of latent variables are compared. Set to a small value (0 .. 5) as necessary (default = 0)", new IntValue(0)));
+      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)));
 
       RegisterEventHandlers();
@@ -201,11 +218,12 @@
       var rows = ProblemData.TrainingIndices.ToArray();
       var targetVars = TargetVariables.CheckedItems.Select(i => i.Value).ToArray();
-      var targetValues = new double[rows.Length,targetVars.Length];
-      
+      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
+      var targetValues = new double[rows.Length, targetVars.Length];
+
       // collect values of all target variables
       var colIdx = 0;
-      foreach(var targetVar in targetVars) {
+      foreach (var targetVar in targetVars) {
         int rowIdx = 0;
-        foreach(var value in problemData.Dataset.GetDoubleValues(targetVar, rows)) {
+        foreach (var value in problemData.Dataset.GetDoubleValues(targetVar, rows)) {
           targetValues[rowIdx, colIdx] = value;
           rowIdx++;
@@ -230,5 +248,6 @@
         alglib.minlbfgscreate(Math.Min(theta.Length, 5), theta, out state);
         alglib.minlbfgssetcond(state, 0.0, 0.0, 0.0, MaximumParameterOptimizationIterations);
-        alglib.minlbfgsoptimize(state, EvaluateObjectiveAndGradient, null, new object[] { trees, targetVars, problemData, nodeIdx, targetValues, rows }); //TODO: create a type
+        alglib.minlbfgsoptimize(state, EvaluateObjectiveAndGradient, null, 
+          new object[] { trees, targetVars, problemData, nodeIdx, targetValues, rows, NumericIntegrationSteps, latentVariables }); //TODO: create a type
         alglib.minlbfgsresults(state, out optTheta, out report);
 
@@ -265,5 +284,6 @@
       double[] grad = new double[optTheta.Length];
       double optQuality = double.NaN;
-      EvaluateObjectiveAndGradient(optTheta, ref optQuality, grad, new object[] { trees, targetVars, problemData, nodeIdx, targetValues, rows });
+      EvaluateObjectiveAndGradient(optTheta, ref optQuality, grad, 
+        new object[] { trees, targetVars, problemData, nodeIdx, targetValues, rows, NumericIntegrationSteps, latentVariables });
       if (double.IsNaN(optQuality) || double.IsInfinity(optQuality)) return 10E6; // return a large value (TODO: be consistent by using NMSE)
 
@@ -279,4 +299,6 @@
       var targetValues = (double[,])((object[])obj)[4];
       var rows = (int[])((object[])obj)[5];
+      var numericIntegrationSteps = (int)((object[])obj)[6];
+      var latentVariables = (string[])((object[])obj)[7];
 
       var predicted = Integrate(
@@ -285,7 +307,8 @@
         problemData.AllowedInputVariables.ToArray(),
         targetVariables,
+        latentVariables,
         rows,
         nodeIdx,                // TODO: is it Ok to use rows here ?
-        x).ToArray();
+        x, numericIntegrationSteps).ToArray();
 
 
@@ -305,8 +328,8 @@
       foreach (var y_pred in predicted) {
         // TODO NMSE to put the same weight on each target regardless of the value range;
-        for(int c = 0;c<y_pred.Length;c++) {
-          
+        for (int c = 0; c < y_pred.Length; c++) {
+
           var y_pred_f = y_pred[c].Item1;
-          var y = targetValues[r,c];
+          var y = targetValues[r, c];
 
           var res = (y - y_pred_f);
@@ -336,5 +359,5 @@
       // TODO extract common functionality from Evaluate and Analyze
       var bestIndividualAndQuality = this.GetBestIndividual(individuals, qualities);
-      var optTheta = ((DoubleArray) bestIndividualAndQuality.Item1["OptTheta"]).ToArray(); // see evaluate
+      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>();
@@ -348,4 +371,5 @@
       var problemData = ProblemData;
       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
 
       var trainingList = new ItemList<DataTable>();
@@ -356,11 +380,13 @@
        problemData.AllowedInputVariables.ToArray(),
        targetVars,
+       latentVariables,
        trainingRows,
        nodeIdx,
-       optTheta).ToArray();
+       optTheta, 
+       NumericIntegrationSteps).ToArray();
 
       for (int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
         var targetVar = targetVars[colIdx];
-        var trainingDataTable = new DataTable(targetVar+ " prediction (training)");
+        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());
@@ -378,7 +404,9 @@
        problemData.AllowedInputVariables.ToArray(),
        targetVars,
+       latentVariables,
        testRows,
        nodeIdx,
-       optTheta).ToArray();
+       optTheta, 
+       NumericIntegrationSteps).ToArray();
 
       for (int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
@@ -400,12 +428,11 @@
     #region interpretation
     private static IEnumerable<Tuple<double, Vector>[]> Integrate(
-      ISymbolicExpressionTree[] trees, IDataset dataset, string[] inputVariables, string[] targetVariables, IEnumerable<int> rows,
-      Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx, double[] parameterValues) {
-
-      int NUM_STEPS = 1;
+      ISymbolicExpressionTree[] trees, IDataset dataset, string[] inputVariables, string[] targetVariables, string[] latentVariables, IEnumerable<int> rows,
+      Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx, double[] parameterValues, int numericIntegrationSteps = 100) {
+
+      int NUM_STEPS = numericIntegrationSteps ;
       double h = 1.0 / NUM_STEPS;
 
       // return first value as stored in the dataset
-
       yield return targetVariables
         .Select(targetVar => Tuple.Create(dataset.GetDoubleValue(targetVar, rows.First()), Vector.Zero))
@@ -422,9 +449,14 @@
         variableValues.Add(varName, Tuple.Create(dataset.GetDoubleValue(varName, t0), Vector.Zero));
       }
+      // add value entries for latent variables which are also integrated
+      foreach(var latentVar in latentVariables) {
+        variableValues.Add(latentVar, Tuple.Create(0.0, Vector.Zero)); // we don't have observations for latent variables -> assume zero as starting value
+      }
+      var calculatedVariables = targetVariables.Concat(latentVariables); // TODO: must conincide with the order of trees in the encoding
 
       foreach (var t in rows.Skip(1)) {
         for (int step = 0; step < NUM_STEPS; step++) {
           var deltaValues = new Dictionary<string, Tuple<double, Vector>>();
-          foreach (var tup in trees.Zip(targetVariables, Tuple.Create)) {
+          foreach (var tup in trees.Zip(calculatedVariables, Tuple.Create)) {
             var tree = tup.Item1;
             var targetVarName = tup.Item2;
@@ -444,4 +476,5 @@
         }
 
+        // only return the target variables for calculation of errors
         yield return targetVariables
           .Select(targetVar => variableValues[targetVar])
@@ -464,5 +497,5 @@
       switch (node.Symbol.Name) {
         case "+": {
-            var l = InterpretRec(node.GetSubtree(0), variableValues, nodeIdx, parameterValues);
+            var l = InterpretRec(node.GetSubtree(0), variableValues, nodeIdx, parameterValues); // TODO capture all parameters into a state type for interpretation
             var r = InterpretRec(node.GetSubtree(1), variableValues, nodeIdx, parameterValues);
 
@@ -520,8 +553,8 @@
      *    |
      *    V
-     * TargetVariables   FunctionSet    MaximumLength
-     *               |   |                 |
-     *               V   V                 |
-     *             Grammar <---------------+ 
+     * TargetVariables   FunctionSet    MaximumLength    NumberOfLatentVariables
+     *               |   |                 |                   |
+     *               V   V                 |                   |
+     *             Grammar <---------------+-------------------
      *                |
      *                V
@@ -539,4 +572,10 @@
 
       MaximumLengthParameter.Value.ValueChanged += MaximumLengthChanged;
+
+      NumberOfLatentVariablesParameter.Value.ValueChanged += NumLatentVariablesChanged;
+    }
+
+    private void NumLatentVariablesChanged(object sender, EventArgs e) {
+      UpdateGrammarAndEncoding();
     }
 
@@ -599,4 +638,7 @@
       return n.Symbol.Name.StartsWith("θ");
     }
+    private static bool IsLatentVariableNode(ISymbolicExpressionTreeNode n) {
+      return n.Symbol.Name.StartsWith("λ");
+    }
 
 
@@ -616,5 +658,8 @@
       var g = CreateGrammar();
       foreach (var targetVar in TargetVariables.CheckedItems) {
-        encoding = encoding.Add(new SymbolicExpressionTreeEncoding(targetVar+"_tree",g, MaximumLength, MaximumLength)); // only limit by length
+        encoding = encoding.Add(new SymbolicExpressionTreeEncoding(targetVar + "_tree", g, MaximumLength, MaximumLength)); // only limit by length
+      }
+      for (int i = 1; i <= NumberOfLatentVariables; i++) {
+        encoding = encoding.Add(new SymbolicExpressionTreeEncoding("λ" + i + "_tree", g, MaximumLength, MaximumLength));
       }
       Encoding = encoding;
@@ -642,4 +687,10 @@
         g.AddTerminalSymbol("θ" + i); // numeric parameter for which the value is optimized using AutoDiff
       }
+
+      // generate symbols for latent variables
+      for (int i = 1; i <= NumberOfLatentVariables; i++) {
+        g.AddTerminalSymbol("λ" + i); // numeric parameter for which the value is optimized using AutoDiff
+      }
+
       return g;
     }
