Line data Source code
1 : !> \file mo_mcmc.f90
2 : !> \brief \copybrief mo_mcmc
3 : !> \details \copydetails mo_mcmc
4 :
5 : !> \brief Monte Carlo Markov Chain sampling.
6 : !> \details This module is Monte Carlo Markov Chain sampling of a posterior parameter distribution.
7 : !> \authors Maren Goehler, Juliane Mai
8 : !> \date Aug 2012
9 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
10 : !! FORCES is released under the LGPLv3+ license \license_note
11 : MODULE mo_mcmc
12 :
13 : USE mo_kind, only : i4, i8, dp
14 : USE mo_xor4096, only : xor4096, xor4096g, get_timeseed, n_save_state
15 : USE mo_append, only : append
16 : USE mo_moment, only : stddev
17 : !$ USE omp_lib, only: OMP_GET_NUM_THREADS
18 : use mo_optimization_utils, only : eval_interface, objective_interface
19 : use mo_message, only : error_message
20 : #ifdef FORCES_WITH_NETCDF
21 : use mo_ncwrite, only : dump_netcdf
22 : #endif
23 :
24 : IMPLICIT NONE
25 :
26 : PUBLIC :: mcmc ! Sample posterior parameter distribution (measurement errors are either known or modeled)
27 : PUBLIC :: mcmc_stddev ! Sample posterior parameter distribution (measurement errors are neither known nor modeled)
28 :
29 : !-----------------------------------------------------------------------------------------------
30 :
31 : !> \brief This module is Monte Carlo Markov Chain sampling of a posterior parameter distribution
32 : !! (measurement errors are either known or modeled).
33 :
34 : !> \details Sample posterior parameter distribution with Metropolis Algorithm.\n
35 : !! This sampling is performed in two steps, i.e. the burn-in phase for adjusting
36 : !! model dependent parameters for the second step which is the proper
37 : !! sampling using the Metropolis Hastings Algorithm.\n\n
38 : !!
39 : !! This sampler does not change the best parameter set, i.e.
40 : !! it cannot be used as an optimiser.\n
41 : !! However, the serial and the parallel version give therefore the bitwise same results.\n
42 : !!
43 : !! <b> 1. Burn in phase: Find the optimal step size </b>
44 : !!
45 : !! \b <i>Purpose:</i>
46 : !!
47 : !! Find optimal stepsize for each parameter such that the
48 : !! acceptance ratio converges to a value around 0.3.
49 : !!
50 : !! <b> <i>Important variables:</i></b>
51 : !!
52 : !! <b> Variable </b> | <b> Description </b>
53 : !! ---------------------- | -----------------------------------------------------------------------
54 : !! burnin_iter | length of markov chain performed to calculate acceptance ratio\n
55 : !! acceptance ratio | ratio between accepted jumps and all trials (LEN)\n
56 : !! acceptance multiplier | stepsize of a parameter is multiplied with this value when jump is accepted (initial : 1.01)
57 : !! rejection multiplier | stepsize of a parameter is multiplied with this value when jump is rejected (initial : 0.99
58 : !! ^ | and will never be changed)
59 : !! stepsize | a new parameter value is chosen based on a uniform distribution
60 : !! ^ | pnew_i = pold_i + Unif(-stepsize_i, stepsize_i) (initial : stepsize_i = 1.0 for all i)
61 : !!
62 : !! \b <i>Algorithm:</i>
63 : !! <ol>
64 : !! <li><i>start a new markov chain of length burnin_iter with initial parameter set is the OPTIMAL one</i>\n
65 : !! - select a set of parameters to change:\n
66 : !! * accurate --> choose one parameter,\n
67 : !! * comput. efficient --> choose all parameters,\n
68 : !! * moderate accurate & efficient --> choose half of the parameters\n
69 : !! - change parameter(s) based on their stepsize\n
70 : !! - decide whether changed parameter set is accepted or rejected:\n
71 : !! * \f$ \mathrm{odds Ratio} = \frac{\mathrm{likelihood}(p_{new})}{\mathrm{likelihood}(p_{old})} \f$\n
72 : !! * random number r = Uniform[0,1]\n
73 : !! * odds Ratio > 0 --> positive accept\n
74 : !! odds Ratio > r --> negative accept\n
75 : !! odds Ratio < r --> reject\n
76 : !! - adapt stepsize of parameters changed:\n
77 : !! * accepted step: stepsize_i = stepsize_i * accptance multiplier\n
78 : !! * rejected step: stepsize_i = stepsize_i * rejection multiplier\n
79 : !! - if step is accepted: for all changed parameter(s) change stepsize\n
80 : !!
81 : !! <li><i>calculate acceptance ratio of the Markov Chain</i>\n
82 : !!
83 : !! <li>adjust acceptance multiplier acc_mult and
84 : !! store good ratios in history list\n
85 : !! - acceptance ratio < 0.23 --> acc_mult = acc_mult * 0.99\n
86 : !! delete history list\n
87 : !! - acceptance ratio > 0.44 --> acc_mult = acc_mult * 1.01\n
88 : !! delete history list\n
89 : !! - 0.23 < acceptance ratio < 0.44\n
90 : !! add acceptance ratio to history list\n
91 : !!
92 : !! <li>check if already 10 values are stored in history list and
93 : !! if they have converged to a value above 0.3\n
94 : !! ( mean above 0.3 and variance less \f$\sqrt{1/12*0.05^2}\f$ = Variance
95 : !! of uniform [acc_ratio +/- 2.5%] )\n
96 : !! - if check is positive abort and save stepsizes\n
97 : !! else goto (1)\n\n
98 : !! </ol>
99 : !!
100 : !! <b> 2. Monte Carlo Markov Chain: Sample posterioir distribution of parameter </b>
101 : !!
102 : !! \b <i>Purpose:</i>
103 : !!
104 : !! use the previous adapted stepsizes and perform ONE monte carlo markov chain\n
105 : !! the accepted parameter sets show the posterior distribution of parameters
106 : !!
107 : !! <b> <i>Important variables:</i></b>
108 : !!
109 : !! <b> Variable </b> | <b> Description </b>
110 : !! ---------------------- | -----------------------------------------------------------------------
111 : !! iter_mcmc | length of the markov chain (>> iter_burnin)\n
112 : !! stepsize | a new parameter value is chosen based on a uniform distribution
113 : !! pnew_i = pold_i + Unif(-stepsize_i, stepsize_i) use stepsizes of the burn-in (1)
114 : !!
115 : !! \b <i>Algorithm:</i>
116 : !! <ol>
117 : !! <li> select a set of parameters to change\n
118 : !! * accurate --> choose one parameter,\n
119 : !! * comput. efficient --> choose all parameters,\n
120 : !! * moderate accurate & efficient --> choose half of the parameters\n
121 : !! <li> change parameter(s) based on their stepsize\n
122 : !! <li> decide whether changed parameter set is accepted or rejected:\n
123 : !! * \f$ \mathrm{odds Ratio} = \frac{\mathrm{likelihood}(p_{new})}{\mathrm{likelihood}(p_{old})} \f$\n
124 : !! * random number r = Uniform[0,1]\n
125 : !! * odds Ratio > 0 --> positive accept\n
126 : !! odds Ratio > r --> negative accept\n
127 : !! odds Ratio < r --> reject\n
128 : !! <li> if step is accepted: save parameter set\n
129 : !! <li> goto (1)\n
130 : !! </ol>
131 : !!
132 : !! \b Example
133 : !!
134 : !! \code{.f90}
135 : !! call mcmc( likelihood, para, rangePar, mcmc_paras, burnin_paras, &
136 : !! seed_in=seeds, printflag_in=printflag, maskpara_in=maskpara, &
137 : !! tmp_file=tmp_file, loglike_in=loglike, &
138 : !! ParaSelectMode_in=ParaSelectMode, &
139 : !! iter_burnin_in=iter_burnin, iter_mcmc_in=iter_mcmc, &
140 : !! chains_in=chains, stepsize_in=stepsize)
141 : !! \endcode
142 : !! See also example in test directory.
143 : !!
144 : !! \b Literature
145 : !!
146 : !! 1. Gelman et. al (1995)
147 : !! _Bayesian Data Analyis_. Chapman & Hall.\n
148 : !! 2. Gelman et. al (1997)
149 : !! _Efficient Metropolis Jumping Rules: Baysian Statistics 5_. pp. 599-607\n
150 : !! 3. Tautenhahn et. al (2012)
151 : !! _Beyond distance-invariant survival in inverse recruitment modeling:
152 : !! A case study in Siberian Pinus sylvestris forests_. Ecological Modelling, 233,
153 : !! 90-103. doi:10.1016/j.ecolmodel.2012.03.009.
154 : !!
155 :
156 : !> \param[in] "real(dp) :: likelihood(x)" Interface Function which calculates likelihood
157 : !! of given parameter set x
158 : !> \param[in] "real(dp) :: para(:)" Inital parameter set (should be GOOD approximation
159 : !! of best parameter set)
160 : !> \param[in] "real(dp) :: rangePar(size(para),2)" Min/max range of parameters
161 : !> \param[in] "integer(i8), optional :: seed_in" User seed to initialise the random number generator \n
162 : !! (default: none --> initialized with timeseed)
163 : !> \param[in] "logical, optional :: printflag_in" Print of output on command line \n
164 : !! (default: .False.)
165 : !> \param[in] "logical, optional :: maskpara_in(size(para))"
166 : !! Parameter will be sampled (.True.) or not (.False.) \n
167 : !! (default: .True.)
168 : !> \param[in] "character(len=*), optional :: tmp_file" filename for temporal data saving: \n
169 : !! every iter_mcmc_in iterations parameter sets are
170 : !! appended to this file \n
171 : !! the number of the chain will be prepended to filename \n
172 : !! output format: netcdf \n
173 : !! (default: no file writing)
174 : !> \param[in] "logical, optional :: loglike_in" true if loglikelihood function is given instead of
175 : !! likelihood function \n
176 : !! (default: .false.)
177 : !> \param[in] "integer(i4), optional :: ParaSelectMode_in"
178 : !! How many parameters will be changed at once?
179 : !! - half of the parameter --> 1_i4
180 : !! - only one parameter --> 2_i4
181 : !! - all parameter --> 3_i4
182 : !! (default: 2_i4)
183 : !> \param[in] "integer(i4), optional :: iter_burnin_in" Length of Markov chains of initial burn-in part \n
184 : !! (default: Max(250, 200*count(maskpara)) )
185 : !> \param[in] "integer(i4), optional :: iter_mcmc_in" Length of Markov chains of proper MCMC part \n
186 : !! (default: 1000 * count(maskpara) )
187 : !> \param[in] "integer(i4), optional :: chains_in" number of parallel mcmc chains \n
188 : !! (default: 5_i4)
189 : !> \param[in] "real(dp), DIMENSION(size(para,1)), optional :: stepsize_in"
190 : !! stepsize for each parameter \n
191 : !! if given burn-in is discarded \n
192 : !! (default: none --> adjusted in burn-in)
193 : !> \param[out] "real(dp), allocatable :: mcmc_paras(:,:)" Parameter sets sampled in proper MCMC part of algorithm
194 : !> \param[out] "real(dp), allocatable :: burnin_paras(:,:)" Parameter sets sampled during burn-in part of algorithm
195 :
196 : !> \note
197 : !> Likelihood has to be defined as a function interface\n
198 : !> The maximal number of parameters is 1000.
199 :
200 : !> \authors Juliane Mai
201 : !> \date Nov 2014
202 :
203 : INTERFACE mcmc
204 : MODULE PROCEDURE mcmc_dp
205 : END INTERFACE mcmc
206 :
207 : !-----------------------------------------------------------------------------------------------
208 :
209 : !> \brief This module is Monte Carlo Markov Chain sampling of a posterior parameter distribution
210 : !! (measurement errors are neither known nor modeled).
211 :
212 : !> \details Sample posterior parameter distribution with Metropolis Algorithm.\n
213 : !! This sampling is performed in two steps, i.e. the burn-in phase for adjusting
214 : !! model dependent parameters for the second step which is the proper
215 : !! sampling using the Metropolis Hastings Algorithm.\n\n
216 :
217 : !> <b> 1. Burn in phase: Find the optimal step size </b>
218 : !!
219 : !! \b <i>Purpose:</i>
220 : !!
221 : !! Find optimal stepsize for each parameter such that the
222 : !! acceptance ratio converges to a value around 0.3.\n
223 : !!
224 : !! <b> <i>Important variables:</i> </b>
225 : !!
226 : !! <b> Variable </b> | <b> Description </b>
227 : !! ---------------------- | -----------------------------------------------------------------------
228 : !! burnin_iter | length of markov chain performed to calculate acceptance ratio\n
229 : !! acceptance ratio | ratio between accepted jumps and all trials (LEN)\n
230 : !! acceptance multiplier | stepsize of a parameter is multiplied with this value when jump is accepted (initial : 1.01)
231 : !! rejection multiplier | stepsize of a parameter is multiplied with this value when jump is rejected (initial : 0.99 and
232 : !! ^ | will never be changed)
233 : !! stepsize | a new parameter value is chosen based on a uniform distribution
234 : !! ^ | pnew_i = pold_i + Unif(-stepsize_i, stepsize_i) (initial : stepsize_i = 1.0 for all i)
235 : !!
236 : !! \b <i>Algorithm:</i>
237 : !! <ol>
238 : !! <li><i>start a new markov chain of length burnin_iter with initial parameter set is the OPTIMAL one</i>\n
239 : !! - select a set of parameters to change:\n
240 : !! * accurate --> choose one parameter,\n
241 : !! * comput. efficient --> choose all parameters,\n
242 : !! * moderate accurate & efficient --> choose half of the parameters\n
243 : !! - change parameter(s) based on their stepsize\n
244 : !! - decide whether changed parameter set is accepted or rejected:\n
245 : !! * \f$ \mathrm{odds Ratio} = \frac{\mathrm{likelihood}(p_{new})}{\mathrm{likelihood}(p_{old})} \f$\n
246 : !! * random number r = Uniform[0,1]\n
247 : !! * odds Ratio > 0 --> positive accept\n
248 : !! odds Ratio > r --> negative accept\n
249 : !! odds Ratio < r --> reject\n
250 : !! - adapt stepsize of parameters changed:\n
251 : !! * accepted step: stepsize_i = stepsize_i * accptance multiplier\n
252 : !! * rejected step: stepsize_i = stepsize_i * rejection multiplier\n
253 : !! - if step is accepted: for all changed parameter(s) change stepsize\n
254 : !!
255 : !! <li><i>calculate acceptance ratio of the Markov Chain</i>\n
256 : !!
257 : !! <li>adjust acceptance multiplier acc_mult and
258 : !! store good ratios in history list\n
259 : !! - acceptance ratio < 0.23 --> acc_mult = acc_mult * 0.99\n
260 : !! delete history list\n
261 : !! - acceptance ratio > 0.44 --> acc_mult = acc_mult * 1.01\n
262 : !! delete history list\n
263 : !! - 0.23 < acceptance ratio < 0.44\n
264 : !! add acceptance ratio to history list\n
265 : !!
266 : !! <li>check if already 10 values are stored in history list and
267 : !! if they have converged to a value above 0.3\n
268 : !! ( mean above 0.3 and variance less \f$\sqrt{1/12*0.05^2}\f$ = Variance
269 : !! of uniform [acc_ratio +/- 2.5%] )\n
270 : !! - if check is positive abort and save stepsizes\n
271 : !! else goto (1)\n\n
272 : !! </ol>
273 : !!
274 : !! <b> 2. Monte Carlo Markov Chain: Sample posterioir distribution of parameter </b>
275 : !!
276 : !! \b <i>Purpose:</i>
277 : !!
278 : !! use the previous adapted stepsizes and perform ONE monte carlo markov chain\n
279 : !! the accepted parameter sets show the posterior distribution of parameters\n
280 : !!
281 : !! <b> <i>Important variables:</i> </b>
282 : !!
283 : !! <b> Variable </b> | <b> Description </b>
284 : !! ---------------------- | -----------------------------------------------------------------------
285 : !! iter_mcmc | length of the markov chain (>> iter_burnin)\n
286 : !! stepsize | a new parameter value is chosen based on a uniform distribution
287 : !! ^ | pnew_i = pold_i + Unif(-stepsize_i, stepsize_i) use stepsizes of the burn-in (1)\n
288 : !!
289 : !! \b <i>Algorithm:</i>
290 : !! <ol>
291 : !! <li> select a set of parameters to change\n
292 : !! * accurate --> choose one parameter,\n
293 : !! * comput. efficient --> choose all parameters,\n
294 : !! * moderate accurate & efficient --> choose half of the parameters\n
295 : !! <li> change parameter(s) based on their stepsize\n
296 : !! <li> decide whether changed parameter set is accepted or rejected:\n
297 : !! * \f$ \mathrm{odds Ratio} = \frac{\mathrm{likelihood}(p_{new})}{\mathrm{likelihood}(p_{old})} \f$\n
298 : !! * random number r = Uniform[0,1]\n
299 : !! * odds Ratio > 0 --> positive accept\n
300 : !! odds Ratio > r --> negative accept\n
301 : !! odds Ratio < r --> reject\n
302 : !! <li> if step is accepted: save parameter set\n
303 : !! <li> goto (1)\n
304 : !! </ol>
305 : !!
306 : !! \b Example
307 : !!
308 : !! \code{.f90}
309 : !! call mcmc( likelihood, para, rangePar, mcmc_paras, burnin_paras, &
310 : !! seed_in=seeds, printflag_in=printflag, maskpara_in=maskpara, &
311 : !! tmp_file=tmp_file, loglike_in=loglike, &
312 : !! ParaSelectMode_in=ParaSelectMode, &
313 : !! iter_burnin_in=iter_burnin, iter_mcmc_in=iter_mcmc, &
314 : !! chains_in=chains, stepsize_in=stepsize)
315 : !! \endcode
316 : !! See also example in test directory.
317 : !!
318 : !! \b Literature
319 : !!
320 : !! 1. Gelman et. al (1995)
321 : !! _Bayesian Data Analyis_. Chapman & Hall.
322 : !! 2. Gelman et. al (1997)
323 : !! _Efficient Metropolis Jumping Rules: Bayesian Statistics 5_. pp. 599-607
324 : !! 3. Tautenhahn et. al (2012)
325 : !! _Beyond distance-invariant survival in inverse recruitment modeling:
326 : !! A case study in Siberian Pinus sylvestris forests_. Ecological Modelling, 233,
327 : !! 90-103. doi:10.1016/j.ecolmodel.2012.03.009.
328 : !!
329 :
330 : !> \param[in] "real(dp) :: likelihood(x,sigma,stddev_new,likeli_new)"
331 : !> Interface Function which calculates likelihood
332 : !> of given parameter set x and given standard deviation sigma
333 : !> and returns optionally the standard deviation stddev_new
334 : !> of the errors using x and
335 : !> likelihood likeli_new using stddev_new
336 :
337 : !> \param[in] "real(dp) :: para(:)" Inital parameter set (should be GOOD approximation
338 : !> of best parameter set)
339 :
340 : !> \param[in] "real(dp) :: rangePar(size(para),2)" Min/max range of parameters
341 : !> \param[in] "integer(i8), optional :: seed_in" User seed to initialise the random number generator \n
342 : !> (default: none --> initialized with timeseed)
343 :
344 : !> \param[in] "logical, optional :: printflag_in" Print of output on command line \n
345 : !> (default: .False.)
346 :
347 : !> \param[in] "logical, optional :: maskpara_in(size(para))"
348 : !> Parameter will be sampled (.True.) or not (.False.) \n
349 : !> (default: .True.)
350 :
351 : !> \param[in] "character(len=*), optional :: tmp_file" filename for temporal data saving: \n
352 : !> every iter_mcmc_in iterations parameter sets are
353 : !> appended to this file \n
354 : !> the number of the chain will be prepended to filename \n
355 : !> output format: netcdf \n
356 : !> (default: no file writing)
357 :
358 : !> \param[in] "logical, optional :: loglike_in" true if loglikelihood function is given instead of
359 : !> likelihood function \n
360 : !> (default: .false.)
361 :
362 : !> \param[in] "integer(i4), optional :: ParaSelectMode_in"
363 : !> How many parameters will be changed at once?
364 : !> - half of the parameter --> 1_i4
365 : !> - only one parameter --> 2_i4
366 : !> - all parameter --> 3_i4
367 : !> (default: 2_i4)
368 :
369 : !> \param[in] "integer(i4), optional :: iter_burnin_in" Length of Markov chains of initial burn-in part \n
370 : !> (default: Max(250, 200*count(maskpara)) )
371 :
372 : !> \param[in] "integer(i4), optional :: iter_mcmc_in" Length of Markov chains of proper MCMC part \n
373 : !> (default: 1000 * count(maskpara) )
374 :
375 : !> \param[in] "integer(i4), optional :: chains_in" number of parallel mcmc chains \n
376 : !> (default: 5_i4)
377 :
378 : !> \param[in] "real(dp), DIMENSION(size(para,1)), optional :: stepsize_in"
379 : !> stepsize for each parameter \n
380 : !> if given burn-in is discarded \n
381 : !> (default: none --> adjusted in burn-in)
382 :
383 : !> \param[out] "real(dp), allocatable :: mcmc_paras(:,:)" Parameter sets sampled in proper MCMC part of algorithm
384 : !> \param[out] "real(dp), allocatable :: burnin_paras(:,:)" Parameter sets sampled during burn-in part of algorithm
385 :
386 : !> \note Restrictions: Likelihood has to be defined as a function interface
387 :
388 : !> \author Maren Goehler
389 : !> \date Sep. 2012
390 : !! - Created using copy of Simulated Annealing:\n
391 : !! - constant temperature T\n
392 : !! - burn-in for stepsize adaption\n
393 : !! - acceptance/ rejection multiplier\n
394 :
395 : !> \author Juliane Mai
396 : !> \date Sep. 2012
397 : !! - Cleaning code and introduce likelihood:\n
398 : !! - likelihood instead of objective function\n
399 : !! - odds ratio\n
400 : !! - convergence criteria for burn-in\n
401 : !! - different modes of parameter selection\n
402 : !! - OpenMP for chains of MCMC\n
403 : !! - optional file for temporal output\n
404 : !> \date Nov 2012
405 : !! - Temporary file writing as NetCDF
406 : !> \date Aug 2013
407 : !! - New likelihood interface to reduce number of function evaluations. Only one seed has to be given.\n
408 : !> \date Nov 2014
409 : !! - add new routine for mcmc:\n
410 : !! - mcmc_stddev, i.e. mcmc with likelihood with "faked sigma".\n
411 : !! - mcmc, i.e. mcmc with "real" likelihood (sigma is given by e.g. error model).\n
412 : !> \date Apr 2015
413 : !! - handling of array-like variables in restart-namelists
414 :
415 : !> \author Mathias Cuntz and Juliane Mai
416 : !> \date Feb 2015
417 : !! - restart
418 :
419 : INTERFACE mcmc_stddev
420 : MODULE PROCEDURE mcmc_stddev_dp
421 : END INTERFACE mcmc_stddev
422 :
423 : !-----------------------------------------------------------------------------------------------
424 :
425 : PRIVATE
426 :
427 : !-----------------------------------------------------------------------------------------------
428 :
429 : CONTAINS
430 :
431 : !-----------------------------------------------------------------------------------------------
432 :
433 4 : SUBROUTINE mcmc_dp(eval, likelihood, para, rangePar, & ! obligatory IN
434 : mcmc_paras, burnin_paras, & ! obligatory OUT
435 2 : seed_in, printflag_in, maskpara_in, & ! optional IN
436 : restart, restart_file, & ! optional IN: if mcmc is restarted and file which contains restart variables
437 : tmp_file, & ! optional IN: filename for temporal output of
438 : ! ! MCMC parasets
439 : loglike_in, & ! optional IN: true if loglikelihood is given
440 : ParaSelectMode_in, & ! optional IN: (=1) half, (=2) one, (=3) all
441 : iter_burnin_in, & ! optional IN: markov length of (1) burn-in
442 : iter_mcmc_in, & ! optional IN: markov length of (2) mcmc
443 : chains_in, & ! optional IN: number of parallel chains of MCMC
444 0 : stepsize_in) ! optional IN: stepsize for each param. (no burn-in)
445 :
446 : IMPLICIT NONE
447 :
448 : procedure(eval_interface), INTENT(IN), POINTER :: eval
449 : procedure(objective_interface), intent(in), pointer :: likelihood
450 :
451 : REAL(DP), DIMENSION(:, :), INTENT(IN) :: rangePar ! range for each parameter
452 : REAL(DP), DIMENSION(:), INTENT(IN) :: para ! initial parameter i
453 : INTEGER(I8), OPTIONAL, INTENT(IN) :: seed_in ! Seeds of random numbers: dim1=chains, dim2=3
454 : ! ! (DEFAULT: Get_timeseed)
455 : LOGICAL, OPTIONAL, INTENT(IN) :: printflag_in ! If command line output is written (.true.)
456 : ! ! (DEFAULT: .false.)
457 : LOGICAL, OPTIONAL, DIMENSION(size(para, 1)), &
458 : INTENT(IN) :: maskpara_in ! true if parameter will be optimized
459 : ! ! false if parameter is discarded in optimization
460 : ! ! (DEFAULT = .true.)
461 : logical, OPTIONAL, INTENT(in) :: restart ! if .true., start from restart file
462 : ! ! (DEFAULT: .false.)
463 : character(len = *), OPTIONAL, INTENT(in) :: restart_file ! restart file name
464 : ! ! (DEFAULT: mo_mcmc.restart)
465 : LOGICAL, OPTIONAL, INTENT(IN) :: loglike_in ! true if loglikelihood is given instead of likelihood
466 : ! ! (DEFAULT: .false.)
467 : INTEGER(I4), OPTIONAL, INTENT(IN) :: ParaSelectMode_in ! how many parameters changed per iteration:
468 : ! ! 1: half of the parameters
469 : ! ! 2: only one (DEFAULT)
470 : ! ! 3: all
471 : INTEGER(I4), OPTIONAL, INTENT(IN) :: iter_burnin_in ! # iterations before checking ac_ratio
472 : ! ! DEFAULT= MIN(250, 20_i4*n)
473 : INTEGER(I4), OPTIONAL, INTENT(IN) :: iter_mcmc_in ! # iterations per chain for posterior sampling
474 : ! ! DEFAULT= 1000_i4*n
475 : INTEGER(I4), OPTIONAL, INTENT(IN) :: chains_in ! number of parallel mcmc chains
476 : ! ! DEFAULT= 5_i4
477 : REAL(DP), OPTIONAL, DIMENSION(size(para, 1)), &
478 : INTENT(IN) :: stepsize_in ! stepsize for each parameter
479 : ! ! if given burn-in is discarded
480 : CHARACTER(len = *), OPTIONAL, INTENT(IN) :: tmp_file ! filename for temporal data saving: every iter_mcmc_in
481 : ! ! iterations parameter sets are appended to this file
482 : ! ! the number of the chain will be prepended to filename
483 : ! ! DEFAULT= no file writing
484 : REAL(DP), DIMENSION(:, :), ALLOCATABLE, INTENT(OUT) :: mcmc_paras
485 : ! ! array to save para values of MCMC runs,
486 : ! ! dim1=sets of all chains, dim2=paras
487 : REAL(DP), DIMENSION(:, :), ALLOCATABLE, INTENT(OUT) :: burnin_paras
488 : ! ! array to save para values of Burn in run
489 :
490 :
491 : ! local variables
492 : INTEGER(I4) :: n ! Number of parameters
493 : LOGICAL :: printflag ! If command line output is written
494 0 : LOGICAL, DIMENSION(size(para, 1)) :: maskpara ! true if parameter will be optimized
495 : LOGICAL, DIMENSION(1000) :: dummy_maskpara ! dummy mask_para (only for namelist)
496 2 : INTEGER(I4), DIMENSION(:), ALLOCATABLE :: truepara ! indexes of parameters to be optimized
497 : LOGICAL :: loglike ! if loglikelihood is given
498 : LOGICAL :: skip_burnin ! if stepsize is given --> burnin is skipped
499 : logical :: itmp_file ! if temporal results wanted
500 : character(len = 1024) :: istmp_file ! local copy of file for temporal results output
501 : logical :: irestart ! if restart wanted
502 : character(len = 1024) :: isrestart_file ! local copy of restart file name
503 :
504 : ! for random numbers
505 2 : INTEGER(I8), dimension(:, :), allocatable :: seeds ! Seeds of random numbers: dim1=chains, dim2=3
506 2 : REAL(DP) :: RN1, RN2, RN3 ! Random numbers
507 2 : integer(I8), dimension(:, :), allocatable :: save_state_1 ! optional argument for restarting RN stream 1:
508 2 : integer(I8), dimension(:, :), allocatable :: save_state_2 ! optional argument for restarting RN stream 2:
509 2 : integer(I8), dimension(:, :), allocatable :: save_state_3 ! optional argument for restarting RN stream 3:
510 : ! ! dim1=chains, dim2=n_save_state
511 :
512 : ! Dummies
513 2 : REAL(DP), DIMENSION(:, :, :), ALLOCATABLE :: tmp
514 : integer(I4) :: idummy
515 : logical :: oddsSwitch1, oddsSwitch2
516 : character(100) :: str
517 : character(200) :: outputfile
518 : integer(i4) :: slash_pos
519 : integer(i4) :: len_filename
520 : character(200) :: filename
521 : character(200) :: path
522 :
523 : ! FOR BURN-IN AND MCMC
524 2 : REAL(DP), DIMENSION(:, :, :), ALLOCATABLE :: mcmc_paras_3d ! array to save para values of MCMC runs,
525 : ! ! dim1=sets, dim2=paras, dim3=chain
526 : integer(I4) :: i
527 4 : integer(I4), dimension(:), allocatable :: Ipos, Ineg
528 : logical :: iStop
529 : INTEGER(I4) :: ParaSelectMode ! how many parameters are changed per jump
530 : INTEGER(I4) :: iter_burnin ! number fo iterations before checking ac_ratio
531 : INTEGER(I4) :: iter_mcmc ! number fo iterations for posterior sampling
532 : INTEGER(I4) :: chains ! number of parallel mcmc chains
533 : INTEGER(I4) :: chain ! counter for chains, 1...chains
534 2 : REAL(DP) :: accMult ! acceptance multiplier for stepsize
535 2 : REAL(DP) :: rejMult ! rejection multiplier for stepsize
536 0 : LOGICAL, DIMENSION(size(para, 1)) :: ChangePara ! logical array to switch if parameter is changed
537 : INTEGER(I4) :: trial ! number of trials for burn-in
538 8 : REAL(DP), DIMENSION(size(para, 1)) :: stepsize ! stepsize adjusted by burn-in and used by mcmc
539 2002 : REAL(DP), DIMENSION(1000) :: dummy_stepsize ! dummy stepsize (only for namelist)
540 8 : REAL(DP), DIMENSION(size(para, 1)) :: paraold ! old parameter set
541 10 : REAL(DP), DIMENSION(size(para, 1)) :: paranew ! new parameter set
542 10 : REAL(DP), DIMENSION(size(para, 1)) :: parabest ! best parameter set overall
543 2002 : REAL(DP), DIMENSION(1000) :: dummy_parabest ! dummy parabest (only for namelist)
544 10 : REAL(DP), DIMENSION(size(para, 1)) :: initial_paraset_mcmc ! best parameterset found in burn-in
545 2002 : REAL(DP), DIMENSION(1000) :: dummy_initial_paraset_mcmc ! dummy initial_paraset_mcmc (only for namelist)
546 2 : REAL(DP) :: likeliold ! likelihood of old parameter set
547 2 : REAL(DP) :: likelinew ! likelihood of new parameter set
548 2 : REAL(DP) :: likelibest ! likelihood of best parameter set overall
549 : INTEGER(I4) :: markov ! counter for markov chain
550 2 : REAL(DP), DIMENSION(:, :), ALLOCATABLE :: burnin_paras_part ! accepted parameter sets of one MC in burn-in
551 2 : REAL(DP) :: oddsRatio ! ratio of likelihoods = likelinew/likeliold
552 2 : REAL(DP), dimension(:), allocatable :: accRatio ! acceptance ratio = (pos/neg) accepted/iter_burnin
553 2 : REAL(DP) :: accratio_stddev ! stddev of accRatios in history
554 2 : REAL(DP), DIMENSION(:), ALLOCATABLE :: history_accratio ! history of 'good' acceptance ratios
555 : INTEGER(I4) :: accratio_n ! number of 'good' acceptance ratios
556 :
557 : ! for checking convergence of MCMC
558 2 : real(dp), allocatable, dimension(:) :: sqrtR
559 2 : real(dp), allocatable, dimension(:) :: vDotJ
560 2 : real(dp), allocatable, dimension(:) :: s2
561 2 : real(dp) :: vDotDot
562 2 : real(dp) :: B
563 2 : real(dp) :: W
564 : integer(i4) :: n_start, n_end, iPar
565 : LOGICAL :: converged ! if MCMC already converged
566 :
567 : ! for OMP
568 : integer(i4) :: n_threads
569 :
570 : namelist /restartnml1/ &
571 : n, printflag, dummy_maskpara, loglike, skip_burnin, &
572 : iStop, ParaSelectMode, iter_burnin, &
573 : iter_mcmc, chains, accMult, rejMult, trial, dummy_stepsize, &
574 : dummy_parabest, likelibest, dummy_initial_paraset_mcmc, &
575 : accratio_stddev, &
576 : accratio_n, vDotDot, B, W, converged, &
577 : n_threads, &
578 : itmp_file, istmp_file
579 :
580 : namelist /restartnml2/ &
581 : truepara, seeds, save_state_1, save_state_2, save_state_3, &
582 : iPos, iNeg, mcmc_paras_3d, &
583 : accRatio, &
584 : sqrtR, vDotJ, s2
585 :
586 : ! CHECKING OPTIONALS
587 :
588 2 : if (present(restart)) then
589 2 : irestart = restart
590 : else
591 : irestart = .false.
592 : end if
593 :
594 2 : if (present(restart_file)) then
595 2 : isrestart_file = restart_file
596 : else
597 0 : isrestart_file = 'mo_mcmc.restart'
598 : end if
599 :
600 2 : if (.not. irestart) then
601 :
602 1 : if (present(tmp_file)) then
603 1 : itmp_file = .true.
604 1 : istmp_file = tmp_file
605 : else
606 0 : itmp_file = .false.
607 0 : istmp_file = ''
608 : end if
609 :
610 1 : if (present(loglike_in)) then
611 1 : loglike = loglike_in
612 : else
613 0 : loglike = .false.
614 : end if
615 :
616 1 : if (present(maskpara_in)) then
617 4 : if (count(maskpara_in) .eq. 0_i4) then
618 0 : stop 'Input argument maskpara: At least one element has to be true'
619 : else
620 4 : maskpara = maskpara_in
621 : end if
622 : else
623 0 : maskpara = .true.
624 : end if
625 :
626 6 : allocate (truepara(count(maskpara)))
627 1 : idummy = 0_i4
628 4 : do i = 1, size(para, 1)
629 4 : if (maskpara(i)) then
630 3 : idummy = idummy + 1_i4
631 3 : truepara(idummy) = i
632 : end if
633 : end do
634 :
635 1 : n = size(truepara, 1)
636 :
637 1 : if (present(ParaSelectMode_in)) then
638 1 : ParaSelectMode = ParaSelectMode_in
639 : else
640 : ! change only one parameter per jump
641 0 : ParaSelectMode = 2_i4
642 : end if
643 :
644 : ! after how many iterations do we compute ac_ratio???
645 1 : if (present(iter_burnin_in)) then
646 0 : if (iter_burnin_in .le. 0_i4) then
647 0 : stop 'Input argument iter_burn_in must be greater than 0!'
648 : else
649 0 : iter_burnin = iter_burnin_in
650 : end if
651 : else
652 1 : iter_burnin = Max(250_i4, 1000_i4 * n)
653 : end if
654 :
655 : ! how many iterations ('jumps') are performed in MCMC
656 : ! iter_mcmc_in is handled later properly (after acceptance ratio of burn_in is known)
657 1 : if (present(iter_mcmc_in)) then
658 0 : if (iter_mcmc_in .le. 0_i4) then
659 0 : stop 'Input argument iter_mcmc must be greater than 0!'
660 : else
661 0 : iter_mcmc = iter_mcmc_in
662 : end if
663 : else
664 1 : iter_mcmc = 1000_i4 * n
665 : end if
666 :
667 1 : if (present(chains_in)) then
668 0 : if (chains_in .lt. 2_i4) then
669 0 : stop 'Input argument chains must be at least 2!'
670 : end if
671 0 : chains = chains_in
672 : else
673 1 : chains = 5_i4
674 : end if
675 :
676 1 : if (present(stepsize_in)) then
677 0 : stepsize = stepsize_in
678 0 : skip_burnin = .true.
679 : else
680 1 : skip_burnin = .false.
681 : end if
682 :
683 1 : if (present(printflag_in)) then
684 1 : printflag = printflag_in
685 : else
686 0 : printflag = .false.
687 : end if
688 :
689 1 : n_threads = 1
690 : !$ write(*,*) '--------------------------------------------------'
691 : !$ write(*,*) ' This program is parallel.'
692 : !$OMP parallel
693 : !$ n_threads = OMP_GET_NUM_THREADS()
694 : !$OMP end parallel
695 : !$ write(*,*) ' ',chains,' MCMC chains will run in ',n_threads,' threads'
696 : !$ write(*,*) '--------------------------------------------------'
697 :
698 1 : if (printflag) then
699 1 : write(*, *) 'Following parameters will be sampled with MCMC: ', truepara
700 : end if
701 :
702 3 : allocate(seeds(chains, 3))
703 3 : allocate(save_state_1(chains, n_save_state))
704 2 : allocate(save_state_2(chains, n_save_state))
705 2 : allocate(save_state_3(chains, n_save_state))
706 6 : allocate(Ipos(chains), Ineg(chains), accRatio(chains))
707 3 : allocate(vDotJ(chains), s2(chains))
708 4 : allocate(sqrtR(size(para)))
709 :
710 6 : Ipos = -9999
711 6 : Ineg = -9999
712 6 : accRatio = -9999.0_dp
713 6 : vDotJ = -9999.0_dp
714 6 : s2 = -9999.0_dp
715 4 : sqrtR = -9999.0_dp
716 :
717 1 : if (present(seed_in)) then
718 4 : seeds(1, :) = (/ seed_in, seed_in + 1000_i8, seed_in + 2000_i8 /)
719 : else
720 : ! Seeds depend on actual time
721 0 : call get_timeseed(seeds(1, :))
722 : end if
723 5 : do chain = 2, chains
724 17 : seeds(chain, :) = seeds(chain - 1_i4, :) + 3000_i8
725 : end do
726 :
727 6 : do chain = 1, chains
728 1325 : call xor4096(seeds(chain, 1), RN1, save_state = save_state_1(chain, :))
729 1325 : call xor4096(seeds(chain, 2), RN2, save_state = save_state_2(chain, :))
730 1326 : call xor4096g(seeds(chain, 3), RN3, save_state = save_state_3(chain, :))
731 : end do
732 19 : seeds = 0_i8
733 :
734 4 : parabest = para
735 :
736 : ! initialize likelihood
737 1 : likelibest = likelihood(parabest, eval)
738 :
739 : !----------------------------------------------------------------------
740 : ! (1) BURN IN
741 : !----------------------------------------------------------------------
742 :
743 1 : if (.not. skip_burnin) then
744 :
745 1 : if (printflag) then
746 1 : write(*, *) ''
747 1 : write(*, *) '--------------------------------------------------'
748 1 : write(*, *) 'Starting Burn-In (iter_burnin = ', iter_burnin, ')'
749 1 : write(*, *) '--------------------------------------------------'
750 1 : write(*, *) ''
751 : end if
752 :
753 : ! INITIALIZATION
754 :
755 : ! probably too large, but large enough to store values of one markovchain
756 4 : allocate(burnin_paras_part(iter_burnin, size(para)))
757 :
758 1 : if (allocated(burnin_paras)) deallocate(burnin_paras)
759 1 : if (allocated(history_accRatio)) deallocate(history_accRatio)
760 :
761 : ! parabestChanged = .false.
762 4 : stepsize = 1.0_dp
763 1 : trial = 1_i4
764 1 : iStop = .false.
765 1 : accMult = 1.01_dp
766 1 : rejMult = 0.99_dp
767 :
768 1 : if (printflag) then
769 1 : write(*, *) ' '
770 1 : write(*, *) 'Start Burn-In with: '
771 1 : write(*, *) ' parabest = ', parabest
772 1 : write(*, *) ' likelihood = ', likelibest
773 1 : write(*, *) ' '
774 : end if
775 :
776 : ! ----------------------------------------------------------------------------------
777 : ! repeats until convergence of acceptance ratio or better parameter set was found
778 : ! ----------------------------------------------------------------------------------
779 14 : convergedBURNIN : do while (.not. iStop)
780 :
781 78 : Ipos = 0_i4 ! positive accepted
782 78 : Ineg = 0_i4 ! negative accepted
783 52 : paraold = parabest
784 13 : likeliold = likelibest
785 :
786 : ! -------------------------------------------------------------------------------
787 : ! do a short-cut MCMC
788 : ! -------------------------------------------------------------------------------
789 39013 : markovchain : do markov = 1, iter_burnin
790 :
791 : ! (A) Generate new parameter set
792 156000 : ChangePara = .false.
793 156000 : paranew = paraold
794 : ! using RN from chain #1
795 : call GenerateNewParameterset_dp(ParaSelectMode, paraold, truepara, rangePar, stepsize, &
796 78000 : save_state_2(1, :), &
797 78000 : save_state_3(1, :), &
798 20787000 : paranew, ChangePara)
799 :
800 : ! (B) new likelihood
801 39000 : likelinew = likelihood(paranew, eval)
802 :
803 39000 : oddsSwitch1 = .false.
804 39000 : if (loglike) then
805 39000 : oddsRatio = likelinew - likeliold
806 39000 : if (oddsRatio .gt. 0.0_dp) oddsSwitch1 = .true.
807 : else
808 0 : oddsRatio = likelinew / likeliold
809 0 : if (oddsRatio .gt. 1.0_dp) oddsSwitch1 = .true.
810 : end if
811 :
812 : ! (C) Accept or Reject?
813 13 : If (oddsSwitch1) then
814 :
815 : ! positive accept
816 7098 : Ipos(1) = Ipos(1) + 1_i4
817 28392 : paraold = paranew
818 7098 : likeliold = likelinew
819 28392 : where (changePara)
820 : stepsize = stepsize * accMult
821 : end where
822 7098 : if (likelinew .gt. likelibest) then
823 0 : parabest = paranew
824 0 : likeliold = likelinew
825 0 : likelibest = likelinew
826 : !
827 0 : write(*, *) ''
828 0 : write(*, *) 'best para changed: ', paranew
829 0 : write(*, *) 'likelihood new: ', likelinew
830 0 : write(*, *) ''
831 : end if
832 28392 : burnin_paras_part(Ipos(1) + Ineg(1), :) = paranew(:)
833 :
834 : else
835 :
836 8454030 : call xor4096(seeds(1, 1), RN1, save_state = save_state_1(1, :))
837 31902 : oddsSwitch2 = .false.
838 31902 : if (loglike) then
839 31902 : if (oddsRatio .lt. -700.0_dp) oddsRatio = -700.0_dp ! to avoid underflow
840 31902 : if (rn1 .lt. exp(oddsRatio)) oddsSwitch2 = .true.
841 : else
842 0 : if (rn1 .lt. oddsRatio) oddsSwitch2 = .true.
843 : end if
844 :
845 : If (oddsSwitch2) then
846 :
847 : ! negative accept
848 7213 : Ineg(1) = Ineg(1) + 1_i4
849 28852 : paraold = paranew
850 7213 : likeliold = likelinew
851 28852 : where (changePara)
852 : stepsize = stepsize * accMult
853 : end where
854 28852 : burnin_paras_part(Ipos(1) + Ineg(1), :) = paranew(:)
855 :
856 : else
857 :
858 : ! reject
859 98756 : where (changePara)
860 : stepsize = stepsize * rejMult
861 : end where
862 :
863 : end if
864 :
865 : end if
866 :
867 : end do markovchain
868 :
869 13 : accRatio(1) = real(Ipos(1) + Ineg(1), dp) / real(iter_burnin, dp)
870 13 : if (printflag) then
871 13 : write(str, '(A,I03,A)') '(A7,I4,A15,F5.3,A17,', size(para, 1), '(E9.2,1X),A1)'
872 13 : write(*, str) 'trial #', trial, ' acc_ratio = ', accRatio(1), ' (stepsize = ', stepsize, ')'
873 : end if
874 :
875 13 : if (Ipos(1) + Ineg(1) .gt. 0_i4) then
876 13 : call append(burnin_paras, burnin_paras_part(1 : Ipos(1) + Ineg(1), :))
877 : end if
878 :
879 : ! adjust acceptance multiplier
880 13 : if (accRatio(1) .lt. 0.234_dp) accMult = accMult * 0.99_dp
881 13 : if (accRatio(1) .gt. 0.441_dp) accMult = accMult * 1.01_dp
882 :
883 : ! store good accRatios in history and delete complete history if bad one appears
884 13 : if (accRatio(1) .lt. 0.234_dp .or. accRatio(1) .gt. 0.441_dp) then
885 3 : if(allocated(history_accRatio)) deallocate(history_accRatio)
886 : else
887 : call append(history_accRatio, accRatio(1))
888 : end if
889 :
890 : ! if in history more than 10 values, check for mean and variance
891 13 : if (allocated(history_accRatio)) then
892 10 : accRatio_n = size(history_accRatio, 1)
893 10 : if (accRatio_n .ge. 10_i4) then
894 1 : idummy = accRatio_n - 9_i4
895 1 : accRatio_stddev = stddev(history_accRatio(idummy : accRatio_n))
896 :
897 : ! Check of Convergence
898 1 : if ((accRatio_stddev .lt. Sqrt(1._dp / 12._dp * 0.05_dp**2))) then
899 1 : iStop = .true.
900 1 : if (printflag) then
901 1 : write(*, *) ''
902 1 : write(*, *) 'STOP BURN-IN with accaptence ratio of ', history_accRatio(accRatio_n)
903 1 : write(*, *) 'final stepsize: ', stepsize
904 1 : write(*, *) 'best parameter set found: ', parabest
905 1 : write(*, *) 'with likelihood: ', likelibest
906 : end if
907 : end if
908 :
909 : end if
910 : end if
911 :
912 13 : trial = trial + 1_i4
913 :
914 : end do convergedBURNIN
915 :
916 : ! end do betterParaFound
917 :
918 : end if ! no burn-in skip
919 : !
920 : ! ------------------------------------
921 : ! start initializing things for MCMC, i.e. initialization and allocation
922 : ! necessary for MCMC restart file
923 : ! ------------------------------------
924 :
925 : ! allocate arrays which will be used later
926 5 : allocate(mcmc_paras_3d(iter_mcmc, size(para), chains))
927 45021 : mcmc_paras_3d = -9999.0_dp
928 :
929 1 : vDotDot = -9999.0_dp
930 1 : B = -9999.0_dp
931 1 : W = -9999.0_dp
932 :
933 : ! just to be sure that all chains start with same initial parameter set
934 : ! in both parallel and sequential mode
935 : ! (although in between a new parabest will be found in chains)
936 4 : initial_paraset_mcmc = parabest
937 :
938 6 : Ipos(:) = 0_i4 ! positive accepted
939 6 : Ineg(:) = 0_i4 ! negative accepted
940 :
941 : ! if all parameters converged: Sqrt(R_i) < 1.1 (see Gelman et. al: Baysian Data Analysis, p. 331ff)
942 1 : converged = .False.
943 :
944 : ! transfer all array-like variables in namelist to fixed-size dummy-arrays
945 1 : dummy_maskpara = .false.
946 4 : dummy_maskpara(1 : size(para, 1)) = maskpara
947 1001 : dummy_stepsize = -9999.0_dp
948 4 : dummy_stepsize(1 : size(para, 1)) = stepsize
949 1001 : dummy_parabest = -9999.0_dp
950 4 : dummy_parabest(1 : size(para, 1)) = parabest
951 1001 : dummy_initial_paraset_mcmc = -9999.0_dp
952 4 : dummy_initial_paraset_mcmc(1 : size(para, 1)) = initial_paraset_mcmc
953 :
954 : ! write restart
955 1 : open(999, file = isrestart_file, status = 'unknown', action = 'write', delim = 'QUOTE')
956 1 : write(999, restartnml1)
957 1 : write(999, restartnml2)
958 1 : close(999)
959 :
960 : else ! --> here starts the restart case
961 :
962 : ! read 1st namelist with allocated/scalar variables
963 1 : open(999, file = isrestart_file, status = 'old', action = 'read', delim = 'QUOTE')
964 1 : read(999, nml = restartnml1)
965 1 : close(999)
966 :
967 : ! transfer all array-like variables in namelist to fixed-size dummy-arrays
968 4 : maskpara = dummy_maskpara(1 : size(para, 1))
969 4 : stepsize = dummy_stepsize(1 : size(para, 1))
970 4 : parabest = dummy_parabest(1 : size(para, 1))
971 4 : initial_paraset_mcmc = dummy_initial_paraset_mcmc(1 : size(para, 1))
972 :
973 : ! allocate global arrays
974 6 : allocate(truepara(count(maskpara)))
975 3 : allocate(seeds(chains, 3))
976 3 : allocate(save_state_1(chains, n_save_state))
977 2 : allocate(save_state_2(chains, n_save_state))
978 2 : allocate(save_state_3(chains, n_save_state))
979 6 : allocate(Ipos(chains), Ineg(chains), accRatio(chains))
980 3 : allocate(vDotJ(chains), s2(chains))
981 4 : allocate(sqrtR(size(para)))
982 1 : if (.not. converged) then
983 0 : if (present(iter_mcmc_in)) then
984 0 : allocate(mcmc_paras_3d(iter_mcmc - iter_mcmc_in, size(para), chains))
985 : else
986 0 : allocate(mcmc_paras_3d(iter_mcmc - 1000_i4 * n, size(para), chains))
987 : end if
988 : else
989 5 : allocate(mcmc_paras_3d(iter_mcmc, size(para), chains))
990 : end if
991 :
992 : end if
993 :
994 : !----------------------------------------------------------------------
995 : ! (2) MCMC
996 : !----------------------------------------------------------------------
997 :
998 2 : if (printflag) then
999 2 : write(*, *) ''
1000 2 : write(*, *) '--------------------------------------------------'
1001 2 : write(*, *) 'Starting MCMC (chains = ', chains, ', iter_mcmc = ', iter_mcmc, ')'
1002 2 : write(*, *) '--------------------------------------------------'
1003 2 : write(*, *) ''
1004 : end if
1005 :
1006 3 : convergedMCMC : do while (.not. converged)
1007 : ! read restart
1008 1 : open(999, file = isrestart_file, status = 'unknown', action = 'read', delim = 'QUOTE')
1009 1 : read(999, restartnml1)
1010 1 : read(999, restartnml2)
1011 1 : close(999)
1012 :
1013 : ! transfer all array-like variables in namelist to fixed-size dummy-arrays
1014 4 : maskpara = dummy_maskpara(1 : size(para, 1))
1015 4 : stepsize = dummy_stepsize(1 : size(para, 1))
1016 4 : parabest = dummy_parabest(1 : size(para, 1))
1017 4 : initial_paraset_mcmc = dummy_initial_paraset_mcmc(1 : size(para, 1))
1018 :
1019 : ! resize mcmc_paras_3d
1020 : ! iter_mcmc was increased --> indicates new length of mcmc_paras_3d
1021 : ! Minval(Ipos+Ineg) --> indicates old length of mcmc_paras_3d
1022 6 : idummy = minval(Ipos + Ineg)
1023 1 : if ((iter_mcmc - idummy) > 0) then
1024 4 : allocate(tmp(idummy, size(para), chains))
1025 21 : tmp(:, :, :) = mcmc_paras_3d(1 : idummy, :, :)
1026 1 : deallocate(mcmc_paras_3d)
1027 5 : allocate(mcmc_paras_3d(iter_mcmc, size(para), chains))
1028 21 : mcmc_paras_3d(1 : idummy, :, :) = tmp(:, :, :)
1029 1 : deallocate(tmp)
1030 : end if
1031 :
1032 : !$OMP parallel default(shared) &
1033 : !$OMP private(chain, paraold, paranew, likeliold, likelinew, oddsSwitch1, oddsSwitch2, RN1, oddsRatio, ChangePara)
1034 : !$OMP do
1035 6 : parallelchain : do chain = 1, chains
1036 :
1037 5 : if (Ipos(chain) + Ineg(chain) .eq. 0_i4) then
1038 20 : paraold = initial_paraset_mcmc ! = parabest of burn-in
1039 : else
1040 0 : paraold = mcmc_paras_3d(Ipos(chain) + Ineg(chain), :, chain)
1041 : end if
1042 5 : likeliold = likelihood(paraold, eval)
1043 :
1044 1 : markovchainMCMC : do
1045 :
1046 : ! (A) Generate new parameter set
1047 191024 : ChangePara = .false.
1048 191024 : paranew = paraold
1049 :
1050 : call GenerateNewParameterset_dp(ParaSelectMode, paraold, truepara, rangePar, stepsize, &
1051 143268 : save_state_2(chain, :), &
1052 143268 : save_state_3(chain, :), &
1053 25549460 : paranew, ChangePara)
1054 :
1055 : ! (B) new likelihood
1056 47756 : likelinew = likelihood(paranew, eval)
1057 47756 : oddsSwitch1 = .false.
1058 47756 : if (loglike) then
1059 47756 : oddsRatio = likelinew - likeliold
1060 47756 : if (oddsRatio .gt. 0.0_dp) oddsSwitch1 = .true.
1061 : else
1062 0 : oddsRatio = likelinew / likeliold
1063 0 : if (oddsRatio .gt. 1.0_dp) oddsSwitch1 = .true.
1064 : end if
1065 :
1066 : ! (C) Accept or Reject?
1067 : If (oddsSwitch1) then
1068 :
1069 : ! positive accept
1070 7508 : Ipos(chain) = Ipos(chain) + 1_i4
1071 30032 : paraold = paranew
1072 7508 : likeliold = likelinew
1073 30032 : mcmc_paras_3d(Ipos(chain) + Ineg(chain), :, chain) = paranew(:)
1074 :
1075 7508 : if (printflag) then
1076 7508 : if ((Ipos(chain) + Ineg(chain) .gt. 0_i4) .and. &
1077 : (mod(Ipos(chain) + Ineg(chain), 10000_i4)) .eq. 0) then
1078 0 : write(*, *) 'Chain ', chain, ': Done ', Ipos(chain) + Ineg(chain), ' samples ...'
1079 : end if
1080 : end if
1081 :
1082 : ! If the following code block is not commented
1083 : ! then mcmc can be used as an optimiser as well
1084 : ! and not 'only' for determination of parameter uncertainties.
1085 : ! However, then the serial and the parallel versions give different results.
1086 : ! if (likelinew .gt. likelibest) then
1087 : ! parabest = paranew
1088 : ! likelibest = likelinew
1089 : ! if (printflag) then
1090 : ! write(*,*) ''
1091 : ! write(*,*) 'best para changed: ',paranew
1092 : ! write(*,*) 'likelihood new: ',likelinew
1093 : ! write(*,*) ''
1094 : ! end if
1095 : ! end if
1096 :
1097 : else
1098 :
1099 10665720 : call xor4096(seeds(chain, 1), RN1, save_state = save_state_1(chain, :))
1100 40248 : oddsSwitch2 = .false.
1101 40248 : if (loglike) then
1102 40248 : if (rn1 .lt. exp(oddsRatio)) oddsSwitch2 = .true.
1103 : else
1104 0 : if (rn1 .lt. oddsRatio) oddsSwitch2 = .true.
1105 : end if
1106 :
1107 : If (oddsSwitch2) then
1108 :
1109 : ! negative accept
1110 7492 : Ineg(chain) = Ineg(chain) + 1_i4
1111 29968 : paraold = paranew
1112 7492 : likeliold = likelinew
1113 29968 : mcmc_paras_3d(Ipos(chain) + Ineg(chain), :, chain) = paranew(:)
1114 :
1115 7492 : if (printflag) then
1116 7492 : if ((Ipos(chain) + Ineg(chain) .gt. 0_i4) .and. &
1117 : (mod(Ipos(chain) + Ineg(chain), 10000_i4)) .eq. 0) then
1118 0 : write(*, *) 'Chain ', chain, ': Done ', Ipos(chain) + Ineg(chain), ' samples ...'
1119 : end if
1120 : end if
1121 :
1122 : end if
1123 :
1124 : end if
1125 :
1126 : ! enough samples were found
1127 47756 : if ((Ipos(chain) + Ineg(chain) .gt. 0_i4) .and. &
1128 5 : (mod(Ipos(chain) + Ineg(chain), iter_mcmc) .eq. 0_i4)) exit
1129 :
1130 : end do markovchainMCMC
1131 :
1132 : end do parallelchain
1133 : !$OMP end do
1134 : !$OMP end parallel
1135 :
1136 : #ifdef FORCES_WITH_NETCDF
1137 : ! write parameter sets to temporal file
1138 1 : if (itmp_file) then
1139 : ! splitting into path and filename
1140 1 : slash_pos = index(istmp_file, '/', .true.)
1141 1 : len_filename = len_trim(istmp_file)
1142 1 : path = istmp_file(1 : slash_pos)
1143 1 : filename = istmp_file(slash_pos + 1 : len_filename)
1144 : !
1145 6 : do chain = 1, chains
1146 5 : write(str, *) chain
1147 5 : write(outputfile, *) trim(adjustl(path)), trim(adjustl(str)), '_', trim(adjustl(filename))
1148 6 : if (present(iter_mcmc_in)) then
1149 0 : allocate(tmp(iter_mcmc_in, size(para, 1), 1))
1150 0 : tmp(:, :, 1) = mcmc_paras_3d(iter_mcmc - iter_mcmc_in + 1_i4 : iter_mcmc, :, chain)
1151 0 : if (iter_mcmc .ne. iter_mcmc_in) then
1152 : ! append
1153 0 : call dump_netcdf(trim(adjustl(outputfile)), tmp, append = .true.)
1154 : else
1155 : ! first time of writing
1156 0 : call dump_netcdf(trim(adjustl(outputfile)), tmp)
1157 : end if
1158 0 : deallocate(tmp)
1159 : else
1160 20 : allocate(tmp(1000_i4 * n, size(para, 1), 1))
1161 45020 : tmp(:, :, 1) = mcmc_paras_3d(iter_mcmc - (1000_i4 * n) + 1_i4 : iter_mcmc, :, chain)
1162 5 : if (iter_mcmc .ne. 1000_i4 * n) then
1163 : ! append
1164 0 : call dump_netcdf(trim(adjustl(outputfile)), tmp, append = .true.)
1165 : else
1166 : ! first time of writing
1167 5 : call dump_netcdf(trim(adjustl(outputfile)), tmp)
1168 : end if
1169 5 : deallocate(tmp)
1170 : end if
1171 : end do
1172 : end if
1173 : #else
1174 : if (itmp_file) &
1175 : call error_message("MCMC: FORCES was compiled without NetCDF support but you set 'tmp_file'")
1176 : #endif
1177 :
1178 : ! test for convergence: Gelman et. al: Baysian Data Analysis, p. 331ff
1179 : ! sqrt(R) = sqrt( ( (n-1)/n * W + 1/n * B ) / W ) < 1.1 for each parameter
1180 : ! n ... last half of the sequence
1181 : ! W ... within sequence variance
1182 : ! B ... between chain variance
1183 :
1184 1 : if (printflag) then
1185 1 : write(*, *) ' '
1186 1 : write(*, *) 'Checking for convergence ....'
1187 : end if
1188 4 : sqrtR = 0.0_dp
1189 6 : n_end = minval(Ipos + Ineg)
1190 1 : n_start = n_end / 2_i4 + 1_i4
1191 :
1192 4 : do iPar = 1, size(truepara)
1193 :
1194 : ! Between chain variances
1195 18 : do chain = 1, chains
1196 15 : vDotJ(chain) = 1.0_dp / real(n_end - n_start + 1_i4, dp) * &
1197 22533 : sum(mcmc_paras_3d(n_start : n_end, truepara(iPar), chain))
1198 : end do
1199 18 : vDotDot = 1.0_dp / real(chains, dp) * sum(vDotJ)
1200 : B = real(n_end - n_start + 1_i4, dp) / real(chains - 1_i4, dp) * &
1201 18 : sum((vDotJ - vDotDot) * (vDotJ - vDotDot))
1202 :
1203 : ! Within chain variances
1204 18 : do chain = 1, chains
1205 15 : s2(chain) = 1.0_dp / real(n_end - n_start, dp) * &
1206 22533 : sum((mcmc_paras_3d(n_start : n_end, truepara(iPar), chain) - vDotJ(chain))**2)
1207 : end do
1208 18 : W = 1.0_dp / real(chains, dp) * Sum(s2)
1209 :
1210 : ! ratio sqrtR
1211 4 : if ((w .lt. tiny(1.0_dp)) .and. (b .lt. tiny(1.0_dp))) then
1212 : ! Mathematica says that this is the limit, if w and b both go to zero
1213 0 : sqrtR(truepara(iPar)) = sqrt(real(n_end - n_start, dp) / real(n_end - n_start + 1_i4, dp))
1214 : else
1215 3 : sqrtR(truepara(iPar)) = real(n_end - n_start, dp) / real(n_end - n_start + 1_i4, dp) * W + &
1216 3 : 1.0_dp / real(n_end - n_start + 1_i4, dp) * B
1217 3 : sqrtR(truepara(iPar)) = sqrt(sqrtR(truepara(iPar)) / W)
1218 : end if
1219 : end do
1220 1 : if (printflag) then
1221 4 : do i = 1, size(para)
1222 4 : if (sqrtR(i) .gt. 1.1_dp) then
1223 0 : write(*, *) ' sqrtR para #', i, ' : ', sqrtR(i), ' <-- FAILED'
1224 : else
1225 3 : write(*, *) ' sqrtR para #', i, ' : ', sqrtR(i)
1226 : end if
1227 : end do
1228 : end if
1229 4 : if (all(sqrtR .lt. 1.1_dp)) then
1230 1 : converged = .true.
1231 1 : if (printflag) then
1232 1 : write(*, *) ' --> converged (all less than 1.1)'
1233 1 : write(*, *) ' '
1234 : end if
1235 : else
1236 : ! increase number of iterations
1237 0 : if (present(iter_mcmc_in)) then
1238 0 : iter_mcmc = iter_mcmc + iter_mcmc_in
1239 : else
1240 0 : iter_mcmc = iter_mcmc + 1000_i4 * n
1241 : end if
1242 :
1243 0 : if (printflag) then
1244 0 : write(*, *) ' --> not converged (not all less than 1.1)'
1245 0 : write(*, *) ' increasing iterations to ', iter_mcmc
1246 0 : write(*, *) ' '
1247 : end if
1248 : end if
1249 :
1250 : ! transfer all array-like variables in namelist to fixed-size dummy-arrays
1251 1 : dummy_maskpara = .false.
1252 4 : dummy_maskpara(1 : size(para, 1)) = maskpara
1253 1001 : dummy_stepsize = -9999.0_dp
1254 4 : dummy_stepsize(1 : size(para, 1)) = stepsize
1255 1001 : dummy_parabest = -9999.0_dp
1256 4 : dummy_parabest(1 : size(para, 1)) = parabest
1257 1001 : dummy_initial_paraset_mcmc = -9999.0_dp
1258 4 : dummy_initial_paraset_mcmc(1 : size(para, 1)) = initial_paraset_mcmc
1259 :
1260 : ! write restart
1261 1 : open(999, file = isrestart_file, status = 'unknown', action = 'write', delim = 'QUOTE')
1262 1 : write(999, restartnml1)
1263 1 : write(999, restartnml2)
1264 3 : close(999)
1265 :
1266 : end do convergedMCMC
1267 :
1268 : ! read restart
1269 2 : open(999, file = isrestart_file, status = 'unknown', action = 'read', delim = 'QUOTE')
1270 2 : read(999, restartnml1)
1271 2 : read(999, restartnml2)
1272 2 : close(999)
1273 :
1274 : ! transfer all array-like variables in namelist to fixed-size dummy-arrays
1275 8 : maskpara = dummy_maskpara(1 : size(para, 1))
1276 8 : stepsize = dummy_stepsize(1 : size(para, 1))
1277 8 : parabest = dummy_parabest(1 : size(para, 1))
1278 8 : initial_paraset_mcmc = dummy_initial_paraset_mcmc(1 : size(para, 1))
1279 :
1280 : ! reshape of mcmc_paras_3d: return only 2d matrix mcmc_paras
1281 8 : allocate(mcmc_paras(size(mcmc_paras_3d, 1) * size(mcmc_paras_3d, 3), size(mcmc_paras_3d, 2)))
1282 12 : do chain = 1, chains
1283 90042 : mcmc_paras((chain - 1_i4) * size(mcmc_paras_3d, 1) + 1_i4 : chain * size(mcmc_paras_3d, 1), :) = mcmc_paras_3d(:, :, chain)
1284 : end do
1285 :
1286 2 : RETURN
1287 6 : END SUBROUTINE mcmc_dp
1288 :
1289 2 : SUBROUTINE mcmc_stddev_dp(eval, likelihood, para, rangePar, & ! obligatory IN
1290 : mcmc_paras, burnin_paras, & ! obligatory OUT
1291 1 : seed_in, printflag_in, maskpara_in, & ! optional IN
1292 : tmp_file, & ! optional IN : filename for temporal output of
1293 : ! ! MCMC parasets
1294 : loglike_in, & ! optional IN : true if loglikelihood is given
1295 : ParaSelectMode_in, & ! optional IN : (=1) half, (=2) one, (=3) all
1296 : iter_burnin_in, & ! optional IN : markov length of (1) burn-in
1297 : iter_mcmc_in, & ! optional IN : markov length of (2) mcmc
1298 : chains_in, & ! optional IN : number of parallel chains of MCMC
1299 0 : stepsize_in) ! optional_IN : stepsize for each param. (no burn-in)
1300 :
1301 : IMPLICIT NONE
1302 :
1303 : procedure(eval_interface), INTENT(IN), POINTER :: eval
1304 : procedure(objective_interface), intent(in), pointer :: likelihood
1305 :
1306 : REAL(DP), DIMENSION(:, :), INTENT(IN) :: rangePar ! range for each parameter
1307 : REAL(DP), DIMENSION(:), INTENT(IN) :: para ! initial parameter i
1308 : INTEGER(I8), OPTIONAL, INTENT(IN) :: seed_in ! Seeds of random numbers: dim1=chains, dim2=3
1309 : ! ! (DEFAULT: Get_timeseed)
1310 : LOGICAL, OPTIONAL, INTENT(IN) :: printflag_in ! If command line output is written (.true.)
1311 : ! ! (DEFAULT: .false.)
1312 : LOGICAL, OPTIONAL, DIMENSION(size(para, 1)), &
1313 : INTENT(IN) :: maskpara_in ! true if parameter will be optimized
1314 : ! ! false if parameter is discarded in optimization
1315 : ! ! DEFAULT = .true.
1316 : LOGICAL, OPTIONAL, INTENT(IN) :: loglike_in ! true if loglikelihood is given instead of likelihood
1317 : ! ! DEFAULT: .false.
1318 : INTEGER(I4), OPTIONAL, INTENT(IN) :: ParaSelectMode_in ! how many parameters changed per iteration:
1319 : ! ! 1: half of the parameters
1320 : ! ! 2: only one (DEFAULT)
1321 : ! ! 3: all
1322 : INTEGER(I4), OPTIONAL, INTENT(IN) :: iter_burnin_in ! # iterations before checking ac_ratio
1323 : ! ! DEFAULT= MIN(250, 20_i4*n)
1324 : INTEGER(I4), OPTIONAL, INTENT(IN) :: iter_mcmc_in ! # iterations per chain for posterior sampling
1325 : ! ! DEFAULT= 1000_i4*n
1326 : INTEGER(I4), OPTIONAL, INTENT(IN) :: chains_in ! number of parallel mcmc chains
1327 : ! ! DEFAULT= 5_i4
1328 : REAL(DP), OPTIONAL, DIMENSION(size(para, 1)), &
1329 : INTENT(IN) :: stepsize_in ! stepsize for each parameter
1330 : ! ! if given burn-in is discarded
1331 : CHARACTER(len = *), OPTIONAL, INTENT(IN) :: tmp_file ! filename for temporal data saving: every iter_mcmc_in
1332 : ! ! iterations parameter sets are appended to this file
1333 : ! ! the number of the chain will be prepended to filename
1334 : ! ! DEFAULT= no file writing
1335 : REAL(DP), DIMENSION(:, :), ALLOCATABLE, INTENT(OUT) :: mcmc_paras
1336 : ! ! array to save para values of MCMC runs,
1337 : ! ! dim1=sets of all chains, dim2=paras
1338 : REAL(DP), DIMENSION(:, :), ALLOCATABLE, INTENT(OUT) :: burnin_paras
1339 : ! ! array to save para values of Burn in run
1340 :
1341 :
1342 : ! local variables
1343 : INTEGER(I4) :: n ! Number of parameters
1344 : LOGICAL :: printflag ! If command line output is written
1345 0 : LOGICAL, DIMENSION(size(para, 1)) :: maskpara ! true if parameter will be optimized
1346 1 : INTEGER(I4), DIMENSION(:), ALLOCATABLE :: truepara ! indexes of parameters to be optimized
1347 : LOGICAL :: loglike ! if loglikelihood is given
1348 :
1349 : ! for random numbers
1350 1 : INTEGER(I8), dimension(:, :), allocatable :: seeds ! Seeds of random numbers: dim1=chains, dim2=3
1351 1 : REAL(DP) :: RN1, RN2, RN3 ! Random numbers
1352 1 : integer(I8), dimension(:, :), allocatable :: save_state_1 ! optional argument for restarting RN stream 1:
1353 1 : integer(I8), dimension(:, :), allocatable :: save_state_2 ! optional argument for restarting RN stream 2:
1354 1 : integer(I8), dimension(:, :), allocatable :: save_state_3 ! optional argument for restarting RN stream 3:
1355 : ! ! dim1=chains, dim2=n_save_state
1356 :
1357 : ! Dummies
1358 1 : REAL(DP), DIMENSION(:, :, :), ALLOCATABLE :: tmp
1359 : integer(I4) :: idummy
1360 : logical :: oddsSwitch1, oddsSwitch2
1361 : character(100) :: str
1362 : character(200) :: outputfile
1363 : integer(i4) :: slash_pos
1364 : integer(i4) :: len_filename
1365 : character(200) :: filename
1366 : character(200) :: path
1367 :
1368 : ! FOR BURN-IN AND MCMC
1369 1 : REAL(DP), DIMENSION(:, :, :), ALLOCATABLE :: mcmc_paras_3d ! array to save para values of MCMC runs,
1370 : ! ! dim1=sets, dim2=paras, dim3=chain
1371 : integer(I4) :: i
1372 2 : integer(I4), dimension(:), allocatable :: Ipos, Ineg
1373 : logical :: iStop
1374 : INTEGER(I4) :: ParaSelectMode ! how many parameters are changed per jump
1375 : INTEGER(I4) :: iter_burnin ! number fo iterations before checking ac_ratio
1376 : INTEGER(I4) :: iter_mcmc ! number fo iterations for posterior sampling
1377 : INTEGER(I4) :: chains ! number of parallel mcmc chains
1378 : INTEGER(I4) :: chain ! counter for chains, 1...chains
1379 1 : REAL(DP) :: accMult ! acceptance multiplier for stepsize
1380 1 : REAL(DP) :: rejMult ! rejection multiplier for stepsize
1381 0 : LOGICAL, DIMENSION(size(para, 1)) :: ChangePara ! logical array to switch if parameter is changed
1382 : INTEGER(I4) :: trial ! number of trials for burn-in
1383 4 : REAL(DP), DIMENSION(size(para, 1)) :: stepsize ! stepsize adjusted by burn-in and used by mcmc
1384 4 : REAL(DP), DIMENSION(size(para, 1)) :: paraold ! old parameter set
1385 5 : REAL(DP), DIMENSION(size(para, 1)) :: paranew ! new parameter set
1386 5 : REAL(DP), DIMENSION(size(para, 1)) :: parabest ! best parameter set overall
1387 5 : REAL(DP), DIMENSION(size(para, 1)) :: initial_paraset_mcmc ! best parameterset found in burn-in
1388 1 : REAL(DP) :: stddev_data ! approximated stddev of data for best paraset found
1389 : LOGICAL :: parabestChanged ! if better parameter set was found during burn-in
1390 1 : REAL(DP) :: likeliold ! likelihood of old parameter set
1391 1 : REAL(DP) :: likelinew ! likelihood of new parameter set
1392 1 : REAL(DP) :: likelibest ! likelihood of best parameter set overall
1393 1 : REAL(DP) :: stddev_new ! standard deviation of errors with current parameter set
1394 1 : REAL(DP) :: likeli_new ! likelihood using stddev_new instead of stddev_data
1395 : INTEGER(I4) :: markov ! counter for markov chain
1396 1 : REAL(DP), DIMENSION(:, :), ALLOCATABLE :: burnin_paras_part ! accepted parameter sets of one MC in burn-in
1397 1 : REAL(DP) :: oddsRatio ! ratio of likelihoods = likelinew/likeliold
1398 1 : REAL(DP), dimension(:), allocatable :: accRatio ! acceptance ratio = (pos/neg) accepted/iter_burnin
1399 1 : REAL(DP) :: accratio_stddev ! stddev of accRatios in history
1400 1 : REAL(DP), DIMENSION(:), ALLOCATABLE :: history_accratio ! history of 'good' acceptance ratios
1401 : INTEGER(I4) :: accratio_n ! number of 'good' acceptance ratios
1402 : !REAL(DP), DIMENSION(:,:), ALLOCATABLE :: history_stepsize ! history of 'good' stepsizes
1403 :
1404 : ! for checking convergence of MCMC
1405 1 : real(dp), allocatable, dimension(:) :: sqrtR
1406 1 : real(dp), allocatable, dimension(:) :: vDotJ
1407 1 : real(dp), allocatable, dimension(:) :: s2
1408 1 : real(dp) :: vDotDot
1409 1 : real(dp) :: B
1410 1 : real(dp) :: W
1411 : integer(i4) :: n_start, n_end, iPar
1412 : LOGICAL :: converged ! if MCMC already converged
1413 :
1414 : ! for OMP
1415 : !$ integer(i4) :: n_threads
1416 :
1417 : ! CHECKING OPTIONALS
1418 :
1419 1 : if (present(loglike_in)) then
1420 1 : loglike = loglike_in
1421 : else
1422 : loglike = .false.
1423 : end if
1424 :
1425 1 : if (present(maskpara_in)) then
1426 4 : if (count(maskpara_in) .eq. 0_i4) then
1427 0 : stop 'Input argument maskpara: At least one element has to be true'
1428 : else
1429 4 : maskpara = maskpara_in
1430 : end if
1431 : else
1432 0 : maskpara = .true.
1433 : end if
1434 :
1435 6 : allocate (truepara(count(maskpara)))
1436 1 : idummy = 0_i4
1437 4 : do i = 1, size(para, 1)
1438 4 : if (maskpara(i)) then
1439 3 : idummy = idummy + 1_i4
1440 3 : truepara(idummy) = i
1441 : end if
1442 : end do
1443 :
1444 1 : n = size(truepara, 1)
1445 :
1446 1 : if (present(ParaSelectMode_in)) then
1447 1 : ParaSelectMode = ParaSelectMode_in
1448 : else
1449 : ! change only one parameter per jump
1450 0 : ParaSelectMode = 2_i4
1451 : end if
1452 :
1453 : ! after how many iterations do we compute ac_ratio???
1454 1 : if (present(iter_burnin_in)) then
1455 0 : if (iter_burnin_in .le. 0_i4) then
1456 0 : stop 'Input argument iter_burn_in must be greater than 0!'
1457 : else
1458 0 : iter_burnin = iter_burnin_in
1459 : end if
1460 : else
1461 1 : iter_burnin = Max(250_i4, 1000_i4 * n)
1462 : end if
1463 :
1464 : ! how many iterations ('jumps') are performed in MCMC
1465 : ! iter_mcmc_in is handled later properly (after acceptance ratio of burn_in is known)
1466 1 : if (present(iter_mcmc_in)) then
1467 0 : if (iter_mcmc_in .le. 0_i4) then
1468 0 : stop 'Input argument iter_mcmc must be greater than 0!'
1469 : else
1470 0 : iter_mcmc = iter_mcmc_in
1471 : end if
1472 : else
1473 1 : iter_mcmc = 1000_i4 * n
1474 : end if
1475 :
1476 1 : if (present(chains_in)) then
1477 0 : if (chains_in .lt. 2_i4) then
1478 0 : stop 'Input argument chains must be at least 2!'
1479 : end if
1480 0 : chains = chains_in
1481 : else
1482 1 : chains = 5_i4
1483 : end if
1484 :
1485 1 : if (present(stepsize_in)) then
1486 0 : stepsize = stepsize_in
1487 : end if
1488 :
1489 1 : if (present(printflag_in)) then
1490 1 : printflag = printflag_in
1491 : else
1492 : printflag = .false.
1493 : end if
1494 :
1495 : !$ write(*,*) '--------------------------------------------------'
1496 : !$ write(*,*) ' This program is parallel.'
1497 : !$OMP parallel
1498 : !$ n_threads = OMP_GET_NUM_THREADS()
1499 : !$OMP end parallel
1500 : !$ write(*,*) ' ',chains,' MCMC chains will run in ',n_threads,' threads'
1501 : !$ write(*,*) '--------------------------------------------------'
1502 :
1503 1 : if (printflag) then
1504 1 : write(*, *) 'Following parameters will be sampled with MCMC: ', truepara
1505 : end if
1506 :
1507 3 : allocate(seeds(chains, 3))
1508 3 : allocate(save_state_1(chains, n_save_state))
1509 2 : allocate(save_state_2(chains, n_save_state))
1510 2 : allocate(save_state_3(chains, n_save_state))
1511 6 : allocate(Ipos(chains), Ineg(chains), accRatio(chains))
1512 3 : allocate(vDotJ(chains), s2(chains))
1513 4 : allocate(sqrtR(size(para)))
1514 :
1515 1 : if (present(seed_in)) then
1516 4 : seeds(1, :) = (/ seed_in, seed_in + 1000_i8, seed_in + 2000_i8 /)
1517 : else
1518 : ! Seeds depend on actual time
1519 0 : call get_timeseed(seeds(1, :))
1520 : end if
1521 5 : do chain = 2, chains
1522 17 : seeds(chain, :) = seeds(chain - 1_i4, :) + 3000_i8
1523 : end do
1524 :
1525 6 : do chain = 1, chains
1526 1325 : call xor4096(seeds(chain, 1), RN1, save_state = save_state_1(chain, :))
1527 1325 : call xor4096(seeds(chain, 2), RN2, save_state = save_state_2(chain, :))
1528 1326 : call xor4096g(seeds(chain, 3), RN3, save_state = save_state_3(chain, :))
1529 : end do
1530 19 : seeds = 0_i8
1531 :
1532 4 : parabest = para
1533 : ! write(*,*) parabest
1534 :
1535 : ! initialize likelihood and sigma
1536 1 : likelibest = likelihood(parabest, eval, 1.0_dp, stddev_new, likeli_new)
1537 1 : likelibest = likeli_new
1538 1 : stddev_data = stddev_new
1539 :
1540 : !----------------------------------------------------------------------
1541 : ! (1) BURN IN
1542 : !----------------------------------------------------------------------
1543 :
1544 1 : if (.not. present(stepsize_in)) then
1545 1 : if (printflag) then
1546 1 : write(*, *) ''
1547 1 : write(*, *) '--------------------------------------------------'
1548 1 : write(*, *) 'Starting Burn-In (iter_burnin = ', iter_burnin, ')'
1549 1 : write(*, *) '--------------------------------------------------'
1550 1 : write(*, *) ''
1551 : end if
1552 :
1553 : ! INITIALIZATION
1554 :
1555 : ! probably too large, but large enough to store values of one markovchain
1556 4 : allocate(burnin_paras_part(iter_burnin, size(para)))
1557 :
1558 1 : parabestChanged = .true.
1559 :
1560 : ! -------------------------------------------------------------------------------------
1561 : ! restarts if a better parameter set was found in between
1562 : ! -------------------------------------------------------------------------------------
1563 2 : betterParaFound : do while (parabestChanged)
1564 :
1565 1 : if (allocated(burnin_paras)) deallocate(burnin_paras)
1566 1 : if (allocated(history_accRatio)) deallocate(history_accRatio)
1567 :
1568 4 : parabestChanged = .false.
1569 4 : stepsize = 1.0_dp
1570 1 : trial = 1_i4
1571 1 : iStop = .false.
1572 1 : accMult = 1.01_dp
1573 1 : rejMult = 0.99_dp
1574 : !stddev_data = stddev_function(parabest)
1575 : !likelibest = likelihood(parabest,stddev_data,stddev_new=stddev_new,likeli_new=likeli_new)
1576 : !likelibest = likeli_new
1577 : !stddev_data = stddev_new
1578 :
1579 1 : if (printflag) then
1580 1 : write(*, *) ' '
1581 1 : write(*, *) 'Restart Burn-In with new approximation of std.dev. of data: '
1582 1 : write(*, *) ' parabest = ', parabest
1583 1 : write(*, *) ' stddev = ', stddev_data
1584 1 : write(*, *) ' likelihood = ', likelibest
1585 1 : write(*, *) ' ssq = ', -2.0_dp * stddev_data**2 * likelibest
1586 1 : write(*, *) ' '
1587 : end if
1588 :
1589 :
1590 : ! ----------------------------------------------------------------------------------
1591 : ! repeats until convergence of acceptance ratio or better parameter set was found
1592 : ! ----------------------------------------------------------------------------------
1593 15 : convergedBURNIN : do while ((.not. iStop) .and. (.not. parabestChanged))
1594 :
1595 78 : Ipos = 0_i4 ! positive accepted
1596 78 : Ineg = 0_i4 ! negative accepted
1597 52 : paraold = parabest
1598 13 : likeliold = likelibest !likelihood(paraold,stddev_data)
1599 :
1600 : ! -------------------------------------------------------------------------------
1601 : ! do a short-cut MCMC
1602 : ! -------------------------------------------------------------------------------
1603 39013 : markovchain : do markov = 1, iter_burnin
1604 :
1605 : ! (A) Generate new parameter set
1606 156000 : ChangePara = .false.
1607 156000 : paranew = paraold
1608 : ! using RN from chain #1
1609 : call GenerateNewParameterset_dp(ParaSelectMode, paraold, truepara, rangePar, stepsize, &
1610 78000 : save_state_2(1, :), &
1611 78000 : save_state_3(1, :), &
1612 20787000 : paranew, ChangePara)
1613 :
1614 : ! (B) new likelihood
1615 39000 : likelinew = likelihood(paranew, eval, stddev_data, stddev_new, likeli_new)
1616 :
1617 39000 : oddsSwitch1 = .false.
1618 39000 : if (loglike) then
1619 39000 : oddsRatio = likelinew - likeliold
1620 39000 : if (oddsRatio .gt. 0.0_dp) oddsSwitch1 = .true.
1621 : else
1622 0 : oddsRatio = likelinew / likeliold
1623 0 : if (oddsRatio .gt. 1.0_dp) oddsSwitch1 = .true.
1624 : end if
1625 : ! write(*,*) 'oddsRatio = ',oddsRatio
1626 :
1627 : ! (C) Accept or Reject?
1628 13 : If (oddsSwitch1) then
1629 :
1630 : ! positive accept
1631 7119 : Ipos(1) = Ipos(1) + 1_i4
1632 28476 : paraold = paranew
1633 7119 : likeliold = likelinew
1634 28476 : where (changePara)
1635 : stepsize = stepsize * accMult
1636 : end where
1637 7119 : if (likelinew .gt. likelibest) then
1638 0 : parabest = paranew
1639 0 : likeliold = likeli_new !JM
1640 : ! Here the sigma is reset!
1641 0 : likelibest = likeli_new
1642 0 : stddev_data = stddev_new
1643 : !
1644 0 : parabestChanged = .true.
1645 0 : write(*, *) ''
1646 0 : write(*, *) 'best para changed: ', paranew
1647 0 : write(*, *) 'likelihood new: ', likelinew
1648 0 : write(*, *) ''
1649 : end if
1650 28476 : burnin_paras_part(Ipos(1) + Ineg(1), :) = paranew(:)
1651 : ! write(*,*) 'positive accept'
1652 :
1653 : else
1654 :
1655 8448465 : call xor4096(seeds(1, 1), RN1, save_state = save_state_1(1, :))
1656 31881 : oddsSwitch2 = .false.
1657 31881 : if (loglike) then
1658 31881 : if (oddsRatio .lt. -700.0_dp) oddsRatio = -700.0_dp ! to avoid underflow
1659 31881 : if (rn1 .lt. exp(oddsRatio)) oddsSwitch2 = .true.
1660 : else
1661 0 : if (rn1 .lt. oddsRatio) oddsSwitch2 = .true.
1662 : end if
1663 :
1664 : If (oddsSwitch2) then
1665 :
1666 : ! negative accept
1667 7191 : Ineg(1) = Ineg(1) + 1_i4
1668 28764 : paraold = paranew
1669 7191 : likeliold = likelinew
1670 : ! stddev_data = stddev_new !JM
1671 28764 : where (changePara)
1672 : stepsize = stepsize * accMult
1673 : end where
1674 28764 : burnin_paras_part(Ipos(1) + Ineg(1), :) = paranew(:)
1675 : ! write(*,*) 'negative accept'
1676 :
1677 : else
1678 :
1679 : ! reject
1680 98760 : where (changePara)
1681 : stepsize = stepsize * rejMult
1682 : end where
1683 : ! write(*,*) 'reject'
1684 :
1685 : end if
1686 :
1687 : end if
1688 :
1689 : ! write(*,*) ''
1690 :
1691 : end do markovchain
1692 :
1693 13 : accRatio(1) = real(Ipos(1) + Ineg(1), dp) / real(iter_burnin, dp)
1694 13 : if (printflag) then
1695 13 : write(str, '(A,I03,A)') '(A7,I4,A15,F5.3,A17,', size(para, 1), '(E9.2,1X),A1)'
1696 13 : write(*, str) 'trial #', trial, ' acc_ratio = ', accRatio(1), ' (stepsize = ', stepsize, ')'
1697 : end if
1698 :
1699 13 : if (Ipos(1) + Ineg(1) .gt. 0_i4) then
1700 13 : call append(burnin_paras, burnin_paras_part(1 : Ipos(1) + Ineg(1), :))
1701 : end if
1702 :
1703 : ! adjust acceptance multiplier
1704 13 : if (accRatio(1) .lt. 0.234_dp) accMult = accMult * 0.99_dp
1705 13 : if (accRatio(1) .gt. 0.441_dp) accMult = accMult * 1.01_dp
1706 :
1707 : ! store good accRatios in history and delete complete history if bad one appears
1708 13 : if (accRatio(1) .lt. 0.234_dp .or. accRatio(1) .gt. 0.441_dp) then
1709 3 : if(allocated(history_accRatio)) deallocate(history_accRatio)
1710 : else
1711 : call append(history_accRatio, accRatio(1))
1712 : end if
1713 :
1714 : ! if in history more than 10 values, check for mean and variance
1715 13 : if (allocated(history_accRatio)) then
1716 10 : accRatio_n = size(history_accRatio, 1)
1717 10 : if (accRatio_n .ge. 10_i4) then
1718 1 : idummy = accRatio_n - 9_i4
1719 1 : accRatio_stddev = stddev(history_accRatio(idummy : accRatio_n))
1720 :
1721 : ! Check of Convergence
1722 1 : if ((accRatio_stddev .lt. Sqrt(1._dp / 12._dp * 0.05_dp**2))) then
1723 1 : iStop = .true.
1724 1 : if (printflag) then
1725 1 : write(*, *) ''
1726 1 : write(*, *) 'STOP BURN-IN with accaptence ratio of ', history_accRatio(accRatio_n)
1727 1 : write(*, *) 'final stepsize: ', stepsize
1728 1 : if (parabestChanged) then
1729 0 : write(*, *) 'better parameter set was found: ', parabest
1730 0 : write(*, *) 'with likelihood: ', likelibest
1731 : end if
1732 : end if
1733 : end if
1734 :
1735 : end if
1736 : end if
1737 :
1738 14 : trial = trial + 1_i4
1739 :
1740 : end do convergedBURNIN
1741 :
1742 : end do betterParaFound
1743 :
1744 : end if
1745 :
1746 : !----------------------------------------------------------------------
1747 : ! (2) MCMC
1748 : !----------------------------------------------------------------------
1749 :
1750 : ! just to be sure that all chains start with same initial parameter set
1751 : ! in both parallel and sequential mode
1752 : ! (although in between a new parabest will be found in chains)
1753 4 : initial_paraset_mcmc = parabest
1754 :
1755 1 : if (printflag) then
1756 1 : write(*, *) ''
1757 1 : write(*, *) '--------------------------------------------------'
1758 1 : write(*, *) 'Starting MCMC (chains = ', chains, ', iter_mcmc = ', iter_mcmc, ')'
1759 1 : write(*, *) '--------------------------------------------------'
1760 1 : write(*, *) ''
1761 : end if
1762 6 : Ipos(:) = 0_i4 ! positive accepted
1763 6 : Ineg(:) = 0_i4 ! negative accepted
1764 :
1765 : ! if all parameters converged: Sqrt(R_i) < 1.1 (see Gelman et. al: Baysian Data Analysis, p. 331ff
1766 : converged = .False.
1767 :
1768 : convergedMCMC : do while (.not. converged)
1769 :
1770 3 : if (.not. allocated(mcmc_paras_3d)) then
1771 : ! first mcmc iterations for all chains
1772 : ! probably too large, but will be resized at the end
1773 5 : allocate(mcmc_paras_3d(iter_mcmc, size(para), chains))
1774 : else
1775 : ! later mcmc iterations for all chains: iter_mcmc was increased
1776 12 : idummy = Minval(Ipos + Ineg)
1777 :
1778 : ! resize mcmc_paras_3d
1779 10 : allocate(tmp(idummy, size(para), chains))
1780 135042 : tmp(:, :, :) = mcmc_paras_3d(1 : idummy, :, :)
1781 2 : deallocate(mcmc_paras_3d)
1782 10 : allocate(mcmc_paras_3d(iter_mcmc, size(para), chains))
1783 135042 : mcmc_paras_3d(1 : idummy, :, :) = tmp(:, :, :)
1784 2 : deallocate(tmp)
1785 : end if
1786 :
1787 : !$OMP parallel default(shared) &
1788 : !$OMP private(chain, paraold, paranew, likeliold, likelinew, oddsSwitch1, oddsSwitch2, RN1, oddsRatio, ChangePara)
1789 : !$OMP do
1790 18 : parallelchain : do chain = 1, chains
1791 :
1792 15 : if (Ipos(chain) + Ineg(chain) .eq. 0_i4) then
1793 20 : paraold = initial_paraset_mcmc ! = parabest of burn-in
1794 : else
1795 40 : paraold = mcmc_paras_3d(Ipos(chain) + Ineg(chain), :, chain)
1796 : end if
1797 15 : likeliold = likelihood(paraold, eval, stddev_data)
1798 :
1799 3 : markovchainMCMC : do
1800 :
1801 : ! (A) Generate new parameter set
1802 572556 : ChangePara = .false.
1803 572556 : paranew = paraold
1804 :
1805 : call GenerateNewParameterset_dp(ParaSelectMode, paraold, truepara, rangePar, stepsize, &
1806 143139 : save_state_2(chain, :), &
1807 143139 : save_state_3(chain, :), &
1808 76006809 : paranew, ChangePara)
1809 :
1810 : ! (B) new likelihood
1811 143139 : likelinew = likelihood(paranew, eval, stddev_data)
1812 143139 : oddsSwitch1 = .false.
1813 143139 : if (loglike) then
1814 143139 : oddsRatio = likelinew - likeliold
1815 143139 : if (oddsRatio .gt. 0.0_dp) oddsSwitch1 = .true.
1816 : else
1817 0 : oddsRatio = likelinew / likeliold
1818 0 : if (oddsRatio .gt. 1.0_dp) oddsSwitch1 = .true.
1819 : end if
1820 :
1821 : ! (C) Accept or Reject?
1822 : If (oddsSwitch1) then
1823 :
1824 : ! positive accept
1825 22561 : Ipos(chain) = Ipos(chain) + 1_i4
1826 90244 : paraold = paranew
1827 22561 : likeliold = likelinew
1828 90244 : mcmc_paras_3d(Ipos(chain) + Ineg(chain), :, chain) = paranew(:)
1829 :
1830 22561 : if (printflag) then
1831 22561 : if ((Ipos(chain) + Ineg(chain) .gt. 0_i4) .and. &
1832 : (mod(Ipos(chain) + Ineg(chain), 10000_i4)) .eq. 0) then
1833 0 : write(*, *) 'Chain ', chain, ': Done ', Ipos(chain) + Ineg(chain), ' samples ...'
1834 : end if
1835 : end if
1836 :
1837 22561 : if (likelinew .gt. likelibest) then
1838 0 : parabest = paranew
1839 0 : likelibest = likelinew
1840 0 : if (printflag) then
1841 0 : write(*, *) ''
1842 0 : write(*, *) 'best para changed: ', paranew
1843 0 : write(*, *) 'likelihood new: ', likelinew
1844 0 : write(*, *) ''
1845 : end if
1846 : end if
1847 :
1848 : else
1849 :
1850 31953170 : call xor4096(seeds(chain, 1), RN1, save_state = save_state_1(chain, :))
1851 120578 : oddsSwitch2 = .false.
1852 120578 : if (loglike) then
1853 120578 : if (rn1 .lt. exp(oddsRatio)) oddsSwitch2 = .true.
1854 : else
1855 0 : if (rn1 .lt. oddsRatio) oddsSwitch2 = .true.
1856 : end if
1857 :
1858 : If (oddsSwitch2) then
1859 :
1860 : ! negative accept
1861 22439 : Ineg(chain) = Ineg(chain) + 1_i4
1862 89756 : paraold = paranew
1863 22439 : likeliold = likelinew
1864 89756 : mcmc_paras_3d(Ipos(chain) + Ineg(chain), :, chain) = paranew(:)
1865 :
1866 22439 : if (printflag) then
1867 22439 : if ((Ipos(chain) + Ineg(chain) .gt. 0_i4) .and. &
1868 : (mod(Ipos(chain) + Ineg(chain), 10000_i4)) .eq. 0) then
1869 0 : write(*, *) 'Chain ', chain, ': Done ', Ipos(chain) + Ineg(chain), ' samples ...'
1870 : end if
1871 : end if
1872 :
1873 : end if
1874 :
1875 : end if
1876 :
1877 : ! enough samples were found
1878 143139 : if ((Ipos(chain) + Ineg(chain) .gt. 0_i4) .and. &
1879 15 : (mod(Ipos(chain) + Ineg(chain), iter_mcmc) .eq. 0_i4)) exit
1880 :
1881 : end do markovchainMCMC
1882 :
1883 : end do parallelchain
1884 : !$OMP end do
1885 : !$OMP end parallel
1886 :
1887 : #ifdef FORCES_WITH_NETCDF
1888 : ! write parameter sets to temporal file
1889 3 : if (present(tmp_file)) then
1890 : ! splitting into path and filename
1891 3 : slash_pos = index(tmp_file, '/', .true.)
1892 3 : len_filename = len_trim(tmp_file)
1893 3 : path = tmp_file(1 : slash_pos)
1894 3 : filename = tmp_file(slash_pos + 1 : len_filename)
1895 : !
1896 18 : do chain = 1, chains
1897 15 : write(str, *) chain
1898 15 : write(outputfile, *) trim(adjustl(path)), trim(adjustl(str)), '_', trim(adjustl(filename))
1899 18 : if (present(iter_mcmc_in)) then
1900 0 : allocate(tmp(iter_mcmc_in, size(para, 1), 1))
1901 0 : tmp(:, :, 1) = mcmc_paras_3d(iter_mcmc - iter_mcmc_in + 1_i4 : iter_mcmc, :, chain)
1902 0 : if (iter_mcmc .ne. iter_mcmc_in) then
1903 : ! append
1904 0 : call dump_netcdf(trim(adjustl(outputfile)), tmp, append = .true.)
1905 : else
1906 : ! first time of writing
1907 0 : call dump_netcdf(trim(adjustl(outputfile)), tmp)
1908 : end if
1909 0 : deallocate(tmp)
1910 : else
1911 60 : allocate(tmp(1000_i4 * n, size(para, 1), 1))
1912 135060 : tmp(:, :, 1) = mcmc_paras_3d(iter_mcmc - (1000_i4 * n) + 1_i4 : iter_mcmc, :, chain)
1913 15 : if (iter_mcmc .ne. 1000_i4 * n) then
1914 : ! append
1915 10 : call dump_netcdf(trim(adjustl(outputfile)), tmp, append = .true.)
1916 : else
1917 : ! first time of writing
1918 5 : call dump_netcdf(trim(adjustl(outputfile)), tmp)
1919 : end if
1920 15 : deallocate(tmp)
1921 : end if
1922 : end do
1923 : end if
1924 : #else
1925 : if (present(tmp_file)) &
1926 : call error_message("MCMC: FORCES was compiled without NetCDF support but you set 'tmp_file'")
1927 : #endif
1928 :
1929 : ! test for convergence: Gelman et. al: Baysian Data Analysis, p. 331ff
1930 : ! sqrt(R) = sqrt( ( (n-1)/n * W + 1/n * B ) / W ) < 1.1 for each parameter
1931 : ! n ... last half of the sequence
1932 : ! W ... within sequence variance
1933 : ! B ... between chain variance
1934 :
1935 3 : if (printflag) then
1936 3 : write(*, *) ' '
1937 3 : write(*, *) 'Checking for convergence ....'
1938 : end if
1939 12 : sqrtR = 0.0_dp
1940 18 : n_end = minval(Ipos + Ineg)
1941 3 : n_start = n_end / 2_i4 + 1_i4
1942 :
1943 12 : do iPar = 1, size(truepara)
1944 :
1945 : ! Between chain variances
1946 54 : do chain = 1, chains
1947 0 : vDotJ(chain) = 1.0_dp / real(n_end - n_start + 1_i4, dp) * &
1948 135054 : sum(mcmc_paras_3d(n_start : n_end, truepara(iPar), chain))
1949 : end do
1950 54 : vDotDot = 1.0_dp / real(chains, dp) * sum(vDotJ)
1951 : B = real(n_end - n_start + 1_i4, dp) / real(chains - 1_i4, dp) * &
1952 54 : sum((vDotJ - vDotDot) * (vDotJ - vDotDot))
1953 :
1954 : ! Within chain variances
1955 54 : do chain = 1, chains
1956 0 : s2(chain) = 1.0_dp / real(n_end - n_start, dp) * &
1957 135054 : sum((mcmc_paras_3d(n_start : n_end, truepara(iPar), chain) - vDotJ(chain))**2)
1958 : end do
1959 54 : W = 1.0_dp / real(chains, dp) * Sum(s2)
1960 :
1961 : ! ratio sqrtR
1962 12 : if ((w .lt. tiny(1.0_dp)) .and. (b .lt. tiny(1.0_dp))) then
1963 : ! Mathematica says that this is the limit, if w and b both go to zero
1964 0 : sqrtR(truepara(iPar)) = sqrt(real(n_end - n_start, dp) / real(n_end - n_start + 1_i4, dp))
1965 : else
1966 9 : sqrtR(truepara(iPar)) = real(n_end - n_start, dp) / real(n_end - n_start + 1_i4, dp) * W + &
1967 9 : 1.0_dp / real(n_end - n_start + 1_i4, dp) * B
1968 9 : sqrtR(truepara(iPar)) = sqrt(sqrtR(truepara(iPar)) / W)
1969 : end if
1970 : end do
1971 3 : if (printflag) then
1972 12 : do i = 1, size(para)
1973 12 : if (sqrtR(i) .gt. 1.1_dp) then
1974 5 : write(*, *) ' sqrtR para #', i, ' : ', sqrtR(i), ' <-- FAILED'
1975 : else
1976 4 : write(*, *) ' sqrtR para #', i, ' : ', sqrtR(i)
1977 : end if
1978 : end do
1979 : end if
1980 6 : if (all(sqrtR .lt. 1.1_dp)) then
1981 1 : converged = .true.
1982 1 : if (printflag) then
1983 1 : write(*, *) ' --> converged (all less than 1.1)'
1984 1 : write(*, *) ' '
1985 : end if
1986 : else
1987 : ! increase number of iterations
1988 2 : if (present(iter_mcmc_in)) then
1989 0 : iter_mcmc = iter_mcmc + iter_mcmc_in
1990 : else
1991 2 : iter_mcmc = iter_mcmc + 1000_i4 * n
1992 : end if
1993 :
1994 2 : if (printflag) then
1995 2 : write(*, *) ' --> not converged (not all less than 1.1)'
1996 2 : write(*, *) ' increasing iterations to ', iter_mcmc
1997 2 : write(*, *) ' '
1998 : end if
1999 : end if
2000 :
2001 : end do convergedMCMC
2002 :
2003 : ! reshape of mcmc_paras_3d: return only 2d matrix mcmc_paras
2004 4 : allocate(mcmc_paras(size(mcmc_paras_3d, 1) * size(mcmc_paras_3d, 3), size(mcmc_paras_3d, 2)))
2005 6 : do chain = 1, chains
2006 135021 : mcmc_paras((chain - 1_i4) * size(mcmc_paras_3d, 1) + 1_i4 : chain * size(mcmc_paras_3d, 1), :) = mcmc_paras_3d(:, :, chain)
2007 : end do
2008 :
2009 1 : RETURN
2010 5 : END SUBROUTINE mcmc_stddev_dp
2011 :
2012 : !***************************************************
2013 : !* PRIVATE FUNCTIONS *
2014 : !***************************************************
2015 :
2016 : !
2017 : real(DP) function parGen_dp(old, dMax, oMin, oMax, RN, inbound)
2018 1 : use mo_kind, only : dp
2019 : implicit none
2020 : !
2021 : REAL(DP), intent(IN) :: old, dMax, oMin, oMax
2022 : REAL(DP), intent(IN) :: RN
2023 : LOGICAL, OPTIONAL, INTENT(OUT) :: inbound
2024 :
2025 : ! local
2026 : LOGICAL :: inbound2
2027 : !
2028 : inbound2 = .true.
2029 : parGen_dp = old + (rn * 2.0_dp * dMax - dMax)
2030 :
2031 : if (parGen_dp < oMin) then
2032 : parGen_dp = oMin
2033 : inbound2 = .false.
2034 : elseif(parGen_dp > oMax)then
2035 : parGen_dp = oMax
2036 : inbound2 = .false.
2037 : end if
2038 :
2039 : if (present(inbound)) inbound = inbound2
2040 :
2041 : end function parGen_dp
2042 :
2043 : ! ************************************************************************
2044 : ! normalized Parameter Change
2045 : ! ************************************************************************
2046 :
2047 268895 : real(DP) function parGenNorm_dp(old, dMax, oMin, oMax, RN, inbound)
2048 : use mo_kind, only : dp
2049 : implicit none
2050 : !
2051 : REAL(DP), intent(IN) :: old, dMax, oMin, oMax
2052 : REAL(DP), intent(IN) :: RN
2053 : LOGICAL, OPTIONAL, INTENT(OUT) :: inbound
2054 :
2055 : ! LOCAL VARIABLES
2056 268895 : REAL(DP) :: RN_scale ! scaled RN
2057 268895 : REAL(DP) :: old_scaled ! for normalization
2058 : LOGICAL :: inbound2
2059 :
2060 268895 : inbound2 = .true.
2061 :
2062 : !(1) normalize
2063 268895 : old_scaled = (old - oMin) / (oMax - oMin)
2064 :
2065 : !(2) scale RN = [-dMax, dMax]
2066 : !RN_scale = (RN * 2.0_dp * dMax) - dMax
2067 : !(2) scale RN distrib. N(0, dMax)
2068 268895 : RN_scale = (RN * dMax)
2069 :
2070 : !(3) generate new parameter value + check inbound
2071 268895 : parGenNorm_dp = old_scaled + RN_scale
2072 : ! Parameter is bounded between Max and Min.
2073 : ! if (parGenNorm_dp < 0.0_dp) then
2074 : ! parGenNorm_dp = 0.0_dp
2075 : ! inbound2 = .false.
2076 : ! elseif(parGenNorm_dp > 1.0_dp)then
2077 : ! parGenNorm_dp = 1.0_dp
2078 : ! inbound2 = .false.
2079 : ! end if
2080 :
2081 268895 : if (present(inbound)) inbound = inbound2
2082 :
2083 : !(4) convert back to original range
2084 268895 : parGenNorm_dp = parGenNorm_dp * (oMax - oMin) + oMin
2085 :
2086 268895 : end function parGenNorm_dp
2087 :
2088 537790 : recursive subroutine GenerateNewParameterset_dp(ParaSelectMode, paraold, truepara, rangePar, stepsize, &
2089 268895 : save_state_2, save_state_3, paranew, ChangePara)
2090 :
2091 : integer(i4), intent(in) :: ParaSelectMode
2092 : real(dp), dimension(:), intent(in) :: paraold
2093 : integer(i4), dimension(:), intent(in) :: truepara
2094 : real(dp), dimension(:, :), intent(in) :: rangePar
2095 : real(dp), dimension(:), intent(in) :: stepsize
2096 : integer(I8), dimension(n_save_state), intent(inout) :: save_state_2 ! optional argument for restarting RN stream 2
2097 : integer(I8), dimension(n_save_state), intent(inout) :: save_state_3 ! optional argument for restarting RN stream 3
2098 : real(dp), dimension(size(paraold)), intent(out) :: paranew
2099 : logical, dimension(size(paraold)), intent(out) :: ChangePara
2100 :
2101 : ! local variables
2102 268895 : REAL(DP) :: RN2, RN3
2103 : logical :: inbound
2104 : integer(i4) :: i, iPar
2105 :
2106 1344475 : paranew = paraold
2107 1075580 : ChangePara = .False.
2108 :
2109 268895 : select case(ParaSelectMode)
2110 : case(1_i4) ! change half of the parameters
2111 0 : do i = 1, size(truepara)
2112 0 : call xor4096(0_i8, RN2, save_state = save_state_2)
2113 0 : if (RN2 > 0.5_dp) then
2114 0 : iPar = truepara(i)
2115 0 : inbound = .false.
2116 0 : call xor4096(0_i8, RN3, save_state = save_state_3)
2117 0 : paranew(iPar) = parGenNorm_dp(paraold(iPar), stepsize(iPar), rangePar(iPar, 1), rangePar(iPar, 2), RN3, inbound)
2118 0 : ChangePara(iPar) = .True.
2119 : end if
2120 : end do
2121 : ! if no parameter was selected, recall and change one
2122 0 : if (count(ChangePara) .eq. 0_i4) then
2123 0 : call xor4096(0_i8, RN2, save_state = save_state_2)
2124 0 : iPar = truepara(int((RN2 * real(size(truepara), dp)) + 1.0_dp, i4))
2125 0 : inbound = .false.
2126 0 : call xor4096(0_i8, RN3, save_state = save_state_3)
2127 0 : paranew(iPar) = parGenNorm_dp(paraold(iPar), stepsize(iPar), rangePar(iPar, 1), rangePar(iPar, 2), RN3, inbound)
2128 0 : ChangePara(iPar) = .True.
2129 : end if
2130 :
2131 : case(2_i4) ! change only one
2132 268895 : call xor4096(0_i8, RN2, save_state = save_state_2)
2133 268895 : iPar = truepara(int((RN2 * real(size(truepara), dp)) + 1.0_dp, i4))
2134 268895 : inbound = .false.
2135 268895 : call xor4096g(0_i8, RN3, save_state = save_state_3)
2136 268895 : paranew(iPar) = parGenNorm_dp(paraold(iPar), stepsize(iPar), rangePar(iPar, 1), rangePar(iPar, 2), RN3, inbound)
2137 : ! write(*,*) 'p_old(',iPar,') = ',paraold(iPar)
2138 : ! write(*,*) 'p_new(',iPar,') = ',paranew(iPar)
2139 268895 : ChangePara(iPar) = .True.
2140 :
2141 : case(3_i4) ! change all
2142 268895 : do i = 1, size(truepara)
2143 0 : iPar = truepara(i)
2144 0 : inbound = .false.
2145 0 : do while(.not. inbound)
2146 0 : call xor4096g(0_i8, RN2, save_state = save_state_3)
2147 0 : paranew(iPar) = parGenNorm_dp(paraold(iPar), stepsize(iPar), rangePar(iPar, 1), rangePar(iPar, 2), RN3, inbound)
2148 : end do
2149 0 : ChangePara(iPar) = .True.
2150 : end do
2151 :
2152 : end select
2153 :
2154 268895 : end subroutine GenerateNewParameterset_dp
2155 :
2156 : END MODULE mo_mcmc
|
