Index: branches/3106_AnalyticContinuedFractionsRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/ContinuedFractionRegression/Algorithm.cs
===================================================================
--- branches/3106_AnalyticContinuedFractionsRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/ContinuedFractionRegression/Algorithm.cs	(revision 17984)
+++ branches/3106_AnalyticContinuedFractionsRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/ContinuedFractionRegression/Algorithm.cs	(revision 17985)
@@ -1,4 +1,5 @@
 ﻿using System;
 using System.Collections.Generic;
+using System.Diagnostics;
 using System.Linq;
 using System.Threading;
@@ -162,8 +163,12 @@
       for (int gen = 1; gen <= numGen && !cancellationToken.IsCancellationRequested; gen++) {
         /* mutate each current solution in the population */
-        var pop_mu = Mutate(pop, mutationRate, rand);
+        var pop_mu = Mutate(pop, mutationRate, rand, trainingData);
         /* generate new population by recombination mechanism */
-        var pop_r = RecombinePopulation(pop_mu, rand, nVars);
-
+        var pop_r = RecombinePopulation(pop_mu, rand, nVars, trainingData);
+
+        // Paper:
+        // A period of individual search operation is performed every generation on all current solutions.
+
+        // Statement by authors:
         // "We ran the Local Search after Mutation and recombination operations. We executed the local-search only on the Current solutions."
         // "We executed the MaintainInvariant() in the following steps:
@@ -175,8 +180,9 @@
         /* local search optimization of current solutions */
         foreach (var agent in pop_r.IterateLevels()) {
-          LocalSearchSimplex(localSearchIterations, localSearchRestarts, localSearchTolerance, evalDelta, agent.current, ref agent.currentObjValue, trainingData, rand);
+          LocalSearchSimplex(localSearchIterations, localSearchRestarts, localSearchTolerance, evalDelta, agent, trainingData, rand);
+          Debug.Assert(agent.pocketObjValue < agent.currentObjValue);
         }
         foreach (var agent in pop_r.IteratePostOrder()) agent.MaintainInvariant(); // post-order to make sure that the root contains the best model 
-
+        foreach (var agent in pop_r.IteratePostOrder()) agent.AssertInvariant();
 
         // for detecting stagnation we track the best objective value since the last reset 
@@ -252,4 +258,7 @@
         agent.MaintainInvariant();
       }
+
+      foreach (var agent in pop.IteratePostOrder()) agent.AssertInvariant();
+
       return pop;
     }
@@ -266,10 +275,13 @@
       root.pocketObjValue = Evaluate(root.pocket, trainingData, Delta);
 
-      foreach (var agent in root.IteratePreOrder()) { agent.MaintainInvariant(); } // Here we push the newly created model down the hierarchy. 
-    }
-
-
-
-    private Agent RecombinePopulation(Agent pop, IRandom rand, int nVars) {
+      foreach (var agent in root.IteratePreOrder()) { agent.MaintainInvariant(); } // Here we use pre-order traversal push the newly created model down the hierarchy. 
+
+      foreach (var agent in root.IteratePostOrder()) agent.AssertInvariant();
+
+    }
+
+
+
+    private Agent RecombinePopulation(Agent pop, IRandom rand, int nVars, double[,] trainingData) {
       var l = pop;
 
@@ -284,13 +296,13 @@
         // Step 4. These steps are executed sequentially from 1 to 4. Similarly, in the
         // recombination of lower-level subpopulations, we will have the new current
-        // (the supporters generated at the previous level) as the leader of the subpopulation.
-        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);
+        // (the supporters generated at the previous level) as the leader of the subpopulation."
+        Recombine(l, s1, SelectRandomOp(rand), rand, nVars, trainingData);
+        Recombine(s3, l, SelectRandomOp(rand), rand, nVars, trainingData);
+        Recombine(s1, s2, SelectRandomOp(rand), rand, nVars, trainingData);
+        Recombine(s2, s3, SelectRandomOp(rand), rand, nVars, trainingData);
 
         // recombination works from top to bottom
         foreach (var child in pop.children) {
-          RecombinePopulation(child, rand, nVars);
+          RecombinePopulation(child, rand, nVars, trainingData);
         }
 
@@ -299,5 +311,7 @@
     }
 
-    private static ContinuedFraction Recombine(ContinuedFraction p1, ContinuedFraction p2, Func<bool[], bool[], bool[]> op, IRandom rand, int nVars) {
+    private ContinuedFraction Recombine(Agent a, Agent b, Func<bool[], bool[], bool[]> op, IRandom rand, int nVars, double[,] trainingData) {
+      ContinuedFraction p1 = a.pocket;
+      ContinuedFraction p2 = b.pocket;
       ContinuedFraction ch = new ContinuedFraction() { h = new Term[p1.h.Length] };
       /* apply a recombination operator chosen uniformly at random on variables of two parents into offspring */
@@ -332,12 +346,11 @@
         ch.h[i].beta = p1.h[i].beta + (rand.NextDouble() * 5 - 1) * (p2.h[i].beta - p1.h[i].beta) / 3.0;
       }
-      // return LocalSearchSimplex(ch, trainingData); // The paper has a local search step here. 
-      // The authors have stated that local search is executed after mutation and recombination
-      // for the current solutions.
-      // Local search and MaintainInvariant is called in the main loop (Alg 1)
+
+      a.current = ch;
+      LocalSearchSimplex(LocalSearchIterations, LocalSearchRestarts, LocalSearchTolerance, Delta, a, trainingData, rand);
       return ch;
     }
 
