source: branches/CMAES/HeuristicLab.Encodings.RealVectorEncoding/3.3/CMAESOperators/CMAInitializer.cs @ 9118

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2013 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 HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Operators;
using HeuristicLab.Optimization;
using HeuristicLab.Parameters;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using System;
using System.Linq;

namespace HeuristicLab.Encodings.RealVectorEncoding {
  [Item("CMAInitializer", "Initializes the covariance matrix and step size variables.")]
  [StorableClass]
  public class CMAInitializer : SingleSuccessorOperator, IRealVectorOperator, ICMAESInitializer, IIterationBasedOperator {

    public Type CMAType {
      get { return typeof(CMAParameters); }
    }

    #region Parameter Properties
    public IValueLookupParameter<IntValue> DimensionParameter {
      get { return (IValueLookupParameter<IntValue>)Parameters["Dimension"]; }
    }

    public IValueLookupParameter<DoubleArray> InitialSigmaParameter {
      get { return (IValueLookupParameter<DoubleArray>)Parameters["InitialSigma"]; }
    }

    public IValueLookupParameter<DoubleArray> MaxSigmaParameter {
      get { return (IValueLookupParameter<DoubleArray>)Parameters["MaxSigma"]; }
    }

    public IValueLookupParameter<DoubleArray> MinSigmaParameter {
      get { return (IValueLookupParameter<DoubleArray>)Parameters["MinSigma"]; }
    }

    public ILookupParameter<IntValue> IterationsParameter {
      get { return (ILookupParameter<IntValue>)Parameters["Iterations"]; }
    }

    public IValueLookupParameter<IntValue> MaximumIterationsParameter {
      get { return (IValueLookupParameter<IntValue>)Parameters["MaximumIterations"]; }
    }

    public IValueLookupParameter<IntValue> InitialIterationsParameter {
      get { return (IValueLookupParameter<IntValue>)Parameters["InitialIterations"]; }
    }

    public ILookupParameter<IntValue> PopulationSizeParameter {
      get { return (ILookupParameter<IntValue>)Parameters["PopulationSize"]; }
    }

    public ILookupParameter<IntValue> MuParameter {
      get { return (ILookupParameter<IntValue>)Parameters["Mu"]; }
    }

    public ILookupParameter<CMAParameters> StrategyParametersParameter {
      get { return (ILookupParameter<CMAParameters>)Parameters["StrategyParameters"]; }
    }
    #endregion

    [StorableConstructor]
    protected CMAInitializer(bool deserializing) : base(deserializing) { }
    protected CMAInitializer(CMAInitializer original, Cloner cloner) : base(original, cloner) { }
    public CMAInitializer()
      : base() {
      Parameters.Add(new ValueLookupParameter<IntValue>("Dimension", "The problem dimension (N)."));
      Parameters.Add(new ValueLookupParameter<DoubleArray>("InitialSigma", "The initial value for Sigma (need to be > 0), can be single dimensioned or an array that should be equal to the size of the vector.", new DoubleArray(new[] { 0.5 })));
      Parameters.Add(new ValueLookupParameter<DoubleArray>("MaxSigma", "The maximum sigma value can be omitted, given as one value for all dimensions or a value for each dimension."));
      Parameters.Add(new ValueLookupParameter<DoubleArray>("MinSigma", "The minimum sigma value can be omitted, given as one value for all dimensions or a value for each dimension."));
      Parameters.Add(new LookupParameter<IntValue>("Iterations", "The current iteration that is being processed."));
      Parameters.Add(new ValueLookupParameter<IntValue>("MaximumIterations", "The maximum number of iterations to be processed."));
      Parameters.Add(new ValueLookupParameter<IntValue>("InitialIterations", "The number of iterations that should be performed using the diagonal covariance matrix only."));
      Parameters.Add(new LookupParameter<IntValue>("PopulationSize", "The population size (lambda)."));
      Parameters.Add(new LookupParameter<IntValue>("Mu", "The number of offspring considered for updating of the strategy parameters."));
      Parameters.Add(new LookupParameter<CMAParameters>("StrategyParameters", "The strategy parameters for real-encoded CMA-ES."));
    }

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

