source: branches/EfficientGlobalOptimization/HeuristicLab.Algorithms.EGO/EfficientGlobalOptimizationAlgorithm.cs @ 14741

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2016 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.Collections.Generic;
using System.Linq;
using System.Threading;
using HeuristicLab.Algorithms.CMAEvolutionStrategy;
using HeuristicLab.Algorithms.DataAnalysis;
using HeuristicLab.Algorithms.EGO;
using HeuristicLab.Analysis;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Encodings.RealVectorEncoding;
using HeuristicLab.Optimization;
using HeuristicLab.Parameters;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Problems.DataAnalysis;
using HeuristicLab.Random;

namespace HeuristicLab.Problems.SurrogateProblem {
  [StorableClass]
  [Creatable(CreatableAttribute.Categories.Algorithms, Priority = 95)]
  [Item("EfficientGlobalOptimizationAlgortihm", "Solves a problem by sequentially learning a model, solving a subproblem on the model and evaluating the best found solution for this subproblem.")]
  public class EfficientGlobalOptimizationAlgorithm : BasicAlgorithm {
    #region Basic-Alg-Essentials
    public override bool SupportsPause => true;
    public override Type ProblemType => typeof(SingleObjectiveBasicProblem<IEncoding>);
    public new SingleObjectiveBasicProblem<IEncoding> Problem
    {
      get { return (SingleObjectiveBasicProblem<IEncoding>)base.Problem; }
      set { base.Problem = value; }
    }
    #endregion

    #region ParameterNames
    private const string GenerationSizeParameterName = "GenerationSize";
    private const string InfillCriterionParameterName = "InfillCriterion";
    private const string InfillOptimizationAlgorithmParameterName = "InfillOptimizationAlgorithm";
    private const string InfillOptimizationRestartsParameterName = "InfillOptimizationRestarts";
    private const string InitialEvaluationsParameterName = "Initial Evaluations";
    private const string MaximumIterationsParameterName = "Maximum Iterations";
    private const string MaximumRuntimeParameterName = "Maximum Runtime";
    private const string RegressionAlgorithmParameterName = "RegressionAlgorithm";
    private const string SeedParameterName = "Seed";
    private const string SetSeedRandomlyParameterName = "SetSeedRandomly";
    #endregion

    #region ResultNames
    private const string BestQualityResultName = "Best Quality";
    private const string BestSolutionResultName = "Best Solution";
    private const string EvaluatedSoultionsResultName = "EvaluatedSolutions";
    private const string IterationsResultName = "Iterations";
    private const string RegressionSolutionResultName = "Model";
    private const string QualitiesChartResultName = "Qualities";
    private const string BestQualitiesRowResultName = "Best Quality";
    private const string CurrentQualitiesRowResultName = "Current Quality";
    private const string WorstQualitiesRowResultName = "Worst Quality";
    #endregion

    #region TransmissionResultNames
    public const string BestInfillSolutionResultName = "BestInfillSolution";
    public const string BestInfillQualityResultName = "BestInfillQuality";
    #endregion

    #region ParameterProperties
    public IFixedValueParameter<IntValue> GenerationSizeParemeter => Parameters[GenerationSizeParameterName] as IFixedValueParameter<IntValue>;
    public IConstrainedValueParameter<IInfillCriterion> InfillCriterionParameter => Parameters[InfillCriterionParameterName] as IConstrainedValueParameter<IInfillCriterion>;
    public IValueParameter<Algorithm> InfillOptimizationAlgorithmParameter => Parameters[InfillOptimizationAlgorithmParameterName] as IValueParameter<Algorithm>;
    public IFixedValueParameter<IntValue> InfillOptimizationRestartsParemeter => Parameters[InfillOptimizationRestartsParameterName] as IFixedValueParameter<IntValue>;
    public IFixedValueParameter<IntValue> InitialEvaluationsParameter => Parameters[InitialEvaluationsParameterName] as IFixedValueParameter<IntValue>;
    public IFixedValueParameter<IntValue> MaximumIterationsParameter => Parameters[MaximumIterationsParameterName] as IFixedValueParameter<IntValue>;
    public IFixedValueParameter<IntValue> MaximumRuntimeParameter => Parameters[MaximumRuntimeParameterName] as IFixedValueParameter<IntValue>;
    public IValueParameter<IDataAnalysisAlgorithm<IRegressionProblem>> RegressionAlgorithmParameter => Parameters[RegressionAlgorithmParameterName] as IValueParameter<IDataAnalysisAlgorithm<IRegressionProblem>>;
    public IFixedValueParameter<IntValue> SeedParameter => Parameters[SeedParameterName] as IFixedValueParameter<IntValue>;
    public FixedValueParameter<BoolValue> SetSeedRandomlyParameter => Parameters[SetSeedRandomlyParameterName] as FixedValueParameter<BoolValue>;
    #endregion

