1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
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
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.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
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...
21 #define No_Standard_RangeError
22 #define No_Standard_OutOfRange
23 #define No_Standard_DimensionError
26 //math_FunctionSetRoot.cxx
28 #include <math_BrentMinimum.hxx>
29 #include <math_Function.hxx>
30 #include <math_FunctionSetRoot.hxx>
31 #include <math_FunctionSetWithDerivatives.hxx>
32 #include <math_Gauss.hxx>
33 #include <math_GaussLeastSquare.hxx>
34 #include <math_IntegerVector.hxx>
35 #include <math_Matrix.hxx>
36 #include <math_SVD.hxx>
37 #include <Precision.hxx>
38 #include <Standard_DimensionError.hxx>
39 #include <StdFail_NotDone.hxx>
41 //===========================================================================
42 // - A partir d une solution de depart, recherche d une direction.( Newton la
43 // plupart du temps, gradient si Newton echoue.
44 // - Recadrage au niveau des bornes avec recalcul de la direction si une
45 // inconnue a une valeur imposee.
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.
52 // On diminue le pas d'avancement ou on change de direction.
54 // Si on depasse le minimum
55 // On fait une interpolation parabolique.
56 // - Si on a progresse sur F
57 // On fait les tests d'arret
59 // On change de direction
60 //============================================================================
61 #define FSR_DEBUG(arg)
62 // Uncomment the following code to have debug output to cout
63 //==========================================================
64 //static Standard_Boolean mydebug = Standard_True;
66 //#define FSR_DEBUG(arg) {if (mydebug) { cout << arg << endl; }}
67 //===========================================================
69 class MyDirFunction : public math_Function
76 math_FunctionSetWithDerivatives *F;
80 MyDirFunction(math_Vector& V1,
84 math_FunctionSetWithDerivatives& f) ;
86 void Initialize(const math_Vector& p0, const math_Vector& dir) const;
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) ;
98 MyDirFunction::MyDirFunction(math_Vector& V1,
102 math_FunctionSetWithDerivatives& f) {
111 void MyDirFunction::Initialize(const math_Vector& p0,
112 const math_Vector& dir) const
119 Standard_Boolean MyDirFunction::Value(const Standard_Real x,
123 for(Standard_Integer i = P->Lower(); i <= P->Upper(); i++) {
125 P->Value(i) = p * x + P0->Value(i);
127 if( F->Value(*P, *FV) )
130 Standard_Real aVal = 0.0;
132 for(Standard_Integer i = FV->Lower(); i <= FV->Upper(); i++)
135 if(aVal <= -1.e+100) // Precision::HalfInfinite() later
136 return Standard_False;
137 else if(aVal >= 1.e+100) // Precision::HalfInfinite() later
138 return Standard_False;
141 fval = 0.5 * (FV->Norm2());
142 return Standard_True;
144 return Standard_False;
147 Standard_Boolean MyDirFunction::Value(const math_Vector& Sol,
154 if( F->Values(Sol, FF, DF) ) {
156 Standard_Real aVal = 0.;
158 for(Standard_Integer i = FF.Lower(); i <= FF.Upper(); i++) {
159 // modified by NIZHNY-MKK Mon Oct 3 17:56:50 2005.BEGIN
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;
168 else if(aVal >= 1.e+100) // Precision::HalfInfinite() later
169 return Standard_False;
170 // modified by NIZHNY-MKK Mon Oct 3 17:57:05 2005.END
174 F2 = 0.5 * (FF.Norm2());
175 GH.TMultiply(DF, FF);
176 for(Standard_Integer i = GH.Lower(); i <= GH.Upper(); i++)
178 if(Precision::IsInfinite((GH.Value(i))))
180 return Standard_False;
184 return Standard_True;
186 return Standard_False;
190 //--------------------------------------------------------------
191 static Standard_Boolean MinimizeDirection(const math_Vector& P0,
192 const math_Vector& P1,
193 const math_Vector& P2,
194 const Standard_Real F1,
196 const math_Vector& Tol,
198 // Purpose : minimisation a partir de 3 points
199 //-------------------------------------------------------
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;
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);
210 if (tol1d > 1.9) return Standard_False; //Pas la peine de se fatiguer
214 math_Vector PP0 = P0 ;
215 math_Vector PP1 = P1 ;
218 invnorme = Delta.Norm();
219 if (invnorme <= Eps) return Standard_False;
220 invnorme = ((Standard_Real) 1) / invnorme;
222 F.Initialize(P1, Delta);
225 FSR_DEBUG (" minimisation dans la direction")
227 cx = (P2-P1).Norm()*invnorme;
229 return Standard_False;
231 math_BrentMinimum Sol(tol1d, 100, tol1d);
232 Sol.Perform(F, ax, bx, cx);
235 tsol = Sol.Location();
236 if (Sol.Minimum() < F1) {
237 Delta.Multiply(tsol);
238 return Standard_True;
241 return Standard_False;
244 //----------------------------------------------------------------------
245 static Standard_Boolean MinimizeDirection(const math_Vector& P,
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,
253 // Purpose: minimisation a partir de 2 points et une derives
254 //----------------------------------------------------------------------
257 if(Precision::IsInfinite(PValue) || Precision::IsInfinite(PDirValue))
259 return Standard_False;
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;
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);
270 if (tol1d > 0.9) return Standard_False;
272 // (1) On realise une premiere interpolation quadratique
273 Standard_Real ax, bx, cx, df1, df2, Delta, tsol, fsol, tsolbis;
274 FSR_DEBUG(" essai d interpolation");
279 if (df1<-Eps && df2>Eps) { // cuvette
280 tsol = - df1 / (df2 - df1);
285 ax = PDirValue - (bx+cx);
287 if (Abs(ax) <= Eps) { // cas lineaire
288 if ((Abs(bx) >= Eps)) tsol = - cx/bx;
291 else { // cas quadratique
292 Delta = bx*bx - 4*ax*cx;
294 // il y a des racines, on prend la plus proche de 0
296 tsol = -(bx + Delta);
297 tsolbis = (Delta - bx);
298 if (Abs(tsolbis) < Abs(tsol)) tsol = tsolbis;
302 // pas ou peu de racine : on "extremise"
308 if (Abs(tsol) >= 1) return Standard_False; //resultat sans interet
310 F.Initialize(P, Dir);
314 good = Standard_True;
316 FSR_DEBUG("t= "<<tsol<<" F = " << fsol << " OldF = "<<PValue);
319 // (2) Si l'on a pas assez progresser on realise une recherche
320 // en bonne et due forme, a partir des inits precedents
321 if ((fsol > 0.2*PValue) && (tol1d < 0.5))
325 ax = tsol; bx = 0.0; cx = 1.0;
328 ax = 0.0; bx = tsol; cx = 1.0;
330 FSR_DEBUG(" minimisation dans la direction");
332 math_BrentMinimum Sol(tol1d, 100, tol1d);
335 Sol.Perform(F, ax, bx, cx);
338 if (Sol.Minimum() <= Result)
340 tsol = Sol.Location();
341 good = Standard_True;
342 Result = Sol.Minimum();
344 // Objective function changes too fast ->
345 // perform additional computations.
346 if (Gradient.Norm2() > 1.0 / Precision::SquareConfusion() &&
348 tsol < cx) // Solution inside of (ax, cx) interval.
350 // First and second part invocation.
351 Sol.Perform(F, ax, (ax + tsol) / 2.0, tsol);
354 if (Sol.Minimum() <= Result)
356 tsol = Sol.Location();
357 good = Standard_True;
358 Result = Sol.Minimum();
362 Sol.Perform(F, tsol, (cx + tsol) / 2.0, cx);
365 if (Sol.Minimum() <= Result)
367 tsol = Sol.Location();
368 good = Standard_True;
369 Result = Sol.Minimum();
373 } // if (Sol.Minimum() <= Result)
374 } // if(Sol.IsDone())
379 // mise a jour du Delta
385 //------------------------------------------------------
386 static void SearchDirection(const math_Matrix& DF,
387 const math_Vector& GH,
388 const math_Vector& FF,
389 Standard_Boolean ChangeDirection,
390 const math_Vector& InvLengthMax,
391 math_Vector& Direction,
395 Standard_Integer Ninc = DF.ColNumber(), Neq = DF.RowNumber();
396 Standard_Real Eps = 1.e-32;
397 if (!ChangeDirection) {
399 for (Standard_Integer i = FF.Lower(); i <= FF.Upper(); i++) {
400 Direction(i) = -FF(i);
402 math_Gauss Solut(DF, 1.e-9);
403 if (Solut.IsDone()) Solut.Solve(Direction);
404 else { // we have to "forget" singular directions.
405 FSR_DEBUG(" Matrice singuliere : On prend SVD");
406 math_SVD SolvebySVD(DF);
407 if (SolvebySVD.IsDone()) SolvebySVD.Solve(-1*FF, Direction);
408 else ChangeDirection = Standard_True;
411 else if (Ninc > Neq) {
413 if (Solut.IsDone()) Solut.Solve(-1*FF, Direction);
414 else ChangeDirection = Standard_True;
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;
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
427 Standard_Real ratio = Abs( Direction(Direction.Lower())
428 *InvLengthMax(Direction.Lower()) );
430 for (i = Direction.Lower()+1; i<=Direction.Upper(); i++) {
431 ratio = Max(ratio, Abs( Direction(i)*InvLengthMax(i)) );
438 if (Dy >= -Eps) { // newton "ne descend pas" on prend le gradient
439 ChangeDirection = Standard_True;
441 if (ChangeDirection) { // On va faire un gradient !
442 for (i = Direction.Lower(); i <= Direction.Upper(); i++) {
443 Direction(i) = - GH(i);
450 //=====================================================================
451 static void SearchDirection(const math_Matrix& DF,
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,
458 const math_Vector& InvLengthMax,
459 math_Vector& Direction,
461 //Purpose : Recherche une direction (et un pas si Newton Fonctionne) le long
463 //=====================================================================
465 Standard_Integer Ninc = DF.ColNumber(), Neq = DF.RowNumber();
466 Standard_Integer i, j, k, Cons = 0;
468 // verification sur les bornes imposees:
470 for (i = 1; i <= Ninc; i++) {
471 if (Constraints(i) != 0) Cons++;
472 // sinon le systeme a resoudre ne change pas.
476 SearchDirection(DF, GH, FF, ChangeDirection, InvLengthMax,
479 else if (Cons == Ninc) { // il n'y a plus rien a faire...
480 for(i = Direction.Lower(); i <= Direction.Upper(); i++) {
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);
491 for (k=1, i = 1; i <= Ninc; i++) {
492 if (Constraints(i) == 0) {
494 MyInvLengthMax(k) = InvLengthMax(i);
495 MyDirection(k) = Direction(i);
496 for (j = 1; j <= Neq; j++) {
497 DF2(j, k) = DF(j, i);
499 k++; //on passe a l'inconnue suivante
503 SearchDirection(DF2, MyGH, FF, ChangeDirection, MyInvLengthMax,
506 // (3) On l'interprete...
507 // Reconstruction de Direction:
508 for (i = 1, k = 1; i <= Ninc; i++) {
509 if (Constraints(i) == 0) {
510 if (!ChangeDirection) {
511 Direction(i) = MyDirection(k);
513 else Direction(i) = - GH(i);
525 //====================================================
526 Standard_Boolean Bounds(const math_Vector& InfBound,
527 const math_Vector& SupBound,
528 const math_Vector& Tol,
530 const math_Vector& SolSave,
531 math_IntegerVector& Constraints,
533 Standard_Boolean& theIsNewSol)
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
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 //======================================================
542 Standard_Boolean Out = Standard_False;
543 Standard_Integer i, Ninc = Sol.Length();
544 Standard_Real monratio = 1;
546 theIsNewSol = Standard_True;
548 // Calcul du ratio de recadrage
549 for (i = 1; i <= Ninc; i++) {
551 Delta(i) = Sol(i) - SolSave(i);
552 if (InfBound(i) == SupBound(i)) {
554 Out = Standard_True; // Ok mais, cela devrait etre eviter
556 else if(Sol(i) < InfBound(i)) {
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) );
563 else if (Sol(i) > SupBound(i)) {
566 // Delta(i) is positive
567 if (Delta(i) > Tol(i))
568 monratio = Min(monratio, (SupBound(i) - SolSave(i))/Delta(i) );
572 if (Out){ // Troncature et derniers recadrage pour blinder (pb numeriques)
573 if (monratio == 0.0) {
574 theIsNewSol = Standard_False;
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);
585 else if (Sol(i) > SupBound(i)) {
586 Sol(i) = SupBound(i);
587 Delta(i) = Sol(i) - SolSave(i);
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)
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),
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)
631 SetTolerance(theTolerance);
634 //=======================================================================
635 //function : math_FunctionSetRoot
636 //purpose : Constructor
637 //=======================================================================
638 math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& theFunction,
639 const Standard_Integer theNbIterations)
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),
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)
667 //=======================================================================
668 //function : ~math_FunctionSetRoot
669 //purpose : Destructor
670 //=======================================================================
671 math_FunctionSetRoot::~math_FunctionSetRoot()
675 //=======================================================================
676 //function : SetTolerance
678 //=======================================================================
679 void math_FunctionSetRoot::SetTolerance(const math_Vector& theTolerance)
681 for (Standard_Integer i = 1; i <= Tol.Length(); ++i)
682 Tol(i) = theTolerance(i);
685 //=======================================================================
688 //=======================================================================
689 void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& theFunction,
690 const math_Vector& theStartingPoint,
691 const Standard_Boolean theStopOnDivergent)
693 Perform(theFunction, theStartingPoint, InfBound, SupBound, theStopOnDivergent);
696 //=======================================================================
699 //=======================================================================
700 void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
701 const math_Vector& StartingPoint,
702 const math_Vector& theInfBound,
703 const math_Vector& theSupBound,
704 Standard_Boolean theStopOnDivergent)
706 Standard_Integer Ninc = F.NbVariables(), Neq = F.NbEquations();
709 (StartingPoint.Length()!= Ninc) ||
710 (theInfBound.Length() != Ninc) ||
711 (theSupBound.Length() != Ninc)) { throw Standard_DimensionError(); }
714 Standard_Boolean ChangeDirection = Standard_False, Sort = Standard_False, isNewSol = Standard_False;
715 Standard_Boolean Good, Verif;
716 Standard_Boolean Stop;
717 const Standard_Real EpsSqrt = 1.e-16, Eps = 1.e-32, Eps2 = 1.e-64, Progres = 0.005;
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
721 math_IntegerVector aConstraints(1, Ninc); // Pour savoir sur quels bord on se trouve
722 for (i = 1; i <= Ninc ; i++) {
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);
728 MyDirFunction F_Dir(Temp1, Temp2, Temp3, Temp4, F);
729 Standard_Integer DescenteIter;
731 Done = Standard_False;
736 myIsDivergent = Standard_False;
737 for (i = 1; i <= Ninc; i++)
739 myIsDivergent |= (Sol(i) < theInfBound(i)) | (theSupBound(i) < Sol(i));
741 if (theStopOnDivergent & myIsDivergent)
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++) {
748 if (Sol(i) <= theInfBound(i)) Sol(i) = theInfBound(i);
749 else if (Sol(i) > theSupBound(i)) Sol(i) = theSupBound(i);
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;
755 if (!theStopOnDivergent || !myIsDivergent)
757 State = F.GetStateNumber();
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...
765 Save(0) = Max (F2, EpsSqrt);
766 Standard_Real aTol_Func = Epsilon(F2);
767 FSR_DEBUG("=== Mode Debug de Function Set Root" << endl);
768 FSR_DEBUG(" F2 Initial = " << F2);
770 if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
771 Done = Standard_False;
772 if (!theStopOnDivergent || !myIsDivergent)
774 Done = Standard_True;
775 State = F.GetStateNumber();
780 for (Kount = 1; Kount <= Itermax; Kount++) {
781 PreviousMinimum = F2;
783 PreviousSolution = Sol;
786 SearchDirection(DF, GH, FF, ChangeDirection, InvLengthMax, DH, Dy);
787 if (Abs(Dy) <= Eps) {
788 Done = Standard_False;
789 if (!theStopOnDivergent || !myIsDivergent)
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();
799 if (ChangeDirection) {
800 Ambda = Ambda2 / Sqrt(Abs(Dy));
801 if (Ambda > 1.0) Ambda = 1.0;
805 Ambda2 = 0.5*Ambda/DH.Norm();
808 for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
809 Sol(i) = Sol(i) + Ambda * DH(i);
812 for (i = 1; i <= Ninc; i++)
814 myIsDivergent |= (Sol(i) < theInfBound(i)) | (theSupBound(i) < Sol(i));
816 if (theStopOnDivergent & myIsDivergent)
821 Sort = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave,
822 aConstraints, Delta, isNewSol);
827 // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
828 if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
829 Done = Standard_False;
830 if (!theStopOnDivergent || !myIsDivergent)
832 State = F.GetStateNumber();
838 FSR_DEBUG("Kount = " << Kount);
839 FSR_DEBUG("Le premier F2 = " << F2);
840 FSR_DEBUG("Direction = " << ChangeDirection);
842 if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
843 Done = Standard_False;
844 if (!theStopOnDivergent || !myIsDivergent)
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();
855 if (Sort || (F2/PreviousMinimum > Progres)) {
857 OldF = PreviousMinimum;
858 Stop = Standard_False;
859 Good = Standard_False;
861 Standard_Boolean Sortbis;
863 // -------------------------------------------------
864 // Traitement standard sans traitement des bords
865 // -------------------------------------------------
866 if (!Sort) { // si l'on est pas sortie on essaye de progresser en avant
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);
874 for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
875 Sol(i) = Sol(i) + Ambda * DH(i);
878 for (i = 1; i <= Ninc; i++)
880 myIsDivergent |= (Sol(i) < theInfBound(i)) | (theSupBound(i) < Sol(i));
882 if (theStopOnDivergent & myIsDivergent)
887 Stop = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave,
888 aConstraints, Delta, isNewSol);
889 FSR_DEBUG(" Augmentation de lambda");
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.
899 Good = MinimizeDirection(SolSave, Delta, OldF, F2, DHSave, GH,
902 else if (ChangeDirection || (DescenteIter>1)
903 || (OldF>PreviousMinimum) ) {
904 // La progression a ete utile, on minimise...
906 Good = MinimizeDirection(PreviousSolution, SolSave, Sol,
907 OldF, Delta, Tol, F_Dir);
916 for (i = 1; i <= Ninc; i++)
918 myIsDivergent |= (Sol(i) < theInfBound(i)) | (theSupBound(i) < Sol(i));
920 if (theStopOnDivergent & myIsDivergent)
925 Sort = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave,
926 aConstraints, Delta, isNewSol);
928 Sort = Standard_False; // On a rejete le point sur la frontiere
930 Stop = Standard_True; // et on sort dans tous les cas...
934 // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
935 if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
936 Done = Standard_False;
937 if (!theStopOnDivergent || !myIsDivergent)
939 State = F.GetStateNumber();
945 if (Abs(Dy) <= Eps) {
948 Done = Standard_False;
949 if (!theStopOnDivergent || !myIsDivergent)
951 Done = Standard_True;
952 ////modified by jgv, 31.08.2011////
953 F.Value(Sol, FF); //update F before GetStateNumber
954 ///////////////////////////////////
955 State = F.GetStateNumber();
959 if (DescenteIter >= 100) {
960 Stop = Standard_True;
963 FSR_DEBUG("--- Sortie du Traitement Standard");
964 FSR_DEBUG(" DescenteIter = "<<DescenteIter << " F2 = " << F2);
966 // ------------------------------------
967 // on passe au traitement des bords
968 // ------------------------------------
970 Stop = (F2 > 1.001*OldF); // Pour ne pas progresser sur le bord
973 while (Sortbis && ((F2<OldF)|| (DescenteIter == 0))
976 // On essaye de progresser sur le bord
979 SearchDirection(DF, GH, FF, aConstraints, Sol,
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) {
985 // Ambda = Ambda2 / Sqrt(Abs(Dy));
986 Ambda = Ambda2 / Sqrt(-Dy);
987 if (Ambda > 1.0) Ambda = 1.0;
991 Ambda2 = 0.5*Ambda/DH.Norm();
994 for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
995 Sol(i) = Sol(i) + Ambda * DH(i);
998 for (i = 1; i <= Ninc; i++)
1000 myIsDivergent |= (Sol(i) < theInfBound(i)) | (theSupBound(i) < Sol(i));
1002 if (theStopOnDivergent & myIsDivergent)
1007 Sortbis = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave,
1008 aConstraints, Delta, 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)
1017 State = F.GetStateNumber();
1023 FSR_DEBUG("--- Iteration au bords : " << DescenteIter);
1024 FSR_DEBUG("--- F2 = " << F2);
1027 Stop = Standard_True;
1030 while((F2/PreviousMinimum > Progres) && (F2<OldF) && (!Stop) ) {
1032 FSR_DEBUG("--- Iteration de descente conditionnel = " << DescenteIter);
1033 if (F2 < OldF && Dy < 0.0) {
1034 // On essaye de progresser dans cette direction.
1037 for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
1038 Sol(i) = Sol(i) + Ambda * DH(i);
1041 for (i = 1; i <= Ninc; i++)
1043 myIsDivergent |= (Sol(i) < theInfBound(i)) | (theSupBound(i) < Sol(i));
1045 if (theStopOnDivergent & myIsDivergent)
1050 Sortbis = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave,
1051 aConstraints, Delta, 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)
1060 State = F.GetStateNumber();
1067 Stop = ((Dy >=0) || (DescenteIter >= 10) || Sortbis);
1069 Stop = ((Dy >=0) || (DescenteIter >= 10));
1071 if (((F2/PreviousMinimum > Progres) &&
1072 (F2>=OldF))||(F2>=PreviousMinimum)) {
1073 // On minimise par Brent
1075 Good = MinimizeDirection(SolSave, Delta, OldF, F2,
1076 DHSave, GH, Tol, F_Dir);
1079 Sort = Standard_False;
1082 Sol = SolSave + Delta;
1084 for (i = 1; i <= Ninc; i++)
1086 myIsDivergent |= (Sol(i) < theInfBound(i)) | (theSupBound(i) < Sol(i));
1088 if (theStopOnDivergent & myIsDivergent)
1093 Sort = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave,
1094 aConstraints, Delta, 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)
1101 State = F.GetStateNumber();
1109 FSR_DEBUG("--- Sortie du Traitement des Bords");
1110 FSR_DEBUG("--- DescenteIter = "<<DescenteIter << " F2 = " << F2);
1114 // ---------------------------------------------
1115 // on passe aux tests d'ARRET
1116 // ---------------------------------------------
1118 // Est ce la solution ?
1119 if (ChangeDirection) Verif = Standard_True;
1120 // Gradient : Il faut eviter de boucler
1122 Verif = Standard_False;
1123 if (Kount > 1) { // Pour accelerer les cas quasi-quadratique
1124 if (Save(Kount-1)<1.e-4*Save(Kount-2)) Verif = Standard_True;
1126 else Verif = (F2 < 1.e-6*Save(0)); //Pour les cas dejas solutions
1129 for(i = Delta.Lower(); i <= Delta.Upper(); i++) {
1130 Delta(i) = PreviousSolution(i) - Sol(i);
1133 if (IsSolutionReached(F)) {
1134 if (PreviousMinimum < F2) {
1137 Done = Standard_False;
1138 if (!theStopOnDivergent || !myIsDivergent)
1140 Done = Standard_True;
1141 ////modified by jgv, 31.08.2011////
1142 F.Value(Sol, FF); //update F before GetStateNumber
1143 ///////////////////////////////////
1144 State = F.GetStateNumber();
1149 //fin du test solution
1151 // Analyse de la progression...
1152 //comparison of current minimum and previous minimum
1153 if ((F2 - PreviousMinimum) <= aTol_Func){
1155 // L'historique est il bon ?
1156 if (F2 >= 0.95*Save(Kount - 5)) {
1157 if (!ChangeDirection) ChangeDirection = Standard_True;
1160 Done = Standard_False;
1161 if (!theStopOnDivergent || !myIsDivergent)
1163 Done = Standard_True;
1164 State = F.GetStateNumber();
1166 return; // si un gain inf a 5% on sort
1169 else ChangeDirection = Standard_False; //Si oui on recommence
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) &&
1177 (F2 > 0.5*PreviousMinimum)) {
1178 ChangeDirection = Standard_True;
1181 // Si l'on ne decide pas de changer de strategie, on verifie,
1182 // si ce n'est dejas fait
1183 if ((!ChangeDirection) && (!Verif)) {
1184 for(i = Delta.Lower(); i <= Delta.Upper(); i++) {
1185 Delta(i) = PreviousSolution(i) - Sol(i);
1187 if (IsSolutionReached(F)) {
1188 Done = Standard_False;
1189 if (!theStopOnDivergent || !myIsDivergent)
1191 Done = Standard_True;
1192 ////modified by jgv, 31.08.2011////
1193 F.Value(Sol, FF); //update F before GetStateNumber
1194 ///////////////////////////////////
1195 State = F.GetStateNumber();
1201 else { // Cas de regression
1202 if (!ChangeDirection) { // On passe au gradient
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)
1210 State = F.GetStateNumber();
1218 if (!theStopOnDivergent || !myIsDivergent)
1220 State = F.GetStateNumber();
1222 return; // y a plus d'issues
1226 if (!theStopOnDivergent || !myIsDivergent)
1228 State = F.GetStateNumber();
1232 //=======================================================================
1235 //=======================================================================
1236 void math_FunctionSetRoot::Dump(Standard_OStream& o) const
1238 o << " math_FunctionSetRoot";
1240 o << " Status = Done\n";
1241 o << " Location value = " << Sol << "\n";
1242 o << " Number of iterations = " << Kount << "\n";
1245 o << "Status = Not Done\n";
1249 //=======================================================================
1252 //=======================================================================
1253 void math_FunctionSetRoot::Root(math_Vector& Root) const
1255 StdFail_NotDone_Raise_if(!Done, " ");
1256 Standard_DimensionError_Raise_if(Root.Length() != Sol.Length(), " ");
1260 //=======================================================================
1261 //function : FunctionSetErrors
1263 //=======================================================================
1264 void math_FunctionSetRoot::FunctionSetErrors(math_Vector& Err) const
1266 StdFail_NotDone_Raise_if(!Done, " ");
1267 Standard_DimensionError_Raise_if(Err.Length() != Sol.Length(), " ");