The ParaMonte Documentation Website
Current view: top level - kernel - ParaMCMCRefinedChain_mod.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: MPI Parallel Kernel - Code Coverage Report Lines: 223 224 99.6 %
Date: 2021-01-08 13:07:16 Functions: 6 6 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       3             : !!!!
       4             : !!!!   MIT License
       5             : !!!!
       6             : !!!!   ParaMonte: plain powerful parallel Monte Carlo library.
       7             : !!!!
       8             : !!!!   Copyright (C) 2012-present, The Computational Data Science Lab
       9             : !!!!
      10             : !!!!   This file is part of the ParaMonte library.
      11             : !!!!
      12             : !!!!   Permission is hereby granted, free of charge, to any person obtaining a
      13             : !!!!   copy of this software and associated documentation files (the "Software"),
      14             : !!!!   to deal in the Software without restriction, including without limitation
      15             : !!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
      16             : !!!!   and/or sell copies of the Software, and to permit persons to whom the
      17             : !!!!   Software is furnished to do so, subject to the following conditions:
      18             : !!!!
      19             : !!!!   The above copyright notice and this permission notice shall be
      20             : !!!!   included in all copies or substantial portions of the Software.
      21             : !!!!
      22             : !!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
      23             : !!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
      24             : !!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
      25             : !!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
      26             : !!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
      27             : !!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
      28             : !!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
      29             : !!!!
      30             : !!!!   ACKNOWLEDGMENT
      31             : !!!!
      32             : !!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
      33             : !!!!   As per the ParaMonte library license agreement terms, if you use any parts of
      34             : !!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
      35             : !!!!   work (education/research/industry/development/...) by citing the ParaMonte
      36             : !!!!   library as described on this page:
      37             : !!!!
      38             : !!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
      39             : !!!!
      40             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      41             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      42             : 
      43             : !>  \brief This module contains the classes and procedures for refining MCMC output chains.
      44             : !>  \author Amir Shahmoradi
      45             : 
      46             : module ParaMCMCRefinedChain_mod
      47             : 
      48             :     use SpecMCMC_SampleRefinementMethod_mod, only: BATCH_MEANS_METHOD_NAME
      49             :     use SpecMCMC_SampleRefinementMethod_mod, only: CUTOFF_AUTOCORR_METHOD_NAME
      50             :     use SpecMCMC_SampleRefinementMethod_mod, only: MAX_CUMSUM_AUTOCORR_METHOD_NAME
      51             :     use ParaMonteChainFileContents_mod, only: Count_type
      52             :     use JaggedArray_mod, only: CharVec_type
      53             :     use Constants_mod, only: IK, RK
      54             :     use Err_mod, only: Err_type
      55             : 
      56             :     implicit none
      57             : 
      58             :     character(*), parameter :: MODULE_NAME = "\paramCMCRefinedChain_mod"
      59             : 
      60             :     !> The `RefinedChain_type` class.
      61             :     type                                    :: RefinedChain_type
      62             :         integer(IK)                         :: ndim  = 0_IK         ! number of sampling variables
      63             :         integer(IK)                         :: numRefinement = 0_IK ! number of refinements, zero if sample size is prescribed by the user
      64             :         type(Count_type)    , allocatable   :: Count(:)             ! compact and verbose counts
      65             :         real(RK)            , allocatable   :: IAC(:,:)             ! size of (ndim,0:numRefinement): The Integrated AutoCorrelation Time at each refinement stage
      66             :         real(RK)            , allocatable   :: LogFuncState(:,:)    ! size of (ndim,Count%compact): LogFuncState is LogFunc + Variables
      67             :         integer(IK)         , allocatable   :: Weight(:)            ! size of (Count%compact): Weight of each state
      68             :         type(CharVec_type)  , allocatable   :: ColHeader(:)         ! refined sample column headers
      69             :         type(Err_type)                      :: Err
      70             :     contains
      71             :         procedure, pass :: get => getRefinedChain
      72             :         procedure, pass :: write => writeRefinedChain
      73             :     end type RefinedChain_type
      74             : 
      75             :    type Method_type
      76             :        logical :: isMax = .false.
      77             :        logical :: isMed = .false.
      78             :        logical :: isMin = .false.
      79             :        logical :: isAvg = .false.
      80             :        logical :: isBatchMeans = .false.
      81             :        logical :: isCutoffAutoCorr = .false.
      82             :        logical :: isViaCompactChain = .true.
      83             :        logical :: isViaVerboseChain = .true.
      84             :        logical :: isMaxCumSumAutoCorr = .false.
      85             :    end type Method_type
      86             : 
      87             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      88             : 
      89             : contains
      90             : 
      91             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      92             : 
      93             :     !> Return the refined Markov chain, given the input Markov chain and its specifications.
      94             :     !> This procedure is a method of the [RefinedChain_type](@ref refinedchain_type) class.
      95             :     !>
      96             :     !> \param[inout] RefinedChain           :   An object of class [RefinedChain_type](@ref refinedchain_type).
      97             :     !> \param[inout] CFC                    :   An object of type [ChainFileContents_type](@ref paramontechainfilecontents_mod::chainfilecontents_type)
      98             :     !!                                          containing the Markov chain.
      99             :     !> \param[out]   Err                    :   An object of class [Err_type](@ref err_mod::err_type) indicating whether any error has occurred or not.
     100             :     !> \param[in]    burninLoc              :   The estimated location of burnin point in the Markov chain (**optional**).
     101             :     !!                                          If not provided, it will be extracted from the components of the input `CFC`.
     102             :     !> \param[in]    refinedChainSize       :   The requested refined sample size (**optional**). If the size of the refined sample is given as input,
     103             :     !!                                          then the requested sample is directly generated based on the input size.
     104             :     !> \param[in]    sampleRefinementCount  :   The maximum number of times the sample can be refined (**optional**, default = `Infinity`).
     105             :     !!                                      :   For example, if set to 1, then only one round of refinement will be performed on the Markov chain.
     106             :     !> \param[in]    sampleRefinementMethod :   The requested method of refining the sample (**optional**, default = "BatchMeans").
     107         585 :     subroutine getRefinedChain  ( RefinedChain              &
     108             :                                 , CFC                       &
     109             :                                 , Err                       &
     110             :                                 , burninLoc                 &
     111             :                                 , refinedChainSize          &
     112             :                                 , sampleRefinementCount     &
     113             :                                 , sampleRefinementMethod    &
     114             :                                 )
     115             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     116             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRefinedChain
     117             : #endif
     118             : 
     119             :         use, intrinsic :: iso_fortran_env, only: output_unit
     120             :         use ParaMonteChainFileContents_mod, only: ChainFileContents_type, NUM_DEF_COL
     121             :         use CrossCorr_mod, only: getBatchMeansIAC, getCumSumIAC, getMaxCumSumIAC
     122             :         use String_mod, only: getLowerCase, replaceStr
     123             :         use Sort_mod, only: getMedian
     124             : 
     125             :         implicit none
     126             : 
     127             :         character(*), parameter :: PROCEDURE_NAME = MODULE_NAME//"@getRefinedChain()"
     128             : 
     129             :         class(RefinedChain_type)    , intent(inout)             :: RefinedChain
     130             :         type(ChainFileContents_type), intent(inout)             :: CFC
     131             :         type(Err_type)              , intent(out)               :: Err
     132             :         character(*)                , intent(in), optional      :: sampleRefinementMethod
     133             :         integer(IK)                 , intent(in), optional      :: burninLoc, refinedChainSize, sampleRefinementCount
     134             :         integer(IK)                                             :: burninLocDefault, i, countCompact, ndimPlusOne
     135             :         integer(IK)                                             :: maxRefinementCurrent, maxRefinementCount
     136         354 :         real(RK)                                                :: integratedAutoCorrTime, ndimPlusOneInverse
     137             :         logical                                                 :: refinedChainSizeIsPresent
     138             :         logical                                                 :: maxRefinementCountIsReached
     139         354 :         character(:)    , allocatable                           :: sampleRefinementMethodDefault, sampleRefinementMethodDefaultLowerCase, srmethod
     140             :         type(Count_type), allocatable                           :: DumCount(:)
     141             :         real(RK)        , allocatable                           :: DumIAC(:,:), DumVec(:)
     142             :         integer(IK)     , allocatable                           :: SampleWeight(:)
     143             :         type(Method_type)                                       :: Method
     144         354 :         type(ChainFileContents_type)                            :: DumCFC
     145             : 
     146         354 :         Err%occurred = .false.
     147             : 
     148             :         ! if the size of the refined sample is given as input, then generate the requested sample straight
     149             : 
     150         354 :         refinedChainSizeIsPresent = present(refinedChainSize)
     151         354 :         if (refinedChainSizeIsPresent) then    ! ignore sampleRefinementCount, even if it is given by the user
     152          85 :             maxRefinementCount = 1_IK
     153             :         else
     154         269 :             maxRefinementCount = 20_IK  ! this is a temporary maximum value, to be increased later if needed
     155         269 :             if (present(sampleRefinementCount)) maxRefinementCount = sampleRefinementCount
     156             :         end if
     157             : 
     158             :         ! this is to avoid memory overflow due to extremely large maxRefinementCount requested by the user
     159             : 
     160         354 :         maxRefinementCurrent = min(2_IK, maxRefinementCount)
     161             : 
     162             :         ! compute ndim and the initial chain size
     163             : 
     164         354 :         RefinedChain%numRefinement = 0_IK
     165             : 
     166         354 :         if (CFC%ndim == 0_IK) then
     167           6 :             RefinedChain%ndim = size(CFC%State(:,1))
     168           6 :             CFC%ndim = RefinedChain%ndim
     169             :         else
     170         348 :             RefinedChain%ndim = CFC%ndim
     171             :         end if
     172             : 
     173             :         ! allocate components
     174             : 
     175         354 :         if (allocated(RefinedChain%IAC)) deallocate(RefinedChain%IAC)
     176         354 :         if (allocated(RefinedChain%Count)) deallocate(RefinedChain%Count)
     177         354 :         allocate(RefinedChain%IAC(0:RefinedChain%ndim,0:maxRefinementCurrent))
     178        1327 :         allocate(RefinedChain%Count(0:maxRefinementCurrent))
     179         354 :         if (CFC%Count%compact == 0_IK) then
     180           6 :             RefinedChain%Count(RefinedChain%numRefinement)%compact = size(CFC%State(1,:))
     181           6 :             CFC%Count%compact = RefinedChain%Count(RefinedChain%numRefinement)%compact
     182             :         else
     183         348 :             RefinedChain%Count(RefinedChain%numRefinement)%compact = CFC%Count%compact
     184             :         end if
     185         354 :         if (CFC%Count%verbose == 0_IK) CFC%Count%verbose = sum(CFC%Weight(1:CFC%Count%compact))
     186         354 :         RefinedChain%Count(RefinedChain%numRefinement)%verbose = CFC%Count%verbose
     187             : 
     188             :         ! define the AIC computation method
     189             : 
     190         354 :         sampleRefinementMethodDefault = BATCH_MEANS_METHOD_NAME; if (present(sampleRefinementMethod)) sampleRefinementMethodDefault = sampleRefinementMethod
     191         354 :         sampleRefinementMethodDefaultLowerCase = getLowerCase(sampleRefinementMethodDefault)
     192             : 
     193         354 :         if (index(sampleRefinementMethodDefaultLowerCase,getLowerCase(BATCH_MEANS_METHOD_NAME))>0) then
     194         294 :             Method%isBatchMeans = .true.
     195         120 :         elseif (index(sampleRefinementMethodDefaultLowerCase,getLowerCase(CUTOFF_AUTOCORR_METHOD_NAME))>0 .or. index(sampleRefinementMethodDefaultLowerCase,"cutoff")>0) then
     196          28 :             Method%isCutoffAutoCorr = .true.
     197          64 :         elseif (index(sampleRefinementMethodDefaultLowerCase,getLowerCase(MAX_CUMSUM_AUTOCORR_METHOD_NAME))>0 .or. index(sampleRefinementMethodDefaultLowerCase,"cumsum")>0) then
     198          32 :             Method%isMaxCumSumAutoCorr = .true.
     199             :         else
     200             :             ! LCOV_EXCL_START
     201             :             Err%occurred = .true.
     202             :             Err%msg = PROCEDURE_NAME // ": Unknown unsupported IAC computation method name: " // sampleRefinementMethodDefault
     203             :             return
     204             :             ! LCOV_EXCL_STOP
     205             :         end if
     206             : 
     207             :         ! define the chain types to use for the AIC computation
     208             : 
     209         354 :         Method%isViaCompactChain = index(sampleRefinementMethodDefaultLowerCase,"compact") > 0
     210         354 :         Method%isViaVerboseChain = index(sampleRefinementMethodDefaultLowerCase,"verbose") > 0
     211         354 :         if (.not.(Method%isViaCompactChain .or. Method%isViaVerboseChain)) then
     212         354 :             Method%isViaCompactChain = .true.
     213         354 :             Method%isViaVerboseChain = .true.
     214             :         end if
     215             : 
     216             :         ! define the statistic to use for the AIC computation
     217             : 
     218         354 :         srmethod = replaceStr(string=sampleRefinementMethodDefaultLowerCase,search=" ",substitute="-")
     219         354 :         Method%isAvg =  index(srmethod,"-avg") > 0 .or. index(srmethod,"-average") > 0
     220         354 :         Method%isMed =  index(srmethod,"-med") > 0 .or. index(srmethod,"-median") > 0
     221         354 :         Method%isMin =  index(srmethod,"-min") > 0 .or. index(srmethod,"-minimum") > 0
     222         354 :         Method%isMax =  index(srmethod,"-max") > 0 .or. index(srmethod,"-maximum") > 0
     223         354 :         if ( Method%isAvg .and. (Method%isMed .or. Method%isMax .or. Method%isMin) ) Err%occurred = .true.
     224         354 :         if ( Method%isMed .and. (Method%isMax .or. Method%isMin) ) Err%occurred = .true.
     225         354 :         if ( Method%isMax .and. Method%isMin ) Err%occurred = .true.
     226         354 :         if (.not.(Method%isAvg .or. Method%isMed .or. Method%isMin .or. Method%isMax)) Method%isAvg = .true. ! default method of AIC summarization.
     227             : 
     228         354 :         if (Method%isAvg) ndimPlusOneInverse = 1._RK / (RefinedChain%ndim + 1_IK)
     229         354 :         if (Method%isMed) ndimPlusOne = RefinedChain%ndim + 1_IK
     230             : 
     231             :         ! assign the column headers
     232             : 
     233         354 :         if (allocated(CFC%ColHeader)) then
     234         768 :             if (allocated(RefinedChain%ColHeader)) deallocate(RefinedChain%ColHeader)
     235        1306 :             allocate(RefinedChain%ColHeader(0:RefinedChain%ndim))
     236        1306 :             do i = 0, RefinedChain%ndim
     237        1306 :                 RefinedChain%ColHeader(i)%record = CFC%ColHeader(i+NUM_DEF_COL)%record
     238             :             end do
     239             :         else
     240             :             ! LCOV_EXCL_START
     241             :             Err%occurred = .true.
     242             :             Err%msg = PROCEDURE_NAME // ": Internal error occurred. CFC%ColHeader is NULL."
     243             :             return
     244             :             ! LCOV_EXCL_STOP
     245             :         end if
     246             : 
     247         354 :         if (present(burninLoc)) then
     248         348 :             burninLocDefault = burninLoc
     249             :         else
     250           6 :             burninLocDefault = CFC%BurninLoc(CFC%Count%compact)
     251             :         end if
     252             : 
     253         354 :         RefinedChain%Count(RefinedChain%numRefinement)%compact = CFC%Count%compact - burninLocDefault + 1
     254         354 :         if (allocated(RefinedChain%LogFuncState)) deallocate(RefinedChain%LogFuncState)
     255         354 :         allocate(RefinedChain%LogFuncState(0:RefinedChain%ndim,RefinedChain%Count(RefinedChain%numRefinement)%compact))
     256      139076 :         RefinedChain%LogFuncState(0                  ,1:RefinedChain%Count(RefinedChain%numRefinement)%compact) = CFC%LogFunc(burninLocDefault:CFC%Count%compact)
     257      384902 :         RefinedChain%LogFuncState(1:RefinedChain%ndim,1:RefinedChain%Count(RefinedChain%numRefinement)%compact) = CFC%State(1:RefinedChain%ndim, burninLocDefault:CFC%Count%compact)
     258             : 
     259             :         ! check if there are more than 1 sample points in the burnin-subtracted CFC
     260             : 
     261         354 :         if (RefinedChain%Count(RefinedChain%numRefinement)%compact==0_IK) then
     262             :             ! LCOV_EXCL_START
     263             :             Err%occurred = .true.
     264             :             Err%msg = PROCEDURE_NAME // ": The size of the refined sample is zero."
     265             :             return
     266             :             ! LCOV_EXCL_STOP
     267         354 :         elseif (RefinedChain%Count(RefinedChain%numRefinement)%compact==1_IK) then
     268           6 :             if (allocated(RefinedChain%Weight)) deallocate(RefinedChain%Weight)
     269           6 :             allocate(RefinedChain%Weight(RefinedChain%Count(RefinedChain%numRefinement)%compact))
     270           6 :             if (refinedChainSizeIsPresent) then
     271           0 :                 RefinedChain%Weight = refinedChainSize
     272             :             else
     273          12 :                 RefinedChain%Weight = 1
     274             :             end if
     275           6 :             if (allocated(RefinedChain%IAC)) deallocate(RefinedChain%IAC)
     276           6 :             if (allocated(RefinedChain%Count)) deallocate(RefinedChain%Count)
     277           6 :             allocate(RefinedChain%IAC(0:RefinedChain%ndim,0:0))
     278          12 :             allocate(RefinedChain%Count(0:0))
     279          30 :             RefinedChain%IAC = 0._RK
     280          12 :             RefinedChain%Count(RefinedChain%numRefinement)%verbose = sum(RefinedChain%Weight)
     281           6 :             return
     282             :         end if
     283             : 
     284      139064 :         RefinedChain%Weight = CFC%Weight(burninLocDefault:CFC%Count%compact) ! Weight is intentionally separately assigned from State here
     285             : 
     286             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     287             : 
     288         491 :         loopRefinement: do
     289             : 
     290        1187 :             countCompact = RefinedChain%Count(RefinedChain%numRefinement)%compact
     291             : 
     292             :             ! set the sample weight
     293             : 
     294        1187 :             if (Method%isViaCompactChain) then
     295         727 :                 if (allocated(SampleWeight)) deallocate(SampleWeight)
     296         460 :             elseif (Method%isViaVerboseChain) then
     297       89936 :                 SampleWeight = RefinedChain%Weight
     298             :             endif
     299             : 
     300             :             ! obtain the IAC for each individual variable
     301             : 
     302        4375 :             do i = 0, RefinedChain%ndim
     303      792268 :                 DumVec = RefinedChain%LogFuncState(i,1:RefinedChain%Count(RefinedChain%numRefinement)%compact)
     304        4375 :                 if (Method%isBatchMeans) then
     305       10568 :                     RefinedChain%IAC(i,RefinedChain%numRefinement) = getBatchMeansIAC   ( np        = countCompact  &
     306             :                                                                                         , Point     = DumVec        & ! LCOV_EXCL_LINE
     307             :                                                                                         , Weight    = SampleWeight  & ! LCOV_EXCL_LINE
     308        2642 :                                                                                         )
     309         546 :                 elseif (Method%isCutoffAutoCorr) then
     310         996 :                     RefinedChain%IAC(i,RefinedChain%numRefinement) = getCumSumIAC       ( np        = countCompact  &
     311             :                                                                                         , Point     = DumVec        & ! LCOV_EXCL_LINE
     312             :                                                                                         , Weight    = SampleWeight  & ! LCOV_EXCL_LINE
     313         249 :                                                                                         )
     314         297 :                 elseif (Method%isMaxCumSumAutoCorr) then
     315        1188 :                     RefinedChain%IAC(i,RefinedChain%numRefinement) = getMaxCumSumIAC    ( np        = countCompact  &
     316             :                                                                                         , Point     = DumVec        & ! LCOV_EXCL_LINE
     317             :                                                                                         , Weight    = SampleWeight  & ! LCOV_EXCL_LINE
     318         297 :                                                                                         )
     319             :                 end if
     320             :             end do
     321             : 
     322        1187 :             if (refinedChainSizeIsPresent) then
     323       83024 :                 integratedAutoCorrTime = real(sum(RefinedChain%Weight),kind=RK) / real(refinedChainSize,kind=RK)
     324             :             else
     325         957 :                 if (Method%isAvg) then
     326        2762 :                     integratedAutoCorrTime = sum( RefinedChain%IAC(0:RefinedChain%ndim,RefinedChain%numRefinement) ) * ndimPlusOneInverse
     327         176 :                 elseif (Method%isMed) then
     328          88 :                     call getMedian(lenArray=ndimPlusOne,Array=RefinedChain%IAC(0:RefinedChain%ndim,RefinedChain%numRefinement),median=integratedAutoCorrTime,Err=Err)
     329          88 :                     if (Err%occurred) then
     330             :                         ! LCOV_EXCL_START
     331             :                         Err%msg = PROCEDURE_NAME//Err%msg
     332             :                         return
     333             :                         ! LCOV_EXCL_STOP
     334             :                     end if
     335          88 :                 elseif (Method%isMax) then
     336         305 :                     integratedAutoCorrTime = maxval( RefinedChain%IAC(0:RefinedChain%ndim,RefinedChain%numRefinement) )
     337          27 :                 elseif (Method%isMin) then
     338         135 :                     integratedAutoCorrTime = minval( RefinedChain%IAC(0:RefinedChain%ndim,RefinedChain%numRefinement) )
     339             :                 end if
     340             :             end if
     341             : 
     342             :             ! so far, we have computed the max IAC of the sample, but no refinement. Refine the sample only if needed.
     343             : 
     344        1187 :             maxRefinementCountIsReached = RefinedChain%numRefinement==maxRefinementCount
     345        1187 :             if (integratedAutoCorrTime<2._RK .or. maxRefinementCountIsReached) then
     346         696 :                 if (Method%isViaCompactChain .and. Method%isViaVerboseChain) then
     347             :                     !if (maxRefinementCountIsReached) maxRefinementCount = maxRefinementCount * 2_IK
     348         348 :                     maxRefinementCount = maxRefinementCount * 2_IK
     349         348 :                     Method%isViaCompactChain = .false.
     350         348 :                     cycle loopRefinement
     351             :                 end if
     352         348 :                 exit loopRefinement
     353             :             end if
     354             : 
     355             :             ! generate the refined sample, dump it in CFC, then put it back into RefinedChain to start over again
     356             : 
     357         491 :             RefinedChain%numRefinement = RefinedChain%numRefinement + 1_IK
     358             : 
     359             :             ! reallocate to bigger array if nedded
     360             : 
     361         491 :             if (RefinedChain%numRefinement>maxRefinementCurrent) then
     362          23 :                 call move_alloc( from = RefinedChain%IAC, to = DumIAC )
     363          23 :                 call move_alloc( from = RefinedChain%Count, to = DumCount )
     364          23 :                 maxRefinementCurrent = min( maxRefinementCurrent*2 , maxRefinementCount )
     365          23 :                 allocate( RefinedChain%IAC( 0:RefinedChain%ndim , 0:maxRefinementCurrent ) )
     366         146 :                 allocate( RefinedChain%Count( 0:maxRefinementCurrent ) )
     367         297 :                 RefinedChain%IAC(0:RefinedChain%ndim,0:RefinedChain%numRefinement-1) = DumIAC
     368          96 :                 RefinedChain%Count(0:RefinedChain%numRefinement-1) = DumCount
     369             :             end if
     370             : 
     371         491 :             if (integratedAutoCorrTime<2._RK) cycle loopRefinement ! no need for refinement. should happen only when transitioning from compact to verbose
     372             : 
     373             :             call refineWeightedSample   ( nd = RefinedChain%ndim                                        & ! LCOV_EXCL_LINE
     374             :                                         , np = countCompact                                             & ! LCOV_EXCL_LINE
     375             :                                         , skip = integratedAutoCorrTime                                 & ! LCOV_EXCL_LINE
     376             :                                         , Sample = RefinedChain%LogFuncState                            & ! LCOV_EXCL_LINE
     377             :                                         , Weight = RefinedChain%Weight                                  & ! LCOV_EXCL_LINE
     378             :                                         , RefinedChain = DumCFC%State                                   & ! LCOV_EXCL_LINE
     379             :                                         , RefinedWeight = DumCFC%Weight                                 & ! LCOV_EXCL_LINE
     380         982 :                                         , PointCount = RefinedChain%Count(RefinedChain%numRefinement)   &
     381             :                                         , refinedChainSize = refinedChainSize                           & ! LCOV_EXCL_LINE
     382         491 :                                         )
     383       71279 :             RefinedChain%Weight        = DumCFC%Weight
     384      266069 :             RefinedChain%LogFuncState  = DumCFC%State
     385             : 
     386             :         end do loopRefinement
     387             : 
     388             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     389             : 
     390         354 :     end subroutine getRefinedChain
     391             : 
     392             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     393             : 
     394             :     !> Return the refined vector of weights of the vector of weights of a weighted Markov chain.
     395             :     !>
     396             :     !> \param[in]   np                  :   The number of elements of the `Weight` vector.
     397             :     !> \param[in]   Weight              :   The input vector of weights.
     398             :     !> \param[in]   skip                :   The size of the jumps that have to be made through the weighted Markov chain.
     399             :     !> \param[in]   refinedChainSize    :   The requested refined sample size (**optional**). If present, then the refined chain (represented by the
     400             :     !>                                  :   vector `Weight`) will be refined such that the resulting refined chain has the size `refinedChainSize`.
     401             :     !>
     402             :     !> \return
     403             :     !> `RefinedWeight` : An array of size `np`, whose elements indicate which points are present in the final refined chain.\n
     404             :     !> Examples:
     405             :     !> ```
     406             :     !>          skip: 1
     407             :     !>        Weight: 5, 0, 1, 3, 1
     408             :     !> RefinedWeight: 5, 0, 1, 3, 1
     409             :     !>          skip: 2
     410             :     !>        Weight: 5, 0, 1, 3, 1
     411             :     !> RefinedWeight: 3, 0, 0, 2, 0
     412             :     !>          skip: 3
     413             :     !>        Weight: 5, 0, 1, 3, 1
     414             :     !> RefinedWeight: 2, 0, 0, 1, 1
     415             :     !> ```
     416         982 :     pure function getRefinedWeight(np,Weight,skip,refinedChainSize) result(RefinedWeight)
     417             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     418             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRefinedWeight
     419             : #endif
     420         354 :         use Math_mod, only: getCumSum
     421             :         use Constants_mod, only: IK, RK
     422             :         implicit none
     423             :         real(RK), intent(in)                :: skip
     424             :         integer(IK), intent(in)             :: np, Weight(np)
     425             :         integer(IK) , intent(in) , optional :: refinedChainSize
     426             :         integer(IK)                         :: RefinedWeight(np)
     427             :         integer(IK)                         :: ip, refinedChainSizeCounter, offset
     428      154641 :         real(RK)                            :: sumSkips, CumSumWeight(np)
     429             :         logical                             :: refinedChainSizeIsPresent
     430         491 :         refinedChainSizeIsPresent = present(refinedChainSize)
     431         491 :         if (refinedChainSizeIsPresent) refinedChainSizeCounter = 0_IK
     432      154150 :         CumSumWeight = real(getCumSum(np,Weight),kind=RK)
     433         491 :         sumSkips = skip
     434         491 :         offset = 1_IK
     435         491 :         ip = offset
     436      154150 :         RefinedWeight = 0_IK
     437      153602 :         loopOverAllSample: do
     438      111189 :             loopOverCurrentSample: do
     439      265282 :                 if (sumSkips>CumSumWeight(ip)) then
     440      154053 :                     exit loopOverCurrentSample
     441      111229 :                 elseif (refinedChainSizeIsPresent) then
     442       11386 :                     if (refinedChainSizeCounter==refinedChainSize) exit loopOverAllSample
     443       11346 :                     refinedChainSizeCounter = refinedChainSizeCounter + 1_IK
     444             :                 end if
     445      111189 :                 RefinedWeight(ip) = RefinedWeight(ip) + 1_IK
     446      111189 :                 sumSkips = sumSkips + skip
     447      111189 :                 cycle loopOverCurrentSample
     448             :             end do loopOverCurrentSample
     449      154053 :             if (ip==np) then
     450         491 :                 if (refinedChainSizeIsPresent) then
     451          60 :                     if (refinedChainSizeCounter<refinedChainSize) then
     452          40 :                         offset = offset + 1_IK
     453          40 :                         if (offset==np) offset = 1_IK
     454          40 :                         ip = offset
     455          40 :                         sumSkips = skip
     456          40 :                         if (offset/=1_IK) sumSkips = sumSkips + CumSumWeight(ip-1)
     457             :                         cycle loopOverAllSample ! LCOV_EXCL_LINE
     458             :                     end if
     459             :                 end if
     460         451 :                 exit loopOverAllSample
     461             :             end if
     462      153562 :             ip = ip + 1_IK
     463             :         end do loopOverAllSample
     464         491 :     end function getRefinedWeight
     465             : 
     466             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     467             : 
     468             :     !> Refined an input weighted sample according to the new requested weights.
     469             :     !>
     470             :     !> \param[in]   nd                  :   The number of dimensions of the input `Sample(0:nd,np)`.
     471             :     !> \param[in]   np                  :   The number of points in the input `Sample(0:nd,np)`.
     472             :     !> \param[in]   skip                :   The jump size with which the input chain has to be refined.
     473             :     !> \param[in]   Sample              :   The input 2-dimensional array of sampled states which has to be refined.
     474             :     !> \param[in]   Weight              :   The weights of the sampled points.
     475             :     !> \param[out]  RefinedChain        :   The refined array.
     476             :     !> \param[out]  RefinedWeight       :   The vector of refined weights corresponding to the output refined array.
     477             :     !> \param[out]  PointCount          :   An object of derived type [Count_type](@ref paramontechainfilecontents_mod::count_type)
     478             :     !>                                      containing the number of points in the refined sample.
     479             :     !> \param[in]   refinedChainSize    :   The requested refined sample size (**optional**). If the size of the refined sample is given as input,
     480             :     !>                                      then the requested sample is directly generated based on the input size.
     481         491 :     pure subroutine refineWeightedSample(nd,np,skip,Sample,Weight,RefinedChain,RefinedWeight,PointCount,refinedChainSize)
     482             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     483             :         !DEC$ ATTRIBUTES DLLEXPORT :: refineWeightedSample
     484             : #endif
     485         491 :         use Constants_mod, only: IK, RK
     486             :         implicit none
     487             :         integer(IK)     , intent(in)                :: nd, np, Weight(np)   !, skip
     488             :         real(RK)        , intent(in)                :: Sample(0:nd,np), skip
     489             :         real(RK)        , intent(out), allocatable  :: RefinedChain(:,:)
     490             :         integer(IK)     , intent(out), allocatable  :: RefinedWeight(:)
     491             :         integer(IK)     , intent(in) , optional     :: refinedChainSize
     492             :         type(Count_type), intent(out)               :: PointCount
     493             :         integer(IK)                                 :: ip, ipRefined, npRefined, UpdatedWeight(np) ! LCOV_EXCL_LINE
     494         491 :         UpdatedWeight = getRefinedWeight(np,Weight,skip,refinedChainSize)
     495      154150 :         npRefined = count(UpdatedWeight>0)
     496         491 :         allocate( RefinedChain(0:nd,npRefined) , RefinedWeight(npRefined) )
     497         491 :         ipRefined = 0
     498         491 :         PointCount%verbose = 0
     499      154150 :         do ip = 1, np
     500      154150 :             if (UpdatedWeight(ip)>0) then
     501       70786 :                 ipRefined = ipRefined + 1
     502      265576 :                 RefinedChain(0:nd,ipRefined) = Sample(0:nd,ip)
     503       70786 :                 RefinedWeight(ipRefined) = UpdatedWeight(ip)
     504       70786 :                 PointCount%verbose = PointCount%verbose + RefinedWeight(ipRefined)
     505             :             end if
     506             :         end do
     507         491 :         PointCount%compact = npRefined
     508         491 :     end subroutine refineWeightedSample
     509             : 
     510             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     511             : 
     512             :     !> Return the best skip size through a Markov chain to refined it to the optimal requested size.
     513             :     !>
     514             :     !> \param[in]   oldSampleSize   :   The original size of the Markov chain.
     515             :     !> \param[in]   newSampleSize   :   The final desired size of the refined sample.
     516             :     !>
     517             :     !> \return
     518             :     !> `skip4NewSampleSize` : The computed skip size.
     519             :     !>
     520             :     !> \warning
     521             :     !> The condition `oldSampleSize >= newSampleSize` must always hold,
     522             :     !> otherwise a negative value for `skip4NewSampleSize` will be returned to indicate the occurrence of an error.
     523          12 :     pure function getSkip4NewSampleSize(oldSampleSize,newSampleSize) result(skip4NewSampleSize)
     524             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     525             :         !DEC$ ATTRIBUTES DLLEXPORT :: getSkip4NewSampleSize
     526             : #endif
     527         491 :         use Constants_mod, only: IK, RK
     528             :         implicit none
     529             :         integer(IK) , intent(in)    :: oldSampleSize,newSampleSize
     530             :         integer(IK)                 :: skip4NewSampleSize, addition, quotient
     531          12 :         if (oldSampleSize < newSampleSize) then
     532           3 :             skip4NewSampleSize = -1_IK
     533           3 :             return
     534             :         end if
     535           9 :         addition = 1
     536           9 :         quotient = oldSampleSize / newSampleSize
     537           9 :         if (mod(oldSampleSize,newSampleSize)==0) addition = 0
     538           9 :         skip4NewSampleSize = quotient + addition
     539          12 :     end function getSkip4NewSampleSize
     540             : 
     541             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     542             : 
     543             :     !> Write the computed refined chain to the specified output file.
     544             :     !>
     545             :     !> \param[in]   RefinedChain                :   An object of class [RefinedChain_type](@ref refinedchain_type)
     546             :     !>                                              containing the refined sample to be written to the output file.
     547             :     !> \param[in]   sampleFileUnit              :   The unit of the file to which the sample must be written.
     548             :     !> \param[in]   sampleFileHeaderFormat      :   The IO format of the header of the sample file.
     549             :     !> \param[in]   sampleFileContentsFormat    :   The IO format of the contents (sampled states) in the sample file.
     550         261 :     subroutine writeRefinedChain(RefinedChain,sampleFileUnit,sampleFileHeaderFormat,sampleFileContentsFormat)
     551             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     552             :         !DEC$ ATTRIBUTES DLLEXPORT :: writeRefinedChain
     553             : #endif
     554             :         implicit none
     555             :         class(RefinedChain_type), intent(in)    :: RefinedChain
     556             :         integer(IK) , intent(in)                :: sampleFileUnit
     557             :         character(*), intent(in)                :: sampleFileHeaderFormat, sampleFileContentsFormat
     558             :         integer(IK)                             :: isample, iweight, i
     559         940 :         write(sampleFileUnit, sampleFileHeaderFormat) (trim(adjustl(RefinedChain%ColHeader(i)%record)),i=0,RefinedChain%ndim)
     560       45324 :         do isample = 1, RefinedChain%Count(RefinedChain%numRefinement)%compact
     561      154708 :             do iweight = 1, RefinedChain%Weight(isample)
     562      154447 :                 write(sampleFileUnit,sampleFileContentsFormat) RefinedChain%LogFuncState(0:RefinedChain%ndim,isample)
     563             :             end do
     564             :         end do
     565         261 :         flush(sampleFileUnit)
     566             :     end subroutine writeRefinedChain ! LCOV_EXCL_LINE
     567             : 
     568             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     569             : 
     570             :     !> Write the computed refined chain to the specified output file.
     571             :     !>
     572             :     !> \param[in]   sampleFilePath      :   The path to the input chain file that must be read.
     573             :     !> \param[in]   delimiter           :   The delimiter used in the file.
     574             :     !> \param[in]   ndim                :   The number of dimensions of the sampled states in the sample file.
     575             :     !>                                      This is basically, the size of the domain of the objective function.
     576             :     !> \param[in]   tenabled            :   An optional input logical value standing for `transpose-enabled`. If `.false.`,
     577             :     !>                                      the input data will be naturally read according to Fortran column-wise data storage
     578             :     !>                                      rule as a matrix of rank `0:nd * 1:np`. If `.false.`, the input sample file will be
     579             :     !>                                      read as a matrix of rank `1:np * 0:nd`. Note that `np` represents the number of rows
     580             :     !>                                      in the files (that is, the number of sampled points, whereas `nd` represents the
     581             :     !>                                      number of columns in the input file (**optional**, default = `.false.`).
     582             :     !>
     583             :     !> \return
     584             :     !> `RefinedChain` : An object of class [RefinedChain_type](@ref refinedchain_type) containing
     585             :     !>                  the sampled states read from the specified input file.
     586         176 :     function readRefinedChain(sampleFilePath,delimiter,ndim,tenabled) result(RefinedChain)
     587             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     588             :         !DEC$ ATTRIBUTES DLLEXPORT :: readRefinedChain
     589             : #endif
     590         261 :         use FileContents_mod, only: getNumRecordInFile
     591             :         use String_mod, only: String_type, num2str
     592             :         implicit none
     593             :         integer(IK) , intent(in)            :: ndim
     594             :         character(*), intent(in)            :: sampleFilePath, delimiter
     595             :         logical     , intent(in), optional  :: tenabled
     596             :         type(RefinedChain_type)             :: RefinedChain
     597             :         character(*), parameter             :: PROCEDURE_NAME = MODULE_NAME//"@readRefinedChain()"
     598             :         integer(IK)                         :: sampleFileUnit, isample, i
     599             :         logical                             :: tenabledDefault
     600         176 :         type(String_type)                   :: Record
     601             : 
     602             :         if (allocated(Record%value)) deallocate(Record%value) ! LCOV_EXCL_LINE
     603         176 :         allocate( character(99999) :: Record%value )
     604             : 
     605         176 :         RefinedChain%ndim = ndim
     606         176 :         RefinedChain%numRefinement = 0_IK
     607         352 :         allocate(RefinedChain%Count(RefinedChain%numRefinement:RefinedChain%numRefinement)) ! allocate just the zeroth level of `RefinedChain`.
     608             : 
     609             :         ! find the number of lines in the sample file
     610             : 
     611         176 :         call getNumRecordInFile(filePath=sampleFilePath, numRecord=RefinedChain%Count(RefinedChain%numRefinement)%verbose, Err=RefinedChain%Err)
     612         176 :         if (RefinedChain%Err%occurred) then
     613             :             ! LCOV_EXCL_START
     614             :             RefinedChain%Err%msg = PROCEDURE_NAME // RefinedChain%Err%msg
     615             :             return
     616             :             ! LCOV_EXCL_STOP
     617             :         end if
     618             : 
     619         176 :         RefinedChain%Count(RefinedChain%numRefinement)%verbose = RefinedChain%Count(RefinedChain%numRefinement)%verbose - 1_IK ! remove header from the count
     620             : 
     621             :         open( newunit = sampleFileUnit &
     622             :             , file = sampleFilePath &
     623             :             , status = "old" &
     624             : #if defined INTEL_COMPILER_ENABLED && defined OS_IS_WINDOWS
     625             :             , SHARED &
     626             : #endif
     627         176 :             )
     628             : 
     629             :         ! read header
     630             : 
     631         176 :         read(sampleFileUnit,"(A)") Record%value
     632        1372 :         Record%Parts = Record%split(Record%value,delimiter,Record%nPart)
     633         176 :         if (Record%nPart/=ndim+1_IK) then
     634             :             ! LCOV_EXCL_START
     635             :             RefinedChain%Err%occurred = .true.
     636             :             RefinedChain%Err%msg = PROCEDURE_NAME // ": The number of header columns ("//num2str(Record%nPart)//") is not equal to ndim + 1: "//num2str(ndim+1_IK)
     637             :             return
     638             :             ! LCOV_EXCL_STOP
     639             :         end if
     640             : 
     641         686 :         allocate(RefinedChain%ColHeader(0:ndim))
     642         686 :         do i = 0, ndim
     643         686 :             RefinedChain%ColHeader(i)%record = Record%Parts(i+1)%record
     644             :         end do
     645             : 
     646             :         ! read contents
     647             : 
     648         176 :         tenabledDefault = .false.
     649         176 :         if (present(tenabled)) tenabledDefault = tenabled
     650             : 
     651         176 :         if (tenabledDefault) then
     652             : 
     653         168 :             allocate( RefinedChain%LogFuncState(RefinedChain%Count(RefinedChain%numRefinement)%verbose, 0:ndim) )
     654      243234 :             do isample = 1, RefinedChain%Count(RefinedChain%numRefinement)%verbose
     655      243066 :                 read(sampleFileUnit, "(A)") Record%value
     656     1213050 :                 Record%Parts = Record%split(trim(adjustl(Record%value)),delimiter)
     657      970155 :                 do i = 0, ndim
     658      969987 :                     read(Record%Parts(i+1)%record,*) RefinedChain%LogFuncState(isample,i)
     659             :                 end do
     660             :             end do
     661             : 
     662             :         else
     663             : 
     664           8 :             allocate( RefinedChain%LogFuncState(0:ndim, RefinedChain%Count(RefinedChain%numRefinement)%verbose) )
     665          88 :             do isample = 1, RefinedChain%Count(RefinedChain%numRefinement)%verbose
     666          80 :                 read(sampleFileUnit, "(A)") Record%value
     667         400 :                 Record%Parts = Record%split(trim(adjustl(Record%value)),delimiter)
     668         328 :                 do i = 0, ndim
     669         320 :                     read(Record%Parts(i+1)%record,*) RefinedChain%LogFuncState(i,isample)
     670             :                 end do
     671             :             end do
     672             : 
     673             :         end if
     674             : 
     675         176 :         close(sampleFileUnit)
     676             : 
     677         862 :     end function readRefinedChain
     678             : 
     679             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     680             : 
     681             : end module ParaMCMCRefinedChain_mod ! LCOV_EXCL_LINE

ParaMonte: Plain Powerful Parallel Monte Carlo Library 
The Computational Data Science Lab
© Copyright 2012 - 2021