freESCMAMutation.txx

Go to the documentation of this file.

/*=========================================================================

  Program:   F.R.E.E. - flexible registration evaluation engine
  Version:   v.1.0.0
  Date:      $Date: 2006/09/01 12:00:00 $
  Module:    $RCSfile: freESCMAMutation.txx,v $
  Language:  C++



  Copyright (c) 2007 Ralf o Floca (Department of Medical Informatics,
  Institute for Medical Biometry and Informatics, University of Heidelberg,
  Germany). All rights reserved.
  See FREECopyright.txt or http://www.mi.med.uni-hd.de/free/copyright.htm
  for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef _freESCMAMutation_txx
#define _freESCMAMutation_txx

#include "freESCMAMutation.h"

#include "vnl/vnl_sample.h"
#include "vnl/algo/vnl_symmetric_eigensystem.h"

namespace FREE
{
namespace ES
{

template <class TIndividual>
CMAMutation<TIndividual>
::CMAMutation()
{
  m_Sigma = 1.0;
  m_Mu_eff = 1;
  m_C_c = 0.8;
  m_C_cov = 0.0;
  m_C_sigma = 0.6;

  m_Mu = 1;

  m_ObjectiveCount = 1;
  m_CurMatrix = MatrixType(1,1,0.0);
  m_CurCentroid = VectorType(1,0.0);
  m_CurMatrixPath = VectorType(1,0.0);
  m_CurStepPath = VectorType(1,0.0);
  m_CurStepLength = 1.0;

  m_SearchPointSelector = 0;
  m_CentroidCombinator = 0;
  m_RandomGenerator = VariateGeneratorType::New();
}

template <class TIndividual>
void
CMAMutation<TIndividual>
::InitializeControlValues_Mean(const long& lN)
{
  m_Mu_cov = 1/m_Mu_eff;
  
  m_C_c = 1.0/lN;
  m_C_cov = (std::min<double>(m_Mu_eff,lN*lN)) / (lN*lN);
  m_C_sigma = 1.0/lN;
  m_D_sigma = 1.0;
};

template <class TIndividual>
void
CMAMutation<TIndividual>
::InitializeControlValues_Hansen(const long& lMu, const long& lN)
{
  if (m_CentroidCombinator.IsNull()) itkExceptionMacro(<<"Error. Cannot initialize control values. No centroid combinator defined.");
  typename CentroidCombinatorType::WeightVectorType weights = m_CentroidCombinator->GetWeights(lMu);

  //calculate Mu_eff
  m_Mu_eff = 0;
  for (typename CentroidCombinatorType::WeightVectorType::iterator pos = weights.begin(); pos != weights.end(); ++pos)
  {
    m_Mu_eff += (*pos) * (*pos);
  }
  m_Mu_eff = 1.0 / m_Mu_eff;
  
  //calculate Mu_cov;
  m_Mu_cov = 1.0 / m_Mu_eff;

  //calculate C_c;
  m_C_c = 4.0 / ( 4.0+lN);

  //calculate C_cov;
  m_C_cov = m_Mu_cov * (2.0/((lN+sqrt(2.0))*(lN+sqrt(2.0))));
  m_C_cov += (1.0-m_Mu_cov)*std::min(1.0,(2*m_Mu_eff - 1.0)/(((lN+2)*(lN+2))+m_Mu_eff));

  //calculate C_sigma;
  m_C_sigma = (m_Mu_eff+2.0) / (lN+m_Mu_eff+3.0);

  //calculate D_sigma;
  m_D_sigma = 1+ 2*std::max<double>(0,sqrt((m_Mu_eff-1.0) / (lN+1.0)) - 1.0) + m_C_sigma;
};


template <class TIndividual>
void
CMAMutation<TIndividual>
::SetSearchPointSelector( SearchPointSelectorType * pSearchPointSelector )
{
  if (m_SearchPointSelector.GetPointer()!=pSearchPointSelector)
  {
    this->Modified();
    m_SearchPointSelector = pSearchPointSelector;
  }
};

template <class TIndividual>
void
CMAMutation<TIndividual>
::SetCentroidCombinator( CentroidCombinatorType * pCentroidCombinator )
{
  if (m_CentroidCombinator.GetPointer()!=pCentroidCombinator)
  {
    this->Modified();
    m_CentroidCombinator = pCentroidCombinator;
  }
};

template <class TIndividual>
void
CMAMutation<TIndividual>
::RegisterStrategicPopulationParameters(PopulationType& population) const
{
  //get any individual of the population to determine the number of objective values
  if (population.Size()<1) itkExceptionMacro(<<"Error. Need at least one individual in the population to determine number of objectiv values in order to properly register strategic population parameters.");
  m_ObjectiveCount = population.ElementAt(0)->ObjectiveParameters().size();

  //no initialze the members properly
  m_CurMatrix = MatrixType(m_ObjectiveCount,m_ObjectiveCount,0.0);
  m_CurMatrix.fill_diagonal(1.0);

  m_CurCentroid = ComputeCentroid(population); //VectorType(m_ObjectiveCount,0.0);

  m_CurMatrixPath = VectorType(m_ObjectiveCount,0.0);
  m_CurStepPath = VectorType(m_ObjectiveCount,0.0);
  
  m_CurStepLength = m_Sigma;

  ConvertMemberToStrategicPopulationParameters(population.StrategicParameters());
};

template <class TIndividual>
void
CMAMutation<TIndividual>
::MutateStrategicPopulationParameters(PopulationType& population) const
{
  ConvertStrategicPopulationParametersToMember(population.StrategicParameters());

  //compute the next strategic paremter (that will be used for the upcoming individuals
  VectorType nextCentroid = ComputeCentroid(population);

  VectorType nextMatrixPath = ComputeMatrixPath(m_CurMatrixPath, m_CurStepLength, m_CurCentroid, nextCentroid);

  MatrixType nextCovMatrix = ComputeMatrix(m_CurMatrix, nextMatrixPath, m_CurStepLength, m_CurCentroid, population);

  VectorType nextStepPath = ComputeStepPath(m_CurStepPath, m_CurMatrix, m_CurStepLength, m_CurCentroid, nextCentroid);

  double nextStepLength = ComputeStepLength(m_CurStepLength, nextStepPath);
  
  //store the newly computed parameters as new current paremeters
  m_CurCentroid = nextCentroid;
  m_CurMatrixPath = nextMatrixPath;
  m_CurMatrix = nextCovMatrix;
  m_CurStepPath = nextStepPath;
  m_CurStepLength = nextStepLength;

  ConvertMemberToStrategicPopulationParameters(population.StrategicParameters());
};

template <class TIndividual>
double
CMAMutation<TIndividual>
::MutateValue(const double& value)
{
  itkExceptionMacro(<<"Functionality not yet implemented.");
  return 0;
};

template <class TIndividual>
void
CMAMutation<TIndividual>
::MutateObjectivParameters(IndividualType* pIndividual, PopulationType* pPopulation) const
{
  if (!pIndividual) itkExceptionMacro(<<"Error. Passed individual is Null.");
  if (!pPopulation) itkExceptionMacro(<<"Error. Passed population is Null.");
  if (pIndividual->ObjectiveParameters().size() != m_ObjectiveCount) itkExceptionMacro(<<"Error. Passed individual has not the expected number of objective parameters. Found: "<<pIndividual->ObjectiveParameters().size()<<"; expected: "<<this->m_ObjectiveCount);
  if (this->m_ObjectiveCount!=this->m_GeneralObjectiveScales.size()) itkExceptionMacro(<<"Error. Number of objective parameters and global scales are not equal. Ensure that scales have been set properly. Objective count:" << this->m_ObjectiveCount <<"; objective scales count: "<< this->m_GeneralObjectiveScales.size());

  //get the strategic parameters
  ConvertStrategicPopulationParametersToMember(pPopulation->StrategicParameters());

  //decompose covariance matrix
  if (!m_CurMatrix.is_finite()) itkExceptionMacro(<<"Error. Current covariance matrix is not finite. Matrix: "<<m_CurMatrix);
  vnl_symmetric_eigensystem< double > eigensystem(m_CurMatrix);

  MatrixType b = eigensystem.V;
  MatrixType d = (eigensystem.D.asMatrix())*(-1);

  //get random vector;
  VectorType z = GetRandomVector(m_ObjectiveCount);

  //compute mutation vector
  VectorType mutationVector = z.pre_multiply(d);
  mutationVector = mutationVector.pre_multiply(b);
  mutationVector = mutationVector*m_CurStepLength;

  //mutate individual
  for (unsigned long index = 0; index<pIndividual->ObjectiveParameters().size(); index++)
  {
    typename IndividualType::ObjectiveParameterType mutatedParameter = (pIndividual->ObjectiveParameters())[index];

    //mutate parameter
    mutatedParameter += mutationVector[index] / this->m_GeneralObjectiveScales[index];
    
    //set new value;
    (pIndividual->ObjectiveParameters())[index] = mutatedParameter;
  }
};

template <class TIndividual>
void
CMAMutation<TIndividual>
::ConvertStrategicPopulationParametersToMember(PopulationStrategicParametersType& parameters) const
{
  typename PopulationType::StrategicParametersType::SelectedParametersType ownParameters =
    parameters.GetStrategicParameters(this->GetNameOfClass());

  //number of parameters stored; 1 marix, 3 vectors and 1 scalar
  long lNeededParameterCount = (m_ObjectiveCount*m_ObjectiveCount)+3*m_ObjectiveCount+1;

  if (ownParameters.size()!=lNeededParameterCount) throwExceptionMacro("Error. Invalid number of strategic parameters (one matrix (n*n), 3 vectors (n) and 1 scalar needed; n: number of objective values). Ensure that strategic parameters are properly registered. Expected parameter count:" << lNeededParameterCount <<"; current parameter count: "<< ownParameters.size());

  long lCurParamIndex = 0;

  //load covariance matrix
  for (MatrixType::iterator pos = m_CurMatrix.begin(); pos != m_CurMatrix.end(); ++pos)
  {
    *pos = (ownParameters[lCurParamIndex])->GetValue();
    lCurParamIndex++;
  }

  //load current centroid
  for (VectorType::iterator pos = m_CurCentroid.begin(); pos != m_CurCentroid.end(); ++pos)
  {
    *pos = (ownParameters[lCurParamIndex])->GetValue();
    lCurParamIndex++;
  }

  //load current evolution path vector for covariance matrix
  for (VectorType::iterator pos = m_CurMatrixPath.begin(); pos != m_CurMatrixPath.end(); ++pos)
  {
    *pos = (ownParameters[lCurParamIndex])->GetValue();
    lCurParamIndex++;
  }

  //load current evolution path vector for step length
  for (VectorType::iterator pos = m_CurStepPath.begin(); pos != m_CurStepPath.end(); ++pos)
  {
    *pos = (ownParameters[lCurParamIndex])->GetValue();
    lCurParamIndex++;
  }

  //load current evolution step length
  m_CurStepLength = (ownParameters[lCurParamIndex])->GetValue();
};

template <class TIndividual>
void
CMAMutation<TIndividual>
::ConvertMemberToStrategicPopulationParameters(PopulationStrategicParametersType& parameters) const
{
  parameters.DeleteSelectedStrategicParameters(parameters.GetStrategicIDs(this->GetNameOfClass()));

  typedef typename IndividualType::StrategicParameterPointer IndividualStrategicParameterPointer;
  //Store covariance matrix
  for (MatrixType::const_iterator pos = m_CurMatrix.begin(); pos != m_CurMatrix.end(); ++pos)
  {
    IndividualStrategicParameterPointer parameter = IndividualType::StrategicParameterType::New();
    parameter->SetValue(*pos);
    parameter->SetHandling(IndividualType::StrategicParameterType::HDynamic);
    parameter->SetOrigin(this->GetNameOfClass());
    parameters.push_back(parameter);
  }

  //Store current centroid
  for (VectorType::const_iterator pos = m_CurCentroid.begin(); pos != m_CurCentroid.end(); ++pos)
  {
    IndividualStrategicParameterPointer parameter = IndividualType::StrategicParameterType::New();
    parameter->SetValue(*pos);
    parameter->SetHandling(IndividualType::StrategicParameterType::HDynamic);
    parameter->SetOrigin(this->GetNameOfClass());
    parameters.push_back(parameter);
  }

  //Store current evolution path vector for covariance matrix
  for (VectorType::const_iterator pos = m_CurMatrixPath.begin(); pos != m_CurMatrixPath.end(); ++pos)
  {
    IndividualStrategicParameterPointer parameter = IndividualType::StrategicParameterType::New();
    parameter->SetValue(*pos);
    parameter->SetHandling(IndividualType::StrategicParameterType::HDynamic);
    parameter->SetOrigin(this->GetNameOfClass());
    parameters.push_back(parameter);
  }

  //Store current evolution path vector for step length control
  for (VectorType::const_iterator pos = m_CurStepPath.begin(); pos != m_CurStepPath.end(); ++pos)
  {
    IndividualStrategicParameterPointer parameter = IndividualType::StrategicParameterType::New();
    parameter->SetValue(*pos);
    parameter->SetHandling(IndividualType::StrategicParameterType::HDynamic);
    parameter->SetOrigin(this->GetNameOfClass());
    parameters.push_back(parameter);
  }

  //Store current step length
  IndividualStrategicParameterPointer parameter = IndividualType::StrategicParameterType::New();
  parameter->SetValue(m_CurStepLength);
  parameter->SetHandling(IndividualType::StrategicParameterType::HDynamic);
  parameter->SetOrigin(this->GetNameOfClass());
  parameters.push_back(parameter);
};

template <class TIndividual>
typename CMAMutation<TIndividual>::VectorType
CMAMutation<TIndividual>::
ComputeCentroid(PopulationType& rPopulation) const
{
  if (m_SearchPointSelector.IsNull()) itkExceptionMacro(<<"Error. Cannot compute centroid. No point selector defined.");
  if (m_CentroidCombinator.IsNull()) itkExceptionMacro(<<"Error. Cannot compute centroid. No centroid combinator defined.");

  typedef typename SearchPointSelectorType::ParentSelectionType PointSelection;
  PointSelection selection = m_SearchPointSelector->Select(rPopulation);

  m_Mu = selection.size();

  IndividualPointer centroidIndividual = m_CentroidCombinator->Recombine(selection);

  VectorType centroid = ConvertObjectivParameterToVector(centroidIndividual.GetPointer());
      
  return centroid;
};

template <class TIndividual>
typename CMAMutation<TIndividual>::VectorType
CMAMutation<TIndividual>::
ComputeMatrixPath(const VectorType& curMatrixPath, double curStepLength, const VectorType& curCentroid, const VectorType& nextCentroid) const
{
  VectorType nextPath = (1 - m_C_c)*curMatrixPath;
  double temp = sqrt(m_C_c*(2-m_C_c)*m_Mu_eff);
  
  VectorType diff = nextCentroid - curCentroid;

  //remove scaling of parameters for internal computation.
  //The rescaling will be done when mutating new individual.
  //pathes, steplength and covariance matrix are without scale
  //TODO: Better option would be to change the estimation of
  //n-dimensional variate in the step length computation to a
  //scaled version.
  for (unsigned long index = 0; index<diff.size(); index++)
  {
    diff[index] = diff[index] * this->m_GeneralObjectiveScales[index];
  }

  nextPath = nextPath + temp * (diff / curStepLength);
  
  return nextPath;
};

template <class TIndividual>
typename CMAMutation<TIndividual>::MatrixType
CMAMutation<TIndividual>::
ComputeMatrix(const MatrixType& curMatrix, const VectorType& nextMatrixPath, double curStepLength, const VectorType& curCentroid, const PopulationType& population) const
{
  //rank one update
  MatrixType rankOneUpdate = outer_product(nextMatrixPath,nextMatrixPath);
  
  //compute rank my update
    //find all new individuals of the last generation
  MatrixType rankMuUpdate(m_ObjectiveCount,m_ObjectiveCount,0.0);

  typedef std::vector<const TIndividual*> SelectedIndividualsType;
  SelectedIndividualsType newIndividuals;
  for (typename PopulationType::const_iterator pos = population.begin(); pos!=population.end(); ++pos)
  {
    if ((*pos)->GetGenerationID()==population.GetGenerationID())
    {
      newIndividuals.push_back((*pos).GetPointer());
    }
  }
    //calculate matrix

  if (newIndividuals.size())
  {
    for (typename SelectedIndividualsType::iterator pos = newIndividuals.begin(); pos!=newIndividuals.end(); ++pos)
    {
      VectorType individual = ConvertObjectivParameterToVector(*pos);
      VectorType variance = curStepLength*(individual - curCentroid);

      //remove scaling of parameters for internal computation.
      //The rescaling will be done when mutating new individual.
      //pathes, steplength and covariance matrix are without scale
      //TODO: Better option would be to change the estimation of
      //n-dimensional variate in the step length computation to a
      //scaled version.
      for (unsigned long index = 0; index<variance.size(); index++)
      {
        variance[index] = variance[index] * this->m_GeneralObjectiveScales[index];
      }

      rankMuUpdate = rankMuUpdate + outer_product(variance, variance);
    }
    rankMuUpdate = rankMuUpdate*(1/newIndividuals.size());
  }

  MatrixType nextMatrix = (1-m_C_cov)*curMatrix + ((m_C_cov*m_Mu_cov)*rankOneUpdate) + (m_C_cov*(1-m_Mu_cov)*rankMuUpdate);

  return nextMatrix;
};

template <class TIndividual>
typename CMAMutation<TIndividual>::VectorType
CMAMutation<TIndividual>::
ComputeStepPath(const VectorType& curStepPath,const MatrixType& curMatrix, double curStepLength, const VectorType& curCentroid, const VectorType& nextCentroid) const
{
  VectorType nextPath = (1 - m_C_sigma)*curStepPath;
  double temp = sqrt(m_C_sigma*(2-m_C_sigma)*m_Mu_eff);

  //compute C^-1/2 to scale the path independently from the direction
  vnl_symmetric_eigensystem< double > eigensystem(curMatrix);

  MatrixType B = eigensystem.V;
  MatrixType D_inv = (eigensystem.D.asMatrix())*(-1);
  MatrixType B_inv = B.transpose();

  VectorType diff = nextCentroid - curCentroid;

  //remove scaling of parameters for internal computation.
  //The rescaling will be done when mutating new individual.
  //pathes, steplength and covariance matrix are without scale
  //TODO: Better option would be to change the estimation of
  //n-dimensional variate in the step length computation to a
  //scaled version.
  for (unsigned long index = 0; index<diff.size(); index++)
  {
    double debugCheck = diff[index] * this->m_GeneralObjectiveScales[index];
    diff[index] = debugCheck;
  }

  nextPath = nextPath + temp * (B * D_inv * B_inv) * (diff / curStepLength);
  
  return nextPath;
};

template <class TIndividual>
double
CMAMutation<TIndividual>::
ComputeStepLength(double curStepLength, const VectorType& nextStepPath) const
{
  double estimatedN = sqrt((double)m_ObjectiveCount)*(1- (1/(4*m_ObjectiveCount)) + (1/(21*m_ObjectiveCount*m_ObjectiveCount)));

    double debugCheck = nextStepPath.magnitude();

  double nextStepLength = curStepLength * exp((m_C_sigma/m_D_sigma)*((nextStepPath.magnitude()/estimatedN) -1));
  //damping value d_sigma (see Hansen) is skipped because it is nearly 1
  return nextStepLength;
};

template <class TIndividual>
typename CMAMutation<TIndividual>::VectorType
CMAMutation<TIndividual>::
ConvertObjectivParameterToVector(const IndividualType* pIndividual) const
{
  if (!pIndividual) itkExceptionMacro(<<"Error. Cannot convert objective individual parameters; passed null pointer.");

  VectorType vector = VectorType(pIndividual->ObjectiveParameters().size(),0.0);

  for (int index = 0; index < vector.size(); index++)
  {
    vector[index] += (pIndividual->ObjectiveParameters())[index];
  }
  return vector;
};

template <class TIndividual>
typename CMAMutation<TIndividual>::VectorType
CMAMutation<TIndividual>::
GetRandomVector(unsigned long lSize) const
{
  VectorType z = VectorType(lSize,0.0);

  for (unsigned long index = 0; index<lSize; index++)
  {
    z[index] = m_RandomGenerator->GetVariate();
  }

  return z;
};


} // end of namespace ES
} // end of namespace FREE

#endif

---