#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2010 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
 *
 * This file is part of HeuristicLab.
 *
 * HeuristicLab is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * HeuristicLab is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
 */
#endregion

using System;
using System.IO;
using System.Linq;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Operators;
using HeuristicLab.Optimization;
using HeuristicLab.Optimization.Operators;
using HeuristicLab.Parameters;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;

namespace HeuristicLab.Analysis.FitnessLandscape {

  [Item("PopulationDistributionAnalyzer", "An operator that analyzes the distribution of fitness values")]
  [StorableClass]
  public class PopulationDistributionAnalyzer : AlgorithmOperator, IAnalyzer {
    public bool EnabledByDefault {
      get { return false; }
    }

    #region Parameters
    public ScopeTreeLookupParameter<DoubleValue> QualityParameter {
      get { return (ScopeTreeLookupParameter<DoubleValue>)Parameters["Quality"]; }
    }
    public ValueLookupParameter<DataTable> FitnessQuantilesParameter {
      get { return (ValueLookupParameter<DataTable>)Parameters["Fitness Quantiles"]; }
    }
    public ValueLookupParameter<DataTable> PopulationDispersionParameter {
      get { return (ValueLookupParameter<DataTable>)Parameters["Population Dispersion"]; }
    }
    public ValueLookupParameter<DataTable> HigherPopulationMomentsParameter {
      get { return (ValueLookupParameter<DataTable>)Parameters["Higher Population Moments"]; }
    }
    public ValueLookupParameter<DataTable> PopulationNormalityParameter {
      get { return (ValueLookupParameter<DataTable>)Parameters["Population Normality"]; }
    }
    public ValueLookupParameter<VariableCollection> ResultsParameter {
      get { return (ValueLookupParameter<VariableCollection>)Parameters["Results"]; }
    }
    public ValueLookupParameter<ResultCollection> PopulationDistributionResultsParameter {
      get { return (ValueLookupParameter<ResultCollection>)Parameters["Population Distribution Results"]; }
    }
    public OptionalValueParameter<StringValue> PopulationLogFileNameParameter {
      get { return (OptionalValueParameter<StringValue>)Parameters["Population Log File Name"]; }
    }
    public IConstrainedValueParameter<IntValue> NQuantilesParameter {
      get { return (IConstrainedValueParameter<IntValue>)Parameters["NQuantiles"]; }
    }
    #endregion

    [StorableConstructor]
    protected PopulationDistributionAnalyzer(bool deserializing) : base(deserializing) { }
    protected PopulationDistributionAnalyzer(PopulationDistributionAnalyzer original, Cloner cloner) : base(original, cloner) { }

    public PopulationDistributionAnalyzer() {
      Parameters.Add(new ScopeTreeLookupParameter<DoubleValue>("Quality", "The quality of the solution"));
      Parameters.Add(new ValueLookupParameter<VariableCollection>("Results", "The collection of all results of this algorithm"));
      Parameters.Add(new ValueLookupParameter<ResultCollection>("Population Distribution Results", "All results from population distribution analysis"));
      Parameters.Add(new ValueLookupParameter<DataTable>("Fitness Quantiles", "Data table with quantiles of the fitness distribution"));
      Parameters.Add(new ValueLookupParameter<DataTable>("Population Dispersion", "Data table dispersion statistics"));
      Parameters.Add(new ValueLookupParameter<DataTable>("Higher Population Moments", "Data table skewness and kurtosis of population's fitness distribution"));
      Parameters.Add(new ValueLookupParameter<DataTable>("Population Normality", "Jarque-Bera Normality Test p-value and 0.05 threshold"));
      Parameters.Add(new OptionalValueParameter<StringValue>("Population Log File Name", "File name of a log file where all population fittness values are logged to"));
      Parameters.Add(new ConstrainedValueParameter<IntValue>("NQuantiles", "Number of quantiles to plot", new ItemSet<IntValue>(Enumerable.Range(1, 50).Select(v => new IntValue(v)))));

      NQuantilesParameter.Value = NQuantilesParameter.ValidValues.Single(v => v.Value == 10);

      var resultsCollector = new ResultsCollector();
      resultsCollector.ResultsParameter.ActualName = PopulationDistributionResultsParameter.Name;
      resultsCollector.CollectedValues.Add(new LookupParameter<DataTable>(FitnessQuantilesParameter.Name));
      resultsCollector.CollectedValues.Add(new LookupParameter<DataTable>(PopulationDispersionParameter.Name));
      resultsCollector.CollectedValues.Add(new LookupParameter<DataTable>(HigherPopulationMomentsParameter.Name));
      resultsCollector.CollectedValues.Add(new LookupParameter<DataTable>(PopulationNormalityParameter.Name));

      var globalResultsCollector = new ResultsCollector();
      globalResultsCollector.CollectedValues.Add(new ValueLookupParameter<ResultCollection>(PopulationDistributionResultsParameter.Name));

      OperatorGraph.InitialOperator = resultsCollector;
      resultsCollector.Successor = globalResultsCollector;
      globalResultsCollector.Successor = null;
    }

