1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
4 // The content of this file is subject to the Open CASCADE Technology Public
5 // License Version 6.5 (the "License"). You may not use the content of this file
6 // except in compliance with the License. Please obtain a copy of the License
7 // at http://www.opencascade.org and read it completely before using this file.
9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
12 // The Original Code and all software distributed under the License is
13 // distributed on an "AS IS" basis, without warranty of any kind, and the
14 // Initial Developer hereby disclaims all such warranties, including without
15 // limitation, any warranties of merchantability, fitness for a particular
16 // purpose or non-infringement. Please see the License for the specific terms
17 // and conditions governing the rights and limitations under the License.
19 // pmn 15/05/97 pas de Gauss avec des pivot trop petit. SVD fait mieux
20 // l'affaire + limitation de la longeur du pas + qq comentaire issus d'EUCLID3
21 // pmn 10/06/97 refonte totale du traitement des bords + ajustement des init
22 // et des tolerances pour brent...
25 #define No_Standard_RangeError
26 #define No_Standard_OutOfRange
27 #define No_Standard_DimensionError
30 //math_FunctionSetRoot.cxx
33 #include <math_FunctionSetRoot.ixx>
34 #include <Standard_DimensionError.hxx>
35 #include <math_Gauss.hxx>
36 #include <math_SVD.hxx>
37 #include <math_GaussLeastSquare.hxx>
38 #include <math_IntegerVector.hxx>
39 #include <math_Function.hxx>
40 #include <math_BrentMinimum.hxx>
41 #include <math_FunctionSetWithDerivatives.hxx>
42 #include <Precision.hxx>
45 //===========================================================================
46 // - A partir d une solution de depart, recherche d une direction.( Newton la
47 // plupart du temps, gradient si Newton echoue.
48 // - Recadrage au niveau des bornes avec recalcul de la direction si une
49 // inconnue a une valeur imposee.
50 // -Si On ne sort pas des bornes
51 // Tant que (On ne progresse pas assez ou on ne change pas de direction)
52 // . Si (Progression encore possible)
53 // Si (On ne sort pas des bornes)
54 // On essaye de progresser dans cette meme direction.
56 // On diminue le pas d'avancement ou on change de direction.
58 // Si on depasse le minimum
59 // On fait une interpolation parabolique.
60 // - Si on a progresse sur F
61 // On fait les tests d'arret
63 // On change de direction
64 //============================================================================
66 #define FSR_DEBUG(arg)
67 // Uncomment the following code to have debug output to cout
69 static Standard_Boolean mydebug = Standard_False;
71 #define FSR_DEBUG(arg) {if (mydebug) { cout << arg << endl; }}
74 class MyDirFunction : public math_Function
81 math_FunctionSetWithDerivatives *F;
85 MyDirFunction(math_Vector& V1,
89 math_FunctionSetWithDerivatives& f) ;
91 void Initialize(const math_Vector& p0, const math_Vector& dir) const;
93 Standard_Boolean Value(const math_Vector& Sol, math_Vector& FF,
94 math_Matrix& DF, math_Vector& GH,
95 Standard_Real& F2, Standard_Real& Gnr1);
96 // Standard_Boolean MyDirFunction::Value(const math_Vector& Sol, math_Vector& FF,
97 // math_Matrix& DF, math_Vector& GH,
98 // Standard_Real& F2, Standard_Real& Gnr1);
99 Standard_Boolean Value(const Standard_Real x, Standard_Real& fval) ;
103 MyDirFunction::MyDirFunction(math_Vector& V1,
107 math_FunctionSetWithDerivatives& f) {
116 void MyDirFunction::Initialize(const math_Vector& p0,
117 const math_Vector& dir) const
124 Standard_Boolean MyDirFunction::Value(const Standard_Real x,
128 for(Standard_Integer i = P->Lower(); i <= P->Upper(); i++) {
130 P->Value(i) = p * x + P0->Value(i);
132 if( F->Value(*P, *FV) ) {
134 Standard_Real aVal = 0.;
136 for(Standard_Integer i = FV->Lower(); i <= FV->Upper(); i++) {
139 if(aVal <= -1.e+100) // Precision::HalfInfinite() later
140 // if(Precision::IsInfinite(Abs(FV->Value(i)))) {
141 // fval = Precision::Infinite();
142 return Standard_False;
144 else if(aVal >= 1.e+100) // Precision::HalfInfinite() later
145 return Standard_False;
148 fval = 0.5 * (FV->Norm2());
149 return Standard_True;
151 return Standard_False;
154 Standard_Boolean MyDirFunction::Value(const math_Vector& Sol,
161 if( F->Values(Sol, FF, DF) ) {
163 Standard_Real aVal = 0.;
165 for(Standard_Integer i = FF.Lower(); i <= FF.Upper(); i++) {
166 // modified by NIZHNY-MKK Mon Oct 3 17:56:50 2005.BEGIN
169 if(aVal <= -1.e+100) // Precision::HalfInfinite() later
170 // if(Precision::IsInfinite(Abs(FF.Value(i)))) {
171 // F2 = Precision::Infinite();
172 // Gnr1 = Precision::Infinite();
173 return Standard_False;
175 else if(aVal >= 1.e+100) // Precision::HalfInfinite() later
176 return Standard_False;
177 // modified by NIZHNY-MKK Mon Oct 3 17:57:05 2005.END
181 F2 = 0.5 * (FF.Norm2());
182 GH.TMultiply(DF, FF);
183 for(Standard_Integer i = GH.Lower(); i <= GH.Upper(); i++)
185 if(Precision::IsInfinite((GH.Value(i))))
187 return Standard_False;
191 return Standard_True;
193 return Standard_False;
197 //--------------------------------------------------------------
198 static Standard_Boolean MinimizeDirection(const math_Vector& P0,
199 const math_Vector& P1,
200 const math_Vector& P2,
201 const Standard_Real F1,
203 const math_Vector& Tol,
205 // Purpose : minimisation a partir de 3 points
206 //-------------------------------------------------------
208 // (1) Evaluation d'un tolerance parametrique 1D
209 Standard_Real tol1d = 2.1 , invnorme, tsol;
210 Standard_Real Eps = 1.e-16;
211 Standard_Real ax, bx, cx;
213 for (Standard_Integer ii =1; ii<=Tol.Length(); ii++) {
214 invnorme = Abs(Delta(ii));
215 if (invnorme>Eps) tol1d = Min(tol1d, Tol(ii) / invnorme);
217 if (tol1d > 1.9) return Standard_False; //Pas la peine de se fatiguer
221 math_Vector PP0 = P0 ;
222 math_Vector PP1 = P1 ;
225 invnorme = Delta.Norm();
226 if (invnorme <= Eps) return Standard_False;
227 invnorme = ((Standard_Real) 1) / invnorme;
229 F.Initialize(P1, Delta);
232 FSR_DEBUG (" minimisation dans la direction")
234 cx = (P2-P1).Norm()*invnorme;
235 if (cx < 1.e-2) return Standard_False;
236 math_BrentMinimum Sol(F, ax, bx, cx, tol1d, 100, tol1d);
238 tsol = Sol.Location();
239 if (Sol.Minimum() < F1) {
240 Delta.Multiply(tsol);
241 return Standard_True;
244 return Standard_False;
247 //----------------------------------------------------------------------
248 static Standard_Boolean MinimizeDirection(const math_Vector& P,
250 const Standard_Real& PValue,
251 const Standard_Real& PDirValue,
252 const math_Vector& Gradient,
253 const math_Vector& DGradient,
254 const math_Vector& Tol,
256 // Purpose: minimisation a partir de 2 points et une derives
257 //----------------------------------------------------------------------
260 // (0) Evaluation d'un tolerance parametrique 1D
261 Standard_Boolean good = Standard_False;
262 Standard_Real Eps = 1.e-20;
263 Standard_Real tol1d = 1.1, Result = PValue, absdir;
265 for (Standard_Integer ii =1; ii<=Tol.Length(); ii++) {
266 absdir = Abs(Dir(ii));
267 if (absdir >Eps) tol1d = Min(tol1d, Tol(ii) / absdir);
269 if (tol1d > 0.9) return Standard_False;
271 // (1) On realise une premiere interpolation quadratique
272 Standard_Real ax, bx, cx, df1, df2, Delta, tsol, fsol, tsolbis;
273 FSR_DEBUG(" essai d interpolation")
278 if (df1<-Eps && df2>Eps) { // cuvette
279 tsol = - df1 / (df2 - df1);
284 ax = PDirValue - (bx+cx);
286 if (Abs(ax) <= Eps) { // cas lineaire
287 if ((Abs(bx) >= Eps)) tsol = - cx/bx;
290 else { // cas quadratique
291 Delta = bx*bx - 4*ax*cx;
293 // il y a des racines, on prend la plus proche de 0
295 tsol = -(bx + Delta);
296 tsolbis = (Delta - bx);
297 if (Abs(tsolbis) < Abs(tsol)) tsol = tsolbis;
301 // pas ou peu de racine : on "extremise"
307 if (Abs(tsol) >= 1) return Standard_False; //resultat sans interet
309 F.Initialize(P, Dir);
313 good = Standard_True;
315 FSR_DEBUG("t= "<<tsol<<" F = " << fsol << " OldF = "<<PValue)
318 // (2) Si l'on a pas assez progresser on realise une recherche
319 // en bonne et due forme, a partir des inits precedents
320 if ((fsol > 0.2*PValue) && (tol1d < 0.5)) {
323 ax = tsol; bx = 0.0; cx = 1.0;
326 ax = 0.0; bx = tsol; cx = 1.0;
328 FSR_DEBUG(" minimisation dans la direction")
329 math_BrentMinimum Sol(F, ax, bx, cx, tol1d, 100, tol1d);
331 if (Sol.Minimum() <= Result) {
332 tsol = Sol.Location();
333 good = Standard_True;
334 FSR_DEBUG("t= "<<tsol<<" F ="<< Sol.Minimum() << " OldF = "<<Result)
339 // mise a jour du Delta
345 //------------------------------------------------------
346 static void SearchDirection(const math_Matrix& DF,
347 const math_Vector& GH,
348 const math_Vector& FF,
349 Standard_Boolean ChangeDirection,
350 const math_Vector& InvLengthMax,
351 math_Vector& Direction,
355 Standard_Integer Ninc = DF.ColNumber(), Neq = DF.RowNumber();
356 Standard_Real Eps = 1.e-32;
357 if (!ChangeDirection) {
359 for (Standard_Integer i = FF.Lower(); i <= FF.Upper(); i++) {
360 Direction(i) = -FF(i);
362 math_Gauss Solut(DF, 1.e-9);
363 if (Solut.IsDone()) Solut.Solve(Direction);
364 else { // we have to "forget" singular directions.
365 FSR_DEBUG(" Matrice singuliere : On prend SVD")
366 math_SVD SolvebySVD(DF);
367 if (SolvebySVD.IsDone()) SolvebySVD.Solve(-1*FF, Direction);
368 else ChangeDirection = Standard_True;
371 else if (Ninc > Neq) {
373 if (Solut.IsDone()) Solut.Solve(-1*FF, Direction);
374 else ChangeDirection = Standard_True;
376 else if (Ninc < Neq) { // Calcul par GaussLeastSquare
377 math_GaussLeastSquare Solut(DF);
378 if (Solut.IsDone()) Solut.Solve(-1*FF, Direction);
379 else ChangeDirection = Standard_True;
382 // Il vaut mieux interdire des directions trops longue
383 // Afin de blinder les cas trop mal conditionne
384 // PMN 12/05/97 Traitement des singularite dans les conges
385 // Sur des surfaces periodiques
387 Standard_Real ratio = Abs( Direction(Direction.Lower())
388 *InvLengthMax(Direction.Lower()) );
390 for (i = Direction.Lower()+1; i<=Direction.Upper(); i++) {
391 ratio = Max(ratio, Abs( Direction(i)*InvLengthMax(i)) );
398 if (Dy >= -Eps) { // newton "ne descend pas" on prend le gradient
399 ChangeDirection = Standard_True;
401 if (ChangeDirection) { // On va faire un gradient !
402 for (i = Direction.Lower(); i <= Direction.Upper(); i++) {
403 Direction(i) = - GH(i);
410 //=====================================================================
411 static void SearchDirection(const math_Matrix& DF,
412 const math_Vector& GH,
413 const math_Vector& FF,
414 const math_IntegerVector& Constraints,
415 // const math_Vector& X, // Le point d'init
416 const math_Vector& , // Le point d'init
417 Standard_Boolean ChangeDirection,
418 const math_Vector& InvLengthMax,
419 math_Vector& Direction,
421 //Purpose : Recherche une direction (et un pas si Newton Fonctionne) le long
423 //=====================================================================
425 Standard_Integer Ninc = DF.ColNumber(), Neq = DF.RowNumber();
426 Standard_Integer i, j, k, Cons = 0;
428 // verification sur les bornes imposees:
430 for (i = 1; i <= Ninc; i++) {
431 if (Constraints(i) != 0) Cons++;
432 // sinon le systeme a resoudre ne change pas.
436 SearchDirection(DF, GH, FF, ChangeDirection, InvLengthMax,
439 else if (Cons == Ninc) { // il n'y a plus rien a faire...
440 for(Standard_Integer i = Direction.Lower(); i <= Direction.Upper(); i++) {
445 else { //(1) Cas general : On definit un sous probleme
446 math_Matrix DF2(1, Neq, 1, Ninc-Cons);
447 math_Vector MyGH(1, Ninc-Cons);
448 math_Vector MyDirection(1, Ninc-Cons);
449 math_Vector MyInvLengthMax(1, Ninc);
451 for (k=1, i = 1; i <= Ninc; i++) {
452 if (Constraints(i) == 0) {
454 MyInvLengthMax(k) = InvLengthMax(i);
455 MyDirection(k) = Direction(i);
456 for (j = 1; j <= Neq; j++) {
457 DF2(j, k) = DF(j, i);
459 k++; //on passe a l'inconnue suivante
463 SearchDirection(DF2, MyGH, FF, ChangeDirection, MyInvLengthMax,
466 // (3) On l'interprete...
467 // Reconstruction de Direction:
468 for (i = 1, k = 1; i <= Ninc; i++) {
469 if (Constraints(i) == 0) {
470 if (!ChangeDirection) {
471 Direction(i) = MyDirection(k);
473 else Direction(i) = - GH(i);
485 //====================================================
486 Standard_Boolean Bounds(const math_Vector& InfBound,
487 const math_Vector& SupBound,
488 const math_Vector& Tol,
490 const math_Vector& SolSave,
491 math_IntegerVector& Constraints,
493 Standard_Boolean& theIsNewSol)
495 // Purpose: Trims an initial solution Sol to be within a domain defined by
496 // InfBound and SupBound. Delta will contain a distance between final Sol and
498 // IsNewSol returns False, if final Sol fully coincides with SolSave, i.e.
499 // if SolSave already lied on a boundary and initial Sol was fully beyond it
500 //======================================================
502 Standard_Boolean Out = Standard_False;
503 Standard_Integer i, Ninc = Sol.Length();
504 Standard_Real monratio = 1;
506 theIsNewSol = Standard_True;
508 // Calcul du ratio de recadrage
509 for (i = 1; i <= Ninc; i++) {
511 Delta(i) = Sol(i) - SolSave(i);
512 if (InfBound(i) == SupBound(i)) {
514 Out = Standard_True; // Ok mais, cela devrait etre eviter
516 else if(Sol(i) < InfBound(i)) {
519 // Delta(i) is negative
520 if (-Delta(i) > Tol(i)) // Afin d'eviter des ratio nulles pour rien
521 monratio = Min(monratio, (InfBound(i) - SolSave(i))/Delta(i) );
523 else if (Sol(i) > SupBound(i)) {
526 // Delta(i) is positive
527 if (Delta(i) > Tol(i))
528 monratio = Min(monratio, (SupBound(i) - SolSave(i))/Delta(i) );
532 if (Out){ // Troncature et derniers recadrage pour blinder (pb numeriques)
533 if (monratio == 0.0) {
534 theIsNewSol = Standard_False;
540 for (i = 1; i <= Ninc; i++) {
541 if(Sol(i) < InfBound(i)) {
542 Sol(i) = InfBound(i);
543 Delta(i) = Sol(i) - SolSave(i);
545 else if (Sol(i) > SupBound(i)) {
546 Sol(i) = SupBound(i);
547 Delta(i) = Sol(i) - SolSave(i);
559 math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
560 const math_Vector& Tolerance,
561 const Standard_Integer NbIterations) :
562 Delta(1, F.NbVariables()),
563 Sol(1, F.NbVariables()),
564 DF(1, F.NbEquations(),
566 Tol(1,F.NbVariables()),
568 InfBound(1, F.NbVariables()),
569 SupBound(1, F.NbVariables()),
571 SolSave(1, F.NbVariables()),
572 GH(1, F.NbVariables()),
573 DH(1, F.NbVariables()),
574 DHSave(1, F.NbVariables()),
575 FF(1, F.NbEquations()),
576 PreviousSolution(1, F.NbVariables()),
577 Save(0, NbIterations),
578 Constraints(1, F.NbVariables()),
579 Temp1(1, F.NbVariables()),
580 Temp2(1, F.NbVariables()),
581 Temp3(1, F.NbVariables()),
582 Temp4(1, F.NbEquations())
585 for (Standard_Integer i = 1; i <= Tol.Length(); i++) {
586 Tol(i) =Tolerance(i);
588 Itermax = NbIterations;
591 math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
592 const Standard_Integer NbIterations) :
593 Delta(1, F.NbVariables()),
594 Sol(1, F.NbVariables()),
595 DF(1, F.NbEquations(),
597 Tol(1, F.NbVariables()),
599 InfBound(1, F.NbVariables()),
600 SupBound(1, F.NbVariables()),
602 SolSave(1, F.NbVariables()),
603 GH(1, F.NbVariables()),
604 DH(1, F.NbVariables()),
605 DHSave(1, F.NbVariables()),
606 FF(1, F.NbEquations()),
607 PreviousSolution(1, F.NbVariables()),
608 Save(0, NbIterations),
609 Constraints(1, F.NbVariables()),
610 Temp1(1, F.NbVariables()),
611 Temp2(1, F.NbVariables()),
612 Temp3(1, F.NbVariables()),
613 Temp4(1, F.NbEquations())
616 Itermax = NbIterations;
621 math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
622 const math_Vector& StartingPoint,
623 const math_Vector& Tolerance,
624 const math_Vector& infBound,
625 const math_Vector& supBound,
626 const Standard_Integer NbIterations,
627 Standard_Boolean theStopOnDivergent) :
628 Delta(1, F.NbVariables()),
629 Sol(1, F.NbVariables()),
630 DF(1, F.NbEquations(),
632 Tol(1,F.NbVariables()),
635 InfBound(1, F.NbVariables()),
636 SupBound(1, F.NbVariables()),
638 SolSave(1, F.NbVariables()),
639 GH(1, F.NbVariables()),
640 DH(1, F.NbVariables()),
641 DHSave(1, F.NbVariables()),
642 FF(1, F.NbEquations()),
643 PreviousSolution(1, F.NbVariables()),
644 Save(0, NbIterations),
645 Constraints(1, F.NbVariables()),
646 Temp1(1, F.NbVariables()),
647 Temp2(1, F.NbVariables()),
648 Temp3(1, F.NbVariables()),
649 Temp4(1, F.NbEquations())
653 for (Standard_Integer i = 1; i <= Tol.Length(); i++) {
654 Tol(i) =Tolerance(i);
656 Itermax = NbIterations;
657 Perform(F, StartingPoint, infBound, supBound, theStopOnDivergent);
660 math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
661 const math_Vector& StartingPoint,
662 const math_Vector& Tolerance,
663 const Standard_Integer NbIterations) :
664 Delta(1, F.NbVariables()),
665 Sol(1, F.NbVariables()),
666 DF(1, F.NbEquations(),
667 1, StartingPoint.Length()),
668 Tol(1,F.NbVariables()),
670 InfBound(1, F.NbVariables()),
671 SupBound(1, F.NbVariables()),
673 SolSave(1, F.NbVariables()),
674 GH(1, F.NbVariables()),
675 DH(1, F.NbVariables()),
676 DHSave(1, F.NbVariables()),
677 FF(1, F.NbEquations()),
678 PreviousSolution(1, F.NbVariables()),
679 Save(0, NbIterations),
680 Constraints(1, F.NbVariables()),
681 Temp1(1, F.NbVariables()),
682 Temp2(1, F.NbVariables()),
683 Temp3(1, F.NbVariables()),
684 Temp4(1, F.NbEquations())
687 for (Standard_Integer i = 1; i <= Tol.Length(); i++) {
688 Tol(i) = Tolerance(i);
690 Itermax = NbIterations;
691 InfBound.Init(RealFirst());
692 SupBound.Init(RealLast());
693 Perform(F, StartingPoint, InfBound, SupBound);
696 void math_FunctionSetRoot::Delete()
699 void math_FunctionSetRoot::SetTolerance(const math_Vector& Tolerance)
701 for (Standard_Integer i = 1; i <= Tol.Length(); i++) {
702 Tol(i) = Tolerance(i);
706 void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
707 const math_Vector& StartingPoint,
708 const math_Vector& InfBound,
709 const math_Vector& SupBound,
710 Standard_Boolean theStopOnDivergent)
712 Standard_Integer Ninc = F.NbVariables(), Neq = F.NbEquations();
715 (StartingPoint.Length()!= Ninc) ||
716 (InfBound.Length() != Ninc) ||
717 (SupBound.Length() != Ninc)) { Standard_DimensionError:: Raise(); }
720 Standard_Boolean ChangeDirection = Standard_False, Sort = Standard_False, isNewSol = Standard_False;
721 Standard_Boolean Good, Verif;
722 Standard_Boolean Stop;
723 const Standard_Real EpsSqrt = 1.e-16, Eps = 1.e-32, Eps2 = 1.e-64, Progres = 0.005;
724 Standard_Real F2, PreviousMinimum, Dy, OldF;
725 Standard_Real Ambda, Ambda2, Gnr1, Oldgr;
726 math_Vector InvLengthMax(1, Ninc); // Pour bloquer les pas a 1/4 du domaine
727 math_IntegerVector Constraints(1, Ninc); // Pour savoir sur quels bord on se trouve
728 for (i = 1; i <= Ninc ; i++) {
729 // modified by NIZHNY-MKK Mon Oct 3 18:03:50 2005
730 // InvLengthMax(i) = 1. / Max(Abs(SupBound(i) - InfBound(i))/4, 1.e-9);
731 InvLengthMax(i) = 1. / Max((SupBound(i) - InfBound(i))/4, 1.e-9);
734 MyDirFunction F_Dir(Temp1, Temp2, Temp3, Temp4, F);
735 Standard_Integer DescenteIter;
737 Done = Standard_False;
742 myIsDivergent = Standard_False;
743 for (i = 1; i <= Ninc; i++)
745 myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
747 if (theStopOnDivergent & myIsDivergent)
751 // Verification de la validite des inconnues par rapport aux bornes.
752 // Recentrage sur les bornes si pas valide.
753 for ( i = 1; i <= Ninc; i++) {
754 if (Sol(i) <= InfBound(i)) Sol(i) = InfBound(i);
755 else if (Sol(i) > SupBound(i)) Sol(i) = SupBound(i);
758 // Calcul de la premiere valeur de F et de son gradient
759 if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
760 Done = Standard_False;
761 if (!theStopOnDivergent || !myIsDivergent)
763 State = F.GetStateNumber();
768 // Le rang 0 de Save ne doit servir q'au test accelarteur en fin de boucle
769 // s'il on est dejas sur la solution, il faut leurer ce test pour eviter
770 // de faire une seconde iteration...
771 Save(0) = Max (F2, EpsSqrt);
772 Standard_Real aTol_Func = Epsilon(F2);
773 FSR_DEBUG("=== Mode Debug de Function Set Root" << endl)
774 FSR_DEBUG(" F2 Initial = " << F2)
776 if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
777 Done = Standard_False;
778 if (!theStopOnDivergent || !myIsDivergent)
780 Done = Standard_True;
781 State = F.GetStateNumber();
786 for (Kount = 1; Kount <= Itermax; Kount++) {
787 PreviousMinimum = F2;
789 PreviousSolution = Sol;
792 SearchDirection(DF, GH, FF, ChangeDirection, InvLengthMax, DH, Dy);
793 if (Abs(Dy) <= Eps) {
794 Done = Standard_False;
795 if (!theStopOnDivergent || !myIsDivergent)
797 Done = Standard_True;
798 ////modified by jgv, 31.08.2011////
799 F.Value(Sol, FF); //update F before GetStateNumber
800 ///////////////////////////////////
801 State = F.GetStateNumber();
805 if (ChangeDirection) {
806 Ambda = Ambda2 / Sqrt(Abs(Dy));
807 if (Ambda > 1.0) Ambda = 1.0;
811 Ambda2 = 0.5*Ambda/DH.Norm();
814 for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
815 Sol(i) = Sol(i) + Ambda * DH(i);
818 for (i = 1; i <= Ninc; i++)
820 myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
822 if (theStopOnDivergent & myIsDivergent)
827 Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
828 Constraints, Delta, isNewSol);
833 // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
834 if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
835 Done = Standard_False;
836 if (!theStopOnDivergent || !myIsDivergent)
838 State = F.GetStateNumber();
844 FSR_DEBUG("Kount = " << Kount)
845 FSR_DEBUG("Le premier F2 = " << F2)
846 FSR_DEBUG("Direction = " << ChangeDirection)
848 if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
849 Done = Standard_False;
850 if (!theStopOnDivergent || !myIsDivergent)
852 Done = Standard_True;
853 ////modified by jgv, 31.08.2011////
854 F.Value(Sol, FF); //update F before GetStateNumber
855 ///////////////////////////////////
856 State = F.GetStateNumber();
861 if (Sort || (F2/PreviousMinimum > Progres)) {
863 OldF = PreviousMinimum;
864 Stop = Standard_False;
865 Good = Standard_False;
867 Standard_Boolean Sortbis;
869 // -------------------------------------------------
870 // Traitement standard sans traitement des bords
871 // -------------------------------------------------
872 if (!Sort) { // si l'on est pas sortie on essaye de progresser en avant
873 while((F2/PreviousMinimum > Progres) && !Stop) {
874 if (F2 < OldF && (Dy < 0.0)) {
875 // On essaye de progresser dans cette direction.
876 FSR_DEBUG(" iteration de descente = " << DescenteIter)
880 for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
881 Sol(i) = Sol(i) + Ambda * DH(i);
884 for (i = 1; i <= Ninc; i++)
886 myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
888 if (theStopOnDivergent & myIsDivergent)
893 Stop = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
894 Constraints, Delta, isNewSol);
895 FSR_DEBUG(" Augmentation de lambda")
899 if ((F2 >= OldF)||(F2 >= PreviousMinimum)) {
900 Good = Standard_False;
901 if (DescenteIter == 0) {
902 // C'est le premier pas qui flanche, on fait une interpolation.
903 // et on minimise si necessaire.
905 Good = MinimizeDirection(SolSave, Delta, OldF, F2, DHSave, GH,
908 else if (ChangeDirection || (DescenteIter>1)
909 || (OldF>PreviousMinimum) ) {
910 // La progression a ete utile, on minimise...
912 Good = MinimizeDirection(PreviousSolution, SolSave, Sol,
913 OldF, Delta, Tol, F_Dir);
922 for (i = 1; i <= Ninc; i++)
924 myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
926 if (theStopOnDivergent & myIsDivergent)
931 Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
932 Constraints, Delta, isNewSol);
934 Sort = Standard_False; // On a rejete le point sur la frontiere
936 Stop = Standard_True; // et on sort dans tous les cas...
940 // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
941 if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
942 Done = Standard_False;
943 if (!theStopOnDivergent || !myIsDivergent)
945 State = F.GetStateNumber();
951 if (Abs(Dy) <= Eps) {
954 Done = Standard_False;
955 if (!theStopOnDivergent || !myIsDivergent)
957 Done = Standard_True;
958 ////modified by jgv, 31.08.2011////
959 F.Value(Sol, FF); //update F before GetStateNumber
960 ///////////////////////////////////
961 State = F.GetStateNumber();
965 if (DescenteIter >= 10) {
966 Stop = Standard_True;
969 FSR_DEBUG("--- Sortie du Traitement Standard")
970 FSR_DEBUG(" DescenteIter = "<<DescenteIter << " F2 = " << F2)
972 // ------------------------------------
973 // on passe au traitement des bords
974 // ------------------------------------
976 Stop = (F2 > 1.001*OldF); // Pour ne pas progresser sur le bord
979 while (Sortbis && ((F2<OldF)|| (DescenteIter == 0))
982 // On essaye de progresser sur le bord
985 SearchDirection(DF, GH, FF, Constraints, Sol,
986 ChangeDirection, InvLengthMax, DH, Dy);
987 FSR_DEBUG(" Conditional Direction = " << ChangeDirection)
988 if (Dy<-Eps) { //Pour eviter des calculs inutiles et des /0...
989 if (ChangeDirection) {
991 // Ambda = Ambda2 / Sqrt(Abs(Dy));
992 Ambda = Ambda2 / Sqrt(-Dy);
993 if (Ambda > 1.0) Ambda = 1.0;
997 Ambda2 = 0.5*Ambda/DH.Norm();
1000 for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
1001 Sol(i) = Sol(i) + Ambda * DH(i);
1004 for (i = 1; i <= Ninc; i++)
1006 myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
1008 if (theStopOnDivergent & myIsDivergent)
1013 Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
1014 Constraints, Delta, isNewSol);
1018 // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
1019 if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
1020 Done = Standard_False;
1021 if (!theStopOnDivergent || !myIsDivergent)
1023 State = F.GetStateNumber();
1029 FSR_DEBUG("--- Iteration au bords : " << DescenteIter)
1030 FSR_DEBUG("--- F2 = " << F2)
1033 Stop = Standard_True;
1036 while((F2/PreviousMinimum > Progres) && (F2<OldF) && (!Stop) ) {
1038 FSR_DEBUG("--- Iteration de descente conditionnel = " << DescenteIter)
1039 if (F2 < OldF && Dy < 0.0) {
1040 // On essaye de progresser dans cette direction.
1043 for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
1044 Sol(i) = Sol(i) + Ambda * DH(i);
1047 for (i = 1; i <= Ninc; i++)
1049 myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
1051 if (theStopOnDivergent & myIsDivergent)
1056 Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
1057 Constraints, Delta, isNewSol);
1061 // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
1062 if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
1063 Done = Standard_False;
1064 if (!theStopOnDivergent || !myIsDivergent)
1066 State = F.GetStateNumber();
1073 Stop = ((Dy >=0) || (DescenteIter >= 10) || Sortbis);
1075 Stop = ((Dy >=0) || (DescenteIter >= 10));
1077 if (((F2/PreviousMinimum > Progres) &&
1078 (F2>=OldF))||(F2>=PreviousMinimum)) {
1079 // On minimise par Brent
1081 Good = MinimizeDirection(SolSave, Delta, OldF, F2,
1082 DHSave, GH, Tol, F_Dir);
1085 Sort = Standard_False;
1088 Sol = SolSave + Delta;
1090 for (i = 1; i <= Ninc; i++)
1092 myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
1094 if (theStopOnDivergent & myIsDivergent)
1099 Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
1100 Constraints, Delta, isNewSol);
1102 // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
1103 if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
1104 Done = Standard_False;
1105 if (!theStopOnDivergent || !myIsDivergent)
1107 State = F.GetStateNumber();
1115 FSR_DEBUG("--- Sortie du Traitement des Bords")
1116 FSR_DEBUG("--- DescenteIter = "<<DescenteIter << " F2 = " << F2)
1120 // ---------------------------------------------
1121 // on passe aux tests d'ARRET
1122 // ---------------------------------------------
1124 // Est ce la solution ?
1125 if (ChangeDirection) Verif = Standard_True;
1126 // Gradient : Il faut eviter de boucler
1128 Verif = Standard_False;
1129 if (Kount > 1) { // Pour accelerer les cas quasi-quadratique
1130 if (Save(Kount-1)<1.e-4*Save(Kount-2)) Verif = Standard_True;
1132 else Verif = (F2 < 1.e-6*Save(0)); //Pour les cas dejas solutions
1135 for(i = Delta.Lower(); i <= Delta.Upper(); i++) {
1136 Delta(i) = PreviousSolution(i) - Sol(i);
1139 if (IsSolutionReached(F)) {
1140 if (PreviousMinimum < F2) {
1143 Done = Standard_False;
1144 if (!theStopOnDivergent || !myIsDivergent)
1146 Done = Standard_True;
1147 ////modified by jgv, 31.08.2011////
1148 F.Value(Sol, FF); //update F before GetStateNumber
1149 ///////////////////////////////////
1150 State = F.GetStateNumber();
1155 //fin du test solution
1157 // Analyse de la progression...
1158 //comparison of current minimum and previous minimum
1159 if ((F2 - PreviousMinimum) <= aTol_Func){
1161 // L'historique est il bon ?
1162 if (F2 >= 0.95*Save(Kount - 5)) {
1163 if (!ChangeDirection) ChangeDirection = Standard_True;
1166 Done = Standard_False;
1167 if (!theStopOnDivergent || !myIsDivergent)
1169 Done = Standard_True;
1170 State = F.GetStateNumber();
1172 return; // si un gain inf a 5% on sort
1175 else ChangeDirection = Standard_False; //Si oui on recommence
1177 else ChangeDirection = Standard_False; //Pas d'historique on continue
1178 // Si le gradient ne diminue pas suffisemment par newton on essaie
1179 // le gradient sauf si f diminue (aussi bizarre que cela puisse
1180 // paraitre avec NEWTON le gradient de f peut augmenter alors que f
1181 // diminue: dans ce cas il faut garder NEWTON)
1182 if ((Gnr1 > 0.9*Oldgr) &&
1183 (F2 > 0.5*PreviousMinimum)) {
1184 ChangeDirection = Standard_True;
1187 // Si l'on ne decide pas de changer de strategie, on verifie,
1188 // si ce n'est dejas fait
1189 if ((!ChangeDirection) && (!Verif)) {
1190 for(i = Delta.Lower(); i <= Delta.Upper(); i++) {
1191 Delta(i) = PreviousSolution(i) - Sol(i);
1193 if (IsSolutionReached(F)) {
1194 Done = Standard_False;
1195 if (!theStopOnDivergent || !myIsDivergent)
1197 Done = Standard_True;
1198 ////modified by jgv, 31.08.2011////
1199 F.Value(Sol, FF); //update F before GetStateNumber
1200 ///////////////////////////////////
1201 State = F.GetStateNumber();
1207 else { // Cas de regression
1208 if (!ChangeDirection) { // On passe au gradient
1209 ChangeDirection = Standard_True;
1210 Sol = PreviousSolution;
1211 // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
1212 if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
1213 Done = Standard_False;
1214 if (!theStopOnDivergent || !myIsDivergent)
1216 State = F.GetStateNumber();
1224 if (!theStopOnDivergent || !myIsDivergent)
1226 State = F.GetStateNumber();
1228 return; // y a plus d'issues
1232 if (!theStopOnDivergent || !myIsDivergent)
1234 State = F.GetStateNumber();
1241 Standard_Boolean math_FunctionSetRoot::IsSolutionReached(math_FunctionSetWithDerivatives& ) {
1242 for(Standard_Integer i = 1; i<= Sol.Length(); i++) {
1243 if(Abs(Delta(i)) > Tol(i)) {return Standard_False;}
1245 return Standard_True;
1249 void math_FunctionSetRoot::Dump(Standard_OStream& o) const {
1250 o<<" math_FunctionSetRoot";
1252 o << " Status = Done\n";
1253 o << " Location value = " << Sol << "\n";
1254 o << " Number of iterations = " << Kount << "\n";
1257 o<<"Status = Not Done\n";
1262 void math_FunctionSetRoot::Root(math_Vector& Root) const{
1263 StdFail_NotDone_Raise_if(!Done, " ");
1264 Standard_DimensionError_Raise_if(Root.Length() != Sol.Length(), " ");
1269 void math_FunctionSetRoot::FunctionSetErrors(math_Vector& Err) const{
1270 StdFail_NotDone_Raise_if(!Done, " ");
1271 Standard_DimensionError_Raise_if(Err.Length() != Sol.Length(), " ");