1 // Created on: 2014-01-20
2 // Created by: Alexaner Malyshev
3 // Copyright (c) 2014-2014 OPEN CASCADE SAS
5 // This file is part of Open CASCADE Technology software library.
7 // This library is free software; you can redistribute it and/or modify it under
8 // the terms of the GNU Lesser General Public License version 2.1 as published
9 // by the Free Software Foundation, with special exception defined in the file
10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
11 // distribution for complete text of the license and disclaimer of any warranty.
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement
16 #include <math_GlobOptMin.hxx>
18 #include <math_BFGS.hxx>
19 #include <math_Matrix.hxx>
20 #include <math_MultipleVarFunctionWithGradient.hxx>
21 #include <math_MultipleVarFunctionWithHessian.hxx>
22 #include <math_NewtonMinimum.hxx>
23 #include <math_Powell.hxx>
24 #include <Standard_Integer.hxx>
25 #include <Standard_Real.hxx>
27 const Handle(Standard_Type)& STANDARD_TYPE(math_GlobOptMin)
29 static Handle(Standard_Type) _atype = new Standard_Type ("math_GlobOptMin", sizeof (math_GlobOptMin));
33 //=======================================================================
34 //function : math_GlobOptMin
35 //purpose : Constructor
36 //=======================================================================
37 math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
38 const math_Vector& theA,
39 const math_Vector& theB,
40 const Standard_Real theC,
41 const Standard_Real theDiscretizationTol,
42 const Standard_Real theSameTol)
43 : myN(theFunc->NbVariables()),
60 for(i = 1; i <= myN; i++)
69 for(i = 1; i <= myN; i++)
71 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
74 myTol = theDiscretizationTol;
75 mySameTol = theSameTol;
77 myDone = Standard_False;
80 //=======================================================================
81 //function : SetGlobalParams
82 //purpose : Set params without memory allocation.
83 //=======================================================================
84 void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
85 const math_Vector& theA,
86 const math_Vector& theB,
87 const Standard_Real theC,
88 const Standard_Real theDiscretizationTol,
89 const Standard_Real theSameTol)
98 for(i = 1; i <= myN; i++)
100 myGlobA(i) = theA(i);
101 myGlobB(i) = theB(i);
107 myTol = theDiscretizationTol;
108 mySameTol = theSameTol;
110 myDone = Standard_False;
113 //=======================================================================
114 //function : SetLocalParams
115 //purpose : Set params without memory allocation.
116 //=======================================================================
117 void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
118 const math_Vector& theLocalB)
125 for(i = 1; i <= myN; i++)
127 myA(i) = theLocalA(i);
128 myB(i) = theLocalB(i);
131 for(i = 1; i <= myN; i++)
133 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
136 myDone = Standard_False;
139 //=======================================================================
141 //purpose : Set algorithm tolerances.
142 //=======================================================================
143 void math_GlobOptMin::SetTol(const Standard_Real theDiscretizationTol,
144 const Standard_Real theSameTol)
146 myTol = theDiscretizationTol;
147 mySameTol = theSameTol;
150 //=======================================================================
152 //purpose : Get algorithm tolerances.
153 //=======================================================================
154 void math_GlobOptMin::GetTol(Standard_Real& theDiscretizationTol,
155 Standard_Real& theSameTol)
157 theDiscretizationTol = myTol;
158 theSameTol = mySameTol;
161 //=======================================================================
162 //function : ~math_GlobOptMin
164 //=======================================================================
165 math_GlobOptMin::~math_GlobOptMin()
169 //=======================================================================
171 //purpose : Compute Global extremum point
172 //=======================================================================
173 // In this algo indexes started from 1, not from 0.
174 void math_GlobOptMin::Perform()
178 // Compute parameters range
179 Standard_Real minLength = RealLast();
180 Standard_Real maxLength = RealFirst();
181 for(i = 1; i <= myN; i++)
183 Standard_Real currentLength = myB(i) - myA(i);
184 if (currentLength < minLength)
185 minLength = currentLength;
186 if (currentLength > maxLength)
187 maxLength = currentLength;
190 myE1 = minLength * myTol / myC;
191 myE2 = maxLength * myTol * 2.0 / myC;
192 myE3 = - maxLength * myTol / 4.0;
194 // Compute start point.
195 math_Vector aPnt(1,myN);
196 for(i = 1; i <= myN; i++)
198 Standard_Real currCentral = (myA(i) + myB(i)) / 2.0;
199 aPnt(i) = currCentral;
202 myFunc->Value(aPnt, myF);
204 math_Vector aExtremumPoint(1,myN);
205 Standard_Real aExtremumValue = RealLast();
206 if (computeLocalExtremum(aPnt, aExtremumValue, aExtremumPoint))
208 // Local Extremum finds better solution than midpoint.
209 if (aExtremumValue < myF)
211 myF = aExtremumValue;
212 aPnt = aExtremumPoint;
217 for(i = 1; i <= myN; i++)
221 computeGlobalExtremum(myN);
223 myDone = Standard_True;
226 //=======================================================================
227 //function : computeLocalExtremum
229 //=======================================================================
230 Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt,
231 Standard_Real& theVal,
232 math_Vector& theOutPnt)
237 if (dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc))
239 math_MultipleVarFunctionWithHessian* myTmp =
240 dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc);
242 math_NewtonMinimum newtonMinimum(*myTmp, thePnt);
243 if (newtonMinimum.IsDone())
245 newtonMinimum.Location(theOutPnt);
246 theVal = newtonMinimum.Minimum();
248 else return Standard_False;
252 if (dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc))
254 math_MultipleVarFunctionWithGradient* myTmp =
255 dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
256 math_BFGS bfgs(*myTmp, thePnt);
259 bfgs.Location(theOutPnt);
260 theVal = bfgs.Minimum();
262 else return Standard_False;
265 // Powell method used.
266 if (dynamic_cast<math_MultipleVarFunction*>(myFunc))
268 math_Matrix m(1, myN, 1, myN, 0.0);
269 for(i = 1; i <= myN; i++)
272 math_Powell powell(*myFunc, thePnt, m, 1e-10);
276 powell.Location(theOutPnt);
277 theVal = powell.Minimum();
279 else return Standard_False;
282 if (isInside(theOutPnt))
283 return Standard_True;
285 return Standard_False;
288 //=======================================================================
289 //function : ComputeGlobalExtremum
291 //=======================================================================
292 void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
295 Standard_Real d; // Functional in moved point.
296 Standard_Real val = RealLast(); // Local extrema computed in moved point.
297 Standard_Real stepBestValue = RealLast();
298 Standard_Real realStep = RealLast();
299 math_Vector stepBestPoint(1, myN);
300 Standard_Boolean isInside = Standard_False;
304 for(myX(j) = myA(j) + myE1; myX(j) < myB(j) + myE1; myX(j) += myV(j))
311 isInside = Standard_False;
312 myFunc->Value(myX, d);
316 isInside = computeLocalExtremum(myX, val, myTmp);
318 stepBestValue = (isInside && (val < d))? val : d;
319 stepBestPoint = (isInside && (val < d))? myTmp : myX;
321 // Solutions are close to each other.
322 if (Abs(stepBestValue - myF) < mySameTol * 0.01)
324 if (!isStored(stepBestPoint))
326 if ((stepBestValue - myF) * myZ > 0.0)
328 for(i = 1; i <= myN; i++)
329 myY.Append(stepBestPoint(i));
334 // New best solution.
335 if ((stepBestValue - myF) * myZ > mySameTol * 0.01)
340 for(i = 1; i <= myN; i++)
341 myY.Append(stepBestPoint(i));
345 realStep = myE2 + Abs(myF - d) / myC;
346 myV(1) = Min(realStep, myMaxV(1));
350 myV(j) = RealLast() / 2.0;
351 computeGlobalExtremum(j - 1);
353 if ((j < myN) && (myV(j + 1) > realStep))
355 if (realStep > myMaxV(j + 1)) // Case of too big step.
356 myV(j + 1) = myMaxV(j + 1);
358 myV(j + 1) = realStep;
363 //=======================================================================
364 //function : IsInside
366 //=======================================================================
367 Standard_Boolean math_GlobOptMin::isInside(const math_Vector& thePnt)
371 for(i = 1; i <= myN; i++)
373 if (thePnt(i) < myGlobA(i) || thePnt(i) > myGlobB(i))
374 return Standard_False;
377 return Standard_True;
379 //=======================================================================
380 //function : IsStored
382 //=======================================================================
383 Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
385 Standard_Integer i,j;
386 Standard_Boolean isSame = Standard_True;
388 for(i = 0; i < mySolCount; i++)
390 isSame = Standard_True;
391 for(j = 1; j <= myN; j++)
393 if ((Abs(thePnt(j) - myY(i * myN + j))) > (myB(j) - myA(j)) * mySameTol)
395 isSame = Standard_False;
399 if (isSame == Standard_True)
400 return Standard_True;
403 return Standard_False;
406 //=======================================================================
407 //function : NbExtrema
409 //=======================================================================
410 Standard_Integer math_GlobOptMin::NbExtrema()
415 //=======================================================================
418 //=======================================================================
419 Standard_Real math_GlobOptMin::GetF()
424 //=======================================================================
427 //=======================================================================
428 Standard_Boolean math_GlobOptMin::isDone()
433 //=======================================================================
436 //=======================================================================
437 void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
441 for(j = 1; j <= myN; j++)
442 theSol(j) = myY((theIndex - 1) * myN + j);