Line data Source code
1 : !> \file mo_anneal.f90
2 : !> \brief \copybrief mo_anneal
3 : !> \details \copydetails mo_anneal
4 :
5 : !> \brief Anneal optimization of cost function.
6 : !> \details Minimization of cost function and temperature finding of minima.
7 : !> \changelog
8 : !! - Juliane Mai, Mar 2012
9 : !! - module implementation
10 : !! - Juliane Mai, May 2012
11 : !! - anneal: sp version
12 : !! - Juliane Mai, May 2012
13 : !! - anneal: documentation
14 : !! - Juliane Mai, May 2012
15 : !! - GetTemperature: sp and dp version
16 : !! - Juliane Mai, Jun 2012
17 : !! - weighted parameter selection
18 : !! - Juliane Mai, Aug 2012
19 : !! - function anneal instead of subroutine
20 : !! - using new module get_timeseed as default for seeding
21 : !! - new optional for minimization or maximization
22 : !! - fixed parameter ranges possible instead of interface range
23 : !! - Juliane Mai, Nov 2012
24 : !! - history of achieved objective function values as optional out only in anneal but not anneal_valid
25 : !! - Juliane Mai, Jan 2013
26 : !! - including DDS features in anneal, i.e. reflection at parameter boundaries,
27 : !! different parameter selection modes (one, all, neighborhood), and
28 : !! different parameter pertubation modes (flexible r=dR (anneal version) or
29 : !! constant r=0.2 (dds version))
30 : !! - remove sp versions
31 : !! - fixed and flexible parameter ranges are now in one function
32 : !! using optional arguments
33 : !! - undef_funcval instead of anneal_valid function
34 : !! - Juliane Mai, Feb 2013
35 : !! - xor4096 optionals combined in save_state
36 : !! - Arya Prasetya, Dec 2021
37 : !! - doxygen documentation anneal and get_temperature
38 : !> \author Juliane Mai
39 : !> \date Mar 2012
40 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
41 : !! FORCES is released under the LGPLv3+ license \license_note
42 : MODULE mo_anneal
43 :
44 : USE mo_kind, ONLY : i4, i8, dp
45 : USE mo_utils, ONLY : le, ge
46 : USE mo_xor4096, ONLY : get_timeseed, xor4096, xor4096g, n_save_state
47 : USE mo_optimization_utils, ONLY : eval_interface, objective_interface
48 :
49 : IMPLICIT NONE
50 :
51 : PUBLIC :: anneal ! Minimization of a cost function via Simaulated Annealing
52 : PUBLIC :: GetTemperature ! Function which returns the optimal initial temperature for
53 : ! ! a given acceptance ratio and initial parameter set
54 :
55 : ! ------------------------------------------------------------------
56 :
57 : !> \brief Optimize cost function with simulated annealing.
58 :
59 : !> \details
60 : !! Optimizes a user provided cost function using the Simulated Annealing strategy.
61 : !!
62 : !! \b Example
63 : !!
64 : !! User defined function `cost_dp` which calculates the cost function value for a
65 : !! parameter set (the interface given below has to be used for this function!).
66 : !!
67 : !! \code{.f90}
68 : !! para = (/ 1.0_dp , 2.0_dp /)
69 : !! prange(1,:) = (/ 0.0_dp, 10.0_dp /)
70 : !! prange(2,:) = (/ 0.0_dp, 50.0_dp /)
71 : !!
72 : !! parabest = anneal(cost_dp, para, prange)
73 : !! \endcode
74 : !!
75 : !! See also test folder for a detailed example, "test/test_mo_anneal".
76 : !!
77 : !! \b Literature
78 : !! 1. S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi.
79 : !! _Optimization by simulated annealing_.
80 : !! Science, 220:671-680, 1983.
81 : !! 2. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller.
82 : !! _Equation of state calculations by fast computing machines_.
83 : !! J. Chem. Phys., 21:1087-1092, June 1953.
84 : !! 3. B. A. Tolson, and C. A. Shoemaker.
85 : !! _Dynamically dimensioned search algorithm for computationally efficient watershed model calibration_.
86 : !! WRR, 43(1), W01413, 2007.
87 : !!
88 : !> \param[in] "INTERFACE :: eval" Interface calculating the eval function at a given point.
89 : !> \param[in] "INTERFACE :: cost" Interface calculating the cost function at a given point.
90 : !> \param[in] "real(dp), dimension(:) :: para" Initial parameter set.
91 : !> \param[in] "real(dp), dimension(size(para),2), optional :: prange"
92 : !! Lower and upper bound per parameter.
93 : !> \param[in] "interface, optional :: prange_func" Interface calculating the
94 : !! feasible range for a parameter
95 : !! at a certain point, if ranges are variable.
96 : !> \param[in] "real(dp), optional :: temp" Initial temperature. \n
97 : !! DEFAULT: Get_Temperature.
98 : !> \param[in] "real(dp), optional :: DT" Geometrical decreement of temperature. \n
99 : !! 0.7<DT<0.999. \n
100 : !! DEFAULT: 0.9_dp.
101 : !> \param[in] "integer(i4), optional :: nITERmax" Maximal number of iterations. \n
102 : !! Will be increased by 10% if stopping criteria of
103 : !! acceptance ratio or epsilon decreement of cost
104 : !! function is not fullfilled. \n
105 : !! DEFAULT: 1000_i4.
106 : !> \param[in] "integer(i4), optional :: LEN" Length of Markov Chain. \n
107 : !! DEFAULT: MAX(250_i4, size(para,1)).
108 : !> \param[in] "integer(i4), optional :: nST" Number of consecutive LEN steps. \n
109 : !! DEFAULT: 5_i4.
110 : !> \param[in] "real(dp), optional :: eps" Stopping criteria of epsilon decreement of
111 : !! cost function \n
112 : !! DEFAULT: 0.0001_dp
113 : !> \param[in] "real(dp), optional :: acc" Stopping criteria for Acceptance Ratio \n
114 : !! acc <= 0.1_dp \n
115 : !! DEFAULT: 0.1_dp
116 : !> \param[in] "integer(i4/i8), dimension(3), optional :: seeds"
117 : !! Seeds of random numbers used for random parameter
118 : !! set generation \n
119 : !! DEFAULT: dependent on current time
120 : !> \param[in] "logical, optional :: printflag"
121 : !! If .true. detailed command line output is written \n
122 : !! DEFAULT: .false.
123 : !> \param[in] "logical, dimension(size(para)), optional :: maskpara"
124 : !! maskpara(i) = .true. --> parameter is optimized \n
125 : !! maskpara(i) = .false. --> parameter is discarded
126 : !! from optimiztaion \n
127 : !! DEFAULT: .true.
128 : !> \param[in] "real(dp), dimension(size(para)), optional :: weight"
129 : !! vector of weights per parameter: \n
130 : !! gives the frequency of parameter to be chosen for
131 : !! optimization (will be scaled to a CDF internally) \n
132 : !! eg. [1,2,1] --> parameter 2 is chosen twice as
133 : !! often as parameter 1 and 2 \n
134 : !! DEFAULT: weight = 1.0_dp
135 : !> \param[in] "integer(i4), optional :: changeParaMode" which and how many param.are changed in one step \n
136 : !! 1 = one parameter \n
137 : !! 2 = all parameter \n
138 : !! 3 = neighborhood parameter \n
139 : !! DEFAULT: 1_i4
140 : !> \param[in] "logical, optional :: reflectionFlag" if new parameter values are Gaussian
141 : !! distributed and reflected (.true.) or
142 : !! uniform in range (.false.) \n
143 : !! DEFAULT: .false.
144 : !> \param[in] "logical, optional :: pertubFlexFlag" if pertubation of Gaussian distributed
145 : !! parameter values is constant
146 : !! at 0.2 (.false.) or
147 : !! depends on dR (.true.) \n
148 : !! DEFAULT: .true.
149 : !> \param[in] "logical, optional :: maxit" maximizing (.true.) or minimizing (.false.)
150 : !! a function \n
151 : !! DEFAULT: .false. (minimization)
152 : !> \param[in] "real(dp), optional :: undef_funcval" objective function value defining invalid
153 : !! model output, e.g. -999.0_dp \n
154 : !> \param[in] "character(len=*) , optional :: tmp_file" file with temporal output
155 : !> \param[out] "real(dp), optional :: funcbest" minimized value of cost function
156 : !> \param[out] "real(dp), dimension(:,:), allocatable, optional :: history"
157 : !! returns a vector of achieved objective
158 : !! after ith model evaluation
159 : !> \retval "real(dp) :: parabest(size(para))" Parameter set minimizing the cost function.
160 :
161 : !> \note
162 : !! - Either fixed parameter range (`prange`) OR flexible parameter range (function interface `prange_func`)
163 : !! has to be given in calling sequence.
164 : !! - Only double precision version available.
165 : !! - If single precision is needed not only dp has to be replaced by sp
166 : !! but also i8 of `save_state` (random number variables) has to be replaced by i4.
167 : !! - ParaChangeMode > 1 is not applied in GetTemperature.
168 : !! - For Temperature estimation always only one single parameter is changed (ParaChangeMode=1)
169 : !! which should give theoretically always the best estimate.
170 : !! - `cost_func` and `prange_func` are user defined functions. See interface definition.
171 :
172 : !> \author Luis Samaniego
173 : !> \date Jan 2000
174 : !> \date Mar 2003
175 : !! - Re-heating
176 :
177 : !> \author Juliane Mai
178 : !> \date Mar 2012
179 : !! - Modular version
180 : !> \date May 2012
181 : !! - sp version
182 : !! - documentation
183 :
184 : ! ------------------------------------------------------------------
185 :
186 : INTERFACE anneal
187 : MODULE PROCEDURE anneal_dp
188 : END INTERFACE anneal
189 :
190 : ! ------------------------------------------------------------------
191 :
192 : !> \brief Find initial temperature for simulated annealing.
193 :
194 : !> \details Determines an initial temperature for Simulated Annealing achieving
195 : !! certain acceptance ratio.
196 : !!
197 : !! \b Example
198 : !!
199 : !! User defined function 'cost_dp' which calculates the cost function value for a
200 : !! parameter set (the interface given below has to be used for this function!).
201 : !!
202 : !! \code{.f90}
203 : !! para = (/ 1.0_dp , 2.0_dp /)
204 : !! acc_goal = 0.95_dp
205 : !! prange(1,:) = (/ 0.0_dp, 10.0_dp /)
206 : !! prange(2,:) = (/ 0.0_dp, 50.0_dp /)
207 : !!
208 : !! temp = GetTemperature(para, cost_dp, acc_goal, prange)
209 : !! \endcode
210 : !!
211 : !! See also test folder for a detailed example, "pf_tests/test_mo_anneal".
212 : !!
213 : !! \b Literature
214 : !! 1. Walid Ben-Ameur.
215 : !! _Compututing the Initial Temperature of Simulated Annealing_.
216 : !! Comput. Opt. and App. (2004).
217 : !!
218 : !> \param[in] "real(dp), dimension(:) :: paraset" Initial (valid) parameter set.
219 : !> \param[in] "INTERFACE :: eval" Interface calculating the eval function at a given point.
220 : !> \param[in] "INTERFACE :: cost" Interface calculating the cost function at a given point.
221 : !> \param[in] "INTERFACE :: prange_func" Interface for functional ranges.
222 : !> \param[in] "real(dp) :: acc_goal" Acceptance Ratio which has to be achieved.
223 : !> \param[in] "real(dp), dimension(size(para),2), optional :: prange"
224 : !! Lower and upper bound per parameter.
225 : !> \param[in] "INTERFACE, optional :: prange_func" Interface calculating the
226 : !! feasible range for a parameter
227 : !! at a certain point, if ranges are variable.
228 : !> \param[in] "integer(i4), optional :: samplesize" Number of iterations the estimation of temperature
229 : !! is based on. \n
230 : !! DEFAULT: Max(20_i4*n,250_i4)
231 : !> \param[in] "logical, dimension(size(para)), optional :: maskpara"
232 : !! maskpara(i) = .true. --> parameter is optimized. \n
233 : !! maskpara(i) = .false. --> parameter is discarded
234 : !! from optimiztaion. \n
235 : !! DEFAULT: .true. .
236 : !> \param[in] "INTEGER(I4/I8), dimension(2), optional :: seeds"
237 : !! Seeds of random numbers used for random parameter
238 : !! set generation. \n
239 : !! DEFAULT: dependent on current time.
240 : !> \param[in] "logical, optional :: printflag"
241 : !! If .true. detailed command line output is written.
242 : !! DEFAULT: .false. .
243 : !> \param[in] "real(dp), dimension(size(para,1)), optional :: weight"
244 : !! Vector of weights per parameter. \n
245 : !! Gives the frequency of parameter to be chosen for
246 : !! optimization (will be scaled to a CDF internally). \n
247 : !! eg. [1,2,1] --> parameter 2 is chosen twice as
248 : !! often as parameter 1 and 2. \n
249 : !! DEFAULT: weight = 1.0_dp.
250 : !> \param[in] "logical, optional :: maxit" Minimizing (.false.) or maximizing (.true.)
251 : !! a function. \n
252 : !! DEFAULT: .false. (minimization).
253 : !> \param[in] "real(dp), optional :: undef_funcval" Objective function value defining invalid
254 : !! model output, e.g. -999.0_dp.
255 :
256 : !> \retval "real(dp) :: temperature" Temperature achieving a certain
257 : !! acceptance ratio in Simulated Annealing.
258 :
259 : !> \note
260 : !! - Either fixed parameter range (`prange`) OR flexible parameter range (function interface `prange_func`)
261 : !! has to be given in calling sequence.
262 : !! - Only double precision version available.
263 : !! - If single precision is needed not only dp has to be replaced by sp
264 : !! but also i8 of `save_state` (random number variables) has to be replaced by i4.
265 : !! - ParaChangeMode > 1 is not applied in GetTemperature.
266 : !! - For Temperature estimation always only one single parameter is changed (ParaChangeMode=1)
267 : !! which should give theoretically always the best estimate.
268 : !! - `cost_dp` and `prange_func` are user defined functions. See interface definition.
269 :
270 : !> \author Juliane Mai
271 : !> \date May 2012
272 :
273 : ! ------------------------------------------------------------------
274 :
275 : INTERFACE GetTemperature
276 : MODULE PROCEDURE GetTemperature_dp
277 : END INTERFACE GetTemperature
278 :
279 : PRIVATE
280 :
281 : INTERFACE Generate_Neighborhood_weight
282 : MODULE PROCEDURE Generate_Neighborhood_weight_dp
283 : END INTERFACE Generate_Neighborhood_weight
284 :
285 : ! ------------------------------------------------------------------
286 :
287 : CONTAINS
288 :
289 4 : FUNCTION anneal_dp(eval, cost, para, & ! obligatory
290 0 : prange, prange_func, & ! optional IN: exactly one of both
291 : temp, Dt, nITERmax, Len, nST, eps, acc, & ! optional IN
292 8 : seeds, printflag, maskpara, weight, & ! optional IN
293 : changeParaMode, reflectionFlag, & ! optional IN
294 : pertubFlexFlag, maxit, undef_funcval, & ! optional IN
295 : tmp_file, & ! optional IN
296 : funcbest, history & ! optional OUT
297 : ) &
298 20 : result(parabest)
299 :
300 : IMPLICIT NONE
301 :
302 : procedure(eval_interface), INTENT(IN), POINTER :: eval
303 : procedure(objective_interface), intent(in), pointer :: cost
304 :
305 : INTERFACE
306 : SUBROUTINE prange_func(paraset, iPar, rangePar)
307 : ! gives the range (min,max) of the parameter iPar at a certain parameter set paraset
308 : use mo_kind
309 : real(dp), dimension(:), INTENT(IN) :: paraset
310 : integer(i4), INTENT(IN) :: iPar
311 : real(dp), dimension(2), INTENT(OUT) :: rangePar
312 : END SUBROUTINE prange_func
313 : END INTERFACE
314 : optional :: prange_func
315 :
316 : real(dp), dimension(:), intent(in) :: para !< initial parameter
317 :
318 : real(dp), optional, dimension(size(para, 1), 2), intent(in) :: prange !< lower and upper limit per parameter
319 : real(dp), optional, intent(in) :: temp !< starting temperature (DEFAULT: Get_Temperature)
320 : real(dp), optional, intent(in) :: Dt !< geometrical decreement, 0.7<DT<0.999 (DEFAULT: 0.9)
321 : integer(i4), optional, intent(in) :: nITERmax !< maximal number of iterations (DEFAULT: 1000)
322 : integer(i4), optional, intent(in) :: Len !< Length of Markov Chain,
323 : !! DEFAULT: max(250, size(para,1))
324 : integer(i4), optional, intent(in) :: nST !< Number of consecutive LEN steps! (DEFAULT: 5)
325 : real(dp), optional, intent(in) :: eps !< epsilon decreement of cost function (DEFAULT: 0.01)
326 : real(dp), optional, intent(in) :: acc !< Acceptance Ratio, <0.1 stopping criteria
327 : !! (DEFAULT: 0.1)
328 : INTEGER(I8), optional, intent(in) :: seeds(3) !< Seeds of random numbers (DEFAULT: Get_timeseed)
329 : logical, optional, intent(in) :: printflag !< If command line output is written (.true.)
330 : !! (DEFAULT: .false.)
331 : logical, optional, dimension(size(para, 1)), intent(in) :: maskpara !< true if parameter will be optimized
332 : !! false if parameter is discarded in optimization
333 : !! (DEFAULT: .true.)
334 : real(dp), optional, dimension(size(para, 1)), intent(in) :: weight !< vector of weights per parameter
335 : !! gives the frequency of parameter to be
336 : !! chosen for optimization (DEFAULT: uniform)
337 : integer(i4), optional, intent(in) :: changeParaMode !< which and how many param. are changed in a step
338 : !! 1 = one parameter 2 = all parameter
339 : !! 3 = neighborhood parameter (DEFAULT: 1_i4)
340 : logical, optional, intent(in) :: reflectionFlag !< if new parameter values are selected normal
341 : !! distributed and reflected (.true.) or
342 : !! uniform in range (.false.) (DEFAULT: .false.)
343 : logical, optional, intent(in) :: pertubFlexFlag !< if pertubation of normal distributed parameter
344 : !! values is constant 0.2 (.false.) or
345 : !! depends on dR (.true.) (DEFAULT: .true.)
346 : logical, optional, intent(in) :: maxit !< Maximization or minimization of function
347 : !! maximization = .true., minimization = .false.
348 : !! (DEFAULT: .false.)
349 : real(dp), optional, intent(in) :: undef_funcval !< objective function value occuring if
350 : !! parameter set leads to invalid model results,
351 : !! e.g. -9999.0_dp! (DEFAULT: not present)
352 : CHARACTER(LEN = *), optional, intent(in) :: tmp_file !< file for temporal output
353 : real(dp), optional, intent(out) :: funcbest !< minimized value of cost function
354 : !! (DEFAULT: not present)
355 : real(dp), optional, dimension(:, :), allocatable, intent(out) :: history !< returns a vector of achieved objective!
356 : !! after ith model evaluation (DEFAULT: not present)
357 :
358 : real(dp), dimension(size(para, 1)) :: parabest !< parameter set minimizing objective
359 :
360 : integer(i4) :: n ! Number of parameters
361 : integer(i4) :: iPar, par ! counter for parameter
362 12 : real(dp), dimension(2) :: iParRange ! parameter's range
363 4 : real(dp) :: T_in ! Temperature
364 4 : real(dp) :: DT_IN ! Temperature decreement
365 : integer(i4) :: nITERmax_in ! maximal number of iterations
366 : integer(i4) :: LEN_IN ! Length of Markov Chain
367 : integer(i4) :: nST_in ! Number of consecutive LEN steps
368 4 : real(dp) :: eps_in ! epsilon decreement of cost function
369 4 : real(dp) :: acc_in ! Acceptance Ratio, <0.1 stopping criteria
370 : INTEGER(I8) :: seeds_in(3) ! Seeds of random numbers
371 : logical :: printflag_in ! If command line output is written
372 : logical :: coststatus ! Checks status of cost function value,
373 : ! ! i.e. is parameter set is feasible
374 0 : logical, dimension(size(para, 1)) :: maskpara_in ! true if parameter will be optimized
375 4 : integer(i4), dimension(:), allocatable :: truepara ! indexes of parameters to be optimized
376 24 : real(dp), dimension(size(para, 1)) :: weight_in ! CDF of parameter chosen for optimization
377 20 : real(dp), dimension(size(para, 1)) :: weightUni ! uniform CDF of parameter chosen for optimization
378 24 : real(dp), dimension(size(para, 1)) :: weightGrad ! linear combination of weight and weightUni
379 : integer(i4) :: changeParaMode_inin ! 1=one, 2=all, 3=neighborhood
380 : logical :: reflectionFlag_inin ! true=gaussian distributed and reflected,
381 : ! ! false=uniform distributed in range
382 : logical :: pertubFlexFlag_inin ! variance of normal distributed parameter values
383 : ! ! .true. = depends on dR, .false. = constant 0.2
384 4 : real(dp) :: maxit_in ! maximization = -1._dp, minimization = 1._dp
385 4 : real(dp) :: costbest ! minimized value of cost function
386 4 : real(dp), dimension(:, :), allocatable :: history_out ! best cost function value after k iterations
387 4 : real(dp), dimension(:, :), allocatable :: history_out_tmp ! helper vector
388 :
389 : type paramLim
390 : real(dp) :: min ! minimum value
391 : real(dp) :: max ! maximum value
392 : real(dp) :: new ! new state value
393 : real(dp) :: old ! old state value
394 : real(dp) :: best ! best value found
395 : real(dp) :: dMult ! sensitivity multiplier for parameter search
396 : end type paramLim
397 4 : type (paramLim), dimension (size(para, 1)), target :: gamma ! Parameter
398 :
399 : ! for random numbers
400 4 : real(dp) :: RN1, RN2, RN3 ! Random numbers
401 : integer(I8), dimension(n_save_state) :: save_state_1 ! optional arguments for restarting RN stream 1
402 : integer(I8), dimension(n_save_state) :: save_state_2 ! optional arguments for restarting RN stream 2
403 : integer(I8), dimension(n_save_state) :: save_state_3 ! optional arguments for restarting RN stream 3
404 : ! for dds parameter selection
405 4 : logical, dimension(size(para, 1)) :: neighborhood ! selected parameter in neighborhood
406 4 : real(dp) :: pertubationR ! neighborhood pertubation size parameter
407 :
408 : ! for SA
409 : integer(i4) :: idummy, i
410 4 : real(dp) :: NormPhi
411 8 : real(dp) :: ac_ratio, pa
412 4 : integer(i4), dimension(:, :), allocatable :: iPos_iNeg_history ! stores iPos and iNeg of nST last Markov Chains
413 4 : real(dp) :: fo, fn, df, fBest
414 8 : real(dp) :: rho, fInc, fbb, dr
415 8 : real(dp) :: T0, DT0
416 : real(dp), parameter :: small = -700._DP
417 : integer(i4) :: j, iter, kk
418 : integer(i4) :: Ipos, Ineg
419 : integer(i4) :: iConL, iConR, iConF
420 : integer(i4) :: iTotalCounter ! includes reheating for final conditions
421 : integer(i4) :: iTotalCounterR ! counter of interations in one reheating
422 : logical :: iStop
423 : logical :: ldummy
424 :
425 : ! CHECKING OPTIONALS
426 :
427 4 : n = size(para, 1)
428 :
429 : ! either range or rangfunc has to be given
430 4 : if (present(prange) .eqv. present(prange_func)) then
431 0 : stop 'anneal: Either range or prange_func has to be given'
432 : end if
433 :
434 4 : if (present(Dt)) then
435 4 : if ((Dt .lt. 0.7_dp) .or. (Dt .gt. 0.999_dp)) then
436 0 : stop 'Input argument DT must lie between 0.7 and 0.999'
437 : else
438 : DT_IN = Dt
439 : end if
440 : else
441 : DT_IN = 0.9_dp
442 : end if
443 :
444 4 : if (present(nITERmax)) then
445 4 : if (nITERmax .lt. 1_I4) then
446 0 : stop 'Input argument nITERmax must be greater 0'
447 : else
448 4 : nITERmax_in = nITERmax
449 : end if
450 : else
451 0 : nITERmax_in = 1000_i4
452 : end if
453 :
454 4 : if (present(Len)) then
455 4 : if (Len .lt. Max(20_i4 * n, 250_i4)) then
456 0 : print*, 'WARNING: Input argument LEN should be greater than Max(250,20*N), N=number of parameters'
457 0 : LEN_IN = Len
458 : else
459 : LEN_IN = Len
460 : end if
461 : else
462 0 : LEN_IN = Max(20_i4 * n, 250_i4)
463 : end if
464 :
465 4 : idummy = nITERmax_in / LEN_IN + 1_i4
466 12 : allocate(history_out(idummy, 2))
467 :
468 4 : if (present(nST)) then
469 4 : if (nST .lt. 0_i4) then
470 0 : stop 'Input argument nST must be greater than 0'
471 : else
472 : nST_in = nST
473 : end if
474 : else
475 : nST_in = 5_i4
476 : end if
477 :
478 12 : allocate(iPos_iNeg_history(nST_in, 2))
479 52 : iPos_iNeg_history = 0_i4
480 :
481 4 : if (present(eps)) then
482 4 : if (le(eps, 0.0_dp)) then
483 0 : stop 'Input argument eps must be greater than 0'
484 : else
485 4 : eps_in = eps
486 : end if
487 : else
488 : eps_in = 0.0001_dp
489 : end if
490 :
491 4 : if (present(acc)) then
492 4 : if (le(acc, 0.0_dp) .or. ge(acc, 1.0_dp)) then
493 0 : stop 'Input argument acc must lie between 0.0 and 1.0'
494 : else
495 4 : acc_in = acc
496 : end if
497 : else
498 : acc_in = 0.1_dp
499 : end if
500 :
501 4 : if (present(seeds)) then
502 4 : seeds_in = seeds
503 : else
504 : ! Seeds depend on actual time
505 0 : call get_timeseed(seeds_in)
506 : end if
507 :
508 4 : if (present(printflag)) then
509 4 : printflag_in = printflag
510 : else
511 0 : printflag_in = .false.
512 : end if
513 :
514 4 : if (present(maskpara)) then
515 20 : if (count(maskpara) .eq. 0_i4) then
516 0 : stop 'Input argument maskpara: At least one element has to be true'
517 : else
518 20 : maskpara_in = maskpara
519 : end if
520 : else
521 0 : maskpara_in = .true.
522 : end if
523 :
524 4 : if (present(maxit)) then
525 4 : ldummy = maxit
526 4 : if (maxit) then
527 : maxit_in = -1._dp
528 : else
529 4 : maxit_in = 1._dp
530 : end if
531 : else
532 0 : ldummy = .false.
533 0 : maxit_in = 1._dp
534 : end if
535 :
536 4 : if (present(temp)) then
537 4 : if ((temp .lt. 0.0_dp)) then
538 0 : stop 'Input argument temp must be greater then zero'
539 : else
540 4 : T_in = temp
541 : end if
542 : else
543 0 : if (present(prange_func)) then
544 0 : if (present(undef_funcval)) then
545 : T_in = GetTemperature(eval, cost, para, 0.95_dp, prange_func = prange_func, &
546 : maskpara = maskpara_in, samplesize = 2_i4 * LEN_IN, &
547 : seeds = seeds_in(1 : 2), printflag = printflag_in, &
548 0 : maxit = ldummy, undef_funcval = undef_funcval)
549 : else
550 : T_in = GetTemperature(eval, cost, para, 0.95_dp, prange_func = prange_func, &
551 : maskpara = maskpara_in, samplesize = 2_i4 * LEN_IN, &
552 : seeds = seeds_in(1 : 2), printflag = printflag_in, &
553 0 : maxit = ldummy)
554 : end if
555 : else
556 0 : if (present(undef_funcval)) then
557 : T_in = GetTemperature(eval, cost, para, 0.95_dp, prange = prange, &
558 : maskpara = maskpara_in, samplesize = 2_i4 * LEN_IN, &
559 : seeds = seeds_in(1 : 2), printflag = printflag_in, &
560 0 : maxit = ldummy, undef_funcval = undef_funcval)
561 : else
562 : T_in = GetTemperature(eval, cost, para, 0.95_dp, prange = prange, &
563 : maskpara = maskpara_in, samplesize = 2_i4 * LEN_IN, &
564 : seeds = seeds_in(1 : 2), printflag = printflag_in, &
565 0 : maxit = ldummy)
566 : end if
567 : end if
568 : end if
569 :
570 4 : if (present(changeParaMode)) then
571 4 : changeParaMode_inin = changeParaMode
572 : else
573 : changeParaMode_inin = 1_i4
574 : end if
575 :
576 4 : if (present(reflectionFlag)) then
577 4 : reflectionFlag_inin = reflectionFlag
578 : else
579 : reflectionFlag_inin = .false.
580 : end if
581 :
582 4 : if (present(pertubFlexFlag)) then
583 4 : pertubFlexFlag_inin = pertubFlexFlag
584 : else
585 : pertubFlexFlag_inin = .true.
586 : end if
587 :
588 : ! Temporal file writing
589 4 : if(present(tmp_file)) then
590 0 : open(unit = 999, file = trim(adjustl(tmp_file)), action = 'write', status = 'unknown')
591 0 : write(999, *) '# settings :: general'
592 0 : write(999, *) '# nIterations iseed'
593 0 : write(999, *) nITERmax_in, seeds_in
594 0 : write(999, *) '# settings :: anneal specific'
595 0 : write(999, *) '# sa_tmp'
596 0 : write(999, *) T_in
597 0 : write(999, *) '# iter bestf (bestx(j),j=1,nopt)'
598 0 : close(999)
599 : end if
600 :
601 : ! INITIALIZATION
602 :
603 28 : allocate (truepara(count(maskpara_in)))
604 4 : idummy = 0_i4
605 20 : do i = 1, N
606 20 : if (maskpara_in(i)) then
607 16 : idummy = idummy + 1_i4
608 16 : truepara(idummy) = i
609 : end if
610 : end do
611 :
612 4 : if (printflag_in) then
613 0 : print*, 'Following parameters will be optimized: ', truepara
614 : end if
615 :
616 20 : weight_in = 0.0_dp
617 20 : weightUni = 0.0_dp
618 4 : if (present(weight)) then
619 56 : where (maskpara_in(:))
620 4 : weight_in(:) = weight(:)
621 : weightUni(:) = 1.0_dp
622 : end where
623 : else
624 0 : where (maskpara_in(:))
625 : weight_in(:) = 1.0_dp
626 : weightUni(:) = 1.0_dp
627 : end where
628 : end if
629 : ! scaling the weights
630 36 : weight_in = weight_in / sum(weight_in)
631 36 : weightUni = weightUni / sum(weightUni)
632 : ! cummulating the weights
633 16 : do i = 2, n
634 12 : weight_in(i) = weight_in(i) + weight_in(i - 1)
635 16 : weightUni(i) = weightUni(i) + weightUni(i - 1)
636 : end do
637 :
638 4 : call xor4096 (seeds_in(1), RN1, save_state = save_state_1)
639 4 : call xor4096 (seeds_in(2), RN2, save_state = save_state_2)
640 4 : call xor4096g(seeds_in(3), RN3, save_state = save_state_3)
641 4 : seeds_in = 0_i8
642 :
643 : ! Start Simulated Annealing routine
644 20 : gamma(:)%dmult = 1.0_dp
645 20 : gamma(:)%new = para(:)
646 20 : gamma(:)%old = para(:)
647 20 : gamma(:)%best = para(:)
648 4 : NormPhi = -9999.9_dp
649 4 : T0 = T_in
650 4 : DT0 = DT_IN
651 :
652 : ! Generate and evaluate the initial solution state
653 20 : fo = cost(gamma(:)%old, eval) * maxit_in
654 4 : if (abs(fo) .lt. tiny(0.0_dp)) fo = 0.0000001_dp * maxit_in
655 :
656 4 : file_write : if (present(tmp_file)) then
657 0 : open(unit = 999, file = trim(adjustl(tmp_file)), action = 'write', position = 'append', recl = (n + 2) * 30)
658 0 : if (.not. ldummy) then
659 0 : write(999, *) '0', fo, gamma(:)%old
660 : else
661 0 : write(999, *) '0', -fo, gamma(:)%old
662 : end if
663 0 : close(999)
664 : end if file_write
665 :
666 : ! initialize counters /var for new SA
667 4 : iConL = 0_i4
668 4 : iConR = 1_i4
669 4 : iConF = 0_i4
670 4 : iTotalCounter = 0_i4
671 4 : iTotalCounterR = 0_i4
672 4 : fn = 0.0_DP
673 4 : ac_ratio = 1.0_DP
674 4 : iStop = .TRUE.
675 4 : iter = 0_i4
676 : ! restoring initial T param.
677 4 : DT_IN = DT0
678 : ! Storing the best solution so far
679 4 : NormPhi = fo * maxit_in
680 4 : fo = fo / NormPhi
681 4 : fBest = 1.0_dp * maxit_in
682 :
683 4 : if (printflag_in) then
684 0 : print '(A15,E15.7,A4,E15.7,A4)', ' start NSe = ', fBest * maxit_in, ' ( ', fBest * normPhi * maxit_in, ' ) '
685 0 : print '(I8, 2I5, 4E15.7)', 1_i4, 0_i4, 0_i4, 1._dp, T_in, fo * normPhi * maxit_in, fBest * normPhi * maxit_in
686 : end if
687 :
688 : ! ****************** Stop Criterium *******************
689 : ! Repeat until the % reduction of nST consecutive Markov
690 : ! chains (LEN) of the objective function (f) <= epsilon
691 : ! ******************************************************
692 5152 : loopTest : do while (iStop)
693 5148 : iter = iter + 1_i4
694 5148 : Ipos = 0_i4
695 5148 : Ineg = 0_i4
696 5148 : fbb = fBest
697 : !LEN_IN = int( real(iter,dp)*sqrt(real(iter,dp)) / nITER + 1.5*real(N,dp),i4)
698 : ! Repeat LEN times with feasible solution
699 5148 : j = 1
700 1292148 : loopLEN : do while (j .le. LEN_IN)
701 :
702 1287000 : iTotalCounterR = iTotalCounterR + 1_i4
703 1287000 : iTotalCounter = iTotalCounter + 1_i4
704 :
705 : ! Generate a random subsequent state and evaluate its objective function
706 : ! (1) Generate new parameter set
707 1287000 : dR = (1.0_DP - real(iTotalCounterR, dp) / real(nITERmax_in, dp))**2.0_DP
708 1287000 : if (.not. (ge(dR, 0.05_DP)) .and. &
709 : (iTotalCounterR <= int(real(nIterMax_in, dp) / 3._dp * 4_dp, i4))) then
710 821168 : dR = 0.05_DP
711 : end if
712 :
713 1287000 : if (pertubFlexFlag_inin) then
714 1287000 : pertubationR = max(dR / 5.0_dp, 0.05_dp)
715 : else
716 0 : pertubationR = 0.2_dp
717 : end if
718 :
719 : ! gradual weights:
720 : ! acc_ratio close to 1 --> weight
721 : ! acc_ratio close to 0 --> weight uniform
722 : ! gradual weighting is linear combination of these two weights
723 6435000 : weightGrad(:) = weightUni(:) + (ac_ratio) * (weight_in(:) - weightUni(:))
724 :
725 1287000 : select case(changeParaMode_inin)
726 : case(1_i4) ! only one parameter is changed
727 : ! (1a) Change one parameter traditionally
728 1287000 : call xor4096(seeds_in(2), RN2, save_state = save_state_2)
729 1287000 : iPar = 1_i4
730 : !
731 3216686 : do while (weightGrad(iPar) .lt. RN2)
732 1929686 : iPar = iPar + 1_i4
733 : end do
734 1287000 : if (present(prange_func)) then
735 6435000 : call prange_func(gamma(:)%old, iPar, iParRange)
736 : else
737 0 : iParRange(1) = prange(iPar, 1)
738 0 : iParRange(2) = prange(iPar, 2)
739 : end if
740 1287000 : gamma(iPar)%min = iParRange(1)
741 1287000 : gamma(iPar)%max = iParRange(2)
742 1287000 : if (reflectionFlag_inin) then
743 0 : call xor4096g(seeds_in(3), RN3, save_state = save_state_3)
744 0 : gamma(iPar)%new = parGen_dds_dp(gamma(iPar)%old, pertubationR, &
745 0 : gamma(iPar)%min, gamma(iPar)%max, RN3)
746 : else
747 1287000 : call xor4096(seeds_in(2), RN2, save_state = save_state_2)
748 1287000 : gamma(iPar)%new = parGen_anneal_dp(gamma(iPar)%old, dR, &
749 2574000 : gamma(iPar)%min, gamma(iPar)%max, RN2)
750 : end if
751 : case(2_i4) ! all parameter are changed
752 0 : do par = 1, size(truepara)
753 0 : iPar = truepara(par)
754 0 : if (present(prange_func)) then
755 0 : call prange_func(gamma(:)%old, iPar, iParRange)
756 : else
757 0 : iParRange(1) = prange(iPar, 1)
758 0 : iParRange(2) = prange(iPar, 2)
759 : end if
760 0 : gamma(iPar)%min = iParRange(1)
761 0 : gamma(iPar)%max = iParRange(2)
762 0 : if (reflectionFlag_inin) then
763 0 : call xor4096g(seeds_in(3), RN3, save_state = save_state_3)
764 0 : gamma(iPar)%new = parGen_dds_dp(gamma(iPar)%old, pertubationR, &
765 0 : gamma(iPar)%min, gamma(iPar)%max, RN3)
766 : else
767 0 : call xor4096(seeds_in(2), RN2, save_state = save_state_2)
768 0 : gamma(iPar)%new = parGen_anneal_dp(gamma(iPar)%old, dR, &
769 0 : gamma(iPar)%min, gamma(iPar)%max, RN2)
770 : end if
771 : end do
772 : case(3_i4) ! parameter in neighborhood are changed
773 : ! Generate new neighborhood
774 : call generate_neighborhood_weight_dp(truepara, weightGrad, save_state_1, &
775 0 : iTotalCounter, nIterMax_in, neighborhood)
776 : !
777 : ! change parameter in neighborhood
778 1287000 : do iPar = 1, n
779 0 : if (neighborhood(iPar)) then
780 : ! find range of parameter
781 0 : if (present(prange_func)) then
782 0 : call prange_func(gamma(:)%old, iPar, iParRange)
783 : else
784 0 : iParRange(1) = prange(iPar, 1)
785 0 : iParRange(2) = prange(iPar, 2)
786 : end if
787 0 : gamma(iPar)%min = iParRange(1)
788 0 : gamma(iPar)%max = iParRange(2)
789 : !
790 0 : if (reflectionFlag_inin) then
791 : ! generate gaussian distributed new parameter value which is reflected if out of bound
792 0 : call xor4096g(seeds_in(3), RN3, save_state = save_state_3)
793 0 : gamma(iPar)%new = parGen_dds_dp(gamma(iPar)%old, pertubationR, &
794 0 : gamma(iPar)%min, gamma(iPar)%max, RN3)
795 : else
796 : ! generate new parameter value uniform distributed in range (no reflection)
797 0 : call xor4096(seeds_in(2), RN2, save_state = save_state_2)
798 0 : gamma(iPar)%new = parGen_anneal_dp(gamma(iPar)%old, dR, &
799 0 : gamma(iPar)%min, gamma(iPar)%max, RN2)
800 : end if
801 : end if
802 : end do
803 : end select
804 :
805 : ! (2) Calculate new objective function value
806 6435000 : fn = cost(gamma(:)%new, eval) * maxit_in
807 1287000 : coststatus = .true.
808 1287000 : if (present(undef_funcval)) then
809 1287000 : if (abs(fn * maxit_in - undef_funcval) .lt. tiny(1.0_dp)) then
810 : coststatus = .false.
811 : end if
812 : end if
813 :
814 5148 : feasible : if (coststatus) then ! feasible parameter set
815 1287000 : fn = fn / normPhi
816 :
817 : ! Change in cost function value
818 1287000 : df = fn - fo
819 :
820 : ! analyze change in the objective function: df
821 1287000 : if (df < 0.0_DP) then
822 :
823 : ! accept the new state
824 86921 : Ipos = Ipos + 1_i4
825 86921 : fo = fn
826 434605 : gamma(:)%old = gamma(:)%new
827 :
828 : ! keep best solution
829 86921 : if (fo < fBest) then
830 400 : fBest = fo
831 400 : gamma(:)%best = gamma(:)%new
832 : end if
833 : else
834 1200079 : if (df > eps_in) then
835 1197051 : rho = -df / T_in
836 1197051 : if (rho < small) then
837 : pa = 0.0_DP
838 : else
839 1197048 : pa = EXP(rho)
840 : end if
841 : !
842 1197051 : call xor4096(seeds_in(1), RN1, save_state = save_state_1)
843 : !
844 1197051 : if (pa > RN1) then
845 : ! accept new state with certain probability
846 85784 : Ineg = Ineg + 1_i4
847 85784 : fo = fn
848 : ! save old state
849 428920 : gamma(:)%old = gamma(:)%new
850 : end if
851 : end if
852 : end if
853 1287000 : j = j + 1
854 : else
855 : ! function value was not valid
856 0 : iTotalCounterR = iTotalCounterR - 1_i4
857 0 : iTotalCounter = iTotalCounter - 1_i4
858 : end if feasible !valid parameter set
859 : end do loopLEN
860 :
861 : ! estimate acceptance ratio
862 5148 : ac_ratio = (real(Ipos, dp) + real(Ineg, dp)) / real(LEN_IN, dp)
863 5148 : if (modulo(iTotalCounter, LEN_IN * nST_in) / LEN_IN .gt. 0_i4) then
864 4120 : idummy = modulo(iTotalCounter, LEN_IN * nST_in) / LEN_IN
865 : else
866 : idummy = nST_in
867 : end if
868 15444 : iPos_iNeg_history(idummy, :) = (/ iPos, iNeg /)
869 :
870 : ! store current best in history vector
871 5148 : history_out(iTotalCounter / LEN_IN, 1) = real(iTotalCounter, dp)
872 5148 : history_out(iTotalCounter / LEN_IN, 2) = maxit_in * fBest * normPhi
873 :
874 5148 : file_write2 : if (present(tmp_file)) then
875 0 : open(unit = 999, file = trim(adjustl(tmp_file)), action = 'write', position = 'append', recl = (n + 2) * 30)
876 0 : write(999, *) iTotalCounter, maxit_in * fBest * normPhi, gamma(:)%best
877 0 : close(999)
878 : end if file_write2
879 :
880 5148 : if (printflag_in) then
881 0 : print '(I8, 2I5, E15.7, 3E15.7)', ITotalCounter, Ipos, Ineg, ac_ratio, T_in, &
882 0 : fo * NormPhi * maxit_in, fBest * NormPhi * maxit_in
883 : end if
884 :
885 : ! Cooling schedule
886 5148 : if (fbb .gt. tiny(1.0_dp)) then
887 5148 : fInc = (fbb - fBest) / fbb
888 : else
889 : fInc = 1.0_dp
890 : end if
891 5148 : if (fInc < 0.00000001_dp) then
892 5086 : iConF = iConF + 1_i4
893 : else
894 : iConF = 0_i4
895 : end if
896 :
897 5148 : if ((ac_ratio < 0.15_dp) .and. (iConF > 5_i4) .and. (iConR <= -3_i4)) then ! - iConR no reheating
898 : ! Re-heating
899 0 : if (printflag_in) then
900 0 : print *, 'Re-heating: ', iConR
901 : end if
902 0 : iConR = iConR + 1_i4
903 0 : T_in = T0 / 2._DP
904 0 : iter = 0_i4 ! for LEN
905 0 : iTotalCounterR = 0_i4 ! for dR
906 : ! start from current best
907 0 : gamma(:)%old = gamma(:)%best
908 : else
909 : ! Update Temperature (geometrical decrement)
910 5148 : if (ac_ratio < 0.4_dp) then
911 : DT_IN = 0.995_dp
912 : else
913 202 : DT_IN = DT0
914 : end if
915 5148 : T_in = T_in * DT_IN
916 : end if
917 :
918 : ! Stop Criteria: consecutive MC with marginal decrements and acceptance ratio
919 5148 : if (fInc < eps_in) then
920 5086 : iConL = iConL + 1_i4
921 : else
922 : iConL = 0_i4
923 : end if
924 66924 : if ((iConL > nST_in) .and. (ac_ratio < acc_in) .and. (iConR > 2_i4) .and. &
925 : (sum(iPos_iNeg_history(:, :)) .lt. 1_i4)) then
926 0 : iStop = .FALSE.
927 : end if
928 : ! a way out maximum
929 66924 : if ((iTotalCounter > nITERmax_in) .and. & !(ac_ratio > acc_in) .and. &
930 : (sum(iPos_iNeg_history(:, :)) .ge. 1_i4)) then
931 :
932 : ! store achieved history so far
933 128 : allocate(history_out_tmp(size(history_out, 1), size(history_out, 2)))
934 55064 : history_out_tmp = history_out
935 32 : deallocate(history_out)
936 :
937 32 : nITERmax_in = max(nITERmax_in + LEN_IN, int(real(nITERmax_in, dp) * 1.10_dp, i4))
938 32 : if (printflag_in) then
939 0 : print *, 'nITERmax changed to =', nITERmax_in
940 : end if
941 :
942 : ! increase lenght of history vecor
943 32 : idummy = nITERmax_in / LEN_IN + 1_i4
944 96 : allocate(history_out(idummy, 2))
945 :
946 27500 : do j = 1, size(history_out_tmp, 1)
947 82436 : history_out(j, :) = history_out_tmp(j, :)
948 : end do
949 :
950 32 : deallocate(history_out_tmp)
951 :
952 : end if
953 : ! STOP condition
954 5152 : if (iTotalCounter > nITERmax_in) then
955 4 : iStop = .FALSE.
956 : end if
957 :
958 : end do loopTest
959 :
960 : ! calculate cost function again (only for check and return values)
961 20 : parabest = gamma(:)%best
962 4 : costbest = cost(parabest, eval) * maxit_in
963 4 : if (present (funcbest)) then
964 4 : funcbest = costbest * maxit_in
965 : end if
966 4 : fo = costbest / NormPhi
967 :
968 4 : if (printflag_in) then
969 0 : print *, ' '
970 0 : print '(A15,E15.7)', ' end cost = ', maxit_in * fBest * normPhi
971 0 : print *, 'end parameter: '
972 0 : do kk = 1, N
973 0 : print '(A10,I3,A3, E15.7)', ' para #', kk, ' = ', gamma(kk)%best
974 : end do
975 :
976 0 : print *, 'Final check: ', (fo - fBest)
977 : end if
978 :
979 4 : if (present(history)) then
980 16 : allocate(history(size(history_out, 1) - 1, size(history_out, 2)))
981 10300 : history(:, :) = history_out(1 : size(history_out, 1) - 1, :)
982 : end if
983 :
984 4 : deallocate(truepara)
985 4 : deallocate(history_out)
986 :
987 4 : return
988 :
989 20 : END FUNCTION anneal_dp
990 :
991 1 : real(dp) function GetTemperature_dp(eval, cost, paraset, acc_goal, & ! obligatory
992 0 : prange, prange_func, & ! optional IN: exactly one of both
993 0 : samplesize, maskpara, seeds, printflag, & ! optional IN
994 0 : weight, maxit, undef_funcval & ! optional IN
995 : )
996 4 : use mo_kind, only : dp, i4, i8
997 : implicit none
998 :
999 : procedure(eval_interface), INTENT(IN), POINTER :: eval
1000 : procedure(objective_interface), intent(in), pointer :: cost
1001 : real(dp), dimension(:), intent(in) :: paraset !< a valid parameter set of the model
1002 : real(dp), intent(in) :: acc_goal !< acceptance ratio to achieve
1003 : real(dp), optional, dimension(size(paraset, 1), 2), intent(in) :: prange !< lower and upper limit per parameter
1004 : integer(i4), optional, intent(in) :: samplesize !< size of random set the acc_estimate is based on.
1005 : !! DEFAULT: Max(250, 20*Number paras)
1006 : logical, optional, dimension(size(paraset, 1)), intent(in) :: maskpara !< true if parameter will be optimized.
1007 : !! false if parameter is discarded in optimization.
1008 : !! DEFAULT: .true.
1009 : integer(i8), optional, dimension(2), intent(in) :: seeds !< Seeds of random numbers. DEFAULT: time dependent.
1010 : logical, optional, intent(in) :: printflag !< .true. if detailed temperature estimation is
1011 : !! printed. DEFAULT: .false.
1012 : real(dp), OPTIONAL, dimension(size(paraset, 1)), INTENT(IN) :: weight !< vector of weights per parameter gives the
1013 : !! frequency of parameter to be chosen for
1014 : !! optimization. DEFAULT: equal weighting
1015 : logical, OPTIONAL, INTENT(IN) :: maxit !< Maxim. or minim. of function
1016 : !! maximization = .true., minimization = .false. .
1017 : !! DEFAULT: .false.
1018 : real(dp), OPTIONAL, INTENT(IN) :: undef_funcval !< objective function value occuring if
1019 : !! parameter set leads to invalid model results,
1020 : !! e.g. -999.0_dp. DEFAULT: not present
1021 :
1022 : INTERFACE
1023 : SUBROUTINE prange_func(paraset, iPar, rangePar)
1024 : ! gives the range (min,max) of the parameter iPar at a certain parameter set paraset
1025 : use mo_kind
1026 : real(dp), dimension(:), INTENT(IN) :: paraset
1027 : integer(i4), INTENT(IN) :: iPar
1028 : real(dp), dimension(2), INTENT(OUT) :: rangePar
1029 : END SUBROUTINE prange_func
1030 : END INTERFACE
1031 : OPTIONAL :: prange_func
1032 : !
1033 : ! Local variables
1034 :
1035 : integer(i4) :: n
1036 : integer(i4) :: samplesize_in
1037 : integer(i4) :: idummy, i, j
1038 : integer(i4) :: iPar
1039 3 : real(dp), dimension(2) :: iParRange
1040 1 : real(dp) :: NormPhi
1041 1 : real(dp) :: fo, fn, dr
1042 1 : real(dp) :: T
1043 :
1044 : logical :: coststatus ! true if model output is valid
1045 0 : logical, dimension(size(paraset, 1)) :: maskpara_in ! true if parameter will be optimized
1046 1 : integer(i4), dimension(:), ALLOCATABLE :: truepara ! indexes of parameters to be optimized
1047 : logical :: printflag_in ! if detailed estimation of temperature is printed
1048 5 : real(dp), dimension(size(paraset, 1)) :: weight_in ! CDF of parameter to chose for optimization
1049 1 : real(dp) :: maxit_in ! Maximization or minimization of function:
1050 : ! ! -1 = maxim, 1 = minim
1051 :
1052 : type paramLim
1053 : real(dp) :: min ! minimum value
1054 : real(dp) :: max ! maximum value
1055 : real(dp) :: new ! new state value
1056 : real(dp) :: old ! old state value
1057 : real(dp) :: best ! best value found
1058 : real(dp) :: dMult ! sensitivity multiplier
1059 : ! ! for parameter search
1060 : end type paramLim
1061 1 : type (paramLim), dimension (size(paraset, 1)), target :: gamma ! Parameter
1062 :
1063 : ! for random numbers
1064 : INTEGER(I8) :: seeds_in(2)
1065 1 : real(dp) :: RN1, RN2 ! Random numbers
1066 : integer(I8), dimension(n_save_state) :: save_state_1 ! optional arguments for restarting RN stream 1
1067 : integer(I8), dimension(n_save_state) :: save_state_2 ! optional arguments for restarting RN stream 2
1068 :
1069 : ! for initial temperature estimate
1070 1 : real(dp) :: acc_estim ! estimate of acceptance probability
1071 : ! depends on temperature
1072 : ! goal: find a temperature such that acc_estim ~ 0.9
1073 1 : real(dp), dimension(:, :), allocatable :: Energy ! dim(LEN,2) cost function values before (:,1) and
1074 : ! ! after (:,2) transition
1075 :
1076 1 : n = size(paraset, 1)
1077 :
1078 : ! either range or rangfunc has to be given
1079 1 : if (present(prange) .eqv. present(prange_func)) then
1080 0 : stop 'anneal: Either range or prange_func has to be given'
1081 : end if
1082 :
1083 1 : if (present(samplesize)) then
1084 1 : if (samplesize .lt. Max(20_i4 * n, 250_i4)) then
1085 : !stop 'Input argument LEN must be greater than Max(250,20*N), N=number of parameters'
1086 0 : print*, 'WARNING (GetTemperature): '
1087 0 : print*, 'Input argument samplesize should be greater than Max(250,20*N), N=number of parameters'
1088 0 : samplesize_in = samplesize
1089 : else
1090 : samplesize_in = samplesize
1091 : end if
1092 : else
1093 0 : samplesize_in = Max(20_i4 * n, 250_i4)
1094 : end if
1095 :
1096 1 : if (present(maskpara)) then
1097 0 : if (count(maskpara) .eq. 0_i4) then
1098 0 : stop 'Input argument maskpara: At least one element has to be true'
1099 : else
1100 0 : maskpara_in = maskpara
1101 : end if
1102 : else
1103 5 : maskpara_in = .true.
1104 : end if
1105 :
1106 1 : if (present(seeds)) then
1107 1 : seeds_in = seeds
1108 : else
1109 : ! Seeds depend on actual time
1110 0 : call get_timeseed(seeds_in)
1111 0 : print*, 'temp: seeds(1)=', seeds_in(1)
1112 : end if
1113 :
1114 1 : if (present(printflag)) then
1115 1 : printflag_in = printflag
1116 : else
1117 : printflag_in = .false.
1118 : end if
1119 :
1120 1 : if (present(maxit)) then
1121 0 : if (maxit) then
1122 : maxit_in = -1._dp
1123 : else
1124 0 : maxit_in = 1._dp
1125 : end if
1126 : else
1127 : maxit_in = 1._dp
1128 : end if
1129 :
1130 3 : allocate(Energy(samplesize_in, 2))
1131 :
1132 7 : allocate (truepara(count(maskpara_in)))
1133 1 : idummy = 0_i4
1134 5 : do i = 1, N
1135 5 : if (maskpara_in(i)) then
1136 4 : idummy = idummy + 1_i4
1137 4 : truepara(idummy) = i
1138 : end if
1139 : end do
1140 :
1141 5 : weight_in = 0.0_dp
1142 1 : if (present(weight)) then
1143 0 : where (maskpara_in(:))
1144 0 : weight_in(:) = weight(:)
1145 : end where
1146 : else
1147 5 : where (maskpara_in(:))
1148 : weight_in(:) = 1.0_dp
1149 : end where
1150 : end if
1151 : ! scaling the weights
1152 9 : weight_in = weight_in / sum(weight_in)
1153 : ! cummulating the weights
1154 4 : do i = 2, n
1155 4 : weight_in(i) = weight_in(i) + weight_in(i - 1)
1156 : end do
1157 :
1158 : ! Setting up the RNG
1159 : ! (1) Seeds depend on actual time or on input seeds
1160 : ! (2) Initialize the streams
1161 1 : call xor4096(seeds_in(1), RN1, save_state = save_state_1)
1162 1 : call xor4096(seeds_in(2), RN2, save_state = save_state_2)
1163 1 : seeds_in = 0_i8
1164 : ! (3) Now ready for calling
1165 :
1166 5 : gamma(:)%dmult = 1.0_dp
1167 5 : gamma(:)%new = paraset(:)
1168 5 : gamma(:)%old = paraset(:)
1169 5 : gamma(:)%best = paraset(:)
1170 1 : NormPhi = -9999.9_dp
1171 :
1172 1 : fo = cost(paraset, eval) * maxit_in
1173 1 : if (abs(fo) .lt. tiny(0.0_dp)) fo = 0.0000001_dp * maxit_in
1174 1 : NormPhi = fo
1175 1 : fo = fo / NormPhi * maxit_in
1176 :
1177 1 : j = 1
1178 501 : loopSamplesize : do while (j .le. samplesize_in)
1179 : ! Generate a random subsequent state and evaluate its objective function
1180 : ! (1) Generate new parameter set
1181 500 : dR = 1.0_DP
1182 :
1183 : ! (1a) Select parameter to be changed
1184 500 : call xor4096(seeds_in(1), RN1, save_state = save_state_1)
1185 500 : iPar = 1_i4
1186 1250 : do while (weight_in(iPar) .lt. RN1)
1187 750 : iPar = iPar + 1_i4
1188 : end do
1189 :
1190 : ! (1b) Generate new value of selected parameter
1191 500 : if (present(prange_func)) then
1192 2500 : call prange_func(gamma(:)%old, iPar, iParRange)
1193 : else
1194 0 : iParRange(1) = prange(iPar, 1)
1195 0 : iParRange(2) = prange(iPar, 2)
1196 : end if
1197 500 : gamma(iPar)%min = iParRange(1)
1198 500 : gamma(iPar)%max = iParRange(2)
1199 500 : call xor4096(seeds_in(2), RN2, save_state = save_state_2)
1200 500 : gamma(iPar)%new = parGen_anneal_dp(gamma(iPar)%old, gamma(iPar)%dMult * dR, &
1201 1000 : gamma(iPar)%min, gamma(iPar)%max, RN2)
1202 : !
1203 : ! (2) Calculate new objective function value and normalize it
1204 2500 : fn = cost(gamma(:)%new, eval) * maxit_in
1205 500 : coststatus = .true.
1206 500 : if (present(undef_funcval)) then
1207 0 : if (abs(fn * maxit_in - undef_funcval) .lt. tiny(1.0_dp)) then
1208 : coststatus = .false.
1209 : end if
1210 : end if
1211 :
1212 1 : feasible : if (coststatus) then ! feasible parameter set
1213 500 : fn = fn / normPhi
1214 : ! Save Energy states of feasible transitions
1215 : ! for adaption of optimal initial temperature
1216 : ! Walid Ben-Ameur: "Comput. the Initial Temperature of Sim. Annealing"
1217 : ! Comput. Opt. and App. 2004
1218 500 : Energy(j, 2) = fn ! E_max_t
1219 500 : Energy(j, 1) = fo ! E_min_t
1220 500 : j = j + 1
1221 : end if feasible !valid parameter set
1222 : end do loopSamplesize
1223 :
1224 : ! estimation of the acceptance probability based on the random set ||<Samplesize>||
1225 : ! only if actual temperature (T) equals initial temperature (temp)
1226 1003 : T = maxval(Energy)
1227 :
1228 1001 : acc_estim = sum(exp(-(Energy(:, 2) / T))) / sum(exp(-(Energy(:, 1) / T)))
1229 1 : if (printflag_in) then
1230 1 : print*, "acc_estimate = ", acc_estim, " ( T = ", T, " )"
1231 : end if
1232 12 : Do While ((acc_estim .lt. 1.0_dp) .and. (abs(acc_estim - acc_goal) .gt. 0.0001_dp))
1233 11 : T = T * (Log(acc_estim) / Log(acc_goal))**(0.5_dp) ! **(1.0/p) with p=1.0
1234 11023 : if (all(T .gt. Energy(:, 1) / 709._dp) .and. all(T .gt. Energy(:, 2) / 709._dp)) then
1235 11011 : acc_estim = sum(exp(-(Energy(:, 2) / T))) / sum(exp(-(Energy(:, 1) / T)))
1236 11 : if (printflag_in) then
1237 11 : print*, "acc_estimate = ", acc_estim, " ( T = ", T, " )"
1238 : end if
1239 : else
1240 0 : T = T / (Log(acc_estim) / Log(acc_goal))**(0.5_dp)
1241 0 : exit
1242 : end if
1243 : end do
1244 1 : GetTemperature_dp = T
1245 :
1246 3 : end function GetTemperature_dp
1247 :
1248 : !***************************************************
1249 : !* PRIVATE FUNCTIONS *
1250 : !***************************************************
1251 :
1252 1287500 : real(dp) function parGen_anneal_dp(old, dMax, oMin, oMax, RN)
1253 1 : use mo_kind, only : dp, i4
1254 : implicit none
1255 :
1256 : real(dp), intent(IN) :: old, dMax, oMin, oMax
1257 : real(dp), intent(IN) :: RN
1258 :
1259 : ! local variables
1260 1287500 : real(dp) :: oi, ox, delta, dMaxScal
1261 : integer(i4) :: iDigit, iszero ! were intent IN before
1262 : !
1263 1287500 : iszero = 1_i4
1264 1287500 : iDigit = 8_i4
1265 :
1266 1287500 : oi = oMin
1267 1287500 : ox = oMax
1268 : ! scaling (new)
1269 1287500 : dMaxScal = dMax * ABS(oMax - oMin)
1270 :
1271 1287500 : if(ge(oi, ox)) then
1272 0 : parGen_anneal_dp = (oi + ox) / 2.0_DP
1273 : else
1274 1287500 : if ((old - dMaxScal) > oi) oi = old - dMaxScal
1275 1287500 : if ((old + dMaxScal) < ox) ox = old + dMaxScal
1276 1287500 : delta = (oi + RN * (ox - oi)) - old
1277 1287500 : if (delta > dMaxScal) delta = dMaxScal
1278 1287500 : if (delta <-dMaxScal) delta = -dMaxScal
1279 1287500 : parGen_anneal_dp = old + dChange_dp(delta, iDigit, isZero)
1280 : end if
1281 :
1282 : ! Parameter is bounded between Max and Min.
1283 : ! Correction from Kumar and Luis
1284 :
1285 1287500 : if (parGen_anneal_dp < oMin) then
1286 : parGen_anneal_dp = oMin
1287 1287480 : elseif(parGen_anneal_dp > oMax)then
1288 116 : parGen_anneal_dp = oMax
1289 : end if
1290 1287500 : end function parGen_anneal_dp
1291 :
1292 0 : real(dp) function parGen_dds_dp(old, perturb, oMin, oMax, RN)
1293 : ! PURPOSE:
1294 : ! Perturb variable using a standard normal random number RN
1295 : ! and reflecting at variable bounds if necessary
1296 : ! Tolson et al. (2007): STEP 4
1297 : !
1298 : implicit none
1299 :
1300 : real(dp), intent(IN) :: old, perturb, oMin, oMax
1301 : real(dp), intent(IN) :: RN
1302 :
1303 : ! generating new value
1304 0 : parGen_dds_dp = old + perturb * (oMax - oMin) * RN
1305 :
1306 : ! reflect one time and set to boundary value
1307 0 : if (parGen_dds_dp .lt. oMin) then
1308 0 : parGen_dds_dp = oMin + (oMin - parGen_dds_dp)
1309 0 : if (parGen_dds_dp .gt. oMax) parGen_dds_dp = oMin
1310 : end if
1311 :
1312 0 : if (parGen_dds_dp .gt. oMax) then
1313 0 : parGen_dds_dp = oMax - (parGen_dds_dp - oMax)
1314 0 : if (parGen_dds_dp .lt. oMin) parGen_dds_dp = oMax
1315 : end if
1316 :
1317 1287500 : end function parGen_dds_dp
1318 :
1319 1287500 : real(dp) function dChange_dp(delta, iDigit, isZero)
1320 0 : use mo_kind, only : i4, i8, dp
1321 : implicit none
1322 :
1323 : integer(i4), intent(IN) :: iDigit, isZero
1324 : real(dp), intent(IN) :: delta
1325 :
1326 : ! local variables
1327 : integer(i4) :: ioszt
1328 : integer(I8) :: iDelta
1329 :
1330 1287500 : ioszt = 10**iDigit
1331 1287500 : iDelta = int(delta * real(ioszt, dp), i8)
1332 1287500 : if(isZero == 1_i4) then
1333 1287500 : if(iDelta == 0_i8) then
1334 0 : if(delta < 0_i4) then
1335 : iDelta = -1_i8
1336 : else
1337 0 : iDelta = 1_i8
1338 : end if
1339 : end if
1340 : end if
1341 1287500 : dChange_dp = real(iDelta, dp) / real(ioszt, dp)
1342 1287500 : end function dChange_dp
1343 :
1344 0 : subroutine generate_neighborhood_weight_dp(truepara, cum_weight, save_state_xor, iTotalCounter, &
1345 0 : nITERmax, neighborhood)
1346 : ! PURPOSE:
1347 : ! generates a new neighborhood
1348 : ! Tolson et al. (2007): STEP 3
1349 : !
1350 : integer(i4), dimension(:), intent(in) :: truepara
1351 : real(dp), dimension(:), intent(in) :: cum_weight
1352 : integer(i8), dimension(n_save_state), intent(inout) :: save_state_xor
1353 : integer(i4), intent(in) :: iTotalCounter
1354 : integer(i4), intent(in) :: nITERmax
1355 : logical, dimension(size(cum_weight)), intent(out) :: neighborhood
1356 :
1357 : ! local variables
1358 : integer(i4) :: ipar, npar
1359 : integer(i4) :: iSize, size_neighbor
1360 0 : real(dp) :: prob
1361 0 : real(dp) :: rn
1362 0 : real(dp) :: weight, norm
1363 0 : real(dp), dimension(size(cum_weight)) :: local_cum_weight
1364 :
1365 0 : prob = 1.0_dp - Log(real(iTotalCounter, dp)) / Log(real(nITERmax, dp))
1366 0 : neighborhood = .false.
1367 :
1368 0 : npar = size(cum_weight)
1369 0 : local_cum_weight = cum_weight
1370 :
1371 : ! How many parameters will be selected for neighborhood?
1372 0 : size_neighbor = 0_i4
1373 0 : do ipar = 1, size(truepara)
1374 0 : call xor4096(0_i8, rn, save_state = save_state_xor)
1375 0 : if (rn < prob) then
1376 0 : size_neighbor = size_neighbor + 1_i4
1377 : end if
1378 : end do
1379 : ! at least one...
1380 0 : size_neighbor = max(1_i4, size_neighbor)
1381 :
1382 : ! Which parameter will be used for neighborhood?
1383 0 : do iSize = 1, size_neighbor
1384 : ! (1) generate RN
1385 0 : call xor4096(0_i8, rn, save_state = save_state_xor)
1386 : !
1387 : ! (2) find location <iPar> in cummulative distribution function
1388 0 : iPar = 1_i4
1389 0 : do while (local_cum_weight(iPar) .lt. rn)
1390 0 : iPar = iPar + 1_i4
1391 : end do
1392 : !
1393 : ! (3) add parameter to neighborhood
1394 0 : neighborhood(iPar) = .true.
1395 : !
1396 : ! (4) recalculate cummulative distribution function
1397 : ! (4.1) Which weight had iPar?
1398 0 : if (iPar .gt. 1_i4) then
1399 0 : weight = local_cum_weight(iPar) - local_cum_weight(iPar - 1)
1400 : else
1401 0 : weight = local_cum_weight(1)
1402 : end if
1403 : ! (4.2) Substract this weight from cummulative array starting at iPar
1404 0 : local_cum_weight(iPar : nPar) = local_cum_weight(iPar : nPar) - weight
1405 : ! (4.3) Renormalize cummulative weight to one
1406 0 : if (count(neighborhood) .lt. size(truepara)) then
1407 0 : norm = 1.0_dp / local_cum_weight(nPar)
1408 : else
1409 : norm = 1.0_dp
1410 : end if
1411 0 : local_cum_weight(:) = local_cum_weight(:) * norm
1412 : !
1413 : end do
1414 :
1415 1287500 : end subroutine generate_neighborhood_weight_dp
1416 :
1417 : END MODULE mo_anneal
|
