Line data Source code
1 : !> \file mo_dds.f90
2 : !> \brief \copybrief mo_dds
3 : !> \details \copydetails mo_dds
4 :
5 : !> \brief Dynamically Dimensioned Search (DDS)
6 : !> \details This module provides routines for Dynamically Dimensioned Search (DDS)
7 : !! of Tolson and Shoemaker (2007). It searches the minimum or maximum of a user-specified function,
8 : !! using an n-dimensional continuous global optimization algorithm (DDS).
9 : !> \authors Bryan Tolson, modified by Rohini Kumar, Matthias Cuntz and Juliane Mai.
10 : !> \date Jul 2012
11 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
12 : !! FORCES is released under the LGPLv3+ license \license_note
13 : module mo_dds
14 :
15 : IMPLICIT NONE
16 :
17 : PRIVATE
18 :
19 : PUBLIC :: DDS ! Dynamically Dimensioned Search (DDS)
20 : PUBLIC :: MDDS ! Modified Dynamically Dimensioned Search (DDS)
21 :
22 : ! ------------------------------------------------------------------
23 :
24 : CONTAINS
25 :
26 : ! ------------------------------------------------------------------
27 :
28 : !> \brief DDS
29 :
30 : !> \details Searches Minimum or Maximum of a user-specified function using
31 : !! Dynamically Dimensioned Search (DDS).\n
32 : !! DDS is an n-dimensional continuous global optimization algorithm.
33 : !! It is coded as a minimizer but one can give maxit=True in a maximization problem,
34 : !! so that the algorithm minimizes the negative of the objective function F=(-1*F).
35 : !!
36 : !! The function to be minimized is the first argument of DDS and must be defined as \n
37 : !! \code{.f90}
38 : !! function func(p)
39 : !! use mo_kind, only: dp
40 : !! implicit none
41 : !! real(dp), dimension(:), intent(in) :: p
42 : !! real(dp) :: func
43 : !! end function func
44 : !! \endcode
45 : !!
46 : !! \b Example
47 : !!
48 : !! \code{.f90}
49 : !! dv_range(:,1) = (/ -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0 /)
50 : !! dv_range(:,2) = (/ 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0 /)
51 : !! dv_ini = (/ -.226265E+01, -.130187E+01, -.151219E+01, 0.133983E+00, 0.988159E+00, &
52 : !! -.495074E+01, -.126574E+02, 0.572684E+00, 0.303864E+01, 0.343031E+01 /)
53 : !! dv_opt = DDS(griewank, dv_ini, dv_range)
54 : !! \endcode
55 : !!
56 : !! See also example in test directory.
57 : !!
58 : !! \b Literature
59 : !! 1. Tolson, B. A., and C. A. Shoemaker (2007)
60 : !! _Dynamically dimensioned search algorithm for computationally efficient watershed
61 : !! model calibration_, Water Resour. Res., 43, W01413, doi:10.1029/2005WR004723.
62 : !!
63 : !> \param[in] "real(dp) :: obj_func(p)" Function on which to search the minimum
64 : !> \param[in] "real(dp) :: pini(:)" inital value of decision variables
65 : !> \param[in] "real(dp) :: prange(size(pini),2)" Min/max range of decision variables
66 : !> \param[in] "real(dp), optional :: r" DDS perturbation parameter\n
67 : !! (default: 0.2)
68 : !> \param[in] "integer(i8), optional :: seed" User seed to initialise the random number generator
69 : !! (default: None)
70 : !> \param[in] "integer(i8), optional :: maxiter" Maximum number of iteration or function evaluation
71 : !! (default: 1000)
72 : !> \param[in] "logical, optional :: maxit" Maximization (.True.) or
73 : !! minimization (.False.) of function
74 : !! (default: .False.)
75 : !> \param[in] "logical, optional :: mask(size(pini))" parameter to be optimized (true or false)
76 : !! (default: .True.)
77 : !> \param[in] "character(len=*) , optional :: tmp_file" file with temporal output
78 : !> \param[out] "real(dp), optional :: funcbest" the best value of the function.
79 : !> \param[out] "real(dp), optional, allocatable :: history(:)" the history of best function values,
80 : !! history(maxiter)=funcbest\n
81 : !! allocatable only to be in correspondance with other optimization routines
82 : !> \retval "real(dp) :: DDS" The parameters of the point which is estimated to minimize the function.
83 :
84 : !> \author Bryan Tolson
85 : !> \date Feb 2007
86 :
87 : !> \author Rohini Kumar
88 : !> \date Feb 2008
89 :
90 : !> \author Matthias Cuntz & Juliane Mai
91 : !> \date Jul 2012
92 : !! - module
93 :
94 : !> \author Juliane Mai
95 : !> \date Aug 2012
96 : !! - optional argument funcbest added
97 : !> \date Nov 2012
98 : !! - masked parameter
99 : !> \date Dec 2012
100 : !! - history output
101 :
102 : !> \author Pallav Kumar Shrestha
103 : !> \date Jun 2019
104 : !! - Added "dds_results.out.current" output file with current iteration (non-updated) parameters
105 :
106 : #ifdef MPI
107 : function DDS(eval, obj_func, pini, prange, r, seed, maxiter, maxit, mask, tmp_file, comm, funcbest, history)
108 : #else
109 4 : function DDS(eval, obj_func, pini, prange, r, seed, maxiter, maxit, mask, tmp_file, funcbest, history)
110 : #endif
111 :
112 : use mo_kind, only : i4, i8, dp
113 : use mo_xor4096, only : xor4096, xor4096g
114 : use mo_optimization_utils, only : eval_interface, objective_interface
115 : #ifdef MPI
116 : use mpi_f08
117 : #endif
118 :
119 : implicit none
120 :
121 : procedure(eval_interface), INTENT(IN), POINTER :: eval
122 : procedure(objective_interface), intent(in), pointer :: obj_func
123 :
124 : real(dp), dimension(:), intent(in) :: pini ! inital value of decision variables
125 : real(dp), dimension(:, :), intent(in) :: prange ! Min/max values of decision variables
126 : real(dp), optional, intent(in) :: r ! DDS perturbation parameter (-> 0.2 by default)
127 : integer(i8), optional, intent(in) :: seed ! User seed to initialise the random number generator
128 : integer(i8), optional, intent(in) :: maxiter ! Maximum number of iteration or function evaluation
129 : logical, optional, intent(in) :: maxit ! Maximization or minimization of function
130 : logical, dimension(:), optional, intent(in) :: mask ! parameter to be optimized (true or false)
131 : character(len = *), optional, intent(in) :: tmp_file ! file for temporal output
132 : #ifdef MPI
133 : type(MPI_Comm), optional, intent(in) :: comm ! MPI communicator
134 : #endif
135 : real(dp), optional, intent(out) :: funcbest ! Best value of the function.
136 : real(dp), dimension(:), optional, intent(out), &
137 : allocatable :: history ! History of objective function values
138 : real(dp), dimension(size(pini)) :: DDS ! Best value of decision variables
139 :
140 : ! Local variables
141 : integer(i4) :: pnum ! Total number of decision variables
142 : integer(i8) :: iseed ! User given seed
143 : integer(i8) :: imaxiter ! Maximum number of iteration or function evaluation
144 1 : real(dp) :: ir ! DDS perturbation parameter
145 1 : real(dp) :: imaxit ! Maximization or minimization of function
146 1 : real(dp) :: of_new, of_best ! intermediate results
147 1 : real(dp) :: Pn, new_value ! intermediate results
148 11 : real(dp), dimension(size(pini)) :: pnew ! Test value of decision variables
149 1 : real(dp) :: ranval ! random value
150 : integer(i8) :: i ! maxiter=i8
151 : integer(i4) :: j, dvn_count, dv ! pnum=i4
152 : integer(i4) :: idummy ! dummy vaiable
153 : integer, dimension(8) :: sdate ! date_and_time return
154 1 : logical, dimension(size(pini)) :: maske ! parameter to be optimized (true or false)
155 1 : integer(i4), dimension(:), allocatable :: truepara ! parameter to be optimized (their indexes)
156 : #ifdef MPI
157 : integer(i4) :: ierror
158 : logical :: do_obj_loop
159 : integer(i4) :: iproc, nproc
160 : #endif
161 :
162 : ! Check input
163 1 : pnum = size(pini)
164 1 : if (size(prange, 1) /= pnum) stop 'Error DDS: size(prange,1) /= size(pini)'
165 1 : if (size(prange, 2) /= 2) stop 'Error DDS: size(prange,2) /= 2'
166 : ! r Perturbation parameter
167 1 : ir = 0.2_dp
168 1 : if (present(r)) ir = r
169 1 : if (ir <= 0.0_dp .or. ir > 1.0_dp) stop 'Error DDS: DDS perturbation parameter (0.0, 1.0]'
170 : ! max. iteration
171 1 : imaxiter = 1000
172 1 : if (present(maxiter)) imaxiter = maxiter
173 1 : if (imaxiter < 6) stop 'Error DDS: max function evals must be minimum 6'
174 : ! history output
175 1 : if (present(history)) then
176 0 : allocate(history(imaxiter))
177 : end if
178 : ! Min or max objective function
179 1 : imaxit = 1.0_dp
180 1 : if (present(maxit)) then
181 1 : if (maxit) imaxit = -1.0_dp
182 : end if
183 : ! Given seed
184 1 : iseed = 0
185 1 : if (present(seed)) iseed = seed
186 1 : iseed = max(iseed, 0_i8)
187 :
188 1 : if (present(mask)) then
189 0 : if (count(mask) .eq. 0_i4) then
190 0 : stop 'Input argument mask: At least one element has to be true'
191 : else
192 0 : maske = mask
193 : end if
194 : else
195 11 : maske = .true.
196 : end if
197 :
198 13 : allocate (truepara(count(maske)))
199 1 : idummy = 0_i4
200 11 : do j = 1, size(pini, 1)
201 11 : if (maske(j)) then
202 10 : idummy = idummy + 1_i4
203 10 : truepara(idummy) = j
204 : end if
205 : end do
206 :
207 : ! Seed random numbers
208 1 : if (iseed == 0) then
209 0 : call date_and_time(values = sdate)
210 : iseed = sdate(1) * 31536000000_i8 + sdate(2) * 2592000000_i8 + sdate(3) * 86400000_i8 + &
211 0 : sdate(5) * 3600000_i8 + sdate(6) * 60000_i8 + sdate(7) * 1000_i8 + sdate(8)
212 0 : call xor4096(iseed, ranval)
213 0 : call xor4096g(iseed, ranval)
214 : else
215 1 : call xor4096(iseed, ranval)
216 1 : call xor4096g(iseed, ranval)
217 : end if
218 :
219 : ! Temporal file writing
220 1 : if(present(tmp_file)) then
221 0 : open(unit = 999, file = trim(adjustl(tmp_file)), action = 'write', status = 'unknown')
222 0 : write(999, *) '# settings :: general'
223 0 : write(999, *) '# nIterations iseed'
224 0 : write(999, *) imaxiter, iseed
225 0 : write(999, *) '# settings :: dds specific'
226 0 : write(999, *) '# dds_r'
227 0 : write(999, *) ir
228 0 : write(999, *) '# iter bestf (bestx(j),j=1,nopt)'
229 0 : close(999)
230 :
231 : ! Second file writing with corresponding parameter set for each iteration rather than the current best
232 0 : open(unit = 998, file = trim(adjustl( tmp_file )) // ".current" , action = 'write', status = 'unknown')
233 0 : write(998, *) '# settings :: general'
234 0 : write(998, *) '# nIterations iseed'
235 0 : write(998, *) imaxiter, iseed
236 0 : write(998, *) '# settings :: dds specific'
237 0 : write(998, *) '# dds_r'
238 0 : write(998, *) ir
239 0 : write(998, *) '# iter bestf (bestx(j),j=1,nopt)'
240 0 : close(998)
241 : end if
242 :
243 : ! Evaluate initial solution and return objective function value
244 : ! and Initialise the other variables (e.g. of_best)
245 : ! imaxit is 1.0 for MIN problems, -1 for MAX problems
246 11 : DDS = pini
247 1 : print *, imaxit
248 : #ifdef MPI
249 : call MPI_Comm_size(comm, nproc, ierror)
250 : do iproc = 1, nproc-1
251 : do_obj_loop = .true.
252 : call MPI_Send(do_obj_loop,1, MPI_LOGICAL,iproc,0,comm,ierror)
253 : end do
254 : #endif
255 1 : of_new = imaxit * obj_func(pini, eval)
256 1 : of_best = of_new
257 1 : if (present(history)) history(1) = of_new
258 :
259 1 : file_write : if (present(tmp_file)) then
260 0 : open(unit = 999, file = trim(adjustl(tmp_file)), action = 'write', position = 'append', recl = (pnum + 2) * 30)
261 : open(unit = 998, file = trim(adjustl( tmp_file )) // ".current" , action = 'write', position = 'append', &
262 0 : recl = (pnum + 2) * 30)
263 0 : if (imaxit .lt. 0.0_dp) then
264 : ! Maximize
265 0 : write(999, *) '0', -of_best, pini
266 0 : write(998, *) '0', -of_best, pini
267 : else
268 : ! Minimize
269 0 : write(999, *) '0', of_best, pini
270 0 : write(998, *) '0', of_best, pini
271 : end if
272 0 : close(999)
273 0 : close(998)
274 : end if file_write
275 :
276 : ! Code below is now the DDS algorithm as presented in Figure 1 of Tolson and Shoemaker (2007)
277 :
278 100000 : do i = 1, imaxiter - 1
279 : ! Determine Decision Variable (DV) selected for perturbation:
280 99999 : Pn = 1.0_dp - log(real(i, dp)) / log(real(imaxiter - 1, dp)) ! probability each DV selected
281 99999 : dvn_count = 0 ! counter for how many DVs selected for perturbation
282 1099989 : pnew = DDS ! define pnew initially as best current solution
283 :
284 : ! Step 3 of Fig 1 of Tolson and Shoemaker (2007)
285 1099989 : do j = 1, size(truepara) !pnum
286 999990 : call xor4096(0_i8, ranval) ! selects next uniform random number in sequence
287 : ! Step 4 of Fig 1 of Tolson and Shoemaker (2007)
288 1099989 : if (ranval < Pn) then ! jth DV selected for perturbation
289 86800 : dvn_count = dvn_count + 1
290 : ! call 1-D perturbation function to get new DV value (new_value)
291 86800 : call neigh_value(DDS(truepara(j)), prange(truepara(j), 1), prange(truepara(j), 2), ir, new_value)
292 86800 : pnew(truepara(j)) = new_value
293 : end if
294 : end do
295 :
296 : ! Step 3 of Fig 1 of Tolson and Shoemaker (2007) in case {N} empty
297 99999 : if (dvn_count == 0) then ! no DVs selected at random, so select one
298 52256 : call xor4096(0_i8, ranval) ! selects next uniform random number in sequence
299 52256 : dv = truepara(int((ranval * real(size(truepara), dp)) + 1.0_dp, i4)) ! index for one DV
300 : ! call 1-D perturbation function to get new DV value (new_value):
301 52256 : call neigh_value(DDS(dv), prange(dv, 1), prange(dv, 2), ir, new_value)
302 52256 : pnew(dv) = new_value ! change relevant DV value in stest
303 : end if
304 :
305 : ! Step 5 of Fig 1 of Tolson and Shoemaker (2007)
306 : ! Evaluate obj function value for test
307 : #ifdef MPI
308 : call MPI_Comm_size(comm, nproc, ierror)
309 : do iproc = 1, nproc-1
310 : do_obj_loop = .true.
311 : call MPI_Send(do_obj_loop,1, MPI_LOGICAL,iproc,0,comm,ierror)
312 : end do
313 : #endif
314 99999 : of_new = imaxit * obj_func(pnew, eval) ! imaxit handles min(=1) and max(=-1) problems
315 : ! update current best solution
316 99999 : if (of_new <= of_best) then
317 45 : of_best = of_new
318 495 : DDS = pnew
319 : end if
320 99999 : if (present(history)) history(i + 1) = of_best
321 :
322 100000 : file_write2 : if (present(tmp_file)) then
323 0 : open(unit = 999, file = trim(adjustl(tmp_file)), action = 'write', position = 'append', recl = (pnum + 2) * 30)
324 : open(unit = 998, file = trim(adjustl( tmp_file )) // ".current" , action = 'write', position = 'append', &
325 0 : recl = (pnum + 2) * 30)
326 0 : if (imaxit .lt. 0.0_dp) then
327 : ! Maximize
328 0 : write(999, *) i, -of_best, dds
329 0 : write(998, *) i, -of_new, pnew
330 : else
331 : ! Minimize
332 0 : write(999, *) i, of_best, dds
333 0 : write(998, *) i, of_new, pnew
334 : end if
335 0 : close(999)
336 0 : close(998)
337 : end if file_write2
338 :
339 : end do
340 1 : if (present(funcbest)) funcbest = of_best
341 : #ifdef MPI
342 : call MPI_Comm_size(comm, nproc, ierror)
343 : do iproc = 1, nproc-1
344 : do_obj_loop = .false.
345 : call MPI_Send(do_obj_loop,1, MPI_LOGICAL,iproc,0,comm,ierror)
346 : end do
347 : #endif
348 :
349 1 : end function DDS
350 :
351 : ! ------------------------------------------------------------------
352 :
353 : !> \brief MDDS
354 :
355 :
356 : !> \details Searches Minimum or Maximum of a user-specified function using the
357 : !! Modified Dynamically Dimensioned Search (DDS).\n
358 : !! DDS is an n-dimensional continuous global optimization algorithm.
359 : !! It is coded as a minimizer but one can give maxit=True in a maximization problem,
360 : !! so that the algorithm minimizes the negative of the objective function F=(-1*F).\n
361 : !! The function to be minimized is the first argument of DDS and must be defined as
362 : !! \code
363 : !! function func(p)
364 : !! use mo_kind, only: dp
365 : !! implicit none
366 : !! real(dp), dimension(:), intent(in) :: p
367 : !! real(dp) :: func
368 : !! end function func
369 : !! \endcode
370 : !!
371 : !! MDDS extents normal DDS by a continuous reduction of the DDS pertubation parameter r from 0.3 to 0.05,
372 : !! and by allowing a larger function value with a certain probablity.
373 :
374 : !> \param[in] "real(dp) :: obj_func(p)" Function on which to search the minimum
375 : !> \param[in] "real(dp) :: pini(:)" inital value of decision variables
376 : !> \param[in] "real(dp) :: prange(size(pini),2)" Min/max range of decision variables
377 : !> \param[in] "integer(i8), optional :: seed" User seed to initialise the random number generator
378 : !! (default: None)
379 : !> \param[in] "integer(i8), optional :: maxiter" Maximum number of iteration or function evaluation
380 : !! (default: 1000)
381 : !> \param[in] "logical, optional :: maxit" Maximization (.True.) or
382 : !! minimization (.False.) of function
383 : !! (default: .False.)
384 : !> \param[in] "logical, optional :: mask(size(pini))" parameter to be optimized (true or false)
385 : !! (default: .True.)
386 : !> \param[in] "character(len=*) , optional :: tmp_file" file with temporal output
387 : !> \param[out] "real(dp), optional :: funcbest" the best value of the function.
388 : !> \param[out] "real(dp), optional, allocatable :: history(:)" the history of best function values,
389 : !! history(maxiter)=funcbest\n
390 : !! allocatable only to be in correspondance with other optimization routines
391 : !> \return real(dp) :: MDDS — The parameters of the point which is estimated to minimize the function.
392 :
393 : !> ## Restrictions
394 : !! None.
395 :
396 : !> ## Example
397 : !!
398 : !! dv_range(:,1) = (/ -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0 /)
399 : !! dv_range(:,2) = (/ 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0 /)
400 : !! dv_ini = (/ -.226265E+01, -.130187E+01, -.151219E+01, 0.133983E+00, 0.988159E+00, &
401 : !! -.495074E+01, -.126574E+02, 0.572684E+00, 0.303864E+01, 0.343031E+01 /)
402 : !! dv_opt = MDDS(griewank, dv_ini, dv_range)
403 : !!
404 : !! See also example in test directory.
405 :
406 : !> ## Literature
407 : !!
408 : !! 1. Tolson, B. A., and C. A. Shoemaker (2007),
409 : !! _Dynamically dimensioned search algorithm for computationally efficient watershed
410 : !! model calibration_, Water Resour. Res., 43, W01413, doi:10.1029/2005WR004723.
411 : !! 2. Huang X-L and Xiong J (2010),
412 : !! _Parameter Optimization of Multi-tank Model with Modified Dynamically Dimensioned
413 : !! Search Algorithm_, Proceedings of the Third International Symposium on Computer
414 : !! Science and Computational Technology(ISCSCT ''10), Jiaozuo, P. R. China,
415 : !! 14-15 August 2010, pp. 283-288\n
416 :
417 : !> \author Written Matthias Cuntz and Juliane Mai
418 : !> \date Aug 2012
419 : ! Modified, Juliane Mai, Nov 2012 - masked parameter
420 : ! Juliane Mai, Dec 2012 - history output
421 :
422 4 : function MDDS(eval, obj_func, pini, prange, seed, maxiter, maxit, mask, tmp_file, funcbest, history)
423 :
424 1 : use mo_kind, only : i4, i8, dp
425 : use mo_xor4096, only : xor4096, xor4096g
426 : use mo_optimization_utils, only : eval_interface, objective_interface
427 :
428 : implicit none
429 :
430 : procedure(eval_interface), INTENT(IN), POINTER :: eval
431 : procedure(objective_interface), intent(in), pointer :: obj_func
432 : real(dp), dimension(:), intent(in) :: pini ! inital value of decision variables
433 : real(dp), dimension(:, :), intent(in) :: prange ! Min/max values of decision variables
434 : integer(i8), optional, intent(in) :: seed ! User seed to initialise the random number generator
435 : integer(i8), optional, intent(in) :: maxiter ! Maximum number of iteration or function evaluation
436 : logical, optional, intent(in) :: maxit ! Maximization or minimization of function
437 : logical, dimension(:), optional, intent(in) :: mask ! parameter to be optimized (true or false)
438 : character(len = *), optional, intent(in) :: tmp_file ! file for temporal output
439 : real(dp), optional, intent(out) :: funcbest ! Best value of the function.
440 : real(dp), dimension(:), optional, intent(out), &
441 : allocatable :: history ! History of objective function values
442 : real(dp), dimension(size(pini)) :: MDDS ! Best value of decision variables
443 :
444 : ! Local variables
445 : integer(i4) :: pnum ! Total number of decision variables
446 : integer(i8) :: iseed ! User given seed
447 : integer(i8) :: imaxiter ! Maximum number of iteration or function evaluation
448 1 : real(dp) :: ir ! MDDS perturbation parameter
449 1 : real(dp) :: imaxit ! Maximization or minimization of function
450 1 : real(dp) :: of_new, of_best ! intermediate results
451 1 : real(dp) :: Pn, new_value ! intermediate results
452 11 : real(dp), dimension(size(pini)) :: pnew ! Test value of decision variables
453 1 : real(dp) :: ranval ! random value
454 : integer(i8) :: i ! maxiter=i8
455 : integer(i4) :: j, dvn_count, dv ! pnum=i4
456 : integer(i4) :: idummy ! dummy vaiable
457 : integer, dimension(8) :: sdate ! date_and_time return
458 1 : logical, dimension(size(pini)) :: maske ! parameter to be optimized (true or false)
459 1 : integer(i4), dimension(:), allocatable :: truepara ! parameter to be optimized (their indexes)
460 :
461 :
462 : ! Check input
463 1 : pnum = size(pini)
464 1 : if (size(prange, 1) /= pnum) stop 'Error MDDS: size(prange,1) /= size(pini)'
465 1 : if (size(prange, 2) /= 2) stop 'Error MDDS: size(prange,2) /= 2'
466 : ! max. iteration
467 1 : imaxiter = 1000
468 1 : if (present(maxiter)) imaxiter = maxiter
469 1 : if (imaxiter < 6) stop 'Error MDDS: max function evals must be minimum 6'
470 : ! history output
471 1 : if (present(history)) then
472 0 : allocate(history(imaxiter))
473 : end if
474 : ! Min or max objective function
475 1 : imaxit = 1.0_dp
476 1 : if (present(maxit)) then
477 1 : if (maxit) imaxit = -1.0_dp
478 : end if
479 : ! Given seed
480 1 : iseed = 0
481 1 : if (present(seed)) iseed = seed
482 1 : iseed = max(iseed, 0_i8)
483 :
484 : ! Seed random numbers
485 1 : if (iseed == 0) then
486 0 : call date_and_time(values = sdate)
487 : iseed = sdate(1) * 31536000000_i8 + sdate(2) * 2592000000_i8 + sdate(3) * 86400000_i8 + &
488 0 : sdate(5) * 3600000_i8 + sdate(6) * 60000_i8 + sdate(7) * 1000_i8 + sdate(8)
489 0 : call xor4096(iseed, ranval)
490 0 : call xor4096g(iseed, ranval)
491 : else
492 1 : call xor4096(iseed, ranval)
493 1 : call xor4096g(iseed, ranval)
494 : end if
495 :
496 : ! Masked parameters
497 1 : if (present(mask)) then
498 0 : if (count(mask) .eq. 0_i4) then
499 0 : stop 'Input argument mask: At least one element has to be true'
500 : else
501 0 : maske = mask
502 : end if
503 : else
504 11 : maske = .true.
505 : end if
506 :
507 13 : allocate (truepara(count(maske)))
508 1 : idummy = 0_i4
509 11 : do j = 1, size(pini, 1)
510 11 : if (maske(j)) then
511 10 : idummy = idummy + 1_i4
512 10 : truepara(idummy) = j
513 : end if
514 : end do
515 :
516 : ! Temporal file writing
517 1 : if(present(tmp_file)) then
518 0 : open(unit = 999, file = trim(adjustl(tmp_file)), action = 'write', status = 'unknown')
519 0 : write(999, *) '# settings :: general'
520 0 : write(999, *) '# nIterations iseed'
521 0 : write(999, *) imaxiter, iseed
522 0 : write(999, *) '# settings :: mdds specific'
523 0 : write(999, *) '# None'
524 0 : write(999, *) ''
525 0 : write(999, *) '# iter bestf (bestx(j),j=1,nopt)'
526 0 : close(999)
527 : end if
528 :
529 : ! Evaluate initial solution and return objective function value
530 : ! and Initialise the other variables (e.g. of_best)
531 : ! imaxit is 1.0 for MIN problems, -1 for MAX problems
532 11 : MDDS = pini
533 1 : of_new = imaxit * obj_func(pini, eval)
534 1 : of_best = of_new
535 1 : if (present(history)) history(1) = of_new
536 :
537 1 : file_write : if (present(tmp_file)) then
538 0 : open(unit = 999, file = trim(adjustl(tmp_file)), action = 'write', position = 'append', recl = (pnum + 2) * 30)
539 0 : if (imaxit .lt. 0.0_dp) then
540 : ! Maximize
541 0 : write(999, *) '0', -of_best, mdds
542 : else
543 : ! Minimize
544 0 : write(999, *) '0', of_best, mdds
545 : end if
546 0 : close(999)
547 : end if file_write
548 :
549 : ! Code below is now the MDDS algorithm as presented in Figure 1 of Tolson and Shoemaker (2007)
550 :
551 100000 : do i = 1, imaxiter - 1
552 : ! Determine Decision Variable (DV) selected for perturbation:
553 99999 : Pn = 1.0_dp - log(real(i, dp)) / log(real(imaxiter - 1, dp)) ! probability each DV selected
554 99999 : dvn_count = 0 ! counter for how many DVs selected for perturbation
555 1099989 : pnew = MDDS ! define pnew initially as best current solution
556 : ! Modifications by Huang et al. (2010)
557 99999 : Pn = max(Pn, 0.05_dp)
558 99999 : ir = max(min(0.3_dp, Pn), 0.05_dp)
559 :
560 : ! Step 3 of Fig 1 of Tolson and Shoemaker (2007)
561 1099989 : do j = 1, pnum
562 999990 : call xor4096(0_i8, ranval) ! selects next uniform random number in sequence
563 : ! Step 4 of Fig 1 of Tolson and Shoemaker (2007)
564 1099989 : if (ranval < Pn) then ! jth DV selected for perturbation
565 98852 : dvn_count = dvn_count + 1
566 : ! call 1-D perturbation function to get new DV value (new_value)
567 98852 : call neigh_value(MDDS(truepara(j)), prange(truepara(j), 1), prange(truepara(j), 2), ir, new_value)
568 98852 : pnew(truepara(j)) = new_value
569 : end if
570 : end do
571 :
572 : ! Step 3 of Fig 1 of Tolson and Shoemaker (2007) in case {N} empty
573 99999 : if (dvn_count == 0) then ! no DVs selected at random, so select one
574 43387 : call xor4096(0_i8, ranval) ! selects next uniform random number in sequence
575 43387 : dv = truepara(int((ranval * real(size(truepara), dp)) + 1.0_dp, i4)) ! index for one DV
576 : ! call 1-D perturbation function to get new DV value (new_value):
577 43387 : call neigh_value(MDDS(dv), prange(dv, 1), prange(dv, 2), ir, new_value)
578 43387 : pnew(dv) = new_value ! change relevant DV value in stest
579 : end if
580 :
581 : ! Step 5 of Fig 1 of Tolson and Shoemaker (2007)
582 : ! Evaluate obj function value for test
583 99999 : of_new = imaxit * obj_func(pnew, eval) ! imaxit handles min(=1) and max(=-1) problems
584 : ! update current best solution
585 99999 : if (of_new <= of_best) then
586 500 : of_best = of_new
587 5500 : MDDS = pnew
588 : else ! Modifications by Huang et al. (2010)
589 99499 : call xor4096(0_i8, ranval)
590 99499 : if (exp(-(of_new - of_best) / of_best) > (1.0_dp - ranval * Pn)) then
591 502 : of_best = of_new
592 5522 : MDDS = pnew
593 : end if
594 : end if
595 99999 : if (present(history)) then
596 0 : if (present(maxit)) then
597 0 : if (maxit) then
598 0 : history(i + 1) = max(history(i), of_best)
599 : else
600 0 : history(i + 1) = min(history(i), of_best)
601 : end if
602 : else
603 0 : history(i + 1) = min(history(i), of_best)
604 : end if
605 : end if
606 :
607 100000 : file_write2 : if (present(tmp_file)) then
608 0 : open(unit = 999, file = trim(adjustl(tmp_file)), action = 'write', position = 'append', recl = (pnum + 2) * 30)
609 0 : if (imaxit .lt. 0.0_dp) then
610 : ! Maximize
611 0 : write(999, *) i, -of_best, mdds
612 : else
613 : ! Minimize
614 0 : write(999, *) i, of_best, mdds
615 : end if
616 0 : close(999)
617 : end if file_write2
618 :
619 : end do
620 1 : if (present(funcbest)) funcbest = of_best
621 :
622 1 : end function MDDS
623 :
624 : ! ------------------------------------------------------------------
625 :
626 : ! Purpose is to generate a neighboring decision variable value for a single
627 : ! decision variable value being perturbed by the DDS optimization algorithm.
628 : ! New DV value respects the upper and lower DV bounds.
629 : ! Coded by Bryan Tolson, Nov 2005.
630 :
631 : ! I/O variable definitions:
632 : ! x_cur - current decision variable (DV) value
633 : ! x_min - min DV value
634 : ! x_max - max DV value
635 : ! r - the neighborhood perturbation factor
636 : ! new_value - new DV variable value (within specified min and max)
637 :
638 281295 : subroutine neigh_value(x_cur, x_min, x_max, r, new_value)
639 :
640 1 : use mo_kind, only : i8, dp
641 : use mo_xor4096, only : xor4096g
642 :
643 : implicit none
644 :
645 : real(dp), intent(in) :: x_cur, x_min, x_max, r
646 : real(dp), intent(out) :: new_value
647 281295 : real(dp) :: x_range
648 281295 : real(dp) :: zvalue
649 :
650 281295 : x_range = x_max - x_min
651 :
652 : ! generate a standard normal random variate (zvalue)
653 281295 : call xor4096g(0_i8, zvalue)
654 :
655 : ! calculate new decision variable value:
656 281295 : new_value = x_cur + zvalue * r * x_range
657 :
658 : ! check new value is within bounds. If not, bounds are reflecting.
659 281295 : if (new_value < x_min) then
660 1777 : new_value = x_min + (x_min - new_value)
661 1777 : if (new_value > x_max) then
662 : ! if reflection goes past x_max then value should be x_min since
663 : ! without reflection the approach goes way past lower bound.
664 : ! This keeps x close to lower bound when x_cur is close to lower bound
665 : ! Practically speaking, this should never happen with r values <0.3.
666 0 : new_value = x_min
667 : end if
668 279518 : else if (new_value > x_max) then
669 1837 : new_value = x_max - (new_value - x_max)
670 1837 : if (new_value < x_min) then
671 : ! if reflection goes past x_min then value should be x_max for same reasons as above.
672 0 : new_value = x_max
673 : end if
674 : end if
675 :
676 281295 : end subroutine neigh_value
677 :
678 : ! ------------------------------------------------------------------
679 :
680 : end module mo_dds
|
