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