Index: /branches/3106_AnalyticContinuedFractionsRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/ContinuedFractionRegression/Algorithm.cs
===================================================================
--- /branches/3106_AnalyticContinuedFractionsRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/ContinuedFractionRegression/Algorithm.cs	(revision 17982)
+++ /branches/3106_AnalyticContinuedFractionsRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/ContinuedFractionRegression/Algorithm.cs	(revision 17983)
@@ -7,5 +7,9 @@
 using HeuristicLab.Core;
 using HeuristicLab.Data;
+using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
+using HeuristicLab.Parameters;
 using HeuristicLab.Problems.DataAnalysis;
+using HeuristicLab.Problems.DataAnalysis.Symbolic;
+using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;
 using HeuristicLab.Random;
 
@@ -15,4 +19,41 @@
   [StorableType("7A375270-EAAF-4AD1-82FF-132318D20E09")]
   public class Algorithm : FixedDataAnalysisAlgorithm<IRegressionProblem> {
+    private const string MutationRateParameterName = "MutationRate";
+    private const string DepthParameterName = "Depth";
+    private const string NumGenerationsParameterName = "Depth";
+
+    #region parameters
+    public IFixedValueParameter<PercentValue> MutationRateParameter => (IFixedValueParameter<PercentValue>)Parameters[MutationRateParameterName];
+    public double MutationRate {
+      get { return MutationRateParameter.Value.Value; }
+      set { MutationRateParameter.Value.Value = value; }
+    }
+    public IFixedValueParameter<IntValue> DepthParameter => (IFixedValueParameter<IntValue>)Parameters[DepthParameterName];
+    public int Depth {
+      get { return DepthParameter.Value.Value; }
+      set { DepthParameter.Value.Value = value; }
+    }
+    public IFixedValueParameter<IntValue> NumGenerationsParameter => (IFixedValueParameter<IntValue>)Parameters[NumGenerationsParameterName];
+    public int NumGenerations {
+      get { return NumGenerationsParameter.Value.Value; }
+      set { NumGenerationsParameter.Value.Value = value; }
+    }
+    #endregion
+
+    // storable ctor
+    [StorableConstructor]
+    public Algorithm(StorableConstructorFlag _) : base(_) { }
+
+    // cloning ctor
+    public Algorithm(Algorithm original, Cloner cloner) : base(original, cloner) { }
+
+
+    // default ctor
+    public Algorithm() : base() {
+      Parameters.Add(new FixedValueParameter<PercentValue>(MutationRateParameterName, "Mutation rate (default 10%)", new PercentValue(0.1)));
+      Parameters.Add(new FixedValueParameter<IntValue>(DepthParameterName, "Depth of the continued fraction representation (default 6)", new IntValue(6)));
+      Parameters.Add(new FixedValueParameter<IntValue>(NumGenerationsParameterName, "The maximum number of generations (default 200)", new IntValue(200)));
+    }
+
     public override IDeepCloneable Clone(Cloner cloner) {
       throw new NotImplementedException();
@@ -25,6 +66,7 @@
         problemData.TrainingIndices);
       var nVars = x.GetLength(1) - 1;
-      var rand = new MersenneTwister(31415);
-      CFRAlgorithm(nVars, depth: 6, 0.10, x, out var best, out var bestObj, rand, numGen: 200, stagnatingGens: 5, cancellationToken);
+      var seed = new System.Random().Next();
+      var rand = new MersenneTwister((uint)seed);
+      CFRAlgorithm(nVars, Depth, MutationRate, x, out var best, out var bestObj, rand, NumGenerations, stagnatingGens: 5, cancellationToken);
     }
 
@@ -47,8 +89,8 @@
         /* local search optimization of current solutions */
         foreach (var agent in pop_r.IterateLevels()) {
-          LocalSearchSimplex(agent.current, ref agent.currentObjValue, trainingData, rand);
-        }
-
-        foreach (var agent in pop_r.IteratePostOrder()) agent.MaintainInvariant(); // Deviates from Alg1 in paper 
+          LocalSearchSimplex(agent.current, ref agent.currentObjValue, trainingData, rand); // CHECK paper states that pocket might also be optimized. Unclear how / when invariants are maintained.
+        }
+
+        foreach (var agent in pop_r.IteratePostOrder()) agent.MaintainInvariant(); // CHECK deviates from Alg1 in paper 
 
         /* replace old population with evolved population */
@@ -56,19 +98,22 @@
 
         /* keep track of the best solution */