    public override IOperation Apply() {
      var N = DimensionParameter.ActualValue.Value;
      var lambda = PopulationSizeParameter.ActualValue.Value;

      var sp = new CMAParameters(N, lambda);

      var s = InitialSigmaParameter.ActualValue;
      if (s == null || s.Length == 0) throw new InvalidOperationException("Initial standard deviation (sigma) must be given.");
      var sigma = s.Max();
      if (sigma <= 0) throw new InvalidOperationException("Initial standard deviation (sigma) must be > 0.");

      var pc = new double[N]; // evolution paths for C
      var ps = new double[N]; // evolution paths for sigma
      var B = new double[N, N]; // B defines the coordinate system
      var D = new double[N]; // diagonal D defines the scaling
      var C = new double[N, N]; // covariance matrix C
      var BDz = new double[N];
      double minSqrtdiagC = int.MaxValue, maxSqrtdiagC = int.MinValue;
      for (int i = 0; i < N; i++) {
        B[i, i] = 1;
        if (s.Length == 1) D[i] = 1;
        else if (s.Length == N) D[i] = s[i] / sigma;
        else throw new InvalidOperationException("Initial standard deviation (sigma) must either contain only one value for all dimension or for every dimension.");
        if (D[i] <= 0) throw new InvalidOperationException("Initial standard deviation (sigma) values must all be > 0.");
        C[i, i] = D[i] * D[i];
        if (Math.Sqrt(C[i, i]) < minSqrtdiagC) minSqrtdiagC = Math.Sqrt(C[i, i]);
        if (Math.Sqrt(C[i, i]) > maxSqrtdiagC) maxSqrtdiagC = Math.Sqrt(C[i, i]);
      }

      // ensure maximal and minimal standard deviations
      var minSigma = MinSigmaParameter.ActualValue;
      if (minSigma != null && minSigma.Length > 0) {
        for (int i = 0; i < N; i++) {
          var d = minSigma[(int)Math.Min(i, minSigma.Length - 1)];
          if (d > sigma * minSqrtdiagC) sigma = d / minSqrtdiagC;
        }
      }
      var maxSigma = MaxSigmaParameter.ActualValue;
      if (maxSigma != null && maxSigma.Length > 0) {
        for (int i = 0; i < N; i++) {
          var d = maxSigma[(int)Math.Min(i, maxSigma.Length - 1)];
          if (d > sigma * maxSqrtdiagC) sigma = d / maxSqrtdiagC;
        }
      }
      // end ensure ...

      // testAndCorrectNumerics
      double fac = 1;
      if (D.Max() < 1e-6)
        fac = 1.0 / D.Max();
      else if (D.Min() > 1e4)
        fac = 1.0 / D.Min();

      if (fac != 1.0) {
        sigma /= fac;
        for (int i = 0; i < N; i++) {
          pc[i] *= fac;
          D[i] *= fac;
          for (int j = 0; j < N; j++)
            C[i, j] *= fac * fac;
        }
      }
      // end testAndCorrectNumerics

      var initialIterations = InitialIterationsParameter.ActualValue;
      if (initialIterations == null) {
        initialIterations = new IntValue(0);
      }

      double maxD = D.Max(), minD = D.Min();
      if (minD == 0) sp.AxisRatio = new DoubleValue(double.PositiveInfinity);
      else sp.AxisRatio = new DoubleValue(maxD / minD);
      sp.PC = new DoubleArray(pc);
      sp.PS = new DoubleArray(ps);
      sp.B = new DoubleMatrix(B);
      sp.D = new DoubleArray(D);
      sp.C = new DoubleMatrix(C);
      sp.BDz = new DoubleArray(BDz);
      sp.Sigma = new DoubleValue(sigma);
      if (maxSigma != null) sp.MaxSigma = (DoubleArray)maxSigma.Clone();
      if (minSigma != null) sp.MinSigma = (DoubleArray)minSigma.Clone();
      sp.InitialIterations = new IntValue(initialIterations.Value);

      MuParameter.ActualValue = new IntValue(sp.Mu.Value);
      StrategyParametersParameter.ActualValue = sp;
      return base.Apply();
    }
  }
}