1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022013 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 {


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  /*double[] eVal;


247  double[,] eVec;


248  var success = alglib.smatrixevd(sp.B, N, 1, false, out eVal, out eVec);


249  sp.B = eVec;


250  sp.D = eVal;*/


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


252 


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


254 


255  // assign D to eigenvalue square roots


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


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


258  DegenerateStateParameter.ActualValue.Value = true;


259  sp.D[i] = 0;


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


261  }


262 


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


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


265  }


266  return base.Apply();


267  }


268 


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


270  bool result = true;


271  // eigendecomposition


272  var offdiag = new double[N];


273  try {


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


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


276  } catch { result = false; }


277 


278  return result;


279  } // eigendecomposition


280 


281 


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


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


284 


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


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


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


288  // Fortran subroutine in EISPACK.


289 


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


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


292  }


293 


294  // Householder reduction to tridiagonal form.


295 


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


297 


298  // Scale to avoid under/overflow.


299 


300  double scale = 0.0;


301  double h = 0.0;


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


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


304  }


305  if (scale == 0.0) {


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


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


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


309  V[i, j] = 0.0;


310  V[j, i] = 0.0;


311  }


312  } else {


313 


314  // Generate Householder vector.


315 


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


317  d[k] /= scale;


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


319  }


320  double f = d[i  1];


321  double g = Math.Sqrt(h);


322  if (f > 0) {


323  g = g;


324  }


325  e[i] = scale * g;


326  h = h  f * g;


327  d[i  1] = f  g;


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


329  e[j] = 0.0;


330  }


331 


332  // Apply similarity transformation to remaining columns.


333 


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


335  f = d[j];


336  V[j, i] = f;


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


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


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


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


341  }


342  e[j] = g;


343  }


344  f = 0.0;


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


346  e[j] /= h;


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


348  }


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


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


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


352  }


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


354  f = d[j];


355  g = e[j];


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


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


358  }


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


360  V[i, j] = 0.0;


361  }


362  }


363  d[i] = h;


364  }


365 


366  // Accumulate transformations.


367 


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


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


370  V[i, i] = 1.0;


371  double h = d[i + 1];


372  if (h != 0.0) {


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


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


375  }


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


377  double g = 0.0;


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


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


380  }


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


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


383  }


384  }


385  }


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


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


388  }


389  }


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


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


392  V[n  1, j] = 0.0;


393  }


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


395  e[0] = 0.0;


396  }


397 


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


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


400 


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


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


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


404  // Fortran subroutine in EISPACK.


405 


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


407  e[i  1] = e[i];


408  }


409  e[n  1] = 0.0;


410 


411  double f = 0.0;


412  double tst1 = 0.0;


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


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


415 


416  // Find small subdiagonal element


417 


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


419  int m = l;


420  while (m < n) {


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


422  break;


423  }


424  m++;


425  }


426 


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


428  // otherwise, iterate.


429 


430  if (m > l) {


431  int iter = 0;


432  do {


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


434 


435  // Compute implicit shift


436 


437  double g = d[l];


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


439  double r = hypot(p, 1.0);


440  if (p < 0) {


441  r = r;


442  }


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


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


445  double dl1 = d[l + 1];


446  double h = g  d[l];


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


448  d[i] = h;


449  }


450  f = f + h;


451 


452  // Implicit QL transformation.


453 


454  p = d[m];


455  double c = 1.0;


456  double c2 = c;


457  double c3 = c;


458  double el1 = e[l + 1];


459  double s = 0.0;


460  double s2 = 0.0;


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


462  c3 = c2;


463  c2 = c;


464  s2 = s;


465  g = c * e[i];


466  h = c * p;


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


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


469  s = e[i] / r;


470  c = p / r;


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


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


473 


474  // Accumulate transformation.


475 


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


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


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


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


480  }


481  }


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


483  e[l] = s * p;


484  d[l] = c * p;


485 


486  // Check for convergence.


487 


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


489  }


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


491  e[l] = 0.0;


492  }


493 


494  // Sort eigenvalues and corresponding vectors.


495 


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


497  int k = i;


498  double p = d[i];


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


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


501  k = j;


502  p = d[j];


503  }


504  }


505  if (k != i) {


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


507  d[i] = p;


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


509  p = V[j, i];


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


511  V[j, k] = p;


512  }


513  }


514  }


515  }


516 


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


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


519  double r = 0;


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


521  r = b / a;


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


523  } else if (b != 0) {


524  r = a / b;


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


526  }


527  return r;


528  }


529  }


530  } 
