muwcommalambdacmafltvecop.cpp

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

    }
    mXmean->resize(lN);
    for(unsigned int i=0; i<lN; ++i) (*mXmean)[i] = (*lMeanFloatVec)[i];

    // Generate lambda children with breeder tree, first build breeder roulette
    RouletteT<unsigned int> lRoulette;
    buildRoulette(lRoulette, ioContext);

    // Keep best individuals if elitism is used
    const unsigned int lElitismKS=mElitismKeepSize->getWrappedValue();
    if(lElitismKS > 0) {
      Individual::Bag lBestInd;
      std::make_heap(ioDeme.begin(), ioDeme.end(), IsLessPointerPredicate());
      for(unsigned int i=0; i<lElitismKS; ++i) {
        lBestInd.push_back(ioDeme[0]);
        std::pop_heap(ioDeme.begin(), ioDeme.end(), IsLessPointerPredicate());
        ioDeme.pop_back();
      }
      ioDeme.resize(0);
      ioDeme.insert(ioDeme.end(), lBestInd.begin(), lBestInd.end());
    }
    else ioDeme.resize(0);

    // Generate the children
    const unsigned int lLambda =
      (unsigned int)std::ceil(mLMRatio->getWrappedValue()*double(lDemeSize));
    Individual::Bag lBagWithMeanInd;
    lBagWithMeanInd.push_back(lMeanInd);
    for(unsigned int i=0; i<lLambda; ++i) {
      unsigned int lIndexBreeder = lRoulette.select(ioContext.getSystem().getRandomizer());
      BreederNode::Handle lSelectedBreeder=getRootNode();
      for(unsigned int j=0; j<lIndexBreeder; ++j)
        lSelectedBreeder=lSelectedBreeder->getNextSibling();
      Beagle_NonNullPointerAssertM(lSelectedBreeder);
      Beagle_NonNullPointerAssertM(lSelectedBreeder->getBreederOp());
      Individual::Handle lBredIndiv =
        lSelectedBreeder->getBreederOp()->breed(lBagWithMeanInd,
                                                lSelectedBreeder->getFirstChild(),
                                                ioContext);
      Beagle_NonNullPointerAssertM(lBredIndiv);
      ioDeme.push_back(lBredIndiv);
    }

    // Check if all individuals have known fitness.
    for(unsigned int i=0; i<ioDeme.size(); ++i) {
      // If there is one invalid fitness, we should exit.
      // Evaluation will be taken elsewhere (we hope), and actual operator will be recalled.
      if((ioDeme[i]->getFitness()==NULL) || (ioDeme[i]->getFitness()->isValid()==false)) return;
    }

  }

  // Keep mu best children
  Beagle_AssertM(ioDeme.size() > lDemeSize);
  std::sort(ioDeme.begin(), ioDeme.end(), IsMorePointerPredicate());
  ioDeme.resize(lDemeSize);

  // Log messages
  Beagle_LogTraceM(
    ioContext.getSystem().getLogger(),
    "replacement-strategy", "Beagle::GA::MuWCommaLambdaCMAFltVecOp",
    "Updating covariance matrix, sigma, and cumulation paths of CMA-ES."
  );

  // Compute some constants
  const double lCC = 4.0 / (double(lN) + 4.0);
  const double lCS = (lMuEff+2.0) / (double(lN)+lMuEff+3.0);
  const double lMuCov = lMuEff;
  double lCCov = (((2.0*lMuEff)-1.0) / ((double(lN)+2.0)*(double(lN)+2.0)) + (2.0*lMuEff));
  lCCov = ((1.0/lMuCov) * (2.0/((double(lN)+1.414)*(double(lN)+1.414)))) +
    ((1.0-(1.0/lMuCov)) * minOf(1.0,lCCov));
  double lDamps = std::sqrt((lMuEff-1.0) / (double(lN)+1.0)) - 1.0;
  lDamps = 1.0 + (2.0*maxOf(0.0,lDamps)) + lCS;
  const double lChiN = std::sqrt(double(lN)) *
    (1.0 - (0.25/double(lN)) + (1.0/(21.0*double(lN)*double(lN))));

  // Compute new xmean
  Vector lXmean_new(lN, 0.0);
  for(unsigned int i=0; i<ioDeme.size(); ++i) {
    GA::FloatVector::Handle lVecI = castHandleT<GA::FloatVector>((*ioDeme[i])[0]);
    for(unsigned int j=0; j<lN; ++j) lXmean_new[j] += (lWeight[i] * (*lVecI)[j]);
  }

  // Compute zmean
  const double lSigma=mSigma->getWrappedValue();
  Vector lZmean(lN,0.0);
  for(unsigned int i=0; i<lN; ++i) lZmean[i] = (lXmean_new[i] - (*mXmean)[i]) / lSigma;
  Matrix lBt;
  mB->transpose(lBt);
  lZmean = lBt * lZmean;
  for(unsigned int i=0; i<lN; ++i) lZmean[i] /= (*mD)[i];

  // Update cumulation paths
  Matrix lBD(lN,lN,0.0);
  for(unsigned int i=0; i<lN; ++i) lBD(i,i) = (*mD)[i];
  lBD = (*mB) * lBD;
  Vector lBZm;
  mB->multiply((Matrix&)lBZm, lZmean);
  lBZm   *= std::sqrt(lCS * (2.0-lCS) * lMuEff);
  (*mPS) *= (1.0-lCS);
  (*mPS) += lBZm;

  double lPSnorm = 0.0;
  for(unsigned int i=0; i<mPS->size(); ++i) lPSnorm += ((*mPS)[i] * (*mPS)[i]);
  lPSnorm = std::sqrt(lPSnorm);
  const double lHLeft =
    lPSnorm / std::sqrt(1.0-std::pow(1.0-lCS, 2.0*double(ioContext.getGeneration()+1)));
  const double lHRight = (1.5+(1.0/(double(lN)-0.5))) * lChiN;
  const bool lHSig = (lHLeft<lHRight);

  (*mPC) *= (1.0-lCC);
  if(lHSig) {
    Vector lBDZm;
    lBD.multiply((Matrix&)lBDZm, lZmean);
    lBDZm *= std::sqrt(lCC * (2.0-lCC) * lMuEff);
    (*mPC) += lBDZm;
  }

  // Adapt step size sigma
  mSigma->getWrappedValue() *= std::exp((lCS/lDamps) * ((lPSnorm/lChiN)-1.0));

  // Adapt covariance matrix C
  Matrix lBDt;
  lBD.transpose(lBDt);
  Matrix lC;
  lBD.multiply(lC, lBDt);

  Matrix lR1Upd;           // Rank one update
  mPC->transpose(lR1Upd);
  lR1Upd = (*mPC) * lR1Upd;
  if(lHSig==false) {
    Matrix lCmod = lC;
    lCmod *= lCC * (2.0 - lCC);
    lR1Upd += lCmod;
  }
  lR1Upd *= (lCCov / lMuCov);

  Matrix lRMuUpd(lN,lN,0.0);      // Rank Mu update
  for(unsigned int i=0; i<ioDeme.size(); ++i) {
    GA::FloatVector::Handle lVecI = castHandleT<GA::FloatVector>((*ioDeme[i])[0]);
    Vector lX_I(lVecI->size());
    for(unsigned int j=0; j<lVecI->size(); ++j) lX_I[j] = (*lVecI)[j];
    lX_I -= (*mXmean);
    lX_I *= (1.0 / lSigma);
    Matrix lX_It;
    lX_I.transpose(lX_It);
    lX_It *= lWeight[i];
    lRMuUpd += (lX_I * lX_It);
  }
  lRMuUpd *= (lCCov * (1.0 - (1.0/lMuCov)));

  lC *= (1.0-lCCov);              // Attenuate old matrix
  lC += lR1Upd;                   // Add rank one update
  lC += lRMuUpd;                  // Add rank mu update

  // Update B and D from C
  for(unsigned int i=1; i<lN; ++i) {
    for(unsigned int j=0; j<i; ++j) lC(i,j) = lC(j,i);  // Enforce symetry
  }
  lC.computeEigens(*mD, *mB);                      // Principal component analysis

  // Handle too large differences in eigenvalues D
  const double lMaxD = mD->getMax();
  double lMinD = (*mD)[0];
  for(unsigned int i=1; i<mD->size(); ++i) {
    lMinD = ((*mD)[i] < lMinD) ? (*mD)[i] : lMinD;
  }
  if(lMaxD > (1e14*lMinD)) (*mD) += ((lMaxD/1e14)-lMinD);

  // Adjust D as standard deviation
  for(unsigned int i=0; i<lN; ++i) (*mD)[i] = std::sqrt((*mD)[i]);

  // Log updated parameters.
  Beagle_LogTraceM(
    ioContext.getSystem().getLogger(),
    "replacement-strategy", "Beagle::GA::MuWCommaLambdaCMAFltVecOp",
    string("CMA-ES updated B matrix (principal components): ")+mB->serialize()
  );
  Beagle_LogTraceM(
    ioContext.getSystem().getLogger(),
    "replacement-strategy", "Beagle::GA::MuWCommaLambdaCMAFltVecOp",
    string("CMA-ES updated D vector (standard deviations): ")+mD->serialize()
  );
  Beagle_LogTraceM(
    ioContext.getSystem().getLogger(),
    "replacement-strategy", "Beagle::GA::MuWCommaLambdaCMAFltVecOp",
    string("CMA-ES updated sigma: ")+mSigma->serialize()
  );
  Beagle_LogVerboseM(
    ioContext.getSystem().getLogger(),
    "replacement-strategy", "Beagle::GA::MuWCommaLambdaCMAFltVecOp",
    string("CMA-ES updated p_c vector (covariance cumulation path): ")+mPC->serialize()
  );
  Beagle_LogVerboseM(
    ioContext.getSystem().getLogger(),
    "replacement-strategy", "Beagle::GA::MuWCommaLambdaCMAFltVecOp",
    string("CMA-ES updated p_s vector (sigma cumulation path): ")+mPS->serialize()
  );
  Beagle_StackTraceEndM("void GA::MuWCommaLambdaCMAFltVecOp::operate(Deme& ioDeme, Context& ioContext)");
}


