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
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 : ! - directory where to write ascii output
75 : character(len = *), intent(in) :: dirConfigOut
76 :
77 : ! - best objective function value obtained during optimization
78 : real(dp), intent(out) :: funcbest
79 :
80 : ! true = parameter will be optimized = parameter(i,4) = 1
81 : ! false = parameter will not be optimized = parameter(i,4) = 0
82 : logical, intent(out), allocatable, dimension(:) :: maskpara
83 :
84 : integer(i4) :: ii
85 :
86 : ! current timer number
87 : integer(i4) :: iTimer
88 :
89 : ! parameter sets sampled during burnin
90 5 : real(dp), allocatable, dimension(:, :) :: burnin_paras
91 :
92 : ! parameter sets sampled during proper mcmc
93 5 : real(dp), allocatable, dimension(:, :) :: mcmc_paras
94 :
95 : ! pre-determined stepsize
96 : ! real(dp), dimension(:), allocatable :: step
97 :
98 : integer(i4) :: npara
99 :
100 : ! global_parameters but includes a and b for likelihood
101 5 : real(dp), allocatable, dimension(:, :) :: local_parameters
102 :
103 : ! maskpara but includes a and b for likelihood
104 5 : logical, allocatable, dimension(:) :: local_maskpara
105 :
106 : ! local seed used for optimization
107 : integer(i8) :: iseed
108 :
109 : ! file for temporal optimization outputs
110 : character(256) :: tFile
111 :
112 : ! file for temporal SCE optimization outputs
113 : character(256) :: pFile
114 :
115 :
116 : ! -------------------------------------------------------------------------
117 : ! START
118 : ! -------------------------------------------------------------------------
119 5 : call message(' Start optimization')
120 5 : iTimer = 1
121 5 : call timer_start(iTimer)
122 :
123 : ! mask parameter which have a FLAG=0 in mhm_parameter.nml
124 : ! maskpara = true : parameter will be optimized
125 : ! maskpara = false : parameter is discarded during optimization
126 5 : npara = size(global_parameters, 1)
127 15 : allocate(maskpara(npara))
128 273 : maskpara = .true.
129 273 : do ii = 1, npara
130 273 : if (nint(global_parameters(ii, 4), i4) .eq. 0_i4) then
131 29 : maskpara(ii) = .false.
132 : end if
133 : end do
134 :
135 : ! add two extra parameter for optimisation of likelihood
136 5 : if (opti_function == 8) then
137 0 : allocate(local_parameters(npara + 2, size(global_parameters, 2)))
138 0 : allocate(local_maskpara(npara + 2))
139 : ! setting step sizes manually
140 : ! allocate(step(npara+2))
141 0 : local_parameters(1 : npara, :) = global_parameters(:, :)
142 0 : local_maskpara(1 : npara) = maskpara(:)
143 0 : local_parameters(npara + 1, 1) = 0.001_dp
144 0 : local_parameters(npara + 1, 2) = 100._dp
145 0 : local_parameters(npara + 1, 3) = 1._dp
146 0 : local_parameters(npara + 1, 4) = 1._dp
147 0 : local_parameters(npara + 1, 5) = 0._dp
148 0 : local_parameters(npara + 2, 1) = 0.001_dp
149 0 : local_parameters(npara + 2, 2) = 10._dp
150 0 : local_parameters(npara + 2, 3) = 0.1_dp
151 0 : local_parameters(npara + 2, 4) = 1._dp
152 0 : local_parameters(npara + 2, 5) = 0._dp
153 0 : local_maskpara(npara + 1 :) = .true.
154 0 : if ((opti_method == 0) .and. (.not. mcmc_opti)) then ! MCMC but only for parameter uncertainties
155 0 : local_parameters(npara + 1, 3) = mcmc_error_params(1)
156 0 : local_parameters(npara + 2, 3) = mcmc_error_params(2)
157 0 : local_maskpara(npara + 1 :) = .false.
158 : end if
159 : else
160 20 : allocate(local_parameters(npara, size(global_parameters, 2)))
161 10 : allocate(local_maskpara(npara))
162 1375 : local_parameters = global_parameters
163 278 : local_maskpara = maskpara
164 : end if
165 :
166 : ! Seed for random numbers in optimisation
167 5 : if (seed .gt. 0_i8) then ! fixed user-defined seed
168 5 : iseed = seed
169 : else ! flexible clock-time seed
170 0 : call get_timeseed(iseed)
171 : end if
172 :
173 0 : select case (opti_method)
174 : case (0)
175 0 : call message(' Use MCMC')
176 :
177 0 : tFile = trim(adjustl(dirConfigOut)) // 'mcmc_tmp_parasets.nc'
178 :
179 : ! setting step sizes manually
180 : ! step=(/ 1.00000000000000, ... , /)
181 :
182 5 : select case (opti_function)
183 : case (8)
184 0 : call message(' Use MCMC')
185 0 : call mcmc(eval, objective, local_parameters(:, 3), local_parameters(:, 1 : 2), mcmc_paras, burnin_paras, &
186 : ParaSelectMode_in = 2_i4, tmp_file = tFile, &
187 : maskpara_in = local_maskpara, &
188 : restart = optimize_restart, restart_file = 'mo_mcmc.restart', &
189 : ! stepsize_in=step, &
190 0 : seed_in = iseed, loglike_in = .true., printflag_in = .true.)
191 : case (4)
192 0 : if (optimize_restart) then
193 0 : call error_message('ERROR: A restart of this optimization method is not implemented yet!')
194 : end if
195 0 : call message(' Use MCMC_STDDEV')
196 0 : call mcmc_stddev(eval, objective, local_parameters(:, 3), local_parameters(:, 1 : 2), mcmc_paras, burnin_paras, &
197 : ParaSelectMode_in = 2_i4, tmp_file = tFile, &
198 : maskpara_in = local_maskpara, &
199 0 : seed_in = iseed, loglike_in = .true., printflag_in = .true.)
200 : case default
201 0 : call error_message("Error objective: This opti_function is either not implemented yet.")
202 : end select
203 :
204 : case (1)
205 5 : call message(' Use DDS')
206 :
207 5 : tFile = trim(adjustl(dirConfigOut)) // 'dds_results.out'
208 :
209 5 : if (optimize_restart) then
210 0 : call error_message('ERROR: A restart of this optimization method is not implemented yet!')
211 : end if
212 : ! use fixed user-defined seed
213 : #ifdef MPI
214 : local_parameters(:, 3) = dds(eval, objective, local_parameters(:, 3), local_parameters(:, 1 : 2), &
215 : maxiter = int(nIterations, i8), r = dds_r, seed = iseed, &
216 : tmp_file = tFile, comm = domainMeta%comMaster, mask = local_maskpara, &
217 : funcbest = funcbest)
218 : #else
219 10 : local_parameters(:, 3) = dds(eval, objective, local_parameters(:, 3), local_parameters(:, 1 : 2), &
220 : maxiter = int(nIterations, i8), r = dds_r, seed = iseed, &
221 : tmp_file = tFile, mask = local_maskpara, &
222 283 : funcbest = funcbest)
223 : #endif
224 :
225 : case (2)
226 0 : call message(' Use Simulated Annealing')
227 :
228 0 : tFile = trim(adjustl(dirConfigOut)) // 'anneal_results.out'
229 :
230 0 : if (optimize_restart) then
231 0 : call error_message('ERROR: A restart of this optimization method is not implemented yet!')
232 : end if
233 :
234 0 : if (sa_temp .gt. 0.0_dp) then
235 : ! use fixed user-defined seed and user-defined initial temperature
236 0 : local_parameters(:, 3) = anneal(eval, objective, local_parameters(:, 3), local_parameters(:, 1 : 2), &
237 : temp = sa_temp, seeds = (/iseed, iseed + 1000_i8, iseed + 2000_i8/), nITERmax = nIterations, &
238 : tmp_file = tFile, maskpara = local_maskpara, &
239 0 : funcbest = funcbest)
240 : else
241 : ! use fixed user-defined seed and adaptive initial temperature
242 0 : local_parameters(:, 3) = anneal(eval, objective, local_parameters(:, 3), local_parameters(:, 1 : 2), &
243 : seeds = (/iseed, iseed + 1000_i8, iseed + 2000_i8/), nITERmax = nIterations, &
244 : tmp_file = tFile, maskpara = local_maskpara, &
245 0 : funcbest = funcbest)
246 : end if
247 : case (3)
248 0 : call message(' Use SCE')
249 :
250 0 : tFile = trim(adjustl(dirConfigOut)) // 'sce_results.out'
251 0 : pFile = trim(adjustl(dirConfigOut)) // 'sce_population.out'
252 :
253 : ! use fixed user-defined seed
254 0 : local_parameters(:, 3) = sce(eval, objective, local_parameters(:, 3), local_parameters(:, 1 : 2), &
255 : mymaxn = int(nIterations, i8), myseed = iseed, myngs = sce_ngs, mynpg = sce_npg, mynps = sce_nps, &
256 : parallel = .false., mymask = local_maskpara, &
257 : restart = optimize_restart, restart_file = 'mo_sce.restart', &
258 : #ifdef MPI
259 : comm = domainMeta%comMaster, &
260 : #endif
261 : tmp_file = tFile, popul_file = pFile, &
262 0 : bestf = funcbest)
263 : case default
264 5 : call error_message('mRM', 'This optimization method is not implemented.')
265 : end select
266 5 : call timer_stop(iTimer)
267 5 : call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
268 :
269 1370 : global_parameters(:, :) = local_parameters(1 : npara, :)
270 273 : maskpara(:) = local_maskpara(1 : npara)
271 :
272 5 : deallocate(local_parameters)
273 5 : deallocate(local_maskpara)
274 :
275 10 : end subroutine optimization
276 :
277 : end module mo_optimization
|