    #region Properties
    public int GenerationSize
    {
      get { return GenerationSizeParemeter.Value.Value; }
    }
    public IInfillCriterion InfillCriterion
    {
      get { return InfillCriterionParameter.Value; }
    }
    public Algorithm InfillOptimizationAlgorithm
    {
      get { return InfillOptimizationAlgorithmParameter.Value; }
    }
    public int InfillOptimizationRestarts
    {
      get { return InfillOptimizationRestartsParemeter.Value.Value; }
    }
    public int InitialEvaluations
    {
      get { return InitialEvaluationsParameter.Value.Value; }
    }
    public int MaximumIterations
    {
      get { return MaximumIterationsParameter.Value.Value; }
    }
    public int MaximumRuntime
    {
      get { return MaximumRuntimeParameter.Value.Value; }
    }
    public IDataAnalysisAlgorithm<IRegressionProblem> RegressionAlgorithm
    {
      get { return RegressionAlgorithmParameter.Value; }
    }
    public int Seed
    {
      get { return SeedParameter.Value.Value; }
    }
    public bool SetSeedRandomly
    {
      get { return SetSeedRandomlyParameter.Value.Value; }
    }
    #endregion

    #region StorableProperties
    [Storable]
    private IRandom Random = new MersenneTwister();
    [Storable]
    private List<Tuple<RealVector, double>> Samples;
    #endregion

    #region ResultsProperties
    private double ResultsBestQuality
    {
      get { return ((DoubleValue)Results[BestQualityResultName].Value).Value; }
      set { ((DoubleValue)Results[BestQualityResultName].Value).Value = value; }
    }
    private RealVector ResultsBestSolution
    {
      get { return (RealVector)Results[BestSolutionResultName].Value; }
      set { Results[BestSolutionResultName].Value = value; }
    }
    private int ResultsEvaluations
    {
      get { return ((IntValue)Results[EvaluatedSoultionsResultName].Value).Value; }
      set { ((IntValue)Results[EvaluatedSoultionsResultName].Value).Value = value; }
    }
    private int ResultsIterations
    {
      get { return ((IntValue)Results[IterationsResultName].Value).Value; }
      set { ((IntValue)Results[IterationsResultName].Value).Value = value; }
    }
    private DataTable ResultsQualities
    {
      get { return (DataTable)Results[QualitiesChartResultName].Value; }
    }
    private DataRow ResultsQualitiesBest
    {
      get { return ResultsQualities.Rows[BestQualitiesRowResultName]; }
    }
    private DataRow ResultsQualitiesWorst
    {
      get { return ResultsQualities.Rows[WorstQualitiesRowResultName]; }
    }
    private DataRow ResultsQualitiesIteration
    {
      get { return ResultsQualities.Rows[CurrentQualitiesRowResultName]; }
    }
    private IRegressionSolution ResultsModel
    {
      get { return (IRegressionSolution)Results[RegressionSolutionResultName].Value; }
      set { Results[RegressionSolutionResultName].Value = value; }
    }
    #endregion

