Line data Source code
1 : !> \file mo_optimization.f90
2 : !> \brief \copybrief mo_optimization
3 : !> \details \copydetails mo_optimization
4 :
5 : !> \brief Wrapper subroutine for optimization against runoff and sm.
6 : !> \details This module provides a wrapper subroutine for optimization of mRM/mHM against runoff or soil moisture.
7 : !> \authors Stephan Thober
8 : !> \date Oct 2015
9 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
10 : !! mHM is released under the LGPLv3+ license \license_note
11 : !> \ingroup f_common
12 : module mo_optimization
13 : use mo_kind, only : i4, i8, dp
14 : use mo_optimization_utils, only : eval_interface, objective_interface, mhm_optimizee
15 : implicit none
16 : private
17 : public :: optimization
18 :
19 : contains
20 : ! ------------------------------------------------------------------
21 :
22 : ! NAME
23 : ! optimization
24 :
25 : ! PURPOSE
26 : !> \brief Wrapper for optimization.
27 :
28 : !> \details This subroutine selects the optimization defined in a namelist,
29 : !> i.e. the global variable \e opti\_method.
30 : !> It return the objective function value for a specific parameter set.
31 :
32 : ! INTENT(IN)
33 : !> \param[in] "procedure(eval_interface) :: eval"
34 : !> \param[in] "procedure(objective_interface) :: objective" - objective function used in the optimization
35 : !> \param[in] "character(len = *) :: dirConfigOut" - directory where to write ascii output
36 :
37 : ! INTENT(OUT)
38 : !> \param[out] "real(dp) :: funcbest" - best objective function value obtained during optimization
39 : !> \param[out] "logical, dimension(:) :: maskpara" true = parameter will be optimized = parameter(i,4) = 1
40 : !> false = parameter will not be optimized = parameter(i,4) = 0
41 :
42 : ! HISTORY
43 : !> \authors Matthias Cuntz, Luis Samaniego, Juliane Mai, Matthias Zink and Stephan Thober
44 :
45 : !> \date Oct 2015
46 :
47 : ! Modifications:
48 :
49 5 : subroutine optimization(eval, objective, dirConfigOut, funcBest, maskpara)
50 :
51 : use mo_anneal, only : anneal
52 : use mo_common_mHM_mRM_variables, only : dds_r, mcmc_error_params, mcmc_opti, nIterations, opti_function, opti_method, &
53 : optimize_restart, sa_temp, sce_ngs, sce_npg, sce_nps, seed
54 : use mo_common_variables, only : global_parameters
55 : #ifdef MPI
56 : use mo_common_variables, only : domainMeta
57 : #endif
58 : use mo_dds, only : dds
59 : use mo_mcmc, only : mcmc, mcmc_stddev
60 : use mo_message, only : message, error_message
61 : use mo_sce, only : sce
62 : use mo_string_utils, only : num2str
63 : use mo_timer, only : timer_get, timer_start, &
64 : timer_stop
65 : use mo_xor4096, only : get_timeseed
66 :
67 : implicit none
68 :
69 : procedure(eval_interface), intent(in), pointer :: eval
70 :
71 : ! - objective function used in the optimization
72 : procedure(objective_interface), intent(in), pointer :: objective
73 :
74 : ! - objective target, using objective and eval
75 : type(mhm_optimizee) :: optimize_target
76 :
77 : ! - directory where to write ascii output
78 : character(len = *), intent(in) :: dirConfigOut
79 :
80 : ! - best objective function value obtained during optimization
81 : real(dp), intent(out) :: funcbest
82 :
83 : ! true = parameter will be optimized = parameter(i,4) = 1
84 : ! false = parameter will not be optimized = parameter(i,4) = 0
85 : logical, intent(out), allocatable, dimension(:) :: maskpara
86 :
87 : integer(i4) :: ii
88 :
89 : ! current timer number
90 : integer(i4) :: iTimer
91 :
92 : ! parameter sets sampled during burnin
93 5 : real(dp), allocatable, dimension(:, :) :: burnin_paras
94 :
95 : ! parameter sets sampled during proper mcmc
96 5 : real(dp), allocatable, dimension(:, :) :: mcmc_paras
97 :
98 : ! pre-determined stepsize
99 : ! real(dp), dimension(:), allocatable :: step
100 :
101 : integer(i4) :: npara
102 :
103 : ! global_parameters but includes a and b for likelihood
104 5 : real(dp), allocatable, dimension(:, :) :: local_parameters
105 :
106 : ! maskpara but includes a and b for likelihood
107 5 : logical, allocatable, dimension(:) :: local_maskpara
108 :
109 : ! local seed used for optimization
110 : integer(i8) :: iseed
111 :
112 : ! file for temporal optimization outputs
113 : character(256) :: tFile
114 :
115 : ! file for temporal SCE optimization outputs
116 : character(256) :: pFile
117 :
118 :
119 : ! -------------------------------------------------------------------------
120 : ! START
121 : ! -------------------------------------------------------------------------
122 5 : call message(' Start optimization')
123 5 : iTimer = 1
124 5 : call timer_start(iTimer)
125 :
126 : ! mask parameter which have a FLAG=0 in mhm_parameter.nml
127 : ! maskpara = true : parameter will be optimized
128 : ! maskpara = false : parameter is discarded during optimization
129 5 : npara = size(global_parameters, 1)
130 15 : allocate(maskpara(npara))
131 273 : maskpara = .true.
132 273 : do ii = 1, npara
133 273 : if (nint(global_parameters(ii, 4), i4) .eq. 0_i4) then
134 29 : maskpara(ii) = .false.
135 : end if
136 : end do
137 :
138 5 : optimize_target%eval_pointer => eval
139 5 : optimize_target%obj_pointer => objective
140 :
141 : ! add two extra parameter for optimisation of likelihood
142 5 : if (opti_function == 8) then
143 0 : allocate(local_parameters(npara + 2, size(global_parameters, 2)))
144 0 : allocate(local_maskpara(npara + 2))
145 : ! setting step sizes manually
146 : ! allocate(step(npara+2))
147 0 : local_parameters(1 : npara, :) = global_parameters(:, :)
148 0 : local_maskpara(1 : npara) = maskpara(:)
149 0 : local_parameters(npara + 1, 1) = 0.001_dp
150 0 : local_parameters(npara + 1, 2) = 100._dp
151 0 : local_parameters(npara + 1, 3) = 1._dp
152 0 : local_parameters(npara + 1, 4) = 1._dp
153 0 : local_parameters(npara + 1, 5) = 0._dp
154 0 : local_parameters(npara + 2, 1) = 0.001_dp
155 0 : local_parameters(npara + 2, 2) = 10._dp
156 0 : local_parameters(npara + 2, 3) = 0.1_dp
157 0 : local_parameters(npara + 2, 4) = 1._dp
158 0 : local_parameters(npara + 2, 5) = 0._dp
159 0 : local_maskpara(npara + 1 :) = .true.
160 0 : if ((opti_method == 0) .and. (.not. mcmc_opti)) then ! MCMC but only for parameter uncertainties
161 0 : local_parameters(npara + 1, 3) = mcmc_error_params(1)
162 0 : local_parameters(npara + 2, 3) = mcmc_error_params(2)
163 0 : local_maskpara(npara + 1 :) = .false.
164 : end if
165 : else
166 20 : allocate(local_parameters(npara, size(global_parameters, 2)))
167 10 : allocate(local_maskpara(npara))
168 1375 : local_parameters = global_parameters
169 278 : local_maskpara = maskpara
170 : end if
171 :
172 : ! Seed for random numbers in optimisation
173 5 : if (seed .gt. 0_i8) then ! fixed user-defined seed
174 5 : iseed = seed
175 : else ! flexible clock-time seed
176 0 : call get_timeseed(iseed)
177 : end if
178 :
179 0 : select case (opti_method)
180 : case (0)
181 0 : call message(' Use MCMC')
182 :
183 0 : tFile = trim(adjustl(dirConfigOut)) // 'mcmc_tmp_parasets.nc'
184 :
185 : ! setting step sizes manually
186 : ! step=(/ 1.00000000000000, ... , /)
187 :
188 5 : select case (opti_function)
189 : case (8)
190 0 : call message(' Use MCMC')
191 0 : call mcmc(optimize_target, local_parameters(:, 3), local_parameters(:, 1 : 2), mcmc_paras, burnin_paras, &
192 : ParaSelectMode_in = 2_i4, tmp_file = tFile, &
193 : maskpara_in = local_maskpara, &
194 : restart = optimize_restart, restart_file = 'mo_mcmc.restart', &
195 : ! stepsize_in=step, &
196 0 : seed_in = iseed, loglike_in = .true., printflag_in = .true.)
197 : case (4)
198 0 : if (optimize_restart) then
199 0 : call error_message('ERROR: A restart of this optimization method is not implemented yet!')
200 : end if
201 0 : call message(' Use MCMC_STDDEV')
202 0 : call mcmc_stddev(optimize_target, local_parameters(:, 3), local_parameters(:, 1 : 2), mcmc_paras, burnin_paras, &
203 : ParaSelectMode_in = 2_i4, tmp_file = tFile, &
204 : maskpara_in = local_maskpara, &
205 0 : seed_in = iseed, loglike_in = .true., printflag_in = .true.)
206 : case default
207 0 : call error_message("Error objective: This opti_function is either not implemented yet.")
208 : end select
209 :
210 : case (1)
211 5 : call message(' Use DDS')
212 :
213 5 : tFile = trim(adjustl(dirConfigOut)) // 'dds_results.out'
214 :
215 5 : if (optimize_restart) then
216 0 : call error_message('ERROR: A restart of this optimization method is not implemented yet!')
217 : end if
218 : ! use fixed user-defined seed
219 : #ifdef MPI
220 : local_parameters(:, 3) = dds(optimize_target, local_parameters(:, 3), local_parameters(:, 1 : 2), &
221 : maxiter = int(nIterations, i8), r = dds_r, seed = iseed, &
222 : tmp_file = tFile, comm = domainMeta%comMaster, mask = local_maskpara, &
223 : funcbest = funcbest)
224 : #else
225 10 : local_parameters(:, 3) = dds(optimize_target, local_parameters(:, 3), local_parameters(:, 1 : 2), &
226 : maxiter = int(nIterations, i8), r = dds_r, seed = iseed, &
227 : tmp_file = tFile, mask = local_maskpara, &
228 283 : funcbest = funcbest)
229 : #endif
230 :
231 : case (2)
232 0 : call message(' Use Simulated Annealing')
233 :
234 0 : tFile = trim(adjustl(dirConfigOut)) // 'anneal_results.out'
235 :
236 0 : if (optimize_restart) then
237 0 : call error_message('ERROR: A restart of this optimization method is not implemented yet!')
238 : end if
239 :
240 0 : if (sa_temp .gt. 0.0_dp) then
241 : ! use fixed user-defined seed and user-defined initial temperature
242 0 : local_parameters(:, 3) = anneal(optimize_target, local_parameters(:, 3), local_parameters(:, 1 : 2), &
243 : temp = sa_temp, seeds = (/iseed, iseed + 1000_i8, iseed + 2000_i8/), nITERmax = nIterations, &
244 : tmp_file = tFile, maskpara = local_maskpara, &
245 0 : funcbest = funcbest)
246 : else
247 : ! use fixed user-defined seed and adaptive initial temperature
248 0 : local_parameters(:, 3) = anneal(optimize_target, local_parameters(:, 3), local_parameters(:, 1 : 2), &
249 : seeds = (/iseed, iseed + 1000_i8, iseed + 2000_i8/), nITERmax = nIterations, &
250 : tmp_file = tFile, maskpara = local_maskpara, &
251 0 : funcbest = funcbest)
252 : end if
253 : case (3)
254 0 : call message(' Use SCE')
255 :
256 0 : tFile = trim(adjustl(dirConfigOut)) // 'sce_results.out'
257 0 : pFile = trim(adjustl(dirConfigOut)) // 'sce_population.out'
258 :
259 : ! use fixed user-defined seed
260 0 : local_parameters(:, 3) = sce(optimize_target, local_parameters(:, 3), local_parameters(:, 1 : 2), &
261 : mymaxn = int(nIterations, i8), myseed = iseed, myngs = sce_ngs, mynpg = sce_npg, mynps = sce_nps, &
262 : parallel = .false., mymask = local_maskpara, &
263 : restart = optimize_restart, restart_file = 'mo_sce.restart', &
264 : #ifdef MPI
265 : comm = domainMeta%comMaster, &
266 : #endif
267 : tmp_file = tFile, popul_file = pFile, &
268 0 : bestf = funcbest)
269 : case default
270 5 : call error_message('mRM', 'This optimization method is not implemented.')
271 : end select
272 5 : call timer_stop(iTimer)
273 5 : call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
274 :
275 1370 : global_parameters(:, :) = local_parameters(1 : npara, :)
276 273 : maskpara(:) = local_maskpara(1 : npara)
277 :
278 5 : deallocate(local_parameters)
279 5 : deallocate(local_maskpara)
280 :
281 10 : end subroutine optimization
282 :
283 : end module mo_optimization
|