    public override IDeepCloneable Clone(Cloner cloner) {
      return new PopulationDistributionAnalyzer(this, cloner);
    }

    public override IOperation Apply() {
      if (PopulationDistributionResultsParameter.ActualValue == null)
        PopulationDistributionResultsParameter.ActualValue = new ResultCollection();
      var qualities = QualityParameter.ActualValue.Select(v => v.Value).ToArray();
      CalculateQuantiles(qualities);
      CalculateDistributionParameters(qualities);
      LogPopulationFitnessValues(qualities);
      return base.Apply();
    }

    private void CalculateQuantiles(double[] qualities) {
      DataTable quantiles = FitnessQuantilesParameter.ActualValue;
      if (quantiles == null) {
        quantiles = new DataTable("Fitness Quantiles");
        quantiles.Description = "The population's fitness quantiles";
        FitnessQuantilesParameter.ActualValue = quantiles;
        for (int i = 0; i <= NQuantilesParameter.Value.Value; i++)
          quantiles.Rows.Add(new DataRow((i * 10).ToString()));
      }
      int n_quantiles = quantiles.Rows.Count;
      for (int i = 0; i < n_quantiles; i++) {
        double v = 0;
        alglib.basestat.samplepercentile(qualities, qualities.Length, 1.0 * i / n_quantiles, ref v);
        quantiles.Rows[(i * 10).ToString()].Values.Add(v);
      }
    }

    private void CalculateDistributionParameters(double[] qualities) {
      DataTable populationDispersion = PopulationDispersionParameter.ActualValue;
      if (populationDispersion == null) {
        populationDispersion = new DataTable("Population Dispersion");
        PopulationDispersionParameter.ActualValue = populationDispersion;
        populationDispersion.Rows.Add(new DataRow("Std. Deviation"));
        populationDispersion.Rows.Add(new DataRow("Mean Difference"));
      }
      DataTable higherPopulationMoments = HigherPopulationMomentsParameter.ActualValue;
      if (higherPopulationMoments == null) {
        higherPopulationMoments = new DataTable("Higher Population Moments");
        HigherPopulationMomentsParameter.ActualValue = higherPopulationMoments;
        higherPopulationMoments.Rows.Add(new DataRow("Skewness"));
        higherPopulationMoments.Rows.Add(new DataRow("Kurtosis"));
      }
      DataTable populationNormality = PopulationNormalityParameter.ActualValue;
      if (populationNormality == null) {
        populationNormality = new DataTable("Population Normality");
        PopulationNormalityParameter.ActualValue = populationNormality;
        populationNormality.Rows.Add(new DataRow("Jarque-Bera P-Value"));
        populationNormality.Rows.Add(new DataRow("0.05"));
      }

      double mean, variance, skewness, kurtosis, p_value;
      mean = variance = skewness = kurtosis = p_value = 0;
      alglib.basestat.samplemoments(qualities, qualities.Length, ref mean, ref variance, ref skewness, ref kurtosis);
      alglib.jarquebera.jarqueberatest(qualities, qualities.Length, ref p_value);
      double mean_difference =
        (from i in Enumerable.Range(0, qualities.Length)
         from j in Enumerable.Range(0, i)
         select Math.Abs(qualities[i] - qualities[j])).Sum()
           * 2 / qualities.Length / (qualities.Length - 1);

      populationDispersion.Rows["Std. Deviation"].Values.Add(Math.Sqrt(variance));
      populationDispersion.Rows["Mean Difference"].Values.Add(mean_difference);
      higherPopulationMoments.Rows["Skewness"].Values.Add(skewness);
      higherPopulationMoments.Rows["Kurtosis"].Values.Add(kurtosis);
      populationNormality.Rows["Jarque-Bera P-Value"].Values.Add(p_value);
      populationNormality.Rows["0.05"].Values.Add(0.05);
    }

    private void LogPopulationFitnessValues(double[] qualities) {
      if (PopulationLogFileNameParameter.ActualValue == null)
        return;
      using (var writer = new StreamWriter(PopulationLogFileNameParameter.Value.Value, true)) {
        foreach (var q in qualities) {
          writer.Write(q);
          writer.Write(";");
        }
        writer.WriteLine();
        writer.Close();
      }
    }
  }
}