-//static const char* sccsid = "@(#)math_FunctionSetRoot.cxx 3.2 95/01/10"; // Do not delete this line. Used by sccs.
+// Copyright (c) 1997-1999 Matra Datavision
+// Copyright (c) 1999-2012 OPEN CASCADE SAS
+//
+// The content of this file is subject to the Open CASCADE Technology Public
+// License Version 6.5 (the "License"). You may not use the content of this file
+// except in compliance with the License. Please obtain a copy of the License
+// at http://www.opencascade.org and read it completely before using this file.
+//
+// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
+// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
+//
+// The Original Code and all software distributed under the License is
+// distributed on an "AS IS" basis, without warranty of any kind, and the
+// Initial Developer hereby disclaims all such warranties, including without
+// limitation, any warranties of merchantability, fitness for a particular
+// purpose or non-infringement. Please see the License for the specific terms
+// and conditions governing the rights and limitations under the License.
+
// pmn 15/05/97 pas de Gauss avec des pivot trop petit. SVD fait mieux
// l'affaire + limitation de la longeur du pas + qq comentaire issus d'EUCLID3
// pmn 10/06/97 refonte totale du traitement des bords + ajustement des init
//===========================================================================
// - A partir d une solution de depart, recherche d une direction.( Newton la
// plupart du temps, gradient si Newton echoue.
-//
// - Recadrage au niveau des bornes avec recalcul de la direction si une
// inconnue a une valeur imposee.
-//
// -Si On ne sort pas des bornes
// Tant que (On ne progresse pas assez ou on ne change pas de direction)
// . Si (Progression encore possible)
// Sinon
// Si on depasse le minimum
// On fait une interpolation parabolique.
-//
// - Si on a progresse sur F
// On fait les tests d'arret
// Sinon
// On change de direction
//============================================================================
+#define FSR_DEBUG(arg)
+// Uncomment the following code to have debug output to cout
+/* * /
+static Standard_Boolean mydebug = Standard_False;
+#undef FSR_DEBUG
+#define FSR_DEBUG(arg) {if (mydebug) { cout << arg << endl; }}
+/* */
+
class MyDirFunction : public math_Function
{
F.Initialize(P1, Delta);
// (2) On minimise
+ FSR_DEBUG (" minimisation dans la direction")
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;
+ FSR_DEBUG(" essai d interpolation")
+
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;
+ FSR_DEBUG("t= "<<tsol<<" F = " << fsol << " OldF = "<<PValue)
}
// (2) Si l'on a pas assez progresser on realise une recherche
else {
ax = 0.0; bx = tsol; cx = 1.0;
}
+ FSR_DEBUG(" minimisation dans la direction")
math_BrentMinimum Sol(F, ax, bx, cx, tol1d, 100, tol1d);
if(Sol.IsDone()) {
if (Sol.Minimum() <= Result) {
tsol = Sol.Location();
good = Standard_True;
+ FSR_DEBUG("t= "<<tsol<<" F ="<< Sol.Minimum() << " OldF = "<<Result)
}
}
}
math_Gauss Solut(DF, 1.e-9);
if (Solut.IsDone()) Solut.Solve(Direction);
else { // we have to "forget" singular directions.
+ FSR_DEBUG(" Matrice singuliere : On prend SVD")
math_SVD SolvebySVD(DF);
if (SolvebySVD.IsDone()) SolvebySVD.Solve(-1*FF, Direction);
else ChangeDirection = Standard_True;
const math_Vector& Tolerance,
const math_Vector& infBound,
const math_Vector& supBound,
- const Standard_Integer NbIterations) :
+ const Standard_Integer NbIterations,
+ Standard_Boolean theStopOnDivergent) :
Delta(1, F.NbVariables()),
Sol(1, F.NbVariables()),
DF(1, F.NbEquations(),
Tol(i) =Tolerance(i);
}
Itermax = NbIterations;
- Perform(F, StartingPoint, infBound, supBound);
+ Perform(F, StartingPoint, infBound, supBound, theStopOnDivergent);
}
math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
const math_Vector& StartingPoint,
const math_Vector& InfBound,
- const math_Vector& SupBound)
+ const math_Vector& SupBound,
+ Standard_Boolean theStopOnDivergent)
{
-
Standard_Integer Ninc = F.NbVariables(), Neq = F.NbEquations();
if ((Neq <= 0) ||
Sol = StartingPoint;
Kount = 0;
+ //
+ myIsDivergent = Standard_False;
+ for (i = 1; i <= Ninc; i++)
+ {
+ myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+ }
+ if (theStopOnDivergent & myIsDivergent)
+ {
+ return;
+ }
// Verification de la validite des inconnues par rapport aux bornes.
// Recentrage sur les bornes si pas valide.
for ( i = 1; i <= Ninc; i++) {
// Calcul de la premiere valeur de F et de son gradient
if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
Done = Standard_False;
- State = F.GetStateNumber();
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ State = F.GetStateNumber();
+ }
return;
}
Ambda2 = Gnr1;
// s'il on est dejas sur la solution, il faut leurer ce test pour eviter
// de faire une seconde iteration...
Save(0) = Max (F2, EpsSqrt);
+ Standard_Real aTol_Func = Epsilon(F2);
+ FSR_DEBUG("=== Mode Debug de Function Set Root" << endl)
+ FSR_DEBUG(" F2 Initial = " << F2)
if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
- Done = Standard_True;
- State = F.GetStateNumber();
+ Done = Standard_False;
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ Done = Standard_True;
+ State = F.GetStateNumber();
+ }
return;
}
SearchDirection(DF, GH, FF, ChangeDirection, InvLengthMax, DH, Dy);
if (Abs(Dy) <= Eps) {
- Done = Standard_True;
- State = F.GetStateNumber();
+ Done = Standard_False;
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ Done = Standard_True;
+ ////modified by jgv, 31.08.2011////
+ F.Value(Sol, FF); //update F before GetStateNumber
+ ///////////////////////////////////
+ State = F.GetStateNumber();
+ }
return;
}
if (ChangeDirection) {
- Ambda = Ambda2 / sqrt(Abs(Dy));
+ Ambda = Ambda2 / Sqrt(Abs(Dy));
if (Ambda > 1.0) Ambda = 1.0;
}
else {
for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
Sol(i) = Sol(i) + Ambda * DH(i);
}
+ //
+ for (i = 1; i <= Ninc; i++)
+ {
+ myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+ }
+ if (theStopOnDivergent & myIsDivergent)
+ {
+ return;
+ }
+ //
Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
Constraints, Delta, 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();
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ State = F.GetStateNumber();
+ }
return;
}
}
+
+ FSR_DEBUG("Kount = " << Kount)
+ FSR_DEBUG("Le premier F2 = " << F2)
+ FSR_DEBUG("Direction = " << ChangeDirection)
if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
- Done = Standard_True;
- State = F.GetStateNumber();
+ Done = Standard_False;
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ Done = Standard_True;
+ ////modified by jgv, 31.08.2011////
+ F.Value(Sol, FF); //update F before GetStateNumber
+ ///////////////////////////////////
+ State = F.GetStateNumber();
+ }
return;
}
while((F2/PreviousMinimum > Progres) && !Stop) {
if (F2 < OldF && (Dy < 0.0)) {
// On essaye de progresser dans cette direction.
+ FSR_DEBUG(" iteration de descente = " << DescenteIter)
DescenteIter++;
SolSave = Sol;
OldF = F2;
for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
Sol(i) = Sol(i) + Ambda * DH(i);
}
+ //
+ for (i = 1; i <= Ninc; i++)
+ {
+ myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+ }
+ if (theStopOnDivergent & myIsDivergent)
+ {
+ return;
+ }
+ //
Stop = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
Constraints, Delta, isNewSol);
+ FSR_DEBUG(" Augmentation de lambda")
Ambda *= 1.7;
}
else {
}
else {
Sol = SolSave+Delta;
+ //
+ for (i = 1; i <= Ninc; i++)
+ {
+ myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+ }
+ if (theStopOnDivergent & myIsDivergent)
+ {
+ return;
+ }
+ //
Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
Constraints, Delta, 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();
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ State = F.GetStateNumber();
+ }
return;
}
}
Dy = GH*DH;
if (Abs(Dy) <= Eps) {
- State = F.GetStateNumber();
- Done = Standard_True;
- if (F2 > OldF) Sol = SolSave;
+ if (F2 > OldF)
+ Sol = SolSave;
+ Done = Standard_False;
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ Done = Standard_True;
+ ////modified by jgv, 31.08.2011////
+ F.Value(Sol, FF); //update F before GetStateNumber
+ ///////////////////////////////////
+ State = F.GetStateNumber();
+ }
return;
}
if (DescenteIter >= 10) {
Stop = Standard_True;
}
}
+ FSR_DEBUG("--- Sortie du Traitement Standard")
+ FSR_DEBUG(" DescenteIter = "<<DescenteIter << " F2 = " << F2)
}
// ------------------------------------
// on passe au traitement des bords
OldF = F2;
SearchDirection(DF, GH, FF, Constraints, Sol,
ChangeDirection, InvLengthMax, DH, Dy);
+ FSR_DEBUG(" Conditional Direction = " << ChangeDirection)
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 {
for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
Sol(i) = Sol(i) + Ambda * DH(i);
}
+ //
+ for (i = 1; i <= Ninc; i++)
+ {
+ myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+ }
+ if (theStopOnDivergent & myIsDivergent)
+ {
+ return;
+ }
+ //
Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
Constraints, Delta, 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();
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ State = F.GetStateNumber();
+ }
return;
}
}
Ambda2 = Gnr1;
+ FSR_DEBUG("--- Iteration au bords : " << DescenteIter)
+ FSR_DEBUG("--- F2 = " << F2)
}
else {
Stop = Standard_True;
while((F2/PreviousMinimum > Progres) && (F2<OldF) && (!Stop) ) {
DescenteIter++;
+ FSR_DEBUG("--- Iteration de descente conditionnel = " << DescenteIter)
if (F2 < OldF && Dy < 0.0) {
// On essaye de progresser dans cette direction.
SolSave = Sol;
for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
Sol(i) = Sol(i) + Ambda * DH(i);
}
+ //
+ for (i = 1; i <= Ninc; i++)
+ {
+ myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+ }
+ if (theStopOnDivergent & myIsDivergent)
+ {
+ return;
+ }
+ //
Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
Constraints, Delta, 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();
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ State = F.GetStateNumber();
+ }
return;
}
}
}
else {
Sol = SolSave + Delta;
+ //
+ for (i = 1; i <= Ninc; i++)
+ {
+ myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+ }
+ if (theStopOnDivergent & myIsDivergent)
+ {
+ 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();
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ State = F.GetStateNumber();
+ }
return;
}
}
}
Dy = GH*DH;
}
+ FSR_DEBUG("--- Sortie du Traitement des Bords")
+ FSR_DEBUG("--- DescenteIter = "<<DescenteIter << " F2 = " << F2)
}
}
for(i = Delta.Lower(); i <= Delta.Upper(); i++) {
Delta(i) = PreviousSolution(i) - Sol(i);
}
+
if (IsSolutionReached(F)) {
if (PreviousMinimum < F2) {
Sol = SolSave;
}
- State = F.GetStateNumber();
- Done = Standard_True;
+ Done = Standard_False;
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ Done = Standard_True;
+ ////modified by jgv, 31.08.2011////
+ F.Value(Sol, FF); //update F before GetStateNumber
+ ///////////////////////////////////
+ State = F.GetStateNumber();
+ }
return;
}
}
//fin du test solution
// Analyse de la progression...
- if (F2 < PreviousMinimum) {
+ //comparison of current minimum and previous minimum
+ if ((F2 - PreviousMinimum) <= aTol_Func){
if (Kount > 5) {
// L'historique est il bon ?
if (F2 >= 0.95*Save(Kount - 5)) {
if (!ChangeDirection) ChangeDirection = Standard_True;
- else return; // si un gain inf a 5% on sort
+ else
+ {
+ Done = Standard_False;
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ Done = Standard_True;
+ State = F.GetStateNumber();
+ }
+ return; // si un gain inf a 5% on sort
+ }
}
else ChangeDirection = Standard_False; //Si oui on recommence
}
Delta(i) = PreviousSolution(i) - Sol(i);
}
if (IsSolutionReached(F)) {
- Done = Standard_True;
- State = F.GetStateNumber();
+ Done = Standard_False;
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ Done = Standard_True;
+ ////modified by jgv, 31.08.2011////
+ F.Value(Sol, FF); //update F before GetStateNumber
+ ///////////////////////////////////
+ State = F.GetStateNumber();
+ }
return;
}
}
// 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();
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ State = F.GetStateNumber();
+ }
return;
}
}
- else return; // y a plus d'issues
+ else
+ {
+
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ State = F.GetStateNumber();
+ }
+ return; // y a plus d'issues
+ }
}
}
+ if (!theStopOnDivergent || !myIsDivergent)
+ {
+ State = F.GetStateNumber();
+ }
}