Line data Source code
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
12 : MODULE 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)) — 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 :
210 : CONTAINS
211 :
212 8 : 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 2 : 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 70 : ) 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 74 : real(dp), dimension(size(pini, 1)) :: bl ! lower bound on parameters
311 4004 : real(dp), dimension(1000) :: dummy_bl ! dummy lower bound (only for namelist)
312 66 : real(dp), dimension(size(pini, 1)) :: bu ! upper bound on parameters
313 4004 : 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 4 : real(dp) :: pcento ! percentage by which the criterion value must change
319 4 : 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 0 : 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 4 : real(dp) :: alpha ! parameter for reflection of points in complex
334 4 : 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 4 : real(dp), dimension(:), allocatable :: history_tmp ! history of best function values after each iteration
341 4 : real(dp), dimension(:), allocatable :: ihistory_tmp ! local history for OpenMP
342 4 : 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 4 : real(dp) :: fpini ! function value at initial point
352 4 : real(dp), dimension(:, :), allocatable :: x ! coordinates of points in the population
353 4 : real(dp), dimension(:), allocatable :: xf ! function values of x
354 70 : real(dp), dimension(size(pini, 1)) :: xx ! coordinates of a single point in x
355 4004 : real(dp), dimension(1000) :: dummy_xx ! dummy xx (only for namelist)
356 4 : real(dp), dimension(:, :), allocatable :: cx ! coordinates of points in a complex
357 4 : real(dp), dimension(:), allocatable :: cf ! function values of cx
358 4 : real(dp), dimension(:, :), allocatable :: s ! coordinates of points in the current sub-complex simplex
359 4 : real(dp), dimension(:), allocatable :: sf ! function values of s
360 4 : real(dp), dimension(:), allocatable :: worstx ! worst point at current shuffling loop
361 4 : real(dp) :: worstf ! function value of worstx(.)
362 4 : real(dp), dimension(:), allocatable :: xnstd ! standard deviation of parameters in the population
363 4 : real(dp) :: gnrng ! normalized geometric mean of parameter ranges
364 4 : integer(i4), dimension(:), allocatable :: lcs ! indices locating position of s(.,.) in x(.,.)
365 4 : real(dp), dimension(:), allocatable :: bound ! length of range of ith variable being optimized
366 4 : 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 4 : 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 70 : real(dp) :: denomi ! for checking improvement of last steps
389 4 : real(dp) :: timeou ! for checking improvement of last steps
390 4 : 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 4 : real(dp) :: rand ! random number
395 4 : integer(i8), dimension(:), allocatable :: iseed ! initial random seed
396 4 : integer(i8), dimension(:, :), allocatable :: save_state_unif ! save state of uniform stream
397 4 : integer(i8), dimension(:, :), allocatable :: save_state_gauss ! save state of gaussian stream
398 4 : 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 4 : real(dp), dimension(:, :), allocatable :: xtmp ! tmp array for complex reduction
402 4 : real(dp), dimension(:), allocatable :: ftmp ! %
403 4 : real(dp) :: large ! for treating NaNs
404 : logical :: ipopul_file_append
405 : integer(i4) :: nonan ! # of non-NaN in history_tmp
406 4 : 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 4 : if (present(restart)) then
426 2 : irestart = restart
427 : else
428 : irestart = .false.
429 : end if
430 :
431 4 : if (present(restart_file)) then
432 4 : isrestart_file = restart_file
433 : else
434 0 : isrestart_file = 'mo_sce.restart'
435 : end if
436 :
437 4 : 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 2 : if (present(parallel)) then
448 2 : parall = parallel
449 : else
450 0 : parall = .false.
451 : end if
452 : #endif
453 :
454 2 : if (parall) then
455 : ! OpenMP or not
456 0 : 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 2 : n_threads = 1
465 : end if
466 :
467 : ! One random number chain per OpenMP thread
468 6 : allocate(rand_tmp(n_threads))
469 4 : allocate(iseed(n_threads))
470 6 : allocate(save_state_unif(n_threads, n_save_state))
471 4 : allocate(save_state_gauss(n_threads, n_save_state))
472 :
473 2 : if (present(mymask)) then
474 2 : if (.not. any(mymask)) stop 'mo_sce: all parameters are masked --> none will be optimized'
475 35 : maskpara = mymask
476 : else
477 0 : maskpara = .true.
478 : end if
479 :
480 : ! number of parameters to optimize
481 35 : nopt = count(maskpara, dim = 1)
482 : ! total number of parameters
483 2 : nn = size(pini, 1)
484 :
485 : ! input checking
486 2 : if (size(prange, dim = 1) .ne. size(pini, 1)) then
487 0 : stop 'mo_sce: prange has not matching rows'
488 : end if
489 2 : if (size(prange, dim = 2) .ne. 2) then
490 0 : stop 'mo_sce: two colums expected for prange'
491 : end if
492 35 : bl(:) = prange(:, 1)
493 35 : bu(:) = prange(:, 2)
494 35 : do ii = 1, nn
495 35 : if(((bu(ii) - bl(ii)) .lt. tiny(1.0_dp)) .and. maskpara(ii)) then
496 0 : write(error_unit, *) 'para #', ii, ' :: range = ( ', bl(ii), ' , ', bu(ii), ' )'
497 0 : stop 'mo_sce: inconsistent or too small parameter range'
498 : end if
499 : end do
500 :
501 : !
502 : ! optionals checking
503 2 : if (present(mymaxn)) then
504 2 : if (mymaxn .lt. 2_i4) stop 'mo_sce: maxn has to be at least 2'
505 2 : maxn = mymaxn
506 : else
507 0 : maxn = 1000_i8
508 : end if
509 2 : if (present(mymaxit)) then
510 2 : maxit = mymaxit
511 : else
512 0 : maxit = .false.
513 : end if
514 2 : if (present(mykstop)) then
515 2 : if (mykstop .lt. 1_i4) stop 'mo_sce: kstop has to be at least 1'
516 2 : kstop = mykstop
517 : else
518 0 : kstop = 10_i4
519 : end if
520 2 : if (present(mypcento)) then
521 2 : if (mypcento .lt. 0_dp) stop 'mo_sce: pcento should be positive'
522 2 : pcento = mypcento
523 : else
524 0 : pcento = 0.0001_dp
525 : end if
526 2 : if (present(mypeps)) then
527 2 : if (mypeps .lt. 0_dp) stop 'mo_sce: peps should be positive'
528 2 : peps = mypeps
529 : else
530 0 : peps = 0.001_dp
531 : end if
532 2 : if (present(myseed)) then
533 2 : if (myseed .lt. 1_i8) stop 'mo_sce: seed should be non-negative'
534 4 : forall(ii = 1 : n_threads) iseed(ii) = myseed + (ii - 1) * 1000_i8
535 : else
536 0 : call get_timeseed(iseed)
537 : end if
538 2 : if (present(myngs)) then
539 2 : if (myngs .lt. 1_i4) stop 'mo_sce: ngs has to be at least 1'
540 2 : ngs = myngs
541 : else
542 0 : ngs = 2_i4
543 : end if
544 2 : if (present(mynpg)) then
545 2 : if (mynpg .lt. 3_i4) stop 'mo_sce: npg has to be at least 3'
546 2 : npg = mynpg
547 : else
548 0 : npg = 2 * nopt + 1
549 : end if
550 2 : if (present(mynps)) then
551 2 : if (mynps .lt. 2_i4) stop 'mo_sce: nps has to be at least 2'
552 2 : nps = mynps
553 : else
554 0 : nps = nopt + 1_i4
555 : end if
556 2 : if (present(mynspl)) then
557 2 : if (mynspl .lt. 3_i4) stop 'mo_sce: nspl has to be at least 3'
558 2 : nspl = mynspl
559 : else
560 0 : nspl = 2 * nopt + 1
561 : end if
562 2 : if (present(mymings)) then
563 2 : if (mymings .lt. 1_i4) stop 'mo_sce: mings has to be at least 1'
564 2 : mings = mymings
565 : else
566 0 : mings = ngs ! no reduction of complexes
567 : end if
568 2 : if (present(myiniflg)) then
569 2 : if ((myiniflg .ne. 1_i4) .and. (myiniflg .ne. 0_i4)) stop 'mo_sce: iniflg has to be 0 or 1'
570 2 : iniflg = myiniflg
571 : else
572 0 : iniflg = 1_i4
573 : end if
574 2 : idot = .false.
575 2 : if (present(myprint)) then
576 2 : if ((myprint .lt. 0_i4) .or. (myprint .gt. 4_i4)) stop 'mo_sce: iprint has to be between 0 and 4'
577 2 : if (myprint > 2_i4) then
578 0 : iprint = myprint - 3
579 0 : idot = .true.
580 : else
581 2 : iprint = myprint
582 : end if
583 : else
584 0 : iprint = 2_i4 ! no printing
585 0 : iprint = 0_i4
586 : end if
587 2 : if (present(myalpha)) then
588 2 : alpha = myalpha
589 : else
590 0 : alpha = 0.8_dp
591 : end if
592 2 : if (present(mybeta)) then
593 2 : beta = mybeta
594 : else
595 0 : beta = 0.45_dp
596 : end if
597 :
598 2 : if (present(popul_file_append)) then
599 0 : ipopul_file_append = popul_file_append
600 : else
601 2 : ipopul_file_append = .false.
602 : end if
603 :
604 2 : if (present(tmp_file)) then
605 2 : itmp_file = .true.
606 2 : istmp_file = tmp_file
607 2 : open(999, file = trim(adjustl(istmp_file)), action = 'write', status = 'unknown')
608 2 : write(999, *) '# settings :: general'
609 2 : write(999, *) '# nIterations seed'
610 2 : write(999, *) maxn, iseed
611 2 : write(999, *) '# settings :: sce specific'
612 2 : write(999, *) '# sce_ngs sce_npg sce_nps'
613 2 : write(999, *) ngs, npg, nps
614 2 : write(999, *) '# nloop icall ngs1 bestf worstf gnrng (bestx(j),j=1,nn)'
615 2 : close(999)
616 : else
617 0 : itmp_file = .false.
618 0 : istmp_file = ''
619 : end if
620 :
621 2 : if (present(popul_file)) then
622 2 : ipopul_file = .true.
623 2 : ispopul_file = popul_file
624 : else
625 0 : ipopul_file = .false.
626 0 : ispopul_file = ''
627 : end if
628 :
629 2 : if (ipopul_file .and. (.not. ipopul_file_append)) then
630 2 : open(999, file = trim(adjustl(ispopul_file)), action = 'write', status = 'unknown')
631 2 : write(999, *) '# xf(i) (x(i,j),j=1,nn)'
632 2 : close(999)
633 : end if
634 :
635 : ! allocation of arrays
636 8 : allocate(x(ngs * npg, nn))
637 6 : allocate(xf(ngs * npg))
638 6 : allocate(worstx(nn))
639 4 : allocate(xnstd(nn))
640 4 : allocate(bound(nn))
641 4 : allocate(unit(nn))
642 6 : allocate(criter(kstop + 1))
643 6 : allocate(xname(nn))
644 6 : allocate(history_tmp(maxn + 3 * ngs * nspl))
645 8 : allocate(xtmp(npg, nn))
646 6 : allocate(ftmp(npg))
647 2 : if (maxit) then
648 : large = -0.5_dp * huge(1.0_dp)
649 : else
650 2 : large = 0.5_dp * huge(1.0_dp)
651 : end if
652 60410 : history_tmp(:) = large
653 24 : criter(:) = large
654 : !
655 : ! initialize variables
656 35 : do ii = 1, nn
657 35 : xname(ii) = compress('p' // num2str(ii, '(I3.3)'))
658 : end do
659 :
660 2 : nloop = 0_i4
661 2 : loop = 0_i4
662 2 : igs = 0_i4
663 : !
664 : ! compute the total number of points in initial population
665 2 : npt = ngs * npg
666 2 : ngs1 = ngs
667 2 : npt1 = npt
668 : !
669 2 : if (iprint .lt. 2) then
670 2 : write(output_unit, *) '=================================================='
671 2 : write(output_unit, *) 'ENTER THE SHUFFLED COMPLEX EVOLUTION GLOBAL SEARCH'
672 2 : write(output_unit, *) '=================================================='
673 : end if
674 : !
675 : ! Print seed
676 2 : if (iprint .lt. 2) then
677 2 : write (*, *) ' Seeds used : ', iseed
678 : end if
679 : ! initialize random number generator: Uniform
680 2 : call xor4096(iseed, rand_tmp, save_state = save_state_unif)
681 : ! initialize random number generator: Gaussian
682 4 : iseed = iseed + 1000_i8
683 : call xor4096g(iseed, rand_tmp, save_state = save_state_gauss)
684 4 : iseed = 0_i8
685 : !
686 : ! compute the bound for parameters being optimized
687 35 : do ii = 1, nn
688 33 : bound(ii) = bu(ii) - bl(ii)
689 35 : 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 2 : 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 2 : if (.not. maxit) then
705 2 : fpini = functn(pini, eval)
706 2 : history_tmp(1) = fpini
707 : else
708 0 : fpini = -functn(pini, eval)
709 0 : history_tmp(1) = -fpini
710 : end if
711 :
712 : ! print the initial point and its criterion value
713 35 : bestx = pini
714 2 : bestf_tmp = fpini
715 2 : if (iprint .lt. 2) then
716 2 : write(output_unit, *) ''
717 2 : write(output_unit, *) ''
718 2 : write(output_unit, *) '*** PRINT THE INITIAL POINT AND ITS CRITERION VALUE ***'
719 2 : call write_best_final()
720 : end if
721 2 : if (itmp_file) then ! initial tmp file
722 2 : open(999, file = trim(adjustl(istmp_file)), action = 'write', position = 'append', recl = (nn + 6) * 30)
723 2 : if (.not. maxit) then
724 2 : write(999, *) 0, 1, ngs1, fpini, fpini, 1.0_dp, pini
725 : else
726 0 : write(999, *) 0, 1, ngs1, -fpini, -fpini, 1.0_dp, pini
727 : end if
728 2 : 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 2 : if (iniflg .eq. 1) then
734 35 : do ii = 1, nn
735 35 : x(1, ii) = pini(ii)
736 : end do
737 2 : xf(1) = fpini
738 : !
739 : ! else, generate a point randomly and set it equal to x(1,.)
740 : else
741 0 : itmp = save_state_unif(1, :)
742 0 : call getpnt(1, bl(1 : nn), bu(1 : nn), unit(1 : nn), pini(1 : nn), maskpara, itmp, xx)
743 0 : save_state_unif(1, :) = itmp
744 0 : do ii = 1, nn
745 0 : x(1, ii) = xx(ii)
746 : end do
747 0 : 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 0 : if (.not. maxit) then
756 0 : xf(1) = functn(xx, eval)
757 : else
758 0 : xf(1) = -functn(xx, eval)
759 : end if
760 : end if
761 : !
762 : ! count function evaluation of the first point
763 2 : icall = 1_i8
764 : ! if (icall .ge. maxn) return
765 : !
766 2 : 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 0 : ithread = 1
774 : !$OMP parallel default(shared) private(ii,jj,ithread,xx)
775 : !$OMP do
776 0 : do ii = 2, npt1
777 : !$ ithread = OMP_GET_THREAD_NUM() + 1
778 0 : itmp = save_state_unif(ithread, :)
779 0 : call getpnt(1, bl(1 : nn), bu(1 : nn), unit(1 : nn), pini(1 : nn), maskpara, itmp, xx)
780 0 : save_state_unif(ithread, :) = itmp
781 0 : do jj = 1, nn
782 0 : x(ii, jj) = xx(jj)
783 : end do
784 0 : 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 0 : if (.not. maxit) then
793 0 : xf(ii) = functn(xx, eval)
794 0 : history_tmp(ii) = xf(ii) ! min(history_tmp(ii-1),xf(ii)) --> will be sorted later
795 : else
796 0 : xf(ii) = -functn(xx, eval)
797 0 : 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 2 : ithread = 1
811 :
812 136 : do ii = 2, npt1
813 134 : call getpnt(1, bl(1 : nn), bu(1 : nn), unit(1 : nn), pini(1 : nn), maskpara, save_state_unif(ithread, :), xx)
814 3803 : do jj = 1, nn
815 3803 : x(ii, jj) = xx(jj)
816 : end do
817 134 : 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 134 : if (.not. maxit) then
826 134 : xf(ii) = functn(xx, eval)
827 134 : icall = icall + 1_i8
828 134 : history_tmp(icall) = xf(ii) ! min(history_tmp(icall-1),xf(ii)) --> will be sorted later
829 : else
830 0 : xf(ii) = -functn(xx, eval)
831 0 : icall = icall + 1_i8
832 0 : history_tmp(icall) = -xf(ii) ! max(history_tmp(icall-1),-xf(ii)) --> will be sorted later
833 : end if
834 :
835 136 : if (icall .ge. maxn) then
836 0 : npt1 = ii
837 0 : if (iprint .lt. 2) then
838 : ! maximum trials reached
839 0 : call write_termination_case(1)
840 0 : call write_best_final()
841 : end if
842 : exit
843 : end if
844 : end do
845 :
846 : end if
847 :
848 2 : icall = int(npt1, i8)
849 : !
850 : ! arrange the points in order of increasing function value
851 2 : if (maxit) then
852 0 : large = minval(xf(1 : npt1))
853 0 : large = merge(0.9_dp * large, 1.1_dp * large, large>0._dp)
854 : else
855 140 : large = maxval(xf(1 : npt1))
856 2 : large = merge(1.1_dp * large, 0.9_dp * large, large>0._dp)
857 : end if
858 274 : 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 138 : nonan = size(pack(history_tmp(1 : npt1), mask = is_finite(history_tmp(1 : npt1))))
862 2 : if (nonan /= npt1) then
863 0 : allocate(htmp(nonan))
864 0 : htmp(1 : nonan) = pack(history_tmp(1 : npt1), mask = is_finite(history_tmp(1 : npt1)))
865 0 : call sort(htmp(1 : nonan))
866 0 : history_tmp(1 : nonan) = htmp(1 : nonan)
867 0 : history_tmp(nonan + 1 : npt1) = special_value(1.0_dp, 'IEEE_QUIET_NAN')
868 0 : deallocate(htmp)
869 : else
870 2 : call sort(history_tmp(1 : npt1))
871 : end if
872 2 : call sort_matrix(x(1 : npt1, 1 : nn), xf(1 : npt1))
873 : !
874 : ! record the best and worst points
875 35 : do ii = 1, nn
876 33 : bestx(ii) = x(1, ii)
877 35 : worstx(ii) = x(npt1, ii)
878 : end do
879 2 : bestf_tmp = xf(1)
880 2 : worstf = xf(npt1)
881 : !
882 : ! compute the parameter range for the initial population
883 2 : 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 2 : if (itmp_file) then
887 2 : call write_best_intermediate(.true.)
888 : end if
889 :
890 : ! write population to file
891 2 : if (ipopul_file) then
892 2 : call write_population(.true.)
893 : end if
894 : !
895 : ! print the results for the initial population
896 2 : print : if (iprint .lt. 2) then
897 : ! write currently best x and xf to screen
898 2 : call write_best_intermediate(.false.)
899 :
900 : ! write population on screen
901 2 : if (iprint .eq. 1) then
902 0 : call write_population(.false.)
903 : end if
904 : end if print
905 : !
906 : ! Maximum number of function evaluations reached?
907 2 : if (icall .ge. maxn) then
908 0 : if (iprint .lt. 2) then
909 : ! maximum trials reached
910 0 : call write_termination_case(1)
911 0 : call write_best_final()
912 : end if
913 0 : 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 0 : return
925 : end if
926 : !
927 : ! Feasible parameter space converged?
928 2 : if (ipcnvg .eq. 1) then
929 0 : if (iprint .lt. 2) then
930 : ! converged because feasible parameter space small
931 0 : call write_termination_case(3)
932 0 : call write_best_final()
933 : end if
934 0 : 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 0 : return
946 : end if
947 : !
948 : ! transfer all array-like variables in namelist to fixed-size dummy-arrays
949 2 : dummy_maskpara = .false.
950 35 : dummy_maskpara(1 : size(pini, 1)) = maskpara
951 2002 : dummy_bl = -9999.0_dp
952 35 : dummy_bl(1 : size(pini, 1)) = bl
953 2002 : dummy_bu = -9999.0_dp
954 35 : dummy_bu(1 : size(pini, 1)) = bu
955 2002 : dummy_xx = -9999.0_dp
956 35 : dummy_xx(1 : size(pini, 1)) = xx
957 : !
958 : ! write restart
959 2 : open(999, file = isrestart_file, status = 'unknown', action = 'write', delim = 'quote')
960 2 : write(999, restartnml1)
961 2 : write(999, restartnml2)
962 6 : close(999)
963 : end if ! restart or not
964 :
965 : if (irestart) then
966 : ! read 1st namelist with allocated/scalar variables
967 2 : open(999, file = isrestart_file, status = 'old', action = 'read', delim = 'quote')
968 2 : read(999, nml = restartnml1)
969 2 : close(999)
970 :
971 : ! transfer all array-like variables in namelist to fixed-size dummy-arrays
972 35 : maskpara = dummy_maskpara(1 : size(pini, 1))
973 35 : bl = dummy_bl(1 : size(pini, 1))
974 35 : bu = dummy_bu(1 : size(pini, 1))
975 35 : xx = dummy_xx(1 : size(pini, 1))
976 :
977 : ! allocate global arrays
978 6 : allocate(rand_tmp(n_threads))
979 4 : allocate(iseed(n_threads))
980 6 : allocate(save_state_unif(n_threads, n_save_state))
981 4 : allocate(save_state_gauss(n_threads, n_save_state))
982 8 : allocate(x(ngs * npg, nn))
983 6 : allocate(xf(ngs * npg))
984 6 : allocate(worstx(nn))
985 4 : allocate(xnstd(nn))
986 4 : allocate(bound(nn))
987 4 : allocate(unit(nn))
988 6 : allocate(criter(kstop + 1))
989 6 : allocate(xname(nn))
990 6 : allocate(history_tmp(maxn + 3 * ngs * nspl))
991 8 : allocate(xtmp(npg, nn))
992 6 : allocate(ftmp(npg))
993 : end if
994 :
995 : !
996 : ! begin the main loop
997 48 : mainloop : do while (icall .lt. maxn)
998 : ! read variables from restart files
999 48 : open(999, file = isrestart_file, status = 'old', action = 'read', delim = 'quote')
1000 48 : read(999, nml = restartnml1)
1001 48 : read(999, nml = restartnml2)
1002 48 : close(999)
1003 : !
1004 : ! transfer all array-like variables in namelist to fixed-size dummy-arrays
1005 1137 : maskpara = dummy_maskpara(1 : size(pini, 1))
1006 1137 : bl = dummy_bl(1 : size(pini, 1))
1007 1137 : bu = dummy_bu(1 : size(pini, 1))
1008 1137 : xx = dummy_xx(1 : size(pini, 1))
1009 : !
1010 48 : nloop = nloop + 1
1011 : !
1012 48 : if (iprint .lt. 2) then
1013 48 : write(output_unit, *) ''
1014 48 : 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 48 : if (parall) then
1021 :
1022 : ! -----------------------------------------------------------------------
1023 : ! Parallel version of complex-loop
1024 : ! -----------------------------------------------------------------------
1025 :
1026 0 : 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 0 : allocate(cx(npg, nn))
1031 0 : allocate(cf(npg))
1032 0 : allocate(lcs(nps))
1033 0 : allocate(s(nps, nn))
1034 0 : allocate(sf(nps))
1035 0 : allocate(ihistory_tmp(maxn + 3 * ngs * nspl))
1036 : !$OMP do
1037 0 : comploop_omp : do igs = 1, ngs1
1038 : !$ ithread = OMP_GET_THREAD_NUM() + 1
1039 : !
1040 : ! assign points into complexes
1041 0 : do k1 = 1, npg
1042 0 : k2 = (k1 - 1) * ngs1 + igs
1043 0 : do jj = 1, nn
1044 0 : cx(k1, jj) = x(k2, jj)
1045 : end do
1046 0 : cf(k1) = xf(k2)
1047 : end do
1048 : !
1049 : ! begin inner loop - random selection of sub-complexes
1050 : ! <alpha> loop from duan(1993)
1051 0 : 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 0 : if (nps .eq. npg) then
1059 0 : do kk = 1, nps
1060 0 : lcs(kk) = kk
1061 : end do
1062 : else
1063 0 : call xor4096(iseed(ithread), rand, save_state = save_state_unif(ithread, :))
1064 0 : lcs(1) = 1 + int(real(npg, dp) + 0.5_dp &
1065 0 : - sqrt((real(npg, dp) + .5_dp)**2 - real(npg, dp) * real(npg + 1, dp) * rand))
1066 0 : do kk = 2, nps
1067 : lpos_ok = .false.
1068 0 : do while (.not. lpos_ok)
1069 0 : lpos_ok = .true.
1070 0 : call xor4096(iseed(ithread), rand, save_state = save_state_unif(ithread, :))
1071 : lpos = 1 + int(real(npg, dp) + 0.5_dp &
1072 0 : - 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 0 : do k1 = 1, kk - 1
1075 0 : if (lpos .eq. lcs(k1)) lpos_ok = (lpos_ok .and. .false.)
1076 : end do
1077 : end do
1078 0 : lcs(kk) = lpos
1079 : end do
1080 : !
1081 : ! arrange the sub-complex in order of increasing function value
1082 0 : call sort(lcs(1 : nps))
1083 : end if
1084 : !
1085 : ! create the sub-complex arrays
1086 0 : do kk = 1, nps
1087 0 : do jj = 1, nn
1088 0 : s(kk, jj) = cx(lcs(kk), jj)
1089 : end do
1090 0 : sf(kk) = cf(lcs(kk))
1091 : end do
1092 : !
1093 : ! remember largest for treating of NaNs
1094 0 : if (maxit) then
1095 0 : large = minval(cf(1 : npg))
1096 0 : large = merge(0.9_dp * large, 1.1_dp * large, large>0._dp)
1097 : else
1098 0 : large = maxval(cf(1 : npg))
1099 0 : 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 0 : icall_merk = icall
1104 0 : iicall = icall
1105 0 : 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 0 : call cce(s(1 : nps, 1 : nn), sf(1 : nps), bl(1 : nn), bu(1 : nn), maskpara, xnstd(1 : nn), &
1111 0 : iicall, maxn, maxit, save_state_gauss(ithread, :), functn, eval, alpha, beta, ihistory_tmp, idot)
1112 : #endif
1113 0 : history_tmp(icall + 1 : icall + (iicall - icall_merk)) = ihistory_tmp(icall_merk + 1 : iicall)
1114 0 : icall = icall + (iicall - icall_merk)
1115 : !
1116 : ! if the sub-complex is accepted, replace the new sub-complex
1117 : ! into the complex
1118 0 : do kk = 1, nps
1119 0 : do jj = 1, nn
1120 0 : cx(lcs(kk), jj) = s(kk, jj)
1121 : end do
1122 0 : cf(lcs(kk)) = sf(kk)
1123 : end do
1124 : !
1125 : ! sort the points
1126 0 : cf(1 : npg) = merge(cf(1 : npg), large, is_finite(cf(1 : npg))) ! NaN and Infinite
1127 0 : 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 0 : do k1 = 1, npg
1136 0 : k2 = (k1 - 1) * ngs1 + igs
1137 0 : do jj = 1, nn
1138 0 : x(k2, jj) = cx(k1, jj)
1139 : end do
1140 0 : 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 0 : deallocate(cx)
1148 0 : deallocate(cf)
1149 0 : deallocate(lcs)
1150 0 : deallocate(s)
1151 0 : deallocate(sf)
1152 0 : deallocate(ihistory_tmp)
1153 : !$OMP end parallel
1154 :
1155 : else
1156 :
1157 : ! -----------------------------------------------------------------------
1158 : ! Non-parallel version of complex-loop
1159 : ! -----------------------------------------------------------------------
1160 :
1161 192 : allocate(cx(npg, nn))
1162 144 : allocate(cf(npg))
1163 144 : allocate(lcs(nps))
1164 192 : allocate(s(nps, nn))
1165 144 : allocate(sf(nps))
1166 :
1167 48 : ithread = 1
1168 :
1169 144 : comploop : do igs = 1, ngs1
1170 : !
1171 : ! assign points into complexes
1172 4548 : do k1 = 1, npg
1173 4452 : k2 = (k1 - 1) * ngs1 + igs
1174 133098 : do jj = 1, nn
1175 133098 : cx(k1, jj) = x(k2, jj)
1176 : end do
1177 4548 : cf(k1) = xf(k2)
1178 : end do
1179 : !
1180 : ! begin inner loop - random selection of sub-complexes
1181 : ! <alpha> loop from duan(1993)
1182 4548 : 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 4452 : if (nps .eq. npg) then
1190 0 : do kk = 1, nps
1191 0 : lcs(kk) = kk
1192 : end do
1193 : else
1194 4452 : call xor4096(iseed(ithread), rand, save_state = save_state_unif(ithread, :))
1195 4452 : lcs(1) = 1 + int(real(npg, dp) + 0.5_dp &
1196 4452 : - sqrt((real(npg, dp) + .5_dp)**2 - real(npg, dp) * real(npg + 1, dp) * rand))
1197 133098 : do kk = 2, nps
1198 : lpos_ok = .false.
1199 333769 : do while (.not. lpos_ok)
1200 205123 : lpos_ok = .true.
1201 205123 : call xor4096(iseed(ithread), rand, save_state = save_state_unif(ithread, :))
1202 : lpos = 1 + int(real(npg, dp) + 0.5_dp &
1203 205123 : - 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 3989114 : do k1 = 1, kk - 1
1206 3860468 : if (lpos .eq. lcs(k1)) lpos_ok = (lpos_ok .and. .false.)
1207 : end do
1208 : end do
1209 133098 : lcs(kk) = lpos
1210 : end do
1211 : !
1212 : ! arrange the sub-complex in order of increasing function value
1213 4452 : call sort(lcs(1 : nps))
1214 : end if
1215 : !
1216 : ! create the sub-complex arrays
1217 137550 : do kk = 1, nps
1218 4106382 : do jj = 1, nn
1219 4106382 : s(kk, jj) = cx(lcs(kk), jj)
1220 : end do
1221 137550 : sf(kk) = cf(lcs(kk))
1222 : end do
1223 : !
1224 : ! remember largest for treating of NaNs
1225 4452 : if (maxit) then
1226 0 : large = minval(cf(1 : npg))
1227 0 : large = merge(0.9_dp * large, 1.1_dp * large, large>0._dp)
1228 : else
1229 270648 : large = maxval(cf(1 : npg))
1230 4452 : 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 17808 : call cce(s(1 : nps, 1 : nn), sf(1 : nps), bl(1 : nn), bu(1 : nn), maskpara, xnstd(1 : nn), &
1239 22260 : 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 137550 : do kk = 1, nps
1245 4106382 : do jj = 1, nn
1246 4106382 : cx(lcs(kk), jj) = s(kk, jj)
1247 : end do
1248 137550 : cf(lcs(kk)) = sf(kk)
1249 : end do
1250 : !
1251 : ! sort the points
1252 527940 : cf(1 : npg) = merge(cf(1 : npg), large, is_finite(cf(1 : npg))) ! NaN and Infinite
1253 4452 : 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 4548 : if (icall .ge. maxn) then
1257 0 : if (iprint .lt. 2) then
1258 : ! maximum trials reached
1259 0 : call write_termination_case(1)
1260 0 : 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 4548 : do k1 = 1, npg
1269 4452 : k2 = (k1 - 1) * ngs1 + igs
1270 133098 : do jj = 1, nn
1271 133098 : x(k2, jj) = cx(k1, jj)
1272 : end do
1273 4548 : xf(k2) = cf(k1)
1274 : end do
1275 144 : if (icall .ge. maxn) then
1276 0 : if (iprint .lt. 2) then
1277 : ! maximum trials reached
1278 0 : call write_termination_case(1)
1279 0 : 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 48 : deallocate(cx)
1288 48 : deallocate(cf)
1289 48 : deallocate(lcs)
1290 48 : deallocate(s)
1291 48 : deallocate(sf)
1292 :
1293 : end if ! end parallel/non-parallel region
1294 : !
1295 : ! re-sort the points
1296 48 : call sort_matrix(x(1 : npt1, 1 : nn), xf(1 : npt1))
1297 : !
1298 : ! record the best and worst points
1299 1137 : do jj = 1, nn
1300 1089 : bestx(jj) = x(1, jj)
1301 1137 : worstx(jj) = x(npt1, jj)
1302 : end do
1303 48 : bestf_tmp = xf(1)
1304 48 : worstf = xf(npt1)
1305 : !
1306 : ! test the population for parameter convergence
1307 48 : 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 48 : if (itmp_file) then
1311 48 : call write_best_intermediate(.true.)
1312 : end if
1313 : !
1314 : ! write population to file
1315 48 : if (ipopul_file) then
1316 48 : call write_population(.true.)
1317 : end if
1318 : !
1319 : ! print the results for current population
1320 48 : if (iprint .lt. 2) then
1321 48 : call write_best_intermediate(.false.)
1322 :
1323 48 : if (iprint .eq. 1) then
1324 0 : 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 48 : if (icall .ge. maxn) then
1332 0 : if (iprint .lt. 2) then
1333 : ! maximum trials reached
1334 0 : call write_termination_case(1)
1335 0 : call write_best_final()
1336 : end if
1337 0 : 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 0 : return
1349 : end if
1350 : !
1351 : ! compute the count on successive loops w/o function improvement
1352 48 : criter(kstop + 1) = bestf_tmp
1353 48 : if (nloop .gt. kstop) then
1354 28 : denomi = 0.5_dp * abs(criter((kstop + 1) - kstop) + criter(kstop + 1))
1355 28 : denomi = max(denomi, 1.0e-15_dp)
1356 28 : timeou = abs(criter((kstop + 1) - kstop) - criter(kstop + 1)) / denomi
1357 28 : if (timeou .lt. pcento) then
1358 2 : if (iprint .lt. 2) then
1359 : ! criterion value has not changed during last loops
1360 2 : call write_termination_case(2)
1361 2 : call write_best_final()
1362 : end if
1363 2 : 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 2 : return
1375 : end if
1376 : end if
1377 506 : criter(1 : kstop) = criter(2 : kstop + 1)
1378 : !
1379 : ! if population is converged into a sufficiently small space
1380 46 : if (ipcnvg .eq. 1) then
1381 2 : if (iprint .lt. 2) then
1382 : ! converged because feasible parameter space small
1383 2 : call write_termination_case(3)
1384 2 : call write_best_final()
1385 : end if
1386 2 : 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 2 : return
1398 : end if
1399 : !
1400 : ! none of the stopping criteria is satisfied, continue search
1401 : !
1402 : ! check for complex number reduction
1403 44 : if (ngs1 .gt. mings) then
1404 0 : ngs2 = ngs1
1405 0 : ngs1 = ngs1 - 1
1406 0 : npt1 = ngs1 * npg
1407 0 : 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 44 : dummy_maskpara = .false.
1412 1067 : dummy_maskpara(1 : size(pini, 1)) = maskpara
1413 44044 : dummy_bl = -9999.0_dp
1414 1067 : dummy_bl(1 : size(pini, 1)) = bl
1415 44044 : dummy_bu = -9999.0_dp
1416 1067 : dummy_bu(1 : size(pini, 1)) = bu
1417 44044 : dummy_xx = -9999.0_dp
1418 1067 : dummy_xx(1 : size(pini, 1)) = xx
1419 : !
1420 : ! write restart
1421 44 : open(999, file = isrestart_file, status = 'unknown', action = 'write', delim = 'quote')
1422 44 : write(999, restartnml1)
1423 44 : write(999, restartnml2)
1424 44 : close(999)
1425 : end do mainloop
1426 :
1427 0 : deallocate(xtmp)
1428 0 : deallocate(ftmp)
1429 :
1430 12 : 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 100 : subroutine write_best_intermediate(to_file)
1442 :
1443 : implicit none
1444 :
1445 : logical, intent(in) :: to_file
1446 : integer(i4) :: ll
1447 :
1448 100 : if (to_file) then
1449 50 : open(999, file = trim(adjustl(istmp_file)), action = 'write', position = 'append', recl = (nn + 6) * 30)
1450 50 : if (.not. maxit) then
1451 1172 : write(999, *) nloop, icall, ngs1, bestf_tmp, worstf, gnrng, (bestx(ll), ll = 1, nn)
1452 : else
1453 0 : write(999, *) nloop, icall, ngs1, -bestf_tmp, -worstf, gnrng, (bestx(ll), ll = 1, nn)
1454 : end if
1455 50 : close(999)
1456 : else
1457 50 : write(format_str1, '(A13,I3,A8)') '( A49, ', nn, '(6X,A4))'
1458 50 : write(format_str2, '(A26,I3,A8)') '(I5,1X,I5,3X,I5,3(E22.14), ', nn, '(G10.3))'
1459 50 : if (nloop == 0) then
1460 2 : write(output_unit, *) ''
1461 2 : write(output_unit, '(A44)') ' *** PRINT THE RESULTS OF THE SCE SEARCH ***'
1462 2 : write(output_unit, *) ''
1463 35 : write(output_unit, format_str1) ' LOOP TRIALS COMPLXS BEST F WORST F PAR RNG ', (xname(ll), ll = 1, nn)
1464 : end if
1465 50 : if (.not. maxit) then
1466 1172 : write(output_unit, format_str2) nloop, icall, ngs1, bestf_tmp, worstf, gnrng, (bestx(ll), ll = 1, nn)
1467 : else
1468 0 : 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 4 : end subroutine write_best_intermediate
1473 :
1474 6 : subroutine write_best_final()
1475 :
1476 : implicit none
1477 :
1478 : integer(i4) :: ll
1479 :
1480 6 : write(format_str1, '(A13,I3,A8)') '( A10, ', nn, '(6X,A4))'
1481 6 : write(format_str2, '(A14,I3,A8)') '(E22.14, ', nn, '(G10.3))'
1482 :
1483 105 : write(output_unit, format_str1) 'CRITERION ', (xname(ll), ll = 1, nn)
1484 6 : if (.not. maxit) then
1485 105 : write(output_unit, format_str2) bestf_tmp, (bestx(ll), ll = 1, nn)
1486 : else
1487 0 : write(output_unit, format_str2) -bestf_tmp, (bestx(ll), ll = 1, nn)
1488 : end if
1489 :
1490 100 : end subroutine write_best_final
1491 :
1492 50 : subroutine write_population(to_file)
1493 :
1494 : logical, intent(in) :: to_file
1495 :
1496 : integer(i4) :: ll, mm
1497 :
1498 50 : if (to_file) then
1499 50 : write(format_str2, '(A13,I3,A9)') '(I4, E22.14, ', nn, '(E22.14))'
1500 50 : open(unit = 999, file = trim(adjustl(ispopul_file)), action = 'write', position = 'append', recl = (nn + 2) * 30)
1501 50 : if (.not. maxit) then
1502 4638 : do mm = 1, npt1
1503 136986 : write(999, *) nloop, xf(mm), (x(mm, ll), ll = 1, nn)
1504 : end do
1505 : else
1506 0 : do mm = 1, npt1
1507 0 : write(999, *) nloop, -xf(mm), (x(mm, ll), ll = 1, nn)
1508 : end do
1509 : end if
1510 50 : close(999)
1511 : else
1512 0 : write(format_str2, '(A12,I3,A8)') '(I4, E22.14, ', nn, '(G10.3))'
1513 0 : write(output_unit, *) ''
1514 0 : write(output_unit, '(A22,I3)') ' POPULATION AT LOOP ', nloop
1515 0 : write(output_unit, '(A27)') '---------------------------'
1516 0 : if (.not. maxit) then
1517 0 : do mm = 1, npt1
1518 0 : write(output_unit, format_str2) nloop, xf(mm), (x(mm, ll), ll = 1, nn)
1519 : end do
1520 : else
1521 0 : do mm = 1, npt1
1522 0 : 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 6 : end subroutine write_population
1528 :
1529 4 : subroutine write_termination_case(case)
1530 :
1531 : integer(i4), intent(in) :: case
1532 :
1533 4 : select case (case)
1534 : case (1) ! maximal number of iterations reached
1535 0 : write(output_unit, *) ''
1536 : write(output_unit, '(A46,A39,I7,A46,I4,A12,I4,A19,I4,A4)') &
1537 0 : '*** OPTIMIZATION SEARCH TERMINATED BECAUSE THE', &
1538 0 : ' LIMIT ON THE MAXIMUM NUMBER OF TRIALS ', maxn, &
1539 0 : ' EXCEEDED. SEARCH WAS STOPPED AT SUB-COMPLEX ', &
1540 0 : loop, ' OF COMPLEX ', igs, ' IN SHUFFLING LOOP ', nloop, ' ***'
1541 : !
1542 : case(2) ! objective not changed during last evolution loops
1543 2 : write(output_unit, *) ''
1544 2 : write(output_unit, '(A72,F8.4,A12,I3,A20)') '*** OPTIMIZATION TERMINATED BECAUSE THE CRITERION VALUE HAS NOT CHANGED ', &
1545 4 : pcento * 100._dp, ' PERCENT IN ', kstop, ' SHUFFLING LOOPS ***'
1546 : !
1547 : case(3) ! complexes converged
1548 2 : write(output_unit, *) ''
1549 2 : write(output_unit, '(A50,A20,F8.4,A34)') '*** OPTIMIZATION TERMINATED BECAUSE THE POPULATION', &
1550 4 : ' HAS CONVERGED INTO ', gnrng * 100._dp, ' PERCENT OF THE FEASIBLE SPACE ***'
1551 : !
1552 : case default
1553 0 : write(error_unit, *) 'This termination case is not implemented!'
1554 0 : stop
1555 : end select
1556 :
1557 4 : write(output_unit, *) ''
1558 4 : write(output_unit, '(A66)') '*** PRINT THE FINAL PARAMETER ESTIMATE AND ITS CRITERION VALUE ***'
1559 :
1560 50 : end subroutine write_termination_case
1561 :
1562 4 : subroutine set_optional()
1563 :
1564 : implicit none
1565 :
1566 4 : if (present(neval)) neval = icall
1567 4 : if (present(bestf) .and. .not. maxit) bestf = bestf_tmp
1568 4 : if (present(bestf) .and. maxit) bestf = -bestf_tmp
1569 4 : if (present(history)) then
1570 12 : allocate(history(icall))
1571 9658 : history(:) = history_tmp(1 : icall)
1572 : end if
1573 :
1574 4 : end subroutine set_optional
1575 :
1576 : end function sce
1577 :
1578 :
1579 100 : 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 1172 : real(dp) :: xsum1 ! sum of all values per parameter
1602 1172 : real(dp) :: xsum2 ! sum of square of all values per parameter
1603 50 : real(dp) :: gsum ! sum of all (range scaled) currently covered parameter ranges:
1604 : ! ! (max-min)/bound
1605 1172 : real(dp), dimension(size(bound, 1)) :: xmin ! minimal value per parameter
1606 1222 : real(dp), dimension(size(bound, 1)) :: xmax ! maximal value per parameter
1607 1222 : real(dp), dimension(size(bound, 1)) :: xmean ! mean value per parameter
1608 : real(dp), parameter :: delta = tiny(1.0_dp)
1609 : !
1610 50 : nn = size(x, 2)
1611 50 : npt1 = size(x, 1)
1612 :
1613 : ! compute maximum, minimum and standard deviation of parameter values
1614 50 : gsum = 0._dp
1615 1172 : do kk = 1, nn
1616 1172 : if (mask(kk)) then
1617 1122 : xmax(kk) = -huge(1.0_dp)
1618 1122 : xmin(kk) = huge(1.0_dp)
1619 1122 : xsum1 = 0._dp
1620 1122 : xsum2 = 0._dp
1621 133470 : do ii = 1, npt1
1622 132348 : xmax(kk) = dmax1(x(ii, kk), xmax(kk))
1623 132348 : xmin(kk) = dmin1(x(ii, kk), xmin(kk))
1624 132348 : xsum1 = xsum1 + x(ii, kk)
1625 133470 : xsum2 = xsum2 + x(ii, kk) * x(ii, kk)
1626 : end do
1627 1122 : xmean(kk) = xsum1 / real(npt1, dp)
1628 1122 : xnstd(kk) = (xsum2 / real(npt1, dp) - xmean(kk) * xmean(kk))
1629 1122 : if (xnstd(kk) .le. delta) then
1630 0 : xnstd(kk) = delta
1631 : end if
1632 1122 : xnstd(kk) = sqrt(xnstd(kk))
1633 1122 : xnstd(kk) = xnstd(kk) / bound(kk)
1634 1122 : gsum = gsum + log(delta + (xmax(kk) - xmin(kk)) / bound(kk))
1635 : end if
1636 : end do
1637 50 : gnrng = exp(gsum / real(nn, dp))
1638 : !
1639 : ! check if normalized standard deviation of parameter is <= eps
1640 50 : ipcnvg = 0_i4
1641 50 : if (gnrng .le. peps) then
1642 2 : ipcnvg = 1_i4
1643 : end if
1644 :
1645 50 : end subroutine parstt
1646 :
1647 0 : 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 50 : 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 0 : n = size(a, dim = 2)
1678 : ! NEW number of complexes
1679 0 : ngs1 = ngs2 - 1
1680 : ! number of points in NEW population
1681 0 : npt = ngs1 * npg
1682 :
1683 0 : do igs = 1, ngs1
1684 0 : do ipg = 1, npg
1685 0 : k1 = (ipg - 1) * ngs2 + igs
1686 0 : k2 = (ipg - 1) * ngs1 + igs
1687 0 : do i = 1, n
1688 0 : b(k2, i) = a(k1, i)
1689 : end do
1690 0 : bf(k2) = af(k1)
1691 : end do
1692 : end do
1693 : !
1694 0 : do j = 1, npt
1695 0 : do i = 1, n
1696 0 : a(j, i) = b(j, i)
1697 : end do
1698 0 : af(j) = bf(j)
1699 : end do
1700 : !
1701 0 : end subroutine comp
1702 :
1703 9004 : subroutine sort_matrix(rb, ra)
1704 : !
1705 0 : 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 4502 : integer(i4), dimension(size(rb, 1)) :: iwk
1718 : !
1719 4502 : n = size(rb, 1)
1720 4502 : m = size(rb, 2)
1721 :
1722 : ! indexes of sorted reference vector
1723 4502 : iwk(:) = sort_index(ra(1 : n))
1724 :
1725 : ! sort reference vector
1726 541668 : ra(1 : n) = ra(iwk)
1727 :
1728 : ! sort each column of second array
1729 134270 : do i = 1, m
1730 16164578 : rb(1 : n, i) = rb(iwk, i)
1731 : end do
1732 :
1733 4502 : end subroutine sort_matrix
1734 :
1735 22806 : 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 4502 : 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 11403 : ibound = 0_i4
1762 11403 : nn = size(x, 1)
1763 :
1764 145177 : do ii = 1, nn
1765 145177 : if (mask(ii)) then
1766 134186 : if ((x(ii) .lt. bl(ii)) .or. (x(ii) .gt. bu(ii))) then
1767 : ! This constraint is violated
1768 412 : ibound = 1_i4
1769 : return
1770 : end if
1771 : end if
1772 : end do
1773 :
1774 11403 : end subroutine chkcst
1775 :
1776 678 : 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 11403 : 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 339 : real(dp) :: rand
1804 : !
1805 339 : nn = size(xi, 1)
1806 6999 : do jj = 1, nn
1807 6999 : if (mask(jj)) then
1808 6660 : ibound = 1
1809 13611 : do while (ibound .eq. 1)
1810 6951 : if (idist .eq. 1) then
1811 3669 : call xor4096(0_i8, rand, save_state = save_state)
1812 3669 : x(jj) = bl(jj) + std(jj) * rand * (bu(jj) - bl(jj))
1813 : end if
1814 6951 : if (idist .eq. 2) then
1815 3282 : call xor4096g(0_i8, rand, save_state = save_state)
1816 3282 : x(jj) = xi(jj) + std(jj) * rand * (bu(jj) - bl(jj))
1817 : end if
1818 : !
1819 : ! Check explicit constraints
1820 : !
1821 41415 : call chkcst((/x(jj)/), (/bl(jj)/), (/bu(jj)/), (/mask(jj)/), ibound)
1822 : end do
1823 : else
1824 0 : x(jj) = xi(jj)
1825 : end if
1826 : end do
1827 :
1828 339 : end subroutine getpnt
1829 :
1830 4452 : subroutine cce(s, sf, bl, bu, maskpara, xnstd, icall, maxn, maxit, save_state_gauss, functn, eval, &
1831 4452 : 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 339 : 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 133098 : real(dp), dimension(size(s, 2)) :: sw ! the worst point of the simplex
1885 133098 : real(dp), dimension(size(s, 2)) :: sb ! the best point of the simplex
1886 137550 : real(dp), dimension(size(s, 2)) :: ce ! the centroid of the simplex excluding wo
1887 137550 : real(dp), dimension(size(s, 2)) :: snew ! new point generated from the simplex
1888 4452 : real(dp) :: fw ! function value of the worst point
1889 4452 : 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 4452 : nps = size(s, 1)
1897 4452 : nn = size(s, 2)
1898 : !
1899 : ! equivalence of variables for readabilty of code
1900 4452 : n = nps
1901 4452 : 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 133098 : do j = 1, m
1908 128646 : sb(j) = s(1, j)
1909 128646 : sw(j) = s(n, j)
1910 128646 : ce(j) = 0.0_dp
1911 3973284 : do i = 1, n - 1
1912 3973284 : ce(j) = ce(j) + s(i, j)
1913 : end do
1914 133098 : ce(j) = ce(j) / real(n - 1, dp)
1915 : end do
1916 4452 : fw = sf(n)
1917 : !
1918 : ! compute the new point snew
1919 : !
1920 : ! first try a reflection step
1921 133098 : do j = 1, m
1922 133098 : if (maskpara(j)) then
1923 128646 : snew(j) = ce(j) + alpha * (ce(j) - sw(j))
1924 : else
1925 0 : snew(j) = s(1, j)
1926 : end if
1927 : end do
1928 : !
1929 : ! check if snew satisfies all constraints
1930 4452 : 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 4452 : if (ibound .eq. 1) then
1937 121 : 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 4452 : 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 4452 : if (.not. maxit) then
1950 4452 : fnew = functn(snew, eval)
1951 4452 : icall = icall + 1_i8
1952 4452 : history(icall) = min(history(icall - 1), fnew)
1953 : else
1954 0 : fnew = -functn(snew, eval)
1955 0 : icall = icall + 1_i8
1956 0 : history(icall) = max(history(icall - 1), -fnew)
1957 : end if
1958 : !
1959 : ! maximum numbers of function evaluations reached
1960 4452 : 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 4452 : if ((fnew .gt. fw) .or. (.not. is_finite(fnew))) then
1965 6417 : do j = 1, m
1966 6417 : if (maskpara(j)) then
1967 6102 : snew(j) = ce(j) - beta * (ce(j) - sw(j))
1968 : else
1969 0 : snew(j) = s(1, j)
1970 : end if
1971 : end do
1972 : !
1973 : ! compute the function value of the contracted point
1974 315 : 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 315 : if (.not. maxit) then
1983 315 : fnew = functn(snew, eval)
1984 315 : icall = icall + 1_i8
1985 315 : history(icall) = min(history(icall - 1), fnew)
1986 : else
1987 0 : fnew = -functn(snew, eval)
1988 0 : icall = icall + 1_i8
1989 0 : history(icall) = max(history(icall - 1), -fnew)
1990 : end if
1991 : !
1992 : ! maximum numbers of function evaluations reached
1993 315 : 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 315 : 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 84 : 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 84 : 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 84 : if (.not. maxit) then
2014 84 : fnew = functn(snew, eval)
2015 84 : icall = icall + 1_i8
2016 84 : history(icall) = min(history(icall - 1), fnew)
2017 : else
2018 0 : fnew = -functn(snew, eval)
2019 0 : icall = icall + 1_i8
2020 0 : history(icall) = max(history(icall - 1), -fnew)
2021 : end if
2022 : !
2023 : ! maximum numbers of function evaluations reached
2024 84 : if (icall .ge. maxn) return
2025 : !
2026 : ! successful mutation
2027 : ! replace the worst point by the new point
2028 336 : do j = 1, m
2029 336 : s(n, j) = snew(j)
2030 : end do
2031 84 : sf(n) = fnew
2032 : end if
2033 : end if
2034 : !
2035 : ! replace the worst point by the new point
2036 133098 : do j = 1, m
2037 133098 : s(n, j) = snew(j)
2038 : end do
2039 4452 : sf(n) = fnew
2040 :
2041 4452 : end subroutine cce
2042 :
2043 : END MODULE mo_sce
|
