28 #ifndef SHARK_ALGORITHMS_DIRECT_SEARCH_VD_CMA_H 29 #define SHARK_ALGORITHMS_DIRECT_SEARCH_VD_CMA_H 55 double chi( std::size_t n ) {
56 return( std::sqrt( static_cast<double>( n ) )*(1. - 1./(4.*n) + 1./(21.*n*n)) );
67 {
return "VDCMA-ES"; }
71 return unsigned( 4. + ::floor( 3. * ::log( static_cast<double>( dimension ) ) ) );
87 double sigma = m_initialSigma;
88 if(m_initialSigma == 0)
89 sigma = 1.0/std::sqrt(
double(p.size()));
108 m_numberOfVariables =
function.numberOfVariables();
111 m_sigma = initialSigma;
114 m_vn.resize(m_numberOfVariables);
115 for(std::size_t i = 0; i != m_numberOfVariables;++i){
116 m_vn(i) =
uni(*mpe_rng,0,1.0/m_numberOfVariables);
122 m_evolutionPathC =
blas::repeat(0.0,m_numberOfVariables);
123 m_evolutionPathSigma =
blas::repeat(0.0,m_numberOfVariables);
126 m_mean = initialSearchPoint;
133 m_weights.resize(m_mu);
134 for (std::size_t i = 0; i < m_mu; i++){
135 m_weights(i) = ::log(mu + 0.5) - ::log(1. + i);
137 m_weights /=
sum(m_weights);
140 m_muEff = 1. /
sum(
sqr(m_weights));
141 m_cSigma = 0.5/(1+std::sqrt(m_numberOfVariables/m_muEff));
142 m_dSigma = 1. + 2. *
std::max(0., std::sqrt((m_muEff-1.)/(m_numberOfVariables+1)) - 1.) + m_cSigma;
144 m_cC = (4. + m_muEff / m_numberOfVariables) / (m_numberOfVariables + 4. + 2 * m_muEff / m_numberOfVariables);
145 double correction = (m_numberOfVariables - 5.0)/6.0;
146 m_c1 = correction*2 / (
sqr(m_numberOfVariables + 1.3) + m_muEff);
147 m_cMu =
std::min(1. - m_c1, correction* 2 * (m_muEff - 2. + 1./m_muEff) / (
sqr(m_numberOfVariables + 2) + m_muEff));
153 std::vector< Individual<RealVector, double, RealVector> > offspring( m_lambda );
156 for( std::size_t i = 0; i < offspring.size(); i++ ) {
157 createSample(offspring[i].searchPoint(),offspring[i].chromosome());
159 penalizingEvaluator(
function, offspring.begin(), offspring.end() );
162 std::vector< Individual<RealVector, double, RealVector> > parents( m_mu );
164 selection(offspring.begin(),offspring.end(),parents.begin(), parents.end());
167 updateStrategyParameters( parents );
191 m_initialSigma = initialSigma;
196 RealVector
const&
mean()
const {
207 return m_evolutionPathC;
212 return m_evolutionPathSigma;
216 std::size_t
mu()
const {
240 RealVector m( m_numberOfVariables, 0. );
241 RealVector z( m_numberOfVariables, 0. );
243 for( std::size_t j = 0; j < offspring.size(); j++ ){
244 noalias(m) += m_weights( j ) * offspring[j].searchPoint();
245 noalias(z) += m_weights( j ) * offspring[j].chromosome();
249 double b=(1/std::sqrt(1+
sqr(m_normv))-1);
253 noalias(m_evolutionPathSigma) = (1. - m_cSigma)*m_evolutionPathSigma + std::sqrt( m_cSigma * (2. - m_cSigma) * m_muEff ) * z;
255 double hSigLHS =
norm_2( m_evolutionPathSigma ) / std::sqrt(1. - pow((1 - m_cSigma), 2.*(m_counter+1)));
256 double hSigRHS = (1.4 + 2 / (m_numberOfVariables+1.)) * chi( m_numberOfVariables );
258 if(hSigLHS < hSigRHS) hSig = 1.;
259 noalias(m_evolutionPathC) = (1. - m_cC ) * m_evolutionPathC + hSig * std::sqrt( m_cC * (2. - m_cC) * m_muEff ) * (m - m_mean) / m_sigma;
268 RealVector meanS(m_numberOfVariables,0.0);
269 RealVector meanT(m_numberOfVariables,0.0);
270 for(std::size_t j = 0; j !=
mu(); ++j){
271 computeSAndTFirst(offspring[j].chromosome(),meanS,meanT,m_cMu*m_weights(j));
273 computeSAndTFirst(m_evolutionPathC/m_D,meanS,meanT,hSig*m_c1);
276 computeSAndTSecond(meanS,meanT);
280 noalias(m_vn) = m_vn*m_normv+meanT/m_normv;
286 m_sigma *= std::exp( (m_cSigma / m_dSigma) * (
norm_2(m_evolutionPathSigma)/ chi( m_numberOfVariables ) - 1.) );
294 void createSample(RealVector& x,RealVector& y)
const{
295 x.resize(m_numberOfVariables);
296 y.resize(m_numberOfVariables);
297 for(std::size_t i = 0; i != m_numberOfVariables; ++i){
298 y(i) =
gauss(*mpe_rng,0,1);
300 double a = std::sqrt(1+
sqr(m_normv))-1;
303 noalias(x) = m_mean+ m_sigma*m_D*y;
309 void computeSAndTFirst(RealVector
const& y, RealVector& s,RealVector& t,
double weight )
const{
310 if(weight == 0)
return;
312 double normv2 =
sqr(m_normv);
313 double gammav = 1+normv2;
315 noalias(s) += weight*(
sqr(y) - (normv2/gammav*yvn)*(y*m_vn)- 1.0);
317 noalias(t) += weight*(yvn*y - 0.5*(
sqr(yvn)+gammav)*m_vn);
321 void computeSAndTSecond(RealVector& s,RealVector& t)
const{
322 RealVector vn2 = m_vn*m_vn;
323 double normv2 =
sqr(m_normv);
324 double gammav = 1+normv2;
326 double alpha =
sqr(normv2)+(2*gammav - std::sqrt(gammav))/
max(vn2);
327 alpha=std::sqrt(alpha);
331 double b=-(1-
sqr(alpha))*
sqr(normv2)/gammav+2*
sqr(alpha);
332 RealVector A= 2.0 - (b+2*
sqr(alpha))*vn2;
333 RealVector invAvn2= vn2/A;
343 std::size_t m_numberOfVariables;
345 std::size_t m_lambda;
347 double m_initialSigma;
357 RealVector m_weights;
359 RealVector m_evolutionPathC;
360 RealVector m_evolutionPathSigma;