Index: branches/HeuristicLab.EvolutionTracking/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Tracking/SchemaDiversification/SchemaEvaluator.cs
===================================================================
--- branches/HeuristicLab.EvolutionTracking/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Tracking/SchemaDiversification/SchemaEvaluator.cs	(revision 15482)
+++ branches/HeuristicLab.EvolutionTracking/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Tracking/SchemaDiversification/SchemaEvaluator.cs	(revision 15906)
@@ -21,4 +21,5 @@
 
 using System;
+using System.Collections.Concurrent;
 using System.Collections.Generic;
 using System.Linq;
@@ -41,5 +42,5 @@
     private const string MinimumSchemaFrequencyParameterName = "MinimumSchemaFrequency";
     private const string MinimumPhenotypicSimilarityParameterName = "MinimumPhenotypicSimilarity";
-    private const string ReplacementRatioParameterName = "ReplacementRatio";
+    private const string MutationRateParameterName = "MutationRate";
     private const string SchemaParameterName = "Schema";
     private const string PopulationSizeParameterName = "PopulationSize";
@@ -52,16 +53,15 @@
     private const string MutatorParameterName = "Mutator";
     private const string CrossoverParameterName = "Crossover";
-    private const string RandomReplacementParameterName = "RandomReplacement";
     private const string NumberOfChangedTreesParameterName = "NumberOfChangedTrees";
     private const string ExecuteInParallelParameterName = "ExecuteInParallel";
     private const string MaxDegreeOfParalellismParameterName = "MaxDegreeOfParallelism";
     private const string ExclusiveMatchingParameterName = "ExclusiveMatching";
-    private const string UseAdaptiveReplacementRatioParameterName = "UseAdaptiveReplacementRatio";
+    private const string UseAdaptiveMutationRateParameterName = "UseAdaptiveMutationRate";
     private const string StrictSchemaMatchingParameterName = "StrictSchemaMatching";
     #endregion
 
     #region parameters
