4bbaf12b |
1 | // Created on: 2014-01-20 |
2 | // Created by: Alexaner Malyshev |
3 | // Copyright (c) 2014-2014 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> |
4bbaf12b |
24 | #include <Standard_Integer.hxx> |
25 | #include <Standard_Real.hxx> |
26 | |
27 | const Handle(Standard_Type)& STANDARD_TYPE(math_GlobOptMin) |
28 | { |
29 | static Handle(Standard_Type) _atype = new Standard_Type ("math_GlobOptMin", sizeof (math_GlobOptMin)); |
30 | return _atype; |
31 | } |
32 | |
33 | //======================================================================= |
34 | //function : math_GlobOptMin |
35 | //purpose : Constructor |
36 | //======================================================================= |
37 | math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc, |
38 | const math_Vector& theA, |
39 | const math_Vector& theB, |
5493d334 |
40 | const Standard_Real theC, |
41 | const Standard_Real theDiscretizationTol, |
42 | const Standard_Real theSameTol) |
4bbaf12b |
43 | : myN(theFunc->NbVariables()), |
44 | myA(1, myN), |
45 | myB(1, myN), |
46 | myGlobA(1, myN), |
47 | myGlobB(1, myN), |
48 | myX(1, myN), |
49 | myTmp(1, myN), |
5493d334 |
50 | myV(1, myN), |
51 | myMaxV(1, myN) |
4bbaf12b |
52 | { |
53 | Standard_Integer i; |
54 | |
55 | myFunc = theFunc; |
56 | myC = theC; |
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 | |
74 | myTol = theDiscretizationTol; |
75 | mySameTol = theSameTol; |
76 | |
4bbaf12b |
77 | myDone = Standard_False; |
78 | } |
79 | |
80 | //======================================================================= |
81 | //function : SetGlobalParams |
82 | //purpose : Set params without memory allocation. |
83 | //======================================================================= |
84 | void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc, |
85 | const math_Vector& theA, |
86 | const math_Vector& theB, |
5493d334 |
87 | const Standard_Real theC, |
88 | const Standard_Real theDiscretizationTol, |
89 | const Standard_Real theSameTol) |
4bbaf12b |
90 | { |
91 | Standard_Integer i; |
92 | |
93 | myFunc = theFunc; |
94 | myC = theC; |
95 | myZ = -1; |
96 | mySolCount = 0; |
97 | |
98 | for(i = 1; i <= myN; i++) |
99 | { |
100 | myGlobA(i) = theA(i); |
101 | myGlobB(i) = theB(i); |
102 | |
103 | myA(i) = theA(i); |
104 | myB(i) = theB(i); |
105 | } |
106 | |
5493d334 |
107 | myTol = theDiscretizationTol; |
108 | mySameTol = theSameTol; |
109 | |
4bbaf12b |
110 | myDone = Standard_False; |
111 | } |
112 | |
113 | //======================================================================= |
114 | //function : SetLocalParams |
115 | //purpose : Set params without memory allocation. |
116 | //======================================================================= |
117 | void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA, |
118 | const math_Vector& theLocalB) |
119 | { |
120 | Standard_Integer i; |
121 | |
122 | myZ = -1; |
123 | mySolCount = 0; |
124 | |
125 | for(i = 1; i <= myN; i++) |
126 | { |
127 | myA(i) = theLocalA(i); |
128 | myB(i) = theLocalB(i); |
129 | } |
130 | |
5493d334 |
131 | for(i = 1; i <= myN; i++) |
132 | { |
133 | myMaxV(i) = (myB(i) - myA(i)) / 3.0; |
134 | } |
135 | |
4bbaf12b |
136 | myDone = Standard_False; |
137 | } |
138 | |
5493d334 |
139 | //======================================================================= |
140 | //function : SetTol |
141 | //purpose : Set algorithm tolerances. |
142 | //======================================================================= |
143 | void math_GlobOptMin::SetTol(const Standard_Real theDiscretizationTol, |
144 | const Standard_Real theSameTol) |
145 | { |
146 | myTol = theDiscretizationTol; |
147 | mySameTol = theSameTol; |
148 | } |
149 | |
150 | //======================================================================= |
151 | //function : GetTol |
152 | //purpose : Get algorithm tolerances. |
153 | //======================================================================= |
154 | void math_GlobOptMin::GetTol(Standard_Real& theDiscretizationTol, |
155 | Standard_Real& theSameTol) |
156 | { |
157 | theDiscretizationTol = myTol; |
158 | theSameTol = mySameTol; |
159 | } |
160 | |
4bbaf12b |
161 | //======================================================================= |
162 | //function : ~math_GlobOptMin |
163 | //purpose : |
164 | //======================================================================= |
165 | math_GlobOptMin::~math_GlobOptMin() |
166 | { |
167 | } |
168 | |
169 | //======================================================================= |
170 | //function : Perform |
171 | //purpose : Compute Global extremum point |
172 | //======================================================================= |
173 | // In this algo indexes started from 1, not from 0. |
174 | void math_GlobOptMin::Perform() |
175 | { |
176 | Standard_Integer i; |
177 | |
797d11c6 |
178 | // Compute initial values for myF, myY, myC. |
179 | computeInitialValues(); |
180 | |
4bbaf12b |
181 | // Compute parameters range |
182 | Standard_Real minLength = RealLast(); |
183 | Standard_Real maxLength = RealFirst(); |
184 | for(i = 1; i <= myN; i++) |
185 | { |
186 | Standard_Real currentLength = myB(i) - myA(i); |
187 | if (currentLength < minLength) |
188 | minLength = currentLength; |
189 | if (currentLength > maxLength) |
190 | maxLength = currentLength; |
191 | } |
192 | |
797d11c6 |
193 | myE1 = minLength * myTol; |
194 | myE2 = maxLength * myTol; |
195 | if (myC > 1.0) |
196 | myE3 = - maxLength * myTol / 4.0; |
197 | else |
198 | myE3 = - maxLength * myTol * myC / 4.0; |
4bbaf12b |
199 | |
200 | computeGlobalExtremum(myN); |
201 | |
202 | myDone = Standard_True; |
4bbaf12b |
203 | } |
204 | |
205 | //======================================================================= |
206 | //function : computeLocalExtremum |
207 | //purpose : |
208 | //======================================================================= |
209 | Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt, |
210 | Standard_Real& theVal, |
211 | math_Vector& theOutPnt) |
212 | { |
213 | Standard_Integer i; |
214 | |
215 | //Newton method |
216 | if (dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc)) |
217 | { |
218 | math_MultipleVarFunctionWithHessian* myTmp = |
219 | dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc); |
220 | |
221 | math_NewtonMinimum newtonMinimum(*myTmp, thePnt); |
222 | if (newtonMinimum.IsDone()) |
223 | { |
224 | newtonMinimum.Location(theOutPnt); |
225 | theVal = newtonMinimum.Minimum(); |
226 | } |
227 | else return Standard_False; |
228 | } else |
229 | |
230 | // BFGS method used. |
231 | if (dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc)) |
232 | { |
233 | math_MultipleVarFunctionWithGradient* myTmp = |
234 | dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc); |
07f1a2e6 |
235 | math_BFGS bfgs(myTmp->NbVariables()); |
236 | bfgs.Perform(*myTmp, thePnt); |
4bbaf12b |
237 | if (bfgs.IsDone()) |
238 | { |
239 | bfgs.Location(theOutPnt); |
240 | theVal = bfgs.Minimum(); |
241 | } |
242 | else return Standard_False; |
243 | } else |
244 | |
245 | // Powell method used. |
246 | if (dynamic_cast<math_MultipleVarFunction*>(myFunc)) |
247 | { |
248 | math_Matrix m(1, myN, 1, myN, 0.0); |
249 | for(i = 1; i <= myN; i++) |
250 | m(1, 1) = 1.0; |
251 | |
252 | math_Powell powell(*myFunc, thePnt, m, 1e-10); |
253 | |
254 | if (powell.IsDone()) |
255 | { |
256 | powell.Location(theOutPnt); |
257 | theVal = powell.Minimum(); |
258 | } |
259 | else return Standard_False; |
260 | } |
261 | |
262 | if (isInside(theOutPnt)) |
263 | return Standard_True; |
264 | else |
265 | return Standard_False; |
266 | } |
267 | |
797d11c6 |
268 | //======================================================================= |
269 | //function : computeInitialValues |
270 | //purpose : |
271 | //======================================================================= |
272 | void math_GlobOptMin::computeInitialValues() |
273 | { |
274 | Standard_Integer i; |
275 | math_Vector aCurrPnt(1, myN); |
276 | math_Vector aBestPnt(1, myN); |
277 | |
278 | Standard_Real aCurrVal = RealLast(); |
279 | Standard_Real aBestVal = RealLast(); |
280 | |
281 | // Check functional value in midpoint, low and upp point border and |
282 | // in each point try to perform local optimization. |
283 | aBestPnt = (myA + myB) * 0.5; |
284 | myFunc->Value(aBestPnt, aBestVal); |
285 | |
286 | for(i = 1; i <= 3; i++) |
287 | { |
288 | aCurrPnt = myA + (myB - myA) * (i - 1) / 2.0; |
289 | |
290 | if(computeLocalExtremum(aCurrPnt, aCurrVal, aCurrPnt)) |
291 | { |
292 | // Local Extremum finds better solution than current point. |
293 | if (aCurrVal < aBestVal) |
294 | { |
295 | aBestVal = aCurrVal; |
296 | aBestPnt = aCurrPnt; |
297 | } |
298 | } |
299 | } |
300 | |
301 | myF = aBestVal; |
302 | myY.Clear(); |
303 | for(i = 1; i <= myN; i++) |
304 | myY.Append(aBestPnt(i)); |
305 | mySolCount++; |
306 | |
307 | // Lipschitz const approximation |
308 | Standard_Real aLipConst = 0.0, aPrevVal; |
309 | Standard_Integer aPntNb = 13; |
310 | myFunc->Value(myA, aPrevVal); |
311 | Standard_Real aStep = (myB - myA).Norm() / aPntNb; |
312 | for(i = 1; i <= aPntNb; i++) |
313 | { |
314 | aCurrPnt = myA + (myB - myA) * i / (aPntNb - 1); |
315 | myFunc->Value(aCurrPnt, aCurrVal); |
316 | |
317 | if(Abs(aCurrVal - aPrevVal) / aStep > aLipConst) |
318 | aLipConst = Abs(aCurrVal - aPrevVal) / aStep; |
319 | |
320 | aPrevVal = aCurrVal; |
321 | } |
322 | aLipConst *= Sqrt(myN); |
323 | |
324 | if (aLipConst < myC * 0.1) |
325 | { |
326 | myC = Max(aLipConst * 0.1, 0.01); |
327 | } |
328 | else if (aLipConst > myC * 10) |
329 | { |
330 | myC = Min(myC * 2, 30.0); |
331 | } |
332 | } |
333 | |
4bbaf12b |
334 | //======================================================================= |
335 | //function : ComputeGlobalExtremum |
336 | //purpose : |
337 | //======================================================================= |
338 | void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j) |
339 | { |
340 | Standard_Integer i; |
341 | Standard_Real d; // Functional in moved point. |
342 | Standard_Real val = RealLast(); // Local extrema computed in moved point. |
343 | Standard_Real stepBestValue = RealLast(); |
5493d334 |
344 | Standard_Real realStep = RealLast(); |
345 | math_Vector stepBestPoint(1, myN); |
4bbaf12b |
346 | Standard_Boolean isInside = Standard_False; |
347 | Standard_Real r; |
348 | |
5493d334 |
349 | |
4bbaf12b |
350 | for(myX(j) = myA(j) + myE1; myX(j) < myB(j) + myE1; myX(j) += myV(j)) |
351 | { |
352 | if (myX(j) > myB(j)) |
353 | myX(j) = myB(j); |
354 | |
355 | if (j == 1) |
356 | { |
357 | isInside = Standard_False; |
358 | myFunc->Value(myX, d); |
359 | r = (d - myF) * myZ; |
360 | if(r > myE3) |
361 | { |
362 | isInside = computeLocalExtremum(myX, val, myTmp); |
363 | } |
364 | stepBestValue = (isInside && (val < d))? val : d; |
365 | stepBestPoint = (isInside && (val < d))? myTmp : myX; |
366 | |
367 | // Solutions are close to each other. |
5493d334 |
368 | if (Abs(stepBestValue - myF) < mySameTol * 0.01) |
4bbaf12b |
369 | { |
370 | if (!isStored(stepBestPoint)) |
371 | { |
372 | if ((stepBestValue - myF) * myZ > 0.0) |
373 | myF = stepBestValue; |
374 | for(i = 1; i <= myN; i++) |
375 | myY.Append(stepBestPoint(i)); |
376 | mySolCount++; |
377 | } |
378 | } |
379 | |
380 | // New best solution. |
5493d334 |
381 | if ((stepBestValue - myF) * myZ > mySameTol * 0.01) |
4bbaf12b |
382 | { |
383 | mySolCount = 0; |
384 | myF = stepBestValue; |
385 | myY.Clear(); |
386 | for(i = 1; i <= myN; i++) |
387 | myY.Append(stepBestPoint(i)); |
388 | mySolCount++; |
389 | } |
390 | |
5493d334 |
391 | realStep = myE2 + Abs(myF - d) / myC; |
392 | myV(1) = Min(realStep, myMaxV(1)); |
4bbaf12b |
393 | } |
394 | else |
395 | { |
396 | myV(j) = RealLast() / 2.0; |
397 | computeGlobalExtremum(j - 1); |
398 | } |
5493d334 |
399 | if ((j < myN) && (myV(j + 1) > realStep)) |
4bbaf12b |
400 | { |
5493d334 |
401 | if (realStep > myMaxV(j + 1)) // Case of too big step. |
402 | myV(j + 1) = myMaxV(j + 1); |
4bbaf12b |
403 | else |
5493d334 |
404 | myV(j + 1) = realStep; |
4bbaf12b |
405 | } |
406 | } |
407 | } |
408 | |
409 | //======================================================================= |
410 | //function : IsInside |
411 | //purpose : |
412 | //======================================================================= |
413 | Standard_Boolean math_GlobOptMin::isInside(const math_Vector& thePnt) |
414 | { |
415 | Standard_Integer i; |
416 | |
417 | for(i = 1; i <= myN; i++) |
418 | { |
419 | if (thePnt(i) < myGlobA(i) || thePnt(i) > myGlobB(i)) |
420 | return Standard_False; |
421 | } |
422 | |
423 | return Standard_True; |
424 | } |
425 | //======================================================================= |
426 | //function : IsStored |
427 | //purpose : |
428 | //======================================================================= |
429 | Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt) |
430 | { |
431 | Standard_Integer i,j; |
432 | Standard_Boolean isSame = Standard_True; |
433 | |
434 | for(i = 0; i < mySolCount; i++) |
435 | { |
436 | isSame = Standard_True; |
437 | for(j = 1; j <= myN; j++) |
438 | { |
5493d334 |
439 | if ((Abs(thePnt(j) - myY(i * myN + j))) > (myB(j) - myA(j)) * mySameTol) |
4bbaf12b |
440 | { |
441 | isSame = Standard_False; |
442 | break; |
443 | } |
444 | } |
445 | if (isSame == Standard_True) |
446 | return Standard_True; |
447 | |
448 | } |
449 | return Standard_False; |
450 | } |
451 | |
452 | //======================================================================= |
453 | //function : NbExtrema |
454 | //purpose : |
455 | //======================================================================= |
456 | Standard_Integer math_GlobOptMin::NbExtrema() |
457 | { |
458 | return mySolCount; |
459 | } |
460 | |
461 | //======================================================================= |
462 | //function : GetF |
463 | //purpose : |
464 | //======================================================================= |
465 | Standard_Real math_GlobOptMin::GetF() |
466 | { |
467 | return myF; |
468 | } |
469 | |
470 | //======================================================================= |
471 | //function : IsDone |
472 | //purpose : |
473 | //======================================================================= |
474 | Standard_Boolean math_GlobOptMin::isDone() |
475 | { |
476 | return myDone; |
477 | } |
478 | |
479 | //======================================================================= |
480 | //function : Points |
481 | //purpose : |
482 | //======================================================================= |
483 | void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol) |
484 | { |
485 | Standard_Integer j; |
486 | |
487 | for(j = 1; j <= myN; j++) |
488 | theSol(j) = myY((theIndex - 1) * myN + j); |
489 | } |