source: branches/CMAES/HeuristicLab.Encodings.RealVectorEncoding/3.3/CMAESOperators/CMAUpdater.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("CMAUpdater", "Updates the covariance matrix and strategy parameters of CMA-ES.")]
  [StorableClass]
  public class CMAUpdater : SingleSuccessorOperator, IRealVectorOperator, ICMAESUpdater, IIterationBasedOperator {

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

    #region Parameter Properties

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

    public IScopeTreeLookupParameter<RealVector> MeansParameter {
      get { return (IScopeTreeLookupParameter<RealVector>)Parameters["Means"]; }
    }

    public IScopeTreeLookupParameter<RealVector> OffspringParameter {
      get { return (IScopeTreeLookupParameter<RealVector>)Parameters["Offspring"]; }
    }

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

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

    [StorableConstructor]
    protected CMAUpdater(bool deserializing) : base(deserializing) { }
    protected CMAUpdater(CMAUpdater original, Cloner cloner) : base(original, cloner) { }
    public CMAUpdater()
      : base() {
      Parameters.Add(new LookupParameter<CMAParameters>("StrategyParameters", "The strategy parameters of CMA-ES."));
      Parameters.Add(new ScopeTreeLookupParameter<RealVector>("Means", "The old and the new mean.", 2));
      Parameters.Add(new ScopeTreeLookupParameter<RealVector>("Offspring", "The created offspring solutions.", 3));
      Parameters.Add(new LookupParameter<IntValue>("Iterations", "The number of iterations passed."));
      Parameters.Add(new ValueLookupParameter<IntValue>("MaximumIterations", "The maximum number of iterations."));
      MeansParameter.ActualName = "RealVector";
      OffspringParameter.ActualName = "RealVector";
    }

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

    public override IOperation Apply() {
      var iterations = IterationsParameter.ActualValue.Value;

      var xold = MeansParameter.ActualValue[0];
      var xmean = MeansParameter.ActualValue[1];
      var offspring = OffspringParameter.ActualValue;

      var sp = StrategyParametersParameter.ActualValue;
      var N = sp.C.Rows;

      for (int i = 0; i < N; i++) {
        sp.BDz[i] = Math.Sqrt(sp.MuEff.Value) * (xmean[i] - xold[i]) / sp.Sigma.Value;
      }

      if (sp.InitialIterations.Value >= iterations) {
        for (int i = 0; i < N; i++) {
          sp.PS[i] = (1 - sp.CS.Value) * sp.PS[i]
                     + Math.Sqrt(sp.CS.Value * (2 - sp.CS.Value)) * sp.BDz[i] / sp.D[i];
        }
      } else {
        var artmp = new double[N];
        for (int i = 0; i < N; i++) {
          var sum = 0.0;
          for (int j = 0; j < N; j++) {
            sum += sp.B[j, i] * sp.BDz[j];
          }
          artmp[i] = sum / sp.D[i];
        }
        for (int i = 0; i < N; i++) {
          var sum = 0.0;
          for (int j = 0; j < N; j++) {
            sum += sp.B[i, j] * artmp[j];
          }
          sp.PS[i] = (1 - sp.CS.Value) * sp.PS[i] + Math.Sqrt(sp.CS.Value * (2 - sp.CS.Value)) * sum;
        }
      }
      var normPS = Math.Sqrt(sp.PS.Select(x => x * x).Sum());
      var hsig = normPS / Math.Sqrt(1 - Math.Pow(1 - sp.CS.Value, 2 * iterations)) / sp.ChiN.Value < 1.4 + 2.0 / (N + 1) ? 1.0 : 0.0;
      for (int i = 0; i < sp.PC.Length; i++) {
        sp.PC[i] = (1 - sp.CC.Value) * sp.PC[i]
                   + hsig * Math.Sqrt(sp.CC.Value * (2 - sp.CC.Value)) * sp.BDz[i];
      }

      if (sp.CCov.Value > 0) {
        if (sp.InitialIterations.Value >= iterations) {
          for (int i = 0; i < N; i++) {
            sp.C[i, i] = (1 - sp.CCovSep.Value) * sp.C[i, i]
                         + sp.CCov.Value * (1 / sp.MuCov.Value)
                         * (sp.PC[i] * sp.PC[i] + (1 - hsig) * sp.CC.Value * (2 - sp.CC.Value) * sp.C[i, i]);
            for (int k = 0; k < sp.Mu.Value; k++) {
              sp.C[i, i] += sp.CCov.Value * (1 - 1 / sp.MuCov.Value) * sp.Weights[k] * (offspring[k][i] - xold[i]) *
                            (offspring[k][i] - xold[i]) / (sp.Sigma.Value * sp.Sigma.Value);
            }
          }
        } else {
          for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
              sp.C[i, j] = (1 - sp.CCov.Value) * sp.C[i, j]
                           + sp.CCov.Value * (1 / sp.MuCov.Value)
                           * (sp.PC[i] * sp.PC[j] + (1 - hsig) * sp.CC.Value * (2 - sp.CC.Value) * sp.C[i, j]);
              for (int k = 0; k < sp.Mu.Value; k++) {
                sp.C[i, j] += sp.CCov.Value * (1 - 1 / sp.MuCov.Value) * sp.Weights[k] * (offspring[k][i] - xold[i]) *
                              (offspring[k][j] - xold[j]) / (sp.Sigma.Value * sp.Sigma.Value);
              }
            }
          }
        }
      }
      sp.Sigma.Value *= Math.Exp((sp.CS.Value / sp.Damps.Value) * (normPS / sp.ChiN.Value - 1));

      double minSqrtdiagC = int.MaxValue, maxSqrtdiagC = int.MinValue;
      for (int i = 0; i < N; i++) {
        if (Math.Sqrt(sp.C[i, i]) < minSqrtdiagC) minSqrtdiagC = Math.Sqrt(sp.C[i, i]);
        if (Math.Sqrt(sp.C[i, i]) > maxSqrtdiagC) maxSqrtdiagC = Math.Sqrt(sp.C[i, i]);
      }

      // ensure maximal and minimal standard deviations
      if (sp.MinSigma != null && sp.MinSigma.Length > 0) {
        for (int i = 0; i < N; i++) {
          var d = sp.MinSigma[Math.Min(i, sp.MinSigma.Length - 1)];
          if (d > sp.Sigma.Value * minSqrtdiagC) sp.Sigma.Value = d / minSqrtdiagC;
        }
      }
      if (sp.MaxSigma != null && sp.MaxSigma.Length > 0) {
        for (int i = 0; i < N; i++) {
          var d = sp.MaxSigma[Math.Min(i, sp.MaxSigma.Length - 1)];
          if (d > sp.Sigma.Value * maxSqrtdiagC) sp.Sigma.Value = d / maxSqrtdiagC;
        }
      }
      // end ensure ...

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

      if (fac != 1.0) {
        sp.Sigma.Value /= fac;
        for (int i = 0; i < N; i++) {
          sp.PC[i] *= fac;
          sp.D[i] *= fac;
          for (int j = 0; j < N; j++)
            sp.C[i, j] *= fac * fac;
        }
      }
      // end testAndCorrectNumerics

      if (sp.InitialIterations.Value >= iterations) {
        for (int i = 0; i < N; i++)
          sp.D[i] = Math.Sqrt(sp.C[i, i]);
      } else {
        var c = new double[sp.C.Rows, sp.C.Columns];
        for (int i = 0; i < sp.C.Rows; i++)
          for (int j = 0; j < sp.C.Columns; j++)
            c[i, j] = sp.C[i, j];
        double[] eVal;
        double[,] eVec;
        if (!alglib.smatrixevd(c, N, 1, false, out eVal, out eVec))
          throw new InvalidOperationException("Unable to perform eigendecomposition of matrix.");
        sp.B = new DoubleMatrix(eVec);

        for (int i = 0; i < sp.D.Length; i++)
          sp.D[i] = Math.Sqrt(eVal[i]);

        if (sp.D.Min() == 0.0) sp.AxisRatio.Value = double.PositiveInfinity;
        else sp.AxisRatio.Value = sp.D.Max() / sp.D.Min();
      }
      return base.Apply();
    }
  }
}