muwcommalambdacmafltvecop.cpp

非常好的进化算法EC 实现平台 可以实现多种算法 GA GP

/*
 *  Open BEAGLE
 *  Copyright (C) 2001-2004 by Christian Gagne and Marc Parizeau
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This library 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
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 *  Contact:
 *  Laboratoire de Vision et Systemes Numeriques
 *  Departement de genie electrique et de genie informatique
 *  Universite Laval, Quebec, Canada, G1K 7P4
 *  http://vision.gel.ulaval.ca
 *
 */

/*!
 *  \file   beagle/GA/src/MuWCommaLambdaCMAFltVecOp.cpp
 *  \brief  Source code of class MuWCommaLambdaCMAFltVecOp.
 *  \author Christian Gagne
 *  \author Marc Parizeau
 *  $Revision: 1.11 $
 *  $Date: 2005/10/04 16:25:09 $
 */

#include "beagle/GA.hpp"

#include <cfloat>
#include <float.h>
#include <cmath>
#include <algorithm>

using namespace Beagle;


/*!
 *  \brief Build CMA-ES (Mu_W+Lambda) replacement strategy operator.
 *  \param inLMRatioName Lamda over Mu parameter name used in the register.
 *  \param inName Name of the CMA-ES (Mu_W+Lambda) operator.
 */
GA::MuWCommaLambdaCMAFltVecOp::MuWCommaLambdaCMAFltVecOp(Beagle::string inLMRatioName, 
														 Beagle::string inName) :
  MuCommaLambdaOp(inLMRatioName,inName)
{ }


/*!
 *  \brief Initialize the operator.
 *  \param ioSystem Reference to the evolutionary system.
 */
