src/engine/Plugin/point.cxx

   1 // --- C++ ---
   2 // --- coding: latin_1 ---
   3 //
   4 //    File
   5 //      creation : 2007-03-30.16.58.42
   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 //      particules du simplex
  16 //
  17 //___________________________________________________________________
  18
  19
  20 #include "point.hxx"
  21
  22 #include <iostream>
  23
  24 #define START   0
  25 #define FIRST   1
  26 #define GOOD    2
  27 #define SOGOOD  3
  28 #define TOOGOOD 4
  29 #define BAD     5
  30 #define TOOBAD  6
  31 #define SOBAD   7
  32 #define BADEST  8
  33
  34 #define SEUIL   1.0e-10
  35
  36 Point::Point(long sz)
  37 {
  38     size = sz;
  39     etat = START;
  40     epsil = SEUIL;
  41     rnd = new Cube(size);
  42     start = (Solution *) NULL;
  43     courant = (std::vector<double> *) NULL;
  44     baryc = (std::vector<double> *) NULL;
  45     minim = (std::vector<double> *) NULL;
  46 }
  47
  48 Point::~Point(void)
  49 {
  50     delete rnd;
  51     if (courant)
  52         delete courant;
  53     if (start) {
  54         delete start;
  55         delete baryc;
  56         delete minim;
  57     }
  58 }
  59
  60 Solution *Point::inform(std::vector<double> &obj)
  61 {
  62     Solution        *res, *tmp;
  63
  64     res = (Solution *) NULL;
  65     switch (etat) {
  66         case START :
  67 #ifdef SHOW
  68             std::cout << "@" ;
  69 #endif
  70             res = new Solution(*courant, obj);
  71             break;
  72         case FIRST :
  73             if (obj[0] > (*start->obj)[0]) {
  74 #ifdef SHOW
  75                 std::cout << "-" ;
  76 #endif
  77                 etat = BAD ;
  78                 delete courant;
  79                 delete &obj;
  80             }
  81             else {
  82 #ifdef SHOW
  83                 std::cout << "+" ;
  84 #endif
  85                 etat = GOOD ;
  86                 tmp = start;
  87                 start = new Solution(*courant, obj);
  88                 delete tmp;
  89             }
  90             break;
  91         case GOOD :
  92             if (obj[0] > (*start->obj)[0]) {
  93 #ifdef SHOW
  94                 std::cout << "P" ;
  95 #endif
  96                 etat = SOGOOD;
  97                 delete courant;
  98                 delete &obj;
  99                 res = start;
 100                 delete baryc;
 101                 delete minim;
 102                 start = (Solution *) NULL;
 103             }
 104             else {
 105 #ifdef SHOW
 106                 std::cout << "p" ;
 107 #endif
 108                 etat = TOOGOOD;
 109                 res = new Solution(*courant, obj);
 110             }
 111             break;
 112         case BAD :
 113             if (obj[0] > (*start->obj)[0]) {
 114 #ifdef SHOW
 115                 std::cout << "M" ;
 116 #endif
 117                 etat = TOOBAD;
 118                 delete courant;
 119                 delete &obj;
 120             }
 121             else {
 122 #ifdef SHOW
 123                 std::cout << "m" ;
 124 #endif
 125                 etat = SOBAD;
 126                 res = new Solution(*courant, obj);
 127             }
 128             break;
 129         case TOOBAD:
 130 #ifdef SHOW
 131             std::cout << "X" ;
 132 #endif
 133             etat = BADEST;
 134             res = new Solution(*courant, obj);
 135             break;
 136         default :
 137             std::cout << "pbl inform" << std::endl ;
 138             break;
 139     }
 140     return res;
 141 }
 142
 143 void Point::mute(Solution &pt, std::vector<double> &bary, std::vector<double> &minm)
 144 {
 145     std::vector<double>     *tmp;
 146
 147     if (start) {
 148         delete baryc;
 149         delete minim;
 150         delete start;
 151     }
 152     start = &pt;
 153     baryc = &bary;
 154     minim = &minm;
 155     etat = FIRST;
 156 }
 157
 158 void Point::reinit(void)
 159 {
 160     etat = START;
 161 }
 162
 163 std::vector<double> *Point::next(void)
 164 {
 165     double      d, dd;
 166     int     i;
 167
 168     switch (etat) {
 169         case START :
 170             courant = rnd->gen();
 171             break;
 172         case FIRST :
 173             courant = symetrique(*start->param, *baryc);
 174             break;
 175         case GOOD :
 176             courant = symetrique(*baryc, *start->param);
 177             break;
 178         case BAD :
 179             courant = milieu(*baryc, *start->param);
 180             break;
 181         case TOOBAD :
 182             courant = milieu(*minim, *start->param);
 183             break;
 184         default:
 185             // raise
 186             std::cout << "pbl next" << std::endl ;
 187     }
 188     if (baryc) {
 189         d=0.0;
 190         for (i=0; i<size; i++) {
 191             dd = (*courant)[0] - (*baryc)[0];
 192             d += dd*dd;
 193         }
 194     }
 195     else
 196         d=1.0;
 197     if (d < epsil) {
 198         //delete courant;
 199         return (std::vector<double>  *) NULL;
 200     }
 201     else
 202         return courant;
 203 }
 204
 205 std::vector<double> *Point::symetrique(std::vector<double> &pt, std::vector<double> &centr)
 206 {
 207     long        i;
 208     double      coef, tmp;
 209     std::vector<double>     *res;
 210
 211     // bounded coef
 212     coef = 1.0;
 213     for (i=0; i<size; i++) {
 214         tmp = (centr[i]-pt[i] > 0.0) ?
 215                 (1.0 - centr[i]) / (centr[i] - pt[i]) :
 216                 centr[i] / (pt[i] - centr[i]) ;
 217         coef = (coef < tmp) ? coef : tmp ;
 218     }
 219     //
 220     res = new std::vector<double>(size);
 221     for (i=0; i<size; i++)
 222         (*res)[i] = centr[i] + coef * (centr[i] - pt[i]);
 223     return res;
 224 }
 225
 226 std::vector<double> *Point::milieu(std::vector<double> &un, std::vector<double> &deux)
 227 {
 228     long        i;
 229     std::vector<double>     *res;
 230
 231     res = new std::vector<double>(size);
 232     for (i=0; i<size; i++)
 233         (*res)[i] = (un[i] + deux[i])/2.0;
 234     return res;
 235 }
 236
 237
 238