Commit | Line | Data |
---|---|---|
b311480e | 1 | // Copyright (c) 1997-1999 Matra Datavision |
973c2be1 | 2 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
b311480e | 3 | // |
973c2be1 | 4 | // This file is part of Open CASCADE Technology software library. |
b311480e | 5 | // |
d5f74e42 | 6 | // This library is free software; you can redistribute it and/or modify it under |
7 | // the terms of the GNU Lesser General Public License version 2.1 as published | |
973c2be1 | 8 | // by the Free Software Foundation, with special exception defined in the file |
9 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT | |
10 | // distribution for complete text of the license and disclaimer of any warranty. | |
b311480e | 11 | // |
973c2be1 | 12 | // Alternatively, this file may be used under the terms of Open CASCADE |
13 | // commercial license or contractual agreement. | |
b311480e | 14 | |
7fd59977 | 15 | // pmn 15/05/97 pas de Gauss avec des pivot trop petit. SVD fait mieux |
16 | // l'affaire + limitation de la longeur du pas + qq comentaire issus d'EUCLID3 | |
17 | // pmn 10/06/97 refonte totale du traitement des bords + ajustement des init | |
18 | // et des tolerances pour brent... | |
19 | ||
0797d9d3 | 20 | //#ifndef OCCT_DEBUG |
7fd59977 | 21 | #define No_Standard_RangeError |
22 | #define No_Standard_OutOfRange | |
23 | #define No_Standard_DimensionError | |
7fd59977 | 24 | |
42cf5bc1 | 25 | //#endif |
7fd59977 | 26 | //math_FunctionSetRoot.cxx |
27 | ||
42cf5bc1 | 28 | #include <math_BrentMinimum.hxx> |
29 | #include <math_Function.hxx> | |
30 | #include <math_FunctionSetRoot.hxx> | |
31 | #include <math_FunctionSetWithDerivatives.hxx> | |
7fd59977 | 32 | #include <math_Gauss.hxx> |
7fd59977 | 33 | #include <math_GaussLeastSquare.hxx> |
34 | #include <math_IntegerVector.hxx> | |
42cf5bc1 | 35 | #include <math_Matrix.hxx> |
36 | #include <math_SVD.hxx> | |
7fd59977 | 37 | #include <Precision.hxx> |
42cf5bc1 | 38 | #include <Standard_DimensionError.hxx> |
39 | #include <StdFail_NotDone.hxx> | |
7fd59977 | 40 | |
41 | //=========================================================================== | |
42 | // - A partir d une solution de depart, recherche d une direction.( Newton la | |
43 | // plupart du temps, gradient si Newton echoue. | |
7fd59977 | 44 | // - Recadrage au niveau des bornes avec recalcul de la direction si une |
45 | // inconnue a une valeur imposee. | |
7fd59977 | 46 | // -Si On ne sort pas des bornes |
47 | // Tant que (On ne progresse pas assez ou on ne change pas de direction) | |
48 | // . Si (Progression encore possible) | |
49 | // Si (On ne sort pas des bornes) | |
50 | // On essaye de progresser dans cette meme direction. | |
51 | // Sinon | |
52 | // On diminue le pas d'avancement ou on change de direction. | |
53 | // Sinon | |
54 | // Si on depasse le minimum | |
55 | // On fait une interpolation parabolique. | |
7fd59977 | 56 | // - Si on a progresse sur F |
57 | // On fait les tests d'arret | |
58 | // Sinon | |
59 | // On change de direction | |
60 | //============================================================================ | |
3e42bd70 J |
61 | #define FSR_DEBUG(arg) |
62 | // Uncomment the following code to have debug output to cout | |
89d8607f | 63 | //========================================================== |
64 | //static Standard_Boolean mydebug = Standard_True; | |
65 | //#undef FSR_DEBUG | |
66 | //#define FSR_DEBUG(arg) {if (mydebug) { cout << arg << endl; }} | |
67 | //=========================================================== | |
fbadd2cc | 68 | |
7fd59977 | 69 | class MyDirFunction : public math_Function |
70 | { | |
71 | ||
89d8607f | 72 | math_Vector *P0; |
73 | math_Vector *Dir; | |
74 | math_Vector *P; | |
75 | math_Vector *FV; | |
76 | math_FunctionSetWithDerivatives *F; | |
7fd59977 | 77 | |
78 | public : | |
79 | ||
89d8607f | 80 | MyDirFunction(math_Vector& V1, |
81 | math_Vector& V2, | |
82 | math_Vector& V3, | |
83 | math_Vector& V4, | |
84 | math_FunctionSetWithDerivatives& f) ; | |
85 | ||
86 | void Initialize(const math_Vector& p0, const math_Vector& dir) const; | |
87 | //For hp : | |
88 | Standard_Boolean Value(const math_Vector& Sol, math_Vector& FF, | |
89 | math_Matrix& DF, math_Vector& GH, | |
90 | Standard_Real& F2, Standard_Real& Gnr1); | |
91 | // Standard_Boolean MyDirFunction::Value(const math_Vector& Sol, math_Vector& FF, | |
92 | // math_Matrix& DF, math_Vector& GH, | |
93 | // Standard_Real& F2, Standard_Real& Gnr1); | |
94 | Standard_Boolean Value(const Standard_Real x, Standard_Real& fval) ; | |
95 | ||
7fd59977 | 96 | }; |
97 | ||
98 | MyDirFunction::MyDirFunction(math_Vector& V1, | |
89d8607f | 99 | math_Vector& V2, |
100 | math_Vector& V3, | |
101 | math_Vector& V4, | |
102 | math_FunctionSetWithDerivatives& f) { | |
103 | ||
104 | P0 = &V1; | |
105 | Dir = &V2; | |
106 | P = &V3; | |
107 | FV = &V4; | |
108 | F = &f; | |
7fd59977 | 109 | } |
110 | ||
111 | void MyDirFunction::Initialize(const math_Vector& p0, | |
89d8607f | 112 | const math_Vector& dir) const |
7fd59977 | 113 | { |
114 | *P0 = p0; | |
115 | *Dir = dir; | |
116 | } | |
117 | ||
118 | ||
119 | Standard_Boolean MyDirFunction::Value(const Standard_Real x, | |
89d8607f | 120 | Standard_Real& fval) |
7fd59977 | 121 | { |
122 | Standard_Real p; | |
123 | for(Standard_Integer i = P->Lower(); i <= P->Upper(); i++) { | |
124 | p = Dir->Value(i); | |
125 | P->Value(i) = p * x + P0->Value(i); | |
126 | } | |
e93e4230 | 127 | if( F->Value(*P, *FV) ) |
128 | { | |
7fd59977 | 129 | |
e93e4230 | 130 | Standard_Real aVal = 0.0; |
7fd59977 | 131 | |
e93e4230 | 132 | for(Standard_Integer i = FV->Lower(); i <= FV->Upper(); i++) |
133 | { | |
7fd59977 | 134 | aVal = FV->Value(i); |
e93e4230 | 135 | if(aVal <= -1.e+100) // Precision::HalfInfinite() later |
136 | return Standard_False; | |
7fd59977 | 137 | else if(aVal >= 1.e+100) // Precision::HalfInfinite() later |
89d8607f | 138 | return Standard_False; |
7fd59977 | 139 | } |
140 | ||
141 | fval = 0.5 * (FV->Norm2()); | |
142 | return Standard_True; | |
143 | } | |
144 | return Standard_False; | |
145 | } | |
146 | ||
147 | Standard_Boolean MyDirFunction::Value(const math_Vector& Sol, | |
89d8607f | 148 | math_Vector& FF, |
149 | math_Matrix& DF, | |
150 | math_Vector& GH, | |
151 | Standard_Real& F2, | |
152 | Standard_Real& Gnr1) | |
7fd59977 | 153 | { |
154 | if( F->Values(Sol, FF, DF) ) { | |
155 | ||
156 | Standard_Real aVal = 0.; | |
157 | ||
158 | for(Standard_Integer i = FF.Lower(); i <= FF.Upper(); i++) { | |
89d8607f | 159 | // modified by NIZHNY-MKK Mon Oct 3 17:56:50 2005.BEGIN |
7fd59977 | 160 | aVal = FF.Value(i); |
161 | if(aVal < 0.) { | |
89d8607f | 162 | if(aVal <= -1.e+100) // Precision::HalfInfinite() later |
163 | // if(Precision::IsInfinite(Abs(FF.Value(i)))) { | |
164 | // F2 = Precision::Infinite(); | |
165 | // Gnr1 = Precision::Infinite(); | |
166 | return Standard_False; | |
7fd59977 | 167 | } |
168 | else if(aVal >= 1.e+100) // Precision::HalfInfinite() later | |
89d8607f | 169 | return Standard_False; |
170 | // modified by NIZHNY-MKK Mon Oct 3 17:57:05 2005.END | |
7fd59977 | 171 | } |
172 | ||
173 | ||
174 | F2 = 0.5 * (FF.Norm2()); | |
175 | GH.TMultiply(DF, FF); | |
12f139fd | 176 | for(Standard_Integer i = GH.Lower(); i <= GH.Upper(); i++) |
177 | { | |
178 | if(Precision::IsInfinite((GH.Value(i)))) | |
179 | { | |
180 | return Standard_False; | |
181 | } | |
182 | } | |
7fd59977 | 183 | Gnr1 = GH.Norm2(); |
184 | return Standard_True; | |
185 | } | |
186 | return Standard_False; | |
187 | } | |
188 | ||
189 | ||
190 | //-------------------------------------------------------------- | |
191 | static Standard_Boolean MinimizeDirection(const math_Vector& P0, | |
89d8607f | 192 | const math_Vector& P1, |
193 | const math_Vector& P2, | |
194 | const Standard_Real F1, | |
195 | math_Vector& Delta, | |
196 | const math_Vector& Tol, | |
197 | MyDirFunction& F) | |
198 | // Purpose : minimisation a partir de 3 points | |
199 | //------------------------------------------------------- | |
7fd59977 | 200 | { |
201 | // (1) Evaluation d'un tolerance parametrique 1D | |
202 | Standard_Real tol1d = 2.1 , invnorme, tsol; | |
203 | Standard_Real Eps = 1.e-16; | |
204 | Standard_Real ax, bx, cx; | |
205 | ||
206 | for (Standard_Integer ii =1; ii<=Tol.Length(); ii++) { | |
207 | invnorme = Abs(Delta(ii)); | |
208 | if (invnorme>Eps) tol1d = Min(tol1d, Tol(ii) / invnorme); | |
209 | } | |
210 | if (tol1d > 1.9) return Standard_False; //Pas la peine de se fatiguer | |
211 | tol1d /= 3; | |
212 | ||
89d8607f | 213 | //JR/Hp : |
7fd59977 | 214 | math_Vector PP0 = P0 ; |
215 | math_Vector PP1 = P1 ; | |
216 | Delta = PP1 - PP0; | |
89d8607f | 217 | // Delta = P1 - P0; |
7fd59977 | 218 | invnorme = Delta.Norm(); |
219 | if (invnorme <= Eps) return Standard_False; | |
220 | invnorme = ((Standard_Real) 1) / invnorme; | |
221 | ||
222 | F.Initialize(P1, Delta); | |
223 | ||
224 | // (2) On minimise | |
3e42bd70 | 225 | FSR_DEBUG (" minimisation dans la direction") |
89d8607f | 226 | ax = -1; bx = 0; |
7fd59977 | 227 | cx = (P2-P1).Norm()*invnorme; |
859a47c3 | 228 | if (cx < 1.e-2) |
229 | return Standard_False; | |
230 | ||
231 | math_BrentMinimum Sol(tol1d, 100, tol1d); | |
232 | Sol.Perform(F, ax, bx, cx); | |
233 | ||
7fd59977 | 234 | if(Sol.IsDone()) { |
235 | tsol = Sol.Location(); | |
236 | if (Sol.Minimum() < F1) { | |
237 | Delta.Multiply(tsol); | |
238 | return Standard_True; | |
239 | } | |
240 | } | |
241 | return Standard_False; | |
242 | } | |
243 | ||
244 | //---------------------------------------------------------------------- | |
245 | static Standard_Boolean MinimizeDirection(const math_Vector& P, | |
89d8607f | 246 | math_Vector& Dir, |
247 | const Standard_Real& PValue, | |
248 | const Standard_Real& PDirValue, | |
249 | const math_Vector& Gradient, | |
250 | const math_Vector& DGradient, | |
251 | const math_Vector& Tol, | |
252 | MyDirFunction& F) | |
253 | // Purpose: minimisation a partir de 2 points et une derives | |
254 | //---------------------------------------------------------------------- | |
7fd59977 | 255 | |
256 | { | |
f4dee9bb | 257 | if(Precision::IsInfinite(PValue) || Precision::IsInfinite(PDirValue)) |
258 | { | |
259 | return Standard_False; | |
260 | } | |
7fd59977 | 261 | // (0) Evaluation d'un tolerance parametrique 1D |
262 | Standard_Boolean good = Standard_False; | |
263 | Standard_Real Eps = 1.e-20; | |
264 | Standard_Real tol1d = 1.1, Result = PValue, absdir; | |
265 | ||
266 | for (Standard_Integer ii =1; ii<=Tol.Length(); ii++) { | |
267 | absdir = Abs(Dir(ii)); | |
268 | if (absdir >Eps) tol1d = Min(tol1d, Tol(ii) / absdir); | |
269 | } | |
270 | if (tol1d > 0.9) return Standard_False; | |
89d8607f | 271 | |
7fd59977 | 272 | // (1) On realise une premiere interpolation quadratique |
273 | Standard_Real ax, bx, cx, df1, df2, Delta, tsol, fsol, tsolbis; | |
89d8607f | 274 | FSR_DEBUG(" essai d interpolation"); |
fbadd2cc | 275 | |
7fd59977 | 276 | df1 = Gradient*Dir; |
277 | df2 = DGradient*Dir; | |
278 | ||
279 | if (df1<-Eps && df2>Eps) { // cuvette | |
280 | tsol = - df1 / (df2 - df1); | |
281 | } | |
282 | else { | |
283 | cx = PValue; | |
284 | bx = df1; | |
285 | ax = PDirValue - (bx+cx); | |
286 | ||
287 | if (Abs(ax) <= Eps) { // cas lineaire | |
288 | if ((Abs(bx) >= Eps)) tsol = - cx/bx; | |
289 | else tsol = 0; | |
290 | } | |
291 | else { // cas quadratique | |
292 | Delta = bx*bx - 4*ax*cx; | |
293 | if (Delta > 1.e-9) { | |
89d8607f | 294 | // il y a des racines, on prend la plus proche de 0 |
295 | Delta = Sqrt(Delta); | |
296 | tsol = -(bx + Delta); | |
297 | tsolbis = (Delta - bx); | |
298 | if (Abs(tsolbis) < Abs(tsol)) tsol = tsolbis; | |
299 | tsol /= 2*ax; | |
7fd59977 | 300 | } |
301 | else { | |
89d8607f | 302 | // pas ou peu de racine : on "extremise" |
303 | tsol = -(0.5*bx)/ax; | |
7fd59977 | 304 | } |
305 | } | |
306 | } | |
307 | ||
308 | if (Abs(tsol) >= 1) return Standard_False; //resultat sans interet | |
309 | ||
310 | F.Initialize(P, Dir); | |
311 | F.Value(tsol, fsol); | |
312 | ||
313 | if (fsol<PValue) { | |
314 | good = Standard_True; | |
315 | Result = fsol; | |
89d8607f | 316 | FSR_DEBUG("t= "<<tsol<<" F = " << fsol << " OldF = "<<PValue); |
7fd59977 | 317 | } |
318 | ||
319 | // (2) Si l'on a pas assez progresser on realise une recherche | |
320 | // en bonne et due forme, a partir des inits precedents | |
e93e4230 | 321 | if ((fsol > 0.2*PValue) && (tol1d < 0.5)) |
322 | { | |
89d8607f | 323 | |
7fd59977 | 324 | if (tsol <0) { |
325 | ax = tsol; bx = 0.0; cx = 1.0; | |
326 | } | |
327 | else { | |
328 | ax = 0.0; bx = tsol; cx = 1.0; | |
329 | } | |
89d8607f | 330 | FSR_DEBUG(" minimisation dans la direction"); |
859a47c3 | 331 | |
332 | math_BrentMinimum Sol(tol1d, 100, tol1d); | |
859a47c3 | 333 | |
e93e4230 | 334 | // Base invocation. |
335 | Sol.Perform(F, ax, bx, cx); | |
336 | if(Sol.IsDone()) | |
337 | { | |
338 | if (Sol.Minimum() <= Result) | |
339 | { | |
3e42bd70 J |
340 | tsol = Sol.Location(); |
341 | good = Standard_True; | |
e93e4230 | 342 | Result = Sol.Minimum(); |
343 | ||
344 | // Objective function changes too fast -> | |
345 | // perform additional computations. | |
346 | if (Gradient.Norm2() > 1.0 / Precision::SquareConfusion() && | |
347 | tsol > ax && | |
348 | tsol < cx) // Solution inside of (ax, cx) interval. | |
349 | { | |
350 | // First and second part invocation. | |
351 | Sol.Perform(F, ax, (ax + tsol) / 2.0, tsol); | |
352 | if(Sol.IsDone()) | |
353 | { | |
354 | if (Sol.Minimum() <= Result) | |
355 | { | |
356 | tsol = Sol.Location(); | |
357 | good = Standard_True; | |
358 | Result = Sol.Minimum(); | |
359 | } | |
360 | } | |
361 | ||
362 | Sol.Perform(F, tsol, (cx + tsol) / 2.0, cx); | |
363 | if(Sol.IsDone()) | |
364 | { | |
365 | if (Sol.Minimum() <= Result) | |
366 | { | |
367 | tsol = Sol.Location(); | |
368 | good = Standard_True; | |
369 | Result = Sol.Minimum(); | |
370 | } | |
371 | } | |
372 | } | |
373 | } // if (Sol.Minimum() <= Result) | |
374 | } // if(Sol.IsDone()) | |
7fd59977 | 375 | } |
e93e4230 | 376 | |
377 | if (good) | |
378 | { | |
7fd59977 | 379 | // mise a jour du Delta |
380 | Dir.Multiply(tsol); | |
381 | } | |
382 | return good; | |
383 | } | |
384 | ||
385 | //------------------------------------------------------ | |
386 | static void SearchDirection(const math_Matrix& DF, | |
89d8607f | 387 | const math_Vector& GH, |
388 | const math_Vector& FF, | |
389 | Standard_Boolean ChangeDirection, | |
390 | const math_Vector& InvLengthMax, | |
391 | math_Vector& Direction, | |
392 | Standard_Real& Dy) | |
7fd59977 | 393 | |
394 | { | |
395 | Standard_Integer Ninc = DF.ColNumber(), Neq = DF.RowNumber(); | |
396 | Standard_Real Eps = 1.e-32; | |
397 | if (!ChangeDirection) { | |
398 | if (Ninc == Neq) { | |
399 | for (Standard_Integer i = FF.Lower(); i <= FF.Upper(); i++) { | |
89d8607f | 400 | Direction(i) = -FF(i); |
7fd59977 | 401 | } |
402 | math_Gauss Solut(DF, 1.e-9); | |
403 | if (Solut.IsDone()) Solut.Solve(Direction); | |
404 | else { // we have to "forget" singular directions. | |
89d8607f | 405 | FSR_DEBUG(" Matrice singuliere : On prend SVD"); |
3e42bd70 | 406 | math_SVD SolvebySVD(DF); |
7fd59977 | 407 | if (SolvebySVD.IsDone()) SolvebySVD.Solve(-1*FF, Direction); |
3e42bd70 J |
408 | else ChangeDirection = Standard_True; |
409 | } | |
7fd59977 | 410 | } |
411 | else if (Ninc > Neq) { | |
412 | math_SVD Solut(DF); | |
413 | if (Solut.IsDone()) Solut.Solve(-1*FF, Direction); | |
414 | else ChangeDirection = Standard_True; | |
415 | } | |
416 | else if (Ninc < Neq) { // Calcul par GaussLeastSquare | |
417 | math_GaussLeastSquare Solut(DF); | |
418 | if (Solut.IsDone()) Solut.Solve(-1*FF, Direction); | |
419 | else ChangeDirection = Standard_True; | |
420 | } | |
421 | } | |
422 | // Il vaut mieux interdire des directions trops longue | |
423 | // Afin de blinder les cas trop mal conditionne | |
424 | // PMN 12/05/97 Traitement des singularite dans les conges | |
425 | // Sur des surfaces periodiques | |
89d8607f | 426 | |
7fd59977 | 427 | Standard_Real ratio = Abs( Direction(Direction.Lower()) |
89d8607f | 428 | *InvLengthMax(Direction.Lower()) ); |
7fd59977 | 429 | Standard_Integer i; |
430 | for (i = Direction.Lower()+1; i<=Direction.Upper(); i++) { | |
431 | ratio = Max(ratio, Abs( Direction(i)*InvLengthMax(i)) ); | |
432 | } | |
433 | if (ratio > 1) { | |
434 | Direction /= ratio; | |
435 | } | |
89d8607f | 436 | |
7fd59977 | 437 | Dy = Direction*GH; |
438 | if (Dy >= -Eps) { // newton "ne descend pas" on prend le gradient | |
439 | ChangeDirection = Standard_True; | |
440 | } | |
441 | if (ChangeDirection) { // On va faire un gradient ! | |
442 | for (i = Direction.Lower(); i <= Direction.Upper(); i++) { | |
443 | Direction(i) = - GH(i); | |
444 | } | |
445 | Dy = - (GH.Norm2()); | |
446 | } | |
447 | } | |
448 | ||
449 | ||
450 | //===================================================================== | |
451 | static void SearchDirection(const math_Matrix& DF, | |
89d8607f | 452 | const math_Vector& GH, |
453 | const math_Vector& FF, | |
454 | const math_IntegerVector& Constraints, | |
455 | // const math_Vector& X, // Le point d'init | |
456 | const math_Vector& , // Le point d'init | |
457 | Standard_Boolean ChangeDirection, | |
7fd59977 | 458 | const math_Vector& InvLengthMax, |
89d8607f | 459 | math_Vector& Direction, |
460 | Standard_Real& Dy) | |
461 | //Purpose : Recherche une direction (et un pas si Newton Fonctionne) le long | |
462 | // d'une frontiere | |
463 | //===================================================================== | |
7fd59977 | 464 | { |
465 | Standard_Integer Ninc = DF.ColNumber(), Neq = DF.RowNumber(); | |
466 | Standard_Integer i, j, k, Cons = 0; | |
467 | ||
468 | // verification sur les bornes imposees: | |
469 | ||
470 | for (i = 1; i <= Ninc; i++) { | |
471 | if (Constraints(i) != 0) Cons++; | |
472 | // sinon le systeme a resoudre ne change pas. | |
473 | } | |
474 | ||
475 | if (Cons == 0) { | |
476 | SearchDirection(DF, GH, FF, ChangeDirection, InvLengthMax, | |
89d8607f | 477 | Direction, Dy); |
7fd59977 | 478 | } |
479 | else if (Cons == Ninc) { // il n'y a plus rien a faire... | |
51740958 | 480 | for(i = Direction.Lower(); i <= Direction.Upper(); i++) { |
89d8607f | 481 | Direction(i) = 0; |
482 | } | |
7fd59977 | 483 | Dy = 0; |
484 | } | |
485 | else { //(1) Cas general : On definit un sous probleme | |
486 | math_Matrix DF2(1, Neq, 1, Ninc-Cons); | |
487 | math_Vector MyGH(1, Ninc-Cons); | |
488 | math_Vector MyDirection(1, Ninc-Cons); | |
489 | math_Vector MyInvLengthMax(1, Ninc); | |
490 | ||
491 | for (k=1, i = 1; i <= Ninc; i++) { | |
492 | if (Constraints(i) == 0) { | |
89d8607f | 493 | MyGH(k) = GH(i); |
494 | MyInvLengthMax(k) = InvLengthMax(i); | |
495 | MyDirection(k) = Direction(i); | |
496 | for (j = 1; j <= Neq; j++) { | |
497 | DF2(j, k) = DF(j, i); | |
498 | } | |
499 | k++; //on passe a l'inconnue suivante | |
500 | } | |
7fd59977 | 501 | } |
502 | //(2) On le resoud | |
503 | SearchDirection(DF2, MyGH, FF, ChangeDirection, MyInvLengthMax, | |
89d8607f | 504 | MyDirection, Dy); |
7fd59977 | 505 | |
506 | // (3) On l'interprete... | |
507 | // Reconstruction de Direction: | |
508 | for (i = 1, k = 1; i <= Ninc; i++) { | |
509 | if (Constraints(i) == 0) { | |
89d8607f | 510 | if (!ChangeDirection) { |
511 | Direction(i) = MyDirection(k); | |
512 | } | |
513 | else Direction(i) = - GH(i); | |
514 | k++; | |
7fd59977 | 515 | } |
516 | else { | |
89d8607f | 517 | Direction(i) = 0.0; |
7fd59977 | 518 | } |
519 | } | |
520 | } | |
521 | } | |
522 | ||
523 | ||
524 | ||
525 | //==================================================== | |
3e42bd70 J |
526 | Standard_Boolean Bounds(const math_Vector& InfBound, |
527 | const math_Vector& SupBound, | |
528 | const math_Vector& Tol, | |
529 | math_Vector& Sol, | |
530 | const math_Vector& SolSave, | |
531 | math_IntegerVector& Constraints, | |
532 | math_Vector& Delta, | |
533 | Standard_Boolean& theIsNewSol) | |
89d8607f | 534 | // |
535 | // Purpose: Trims an initial solution Sol to be within a domain defined by | |
536 | // InfBound and SupBound. Delta will contain a distance between final Sol and | |
537 | // SolSave. | |
538 | // IsNewSol returns False, if final Sol fully coincides with SolSave, i.e. | |
539 | // if SolSave already lied on a boundary and initial Sol was fully beyond it | |
540 | //====================================================== | |
7fd59977 | 541 | { |
542 | Standard_Boolean Out = Standard_False; | |
543 | Standard_Integer i, Ninc = Sol.Length(); | |
544 | Standard_Real monratio = 1; | |
89d8607f | 545 | |
3e42bd70 J |
546 | theIsNewSol = Standard_True; |
547 | ||
7fd59977 | 548 | // Calcul du ratio de recadrage |
549 | for (i = 1; i <= Ninc; i++) { | |
550 | Constraints(i) = 0; | |
551 | Delta(i) = Sol(i) - SolSave(i); | |
552 | if (InfBound(i) == SupBound(i)) { | |
553 | Constraints(i) = 1; | |
554 | Out = Standard_True; // Ok mais, cela devrait etre eviter | |
555 | } | |
556 | else if(Sol(i) < InfBound(i)) { | |
557 | Constraints(i) = 1; | |
558 | Out = Standard_True; | |
3e42bd70 J |
559 | // Delta(i) is negative |
560 | if (-Delta(i) > Tol(i)) // Afin d'eviter des ratio nulles pour rien | |
561 | monratio = Min(monratio, (InfBound(i) - SolSave(i))/Delta(i) ); | |
7fd59977 | 562 | } |
563 | else if (Sol(i) > SupBound(i)) { | |
564 | Constraints(i) = 1; | |
565 | Out = Standard_True; | |
3e42bd70 J |
566 | // Delta(i) is positive |
567 | if (Delta(i) > Tol(i)) | |
568 | monratio = Min(monratio, (SupBound(i) - SolSave(i))/Delta(i) ); | |
7fd59977 | 569 | } |
570 | } | |
571 | ||
572 | if (Out){ // Troncature et derniers recadrage pour blinder (pb numeriques) | |
3e42bd70 J |
573 | if (monratio == 0.0) { |
574 | theIsNewSol = Standard_False; | |
575 | Sol = SolSave; | |
576 | Delta.Init (0.0); | |
577 | } else { | |
578 | Delta *= monratio; | |
579 | Sol = SolSave+Delta; | |
580 | for (i = 1; i <= Ninc; i++) { | |
581 | if(Sol(i) < InfBound(i)) { | |
582 | Sol(i) = InfBound(i); | |
583 | Delta(i) = Sol(i) - SolSave(i); | |
584 | } | |
585 | else if (Sol(i) > SupBound(i)) { | |
586 | Sol(i) = SupBound(i); | |
587 | Delta(i) = Sol(i) - SolSave(i); | |
588 | } | |
7fd59977 | 589 | } |
590 | } | |
591 | } | |
592 | return Out; | |
593 | } | |
594 | ||
595 | ||
596 | ||
597 | ||
859a47c3 | 598 | //======================================================================= |
599 | //function : math_FunctionSetRoot | |
600 | //purpose : Constructor | |
601 | //======================================================================= | |
602 | math_FunctionSetRoot::math_FunctionSetRoot( | |
603 | math_FunctionSetWithDerivatives& theFunction, | |
604 | const math_Vector& theTolerance, | |
605 | const Standard_Integer theNbIterations) | |
606 | ||
607 | : Delta(1, theFunction.NbVariables()), | |
608 | Sol (1, theFunction.NbVariables()), | |
609 | DF (1, theFunction.NbEquations() , 1, theFunction.NbVariables()), | |
610 | Tol (1, theFunction.NbVariables()), | |
611 | Done (Standard_False), | |
612 | Kount (0), | |
613 | State (0), | |
614 | Itermax (theNbIterations), | |
615 | InfBound(1, theFunction.NbVariables(), RealFirst()), | |
616 | SupBound(1, theFunction.NbVariables(), RealLast ()), | |
617 | SolSave (1, theFunction.NbVariables()), | |
618 | GH (1, theFunction.NbVariables()), | |
619 | DH (1, theFunction.NbVariables()), | |
620 | DHSave (1, theFunction.NbVariables()), | |
621 | FF (1, theFunction.NbEquations()), | |
622 | PreviousSolution(1, theFunction.NbVariables()), | |
623 | Save (0, theNbIterations), | |
624 | Constraints(1, theFunction.NbVariables()), | |
625 | Temp1 (1, theFunction.NbVariables()), | |
626 | Temp2 (1, theFunction.NbVariables()), | |
627 | Temp3 (1, theFunction.NbVariables()), | |
628 | Temp4 (1, theFunction.NbEquations()), | |
629 | myIsDivergent(Standard_False) | |
7fd59977 | 630 | { |
859a47c3 | 631 | SetTolerance(theTolerance); |
7fd59977 | 632 | } |
633 | ||
859a47c3 | 634 | //======================================================================= |
635 | //function : math_FunctionSetRoot | |
636 | //purpose : Constructor | |
637 | //======================================================================= | |
638 | math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& theFunction, | |
639 | const Standard_Integer theNbIterations) | |
640 | ||
641 | : Delta(1, theFunction.NbVariables()), | |
642 | Sol (1, theFunction.NbVariables()), | |
643 | DF (1, theFunction.NbEquations() , 1, theFunction.NbVariables()), | |
644 | Tol (1, theFunction.NbVariables()), | |
645 | Done (Standard_False), | |
646 | Kount (0), | |
647 | State (0), | |
648 | Itermax (theNbIterations), | |
649 | InfBound(1, theFunction.NbVariables(), RealFirst()), | |
650 | SupBound(1, theFunction.NbVariables(), RealLast ()), | |
651 | SolSave (1, theFunction.NbVariables()), | |
652 | GH (1, theFunction.NbVariables()), | |
653 | DH (1, theFunction.NbVariables()), | |
654 | DHSave (1, theFunction.NbVariables()), | |
655 | FF (1, theFunction.NbEquations()), | |
656 | PreviousSolution(1, theFunction.NbVariables()), | |
657 | Save (0, theNbIterations), | |
658 | Constraints(1, theFunction.NbVariables()), | |
659 | Temp1 (1, theFunction.NbVariables()), | |
660 | Temp2 (1, theFunction.NbVariables()), | |
661 | Temp3 (1, theFunction.NbVariables()), | |
662 | Temp4 (1, theFunction.NbEquations()), | |
663 | myIsDivergent(Standard_False) | |
7fd59977 | 664 | { |
7fd59977 | 665 | } |
666 | ||
859a47c3 | 667 | //======================================================================= |
668 | //function : ~math_FunctionSetRoot | |
669 | //purpose : Destructor | |
670 | //======================================================================= | |
671 | math_FunctionSetRoot::~math_FunctionSetRoot() | |
89d8607f | 672 | { |
7fd59977 | 673 | } |
674 | ||
859a47c3 | 675 | //======================================================================= |
676 | //function : SetTolerance | |
677 | //purpose : | |
678 | //======================================================================= | |
679 | void math_FunctionSetRoot::SetTolerance(const math_Vector& theTolerance) | |
6da30ff1 | 680 | { |
859a47c3 | 681 | for (Standard_Integer i = 1; i <= Tol.Length(); ++i) |
682 | Tol(i) = theTolerance(i); | |
6da30ff1 | 683 | } |
7fd59977 | 684 | |
859a47c3 | 685 | //======================================================================= |
686 | //function : Perform | |
687 | //purpose : | |
688 | //======================================================================= | |
689 | void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& theFunction, | |
690 | const math_Vector& theStartingPoint, | |
691 | const Standard_Boolean theStopOnDivergent) | |
7fd59977 | 692 | { |
859a47c3 | 693 | Perform(theFunction, theStartingPoint, InfBound, SupBound, theStopOnDivergent); |
7fd59977 | 694 | } |
695 | ||
859a47c3 | 696 | //======================================================================= |
697 | //function : Perform | |
698 | //purpose : | |
699 | //======================================================================= | |
7fd59977 | 700 | void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F, |
89d8607f | 701 | const math_Vector& StartingPoint, |
75259fc5 | 702 | const math_Vector& theInfBound, |
703 | const math_Vector& theSupBound, | |
89d8607f | 704 | Standard_Boolean theStopOnDivergent) |
7fd59977 | 705 | { |
7fd59977 | 706 | Standard_Integer Ninc = F.NbVariables(), Neq = F.NbEquations(); |
89d8607f | 707 | |
7fd59977 | 708 | if ((Neq <= 0) || |
89d8607f | 709 | (StartingPoint.Length()!= Ninc) || |
75259fc5 | 710 | (theInfBound.Length() != Ninc) || |
711 | (theSupBound.Length() != Ninc)) { Standard_DimensionError:: Raise(); } | |
7fd59977 | 712 | |
713 | Standard_Integer i; | |
3e42bd70 | 714 | Standard_Boolean ChangeDirection = Standard_False, Sort = Standard_False, isNewSol = Standard_False; |
7fd59977 | 715 | Standard_Boolean Good, Verif; |
716 | Standard_Boolean Stop; | |
3e42bd70 | 717 | const Standard_Real EpsSqrt = 1.e-16, Eps = 1.e-32, Eps2 = 1.e-64, Progres = 0.005; |
7fd59977 | 718 | Standard_Real F2, PreviousMinimum, Dy, OldF; |
719 | Standard_Real Ambda, Ambda2, Gnr1, Oldgr; | |
720 | math_Vector InvLengthMax(1, Ninc); // Pour bloquer les pas a 1/4 du domaine | |
75259fc5 | 721 | math_IntegerVector aConstraints(1, Ninc); // Pour savoir sur quels bord on se trouve |
7fd59977 | 722 | for (i = 1; i <= Ninc ; i++) { |
f4dee9bb | 723 | const Standard_Real aSupBound = Min (theSupBound(i), Precision::Infinite()); |
724 | const Standard_Real anInfBound = Max (theInfBound(i), -Precision::Infinite()); | |
725 | InvLengthMax(i) = 1. / Max((aSupBound - anInfBound)/4, 1.e-9); | |
89d8607f | 726 | } |
7fd59977 | 727 | |
728 | MyDirFunction F_Dir(Temp1, Temp2, Temp3, Temp4, F); | |
729 | Standard_Integer DescenteIter; | |
730 | ||
731 | Done = Standard_False; | |
732 | Sol = StartingPoint; | |
733 | Kount = 0; | |
734 | ||
b659a6dc | 735 | // |
736 | myIsDivergent = Standard_False; | |
737 | for (i = 1; i <= Ninc; i++) | |
738 | { | |
75259fc5 | 739 | myIsDivergent |= (Sol(i) < theInfBound(i)) | (theSupBound(i) < Sol(i)); |
b659a6dc | 740 | } |
741 | if (theStopOnDivergent & myIsDivergent) | |
742 | { | |
743 | return; | |
744 | } | |
7fd59977 | 745 | // Verification de la validite des inconnues par rapport aux bornes. |
746 | // Recentrage sur les bornes si pas valide. | |
747 | for ( i = 1; i <= Ninc; i++) { | |
75259fc5 | 748 | if (Sol(i) <= theInfBound(i)) Sol(i) = theInfBound(i); |
749 | else if (Sol(i) > theSupBound(i)) Sol(i) = theSupBound(i); | |
7fd59977 | 750 | } |
751 | ||
752 | // Calcul de la premiere valeur de F et de son gradient | |
753 | if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) { | |
754 | Done = Standard_False; | |
b659a6dc | 755 | if (!theStopOnDivergent || !myIsDivergent) |
756 | { | |
757 | State = F.GetStateNumber(); | |
758 | } | |
7fd59977 | 759 | return; |
760 | } | |
761 | Ambda2 = Gnr1; | |
762 | // Le rang 0 de Save ne doit servir q'au test accelarteur en fin de boucle | |
763 | // s'il on est dejas sur la solution, il faut leurer ce test pour eviter | |
764 | // de faire une seconde iteration... | |
3e42bd70 | 765 | Save(0) = Max (F2, EpsSqrt); |
5368adff | 766 | Standard_Real aTol_Func = Epsilon(F2); |
89d8607f | 767 | FSR_DEBUG("=== Mode Debug de Function Set Root" << endl); |
768 | FSR_DEBUG(" F2 Initial = " << F2); | |
7fd59977 | 769 | |
3e42bd70 | 770 | if ((F2 <= Eps) || (Gnr1 <= Eps2)) { |
b659a6dc | 771 | Done = Standard_False; |
772 | if (!theStopOnDivergent || !myIsDivergent) | |
773 | { | |
774 | Done = Standard_True; | |
775 | State = F.GetStateNumber(); | |
776 | } | |
7fd59977 | 777 | return; |
778 | } | |
779 | ||
780 | for (Kount = 1; Kount <= Itermax; Kount++) { | |
781 | PreviousMinimum = F2; | |
782 | Oldgr = Gnr1; | |
783 | PreviousSolution = Sol; | |
784 | SolSave = Sol; | |
785 | ||
786 | SearchDirection(DF, GH, FF, ChangeDirection, InvLengthMax, DH, Dy); | |
787 | if (Abs(Dy) <= Eps) { | |
b659a6dc | 788 | Done = Standard_False; |
789 | if (!theStopOnDivergent || !myIsDivergent) | |
790 | { | |
791 | Done = Standard_True; | |
792 | ////modified by jgv, 31.08.2011//// | |
793 | F.Value(Sol, FF); //update F before GetStateNumber | |
794 | /////////////////////////////////// | |
795 | State = F.GetStateNumber(); | |
796 | } | |
7fd59977 | 797 | return; |
798 | } | |
799 | if (ChangeDirection) { | |
3e42bd70 | 800 | Ambda = Ambda2 / Sqrt(Abs(Dy)); |
7fd59977 | 801 | if (Ambda > 1.0) Ambda = 1.0; |
802 | } | |
803 | else { | |
804 | Ambda = 1.0; | |
805 | Ambda2 = 0.5*Ambda/DH.Norm(); | |
806 | } | |
807 | ||
808 | for( i = Sol.Lower(); i <= Sol.Upper(); i++) { | |
809 | Sol(i) = Sol(i) + Ambda * DH(i); | |
810 | } | |
b659a6dc | 811 | // |
812 | for (i = 1; i <= Ninc; i++) | |
813 | { | |
75259fc5 | 814 | myIsDivergent |= (Sol(i) < theInfBound(i)) | (theSupBound(i) < Sol(i)); |
b659a6dc | 815 | } |
816 | if (theStopOnDivergent & myIsDivergent) | |
817 | { | |
818 | return; | |
819 | } | |
820 | // | |
75259fc5 | 821 | Sort = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave, |
822 | aConstraints, Delta, isNewSol); | |
89d8607f | 823 | |
7fd59977 | 824 | |
7fd59977 | 825 | DHSave = GH; |
3e42bd70 | 826 | if (isNewSol) { |
89d8607f | 827 | // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1); |
3e42bd70 J |
828 | if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) { |
829 | Done = Standard_False; | |
b659a6dc | 830 | if (!theStopOnDivergent || !myIsDivergent) |
831 | { | |
832 | State = F.GetStateNumber(); | |
833 | } | |
3e42bd70 J |
834 | return; |
835 | } | |
7fd59977 | 836 | } |
89d8607f | 837 | |
838 | FSR_DEBUG("Kount = " << Kount); | |
839 | FSR_DEBUG("Le premier F2 = " << F2); | |
840 | FSR_DEBUG("Direction = " << ChangeDirection); | |
7fd59977 | 841 | |
3e42bd70 | 842 | if ((F2 <= Eps) || (Gnr1 <= Eps2)) { |
b659a6dc | 843 | Done = Standard_False; |
844 | if (!theStopOnDivergent || !myIsDivergent) | |
845 | { | |
846 | Done = Standard_True; | |
847 | ////modified by jgv, 31.08.2011//// | |
848 | F.Value(Sol, FF); //update F before GetStateNumber | |
849 | /////////////////////////////////// | |
850 | State = F.GetStateNumber(); | |
851 | } | |
7fd59977 | 852 | return; |
853 | } | |
854 | ||
855 | if (Sort || (F2/PreviousMinimum > Progres)) { | |
856 | Dy = GH*DH; | |
857 | OldF = PreviousMinimum; | |
858 | Stop = Standard_False; | |
859 | Good = Standard_False; | |
860 | DescenteIter = 0; | |
861 | Standard_Boolean Sortbis; | |
862 | ||
863 | // ------------------------------------------------- | |
864 | // Traitement standard sans traitement des bords | |
865 | // ------------------------------------------------- | |
866 | if (!Sort) { // si l'on est pas sortie on essaye de progresser en avant | |
89d8607f | 867 | while((F2/PreviousMinimum > Progres) && !Stop) { |
868 | if (F2 < OldF && (Dy < 0.0)) { | |
869 | // On essaye de progresser dans cette direction. | |
870 | FSR_DEBUG(" iteration de descente = " << DescenteIter); | |
871 | DescenteIter++; | |
872 | SolSave = Sol; | |
873 | OldF = F2; | |
874 | for( i = Sol.Lower(); i <= Sol.Upper(); i++) { | |
875 | Sol(i) = Sol(i) + Ambda * DH(i); | |
876 | } | |
877 | // | |
878 | for (i = 1; i <= Ninc; i++) | |
879 | { | |
75259fc5 | 880 | myIsDivergent |= (Sol(i) < theInfBound(i)) | (theSupBound(i) < Sol(i)); |
89d8607f | 881 | } |
882 | if (theStopOnDivergent & myIsDivergent) | |
883 | { | |
884 | return; | |
885 | } | |
886 | // | |
75259fc5 | 887 | Stop = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave, |
888 | aConstraints, Delta, isNewSol); | |
89d8607f | 889 | FSR_DEBUG(" Augmentation de lambda"); |
890 | Ambda *= 1.7; | |
b659a6dc | 891 | } |
89d8607f | 892 | else { |
893 | if ((F2 >= OldF)||(F2 >= PreviousMinimum)) { | |
894 | Good = Standard_False; | |
895 | if (DescenteIter == 0) { | |
896 | // C'est le premier pas qui flanche, on fait une interpolation. | |
897 | // et on minimise si necessaire. | |
898 | DescenteIter++; | |
899 | Good = MinimizeDirection(SolSave, Delta, OldF, F2, DHSave, GH, | |
900 | Tol, F_Dir); | |
901 | } | |
902 | else if (ChangeDirection || (DescenteIter>1) | |
903 | || (OldF>PreviousMinimum) ) { | |
904 | // La progression a ete utile, on minimise... | |
905 | DescenteIter++; | |
906 | Good = MinimizeDirection(PreviousSolution, SolSave, Sol, | |
907 | OldF, Delta, Tol, F_Dir); | |
908 | } | |
909 | if (!Good) { | |
910 | Sol = SolSave; | |
911 | F2 = OldF; | |
912 | } | |
913 | else { | |
914 | Sol = SolSave+Delta; | |
915 | // | |
916 | for (i = 1; i <= Ninc; i++) | |
917 | { | |
75259fc5 | 918 | myIsDivergent |= (Sol(i) < theInfBound(i)) | (theSupBound(i) < Sol(i)); |
89d8607f | 919 | } |
920 | if (theStopOnDivergent & myIsDivergent) | |
921 | { | |
922 | return; | |
923 | } | |
924 | // | |
75259fc5 | 925 | Sort = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave, |
926 | aConstraints, Delta, isNewSol); | |
89d8607f | 927 | } |
928 | Sort = Standard_False; // On a rejete le point sur la frontiere | |
929 | } | |
930 | Stop = Standard_True; // et on sort dans tous les cas... | |
b659a6dc | 931 | } |
89d8607f | 932 | DHSave = GH; |
3e42bd70 | 933 | if (isNewSol) { |
89d8607f | 934 | // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1); |
3e42bd70 J |
935 | if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) { |
936 | Done = Standard_False; | |
b659a6dc | 937 | if (!theStopOnDivergent || !myIsDivergent) |
938 | { | |
939 | State = F.GetStateNumber(); | |
940 | } | |
3e42bd70 J |
941 | return; |
942 | } | |
943 | } | |
89d8607f | 944 | Dy = GH*DH; |
945 | if (Abs(Dy) <= Eps) { | |
946 | if (F2 > OldF) | |
3e42bd70 | 947 | Sol = SolSave; |
89d8607f | 948 | Done = Standard_False; |
949 | if (!theStopOnDivergent || !myIsDivergent) | |
950 | { | |
951 | Done = Standard_True; | |
b659a6dc | 952 | ////modified by jgv, 31.08.2011//// |
953 | F.Value(Sol, FF); //update F before GetStateNumber | |
954 | /////////////////////////////////// | |
89d8607f | 955 | State = F.GetStateNumber(); |
956 | } | |
957 | return; | |
958 | } | |
959 | if (DescenteIter >= 100) { | |
960 | Stop = Standard_True; | |
961 | } | |
962 | } | |
963 | FSR_DEBUG("--- Sortie du Traitement Standard"); | |
964 | FSR_DEBUG(" DescenteIter = "<<DescenteIter << " F2 = " << F2); | |
7fd59977 | 965 | } |
966 | // ------------------------------------ | |
967 | // on passe au traitement des bords | |
968 | // ------------------------------------ | |
969 | if (Sort) { | |
89d8607f | 970 | Stop = (F2 > 1.001*OldF); // Pour ne pas progresser sur le bord |
971 | Sortbis = Sort; | |
972 | DescenteIter = 0; | |
973 | while (Sortbis && ((F2<OldF)|| (DescenteIter == 0)) | |
974 | && (!Stop)) { | |
975 | DescenteIter++; | |
976 | // On essaye de progresser sur le bord | |
977 | SolSave = Sol; | |
978 | OldF = F2; | |
75259fc5 | 979 | SearchDirection(DF, GH, FF, aConstraints, Sol, |
89d8607f | 980 | ChangeDirection, InvLengthMax, DH, Dy); |
981 | FSR_DEBUG(" Conditional Direction = " << ChangeDirection); | |
982 | if (Dy<-Eps) { //Pour eviter des calculs inutiles et des /0... | |
983 | if (ChangeDirection) { | |
984 | ||
985 | // Ambda = Ambda2 / Sqrt(Abs(Dy)); | |
986 | Ambda = Ambda2 / Sqrt(-Dy); | |
987 | if (Ambda > 1.0) Ambda = 1.0; | |
988 | } | |
989 | else { | |
990 | Ambda = 1.0; | |
991 | Ambda2 = 0.5*Ambda/DH.Norm(); | |
992 | } | |
993 | ||
994 | for( i = Sol.Lower(); i <= Sol.Upper(); i++) { | |
995 | Sol(i) = Sol(i) + Ambda * DH(i); | |
996 | } | |
997 | // | |
998 | for (i = 1; i <= Ninc; i++) | |
999 | { | |
75259fc5 | 1000 | myIsDivergent |= (Sol(i) < theInfBound(i)) | (theSupBound(i) < Sol(i)); |
89d8607f | 1001 | } |
1002 | if (theStopOnDivergent & myIsDivergent) | |
1003 | { | |
3e42bd70 J |
1004 | return; |
1005 | } | |
89d8607f | 1006 | // |
75259fc5 | 1007 | Sortbis = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave, |
1008 | aConstraints, Delta, isNewSol); | |
89d8607f | 1009 | |
1010 | DHSave = GH; | |
1011 | if (isNewSol) { | |
1012 | // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1); | |
1013 | if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) { | |
1014 | Done = Standard_False; | |
1015 | if (!theStopOnDivergent || !myIsDivergent) | |
1016 | { | |
1017 | State = F.GetStateNumber(); | |
1018 | } | |
1019 | return; | |
b659a6dc | 1020 | } |
3e42bd70 | 1021 | } |
89d8607f | 1022 | Ambda2 = Gnr1; |
1023 | FSR_DEBUG("--- Iteration au bords : " << DescenteIter); | |
1024 | FSR_DEBUG("--- F2 = " << F2); | |
3e42bd70 | 1025 | } |
89d8607f | 1026 | else { |
1027 | Stop = Standard_True; | |
1028 | } | |
1029 | ||
1030 | while((F2/PreviousMinimum > Progres) && (F2<OldF) && (!Stop) ) { | |
1031 | DescenteIter++; | |
1032 | FSR_DEBUG("--- Iteration de descente conditionnel = " << DescenteIter); | |
1033 | if (F2 < OldF && Dy < 0.0) { | |
1034 | // On essaye de progresser dans cette direction. | |
1035 | SolSave = Sol; | |
1036 | OldF = F2; | |
1037 | for( i = Sol.Lower(); i <= Sol.Upper(); i++) { | |
1038 | Sol(i) = Sol(i) + Ambda * DH(i); | |
1039 | } | |
1040 | // | |
1041 | for (i = 1; i <= Ninc; i++) | |
b659a6dc | 1042 | { |
75259fc5 | 1043 | myIsDivergent |= (Sol(i) < theInfBound(i)) | (theSupBound(i) < Sol(i)); |
b659a6dc | 1044 | } |
89d8607f | 1045 | if (theStopOnDivergent & myIsDivergent) |
1046 | { | |
1047 | return; | |
1048 | } | |
1049 | // | |
75259fc5 | 1050 | Sortbis = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave, |
1051 | aConstraints, Delta, isNewSol); | |
89d8607f | 1052 | } |
1053 | DHSave = GH; | |
1054 | if (isNewSol) { | |
1055 | // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1); | |
1056 | if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) { | |
1057 | Done = Standard_False; | |
1058 | if (!theStopOnDivergent || !myIsDivergent) | |
1059 | { | |
1060 | State = F.GetStateNumber(); | |
1061 | } | |
1062 | return; | |
1063 | } | |
1064 | } | |
1065 | Ambda2 = Gnr1; | |
1066 | Dy = GH*DH; | |
1067 | Stop = ((Dy >=0) || (DescenteIter >= 10) || Sortbis); | |
1068 | } | |
1069 | Stop = ((Dy >=0) || (DescenteIter >= 10)); | |
1070 | } | |
1071 | if (((F2/PreviousMinimum > Progres) && | |
1072 | (F2>=OldF))||(F2>=PreviousMinimum)) { | |
1073 | // On minimise par Brent | |
1074 | DescenteIter++; | |
1075 | Good = MinimizeDirection(SolSave, Delta, OldF, F2, | |
1076 | DHSave, GH, Tol, F_Dir); | |
1077 | if (!Good) { | |
1078 | Sol = SolSave; | |
1079 | Sort = Standard_False; | |
1080 | } | |
1081 | else { | |
1082 | Sol = SolSave + Delta; | |
1083 | // | |
1084 | for (i = 1; i <= Ninc; i++) | |
1085 | { | |
75259fc5 | 1086 | myIsDivergent |= (Sol(i) < theInfBound(i)) | (theSupBound(i) < Sol(i)); |
89d8607f | 1087 | } |
1088 | if (theStopOnDivergent & myIsDivergent) | |
1089 | { | |
3e42bd70 J |
1090 | return; |
1091 | } | |
89d8607f | 1092 | // |
75259fc5 | 1093 | Sort = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave, |
1094 | aConstraints, Delta, isNewSol); | |
89d8607f | 1095 | if (isNewSol) { |
1096 | // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1); | |
1097 | if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) { | |
1098 | Done = Standard_False; | |
1099 | if (!theStopOnDivergent || !myIsDivergent) | |
1100 | { | |
1101 | State = F.GetStateNumber(); | |
1102 | } | |
1103 | return; | |
1104 | } | |
1105 | } | |
1106 | } | |
1107 | Dy = GH*DH; | |
1108 | } | |
1109 | FSR_DEBUG("--- Sortie du Traitement des Bords"); | |
1110 | FSR_DEBUG("--- DescenteIter = "<<DescenteIter << " F2 = " << F2); | |
7fd59977 | 1111 | } |
1112 | } | |
1113 | ||
1114 | // --------------------------------------------- | |
1115 | // on passe aux tests d'ARRET | |
1116 | // --------------------------------------------- | |
1117 | Save(Kount) = F2; | |
1118 | // Est ce la solution ? | |
1119 | if (ChangeDirection) Verif = Standard_True; | |
89d8607f | 1120 | // Gradient : Il faut eviter de boucler |
7fd59977 | 1121 | else { |
1122 | Verif = Standard_False; | |
1123 | if (Kount > 1) { // Pour accelerer les cas quasi-quadratique | |
89d8607f | 1124 | if (Save(Kount-1)<1.e-4*Save(Kount-2)) Verif = Standard_True; |
7fd59977 | 1125 | } |
1126 | else Verif = (F2 < 1.e-6*Save(0)); //Pour les cas dejas solutions | |
1127 | } | |
1128 | if (Verif) { | |
1129 | for(i = Delta.Lower(); i <= Delta.Upper(); i++) { | |
89d8607f | 1130 | Delta(i) = PreviousSolution(i) - Sol(i); |
7fd59977 | 1131 | } |
89d8607f | 1132 | |
7fd59977 | 1133 | if (IsSolutionReached(F)) { |
89d8607f | 1134 | if (PreviousMinimum < F2) { |
1135 | Sol = SolSave; | |
1136 | } | |
1137 | Done = Standard_False; | |
1138 | if (!theStopOnDivergent || !myIsDivergent) | |
1139 | { | |
1140 | Done = Standard_True; | |
b659a6dc | 1141 | ////modified by jgv, 31.08.2011//// |
1142 | F.Value(Sol, FF); //update F before GetStateNumber | |
1143 | /////////////////////////////////// | |
89d8607f | 1144 | State = F.GetStateNumber(); |
1145 | } | |
1146 | return; | |
7fd59977 | 1147 | } |
1148 | } | |
1149 | //fin du test solution | |
89d8607f | 1150 | |
7fd59977 | 1151 | // Analyse de la progression... |
5368adff | 1152 | //comparison of current minimum and previous minimum |
1153 | if ((F2 - PreviousMinimum) <= aTol_Func){ | |
7fd59977 | 1154 | if (Kount > 5) { |
89d8607f | 1155 | // L'historique est il bon ? |
1156 | if (F2 >= 0.95*Save(Kount - 5)) { | |
1157 | if (!ChangeDirection) ChangeDirection = Standard_True; | |
1158 | else | |
1159 | { | |
1160 | Done = Standard_False; | |
1161 | if (!theStopOnDivergent || !myIsDivergent) | |
1162 | { | |
1163 | Done = Standard_True; | |
1164 | State = F.GetStateNumber(); | |
1165 | } | |
1166 | return; // si un gain inf a 5% on sort | |
1167 | } | |
1168 | } | |
1169 | else ChangeDirection = Standard_False; //Si oui on recommence | |
7fd59977 | 1170 | } |
1171 | else ChangeDirection = Standard_False; //Pas d'historique on continue | |
1172 | // Si le gradient ne diminue pas suffisemment par newton on essaie | |
1173 | // le gradient sauf si f diminue (aussi bizarre que cela puisse | |
1174 | // paraitre avec NEWTON le gradient de f peut augmenter alors que f | |
1175 | // diminue: dans ce cas il faut garder NEWTON) | |
1176 | if ((Gnr1 > 0.9*Oldgr) && | |
89d8607f | 1177 | (F2 > 0.5*PreviousMinimum)) { |
1178 | ChangeDirection = Standard_True; | |
7fd59977 | 1179 | } |
1180 | ||
1181 | // Si l'on ne decide pas de changer de strategie, on verifie, | |
1182 | // si ce n'est dejas fait | |
1183 | if ((!ChangeDirection) && (!Verif)) { | |
89d8607f | 1184 | for(i = Delta.Lower(); i <= Delta.Upper(); i++) { |
1185 | Delta(i) = PreviousSolution(i) - Sol(i); | |
1186 | } | |
1187 | if (IsSolutionReached(F)) { | |
1188 | Done = Standard_False; | |
1189 | if (!theStopOnDivergent || !myIsDivergent) | |
1190 | { | |
1191 | Done = Standard_True; | |
b659a6dc | 1192 | ////modified by jgv, 31.08.2011//// |
1193 | F.Value(Sol, FF); //update F before GetStateNumber | |
1194 | /////////////////////////////////// | |
89d8607f | 1195 | State = F.GetStateNumber(); |
1196 | } | |
1197 | return; | |
1198 | } | |
7fd59977 | 1199 | } |
1200 | } | |
1201 | else { // Cas de regression | |
1202 | if (!ChangeDirection) { // On passe au gradient | |
89d8607f | 1203 | ChangeDirection = Standard_True; |
1204 | Sol = PreviousSolution; | |
1205 | // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1); | |
1206 | if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) { | |
1207 | Done = Standard_False; | |
1208 | if (!theStopOnDivergent || !myIsDivergent) | |
1209 | { | |
1210 | State = F.GetStateNumber(); | |
1211 | } | |
1212 | return; | |
1213 | } | |
7fd59977 | 1214 | } |
5368adff | 1215 | else |
1216 | { | |
89d8607f | 1217 | |
b659a6dc | 1218 | if (!theStopOnDivergent || !myIsDivergent) |
1219 | { | |
1220 | State = F.GetStateNumber(); | |
1221 | } | |
5368adff | 1222 | return; // y a plus d'issues |
1223 | } | |
7fd59977 | 1224 | } |
1225 | } | |
b659a6dc | 1226 | if (!theStopOnDivergent || !myIsDivergent) |
1227 | { | |
1228 | State = F.GetStateNumber(); | |
1229 | } | |
7fd59977 | 1230 | } |
1231 | ||
859a47c3 | 1232 | //======================================================================= |
1233 | //function : Dump | |
1234 | //purpose : | |
1235 | //======================================================================= | |
1236 | void math_FunctionSetRoot::Dump(Standard_OStream& o) const | |
1237 | { | |
1238 | o << " math_FunctionSetRoot"; | |
89d8607f | 1239 | if (Done) { |
1240 | o << " Status = Done\n"; | |
1241 | o << " Location value = " << Sol << "\n"; | |
1242 | o << " Number of iterations = " << Kount << "\n"; | |
1243 | } | |
1244 | else { | |
859a47c3 | 1245 | o << "Status = Not Done\n"; |
89d8607f | 1246 | } |
7fd59977 | 1247 | } |
1248 | ||
859a47c3 | 1249 | //======================================================================= |
1250 | //function : Root | |
1251 | //purpose : | |
1252 | //======================================================================= | |
1253 | void math_FunctionSetRoot::Root(math_Vector& Root) const | |
1254 | { | |
7fd59977 | 1255 | StdFail_NotDone_Raise_if(!Done, " "); |
1256 | Standard_DimensionError_Raise_if(Root.Length() != Sol.Length(), " "); | |
1257 | Root = Sol; | |
1258 | } | |
1259 | ||
859a47c3 | 1260 | //======================================================================= |
1261 | //function : FunctionSetErrors | |
1262 | //purpose : | |
1263 | //======================================================================= | |
1264 | void math_FunctionSetRoot::FunctionSetErrors(math_Vector& Err) const | |
1265 | { | |
7fd59977 | 1266 | StdFail_NotDone_Raise_if(!Done, " "); |
1267 | Standard_DimensionError_Raise_if(Err.Length() != Sol.Length(), " "); | |
1268 | Err = Delta; | |
1269 | } |