diff --git a/src/mo_anneal.f90 b/src/mo_anneal.f90 index d15abec6b235dea4e106d56749ed2ab116480461..7b1ab73a63c20b5790f12a30fb592c80b515da02 100644 --- a/src/mo_anneal.f90 +++ b/src/mo_anneal.f90 @@ -44,7 +44,7 @@ MODULE mo_anneal USE mo_kind, ONLY : i4, i8, dp USE mo_utils, ONLY : le, ge USE mo_xor4096, ONLY : get_timeseed, xor4096, xor4096g, n_save_state - USE mo_optimization_utils, ONLY : eval_interface, objective_interface + USE mo_optimization_utils, ONLY : optimizee IMPLICIT NONE @@ -61,15 +61,22 @@ MODULE mo_anneal !! !! \b Example !! - !! User defined function `cost_dp` which calculates the cost function value for a + !! User defined extended type of optimizee `objective` which calculates the cost function value for a !! parameter set (the interface given below has to be used for this function!). !! !! \code{.f90} + !! use mo_optimization_utils, only: optimizee + !! type, extends(optimizee) :: your_optimizee + !! contains + !! ! to be implemented with signature: your_evaluate(self, parameters, sigma, stddev_new, likeli_new) + !! procedure :: evaluate => your_evaluate + !! end type function_optimizee + !! type(your_optimizee) :: objective !! para = (/ 1.0_dp , 2.0_dp /) !! prange(1,:) = (/ 0.0_dp, 10.0_dp /) !! prange(2,:) = (/ 0.0_dp, 50.0_dp /) !! - !! parabest = anneal(cost_dp, para, prange) + !! parabest = anneal(objective, para, prange) !! \endcode !! !! See also test folder for a detailed example, "test/test_mo_anneal". @@ -85,8 +92,7 @@ MODULE mo_anneal !! _Dynamically dimensioned search algorithm for computationally efficient watershed model calibration_. !! WRR, 43(1), W01413, 2007. !! - !> \param[in] "INTERFACE :: eval" Interface calculating the eval function at a given point. - !> \param[in] "INTERFACE :: cost" Interface calculating the cost function at a given point. + !> \param[inout] "class(optimizee) :: objective" cost function to search the optimum !> \param[in] "real(dp), dimension(:) :: para" Initial parameter set. !> \param[in] "real(dp), dimension(size(para),2), optional :: prange" !! Lower and upper bound per parameter. @@ -216,9 +222,7 @@ MODULE mo_anneal !! Comput. Opt. and App. (2004). !! !> \param[in] "real(dp), dimension(:) :: paraset" Initial (valid) parameter set. - !> \param[in] "INTERFACE :: eval" Interface calculating the eval function at a given point. - !> \param[in] "INTERFACE :: cost" Interface calculating the cost function at a given point. - !> \param[in] "INTERFACE :: prange_func" Interface for functional ranges. + !> \param[inout] "class(optimizee) :: objective" cost function to search the optimum !> \param[in] "real(dp) :: acc_goal" Acceptance Ratio which has to be achieved. !> \param[in] "real(dp), dimension(size(para),2), optional :: prange" !! Lower and upper bound per parameter. @@ -286,7 +290,7 @@ MODULE mo_anneal CONTAINS - FUNCTION anneal_dp(eval, cost, para, & ! obligatory + FUNCTION anneal_dp(objective, para, & ! obligatory prange, prange_func, & ! optional IN: exactly one of both temp, Dt, nITERmax, Len, nST, eps, acc, & ! optional IN seeds, printflag, maskpara, weight, & ! optional IN @@ -298,9 +302,7 @@ CONTAINS result(parabest) IMPLICIT NONE - - procedure(eval_interface), INTENT(IN), POINTER :: eval - procedure(objective_interface), intent(in), pointer :: cost + class(optimizee), intent(inout) :: objective INTERFACE SUBROUTINE prange_func(paraset, iPar, rangePar) @@ -542,24 +544,24 @@ CONTAINS else if (present(prange_func)) then if (present(undef_funcval)) then - T_in = GetTemperature(eval, cost, para, 0.95_dp, prange_func = prange_func, & + T_in = GetTemperature(objective, para, 0.95_dp, prange_func = prange_func, & maskpara = maskpara_in, samplesize = 2_i4 * LEN_IN, & seeds = seeds_in(1 : 2), printflag = printflag_in, & maxit = ldummy, undef_funcval = undef_funcval) else - T_in = GetTemperature(eval, cost, para, 0.95_dp, prange_func = prange_func, & + T_in = GetTemperature(objective, para, 0.95_dp, prange_func = prange_func, & maskpara = maskpara_in, samplesize = 2_i4 * LEN_IN, & seeds = seeds_in(1 : 2), printflag = printflag_in, & maxit = ldummy) end if else if (present(undef_funcval)) then - T_in = GetTemperature(eval, cost, para, 0.95_dp, prange = prange, & + T_in = GetTemperature(objective, para, 0.95_dp, prange = prange, & maskpara = maskpara_in, samplesize = 2_i4 * LEN_IN, & seeds = seeds_in(1 : 2), printflag = printflag_in, & maxit = ldummy, undef_funcval = undef_funcval) else - T_in = GetTemperature(eval, cost, para, 0.95_dp, prange = prange, & + T_in = GetTemperature(objective, para, 0.95_dp, prange = prange, & maskpara = maskpara_in, samplesize = 2_i4 * LEN_IN, & seeds = seeds_in(1 : 2), printflag = printflag_in, & maxit = ldummy) @@ -650,7 +652,7 @@ CONTAINS DT0 = DT_IN ! Generate and evaluate the initial solution state - fo = cost(gamma(:)%old, eval) * maxit_in + fo = objective%evaluate(gamma(:)%old) * maxit_in if (abs(fo) .lt. tiny(0.0_dp)) fo = 0.0000001_dp * maxit_in file_write : if (present(tmp_file)) then @@ -803,7 +805,7 @@ CONTAINS end select ! (2) Calculate new objective function value - fn = cost(gamma(:)%new, eval) * maxit_in + fn = objective%evaluate(gamma(:)%new) * maxit_in coststatus = .true. if (present(undef_funcval)) then if (abs(fn * maxit_in - undef_funcval) .lt. tiny(1.0_dp)) then @@ -959,7 +961,7 @@ CONTAINS ! calculate cost function again (only for check and return values) parabest = gamma(:)%best - costbest = cost(parabest, eval) * maxit_in + costbest = objective%evaluate(parabest) * maxit_in if (present (funcbest)) then funcbest = costbest * maxit_in end if @@ -988,7 +990,7 @@ CONTAINS END FUNCTION anneal_dp - real(dp) function GetTemperature_dp(eval, cost, paraset, acc_goal, & ! obligatory + real(dp) function GetTemperature_dp(objective, paraset, acc_goal, & ! obligatory prange, prange_func, & ! optional IN: exactly one of both samplesize, maskpara, seeds, printflag, & ! optional IN weight, maxit, undef_funcval & ! optional IN @@ -996,8 +998,7 @@ CONTAINS use mo_kind, only : dp, i4, i8 implicit none - procedure(eval_interface), INTENT(IN), POINTER :: eval - procedure(objective_interface), intent(in), pointer :: cost + class(optimizee), intent(inout) :: objective real(dp), dimension(:), intent(in) :: paraset !< a valid parameter set of the model real(dp), intent(in) :: acc_goal !< acceptance ratio to achieve real(dp), optional, dimension(size(paraset, 1), 2), intent(in) :: prange !< lower and upper limit per parameter @@ -1169,7 +1170,7 @@ CONTAINS gamma(:)%best = paraset(:) NormPhi = -9999.9_dp - fo = cost(paraset, eval) * maxit_in + fo = objective%evaluate(paraset) * maxit_in if (abs(fo) .lt. tiny(0.0_dp)) fo = 0.0000001_dp * maxit_in NormPhi = fo fo = fo / NormPhi * maxit_in @@ -1201,7 +1202,7 @@ CONTAINS gamma(iPar)%min, gamma(iPar)%max, RN2) ! ! (2) Calculate new objective function value and normalize it - fn = cost(gamma(:)%new, eval) * maxit_in + fn = objective%evaluate(gamma(:)%new) * maxit_in coststatus = .true. if (present(undef_funcval)) then if (abs(fn * maxit_in - undef_funcval) .lt. tiny(1.0_dp)) then diff --git a/src/mo_dds.F90 b/src/mo_dds.F90 index 8757c69bf1f067813e3b96e6015e1816107f513b..1b0ff39408fe19f3ff3ebbb86a1611d271009a7e 100644 --- a/src/mo_dds.F90 +++ b/src/mo_dds.F90 @@ -12,6 +12,8 @@ !! FORCES is released under the LGPLv3+ license \license_note module mo_dds + use mo_optimization_utils, only : optimizee + IMPLICIT NONE PRIVATE @@ -35,22 +37,25 @@ CONTAINS !! !! The function to be minimized is the first argument of DDS and must be defined as \n !! \code{.f90} - !! function func(p) - !! use mo_kind, only: dp - !! implicit none - !! real(dp), dimension(:), intent(in) :: p - !! real(dp) :: func - !! end function func + !! use mo_optimization_utils, only: optimizee + !! type, extends(optimizee) :: your_optimizee + !! contains + !! ! to be implemented with signature: your_evaluate(self, parameters, sigma, stddev_new, likeli_new) + !! procedure :: evaluate => your_evaluate + !! end type function_optimizee !! \endcode !! !! \b Example !! !! \code{.f90} + !! use mo_opt_functions, only: griewank + !! use mo_optimization_utils, only: function_optimizee !! dv_range(:,1) = (/ -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0 /) !! dv_range(:,2) = (/ 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0 /) !! dv_ini = (/ -.226265E+01, -.130187E+01, -.151219E+01, 0.133983E+00, 0.988159E+00, & !! -.495074E+01, -.126574E+02, 0.572684E+00, 0.303864E+01, 0.343031E+01 /) - !! dv_opt = DDS(griewank, dv_ini, dv_range) + !! objective%func_pointer => griewank + !! dv_opt = DDS(objective, dv_ini, dv_range) !! \endcode !! !! See also example in test directory. @@ -60,7 +65,7 @@ CONTAINS !! _Dynamically dimensioned search algorithm for computationally efficient watershed !! model calibration_, Water Resour. Res., 43, W01413, doi:10.1029/2005WR004723. !! - !> \param[in] "real(dp) :: obj_func(p)" Function on which to search the minimum + !> \param[inout] "class(optimizee) :: objective" Objective to search the optimum !> \param[in] "real(dp) :: pini(:)" inital value of decision variables !> \param[in] "real(dp) :: prange(size(pini),2)" Min/max range of decision variables !> \param[in] "real(dp), optional :: r" DDS perturbation parameter\n @@ -104,22 +109,19 @@ CONTAINS !! - Added "dds_results.out.current" output file with current iteration (non-updated) parameters #ifdef MPI - function DDS(eval, obj_func, pini, prange, r, seed, maxiter, maxit, mask, tmp_file, comm, funcbest, history) + function DDS(objective, pini, prange, r, seed, maxiter, maxit, mask, tmp_file, comm, funcbest, history) #else - function DDS(eval, obj_func, pini, prange, r, seed, maxiter, maxit, mask, tmp_file, funcbest, history) + function DDS(objective, pini, prange, r, seed, maxiter, maxit, mask, tmp_file, funcbest, history) #endif use mo_kind, only : i4, i8, dp use mo_xor4096, only : xor4096, xor4096g - use mo_optimization_utils, only : eval_interface, objective_interface #ifdef MPI use mpi_f08 #endif implicit none - - procedure(eval_interface), INTENT(IN), POINTER :: eval - procedure(objective_interface), intent(in), pointer :: obj_func + class(optimizee), intent(inout) :: objective real(dp), dimension(:), intent(in) :: pini ! inital value of decision variables real(dp), dimension(:, :), intent(in) :: prange ! Min/max values of decision variables @@ -252,7 +254,7 @@ CONTAINS call MPI_Send(do_obj_loop,1, MPI_LOGICAL,iproc,0,comm,ierror) end do #endif - of_new = imaxit * obj_func(pini, eval) + of_new = imaxit * objective%evaluate(pini) of_best = of_new if (present(history)) history(1) = of_new @@ -311,7 +313,7 @@ CONTAINS call MPI_Send(do_obj_loop,1, MPI_LOGICAL,iproc,0,comm,ierror) end do #endif - of_new = imaxit * obj_func(pnew, eval) ! imaxit handles min(=1) and max(=-1) problems + of_new = imaxit * objective%evaluate(pnew) ! imaxit handles min(=1) and max(=-1) problems ! update current best solution if (of_new <= of_best) then of_best = of_new @@ -359,19 +361,19 @@ CONTAINS !! It is coded as a minimizer but one can give maxit=True in a maximization problem, !! so that the algorithm minimizes the negative of the objective function F=(-1*F).\n !! The function to be minimized is the first argument of DDS and must be defined as - !! \code - !! function func(p) - !! use mo_kind, only: dp - !! implicit none - !! real(dp), dimension(:), intent(in) :: p - !! real(dp) :: func - !! end function func + !! \code{.f90} + !! use mo_optimization_utils, only: optimizee + !! type, extends(optimizee) :: your_optimizee + !! contains + !! ! to be implemented with signature: your_evaluate(self, parameters, sigma, stddev_new, likeli_new) + !! procedure :: evaluate => your_evaluate + !! end type function_optimizee !! \endcode !! !! MDDS extents normal DDS by a continuous reduction of the DDS pertubation parameter r from 0.3 to 0.05, !! and by allowing a larger function value with a certain probablity. - !> \param[in] "real(dp) :: obj_func(p)" Function on which to search the minimum + !> \param[inout] "class(optimizee) :: objective" Objective to search the optimum !> \param[in] "real(dp) :: pini(:)" inital value of decision variables !> \param[in] "real(dp) :: prange(size(pini),2)" Min/max range of decision variables !> \param[in] "integer(i8), optional :: seed" User seed to initialise the random number generator @@ -395,11 +397,16 @@ CONTAINS !> ## Example !! - !! dv_range(:,1) = (/ -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0 /) - !! dv_range(:,2) = (/ 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0 /) - !! dv_ini = (/ -.226265E+01, -.130187E+01, -.151219E+01, 0.133983E+00, 0.988159E+00, & - !! -.495074E+01, -.126574E+02, 0.572684E+00, 0.303864E+01, 0.343031E+01 /) - !! dv_opt = MDDS(griewank, dv_ini, dv_range) + !! \code{.f90} + !! use mo_opt_functions, only: griewank + !! use mo_optimization_utils, only: function_optimizee + !! dv_range(:,1) = (/ -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0 /) + !! dv_range(:,2) = (/ 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0 /) + !! dv_ini = (/ -.226265E+01, -.130187E+01, -.151219E+01, 0.133983E+00, 0.988159E+00, & + !! -.495074E+01, -.126574E+02, 0.572684E+00, 0.303864E+01, 0.343031E+01 /) + !! objective%func_pointer => griewank + !! dv_opt = MDDS(objective, dv_ini, dv_range) + !! \endcode !! !! See also example in test directory. @@ -419,16 +426,14 @@ CONTAINS ! Modified, Juliane Mai, Nov 2012 - masked parameter ! Juliane Mai, Dec 2012 - history output - function MDDS(eval, obj_func, pini, prange, seed, maxiter, maxit, mask, tmp_file, funcbest, history) + function MDDS(objective, pini, prange, seed, maxiter, maxit, mask, tmp_file, funcbest, history) use mo_kind, only : i4, i8, dp use mo_xor4096, only : xor4096, xor4096g - use mo_optimization_utils, only : eval_interface, objective_interface implicit none - procedure(eval_interface), INTENT(IN), POINTER :: eval - procedure(objective_interface), intent(in), pointer :: obj_func + class(optimizee), intent(inout) :: objective real(dp), dimension(:), intent(in) :: pini ! inital value of decision variables real(dp), dimension(:, :), intent(in) :: prange ! Min/max values of decision variables integer(i8), optional, intent(in) :: seed ! User seed to initialise the random number generator @@ -530,7 +535,7 @@ CONTAINS ! and Initialise the other variables (e.g. of_best) ! imaxit is 1.0 for MIN problems, -1 for MAX problems MDDS = pini - of_new = imaxit * obj_func(pini, eval) + of_new = imaxit * objective%evaluate(pini) of_best = of_new if (present(history)) history(1) = of_new @@ -580,7 +585,7 @@ CONTAINS ! Step 5 of Fig 1 of Tolson and Shoemaker (2007) ! Evaluate obj function value for test - of_new = imaxit * obj_func(pnew, eval) ! imaxit handles min(=1) and max(=-1) problems + of_new = imaxit * objective%evaluate(pnew) ! imaxit handles min(=1) and max(=-1) problems ! update current best solution if (of_new <= of_best) then of_best = of_new diff --git a/src/mo_mcmc.F90 b/src/mo_mcmc.F90 index ab4e12a1dbb784353b24104a95a88abc0ab6ca4e..5a5c522786699410b5c6d253418a05ec2d115a22 100644 --- a/src/mo_mcmc.F90 +++ b/src/mo_mcmc.F90 @@ -15,7 +15,7 @@ MODULE mo_mcmc USE mo_append, only : append USE mo_moment, only : stddev !$ USE omp_lib, only: OMP_GET_NUM_THREADS - use mo_optimization_utils, only : eval_interface, objective_interface + use mo_optimization_utils, only : optimizee use mo_message, only : error_message #ifdef FORCES_WITH_NETCDF use mo_ncwrite, only : dump_netcdf @@ -132,7 +132,7 @@ MODULE mo_mcmc !! \b Example !! !! \code{.f90} - !! call mcmc( likelihood, para, rangePar, mcmc_paras, burnin_paras, & + !! call mcmc( objective, para, rangePar, mcmc_paras, burnin_paras, & !! seed_in=seeds, printflag_in=printflag, maskpara_in=maskpara, & !! tmp_file=tmp_file, loglike_in=loglike, & !! ParaSelectMode_in=ParaSelectMode, & @@ -153,8 +153,7 @@ MODULE mo_mcmc !! 90-103. doi:10.1016/j.ecolmodel.2012.03.009. !! - !> \param[in] "real(dp) :: likelihood(x)" Interface Function which calculates likelihood - !! of given parameter set x + !> \param[inout] "class(optimizee) :: objective" Likelihood objective to search the optimum !> \param[in] "real(dp) :: para(:)" Inital parameter set (should be GOOD approximation !! of best parameter set) !> \param[in] "real(dp) :: rangePar(size(para),2)" Min/max range of parameters @@ -194,7 +193,7 @@ MODULE mo_mcmc !> \param[out] "real(dp), allocatable :: burnin_paras(:,:)" Parameter sets sampled during burn-in part of algorithm !> \note - !> Likelihood has to be defined as a function interface\n + !> Likelihood objective has to be implented as an extended type of optimizee\n !> The maximal number of parameters is 1000. !> \authors Juliane Mai @@ -306,7 +305,7 @@ MODULE mo_mcmc !! \b Example !! !! \code{.f90} - !! call mcmc( likelihood, para, rangePar, mcmc_paras, burnin_paras, & + !! call mcmc( objective, para, rangePar, mcmc_paras, burnin_paras, & !! seed_in=seeds, printflag_in=printflag, maskpara_in=maskpara, & !! tmp_file=tmp_file, loglike_in=loglike, & !! ParaSelectMode_in=ParaSelectMode, & @@ -327,8 +326,8 @@ MODULE mo_mcmc !! 90-103. doi:10.1016/j.ecolmodel.2012.03.009. !! - !> \param[in] "real(dp) :: likelihood(x,sigma,stddev_new,likeli_new)" - !> Interface Function which calculates likelihood + !> \param[inout] "class(optimizee) :: objective" Likelihood objective to search the optimum, with + !> type-bound procedure 'evaluate' which calculates likelihood !> of given parameter set x and given standard deviation sigma !> and returns optionally the standard deviation stddev_new !> of the errors using x and @@ -430,7 +429,7 @@ CONTAINS !----------------------------------------------------------------------------------------------- - SUBROUTINE mcmc_dp(eval, likelihood, para, rangePar, & ! obligatory IN + SUBROUTINE mcmc_dp(objective, para, rangePar, & ! obligatory IN mcmc_paras, burnin_paras, & ! obligatory OUT seed_in, printflag_in, maskpara_in, & ! optional IN restart, restart_file, & ! optional IN: if mcmc is restarted and file which contains restart variables @@ -445,8 +444,7 @@ CONTAINS IMPLICIT NONE - procedure(eval_interface), INTENT(IN), POINTER :: eval - procedure(objective_interface), intent(in), pointer :: likelihood + class(optimizee), intent(inout) :: objective REAL(DP), DIMENSION(:, :), INTENT(IN) :: rangePar ! range for each parameter REAL(DP), DIMENSION(:), INTENT(IN) :: para ! initial parameter i @@ -734,7 +732,7 @@ CONTAINS parabest = para ! initialize likelihood - likelibest = likelihood(parabest, eval) + likelibest = objective%evaluate(parabest) !---------------------------------------------------------------------- ! (1) BURN IN @@ -798,7 +796,7 @@ CONTAINS paranew, ChangePara) ! (B) new likelihood - likelinew = likelihood(paranew, eval) + likelinew = objective%evaluate(paranew) oddsSwitch1 = .false. if (loglike) then @@ -1039,7 +1037,7 @@ CONTAINS else paraold = mcmc_paras_3d(Ipos(chain) + Ineg(chain), :, chain) end if - likeliold = likelihood(paraold, eval) + likeliold = objective%evaluate(paraold) markovchainMCMC : do @@ -1053,7 +1051,7 @@ CONTAINS paranew, ChangePara) ! (B) new likelihood - likelinew = likelihood(paranew, eval) + likelinew = objective%evaluate(paranew) oddsSwitch1 = .false. if (loglike) then oddsRatio = likelinew - likeliold @@ -1286,7 +1284,7 @@ CONTAINS RETURN END SUBROUTINE mcmc_dp - SUBROUTINE mcmc_stddev_dp(eval, likelihood, para, rangePar, & ! obligatory IN + SUBROUTINE mcmc_stddev_dp(objective, para, rangePar, & ! obligatory IN mcmc_paras, burnin_paras, & ! obligatory OUT seed_in, printflag_in, maskpara_in, & ! optional IN tmp_file, & ! optional IN : filename for temporal output of @@ -1300,8 +1298,7 @@ CONTAINS IMPLICIT NONE - procedure(eval_interface), INTENT(IN), POINTER :: eval - procedure(objective_interface), intent(in), pointer :: likelihood + class(optimizee), intent(inout) :: objective REAL(DP), DIMENSION(:, :), INTENT(IN) :: rangePar ! range for each parameter REAL(DP), DIMENSION(:), INTENT(IN) :: para ! initial parameter i @@ -1533,7 +1530,7 @@ CONTAINS ! write(*,*) parabest ! initialize likelihood and sigma - likelibest = likelihood(parabest, eval, 1.0_dp, stddev_new, likeli_new) + likelibest = objective%evaluate(parabest, 1.0_dp, stddev_new, likeli_new) likelibest = likeli_new stddev_data = stddev_new @@ -1612,7 +1609,7 @@ CONTAINS paranew, ChangePara) ! (B) new likelihood - likelinew = likelihood(paranew, eval, stddev_data, stddev_new, likeli_new) + likelinew = objective%evaluate(paranew, stddev_data, stddev_new, likeli_new) oddsSwitch1 = .false. if (loglike) then @@ -1794,7 +1791,7 @@ CONTAINS else paraold = mcmc_paras_3d(Ipos(chain) + Ineg(chain), :, chain) end if - likeliold = likelihood(paraold, eval, stddev_data) + likeliold = objective%evaluate(paraold, stddev_data) markovchainMCMC : do @@ -1808,7 +1805,7 @@ CONTAINS paranew, ChangePara) ! (B) new likelihood - likelinew = likelihood(paranew, eval, stddev_data) + likelinew = objective%evaluate(paranew, stddev_data) oddsSwitch1 = .false. if (loglike) then oddsRatio = likelinew - likeliold diff --git a/src/mo_optimization_utils.f90 b/src/mo_optimization_utils.f90 index 57d1939e9fba20494cd08abe941bb5fc234319dc..d7e48755af6cdbcb14282bbad2854c49cfbf9528 100644 --- a/src/mo_optimization_utils.f90 +++ b/src/mo_optimization_utils.f90 @@ -8,6 +8,7 @@ module mo_optimization_utils use mo_kind, only : dp + use mo_message, only : error_message implicit none @@ -24,7 +25,7 @@ module mo_optimization_utils abstract interface subroutine eval_interface(config, opti_sim) use mo_optimization_types, only : config_t, opti_sim_t - type(config_t), intent(in) :: config + type(config_t), intent(in) :: config type(opti_sim_t), dimension(:), pointer, optional, intent(inout) :: opti_sim end subroutine end interface @@ -33,7 +34,7 @@ module mo_optimization_utils !> \details The optional arguments are motivated by likelihood objective functions. ! ToDo: eval optional interface - function objective_interface (parameterset, eval, arg1, arg2, arg3) + function objective_interface(parameterset, eval, arg1, arg2, arg3) use mo_kind, only : dp import eval_interface real(dp), intent(in), dimension(:) :: parameterset !< parameter set @@ -46,4 +47,134 @@ module mo_optimization_utils end function objective_interface end interface + + !> \brief abstract type 'optimizee' to be used by optimizers + type, abstract :: optimizee + contains + !> \brief evaluate the optimizee + procedure(evaluate_interface), deferred :: evaluate + end type optimizee + + !> \brief Abstract interface for deferred procedures + abstract interface + function evaluate_interface(self, parameters, sigma, stddev_new, likeli_new) result(value) + use mo_kind, only : dp + import :: optimizee + class(optimizee), intent(inout) :: self + real(dp), intent(in), dimension(:) :: parameters !< parameter set + real(dp), intent(in), optional :: sigma !< likelihood: standard deviation of data + real(dp), intent(out), optional :: stddev_new !< likelihood: standard deviation of errors using parameters + real(dp), intent(out), optional :: likeli_new !< likelihood: likelihood using stddev_new + real(dp) :: value !< output value + end function evaluate_interface + end interface + + !> \brief Optimizee for a simple function: f(vec) + type, extends(optimizee) :: function_optimizee + procedure(func_interface), pointer, nopass :: func_pointer => null() !< Pointer to the function + contains + procedure :: evaluate => evaluate_function + end type function_optimizee + + !> \brief Abstract interface for the function pointer + abstract interface + function func_interface(parameters) result(value) + use mo_kind, only : dp + real(dp), dimension(:), intent(in) :: parameters !< parameter set + real(dp) :: value !< output value + end function func_interface + end interface + + !> \brief Optimizee for a likelihood function: f(vec, sigma, stddev_new, likeli_new) + type, extends(optimizee) :: likelihood_optimizee + procedure(likelihood_interface), pointer, nopass :: likelihood_pointer => null() !< Pointer to the likelihood + contains + procedure :: evaluate => evaluate_likelihood + end type likelihood_optimizee + + !> \brief Abstract interface for the likelihood pointer + abstract interface + function likelihood_interface(parameters, sigma, stddev_new, likeli_new) result(value) + use mo_kind, only : dp + real(dp), dimension(:), intent(in) :: parameters !< parameter set + real(dp), intent(in), optional :: sigma !< standard deviation of data + real(dp), intent(out), optional :: stddev_new !< standard deviation of errors using parameters + real(dp), intent(out), optional :: likeli_new !< likelihood using stddev_new, + real(dp) :: value !< output value + end function likelihood_interface + end interface + + !> \brief Optimizee for a eval-objective pair + type, extends(optimizee) :: eval_optimizee + procedure(eval_interface), pointer, nopass :: eval_pointer => null() !< Pointer to the eval + procedure(objective_interface), pointer, nopass :: obj_pointer => null() !< Pointer to the objective + contains + procedure :: evaluate => evaluate_obj_eval + end type eval_optimizee + +contains + + !> \brief Implementation of the evaluate procedure for a simple function + function evaluate_function(self, parameters, sigma, stddev_new, likeli_new) result(value) + class(function_optimizee), intent(inout) :: self + real(dp), dimension(:), intent(in) :: parameters !< parameter set + real(dp), intent(in), optional :: sigma !< likelihood: standard deviation of data + real(dp), intent(out), optional :: stddev_new !< likelihood: standard deviation of errors using parameters + real(dp), intent(out), optional :: likeli_new !< likelihood: likelihood using stddev_new + real(dp) :: value !< output value + + ! Ensure the function pointer is set + if (.not. associated(self%func_pointer)) then + call error_message("Function pointer is not set in function_optimizee!") + end if + + ! Check optional arguments + if (present(sigma) .or. present(stddev_new) .or. present(likeli_new)) then + call error_message("function_optimizee doesn't support 'sigma', 'stddev_new' or 'likeli_new'!") + end if + + ! Call the function pointer + value = self%func_pointer(parameters) + end function evaluate_function + + !> \brief Implementation of the evaluate procedure for a likelihood function + function evaluate_likelihood(self, parameters, sigma, stddev_new, likeli_new) result(value) + class(likelihood_optimizee), intent(inout) :: self + real(DP), dimension(:), intent(in) :: parameters + real(DP), intent(in), optional :: sigma + real(DP), intent(out), optional :: stddev_new + real(DP), intent(out), optional :: likeli_new + real(DP) :: value + + ! Ensure the likelihood function pointer is set + if (.not. associated(self%likelihood_pointer)) then + call error_message("Likelihood function pointer is not set in likelihood_optimizee!") + end if + + ! Call the likelihood function pointer + value = self%likelihood_pointer(parameters, sigma, stddev_new, likeli_new) + end function evaluate_likelihood + + !> \brief Implementation of the evaluate procedure for a eval-objective pair + function evaluate_obj_eval(self, parameters, sigma, stddev_new, likeli_new) result(value) + class(eval_optimizee), intent(inout) :: self + real(DP), dimension(:), intent(in) :: parameters + real(DP), intent(in), optional :: sigma + real(DP), intent(out), optional :: stddev_new + real(DP), intent(out), optional :: likeli_new + real(DP) :: value + + ! Ensure the eval function pointer is set + if (.not. associated(self%eval_pointer)) then + call error_message("Eval function pointer is not set in eval_optimizee!") + end if + ! Ensure the objective function pointer is set + if (.not. associated(self%obj_pointer)) then + call error_message("Objective function pointer is not set in eval_optimizee!") + end if + + ! Call the objective function pointer + value = self%obj_pointer(parameters, self%eval_pointer, sigma, stddev_new, likeli_new) + end function evaluate_obj_eval + end module mo_optimization_utils diff --git a/src/mo_sce.F90 b/src/mo_sce.F90 index 484d0640ad1729e28cda3658dfe0d77c3204bca6..6bd015881c23fc8e8695ab8cf5ef1a8c2beb7407 100644 --- a/src/mo_sce.F90 +++ b/src/mo_sce.F90 @@ -10,6 +10,7 @@ !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved. !! FORCES is released under the LGPLv3+ license \license_note MODULE mo_sce + use mo_optimization_utils, only : optimizee IMPLICIT NONE @@ -48,26 +49,27 @@ MODULE mo_sce !! Rainfall-runoff Models' by Duan, Q., S. Sorooshian, and V.K. Gupta, !! Water Resources Research, Vol 28(4), pp.1015-1031, 1992.\n !! \n - !! The function to be minimized is the first argument of DDS and must be defined as \n + !! The function to be minimized is the first argument of SCE and must be defined as \n !! \code{.f90} - !! function functn(p) - !! use mo_kind, only: dp - !! implicit none - !! real(dp), dimension(:), intent(in) :: p - !! real(dp) :: functn - !! end function functn + !! use mo_optimization_utils, only: optimizee + !! type, extends(optimizee) :: your_optimizee + !! contains + !! ! to be implemented with signature: your_evaluate(self, parameters, sigma, stddev_new, likeli_new) + !! procedure :: evaluate => your_evaluate + !! end type function_optimizee !! \endcode !! !! \b Example !! \code{.f90} !! use mo_opt_functions, only: griewank - !! + !! use mo_optimization_utils, only: function_optimizee + !! type(function_optimizee) :: objective !! prange(:,1) = (/ -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0 /) !! prange(:,2) = (/ 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0 /) !! pini = (/ -.226265E+01, -.130187E+01, -.151219E+01, 0.133983E+00, 0.988159E+00, & !! -.495074E+01, -.126574E+02, 0.572684E+00, 0.303864E+01, 0.343031E+01 /) - !! - !! opt = sce(griewank, pini, prange) + !! objective%func_pointer => griewank + !! opt = sce(objective, pini, prange) !! \endcode !! -> see also example in test directory !! @@ -76,12 +78,11 @@ MODULE mo_sce !! Effective and Efficient Global Optimization for Conceptual Rainfall-runoff Models !! Water Resources Research, Vol 28(4), pp.1015-1031, 1992. !! - !> \param[in] "real(dp) :: functn(p)" Function on which to search the optimum !> \param[in] "real(dp) :: pini(:)" inital value of decision variables !> \param[in] "real(dp) :: prange(size(pini),2)" Min/max range of decision variables ! ! INTENT(INOUT) - ! None + !> \param[inout] "class(optimizee) :: objective" Objective to search the optimum ! INTENT(OUT) ! None @@ -209,7 +210,7 @@ MODULE mo_sce CONTAINS - function sce(eval, functn, pini, prange, & ! IN + function sce(objective, pini, prange, & ! IN mymaxn, mymaxit, mykstop, mypcento, mypeps, myseed, & ! Optional IN myngs, mynpg, mynps, mynspl, mymings, myiniflg, myprint, & ! Optional IN mymask, myalpha, mybeta, & ! Optional IN @@ -229,15 +230,13 @@ CONTAINS !$ use omp_lib, only: OMP_GET_THREAD_NUM, OMP_GET_NUM_THREADS use iso_fortran_env, only : output_unit, error_unit use mo_utils, only : is_finite, special_value - use mo_optimization_utils, only : eval_interface, objective_interface #ifdef MPI use mpi_f08 #endif implicit none ! - procedure(eval_interface), INTENT(IN), POINTER :: eval - procedure(objective_interface), intent(in), pointer :: functn + class(optimizee), INTENT(inout) :: objective real(dp), dimension(:), intent(in) :: pini ! initial parameter set real(dp), dimension(:, :), intent(in) :: prange ! lower and upper bound on parameters ! @@ -702,10 +701,10 @@ CONTAINS end do #endif if (.not. maxit) then - fpini = functn(pini, eval) + fpini = objective%evaluate(pini) history_tmp(1) = fpini else - fpini = -functn(pini, eval) + fpini = -objective%evaluate(pini) history_tmp(1) = -fpini end if @@ -753,9 +752,9 @@ CONTAINS end do #endif if (.not. maxit) then - xf(1) = functn(xx, eval) + xf(1) = objective%evaluate(xx) else - xf(1) = -functn(xx, eval) + xf(1) = -objective%evaluate(xx) end if end if ! @@ -790,10 +789,10 @@ CONTAINS end do #endif if (.not. maxit) then - xf(ii) = functn(xx, eval) + xf(ii) = objective%evaluate(xx) history_tmp(ii) = xf(ii) ! min(history_tmp(ii-1),xf(ii)) --> will be sorted later else - xf(ii) = -functn(xx, eval) + xf(ii) = -objective%evaluate(xx) history_tmp(ii) = -xf(ii) ! max(history_tmp(ii-1),-xf(ii)) --> will be sorted later end if end do @@ -823,11 +822,11 @@ CONTAINS end do #endif if (.not. maxit) then - xf(ii) = functn(xx, eval) + xf(ii) = objective%evaluate(xx) icall = icall + 1_i8 history_tmp(icall) = xf(ii) ! min(history_tmp(icall-1),xf(ii)) --> will be sorted later else - xf(ii) = -functn(xx, eval) + xf(ii) = -objective%evaluate(xx) icall = icall + 1_i8 history_tmp(icall) = -xf(ii) ! max(history_tmp(icall-1),-xf(ii)) --> will be sorted later end if @@ -1105,10 +1104,10 @@ CONTAINS ihistory_tmp = history_tmp #ifdef MPI call cce(s(1 : nps, 1 : nn), sf(1 : nps), bl(1 : nn), bu(1 : nn), maskpara, xnstd(1 : nn), & - iicall, maxn, maxit, save_state_gauss(ithread, :), functn, eval, alpha, beta, ihistory_tmp, idot, comm) + iicall, maxn, maxit, save_state_gauss(ithread, :), objective, alpha, beta, ihistory_tmp, idot, comm) #else call cce(s(1 : nps, 1 : nn), sf(1 : nps), bl(1 : nn), bu(1 : nn), maskpara, xnstd(1 : nn), & - iicall, maxn, maxit, save_state_gauss(ithread, :), functn, eval, alpha, beta, ihistory_tmp, idot) + iicall, maxn, maxit, save_state_gauss(ithread, :), objective, alpha, beta, ihistory_tmp, idot) #endif history_tmp(icall + 1 : icall + (iicall - icall_merk)) = ihistory_tmp(icall_merk + 1 : iicall) icall = icall + (iicall - icall_merk) @@ -1233,10 +1232,10 @@ CONTAINS ! use the sub-complex to generate new point(s) #ifdef MPI call cce(s(1 : nps, 1 : nn), sf(1 : nps), bl(1 : nn), bu(1 : nn), maskpara, xnstd(1 : nn), & - icall, maxn, maxit, save_state_gauss(ithread, :), functn, eval, alpha, beta, history_tmp, idot, comm) + icall, maxn, maxit, save_state_gauss(ithread, :), objective, alpha, beta, history_tmp, idot, comm) #else call cce(s(1 : nps, 1 : nn), sf(1 : nps), bl(1 : nn), bu(1 : nn), maskpara, xnstd(1 : nn), & - icall, maxn, maxit, save_state_gauss(ithread, :), functn, eval, alpha, beta, history_tmp, idot) + icall, maxn, maxit, save_state_gauss(ithread, :), objective, alpha, beta, history_tmp, idot) #endif ! ! if the sub-complex is accepted, replace the new sub-complex @@ -1827,7 +1826,7 @@ CONTAINS end subroutine getpnt - subroutine cce(s, sf, bl, bu, maskpara, xnstd, icall, maxn, maxit, save_state_gauss, functn, eval, & + subroutine cce(s, sf, bl, bu, maskpara, xnstd, icall, maxn, maxit, save_state_gauss, objective, & alpha, beta, history, idot & #ifdef MPI , comm & @@ -1843,7 +1842,6 @@ CONTAINS use mo_xor4096, only : n_save_state use iso_fortran_env, only : output_unit use mo_utils, only : is_finite - use mo_optimization_utils, only : eval_interface, objective_interface #ifdef MPI use mpi_f08 #endif @@ -1862,8 +1860,7 @@ CONTAINS integer(i8), intent(in) :: maxn ! maximal number of function evaluations allowed logical, intent(in) :: maxit ! minimization (true) or maximization (false) integer(i8), dimension(n_save_state), intent(inout) :: save_state_gauss ! seed for random number generation - procedure(objective_interface), intent(in), pointer :: functn - procedure(eval_interface), INTENT(IN), POINTER :: eval + class(optimizee), intent(inout) :: objective real(dp), intent(in) :: alpha ! parameter for reflection steps real(dp), intent(in) :: beta ! parameter for contraction steps real(dp), dimension(maxn), intent(inout) :: history ! history of best function value @@ -1947,11 +1944,11 @@ CONTAINS end do #endif if (.not. maxit) then - fnew = functn(snew, eval) + fnew = objective%evaluate(snew) icall = icall + 1_i8 history(icall) = min(history(icall - 1), fnew) else - fnew = -functn(snew, eval) + fnew = -objective%evaluate(snew) icall = icall + 1_i8 history(icall) = max(history(icall - 1), -fnew) end if @@ -1980,11 +1977,11 @@ CONTAINS end do #endif if (.not. maxit) then - fnew = functn(snew, eval) + fnew = objective%evaluate(snew) icall = icall + 1_i8 history(icall) = min(history(icall - 1), fnew) else - fnew = -functn(snew, eval) + fnew = -objective%evaluate(snew) icall = icall + 1_i8 history(icall) = max(history(icall - 1), -fnew) end if @@ -2011,11 +2008,11 @@ CONTAINS end do #endif if (.not. maxit) then - fnew = functn(snew, eval) + fnew = objective%evaluate(snew) icall = icall + 1_i8 history(icall) = min(history(icall - 1), fnew) else - fnew = -functn(snew, eval) + fnew = -objective%evaluate(snew) icall = icall + 1_i8 history(icall) = max(history(icall - 1), -fnew) end if diff --git a/src/pf_tests/test_mo_anneal.pf b/src/pf_tests/test_mo_anneal.pf index 38407716071630ab2c49806b8b1d1f1da08d1e1a..ba8db78226e41cff4c2a9391f5677f2a99d709dc 100644 --- a/src/pf_tests/test_mo_anneal.pf +++ b/src/pf_tests/test_mo_anneal.pf @@ -4,10 +4,9 @@ module test_mo_anneal use mo_kind, only: dp, i4, i8 use mo_anneal, only: anneal !, anneal_valid use mo_anneal, only: GetTemperature !, GetTemperature_valid - use mo_cost, only: range_dp, cost_objective !, cost_valid_dp + use mo_cost, only: range_dp, cost_dp use mo_xor4096, only: get_timeseed - use mo_opt_functions, only: eval_dummy - use mo_optimization_utils, only: eval_interface, objective_interface + use mo_optimization_utils, only: function_optimizee use mo_message, only: error_message implicit none @@ -30,8 +29,7 @@ module test_mo_anneal real(dp) :: Tstart, Tend real(dp), dimension(:,:), allocatable :: history - procedure(eval_interface), pointer :: eval - procedure(objective_interface), pointer :: obj_func + type(function_optimizee) :: objective ! time dependent seeds call get_timeseed(seeds) @@ -44,16 +42,15 @@ module test_mo_anneal para(3) = 0.5_dp para(4) = 0.4_dp - eval => eval_dummy - obj_func => cost_objective + objective%func_pointer => cost_dp ! Initialization print*, '-----------------------------------' print*, ' INITIALIZATION ' print*, '-----------------------------------' print*, 'Initial parameter set: ',para - print*, 'Initial cost function value: ',cost_objective(para(:), eval) - costbestAll = cost_objective(para(:), eval) + print*, 'Initial cost function value: ',objective%evaluate(para(:)) + costbestAll = objective%evaluate(para(:)) parabestAll = para(:) print*, '-----------------------------------' print*, ' INITIAL TEMPERATURE ' @@ -62,7 +59,7 @@ module test_mo_anneal seeds(1) = 854_i8 seeds(2) = seeds(1) + 1000_i8 print*, 'seeds used: ', seeds(1:2) - temperature = GetTemperature( eval, obj_func, para, 0.95_dp, & + temperature = GetTemperature(objective, para, 0.95_dp, & ! optionals prange_func=range_dp, & samplesize=500_i4, & @@ -82,7 +79,7 @@ module test_mo_anneal print*, 'seeds used: ', seeds(1:3) ! call cpu_time(Tstart) - parabest = anneal(eval, obj_func, para, & + parabest = anneal(objective, para, & ! optionals prange_func=range_dp, & maxit=.false., & @@ -122,9 +119,9 @@ module test_mo_anneal print*, ' parabest = ', parabestAll ! Is program running properly? costbestAll = 3.1142480812726376E-02 - + @assertEqual(costbestAll, 3.1142480812726376E-02_dp, tolerance=t) - + end subroutine test_anneal_dp end module test_mo_anneal diff --git a/src/pf_tests/test_mo_dds.pf b/src/pf_tests/test_mo_dds.pf index e6986dcae3813b3f59d9c458c328aaa27a6221fa..f794c581c6fe9f89af395c1b3bc1c656230cc517 100644 --- a/src/pf_tests/test_mo_dds.pf +++ b/src/pf_tests/test_mo_dds.pf @@ -1,12 +1,12 @@ module test_mo_dds - + use funit use mo_kind, only: i4, i8, dp use mo_dds, only: dds, mdds - use mo_opt_functions, only: eval_dummy, griewank_objective - use mo_optimization_utils, only: eval_interface, objective_interface + use mo_opt_functions, only: griewank + use mo_optimization_utils, only: function_optimizee use mo_message, only: error_message - + implicit none integer(i4), parameter :: np = 10 @@ -24,9 +24,8 @@ module test_mo_dds ! Output of DDS real(dp), dimension(np) :: dv_opt ! Best value of decision variables - procedure(eval_interface), pointer :: eval - procedure(objective_interface), pointer :: obj_func - + type(function_optimizee) :: objective + contains @test @@ -34,7 +33,7 @@ contains real(dp), dimension(10) :: dds_should real(dp) :: t = 1.E-4 - + write(*,*) '--------------------------------------------------' write(*,*) ' Dynamically Dimensioned Search (DDS) Algorithm ' write(*,*) ' version 1.1 - by Bryan Tolson ' @@ -48,10 +47,9 @@ contains nIterMax = 100000_i8 to_max = .False. seed = 123456789_i8 - eval => eval_dummy - obj_func => griewank_objective - dv_opt = DDS(eval, obj_func, dv_ini, dv_range, r=r_val, seed=seed, maxiter=nIterMax, maxit=to_max) + objective%func_pointer => griewank + dv_opt = DDS(objective, dv_ini, dv_range, r=r_val, seed=seed, maxiter=nIterMax, maxit=to_max) write(*,*) '' write(*,*) 'DDS' @@ -61,12 +59,12 @@ contains dds_should = (/ 3.0825, 0.0629, 0.0362, -0.0063, 0.0026, -7.6358, 0.0718, -0.0528, -0.1008, 0.0037 /) @assertEqual(dv_opt, dds_should, tolerance = t) - + end subroutine test_dds @test subroutine test_mdds() - + real(dp), dimension(10) :: mdds_should real(dp) :: t = 1.E-4 @@ -78,7 +76,8 @@ contains nIterMax = 100000_i8 to_max = .False. seed = 123456789_i8 - dv_opt = MDDS(eval, obj_func, dv_ini, dv_range, seed=seed, maxiter=nIterMax, maxit=to_max) + objective%func_pointer => griewank + dv_opt = MDDS(objective, dv_ini, dv_range, seed=seed, maxiter=nIterMax, maxit=to_max) write(*,*) '' write(*,*) 'MDDS' @@ -88,7 +87,7 @@ contains mdds_should = (/ 3.1209, 0.0577, -0.0168, -6.3592, -0.0393, 0.0191, -8.2359, 0.0158, -9.4605, -0.0736 /) @assertEqual(dv_opt, mdds_should, tolerance = t) - + end subroutine test_mdds - -end module test_mo_dds \ No newline at end of file + +end module test_mo_dds diff --git a/src/pf_tests/test_mo_mcmc.pf b/src/pf_tests/test_mo_mcmc.pf index 238cad7e70020f03417925e660d7875f20f15260..a723f11efe1a172ba55201237cde7883a61fbb89 100644 --- a/src/pf_tests/test_mo_mcmc.pf +++ b/src/pf_tests/test_mo_mcmc.pf @@ -5,7 +5,7 @@ module test_mo_mcmc use mo_likelihood, only: setmeas, loglikelihood_dp, loglikelihood_stddev_dp, model_dp use mo_mcmc, only: mcmc, mcmc_stddev use mo_moment, only: mean, stddev - use mo_optimization_utils, only: eval_interface, objective_interface + use mo_optimization_utils, only: eval_interface, objective_interface, eval_optimizee use mo_message, only: error_message implicit none @@ -13,6 +13,7 @@ module test_mo_mcmc ! for running mcmc procedure(eval_interface), pointer :: eval_func procedure(objective_interface), pointer :: likelihood + type(eval_optimizee) :: objective real(dp) :: p real(dp) :: likeli_new @@ -46,8 +47,9 @@ contains write(*,*) ' (A1) "real" likelihood (sigma is an error model_dp or given) --> e.g. loglikelihood of mo_likelihood' write(*,*) ' full run ' write(*,*) '---------------------------------------------------------------------------------------------' - eval_func => model_dp - p = loglikelihood_dp(parabest, eval_func) + objective%eval_pointer => model_dp + objective%obj_pointer => loglikelihood_dp + p = objective%evaluate(parabest) write(*,*) 'log-likelihood = ',p ! initializing the ranges of parameters @@ -64,8 +66,7 @@ contains ! (1) Burn-in will be performed to optimize settings for MCMC ! (2) posterior distribution of the parameters at the minimum (best parameterset) ! will be sampled by MCMC - likelihood => loglikelihood_dp - call mcmc(eval_func, likelihood, parabest, rangePar, mcmc_paras, burnin_paras, & + call mcmc(objective, parabest, rangePar, mcmc_paras, burnin_paras, & ParaSelectMode_in=2_i4,tmp_file='A_make_check_test_file', & restart=.false., restart_file='restart_make_check_test_file', & seed_in=seed, loglike_in=.true., maskpara_in=maskpara, printflag_in=.true.) @@ -114,7 +115,9 @@ contains write(*,*) '---------------------------------------------------------------------------------------------' ! starting MCMC: ! (1) starting from restart file - call mcmc(eval_func, likelihood, parabest, rangePar, mcmc_paras, burnin_paras, & + objective%eval_pointer => model_dp + objective%obj_pointer => loglikelihood_dp + call mcmc(objective, parabest, rangePar, mcmc_paras, burnin_paras, & ParaSelectMode_in=2_i4,tmp_file='A_make_check_test_file', & restart=.true., restart_file='restart_make_check_test_file', & seed_in=seed, loglike_in=.true., maskpara_in=maskpara, printflag_in=.true.) @@ -165,7 +168,9 @@ contains ! sigma of errors is unknown --> initial guess e.g. 1.0_dp ! stddev_new = standard deviation of errors using paraset ! likeli_new = likelihood using stddev_new - p = loglikelihood_stddev_dp(parabest,eval_func,1.0_dp,likeli_new=likeli_new,stddev_new=stddev_new) + objective%eval_pointer => model_dp + objective%obj_pointer => loglikelihood_stddev_dp + p = objective%evaluate(parabest,1.0_dp,likeli_new=likeli_new,stddev_new=stddev_new) write(*,*) 'guessed log-likelihood = ',p write(*,*) 'best log-likelihood = ',likeli_new @@ -183,8 +188,7 @@ contains ! (1) Burn-in will be performed to optimize settings for MCMC ! (2) posterior distribution of the parameters at the minimum (best parameterset) ! will be sampled by MCMC - likelihood => loglikelihood_stddev_dp - call mcmc_stddev(eval_func, likelihood, parabest, rangePar, mcmc_paras, burnin_paras, & + call mcmc_stddev(objective, parabest, rangePar, mcmc_paras, burnin_paras, & ParaSelectMode_in=2_i4,tmp_file='B_make_check_test_file', & seed_in=seed, loglike_in=.true., maskpara_in=maskpara, printflag_in=.true.) @@ -221,5 +225,5 @@ contains deallocate(burnin_paras) end subroutine test_mcmc_stddev - + end module test_mo_mcmc \ No newline at end of file diff --git a/src/pf_tests/test_mo_sce.pf b/src/pf_tests/test_mo_sce.pf index fa8d98b6970c6ec5ae4b86db5f8655a0e09fee66..9a12259e47512ddfd632f2ba36c4611f836196d3 100644 --- a/src/pf_tests/test_mo_sce.pf +++ b/src/pf_tests/test_mo_sce.pf @@ -1,12 +1,12 @@ module test_mo_sce - + use funit use mo_kind, only: i4, i8, dp use mo_sce, only: sce - use mo_opt_functions, only: eval_dummy, ackley_objective, griewank_objective - use mo_optimization_utils, only: eval_interface, objective_interface + use mo_opt_functions, only: ackley, griewank + use mo_optimization_utils, only: function_optimizee use mo_message, only: error_message - + implicit none real(dp), dimension(:), allocatable :: pini @@ -21,9 +21,8 @@ module test_mo_sce integer(i4) :: n logical :: parallel - procedure(eval_interface), pointer :: eval - procedure(objective_interface), pointer :: obj_func - + type(function_optimizee) :: objective + contains @test @@ -31,7 +30,7 @@ contains real(dp) :: t = 1.e-7_dp integer(i4) :: iPar - + ! Dimension of test function = number of parameters n = 30 @@ -51,10 +50,9 @@ contains mask(:) = .true. pini(:) = 0.5_dp parallel=.false. - eval => eval_dummy - obj_func => ackley_objective - opt = sce(eval, & ! Mandatory IN: eval function - obj_func, & ! Mandatory IN: Objective function + objective%func_pointer => ackley + opt = sce( & + objective, & ! Mandatory INOUT: objective pini, & ! Mandatory IN: initial guess prange, & ! Mandatory IN: range for each parameter (min, max) mymaxn=30000_i8, & ! Optional IN: maximal number of function evaluations @@ -97,7 +95,7 @@ contains write(*,*) '' write(*,*) 'number of function evaluations: neval = ', neval - write(*,*) 'best function value found: bestf = ', obj_func(opt, eval) + write(*,*) 'best function value found: bestf = ', objective%evaluate(opt) write(*,*) 'global minimal function value: optf = ', 0.0_dp write(*,*) '' @@ -116,13 +114,13 @@ contains end if ! Check restart - + bestf = 0. neval = 0 history = 0. - opt = sce(eval, & ! Mandatory IN : Objective function - obj_func, & ! Mandatory IN: Objective function + opt = sce( & + objective, & ! Mandatory INOUT: objective pini, & ! Mandatory IN : initial guess prange, & ! Mandatory IN : range for each parameter (min, max) restart=.true., & ! Do the restart @@ -134,7 +132,7 @@ contains write(*,*) '' write(*,*) 'number of function evaluations: neval = ', neval - write(*,*) 'best function value found: bestf = ', obj_func(opt, eval) + write(*,*) 'best function value found: bestf = ', objective%evaluate(opt) write(*,*) 'global minimal function value: optf = ', 0.0_dp write(*,*) '' @@ -159,7 +157,7 @@ contains real(dp) :: t = 1.e-7_dp integer(i4) :: iPar - + ! Dimension of test function = number of parameters n = 3 @@ -179,10 +177,9 @@ contains mask(:) = .true. pini(:) = 0.5_dp parallel=.false. - eval => eval_dummy - obj_func => griewank_objective - opt = sce(eval, & ! Mandatory IN: eval function - obj_func, & ! Mandatory IN: Objective function + objective%func_pointer => griewank + opt = sce( & + objective, & ! Mandatory INOUT: objective pini, & ! Mandatory IN: initial guess prange, & ! Mandatory IN: range for each parameter (min, max) mymaxn=30000_i8, & ! Optional IN: maximal number of function evaluations @@ -225,7 +222,7 @@ contains write(*,*) '' write(*,*) 'number of function evaluations: neval = ', neval - write(*,*) 'best function value found: bestf = ', obj_func(opt, eval) + write(*,*) 'best function value found: bestf = ', objective%evaluate(opt) write(*,*) 'global minimal function value: optf = ', 0.0_dp write(*,*) '' @@ -244,13 +241,13 @@ contains end if ! Check restart - + bestf = 0. neval = 0 history = 0. - opt = sce(eval, & ! Mandatory IN : Objective function - obj_func, & ! Mandatory IN: Objective function + opt = sce( & + objective, & ! Mandatory INOUT: objective pini, & ! Mandatory IN : initial guess prange, & ! Mandatory IN : range for each parameter (min, max) restart=.true., & ! Do the restart @@ -262,7 +259,7 @@ contains write(*,*) '' write(*,*) 'number of function evaluations: neval = ', neval - write(*,*) 'best function value found: bestf = ', obj_func(opt, eval) + write(*,*) 'best function value found: bestf = ', objective%evaluate(opt) write(*,*) 'global minimal function value: optf = ', 0.0_dp write(*,*) '' @@ -281,5 +278,5 @@ contains deallocate(pini, prange, opt, mask) end subroutine test_sce_griewank - -end module test_mo_sce \ No newline at end of file + +end module test_mo_sce