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.
16 #define No_Standard_RangeError
17 #define No_Standard_OutOfRange
18 #define No_Standard_DimensionError
22 #include <math_BFGS.hxx>
23 #include <math_BracketMinimum.hxx>
24 #include <math_BrentMinimum.hxx>
25 #include <math_FunctionWithDerivative.hxx>
26 #include <math_Matrix.hxx>
27 #include <math_MultipleVarFunctionWithGradient.hxx>
28 #include <Precision.hxx>
30 // l'utilisation de math_BrentMinumim pur trouver un minimum dans une direction
31 // donnee n'est pas du tout optimale. voir peut etre interpolation cubique
32 // classique et aussi essayer "recherche unidimensionnelle economique"
33 // PROGRAMMATION MATHEMATIQUE (theorie et algorithmes) tome1 page 82.
35 // Target function for 1D problem, point and direction are known.
36 class DirFunction : public math_FunctionWithDerivative
43 math_MultipleVarFunctionWithGradient *F;
48 DirFunction(math_Vector& V1,
52 math_MultipleVarFunctionWithGradient& f)
60 //! Sets point and direction.
61 void Initialize(const math_Vector& p0,
62 const math_Vector& dir) const
68 void TheGradient(math_Vector& Grad)
73 virtual Standard_Boolean Value(const Standard_Real x,
80 return F->Value(*P, fval);
83 virtual Standard_Boolean Values(const Standard_Real x,
91 if (F->Values(*P, fval, *G))
93 D = (*G).Multiplied(*Dir);
97 return Standard_False;
99 virtual Standard_Boolean Derivative(const Standard_Real x,
107 if (F->Values(*P, fval, *G))
109 D = (*G).Multiplied(*Dir);
110 return Standard_True;
113 return Standard_False;
119 //=============================================================================
120 //function : ComputeInitScale
121 //purpose : Compute the appropriate initial value of scale factor to apply
122 // to the direction to approach to the minimum of the function
123 //=============================================================================
124 static Standard_Boolean ComputeInitScale(const Standard_Real theF0,
125 const math_Vector& theDir,
126 const math_Vector& theGr,
127 Standard_Real& theScale)
129 const Standard_Real dy1 = theGr * theDir;
130 if (Abs(dy1) < RealSmall())
131 return Standard_False;
133 const Standard_Real aHnr1 = theDir.Norm2();
134 const Standard_Real alfa = 0.7*(-theF0) / dy1;
135 theScale = 0.015 / Sqrt(aHnr1);
139 return Standard_True;
142 //=============================================================================
143 //function : ComputeMinMaxScale
144 //purpose : For a given point and direction, and bounding box,
145 // find min and max scale factors with which the point reaches borders
146 // if we apply translation Point+Dir*Scale.
147 //return : True if found, False if point is out of bounds.
148 //=============================================================================
149 static Standard_Boolean ComputeMinMaxScale(const math_Vector& thePoint,
150 const math_Vector& theDir,
151 const math_Vector& theLeft,
152 const math_Vector& theRight,
153 Standard_Real& theMinScale,
154 Standard_Real& theMaxScale)
156 for (Standard_Integer anIdx = 1; anIdx <= theLeft.Upper(); anIdx++)
158 const Standard_Real aLeft = theLeft(anIdx) - thePoint(anIdx);
159 const Standard_Real aRight = theRight(anIdx) - thePoint(anIdx);
160 if (Abs(theDir(anIdx)) > RealSmall())
162 // Use PConfusion to get off a little from the bounds to prevent
163 // possible refuse in Value function.
164 const Standard_Real aLScale = (aLeft + Precision::PConfusion()) / theDir(anIdx);
165 const Standard_Real aRScale = (aRight - Precision::PConfusion()) / theDir(anIdx);
166 if (Abs(aLeft) < Precision::PConfusion())
168 // Point is on the left border.
169 theMaxScale = Min(theMaxScale, Max(0., aRScale));
170 theMinScale = Max(theMinScale, Min(0., aRScale));
172 else if (Abs(aRight) < Precision::PConfusion())
174 // Point is on the right border.
175 theMaxScale = Min(theMaxScale, Max(0., aLScale));
176 theMinScale = Max(theMinScale, Min(0., aLScale));
178 else if (aLeft * aRight < 0)
180 // Point is inside allowed range.
181 theMaxScale = Min(theMaxScale, Max(aLScale, aRScale));
182 theMinScale = Max(theMinScale, Min(aLScale, aRScale));
185 // point is out of bounds
186 return Standard_False;
190 // Direction is parallel to the border.
191 // Check that the point is not out of bounds
192 if (aLeft > Precision::PConfusion() ||
193 aRight < -Precision::PConfusion())
195 return Standard_False;
199 return Standard_True;
202 //=============================================================================
203 //function : MinimizeDirection
204 //purpose : Solves 1D minimization problem when point and directions
206 //=============================================================================
207 static Standard_Boolean MinimizeDirection(math_Vector& P,
211 Standard_Real& Result,
213 Standard_Boolean isBounds,
214 const math_Vector& theLeft,
215 const math_Vector& theRight)
217 Standard_Real lambda;
218 if (!ComputeInitScale(F0, Dir, Gr, lambda))
219 return Standard_False;
221 // by default the scaling range is unlimited
222 Standard_Real aMinLambda = -Precision::Infinite();
223 Standard_Real aMaxLambda = Precision::Infinite();
226 // limit the scaling range taking into account the bounds
227 if (!ComputeMinMaxScale(P, Dir, theLeft, theRight, aMinLambda, aMaxLambda))
228 return Standard_False;
230 if (aMinLambda > -Precision::PConfusion() && aMaxLambda < Precision::PConfusion())
232 // Point is on the border and the direction shows outside.
233 // Make direction to go along the border
234 for (Standard_Integer anIdx = 1; anIdx <= theLeft.Upper(); anIdx++)
236 if ((Abs(P(anIdx) - theRight(anIdx)) < Precision::PConfusion() && Dir(anIdx) > 0.0) ||
237 (Abs(P(anIdx) - theLeft(anIdx)) < Precision::PConfusion() && Dir(anIdx) < 0.0))
243 // re-compute scale values with new direction
244 if (!ComputeInitScale(F0, Dir, Gr, lambda))
245 return Standard_False;
246 if (!ComputeMinMaxScale(P, Dir, theLeft, theRight, aMinLambda, aMaxLambda))
247 return Standard_False;
249 lambda = Min(lambda, aMaxLambda);
250 lambda = Max(lambda, aMinLambda);
253 F.Initialize(P, Dir);
255 if (!F.Value(lambda, F1))
256 return Standard_False;
258 math_BracketMinimum Bracket(0.0, lambda);
260 Bracket.SetLimits(aMinLambda, aMaxLambda);
264 if (Bracket.IsDone())
266 // find minimum inside the bracket
267 Standard_Real ax, xx, bx, Fax, Fxx, Fbx;
268 Bracket.Values(ax, xx, bx);
269 Bracket.FunctionValues(Fax, Fxx, Fbx);
271 Standard_Integer niter = 100;
272 Standard_Real tol = 1.e-03;
273 math_BrentMinimum Sol(tol, Fxx, niter, 1.e-08);
274 Sol.Perform(F, ax, xx, bx);
277 Standard_Real Scale = Sol.Location();
278 Result = Sol.Minimum();
281 return Standard_True;
286 // Bracket definition is failure. If the bounds are defined then
287 // set current point to intersection with bounds
288 Standard_Real aFMin, aFMax;
289 if (!F.Value(aMinLambda, aFMin))
290 return Standard_False;
291 if (!F.Value(aMaxLambda, aFMax))
292 return Standard_False;
293 Standard_Real aBestLambda;
296 aBestLambda = aMinLambda;
301 aBestLambda = aMaxLambda;
304 Dir.Multiply(aBestLambda);
306 return Standard_True;
308 return Standard_False;
311 //=============================================================================
313 //purpose : Performs minimization problem using BFGS method.
314 //=============================================================================
315 void math_BFGS::Perform(math_MultipleVarFunctionWithGradient& F,
316 const math_Vector& StartingPoint)
318 const Standard_Integer n = TheLocation.Length();
319 Standard_Boolean Good = Standard_True;
320 Standard_Integer j, i;
321 Standard_Real fae, fad, fac;
323 math_Vector xi(1, n), dg(1, n), hdg(1, n);
324 math_Matrix hessin(1, n, 1, n);
327 math_Vector Temp1(1, n);
328 math_Vector Temp2(1, n);
329 math_Vector Temp3(1, n);
330 math_Vector Temp4(1, n);
331 DirFunction F_Dir(Temp1, Temp2, Temp3, Temp4, F);
333 TheLocation = StartingPoint;
334 Good = F.Values(TheLocation, PreviousMinimum, TheGradient);
337 Done = Standard_False;
338 TheStatus = math_FunctionError;
341 for (i = 1; i <= n; i++)
344 xi(i) = -TheGradient(i);
348 for (nbiter = 1; nbiter <= Itermax; nbiter++)
350 TheMinimum = PreviousMinimum;
351 const Standard_Boolean IsGood = MinimizeDirection(TheLocation, TheMinimum, TheGradient,
352 xi, TheMinimum, F_Dir, myIsBoundsDefined,
355 if (IsSolutionReached(F))
357 Done = Standard_True;
364 Done = Standard_False;
365 TheStatus = math_DirectionSearchError;
368 PreviousMinimum = TheMinimum;
372 Good = F.Values(TheLocation, TheMinimum, TheGradient);
375 Done = Standard_False;
376 TheStatus = math_FunctionError;
380 for (i = 1; i <= n; i++)
381 dg(i) = TheGradient(i) - dg(i);
383 for (i = 1; i <= n; i++)
386 for (j = 1; j <= n; j++)
387 hdg(i) += hessin(i, j) * dg(j);
391 for (i = 1; i <= n; i++)
393 fac += dg(i) * xi(i);
394 fae += dg(i) * hdg(i);
399 for (i = 1; i <= n; i++)
400 dg(i) = fac * xi(i) - fad * hdg(i);
402 for (i = 1; i <= n; i++)
404 for (j = 1; j <= n; j++)
405 hessin(i, j) += fac * xi(i) * xi(j) - fad * hdg(i) * hdg(j) + fae * dg(i) * dg(j);
408 for (i = 1; i <= n; i++)
411 for (j = 1; j <= n; j++)
412 xi(i) -= hessin(i, j) * TheGradient(j);
415 Done = Standard_False;
416 TheStatus = math_TooManyIterations;
420 //=============================================================================
421 //function : IsSolutionReached
422 //purpose : Checks whether solution reached or not.
423 //=============================================================================
424 Standard_Boolean math_BFGS::IsSolutionReached(math_MultipleVarFunctionWithGradient&) const
427 return 2.0 * fabs(TheMinimum - PreviousMinimum) <=
428 XTol * (fabs(TheMinimum) + fabs(PreviousMinimum) + EPSZ);
431 //=============================================================================
432 //function : math_BFGS
433 //purpose : Constructor.
434 //=============================================================================
435 math_BFGS::math_BFGS(const Standard_Integer NbVariables,
436 const Standard_Real Tolerance,
437 const Standard_Integer NbIterations,
438 const Standard_Real ZEPS)
439 : TheStatus(math_OK),
440 TheLocation(1, NbVariables),
441 TheGradient(1, NbVariables),
447 myIsBoundsDefined(Standard_False),
448 myLeft(1, NbVariables, 0.0),
449 myRight(1, NbVariables, 0.0),
450 Done(Standard_False),
451 Itermax(NbIterations)
455 //=============================================================================
456 //function : ~math_BFGS
457 //purpose : Destructor.
458 //=============================================================================
459 math_BFGS::~math_BFGS()
463 //=============================================================================
465 //purpose : Prints dump.
466 //=============================================================================
467 void math_BFGS::Dump(Standard_OStream& o) const
470 o << "math_BFGS resolution: ";
473 o << " Status = Done \n";
474 o << " Location Vector = " << Location() << "\n";
475 o << " Minimum value = " << Minimum() << "\n";
476 o << " Number of iterations = " << NbIterations() << "\n";
479 o << " Status = not Done because " << (Standard_Integer)TheStatus << "\n";
482 //=============================================================================
483 //function : SetBoundary
484 //purpose : Set boundaries for conditional optimization
485 //=============================================================================
486 void math_BFGS::SetBoundary(const math_Vector& theLeftBorder,
487 const math_Vector& theRightBorder)
489 myLeft = theLeftBorder;
490 myRight = theRightBorder;
491 myIsBoundsDefined = Standard_True;