Salome HOME
updated copyright message
[modules/yacs.git] / src / engine / Plugin / point.cxx
1 // Copyright (C) 2006-2023  CEA/DEN, EDF R&D
2 //
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
7 //
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 // Lesser General Public License for more details.
12 //
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
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