source: branches/MOCMAEvolutionStrategy/HeuristicLab.Algorithms.MOCMAEvolutionStrategy/3.3/MocmaesIndividual.cs @ 14498

using System;
using System.Linq;
using HeuristicLab.Algorithms.MOCMAEvolutionStrategy;
using HeuristicLab.Common;
using HeuristicLab.Encodings.RealVectorEncoding;
using HeuristicLab.Random;

namespace HeuristicLab.Algorithms.CMAEvolutionStrategy {
  internal class MOCMAESIndividual {
    public const int Success = 1;
    public const int NoSuccess = 2;
    public const int Failure = 3;

    //MOCMAES-Params
    private readonly MOCMAESParameters strategy;

    //Chromosome
    public RealVector Mean { get; private set; }
    private double sigma;//stepsize
    private RealVector evolutionPath; // pc
    private RealVector lastStep;
    private RealVector lastZ;
    private double[,] lowerCholesky;


    //Phenotype
    public double[] Fitness { get; set; }
    public double[] PenalizedFitness { get; set; }
    public bool Selected { get; set; }
    public double Rank { get; set; }
    public double SuccessProbability { get; set; }

    /// <summary>
    /// 
    /// </summary>
    /// <param name="mean">has to be 0-vector with correct lenght</param>
    /// <param name="pSucc">has to be ptargetsucc</param>
    /// <param name="sigma">initialSigma</param>
    /// <param name="pc">has to be 0-vector with correct lenght</param>
    /// <param name="c">has to be a symmetric positive definit Covariance matrix</param>
    public MOCMAESIndividual(RealVector mean, double pSucc, double sigma, RealVector pc, double[,] c, MOCMAESParameters strategy) {
      Mean = mean;
      SuccessProbability = pSucc;
      this.sigma = sigma;
      evolutionPath = pc;
      CholeskyDecomposition(c);
      Selected = true;
      this.strategy = strategy;
    }

    private void CholeskyDecomposition(double[,] c) {
      if (!alglib.spdmatrixcholesky(ref c, c.GetLength(0), false))
        throw new ArgumentException("Covariancematrix is not symmetric positiv definit");
      lowerCholesky = (double[,])c.Clone();
    }

    public MOCMAESIndividual(MOCMAESIndividual other) {
      SuccessProbability = other.SuccessProbability;
      sigma = other.sigma;
      evolutionPath = (RealVector)other.evolutionPath.Clone();
      Mean = (RealVector)other.Mean.Clone();
      lowerCholesky = (double[,])other.lowerCholesky.Clone();
      Selected = true;
      strategy = other.strategy;
    }

    public void UpdateEvolutionPath(double learningRate, double updateWeight) {
      updateWeight = Math.Sqrt(updateWeight);
      for (var i = 0; i < evolutionPath.Length; i++) {
        evolutionPath[i] *= learningRate;
        evolutionPath[i] += updateWeight * lastStep[i];
      }
    }

    public double GetNormSqr() {
      return lastZ.Sum(d => d * d);
    }

    private void CholeskyUpdate(RealVector v, double alpha, double beta) {
      var n = v.Length;
      var temp = new double[n];
      for (var i = 0; i < n; i++) temp[i] = v[i];
      double betaPrime = 1;
      var a = Math.Sqrt(alpha);
      for (var j = 0; j < n; j++) {
        var ljj = a * lowerCholesky[j, j];
        var dj = ljj * ljj;
        var wj = temp[j];
        var swj2 = beta * wj * wj;
        var gamma = dj * betaPrime + swj2;
        var x1 = dj + swj2 / betaPrime;
        if (x1 < 0.0) throw new ArgumentException("Update makes Covariancematrix indefinite");
        var nLjj = Math.Sqrt(x1);
        lowerCholesky[j, j] = nLjj;
        betaPrime += swj2 / dj;
        if (j + 1 >= n) continue;
        for (var i = j + 1; i < n; i++) lowerCholesky[i, j] *= a;
        for (var i = j + 1; i < n; i++) temp[i] = wj / ljj * lowerCholesky[i, j];
        if (gamma.IsAlmost(0)) continue;
        for (var i = j + 1; i < n; i++) lowerCholesky[i, j] *= nLjj / ljj;
        for (var i = j + 1; i < n; i++) lowerCholesky[i, j] += nLjj * beta * wj / gamma * temp[i];

      }

    }

    public void Mutate(NormalDistributedRandom gauss) {
      //sampling a random z from N(0,I) where I is the Identity matrix;
      lastZ = new RealVector(Mean.Length);
      var n = lastZ.Length;
      for (var i = 0; i < n; i++) lastZ[i] = gauss.NextDouble();
      //Matrixmultiplication: lastStep = lowerCholesky * lastZ;
      lastStep = new RealVector(Mean.Length);
      for (var i = 0; i < n; i++) {
        double sum = 0;
        for (var j = 0; j <= i; j++) sum += lowerCholesky[i, j] * lastZ[j];
        lastStep[i] = sum;
      }
      //add the step to x weighted by stepsize;
      for (var i = 0; i < Mean.Length; i++) Mean[i] += sigma * lastStep[i];
    }

    public void UpdateAsParent(int v) {
      SuccessProbability = (1 - strategy.StepSizeLearningRate) * SuccessProbability + strategy.StepSizeLearningRate * (v == Success ? 1 : 0);
      sigma *= Math.Exp(1 / strategy.StepSizeDampeningFactor * (SuccessProbability - strategy.TargetSuccessProbability) / (1 - strategy.TargetSuccessProbability));
      if (v != Failure) return;
      if (SuccessProbability < strategy.SuccessThreshold) {
        var stepNormSqr = GetNormSqr();
        var rate = strategy.CovarianceMatrixUnlearningRate;
        if (stepNormSqr > 1 && 1 < strategy.CovarianceMatrixUnlearningRate * (2 * stepNormSqr - 1)) {
          rate = 1 / (2 * stepNormSqr - 1);
        }
        CholeskyUpdate(lastStep, 1 + rate, -rate);
      } else {
        RoundUpdate();
      }

    }

    public void UpdateAsOffspring() {
      SuccessProbability = (1 - strategy.StepSizeLearningRate) * SuccessProbability + strategy.StepSizeLearningRate;
      sigma *= Math.Exp(1 / strategy.StepSizeDampeningFactor * (SuccessProbability - strategy.TargetSuccessProbability) / (1 - strategy.TargetSuccessProbability));
      var evolutionpathUpdateWeight = strategy.EvolutionPathLearningRate * (2.0 - strategy.EvolutionPathLearningRate);
      if (SuccessProbability < strategy.SuccessThreshold) {
        UpdateEvolutionPath(1 - strategy.EvolutionPathLearningRate, evolutionpathUpdateWeight);
        CholeskyUpdate(evolutionPath, 1 - strategy.CovarianceMatrixLearningRate, strategy.CovarianceMatrixLearningRate);
      } else {
        RoundUpdate();
      }
    }

    private void RoundUpdate() {
      var evolutionPathUpdateWeight = strategy.EvolutionPathLearningRate * (2.0 - strategy.EvolutionPathLearningRate);
      UpdateEvolutionPath(1 - strategy.EvolutionPathLearningRate, 0);
      CholeskyUpdate(evolutionPath, 1 - strategy.CovarianceMatrixLearningRate + evolutionPathUpdateWeight, strategy.CovarianceMatrixLearningRate);
    }

  }

}