/*!
 *  \brief Post-initialization hook.
 *  \param ioSystem Reference to the evolutionary system.
 *  \throw Beagle::ValidationException If a parameter is invalid.
 */
void GA::MuWCommaLambdaCMAFltVecOp::postInit(System& ioSystem)
{
  Beagle_StackTraceBeginM();
  MuCommaLambdaOp::postInit(ioSystem);

  UInt::Handle lFloatVectorSize;
  if(ioSystem.getRegister().isRegistered("ga.init.vectorsize")) {
    lFloatVectorSize = castHandleT<UInt>(ioSystem.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();

  if((mB->getRows()==0) && (mB->getCols()==0)) mB->setIdentity(lN);
  else if((mB->getRows()!=mB->getCols()) || (mB->getRows()!=lN)) {
    std::ostringstream lOSS;
    lOSS << "Matrix B given by parameter \"ga.cmaes.b\" must be square ";
    lOSS << "and with a number of rows (and columns) equals to the size ";
    lOSS << "of the float vector individuals.";
    throw ValidationException(lOSS.str().c_str());
  }

  if(mD->size()==0) {
    mD->resize(lN);
    for(unsigned int i=0; i<lN; ++i) (*mD)[i] = 1.0;
  }
  else if(mD->size()!=lN) {
    std::ostringstream lOSS;
    lOSS << "Vector D given by parameter \"ga.cmaes.d\" must have a size ";
    lOSS << "equals to the size of the float vector individuals.";
    throw ValidationException(lOSS.str().c_str());
  }

  if(mPC->size()==0) {
    mPC->resize(lN);
    for(unsigned int i=0; i<lN; ++i) (*mPC)[i] = 0.0;
  }
  else if(mPC->getRows()!=lN) {
    std::ostringstream lOSS;
    lOSS << "Cumulation path vector p_c given by parameter \"ga.cmaes.pc\" ";
    lOSS << "must have the the size of the float vector individuals.";
    throw ValidationException(lOSS.str().c_str());
  }

  if(mPS->size()==0) {
    mPS->resize(lN);
    for(unsigned int i=0; i<lN; ++i) (*mPS)[i] = 0.0;
  }
  else if(mPS->getRows()!=lN) {
    std::ostringstream lOSS;
    lOSS << "Cumulation path vector p_s given by parameter \"ga.cmaes.ps\" ";
    lOSS << "must have the the size of the float vector individuals.";
    throw ValidationException(lOSS.str().c_str());
  }
  Beagle_StackTraceEndM("void GA::MuWCommaLambdaCMAFltVecOp::postInit(System& ioSystem)");
}