source: branches/EfficientGlobalOptimization/HeuristicLab.Algorithms.EGO/InfillCriteria/RobustImprovement.cs @ 14818

#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.Parameters;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Problems.DataAnalysis;

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

  [StorableClass]
  [Item("RobustImprovementMeassure", "Adding or Subtracting the variance * factor to the model estimation")]
  public class RobustImprovement : InfillCriterionBase {

    #region ParameterNames
    private const string KParameterName = "NearestNeighbours";
    #endregion

    #region ParameterProperties
    public IFixedValueParameter<IntValue> KParameter => Parameters[KParameterName] as IFixedValueParameter<IntValue>;

    #endregion

    #region Properties
    private int K => KParameter.Value.Value;

    [Storable]
    private double MaxSolutionDist;

    [Storable]
    //TODO use VP-Tree instead of array
    private RealVector[] Data;
    #endregion

    #region HL-Constructors, Serialization and Cloning 
    [StorableConstructor]
    private RobustImprovement(bool deserializing) : base(deserializing) { }

    private RobustImprovement(RobustImprovement original, Cloner cloner) : base(original, cloner) {
      MaxSolutionDist = original.MaxSolutionDist;
      Data = original.Data != null ? original.Data.Select(cloner.Clone).ToArray() : null;
    }
    public RobustImprovement() {
      Parameters.Add(new FixedValueParameter<IntValue>(KParameterName, "A value larger than 0 indicating how many nearestNeighbours shall be used to determine the RI meassure", new IntValue(3)));
    }
    public override IDeepCloneable Clone(Cloner cloner) {
      return new RobustImprovement(this, cloner);
    }
    #endregion


    public override double Evaluate(RealVector vector) {
      List<RealVector> nearestNeighbours;
      List<double> distances;
      Search(vector, K, out nearestNeighbours, out distances);
      var distVectors = nearestNeighbours.Select(x => Minus(x, vector)).ToList();
      var sum = 0.0;
      var wsum = 1.0; //weights for angular distance
      var used = new HashSet<RealVector>();
      foreach (var distVector in distVectors) {
        var d = Math.Pow(distances[used.Count], 0.5);
        if (used.Count == 0) {
          sum += d;
        } else {
          var w = used.Select(x => Angular(distVector, x)).Min();
          sum += w * d;
          wsum += w;
        }
        used.Add(distVector);
      }
      sum /= wsum * MaxSolutionDist; //normalize
      return sum;
    }
    public override bool Maximization() {
      return ExpensiveMaximization;
    }
    protected override void Initialize() {
      var model = RegressionSolution.Model as IConfidenceRegressionModel;
      if (model == null) throw new ArgumentException("can not calculate EI without confidence measure");
      Data = new RealVector[RegressionSolution.ProblemData.Dataset.Rows];
      for (var i = 0; i < Data.Length; i++) {
        Data[i] = new RealVector(Encoding.Length);
        for (var j = 0; j < Encoding.Length; j++)
          Data[i][j] = RegressionSolution.ProblemData.Dataset.GetDoubleValue(i, j);
      }

      var maxSolution = new double[Encoding.Length];
      var minSolution = new double[Encoding.Length];
      for (var i = 0; i < Encoding.Length; i++) {
        var j = i % Encoding.Bounds.Rows;
        maxSolution[i] = Encoding.Bounds[j, 1];
        minSolution[i] = Encoding.Bounds[j, 0];
      }
      MaxSolutionDist = Euclidian(maxSolution, minSolution) / Data.Length;
    }

    #region Helpers
    private static double Euclidian(IEnumerable<double> a, IEnumerable<double> b) {
      return Math.Sqrt(a.Zip(b, (d, d1) => d - d1).Sum(d => d * d));
    }
    private static double Angular(RealVector a, RealVector b) {
      var innerProduct = a.Zip(b, (x, y) => x * y).Sum();
      var res = Math.Acos(innerProduct / (Norm(a) * Norm(b))) / Math.PI;
      return double.IsNaN(res) ? 0 : res;
    }
    private static double Norm(IEnumerable<double> a) {
      return Math.Sqrt(a.Sum(d => d * d));
    }
    private static RealVector Minus(RealVector a, RealVector b) {
      return new RealVector(a.Zip(b, (d, d1) => d - d1).ToArray());
    }

    private void Search(RealVector vector, int k, out List<RealVector> nearestNeighbours, out List<double> distances) {
      var neighbours = new SortedList<double, RealVector>(new DuplicateKeyComparer<double>());
      foreach (var n in Data) neighbours.Add(Euclidian(n, vector), n);
      nearestNeighbours = new List<RealVector>();

      distances = new List<double>();
      foreach (var entry in neighbours) {
        nearestNeighbours.Add(entry.Value);
        distances.Add(entry.Key);
        if (distances.Count == k) break;
      }
    }
    #endregion

    public class DuplicateKeyComparer<TKey> : IComparer<TKey> where TKey : IComparable {
      public int Compare(TKey x, TKey y) {
        var result = x.CompareTo(y);
        return result == 0 ? 1 : result;
      }
    }
  }
}