void GA::MuWCommaLambdaCMAFltVecOp::initialize(System& ioSystem)
{
  Beagle_StackTraceBeginM();
  MuCommaLambdaOp::initialize(ioSystem);

  if(ioSystem.getRegister().isRegistered("ga.cmaes.b")) {
    mB = castHandleT<Matrix>(ioSystem.getRegister()["ga.cmaes.b"]);
  } else {
    mB = new Matrix;
    Register::Description lDescription(
      "CMA-ES B matrix",
      "Matrix",
      "",
      "CMA-ES B matrix containing the covariance matrix eigenvectors."
    );
    ioSystem.getRegister().addEntry("ga.cmaes.b", mB, lDescription);
  }

  if(ioSystem.getRegister().isRegistered("ga.cmaes.d")) {
    mD = castHandleT<Vector>(ioSystem.getRegister()["ga.cmaes.d"]);
  } else {
    mD = new Vector;
    Register::Description lDescription(
      "CMA-ES D vector",
      "Vector",
      "",
      "CMA-ES D vector containing the square root of the covariance matrix eigenvalues."
    );
    ioSystem.getRegister().addEntry("ga.cmaes.d", mD, lDescription);
  }

  if(ioSystem.getRegister().isRegistered("ga.cmaes.pc")) {
    mPC = castHandleT<Vector>(ioSystem.getRegister()["ga.cmaes.pc"]);
  } else {
    mPC = new Vector;
    Register::Description lDescription(
      "CMA-ES p_c vector",
      "Vector",
      "",
      "CMA-ES p_c vector containing the covariance matrix evolution path."
    );
    ioSystem.getRegister().addEntry("ga.cmaes.pc", mPC, lDescription);
  }

  if(ioSystem.getRegister().isRegistered("ga.cmaes.ps")) {
    mPS = castHandleT<Vector>(ioSystem.getRegister()["ga.cmaes.ps"]);
  } else {
    mPS = new Vector;
    Register::Description lDescription(
      "CMA-ES p_s vector",
      "Vector",
      "",
      "CMA-ES p_s vector containing the sigma value evolution path."
    );
    ioSystem.getRegister().addEntry("ga.cmaes.ps", mPS, lDescription);
  }

  if(ioSystem.getRegister().isRegistered("ga.cmaes.sigma")) {
    mSigma = castHandleT<Double>(ioSystem.getRegister()["ga.cmaes.sigma"]);
  } else {
    mSigma = new Double(0.5);
    Register::Description lDescription(
      "CMA-ES sigma value",
      "Double",
      "0.5",
      "CMA-ES sigma value moduling the mutation step size."
    );
    ioSystem.getRegister().addEntry("ga.cmaes.sigma", mSigma, lDescription);
  }

  if(ioSystem.getRegister().isRegistered("ga.cmaes.xmean")) {
    mXmean = castHandleT<Vector>(ioSystem.getRegister()["ga.cmaes.xmean"]);
  } else {
    mXmean = new Vector;
    Register::Description lDescription(
      "CMA-ES mean vector",
      "Vector",
      "",
      "CMA-ES mean vector containing the average individual."
    );
    ioSystem.getRegister().addEntry("ga.cmaes.xmean", mXmean, lDescription);
  }

  if(ioSystem.getRegister().isRegistered("ga.float.maxvalue")) {
    mMaxValue = castHandleT<DoubleArray>(ioSystem.getRegister()["ga.float.maxvalue"]);
  } else {
    mMaxValue = new DoubleArray(1,DBL_MAX);
    std::ostringstream lOSS;
    lOSS << "Maximum values assigned to vector's floats. ";
    lOSS << "Value can be a scalar, which limit the value for all float ";
    lOSS << "vector parameters, or a vector which limit the value for the parameters ";
    lOSS << "individually. If the maximum value is smaller than the ";
    lOSS << "float vector size, the limit used for the last values of the float vector ";
    lOSS << "is equal to the last value of the maximum value vector.";
    Register::Description lDescription(
      "Maximum vector values",
      "DoubleArray",
      dbl2str(DBL_MAX),
      lOSS.str().c_str()
    );
    ioSystem.getRegister().addEntry("ga.float.maxvalue", mMaxValue, lDescription);
  }

  if(ioSystem.getRegister().isRegistered("ga.float.minvalue")) {
    mMinValue = castHandleT<DoubleArray>(ioSystem.getRegister()["ga.float.minvalue"]);
  } else {
    mMinValue = new DoubleArray(1,DBL_MIN);
    std::ostringstream lOSS;
    lOSS << "Minimum  values assigned to vector's floats. ";
    lOSS << "Value can be a scalar, which limit the value for all float ";
    lOSS << "vector parameters, or a vector which limit the value for the parameters ";
    lOSS << "individually. If the minimum value is smaller than the ";
    lOSS << "float vector size, the limit used for the last values of the float vector ";
    lOSS << "is equal to the last value of the minimum value vector.";
    Register::Description lDescription(
      "Minimum values",
      "DoubleArray",
      dbl2str(DBL_MIN),
      lOSS.str().c_str()
    );
    ioSystem.getRegister().addEntry("ga.float.minvalue", mMinValue, lDescription);
  }
  Beagle_StackTraceEndM("void GA::MuWCommaLambdaCMAFltVecOp::initialize(System& ioSystem)");
}


/*!
 *  \brief Apply the CMA-ES (Mu_W+Lambda) replacement strategy operation on a deme.
 *  \param ioDeme Reference to the deme on which the operation takes place.
 *  \param ioContext Evolutionary context of the operation.
 *  \throw Beagle::ValidationException If a parameter is missing or have a bad value.
 *  \throw Beagle::AssertException If an invalid condition appears.
 *
 *  This routine is mostly inspired from matlab routine purecmaes.m and the document
 *  "The CMA Evolution Strategy: A Tutorial" (October 15, 2004), both
 *  from Nikolaus Hansen.
 *  See http://www.bionik.tu-berlin.de/user/niko/cmaes_inmatlab.html
 */