-    public ILookupParameter<BoolValue> UseAdaptiveReplacementRatioParameter {
-      get { return (ILookupParameter<BoolValue>)Parameters[UseAdaptiveReplacementRatioParameterName]; }
+    public ILookupParameter<BoolValue> UseAdaptiveMutationRateParameter {
+      get { return (ILookupParameter<BoolValue>)Parameters[UseAdaptiveMutationRateParameterName]; }
     }
     public ILookupParameter<BoolValue> ExclusiveMatchingParameter {
@@ -89,7 +89,4 @@
       get { return (ILookupParameter<BoolValue>)Parameters[ApplyLinearScalingParameterName]; }
     }
-    public ILookupParameter<BoolValue> RandomReplacementParameter {
-      get { return (ILookupParameter<BoolValue>)Parameters[RandomReplacementParameterName]; }
-    }
     public ILookupParameter<ISymbolicExpressionTreeCrossover> CrossoverParameter {
       get { return (ILookupParameter<ISymbolicExpressionTreeCrossover>)Parameters[CrossoverParameterName]; }
@@ -107,6 +104,6 @@
       get { return (ILookupParameter<PercentValue>)Parameters[MinimumSchemaFrequencyParameterName]; }
     }
-    public ILookupParameter<PercentValue> ReplacementRatioParameter {
-      get { return (ILookupParameter<PercentValue>)Parameters[ReplacementRatioParameterName]; }
+    public ILookupParameter<PercentValue> MutationRateParameter {
+      get { return (ILookupParameter<PercentValue>)Parameters[MutationRateParameterName]; }
     }
     public ILookupParameter<PercentValue> MinimumPhenotypicSimilarityParameter {
@@ -123,7 +120,6 @@
     #region parameter properties
     public PercentValue MinimumSchemaFrequency { get { return MinimumSchemaFrequencyParameter.ActualValue; } }
-    public PercentValue ReplacementRatio { get { return ReplacementRatioParameter.ActualValue; } }
+    public PercentValue MutationRate { get { return MutationRateParameter.ActualValue; } }
     public PercentValue MinimumPhenotypicSimilarity { get { return MinimumPhenotypicSimilarityParameter.ActualValue; } }
-    public BoolValue RandomReplacement { get { return RandomReplacementParameter.ActualValue; } }
     public IntValue NumberOfChangedTrees { get { return NumberOfChangedTreesParameter.ActualValue; } }
     #endregion
@@ -135,5 +131,5 @@
 
     [Storable]
-    public string ReplacementRule { get; set; }
+    public string MutationRateUpdateRule { get; set; }
 
     [Storable]
@@ -171,5 +167,5 @@
       Parameters.Add(new LookupParameter<ISymbolicExpressionTree>(SchemaParameterName, "The current schema to be evaluated"));
       Parameters.Add(new LookupParameter<PercentValue>(MinimumSchemaFrequencyParameterName));
-      Parameters.Add(new LookupParameter<PercentValue>(ReplacementRatioParameterName));
+      Parameters.Add(new LookupParameter<PercentValue>(MutationRateParameterName));
       Parameters.Add(new LookupParameter<PercentValue>(MinimumPhenotypicSimilarityParameterName));
       Parameters.Add(new LookupParameter<ISymbolicExpressionTree>(PopulationSizeParameterName));
@@ -183,10 +179,9 @@
       Parameters.Add(new LookupParameter<ISymbolicExpressionTreeManipulator>(MutatorParameterName));
       Parameters.Add(new LookupParameter<ISymbolicExpressionTreeCrossover>(CrossoverParameterName));
-      Parameters.Add(new LookupParameter<BoolValue>(RandomReplacementParameterName));
       Parameters.Add(new LookupParameter<IntValue>(NumberOfChangedTreesParameterName));
       Parameters.Add(new LookupParameter<BoolValue>(ExecuteInParallelParameterName));
       Parameters.Add(new LookupParameter<IntValue>(MaxDegreeOfParalellismParameterName));
       Parameters.Add(new LookupParameter<BoolValue>(ExclusiveMatchingParameterName));
-      Parameters.Add(new LookupParameter<BoolValue>(UseAdaptiveReplacementRatioParameterName));
+      Parameters.Add(new LookupParameter<BoolValue>(UseAdaptiveMutationRateParameterName));
       #endregion
     }
@@ -218,20 +213,24 @@
 
       var individuals = ExecutionContext.Scope.SubScopes; // the scopes represent the individuals
-      var trees = individuals.Select(x => (ISymbolicExpressionTree)x.Variables["SymbolicExpressionTree"].Value).ToList();
+      var trees = new ISymbolicExpressionTree[individuals.Count];
+      var qualities = new double[individuals.Count];
+
+      for (int i = 0; i < individuals.Count; ++i) {
+        trees[i] = (ISymbolicExpressionTree)individuals[i].Variables["SymbolicExpressionTree"].Value;
+        qualities[i] = ((DoubleValue)individuals[i].Variables["Quality"].Value).Value;
+      }
 
       var random = RandomParameter.ActualValue;
-
-      var s = Schema;
-      var sRoot = s.Root.GetSubtree(0).GetSubtree(0);
-      int countThreshold = (int)Math.Max(2, Math.Round(MinimumSchemaFrequency.Value * individuals.Count));
+      var sRoot = Schema.Root.GetSubtree(0).GetSubtree(0);
 
       // first apply the length and root equality checks in order to filter the individuals
       var exclusiveMatching = ExclusiveMatchingParameter.ActualValue.Value;
       var filtered = new List<int>();
-      for (int i = 0; i < individuals.Count; ++i) {
+
+      for (int i = 0; i < trees.Length; ++i) {
         if (exclusiveMatching && individuals[i].Variables.ContainsKey("AlreadyMatched")) continue;
         var t = trees[i];
         var tRoot = t.Root.GetSubtree(0).GetSubtree(0);
-        if (t.Length < s.Length || !qm.EqualityComparer.Equals(tRoot, sRoot)) continue;
+        if (t.Length < Schema.Length || !qm.EqualityComparer.Equals(tRoot, sRoot)) continue;
         filtered.Add(i);
       }
@@ -239,4 +238,5 @@
       // if we don't have enough filtered individuals, then we are done
       // if the schema exceeds the minimum frequency, it gets processed further
+      int countThreshold = (int)Math.Max(2, Math.Round(MinimumSchemaFrequency.Value * individuals.Count));
       if (filtered.Count < countThreshold) {
         return base.Apply();
@@ -244,41 +244,41 @@
 
       // check if the filtered individuals match the schema
+      var matching = new List<int>();
       var sNodes = QueryMatch.InitializePostOrder(sRoot);
-      var matchingIndividuals = new ScopeList();
+
       bool executeInParallel = ExecuteInParallelParameter.ActualValue.Value;
       int maxDegreeOfParallelism = MaxDegreeOfParallelismParameter.ActualValue.Value;
 
       if (executeInParallel) {
-        var matching = new bool[filtered.Count];
-
-        Parallel.For(0, filtered.Count, new ParallelOptions { MaxDegreeOfParallelism = maxDegreeOfParallelism },
-          i => {
-            var index = filtered[i];
-            var t = trees[index];
-            var tNodes = QueryMatch.InitializePostOrder(t.Root.GetSubtree(0).GetSubtree(0));
-            if (qm.Match(tNodes, sNodes)) {
-              matching[i] = true;
-            }
-          });
-
-        matchingIndividuals.AddRange(filtered.Where((x, i) => matching[i]).Select(x => individuals[x]));
+        var partitioner = Partitioner.Create(0, filtered.Count);
+        var po = new ParallelOptions { MaxDegreeOfParallelism = maxDegreeOfParallelism };
+
+        Parallel.ForEach(partitioner, po, (range, loop) => {
+          var partial = new List<int>();
+          for (int i = range.Item1; i < range.Item2; ++i) {
+            var idx = filtered[i];
+            var tRoot = trees[idx].Root.GetSubtree(0).GetSubtree(0);
+            var tNodes = QueryMatch.InitializePostOrder(tRoot);
+
+            if (qm.Match(tNodes, sNodes)) { partial.Add(idx); }
+          }
+          lock (matching) { matching.AddRange(partial); }
+        });
       } else {
         for (int i = 0; i < filtered.Count; ++i) {
-          // break early if it becomes impossible to reach the minimum threshold
-          if (matchingIndividuals.Count + filtered.Count - i < countThreshold)
-            break;
-
           var index = filtered[i];
           var tRoot = trees[index].Root.GetSubtree(0).GetSubtree(0);
           var tNodes = QueryMatch.InitializePostOrder(tRoot);
-          if (qm.Match(tNodes, sNodes))
-            matchingIndividuals.Add(individuals[index]);
+          if (qm.Match(tNodes, sNodes)) {
+            matching.Add(index);
+          }
         }
       }
-
-      if (matchingIndividuals.Count < countThreshold) {
+      if (matching.Count < countThreshold) {
         return base.Apply();
       }
 
+      matching.Sort((a, b) => qualities[a].CompareTo(qualities[b])); // sort by ascending quality
+      var matchingIndividuals = matching.Select(x => individuals[x]).ToArray(); // fittest individual will be last in the array
       var similarity = CalculateSimilarity(matchingIndividuals, calculator, executeInParallel, maxDegreeOfParallelism);
       if (similarity < MinimumPhenotypicSimilarity.Value) {
@@ -286,35 +286,39 @@
       }
 
-      double replacementRatio;
-      var adaptiveReplacementRatio = UseAdaptiveReplacementRatioParameter.ActualValue.Value;
-
-      if (adaptiveReplacementRatio) {
-        var r = (double)matchingIndividuals.Count / individuals.Count;
-        replacementRatio = CalculateReplacementRatio(r);
+      double mutationRate;
+      var useAdaptiveMutationRate = UseAdaptiveMutationRateParameter.ActualValue.Value;
+
+      if (useAdaptiveMutationRate) {
+        var r = (double)matchingIndividuals.Length / individuals.Count;
+        mutationRate = CalculateMutationRate(r);
       } else {
-        replacementRatio = ReplacementRatio.Value;
-      }
-
-      int n = (int)Math.Floor(matchingIndividuals.Count * replacementRatio);
-      var individualsToReplace = RandomReplacement.Value ? matchingIndividuals.SampleRandomWithoutRepetition(random, n).ToList()
-                                                         : matchingIndividuals.OrderBy(x => (DoubleValue)x.Variables["Quality"].Value).Take(n).ToList();
-      var mutationOc = new OperationCollection { Parallel = false }; // cannot be parallel due to the before/after operators which insert vertices in the genealogy graph
-      var updateEstimatedValues = new OperationCollection { Parallel = true }; // evaluation should be done in parallel when possible
-      foreach (var ind in individualsToReplace) {
-        var mutatorOp = ExecutionContext.CreateChildOperation(SchemaManipulator, ind);
-        var updateOp = ExecutionContext.CreateChildOperation(updateQualityOperator, ind);
-        mutationOc.Add(mutatorOp);
-        updateEstimatedValues.Add(updateOp);
+        mutationRate = MutationRate.Value;
+      }
+
+      var mutations = new OperationCollection { Parallel = false }; // cannot be parallel due to the before/after operators which insert vertices in the genealogy graph
+      var updates = new OperationCollection { Parallel = true }; // evaluation should be done in parallel when possible
+
+      // use length - 1 because we don't want to mutate the best individual in each schema group (which could also be the overall elite)
+      for (int i = 0; i < matchingIndividuals.Length - 1; ++i) {
+        if (random.NextDouble() > mutationRate) continue;
+
+        var ind = matchingIndividuals[i];
+
+        var mutate = ExecutionContext.CreateChildOperation(SchemaManipulator, ind);
+        var update = ExecutionContext.CreateChildOperation(updateQualityOperator, ind);
+
+        mutations.Add(mutate);
+        updates.Add(update);
+
         if (exclusiveMatching)
           ind.Variables.Add(new Core.Variable("AlreadyMatched"));
       }
-
-      NumberOfChangedTrees.Value += individualsToReplace.Count; // a lock is not necessary here because the SchemaEvaluators cannot be executed in parallel
-
-      return new OperationCollection(mutationOc, updateEstimatedValues, base.Apply());
-    }
-
-    private double CalculateReplacementRatio(double r) {
-      switch (ReplacementRule) {
+      NumberOfChangedTrees.Value += mutations.Count;
+
+      return new OperationCollection(mutations, updates, base.Apply());
+    }
+
+    private double CalculateMutationRate(double r) {
+      switch (MutationRateUpdateRule) {
         case "f(x) = x": {
             return r;
@@ -340,19 +344,34 @@
     }
 
-    public static double CalculateSimilarity(ScopeList individuals, ISolutionSimilarityCalculator calculator, bool parallel = false, int nThreads = -1) {
+    public static double CalculateSimilarity(IScope[] individuals, ISolutionSimilarityCalculator calculator, bool parallel = false, int nThreads = -1) {
+      if (individuals.Length < 2)
+        return double.NaN;
+
       double similarity = 0;
-      int count = individuals.Count;
-      if (count < 2) return double.NaN;
+      int count = individuals.Length;
+      int n = count * (count - 1) / 2;
+
       if (parallel) {
-        var parallelOptions = new ParallelOptions { MaxDegreeOfParallelism = nThreads };
-        var simArray = new double[count - 1];
-        Parallel.For(0, count - 1, parallelOptions, i => {
-          double innerSim = 0;
+        var ii = new int[n];
+        var jj = new int[n];
+        int k = 0;
+        for (int i = 0; i < count - 1; ++i)
           for (int j = i + 1; j < count; ++j) {
-            innerSim += calculator.CalculateSolutionSimilarity(individuals[i], individuals[j]);
-          }
-          simArray[i] = innerSim;
+            ii[k] = i;
+            jj[k] = j;
+            ++k;
+          }
+        var po = new ParallelOptions { MaxDegreeOfParallelism = nThreads };
+        var partitioner = Partitioner.Create(0, n);
+        var locker = new object();
+        Parallel.ForEach(partitioner, new ParallelOptions { MaxDegreeOfParallelism = 4 }, (range, loop) => {
+          var partial = 0d;
+          for (int idx = range.Item1; idx < range.Item2; ++idx) {
+            int i = ii[idx];
+            int j = jj[idx];
+            partial += calculator.CalculateSolutionSimilarity(individuals[i], individuals[j]);
+          }
+          lock (locker) { similarity += partial; }
         });
-        similarity = simArray.Sum();
       } else {
         for (int i = 0; i < count - 1; ++i) {
@@ -362,5 +381,5 @@
         }
       }
-      return similarity / (count * (count - 1) / 2.0);
+      return similarity / n;
     }
   }