-        if (bestObj > pop.pocketObjValue) {
+        if (bestObj > pop.pocketObjValue) { // CHECK: comparison obviously wrong in the paper
           best = pop.pocket;
           bestObj = pop.pocketObjValue;
           bestObjGen = gen;
-          Results.AddOrUpdateResult("MSE (best)", new DoubleValue(bestObj));
+          // Results.AddOrUpdateResult("MSE (best)", new DoubleValue(bestObj));
+          // Results.AddOrUpdateResult("Solution", CreateSymbolicRegressionSolution(best, Problem.ProblemData));
         }
 
 
         if (gen > bestObjGen + stagnatingGens) {
-          bestObjGen = gen; // wait at least stagnatingGens until resetting again
-          // Reset(pop, nVars, depth, rand, trainingData);
-          InitialPopulation(nVars, depth, rand, trainingData);
-        }
-      }
-    }
+          bestObjGen = gen; // CHECK: unspecified in the paper: wait at least stagnatingGens until resetting again
+          Reset(pop, nVars, depth, rand, trainingData);
+          // InitialPopulation(nVars, depth, rand, trainingData); CHECK reset is not specified in the paper
+        }
+      }
+    }
+
+
 
     private Agent InitialPopulation(int nVars, int depth, IRandom rand, double[,] trainingData) {
@@ -116,16 +161,26 @@
     private Agent RecombinePopulation(Agent pop, IRandom rand, int nVars) {
       var l = pop;
+
       if (pop.children.Count > 0) {
         var s1 = pop.children[0];
         var s2 = pop.children[1];
         var s3 = pop.children[2];
-        l.current = Recombine(l.pocket, s1.current, SelectRandomOp(rand), rand, nVars);
-        s3.current = Recombine(s3.pocket, l.current, SelectRandomOp(rand), rand, nVars);
-        s1.current = Recombine(s1.pocket, s2.current, SelectRandomOp(rand), rand, nVars);
-        s2.current = Recombine(s2.pocket, s3.current, SelectRandomOp(rand), rand, nVars);
-      }
-
-      foreach (var child in pop.children) {
-        RecombinePopulation(child, rand, nVars);
+
+        // CHECK Deviates from paper (recombine all models in the current pop before updating the population)
+        var l_current = Recombine(l.pocket, s1.current, SelectRandomOp(rand), rand, nVars);
+        var s3_current = Recombine(s3.pocket, l.current, SelectRandomOp(rand), rand, nVars);
+        var s1_current = Recombine(s1.pocket, s2.current, SelectRandomOp(rand), rand, nVars);
+        var s2_current = Recombine(s2.pocket, s3.current, SelectRandomOp(rand), rand, nVars);
+
+        // recombination works from top to bottom
+        // CHECK do we use the new current solutions (s1_current .. s3_current) already in the next levels?
+        foreach (var child in pop.children) {
+          RecombinePopulation(child, rand, nVars);
+        }
+
+        l.current = l_current;
+        s3.current = s3_current;
+        s1.current = s1_current;
+        s2.current = s2_current;
       }
       return pop;
@@ -158,5 +213,5 @@
     private static ContinuedFraction Recombine(ContinuedFraction p1, ContinuedFraction p2, Func<bool[], bool[], bool[]> op, IRandom rand, int nVars) {
       ContinuedFraction ch = new ContinuedFraction() { h = new Term[p1.h.Length] };
-      /* apply a recombination operator chosen uniformly at random on variable sof two parents into offspring */
+      /* apply a recombination operator chosen uniformly at random on variables of two parents into offspring */
       ch.vars = op(p1.vars, p2.vars);
 
@@ -168,7 +223,7 @@
         /* recombine coefficient values for variables */
         var coefx = new double[nVars];
-        var varsx = new bool[nVars]; // TODO: deviates from paper -> check
+        var varsx = new bool[nVars]; // CHECK: deviates from paper, probably forgotten in the pseudo-code
         for (int vi = 1; vi < nVars; vi++) {
-          if (ch.vars[vi]) {
+          if (ch.vars[vi]) {  // CHECK: paper uses featAt()
             if (varsa[vi] && varsb[vi]) {
               coefx[vi] = coefa[vi] + (rand.NextDouble() * 5 - 1) * (coefb[vi] - coefa[vi]) / 3.0;
@@ -190,5 +245,5 @@
       }
       /* update current solution and apply local search */
-      // return LocalSearchSimplex(ch, trainingData); // Deviates from paper because Alg1 also has LocalSearch after Recombination
+      // return LocalSearchSimplex(ch, trainingData); // CHECK: Deviates from paper because Alg1 also has LocalSearch after Recombination
       return ch;
     }
@@ -220,6 +275,6 @@
         /* Case 1: cfrac variable is turned ON: Turn OFF the variable, and either 'Remove' or 
          * 'Remember' the coefficient value at random */
-        if (cfrac.vars[vIdx]) {
-          h.vars[vIdx] = false;
+        if (cfrac.vars[vIdx]) {  // CHECK: paper uses varAt()
+          h.vars[vIdx] = false;  // CHECK: paper uses varAt()
           h.coef[vIdx] = coinToss(0, h.coef[vIdx]);
         } else {
@@ -227,5 +282,5 @@
            * or 'Replace' the coefficient with a random value between -3 and 3 at random */
           if (!h.vars[vIdx]) {
-            h.vars[vIdx] = true;
+            h.vars[vIdx] = true;  // CHECK: paper uses varAt()
             h.coef[vIdx] = coinToss(0, rand.NextDouble() * 6 - 3);
           }
@@ -233,12 +288,13 @@
       }
       /* toggle the randomly selected variable */
-      cfrac.vars[vIdx] = !cfrac.vars[vIdx];
+      cfrac.vars[vIdx] = !cfrac.vars[vIdx];  // CHECK: paper uses varAt()
     }
 
     private void ModifyVariable(ContinuedFraction cfrac, IRandom rand) {
       /* randomly select a variable which is turned ON */
-      var candVars = cfrac.vars.Count(vi => vi);
-      if (candVars == 0) return; // no variable active
-      var vIdx = rand.Next(candVars);
+      var candVars = new List<int>();
+      for (int i = 0; i < cfrac.vars.Length; i++) if (cfrac.vars[i]) candVars.Add(i);  // CHECK: paper uses varAt()
+      if (candVars.Count == 0) return; // no variable active
+      var vIdx = candVars[rand.Next(candVars.Count)];
 
       /* randomly select a term (h) of continued fraction */
@@ -246,5 +302,5 @@
 
       /* modify the coefficient value*/
-      if (h.vars[vIdx]) {
+      if (h.vars[vIdx]) {  // CHECK: paper uses varAt()
         h.coef[vIdx] = 0.0;
       } else {
@@ -252,5 +308,5 @@
       }
       /* Toggle the randomly selected variable */
-      h.vars[vIdx] = !h.vars[vIdx];
+      h.vars[vIdx] = !h.vars[vIdx]; // CHECK: paper uses varAt()
     }
 
@@ -268,5 +324,5 @@
         sum += res * res;
       }
-      var delta = 0.1; // TODO
+      var delta = 0.1;
       return sum / trainingData.GetLength(0) * (1 + delta * cfrac.vars.Count(vi => vi));
     }
@@ -281,4 +337,6 @@
         res = numerator / denom;
       }
+      var h0 = cfrac.h[0];
+      res += h0.beta + dot(h0.vars, h0.coef, dataPoint);
       return res;
     }
@@ -329,6 +387,4 @@
 
         var newQuality = Evaluate(ch, trainingData);
-
-        // TODO: optionally use regularization (ridge / LASSO)
 
         if (newQuality < bestQuality) {
@@ -377,4 +433,64 @@
       }
     }
+
+    Symbol addSy = new Addition();
+    Symbol mulSy = new Multiplication();
+    Symbol divSy = new Division();
+    Symbol startSy = new StartSymbol();
+    Symbol progSy = new ProgramRootSymbol();
+    Symbol constSy = new Constant();
+    Symbol varSy = new Problems.DataAnalysis.Symbolic.Variable();
+
+    private ISymbolicRegressionSolution CreateSymbolicRegressionSolution(ContinuedFraction cfrac, IRegressionProblemData problemData) {
+      var variables = problemData.AllowedInputVariables.ToArray();
+      ISymbolicExpressionTreeNode res = null;
+      for (int i = cfrac.h.Length - 1; i > 1; i -= 2) {
+        var hi = cfrac.h[i];
+        var hi1 = cfrac.h[i - 1];
+        var denom = CreateLinearCombination(hi.vars, hi.coef, variables, hi.beta);
+        if (res != null) {
+          denom.AddSubtree(res);
+        }
+
+        var numerator = CreateLinearCombination(hi1.vars, hi1.coef, variables, hi1.beta);
+
+        res = divSy.CreateTreeNode();
+        res.AddSubtree(numerator);
+        res.AddSubtree(denom);
+      }
+
+      var h0 = cfrac.h[0];
+      var h0Term = CreateLinearCombination(h0.vars, h0.coef, variables, h0.beta);
+      h0Term.AddSubtree(res);
+
+      var progRoot = progSy.CreateTreeNode();
+      var start = startSy.CreateTreeNode();
+      progRoot.AddSubtree(start);
+      start.AddSubtree(h0Term);
+
+      var model = new SymbolicRegressionModel(problemData.TargetVariable, new SymbolicExpressionTree(progRoot), new SymbolicDataAnalysisExpressionTreeBatchInterpreter());
+      var sol = new SymbolicRegressionSolution(model, (IRegressionProblemData)problemData.Clone());
+      return sol;
+    }
+
+    private ISymbolicExpressionTreeNode CreateLinearCombination(bool[] vars, double[] coef, string[] variables, double beta) {
+      var sum = addSy.CreateTreeNode();
+      for (int i = 0; i < vars.Length; i++) {
+        if (vars[i]) {
+          var varNode = (VariableTreeNode)varSy.CreateTreeNode();
+          varNode.Weight = coef[i];
+          varNode.VariableName = variables[i];
+          sum.AddSubtree(varNode);
+        }
+      }
+      sum.AddSubtree(CreateConstant(beta));
+      return sum;
+    }
+
+    private ISymbolicExpressionTreeNode CreateConstant(double value) {
+      var constNode = (ConstantTreeNode)constSy.CreateTreeNode();
+      constNode.Value = value;
+      return constNode;
+    }
   }
 
