Index: branches/PerformanceComparison/HeuristicLab.OptimizationExpertSystem.Common/3.3/ProblemCharacteristicAnalysis/QAP/QAPDirectedWalk.cs
===================================================================
--- branches/PerformanceComparison/HeuristicLab.OptimizationExpertSystem.Common/3.3/ProblemCharacteristicAnalysis/QAP/QAPBestDirectedWalk.cs	(revision 13861)
+++ branches/PerformanceComparison/HeuristicLab.OptimizationExpertSystem.Common/3.3/ProblemCharacteristicAnalysis/QAP/QAPDirectedWalk.cs	(revision 14429)
@@ -85,91 +85,88 @@
 
     public override IEnumerable<IResult> Calculate() {
+      IRandom random = Seed.HasValue ? new MersenneTwister((uint)Seed.Value) : new MersenneTwister();
+      var qap = (QuadraticAssignmentProblem)Problem;
       var pathCount = Paths;
-      var random = Seed.HasValue ? new MersenneTwister((uint)Seed.Value) : new MersenneTwister();
-      var bestImprovement = BestImprovement;
-      var qap = (QuadraticAssignmentProblem)Problem;
-      var N = qap.Weights.Rows;
-      var start = new Permutation(PermutationTypes.Absolute, N, random);
-      var end = start;
-      List<double> ips = new List<double>(), ups = new List<double>(), sd2 = new List<double>();
-      var pathsD1 = new List<List<Tuple<double, double>>>();
-      for (var i = 0; i < pathCount; i++) {
-        var trajD1 = new List<Tuple<double, double>>();
-        pathsD1.Add(trajD1);
-
-        if ((i + 1) % N == 0) {
-          start = new Permutation(PermutationTypes.Absolute, N, random);
-          end = start;
-        }
-        end = (Permutation)end.Clone();
-        Rot1(end);
-
-        var hist = new Tuple<Permutation, double>[5];
-        var walkDist = 0.0;
-        var prevVal = double.NaN;
-        var sumD2 = 0.0;
-        var inflectionPoints = 0;
-        var undulationPoints = 0;
-        var countPoints = 0;
-        var counter = 0;
-        var path = bestImprovement ? BestImprovementWalk(qap, start, QAPEvaluator.Apply(start, qap.Weights, qap.Distances), end)
-                                  : FirstImprovementWalk(qap, start, QAPEvaluator.Apply(start, qap.Weights, qap.Distances), end, random);
-        foreach (var next in path) {
-          if (hist[0] != null) {
-            var dist = Dist(next.Item1, hist[0].Item1);
-            walkDist += dist;
+
+      var perm = new Permutation(PermutationTypes.Absolute, qap.Weights.Rows, random);
+      var permutations = new List<Permutation> { perm };
+      while (permutations.Count < pathCount + 1) {
+        perm = (Permutation)permutations.Last().Clone();
+        BiasedShuffle(perm, random);
+        permutations.Add(perm);
+      }
+
+      var trajectories = Run(random, (QuadraticAssignmentProblem)Problem, permutations, BestImprovement).ToList();
+      var firstDerivatives = trajectories.Select(path => ApproximateDerivative(path).ToList()).ToList();
+      var secondDerivatives = firstDerivatives.Select(d1 => ApproximateDerivative(d1).ToList()).ToList();
+      
+      var props = GetCharacteristics(trajectories, firstDerivatives, secondDerivatives).ToDictionary(x => x.Item1, x => x.Item2);
+      foreach (var chara in characteristics.CheckedItems.Select(x => x.Value.Value)) {
+        if (chara == "Swap2.Sharpness") yield return new Result("Swap2.Sharpness", new DoubleValue(props["Sharpness"]));
+        if (chara == "Swap2.Bumpiness") yield return new Result("Swap2.Bumpiness", new DoubleValue(props["Bumpiness"]));
+        if (chara == "Swap2.Flatness") yield return new Result("Swap2.Flatness", new DoubleValue(props["Flatness"]));
+        if (chara == "Swap2.Steadiness") yield return new Result("Swap2.Steadiness", new DoubleValue(props["Steadiness"]));
+      }
+    }
+
+    public static IEnumerable<List<Tuple<Permutation, double>>> Run(IRandom random, QuadraticAssignmentProblem qap, IEnumerable<Permutation> permutations, bool bestImprovement = true) {
+      var iter = permutations.GetEnumerator();
+      if (!iter.MoveNext()) yield break;
+
+      var start = iter.Current;
+      while (iter.MoveNext()) {
+        var end = iter.Current;
+
+        var walk = (bestImprovement ? BestDirectedWalk(qap, start, end) : FirstDirectedWalk(random, qap, start, end)).ToList();
+        var max = walk.Max(x => x.Item2);
+        var min = walk.Min(x => x.Item2);
+        if (max > min)
+          yield return walk.Select(x => Tuple.Create(x.Item1, (x.Item2 - min) / (max - min))).ToList();
+        else yield return walk.Select(x => Tuple.Create(x.Item1, 0.0)).ToList();
+        start = end;
+      } // end paths
+    }
+
+    private IEnumerable<Tuple<string, double>> GetCharacteristics(List<List<Tuple<Permutation, double>>> f, List<List<Tuple<Permutation, double>>> f1, List<List<Tuple<Permutation, double>>> f2) {
+      var sharpness = f2.Average(x => Area(x));
+      var bumpiness = 0.0;
+      var flatness = 0.0;
+      var downPointing = f1.Where(x => x.Min(y => y.Item2) < 0).ToList();
+
+      var steadiness = 0.0;
+      foreach (var path in downPointing) {
+        steadiness += ComBelowZero(path);
+      }
+      if (downPointing.Count > 0) steadiness /= downPointing.Count;
+
+      for (var p = 0; p < f2.Count; p++) {
+        var bump = 0;
+        var flat = 0;
+        for (var i = 0; i < f2[p].Count - 1; i++) {
+          if ((f2[p][i].Item2 > 0 && f2[p][i + 1].Item2 < 0) || (f2[p][i].Item2 < 0 && f2[p][i + 1].Item2 > 0)) {
+            bump++;
+          } else if (f2[p][i].Item2 == 0) {
+            flat++;
           }
-          // from the past 5 values we can calculate the 2nd derivative
-          // first derivative in point 2 as differential between points 1 and 3
-          // first derivative in point 4 as differential between points 3 and 5
-          // second derivative in point 3 as differential between the first derivatives in points 2 and 4
-          hist[4] = hist[3];
-          hist[3] = hist[2];
-          hist[2] = hist[1];
-          hist[1] = hist[0];
-          hist[0] = next;
-          counter++;
-
-          if (counter < 3) continue;
-          var grad1 = (hist[0].Item2 - hist[2].Item2) / Dist(hist[0].Item1, hist[2].Item1);
-
-          if (!double.IsNaN(grad1) && !double.IsInfinity(grad1))
-            trajD1.Add(Tuple.Create(walkDist, grad1));
-
-          if (counter < 5) continue;
-          countPoints++;
-          var grad2 = (hist[2].Item2 - hist[4].Item2) / Dist(hist[2].Item1, hist[4].Item1);
-          var dgrad = (grad1 - grad2) / Dist(hist[1].Item1, hist[3].Item1);
-
-          if (double.IsNaN(dgrad) || double.IsInfinity(dgrad)) continue;
-          if (!double.IsNaN(prevVal)) {
-            if (prevVal < 0 && dgrad > 0 || prevVal > 0 && dgrad < 0)
-              inflectionPoints++;
-            else if (prevVal.IsAlmost(0) && dgrad.IsAlmost(0)
-                || prevVal.IsAlmost(0) && !dgrad.IsAlmost(0)
-                || !prevVal.IsAlmost(0) && dgrad.IsAlmost(0))
-              undulationPoints++;
-          }
-          sumD2 += Math.Abs(grad1 - grad2);
-          prevVal = dgrad;
-        }
-        start = end;
-        ips.Add(inflectionPoints / (double)countPoints);
-        ups.Add(undulationPoints / (double)countPoints);
-        sd2.Add(sumD2 / walkDist);
-      } // end paths
-      var avgZero = pathsD1.Select(path => path.SkipWhile(v => v.Item2 < 0).First().Item1 / path.Last().Item1).Median();
-      
-      foreach (var chara in characteristics.CheckedItems.Select(x => x.Value.Value)) {
-        if (chara == "Swap2.Sharpness") yield return new Result("Swap2.Sharpness", new DoubleValue(sd2.Average()));
-        if (chara == "Swap2.Bumpiness") yield return new Result("Swap2.Bumpiness", new DoubleValue(ips.Average()));
-        if (chara == "Swap2.Flatness") yield return new Result("Swap2.Flatness", new DoubleValue(ups.Average()));
-        if (chara == "Swap2.Steadiness") yield return new Result("Swap2.Steadiness", new DoubleValue(avgZero));
-      }
-    }
-
-    public IEnumerable<Tuple<Permutation, double>> BestImprovementWalk(QuadraticAssignmentProblem qap, Permutation start, double fitness, Permutation end) {
+        }
+        bumpiness += bump / (f2[p].Count - 1.0);
+        flatness += flat / (f2[p].Count - 1.0);
+      }
+      bumpiness /= f2.Count;
+      flatness /= f2.Count;
+      return new[] {
+      Tuple.Create("Sharpness", sharpness),
+      Tuple.Create("Bumpiness", bumpiness),
+      Tuple.Create("Flatness", flatness),
+      Tuple.Create("Steadiness", steadiness)
+    };
+    }
+
+    public static IEnumerable<Tuple<Permutation, double>> BestDirectedWalk(QuadraticAssignmentProblem qap, Permutation start, Permutation end) {
       var N = qap.Weights.Rows;
       var sol = start;
+      var fitness = QAPEvaluator.Apply(sol, qap.Weights, qap.Distances);
+      yield return Tuple.Create(sol, fitness);
+
       var invSol = GetInverse(sol);
       // we require at most N-1 steps to move from one permutation to another
@@ -197,7 +194,10 @@
     }
 
-    public IEnumerable<Tuple<Permutation, double>> FirstImprovementWalk(QuadraticAssignmentProblem qap, Permutation start, double fitness, Permutation end, IRandom random) {
+    public static IEnumerable<Tuple<Permutation, double>> FirstDirectedWalk(IRandom random, QuadraticAssignmentProblem qap, Permutation start, Permutation end) {
       var N = qap.Weights.Rows;
       var sol = start;
+      var fitness = QAPEvaluator.Apply(sol, qap.Weights, qap.Distances);
+      yield return Tuple.Create(sol, fitness);
+
       var invSol = GetInverse(sol);
       // randomize the order in which improvements are tried
@@ -229,18 +229,96 @@
     }
 
+    private static double Area(IEnumerable<Tuple<Permutation, double>> path) {
+      var iter = path.GetEnumerator();
+      if (!iter.MoveNext()) return 0.0;
+      var area = 0.0;
+      var prev = iter.Current;
+      while (iter.MoveNext()) {
+        area += TrapezoidArea(prev, iter.Current);
+        prev = iter.Current;
+      }
+      return area;
+    }
+
+    private static double TrapezoidArea(Tuple<Permutation, double> a, Tuple<Permutation, double> b) {
+      var area = 0.0;
+      var dist = Dist(a.Item1, b.Item1);
+      if ((a.Item2 <= 0 && b.Item2 <= 0) || (a.Item2 >= 0 && b.Item2 >= 0))
+        area += dist * (Math.Abs(a.Item2) + Math.Abs(b.Item2)) / 2.0;
+      else {
+        var k = (b.Item2 - a.Item2) / dist;
+        var d = a.Item2;
+        var x = -d / k;
+        area += Math.Abs(x * a.Item2 / 2.0);
+        area += Math.Abs((dist - x) * b.Item2 / 2.0);
+      }
+      return area;
+    }
+
+    private static double ComBelowZero(IEnumerable<Tuple<Permutation, double>> path) {
+      var area = 0.0;
+      var com = 0.0;
+      var nwalkDist = 0.0;
+      Tuple<Permutation, double> prev = null;
+      var iter = path.GetEnumerator();
+      while (iter.MoveNext()) {
+        var c = iter.Current;
+        if (prev != null) {
+          var ndist = Dist(prev.Item1, c.Item1) / (double)c.Item1.Length;
+          nwalkDist += ndist;
+          if (prev.Item2 < 0 || c.Item2 < 0) {
+            var a = TrapezoidArea(prev, c) / (double)c.Item1.Length;
+            area += a;
+            com += (nwalkDist - (ndist / 2.0)) * a;
+          }
+        }
+        prev = c;
+      }
+      return com / area;
+    }
+
+    private static IEnumerable<Tuple<Permutation, double>> ApproximateDerivative(IEnumerable<Tuple<Permutation, double>> data) {
+      Tuple<Permutation, double> prev = null, prev2 = null;
+      foreach (var d in data) {
+        if (prev == null) {
+          prev = d;
+          continue;
+        }
+        if (prev2 == null) {
+          prev2 = prev;
+          prev = d;
+          continue;
+        }
+        var dist = Dist(prev2.Item1, d.Item1);
+        yield return Tuple.Create(prev.Item1, (d.Item2 - prev2.Item2) / (double)dist);
+        prev2 = prev;
+        prev = d;
+      }
+    }
+
     private static double Dist(Permutation a, Permutation b) {
-      return a.Where((t, i) => t != b[i]).Count();
+      var dist = 0;
+      for (var i = 0; i < a.Length; i++)
+        if (a[i] != b[i]) dist++;
+      return dist;
     }
 
     private static int[] GetInverse(Permutation p) {
       var inv = new int[p.Length];
-      for (var i = 0; i < p.Length; i++) inv[p[i]] = i;
+      for (var i = 0; i < p.Length; i++) {
+        inv[p[i]] = i;
+      }
       return inv;
     }
 
-    private static void Rot1(Permutation p) {
-      var first = p[0];
-      for (var i = 0; i < p.Length - 1; i++) p[i] = p[i + 1];
-      p[p.Length - 1] = first;
+    // permutation must be strictly different in every position
+    private static void BiasedShuffle(Permutation p, IRandom random) {
+      for (var i = p.Length - 1; i > 0; i--) {
+        // Swap element "i" with a random earlier element (excluding itself)
+        var swapIndex = random.Next(i);
+        var h = p[swapIndex];
+        p[swapIndex] = p[i];
+        p[i] = h;
+      }
     }
   }
