OpenFPM_pdata  1.1.0
Project that contain the implementation of distributed structures
 All Data Structures Namespaces Functions Variables Typedefs Enumerations Friends Pages
main.cpp
1 
144 
146 #define EIGEN_USE_LAPACKE
147 #include "Vector/vector_dist.hpp"
148 #include "DMatrix/EMatrix.hpp"
149 #include <Eigen/Eigenvalues>
150 #include <Eigen/Jacobi>
151 #include <limits>
152 #include "Vector/vector_dist.hpp"
153 #include <f15_cec_fun.hpp>
154 #include <boost/math/special_functions/sign.hpp>
155 
157 
159 
160 // PARAMETERS
161 constexpr int dim = 10;
162 // when you set dim set also lambda to to 4+std::floor(3*log(dim))
163 constexpr int lambda = 7;
164 constexpr int mu = lambda/2;
165 double psoWeight = 0.7;
166 // number of cma-step before pso step
167 int N_pso = 200;
168 double stopTolX = 2e-11;
169 double stopTolUpX = 2000.0;
170 int restart_cma = 1;
171 size_t max_fun_eval = 30000000;
172 constexpr int hist_size = 21;
173 
174 // Convenient global variables (Their value is set after)
175 double mu_eff = 1.0;
176 double cs = 1.0;
177 double cc = 1.0;
178 double ccov = 1.0;
179 double chiN;
180 double d_amps = 1.0;
181 double stop_fitness = 1.0;
182 int eigeneval = 0;
183 double t_c = 0.1;
184 double b = 0.1;
185 
187 
189 
191 
192 constexpr int sigma = 0;
193 constexpr int Cov_m = 1;
194 constexpr int B = 2;
195 constexpr int D = 3;
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;
209 
210 typedef vector_dist<dim,double, aggregate<double,
211  EMatrixXd,
212  EMatrixXd,
213  Eigen::DiagonalMatrix<double,Eigen::Dynamic>,
214  EVectorXd[lambda],
215  EVectorXd,
216  EVectorXd,
217  int[lambda],
218  int,
219  double [hist_size],
220  double [dim],
221  double,
222  EVectorXd,
223  int,
224  bool,
225  EVectorXd,
226  EVectorXd> > particle_type;
227 
229 
244 
246 double generateGaussianNoise(double mu, double sigma)
247 {
248  static const double epsilon = std::numeric_limits<double>::min();
249  static const double two_pi = 2.0*3.14159265358979323846;
250 
251  thread_local double z1;
252  thread_local double generate;
253  generate = !generate;
254 
255  if (!generate)
256  {return z1 * sigma + mu;}
257 
258  double u1, u2;
259  do
260  {
261  u1 = rand() * (1.0 / RAND_MAX);
262  u2 = rand() * (1.0 / RAND_MAX);
263  }
264  while ( u1 <= epsilon );
265 
266  double z0;
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;
270 }
271 
272 template<unsigned int dim>
273 EVectorXd generateGaussianVector()
274 {
275  EVectorXd tmp;
276  tmp.resize(dim);
277 
278  for (size_t i = 0 ; i < dim ; i++)
279  {
280  tmp(i) = generateGaussianNoise(0,1);
281  }
282 
283  return tmp;
284 }
285 
287 
288 template<unsigned int dim>
289 void fill_vector(double (& f)[dim], EVectorXd & ev)
290 {
291  for (size_t i = 0 ; i < dim ; i++)
292  {ev(i) = f[i];}
293 }
294 
295 void fill_vector(const double * f, EVectorXd & ev)
296 {
297  for (size_t i = 0 ; i < ev.size() ; i++)
298  {ev(i) = f[i];}
299 }
300 
301 struct fun_index
302 {
303  double f;
304  int id;
305 
306  bool operator<(const fun_index & tmp) const
307  {
308  return f < tmp.f;
309  }
310 };
311 
312 double wm[mu];
313 
314 void init_weight()
315 {
316  for (size_t i = 0 ; i < mu ; i++)
317  {wm[i] = log(double(mu)+1.0) - log(double(i)+1.0);}
318 
319  double tot = 0.0;
320 
321  for (size_t i = 0 ; i < mu ; i++)
322  {tot += wm[i];}
323 
324  double sum = 0.0;
325  double sum2 = 0.0;
326 
327  for (size_t i = 0 ; i < mu ; i++)
328  {
329  wm[i] /= tot;
330  sum += wm[i];
331  sum2 += wm[i]*wm[i];
332  }
333 
334  // also set mu_eff
335  mu_eff=sum*sum/sum2;
336 
337 }
338 
339 double weight_sample(int i)
340 {
341  return wm[i];
342 }
343 
356 
358 void create_rotmat(EVectorXd & S,EVectorXd & T, EMatrixXd & R)
359 {
360  EVectorXd S_work(dim);
361  EVectorXd T_work(dim);
362  EVectorXd S_sup(dim);
363  EVectorXd T_sup(dim);
364 
365  EMatrixXd R_tar(dim,dim);
366  EMatrixXd R_tmp(dim,dim);
367  EMatrixXd R_sup(dim,dim);
368  double G_S,G_C;
369  EMatrixXd S_tmp(2,2);
370  EMatrixXd T_tmp(2,2);
371  int p,q,i;
372 
373  S_work = S;
374  T_work = T;
375 
376  R.setIdentity();
377  R_tar = R;
378  R_tmp = R;
379 
380  for (p = dim - 2; p >= 0 ; p -= 1)
381  {
382 
383  for (q = dim - 1 ; q >= p+1 ; q-= 1)
384  {
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);
389 
390  // Perform Givens Rotation on start vector
391 
392  Eigen::JacobiRotation<double> G;
393  double z;
394  G.makeGivens(S_tmp(0), S_tmp(1),&z);
395 
396  // Check direction of rotation
397  double sign = 1.0;
398  if (z < 0.0)
399  {sign = -1.0;}
400 
401  // Build a Rotation Matrix out of G_C and G_S
402  R_tmp.setIdentity();
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();
407 
408  // Rotate start vector and update R
409  // S_work = R_tmp*S_work
410 
411  S_work = R_tmp*S_work;
412  // R = R_tmp*R
413  R = R_tmp*R;
414 
415  // Perform Givens Rotation on target vector
416 
417  G.makeGivens(T_tmp(0), T_tmp(1),&z);
418 
419  sign = 1.0;
420  if (z < 0.0)
421  {sign = -1.0;}
422 
423  R_tmp.setIdentity();
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();
428 
429  // Rotate target vector and update R_tar
430 
431  T_work = R_tmp*T_work;
432  R_tar = R_tmp*R_tar;
433  }
434  }
435 
436  R = R_tar.transpose()*R;
437 
438  // Check the rotation
439 
440  EVectorXd Check(dim);
441  Check = R*S;
442 }
443 
445 
460 
462 void updatePso(openfpm::vector<double> & best_sol,
463  double sigma,
464  EVectorXd & xmean,
465  EVectorXd & xold,
466  EMatrixXd & B,
467  Eigen::DiagonalMatrix<double,Eigen::Dynamic> & D,
468  EMatrixXd & C_pso)
469 {
470  EVectorXd best_sol_ei(dim);
471 
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);
477  EVectorXd bias(dim);
478  bias.setZero();
479 
480  // Rotation Matrix
481  EMatrixXd R(dim,dim);
482 
483  if (gb_vec_length > 0.0)
484  {
485  if(sigma < gb_vec_length)
486  {
487  if(sigma/gb_vec_length <= t_c*gb_vec_length)
488  {bias = 0.5*gb_vec;}
489  else
490  {bias = sigma*gb_vec/gb_vec_length;}
491  }
492  else
493  {bias.setZero();}
494  }
495 
496  xmean = xmean + bias;
497 
498  if (psoWeight < 1.0)
499  {
500  EMatrixXd B_rot(dim,dim);
501  Eigen::DiagonalMatrix<double,Eigen::Dynamic> D_square(dim);
502 
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);}
507 
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();
511 
512  EMatrixXd trUp = C_pso.triangularView<Eigen::Upper>();
513  EMatrixXd trDw = C_pso.triangularView<Eigen::StrictlyUpper>();
514  C_pso = trUp + trDw.transpose();
515  }
516 }
517 
519 
536 
538 void broadcast_best_solution(particle_type & vd,
539  openfpm::vector<double> & best_sol,
540  double & best,
541  double best_sample,
542  openfpm::vector<double> & best_sample_sol)
543 {
544  best_sol.resize(dim);
545  auto & v_cl = create_vcluster();
546 
547  double best_old = best_sample;
548  v_cl.min(best_sample);
549  v_cl.execute();
550 
551  // The old solution remain the best
552  if (best < best_sample)
553  {return;}
554 
555  best = best_sample;
556 
557  size_t rank;
558  if (best_old == best_sample)
559  {
560  rank = v_cl.getProcessUnitID();
561 
562  // we own the minimum and we decide who broad cast
563  v_cl.min(rank);
564  v_cl.execute();
565 
566  if (rank == v_cl.getProcessUnitID())
567  {
568  for (size_t i = 0 ; i < dim ; i++)
569  {best_sol.get(i) = best_sample_sol.get(i);}
570  }
571  }
572  else
573  {
574  rank = std::numeric_limits<size_t>::max();
575 
576  // we do not own decide who broad cast
577  v_cl.min(rank);
578  v_cl.execute();
579  }
580 
581  // now we broad cast the best solution across processors
582 
583  v_cl.Bcast(best_sol,rank);
584  v_cl.execute();
585 }
586 
588 
608 
610 void cmaes_myprctile(openfpm::vector<fun_index> & f_obj, double (& perc)[2], double (& res)[2])
611 {
612  double sar[lambda];
613  double availablepercentiles[lambda];
614  int idx[hist_size];
615  int i,k;
616 
617  for (size_t i = 0 ; i < lambda ; i++)
618  {
619  availablepercentiles[i] = 0.0;
620  sar[i] = f_obj.get(i).f;
621  }
622  std::sort(&sar[0],&sar[lambda]);
623 
624  for (size_t i = 0 ; i < 2 ; i++)
625  {
626  if (perc[i] <= (100.0*0.5/lambda))
627  {res[i] = sar[0];}
628  else if (perc[i] >= (100.0*(lambda-0.5)/lambda) )
629  {res[i] = sar[lambda-1];}
630  else
631  {
632  for (size_t j = 0 ; j < lambda ; j++)
633  {availablepercentiles[j] = 100.0 * ((double(j)+1.0)-0.5) / lambda;}
634 
635  for (k = 0 ; k < lambda ; k++)
636  {if(availablepercentiles[k] >= perc[i]) {break;}}
637  k-=1;
638 
639  res[i] = sar[k] + (sar[k+1]-sar[k]) * (perc[i]
640  -availablepercentiles[k]) / (availablepercentiles[k+1] - availablepercentiles[k]);
641  }
642  }
643 }
644 
646 
647 double maxval(double (& buf)[hist_size], bool (& mask)[hist_size])
648 {
649  double max = 0.0;
650  for (size_t i = 0 ; i < hist_size ; i++)
651  {
652  if (buf[i] > max && mask[i] == true)
653  {max = buf[i];}
654  }
655 
656  return max;
657 }
658 
659 double minval(double (& buf)[hist_size], bool (& mask)[hist_size])
660 {
661  double min = std::numeric_limits<double>::max();
662  for (size_t i = 0 ; i < hist_size ; i++)
663  {
664  if (buf[i] < min && mask[i] == true)
665  {min = buf[i];}
666  }
667 
668  return min;
669 }
670 
672 
673 void cmaes_intobounds(EVectorXd & x, EVectorXd & xout,bool (& idx)[dim], bool & idx_any)
674 {
675  idx_any = false;
676  for (size_t i = 0; i < dim ; i++)
677  {
678  if(x(i) < -5.0)
679  {
680  xout(i) = -5.0;
681  idx[i] = true;
682  idx_any = true;
683  }
684  else if (x(i) > 5.0)
685  {
686  xout(i) = 5.0;
687  idx[i] = true;
688  idx_any = true;
689  }
690  else
691  {
692  xout(i) = x(i);
693  idx[i] = false;
694  }
695  }
696 }
697 
699 
701 
702 void cmaes_handlebounds(openfpm::vector<fun_index> & f_obj,
703  double sigma,
704  double & validfit,
705  EVectorXd (& arxvalid)[lambda],
706  EVectorXd (& arx)[lambda],
707  EMatrixXd & C,
708  EVectorXd & xmean,
709  EVectorXd & xold,
710  double (& weight)[dim],
711  double (& fithist)[hist_size],
712  bool & iniphase,
713  double & validfitval,
714  double mu_eff,
715  int step,
716  int last_restart)
717 {
718  double val[2];
719  double value;
720  double diag[dim];
721  double meandiag;
722  int i,k,maxI;
723  bool mask[hist_size];
724  bool idx[dim];
725  EVectorXd tx(dim);
726  int dfitidx[hist_size];
727  double dfitsort[hist_size];
728  double prct[2] = {25.0,75.0};
729  bool idx_any;
730 
731  for (size_t i = 0 ; i < hist_size ; i++)
732  {
733  dfitsort[i] = 0.0;
734  dfitidx[i] = 0;
735 
736  if (fithist[i] > 0.0)
737  {mask[i] = true;}
738  else
739  {mask[i] = false;}
740  }
741 
742  for (size_t i = 0 ; i < dim ; i++)
743  {diag[i] = C(i,i);}
744 
745  maxI = 0;
746 
747  meandiag = C.trace()/dim;
748 
749  cmaes_myprctile(f_obj, prct, val);
750  value = (val[1] - val[0]) / dim / meandiag / (sigma*sigma);
751 
752  if (value >= std::numeric_limits<double>::max())
753  {
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);
757  }
758  else if(value == 0.0)
759  {
760  value = minval(fithist,mask);
761  }
762  else if (validfit == 0.0)
763  {
764  for (size_t i = 0 ; i < hist_size ; i++)
765  {fithist[i] = -1.0;}
766  validfit = 1;
767  }
768 
769  for (size_t i = 0; i < hist_size ; i++)
770  {
771  if(fithist[i] < 0.0)
772  {
773  fithist[i] = value;
774  maxI = i;
775  break;
776  }
777  else if(i == hist_size-1)
778  {
779  for (size_t k = 0 ; k < hist_size-1 ; k++)
780  {fithist[k] = fithist[k+1];}
781  fithist[i] = value;
782  maxI = i;
783  }
784  }
785 
786  cmaes_intobounds(xmean,tx,idx,idx_any);
787 
788  if (iniphase)
789  {
790  if (idx_any)
791  {
792  if(maxI == 0)
793  {value = fithist[0];}
794  else
795  {
796  openfpm::vector<fun_index> fitsort(maxI+1);
797  for (size_t i = 0 ; i <= maxI; i++)
798  {
799  fitsort.get(i).f = fithist[i];
800  fitsort.get(i).id = i;
801  }
802 
803  fitsort.sort();
804  for (size_t k = 0; k <= maxI ; k++)
805  {fitsort.get(k).f = fithist[fitsort.get(k).id];}
806 
807  if ((maxI+1) % 2 == 0)
808  {value = (fitsort.get(maxI/2).f+fitsort.get(maxI/2+1).f)/2.0;}
809  else
810  {value = fitsort.get(maxI/2).f;}
811  }
812  for (size_t i = 0 ; i < dim ; i++)
813  {
814  diag[i] = diag[i]/meandiag;
815  weight[i] = 2.0002 * value / diag[i];
816  }
817  if (validfitval == 1.0 && step-last_restart > 2)
818  {
819  iniphase = false;
820  }
821  }
822  }
823 
824  if(idx_any)
825  {
826  tx = xmean - tx;
827  for(size_t i = 0 ; i < dim ; i++)
828  {
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)))) );
831  }
832  for (size_t i = 0 ; i < dim ; i++)
833  {
834  if (idx[i] == true)
835  {
836  weight[i] = pow(1.2,(std::max(1.0,mu_eff/10.0/dim)))*weight[i];
837  }
838  }
839  }
840  double arpenalty[lambda];
841  for (size_t i = 0 ; i < lambda ; i++)
842  {
843  arpenalty[i] = 0.0;
844  for (size_t j = 0 ; j < dim ; j++)
845  {
846  arpenalty[i] += weight[j] * (arxvalid[i](j) - arx[i](j))*(arxvalid[i](j) - arx[i](j));
847  }
848  f_obj.get(i).f += arpenalty[i];
849  }
850 // fitness%sel = fitness%raw + bnd%arpenalty;
851 }
852 
854 
855 double adjust_sigma(double sigma, EMatrixXd & C)
856 {
857  for (size_t i = 0 ; i < dim ; i++)
858  {
859  if (sigma*sqrt(C(i,i)) > 5.0)
860  {sigma = 5.0/sqrt(C(i,i));}
861  }
862 
863  return sigma;
864 }
865 
883 
885 void cma_step(particle_type & vd, int step, double & best,
886  int & best_i, openfpm::vector<double> & best_sol,
887  size_t & fun_eval)
888 {
889  size_t fe = 0;
890  EVectorXd xmean(dim);
891  EVectorXd mean_z(dim);
892  EVectorXd arxvalid[lambda];
893  EVectorXd arx[lambda];
894 
895  for (size_t i = 0 ; i < lambda ; i++)
896  {
897  arx[i].resize(dim);
898  arxvalid[i].resize(dim);
899  }
900 
901  double best_sample = std::numeric_limits<double>::max();
902  openfpm::vector<double> best_sample_sol(dim);
903 
904  openfpm::vector<fun_index> f_obj(lambda);
905 
906  int counteval = step*lambda;
907 
908  auto it = vd.getDomainIterator();
909  while (it.isNext())
910  {
911  auto p = it.get();
912 
913  if (vd.getProp<stop>(p) == true)
914  {++it;continue;}
915 
916  EVectorXd (& arz)[lambda] = vd.getProp<Zeta>(p);
917 
918  // fill the mean vector;
919 
920  fill_vector(vd.getPos(p),xmean);
921 
922  for (size_t j = 0 ; j < lambda ; j++)
923  {
924  vd.getProp<Zeta>(p)[j] = generateGaussianVector<dim>();
925  arx[j] = xmean + vd.getProp<sigma>(p)*vd.getProp<B>(p)*vd.getProp<D>(p)*vd.getProp<Zeta>(p)[j];
926 
927  // sample point has to be inside -5.0 and 5.0
928  for (size_t i = 0 ; i < dim ; i++)
929  {
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;}
934  else
935  {arxvalid[j](i) = arx[j](i);}
936  }
937 
938  f_obj.get(j).f = hybrid_composition<dim>(arxvalid[j]);
939  f_obj.get(j).id = j;
940  fe++;
941 
942  // Get the best ever
943  if (f_obj.get(j).f < best_sample)
944  {
945  best_sample = f_obj.get(j).f;
946 
947  // Copy the new mean as position of the particle
948  for (size_t i = 0 ; i < dim ; i++)
949  {best_sample_sol.get(i) = arxvalid[j](i);}
950  }
951  }
952 
953  // Add penalities for out of bound points
954  cmaes_handlebounds(f_obj,vd.getProp<sigma>(p),
955  vd.getProp<validfit>(p),arxvalid,
956  arx,vd.getProp<Cov_m>(p),
957  xmean,vd.getProp<xold>(p),vd.getProp<weight>(p),
958  vd.getProp<fithist>(p),vd.getProp<iniphase>(p),
959  vd.getProp<validfit>(p),mu_eff,
960  step,vd.getProp<last_restart>(p));
961 
962  f_obj.sort();
963 
964  for (size_t j = 0 ; j < lambda ; j++)
965  {vd.getProp<ord>(p)[j] = f_obj.get(j).id;}
966 
967  vd.getProp<xold>(p) = xmean;
968 
969  // Calculate weighted mean
970 
971  xmean.setZero();
972  mean_z.setZero();
973  for (size_t j = 0 ; j < mu ; j++)
974  {
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]];
977  }
978 
979  vd.getProp<xmean_st>(p) = xmean;
980  vd.getProp<meanz_st>(p) = mean_z;
981 
982  ++it;
983  }
984 
985  // Find the best point across processors
986  broadcast_best_solution(vd,best_sol,best,best_sample,best_sample_sol);
987 
988  // bool calculate B and D
989  bool calc_bd = counteval - eigeneval > lambda/(ccov)/dim/10;
990  if (calc_bd == true)
991  {eigeneval = counteval;}
992 
993  auto it2 = vd.getDomainIterator();
994  while (it2.isNext())
995  {
996  auto p = it2.get();
997 
998  if (vd.getProp<stop>(p) == true)
999  {++it2;continue;}
1000 
1001  xmean = vd.getProp<xmean_st>(p);
1002  mean_z = vd.getProp<meanz_st>(p);
1003 
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;
1005 
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);
1007 
1008  vd.getProp<path_c>(p) = (1-cc)*vd.getProp<path_c>(p) + hsig * sqrt(cc*(2-cc)*mu_eff)*(vd.getProp<B>(p)*vd.getProp<D>(p)*mean_z);
1009 
1010  if (step % N_pso == 0)
1011  {
1012  EMatrixXd C_pso(dim,dim);
1013  updatePso(best_sol,vd.getProp<sigma>(p),xmean,vd.getProp<xold>(p),vd.getProp<B>(p),vd.getProp<D>(p),C_pso);
1014 
1015  // Adapt covariance matrix C
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());
1018 
1019  for (size_t i = 0 ; i < mu ; i++)
1020  {vd.getProp<Cov_m>(p) += ccov*(1.0-1.0/mu_eff)*(vd.getProp<B>(p)*vd.getProp<D>(p)*vd.getProp<Zeta>(p)[vd.getProp<ord>(p)[i]])*weight_sample(i)*
1021  (vd.getProp<B>(p)*vd.getProp<D>(p)*vd.getProp<Zeta>(p)[vd.getProp<ord>(p)[i]]).transpose();
1022  }
1023 
1024  vd.getProp<Cov_m>(p) = psoWeight*vd.getProp<Cov_m>(p) + (1.0 - psoWeight)*C_pso;
1025  }
1026  else
1027  {
1028  // Adapt covariance matrix C
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());
1031 
1032  for (size_t i = 0 ; i < mu ; i++)
1033  {vd.getProp<Cov_m>(p) += ccov*(1.0-1.0/mu_eff)*(vd.getProp<B>(p)*vd.getProp<D>(p)*vd.getProp<Zeta>(p)[vd.getProp<ord>(p)[i]])*weight_sample(i)*
1034  (vd.getProp<B>(p)*vd.getProp<D>(p)*vd.getProp<Zeta>(p)[vd.getProp<ord>(p)[i]]).transpose();
1035  }
1036  }
1037 
1038  // Numeric error
1039 
1040  double smaller = std::numeric_limits<double>::max();
1041  for (size_t i = 0 ; i < dim ; i++)
1042  {
1043  if (vd.getProp<sigma>(p)*sqrt(vd.getProp<D>(p).diagonal()[i]) > 5.0)
1044  {
1045  if (smaller > 5.0/sqrt(vd.getProp<D>(p).diagonal()[i]))
1046  {smaller = 5.0/sqrt(vd.getProp<D>(p).diagonal()[i]);}
1047  }
1048  }
1049  if (smaller != std::numeric_limits<double>::max())
1050  {vd.getProp<sigma>(p) = smaller;}
1051 
1052  //Adapt step-size sigma
1053  vd.getProp<sigma>(p) = vd.getProp<sigma>(p)*exp((cs/d_amps)*(vd.getProp<path_s>(p).norm()/chiN - 1));
1054 
1055  // Update B and D from C
1056 
1057  if (calc_bd)
1058  {
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();
1062 
1063  // Eigen decomposition
1064  Eigen::SelfAdjointEigenSolver<EMatrixXd> eig_solver;
1065 
1066  eig_solver.compute(vd.getProp<Cov_m>(p));
1067 
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();
1071 
1072  // Make first component always positive
1073  for (size_t i = 0 ; i < dim ; i++)
1074  {
1075  if (vd.getProp<B>(p)(0,i) < 0)
1076  {vd.getProp<B>(p).col(i) = - vd.getProp<B>(p).col(i);}
1077  }
1078 
1079  EMatrixXd tmp = vd.getProp<B>(p).transpose();
1080  }
1081 
1082  // Copy the new mean as position of the particle
1083  for (size_t i = 0 ; i < dim ; i++)
1084  {vd.getPos(p)[i] = xmean(i);}
1085 
1086  vd.getProp<sigma>(p) = adjust_sigma(vd.getProp<sigma>(p),vd.getProp<Cov_m>(p));
1087 
1088  // Stop conditions
1089  bool stop_tol = true;
1090  bool stop_tolX = true;
1091  for (size_t i = 0 ; i < dim ; i++)
1092  {
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;
1095  }
1096 
1097  vd.getProp<stop>(p) = stop_tol | stop_tolX;
1098 
1099  // Escape flat fitness, or better terminate?
1100  if (f_obj.get(0).f == f_obj.get(std::ceil(0.7*lambda)).f )
1101  {
1102  vd.getProp<sigma>(p) = vd.getProp<sigma>(p)*exp(0.2+cs/d_amps);
1103  std::cout << "warning: flat fitness, consider reformulating the objective";
1104 
1105  // Stop it
1106  vd.getProp<stop>(p) = true;
1107  }
1108 
1109  if (vd.getProp<stop>(p) == true)
1110  {std::cout << "Stopped" << std::endl;}
1111 
1112  if (restart_cma && vd.getProp<stop>(p) == true)
1113  {
1114  std::cout << "------- Restart #" << std::endl;
1115 
1116  std::cout << "---------------------------------" << std::endl;
1117  std::cout << "Best: " << best << " " << fun_eval << std::endl;
1118  std::cout << "---------------------------------" << std::endl;
1119 
1120  vd.getProp<last_restart>(p) = step;
1121  vd.getProp<xold>(p).setZero();
1122 
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();
1127  vd.getProp<Cov_m>(p) = vd.getProp<B>(p)*vd.getProp<D>(p)*vd.getProp<D>(p)*vd.getProp<B>(p);
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);
1132  vd.getProp<stop>(p) = false;
1133  vd.getProp<iniphase>(p) = true;
1134  vd.getProp<last_restart>(p) = 0;
1135  vd.getProp<sigma>(p) = 2.0;
1136 
1137  // a different point in space
1138  for (size_t i = 0 ; i < dim ; i++)
1139  {
1140  // we define x, assign a random position between 0.0 and 1.0
1141  vd.getPos(p)[i] = 10.0*(double)rand() / RAND_MAX - 5.0;
1142  }
1143 
1144  // Initialize the bound history
1145 
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;
1150  }
1151 
1152  ++it2;
1153  }
1154 
1155  auto & v_cl = create_vcluster();
1156  v_cl.sum(fe);
1157  v_cl.execute();
1158 
1159  fun_eval += fe;
1160 }
1161 
1163 
1177 int main(int argc, char* argv[])
1178 {
1179  // initialize the library
1180  openfpm_init(&argc,&argv);
1181 
1182  auto & v_cl = create_vcluster();
1183 
1184  // Here we define our domain a 2D box with internals from 0 to 1.0 for x and y
1185  Box<dim,double> domain;
1186 
1187  for (size_t i = 0 ; i < dim ; i++)
1188  {
1189  domain.setLow(i,0.0);
1190  domain.setHigh(i,1.0);
1191  }
1192 
1193  // Here we define the boundary conditions of our problem
1194  size_t bc[dim];
1195  for (size_t i = 0 ; i < dim ; i++)
1196  {bc[i] = NON_PERIODIC;};
1197 
1198  prepare_f15<dim>();
1199 
1200  // extended boundary around the domain, and the processor domain
1201  Ghost<dim,double> g(0.0);
1202 
1203  particle_type vd(16,domain,bc,g);
1204 
1205  // Initialize constants
1206 
1207  stop_fitness = 1e-10;
1208  size_t stopeval = 1e3*dim*dim;
1209 
1210  // Strategy parameter setting: Selection
1211  init_weight();
1212 
1213  // Strategy parameter setting: Adaptation
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;
1219 
1220  chiN = sqrt(dim)*(1.0-1.0/(4.0*dim)+1.0/(21.0*dim*dim));
1221 
1223 
1224 
1225  // initialize the srand
1226  int seed = 24756*v_cl.rank()*v_cl.rank() + time(NULL);
1227  srand(seed);
1228 
1229  auto it = vd.getDomainIterator();
1230 
1231  while (it.isNext())
1232  {
1233  auto p = it.get();
1234 
1235  for (size_t i = 0 ; i < dim ; i++)
1236  {
1237  // we define x, assign a random position between 0.0 and 1.0
1238  vd.getPos(p)[i] = 10.0*(double)rand() / RAND_MAX - 5.0;
1239  }
1240 
1241  vd.getProp<sigma>(p) = 2.0;
1242 
1243  // Initialize the covariant Matrix,B and D to identity
1244 
1245  vd.getProp<D>(p).resize(dim);
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();
1250  vd.getProp<Cov_m>(p) = vd.getProp<B>(p)*vd.getProp<D>(p)*vd.getProp<D>(p)*vd.getProp<B>(p);
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);
1255  vd.getProp<stop>(p) = false;
1256  vd.getProp<iniphase>(p) = true;
1257  vd.getProp<last_restart>(p) = 0;
1258 
1259  // Initialize the bound history
1260 
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;
1265 
1266  // next particle
1267  ++it;
1268  }
1269 
1270  if (v_cl.rank() == 0)
1271  {std::cout << "Starting PS-CMA-ES" << std::endl;}
1272 
1273  double best = 0.0;
1274  int best_i = 0;
1275 
1276  best = std::numeric_limits<double>::max();
1277  openfpm::vector<double> best_sol(dim);
1278  // now do several iteration
1279 
1280  int stop_cond = 0;
1281  size_t fun_eval = 0;
1282  int i = 0;
1283  while (fun_eval < max_fun_eval && best > 120.000001)
1284  {
1285  // sample offspring
1286  cma_step(vd,i+1,best,best_i,best_sol,fun_eval);
1287 
1288  i++;
1289  }
1290 
1291  if (v_cl.rank() == 0)
1292  {
1293  std::cout << "Best solution: " << best << " with " << fun_eval << std::endl;
1294  std::cout << "at: " << std::endl;
1295 
1296  for (size_t i = 0 ; i < best_sol.size() ; i++)
1297  {
1298  std::cout << best_sol.get(i) << " ";
1299  }
1300  }
1301 
1302  openfpm_finalize();
1303 
1305 
1315 }
Derivative second order on h (spacing)
Definition: Derivative.hpp:28
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
Definition: Box.hpp:467
size_t size()
Stub size.
Definition: map_vector.hpp:70
Definition: Ghost.hpp:39
void setLow(int i, T val)
set the low interval of the box
Definition: Box.hpp:456
This class represent an N-dimensional box.
Definition: Box.hpp:56
vector_dist_iterator getDomainIterator() const
Get an iterator that traverse the particles in the domain.
Distributed vector.
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...
Definition: aggregate.hpp:81
It model an expression expr1 + ... exprn.
Definition: sum.hpp:92
Implementation of 1-D std::vector like structure.
Definition: map_vector.hpp:61