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