0.6.2-dev0
FORCES
FORtran lib for Comp. Env. Sys.
Loading...
Searching...
No Matches
mo_sce.F90
Go to the documentation of this file.
1!> \file mo_sce.f90
2!> \brief \copybrief mo_sce
3!> \details \copydetails mo_sce
4
5!> \brief Shuffled Complex Evolution optimization algorithm.
6!> \details Optimization algorithm using Shuffled Complex Evolution strategy.
7!! Original version 2.1 of Qingyun Duan (1992) rewritten in Fortran 90.
8!> \authors Juliane Mai
9!> \date Feb 2013
10!> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
11!! FORCES is released under the LGPLv3+ license \license_note
12MODULE mo_sce
13
14 IMPLICIT NONE
15
16 PUBLIC :: sce ! sce optimization
17
18 ! ------------------------------------------------------------------
19
20 ! NAME
21 ! sce
22
23 ! PURPOSE
24 !> \brief Shuffled Complex Evolution (SCE) algorithm for global optimization.
25 !
26 !> \details Shuffled Complex Evolution method for global optimization -- version 2.1\n
27 !! \n
28 !! by Qingyun Duan\n
29 !! Department of Hydrology & Water Resources\n
30 !! University of Arizona, Tucson, AZ 85721\n
31 !! (602) 621-9360, email: duan@hwr.arizona.edu\n
32 !! \n
33 !! Written by Qingyun Duan, Oct 1990.\n
34 !! Revised by Qingyun Duan, Aug 1991.\n
35 !! Revised by Qingyun Duan, Apr 1992.\n
36 !! \n
37 !! Re-written by Juliane Mai, Feb 2013.\n
38 !! \n
39 !! \b Statement by Qingyun Duan:\n
40 !! \n
41 !! This general purpose global optimization program is developed at
42 !! the Department of Hydrology & Water Resources of the University
43 !! of Arizona. Further information regarding the SCE-UA method can
44 !! be obtained from Dr. Q. Duan, Dr. S. Sorooshian or Dr. V.K. Gupta
45 !! at the address and phone number listed above. We request all
46 !! users of this program make proper reference to the paper entitled
47 !! 'Effective and Efficient Global Optimization for Conceptual
48 !! Rainfall-runoff Models' by Duan, Q., S. Sorooshian, and V.K. Gupta,
49 !! Water Resources Research, Vol 28(4), pp.1015-1031, 1992.\n
50 !! \n
51 !! The function to be minimized is the first argument of DDS and must be defined as \n
52 !! \code{.f90}
53 !! function functn(p)
54 !! use mo_kind, only: dp
55 !! implicit none
56 !! real(dp), dimension(:), intent(in) :: p
57 !! real(dp) :: functn
58 !! end function functn
59 !! \endcode
60 !!
61 !! \b Example
62 !! \code{.f90}
63 !! use mo_opt_functions, only: griewank
64 !!
65 !! prange(:,1) = (/ -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0 /)
66 !! prange(:,2) = (/ 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0 /)
67 !! pini = (/ -.226265E+01, -.130187E+01, -.151219E+01, 0.133983E+00, 0.988159E+00, &
68 !! -.495074E+01, -.126574E+02, 0.572684E+00, 0.303864E+01, 0.343031E+01 /)
69 !!
70 !! opt = sce(griewank, pini, prange)
71 !! \endcode
72 !! -> see also example in test directory
73 !!
74 !! \b Literature
75 !! 1. Duan, Q., S. Sorooshian, and V.K. Gupta -
76 !! Effective and Efficient Global Optimization for Conceptual Rainfall-runoff Models
77 !! Water Resources Research, Vol 28(4), pp.1015-1031, 1992.
78 !!
79 !> \param[in] "real(dp) :: functn(p)" Function on which to search the optimum
80 !> \param[in] "real(dp) :: pini(:)" inital value of decision variables
81 !> \param[in] "real(dp) :: prange(size(pini),2)" Min/max range of decision variables
82 !
83 ! INTENT(INOUT)
84 ! None
85
86 ! INTENT(OUT)
87 ! None
88 !
89 ! INTENT(IN), OPTIONAL
90 !> \param[in] "integer(i8), optional :: mymaxn" max no. of trials allowed before optimization is terminated\n
91 !> DEFAULT: 1000_i8
92 !> \param[in] "logical, optional :: mymaxit" maximization (.true.) or minimization (.false.) of function\n
93 !> DEFAULT: false
94 !> \param[in] "integer(i4), optional :: mykstop" number of shuffling loops in which the criterion value must
95 !> change by given percentage before optimiz. is terminated\n
96 !> DEFAULT: 10_i4
97 !> \param[in] "real(dp), optional :: mypcento" percentage by which the criterion value must change in
98 !> given number of shuffling loops\n
99 !> DEFAULT: 0.0001_dp
100 !> \param[in] "real(dp), optional :: mypeps" optimization is terminated if volume of complex has
101 !> converged to given percentage of feasible space\n
102 !> DEFAULT: 0.001_dp
103 !> \param[in] "integer(i8), optional :: myseed" initial random seed\n
104 !> DEFAULT: get_timeseed
105 !> \param[in] "integer(i4), optional :: myngs" number of complexes in the initial population\n
106 !> DEFAULT: 2_i4
107 !> \param[in] "integer(i4), optional :: mynpg" number of points in each complex\n
108 !> DEFAULT: 2*n+1
109 !> \param[in] "integer(i4), optional :: mynps" number of points in a sub-complex\n
110 !> DEFAULT: n+1
111 !> \param[in] "integer(i4), optional :: mynspl" number of evolution steps allowed for each complex before
112 !> complex shuffling\n
113 !> DEFAULT: 2*n+1
114 !> \param[in] "integer(i4), optional :: mymings" minimum number of complexes required, if the number of
115 !> complexes is allowed to reduce as the
116 !> optimization proceeds\n
117 !> DEFAULT: ngs = number of complexes in initial population
118 !> \param[in] "integer(i4), optional :: myiniflg" flag on whether to include the initial point in population\n
119 !> 0, not included\n
120 !> 1, included (DEFAULT)
121 !> \param[in] "integer(i4), optional :: myprint" flag for controlling print-out after each shuffling loop\n
122 !> 0, print information on the best point of the population\n
123 !> 1, print information on every point of the population\n
124 !> 2, no printing (DEFAULT)
125 !> 3, same as 0 but print progress '.' on every function call
126 !> 4, same as 1 but print progress '.' on every function call
127 !> \param[in] "logical, optional :: mymask(size(pini))"
128 !> parameter included in optimization (true) or discarded (false)
129 !> DEFAULT: .true.
130 !> \param[in] "real(dp), optional :: myalpha" parameter for reflection of points in complex\n
131 !> DEFAULT: 0.8_dp
132 !> \param[in] "real(dp), optional :: mybeta" parameter for contraction of points in complex\n
133 !> DEFAULT: 0.45_dp
134 !> \param[in] "character(len=*), optional :: tmp_file" if given: write results after each evolution loop
135 !> to temporal output file of that name\n
136 !> \# of headlines: 7\n
137 !> format: '\# nloop icall ngs1 bestf worstf ... \n
138 !> ... gnrng (bestx(j),j=1,nn)'
139 !> \param[in] "character(len=*), optional :: popul_file" if given: write whole population to file of that name\n
140 !> \# of headlines: 1 \n
141 !> format: \#_evolution_loop, xf(i), (x(i,j),j=1,nn)\n
142 !> total number of lines written <= neval <= mymaxn\n
143 !> \param[in] "logical, optional :: popul_file_append" if true, append to existing population file (default: false)\n
144 !> \param[in] "logical, optional :: parallel" sce runs in parallel (true) or not (false)
145 !> parallel sce should only be used if model/ objective
146 !> is not parallel
147 !> DEFAULT: .false.
148 !> \param[in] "logical, optional :: restart" if .true.: restart former sce run from restart_file
149 !> \param[in] "character(len=*), optional :: restart_file" file name for read/write of restart file
150 !> (default: mo_sce.restart)
151 !
152 ! INTENT(INOUT), OPTIONAL
153 ! None
154 !
155 ! INTENT(OUT), OPTIONAL
156 !> \param[out] "real(dp), optional :: bestf" the best value of the function.
157 !> \param[out] "integer(i8), optional :: neval" number of function evaluations needed.
158 !> \param[out] "real(dp), optional, allocatable :: history(:)" the history of best function values,
159 !> history(neval)=bestf
160 !
161 ! RETURN
162 !> \return real(dp) :: bestx(size(pini)) &mdash; The parameters of the point which is estimated
163 !! to minimize/maximize the function.
164
165 ! RESTRICTIONS
166 !> \note Maximal number of parameters is 1000.\n
167 !! SCE is OpenMP enabled on the loop over the complexes.\n
168 !! OMP_NUM_THREADS > 1 does not give reproducible results even when seeded!
169 !!
170 ! HISTORY
171 !> \author Juliane Mai, Matthias Cuntz
172 !> \date Feb 2013
173 !! - first version
174 !> \date Jul 2013
175 !! - OpenMP
176 !! - NaN and Inf in objective function
177 !> \date Oct 2013
178 !! - added peps as optional argument
179 !! - allow for masked parameters
180 !! - write population to file --> popul_file
181 !! - write intermediate results to file --> tmp_file
182 !! - flag parallel introduced
183 !> \date Nov 2013
184 !! - progress dots (MC)
185 !! - use iso_fortran_env
186 !! - treat functn=NaN as worse function value in cce
187 !> \date May 2014
188 !! - sort -> orderpack (MC)
189 !! - popul_file_append (MC)
190 !! - sort with NaNs (MC)
191 !> \date Feb 2015
192 !! - restart (MC, JM)
193 !> \date Mar 2015
194 !! - use is_finite and special_value from mo_utils (MC)
195 !> \date Apr 2015
196 !! - handling of array-like variables in restart-namelists (JM)
197
198 ! ------------------------------------------------------------------
199
200 PRIVATE
201
202 ! chkcst ! check if a trial point satisfies all constraints.
203 ! comp ! reduces the number of complexes and points in a population
204 ! parstt ! checks for parameter convergence
205 ! getpnt ! generated a new point within its range
206 ! cce ! generates new point(s) from sub-complex
207
208 ! ------------------------------------------------------------------
209
210CONTAINS
211
212 function sce(eval, functn, pini, prange, & ! IN
213 mymaxn, mymaxit, mykstop, mypcento, mypeps, myseed, & ! Optional IN
214 myngs, mynpg, mynps, mynspl, mymings, myiniflg, myprint, & ! Optional IN
215 mymask, myalpha, mybeta, & ! Optional IN
216 tmp_file, popul_file, & ! Optional IN
217 popul_file_append, & ! Optional IN
218 parallel, restart, restart_file, & ! OPTIONAL IN
219#ifdef MPI
220 comm, &
221#endif
222 bestf, neval, history & ! Optional OUT
223 ) result(bestx)
224
225 use mo_kind, only : i4, i8, dp
226 use mo_orderpack, only : sort
227 use mo_string_utils, only : num2str, compress
228 use mo_xor4096, only : get_timeseed, n_save_state, xor4096, xor4096g
229 !$ use omp_lib, only: OMP_GET_THREAD_NUM, OMP_GET_NUM_THREADS
230 use iso_fortran_env, only : output_unit, error_unit
231 use mo_utils, only : is_finite, special_value
232 use mo_optimization_utils, only : eval_interface, objective_interface
233#ifdef MPI
234 use mpi_f08
235#endif
236
237 implicit none
238 !
239 procedure(eval_interface), INTENT(IN), POINTER :: eval
240 procedure(objective_interface), intent(in), pointer :: functn
241 real(dp), dimension(:), intent(in) :: pini ! initial parameter set
242 real(dp), dimension(:, :), intent(in) :: prange ! lower and upper bound on parameters
243 !
244 ! algorithmic control & convergence parameter
245 integer(i8), optional, intent(in) :: mymaxn ! max no. of trials allowed before optimization is terminated
246 ! ! DEFAULT: 1000_i8
247 logical, optional, intent(in) :: mymaxit ! maximization (.true.) or minimization (.false.) of function
248 ! ! DEFAULT: false
249 integer(i4), optional, intent(in) :: mykstop ! number of shuffling loops in which the criterion value must
250 ! ! change by given percentage before optimiz. is terminated
251 ! ! DEFAULT: 10_i4
252 real(dp), optional, intent(in) :: mypcento ! percentage by which the criterion value must change in
253 ! ! given number of shuffling loops
254 ! ! DEFAULT: 0.0001_dp
255 real(dp), optional, intent(in) :: mypeps ! optimization is terminated if volume of complex has
256 ! ! converged to given percentage of feasible space
257 ! ! DEFAULT: 0.001_dp
258 integer(i8), optional, intent(in) :: myseed ! initial random seed
259 ! ! DEFAULT: get_timeseed
260 integer(i4), optional, intent(in) :: myngs ! number of complexes in the initial population
261 ! ! DEFAULT: 2_i4
262 integer(i4), optional, intent(in) :: mynpg ! number of points in each complex
263 ! ! DEFAULT: 2*n+1
264 integer(i4), optional, intent(in) :: mynps ! number of points in a sub-complex
265 ! ! DEFAULT: n+1
266 integer(i4), optional, intent(in) :: mynspl ! number of evolution steps allowed for each complex before
267 ! ! complex shuffling
268 ! ! DEFAULT: 2*n+1
269 integer(i4), optional, intent(in) :: mymings ! minimum number of complexes required, if the number of
270 ! ! complexes is allowed to reduce as the
271 ! ! optimization proceeds
272 ! ! DEFAULT: ngs = number of complexes in initial population
273 integer(i4), optional, intent(in) :: myiniflg ! flag on whether to include the initial point in population
274 ! ! 0, not included
275 ! ! 1, included (DEFAULT)
276 integer(i4), optional, intent(in) :: myprint ! flag for controlling print-out after each shuffling loop
277 ! ! 0, print information on the best point of the population
278 ! ! 1, print information on every point of the population
279 ! ! 2, no printing (DEFAULT)
280 logical, optional, intent(in), dimension(size(pini, 1)) :: mymask ! parameter included in optimization(true) or discarded(false)
281 ! ! DEFAULT: .true.
282 real(dp), optional, intent(in) :: myalpha ! parameter for reflection of points in complex
283 ! ! DEFAULT: 0.8_dp
284 real(dp), optional, intent(in) :: mybeta ! parameter for contraction of points in complex
285 ! ! DEFAULT: 0.45_dp
286 character(len = *), optional, intent(in) :: tmp_file ! file for temporal output: write results after evolution loop
287 ! ! # of headlines: 7
288 ! ! format: '# nloop icall ngs1 bestf worstf ...
289 ! ! ... gnrng (bestx(j),j=1,nn)'
290 character(len = *), optional, intent(in) :: popul_file ! file for temporal output: writes whole population
291 ! ! # of headlines: 1
292 ! ! format: #_evolution_loop, xf(i), (x(i,j),j=1,nn)
293 ! ! total number of lines written <= neval <= mymaxn
294 logical, optional, intent(in) :: popul_file_append ! if true, append to popul_file
295 logical, optional, intent(in) :: parallel ! sce runs in parallel (true) or not (false)
296 ! ! parallel sce should only be used if model/ objective
297 ! ! is not parallel
298 ! ! DEAFULT: .false.
299 logical, optional, intent(in) :: restart ! if .true., start from restart file
300 character(len = *), optional, intent(in) :: restart_file ! restart file name (default: mo_sce.restart)
301#ifdef MPI
302 type(mpi_comm), optional, intent(in) :: comm ! MPI communicator
303#endif
304 real(dp), optional, intent(out) :: bestf ! function value of bestx(.)
305 integer(i8), optional, intent(out) :: neval ! number of function evaluations
306 real(dp), optional, dimension(:), allocatable, intent(out) :: history ! history of best function values after each iteration
307 real(dp), dimension(size(pini, 1)) :: bestx ! best point at current shuffling loop (is RETURN)
308 !
309 ! optionals transfer
310 real(dp), dimension(size(pini, 1)) :: bl ! lower bound on parameters
311 real(dp), dimension(1000) :: dummy_bl ! dummy lower bound (only for namelist)
312 real(dp), dimension(size(pini, 1)) :: bu ! upper bound on parameters
313 real(dp), dimension(1000) :: dummy_bu ! dummy upper bound (only for namelist)
314 integer(i8) :: maxn ! max no. of trials allowed before optimization is terminated
315 logical :: maxit ! minimization (false) or maximization (true)
316 integer(i4) :: kstop ! number of shuffling loops in which the criterion value
317 ! ! must change
318 real(dp) :: pcento ! percentage by which the criterion value must change
319 real(dp) :: peps ! optimization is terminated if volume of complex has
320 ! ! converged to given percentage of feasible space
321 integer(i4) :: ngs ! number of complexes in the initial population
322 integer(i4) :: npg ! number of points in each complex
323 integer(i4) :: nps ! number of points in a sub-complex
324 integer(i4) :: nspl ! number of evolution steps allowed for each complex before
325 ! ! complex shuffling
326 integer(i4) :: mings ! minimum number of complexes required
327 integer(i4) :: iniflg ! flag on whether to include the initial point in population
328 integer(i4) :: iprint ! flag for controlling print-out after each shuffling loop
329 logical :: idot ! controls progress report .
330 logical, dimension(size(pini, 1)) :: maskpara ! mask(i) = .true. --> parameter i will be optimized
331 ! ! mask(i) = .false. --> parameter i will not be optimized
332 logical, dimension(1000) :: dummy_maskpara ! dummy maskpara (only for namelist)
333 real(dp) :: alpha ! parameter for reflection of points in complex
334 real(dp) :: beta ! parameter for contraction of points in complex
335 logical :: itmp_file ! if temporal results wanted
336 character(len = 1024) :: istmp_file ! local copy of file for temporal results output
337 logical :: ipopul_file ! if temporal population wanted
338 character(len = 1024) :: ispopul_file ! local copy of file for temporal population output
339
340 real(dp), dimension(:), allocatable :: history_tmp ! history of best function values after each iteration
341 real(dp), dimension(:), allocatable :: ihistory_tmp ! local history for OpenMP
342 real(dp) :: bestf_tmp ! function value of bestx(.)
343 logical :: parall ! if sce is used parallel or not
344 logical :: irestart ! if restart wanted
345 character(len = 1024) :: isrestart_file ! local copy of restart file name
346 !
347 ! local variables
348 integer(i4) :: nopt ! number of parameters to be optimized
349 integer(i4) :: nn ! total number of parameters
350 integer(i4) :: npt ! total number of points in initial population (npt=ngs*npg)
351 real(dp) :: fpini ! function value at initial point
352 real(dp), dimension(:, :), allocatable :: x ! coordinates of points in the population
353 real(dp), dimension(:), allocatable :: xf ! function values of x
354 real(dp), dimension(size(pini, 1)) :: xx ! coordinates of a single point in x
355 real(dp), dimension(1000) :: dummy_xx ! dummy xx (only for namelist)
356 real(dp), dimension(:, :), allocatable :: cx ! coordinates of points in a complex
357 real(dp), dimension(:), allocatable :: cf ! function values of cx
358 real(dp), dimension(:, :), allocatable :: s ! coordinates of points in the current sub-complex simplex
359 real(dp), dimension(:), allocatable :: sf ! function values of s
360 real(dp), dimension(:), allocatable :: worstx ! worst point at current shuffling loop
361 real(dp) :: worstf ! function value of worstx(.)
362 real(dp), dimension(:), allocatable :: xnstd ! standard deviation of parameters in the population
363 real(dp) :: gnrng ! normalized geometric mean of parameter ranges
364 integer(i4), dimension(:), allocatable :: lcs ! indices locating position of s(.,.) in x(.,.)
365 real(dp), dimension(:), allocatable :: bound ! length of range of ith variable being optimized
366 real(dp), dimension(:), allocatable :: unit ! standard deviation of initial parameter distributions
367 integer(i4) :: ngs1 ! number of complexes in current population
368 integer(i4) :: ngs2 ! number of complexes in last population
369 real(dp), dimension(:), allocatable :: criter ! vector containing the best criterion values of the last
370 ! ! (kstop+1) shuffling loops
371 integer(i4) :: ipcnvg ! flag indicating whether parameter convergence is reached
372 ! ! (i.e., check if gnrng is less than 0.001)
373 ! ! 0, parameter convergence not satisfied
374 ! ! 1, parameter convergence satisfied
375 integer(i4) :: nloop
376 integer(i4) :: loop
377 integer(i4) :: ii
378 integer(i4) :: jj
379 integer(i4) :: kk
380 integer(i4) :: lpos
381 logical :: lpos_ok ! for selction of points based on triangular
382 ! ! probability distribution
383 integer(i4) :: npt1
384 integer(i8) :: icall ! counter for function evaluations
385 integer(i8) :: icall_merk, iicall ! local icall for OpenMP
386 integer(i4) :: igs
387 integer(i4) :: k1, k2
388 real(dp) :: denomi ! for checking improvement of last steps
389 real(dp) :: timeou ! for checking improvement of last steps
390 character(4), dimension(:), allocatable :: xname ! parameter names: "p1", "p2", "p3", ...
391 character(512) :: format_str1 ! format string
392 character(512) :: format_str2 ! format string
393 ! for random numbers
394 real(dp) :: rand ! random number
395 integer(i8), dimension(:), allocatable :: iseed ! initial random seed
396 integer(i8), dimension(:, :), allocatable :: save_state_unif ! save state of uniform stream
397 integer(i8), dimension(:, :), allocatable :: save_state_gauss ! save state of gaussian stream
398 real(dp), dimension(:), allocatable :: rand_tmp ! random number
399 integer(i4) :: ithread, n_threads ! OMP or not
400 integer(i8), dimension(n_save_state) :: itmp
401 real(dp), dimension(:, :), allocatable :: xtmp ! tmp array for complex reduction
402 real(dp), dimension(:), allocatable :: ftmp ! %
403 real(dp) :: large ! for treating NaNs
404 logical :: ipopul_file_append
405 integer(i4) :: nonan ! # of non-NaN in history_tmp
406 real(dp), dimension(:), allocatable :: htmp ! tmp storage for history_tmp
407#ifdef MPI
408 integer(i4) :: ierror
409 logical :: do_obj_loop
410 integer(i4) :: iproc, nproc
411#endif
412
413 namelist /restartnml1/ &
414 bestx, dummy_bl, dummy_bu, maxn, maxit, kstop, pcento, peps, ngs, npg, nps, nspl, &
415 mings, iniflg, iprint, idot, dummy_maskpara, alpha, beta, bestf_tmp, &
416 parall, nopt, nn, npt, fpini, dummy_xx, worstf, gnrng, &
417 ngs1, ngs2, nloop, npt1, icall, &
418 n_threads, ipopul_file_append, itmp_file, istmp_file, &
419 ipopul_file, ispopul_file
420 namelist /restartnml2/ &
421 history_tmp, x, xf, worstx, xnstd, &
422 bound, unit, criter, xname, iseed, save_state_unif, &
423 save_state_gauss
424
425 if (present(restart)) then
426 irestart = restart
427 else
428 irestart = .false.
429 end if
430
431 if (present(restart_file)) then
432 isrestart_file = restart_file
433 else
434 isrestart_file = 'mo_sce.restart'
435 end if
436
437 if (.not. irestart) then
438
439#ifdef MPI
440 parall = .false.
441 if (present(parallel)) then
442 if (parallel) then
443 write(*, *) 'WARNING: sce: openMP parallelization with MPI enabled is no implemented yet'
444 end if
445 end if
446#else
447 if (present(parallel)) then
448 parall = parallel
449 else
450 parall = .false.
451 end if
452#endif
453
454 if (parall) then
455 ! OpenMP or not
456 n_threads = 1
457 !$ write(output_unit,*) '--------------------------------------------------'
458 !$OMP parallel
459 !$ n_threads = OMP_GET_NUM_THREADS()
460 !$OMP end parallel
461 !$ write(output_unit,*) ' SCE is parellel with ',n_threads,' threads'
462 !$ write(output_unit,*) '--------------------------------------------------'
463 else
464 n_threads = 1
465 end if
466
467 ! One random number chain per OpenMP thread
468 allocate(rand_tmp(n_threads))
469 allocate(iseed(n_threads))
470 allocate(save_state_unif(n_threads, n_save_state))
471 allocate(save_state_gauss(n_threads, n_save_state))
472
473 if (present(mymask)) then
474 if (.not. any(mymask)) stop 'mo_sce: all parameters are masked --> none will be optimized'
475 maskpara = mymask
476 else
477 maskpara = .true.
478 end if
479
480 ! number of parameters to optimize
481 nopt = count(maskpara, dim = 1)
482 ! total number of parameters
483 nn = size(pini, 1)
484
485 ! input checking
486 if (size(prange, dim = 1) .ne. size(pini, 1)) then
487 stop 'mo_sce: prange has not matching rows'
488 end if
489 if (size(prange, dim = 2) .ne. 2) then
490 stop 'mo_sce: two colums expected for prange'
491 end if
492 bl(:) = prange(:, 1)
493 bu(:) = prange(:, 2)
494 do ii = 1, nn
495 if(((bu(ii) - bl(ii)) .lt. tiny(1.0_dp)) .and. maskpara(ii)) then
496 write(error_unit, *) 'para #', ii, ' :: range = ( ', bl(ii), ' , ', bu(ii), ' )'
497 stop 'mo_sce: inconsistent or too small parameter range'
498 end if
499 end do
500
501 !
502 ! optionals checking
503 if (present(mymaxn)) then
504 if (mymaxn .lt. 2_i4) stop 'mo_sce: maxn has to be at least 2'
505 maxn = mymaxn
506 else
507 maxn = 1000_i8
508 end if
509 if (present(mymaxit)) then
510 maxit = mymaxit
511 else
512 maxit = .false.
513 end if
514 if (present(mykstop)) then
515 if (mykstop .lt. 1_i4) stop 'mo_sce: kstop has to be at least 1'
516 kstop = mykstop
517 else
518 kstop = 10_i4
519 end if
520 if (present(mypcento)) then
521 if (mypcento .lt. 0_dp) stop 'mo_sce: pcento should be positive'
522 pcento = mypcento
523 else
524 pcento = 0.0001_dp
525 end if
526 if (present(mypeps)) then
527 if (mypeps .lt. 0_dp) stop 'mo_sce: peps should be positive'
528 peps = mypeps
529 else
530 peps = 0.001_dp
531 end if
532 if (present(myseed)) then
533 if (myseed .lt. 1_i8) stop 'mo_sce: seed should be non-negative'
534 forall(ii = 1 : n_threads) iseed(ii) = myseed + (ii - 1) * 1000_i8
535 else
536 call get_timeseed(iseed)
537 end if
538 if (present(myngs)) then
539 if (myngs .lt. 1_i4) stop 'mo_sce: ngs has to be at least 1'
540 ngs = myngs
541 else
542 ngs = 2_i4
543 end if
544 if (present(mynpg)) then
545 if (mynpg .lt. 3_i4) stop 'mo_sce: npg has to be at least 3'
546 npg = mynpg
547 else
548 npg = 2 * nopt + 1
549 end if
550 if (present(mynps)) then
551 if (mynps .lt. 2_i4) stop 'mo_sce: nps has to be at least 2'
552 nps = mynps
553 else
554 nps = nopt + 1_i4
555 end if
556 if (present(mynspl)) then
557 if (mynspl .lt. 3_i4) stop 'mo_sce: nspl has to be at least 3'
558 nspl = mynspl
559 else
560 nspl = 2 * nopt + 1
561 end if
562 if (present(mymings)) then
563 if (mymings .lt. 1_i4) stop 'mo_sce: mings has to be at least 1'
564 mings = mymings
565 else
566 mings = ngs ! no reduction of complexes
567 end if
568 if (present(myiniflg)) then
569 if ((myiniflg .ne. 1_i4) .and. (myiniflg .ne. 0_i4)) stop 'mo_sce: iniflg has to be 0 or 1'
570 iniflg = myiniflg
571 else
572 iniflg = 1_i4
573 end if
574 idot = .false.
575 if (present(myprint)) then
576 if ((myprint .lt. 0_i4) .or. (myprint .gt. 4_i4)) stop 'mo_sce: iprint has to be between 0 and 4'
577 if (myprint > 2_i4) then
578 iprint = myprint - 3
579 idot = .true.
580 else
581 iprint = myprint
582 end if
583 else
584 iprint = 2_i4 ! no printing
585 iprint = 0_i4
586 end if
587 if (present(myalpha)) then
588 alpha = myalpha
589 else
590 alpha = 0.8_dp
591 end if
592 if (present(mybeta)) then
593 beta = mybeta
594 else
595 beta = 0.45_dp
596 end if
597
598 if (present(popul_file_append)) then
599 ipopul_file_append = popul_file_append
600 else
601 ipopul_file_append = .false.
602 end if
603
604 if (present(tmp_file)) then
605 itmp_file = .true.
606 istmp_file = tmp_file
607 open(999, file = trim(adjustl(istmp_file)), action = 'write', status = 'unknown')
608 write(999, *) '# settings :: general'
609 write(999, *) '# nIterations seed'
610 write(999, *) maxn, iseed
611 write(999, *) '# settings :: sce specific'
612 write(999, *) '# sce_ngs sce_npg sce_nps'
613 write(999, *) ngs, npg, nps
614 write(999, *) '# nloop icall ngs1 bestf worstf gnrng (bestx(j),j=1,nn)'
615 close(999)
616 else
617 itmp_file = .false.
618 istmp_file = ''
619 end if
620
621 if (present(popul_file)) then
622 ipopul_file = .true.
623 ispopul_file = popul_file
624 else
625 ipopul_file = .false.
626 ispopul_file = ''
627 end if
628
629 if (ipopul_file .and. (.not. ipopul_file_append)) then
630 open(999, file = trim(adjustl(ispopul_file)), action = 'write', status = 'unknown')
631 write(999, *) '# xf(i) (x(i,j),j=1,nn)'
632 close(999)
633 end if
634
635 ! allocation of arrays
636 allocate(x(ngs * npg, nn))
637 allocate(xf(ngs * npg))
638 allocate(worstx(nn))
639 allocate(xnstd(nn))
640 allocate(bound(nn))
641 allocate(unit(nn))
642 allocate(criter(kstop + 1))
643 allocate(xname(nn))
644 allocate(history_tmp(maxn + 3 * ngs * nspl))
645 allocate(xtmp(npg, nn))
646 allocate(ftmp(npg))
647 if (maxit) then
648 large = -0.5_dp * huge(1.0_dp)
649 else
650 large = 0.5_dp * huge(1.0_dp)
651 end if
652 history_tmp(:) = large
653 criter(:) = large
654 !
655 ! initialize variables
656 do ii = 1, nn
657 xname(ii) = compress('p' // num2str(ii, '(I3.3)'))
658 end do
659
660 nloop = 0_i4
661 loop = 0_i4
662 igs = 0_i4
663 !
664 ! compute the total number of points in initial population
665 npt = ngs * npg
666 ngs1 = ngs
667 npt1 = npt
668 !
669 if (iprint .lt. 2) then
670 write(output_unit, *) '=================================================='
671 write(output_unit, *) 'ENTER THE SHUFFLED COMPLEX EVOLUTION GLOBAL SEARCH'
672 write(output_unit, *) '=================================================='
673 end if
674 !
675 ! Print seed
676 if (iprint .lt. 2) then
677 write (*, *) ' Seeds used : ', iseed
678 end if
679 ! initialize random number generator: Uniform
680 call xor4096(iseed, rand_tmp, save_state = save_state_unif)
681 ! initialize random number generator: Gaussian
682 iseed = iseed + 1000_i8
683 call xor4096g(iseed, rand_tmp, save_state = save_state_gauss)
684 iseed = 0_i8
685 !
686 ! compute the bound for parameters being optimized
687 do ii = 1, nn
688 bound(ii) = bu(ii) - bl(ii)
689 unit(ii) = 1.0_dp
690 end do
691
692 !--------------------------------------------------
693 ! compute the function value of the initial point
694 !--------------------------------------------------
695 ! function evaluation will be counted later...
696 if (idot) write(output_unit, '(A1)') '.'
697#ifdef MPI
698 call mpi_comm_size(comm, nproc, ierror)
699 do iproc = 1, nproc-1
700 do_obj_loop = .true.
701 call mpi_send(do_obj_loop,1, mpi_logical,iproc,0,comm,ierror)
702 end do
703#endif
704 if (.not. maxit) then
705 fpini = functn(pini, eval)
706 history_tmp(1) = fpini
707 else
708 fpini = -functn(pini, eval)
709 history_tmp(1) = -fpini
710 end if
711
712 ! print the initial point and its criterion value
713 bestx = pini
714 bestf_tmp = fpini
715 if (iprint .lt. 2) then
716 write(output_unit, *) ''
717 write(output_unit, *) ''
718 write(output_unit, *) '*** PRINT THE INITIAL POINT AND ITS CRITERION VALUE ***'
719 call write_best_final()
720 end if
721 if (itmp_file) then ! initial tmp file
722 open(999, file = trim(adjustl(istmp_file)), action = 'write', position = 'append', recl = (nn + 6) * 30)
723 if (.not. maxit) then
724 write(999, *) 0, 1, ngs1, fpini, fpini, 1.0_dp, pini
725 else
726 write(999, *) 0, 1, ngs1, -fpini, -fpini, 1.0_dp, pini
727 end if
728 close(999)
729 end if
730 !
731 ! generate an initial set of npt1 points in the parameter space
732 ! if iniflg is equal to 1, set x(1,.) to initial point pini(.)
733 if (iniflg .eq. 1) then
734 do ii = 1, nn
735 x(1, ii) = pini(ii)
736 end do
737 xf(1) = fpini
738 !
739 ! else, generate a point randomly and set it equal to x(1,.)
740 else
741 itmp = save_state_unif(1, :)
742 call getpnt(1, bl(1 : nn), bu(1 : nn), unit(1 : nn), pini(1 : nn), maskpara, itmp, xx)
743 save_state_unif(1, :) = itmp
744 do ii = 1, nn
745 x(1, ii) = xx(ii)
746 end do
747 if (idot) write(output_unit, '(A1)') '.'
748#ifdef MPI
749 call mpi_comm_size(comm, nproc, ierror)
750 do iproc = 1, nproc-1
751 do_obj_loop = .true.
752 call mpi_send(do_obj_loop,1, mpi_logical,iproc,0,comm,ierror)
753 end do
754#endif
755 if (.not. maxit) then
756 xf(1) = functn(xx, eval)
757 else
758 xf(1) = -functn(xx, eval)
759 end if
760 end if
761 !
762 ! count function evaluation of the first point
763 icall = 1_i8
764 ! if (icall .ge. maxn) return
765 !
766 if (parall) then
767
768 ! -----------------------------------------------------------------------
769 ! Parallel version of complex-loop
770 ! -----------------------------------------------------------------------
771 ! generate npt1-1 random points distributed uniformly in the parameter
772 ! space, and compute the corresponding function values
773 ithread = 1
774 !$OMP parallel default(shared) private(ii,jj,ithread,xx)
775 !$OMP do
776 do ii = 2, npt1
777 !$ ithread = OMP_GET_THREAD_NUM() + 1
778 itmp = save_state_unif(ithread, :)
779 call getpnt(1, bl(1 : nn), bu(1 : nn), unit(1 : nn), pini(1 : nn), maskpara, itmp, xx)
780 save_state_unif(ithread, :) = itmp
781 do jj = 1, nn
782 x(ii, jj) = xx(jj)
783 end do
784 if (idot) write(output_unit, '(A1)') '.'
785#ifdef MPI
786 call mpi_comm_size(comm, nproc, ierror)
787 do iproc = 1, nproc-1
788 do_obj_loop = .true.
789 call mpi_send(do_obj_loop,1, mpi_logical,iproc,0,comm,ierror)
790 end do
791#endif
792 if (.not. maxit) then
793 xf(ii) = functn(xx, eval)
794 history_tmp(ii) = xf(ii) ! min(history_tmp(ii-1),xf(ii)) --> will be sorted later
795 else
796 xf(ii) = -functn(xx, eval)
797 history_tmp(ii) = -xf(ii) ! max(history_tmp(ii-1),-xf(ii)) --> will be sorted later
798 end if
799 end do
800 !$OMP end do
801 !$OMP end parallel
802
803 else
804
805 ! -----------------------------------------------------------------------
806 ! Non-Parallel version of complex-loop
807 ! -----------------------------------------------------------------------
808 ! generate npt1-1 random points distributed uniformly in the parameter
809 ! space, and compute the corresponding function values
810 ithread = 1
811
812 do ii = 2, npt1
813 call getpnt(1, bl(1 : nn), bu(1 : nn), unit(1 : nn), pini(1 : nn), maskpara, save_state_unif(ithread, :), xx)
814 do jj = 1, nn
815 x(ii, jj) = xx(jj)
816 end do
817 if (idot) write(output_unit, '(A1)') '.'
818#ifdef MPI
819 call mpi_comm_size(comm, nproc, ierror)
820 do iproc = 1, nproc-1
821 do_obj_loop = .true.
822 call mpi_send(do_obj_loop,1, mpi_logical,iproc,0,comm,ierror)
823 end do
824#endif
825 if (.not. maxit) then
826 xf(ii) = functn(xx, eval)
827 icall = icall + 1_i8
828 history_tmp(icall) = xf(ii) ! min(history_tmp(icall-1),xf(ii)) --> will be sorted later
829 else
830 xf(ii) = -functn(xx, eval)
831 icall = icall + 1_i8
832 history_tmp(icall) = -xf(ii) ! max(history_tmp(icall-1),-xf(ii)) --> will be sorted later
833 end if
834
835 if (icall .ge. maxn) then
836 npt1 = ii
837 if (iprint .lt. 2) then
838 ! maximum trials reached
839 call write_termination_case(1)
840 call write_best_final()
841 end if
842 exit
843 end if
844 end do
845
846 end if
847
848 icall = int(npt1, i8)
849 !
850 ! arrange the points in order of increasing function value
851 if (maxit) then
852 large = minval(xf(1 : npt1))
853 large = merge(0.9_dp * large, 1.1_dp * large, large>0._dp)
854 else
855 large = maxval(xf(1 : npt1))
856 large = merge(1.1_dp * large, 0.9_dp * large, large>0._dp)
857 end if
858 xf(1 : npt1) = merge(xf(1 : npt1), large, is_finite(xf(1 : npt1))) ! NaN and Infinite
859 ! sort does not work with NaNs
860 ! -> get history_tmp w/o NaN, sort it, and set the rest to NaN
861 nonan = size(pack(history_tmp(1 : npt1), mask = is_finite(history_tmp(1 : npt1))))
862 if (nonan /= npt1) then
863 allocate(htmp(nonan))
864 htmp(1 : nonan) = pack(history_tmp(1 : npt1), mask = is_finite(history_tmp(1 : npt1)))
865 call sort(htmp(1 : nonan))
866 history_tmp(1 : nonan) = htmp(1 : nonan)
867 history_tmp(nonan + 1 : npt1) = special_value(1.0_dp, 'IEEE_QUIET_NAN')
868 deallocate(htmp)
869 else
870 call sort(history_tmp(1 : npt1))
871 end if
872 call sort_matrix(x(1 : npt1, 1 : nn), xf(1 : npt1))
873 !
874 ! record the best and worst points
875 do ii = 1, nn
876 bestx(ii) = x(1, ii)
877 worstx(ii) = x(npt1, ii)
878 end do
879 bestf_tmp = xf(1)
880 worstf = xf(npt1)
881 !
882 ! compute the parameter range for the initial population
883 call parstt(x(1 : npt1, 1 : nn), bound, peps, maskpara, xnstd, gnrng, ipcnvg)
884 !
885 ! write currently best x and xf to temporal file
886 if (itmp_file) then
887 call write_best_intermediate(.true.)
888 end if
889
890 ! write population to file
891 if (ipopul_file) then
892 call write_population(.true.)
893 end if
894 !
895 ! print the results for the initial population
896 print : if (iprint .lt. 2) then
897 ! write currently best x and xf to screen
898 call write_best_intermediate(.false.)
899
900 ! write population on screen
901 if (iprint .eq. 1) then
902 call write_population(.false.)
903 end if
904 end if print
905 !
906 ! Maximum number of function evaluations reached?
907 if (icall .ge. maxn) then
908 if (iprint .lt. 2) then
909 ! maximum trials reached
910 call write_termination_case(1)
911 call write_best_final()
912 end if
913 call set_optional()
914 ! -----------------------
915 ! Abort subroutine
916 ! -----------------------
917#ifdef MPI
918 call mpi_comm_size(comm, nproc, ierror)
919 do iproc = 1, nproc-1
920 do_obj_loop = .false.
921 call mpi_send(do_obj_loop,1, mpi_logical,iproc,0,comm,ierror)
922 end do
923#endif
924 return
925 end if
926 !
927 ! Feasible parameter space converged?
928 if (ipcnvg .eq. 1) then
929 if (iprint .lt. 2) then
930 ! converged because feasible parameter space small
931 call write_termination_case(3)
932 call write_best_final()
933 end if
934 call set_optional()
935 ! -----------------------
936 ! Abort subroutine
937 ! -----------------------
938#ifdef MPI
939 call mpi_comm_size(comm, nproc, ierror)
940 do iproc = 1, nproc-1
941 do_obj_loop = .false.
942 call mpi_send(do_obj_loop,1, mpi_logical,iproc,0,comm,ierror)
943 end do
944#endif
945 return
946 end if
947 !
948 ! transfer all array-like variables in namelist to fixed-size dummy-arrays
949 dummy_maskpara = .false.
950 dummy_maskpara(1 : size(pini, 1)) = maskpara
951 dummy_bl = -9999.0_dp
952 dummy_bl(1 : size(pini, 1)) = bl
953 dummy_bu = -9999.0_dp
954 dummy_bu(1 : size(pini, 1)) = bu
955 dummy_xx = -9999.0_dp
956 dummy_xx(1 : size(pini, 1)) = xx
957 !
958 ! write restart
959 open(999, file = isrestart_file, status = 'unknown', action = 'write', delim = 'quote')
960 write(999, restartnml1)
961 write(999, restartnml2)
962 close(999)
963 end if ! restart or not
964
965 if (irestart) then
966 ! read 1st namelist with allocated/scalar variables
967 open(999, file = isrestart_file, status = 'old', action = 'read', delim = 'quote')
968 read(999, nml = restartnml1)
969 close(999)
970
971 ! transfer all array-like variables in namelist to fixed-size dummy-arrays
972 maskpara = dummy_maskpara(1 : size(pini, 1))
973 bl = dummy_bl(1 : size(pini, 1))
974 bu = dummy_bu(1 : size(pini, 1))
975 xx = dummy_xx(1 : size(pini, 1))
976
977 ! allocate global arrays
978 allocate(rand_tmp(n_threads))
979 allocate(iseed(n_threads))
980 allocate(save_state_unif(n_threads, n_save_state))
981 allocate(save_state_gauss(n_threads, n_save_state))
982 allocate(x(ngs * npg, nn))
983 allocate(xf(ngs * npg))
984 allocate(worstx(nn))
985 allocate(xnstd(nn))
986 allocate(bound(nn))
987 allocate(unit(nn))
988 allocate(criter(kstop + 1))
989 allocate(xname(nn))
990 allocate(history_tmp(maxn + 3 * ngs * nspl))
991 allocate(xtmp(npg, nn))
992 allocate(ftmp(npg))
993 end if
994
995 !
996 ! begin the main loop
997 mainloop : do while (icall .lt. maxn)
998 ! read variables from restart files
999 open(999, file = isrestart_file, status = 'old', action = 'read', delim = 'quote')
1000 read(999, nml = restartnml1)
1001 read(999, nml = restartnml2)
1002 close(999)
1003 !
1004 ! transfer all array-like variables in namelist to fixed-size dummy-arrays
1005 maskpara = dummy_maskpara(1 : size(pini, 1))
1006 bl = dummy_bl(1 : size(pini, 1))
1007 bu = dummy_bu(1 : size(pini, 1))
1008 xx = dummy_xx(1 : size(pini, 1))
1009 !
1010 nloop = nloop + 1
1011 !
1012 if (iprint .lt. 2) then
1013 write(output_unit, *) ''
1014 write(output_unit, '(A28,I4)') ' *** Evolution Loop Number ', nloop
1015 end if
1016 !
1017 ! begin loop on complexes
1018 ! <beta> loop from duan(1993)
1019
1020 if (parall) then
1021
1022 ! -----------------------------------------------------------------------
1023 ! Parallel version of complex-loop
1024 ! -----------------------------------------------------------------------
1025
1026 ithread = 1
1027 !$OMP parallel default(shared) &
1028 !$OMP private(igs, loop, ithread, kk, k1, k2, jj, lpos_ok, lpos, rand, cx, cf, lcs, s, sf) &
1029 !$OMP private(icall_merk, iicall, ihistory_tmp, large)
1030 allocate(cx(npg, nn))
1031 allocate(cf(npg))
1032 allocate(lcs(nps))
1033 allocate(s(nps, nn))
1034 allocate(sf(nps))
1035 allocate(ihistory_tmp(maxn + 3 * ngs * nspl))
1036 !$OMP do
1037 comploop_omp : do igs = 1, ngs1
1038 !$ ithread = OMP_GET_THREAD_NUM() + 1
1039 !
1040 ! assign points into complexes
1041 do k1 = 1, npg
1042 k2 = (k1 - 1) * ngs1 + igs
1043 do jj = 1, nn
1044 cx(k1, jj) = x(k2, jj)
1045 end do
1046 cf(k1) = xf(k2)
1047 end do
1048 !
1049 ! begin inner loop - random selection of sub-complexes
1050 ! <alpha> loop from duan(1993)
1051 subcomploop_omp : do loop = 1, nspl
1052 !
1053 ! choose a sub-complex (nps points) according to a linear
1054 ! probability distribution
1055 !
1056 ! if number of points in subcomplex (nps) = number of points in complex (npg)
1057 ! --> select all
1058 if (nps .eq. npg) then
1059 do kk = 1, nps
1060 lcs(kk) = kk
1061 end do
1062 else
1063 call xor4096(iseed(ithread), rand, save_state = save_state_unif(ithread, :))
1064 lcs(1) = 1 + int(real(npg, dp) + 0.5_dp &
1065 - sqrt((real(npg, dp) + .5_dp)**2 - real(npg, dp) * real(npg + 1, dp) * rand))
1066 do kk = 2, nps
1067 lpos_ok = .false.
1068 do while (.not. lpos_ok)
1069 lpos_ok = .true.
1070 call xor4096(iseed(ithread), rand, save_state = save_state_unif(ithread, :))
1071 lpos = 1 + int(real(npg, dp) + 0.5_dp &
1072 - sqrt((real(npg, dp) + .5_dp)**2 - real(npg, dp) * real(npg + 1, dp) * rand))
1073 ! test if point already chosen: returns lpos_ok=false if any point already exists
1074 do k1 = 1, kk - 1
1075 if (lpos .eq. lcs(k1)) lpos_ok = (lpos_ok .and. .false.)
1076 end do
1077 end do
1078 lcs(kk) = lpos
1079 end do
1080 !
1081 ! arrange the sub-complex in order of increasing function value
1082 call sort(lcs(1 : nps))
1083 end if
1084 !
1085 ! create the sub-complex arrays
1086 do kk = 1, nps
1087 do jj = 1, nn
1088 s(kk, jj) = cx(lcs(kk), jj)
1089 end do
1090 sf(kk) = cf(lcs(kk))
1091 end do
1092 !
1093 ! remember largest for treating of NaNs
1094 if (maxit) then
1095 large = minval(cf(1 : npg))
1096 large = merge(0.9_dp * large, 1.1_dp * large, large>0._dp)
1097 else
1098 large = maxval(cf(1 : npg))
1099 large = merge(1.1_dp * large, 0.9_dp * large, large>0._dp)
1100 end if
1101 !
1102 ! use the sub-complex to generate new point(s)
1103 icall_merk = icall
1104 iicall = icall
1105 ihistory_tmp = history_tmp
1106#ifdef MPI
1107 call cce(s(1 : nps, 1 : nn), sf(1 : nps), bl(1 : nn), bu(1 : nn), maskpara, xnstd(1 : nn), &
1108 iicall, maxn, maxit, save_state_gauss(ithread, :), functn, eval, alpha, beta, ihistory_tmp, idot, comm)
1109#else
1110 call cce(s(1 : nps, 1 : nn), sf(1 : nps), bl(1 : nn), bu(1 : nn), maskpara, xnstd(1 : nn), &
1111 iicall, maxn, maxit, save_state_gauss(ithread, :), functn, eval, alpha, beta, ihistory_tmp, idot)
1112#endif
1113 history_tmp(icall + 1 : icall + (iicall - icall_merk)) = ihistory_tmp(icall_merk + 1 : iicall)
1114 icall = icall + (iicall - icall_merk)
1115 !
1116 ! if the sub-complex is accepted, replace the new sub-complex
1117 ! into the complex
1118 do kk = 1, nps
1119 do jj = 1, nn
1120 cx(lcs(kk), jj) = s(kk, jj)
1121 end do
1122 cf(lcs(kk)) = sf(kk)
1123 end do
1124 !
1125 ! sort the points
1126 cf(1 : npg) = merge(cf(1 : npg), large, is_finite(cf(1 : npg))) ! NaN and Infinite
1127 call sort_matrix(cx(1 : npg, 1 : nn), cf(1 : npg))
1128 !
1129 ! ! if maximum number of runs exceeded, break out of the loop
1130 ! if (icall .ge. maxn) exit
1131 !
1132 end do subcomploop_omp ! <alpha loop>
1133 !
1134 ! replace the new complex into original array x(.,.)
1135 do k1 = 1, npg
1136 k2 = (k1 - 1) * ngs1 + igs
1137 do jj = 1, nn
1138 x(k2, jj) = cx(k1, jj)
1139 end do
1140 xf(k2) = cf(k1)
1141 end do
1142 ! if (icall .ge. maxn) exit
1143 !
1144 ! end loop on complexes
1145 end do comploop_omp ! <beta loop>
1146 !$OMP end do
1147 deallocate(cx)
1148 deallocate(cf)
1149 deallocate(lcs)
1150 deallocate(s)
1151 deallocate(sf)
1152 deallocate(ihistory_tmp)
1153 !$OMP end parallel
1154
1155 else
1156
1157 ! -----------------------------------------------------------------------
1158 ! Non-parallel version of complex-loop
1159 ! -----------------------------------------------------------------------
1160
1161 allocate(cx(npg, nn))
1162 allocate(cf(npg))
1163 allocate(lcs(nps))
1164 allocate(s(nps, nn))
1165 allocate(sf(nps))
1166
1167 ithread = 1
1168
1169 comploop : do igs = 1, ngs1
1170 !
1171 ! assign points into complexes
1172 do k1 = 1, npg
1173 k2 = (k1 - 1) * ngs1 + igs
1174 do jj = 1, nn
1175 cx(k1, jj) = x(k2, jj)
1176 end do
1177 cf(k1) = xf(k2)
1178 end do
1179 !
1180 ! begin inner loop - random selection of sub-complexes
1181 ! <alpha> loop from duan(1993)
1182 subcomploop : do loop = 1, nspl
1183 !
1184 ! choose a sub-complex (nps points) according to a linear
1185 ! probability distribution
1186 !
1187 ! if number of points in subcomplex (nps) = number of points in complex (npg)
1188 ! --> select all
1189 if (nps .eq. npg) then
1190 do kk = 1, nps
1191 lcs(kk) = kk
1192 end do
1193 else
1194 call xor4096(iseed(ithread), rand, save_state = save_state_unif(ithread, :))
1195 lcs(1) = 1 + int(real(npg, dp) + 0.5_dp &
1196 - sqrt((real(npg, dp) + .5_dp)**2 - real(npg, dp) * real(npg + 1, dp) * rand))
1197 do kk = 2, nps
1198 lpos_ok = .false.
1199 do while (.not. lpos_ok)
1200 lpos_ok = .true.
1201 call xor4096(iseed(ithread), rand, save_state = save_state_unif(ithread, :))
1202 lpos = 1 + int(real(npg, dp) + 0.5_dp &
1203 - sqrt((real(npg, dp) + .5_dp)**2 - real(npg, dp) * real(npg + 1, dp) * rand))
1204 ! test if point already chosen: returns lpos_ok=false if any point already exists
1205 do k1 = 1, kk - 1
1206 if (lpos .eq. lcs(k1)) lpos_ok = (lpos_ok .and. .false.)
1207 end do
1208 end do
1209 lcs(kk) = lpos
1210 end do
1211 !
1212 ! arrange the sub-complex in order of increasing function value
1213 call sort(lcs(1 : nps))
1214 end if
1215 !
1216 ! create the sub-complex arrays
1217 do kk = 1, nps
1218 do jj = 1, nn
1219 s(kk, jj) = cx(lcs(kk), jj)
1220 end do
1221 sf(kk) = cf(lcs(kk))
1222 end do
1223 !
1224 ! remember largest for treating of NaNs
1225 if (maxit) then
1226 large = minval(cf(1 : npg))
1227 large = merge(0.9_dp * large, 1.1_dp * large, large>0._dp)
1228 else
1229 large = maxval(cf(1 : npg))
1230 large = merge(1.1_dp * large, 0.9_dp * large, large>0._dp)
1231 end if
1232 !
1233 ! use the sub-complex to generate new point(s)
1234#ifdef MPI
1235 call cce(s(1 : nps, 1 : nn), sf(1 : nps), bl(1 : nn), bu(1 : nn), maskpara, xnstd(1 : nn), &
1236 icall, maxn, maxit, save_state_gauss(ithread, :), functn, eval, alpha, beta, history_tmp, idot, comm)
1237#else
1238 call cce(s(1 : nps, 1 : nn), sf(1 : nps), bl(1 : nn), bu(1 : nn), maskpara, xnstd(1 : nn), &
1239 icall, maxn, maxit, save_state_gauss(ithread, :), functn, eval, alpha, beta, history_tmp, idot)
1240#endif
1241 !
1242 ! if the sub-complex is accepted, replace the new sub-complex
1243 ! into the complex
1244 do kk = 1, nps
1245 do jj = 1, nn
1246 cx(lcs(kk), jj) = s(kk, jj)
1247 end do
1248 cf(lcs(kk)) = sf(kk)
1249 end do
1250 !
1251 ! sort the points
1252 cf(1 : npg) = merge(cf(1 : npg), large, is_finite(cf(1 : npg))) ! NaN and Infinite
1253 call sort_matrix(cx(1 : npg, 1 : nn), cf(1 : npg))
1254 !
1255 ! if maximum number of runs exceeded, break out of the loop
1256 if (icall .ge. maxn) then
1257 if (iprint .lt. 2) then
1258 ! maximum trials reached
1259 call write_termination_case(1)
1260 call write_best_final()
1261 end if
1262 exit
1263 end if
1264 !
1265 end do subcomploop ! <alpha loop>
1266 !
1267 ! replace the new complex into original array x(.,.)
1268 do k1 = 1, npg
1269 k2 = (k1 - 1) * ngs1 + igs
1270 do jj = 1, nn
1271 x(k2, jj) = cx(k1, jj)
1272 end do
1273 xf(k2) = cf(k1)
1274 end do
1275 if (icall .ge. maxn) then
1276 if (iprint .lt. 2) then
1277 ! maximum trials reached
1278 call write_termination_case(1)
1279 call write_best_final()
1280 end if
1281 exit
1282 end if
1283 !
1284 ! end loop on complexes
1285 end do comploop ! <beta loop>
1286
1287 deallocate(cx)
1288 deallocate(cf)
1289 deallocate(lcs)
1290 deallocate(s)
1291 deallocate(sf)
1292
1293 end if ! end parallel/non-parallel region
1294 !
1295 ! re-sort the points
1296 call sort_matrix(x(1 : npt1, 1 : nn), xf(1 : npt1))
1297 !
1298 ! record the best and worst points
1299 do jj = 1, nn
1300 bestx(jj) = x(1, jj)
1301 worstx(jj) = x(npt1, jj)
1302 end do
1303 bestf_tmp = xf(1)
1304 worstf = xf(npt1)
1305 !
1306 ! test the population for parameter convergence
1307 call parstt(x(1 : npt1, 1 : nn), bound, peps, maskpara, xnstd, gnrng, ipcnvg)
1308 !
1309 ! write currently best x and xf to temporal file
1310 if (itmp_file) then
1311 call write_best_intermediate(.true.)
1312 end if
1313 !
1314 ! write population to file
1315 if (ipopul_file) then
1316 call write_population(.true.)
1317 end if
1318 !
1319 ! print the results for current population
1320 if (iprint .lt. 2) then
1321 call write_best_intermediate(.false.)
1322
1323 if (iprint .eq. 1) then
1324 call write_population(.false.)
1325 end if
1326
1327 end if
1328 !
1329 ! test if maximum number of function evaluations exceeded
1330 ! Maximum number of function evaluations reached?
1331 if (icall .ge. maxn) then
1332 if (iprint .lt. 2) then
1333 ! maximum trials reached
1334 call write_termination_case(1)
1335 call write_best_final()
1336 end if
1337 call set_optional()
1338 ! -----------------------
1339 ! Abort subroutine
1340 ! -----------------------
1341#ifdef MPI
1342 call mpi_comm_size(comm, nproc, ierror)
1343 do iproc = 1, nproc-1
1344 do_obj_loop = .false.
1345 call mpi_send(do_obj_loop,1, mpi_logical,iproc,0,comm,ierror)
1346 end do
1347#endif
1348 return
1349 end if
1350 !
1351 ! compute the count on successive loops w/o function improvement
1352 criter(kstop + 1) = bestf_tmp
1353 if (nloop .gt. kstop) then
1354 denomi = 0.5_dp * abs(criter((kstop + 1) - kstop) + criter(kstop + 1))
1355 denomi = max(denomi, 1.0e-15_dp)
1356 timeou = abs(criter((kstop + 1) - kstop) - criter(kstop + 1)) / denomi
1357 if (timeou .lt. pcento) then
1358 if (iprint .lt. 2) then
1359 ! criterion value has not changed during last loops
1360 call write_termination_case(2)
1361 call write_best_final()
1362 end if
1363 call set_optional()
1364 ! -----------------------
1365 ! Abort subroutine
1366 ! -----------------------
1367#ifdef MPI
1368 call mpi_comm_size(comm, nproc, ierror)
1369 do iproc = 1, nproc-1
1370 do_obj_loop = .false.
1371 call mpi_send(do_obj_loop,1, mpi_logical,iproc,0,comm,ierror)
1372 end do
1373#endif
1374 return
1375 end if
1376 end if
1377 criter(1 : kstop) = criter(2 : kstop + 1)
1378 !
1379 ! if population is converged into a sufficiently small space
1380 if (ipcnvg .eq. 1) then
1381 if (iprint .lt. 2) then
1382 ! converged because feasible parameter space small
1383 call write_termination_case(3)
1384 call write_best_final()
1385 end if
1386 call set_optional()
1387 ! -----------------------
1388 ! Abort subroutine
1389 ! -----------------------
1390#ifdef MPI
1391 call mpi_comm_size(comm, nproc, ierror)
1392 do iproc = 1, nproc-1
1393 do_obj_loop = .false.
1394 call mpi_send(do_obj_loop,1, mpi_logical,iproc,0,comm,ierror)
1395 end do
1396#endif
1397 return
1398 end if
1399 !
1400 ! none of the stopping criteria is satisfied, continue search
1401 !
1402 ! check for complex number reduction
1403 if (ngs1 .gt. mings) then
1404 ngs2 = ngs1
1405 ngs1 = ngs1 - 1
1406 npt1 = ngs1 * npg
1407 call comp(ngs2, npg, x(1 : ngs2 * npg, 1 : nn), xf(1 : ngs2 * npg), xtmp(1 : ngs2 * npg, 1 : nn), ftmp(1 : ngs2 * npg))
1408 end if
1409 !
1410 ! transfer all array-like variables in namelist to fixed-size dummy-arrays
1411 dummy_maskpara = .false.
1412 dummy_maskpara(1 : size(pini, 1)) = maskpara
1413 dummy_bl = -9999.0_dp
1414 dummy_bl(1 : size(pini, 1)) = bl
1415 dummy_bu = -9999.0_dp
1416 dummy_bu(1 : size(pini, 1)) = bu
1417 dummy_xx = -9999.0_dp
1418 dummy_xx(1 : size(pini, 1)) = xx
1419 !
1420 ! write restart
1421 open(999, file = isrestart_file, status = 'unknown', action = 'write', delim = 'quote')
1422 write(999, restartnml1)
1423 write(999, restartnml2)
1424 close(999)
1425 end do mainloop
1426
1427 deallocate(xtmp)
1428 deallocate(ftmp)
1429
1430 call set_optional()
1431
1432#ifdef MPI
1433 call mpi_comm_size(comm, nproc, ierror)
1434 do iproc = 1, nproc-1
1435 do_obj_loop = .false.
1436 call mpi_send(do_obj_loop,1, mpi_logical,iproc,0,comm,ierror)
1437 end do
1438#endif
1439 contains
1440
1441 subroutine write_best_intermediate(to_file)
1442
1443 implicit none
1444
1445 logical, intent(in) :: to_file
1446 integer(i4) :: ll
1447
1448 if (to_file) then
1449 open(999, file = trim(adjustl(istmp_file)), action = 'write', position = 'append', recl = (nn + 6) * 30)
1450 if (.not. maxit) then
1451 write(999, *) nloop, icall, ngs1, bestf_tmp, worstf, gnrng, (bestx(ll), ll = 1, nn)
1452 else
1453 write(999, *) nloop, icall, ngs1, -bestf_tmp, -worstf, gnrng, (bestx(ll), ll = 1, nn)
1454 end if
1455 close(999)
1456 else
1457 write(format_str1, '(A13,I3,A8)') '( A49, ', nn, '(6X,A4))'
1458 write(format_str2, '(A26,I3,A8)') '(I5,1X,I5,3X,I5,3(E22.14), ', nn, '(G10.3))'
1459 if (nloop == 0) then
1460 write(output_unit, *) ''
1461 write(output_unit, '(A44)') ' *** PRINT THE RESULTS OF THE SCE SEARCH ***'
1462 write(output_unit, *) ''
1463 write(output_unit, format_str1) ' LOOP TRIALS COMPLXS BEST F WORST F PAR RNG ', (xname(ll), ll = 1, nn)
1464 end if
1465 if (.not. maxit) then
1466 write(output_unit, format_str2) nloop, icall, ngs1, bestf_tmp, worstf, gnrng, (bestx(ll), ll = 1, nn)
1467 else
1468 write(output_unit, format_str2) nloop, icall, ngs1, -bestf_tmp, -worstf, gnrng, (bestx(ll), ll = 1, nn)
1469 end if
1470 end if
1471
1472 end subroutine write_best_intermediate
1473
1474 subroutine write_best_final()
1475
1476 implicit none
1477
1478 integer(i4) :: ll
1479
1480 write(format_str1, '(A13,I3,A8)') '( A10, ', nn, '(6X,A4))'
1481 write(format_str2, '(A14,I3,A8)') '(E22.14, ', nn, '(G10.3))'
1482
1483 write(output_unit, format_str1) 'CRITERION ', (xname(ll), ll = 1, nn)
1484 if (.not. maxit) then
1485 write(output_unit, format_str2) bestf_tmp, (bestx(ll), ll = 1, nn)
1486 else
1487 write(output_unit, format_str2) -bestf_tmp, (bestx(ll), ll = 1, nn)
1488 end if
1489
1490 end subroutine write_best_final
1491
1492 subroutine write_population(to_file)
1493
1494 logical, intent(in) :: to_file
1495
1496 integer(i4) :: ll, mm
1497
1498 if (to_file) then
1499 write(format_str2, '(A13,I3,A9)') '(I4, E22.14, ', nn, '(E22.14))'
1500 open(unit = 999, file = trim(adjustl(ispopul_file)), action = 'write', position = 'append', recl = (nn + 2) * 30)
1501 if (.not. maxit) then
1502 do mm = 1, npt1
1503 write(999, *) nloop, xf(mm), (x(mm, ll), ll = 1, nn)
1504 end do
1505 else
1506 do mm = 1, npt1
1507 write(999, *) nloop, -xf(mm), (x(mm, ll), ll = 1, nn)
1508 end do
1509 end if
1510 close(999)
1511 else
1512 write(format_str2, '(A12,I3,A8)') '(I4, E22.14, ', nn, '(G10.3))'
1513 write(output_unit, *) ''
1514 write(output_unit, '(A22,I3)') ' POPULATION AT LOOP ', nloop
1515 write(output_unit, '(A27)') '---------------------------'
1516 if (.not. maxit) then
1517 do mm = 1, npt1
1518 write(output_unit, format_str2) nloop, xf(mm), (x(mm, ll), ll = 1, nn)
1519 end do
1520 else
1521 do mm = 1, npt1
1522 write(output_unit, format_str2) nloop, -xf(mm), (x(mm, ll), ll = 1, nn)
1523 end do
1524 end if
1525 end if
1526
1527 end subroutine write_population
1528
1529 subroutine write_termination_case(case)
1530
1531 integer(i4), intent(in) :: case
1532
1533 select case (case)
1534 case (1) ! maximal number of iterations reached
1535 write(output_unit, *) ''
1536 write(output_unit, '(A46,A39,I7,A46,I4,A12,I4,A19,I4,A4)') &
1537 '*** OPTIMIZATION SEARCH TERMINATED BECAUSE THE', &
1538 ' LIMIT ON THE MAXIMUM NUMBER OF TRIALS ', maxn, &
1539 ' EXCEEDED. SEARCH WAS STOPPED AT SUB-COMPLEX ', &
1540 loop, ' OF COMPLEX ', igs, ' IN SHUFFLING LOOP ', nloop, ' ***'
1541 !
1542 case(2) ! objective not changed during last evolution loops
1543 write(output_unit, *) ''
1544 write(output_unit, '(A72,F8.4,A12,I3,A20)') '*** OPTIMIZATION TERMINATED BECAUSE THE CRITERION VALUE HAS NOT CHANGED ', &
1545 pcento * 100._dp, ' PERCENT IN ', kstop, ' SHUFFLING LOOPS ***'
1546 !
1547 case(3) ! complexes converged
1548 write(output_unit, *) ''
1549 write(output_unit, '(A50,A20,F8.4,A34)') '*** OPTIMIZATION TERMINATED BECAUSE THE POPULATION', &
1550 ' HAS CONVERGED INTO ', gnrng * 100._dp, ' PERCENT OF THE FEASIBLE SPACE ***'
1551 !
1552 case default
1553 write(error_unit, *) 'This termination case is not implemented!'
1554 stop
1555 end select
1556
1557 write(output_unit, *) ''
1558 write(output_unit, '(A66)') '*** PRINT THE FINAL PARAMETER ESTIMATE AND ITS CRITERION VALUE ***'
1559
1560 end subroutine write_termination_case
1561
1562 subroutine set_optional()
1563
1564 implicit none
1565
1566 if (present(neval)) neval = icall
1567 if (present(bestf) .and. .not. maxit) bestf = bestf_tmp
1568 if (present(bestf) .and. maxit) bestf = -bestf_tmp
1569 if (present(history)) then
1570 allocate(history(icall))
1571 history(:) = history_tmp(1 : icall)
1572 end if
1573
1574 end subroutine set_optional
1575
1576 end function sce
1577
1578
1579 subroutine parstt(x, bound, peps, mask, xnstd, gnrng, ipcnvg)
1580 !
1581 ! subroutine checking for parameter convergence
1582 !
1583 use mo_kind, only : i4, dp
1584
1585 implicit none
1586
1587 real(dp), dimension(:, :), intent(in) :: x ! points in population, cols=nn, rows=npt1 (!)
1588 real(dp), dimension(:), intent(in) :: bound ! difference of upper and lower limit per parameter
1589 real(dp), intent(in) :: peps ! optimization is terminated if volume of complex has
1590 ! ! converged to given percentage of feasible space
1591 logical, dimension(:), intent(in) :: mask ! mask of parameters
1592 real(dp), dimension(size(bound, 1)), intent(out) :: xnstd ! std. deviation of points in population per parameter
1593 real(dp), intent(out) :: gnrng ! fraction of feasible space covered by complexes
1594 integer(i4), intent(out) :: ipcnvg ! 1 : population converged into sufficiently small space
1595 ! ! 0 : not converged
1596 !
1597 ! local variables
1598 integer(i4) :: nn ! number of parameters
1599 integer(i4) :: npt1 ! number of points in current population
1600 integer(i4) :: ii, kk
1601 real(dp) :: xsum1 ! sum of all values per parameter
1602 real(dp) :: xsum2 ! sum of square of all values per parameter
1603 real(dp) :: gsum ! sum of all (range scaled) currently covered parameter ranges:
1604 ! ! (max-min)/bound
1605 real(dp), dimension(size(bound, 1)) :: xmin ! minimal value per parameter
1606 real(dp), dimension(size(bound, 1)) :: xmax ! maximal value per parameter
1607 real(dp), dimension(size(bound, 1)) :: xmean ! mean value per parameter
1608 real(dp), parameter :: delta = tiny(1.0_dp)
1609 !
1610 nn = size(x, 2)
1611 npt1 = size(x, 1)
1612
1613 ! compute maximum, minimum and standard deviation of parameter values
1614 gsum = 0._dp
1615 do kk = 1, nn
1616 if (mask(kk)) then
1617 xmax(kk) = -huge(1.0_dp)
1618 xmin(kk) = huge(1.0_dp)
1619 xsum1 = 0._dp
1620 xsum2 = 0._dp
1621 do ii = 1, npt1
1622 xmax(kk) = dmax1(x(ii, kk), xmax(kk))
1623 xmin(kk) = dmin1(x(ii, kk), xmin(kk))
1624 xsum1 = xsum1 + x(ii, kk)
1625 xsum2 = xsum2 + x(ii, kk) * x(ii, kk)
1626 end do
1627 xmean(kk) = xsum1 / real(npt1, dp)
1628 xnstd(kk) = (xsum2 / real(npt1, dp) - xmean(kk) * xmean(kk))
1629 if (xnstd(kk) .le. delta) then
1630 xnstd(kk) = delta
1631 end if
1632 xnstd(kk) = sqrt(xnstd(kk))
1633 xnstd(kk) = xnstd(kk) / bound(kk)
1634 gsum = gsum + log(delta + (xmax(kk) - xmin(kk)) / bound(kk))
1635 end if
1636 end do
1637 gnrng = exp(gsum / real(nn, dp))
1638 !
1639 ! check if normalized standard deviation of parameter is <= eps
1640 ipcnvg = 0_i4
1641 if (gnrng .le. peps) then
1642 ipcnvg = 1_i4
1643 end if
1644
1645 end subroutine parstt
1646
1647 subroutine comp(ngs2, npg, a, af, b, bf)
1648 !
1649 ! This subroutine reduce
1650 ! input matrix a(n,ngs2*npg) to matrix b(n,ngs1*npg) and
1651 ! input vector af(ngs2*npg) to vector bf(ngs1*npg)
1652 !
1653 ! Needed for complex reduction:
1654 ! Number of complexes decreases from ngs2 to ngs1
1655 ! Therefore, number of points in population decreases (npt+npg) to npt
1656 !
1657 use mo_kind, only : i4, dp
1658
1659 implicit none
1660
1661 integer(i4), intent(in) :: ngs2 ! OLD number of complexes = ngs1+1
1662 integer(i4), intent(in) :: npg ! number of points in each complex
1663 real(dp), dimension(:, :), intent(inout) :: a ! OLD points, ncols=n, nrows=ngs2*npg
1664 real(dp), dimension(:), intent(inout) :: af ! OLD function values, nrows=ngs2*npg
1665 real(dp), dimension(size(a, 1), size(a, 2)), intent(out) :: b ! NEW points, ncols=n, nrows=ngs1*npg=npt
1666 real(dp), dimension(size(af)), intent(out) :: bf ! NEW function values, nrows=ngs1*npg=npt
1667 !
1668 ! local variables
1669 integer(i4) :: n ! cols: number of parameters: nopt
1670 integer(i4) :: npt ! rows: number of points in NEW population: npt1
1671 integer(i4) :: ngs1 ! NEW number of complexes
1672 integer(i4) :: i, j
1673 integer(i4) :: igs ! current new complex
1674 integer(i4) :: ipg ! current new point in complex igs
1675 integer(i4) :: k1, k2
1676
1677 n = size(a, dim = 2)
1678 ! NEW number of complexes
1679 ngs1 = ngs2 - 1
1680 ! number of points in NEW population
1681 npt = ngs1 * npg
1682
1683 do igs = 1, ngs1
1684 do ipg = 1, npg
1685 k1 = (ipg - 1) * ngs2 + igs
1686 k2 = (ipg - 1) * ngs1 + igs
1687 do i = 1, n
1688 b(k2, i) = a(k1, i)
1689 end do
1690 bf(k2) = af(k1)
1691 end do
1692 end do
1693 !
1694 do j = 1, npt
1695 do i = 1, n
1696 a(j, i) = b(j, i)
1697 end do
1698 af(j) = bf(j)
1699 end do
1700 !
1701 end subroutine comp
1702
1703 subroutine sort_matrix(rb, ra)
1704 !
1705 use mo_kind, only : i4, dp
1706 use mo_orderpack, only : sort_index
1707
1708 implicit none
1709
1710 real(dp), dimension(:), intent(inout) :: ra ! rows: number of points in population --> npt1
1711 real(dp), dimension(:, :), intent(inout) :: rb ! rows: number of points in population --> npt1
1712 ! ! cols: number of parameters --> nn
1713 ! local variables
1714 integer(i4) :: n ! number of points in population --> npt1
1715 integer(i4) :: m ! number of parameters --> nn
1716 integer(i4) :: i
1717 integer(i4), dimension(size(rb, 1)) :: iwk
1718 !
1719 n = size(rb, 1)
1720 m = size(rb, 2)
1721
1722 ! indexes of sorted reference vector
1723 iwk(:) = sort_index(ra(1 : n))
1724
1725 ! sort reference vector
1726 ra(1 : n) = ra(iwk)
1727
1728 ! sort each column of second array
1729 do i = 1, m
1730 rb(1 : n, i) = rb(iwk, i)
1731 end do
1732
1733 end subroutine sort_matrix
1734
1735 subroutine chkcst(x, bl, bu, mask, ibound)
1736 !
1737 ! This subroutine check if the trial point satisfies all
1738 ! constraints.
1739 !
1740 ! ibound - violation indicator
1741 ! = 0 no violation
1742 ! = 1 violation
1743 ! nn = number of optimizing variables
1744 ! ii = the ii'th variable of the arrays x, bl, and bu
1745 !
1746 ! Note: removed checking for implicit constraints (was anyway empty)
1747 !
1748 use mo_kind, only : i4, dp
1749 implicit none
1750
1751 real(dp), dimension(:), intent(in) :: x ! trial point dim=(nn)
1752 real(dp), dimension(:), intent(in) :: bl ! lower bound dim=(nn)
1753 real(dp), dimension(:), intent(in) :: bu ! upper bound dim=(nn)
1754 logical, dimension(:), intent(in) :: mask ! parameter mask
1755 integer(i4), intent(out) :: ibound ! violation indicator
1756
1757 ! local variables
1758 integer(i4) :: ii
1759 integer(i4) :: nn
1760 !
1761 ibound = 0_i4
1762 nn = size(x, 1)
1763
1764 do ii = 1, nn
1765 if (mask(ii)) then
1766 if ((x(ii) .lt. bl(ii)) .or. (x(ii) .gt. bu(ii))) then
1767 ! This constraint is violated
1768 ibound = 1_i4
1769 return
1770 end if
1771 end if
1772 end do
1773
1774 end subroutine chkcst
1775
1776 subroutine getpnt(idist, bl, bu, std, xi, mask, save_state, x)
1777 !
1778 ! This subroutine generates a new point within feasible region
1779 !
1780 ! Note: checking of implicit constraints removed
1781 !
1782 use mo_kind, only : i4, i8, dp
1783 use mo_xor4096, only : n_save_state, xor4096, xor4096g
1784
1785 implicit none
1786
1787 integer(i4), intent(in) :: idist ! idist = probability flag
1788 ! ! = 1 - uniform distribution
1789 ! ! = 2 - Gaussian distribution
1790 real(dp), dimension(:), intent(in) :: bl ! lower bound
1791 real(dp), dimension(:), intent(in) :: bu ! upper bound
1792 real(dp), dimension(:), intent(in) :: std ! standard deviation of probability distribution
1793 real(dp), dimension(:), intent(in) :: xi ! focal point
1794 logical, dimension(:), intent(in) :: mask ! mask of parameters
1795 integer(i8), dimension(n_save_state), intent(inout) :: save_state ! save state of random number stream
1796 ! ! --> stream has to be according to idist
1797 real(dp), dimension(size(xi, 1)), intent(out) :: x ! new point
1798 !
1799 ! local variables
1800 integer(i4) :: nn ! number of parameters
1801 integer(i4) :: jj
1802 integer(i4) :: ibound ! 0=point in bound, 1=point out of bound
1803 real(dp) :: rand
1804 !
1805 nn = size(xi, 1)
1806 do jj = 1, nn
1807 if (mask(jj)) then
1808 ibound = 1
1809 do while (ibound .eq. 1)
1810 if (idist .eq. 1) then
1811 call xor4096(0_i8, rand, save_state = save_state)
1812 x(jj) = bl(jj) + std(jj) * rand * (bu(jj) - bl(jj))
1813 end if
1814 if (idist .eq. 2) then
1815 call xor4096g(0_i8, rand, save_state = save_state)
1816 x(jj) = xi(jj) + std(jj) * rand * (bu(jj) - bl(jj))
1817 end if
1818 !
1819 ! Check explicit constraints
1820 !
1821 call chkcst((/x(jj)/), (/bl(jj)/), (/bu(jj)/), (/mask(jj)/), ibound)
1822 end do
1823 else
1824 x(jj) = xi(jj)
1825 end if
1826 end do
1827
1828 end subroutine getpnt
1829
1830 subroutine cce(s, sf, bl, bu, maskpara, xnstd, icall, maxn, maxit, save_state_gauss, functn, eval, &
1831 alpha, beta, history, idot &
1832#ifdef MPI
1833 , comm &
1834#endif
1835 )
1836 !
1837 ! algorithm generate a new point(s) from a sub-complex
1838 !
1839 ! Note: new intent IN variables for flexible reflection & contraction: alpha, beta
1840 ! new intent INOUT variable for history of objective function values
1841 !
1842 use mo_kind, only : i4, i8, dp
1843 use mo_xor4096, only : n_save_state
1844 use iso_fortran_env, only : output_unit
1845 use mo_utils, only : is_finite
1846 use mo_optimization_utils, only : eval_interface, objective_interface
1847#ifdef MPI
1848 use mpi_f08
1849#endif
1850
1851 implicit none
1852
1853 real(dp), dimension(:, :), intent(inout) :: s ! points in sub-complex, cols=nn, rows=nps
1854 real(dp), dimension(:), intent(inout) :: sf ! objective function value of points in
1855 ! ! sub-complex, rows=nps
1856 real(dp), dimension(:), intent(in) :: bl ! lower bound per parameter
1857 real(dp), dimension(:), intent(in) :: bu ! upper bound per parameter
1858 logical, dimension(:), intent(in) :: maskpara ! mask of parameters
1859 real(dp), dimension(:), intent(in) :: xnstd ! standard deviation of points in
1860 ! ! sub-complex per parameter
1861 integer(i8), intent(inout) :: icall ! number of function evaluations
1862 integer(i8), intent(in) :: maxn ! maximal number of function evaluations allowed
1863 logical, intent(in) :: maxit ! minimization (true) or maximization (false)
1864 integer(i8), dimension(n_save_state), intent(inout) :: save_state_gauss ! seed for random number generation
1865 procedure(objective_interface), intent(in), pointer :: functn
1866 procedure(eval_interface), INTENT(IN), POINTER :: eval
1867 real(dp), intent(in) :: alpha ! parameter for reflection steps
1868 real(dp), intent(in) :: beta ! parameter for contraction steps
1869 real(dp), dimension(maxn), intent(inout) :: history ! history of best function value
1870 logical, intent(in) :: idot ! if true: progress report with '.'
1871#ifdef MPI
1872 type(mpi_comm), optional, intent(in) :: comm ! MPI communicator
1873#endif
1874 !
1875 ! local variables
1876 integer(i4) :: nn ! number of parameters
1877 integer(i4) :: nps ! number of points in a sub-complex
1878 integer(i4) :: i, j
1879 integer(i4) :: n ! nps = number of points in sub-complex
1880 integer(i4) :: m ! nn = number of parameters
1881 integer(i4) :: ibound ! ibound = flag indicating if constraints are violated
1882 ! ! = 1 yes
1883 ! ! = 0 no
1884 real(dp), dimension(size(s, 2)) :: sw ! the worst point of the simplex
1885 real(dp), dimension(size(s, 2)) :: sb ! the best point of the simplex
1886 real(dp), dimension(size(s, 2)) :: ce ! the centroid of the simplex excluding wo
1887 real(dp), dimension(size(s, 2)) :: snew ! new point generated from the simplex
1888 real(dp) :: fw ! function value of the worst point
1889 real(dp) :: fnew ! function value of the new point
1890#ifdef MPI
1891 integer(i4) :: ierror
1892 logical :: do_obj_loop
1893 integer(i4) :: iproc, nproc
1894#endif
1895
1896 nps = size(s, 1)
1897 nn = size(s, 2)
1898 !
1899 ! equivalence of variables for readabilty of code
1900 n = nps
1901 m = nn
1902 !
1903 ! identify the worst point wo of the sub-complex s
1904 ! compute the centroid ce of the remaining points
1905 ! compute step, the vector between wo and ce
1906 ! identify the worst function value fw
1907 do j = 1, m
1908 sb(j) = s(1, j)
1909 sw(j) = s(n, j)
1910 ce(j) = 0.0_dp
1911 do i = 1, n - 1
1912 ce(j) = ce(j) + s(i, j)
1913 end do
1914 ce(j) = ce(j) / real(n - 1, dp)
1915 end do
1916 fw = sf(n)
1917 !
1918 ! compute the new point snew
1919 !
1920 ! first try a reflection step
1921 do j = 1, m
1922 if (maskpara(j)) then
1923 snew(j) = ce(j) + alpha * (ce(j) - sw(j))
1924 else
1925 snew(j) = s(1, j)
1926 end if
1927 end do
1928 !
1929 ! check if snew satisfies all constraints
1930 call chkcst(snew(1 : nn), bl(1 : nn), bu(1 : nn), maskpara(1 : nn), ibound)
1931 !
1932 ! snew is outside the bound,
1933 ! choose a point at random within feasible region according to
1934 ! a normal distribution with best point of the sub-complex
1935 ! as mean and standard deviation of the population as std
1936 if (ibound .eq. 1) then
1937 call getpnt(2, bl(1 : nn), bu(1 : nn), xnstd(1 : nn), sb(1 : nn), maskpara(1 : nn), save_state_gauss, snew)
1938 end if
1939 !
1940 ! compute the function value at snew
1941 if (idot) write(output_unit, '(A1)') '.'
1942#ifdef MPI
1943 call mpi_comm_size(comm, nproc, ierror)
1944 do iproc = 1, nproc-1
1945 do_obj_loop = .true.
1946 call mpi_send(do_obj_loop,1, mpi_logical,iproc,0,comm,ierror)
1947 end do
1948#endif
1949 if (.not. maxit) then
1950 fnew = functn(snew, eval)
1951 icall = icall + 1_i8
1952 history(icall) = min(history(icall - 1), fnew)
1953 else
1954 fnew = -functn(snew, eval)
1955 icall = icall + 1_i8
1956 history(icall) = max(history(icall - 1), -fnew)
1957 end if
1958 !
1959 ! maximum numbers of function evaluations reached
1960 if (icall .ge. maxn) return
1961 !
1962 ! compare fnew with the worst function value fw
1963 ! fnew is greater than fw, so try a contraction step
1964 if ((fnew .gt. fw) .or. (.not. is_finite(fnew))) then
1965 do j = 1, m
1966 if (maskpara(j)) then
1967 snew(j) = ce(j) - beta * (ce(j) - sw(j))
1968 else
1969 snew(j) = s(1, j)
1970 end if
1971 end do
1972 !
1973 ! compute the function value of the contracted point
1974 if (idot) write(output_unit, '(A1)') '.'
1975#ifdef MPI
1976 call mpi_comm_size(comm, nproc, ierror)
1977 do iproc = 1, nproc-1
1978 do_obj_loop = .true.
1979 call mpi_send(do_obj_loop,1, mpi_logical,iproc,0,comm,ierror)
1980 end do
1981#endif
1982 if (.not. maxit) then
1983 fnew = functn(snew, eval)
1984 icall = icall + 1_i8
1985 history(icall) = min(history(icall - 1), fnew)
1986 else
1987 fnew = -functn(snew, eval)
1988 icall = icall + 1_i8
1989 history(icall) = max(history(icall - 1), -fnew)
1990 end if
1991 !
1992 ! maximum numbers of function evaluations reached
1993 if (icall .ge. maxn) return
1994 !
1995 ! compare fnew to the worst value fw
1996 ! if fnew is less than or equal to fw, then accept the point and return
1997 if ((fnew .gt. fw) .or. (.not. is_finite(fnew))) then
1998 !
1999 ! if both reflection and contraction fail, choose another point
2000 ! according to a normal distribution with best point of the sub-complex
2001 ! as mean and standard deviation of the population as std
2002 call getpnt(2, bl(1 : nn), bu(1 : nn), xnstd(1 : nn), sb(1 : nn), maskpara(1 : nn), save_state_gauss, snew)
2003 !
2004 ! compute the function value at the random point
2005 if (idot) write(output_unit, '(A1)') '.'
2006#ifdef MPI
2007 call mpi_comm_size(comm, nproc, ierror)
2008 do iproc = 1, nproc-1
2009 do_obj_loop = .true.
2010 call mpi_send(do_obj_loop,1, mpi_logical,iproc,0,comm,ierror)
2011 end do
2012#endif
2013 if (.not. maxit) then
2014 fnew = functn(snew, eval)
2015 icall = icall + 1_i8
2016 history(icall) = min(history(icall - 1), fnew)
2017 else
2018 fnew = -functn(snew, eval)
2019 icall = icall + 1_i8
2020 history(icall) = max(history(icall - 1), -fnew)
2021 end if
2022 !
2023 ! maximum numbers of function evaluations reached
2024 if (icall .ge. maxn) return
2025 !
2026 ! successful mutation
2027 ! replace the worst point by the new point
2028 do j = 1, m
2029 s(n, j) = snew(j)
2030 end do
2031 sf(n) = fnew
2032 end if
2033 end if
2034 !
2035 ! replace the worst point by the new point
2036 do j = 1, m
2037 s(n, j) = snew(j)
2038 end do
2039 sf(n) = fnew
2040
2041 end subroutine cce
2042
2043END MODULE mo_sce
mo_optimization_utils::eval_interface
Interface for evaluation routine.
Definition mo_optimization_utils.f90:16
mo_optimization_utils::objective_interface
Interface for objective function.
Definition mo_optimization_utils.f90:38
mo_orderpack::sort_index
Gives the indeces that would sort an array in ascending order.
Definition mo_orderpack.f90:319
mo_orderpack::sort
Sorts the input array in ascending order.
Definition mo_orderpack.f90:312
mo_string_utils::num2str
Convert to string.
Definition mo_string_utils.f90:83
mo_utils::is_finite
.true. if not IEEE Inf.
Definition mo_utils.F90:324
mo_utils::special_value
Special IEEE values.
Definition mo_utils.F90:496
mo_xor4096::get_timeseed
Generate seed for xor4096.
Definition mo_xor4096.f90:72
mo_xor4096::xor4096
Uniform XOR4096-based RNG.
Definition mo_xor4096.f90:200
mo_xor4096::xor4096g
Gaussian XOR4096-based RNG.
Definition mo_xor4096.f90:287
mo_kind
Define number representations.
Definition mo_kind.F90:17
mo_kind::i4
integer, parameter i4
4 Byte Integer Kind
Definition mo_kind.F90:40
mo_kind::i8
integer, parameter i8
8 Byte Integer Kind
Definition mo_kind.F90:42
mo_kind::dp
integer, parameter dp
Double Precision Real Kind.
Definition mo_kind.F90:46
mo_optimization_utils
Utility functions, such as interface definitions, for optimization routines.
Definition mo_optimization_utils.f90:8
mo_orderpack
Sort and ranking routines.
Definition mo_orderpack.f90:274
mo_sce
Shuffled Complex Evolution optimization algorithm.
Definition mo_sce.F90:12
mo_sce::sce
real(dp) function, dimension(size(pini, 1)), public sce(eval, functn, pini, prange, mymaxn, mymaxit, mykstop, mypcento, mypeps, myseed, myngs, mynpg, mynps, mynspl, mymings, myiniflg, myprint, mymask, myalpha, mybeta, tmp_file, popul_file, popul_file_append, parallel, restart, restart_file, bestf, neval, history)
Shuffled Complex Evolution (SCE) algorithm for global optimization.
Definition mo_sce.F90:224
mo_string_utils
String utilities.
Definition mo_string_utils.f90:11
mo_string_utils::compress
character(len(whitespaces)) function, public compress(whitespaces, n)
Remove white spaces.
Definition mo_string_utils.f90:157
mo_utils
General utilities for the CHS library.
Definition mo_utils.F90:20
mo_xor4096
XOR4096-based random number generator.
Definition mo_xor4096.f90:11
mo_xor4096::n_save_state
integer(i4), parameter, public n_save_state
Dimension of vector saving the state of a stream.
Definition mo_xor4096.f90:296