mirror of
https://git.dev.opencascade.org/repos/occt.git
synced 2025-04-26 10:19:45 +03:00
a5a7b3185b
Update empty method guards to new style with regex (see PR). Used clang-format 18.1.8. New actions to validate code formatting is added. Update .clang-format with disabling of include sorting. It is temporary changes, then include will be sorted. Apply formatting for /src and /tools folder. The files with .hxx,.cxx,.lxx,.h,.pxx,.hpp,*.cpp extensions.
478 lines
15 KiB
C++
478 lines
15 KiB
C++
// Copyright (c) 1997-1999 Matra Datavision
|
|
// Copyright (c) 1999-2014 OPEN CASCADE SAS
|
|
//
|
|
// This file is part of Open CASCADE Technology software library.
|
|
//
|
|
// This library is free software; you can redistribute it and/or modify it under
|
|
// the terms of the GNU Lesser General Public License version 2.1 as published
|
|
// by the Free Software Foundation, with special exception defined in the file
|
|
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
|
|
// distribution for complete text of the license and disclaimer of any warranty.
|
|
//
|
|
// Alternatively, this file may be used under the terms of Open CASCADE
|
|
// commercial license or contractual agreement.
|
|
|
|
// #ifndef OCCT_DEBUG
|
|
#define No_Standard_RangeError
|
|
#define No_Standard_OutOfRange
|
|
#define No_Standard_DimensionError
|
|
|
|
// #endif
|
|
|
|
#include <math_BFGS.hxx>
|
|
#include <math_BracketMinimum.hxx>
|
|
#include <math_BrentMinimum.hxx>
|
|
#include <math_FunctionWithDerivative.hxx>
|
|
#include <math_Matrix.hxx>
|
|
#include <math_MultipleVarFunctionWithGradient.hxx>
|
|
#include <Precision.hxx>
|
|
|
|
// l'utilisation de math_BrentMinumim pur trouver un minimum dans une direction
|
|
// donnee n'est pas du tout optimale. voir peut etre interpolation cubique
|
|
// classique et aussi essayer "recherche unidimensionnelle economique"
|
|
// PROGRAMMATION MATHEMATIQUE (theorie et algorithmes) tome1 page 82.
|
|
|
|
// Target function for 1D problem, point and direction are known.
|
|
class DirFunction : public math_FunctionWithDerivative
|
|
{
|
|
|
|
math_Vector* P0;
|
|
math_Vector* Dir;
|
|
math_Vector* P;
|
|
math_Vector* G;
|
|
math_MultipleVarFunctionWithGradient* F;
|
|
|
|
public:
|
|
//! Ctor.
|
|
DirFunction(math_Vector& V1,
|
|
math_Vector& V2,
|
|
math_Vector& V3,
|
|
math_Vector& V4,
|
|
math_MultipleVarFunctionWithGradient& f)
|
|
: P0(&V1),
|
|
Dir(&V2),
|
|
P(&V3),
|
|
G(&V4),
|
|
F(&f)
|
|
{
|
|
}
|
|
|
|
//! Sets point and direction.
|
|
void Initialize(const math_Vector& p0, const math_Vector& dir) const
|
|
{
|
|
*P0 = p0;
|
|
*Dir = dir;
|
|
}
|
|
|
|
void TheGradient(math_Vector& Grad) { Grad = *G; }
|
|
|
|
virtual Standard_Boolean Value(const Standard_Real x, Standard_Real& fval)
|
|
{
|
|
*P = *Dir;
|
|
P->Multiply(x);
|
|
P->Add(*P0);
|
|
fval = 0.;
|
|
return F->Value(*P, fval);
|
|
}
|
|
|
|
virtual Standard_Boolean Values(const Standard_Real x, Standard_Real& fval, Standard_Real& D)
|
|
{
|
|
*P = *Dir;
|
|
P->Multiply(x);
|
|
P->Add(*P0);
|
|
fval = D = 0.;
|
|
if (F->Values(*P, fval, *G))
|
|
{
|
|
D = (*G).Multiplied(*Dir);
|
|
return Standard_True;
|
|
}
|
|
|
|
return Standard_False;
|
|
}
|
|
|
|
virtual Standard_Boolean Derivative(const Standard_Real x, Standard_Real& D)
|
|
{
|
|
*P = *Dir;
|
|
P->Multiply(x);
|
|
P->Add(*P0);
|
|
Standard_Real fval;
|
|
D = 0.;
|
|
if (F->Values(*P, fval, *G))
|
|
{
|
|
D = (*G).Multiplied(*Dir);
|
|
return Standard_True;
|
|
}
|
|
|
|
return Standard_False;
|
|
}
|
|
};
|
|
|
|
//=============================================================================
|
|
// function : ComputeInitScale
|
|
// purpose : Compute the appropriate initial value of scale factor to apply
|
|
// to the direction to approach to the minimum of the function
|
|
//=============================================================================
|
|
static Standard_Boolean ComputeInitScale(const Standard_Real theF0,
|
|
const math_Vector& theDir,
|
|
const math_Vector& theGr,
|
|
Standard_Real& theScale)
|
|
{
|
|
const Standard_Real dy1 = theGr * theDir;
|
|
if (Abs(dy1) < RealSmall())
|
|
return Standard_False;
|
|
|
|
const Standard_Real aHnr1 = theDir.Norm2();
|
|
const Standard_Real alfa = 0.7 * (-theF0) / dy1;
|
|
theScale = 0.015 / Sqrt(aHnr1);
|
|
if (theScale > alfa)
|
|
theScale = alfa;
|
|
|
|
return Standard_True;
|
|
}
|
|
|
|
//=============================================================================
|
|
// function : ComputeMinMaxScale
|
|
// purpose : For a given point and direction, and bounding box,
|
|
// find min and max scale factors with which the point reaches borders
|
|
// if we apply translation Point+Dir*Scale.
|
|
// return : True if found, False if point is out of bounds.
|
|
//=============================================================================
|
|
static Standard_Boolean ComputeMinMaxScale(const math_Vector& thePoint,
|
|
const math_Vector& theDir,
|
|
const math_Vector& theLeft,
|
|
const math_Vector& theRight,
|
|
Standard_Real& theMinScale,
|
|
Standard_Real& theMaxScale)
|
|
{
|
|
for (Standard_Integer anIdx = 1; anIdx <= theLeft.Upper(); anIdx++)
|
|
{
|
|
const Standard_Real aLeft = theLeft(anIdx) - thePoint(anIdx);
|
|
const Standard_Real aRight = theRight(anIdx) - thePoint(anIdx);
|
|
if (Abs(theDir(anIdx)) > RealSmall())
|
|
{
|
|
// Use PConfusion to get off a little from the bounds to prevent
|
|
// possible refuse in Value function.
|
|
const Standard_Real aLScale = (aLeft + Precision::PConfusion()) / theDir(anIdx);
|
|
const Standard_Real aRScale = (aRight - Precision::PConfusion()) / theDir(anIdx);
|
|
if (Abs(aLeft) < Precision::PConfusion())
|
|
{
|
|
// Point is on the left border.
|
|
theMaxScale = Min(theMaxScale, Max(0., aRScale));
|
|
theMinScale = Max(theMinScale, Min(0., aRScale));
|
|
}
|
|
else if (Abs(aRight) < Precision::PConfusion())
|
|
{
|
|
// Point is on the right border.
|
|
theMaxScale = Min(theMaxScale, Max(0., aLScale));
|
|
theMinScale = Max(theMinScale, Min(0., aLScale));
|
|
}
|
|
else if (aLeft * aRight < 0)
|
|
{
|
|
// Point is inside allowed range.
|
|
theMaxScale = Min(theMaxScale, Max(aLScale, aRScale));
|
|
theMinScale = Max(theMinScale, Min(aLScale, aRScale));
|
|
}
|
|
else
|
|
// point is out of bounds
|
|
return Standard_False;
|
|
}
|
|
else
|
|
{
|
|
// Direction is parallel to the border.
|
|
// Check that the point is not out of bounds
|
|
if (aLeft > Precision::PConfusion() || aRight < -Precision::PConfusion())
|
|
{
|
|
return Standard_False;
|
|
}
|
|
}
|
|
}
|
|
return Standard_True;
|
|
}
|
|
|
|
//=============================================================================
|
|
// function : MinimizeDirection
|
|
// purpose : Solves 1D minimization problem when point and directions
|
|
// are known.
|
|
//=============================================================================
|
|
static Standard_Boolean MinimizeDirection(math_Vector& P,
|
|
Standard_Real F0,
|
|
math_Vector& Gr,
|
|
math_Vector& Dir,
|
|
Standard_Real& Result,
|
|
DirFunction& F,
|
|
Standard_Boolean isBounds,
|
|
const math_Vector& theLeft,
|
|
const math_Vector& theRight)
|
|
{
|
|
Standard_Real lambda;
|
|
if (!ComputeInitScale(F0, Dir, Gr, lambda))
|
|
return Standard_False;
|
|
|
|
// by default the scaling range is unlimited
|
|
Standard_Real aMinLambda = -Precision::Infinite();
|
|
Standard_Real aMaxLambda = Precision::Infinite();
|
|
if (isBounds)
|
|
{
|
|
// limit the scaling range taking into account the bounds
|
|
if (!ComputeMinMaxScale(P, Dir, theLeft, theRight, aMinLambda, aMaxLambda))
|
|
return Standard_False;
|
|
|
|
if (aMinLambda > -Precision::PConfusion() && aMaxLambda < Precision::PConfusion())
|
|
{
|
|
// Point is on the border and the direction shows outside.
|
|
// Make direction to go along the border
|
|
for (Standard_Integer anIdx = 1; anIdx <= theLeft.Upper(); anIdx++)
|
|
{
|
|
if ((Abs(P(anIdx) - theRight(anIdx)) < Precision::PConfusion() && Dir(anIdx) > 0.0)
|
|
|| (Abs(P(anIdx) - theLeft(anIdx)) < Precision::PConfusion() && Dir(anIdx) < 0.0))
|
|
{
|
|
Dir(anIdx) = 0.0;
|
|
}
|
|
}
|
|
|
|
// re-compute scale values with new direction
|
|
if (!ComputeInitScale(F0, Dir, Gr, lambda))
|
|
return Standard_False;
|
|
if (!ComputeMinMaxScale(P, Dir, theLeft, theRight, aMinLambda, aMaxLambda))
|
|
return Standard_False;
|
|
}
|
|
lambda = Min(lambda, aMaxLambda);
|
|
lambda = Max(lambda, aMinLambda);
|
|
}
|
|
|
|
F.Initialize(P, Dir);
|
|
Standard_Real F1;
|
|
if (!F.Value(lambda, F1))
|
|
return Standard_False;
|
|
|
|
math_BracketMinimum Bracket(0.0, lambda);
|
|
if (isBounds)
|
|
Bracket.SetLimits(aMinLambda, aMaxLambda);
|
|
Bracket.SetFA(F0);
|
|
Bracket.SetFB(F1);
|
|
Bracket.Perform(F);
|
|
if (Bracket.IsDone())
|
|
{
|
|
// find minimum inside the bracket
|
|
Standard_Real ax, xx, bx, Fax, Fxx, Fbx;
|
|
Bracket.Values(ax, xx, bx);
|
|
Bracket.FunctionValues(Fax, Fxx, Fbx);
|
|
|
|
Standard_Integer niter = 100;
|
|
Standard_Real tol = 1.e-03;
|
|
math_BrentMinimum Sol(tol, Fxx, niter, 1.e-08);
|
|
Sol.Perform(F, ax, xx, bx);
|
|
if (Sol.IsDone())
|
|
{
|
|
Standard_Real Scale = Sol.Location();
|
|
Result = Sol.Minimum();
|
|
Dir.Multiply(Scale);
|
|
P.Add(Dir);
|
|
return Standard_True;
|
|
}
|
|
}
|
|
else if (isBounds)
|
|
{
|
|
// Bracket definition is failure. If the bounds are defined then
|
|
// set current point to intersection with bounds
|
|
Standard_Real aFMin, aFMax;
|
|
if (!F.Value(aMinLambda, aFMin))
|
|
return Standard_False;
|
|
if (!F.Value(aMaxLambda, aFMax))
|
|
return Standard_False;
|
|
Standard_Real aBestLambda;
|
|
if (aFMin < aFMax)
|
|
{
|
|
aBestLambda = aMinLambda;
|
|
Result = aFMin;
|
|
}
|
|
else
|
|
{
|
|
aBestLambda = aMaxLambda;
|
|
Result = aFMax;
|
|
}
|
|
Dir.Multiply(aBestLambda);
|
|
P.Add(Dir);
|
|
return Standard_True;
|
|
}
|
|
return Standard_False;
|
|
}
|
|
|
|
//=============================================================================
|
|
// function : Perform
|
|
// purpose : Performs minimization problem using BFGS method.
|
|
//=============================================================================
|
|
void math_BFGS::Perform(math_MultipleVarFunctionWithGradient& F, const math_Vector& StartingPoint)
|
|
{
|
|
const Standard_Integer n = TheLocation.Length();
|
|
Standard_Boolean Good = Standard_True;
|
|
Standard_Integer j, i;
|
|
Standard_Real fae, fad, fac;
|
|
|
|
math_Vector xi(1, n), dg(1, n), hdg(1, n);
|
|
math_Matrix hessin(1, n, 1, n);
|
|
hessin.Init(0.0);
|
|
|
|
math_Vector Temp1(1, n);
|
|
math_Vector Temp2(1, n);
|
|
math_Vector Temp3(1, n);
|
|
math_Vector Temp4(1, n);
|
|
DirFunction F_Dir(Temp1, Temp2, Temp3, Temp4, F);
|
|
|
|
TheLocation = StartingPoint;
|
|
Good = F.Values(TheLocation, PreviousMinimum, TheGradient);
|
|
if (!Good)
|
|
{
|
|
Done = Standard_False;
|
|
TheStatus = math_FunctionError;
|
|
return;
|
|
}
|
|
for (i = 1; i <= n; i++)
|
|
{
|
|
hessin(i, i) = 1.0;
|
|
xi(i) = -TheGradient(i);
|
|
}
|
|
|
|
for (nbiter = 1; nbiter <= Itermax; nbiter++)
|
|
{
|
|
TheMinimum = PreviousMinimum;
|
|
const Standard_Boolean IsGood = MinimizeDirection(TheLocation,
|
|
TheMinimum,
|
|
TheGradient,
|
|
xi,
|
|
TheMinimum,
|
|
F_Dir,
|
|
myIsBoundsDefined,
|
|
myLeft,
|
|
myRight);
|
|
|
|
if (IsSolutionReached(F))
|
|
{
|
|
Done = Standard_True;
|
|
TheStatus = math_OK;
|
|
return;
|
|
}
|
|
|
|
if (!IsGood)
|
|
{
|
|
Done = Standard_False;
|
|
TheStatus = math_DirectionSearchError;
|
|
return;
|
|
}
|
|
PreviousMinimum = TheMinimum;
|
|
|
|
dg = TheGradient;
|
|
|
|
Good = F.Values(TheLocation, TheMinimum, TheGradient);
|
|
if (!Good)
|
|
{
|
|
Done = Standard_False;
|
|
TheStatus = math_FunctionError;
|
|
return;
|
|
}
|
|
|
|
for (i = 1; i <= n; i++)
|
|
dg(i) = TheGradient(i) - dg(i);
|
|
|
|
for (i = 1; i <= n; i++)
|
|
{
|
|
hdg(i) = 0.0;
|
|
for (j = 1; j <= n; j++)
|
|
hdg(i) += hessin(i, j) * dg(j);
|
|
}
|
|
|
|
fac = fae = 0.0;
|
|
for (i = 1; i <= n; i++)
|
|
{
|
|
fac += dg(i) * xi(i);
|
|
fae += dg(i) * hdg(i);
|
|
}
|
|
fac = 1.0 / fac;
|
|
fad = 1.0 / fae;
|
|
|
|
for (i = 1; i <= n; i++)
|
|
dg(i) = fac * xi(i) - fad * hdg(i);
|
|
|
|
for (i = 1; i <= n; i++)
|
|
{
|
|
for (j = 1; j <= n; j++)
|
|
hessin(i, j) += fac * xi(i) * xi(j) - fad * hdg(i) * hdg(j) + fae * dg(i) * dg(j);
|
|
}
|
|
|
|
for (i = 1; i <= n; i++)
|
|
{
|
|
xi(i) = 0.0;
|
|
for (j = 1; j <= n; j++)
|
|
xi(i) -= hessin(i, j) * TheGradient(j);
|
|
}
|
|
}
|
|
Done = Standard_False;
|
|
TheStatus = math_TooManyIterations;
|
|
return;
|
|
}
|
|
|
|
//=============================================================================
|
|
// function : IsSolutionReached
|
|
// purpose : Checks whether solution reached or not.
|
|
//=============================================================================
|
|
Standard_Boolean math_BFGS::IsSolutionReached(math_MultipleVarFunctionWithGradient&) const
|
|
{
|
|
|
|
return 2.0 * fabs(TheMinimum - PreviousMinimum)
|
|
<= XTol * (fabs(TheMinimum) + fabs(PreviousMinimum) + EPSZ);
|
|
}
|
|
|
|
//=================================================================================================
|
|
|
|
math_BFGS::math_BFGS(const Standard_Integer NbVariables,
|
|
const Standard_Real Tolerance,
|
|
const Standard_Integer NbIterations,
|
|
const Standard_Real ZEPS)
|
|
: TheStatus(math_OK),
|
|
TheLocation(1, NbVariables),
|
|
TheGradient(1, NbVariables),
|
|
PreviousMinimum(0.),
|
|
TheMinimum(0.),
|
|
XTol(Tolerance),
|
|
EPSZ(ZEPS),
|
|
nbiter(0),
|
|
myIsBoundsDefined(Standard_False),
|
|
myLeft(1, NbVariables, 0.0),
|
|
myRight(1, NbVariables, 0.0),
|
|
Done(Standard_False),
|
|
Itermax(NbIterations)
|
|
{
|
|
}
|
|
|
|
//=================================================================================================
|
|
|
|
math_BFGS::~math_BFGS() {}
|
|
|
|
//=================================================================================================
|
|
|
|
void math_BFGS::Dump(Standard_OStream& o) const
|
|
{
|
|
|
|
o << "math_BFGS resolution: ";
|
|
if (Done)
|
|
{
|
|
o << " Status = Done \n";
|
|
o << " Location Vector = " << Location() << "\n";
|
|
o << " Minimum value = " << Minimum() << "\n";
|
|
o << " Number of iterations = " << NbIterations() << "\n";
|
|
}
|
|
else
|
|
o << " Status = not Done because " << (Standard_Integer)TheStatus << "\n";
|
|
}
|
|
|
|
//=============================================================================
|
|
// function : SetBoundary
|
|
// purpose : Set boundaries for conditional optimization
|
|
//=============================================================================
|
|
void math_BFGS::SetBoundary(const math_Vector& theLeftBorder, const math_Vector& theRightBorder)
|
|
{
|
|
myLeft = theLeftBorder;
|
|
myRight = theRightBorder;
|
|
myIsBoundsDefined = Standard_True;
|
|
}
|