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