Line data Source code
1 : !> \file mo_nelmin.f90
2 : !> \brief \copybrief mo_nelmin
3 : !> \details \copydetails mo_nelmin
4 :
5 : !> \brief Nelder-Mead algorithm
6 : !> \details This module provides NELMIN, which minimizes a function using the Nelder-Mead algorithm
7 : !! with the Applied Statistics algorithms No. 047.
8 : !! Original FORTRAN77 version by R ONeill.
9 : !! FORTRAN90 version by John Burkardt.
10 : !> \author Matthias Cuntz
11 : !> \date 2012
12 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
13 : !! FORCES is released under the LGPLv3+ license \license_note
14 : MODULE mo_nelmin
15 :
16 : USE mo_kind, ONLY: i4, sp, dp
17 : USE mo_utils, ONLY: ne
18 :
19 : IMPLICIT NONE
20 :
21 : PUBLIC :: nelmin ! minimizes a function using the Nelder-Mead algorithm
22 : PUBLIC :: nelminxy ! same as nelmin but pass x and y to function
23 : PUBLIC :: nelminrange ! same as nelmin but with given range of parameter
24 :
25 : ! ------------------------------------------------------------------
26 : !> \brief Minimizes a user-specified function using the Nelder-Mead algorithm.
27 : !> \details Simplex function minimisation procedure due to Nelder and Mead (1965),
28 : !! as implemented by O'Neill(1971, Appl.Statist. 20, 338-45), with
29 : !! subsequent comments by Chambers+Ertel(1974, 23, 250-1), Benyon(1976,
30 : !! 25, 97) and Hill(1978, 27, 380-2)
31 : !!
32 : !! The function to be minimized is the first argument of nelmin and must be defined as
33 : !! \code{.f90}
34 : !! FUNCTION func(p)
35 : !! USE mo_kind, ONLY: dp
36 : !! IMPLICIT NONE
37 : !! REAL(dp), DIMENSION(:), INTENT(IN) :: p
38 : !! REAL(dp) :: func
39 : !! END FUNCTION func
40 : !! \endcode
41 : !!
42 : !! and for nelminxy
43 : !! \code{.f90}
44 : !! FUNCTION func(p,x,y)
45 : !! USE mo_kind, ONLY: dp
46 : !! IMPLICIT NONE
47 : !! REAL(dp), DIMENSION(:), INTENT(IN) :: p
48 : !! REAL(dp), DIMENSION(:), INTENT(IN) :: x
49 : !! REAL(dp), DIMENSION(:), INTENT(IN) :: y
50 : !! REAL(dp) :: func
51 : !! END FUNCTION func
52 : !! \endcode
53 : !!
54 : !! This routine does not include a termination test using the
55 : !! fitting of a quadratic surface.
56 : !!
57 : !! \b Calling sequence
58 : !! \code{.f90}
59 : !! pmin = nelmin(func, pstart, funcmin, varmin, step, konvge, maxeval, &
60 : !! neval, numrestart, ierror)
61 : !! pmin = nelminxy(func, pstart, xx, yy, funcmin, varmin, step, konvge, maxeval, &
62 : !! neval, numrestart, ierror)
63 : !! pmin = nelmin(func, pstart, prange, funcmin, varmin, step, konvge, maxeval, &
64 : !! neval, numrestart, ierror)
65 : !! \endcode
66 : !!
67 : !! Original FORTRAN77 version by R ONeill.
68 : !! FORTRAN90 version by John Burkardt.
69 : !!
70 : !! \b Example
71 : !! \code{.f90}
72 : !! pstart(1:n) = (/ -1.2_dp, 1.0_dp /)
73 : !! step(1:n) = (/ 1.0_dp, 1.0_dp /)
74 : !! pmin = nelmin(rosenbrock, pstart, step=step, ierror=ierror)
75 : !! \endcode
76 : !!
77 : !! -> see also example in test directory
78 : !!
79 : !! \b Literature
80 : !! 1. Simplex function minimisation procedure due to Nelder and Mead (1965),
81 : !! as implemented by O'Neill(1971, Appl.Statist. 20, 338-45), with
82 : !! subsequent comments by Chambers+Ertel(1974, 23, 250-1), Benyon(1976,
83 : !! 25, 97) and Hill(1978, 27, 380-2)
84 : !! 2. John Nelder, Roger Mead, A simplex method for function minimization,
85 : !! Computer Journal, Volume 7, 1965, pages 308-313.
86 : !! 3. R ONeill, Algorithm AS 47: Function Minimization Using a Simplex Procedure,
87 : !! Applied Statistics, Volume 20, Number 3, 1971, pages 338-345.
88 : !!
89 : !> \param[in] "real(dp) :: func(p,xx,yy)" Function on which to search the minimum
90 : !> \param[in] "real(dp) :: pstart(:)" Starting point for the iteration.
91 : !> \param[in] "real(dp) :: xx" (ONLY nelminxy) First values to pass as function arguments
92 : !> \param[in] "real(dp) :: yy" (ONLY nelminxy) Second values to pass as function arguments
93 : !> \param[in] "real(dp) :: prange(:,2)" (ONLY nelminrange) Range of parameters (upper and lower bound).
94 : !> \param[in] "real(dp), optional :: varmin" the terminating limit for the variance
95 : !! of the function values. varmin>0 is required.
96 : !> \param[in] "real(dp), optional :: step(:)" determines the size and shape of the initial simplex.
97 : !! The relative magnitudes of its elements should reflect
98 : !! the units of the variables. size(step)=size(start)
99 : !> \param[in] "integer(i4), optional :: konvge" the convergence check is carried out every konvge iterations.
100 : !! konvge>0 is required.
101 : !> \param[in] "integer(i4), optional :: maxeval" the maximum number of function evaluations.
102 : !! default: 1000
103 : !> \param[out] "real(dp), optional :: funcmin" the minimum value of the function.
104 : !> \param[out] "integer(i4), optional :: neval" the number of function evaluations used.
105 : !> \param[out] "integer(i4), optional :: numrestart" the number of restarts.
106 : !> \param[out] "integer(i4), optional :: ierror" error indicator.
107 : !! 0: no errors detected.
108 : !! 1: varmin or konvge have an illegal value.
109 : !! 2: terminated because maxeval was exceeded without convergence.
110 : !> \param[out] "real(dp), optional, allocatable :: history(:)" the history of best function values, history(neval)=funcmin
111 : !> \return real(sp/dp) :: nelmin(size(start)) — coordinates of the point which is estimated to minimize the function.
112 :
113 : !> \author Matthias Cuntz
114 : !> \date Jul 2012
115 : !! - i4, dp, intent
116 : !! - function, optional
117 : !! - nelminxy
118 : !> \author Juliane Mai
119 : !> \date Aug 2012
120 : !! - nelminrange
121 : !> \date Dec 2012
122 : !! - history output
123 : INTERFACE nelminrange
124 : MODULE PROCEDURE nelminrange_dp, nelminrange_sp
125 : END INTERFACE nelminrange
126 :
127 : ! ------------------------------------------------------------------
128 :
129 : PRIVATE
130 :
131 : ! ------------------------------------------------------------------
132 :
133 : CONTAINS
134 :
135 : !> \copydoc nelminrange
136 19 : function nelmin(func, pstart, varmin, step, konvge, maxeval, &
137 : funcmin, neval, numrestart, ierror, history)
138 :
139 : implicit none
140 :
141 : INTERFACE
142 : FUNCTION func(pp)
143 : USE mo_kind, ONLY: dp
144 : IMPLICIT NONE
145 : REAL(dp), DIMENSION(:), INTENT(IN) :: pp
146 : REAL(dp) :: func
147 : END FUNCTION func
148 : END INTERFACE
149 : real(dp), intent(IN) :: pstart(:)
150 : real(dp), optional, intent(IN) :: varmin
151 : real(dp), optional, intent(IN) :: step(:)
152 : integer(i4), optional, intent(IN) :: konvge
153 : integer(i4), optional, intent(IN) :: maxeval
154 : real(dp), optional, intent(OUT) :: funcmin
155 : integer(i4), optional, intent(OUT) :: neval
156 : integer(i4), optional, intent(OUT) :: numrestart
157 : integer(i4), optional, intent(OUT) :: ierror
158 : real(dp), optional, intent(OUT), allocatable :: history(:) ! History of objective function values
159 : real(dp) :: nelmin(size(pstart))
160 :
161 : real(dp), parameter :: ccoeff = 0.5_dp
162 : real(dp), parameter :: ecoeff = 2.0_dp
163 : real(dp), parameter :: rcoeff = 1.0_dp
164 : real(dp), parameter :: eps = 0.001_dp
165 18 : real(dp) :: ipstart(size(pstart))
166 18 : real(dp) :: ifuncmin
167 18 : real(dp) :: ivarmin
168 24 : real(dp) :: istep(size(pstart))
169 : integer(i4) :: ikonvge
170 : integer(i4) :: imaxeval
171 : integer(i4) :: ineval
172 : integer(i4) :: inumrestart
173 : integer(i4) :: iierror
174 : integer(i4) :: n, nn
175 : integer(i4) :: i
176 : integer(i4) :: ihi
177 : integer(i4) :: ilo
178 : integer(i4) :: j
179 : integer(i4) :: jcount
180 : integer(i4) :: l
181 18 : real(dp) :: del
182 66 : real(dp) :: p(size(pstart),size(pstart)+1)
183 18 : real(dp) :: p2star(size(pstart))
184 24 : real(dp) :: pbar(size(pstart))
185 24 : real(dp) :: pstar(size(pstart))
186 6 : real(dp) :: rq
187 6 : real(dp) :: x
188 30 : real(dp) :: y(size(pstart)+1)
189 6 : real(dp) :: y2star
190 6 : real(dp) :: ylo
191 6 : real(dp) :: ystar
192 6 : real(dp) :: z
193 18 : real(dp) :: dn, dnn
194 36 : real(dp) :: p0(size(pstart)), y0
195 6 : real(dp), allocatable :: history_tmp(:)
196 : !
197 : ! Defaults
198 : !
199 18 : nelmin(:) = 0.
200 6 : if (present(varmin)) then
201 3 : if (varmin <= 0.0_dp) stop 'Error nelmin: varmin<0'
202 : ivarmin = varmin
203 : else
204 : ivarmin = 1.0e-9_dp
205 : endif
206 : ! maximal number of function evaluations
207 6 : if (present(maxeval)) then
208 4 : if (maxeval <= 1) stop 'Error nelmin: maxeval<=1'
209 : imaxeval = maxeval
210 : else
211 : imaxeval = 1000
212 : endif
213 : ! history output
214 6 : if (present(history)) then
215 : ! worst case length
216 3 : allocate(history_tmp(imaxeval+3*size(ipstart)+1))
217 : end if
218 6 : if (present(konvge)) then
219 4 : ikonvge = konvge
220 : else
221 2 : ikonvge = imaxeval / 10
222 : endif
223 6 : if (ikonvge < 1) stop 'Error nelmin: konvg<1'
224 : !
225 6 : if (present(step)) then
226 6 : istep = step
227 : else ! if not given, deviate initial by 1%
228 4 : y0 = func(pstart)
229 12 : do i=1, size(pstart)
230 24 : p0 = pstart
231 8 : if (ne(p0(i),0.0_dp)) then
232 8 : p0(i) = 1.01_dp*p0(i)
233 : else
234 0 : p0(i) = 0.01_dp
235 : endif
236 12 : istep(i) = abs(func(p0) - y0)
237 : enddo
238 : endif
239 : !
240 : ! Check the input parameters.
241 : !
242 : !
243 : ! Initialization.
244 : !
245 18 : ipstart = pstart
246 6 : n = size(ipstart)
247 6 : ineval = 0_i4
248 6 : inumrestart = 0_i4
249 6 : jcount = ikonvge
250 6 : dn = real(n, dp)
251 6 : nn = n + 1
252 6 : dnn = real(nn, dp)
253 6 : del = 1.0_dp
254 6 : rq = ivarmin * dn
255 : !
256 : ! Initial or restarted loop.
257 : !
258 0 : do
259 18 : p(1:n,nn) = ipstart(1:n)
260 6 : y(nn) = func(ipstart)
261 6 : ineval = ineval + 1
262 6 : if (present(history)) history_tmp(ineval) = y(nn)
263 : !
264 : ! Define the initial simplex.
265 : !
266 18 : do j = 1, n
267 12 : x = ipstart(j)
268 12 : ipstart(j) = ipstart(j) + istep(j) * del
269 36 : p(1:n,j) = ipstart(1:n)
270 12 : y(j) = func(ipstart)
271 12 : ineval = ineval + 1
272 12 : ipstart(j) = x
273 18 : if (present(history)) history_tmp(ineval) = Min( y(j), history_tmp(ineval-1) )
274 : end do
275 : !
276 : ! Find highest and lowest Y values. FUNCMIN = Y(IHI) indicates
277 : ! the vertex of the simplex to be replaced.
278 : !
279 24 : ilo = minloc(y(1:n+1), 1)
280 6 : ylo = y(ilo)
281 : !
282 : ! Inner loop.
283 : !
284 591 : do while ( ineval < imaxeval )
285 : !
286 : ! FUNCMIN is, of course, the HIGHEST value???
287 : !
288 2364 : ihi = maxloc(y(1:n+1), 1)
289 591 : ifuncmin = y(ihi)
290 : !
291 : ! Calculate PBAR, the centroid of the simplex vertices
292 : ! excepting the vertex with Y value FUNCMIN.
293 : !
294 1773 : do i = 1, n
295 5319 : pbar(i) = ( sum(p(i,1:n+1)) - p(i,ihi) ) / dn
296 : end do
297 : !
298 : ! Reflection through the centroid.
299 : !
300 1773 : pstar(1:n) = pbar(1:n) + rcoeff * ( pbar(1:n) - p(1:n,ihi) )
301 591 : ystar = func(pstar)
302 591 : ineval = ineval + 1
303 591 : if (present(history)) history_tmp(ineval) = Min( ystar, history_tmp(ineval-1) )
304 : !
305 : ! Successful reflection, so extension.
306 : !
307 591 : if ( ystar < ylo ) then
308 :
309 651 : p2star(1:n) = pbar(1:n) + ecoeff * ( pstar(1:n) - pbar(1:n) )
310 217 : y2star = func(p2star)
311 217 : ineval = ineval + 1
312 217 : if (present(history)) history_tmp(ineval) = Min( y2star, history_tmp(ineval-1) )
313 : !
314 : ! Retain extension or contraction.
315 : !
316 217 : if ( ystar < y2star ) then
317 333 : p(1:n,ihi) = pstar(1:n)
318 111 : y(ihi) = ystar
319 : else
320 318 : p(1:n,ihi) = p2star(1:n)
321 106 : y(ihi) = y2star
322 : end if
323 : !
324 : ! No extension.
325 : !
326 : else
327 : l = 0
328 1496 : do i = 1, nn
329 1496 : if ( ystar < y(i) ) l = l + 1
330 : end do
331 :
332 374 : if ( 1 < l ) then
333 237 : p(1:n,ihi) = pstar(1:n)
334 79 : y(ihi) = ystar
335 : !
336 : ! Contraction on the Y(IHI) side of the centroid.
337 : !
338 295 : else if ( l == 0 ) then
339 735 : p2star(1:n) = pbar(1:n) + ccoeff * ( p(1:n,ihi) - pbar(1:n) )
340 245 : y2star = func(p2star)
341 245 : ineval = ineval + 1
342 245 : if (present(history)) history_tmp(ineval) = Min( y2star, history_tmp(ineval-1) )
343 : !
344 : ! Contract the whole simplex.
345 : !
346 245 : if ( y(ihi) < y2star ) then
347 4 : do j = 1, n + 1
348 9 : p(1:n,j) = ( p(1:n,j) + p(1:n,ilo) ) * 0.5_dp
349 9 : nelmin(1:n) = p(1:n,j)
350 3 : y(j) = func(nelmin)
351 3 : ineval = ineval + 1
352 4 : if (present(history)) history_tmp(ineval) = Min( y(j), history_tmp(ineval-1) )
353 : end do
354 4 : ilo = minloc(y(1:n+1), 1)
355 1 : ylo = y(ilo)
356 :
357 1 : cycle
358 : !
359 : ! Retain contraction.
360 : !
361 : else
362 732 : p(1:n,ihi) = p2star(1:n)
363 244 : y(ihi) = y2star
364 : end if
365 : !
366 : ! Contraction on the reflection side of the centroid.
367 : !
368 50 : else if ( l == 1 ) then
369 :
370 150 : p2star(1:n) = pbar(1:n) + ccoeff * ( pstar(1:n) - pbar(1:n) )
371 50 : y2star = func(p2star)
372 50 : ineval = ineval + 1
373 50 : if (present(history)) history_tmp(ineval) = Min( y2star, history_tmp(ineval-1) )
374 : !
375 : ! Retain reflection?
376 : !
377 50 : if ( y2star <= ystar ) then
378 150 : p(1:n,ihi) = p2star(1:n)
379 50 : y(ihi) = y2star
380 : else
381 0 : p(1:n,ihi) = pstar(1:n)
382 0 : y(ihi) = ystar
383 : end if
384 :
385 : end if
386 :
387 : end if
388 : !
389 : ! Check if YLO improved.
390 : !
391 590 : if ( y(ihi) < ylo ) then
392 319 : ylo = y(ihi)
393 319 : ilo = ihi
394 : end if
395 :
396 590 : jcount = jcount - 1
397 :
398 590 : if ( 0 < jcount ) cycle
399 : !
400 : ! Check to see if minimum reached.
401 : !
402 38 : if ( ineval <= imaxeval ) then
403 128 : jcount = ikonvge
404 128 : x = sum ( y(1:n+1) ) / dnn
405 128 : z = sum ( ( y(1:n+1) - x )**2 )
406 32 : if ( z <= rq ) exit
407 : end if
408 :
409 : end do
410 : !
411 : ! Factorial tests to check that FUNCMIN is a local minimum.
412 : !
413 18 : nelmin(1:n) = p(1:n,ilo)
414 6 : ifuncmin = y(ilo)
415 :
416 6 : if ( imaxeval < ineval ) then
417 : iierror = 2
418 : exit
419 : end if
420 :
421 18 : iierror = 0
422 :
423 18 : do i = 1, n
424 12 : del = istep(i) * eps
425 12 : nelmin(i) = nelmin(i) + del
426 12 : z = func(nelmin)
427 12 : ineval = ineval + 1
428 12 : if (present(history)) history_tmp(ineval) = Min( z, history_tmp(ineval-1) )
429 12 : if ( z < ifuncmin ) then
430 : iierror = 2
431 : exit
432 : end if
433 12 : nelmin(i) = nelmin(i) - del - del
434 12 : z = func(nelmin)
435 12 : ineval = ineval + 1
436 12 : if (present(history)) history_tmp(ineval) = Min( z, history_tmp(ineval-1) )
437 12 : if ( z < ifuncmin ) then
438 : iierror = 2
439 : exit
440 : end if
441 18 : nelmin(i) = nelmin(i) + del
442 : end do
443 :
444 6 : if ( iierror == 0 ) exit
445 : !
446 : ! Restart the procedure.
447 : !
448 0 : ipstart(1:n) = nelmin(1:n)
449 0 : del = eps
450 0 : inumrestart = inumrestart + 1
451 :
452 : end do
453 :
454 6 : if (present(funcmin)) then
455 1 : funcmin = ifuncmin
456 : endif
457 6 : if (present(neval)) then
458 1 : neval = ineval
459 : endif
460 6 : if (present(numrestart)) then
461 1 : numrestart = inumrestart
462 : endif
463 6 : if (present(ierror)) then
464 1 : ierror = iierror
465 : endif
466 6 : if (present(history)) then
467 3 : allocate(history(ineval))
468 158 : history(:) = history_tmp(1:ineval)
469 : end if
470 :
471 6 : end function nelmin
472 :
473 :
474 : !> \copydoc nelminrange
475 10 : function nelminxy(func, pstart, xx, yy, varmin, step, konvge, maxeval, &
476 : funcmin, neval, numrestart, ierror, history)
477 :
478 : implicit none
479 :
480 : INTERFACE
481 : FUNCTION func(pp, xxx, yyy)
482 : USE mo_kind, ONLY: dp
483 : IMPLICIT NONE
484 : REAL(dp), DIMENSION(:), INTENT(IN) :: pp
485 : REAL(dp), DIMENSION(:), INTENT(IN) :: xxx
486 : REAL(dp), DIMENSION(:), INTENT(IN) :: yyy
487 : REAL(dp) :: func
488 : END FUNCTION func
489 : END INTERFACE
490 : real(dp), intent(IN) :: pstart(:)
491 : real(dp), intent(IN) :: xx(:)
492 : real(dp), intent(IN) :: yy(:)
493 : real(dp), optional, intent(OUT) :: funcmin
494 : real(dp), optional, intent(IN) :: varmin
495 : real(dp), optional, intent(IN) :: step(:)
496 : integer(i4), optional, intent(IN) :: konvge
497 : integer(i4), optional, intent(IN) :: maxeval
498 : integer(i4), optional, intent(OUT) :: neval
499 : integer(i4), optional, intent(OUT) :: numrestart
500 : integer(i4), optional, intent(OUT) :: ierror
501 : real(dp), optional, intent(OUT), allocatable :: history(:) ! History of objective function values
502 : real(dp) :: nelminxy(size(pstart))
503 :
504 : real(dp), parameter :: ccoeff = 0.5_dp
505 : real(dp), parameter :: ecoeff = 2.0_dp
506 : real(dp), parameter :: rcoeff = 1.0_dp
507 : real(dp), parameter :: eps = 0.001_dp
508 3 : real(dp) :: ipstart(size(pstart))
509 3 : real(dp) :: ifuncmin
510 3 : real(dp) :: ivarmin
511 4 : real(dp) :: istep(size(pstart))
512 : integer(i4) :: ikonvge
513 : integer(i4) :: imaxeval
514 : integer(i4) :: ineval
515 : integer(i4) :: inumrestart
516 : integer(i4) :: iierror
517 : integer(i4) :: n, nn
518 : integer(i4) :: i
519 : integer(i4) :: ihi
520 : integer(i4) :: ilo
521 : integer(i4) :: j
522 : integer(i4) :: jcount
523 : integer(i4) :: l
524 3 : real(dp) :: del
525 11 : real(dp) :: p(size(pstart),size(pstart)+1)
526 3 : real(dp) :: p2star(size(pstart))
527 4 : real(dp) :: pbar(size(pstart))
528 4 : real(dp) :: pstar(size(pstart))
529 1 : real(dp) :: rq
530 1 : real(dp) :: x
531 5 : real(dp) :: y(size(pstart)+1)
532 1 : real(dp) :: y2star
533 1 : real(dp) :: ylo
534 1 : real(dp) :: ystar
535 1 : real(dp) :: z
536 3 : real(dp) :: dn, dnn
537 6 : real(dp) :: p0(size(pstart)), y0
538 1 : real(dp), allocatable :: history_tmp(:)
539 : !
540 : ! Defaults
541 : !
542 3 : nelminxy(:) = 0.
543 1 : if (present(varmin)) then
544 0 : if (varmin <= 0.0_dp) stop 'Error nelminxy: varmin<0'
545 : ivarmin = varmin
546 : else
547 : ivarmin = 1.0e-9_dp
548 : endif
549 : ! maximal number of function evaluations
550 1 : if (present(maxeval)) then
551 0 : if (maxeval <= 1) stop 'Error nelminxy: maxeval<=1'
552 : imaxeval = maxeval
553 : else
554 : imaxeval = 1000
555 : endif
556 : ! history output
557 1 : if (present(history)) then
558 : ! worst case length
559 3 : allocate(history_tmp(imaxeval+3*size(ipstart)+1))
560 : end if
561 1 : if (present(konvge)) then
562 0 : ikonvge = konvge
563 : else
564 1 : ikonvge = imaxeval / 10
565 : endif
566 1 : if (ikonvge < 1) stop 'Error nelminxy: konvg<1'
567 1 : if (size(xx) /= size(yy)) stop 'Error nelminxy: size(xx) /= size(yy)'
568 : !
569 1 : if (present(step)) then
570 0 : istep = step
571 : else ! if not given, deviate initial by 1%
572 1 : y0 = func(pstart, xx, yy)
573 3 : do i=1, size(pstart)
574 6 : p0 = pstart
575 2 : if (ne(p0(i),0.0_dp)) then
576 2 : p0(i) = 1.01_dp*p0(i)
577 : else
578 0 : p0(i) = 0.01_dp
579 : endif
580 3 : istep(i) = abs(func(p0, xx, yy) - y0)
581 : enddo
582 : endif
583 : !
584 : ! Check the input parameters.
585 : !
586 : !
587 : ! Initialization.
588 : !
589 3 : ipstart = pstart
590 1 : n = size(ipstart)
591 1 : ineval = 0_i4
592 1 : inumrestart = 0_i4
593 1 : jcount = ikonvge
594 1 : dn = real(n, dp)
595 1 : nn = n + 1
596 1 : dnn = real(nn, dp)
597 1 : del = 1.0_dp
598 1 : rq = ivarmin * dn
599 : !
600 : ! Initial or restarted loop.
601 : !
602 0 : do
603 3 : p(1:n,nn) = ipstart(1:n)
604 1 : y(nn) = func(ipstart, xx, yy)
605 1 : ineval = ineval + 1
606 1 : if (present(history)) history_tmp(ineval) = y(nn)
607 : !
608 : ! Define the initial simplex.
609 : !
610 3 : do j = 1, n
611 2 : x = ipstart(j)
612 2 : ipstart(j) = ipstart(j) + istep(j) * del
613 6 : p(1:n,j) = ipstart(1:n)
614 2 : y(j) = func(ipstart, xx, yy)
615 2 : ineval = ineval + 1
616 2 : ipstart(j) = x
617 3 : if (present(history)) history_tmp(ineval) = Min( y(j), history_tmp(ineval-1) )
618 : end do
619 : !
620 : ! Find highest and lowest Y values. FUNCMIN = Y(IHI) indicates
621 : ! the vertex of the simplex to be replaced.
622 : !
623 4 : ilo = minloc(y(1:n+1), 1)
624 1 : ylo = y(ilo)
625 : !
626 : ! Inner loop.
627 : !
628 201 : do while ( ineval < imaxeval )
629 : !
630 : ! FUNCMIN is, of course, the HIGHEST value???
631 : !
632 804 : ihi = maxloc(y(1:n+1), 1)
633 201 : ifuncmin = y(ihi)
634 : !
635 : ! Calculate PBAR, the centroid of the simplex vertices
636 : ! excepting the vertex with Y value FUNCMIN.
637 : !
638 603 : do i = 1, n
639 1809 : pbar(i) = ( sum(p(i,1:n+1)) - p(i,ihi) ) / dn
640 : end do
641 : !
642 : ! Reflection through the centroid.
643 : !
644 603 : pstar(1:n) = pbar(1:n) + rcoeff * ( pbar(1:n) - p(1:n,ihi) )
645 201 : ystar = func(pstar, xx, yy)
646 201 : ineval = ineval + 1
647 201 : if (present(history)) history_tmp(ineval) = Min( ystar, history_tmp(ineval-1) )
648 : !
649 : ! Successful reflection, so extension.
650 : !
651 201 : if ( ystar < ylo ) then
652 :
653 240 : p2star(1:n) = pbar(1:n) + ecoeff * ( pstar(1:n) - pbar(1:n) )
654 80 : y2star = func(p2star, xx, yy)
655 80 : ineval = ineval + 1
656 80 : if (present(history)) history_tmp(ineval) = Min( y2star, history_tmp(ineval-1) )
657 : !
658 : ! Retain extension or contraction.
659 : !
660 80 : if ( ystar < y2star ) then
661 111 : p(1:n,ihi) = pstar(1:n)
662 37 : y(ihi) = ystar
663 : else
664 129 : p(1:n,ihi) = p2star(1:n)
665 43 : y(ihi) = y2star
666 : end if
667 : !
668 : ! No extension.
669 : !
670 : else
671 : l = 0
672 484 : do i = 1, nn
673 484 : if ( ystar < y(i) ) l = l + 1
674 : end do
675 :
676 121 : if ( 1 < l ) then
677 87 : p(1:n,ihi) = pstar(1:n)
678 29 : y(ihi) = ystar
679 : !
680 : ! Contraction on the Y(IHI) side of the centroid.
681 : !
682 92 : else if ( l == 0 ) then
683 237 : p2star(1:n) = pbar(1:n) + ccoeff * ( p(1:n,ihi) - pbar(1:n) )
684 79 : y2star = func(p2star, xx, yy)
685 79 : ineval = ineval + 1
686 79 : if (present(history)) history_tmp(ineval) = Min( y2star, history_tmp(ineval-1) )
687 : !
688 : ! Contract the whole simplex.
689 : !
690 79 : if ( y(ihi) < y2star ) then
691 4 : do j = 1, n + 1
692 9 : p(1:n,j) = ( p(1:n,j) + p(1:n,ilo) ) * 0.5_dp
693 9 : nelminxy(1:n) = p(1:n,j)
694 3 : y(j) = func(nelminxy, xx, yy)
695 3 : ineval = ineval + 1
696 4 : if (present(history)) history_tmp(ineval) = Min( y(j), history_tmp(ineval-1) )
697 : end do
698 4 : ilo = minloc(y(1:n+1), 1)
699 1 : ylo = y(ilo)
700 :
701 1 : cycle
702 : !
703 : ! Retain contraction.
704 : !
705 : else
706 234 : p(1:n,ihi) = p2star(1:n)
707 78 : y(ihi) = y2star
708 : end if
709 : !
710 : ! Contraction on the reflection side of the centroid.
711 : !
712 13 : else if ( l == 1 ) then
713 :
714 39 : p2star(1:n) = pbar(1:n) + ccoeff * ( pstar(1:n) - pbar(1:n) )
715 13 : y2star = func(p2star, xx, yy)
716 13 : ineval = ineval + 1
717 13 : if (present(history)) history_tmp(ineval) = Min( y2star, history_tmp(ineval-1) )
718 : !
719 : ! Retain reflection?
720 : !
721 13 : if ( y2star <= ystar ) then
722 39 : p(1:n,ihi) = p2star(1:n)
723 13 : y(ihi) = y2star
724 : else
725 0 : p(1:n,ihi) = pstar(1:n)
726 0 : y(ihi) = ystar
727 : end if
728 :
729 : end if
730 :
731 : end if
732 : !
733 : ! Check if YLO improved.
734 : !
735 200 : if ( y(ihi) < ylo ) then
736 97 : ylo = y(ihi)
737 97 : ilo = ihi
738 : end if
739 :
740 200 : jcount = jcount - 1
741 :
742 200 : if ( 0 < jcount ) cycle
743 : !
744 : ! Check to see if minimum reached.
745 : !
746 3 : if ( ineval <= imaxeval ) then
747 8 : jcount = ikonvge
748 8 : x = sum ( y(1:n+1) ) / dnn
749 8 : z = sum ( ( y(1:n+1) - x )**2 )
750 2 : if ( z <= rq ) exit
751 : end if
752 :
753 : end do
754 : !
755 : ! Factorial tests to check that FUNCMIN is a local minimum.
756 : !
757 3 : nelminxy(1:n) = p(1:n,ilo)
758 1 : ifuncmin = y(ilo)
759 :
760 1 : if ( imaxeval < ineval ) then
761 : iierror = 2
762 : exit
763 : end if
764 :
765 3 : iierror = 0
766 :
767 3 : do i = 1, n
768 2 : del = istep(i) * eps
769 2 : nelminxy(i) = nelminxy(i) + del
770 2 : z = func(nelminxy, xx, yy)
771 2 : ineval = ineval + 1
772 2 : if (present(history)) history_tmp(ineval) = Min( z, history_tmp(ineval-1) )
773 2 : if ( z < ifuncmin ) then
774 : iierror = 2
775 : exit
776 : end if
777 2 : nelminxy(i) = nelminxy(i) - del - del
778 2 : z = func(nelminxy, xx, yy)
779 2 : ineval = ineval + 1
780 2 : if (present(history)) history_tmp(ineval) = Min( z, history_tmp(ineval-1) )
781 2 : if ( z < ifuncmin ) then
782 : iierror = 2
783 : exit
784 : end if
785 3 : nelminxy(i) = nelminxy(i) + del
786 : end do
787 :
788 1 : if ( iierror == 0 ) exit
789 : !
790 : ! Restart the procedure.
791 : !
792 0 : ipstart(1:n) = nelminxy(1:n)
793 0 : del = eps
794 0 : inumrestart = inumrestart + 1
795 :
796 : end do
797 :
798 1 : if (present(funcmin)) then
799 0 : funcmin = ifuncmin
800 : endif
801 1 : if (present(neval)) then
802 1 : neval = ineval
803 : endif
804 1 : if (present(numrestart)) then
805 0 : numrestart = inumrestart
806 : endif
807 1 : if (present(ierror)) then
808 0 : ierror = iierror
809 : endif
810 1 : if (present(history)) then
811 3 : allocate(history(ineval))
812 384 : history(:) = history_tmp(1:ineval)
813 : end if
814 :
815 1 : end function nelminxy
816 :
817 13 : function nelminrange_dp(func, pstart, prange, varmin, step, konvge, maxeval, &
818 : funcmin, neval, numrestart, ierror, history)
819 :
820 : implicit none
821 :
822 : INTERFACE
823 : FUNCTION func(pp)
824 : USE mo_kind, ONLY: dp
825 : IMPLICIT NONE
826 : REAL(dp), DIMENSION(:), INTENT(IN) :: pp
827 : REAL(dp) :: func
828 : END FUNCTION func
829 : END INTERFACE
830 : real(dp), intent(IN) :: pstart(:)
831 : real(dp), intent(IN) :: prange(:,:)
832 : real(dp), optional, intent(IN) :: varmin
833 : real(dp), optional, intent(IN) :: step(:)
834 : integer(i4), optional, intent(IN) :: konvge
835 : integer(i4), optional, intent(IN) :: maxeval
836 : real(dp), optional, intent(OUT) :: funcmin
837 : integer(i4), optional, intent(OUT) :: neval
838 : integer(i4), optional, intent(OUT) :: numrestart
839 : integer(i4), optional, intent(OUT) :: ierror
840 : real(dp), optional, intent(OUT), allocatable :: history(:) ! History of objective function values
841 : real(dp) :: nelminrange_dp(size(pstart))
842 :
843 : real(dp), parameter :: ccoeff = 0.5_dp
844 : real(dp), parameter :: ecoeff = 2.0_dp
845 : real(dp), parameter :: rcoeff = 1.0_dp
846 : real(dp), parameter :: eps = 0.001_dp
847 10 : real(dp) :: ipstart(size(pstart))
848 4 : real(dp) :: ifuncmin
849 10 : real(dp) :: ivarmin
850 14 : real(dp) :: istep(size(pstart))
851 : integer(i4) :: ikonvge
852 : integer(i4) :: imaxeval
853 : integer(i4) :: ineval
854 : integer(i4) :: inumrestart
855 : integer(i4) :: iierror
856 : integer(i4) :: n, nn
857 : integer(i4) :: i
858 : integer(i4) :: ihi
859 : integer(i4) :: ilo
860 : integer(i4) :: j
861 : integer(i4) :: jcount
862 : integer(i4) :: l
863 4 : real(dp) :: del
864 34 : real(dp) :: p(size(pstart),size(pstart)+1)
865 10 : real(dp) :: p2star(size(pstart))
866 14 : real(dp) :: pbar(size(pstart))
867 14 : real(dp) :: pstar(size(pstart))
868 4 : real(dp) :: rq
869 4 : real(dp) :: x
870 18 : real(dp) :: y(size(pstart)+1)
871 4 : real(dp) :: y2star
872 4 : real(dp) :: ylo
873 4 : real(dp) :: ystar
874 4 : real(dp) :: z
875 4 : real(dp) :: dn, dnn
876 22 : real(dp) :: p0(size(pstart)), y0
877 4 : real(dp), allocatable :: history_tmp(:)
878 : !
879 : ! Defaults
880 : !
881 10 : nelminrange_dp(:) = 0.5_dp * ( prange(:,1) + prange(:,2) )
882 :
883 4 : if (present(varmin)) then
884 0 : if (varmin <= 0.0_dp) stop 'Error nelmin: varmin<0'
885 : ivarmin = varmin
886 : else
887 : ivarmin = 1.0e-9_dp
888 : endif
889 : ! maximal number of function evaluations
890 4 : if (present(maxeval)) then
891 0 : if (maxeval <= 1) stop 'Error nelmin: maxeval<=1'
892 : imaxeval = maxeval
893 : else
894 : imaxeval = 1000
895 : endif
896 : ! history output
897 4 : if (present(history)) then
898 : ! worst case length
899 0 : allocate(history_tmp(imaxeval+3*size(ipstart)+1))
900 : end if
901 4 : if (present(konvge)) then
902 0 : ikonvge = konvge
903 : else
904 4 : ikonvge = imaxeval / 10
905 : endif
906 4 : if (ikonvge < 1) stop 'Error nelmin: konvg<1'
907 : !
908 4 : if (present(step)) then
909 0 : istep = step
910 : else ! if not given, deviate initial by 1%
911 4 : y0 = func(pstart)
912 10 : do i=1, size(pstart)
913 16 : p0 = pstart
914 6 : if (ne(p0(i),0.0_dp)) then
915 : ! bound to range
916 6 : p0(i) = min( prange(i,2) , max( prange(i,1) , 1.01_dp*p0(i) ) )
917 : else
918 : ! bound to range
919 0 : p0(i) = min( prange(i,2) , max( prange(i,1) , 0.01_dp ) )
920 : endif
921 10 : istep(i) = abs(func(p0) - y0)
922 : enddo
923 : endif
924 : !
925 : ! Check the input parameters.
926 : !
927 4 : if (size(prange,2) .ne. 2_i4) stop 'nelminrange_dp: range has to be array of size (:,2)'
928 4 : if (size(prange,1) .ne. size(ipstart)) stop 'nelminrange_dp: range has to be given for each dimension'
929 20 : if (.not. (all(prange(:,1) .le. pstart(:)) .and. all(prange(:,2) .ge. pstart(:)) )) then
930 0 : stop 'nelminrange_dp: starting point is not in range'
931 : end if
932 : !
933 : ! Initialization.
934 : !
935 10 : ipstart = pstart
936 4 : n = size(ipstart)
937 4 : ineval = 0_i4
938 4 : inumrestart = 0_i4
939 4 : jcount = ikonvge
940 4 : dn = real(n, dp)
941 4 : nn = n + 1
942 4 : dnn = real(nn, dp)
943 4 : del = 1.0_dp
944 4 : rq = ivarmin * dn
945 : !
946 : ! Initial or restarted loop.
947 : !
948 0 : do
949 10 : p(1:n,nn) = ipstart(1:n)
950 4 : y(nn) = func(ipstart)
951 4 : ineval = ineval + 1
952 4 : if (present(history)) history_tmp(ineval) = y(nn)
953 : !
954 : ! Define the initial simplex.
955 : !
956 10 : do j = 1, n
957 6 : x = ipstart(j)
958 : ! bound to range
959 6 : ipstart(j) = min( prange(j,2) , max( prange(j,1) , ipstart(j) + istep(j) * del ) )
960 16 : p(1:n,j) = ipstart(1:n)
961 6 : y(j) = func(ipstart)
962 6 : ineval = ineval + 1
963 6 : ipstart(j) = x
964 10 : if (present(history)) history_tmp(ineval) = Min( y(j), history_tmp(ineval-1) )
965 : end do
966 : !
967 : ! Find highest and lowest Y values. FUNCMIN = Y(IHI) indicates
968 : ! the vertex of the simplex to be replaced.
969 : !
970 14 : ilo = minloc(y(1:n+1), 1)
971 4 : ylo = y(ilo)
972 : !
973 : ! Inner loop.
974 : !
975 505 : do while ( ineval < imaxeval )
976 : !
977 : ! FUNCMIN is, of course, the HIGHEST value???
978 : !
979 1735 : ihi = maxloc(y(1:n+1), 1)
980 505 : ifuncmin = y(ihi)
981 : !
982 : ! Calculate PBAR, the centroid of the simplex vertices
983 : ! excepting the vertex with Y value FUNCMIN.
984 : !
985 1230 : do i = 1, n
986 3120 : pbar(i) = ( sum(p(i,1:n+1)) - p(i,ihi) ) / dn
987 : end do
988 : !
989 : ! Reflection through the centroid.
990 : !
991 : ! bound to range
992 1230 : pstar(1:n) = min( prange(1:n,2) , max( prange(1:n,1) , pbar(1:n) + rcoeff * ( pbar(1:n) - p(1:n,ihi) ) ) )
993 505 : ystar = func(pstar)
994 505 : ineval = ineval + 1
995 505 : if (present(history)) history_tmp(ineval) = Min( ystar, history_tmp(ineval-1) )
996 : !
997 : ! Successful reflection, so extension.
998 : !
999 505 : if ( ystar < ylo ) then
1000 :
1001 : ! bound to range
1002 220 : p2star(1:n) = min( prange(1:n,2) , max( prange(1:n,1) , pbar(1:n) + ecoeff * ( pstar(1:n) - pbar(1:n) ) ) )
1003 78 : y2star = func(p2star)
1004 78 : ineval = ineval + 1
1005 78 : if (present(history)) history_tmp(ineval) = Min( y2star, history_tmp(ineval-1) )
1006 : !
1007 : ! Retain extension or contraction.
1008 : !
1009 78 : if ( ystar < y2star ) then
1010 50 : p(1:n,ihi) = pstar(1:n)
1011 19 : y(ihi) = ystar
1012 : else
1013 170 : p(1:n,ihi) = p2star(1:n)
1014 59 : y(ihi) = y2star
1015 : end if
1016 : !
1017 : ! No extension.
1018 : !
1019 : else
1020 : l = 0
1021 1437 : do i = 1, nn
1022 1437 : if ( ystar < y(i) ) l = l + 1
1023 : end do
1024 :
1025 427 : if ( 1 < l ) then
1026 84 : p(1:n,ihi) = pstar(1:n)
1027 28 : y(ihi) = ystar
1028 : !
1029 : ! Contraction on the Y(IHI) side of the centroid.
1030 : !
1031 399 : else if ( l == 0 ) then
1032 : ! bound to range
1033 692 : p2star(1:n) = min( prange(1:n,2) , max( prange(1:n,1) , pbar(1:n) + ccoeff * ( p(1:n,ihi) - pbar(1:n) ) ) )
1034 295 : y2star = func(p2star)
1035 295 : ineval = ineval + 1
1036 295 : if (present(history)) history_tmp(ineval) = Min( y2star, history_tmp(ineval-1) )
1037 : !
1038 : ! Contract the whole simplex.
1039 : !
1040 295 : if ( y(ihi) < y2star ) then
1041 335 : do j = 1, n + 1
1042 : ! bound to range
1043 520 : p(1:n,j) = min( prange(1:n,2) , max( prange(1:n,1) , ( p(1:n,j) + p(1:n,ilo) ) * 0.5_dp ) )
1044 520 : nelminrange_dp(1:n) = p(1:n,j)
1045 230 : y(j) = func(nelminrange_dp)
1046 230 : ineval = ineval + 1
1047 335 : if (present(history)) history_tmp(ineval) = Min( y(j), history_tmp(ineval-1) )
1048 : end do
1049 335 : ilo = minloc(y(1:n+1), 1)
1050 105 : ylo = y(ilo)
1051 :
1052 105 : cycle
1053 : !
1054 : ! Retain contraction.
1055 : !
1056 : else
1057 462 : p(1:n,ihi) = p2star(1:n)
1058 190 : y(ihi) = y2star
1059 : end if
1060 : !
1061 : ! Contraction on the reflection side of the centroid.
1062 : !
1063 104 : else if ( l == 1 ) then
1064 :
1065 : ! bound to range
1066 234 : p2star(1:n) = min( prange(1:n,2) , max( prange(1:n,1) , pbar(1:n) + ccoeff * ( pstar(1:n) - pbar(1:n) ) ) )
1067 104 : y2star = func(p2star)
1068 104 : ineval = ineval + 1
1069 104 : if (present(history)) history_tmp(ineval) = Min( y2star, history_tmp(ineval-1) )
1070 : !
1071 : ! Retain reflection?
1072 : !
1073 104 : if ( y2star <= ystar ) then
1074 78 : p(1:n,ihi) = p2star(1:n)
1075 29 : y(ihi) = y2star
1076 : else
1077 156 : p(1:n,ihi) = pstar(1:n)
1078 75 : y(ihi) = ystar
1079 : end if
1080 :
1081 : end if
1082 :
1083 : end if
1084 : !
1085 : ! Check if YLO improved.
1086 : !
1087 400 : if ( y(ihi) < ylo ) then
1088 138 : ylo = y(ihi)
1089 138 : ilo = ihi
1090 : end if
1091 :
1092 400 : jcount = jcount - 1
1093 :
1094 400 : if ( 0 < jcount ) cycle
1095 : !
1096 : ! Check to see if minimum reached.
1097 : !
1098 8 : if ( ineval <= imaxeval ) then
1099 14 : jcount = ikonvge
1100 14 : x = sum ( y(1:n+1) ) / dnn
1101 14 : z = sum ( ( y(1:n+1) - x )**2 )
1102 4 : if ( z <= rq ) exit
1103 : end if
1104 :
1105 : end do
1106 : !
1107 : ! Factorial tests to check that FUNCMIN is a local minimum.
1108 : !
1109 : ! bound to range
1110 10 : nelminrange_dp(1:n) = min( prange(1:n,2) , max( prange(1:n,1) , p(1:n,ilo) ) )
1111 4 : ifuncmin = y(ilo)
1112 :
1113 4 : if ( imaxeval < ineval ) then
1114 : iierror = 2
1115 : exit
1116 : end if
1117 :
1118 10 : iierror = 0
1119 :
1120 10 : do i = 1, n
1121 6 : del = istep(i) * eps
1122 : ! bound to range: probably not necessary
1123 6 : nelminrange_dp(i) = min( prange(i,2) , max( prange(i,1) , nelminrange_dp(i) + del ) )
1124 6 : z = func(nelminrange_dp)
1125 6 : ineval = ineval + 1
1126 6 : if (present(history)) history_tmp(ineval) = Min( z, history_tmp(ineval-1) )
1127 6 : if ( z < ifuncmin ) then
1128 : iierror = 2
1129 : exit
1130 : end if
1131 : ! bound to range: probably not necessary
1132 6 : nelminrange_dp(i) = min( prange(i,2) , max( prange(i,1) , nelminrange_dp(i) - del - del ) )
1133 6 : z = func(nelminrange_dp)
1134 6 : ineval = ineval + 1
1135 6 : if (present(history)) history_tmp(ineval) = Min( z, history_tmp(ineval-1) )
1136 6 : if ( z < ifuncmin ) then
1137 : iierror = 2
1138 : exit
1139 : end if
1140 : ! bound to range: probably not necessary
1141 10 : nelminrange_dp(i) = min( prange(i,2) , max( prange(i,1) , nelminrange_dp(i) + del ) )
1142 : end do
1143 :
1144 4 : if ( iierror == 0 ) exit
1145 : !
1146 : ! Restart the procedure.
1147 : !
1148 : ! bound to range: probably not necessary
1149 0 : ipstart(1:n) = min( prange(1:n,2) , max( prange(1:n,1) , nelminrange_dp(1:n) ) )
1150 0 : del = eps
1151 0 : inumrestart = inumrestart + 1
1152 :
1153 : end do
1154 :
1155 4 : if (present(funcmin)) then
1156 0 : funcmin = ifuncmin
1157 : endif
1158 4 : if (present(neval)) then
1159 0 : neval = ineval
1160 : endif
1161 4 : if (present(numrestart)) then
1162 0 : numrestart = inumrestart
1163 : endif
1164 4 : if (present(ierror)) then
1165 0 : ierror = iierror
1166 : endif
1167 4 : if (present(history)) then
1168 0 : allocate(history(ineval))
1169 0 : history(:) = history_tmp(1:ineval)
1170 : end if
1171 :
1172 4 : end function nelminrange_dp
1173 :
1174 4 : function nelminrange_sp(func, pstart, prange, varmin, step, konvge, maxeval, &
1175 : funcmin, neval, numrestart, ierror, history)
1176 :
1177 : implicit none
1178 :
1179 : INTERFACE
1180 : FUNCTION func(pp)
1181 : USE mo_kind, ONLY: sp
1182 : IMPLICIT NONE
1183 : REAL(sp), DIMENSION(:), INTENT(IN) :: pp
1184 : REAL(sp) :: func
1185 : END FUNCTION func
1186 : END INTERFACE
1187 : real(sp), intent(IN) :: pstart(:)
1188 : real(sp), intent(IN) :: prange(:,:)
1189 : real(sp), optional, intent(IN) :: varmin
1190 : real(sp), optional, intent(IN) :: step(:)
1191 : integer(i4), optional, intent(IN) :: konvge
1192 : integer(i4), optional, intent(IN) :: maxeval
1193 : real(sp), optional, intent(OUT) :: funcmin
1194 : integer(i4), optional, intent(OUT) :: neval
1195 : integer(i4), optional, intent(OUT) :: numrestart
1196 : integer(i4), optional, intent(OUT) :: ierror
1197 : real(sp), optional, intent(OUT), allocatable :: history(:) ! History of objective function values
1198 : real(sp) :: nelminrange_sp(size(pstart))
1199 :
1200 : real(sp), parameter :: ccoeff = 0.5_sp
1201 : real(sp), parameter :: ecoeff = 2.0_sp
1202 : real(sp), parameter :: rcoeff = 1.0_sp
1203 : real(sp), parameter :: eps = 0.001_sp
1204 0 : real(sp) :: ipstart(size(pstart))
1205 0 : real(sp) :: ifuncmin
1206 0 : real(sp) :: ivarmin
1207 0 : real(sp) :: istep(size(pstart))
1208 : integer(i4) :: ikonvge
1209 : integer(i4) :: imaxeval
1210 : integer(i4) :: ineval
1211 : integer(i4) :: inumrestart
1212 : integer(i4) :: iierror
1213 : integer(i4) :: n, nn
1214 : integer(i4) :: i
1215 : integer(i4) :: ihi
1216 : integer(i4) :: ilo
1217 : integer(i4) :: j
1218 : integer(i4) :: jcount
1219 : integer(i4) :: l
1220 0 : real(sp) :: del
1221 0 : real(sp) :: p(size(pstart),size(pstart)+1)
1222 0 : real(sp) :: p2star(size(pstart))
1223 0 : real(sp) :: pbar(size(pstart))
1224 0 : real(sp) :: pstar(size(pstart))
1225 0 : real(sp) :: rq
1226 0 : real(sp) :: x
1227 0 : real(sp) :: y(size(pstart)+1)
1228 0 : real(sp) :: y2star
1229 0 : real(sp) :: ylo
1230 0 : real(sp) :: ystar
1231 0 : real(sp) :: z
1232 0 : real(sp) :: dn, dnn
1233 0 : real(sp) :: p0(size(pstart)), y0
1234 0 : real(sp), allocatable :: history_tmp(:)
1235 : !
1236 : ! Defaults
1237 : !
1238 0 : nelminrange_sp(:) = 0.5_sp * ( prange(:,1) + prange(:,2) )
1239 :
1240 0 : if (present(varmin)) then
1241 0 : if (varmin <= 0.0_sp) stop 'Error nelmin: varmin<0'
1242 : ivarmin = varmin
1243 : else
1244 : ivarmin = 1.0e-9_sp
1245 : endif
1246 : ! maximal number of function evaluations
1247 0 : if (present(maxeval)) then
1248 0 : if (maxeval <= 1) stop 'Error nelmin: maxeval<=1'
1249 : imaxeval = maxeval
1250 : else
1251 : imaxeval = 1000
1252 : endif
1253 : ! history output
1254 0 : if (present(history)) then
1255 : ! worst case length
1256 0 : allocate(history_tmp(imaxeval+3*size(ipstart)+1))
1257 : end if
1258 0 : if (present(konvge)) then
1259 0 : ikonvge = konvge
1260 : else
1261 0 : ikonvge = imaxeval / 10
1262 : endif
1263 0 : if (ikonvge < 1) stop 'Error nelmin: konvg<1'
1264 : !
1265 0 : if (present(step)) then
1266 0 : istep = step
1267 : else ! if not given, deviate initial by 1%
1268 0 : y0 = func(pstart)
1269 0 : do i=1, size(pstart)
1270 0 : p0 = pstart
1271 0 : if (ne(p0(i),0.0_sp)) then
1272 : ! bound to range
1273 0 : p0(i) = min( prange(i,2) , max( prange(i,1) , 1.01_sp*p0(i) ) )
1274 : else
1275 : ! bound to range
1276 0 : p0(i) = min( prange(i,2) , max( prange(i,1) , 0.01_sp ) )
1277 : endif
1278 0 : istep(i) = abs(func(p0) - y0)
1279 : enddo
1280 : endif
1281 : !
1282 : ! Check the input parameters.
1283 : !
1284 0 : if (size(prange,2) .ne. 2_i4) stop 'nelminrange_sp: range has to be array of size (:,2)'
1285 0 : if (size(prange,1) .ne. size(ipstart)) stop 'nelminrange_sp: range has to be given for each dimension'
1286 0 : if (.not. (all(prange(:,1) .le. pstart(:)) .and. all(prange(:,2) .ge. pstart(:)) )) then
1287 0 : stop 'nelminrange_sp: starting point is not in range'
1288 : end if
1289 : !
1290 : ! Initialization.
1291 : !
1292 0 : ipstart = pstart
1293 0 : n = size(ipstart)
1294 0 : ineval = 0_i4
1295 0 : inumrestart = 0_i4
1296 0 : jcount = ikonvge
1297 0 : dn = real(n, sp)
1298 0 : nn = n + 1
1299 0 : dnn = real(nn, sp)
1300 0 : del = 1.0_sp
1301 0 : rq = ivarmin * dn
1302 : !
1303 : ! Initial or restarted loop.
1304 : !
1305 0 : do
1306 0 : p(1:n,nn) = ipstart(1:n)
1307 0 : y(nn) = func(ipstart)
1308 0 : ineval = ineval + 1
1309 0 : if (present(history)) history_tmp(ineval) = y(nn)
1310 : !
1311 : ! Define the initial simplex.
1312 : !
1313 0 : do j = 1, n
1314 0 : x = ipstart(j)
1315 : ! bound to range
1316 0 : ipstart(j) = min( prange(j,2) , max( prange(j,1) , ipstart(j) + istep(j) * del ) )
1317 0 : p(1:n,j) = ipstart(1:n)
1318 0 : y(j) = func(ipstart)
1319 0 : ineval = ineval + 1
1320 0 : ipstart(j) = x
1321 0 : if (present(history)) history_tmp(ineval) = Min( y(j), history_tmp(ineval-1) )
1322 : end do
1323 : !
1324 : ! Find highest and lowest Y values. FUNCMIN = Y(IHI) indicates
1325 : ! the vertex of the simplex to be replaced.
1326 : !
1327 0 : ilo = minloc(y(1:n+1), 1)
1328 0 : ylo = y(ilo)
1329 : !
1330 : ! Inner loop.
1331 : !
1332 0 : do while ( ineval < imaxeval )
1333 : !
1334 : ! FUNCMIN is, of course, the HIGHEST value???
1335 : !
1336 0 : ihi = maxloc(y(1:n+1), 1)
1337 0 : ifuncmin = y(ihi)
1338 : !
1339 : ! Calculate PBAR, the centroid of the simplex vertices
1340 : ! excepting the vertex with Y value FUNCMIN.
1341 : !
1342 0 : do i = 1, n
1343 0 : pbar(i) = ( sum(p(i,1:n+1)) - p(i,ihi) ) / dn
1344 : end do
1345 : !
1346 : ! Reflection through the centroid.
1347 : !
1348 : ! bound to range
1349 0 : pstar(1:n) = min( prange(1:n,2) , max( prange(1:n,1) , pbar(1:n) + rcoeff * ( pbar(1:n) - p(1:n,ihi) ) ) )
1350 0 : ystar = func(pstar)
1351 0 : ineval = ineval + 1
1352 0 : if (present(history)) history_tmp(ineval) = Min( ystar, history_tmp(ineval-1) )
1353 : !
1354 : ! Successful reflection, so extension.
1355 : !
1356 0 : if ( ystar < ylo ) then
1357 :
1358 : ! bound to range
1359 0 : p2star(1:n) = min( prange(1:n,2) , max( prange(1:n,1) , pbar(1:n) + ecoeff * ( pstar(1:n) - pbar(1:n) ) ) )
1360 0 : y2star = func(p2star)
1361 0 : ineval = ineval + 1
1362 0 : if (present(history)) history_tmp(ineval) = Min( y2star, history_tmp(ineval-1) )
1363 : !
1364 : ! Retain extension or contraction.
1365 : !
1366 0 : if ( ystar < y2star ) then
1367 0 : p(1:n,ihi) = pstar(1:n)
1368 0 : y(ihi) = ystar
1369 : else
1370 0 : p(1:n,ihi) = p2star(1:n)
1371 0 : y(ihi) = y2star
1372 : end if
1373 : !
1374 : ! No extension.
1375 : !
1376 : else
1377 : l = 0
1378 0 : do i = 1, nn
1379 0 : if ( ystar < y(i) ) l = l + 1
1380 : end do
1381 :
1382 0 : if ( 1 < l ) then
1383 0 : p(1:n,ihi) = pstar(1:n)
1384 0 : y(ihi) = ystar
1385 : !
1386 : ! Contraction on the Y(IHI) side of the centroid.
1387 : !
1388 0 : else if ( l == 0 ) then
1389 : ! bound to range
1390 0 : p2star(1:n) = min( prange(1:n,2) , max( prange(1:n,1) , pbar(1:n) + ccoeff * ( p(1:n,ihi) - pbar(1:n) ) ) )
1391 0 : y2star = func(p2star)
1392 0 : ineval = ineval + 1
1393 0 : if (present(history)) history_tmp(ineval) = Min( y2star, history_tmp(ineval-1) )
1394 : !
1395 : ! Contract the whole simplex.
1396 : !
1397 0 : if ( y(ihi) < y2star ) then
1398 0 : do j = 1, n + 1
1399 : ! bound to range
1400 0 : p(1:n,j) = min( prange(1:n,2) , max( prange(1:n,1) , ( p(1:n,j) + p(1:n,ilo) ) * 0.5_sp ) )
1401 0 : nelminrange_sp(1:n) = p(1:n,j)
1402 0 : y(j) = func(nelminrange_sp)
1403 0 : ineval = ineval + 1
1404 0 : if (present(history)) history_tmp(ineval) = Min( y(j), history_tmp(ineval-1) )
1405 : end do
1406 0 : ilo = minloc(y(1:n+1), 1)
1407 0 : ylo = y(ilo)
1408 :
1409 0 : cycle
1410 : !
1411 : ! Retain contraction.
1412 : !
1413 : else
1414 0 : p(1:n,ihi) = p2star(1:n)
1415 0 : y(ihi) = y2star
1416 : end if
1417 : !
1418 : ! Contraction on the reflection side of the centroid.
1419 : !
1420 0 : else if ( l == 1 ) then
1421 :
1422 : ! bound to range
1423 0 : p2star(1:n) = min( prange(1:n,2) , max( prange(1:n,1) , pbar(1:n) + ccoeff * ( pstar(1:n) - pbar(1:n) ) ) )
1424 0 : y2star = func(p2star)
1425 0 : ineval = ineval + 1
1426 0 : if (present(history)) history_tmp(ineval) = Min( y2star, history_tmp(ineval-1) )
1427 : !
1428 : ! Retain reflection?
1429 : !
1430 0 : if ( y2star <= ystar ) then
1431 0 : p(1:n,ihi) = p2star(1:n)
1432 0 : y(ihi) = y2star
1433 : else
1434 0 : p(1:n,ihi) = pstar(1:n)
1435 0 : y(ihi) = ystar
1436 : end if
1437 :
1438 : end if
1439 :
1440 : end if
1441 : !
1442 : ! Check if YLO improved.
1443 : !
1444 0 : if ( y(ihi) < ylo ) then
1445 0 : ylo = y(ihi)
1446 0 : ilo = ihi
1447 : end if
1448 :
1449 0 : jcount = jcount - 1
1450 :
1451 0 : if ( 0 < jcount ) cycle
1452 : !
1453 : ! Check to see if minimum reached.
1454 : !
1455 0 : if ( ineval <= imaxeval ) then
1456 0 : jcount = ikonvge
1457 0 : x = sum ( y(1:n+1) ) / dnn
1458 0 : z = sum ( ( y(1:n+1) - x )**2 )
1459 0 : if ( z <= rq ) exit
1460 : end if
1461 :
1462 : end do
1463 : !
1464 : ! Factorial tests to check that FUNCMIN is a local minimum.
1465 : !
1466 : ! bound to range
1467 0 : nelminrange_sp(1:n) = min( prange(1:n,2) , max( prange(1:n,1) , p(1:n,ilo) ) )
1468 0 : ifuncmin = y(ilo)
1469 :
1470 0 : if ( imaxeval < ineval ) then
1471 : iierror = 2
1472 : exit
1473 : end if
1474 :
1475 0 : iierror = 0
1476 :
1477 0 : do i = 1, n
1478 0 : del = istep(i) * eps
1479 : ! bound to range: probably not necessary
1480 0 : nelminrange_sp(i) = min( prange(i,2) , max( prange(i,1) , nelminrange_sp(i) + del ) )
1481 0 : z = func(nelminrange_sp)
1482 0 : ineval = ineval + 1
1483 0 : if (present(history)) history_tmp(ineval) = Min( z, history_tmp(ineval-1) )
1484 0 : if ( z < ifuncmin ) then
1485 : iierror = 2
1486 : exit
1487 : end if
1488 : ! bound to range: probably not necessary
1489 0 : nelminrange_sp(i) = min( prange(i,2) , max( prange(i,1) , nelminrange_sp(i) - del - del ) )
1490 0 : z = func(nelminrange_sp)
1491 0 : ineval = ineval + 1
1492 0 : if (present(history)) history_tmp(ineval) = Min( z, history_tmp(ineval-1) )
1493 0 : if ( z < ifuncmin ) then
1494 : iierror = 2
1495 : exit
1496 : end if
1497 : ! bound to range: probably not necessary
1498 0 : nelminrange_sp(i) = min( prange(i,2) , max( prange(i,1) , nelminrange_sp(i) + del ) )
1499 : end do
1500 :
1501 0 : if ( iierror == 0 ) exit
1502 : !
1503 : ! Restart the procedure.
1504 : !
1505 : ! bound to range: probably not necessary
1506 0 : ipstart(1:n) = min( prange(1:n,2) , max( prange(1:n,1) , nelminrange_sp(1:n) ) )
1507 0 : del = eps
1508 0 : inumrestart = inumrestart + 1
1509 :
1510 : end do
1511 :
1512 0 : if (present(funcmin)) then
1513 0 : funcmin = ifuncmin
1514 : endif
1515 0 : if (present(neval)) then
1516 0 : neval = ineval
1517 : endif
1518 0 : if (present(numrestart)) then
1519 0 : numrestart = inumrestart
1520 : endif
1521 0 : if (present(ierror)) then
1522 0 : ierror = iierror
1523 : endif
1524 0 : if (present(history)) then
1525 0 : allocate(history(ineval))
1526 0 : history(:) = history_tmp(1:ineval)
1527 : end if
1528 :
1529 0 : end function nelminrange_sp
1530 :
1531 :
1532 : ! ------------------------------------------------------------------
1533 :
1534 : END MODULE mo_nelmin
|