    #region HLConstructors
    [StorableConstructor]
    protected EfficientGlobalOptimizationAlgorithm(bool deserializing) : base(deserializing) { }
    [StorableHook(HookType.AfterDeserialization)]
    private void AfterDeseialization() {
      RegisterEventhandlers();
    }
    protected EfficientGlobalOptimizationAlgorithm(EfficientGlobalOptimizationAlgorithm original, Cloner cloner)
      : base(original, cloner) {
      Random = cloner.Clone(Random);
      if (original.Samples != null) Samples = original.Samples.Select(x => new Tuple<RealVector, double>(cloner.Clone(x.Item1), x.Item2)).ToList();
      RegisterEventhandlers();
    }
    public override IDeepCloneable Clone(Cloner cloner) { return new EfficientGlobalOptimizationAlgorithm(this, cloner); }
    public EfficientGlobalOptimizationAlgorithm() {
      var cmaes = new CMAEvolutionStrategy {
        MaximumGenerations = 300,
        PopulationSize = 50
      };
      var model = new GaussianProcessRegression {
        Problem = new RegressionProblem()
      };
      model.CovarianceFunctionParameter.Value = new CovarianceRationalQuadraticIso();
      Parameters.Add(new FixedValueParameter<IntValue>(MaximumIterationsParameterName, "", new IntValue(int.MaxValue)));
      Parameters.Add(new FixedValueParameter<IntValue>(InitialEvaluationsParameterName, "", new IntValue(10)));
      Parameters.Add(new FixedValueParameter<IntValue>(MaximumRuntimeParameterName, "The maximum runtime in seconds after which the algorithm stops. Use -1 to specify no limit for the runtime", new IntValue(3600)));
      Parameters.Add(new FixedValueParameter<IntValue>(SeedParameterName, "The random seed used to initialize the new pseudo random number generator.", new IntValue(0)));
      Parameters.Add(new FixedValueParameter<BoolValue>(SetSeedRandomlyParameterName, "True if the random seed should be set to a random value, otherwise false.", new BoolValue(true)));
      Parameters.Add(new ValueParameter<IDataAnalysisAlgorithm<IRegressionProblem>>(RegressionAlgorithmParameterName, "The model used to approximate the problem", model));
      Parameters.Add(new ValueParameter<Algorithm>(InfillOptimizationAlgorithmParameterName, "The algorithm used to solve the expected improvement subproblem", cmaes));
      Parameters.Add(new FixedValueParameter<IntValue>(InfillOptimizationRestartsParameterName, "Number of restarts of the SubAlgortihm to avoid local optima", new IntValue(1)));
      Parameters.Add(new FixedValueParameter<IntValue>(GenerationSizeParameterName, "Number points that are sampled every iteration (stadard EGO: 1)", new IntValue(1)));
      Parameters.Add(new ConstrainedValueParameter<IInfillCriterion>(InfillCriterionParameterName, "Decision what value should decide the next sample"));
      InfillCriterionParameter.ValidValues.Add(new ExpectedImprovement());
      InfillCriterionParameter.ValidValues.Add(new ExpectedQuality());
      InfillCriterionParameter.ValidValues.Add(new ConfidenceBound());
      SetInfillProblem();
      RegisterEventhandlers();
    }
    #endregion

    protected override void Initialize(CancellationToken cancellationToken) {
      base.Initialize(cancellationToken);
      //encoding
      var enc = Problem.Encoding as RealVectorEncoding;
      if (enc == null) throw new ArgumentException("The EGO algorithm can only be applied to RealVectorEncodings");

      //random
      if (SetSeedRandomly) SeedParameter.Value.Value = new System.Random().Next();
      Random.Reset(Seed);
      Samples = new List<Tuple<RealVector, double>>();

      //results
      Results.Add(new Result(IterationsResultName, new IntValue(0)));
      Results.Add(new Result(EvaluatedSoultionsResultName, new IntValue(0)));
      Results.Add(new Result(BestSolutionResultName, new RealVector(1)));
      Results.Add(new Result(BestQualityResultName, new DoubleValue(Problem.Maximization ? double.MinValue : double.MaxValue)));
      Results.Add(new Result(RegressionSolutionResultName, typeof(IRegressionSolution)));
      var table = new DataTable(QualitiesChartResultName);
      table.Rows.Add(new DataRow(BestQualitiesRowResultName));
      table.Rows.Add(new DataRow(WorstQualitiesRowResultName));
      table.Rows.Add(new DataRow(CurrentQualitiesRowResultName));
      Results.Add(new Result(QualitiesChartResultName, table));

      ResultsQualities.Rows.Add(new DataRow("DEBUG:ModelBuildingIterations"));

      //initial samples
      var points = EgoUtilities.GetUniformRandomDesign(InitialEvaluations, enc.Length, enc.Bounds, Random);
      foreach (var t in points) {
        Samples.Add(Evaluate(t));
        cancellationToken.ThrowIfCancellationRequested();
      }

      Analyze();
    }

