source: branches/EfficientGlobalOptimization/HeuristicLab.Algorithms.EGO/SamplingMethods/LatinHyperCubeDesign.cs @ 15870

#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 HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Encodings.RealVectorEncoding;
using HeuristicLab.Optimization;
using HeuristicLab.Parameters;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;

// ReSharper disable once CheckNamespace
namespace HeuristicLab.Algorithms.EGO {

  [StorableClass]
  [Item("LatinHyperCubeDesign", "A latin hypercube sampling strategy for real valued optimization")]
  public class LatinHyperCubeDesign : ParameterizedNamedItem, IInitialSampling {

    public const string DesigningAlgorithmParamterName = "DesignningAlgorithm";
    public IValueParameter<IAlgorithm> DesigningAlgorithmParameter => Parameters[DesigningAlgorithmParamterName] as IValueParameter<IAlgorithm>;
    private IAlgorithm DesigningAlgorithm => DesigningAlgorithmParameter.Value;

    [StorableConstructor]
    private LatinHyperCubeDesign(bool deserializing) : base(deserializing) { }
    private LatinHyperCubeDesign(LatinHyperCubeDesign original, Cloner cloner) : base(original, cloner) { }
    public LatinHyperCubeDesign() {
      var es = new EvolutionStrategy.EvolutionStrategy {
        PopulationSize = { Value = 50 },
        Children = { Value = 100 },
        MaximumGenerations = { Value = 300 }
      };
      Parameters.Add(new ValueParameter<IAlgorithm>(DesigningAlgorithmParamterName, "The Algorithm used for optimizing the Latin Hypercube (minimax) criterion.", es));
    }
    public override IDeepCloneable Clone(Cloner cloner) {
      return new LatinHyperCubeDesign(this, cloner);
    }

    public RealVector[] GetSamples(int noSamples, RealVector[] existingSamples, RealVectorEncoding encoding, IRandom random) {
      var lhsprob = new LHSProblem(noSamples, existingSamples, encoding);
      DesigningAlgorithm.Problem = lhsprob;
      EgoUtilities.SyncRunSubAlgorithm(DesigningAlgorithm, random.Next());
      return lhsprob.BestSolution;
    }
  }

  [StorableClass]
  internal class LHSProblem : SingleObjectiveBasicProblem<RealVectorEncoding> {
    [Storable]
    private int NoVectors;
    [Storable]
    private int NoDimensions;
    [Storable]
    private RealVector[] ExistingSamples;

    [Storable]
    public double BestQuality = double.MinValue;
    [Storable]
    public RealVector[] BestSolution;

    [StorableConstructor]
    private LHSProblem(bool deserializing) : base(deserializing) { }

    private LHSProblem(LHSProblem original, Cloner cloner) : base(original, cloner) {
      NoDimensions = original.NoDimensions;
      NoVectors = original.NoVectors;
      ExistingSamples = original.ExistingSamples?.Select(cloner.Clone).ToArray();
    }

    public LHSProblem(int novectors, RealVector[] existingSamples, RealVectorEncoding encoding) {
      NoVectors = novectors;
      NoDimensions = encoding.Length;
      ExistingSamples = existingSamples;
      Encoding.Length = novectors * encoding.Length;
      var b = new DoubleMatrix(Encoding.Length, 2);
      for (int i = 0; i < novectors; i++) {
        for (int j = 0; j < encoding.Length; j++) {
          var k = j % encoding.Bounds.Rows;
          b[i * encoding.Length + j, 0] = encoding.Bounds[k, 0];
          b[i * encoding.Length + j, 1] = encoding.Bounds[k, 1];
        }
      }
      Encoding.Bounds = b;

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

    public override double Evaluate(Individual individual, IRandom random) {
      return Enumerable.Range(0, NoVectors).Select(x => MinDistance(x, individual.RealVector())).Min();
    }

    //returns -1 if no distance can be calculated
    private double MinDistance(int pointNumber, RealVector design) {
      var min = double.MaxValue;
      if (ExistingSamples != null && ExistingSamples.Length != 0) min = ExistingSamples.Select(x => Distance(x, ExtractPoint(pointNumber, design))).Min();
      if (pointNumber == 0) return min;
      var min2 = Enumerable.Range(0, pointNumber).Select(i => Distance(ExtractPoint(i, design), ExtractPoint(pointNumber, design))).Min();
      return min < 0 ? min2 : Math.Min(min, min2);
    }
    private IEnumerable<double> ExtractPoint(int pointNumber, RealVector design) {
      return design.Skip(NoDimensions * pointNumber).Take(NoDimensions);
    }
    private static double Distance(IEnumerable<double> a, IEnumerable<double> b) {
      return a.Zip(b, (d, d1) => d - d1).Sum(d => d * d);
    }

    public override void Analyze(Individual[] individuals, double[] qualities, ResultCollection results, IRandom random) {
      base.Analyze(individuals, qualities, results, random);
      var i = qualities.ArgMin(x => x);
      if (BestQuality > qualities[i]) return;
      var r = individuals[i].RealVector();
      BestSolution = Enumerable.Range(0, NoVectors).Select(x => new RealVector(ExtractPoint(x, r).ToArray())).ToArray();
    }

    public override bool Maximization => true;
  }




}