source: branches/1614_GeneralizedQAP/HeuristicLab.Analysis.FitnessLandscape/3.3/ProblemCharacteristicAnalysis/CurveAnalysis.cs @ 15713

using System;
using System.Collections.Generic;
using System.Linq;

namespace HeuristicLab.Analysis.FitnessLandscape {
  public static class CurveAnalysis<T> {
    public static CurveAnalysisResult GetCharacteristics(List<List<Tuple<T, double>>> trajectories, Func<T, T, double> distFunc) {
      trajectories = trajectories.Where(x => x.Count > 5).ToList();
      var f1 = trajectories.Select(path => ApproximateDerivative(path, distFunc).ToList()).ToList();
      var f2 = f1.Select(d1 => ApproximateDerivative(d1, distFunc).ToList()).ToList();

      var sharpness = f1.Average(x => x.Average(y => Math.Abs(y.Item2)));//f2.Average(x => Area(x, distFunc));
      var bumpiness = 0.0;
      var flatness = 0.0;

      for (var p = 0; p < f2.Count; p++) {
        if (f2[p].Count <= 2) continue;
        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++;
          }
        }
        bumpiness += bump / (f2[p].Count - 1.0);
        flatness += flat / (f2[p].Count - 1.0);
      }
      bumpiness /= f2.Count;
      flatness /= f2.Count;
      return new CurveAnalysisResult(sharpness, bumpiness, flatness);
    }

    private static double Area(IEnumerable<Tuple<T, double>> path, Func<T, T, double> distFunc) {
      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, distFunc);
        prev = iter.Current;
      }
      return area;
    }

    private static double TrapezoidArea(Tuple<T, double> a, Tuple<T, double> b, Func<T, T, double> distFunc) {
      var area = 0.0;
      var dist = distFunc(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 IEnumerable<Tuple<T, double>> ApproximateDerivative(IEnumerable<Tuple<T, double>> data, Func<T, T, double> distFunc) {
      Tuple<T, 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 = distFunc(prev2.Item1, d.Item1);
        yield return Tuple.Create(prev.Item1, (d.Item2 - prev2.Item2) / dist);
        prev2 = prev;
        prev = d;
      }
    }
  }

  public class CurveAnalysisResult {
    public double Sharpness { get; private set; }
    public double Bumpiness { get; private set; }
    public double Flatness { get; private set; }

    public string[] GetNames() {
      return new[] { "Sharpness", "Bumpiness", "Flatness" };
    }

    public double[] GetValues() {
      return new[] { Sharpness, Bumpiness, Flatness };
    }

    public CurveAnalysisResult(double sharpness, double bumpiness, double flatness) {
      Sharpness = sharpness;
      Bumpiness = bumpiness;
      Flatness = flatness;
    }
  }
}