LCOV - code coverage report
Current view: top level - common - mo_optimization.F90 (source / functions) Hit Total Coverage
Test: mHM coverage Lines: 35 82 42.7 %
Date: 2024-04-30 08:53:32 Functions: 1 1 100.0 %

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

Generated by: LCOV version 1.16