0026484: BRepExtrema_DistShapeShape hangs
[occt.git] / src / math / math_GlobOptMin.cxx
CommitLineData
4bbaf12b 1// Created on: 2014-01-20
2// Created by: Alexaner Malyshev
4b65fc77 3// Copyright (c) 2014-2015 OPEN CASCADE SAS
4bbaf12b 4//
5// This file is part of Open CASCADE Technology software library.
6//
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.
12//
13// Alternatively, this file may be used under the terms of Open CASCADE
14// commercial license or contractual agreement
15
16#include <math_GlobOptMin.hxx>
17
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>
4bbaf12b 24#include <Standard_Integer.hxx>
25#include <Standard_Real.hxx>
e8746a26 26#include <Precision.hxx>
4bbaf12b 27
4bbaf12b 28
29//=======================================================================
30//function : math_GlobOptMin
31//purpose : Constructor
32//=======================================================================
33math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
34 const math_Vector& theA,
35 const math_Vector& theB,
5493d334 36 const Standard_Real theC,
37 const Standard_Real theDiscretizationTol,
38 const Standard_Real theSameTol)
4bbaf12b 39: myN(theFunc->NbVariables()),
40 myA(1, myN),
41 myB(1, myN),
42 myGlobA(1, myN),
43 myGlobB(1, myN),
44 myX(1, myN),
45 myTmp(1, myN),
5493d334 46 myV(1, myN),
3f733bb1 47 myMaxV(1, myN),
4b65fc77 48 myExpandCoeff(1, myN),
49 myCellSize(0, myN - 1),
50 myFilter(theFunc->NbVariables())
4bbaf12b 51{
52 Standard_Integer i;
53
54 myFunc = theFunc;
55 myC = theC;
78e7cada 56 myIsFindSingleSolution = Standard_False;
4bbaf12b 57 myZ = -1;
58 mySolCount = 0;
59
60 for(i = 1; i <= myN; i++)
61 {
62 myGlobA(i) = theA(i);
63 myGlobB(i) = theB(i);
64
65 myA(i) = theA(i);
66 myB(i) = theB(i);
67 }
68
5493d334 69 for(i = 1; i <= myN; i++)
70 {
71 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
72 }
73
3f733bb1 74 myExpandCoeff(1) = 1.0;
75 for(i = 2; i <= myN; i++)
76 {
77 myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
78 }
79
5493d334 80 myTol = theDiscretizationTol;
81 mySameTol = theSameTol;
82
4b65fc77 83 const Standard_Integer aMaxSquareSearchSol = 200;
84 Standard_Integer aSolNb = Standard_Integer(Pow(3.0, Standard_Real(myN)));
85 myMinCellFilterSol = Max(2 * aSolNb, aMaxSquareSearchSol);
86 initCellSize();
87
4bbaf12b 88 myDone = Standard_False;
89}
90
91//=======================================================================
92//function : SetGlobalParams
93//purpose : Set params without memory allocation.
94//=======================================================================
95void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
96 const math_Vector& theA,
97 const math_Vector& theB,
5493d334 98 const Standard_Real theC,
99 const Standard_Real theDiscretizationTol,
100 const Standard_Real theSameTol)
4bbaf12b 101{
102 Standard_Integer i;
103
104 myFunc = theFunc;
105 myC = theC;
106 myZ = -1;
107 mySolCount = 0;
108
109 for(i = 1; i <= myN; i++)
110 {
111 myGlobA(i) = theA(i);
112 myGlobB(i) = theB(i);
113
114 myA(i) = theA(i);
115 myB(i) = theB(i);
116 }
117
3f733bb1 118 for(i = 1; i <= myN; i++)
119 {
120 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
121 }
122
123 myExpandCoeff(1) = 1.0;
124 for(i = 2; i <= myN; i++)
125 {
126 myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
127 }
128
5493d334 129 myTol = theDiscretizationTol;
130 mySameTol = theSameTol;
131
4b65fc77 132 initCellSize();
133
4bbaf12b 134 myDone = Standard_False;
135}
136
137//=======================================================================
138//function : SetLocalParams
139//purpose : Set params without memory allocation.
140//=======================================================================
141void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
142 const math_Vector& theLocalB)
143{
144 Standard_Integer i;
145
146 myZ = -1;
147 mySolCount = 0;
148
149 for(i = 1; i <= myN; i++)
150 {
151 myA(i) = theLocalA(i);
152 myB(i) = theLocalB(i);
153 }
154
5493d334 155 for(i = 1; i <= myN; i++)
156 {
157 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
158 }
159
3f733bb1 160 myExpandCoeff(1) = 1.0;
161 for(i = 2; i <= myN; i++)
162 {
163 myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
164 }
165
4bbaf12b 166 myDone = Standard_False;
167}
168
169//=======================================================================
5493d334 170//function : SetTol
171//purpose : Set algorithm tolerances.
172//=======================================================================
173void math_GlobOptMin::SetTol(const Standard_Real theDiscretizationTol,
174 const Standard_Real theSameTol)
175{
176 myTol = theDiscretizationTol;
177 mySameTol = theSameTol;
178}
179
180//=======================================================================
181//function : GetTol
182//purpose : Get algorithm tolerances.
183//=======================================================================
184void math_GlobOptMin::GetTol(Standard_Real& theDiscretizationTol,
185 Standard_Real& theSameTol)
186{
187 theDiscretizationTol = myTol;
188 theSameTol = mySameTol;
189}
190
191//=======================================================================
4bbaf12b 192//function : ~math_GlobOptMin
193//purpose :
194//=======================================================================
195math_GlobOptMin::~math_GlobOptMin()
196{
197}
198
199//=======================================================================
200//function : Perform
201//purpose : Compute Global extremum point
202//=======================================================================
203// In this algo indexes started from 1, not from 0.
78e7cada 204void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution)
4bbaf12b 205{
206 Standard_Integer i;
207
208 // Compute parameters range
209 Standard_Real minLength = RealLast();
210 Standard_Real maxLength = RealFirst();
211 for(i = 1; i <= myN; i++)
212 {
213 Standard_Real currentLength = myB(i) - myA(i);
214 if (currentLength < minLength)
215 minLength = currentLength;
216 if (currentLength > maxLength)
217 maxLength = currentLength;
218 }
219
e8746a26 220 if (minLength < Precision::PConfusion())
221 {
222 #ifdef OCCT_DEBUG
223 cout << "math_GlobOptMin::Perform(): Degenerated parameters space" << endl;
224 #endif
225
226 return;
227 }
228
229 // Compute initial values for myF, myY, myC.
230 computeInitialValues();
231
797d11c6 232 myE1 = minLength * myTol;
233 myE2 = maxLength * myTol;
78e7cada 234
235 myIsFindSingleSolution = isFindSingleSolution;
236 if (isFindSingleSolution)
237 {
238 // Run local optimization
239 // if current value better than optimal.
240 myE3 = 0.0;
241 }
797d11c6 242 else
78e7cada 243 {
244 if (myC > 1.0)
245 myE3 = - maxLength * myTol / 4.0;
246 else
247 myE3 = - maxLength * myTol * myC / 4.0;
248 }
4bbaf12b 249
4b65fc77 250 isFirstCellFilterInvoke = Standard_True;
4bbaf12b 251 computeGlobalExtremum(myN);
252
253 myDone = Standard_True;
4bbaf12b 254}
255
256//=======================================================================
257//function : computeLocalExtremum
258//purpose :
259//=======================================================================
260Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt,
261 Standard_Real& theVal,
262 math_Vector& theOutPnt)
263{
264 Standard_Integer i;
265
266 //Newton method
267 if (dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc))
268 {
269 math_MultipleVarFunctionWithHessian* myTmp =
270 dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc);
859a47c3 271 math_NewtonMinimum newtonMinimum(*myTmp);
91806b90 272 newtonMinimum.SetBoundary(myGlobA, myGlobB);
859a47c3 273 newtonMinimum.Perform(*myTmp, thePnt);
274
4bbaf12b 275 if (newtonMinimum.IsDone())
276 {
277 newtonMinimum.Location(theOutPnt);
278 theVal = newtonMinimum.Minimum();
279 }
280 else return Standard_False;
281 } else
282
283 // BFGS method used.
284 if (dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc))
285 {
286 math_MultipleVarFunctionWithGradient* myTmp =
287 dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
07f1a2e6 288 math_BFGS bfgs(myTmp->NbVariables());
289 bfgs.Perform(*myTmp, thePnt);
4bbaf12b 290 if (bfgs.IsDone())
291 {
292 bfgs.Location(theOutPnt);
293 theVal = bfgs.Minimum();
294 }
295 else return Standard_False;
296 } else
297
298 // Powell method used.
299 if (dynamic_cast<math_MultipleVarFunction*>(myFunc))
300 {
301 math_Matrix m(1, myN, 1, myN, 0.0);
302 for(i = 1; i <= myN; i++)
303 m(1, 1) = 1.0;
304
859a47c3 305 math_Powell powell(*myFunc, 1e-10);
306 powell.Perform(*myFunc, thePnt, m);
4bbaf12b 307
308 if (powell.IsDone())
309 {
310 powell.Location(theOutPnt);
311 theVal = powell.Minimum();
312 }
313 else return Standard_False;
314 }
315
316 if (isInside(theOutPnt))
317 return Standard_True;
318 else
319 return Standard_False;
320}
321
322//=======================================================================
797d11c6 323//function : computeInitialValues
324//purpose :
325//=======================================================================
326void math_GlobOptMin::computeInitialValues()
327{
328 Standard_Integer i;
329 math_Vector aCurrPnt(1, myN);
330 math_Vector aBestPnt(1, myN);
e8746a26 331 math_Vector aParamStep(1, myN);
797d11c6 332 Standard_Real aCurrVal = RealLast();
333 Standard_Real aBestVal = RealLast();
334
335 // Check functional value in midpoint, low and upp point border and
336 // in each point try to perform local optimization.
337 aBestPnt = (myA + myB) * 0.5;
338 myFunc->Value(aBestPnt, aBestVal);
339
340 for(i = 1; i <= 3; i++)
341 {
342 aCurrPnt = myA + (myB - myA) * (i - 1) / 2.0;
343
344 if(computeLocalExtremum(aCurrPnt, aCurrVal, aCurrPnt))
345 {
346 // Local Extremum finds better solution than current point.
347 if (aCurrVal < aBestVal)
348 {
349 aBestVal = aCurrVal;
350 aBestPnt = aCurrPnt;
351 }
352 }
353 }
354
355 myF = aBestVal;
356 myY.Clear();
357 for(i = 1; i <= myN; i++)
358 myY.Append(aBestPnt(i));
359 mySolCount++;
360
361 // Lipschitz const approximation
e8746a26 362 Standard_Real aLipConst = 0.0, aPrevValDiag, aPrevValProj;
797d11c6 363 Standard_Integer aPntNb = 13;
e8746a26 364 myFunc->Value(myA, aPrevValDiag);
365 aPrevValProj = aPrevValDiag;
797d11c6 366 Standard_Real aStep = (myB - myA).Norm() / aPntNb;
e8746a26 367 aParamStep = (myB - myA) / aPntNb;
797d11c6 368 for(i = 1; i <= aPntNb; i++)
369 {
e8746a26 370 aCurrPnt = myA + aParamStep * i;
797d11c6 371
e8746a26 372 // Walk over diagonal.
373 myFunc->Value(aCurrPnt, aCurrVal);
374 aLipConst = Max (Abs(aCurrVal - aPrevValDiag), aLipConst);
375 aPrevValDiag = aCurrVal;
797d11c6 376
e8746a26 377 // Walk over diag in projected space aPnt(1) = myA(1) = const.
378 aCurrPnt(1) = myA(1);
379 myFunc->Value(aCurrPnt, aCurrVal);
380 aLipConst = Max (Abs(aCurrVal - aPrevValProj), aLipConst);
381 aPrevValProj = aCurrVal;
797d11c6 382 }
e8746a26 383
384 aLipConst *= Sqrt(myN) / aStep;
797d11c6 385
386 if (aLipConst < myC * 0.1)
387 {
388 myC = Max(aLipConst * 0.1, 0.01);
389 }
390 else if (aLipConst > myC * 10)
391 {
392 myC = Min(myC * 2, 30.0);
393 }
394}
395
396//=======================================================================
4bbaf12b 397//function : ComputeGlobalExtremum
398//purpose :
399//=======================================================================
400void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
401{
402 Standard_Integer i;
403 Standard_Real d; // Functional in moved point.
404 Standard_Real val = RealLast(); // Local extrema computed in moved point.
3f733bb1 405 Standard_Real aStepBestValue = RealLast();
406 Standard_Real aRealStep = 0.0;
407 math_Vector aStepBestPoint(1, myN);
4bbaf12b 408 Standard_Boolean isInside = Standard_False;
409 Standard_Real r;
debc95ee 410 Standard_Boolean isReached = Standard_False;
4bbaf12b 411
debc95ee 412 for(myX(j) = myA(j) + myE1;
413 (myX(j) < myB(j) + myE1) && (!isReached);
414 myX(j) += myV(j))
4bbaf12b 415 {
416 if (myX(j) > myB(j))
debc95ee 417 {
4bbaf12b 418 myX(j) = myB(j);
debc95ee 419 isReached = Standard_True;
420 }
4bbaf12b 421
422 if (j == 1)
423 {
424 isInside = Standard_False;
425 myFunc->Value(myX, d);
3f733bb1 426 r = (d + myZ * myC * aRealStep - myF) * myZ;
4bbaf12b 427 if(r > myE3)
428 {
429 isInside = computeLocalExtremum(myX, val, myTmp);
430 }
3f733bb1 431 aStepBestValue = (isInside && (val < d))? val : d;
432 aStepBestPoint = (isInside && (val < d))? myTmp : myX;
4bbaf12b 433
78e7cada 434 // Solutions are close to each other
435 // and it is allowed to have more than one solution.
436 if (Abs(aStepBestValue - myF) < mySameTol * 0.01 &&
437 !myIsFindSingleSolution)
4bbaf12b 438 {
3f733bb1 439 if (!isStored(aStepBestPoint))
4bbaf12b 440 {
3f733bb1 441 if ((aStepBestValue - myF) * myZ > 0.0)
442 myF = aStepBestValue;
4bbaf12b 443 for(i = 1; i <= myN; i++)
3f733bb1 444 myY.Append(aStepBestPoint(i));
4bbaf12b 445 mySolCount++;
446 }
447 }
448
78e7cada 449 // New best solution:
450 // new point is out of (mySameTol * 0.01) surrounding or
451 // new point is better than old + single point search.
452 Standard_Real aFunctionalDelta = (aStepBestValue - myF) * myZ;
453 if (aFunctionalDelta > mySameTol * 0.01 ||
454 (aFunctionalDelta > 0.0 && myIsFindSingleSolution))
4bbaf12b 455 {
456 mySolCount = 0;
3f733bb1 457 myF = aStepBestValue;
4bbaf12b 458 myY.Clear();
459 for(i = 1; i <= myN; i++)
3f733bb1 460 myY.Append(aStepBestPoint(i));
4bbaf12b 461 mySolCount++;
4b65fc77 462
463 isFirstCellFilterInvoke = Standard_True;
4bbaf12b 464 }
465
3f733bb1 466 aRealStep = myE2 + Abs(myF - d) / myC;
467 myV(1) = Min(aRealStep, myMaxV(1));
4bbaf12b 468 }
469 else
470 {
471 myV(j) = RealLast() / 2.0;
472 computeGlobalExtremum(j - 1);
3f733bb1 473
474 // Nullify steps on lower dimensions.
475 for(i = 1; i < j; i++)
476 myV(i) = 0.0;
4bbaf12b 477 }
3f733bb1 478 // Compute step in (j + 1) dimension according to scale.
479 if (j < myN)
4bbaf12b 480 {
3f733bb1 481 Standard_Real aUpperDimStep = myV(j) * myExpandCoeff(j + 1);
482 if (myV(j + 1) > aUpperDimStep)
483 {
484 if (aUpperDimStep > myMaxV(j + 1)) // Case of too big step.
485 myV(j + 1) = myMaxV(j + 1);
486 else
487 myV(j + 1) = aUpperDimStep;
488 }
4bbaf12b 489 }
490 }
491}
492
493//=======================================================================
494//function : IsInside
495//purpose :
496//=======================================================================
497Standard_Boolean math_GlobOptMin::isInside(const math_Vector& thePnt)
498{
499 Standard_Integer i;
500
501 for(i = 1; i <= myN; i++)
502 {
503 if (thePnt(i) < myGlobA(i) || thePnt(i) > myGlobB(i))
504 return Standard_False;
505 }
506
507 return Standard_True;
508}
509//=======================================================================
510//function : IsStored
511//purpose :
512//=======================================================================
513Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
514{
515 Standard_Integer i,j;
516 Standard_Boolean isSame = Standard_True;
20a216fe 517 math_Vector aTol(1, myN);
518 aTol = (myB - myA) * mySameTol;
4bbaf12b 519
4b65fc77 520 // C1 * n^2 = C2 * 3^dim * n
521 if (mySolCount < myMinCellFilterSol)
4bbaf12b 522 {
4b65fc77 523 for(i = 0; i < mySolCount; i++)
4bbaf12b 524 {
4b65fc77 525 isSame = Standard_True;
526 for(j = 1; j <= myN; j++)
4bbaf12b 527 {
4b65fc77 528 if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
529 {
530 isSame = Standard_False;
531 break;
532 }
4bbaf12b 533 }
4b65fc77 534 if (isSame == Standard_True)
535 return Standard_True;
4bbaf12b 536 }
4b65fc77 537 }
538 else
539 {
50bc8f96 540 NCollection_CellFilter_Inspector anInspector(myN, Precision::PConfusion());
4b65fc77 541 if (isFirstCellFilterInvoke)
542 {
543 myFilter.Reset(myCellSize);
4bbaf12b 544
4b65fc77 545 // Copy initial data into cell filter.
546 for(Standard_Integer aSolIdx = 0; aSolIdx < mySolCount; aSolIdx++)
547 {
548 math_Vector aVec(1, myN);
549 for(Standard_Integer aSolDim = 1; aSolDim <= myN; aSolDim++)
550 aVec(aSolDim) = myY(aSolIdx * myN + aSolDim);
551
552 myFilter.Add(aVec, aVec);
553 }
554 }
555
556 isFirstCellFilterInvoke = Standard_False;
557
558 math_Vector aLow(1, myN), anUp(1, myN);
559 anInspector.Shift(thePnt, myCellSize, aLow, anUp);
560
561 anInspector.ClearFind();
562 anInspector.SetCurrent(thePnt);
563 myFilter.Inspect(aLow, anUp, anInspector);
564 if (!anInspector.isFind())
565 {
566 // Point is out of close cells, add new one.
567 myFilter.Add(thePnt, thePnt);
568 }
4bbaf12b 569 }
570 return Standard_False;
571}
572
573//=======================================================================
574//function : NbExtrema
575//purpose :
576//=======================================================================
577Standard_Integer math_GlobOptMin::NbExtrema()
578{
579 return mySolCount;
580}
581
582//=======================================================================
583//function : GetF
584//purpose :
585//=======================================================================
586Standard_Real math_GlobOptMin::GetF()
587{
588 return myF;
589}
590
591//=======================================================================
592//function : IsDone
593//purpose :
594//=======================================================================
595Standard_Boolean math_GlobOptMin::isDone()
596{
597 return myDone;
598}
599
600//=======================================================================
601//function : Points
602//purpose :
603//=======================================================================
604void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
605{
606 Standard_Integer j;
607
608 for(j = 1; j <= myN; j++)
609 theSol(j) = myY((theIndex - 1) * myN + j);
610}
4b65fc77 611
612//=======================================================================
613//function : initCellSize
614//purpose :
615//=======================================================================
616void math_GlobOptMin::initCellSize()
617{
618 for(Standard_Integer anIdx = 1; anIdx <= myN; anIdx++)
619 {
620 myCellSize(anIdx - 1) = (myGlobB(anIdx) - myGlobA(anIdx))
621 * Precision::PConfusion() / (2.0 * Sqrt(2.0));
622 }
623}