    protected override void Run(CancellationToken cancellationToken) {
      for (ResultsIterations = 0; ResultsIterations < MaximumIterations; ResultsIterations++) {
        try {
          ResultsModel = BuildModel();
          cancellationToken.ThrowIfCancellationRequested();
          for (var i = 0; i < GenerationSize; i++) {
            var samplepoint = OptimizeInfillProblem();
            var sample = Evaluate(samplepoint);
            Samples.Add(sample);
            cancellationToken.ThrowIfCancellationRequested();
          }

        }
        finally {
          Analyze();
        }
      }
    }

    #region Eventhandling
    private void RegisterEventhandlers() {
      DeregisterEventhandlers();
      RegressionAlgorithmParameter.ValueChanged += OnModelAlgorithmChanged;
      InfillOptimizationAlgorithmParameter.ValueChanged += OnInfillOptimizationAlgorithmChanged;
      InfillOptimizationAlgorithm.ProblemChanged += InfillOptimizationProblemChanged;
      InfillCriterionParameter.ValueChanged += InfillCriterionChanged;

    }
    private void DeregisterEventhandlers() {
      RegressionAlgorithmParameter.ValueChanged -= OnModelAlgorithmChanged;
      InfillOptimizationAlgorithmParameter.ValueChanged -= OnInfillOptimizationAlgorithmChanged;
      InfillOptimizationAlgorithm.ProblemChanged -= InfillOptimizationProblemChanged;
      InfillCriterionParameter.ValueChanged -= InfillCriterionChanged;
    }
    private void OnInfillOptimizationAlgorithmChanged(object sender, EventArgs args) {
      SetInfillProblem();
      InfillOptimizationAlgorithm.ProblemChanged -= InfillOptimizationProblemChanged;
      InfillOptimizationAlgorithm.ProblemChanged += InfillOptimizationProblemChanged;
    }
    private void InfillOptimizationProblemChanged(object sender, EventArgs e) {
      InfillOptimizationAlgorithm.ProblemChanged -= InfillOptimizationProblemChanged;
      SetInfillProblem();
      InfillOptimizationAlgorithm.ProblemChanged += InfillOptimizationProblemChanged;
    }
    private void InfillCriterionChanged(object sender, EventArgs e) {
      var infillProblem = InfillOptimizationAlgorithm.Problem as InfillProblem;
      if (infillProblem == null) throw new ArgumentException("InfillOptimizationAlgorithm has no InfillProblem. Troubles with Eventhandling?");
      infillProblem.InfillCriterion = InfillCriterion;
    }
    private void OnModelAlgorithmChanged(object sender, EventArgs args) {
      RegressionAlgorithm.Problem = new RegressionProblem();
    }
    protected override void OnExecutionTimeChanged() {
      base.OnExecutionTimeChanged();
      if (CancellationTokenSource == null) return;
      if (MaximumRuntime == -1) return;
      if (ExecutionTime.TotalSeconds > MaximumRuntime) CancellationTokenSource.Cancel();
    }
    protected override void OnPaused() {
      base.OnPaused();
      if (InfillOptimizationAlgorithm.ExecutionState == ExecutionState.Started) InfillOptimizationAlgorithm.Pause();
      if (RegressionAlgorithm.ExecutionState == ExecutionState.Started) RegressionAlgorithm.Pause();

    }
    protected override void OnStopped() {
      base.OnStopped();
      if (InfillOptimizationAlgorithm.ExecutionState != ExecutionState.Stopped) InfillOptimizationAlgorithm.Stop();
      if (RegressionAlgorithm.ExecutionState != ExecutionState.Stopped) RegressionAlgorithm.Stop();
    }
    protected override void OnProblemChanged() {
      base.OnProblemChanged();
      var infillProblem = InfillOptimizationAlgorithm.Problem as InfillProblem;
      if (infillProblem == null) throw new ArgumentException("InfillOptimizationAlgorithm has no InfillProblem. Troubles with Eventhandling?");
      infillProblem.Problem = Problem;
    }
    #endregion

