1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022015 Heuristic and Evolutionary Algorithms Laboratory (HEAL)


4  *


5  * This file is part of HeuristicLab.


6  *


7  * HeuristicLab is free software: you can redistribute it and/or modify


8  * it under the terms of the GNU General Public License as published by


9  * the Free Software Foundation, either version 3 of the License, or


10  * (at your option) any later version.


11  *


12  * HeuristicLab is distributed in the hope that it will be useful,


13  * but WITHOUT ANY WARRANTY; without even the implied warranty of


14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


15  * GNU General Public License for more details.


16  *


17  * You should have received a copy of the GNU General Public License


18  * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.


19  */


20  #endregion


21 


22  using HeuristicLab.Common;


23  using HeuristicLab.Core;


24  using HeuristicLab.Data;


25  using HeuristicLab.Encodings.RealVectorEncoding;


26  using HeuristicLab.Operators;


27  using HeuristicLab.Optimization;


28  using HeuristicLab.Parameters;


29  using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;


30  using System;


31  using System.Linq;


32 


33  namespace HeuristicLab.Algorithms.CMAEvolutionStrategy {


34  [Item("CMAUpdater", "Updates the covariance matrix and strategy parameters of CMAES.")]


35  [StorableClass]


36  public class CMAUpdater : SingleSuccessorOperator, ICMAUpdater, IIterationBasedOperator, ISingleObjectiveOperator {


37 


38  public Type CMAType {


39  get { return typeof(CMAParameters); }


40  }


41 


42  #region Parameter Properties


43  public ILookupParameter<CMAParameters> StrategyParametersParameter {


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


45  }


46 


47  public ILookupParameter<RealVector> MeanParameter {


48  get { return (ILookupParameter<RealVector>)Parameters["Mean"]; }


49  }


50 


51  public ILookupParameter<RealVector> OldMeanParameter {


52  get { return (ILookupParameter<RealVector>)Parameters["OldMean"]; }


53  }


54 


55  public IScopeTreeLookupParameter<RealVector> OffspringParameter {


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


57  }


58 


59  public IScopeTreeLookupParameter<DoubleValue> QualityParameter {


60  get { return (IScopeTreeLookupParameter<DoubleValue>)Parameters["Quality"]; }


61  }


62 


63  public ILookupParameter<IntValue> IterationsParameter {


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


65  }


66 


67  public IValueLookupParameter<IntValue> MaximumIterationsParameter {


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


69  }


70 


71  public IValueLookupParameter<IntValue> MaximumEvaluatedSolutionsParameter {


72  get { return (IValueLookupParameter<IntValue>)Parameters["MaximumEvaluatedSolutions"]; }


73  }


74 


75  public ILookupParameter<BoolValue> DegenerateStateParameter {


76  get { return (ILookupParameter<BoolValue>)Parameters["DegenerateState"]; }


77  }


78  #endregion


79 


80  [StorableConstructor]


81  protected CMAUpdater(bool deserializing) : base(deserializing) { }


82  protected CMAUpdater(CMAUpdater original, Cloner cloner) : base(original, cloner) { }


83  public CMAUpdater()


84  : base() {


85  Parameters.Add(new LookupParameter<CMAParameters>("StrategyParameters", "The strategy parameters of CMAES."));


86  Parameters.Add(new LookupParameter<RealVector>("Mean", "The new mean."));


87  Parameters.Add(new LookupParameter<RealVector>("OldMean", "The old mean."));


88  Parameters.Add(new ScopeTreeLookupParameter<RealVector>("Offspring", "The created offspring solutions."));


89  Parameters.Add(new ScopeTreeLookupParameter<DoubleValue>("Quality", "The quality of the offspring."));


90  Parameters.Add(new LookupParameter<IntValue>("Iterations", "The number of iterations passed."));


91  Parameters.Add(new ValueLookupParameter<IntValue>("MaximumIterations", "The maximum number of iterations."));


92  Parameters.Add(new ValueLookupParameter<IntValue>("MaximumEvaluatedSolutions", "The maximum number of evaluated solutions."));


93  Parameters.Add(new LookupParameter<BoolValue>("DegenerateState", "Whether the algorithm state has degenerated and should be terminated."));


94  MeanParameter.ActualName = "XMean";


95  OldMeanParameter.ActualName = "XOld";


96  OffspringParameter.ActualName = "RealVector";


97  }


98 


99  public override IDeepCloneable Clone(Cloner cloner) {


100  return new CMAUpdater(this, cloner);


101  }


102 


103  public override IOperation Apply() {


104  var iterations = IterationsParameter.ActualValue.Value;


105 


106  var xold = OldMeanParameter.ActualValue;


107  var xmean = MeanParameter.ActualValue;


108  var offspring = OffspringParameter.ActualValue;


109  var quality = QualityParameter.ActualValue;


110  var lambda = offspring.Length;


111 


112  var N = xmean.Length;


113  var sp = StrategyParametersParameter.ActualValue;


114 


115  #region Initialize default values for strategy parameter adjustment


116  if (sp.ChiN == 0) sp.ChiN = Math.Sqrt(N) * (1.0  1.0 / (4.0 * N) + 1.0 / (21.0 * N * N));


117  if (sp.MuEff == 0) sp.MuEff = sp.Weights.Sum() * sp.Weights.Sum() / sp.Weights.Sum(x => x * x);


118  if (sp.CS == 0) sp.CS = (sp.MuEff + 2) / (N + sp.MuEff + 3);


119  if (sp.Damps == 0) {


120  var maxIterations = MaximumIterationsParameter.ActualValue.Value;


121  var maxEvals = MaximumEvaluatedSolutionsParameter.ActualValue.Value;


122  sp.Damps = 2 * Math.Max(0, Math.Sqrt((sp.MuEff  1) / (N + 1))  1)


123  * Math.Max(0.3, 1  N / (1e6 + Math.Min(maxIterations, maxEvals / lambda))) + sp.CS + 1;


124  }


125  if (sp.CC == 0) sp.CC = 4.0 / (N + 4);


126  if (sp.MuCov == 0) sp.MuCov = sp.MuEff;


127  if (sp.CCov == 0) sp.CCov = 2.0 / ((N + 1.41) * (N + 1.41) * sp.MuCov)


128  + (1  (1.0 / sp.MuCov)) * Math.Min(1, (2 * sp.MuEff  1) / (sp.MuEff + (N + 2) * (N + 2)));


129  if (sp.CCovSep == 0) sp.CCovSep = Math.Min(1, sp.CCov * (N + 1.5) / 3);


130  #endregion


131 


132  sp.QualityHistory.Enqueue(quality[0].Value);


133  while (sp.QualityHistory.Count > sp.QualityHistorySize && sp.QualityHistorySize >= 0)


134  sp.QualityHistory.Dequeue();


135 


136  for (int i = 0; i < N; i++) {


137  sp.BDz[i] = Math.Sqrt(sp.MuEff) * (xmean[i]  xold[i]) / sp.Sigma;


138  }


139 


140  if (sp.InitialIterations >= iterations) {


141  for (int i = 0; i < N; i++) {


142  sp.PS[i] = (1  sp.CS) * sp.PS[i]


143  + Math.Sqrt(sp.CS * (2  sp.CS)) * sp.BDz[i] / sp.D[i];


144  }


145  } else {


146  var artmp = new double[N];


147  for (int i = 0; i < N; i++) {


148  var sum = 0.0;


149  for (int j = 0; j < N; j++) {


150  sum += sp.B[j, i] * sp.BDz[j];


151  }


152  artmp[i] = sum / sp.D[i];


153  }


154  for (int i = 0; i < N; i++) {


155  var sum = 0.0;


156  for (int j = 0; j < N; j++) {


157  sum += sp.B[i, j] * artmp[j];


158  }


159  sp.PS[i] = (1  sp.CS) * sp.PS[i] + Math.Sqrt(sp.CS * (2  sp.CS)) * sum;


160  }


161  }


162  var normPS = Math.Sqrt(sp.PS.Select(x => x * x).Sum());


163  var hsig = normPS / Math.Sqrt(1  Math.Pow(1  sp.CS, 2 * iterations)) / sp.ChiN < 1.4 + 2.0 / (N + 1) ? 1.0 : 0.0;


164  for (int i = 0; i < sp.PC.Length; i++) {


165  sp.PC[i] = (1  sp.CC) * sp.PC[i]


166  + hsig * Math.Sqrt(sp.CC * (2  sp.CC)) * sp.BDz[i];


167  }


168 


169  if (sp.CCov > 0) {


170  if (sp.InitialIterations >= iterations) {


171  for (int i = 0; i < N; i++) {


172  sp.C[i, i] = (1  sp.CCovSep) * sp.C[i, i]


173  + sp.CCov * (1 / sp.MuCov)


174  * (sp.PC[i] * sp.PC[i] + (1  hsig) * sp.CC * (2  sp.CC) * sp.C[i, i]);


175  for (int k = 0; k < sp.Mu; k++) {


176  sp.C[i, i] += sp.CCov * (1  1 / sp.MuCov) * sp.Weights[k] * (offspring[k][i]  xold[i]) *


177  (offspring[k][i]  xold[i]) / (sp.Sigma * sp.Sigma);


178  }


179  }


180  } else {


181  for (int i = 0; i < N; i++) {


182  for (int j = 0; j < N; j++) {


183  sp.C[i, j] = (1  sp.CCov) * sp.C[i, j]


184  + sp.CCov * (1 / sp.MuCov)


185  * (sp.PC[i] * sp.PC[j] + (1  hsig) * sp.CC * (2  sp.CC) * sp.C[i, j]);


186  for (int k = 0; k < sp.Mu; k++) {


187  sp.C[i, j] += sp.CCov * (1  1 / sp.MuCov) * sp.Weights[k] * (offspring[k][i]  xold[i]) *


188  (offspring[k][j]  xold[j]) / (sp.Sigma * sp.Sigma);


189  }


190  }


191  }


192  }


193  }


194  sp.Sigma *= Math.Exp((sp.CS / sp.Damps) * (normPS / sp.ChiN  1));


195 


196  double minSqrtdiagC = int.MaxValue, maxSqrtdiagC = int.MinValue;


197  for (int i = 0; i < N; i++) {


198  if (Math.Sqrt(sp.C[i, i]) < minSqrtdiagC) minSqrtdiagC = Math.Sqrt(sp.C[i, i]);


199  if (Math.Sqrt(sp.C[i, i]) > maxSqrtdiagC) maxSqrtdiagC = Math.Sqrt(sp.C[i, i]);


200  }


201 


202  // ensure maximal and minimal standard deviations


203  if (sp.SigmaBounds != null && sp.SigmaBounds.GetLength(0) > 0) {


204  for (int i = 0; i < N; i++) {


205  var d = sp.SigmaBounds[Math.Min(i, sp.SigmaBounds.GetLength(0)  1), 0];


206  if (d > sp.Sigma * minSqrtdiagC) sp.Sigma = d / minSqrtdiagC;


207  }


208  for (int i = 0; i < N; i++) {


209  var d = sp.SigmaBounds[Math.Min(i, sp.SigmaBounds.GetLength(0)  1), 1];


210  if (d > sp.Sigma * maxSqrtdiagC) sp.Sigma = d / maxSqrtdiagC;


211  }


212  }


213  // end ensure ...


214 


215  // testAndCorrectNumerics


216  double fac = 1;


217  if (sp.D.Max() < 1e6)


218  fac = 1.0 / sp.D.Max();


219  else if (sp.D.Min() > 1e4)


220  fac = 1.0 / sp.D.Min();


221 


222  if (fac != 1.0) {


223  sp.Sigma /= fac;


224  for (int i = 0; i < N; i++) {


225  sp.PC[i] *= fac;


226  sp.D[i] *= fac;


227  for (int j = 0; j < N; j++)


228  sp.C[i, j] *= fac * fac;


229  }


230  }


231  // end testAndCorrectNumerics


232 


233 


234  if (sp.InitialIterations >= iterations) {


235  for (int i = 0; i < N; i++)


236  sp.D[i] = Math.Sqrt(sp.C[i, i]);


237  DegenerateStateParameter.ActualValue = new BoolValue(false);


238  } else {


239 


240  // set B < C


241  for (int i = 0; i < N; i++) {


242  for (int j = 0; j < N; j++) {


243  sp.B[i, j] = sp.C[i, j];


244  }


245  }


246  var success = Eigendecomposition(N, sp.B, sp.D);


247 


248  DegenerateStateParameter.ActualValue = new BoolValue(!success);


249 


250  // assign D to eigenvalue square roots


251  for (int i = 0; i < N; i++) {


252  if (sp.D[i] < 0) { // numerical problem?


253  DegenerateStateParameter.ActualValue.Value = true;


254  sp.D[i] = 0;


255  } else sp.D[i] = Math.Sqrt(sp.D[i]);


256  }


257 


258  if (sp.D.Min() == 0.0) sp.AxisRatio = double.PositiveInfinity;


259  else sp.AxisRatio = sp.D.Max() / sp.D.Min();


260  }


261  return base.Apply();


262  }


263 


264  private bool Eigendecomposition(int N, double[,] B, double[] diagD) {


265  bool result = true;


266  // eigendecomposition


267  var offdiag = new double[N];


268  try {


269  tred2(N, B, diagD, offdiag);


270  tql2(N, diagD, offdiag, B);


271  } catch { result = false; }


272 


273  return result;


274  } // eigendecomposition


275 


276 


277  // Symmetric Householder reduction to tridiagonal form, taken from JAMA package.


278  private void tred2(int n, double[,] V, double[] d, double[] e) {


279 


280  // This is derived from the Algol procedures tred2 by


281  // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for


282  // Auto. Comp., Vol.iiLinear Algebra, and the corresponding


283  // Fortran subroutine in EISPACK.


284 


285  for (int j = 0; j < n; j++) {


286  d[j] = V[n  1, j];


287  }


288 


289  // Householder reduction to tridiagonal form.


290 


291  for (int i = n  1; i > 0; i) {


292 


293  // Scale to avoid under/overflow.


294 


295  double scale = 0.0;


296  double h = 0.0;


297  for (int k = 0; k < i; k++) {


298  scale = scale + Math.Abs(d[k]);


299  }


300  if (scale == 0.0) {


301  e[i] = d[i  1];


302  for (int j = 0; j < i; j++) {


303  d[j] = V[i  1, j];


304  V[i, j] = 0.0;


305  V[j, i] = 0.0;


306  }


307  } else {


308 


309  // Generate Householder vector.


310 


311  for (int k = 0; k < i; k++) {


312  d[k] /= scale;


313  h += d[k] * d[k];


314  }


315  double f = d[i  1];


316  double g = Math.Sqrt(h);


317  if (f > 0) {


318  g = g;


319  }


320  e[i] = scale * g;


321  h = h  f * g;


322  d[i  1] = f  g;


323  for (int j = 0; j < i; j++) {


324  e[j] = 0.0;


325  }


326 


327  // Apply similarity transformation to remaining columns.


328 


329  for (int j = 0; j < i; j++) {


330  f = d[j];


331  V[j, i] = f;


332  g = e[j] + V[j, j] * f;


333  for (int k = j + 1; k <= i  1; k++) {


334  g += V[k, j] * d[k];


335  e[k] += V[k, j] * f;


336  }


337  e[j] = g;


338  }


339  f = 0.0;


340  for (int j = 0; j < i; j++) {


341  e[j] /= h;


342  f += e[j] * d[j];


343  }


344  double hh = f / (h + h);


345  for (int j = 0; j < i; j++) {


346  e[j] = hh * d[j];


347  }


348  for (int j = 0; j < i; j++) {


349  f = d[j];


350  g = e[j];


351  for (int k = j; k <= i  1; k++) {


352  V[k, j] = (f * e[k] + g * d[k]);


353  }


354  d[j] = V[i  1, j];


355  V[i, j] = 0.0;


356  }


357  }


358  d[i] = h;


359  }


360 


361  // Accumulate transformations.


362 


363  for (int i = 0; i < n  1; i++) {


364  V[n  1, i] = V[i, i];


365  V[i, i] = 1.0;


366  double h = d[i + 1];


367  if (h != 0.0) {


368  for (int k = 0; k <= i; k++) {


369  d[k] = V[k, i + 1] / h;


370  }


371  for (int j = 0; j <= i; j++) {


372  double g = 0.0;


373  for (int k = 0; k <= i; k++) {


374  g += V[k, i + 1] * V[k, j];


375  }


376  for (int k = 0; k <= i; k++) {


377  V[k, j] = g * d[k];


378  }


379  }


380  }


381  for (int k = 0; k <= i; k++) {


382  V[k, i + 1] = 0.0;


383  }


384  }


385  for (int j = 0; j < n; j++) {


386  d[j] = V[n  1, j];


387  V[n  1, j] = 0.0;


388  }


389  V[n  1, n  1] = 1.0;


390  e[0] = 0.0;


391  }


392 


393  // Symmetric tridiagonal QL algorithm, taken from JAMA package.


394  private void tql2(int n, double[] d, double[] e, double[,] V) {


395 


396  // This is derived from the Algol procedures tql2, by


397  // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for


398  // Auto. Comp., Vol.iiLinear Algebra, and the corresponding


399  // Fortran subroutine in EISPACK.


400 


401  for (int i = 1; i < n; i++) {


402  e[i  1] = e[i];


403  }


404  e[n  1] = 0.0;


405 


406  double f = 0.0;


407  double tst1 = 0.0;


408  double eps = Math.Pow(2.0, 52.0);


409  for (int l = 0; l < n; l++) {


410 


411  // Find small subdiagonal element


412 


413  tst1 = Math.Max(tst1, Math.Abs(d[l]) + Math.Abs(e[l]));


414  int m = l;


415  while (m < n) {


416  if (Math.Abs(e[m]) <= eps * tst1) {


417  break;


418  }


419  m++;


420  }


421 


422  // If m == l, d[l] is an eigenvalue,


423  // otherwise, iterate.


424 


425  if (m > l) {


426  int iter = 0;


427  do {


428  iter = iter + 1; // (Could check iteration count here.)


429 


430  // Compute implicit shift


431 


432  double g = d[l];


433  double p = (d[l + 1]  g) / (2.0 * e[l]);


434  double r = hypot(p, 1.0);


435  if (p < 0) {


436  r = r;


437  }


438  d[l] = e[l] / (p + r);


439  d[l + 1] = e[l] * (p + r);


440  double dl1 = d[l + 1];


441  double h = g  d[l];


442  for (int i = l + 2; i < n; i++) {


443  d[i] = h;


444  }


445  f = f + h;


446 


447  // Implicit QL transformation.


448 


449  p = d[m];


450  double c = 1.0;


451  double c2 = c;


452  double c3 = c;


453  double el1 = e[l + 1];


454  double s = 0.0;


455  double s2 = 0.0;


456  for (int i = m  1; i >= l; i) {


457  c3 = c2;


458  c2 = c;


459  s2 = s;


460  g = c * e[i];


461  h = c * p;


462  r = hypot(p, e[i]);


463  e[i + 1] = s * r;


464  s = e[i] / r;


465  c = p / r;


466  p = c * d[i]  s * g;


467  d[i + 1] = h + s * (c * g + s * d[i]);


468 


469  // Accumulate transformation.


470 


471  for (int k = 0; k < n; k++) {


472  h = V[k, i + 1];


473  V[k, i + 1] = s * V[k, i] + c * h;


474  V[k, i] = c * V[k, i]  s * h;


475  }


476  }


477  p = s * s2 * c3 * el1 * e[l] / dl1;


478  e[l] = s * p;


479  d[l] = c * p;


480 


481  // Check for convergence.


482 


483  } while (Math.Abs(e[l]) > eps * tst1);


484  }


485  d[l] = d[l] + f;


486  e[l] = 0.0;


487  }


488 


489  // Sort eigenvalues and corresponding vectors.


490 


491  for (int i = 0; i < n  1; i++) {


492  int k = i;


493  double p = d[i];


494  for (int j = i + 1; j < n; j++) {


495  if (d[j] < p) { // NH find smallest k>i


496  k = j;


497  p = d[j];


498  }


499  }


500  if (k != i) {


501  d[k] = d[i]; // swap k and i


502  d[i] = p;


503  for (int j = 0; j < n; j++) {


504  p = V[j, i];


505  V[j, i] = V[j, k];


506  V[j, k] = p;


507  }


508  }


509  }


510  }


511 


512  /** sqrt(a^2 + b^2) without under/overflow. **/


513  private double hypot(double a, double b) {


514  double r = 0;


515  if (Math.Abs(a) > Math.Abs(b)) {


516  r = b / a;


517  r = Math.Abs(a) * Math.Sqrt(1 + r * r);


518  } else if (b != 0) {


519  r = a / b;


520  r = Math.Abs(b) * Math.Sqrt(1 + r * r);


521  }


522  return r;


523  }


524  }


525  } 
