src/engine/Plugin/simplex.cxx

   1 // --- C++ ---
   2 // --- coding: latin_1 ---
   3 //
   4 //    File
   5 //      creation : 2007-03-30.13.44.14
   6 //      revision : $Id$
   7 //
   8 //    Copyright © 2007 Commissariat à l'Energie Atomique
   9 //      par Gilles ARNAUD (DM2S/SFME/LETR)
  10 //        C.E. Saclay; Bat 454; 91191 GIF/YVETTE CEDEX; France
  11 //        Tel: 01 69 08 38 86; Fax : 33 1 69 08 85 68
  12 //        Gilles.Arnaud@cea.fr
  13 //
  14 //    Object
  15 //      version distribue du simplex de nelder mead
  16 //
  17 //___________________________________________________________________
  18
  19
  20 #include "simplex.hxx"
  21
  22 #include <iostream>
  23
  24 #include "aleas.hxx"
  25
  26 Simplex::Simplex(long sz, long nbgen, Maestro &maest)
  27 {
  28     size = sz;
  29     nbin = nbgen;
  30     calc = &maest;
  31     maxe = 40000;
  32     nbeval = 0;
  33 }
  34
  35 Simplex::~Simplex(void)
  36 {
  37     int     i;
  38
  39     if (work.size() == size) {
  40         for (i=0; i<size; i++)
  41             delete work[i];
  42         for (i=0; i<nbin; i++)
  43             delete simplx[i];
  44     }
  45 }
  46
  47 void Simplex::setStop(long max)
  48 {
  49     maxe = max;
  50 }
  51
  52 void Simplex::start(void)
  53 {
  54     long    i;
  55
  56     work.resize(size);
  57     for (i=0; i<size; i++) {
  58         work[i] = new Point(nbin);
  59         calc->put(i, *(work[i]->next()));
  60     }
  61     nbeval = maxe;
  62 }
  63
  64
  65 int Simplex::next(void)
  66 {
  67     long        id;
  68     Solution    *pt, *pnt;
  69     std::vector<double>     *next, *res;
  70
  71     res = calc->get(&id);
  72     nbeval--;
  73     pt = work[id]->inform(*res);
  74     if (pt) {
  75         pnt = add(pt);
  76         if (pnt)
  77             work[id]->mute(*pnt, *barycentre(), *minimum());
  78         else
  79             work[id]->reinit();
  80     }
  81     next = work[id]->next();
  82     if (nbeval > size)
  83         if (next)
  84             calc->put(id, *next);
  85         //else
  86         //    nbeval = size - 1;
  87     return (nbeval > 0);
  88 }
  89
  90 void Simplex::finish(void)
  91 {
  92     std::cout << maxe - nbeval << std::endl;
  93     return;
  94 }
  95
  96 Solution *Simplex::solution(void)
  97 {
  98     return simplx[0];
  99 }
 100
 101 Solution *Simplex::add(Solution *sol)
 102 {
 103     int     i;
 104     Solution    *ret, *swp;
 105
 106     if (simplx.size() < nbin) {
 107         ret = (Solution *) NULL;
 108         i = simplx.size();
 109         simplx.push_back(sol);
 110         while (i && simplx[i]->obj[0] < simplx[i-1]->obj[0]) {
 111             // TODO swap interne
 112             swp = simplx[i];
 113             simplx[i] = simplx[i-1];
 114             simplx[i-1] = swp;
 115             i--;
 116         }
 117     }
 118     else if (sol->obj[0] > simplx[nbin-1]->obj[0])
 119         ret = sol;
 120     else {
 121         i = nbin -1;
 122         ret = simplx[i];
 123         simplx[i] = sol;
 124         while (i && simplx[i]->obj[0] < simplx[i-1]->obj[0]) {
 125             swp = simplx[i];
 126             simplx[i] = simplx[i-1];
 127             simplx[i-1] = swp;
 128             i--;
 129         }
 130     }
 131     return ret;
 132 }
 133
 134 std::vector<double> *Simplex::minimum(void)
 135 {
 136     return new std::vector<double>(*simplx[0]->param);
 137 }
 138
 139 std::vector<double> *Simplex::barycentre(void)
 140 {
 141     int     i,j;
 142     std::vector<double>     *ret;
 143
 144     ret = new std::vector<double>(nbin);
 145     for (i=0; i<nbin; i++) {
 146         (*ret)[i] = 0.0;
 147         for (j=0; j<nbin; j++)
 148             (*ret)[i] += (*simplx[j]->param)[i];
 149         (*ret)[i] /= nbin;
 150     }
 151     return ret;
 152 }