source: addons/HeuristicLab.FitnessLandscapeAnalysis/HeuristicLab.Analysis.FitnessLandscape/Analysis/RuggednessCalculator.cs

using System;
using System.Linq;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Operators;
using HeuristicLab.Parameters;
using System.Collections.Generic;
using HEAL.Attic;

namespace HeuristicLab.Analysis.FitnessLandscape.Analysis {
  [Item("Ruggedness Calculator", "Calculates ruggedness descriptors froma given quality trail.")]
  [StorableType("1739820A-DC44-4371-A912-E44C1DD46940")]
  public class RuggednessCalculator : SingleSuccessorOperator {

    #region Parameters
    public LookupParameter<DataTable> QualityTrailParameter {
      get { return (LookupParameter<DataTable>)Parameters["QualityTrail"]; }
    }
    public LookupParameter<IntValue> CorrelationLengthParameter {
      get { return (LookupParameter<IntValue>)Parameters["CorrelationLength"]; }
    }
    public LookupParameter<DoubleArray> AutoCorrelationParameter {
      get { return (LookupParameter<DoubleArray>)Parameters["AutoCorrelation"]; }
    }
    #endregion


    #region Constructors & Cloning
    [StorableConstructor]
    protected RuggednessCalculator(StorableConstructorFlag _) : base(_) { }
    protected RuggednessCalculator(RuggednessCalculator original, Cloner cloner) : base(original, cloner) { }
    public RuggednessCalculator() {
      Parameters.Add(new LookupParameter<DataTable>("QualityTrail", "Historical values of walk qualities"));
      Parameters.Add(new LookupParameter<IntValue>("CorrelationLength", "Average maximum distances between correlated quality values."));
      Parameters.Add(new LookupParameter<DoubleArray>("AutoCorrelation", "AutoCorrelation"));
    }
    public override IDeepCloneable Clone(Cloner cloner) {
      return new RuggednessCalculator(this, cloner);
    }
    #endregion

    public override IOperation Apply() {
      double[] qualities = QualityTrailParameter.ActualValue.Rows.First().Values.ToArray();
      double[] autocorrelation;
      CorrelationLengthParameter.ActualValue = new IntValue(CalculateCorrelationLength(qualities, out autocorrelation));
      AutoCorrelationParameter.ActualValue = new DoubleArray(autocorrelation);
      return base.Apply();
    }

    public static int CalculateCorrelationLength(double[] qualities, out double[] acf) {
      double[] correlations = new double[qualities.Length];
      alglib.corr.corrr1dcircular(qualities, qualities.Length, qualities, qualities.Length, ref correlations);
      double mean = 0;
      double variance = 0;
      double skewness = 0;
      double kurtosis = 0;
      alglib.basestat.samplemoments(qualities, qualities.Length, ref mean, ref variance, ref skewness, ref kurtosis);
      List<double> autocorrelation = new List<double>() { 1.0 };
      double bound = Math.Min(2 / Math.Sqrt(qualities.Length), 1.0);
      int correlationLength = 1;
      for (; correlationLength < qualities.Length; correlationLength++) {
        double value = correlations[correlationLength]/qualities.Length - mean * mean;
        if (variance > 0)
          value = Math.Max(Math.Min(value/variance, 1.0), -1.0);
        else
          value = 1;
        autocorrelation.Add(value);
        if (Math.Abs(value) < bound) break;
      }
      acf = autocorrelation.ToArray();
      return correlationLength-1;
    }

    public static bool AnyGreaterOne(double[] values) {
      return values.Any(d => d > 1);
    }
  }
}