Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This include file contains procedure implementations of [pm_optimization](@ref pm_optimization).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Thursday 1:45 AM, August 22, 2019, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : #define SWAP(X,Y,TEMP) TEMP = X; X = Y; Y = TEMP
28 :
29 : !%%%%%%%%%%%%%%%%%%%
30 : #if isBracketMax_ENABLED
31 : !%%%%%%%%%%%%%%%%%%%
32 :
33 8 : bracketed = xlow <= xmax .and. xmax <= xupp .and. flow <= fmax .and. fmax >= fupp
34 :
35 : !%%%%%%%%%%%%%%%%%%%
36 : #elif isBracketMin_ENABLED
37 : !%%%%%%%%%%%%%%%%%%%
38 :
39 8 : bracketed = xlow <= xmin .and. xmin <= xupp .and. flow >= fmin .and. fmin <= fupp
40 :
41 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42 : #elif setBracketMax_ENABLED || setBracketMin_ENABLED
43 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44 :
45 : #if setBracketMax_ENABLED
46 : #define COMPARES_WITH >
47 : #define xmin xmax
48 : #define fmin fmax
49 : #elif setBracketMin_ENABLED
50 : #define COMPARES_WITH <
51 : #else
52 : #error "Unrecognized interface."
53 : #endif
54 : real(RKC) :: lowf, uppf, fnew, xnew, ulim, q, r
55 : real(RKC) , parameter :: GLIMIT = 100._RKC, GOLDEN = real(GOLDEN_RATIO, RKC)
56 : real(RKC) , parameter :: SMALL = sqrt(epsilon(1._RKC))**3 ! Protects against trying to achieve fractional accuracy for a minimum whose value is exactly zero.
57 132 : CHECK_ASSERTION(__LINE__, xlow <= xupp, SK_"@setBracketMin(): The condition `xlow < xupp` must hold. xlow, xupp = "//getStr([xlow, xupp]))
58 : ! Swap `xlow` and `xupp` if needed, to ensure downhill direction from `xlow` to xupp`.
59 44 : lowf = getFunc(xlow)
60 44 : uppf = getFunc(xupp)
61 44 : if (lowf COMPARES_WITH uppf) then
62 16 : ulim = xlow
63 16 : xlow = xupp
64 16 : xupp = ulim
65 : ulim = lowf
66 : lowf = uppf
67 : uppf = ulim
68 : end if
69 44 : xmin = xupp + GOLDEN * (xupp - xlow)
70 44 : fmin = getFunc(xmin)
71 51 : loopRefine: do niter = 1, niter
72 51 : if (uppf COMPARES_WITH fmin) exit loopRefine
73 19 : r = (xupp - xlow) * (uppf - fmin)
74 19 : q = (xupp - xmin) * (uppf - lowf)
75 19 : xnew = xupp - ((xupp - xmin) * q - (xupp - xlow) * r) / (2 * sign(max(abs(q - r), SMALL), q - r))
76 19 : ulim = xupp + GLIMIT * (xmin - xupp)
77 19 : if (0._RKC < (xupp - xnew) * (xnew - xmin)) then
78 12 : fnew = getFunc(xnew)
79 12 : if (fnew COMPARES_WITH fmin) then
80 12 : xlow = xupp
81 : lowf = uppf
82 12 : xupp = xnew
83 : uppf = fnew
84 12 : exit loopRefine
85 0 : else if (uppf COMPARES_WITH fnew) then
86 0 : xmin = xnew
87 0 : fmin = fnew
88 0 : exit loopRefine
89 : end if
90 0 : xnew = xmin + GOLDEN * (xmin - xupp)
91 0 : fnew = getFunc(xnew)
92 7 : else if (0._RKC < (xmin - xnew) * (xnew - ulim)) then
93 7 : fnew = getFunc(xnew)
94 7 : if (fnew COMPARES_WITH fmin) then
95 7 : xupp = xmin
96 7 : xmin = xnew
97 7 : xnew = xmin + GOLDEN * (xmin - xupp)
98 : ! shift.
99 : uppf = fmin
100 7 : fmin = fnew
101 7 : fnew = getFunc(xnew)
102 : end if
103 0 : else if ((xnew - ulim) * (ulim - xmin) < 0._RKC) then
104 0 : xnew = xmin + GOLDEN * (xmin - xupp)
105 0 : fnew = getFunc(xnew)
106 : else
107 0 : xnew = ulim
108 0 : fnew = getFunc(xnew)
109 : end if
110 : ! shift.
111 7 : xlow = xupp
112 7 : xupp = xmin
113 7 : xmin = xnew
114 : ! shift.
115 : lowf = uppf
116 7 : uppf = fmin
117 39 : fmin = fnew
118 : end do loopRefine
119 : ! final swap.
120 44 : SWAP(xmin,xupp,ulim)
121 44 : SWAP(fmin,uppf,ulim)
122 44 : if (xupp < xlow) then
123 16 : SWAP(xlow,xupp,ulim)
124 : end if
125 44 : if (present(flow)) flow = lowf
126 44 : if (present(fupp)) fupp = uppf
127 : #undef xmin
128 : #undef fmin
129 :
130 : !%%%%%%%%%%%%%%%%%%
131 : #elif getMinBrent_ENABLED
132 : !%%%%%%%%%%%%%%%%%%
133 :
134 : integer(IK) :: retin, niter_init
135 : real(RKC) :: lowx, uppx, minf, lot
136 : integer(IK), parameter :: MAXITER = int(100 * precision(xmin) / 53._RKC, IK)
137 4 : if (present(xlow)) then
138 1 : lowx = xlow
139 : else
140 3 : lowx = 0.1_RKC
141 : end if
142 4 : if (present(xupp)) then
143 1 : uppx = xupp
144 : else
145 3 : uppx = 0.9_RKC
146 : end if
147 4 : if (present(tol)) then
148 1 : lot = tol
149 : else
150 : lot = sqrt(epsilon(xmin))
151 : end if
152 4 : if (present(niter)) then
153 3 : niter_init = niter
154 : else
155 : niter_init = MAXITER
156 : end if
157 4 : if (present(xlow) .and. present(xupp)) then
158 1 : xmin = .5_RKC * (xupp - xlow)
159 1 : minf = getFunc(xmin)
160 : else
161 3 : retin = niter_init
162 3 : call setBracketMin(getFunc, retin, xmin, lowx, uppx, minf)
163 3 : if (niter_init < retin) then
164 0 : if (present(niter)) then
165 0 : niter = retin
166 0 : return
167 : else
168 0 : error stop "@getMinBrent(): Bracketing failed. The specified bracket `[xlow, xupp]` may be too wide given the specified `niter` or the input function is concave without minimum."
169 : end if
170 : end if
171 : end if
172 4 : retin = niter_init
173 4 : call setMinBrent(getFunc, xmin, lowx, uppx, minf, lot, retin)
174 4 : if (niter_init < retin) then
175 0 : if (present(niter)) then
176 0 : niter = retin
177 0 : return
178 : else
179 0 : error stop "@getMinBrent(): The Brent minimizer failed to converge. Specifying a larger value for `niter` might resolve the error."
180 : end if
181 : end if
182 4 : if (present(fmin)) fmin = minf
183 :
184 : !%%%%%%%%%%%%%%%%%%
185 : #elif setMinBrent_ENABLED
186 : !%%%%%%%%%%%%%%%%%%
187 :
188 : real(RKC) , parameter :: SMALL = sqrt(epsilon(1._RKC))**3 ! Protects against trying to achieve fractional accuracy for a minimum whose value is exactly zero.
189 : real(RKC) , parameter :: GOLEN_MEAN = .5_RKC * (3._RKC - sqrt(5._RKC)) ! Golden Section switch criterion
190 : real(RKC) :: xz, xw, xold, fold, fv, fw, p, q, r, d, etemp, relerr, abstol, tolby2, center
191 :
192 32 : CHECK_ASSERTION(__LINE__, 0 < tol, SK_"@setMinBrent(): The condition `0 < tol` must hold. tol = "//getStr(tol))
193 32 : CHECK_ASSERTION(__LINE__, 0 < niter, SK_"@setMinBrent(): The condition `0 < niter` must hold. niter = "//getStr(niter))
194 128 : CHECK_ASSERTION(__LINE__, xlow <= xmin .and. xmin <= xupp, SK_"@setMinBrent(): The condition `xlow < xmin .and. xmin < xupp` must hold. xlow, xmin, xupp = "//getStr([xlow, xmin, xupp]))
195 225 : CHECK_ASSERTION(__LINE__, getFunc(xlow) >= fmin .or. fmin <= getFunc(xupp), SK_"@setMinBrent(): The condition `getFunc(xlow) >= fmin .or. fmin =< getFunc(xupp)` must hold. xlow, xmin, xupp, flow, fmin, fupp = "//getStr([xlow, xmin, xupp, getFunc(xlow), fmin, getFunc(xupp)]))
196 :
197 32 : xw = xmin
198 : xmin = xmin
199 : xold = xmin
200 : relerr = 0._RKC
201 : !fmin = getFunc(xmin)
202 32 : fv = fmin
203 : fw = fmin
204 109 : do niter = 1, niter
205 85 : center = .5_RKC * (xlow + xupp)
206 85 : abstol = tol * abs(xmin) + SMALL
207 85 : tolby2 = 2._RKC * abstol
208 85 : if (abs(xmin - center) <= (tolby2 - 0.5_RKC * (xupp - xlow))) return
209 77 : if (abstol < abs(relerr)) then
210 45 : r = (xmin - xw) * (fmin - fv)
211 45 : q = (xmin - xold) * (fmin - fw)
212 45 : p = (xmin - xold) * q - (xmin - xw) * r
213 45 : q = 2.0_RKC * (q - r)
214 45 : if (0._RKC < q) p = - p
215 45 : q = abs(q)
216 : etemp = relerr
217 : relerr = d
218 45 : if (abs(0.5_RKC * q * etemp) <= abs(p) .or. p <= q * (xlow - xmin) .or. q * (xupp - xmin) <= p) then
219 18 : if (xmin < center) then
220 8 : relerr = xupp - xmin
221 : else
222 10 : relerr = xlow - xmin
223 : end if
224 18 : d = GOLEN_MEAN * relerr
225 : else
226 27 : d = p / q
227 27 : xz = xmin + d
228 27 : if (xz - xlow < tolby2 .or. xupp - xz < tolby2) d = sign(abstol, center - xmin)
229 : end if
230 : else
231 32 : if (xmin < center) then
232 15 : relerr = xupp - xmin
233 : else
234 17 : relerr = xlow - xmin
235 : end if
236 32 : d = GOLEN_MEAN * relerr
237 : end if
238 77 : if (abs(d) < abstol) then
239 6 : xz = xmin + sign(abstol, d)
240 : else
241 71 : xz = xmin + d
242 : end if
243 77 : fold = getFunc(xz)
244 101 : if (fmin < fold) then
245 63 : if (xz < xmin) then
246 34 : xlow = xz
247 : else
248 29 : xupp = xz
249 : end if
250 63 : if (fold <= fw .or. xw == xmin) then
251 : xold = xw
252 : fv = fw
253 : xw = xz
254 : fw = fold
255 5 : else if (fold <= fv .or. xold == xmin .or. xold == xw) then
256 : xold = xz
257 : fv = fold
258 : end if
259 : else
260 14 : if (xz < xmin) then
261 8 : xupp = xmin
262 : else
263 6 : xlow = xmin
264 : end if
265 : ! shift.
266 : xold = xw
267 : xw = xmin
268 14 : xmin = xz
269 : ! shift.
270 : fv = fw
271 : fw = fmin
272 14 : fmin = fold
273 : end if
274 : end do
275 :
276 : !%%%%%%%%%%%%%%%%%%%%%%%%
277 : #elif isFailedMinPowell_ENABLED
278 : !%%%%%%%%%%%%%%%%%%%%%%%%
279 :
280 : integer(IK) :: retin
281 4 : real(RKC) :: dirset(size(xmin, 1, IK), size(xmin, 1, IK)), lot, minf
282 : integer(IK) , parameter :: MAXITER = int(100 * precision(xmin) / 53._RKC, IK)
283 2 : call setMatInit(dirset, uppLowDia, vupp = 0._RKC, vlow = 0._RKC, vdia = 1._RKC, nrow = size(xmin, 1, IK), ncol = size(xmin, 1, IK), roff = 0_IK, coff = 0_IK)
284 2 : if (present(tol)) then
285 1 : lot = tol
286 : else
287 1 : lot = sqrt(epsilon(xmin))
288 : end if
289 2 : if (present(niter)) then
290 1 : retin = niter
291 : else
292 1 : retin = MAXITER
293 : end if
294 2 : if (present(fmin)) then
295 1 : minf = fmin
296 : else
297 1 : minf = getFunc(xmin)
298 : end if
299 2 : call setMinPowell(getFunc, xmin, minf, dirset, lot, retin)
300 2 : if (present(fmin)) fmin = minf
301 2 : failed = MAXITER < retin
302 :
303 : !%%%%%%%%%%%%%%%%%%%
304 : #elif setMinPowell_ENABLED
305 : !%%%%%%%%%%%%%%%%%%%
306 :
307 : integer(IK) :: ndim, idim, ibig, retin
308 6 : real(RKC) :: diff, fold, fnew, t, xold(size(xmin, 1, IK)), xnew(size(xmin, 1, IK)), xdir(size(xmin, 1, IK)), xtmp(size(xmin, 1, IK))
309 : real(RKC) , parameter :: SMALL = sqrt(epsilon(1._RKC))**3 ! Protects against trying to achieve fractional accuracy for a minimum whose value is exactly zero.
310 3 : ndim = size(xmin, 1, IK)
311 3 : CHECK_ASSERTION(__LINE__, 0 < tol, SK_"@setMinPowell(): The condition `0 < tol` must hold. tol = "//getStr(tol))
312 3 : CHECK_ASSERTION(__LINE__, 0 < niter, SK_"@setMinPowell(): The condition `0 < niter` must hold. niter = "//getStr(niter))
313 24 : CHECK_ASSERTION(__LINE__, all(ndim == shape(dirset, IK)), SK_"@setMinPowell(): The condition `all(size(xmin) == shape(dirset))` must hold. size(xmin), shape(dirset) = "//getStr([size(xmin, 1, IK), shape(dirset, IK)]))
314 9 : CHECK_ASSERTION(__LINE__, abs(fmin - getFunc(xmin)) < 100 * epsilon(0._RKC), SK_"@setMinPowell(): The condition `fmin - getFunc(xmin)` must hold. fmin, getFunc(xmin) = "//getStr([fmin, getFunc(xmin)]))
315 3 : retin = niter
316 15 : xold = xmin
317 : do
318 6 : fold = fmin
319 : ibig = 0_IK
320 : diff = 0._RKC
321 30 : do idim = 1, ndim
322 120 : xdir = dirset(1 : ndim, idim)
323 24 : fnew = fmin
324 24 : call setMinDir() ! minimize along the specified direction.
325 24 : if (retin < niter) return ! error occurred.
326 30 : if (diff < fnew - fmin) then
327 : diff = fnew - fmin
328 : ibig = idim
329 : end if
330 : end do
331 6 : if (2 * (fold - fmin) <= tol * (abs(fold) + abs(fmin)) + SMALL) return
332 : !if (niter < iter) return
333 15 : xnew = 2 * xmin - xold
334 15 : xdir = xmin - xold
335 15 : xold = xmin
336 3 : fnew = getFunc(xnew)
337 3 : CHECK_ASSERTION(__LINE__, ibig /= 0_IK, SK_"@setMinPowell(): Internal error occurred: ibig = 0")
338 3 : if (fold < fnew) then
339 3 : t = 2 * (fold - 2 * fmin + fnew) * (fold - fmin - diff)**2 - diff * (fold - fnew)**2
340 3 : if (t < 0._RKC) then
341 0 : call setMinDir() ! minimize along the specified direction.
342 0 : if (retin < niter) return ! error occurred.
343 0 : dirset(1 : ndim, ibig) = dirset(1 : ndim, ndim)
344 0 : dirset(1 : ndim, ndim) = xdir
345 : end if
346 : end if
347 : end do
348 :
349 : contains
350 :
351 189 : function getFuncDir(x) result(func)
352 : real(RKC), intent(in) :: x
353 : real(RKC) :: func
354 945 : xtmp = xmin + x * xdir
355 189 : func = getFunc(xtmp)
356 189 : end function
357 :
358 24 : subroutine setMinDir()
359 : integer(IK) :: jdim
360 : real(RKC) :: minx, xlow, xupp
361 : ! \todo: There must be a better initial bracketing guess at each iteration.
362 24 : xlow = 0._RKC
363 24 : xupp = 1._RKC
364 24 : niter = 1000
365 24 : call setBracketMin(getFuncDir, niter, minx, xlow, xupp, fmin)
366 24 : if (1000 < niter) then
367 0 : niter = retin + 1
368 : else
369 24 : call setMinBrent(getFuncDir, minx, xlow, xupp, fmin, tol, niter)
370 24 : do concurrent(jdim = 1 : ndim)
371 96 : xdir(jdim) = xdir(jdim) * minx
372 120 : xmin(jdim) = xmin(jdim) + xdir(jdim)
373 : end do
374 : end if
375 24 : end subroutine
376 :
377 : #else
378 : !%%%%%%%%%%%%%%%%%%%%%%%%
379 : #error "Unrecognized interface."
380 : !%%%%%%%%%%%%%%%%%%%%%%%%
381 : #endif
382 : #undef COMPARES_WITH
383 : #undef SWAP
|