-    private Agent Mutate(Agent pop, double mutationRate, IRandom rand) {
+    private Agent Mutate(Agent pop, double mutationRate, IRandom rand, double[,] trainingData) {
       foreach (var agent in pop.IterateLevels()) {
         if (rand.NextDouble() < mutationRate) {
@@ -347,4 +360,9 @@
           else
             ModifyVariable(agent.current, rand); // soft mutation
+
+          // Finally, the local search operation is executed on the mutated solution in order to optimize
+          // non-zero coefficients. We do not apply mutation on pocket solutions because we consider them as a "collective memory"
+          // of good models visited in the past. They influence the search process via recombination only.
+          LocalSearchSimplex(LocalSearchIterations, LocalSearchRestarts, LocalSearchTolerance, Delta, agent, trainingData, rand);
         }
       }
@@ -365,6 +383,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]) {  // CHECK: paper uses varAt()
-          h.vars[vIdx] = false;  // CHECK: paper uses varAt()
+        if (cfrac.vars[vIdx]) {
+          h.vars[vIdx] = false;
           h.coef[vIdx] = coinToss(0, h.coef[vIdx]);
         } else {
@@ -372,5 +390,5 @@
            * or 'Replace' the coefficient with a random value between -3 and 3 at random */
           if (!h.vars[vIdx]) {
-            h.vars[vIdx] = true;  // CHECK: paper uses varAt()
+            h.vars[vIdx] = true;
             h.coef[vIdx] = coinToss(0, rand.NextDouble() * 6 - 3);
           }
@@ -378,5 +396,5 @@
       }
       /* toggle the randomly selected variable */
-      cfrac.vars[vIdx] = !cfrac.vars[vIdx];  // CHECK: paper uses varAt()
+      cfrac.vars[vIdx] = !cfrac.vars[vIdx];
     }
 
@@ -384,5 +402,5 @@
       /* randomly select a variable which is turned ON */
       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()
+      for (int i = 0; i < cfrac.vars.Length; i++) if (cfrac.vars[i]) candVars.Add(i);
       if (candVars.Count == 0) return; // no variable active
       var vIdx = candVars[rand.Next(candVars.Count)];
@@ -392,5 +410,5 @@
 
       /* modify the coefficient value */
-      if (h.vars[vIdx]) {  // CHECK: paper uses varAt()
+      if (h.vars[vIdx]) {
         h.coef[vIdx] = 0.0;
       } else {
@@ -398,5 +416,5 @@
       }
       /* Toggle the randomly selected variable */
-      h.vars[vIdx] = !h.vars[vIdx]; // CHECK: paper uses varAt()
+      h.vars[vIdx] = !h.vars[vIdx];
     }
 
@@ -465,5 +483,5 @@
 
 
-    private static void LocalSearchSimplex(int iterations, int restarts, double tolerance, double delta, ContinuedFraction ch, ref double quality, double[,] trainingData, IRandom rand) {
+    private static void LocalSearchSimplex(int iterations, int restarts, double tolerance, double delta, Agent a, double[,] trainingData, IRandom rand) {
       double uniformPeturbation = 1.0;
       int maxEvals = iterations;
@@ -472,5 +490,6 @@
       int numSelectedRows = numRows / 5; // 20% of the training samples
 
-      quality = Evaluate(ch, trainingData, delta); // get quality with origial coefficients
+      var ch = a.current;
+      var quality = Evaluate(ch, trainingData, delta); // get quality with original coefficients
 
       double[] origCoeff = ExtractCoeff(ch);
@@ -507,6 +526,17 @@
       } // reps
 
-      SetCoeff(ch, bestCoeff);
-      quality = bestQuality;
+      SetCoeff(a.current, bestCoeff);
+      a.currentObjValue = bestQuality;
+
+      // Unclear what the following means exactly. 
+      // 
+      // "We remind again that
+      // each solution corresponds to a single model, this means that if a current model becomes better than its corresponding
+      // pocket model (in terms of the guiding function of the solution), then an individual search optimization step is also
+      // performed on the pocket solution/ model before we swap it with the current solution/ model. Individual search can then
+      // make a current model better than the pocket model (again, according to the guiding function), and in that case they
+      // switch positions within the agent that contains both of them."
+
+      a.MaintainPocketCurrentInvariant();
     }
 
@@ -531,5 +561,5 @@
         coeff.Add(hi.beta);
         for (int vIdx = 0; vIdx < hi.vars.Length; vIdx++) {
-          if (hi.vars[vIdx]) coeff.Add(hi.coef[vIdx]);
+          if (hi.vars[vIdx] && hi.coef[vIdx] != 0) coeff.Add(hi.coef[vIdx]); // paper states twice (for mutation and recombination) that non-zero coefficients are optimized
         }
       }
@@ -542,5 +572,5 @@
         hi.beta = curCoeff[k++];
         for (int vIdx = 0; vIdx < hi.vars.Length; vIdx++) {
-          if (hi.vars[vIdx]) hi.coef[vIdx] = curCoeff[k++];
+          if (hi.vars[vIdx] && hi.coef[vIdx] != 0) hi.coef[vIdx] = curCoeff[k++]; // paper states twice (for mutation and recombination) that non-zero coefficients are optimized
         }
       }
