OpenCascade/occt
1
0
Fork 0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-04-09 18:50:54 +03:00
Code Packages Releases Activity

Rollback integration OCC22567 Speed up of math_FunctionSetRoot (used in Extrema)

Browse Source
This commit is contained in:
RLN and KGV 2011-07-22 07:59:36 +00:00 committed by bugmaster
parent e33e7e78af
commit fbadd2ccce
1 changed files with 124 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

226
src/math/math_FunctionSetRoot.cxx
View File

@ -49,6 +49,8 @@
// On change de direction // On change de direction
//============================================================================ //============================================================================
static Standard_Boolean mydebug = Standard_False;
class MyDirFunction : public math_Function class MyDirFunction : public math_Function
{ {
@ -200,6 +202,7 @@ static Standard_Boolean MinimizeDirection(const math_Vector& P0,
F.Initialize(P1, Delta); F.Initialize(P1, Delta);
// (2) On minimise // (2) On minimise
if (mydebug)cout << " minimisation dans la direction" << endl;
ax = -1; bx = 0; ax = -1; bx = 0;
cx = (P2-P1).Norm()*invnorme; cx = (P2-P1).Norm()*invnorme;
if (cx < 1.e-2) return Standard_False; if (cx < 1.e-2) return Standard_False;
@ -240,6 +243,8 @@ static Standard_Boolean MinimizeDirection(const math_Vector& P,
// (1) On realise une premiere interpolation quadratique // (1) On realise une premiere interpolation quadratique
Standard_Real ax, bx, cx, df1, df2, Delta, tsol, fsol, tsolbis; Standard_Real ax, bx, cx, df1, df2, Delta, tsol, fsol, tsolbis;
if (mydebug) { cout << " essai d interpolation" << endl;}
df1 = Gradient*Dir; df1 = Gradient*Dir;
df2 = DGradient*Dir; df2 = DGradient*Dir;
@ -259,7 +264,7 @@ static Standard_Boolean MinimizeDirection(const math_Vector& P,
Delta = bx*bx - 4*ax*cx; Delta = bx*bx - 4*ax*cx;
if (Delta > 1.e-9) { if (Delta > 1.e-9) {
// il y a des racines, on prend la plus proche de 0 // il y a des racines, on prend la plus proche de 0
Delta = sqrt(Delta); Delta = Sqrt(Delta);
tsol = -(bx + Delta); tsol = -(bx + Delta);
tsolbis = (Delta - bx); tsolbis = (Delta - bx);
if (Abs(tsolbis) < Abs(tsol)) tsol = tsolbis; if (Abs(tsolbis) < Abs(tsol)) tsol = tsolbis;
@ -280,6 +285,7 @@ static Standard_Boolean MinimizeDirection(const math_Vector& P,
if (fsol<PValue) { if (fsol<PValue) {
good = Standard_True; good = Standard_True;
Result = fsol; Result = fsol;
if (mydebug) cout << "t= "<<tsol<<" F = " << fsol << " OldF = "<<PValue <<endl;
} }
// (2) Si l'on a pas assez progresser on realise une recherche // (2) Si l'on a pas assez progresser on realise une recherche
@ -292,11 +298,14 @@ static Standard_Boolean MinimizeDirection(const math_Vector& P,
else { else {
ax = 0.0; bx = tsol; cx = 1.0; ax = 0.0; bx = tsol; cx = 1.0;
} }
if (mydebug) cout << " minimisation dans la direction" << endl;
math_BrentMinimum Sol(F, ax, bx, cx, tol1d, 100, tol1d); math_BrentMinimum Sol(F, ax, bx, cx, tol1d, 100, tol1d);
if(Sol.IsDone()) { if(Sol.IsDone()) {
if (Sol.Minimum() <= Result) { if (Sol.Minimum() <= Result) {
tsol = Sol.Location(); tsol = Sol.Location();
good = Standard_True; good = Standard_True;
if (mydebug) cout << "t= "<<tsol<<" F ="<< Sol.Minimum()
<< " OldF = "<<Result <<endl;
} }
} }
} }
@ -327,10 +336,11 @@ static void SearchDirection(const math_Matrix& DF,
math_Gauss Solut(DF, 1.e-9); math_Gauss Solut(DF, 1.e-9);
if (Solut.IsDone()) Solut.Solve(Direction); if (Solut.IsDone()) Solut.Solve(Direction);
else { // we have to "forget" singular directions. else { // we have to "forget" singular directions.
math_SVD SolvebySVD(DF); if (mydebug) cout << " Matrice singuliere : On prend SVD" << endl;
math_SVD SolvebySVD(DF);
if (SolvebySVD.IsDone()) SolvebySVD.Solve(-1*FF, Direction); if (SolvebySVD.IsDone()) SolvebySVD.Solve(-1*FF, Direction);
else ChangeDirection = Standard_True; else ChangeDirection = Standard_True;
} }
} }
else if (Ninc > Neq) { else if (Ninc > Neq) {
math_SVD Solut(DF); math_SVD Solut(DF);
@ -447,28 +457,22 @@ static void SearchDirection(const math_Matrix& DF,
//==================================================== //====================================================
Standard_Boolean Bounds(const math_Vector& InfBound, Standard_Boolean Bounds(const math_Vector& InfBound,
const math_Vector& SupBound, const math_Vector& SupBound,
const math_Vector& Tol, const math_Vector& Tol,
math_Vector& Sol, math_Vector& Sol,
const math_Vector& SolSave, const math_Vector& SolSave,
math_IntegerVector& Constraints, math_IntegerVector& Constraints,
math_Vector& Delta, math_Vector& Delta)
Standard_Boolean& theIsNewSol)
// //
// Purpose: Trims an initial solution Sol to be within a domain defined by // Purpose : Troncate un pas d'optimisation pour rester
// InfBound and SupBound. Delta will contain a distance between final Sol and // dans le domaine, Delta donne le pas final
// SolSave.
// IsNewSol returns False, if final Sol fully coincides with SolSave, i.e.
// if SolSave already lied on a boundary and initial Sol was fully beyond it
//====================================================== //======================================================
{ {
Standard_Boolean Out = Standard_False; Standard_Boolean Out = Standard_False;
Standard_Integer i, Ninc = Sol.Length(); Standard_Integer i, Ninc = Sol.Length();
Standard_Real monratio = 1; Standard_Real monratio = 1;
theIsNewSol = Standard_True;
// Calcul du ratio de recadrage // Calcul du ratio de recadrage
for (i = 1; i <= Ninc; i++) { for (i = 1; i <= Ninc; i++) {
Constraints(i) = 0; Constraints(i) = 0;
@ -480,36 +484,28 @@ Standard_Boolean Bounds(const math_Vector& InfBound,
else if(Sol(i) < InfBound(i)) { else if(Sol(i) < InfBound(i)) {
Constraints(i) = 1; Constraints(i) = 1;
Out = Standard_True; Out = Standard_True;
// Delta(i) is negative if (Abs(Delta(i)) > Tol(i)) // Afin d'eviter des ratio nulles pour rien
if (-Delta(i) > Tol(i)) // Afin d'eviter des ratio nulles pour rien monratio = Min(monratio, Abs( (SolSave(i)-InfBound(i))/Delta(i)) );
monratio = Min(monratio, (InfBound(i) - SolSave(i))/Delta(i) );
} }
else if (Sol(i) > SupBound(i)) { else if (Sol(i) > SupBound(i)) {
Constraints(i) = 1; Constraints(i) = 1;
Out = Standard_True; Out = Standard_True;
// Delta(i) is positive if (Abs(Delta(i)) > Tol(i))
if (Delta(i) > Tol(i)) monratio = Min(monratio, Abs( (SolSave(i)-SupBound(i))/Delta(i)) );
monratio = Min(monratio, (SupBound(i) - SolSave(i))/Delta(i) );
} }
} }
if (Out){ // Troncature et derniers recadrage pour blinder (pb numeriques) if (Out){ // Troncature et derniers recadrage pour blinder (pb numeriques)
if (monratio == 0.0) { Delta *= monratio;
theIsNewSol = Standard_False; Sol = SolSave+Delta;
Sol = SolSave; for (i = 1; i <= Ninc; i++) {
Delta.Init (0.0); if(Sol(i) < InfBound(i)) {
} else { Sol(i) = InfBound(i);
Delta *= monratio; Delta(i) = Sol(i) - SolSave(i);
Sol = SolSave+Delta; }
for (i = 1; i <= Ninc; i++) { else if (Sol(i) > SupBound(i)) {
if(Sol(i) < InfBound(i)) { Sol(i) = SupBound(i);
Sol(i) = InfBound(i); Delta(i) = Sol(i) - SolSave(i);
Delta(i) = Sol(i) - SolSave(i);
}
else if (Sol(i) > SupBound(i)) {
Sol(i) = SupBound(i);
Delta(i) = Sol(i) - SolSave(i);
}
} }
} }
} }
@ -680,10 +676,10 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
(SupBound.Length() != Ninc)) { Standard_DimensionError:: Raise(); } (SupBound.Length() != Ninc)) { Standard_DimensionError:: Raise(); }
Standard_Integer i; Standard_Integer i;
Standard_Boolean ChangeDirection = Standard_False, Sort = Standard_False, isNewSol = Standard_False; Standard_Boolean ChangeDirection = Standard_False, Sort = Standard_False;
Standard_Boolean Good, Verif; Standard_Boolean Good, Verif;
Standard_Boolean Stop; Standard_Boolean Stop;
const Standard_Real EpsSqrt = 1.e-16, Eps = 1.e-32, Eps2 = 1.e-64, Progres = 0.005; Standard_Real Eps = 1.e-32, Progres = 0.005;
Standard_Real F2, PreviousMinimum, Dy, OldF; Standard_Real F2, PreviousMinimum, Dy, OldF;
Standard_Real Ambda, Ambda2, Gnr1, Oldgr; Standard_Real Ambda, Ambda2, Gnr1, Oldgr;
math_Vector InvLengthMax(1, Ninc); // Pour bloquer les pas a 1/4 du domaine math_Vector InvLengthMax(1, Ninc); // Pour bloquer les pas a 1/4 du domaine
@ -718,9 +714,14 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
// Le rang 0 de Save ne doit servir q'au test accelarteur en fin de boucle // Le rang 0 de Save ne doit servir q'au test accelarteur en fin de boucle
// s'il on est dejas sur la solution, il faut leurer ce test pour eviter // s'il on est dejas sur la solution, il faut leurer ce test pour eviter
// de faire une seconde iteration... // de faire une seconde iteration...
Save(0) = Max (F2, EpsSqrt); Save(0) = Max (F2, Sqrt(Eps));
if ((F2 <= Eps) || (Gnr1 <= Eps2)) { if (mydebug) {
cout << "=== Mode Debug de Function Set Root"<<endl;
cout << " F2 Initial = " << F2 << endl;
}
if ((F2 <= Eps) || (Sqrt(Gnr1) <= Eps)) {
Done = Standard_True; Done = Standard_True;
State = F.GetStateNumber(); State = F.GetStateNumber();
return; return;
@ -739,7 +740,7 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
return; return;
} }
if (ChangeDirection) { if (ChangeDirection) {
Ambda = Ambda2 / sqrt(Abs(Dy)); Ambda = Ambda2/Sqrt(Abs(Dy));
if (Ambda > 1.0) Ambda = 1.0; if (Ambda > 1.0) Ambda = 1.0;
} }
else { else {
@ -751,20 +752,26 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
Sol(i) = Sol(i) + Ambda * DH(i); Sol(i) = Sol(i) + Ambda * DH(i);
} }
Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave, Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
Constraints, Delta, isNewSol); Constraints, Delta);
DHSave = GH; DHSave = GH;
if (isNewSol) { // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1); if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) { Done = Standard_False;
Done = Standard_False; State = F.GetStateNumber();
State = F.GetStateNumber(); return;
return;
}
} }
if ((F2 <= Eps) || (Gnr1 <= Eps2)) { #if DEB
if (mydebug) {
cout << "Kount = " << Kount << endl;
cout << "Le premier F2 = " << F2 << endl;
cout << "Direction = " << ChangeDirection << endl;
}
#endif
if ((F2 <= Eps) || (Sqrt(Gnr1) <= Eps)) {
Done = Standard_True; Done = Standard_True;
State = F.GetStateNumber(); State = F.GetStateNumber();
return; return;
@ -785,14 +792,16 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
while((F2/PreviousMinimum > Progres) && !Stop) { while((F2/PreviousMinimum > Progres) && !Stop) {
if (F2 < OldF && (Dy < 0.0)) { if (F2 < OldF && (Dy < 0.0)) {
// On essaye de progresser dans cette direction. // On essaye de progresser dans cette direction.
if (mydebug) cout << " iteration de descente = " << DescenteIter<<endl;
DescenteIter++; DescenteIter++;
SolSave = Sol; SolSave = Sol;
OldF = F2; OldF = F2;
for( i = Sol.Lower(); i <= Sol.Upper(); i++) { for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
Sol(i) = Sol(i) + Ambda * DH(i); Sol(i) = Sol(i) + Ambda * DH(i);
} }
Stop = Bounds(InfBound, SupBound, Tol, Sol, SolSave, Stop = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
Constraints, Delta, isNewSol); Constraints, Delta);
if (mydebug) { cout << " Augmentation de lambda" << endl;}
Ambda *= 1.7; Ambda *= 1.7;
} }
else { else {
@ -817,23 +826,21 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
F2 = OldF; F2 = OldF;
} }
else { else {
Sol = SolSave+Delta; Sol = SolSave+Delta;
Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
Constraints, Delta, isNewSol);
} }
Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
Constraints, Delta);
Sort = Standard_False; // On a rejete le point sur la frontiere Sort = Standard_False; // On a rejete le point sur la frontiere
} }
Stop = Standard_True; // et on sort dans tous les cas... Stop = Standard_True; // et on sort dans tous les cas...
} }
DHSave = GH; DHSave = GH;
if (isNewSol) { // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1); if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) { Done = Standard_False;
Done = Standard_False; State = F.GetStateNumber();
State = F.GetStateNumber(); return;
return; }
}
}
Dy = GH*DH; Dy = GH*DH;
if (Abs(Dy) <= Eps) { if (Abs(Dy) <= Eps) {
State = F.GetStateNumber(); State = F.GetStateNumber();
@ -845,6 +852,10 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
Stop = Standard_True; Stop = Standard_True;
} }
} }
if (mydebug) {
cout << "--- Sortie du Traitement Standard"<<endl;
cout << " DescenteIter = "<<DescenteIter << " F2 = " << F2 << endl;
}
} }
// ------------------------------------ // ------------------------------------
// on passe au traitement des bords // on passe au traitement des bords
@ -861,11 +872,12 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
OldF = F2; OldF = F2;
SearchDirection(DF, GH, FF, Constraints, Sol, SearchDirection(DF, GH, FF, Constraints, Sol,
ChangeDirection, InvLengthMax, DH, Dy); ChangeDirection, InvLengthMax, DH, Dy);
if (mydebug) { cout << " Conditional Direction = " << ChangeDirection << endl;}
if (Dy<-Eps) { //Pour eviter des calculs inutiles et des /0... if (Dy<-Eps) { //Pour eviter des calculs inutiles et des /0...
if (ChangeDirection) { if (ChangeDirection) {
// Ambda = Ambda2 / sqrt(Abs(Dy)); // Ambda = Ambda2/Sqrt(Abs(Dy));
Ambda = Ambda2 / sqrt(-Dy); Ambda = Ambda2/Sqrt(-Dy);
if (Ambda > 1.0) Ambda = 1.0; if (Ambda > 1.0) Ambda = 1.0;
} }
else { else {
@ -877,18 +889,20 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
Sol(i) = Sol(i) + Ambda * DH(i); Sol(i) = Sol(i) + Ambda * DH(i);
} }
Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave, Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
Constraints, Delta, isNewSol); Constraints, Delta);
DHSave = GH; DHSave = GH;
if (isNewSol) { // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1); if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) { Done = Standard_False;
Done = Standard_False; State = F.GetStateNumber();
State = F.GetStateNumber(); return;
return; }
}
}
Ambda2 = Gnr1; Ambda2 = Gnr1;
if (mydebug) {
cout << "--- Iteration au bords : " << DescenteIter << endl;
cout << "--- F2 = " << F2 << endl;
}
} }
else { else {
Stop = Standard_True; Stop = Standard_True;
@ -896,6 +910,8 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
while((F2/PreviousMinimum > Progres) && (F2<OldF) && (!Stop) ) { while((F2/PreviousMinimum > Progres) && (F2<OldF) && (!Stop) ) {
DescenteIter++; DescenteIter++;
if (mydebug) cout << "--- Iteration de descente conditionnel = "
<< DescenteIter<<endl;
if (F2 < OldF && Dy < 0.0) { if (F2 < OldF && Dy < 0.0) {
// On essaye de progresser dans cette direction. // On essaye de progresser dans cette direction.
SolSave = Sol; SolSave = Sol;
@ -904,17 +920,15 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
Sol(i) = Sol(i) + Ambda * DH(i); Sol(i) = Sol(i) + Ambda * DH(i);
} }
Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave, Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
Constraints, Delta, isNewSol); Constraints, Delta );
} }
DHSave = GH; DHSave = GH;
if (isNewSol) { // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1); if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) { Done = Standard_False;
Done = Standard_False; State = F.GetStateNumber();
State = F.GetStateNumber(); return;
return; }
}
}
Ambda2 = Gnr1; Ambda2 = Gnr1;
Dy = GH*DH; Dy = GH*DH;
Stop = ((Dy >=0) || (DescenteIter >= 10) || Sortbis); Stop = ((Dy >=0) || (DescenteIter >= 10) || Sortbis);
@ -929,23 +943,24 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
DHSave, GH, Tol, F_Dir); DHSave, GH, Tol, F_Dir);
if (!Good) { if (!Good) {
Sol = SolSave; Sol = SolSave;
Sort = Standard_False;
} }
else { else {
Sol = SolSave + Delta; Sol = SolSave + Delta;
Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave, }
Constraints, Delta, isNewSol); Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
if (isNewSol) { Constraints, Delta);
// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1); // F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) { if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
Done = Standard_False; Done = Standard_False;
State = F.GetStateNumber(); State = F.GetStateNumber();
return; return;
}
}
} }
Dy = GH*DH; Dy = GH*DH;
} }
if (mydebug) {
cout << "--- Sortie du Traitement des Bords"<<endl;
cout << "--- DescenteIter = "<<DescenteIter << " F2 = " << F2 << endl;
}
} }
} }
@ -1063,3 +1078,10 @@ void math_FunctionSetRoot::FunctionSetErrors(math_Vector& Err) const{
Standard_DimensionError_Raise_if(Err.Length() != Sol.Length(), " "); Standard_DimensionError_Raise_if(Err.Length() != Sol.Length(), " ");
Err = Delta; Err = Delta;
} }