99e968d3d2660586b0d0a98bacf33d10f8d29382
[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   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
69   for(i = 1; i <= myN; i++)
70   {
71     myMaxV(i) = (myB(i) - myA(i)) / 3.0;
72   }
73
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
80   myTol = theDiscretizationTol;
81   mySameTol = theSameTol;
82
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
88   myDone = Standard_False;
89 }
90
91 //=======================================================================
92 //function : SetGlobalParams
93 //purpose  : Set params without memory allocation.
94 //=======================================================================
95 void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
96                                       const math_Vector& theA,
97                                       const math_Vector& theB,
98                                       const Standard_Real theC,
99                                       const Standard_Real theDiscretizationTol,
100                                       const Standard_Real theSameTol)
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
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
129   myTol = theDiscretizationTol;
130   mySameTol = theSameTol;
131
132   initCellSize();
133
134   myDone = Standard_False;
135 }
136
137 //=======================================================================
138 //function : SetLocalParams
139 //purpose  : Set params without memory allocation.
140 //=======================================================================
141 void 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
155   for(i = 1; i <= myN; i++)
156   {
157     myMaxV(i) = (myB(i) - myA(i)) / 3.0;
158   }
159
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
166   myDone = Standard_False;
167 }
168
169 //=======================================================================
170 //function : SetTol
171 //purpose  : Set algorithm tolerances.
172 //=======================================================================
173 void 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 //=======================================================================
184 void math_GlobOptMin::GetTol(Standard_Real& theDiscretizationTol,
185                              Standard_Real& theSameTol)
186 {
187   theDiscretizationTol = myTol;
188   theSameTol = mySameTol;
189 }
190
191 //=======================================================================
192 //function : ~math_GlobOptMin
193 //purpose  : 
194 //=======================================================================
195 math_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.
204 void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution)
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
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
232   myE1 = minLength * myTol;
233   myE2 = maxLength * myTol;
234
235   myIsFindSingleSolution = isFindSingleSolution;
236   if (isFindSingleSolution)
237   {
238     // Run local optimization 
239     // if current value better than optimal.
240     myE3 = 0.0;
241   }
242   else
243   {
244     if (myC > 1.0)
245       myE3 = - maxLength * myTol / 4.0;
246     else
247       myE3 = - maxLength * myTol * myC / 4.0;
248   }
249
250   isFirstCellFilterInvoke = Standard_True;
251   computeGlobalExtremum(myN);
252
253   myDone = Standard_True;
254 }
255
256 //=======================================================================
257 //function : computeLocalExtremum
258 //purpose  :
259 //=======================================================================
260 Standard_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* aTmp = 
270       dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc);
271     math_NewtonMinimum newtonMinimum(*aTmp);
272     newtonMinimum.SetBoundary(myGlobA, myGlobB);
273     newtonMinimum.Perform(*aTmp, thePnt);
274
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* aTmp =
287       dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
288     math_BFGS bfgs(aTmp->NbVariables());
289     bfgs.Perform(*aTmp, thePnt);
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
305     math_Powell powell(*myFunc, 1e-10);
306     powell.Perform(*myFunc, thePnt, m);
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 //=======================================================================
323 //function : computeInitialValues
324 //purpose  : 
325 //=======================================================================
326 void math_GlobOptMin::computeInitialValues()
327 {
328   Standard_Integer i;
329   math_Vector aCurrPnt(1, myN);
330   math_Vector aBestPnt(1, myN);
331   math_Vector aParamStep(1, myN);
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
362   Standard_Real aLipConst = 0.0, aPrevValDiag, aPrevValProj;
363   Standard_Integer aPntNb = 13;
364   myFunc->Value(myA, aPrevValDiag);
365   aPrevValProj = aPrevValDiag;
366   Standard_Real aStep = (myB - myA).Norm() / aPntNb;
367   aParamStep = (myB - myA) / aPntNb;
368   for(i = 1; i <= aPntNb; i++)
369   {
370     aCurrPnt = myA + aParamStep * i;
371
372     // Walk over diagonal.
373     myFunc->Value(aCurrPnt, aCurrVal);
374     aLipConst = Max (Abs(aCurrVal - aPrevValDiag), aLipConst);
375     aPrevValDiag = aCurrVal;
376
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;
382   }
383
384   aLipConst *= Sqrt(myN) / aStep;
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 //=======================================================================
397 //function : ComputeGlobalExtremum
398 //purpose  :
399 //=======================================================================
400 void 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.
405   Standard_Real aStepBestValue = RealLast();
406   Standard_Real aRealStep = 0.0;
407   math_Vector aStepBestPoint(1, myN);
408   Standard_Boolean isInside = Standard_False;
409   Standard_Real r;
410   Standard_Boolean isReached = Standard_False;
411
412   for(myX(j) = myA(j) + myE1; 
413      (myX(j) < myB(j) + myE1) && (!isReached);
414       myX(j) += myV(j))
415   {
416     if (myX(j) > myB(j))
417     {
418       myX(j) = myB(j);
419       isReached = Standard_True;
420     }
421
422     if (j == 1)
423     {
424       isInside = Standard_False;
425       myFunc->Value(myX, d);
426       r = (d + myZ * myC * aRealStep - 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       aRealStep = myE2 + Abs(myF - d) / myC;
467       myV(1) = Min(aRealStep, myMaxV(1));
468     }
469     else
470     {
471       myV(j) = RealLast() / 2.0;
472       computeGlobalExtremum(j - 1);
473
474       // Nullify steps on lower dimensions.
475       for(i = 1; i < j; i++)
476         myV(i) = 0.0;
477     }
478     // Compute step in (j + 1) dimension according to scale.
479     if (j < myN)
480     {
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       }
489     }
490   }
491 }
492
493 //=======================================================================
494 //function : IsInside
495 //purpose  :
496 //=======================================================================
497 Standard_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 //=======================================================================
513 Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
514 {
515   Standard_Integer i,j;
516   Standard_Boolean isSame = Standard_True;
517   math_Vector aTol(1, myN);
518   aTol = (myB -  myA) * mySameTol;
519
520   // C1 * n^2 = C2 * 3^dim * n
521   if (mySolCount < myMinCellFilterSol)
522   {
523     for(i = 0; i < mySolCount; i++)
524     {
525       isSame = Standard_True;
526       for(j = 1; j <= myN; j++)
527       {
528         if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
529         {
530           isSame = Standard_False;
531           break;
532         }
533       }
534       if (isSame == Standard_True)
535         return Standard_True;
536     }
537   }
538   else
539   {
540     NCollection_CellFilter_Inspector anInspector(myN, Precision::PConfusion());
541     if (isFirstCellFilterInvoke)
542     {
543       myFilter.Reset(myCellSize);
544
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     }
569   }
570   return Standard_False;
571 }
572
573 //=======================================================================
574 //function : NbExtrema
575 //purpose  :
576 //=======================================================================
577 Standard_Integer math_GlobOptMin::NbExtrema()
578 {
579   return mySolCount;
580 }
581
582 //=======================================================================
583 //function : GetF
584 //purpose  :
585 //=======================================================================
586 Standard_Real math_GlobOptMin::GetF()
587 {
588   return myF;
589 }
590
591 //=======================================================================
592 //function : IsDone
593 //purpose  :
594 //=======================================================================
595 Standard_Boolean math_GlobOptMin::isDone()
596 {
597   return myDone;
598 }
599
600 //=======================================================================
601 //function : Points
602 //purpose  :
603 //=======================================================================
604 void 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 }
611
612 //=======================================================================
613 //function : initCellSize
614 //purpose  :
615 //=======================================================================
616 void 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 }