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