146 #define EIGEN_USE_LAPACKE
147 #include "Vector/vector_dist.hpp"
148 #include "DMatrix/EMatrix.hpp"
149 #include <Eigen/Eigenvalues>
150 #include <Eigen/Jacobi>
152 #include "Vector/vector_dist.hpp"
153 #include <f15_cec_fun.hpp>
154 #include <boost/math/special_functions/sign.hpp>
161 constexpr
int dim = 10;
163 constexpr
int lambda = 7;
164 constexpr
int mu = lambda/2;
165 double psoWeight = 0.7;
168 double stopTolX = 2e-11;
169 double stopTolUpX = 2000.0;
171 size_t max_fun_eval = 30000000;
172 constexpr
int hist_size = 21;
181 double stop_fitness = 1.0;
192 constexpr
int sigma = 0;
193 constexpr
int Cov_m = 1;
196 constexpr
int Zeta = 4;
197 constexpr
int path_s = 5;
198 constexpr
int path_c = 6;
199 constexpr
int ord = 7;
200 constexpr
int stop = 8;
201 constexpr
int fithist = 9;
202 constexpr
int weight = 10;
203 constexpr
int validfit = 11;
204 constexpr
int xold = 12;
205 constexpr
int last_restart = 13;
206 constexpr
int iniphase = 14;
207 constexpr
int xmean_st = 15;
208 constexpr
int meanz_st = 16;
213 Eigen::DiagonalMatrix<double,Eigen::Dynamic>,
246 double generateGaussianNoise(
double mu,
double sigma)
248 static const double epsilon = std::numeric_limits<double>::min();
249 static const double two_pi = 2.0*3.14159265358979323846;
251 thread_local
double z1;
252 thread_local
double generate;
253 generate = !generate;
256 {
return z1 * sigma + mu;}
261 u1 = rand() * (1.0 / RAND_MAX);
262 u2 = rand() * (1.0 / RAND_MAX);
264 while ( u1 <= epsilon );
267 z0 = sqrt(-2.0 * log(u2)) * cos(two_pi * u1);
268 z1 = sqrt(-2.0 * log(u2)) * sin(two_pi * u1);
269 return z0 * sigma + mu;
272 template<
unsigned int dim>
273 EVectorXd generateGaussianVector()
278 for (
size_t i = 0 ; i < dim ; i++)
280 tmp(i) = generateGaussianNoise(0,1);
288 template<
unsigned int dim>
289 void fill_vector(
double (& f)[dim], EVectorXd & ev)
291 for (
size_t i = 0 ; i < dim ; i++)
295 void fill_vector(
const double * f, EVectorXd & ev)
297 for (
size_t i = 0 ; i < ev.size() ; i++)
306 bool operator<(
const fun_index & tmp)
const
316 for (
size_t i = 0 ; i < mu ; i++)
317 {wm[i] = log(
double(mu)+1.0) - log(
double(i)+1.0);}
321 for (
size_t i = 0 ; i < mu ; i++)
327 for (
size_t i = 0 ; i < mu ; i++)
339 double weight_sample(
int i)
358 void create_rotmat(EVectorXd & S,EVectorXd & T, EMatrixXd & R)
360 EVectorXd S_work(dim);
361 EVectorXd T_work(dim);
362 EVectorXd S_sup(dim);
363 EVectorXd T_sup(dim);
365 EMatrixXd R_tar(dim,dim);
366 EMatrixXd R_tmp(dim,dim);
367 EMatrixXd R_sup(dim,dim);
369 EMatrixXd S_tmp(2,2);
370 EMatrixXd T_tmp(2,2);
380 for (p = dim - 2; p >= 0 ; p -= 1)
383 for (q = dim - 1 ; q >= p+1 ; q-= 1)
385 T_tmp(0) = T_work(p);
386 T_tmp(1) = T_work(q);
387 S_tmp(0) = S_work(p);
388 S_tmp(1) = S_work(q);
392 Eigen::JacobiRotation<double> G;
394 G.makeGivens(S_tmp(0), S_tmp(1),&z);
403 R_tmp(p,p) = sign*G.c();
404 R_tmp(q,q) = sign*G.c();
405 R_tmp(p,q) = sign*-G.s();
406 R_tmp(q,p) = sign*G.s();
411 S_work = R_tmp*S_work;
417 G.makeGivens(T_tmp(0), T_tmp(1),&z);
424 R_tmp(p,p) = sign*G.c();
425 R_tmp(q,q) = sign*G.c();
426 R_tmp(p,q) = sign*-G.s();
427 R_tmp(q,p) = sign*G.s();
431 T_work = R_tmp*T_work;
436 R = R_tar.transpose()*R;
440 EVectorXd Check(dim);
467 Eigen::DiagonalMatrix<double,Eigen::Dynamic> &
D,
470 EVectorXd best_sol_ei(dim);
472 double bias_weight = psoWeight;
473 fill_vector(&best_sol.get(0),best_sol_ei);
474 EVectorXd gb_vec = best_sol_ei-xmean;
475 double gb_vec_length = sqrt(gb_vec.transpose() * gb_vec);
476 EVectorXd b_main = B.col(dim-1);
481 EMatrixXd R(dim,dim);
483 if (gb_vec_length > 0.0)
485 if(sigma < gb_vec_length)
487 if(sigma/gb_vec_length <= t_c*gb_vec_length)
490 {bias = sigma*gb_vec/gb_vec_length;}
496 xmean = xmean + bias;
500 EMatrixXd B_rot(dim,dim);
501 Eigen::DiagonalMatrix<double,Eigen::Dynamic> D_square(dim);
503 EVectorXd gb_vec_old = best_sol_ei - xold;
504 create_rotmat(b_main,gb_vec_old,R);
505 for (
size_t i = 0 ; i < dim ; i++)
506 {B_rot.col(i) = R*B.col(i);}
508 for (
size_t i = 0 ; i < dim ; i++)
509 {D_square.diagonal()[i] = D.diagonal()[i] * D.diagonal()[i];}
510 C_pso = B_rot * D_square * B_rot.transpose();
512 EMatrixXd trUp = C_pso.triangularView<Eigen::Upper>();
513 EMatrixXd trDw = C_pso.triangularView<Eigen::StrictlyUpper>();
514 C_pso = trUp + trDw.transpose();
544 best_sol.resize(dim);
545 auto & v_cl = create_vcluster();
547 double best_old = best_sample;
548 v_cl.min(best_sample);
552 if (best < best_sample)
558 if (best_old == best_sample)
560 rank = v_cl.getProcessUnitID();
566 if (rank == v_cl.getProcessUnitID())
568 for (
size_t i = 0 ; i < dim ; i++)
569 {best_sol.get(i) = best_sample_sol.get(i);}
574 rank = std::numeric_limits<size_t>::max();
583 v_cl.Bcast(best_sol,rank);
613 double availablepercentiles[lambda];
617 for (
size_t i = 0 ; i < lambda ; i++)
619 availablepercentiles[i] = 0.0;
620 sar[i] = f_obj.get(i).f;
622 std::sort(&sar[0],&sar[lambda]);
624 for (
size_t i = 0 ; i < 2 ; i++)
626 if (perc[i] <= (100.0*0.5/lambda))
628 else if (perc[i] >= (100.0*(lambda-0.5)/lambda) )
629 {res[i] = sar[lambda-1];}
632 for (
size_t j = 0 ; j < lambda ; j++)
633 {availablepercentiles[j] = 100.0 * ((double(j)+1.0)-0.5) / lambda;}
635 for (k = 0 ; k < lambda ; k++)
636 {
if(availablepercentiles[k] >= perc[i]) {
break;}}
639 res[i] = sar[k] + (sar[k+1]-sar[k]) * (perc[i]
640 -availablepercentiles[k]) / (availablepercentiles[k+1] - availablepercentiles[k]);
647 double maxval(
double (& buf)[hist_size],
bool (& mask)[hist_size])
650 for (
size_t i = 0 ; i < hist_size ; i++)
652 if (buf[i] > max && mask[i] ==
true)
659 double minval(
double (& buf)[hist_size],
bool (& mask)[hist_size])
661 double min = std::numeric_limits<double>::max();
662 for (
size_t i = 0 ; i < hist_size ; i++)
664 if (buf[i] < min && mask[i] ==
true)
673 void cmaes_intobounds(EVectorXd & x, EVectorXd & xout,
bool (& idx)[dim],
bool & idx_any)
676 for (
size_t i = 0; i < dim ; i++)
705 EVectorXd (& arxvalid)[lambda],
706 EVectorXd (& arx)[lambda],
710 double (& weight)[dim],
711 double (& fithist)[hist_size],
713 double & validfitval,
723 bool mask[hist_size];
726 int dfitidx[hist_size];
727 double dfitsort[hist_size];
728 double prct[2] = {25.0,75.0};
731 for (
size_t i = 0 ; i < hist_size ; i++)
736 if (fithist[i] > 0.0)
742 for (
size_t i = 0 ; i < dim ; i++)
747 meandiag = C.trace()/dim;
749 cmaes_myprctile(f_obj, prct, val);
750 value = (val[1] - val[0]) / dim / meandiag / (sigma*sigma);
752 if (value >= std::numeric_limits<double>::max())
754 auto & v_cl = create_vcluster();
755 std::cout <<
"Process " << v_cl.rank() <<
" warning: Non-finite fitness range" << std::endl;
756 value = maxval(fithist,mask);
758 else if(value == 0.0)
760 value = minval(fithist,mask);
762 else if (validfit == 0.0)
764 for (
size_t i = 0 ; i < hist_size ; i++)
769 for (
size_t i = 0; i < hist_size ; i++)
777 else if(i == hist_size-1)
779 for (
size_t k = 0 ; k < hist_size-1 ; k++)
780 {fithist[k] = fithist[k+1];}
786 cmaes_intobounds(xmean,tx,idx,idx_any);
793 {value = fithist[0];}
797 for (
size_t i = 0 ; i <= maxI; i++)
799 fitsort.get(i).f = fithist[i];
800 fitsort.get(i).id = i;
804 for (
size_t k = 0; k <= maxI ; k++)
805 {fitsort.get(k).f = fithist[fitsort.get(k).id];}
807 if ((maxI+1) % 2 == 0)
808 {value = (fitsort.get(maxI/2).f+fitsort.get(maxI/2+1).f)/2.0;}
810 {value = fitsort.get(maxI/2).f;}
812 for (
size_t i = 0 ; i < dim ; i++)
814 diag[i] = diag[i]/meandiag;
815 weight[i] = 2.0002 * value / diag[i];
817 if (validfitval == 1.0 && step-last_restart > 2)
827 for(
size_t i = 0 ; i < dim ; i++)
829 idx[i] = (idx[i] && (fabs(tx(i)) > 3.0*std::max(1.0,sqrt(dim)/mu_eff) * sigma * sqrt(diag[i])));
830 idx[i] = (idx[i] && (std::copysign(1.0,tx(i)) == std::copysign(1.0,(xmean(i)-xold(i)))) );
832 for (
size_t i = 0 ; i < dim ; i++)
836 weight[i] = pow(1.2,(std::max(1.0,mu_eff/10.0/dim)))*weight[i];
840 double arpenalty[lambda];
841 for (
size_t i = 0 ; i < lambda ; i++)
844 for (
size_t j = 0 ; j < dim ; j++)
846 arpenalty[i] += weight[j] * (arxvalid[i](j) - arx[i](j))*(arxvalid[i](j) - arx[i](j));
848 f_obj.get(i).f += arpenalty[i];
855 double adjust_sigma(
double sigma, EMatrixXd & C)
857 for (
size_t i = 0 ; i < dim ; i++)
859 if (sigma*sqrt(C(i,i)) > 5.0)
860 {sigma = 5.0/sqrt(C(i,i));}
890 EVectorXd xmean(dim);
891 EVectorXd mean_z(dim);
892 EVectorXd arxvalid[lambda];
893 EVectorXd arx[lambda];
895 for (
size_t i = 0 ; i < lambda ; i++)
898 arxvalid[i].resize(dim);
901 double best_sample = std::numeric_limits<double>::max();
906 int counteval = step*lambda;
913 if (vd.
getProp<stop>(p) ==
true)
916 EVectorXd (& arz)[lambda] = vd.
getProp<Zeta>(p);
920 fill_vector(vd.
getPos(p),xmean);
922 for (
size_t j = 0 ; j < lambda ; j++)
924 vd.
getProp<Zeta>(p)[j] = generateGaussianVector<dim>();
928 for (
size_t i = 0 ; i < dim ; i++)
930 if (arx[j](i) < -5.0)
931 {arxvalid[j](i) = -5.0;}
932 else if (arx[j](i) > 5.0)
933 {arxvalid[j](i) = 5.0;}
935 {arxvalid[j](i) = arx[j](i);}
938 f_obj.get(j).f = hybrid_composition<dim>(arxvalid[j]);
943 if (f_obj.get(j).f < best_sample)
945 best_sample = f_obj.get(j).f;
948 for (
size_t i = 0 ; i < dim ; i++)
949 {best_sample_sol.get(i) = arxvalid[j](i);}
954 cmaes_handlebounds(f_obj,vd.
getProp<sigma>(p),
955 vd.
getProp<validfit>(p),arxvalid,
959 vd.
getProp<validfit>(p),mu_eff,
960 step,vd.
getProp<last_restart>(p));
964 for (
size_t j = 0 ; j < lambda ; j++)
965 {vd.
getProp<ord>(p)[j] = f_obj.get(j).id;}
973 for (
size_t j = 0 ; j < mu ; j++)
975 xmean += weight_sample(j)*arx[vd.
getProp<ord>(p)[j]];
976 mean_z += weight_sample(j)*vd.
getProp<Zeta>(p)[vd.
getProp<ord>(p)[j]];
979 vd.
getProp<xmean_st>(p) = xmean;
980 vd.
getProp<meanz_st>(p) = mean_z;
986 broadcast_best_solution(vd,best_sol,best,best_sample,best_sample_sol);
989 bool calc_bd = counteval - eigeneval > lambda/(ccov)/dim/10;
991 {eigeneval = counteval;}
998 if (vd.
getProp<stop>(p) ==
true)
1001 xmean = vd.
getProp<xmean_st>(p);
1002 mean_z = vd.
getProp<meanz_st>(p);
1004 vd.
getProp<path_s>(p) = vd.
getProp<path_s>(p)*(1.0 - cs) + sqrt(cs*(2.0-cs)*mu_eff)*vd.
getProp<B>(p)*mean_z;
1006 double hsig = vd.
getProp<path_s>(p).norm()/sqrt(1.0-pow((1.0-cs),(2.0*
double((step-vd.
getProp<last_restart>(p))))))/chiN < 1.4 + 2.0/(dim+1);
1010 if (step % N_pso == 0)
1012 EMatrixXd C_pso(dim,dim);
1016 vd.
getProp<Cov_m>(p) = (1.0-ccov+(1.0-hsig)*ccov*cc*(2.0-cc)/mu_eff)*vd.
getProp<Cov_m>(p) +
1017 ccov*(1.0/mu_eff)*(vd.
getProp<path_c>(p)*vd.
getProp<path_c>(p).transpose());
1019 for (
size_t i = 0 ; i < mu ; i++)
1024 vd.
getProp<Cov_m>(p) = psoWeight*vd.
getProp<Cov_m>(p) + (1.0 - psoWeight)*C_pso;
1029 vd.
getProp<Cov_m>(p) = (1.0-ccov+(1.0-hsig)*ccov*cc*(2.0-cc)/mu_eff)*vd.
getProp<Cov_m>(p) +
1030 ccov*(1.0/mu_eff)*(vd.
getProp<path_c>(p)*vd.
getProp<path_c>(p).transpose());
1032 for (
size_t i = 0 ; i < mu ; i++)
1040 double smaller = std::numeric_limits<double>::max();
1041 for (
size_t i = 0 ; i < dim ; i++)
1043 if (vd.
getProp<sigma>(p)*sqrt(vd.
getProp<D>(p).diagonal()[i]) > 5.0)
1045 if (smaller > 5.0/sqrt(vd.
getProp<D>(p).diagonal()[i]))
1046 {smaller = 5.0/sqrt(vd.
getProp<D>(p).diagonal()[i]);}
1049 if (smaller != std::numeric_limits<double>::max())
1050 {vd.
getProp<sigma>(p) = smaller;}
1053 vd.
getProp<sigma>(p) = vd.
getProp<sigma>(p)*exp((cs/d_amps)*(vd.
getProp<path_s>(p).norm()/chiN - 1));
1059 EMatrixXd trUp = vd.
getProp<Cov_m>(p).triangularView<Eigen::Upper>();
1060 EMatrixXd trDw = vd.
getProp<Cov_m>(p).triangularView<Eigen::StrictlyUpper>();
1061 vd.
getProp<Cov_m>(p) = trUp + trDw.transpose();
1064 Eigen::SelfAdjointEigenSolver<EMatrixXd> eig_solver;
1066 eig_solver.compute(vd.
getProp<Cov_m>(p));
1068 for (
size_t i = 0 ; i < eig_solver.eigenvalues().size() ; i++)
1069 {vd.
getProp<D>(p).diagonal()[i] = sqrt(eig_solver.eigenvalues()[i]);}
1070 vd.
getProp<B>(p) = eig_solver.eigenvectors();
1073 for (
size_t i = 0 ; i < dim ; i++)
1075 if (vd.
getProp<B>(p)(0,i) < 0)
1079 EMatrixXd tmp = vd.
getProp<B>(p).transpose();
1083 for (
size_t i = 0 ; i < dim ; i++)
1084 {vd.
getPos(p)[i] = xmean(i);}
1089 bool stop_tol =
true;
1090 bool stop_tolX =
true;
1091 for (
size_t i = 0 ; i < dim ; i++)
1093 stop_tol &= (vd.
getProp<sigma>(p)*std::max(fabs(vd.
getProp<path_c>(p)(i)),sqrt(vd.
getProp<Cov_m>(p)(i,i)))) < stopTolX;
1094 stop_tolX &= vd.
getProp<sigma>(p)*sqrt(vd.
getProp<D>(p).diagonal()[i]) > stopTolUpX;
1097 vd.
getProp<stop>(p) = stop_tol | stop_tolX;
1100 if (f_obj.get(0).f == f_obj.get(std::ceil(0.7*lambda)).f )
1103 std::cout <<
"warning: flat fitness, consider reformulating the objective";
1109 if (vd.
getProp<stop>(p) ==
true)
1110 {std::cout <<
"Stopped" << std::endl;}
1112 if (restart_cma && vd.
getProp<stop>(p) ==
true)
1114 std::cout <<
"------- Restart #" << std::endl;
1116 std::cout <<
"---------------------------------" << std::endl;
1117 std::cout <<
"Best: " << best <<
" " << fun_eval << std::endl;
1118 std::cout <<
"---------------------------------" << std::endl;
1120 vd.
getProp<last_restart>(p) = step;
1121 vd.
getProp<xold>(p).setZero();
1123 for (
size_t i = 0 ; i < vd.
getProp<D>(p).diagonal().size() ; i++)
1124 {vd.
getProp<D>(p).diagonal()[i] = 1.0;}
1125 vd.
getProp<B>(p).resize(dim,dim);
1126 vd.
getProp<B>(p).setIdentity();
1128 vd.
getProp<path_s>(p).resize(dim);
1129 vd.
getProp<path_s>(p).setZero(dim);
1130 vd.
getProp<path_c>(p).resize(dim);
1131 vd.
getProp<path_c>(p).setZero(dim);
1133 vd.
getProp<iniphase>(p) =
true;
1134 vd.
getProp<last_restart>(p) = 0;
1138 for (
size_t i = 0 ; i < dim ; i++)
1141 vd.
getPos(p)[i] = 10.0*(double)rand() / RAND_MAX - 5.0;
1146 for (
size_t i = 0 ; i < hist_size ; i++)
1147 {vd.
getProp<fithist>(p)[i] = -1.0;}
1148 vd.
getProp<fithist>(p)[0] = 1.0;
1149 vd.
getProp<validfit>(p) = 0.0;
1155 auto & v_cl = create_vcluster();
1177 int main(
int argc,
char* argv[])
1180 openfpm_init(&argc,&argv);
1182 auto & v_cl = create_vcluster();
1187 for (
size_t i = 0 ; i < dim ; i++)
1195 for (
size_t i = 0 ; i < dim ; i++)
1196 {bc[i] = NON_PERIODIC;};
1207 stop_fitness = 1e-10;
1208 size_t stopeval = 1e3*dim*dim;
1214 cc = 4.0 / (dim+4.0);
1215 cs = (mu_eff+2.0) / (double(dim)+mu_eff+3.0);
1216 ccov = (1.0/mu_eff) * 2.0/((dim+1.41)*(dim+1.41)) +
1217 (1.0 - 1.0/mu_eff)* std::min(1.0,(2.0*mu_eff-1.0)/((dim+2.0)*(dim+2.0) + mu_eff));
1218 d_amps = 1 + 2*std::max(0.0, sqrt((mu_eff-1.0)/(dim+1))-1) + cs;
1220 chiN = sqrt(dim)*(1.0-1.0/(4.0*dim)+1.0/(21.0*dim*dim));
1226 int seed = 24756*v_cl.rank()*v_cl.rank() + time(NULL);
1235 for (
size_t i = 0 ; i < dim ; i++)
1238 vd.
getPos(p)[i] = 10.0*(double)rand() / RAND_MAX - 5.0;
1246 for (
size_t i = 0 ; i < vd.
getProp<D>(p).diagonal().size() ; i++)
1247 {vd.
getProp<D>(p).diagonal()[i] = 1.0;}
1248 vd.
getProp<B>(p).resize(dim,dim);
1249 vd.
getProp<B>(p).setIdentity();
1251 vd.
getProp<path_s>(p).resize(dim);
1252 vd.
getProp<path_s>(p).setZero(dim);
1253 vd.
getProp<path_c>(p).resize(dim);
1254 vd.
getProp<path_c>(p).setZero(dim);
1256 vd.
getProp<iniphase>(p) =
true;
1257 vd.
getProp<last_restart>(p) = 0;
1261 for (
size_t i = 0 ; i < hist_size ; i++)
1262 {vd.
getProp<fithist>(p)[i] = -1.0;}
1263 vd.
getProp<fithist>(p)[0] = 1.0;
1264 vd.
getProp<validfit>(p) = 0.0;
1270 if (v_cl.rank() == 0)
1271 {std::cout <<
"Starting PS-CMA-ES" << std::endl;}
1276 best = std::numeric_limits<double>::max();
1281 size_t fun_eval = 0;
1283 while (fun_eval < max_fun_eval && best > 120.000001)
1286 cma_step(vd,i+1,best,best_i,best_sol,fun_eval);
1291 if (v_cl.rank() == 0)
1293 std::cout <<
"Best solution: " << best <<
" with " << fun_eval << std::endl;
1294 std::cout <<
"at: " << std::endl;
1296 for (
size_t i = 0 ; i < best_sol.
size() ; i++)
1298 std::cout << best_sol.get(i) <<
" ";
Derivative second order on h (spacing)
auto getProp(vect_dist_key_dx vec_key) -> decltype(v_prp.template get< id >(vec_key.getKey()))
Get the property of an element.
auto getPos(vect_dist_key_dx vec_key) -> decltype(v_pos.template get< 0 >(vec_key.getKey()))
Get the position of an element.
void setHigh(int i, T val)
set the high interval of the box
void setLow(int i, T val)
set the low interval of the box
This class represent an N-dimensional box.
vector_dist_iterator getDomainIterator() const
Get an iterator that traverse the particles in the domain.
vect_dist_key_dx get()
Get the actual key.
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
It model an expression expr1 + ... exprn.
Implementation of 1-D std::vector like structure.