void GA::MuWCommaLambdaCMAFltVecOp::operate(Deme& ioDeme, Context& ioContext)
{
  Beagle_StackTraceBeginM();
  // Get real popsize and size of float vectors from register.
  UIntArray::Handle lPopSize;
  if(ioContext.getSystem().getRegister().isRegistered("ec.pop.size")) {
    lPopSize = castHandleT<UIntArray>(ioContext.getSystem().getRegister()["ec.pop.size"]);
  }
  else {
    std::ostringstream lOSS;
    lOSS << "Population size parameter \"ec.pop.size\" is not found in register!";
    throw ValidationException(lOSS.str().c_str());
  }
  const unsigned int lDemeSize = (*lPopSize)[ioContext.getDemeIndex()];
  UInt::Handle lFloatVectorSize;
  if(ioContext.getSystem().getRegister().isRegistered("ga.init.vectorsize")) {
    lFloatVectorSize = castHandleT<UInt>(ioContext.getSystem().getRegister()["ga.init.vectorsize"]);
  }
  else {
    std::ostringstream lOSS;
    lOSS << "GA::MuWCommaLambdaCMAFltVecOp must be used in fixed-lenght float vector ";
    lOSS << "individuals. Parameter \"ga.init.vectorsize\" is not in register, ";
    lOSS << "while it is needed to set initial size of the different CMA-ES matrices ";
    lOSS << "and vectors.";
    throw ValidationException(lOSS.str().c_str());
  }
  const unsigned int lN=lFloatVectorSize->getWrappedValue();

  // Compute weights and effective mu  
  Vector lWeight(lDemeSize, std::log(double(lDemeSize+1)));
  double lSumWeight=0.0;
  double lMuEff=0.0;
  for(unsigned int i=0; i<lWeight.size(); ++i) {
    lWeight[i] -= std::log(double(i+1));
    lSumWeight += lWeight[i];
  }
  for(unsigned int i=0; i<lWeight.size(); ++i) {
    lWeight[i] /= lSumWeight;
    lMuEff += (lWeight[i] * lWeight[i]);
  }
  lMuEff = 1.0 / lMuEff;

  // If the replacement strategy possess a breeder tree
  if(getRootNode()!=NULL) {

    // Check parameters and log some information
    Beagle_AssertM(ioDeme.size()==lDemeSize);
    Beagle_NonNullPointerAssertM(mElitismKeepSize);
    Beagle_ValidateParameterM(mLMRatio->getWrappedValue() >= 1.0,
                              mLMRatioName,
                              "The LM ratio must be higher or equal to 1.0.");
    Beagle_ValidateParameterM(mElitismKeepSize->getWrappedValue() <= ioDeme.size(),
                              "ec.elite.keepsize",
                              "The elistism keepsize must be less than the deme size!");

    Beagle_LogTraceM(
      ioContext.getSystem().getLogger(),
      "replacement-strategy", "Beagle::GA::MuWCommaLambdaCMAFltVecOp",
      string("Using CMA-ES (mu_w,lambda) replacement strategy to process the ")+
      uint2ordinal(ioContext.getDemeIndex()+1)+" deme"
    );
    Beagle_LogObjectM(
      ioContext.getSystem().getLogger(),
      Logger::eTrace,
      "replacement-strategy", "Beagle::GA::MuWCommaLambdaCMAFltVecOp",
      (*this)
    );

    // Create weighted mean individual.
    std::sort(ioDeme.begin(), ioDeme.end(), IsMorePointerPredicate());
    Individual::Handle lMeanInd = castHandleT<Individual>(ioDeme.getTypeAlloc()->allocate());
    lMeanInd->resize(1);
    GA::FloatVector::Handle lMeanFloatVec = castHandleT<GA::FloatVector>((*lMeanInd)[0]);
    lMeanFloatVec->resize(lN, 0.0);
    for(unsigned int i=0; i<ioDeme.size(); ++i) {
      Beagle_AssertM(ioDeme[i]->size()==1);
      GA::FloatVector::Handle lVecI = castHandleT<GA::FloatVector>((*ioDeme[i])[0]);
      Beagle_AssertM(lVecI->size()==lN);
      for(unsigned int j=0; j<lN; ++j) (*lMeanFloatVec)[j] += (lWeight[i] * (*lVecI)[j]);