// On change de direction
//============================================================================
-static Standard_Boolean mydebug = Standard_False;
-
class MyDirFunction : public math_Function
{
F.Initialize(P1, Delta);
// (2) On minimise
- if (mydebug)cout << " minimisation dans la direction" << endl;
ax = -1; bx = 0;
cx = (P2-P1).Norm()*invnorme;
if (cx < 1.e-2) return Standard_False;
// (1) On realise une premiere interpolation quadratique
Standard_Real ax, bx, cx, df1, df2, Delta, tsol, fsol, tsolbis;
- if (mydebug) { cout << " essai d interpolation" << endl;}
-
df1 = Gradient*Dir;
df2 = DGradient*Dir;
Delta = bx*bx - 4*ax*cx;
if (Delta > 1.e-9) {
// il y a des racines, on prend la plus proche de 0
- Delta = Sqrt(Delta);
+ Delta = sqrt(Delta);
tsol = -(bx + Delta);
tsolbis = (Delta - bx);
if (Abs(tsolbis) < Abs(tsol)) tsol = tsolbis;
if (fsol<PValue) {
good = Standard_True;
Result = fsol;
- if (mydebug) cout << "t= "<<tsol<<" F = " << fsol << " OldF = "<<PValue <<endl;
}
// (2) Si l'on a pas assez progresser on realise une recherche
else {
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);
if(Sol.IsDone()) {
if (Sol.Minimum() <= Result) {
- tsol = Sol.Location();
- good = Standard_True;
- if (mydebug) cout << "t= "<<tsol<<" F ="<< Sol.Minimum()
- << " OldF = "<<Result <<endl;
+ tsol = Sol.Location();
+ good = Standard_True;
}
}
}
math_Gauss Solut(DF, 1.e-9);
if (Solut.IsDone()) Solut.Solve(Direction);
else { // we have to "forget" singular directions.
- if (mydebug) cout << " Matrice singuliere : On prend SVD" << endl;
- math_SVD SolvebySVD(DF);
+ math_SVD SolvebySVD(DF);
if (SolvebySVD.IsDone()) SolvebySVD.Solve(-1*FF, Direction);
- else ChangeDirection = Standard_True;
- }
+ else ChangeDirection = Standard_True;
+ }
}
else if (Ninc > Neq) {
math_SVD Solut(DF);
//====================================================
-Standard_Boolean Bounds(const math_Vector& InfBound,
- const math_Vector& SupBound,
- const math_Vector& Tol,
- math_Vector& Sol,
- const math_Vector& SolSave,
- math_IntegerVector& Constraints,
- math_Vector& Delta)
+Standard_Boolean Bounds(const math_Vector& InfBound,
+ const math_Vector& SupBound,
+ const math_Vector& Tol,
+ math_Vector& Sol,
+ const math_Vector& SolSave,
+ math_IntegerVector& Constraints,
+ math_Vector& Delta,
+ Standard_Boolean& theIsNewSol)
//
-// Purpose : Troncate un pas d'optimisation pour rester
-// dans le domaine, Delta donne le pas final
+// Purpose: Trims an initial solution Sol to be within a domain defined by
+// InfBound and SupBound. Delta will contain a distance between final Sol and
+// 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_Integer i, Ninc = Sol.Length();
Standard_Real monratio = 1;
+ theIsNewSol = Standard_True;
+
// Calcul du ratio de recadrage
for (i = 1; i <= Ninc; i++) {
Constraints(i) = 0;
else if(Sol(i) < InfBound(i)) {
Constraints(i) = 1;
Out = Standard_True;
- if (Abs(Delta(i)) > Tol(i)) // Afin d'eviter des ratio nulles pour rien
- monratio = Min(monratio, Abs( (SolSave(i)-InfBound(i))/Delta(i)) );
+ // Delta(i) is negative
+ if (-Delta(i) > Tol(i)) // Afin d'eviter des ratio nulles pour rien
+ monratio = Min(monratio, (InfBound(i) - SolSave(i))/Delta(i) );
}
else if (Sol(i) > SupBound(i)) {
Constraints(i) = 1;
Out = Standard_True;
- if (Abs(Delta(i)) > Tol(i))
- monratio = Min(monratio, Abs( (SolSave(i)-SupBound(i))/Delta(i)) );
+ // Delta(i) is positive
+ if (Delta(i) > Tol(i))
+ monratio = Min(monratio, (SupBound(i) - SolSave(i))/Delta(i) );
}
}
if (Out){ // Troncature et derniers recadrage pour blinder (pb numeriques)
- Delta *= monratio;
- Sol = SolSave+Delta;
- for (i = 1; i <= Ninc; i++) {
- if(Sol(i) < InfBound(i)) {
- Sol(i) = InfBound(i);
- Delta(i) = Sol(i) - SolSave(i);
- }
- else if (Sol(i) > SupBound(i)) {
- Sol(i) = SupBound(i);
- Delta(i) = Sol(i) - SolSave(i);
+ if (monratio == 0.0) {
+ theIsNewSol = Standard_False;
+ Sol = SolSave;
+ Delta.Init (0.0);
+ } else {
+ Delta *= monratio;
+ Sol = SolSave+Delta;
+ for (i = 1; i <= Ninc; i++) {
+ if(Sol(i) < InfBound(i)) {
+ Sol(i) = InfBound(i);
+ Delta(i) = Sol(i) - SolSave(i);
+ }
+ else if (Sol(i) > SupBound(i)) {
+ Sol(i) = SupBound(i);
+ Delta(i) = Sol(i) - SolSave(i);
+ }
}
}
}
(SupBound.Length() != Ninc)) { Standard_DimensionError:: Raise(); }
Standard_Integer i;
- Standard_Boolean ChangeDirection = Standard_False, Sort = Standard_False;
+ Standard_Boolean ChangeDirection = Standard_False, Sort = Standard_False, isNewSol = Standard_False;
Standard_Boolean Good, Verif;
Standard_Boolean Stop;
- Standard_Real Eps = 1.e-32, Progres = 0.005;
+ const Standard_Real EpsSqrt = 1.e-16, Eps = 1.e-32, Eps2 = 1.e-64, Progres = 0.005;
Standard_Real F2, PreviousMinimum, Dy, OldF;
Standard_Real Ambda, Ambda2, Gnr1, Oldgr;
math_Vector InvLengthMax(1, Ninc); // Pour bloquer les pas a 1/4 du domaine
// 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
// de faire une seconde iteration...
- Save(0) = Max (F2, Sqrt(Eps));
-
- if (mydebug) {
- cout << "=== Mode Debug de Function Set Root"<<endl;
- cout << " F2 Initial = " << F2 << endl;
- }
+ Save(0) = Max (F2, EpsSqrt);
- if ((F2 <= Eps) || (Sqrt(Gnr1) <= Eps)) {
+ if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
Done = Standard_True;
State = F.GetStateNumber();
return;
return;
}
if (ChangeDirection) {
- Ambda = Ambda2/Sqrt(Abs(Dy));
+ Ambda = Ambda2 / sqrt(Abs(Dy));
if (Ambda > 1.0) Ambda = 1.0;
}
else {
Sol(i) = Sol(i) + Ambda * DH(i);
}
Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
- Constraints, Delta);
+ Constraints, Delta, isNewSol);
DHSave = GH;
-// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
- if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
- Done = Standard_False;
- State = F.GetStateNumber();
- return;
- }
-
-#if DEB
- if (mydebug) {
- cout << "Kount = " << Kount << endl;
- cout << "Le premier F2 = " << F2 << endl;
- cout << "Direction = " << ChangeDirection << endl;
+ if (isNewSol) {
+// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+ if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
+ Done = Standard_False;
+ State = F.GetStateNumber();
+ return;
+ }
}
-#endif
- if ((F2 <= Eps) || (Sqrt(Gnr1) <= Eps)) {
+ if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
Done = Standard_True;
State = F.GetStateNumber();
return;
while((F2/PreviousMinimum > Progres) && !Stop) {
if (F2 < OldF && (Dy < 0.0)) {
// On essaye de progresser dans cette direction.
- if (mydebug) cout << " iteration de descente = " << DescenteIter<<endl;
DescenteIter++;
SolSave = Sol;
OldF = F2;
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,
- Constraints, Delta);
- if (mydebug) { cout << " Augmentation de lambda" << endl;}
+ Constraints, Delta, isNewSol);
Ambda *= 1.7;
}
else {
F2 = OldF;
}
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
}
Stop = Standard_True; // et on sort dans tous les cas...
}
DHSave = GH;
-// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
- if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
- Done = Standard_False;
- State = F.GetStateNumber();
- return;
- }
+ if (isNewSol) {
+// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+ if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
+ Done = Standard_False;
+ State = F.GetStateNumber();
+ return;
+ }
+ }
Dy = GH*DH;
if (Abs(Dy) <= Eps) {
State = F.GetStateNumber();
Stop = Standard_True;
}
}
- if (mydebug) {
- cout << "--- Sortie du Traitement Standard"<<endl;
- cout << " DescenteIter = "<<DescenteIter << " F2 = " << F2 << endl;
- }
}
// ------------------------------------
// on passe au traitement des bords
OldF = F2;
SearchDirection(DF, GH, FF, Constraints, Sol,
ChangeDirection, InvLengthMax, DH, Dy);
- if (mydebug) { cout << " Conditional Direction = " << ChangeDirection << endl;}
if (Dy<-Eps) { //Pour eviter des calculs inutiles et des /0...
if (ChangeDirection) {
-// Ambda = Ambda2/Sqrt(Abs(Dy));
- Ambda = Ambda2/Sqrt(-Dy);
+// Ambda = Ambda2 / sqrt(Abs(Dy));
+ Ambda = Ambda2 / sqrt(-Dy);
if (Ambda > 1.0) Ambda = 1.0;
}
else {
Sol(i) = Sol(i) + Ambda * DH(i);
}
Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
- Constraints, Delta);
+ Constraints, Delta, isNewSol);
DHSave = GH;
-// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
- if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
- Done = Standard_False;
- State = F.GetStateNumber();
- return;
- }
+ if (isNewSol) {
+// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+ if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
+ Done = Standard_False;
+ State = F.GetStateNumber();
+ return;
+ }
+ }
Ambda2 = Gnr1;
- if (mydebug) {
- cout << "--- Iteration au bords : " << DescenteIter << endl;
- cout << "--- F2 = " << F2 << endl;
- }
}
else {
Stop = Standard_True;
while((F2/PreviousMinimum > Progres) && (F2<OldF) && (!Stop) ) {
DescenteIter++;
- if (mydebug) cout << "--- Iteration de descente conditionnel = "
- << DescenteIter<<endl;
if (F2 < OldF && Dy < 0.0) {
// On essaye de progresser dans cette direction.
SolSave = Sol;
Sol(i) = Sol(i) + Ambda * DH(i);
}
Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
- Constraints, Delta );
+ Constraints, Delta, isNewSol);
}
DHSave = GH;
-// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
- if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
- Done = Standard_False;
- State = F.GetStateNumber();
- return;
- }
+ if (isNewSol) {
+// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+ if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
+ Done = Standard_False;
+ State = F.GetStateNumber();
+ return;
+ }
+ }
Ambda2 = Gnr1;
Dy = GH*DH;
Stop = ((Dy >=0) || (DescenteIter >= 10) || Sortbis);
DHSave, GH, Tol, F_Dir);
if (!Good) {
Sol = SolSave;
+ Sort = Standard_False;
}
else {
Sol = SolSave + Delta;
- }
- Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
- Constraints, Delta);
-// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
- if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
- Done = Standard_False;
- State = F.GetStateNumber();
- return;
+ Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
+ Constraints, Delta, isNewSol);
+ if (isNewSol) {
+// F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+ if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
+ Done = Standard_False;
+ State = F.GetStateNumber();
+ return;
+ }
+ }
}
Dy = GH*DH;
}
- if (mydebug) {
- cout << "--- Sortie du Traitement des Bords"<<endl;
- cout << "--- DescenteIter = "<<DescenteIter << " F2 = " << F2 << endl;
- }
}
}
Standard_DimensionError_Raise_if(Err.Length() != Sol.Length(), " ");
Err = Delta;
}
-
-
-
-
-
-
-