    #region helpers
    private void SetInfillProblem() {
      var infillProblem = new InfillProblem {
        InfillCriterion = InfillCriterion,
        Problem = Problem
      };
      InfillOptimizationAlgorithm.Problem = infillProblem;
    }
    private IRegressionSolution BuildModel() {
      var dataset = EgoUtilities.GetDataSet(Samples);
      var problemdata = new RegressionProblemData(dataset, dataset.VariableNames.Where(x => !x.Equals("output")), "output");
      problemdata.TrainingPartition.Start = 0;
      problemdata.TrainingPartition.End = Samples.Count;
      problemdata.TestPartition.Start = Samples.Count;
      problemdata.TestPartition.End = Samples.Count;

      //train 
      var problem = (RegressionProblem)RegressionAlgorithm.Problem;
      problem.ProblemDataParameter.Value = problemdata;
      var i = 0;
      IRegressionSolution solution = null;
      double r2 = 0;
      while ((solution == null || RegressionAlgorithm is GaussianProcessRegression && r2 < 0.95) && i++ < 100) {  //TODO: ask why GP degenerates to NaN so often
        var results = EgoUtilities.SyncRunSubAlgorithm(RegressionAlgorithm, Random.Next(int.MaxValue));
        solution = results.Select(x => x.Value).OfType<IRegressionSolution>().SingleOrDefault();
        r2 = solution?.TrainingRSquared ?? 0;
      }
      ResultsQualities.Rows["DEBUG:ModelBuildingIterations"].Values.Add(i);
      if (solution == null) throw new ArgumentException("The Algorithm did not return a Model");
      return solution;
    }
    private RealVector OptimizeInfillProblem() {
      //parameterize and check InfillProblem
      var infillProblem = InfillOptimizationAlgorithm.Problem as InfillProblem;
      if (infillProblem == null) throw new ArgumentException("InfillOptimizationAlgorithm does not have InfillProblem. Problem with Eventhandling?");
      if (infillProblem.InfillCriterion != InfillCriterion) throw new ArgumentException("InfillCiriterion for Problem is not correct. Problem with Eventhandling?");
      if (infillProblem.Problem != Problem) throw new ArgumentException("Expensive real problem is not correctly set in InfillProblem. Problem with Eventhandling?");
      infillProblem.RegressionSolution = ResultsModel;

      RealVector bestVector = null;
      var bestValue = infillProblem.Maximization ? double.NegativeInfinity : double.PositiveInfinity;

      for (var i = 0; i < InfillOptimizationRestarts; i++) {
        //optimize
        var res = EgoUtilities.SyncRunSubAlgorithm(InfillOptimizationAlgorithm, Random.Next(int.MaxValue));

        //extract results
        if (!res.ContainsKey(BestInfillSolutionResultName)) throw new ArgumentException("The InfillOptimizationAlgorithm did not return a best solution");
        var v = res[BestInfillSolutionResultName].Value as RealVector;
        if (!res.ContainsKey(BestInfillQualityResultName)) throw new ArgumentException("The InfillOptimizationAlgorithm did not return a best quality");
        var d = res[BestInfillQualityResultName].Value as DoubleValue;
        if (d == null || v == null) throw new ArgumentException("The InfillOptimizationAlgorithm did not return the expected result types");

        //check for improvement
        if (infillProblem.Maximization != d.Value > bestValue) continue;
        bestValue = d.Value;
        bestVector = v;
      }

      InfillOptimizationAlgorithm.Runs.Clear();
      return bestVector;
    }
    private Tuple<RealVector, double> Evaluate(RealVector point) {
      var scope = new Scope();
      scope.Variables.Add(new Variable(Problem.Encoding.Name, point));
      var ind = new SingleEncodingIndividual(Problem.Encoding, scope);
      var q = Problem.Evaluate(ind, Random);
      return new Tuple<RealVector, double>(point, q);
    }
    private void Analyze() {
      ResultsEvaluations = Samples.Count;
      var max = Samples.ArgMax(x => x.Item2);
      var min = Samples.ArgMin(x => x.Item2);
      var best = Samples[Problem.Maximization ? max : min];
      ResultsBestQuality = best.Item2;
      ResultsBestSolution = best.Item1;
      ResultsQualitiesBest.Values.Add(ResultsBestQuality);
      ResultsQualitiesIteration.Values.Add(Samples[Samples.Count - 1].Item2);
      ResultsQualitiesWorst.Values.Add(Samples[Problem.Maximization ? min : max].Item2);
    }
    #endregion
  }
}