The ParaMonte Documentation Website
Current view: top level - kernel - ParaDXXXProposal_mod.inc.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: MPI Parallel Kernel - Code Coverage Report Lines: 196 196 100.0 %
Date: 2021-01-08 13:07:16 Functions: 32 32 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
      44             : !> This file implements the body of the Proposal modules of the `ParaDRAM` and `ParaDISE` samplers.
      45             : !>
      46             : !> \remark
      47             : !> This module requires preprocessing, prior to compilation.
      48             : !>
      49             : !> \author Amir Shahmoradi
      50             : 
      51             : #if defined UNIFORM
      52             : 
      53             : #define GET_RANDOM_PROPOSAL getRandMVU
      54             : 
      55             : #elif defined NORMAL
      56             : 
      57             : #define GET_RANDOM_PROPOSAL getRandMVN
      58             : 
      59             : #else
      60             : 
      61             : #error "Unknown Proposal model in ParaDXXXProposal_mod.inc.f90"
      62             : 
      63             : #endif
      64             : 
      65             : 
      66             : #if defined PARADRAM
      67             : 
      68             : #define ParaDXXX ParaDRAM
      69             :     use ParaDRAMProposalAbstract_mod, only: ProposalAbstract_type, ProposalErr
      70             : 
      71             : #elif defined PARADISE
      72             : 
      73             : #define ParaDXXX ParaDISE
      74             :     use ParaDISEProposalAbstract_mod, only: ProposalAbstract_type, ProposalErr
      75             : 
      76             : #endif
      77             : 
      78             :     use ParaMonte_mod, only: Image_type
      79             :     use Constants_mod, only: IK, RK, PMSM
      80             :     use String_mod, only: IntStr_type
      81             : 
      82             :     implicit none
      83             : 
      84             :     private
      85             :     public :: Proposal_type
      86             : 
      87             : #if defined NORMAL
      88             :     character(*), parameter         :: MODULE_NAME = "@"//PMSM%ParaDXXX//"ProposalNormal_mod"
      89             : #elif defined UNIFORM
      90             :     character(*), parameter         :: MODULE_NAME = "@"//PMSM%ParaDXXX//"ProposalUniform_mod"
      91             : #endif
      92             : 
      93             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      94             : 
      95             :     type                            :: AccRate_type
      96             :         real(RK)                    :: sumUpToLastUpdate
      97             :         real(RK)                    :: target
      98             :     end type AccRate_type
      99             : 
     100             :     !> The `Proposal_type` class.
     101             :     type, extends(ProposalAbstract_type) :: Proposal_type
     102             :        !type(AccRate_type)          :: AccRate
     103             :     contains
     104             :         procedure   , nopass        :: getNew
     105             : #if defined PARADISE
     106             :         procedure   , nopass        :: getLogProb
     107             : #endif
     108             :         procedure   , nopass        :: doAdaptation
     109             :        !procedure   , nopass        :: readRestartFile
     110             :        !procedure   , nopass        :: writeRestartFile
     111             : #if defined CAF_ENABLED || defined MPI_ENABLED
     112             :         procedure   , nopass        :: bcastAdaptation
     113             : #endif
     114             :     end type Proposal_type
     115             : 
     116             :     interface Proposal_type
     117             :         module procedure :: constructProposalSymmetric
     118             :     end interface Proposal_type
     119             : 
     120             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     121             : 
     122             :     ! Covariance Matrix of the proposal distribution. Last index belongs to delayed rejection
     123             : 
     124             : #if defined CAF_ENABLED
     125             :     real(RK)        , save  , allocatable   :: comv_CholDiagLower(:,:,:)[:]
     126             : #else
     127             :     real(RK)        , save  , allocatable   :: comv_CholDiagLower(:,:,:)
     128             : #endif
     129             :     real(RK)        , save  , allocatable   :: mv_logSqrtDetInvCovMat(:)
     130             :     real(RK)        , save  , allocatable   :: mv_InvCovMat(:,:,:)
     131             : 
     132             : #if defined MPI_ENABLED
     133             :     integer(IK)     , save                  :: mc_ndimSqPlusNdim
     134             : #endif
     135             :     type(Image_type), save                  :: mc_Image
     136             :     integer(IK)     , save                  :: mc_ndim
     137             :     integer(IK)     , save                  :: mc_logFileUnit
     138             :     integer(IK)     , save                  :: mc_restartFileUnit
     139             :     logical         , save                  :: mc_scalingRequested
     140             :     real(RK)        , save                  :: mc_defaultScaleFactorSq
     141             :     integer(IK)     , save                  :: mc_DelayedRejectionCount
     142             :     integer(IK)     , save                  :: mc_MaxNumDomainCheckToWarn
     143             :     integer(IK)     , save                  :: mc_MaxNumDomainCheckToStop
     144             :     logical         , save                  :: mc_delayedRejectionRequested
     145             :     real(RK)        , save                  :: mc_ndimInverse
     146             :     real(RK)        , save                  :: mc_targetAcceptanceRate
     147             :     real(RK)        , save                  :: mc_TargetAcceptanceRateLimit(2)
     148             :    !real(RK)        , save                  :: mc_maxScaleFactor != 2._RK
     149             :    !real(RK)        , save                  :: mc_maxScaleFactorSq != mc_maxScaleFactor**2
     150             :     real(RK)        , save  , allocatable   :: mc_DelayedRejectionScaleFactorVec(:)
     151             :     real(RK)        , save  , allocatable   :: mc_DomainLowerLimitVec(:)
     152             :     real(RK)        , save  , allocatable   :: mc_DomainUpperLimitVec(:)
     153             :     logical         , save  , allocatable   :: mc_isAsciiRestartFileFormat
     154             :     logical         , save  , allocatable   :: mc_isBinaryRestartFileFormat
     155             :     character(:)    , save  , allocatable   :: mc_MaxNumDomainCheckToWarnMsg
     156             :     character(:)    , save  , allocatable   :: mc_MaxNumDomainCheckToStopMsg
     157             :     character(:)    , save  , allocatable   :: mc_negativeTotalVariationMsg
     158             :     character(:)    , save  , allocatable   :: mc_restartFileFormat
     159             :     character(:)    , save  , allocatable   :: mc_methodBrand
     160             :     character(:)    , save  , allocatable   :: mc_methodName
     161             : #if defined UNIFORM
     162             :     real(RK)        , save  , allocatable   :: mc_negLogVolUnitBall
     163             : #endif
     164             : 
     165             :     ! the following had to be defined globally for the sake of restart file generation
     166             : 
     167             :     real(RK)        , save  , allocatable   :: mv_MeanOld_save(:)
     168             :     real(RK)        , save                  :: mv_logSqrtDetOld_save
     169             :     real(RK)        , save                  :: mv_adaptiveScaleFactorSq_save    ! = 1._RK
     170             :     integer(IK)     , save                  :: mv_sampleSizeOld_save            ! = 0_IK
     171             : 
     172             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     173             : 
     174             : contains
     175             : 
     176             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     177             : 
     178             :     !> This is the constructor of the [ParaDXXXProposalNormal_type](@ref paradxxxproposalnormal_type)
     179             :     !> or [ParaDXXXProposalUniform_type](@ref paradxxxproposaluniform_type) classes.\n
     180             :     !> Return an object of [ParaDXXXProposalNormal_type](@ref paradxxxproposalnormal_type)
     181             :     !> or [ParaDXXXProposalUniform_type](@ref paradxxxproposaluniform_type) classes given the input simulation characteristics.
     182             :     !>
     183             :     !> \remark
     184             :     !> This interface madness is a result of the internal compiler bug in GFortran as of Jan 2020, which diagnoses a `ParaDXXX_type`
     185             :     !> argument as circular dependency due to this constructor appearing in the type-bound setup procedure of `ParaDXXX_type`.
     186             :     !> Intel does not complain. Until GFortran comes up with a fix, we have to live with this interface.
     187         687 :     function constructProposalSymmetric ( ndim &
     188             :                                         , SpecBase &
     189             :                                         , SpecMCMC &
     190             :                                         , SpecDRAM &
     191             :                                         , Image &
     192             :                                         , name &
     193             :                                         , brand &
     194             :                                         , LogFile &
     195             :                                         , RestartFile &
     196             :                                         , isFreshRun &
     197             :                                         ) result(Proposal)
     198             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     199             :         !DEC$ ATTRIBUTES DLLEXPORT :: constructProposalSymmetric
     200             : #endif
     201             :         use Constants_mod, only: IK, RK, NULL_RK
     202             :         use ParaMonte_mod, only: Image_type
     203             :         use ParaMonte_mod, only: LogFile_type
     204             :         use ParaMonte_mod, only: RestartFile_type
     205             :         use SpecBase_mod, only: SpecBase_type
     206             :         use SpecMCMC_mod, only: SpecMCMC_type
     207             :         use SpecDRAM_mod, only: SpecDRAM_type
     208             :         use Matrix_mod, only: getCholeskyFactor
     209             :         use String_mod, only: num2str
     210             :         use Err_mod, only: abort
     211             : #if defined UNIFORM
     212             :         use Math_mod, only: getLogVolUnitBall
     213             : #endif
     214             :         implicit none
     215             : 
     216             :         integer(IK)             , intent(in)    :: ndim
     217             :         type(SpecBase_type)     , intent(in)    :: SpecBase
     218             :         type(SpecMCMC_type)     , intent(in)    :: SpecMCMC
     219             :         type(SpecDRAM_type)     , intent(in)    :: SpecDRAM
     220             :         type(Image_type)        , intent(in)    :: Image
     221             :         character(*)            , intent(in)    :: name
     222             :         character(*)            , intent(in)    :: brand
     223             :         type(LogFile_type)      , intent(in)    :: LogFile
     224             :         type(RestartFile_type)  , intent(in)    :: RestartFile
     225             :         logical                                 :: isFreshRun
     226             : 
     227             :         type(Proposal_type)                     :: Proposal
     228             : 
     229             :         character(*), parameter                 :: PROCEDURE_NAME = MODULE_NAME // "@constructProposalSymmetric()"
     230             :         integer                                 :: i, j
     231             : 
     232         687 :         ProposalErr%occurred = .false.
     233             : 
     234             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     235             :         ! setup sampler update global save variables
     236             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     237             : 
     238         687 :         if (allocated(mv_MeanOld_save)) deallocate(mv_MeanOld_save); allocate(mv_MeanOld_save(ndim))
     239        1749 :         mv_MeanOld_save(1:ndim) = SpecMCMC%StartPointVec%Val
     240         687 :         mv_logSqrtDetOld_save   = NULL_RK
     241         687 :         mv_sampleSizeOld_save   = 1_IK
     242         687 :         mv_adaptiveScaleFactorSq_save = 1._RK
     243             : 
     244             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     245             :         ! setup general proposal specifications
     246             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     247             : 
     248             : #if defined MPI_ENABLED
     249         687 :         mc_ndimSqPlusNdim                   = ndim*(ndim+1_IK)
     250             : #endif
     251         687 :         mc_ndim                             = ndim
     252        2319 :         mc_DomainLowerLimitVec              = SpecBase%DomainLowerLimitVec%Val
     253        2319 :         mc_DomainUpperLimitVec              = SpecBase%DomainUpperLimitVec%Val
     254        1473 :         mc_DelayedRejectionScaleFactorVec   = SpecDRAM%delayedRejectionScaleFactorVec%Val
     255             :        !mc_isNormal                         = SpecMCMC%ProposalModel%isNormal
     256         687 :         mc_Image                            = Image
     257         687 :         mc_methodName                       = name
     258         687 :         mc_methodBrand                      = brand
     259         687 :         mc_logFileUnit                      = LogFile%unit
     260         687 :         mc_restartFileUnit                  = RestartFile%unit
     261         687 :         mc_restartFileFormat                = RestartFile%format
     262         687 :         mc_isBinaryRestartFileFormat        = SpecBase%RestartFileFormat%isBinary
     263         687 :         mc_isAsciiRestartFileFormat         = SpecBase%RestartFileFormat%isAscii
     264         687 :         mc_defaultScaleFactorSq             = SpecMCMC%ScaleFactor%val**2
     265             :        !Proposal%AccRate%sumUpToLastUpdate  = 0._RK
     266         687 :         mc_maxNumDomainCheckToWarn          = SpecBase%MaxNumDomainCheckToWarn%val
     267         687 :         mc_maxNumDomainCheckToStop          = SpecBase%MaxNumDomainCheckToStop%val
     268         687 :         mc_delayedRejectionCount            = SpecDRAM%DelayedRejectionCount%val
     269         687 :         mc_delayedRejectionRequested        = mc_DelayedRejectionCount > 0_IK
     270         687 :         mc_scalingRequested                 = SpecBase%TargetAcceptanceRate%scalingRequested
     271         687 :         mc_ndimInverse                      = 1._RK/real(ndim,kind=RK)
     272        2061 :         mc_TargetAcceptanceRateLimit        = SpecBase%TargetAcceptanceRate%Val
     273         687 :         mc_targetAcceptanceRate             = 0.234_RK
     274         795 :         if (mc_scalingRequested) mc_targetAcceptanceRate = sum(mc_TargetAcceptanceRateLimit) / 2._RK
     275             :        !mc_maxScaleFactorSq                 = 4._RK**mc_ndimInverse
     276             :        !mc_maxScaleFactor                   = sqrt(mc_maxScaleFactorSq)
     277             : 
     278             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     279             :         ! setup ProposalSymmetric specifications
     280             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     281             : 
     282             :         ! setup covariance matrix
     283             : 
     284         687 :         if (allocated(mv_InvCovMat)) deallocate(mv_InvCovMat)
     285         687 :         allocate( mv_InvCovMat(ndim,0:ndim,0:mc_DelayedRejectionCount) )
     286         687 :         if (allocated(mv_logSqrtDetInvCovMat)) deallocate(mv_logSqrtDetInvCovMat)
     287         687 :         allocate( mv_logSqrtDetInvCovMat(0:mc_DelayedRejectionCount) )
     288             : 
     289         687 :         if (allocated(comv_CholDiagLower)) deallocate(comv_CholDiagLower)
     290             : #if defined CAF_ENABLED
     291             :         ! on the second dimension, the zeroth index refers to the Diagonal elements of the Cholesky lower triangular matrix
     292             :         ! This rearrangement was done for more efficient communication of the matrix across processes.
     293             :         allocate( comv_CholDiagLower(ndim,0:ndim,0:mc_DelayedRejectionCount)[*] )
     294             : #else
     295         687 :         allocate( comv_CholDiagLower(ndim,0:ndim,0:mc_DelayedRejectionCount) )
     296             : #endif
     297             : 
     298        3561 :         comv_CholDiagLower(1:ndim,1:ndim,0) = SpecMCMC%ProposalStartCovMat%Val
     299             : 
     300             :         ! Now scale the covariance matrix
     301             : 
     302        1749 :         do j = 1, ndim
     303        3186 :             do i = 1, j
     304        2499 :                 comv_CholDiagLower(i,j,0) = comv_CholDiagLower(i,j,0) * mc_defaultScaleFactorSq
     305             :             end do
     306             :         end do
     307             : 
     308             :         ! Now get the Cholesky Factor of the Covariance Matrix. Lower comv_CholDiagLower will be the CholFac
     309             : 
     310         687 :         call getCholeskyFactor( ndim, comv_CholDiagLower(:,1:ndim,0), comv_CholDiagLower(1:ndim,0,0) ) ! The `:` instead of `1:ndim` avoids temporary array creation.
     311         687 :         if (comv_CholDiagLower(1,0,0)<0._RK) then
     312             :         ! LCOV_EXCL_START
     313             :             ProposalErr%msg = mc_Image%name // PROCEDURE_NAME // ": Singular input covariance matrix by user was detected. This is strange.\nCovariance matrix lower triangle:"
     314             :             do j = 1, ndim
     315             :                 do i = 1, j
     316             :                     ProposalErr%msg = ProposalErr%msg // "\n" // num2str(comv_CholDiagLower(1:i,j,0))
     317             :                 end do
     318             :             end do
     319             :             ProposalErr%occurred = .true.
     320             :             call abort( Err = ProposalErr, prefix = mc_methodBrand, newline = "\n", outputUnit = mc_logFileUnit )
     321             :             return
     322             :         ! LCOV_EXCL_STOP
     323             :         end if
     324             : 
     325         687 :         if (mc_delayedRejectionRequested) call updateDelRejCholDiagLower()
     326         687 :         call getInvCovMat()
     327        1749 :         mv_logSqrtDetOld_save = sum(log( comv_CholDiagLower(1:ndim,0,0) ))
     328             : 
     329             :         ! Scale the higher-stage delayed-rejection Cholesky Lower matrices
     330             : 
     331         687 :         if (mc_delayedRejectionRequested) call updateDelRejCholDiagLower()
     332             : 
     333             :         ! This will be used for Domain boundary checking during the simulation
     334             : 
     335         687 :         mc_negativeTotalVariationMsg = MODULE_NAME//"@doAdaptation(): Non-positive adaptation measure detected, possibly due to round-off error: "
     336             :         mc_MaxNumDomainCheckToWarnMsg = MODULE_NAME//"@getNew(): "//num2str(mc_MaxNumDomainCheckToWarn)//&
     337         687 :                                         " proposals were drawn out of the objective function's Domain without any acceptance."
     338             :         mc_MaxNumDomainCheckToStopMsg = MODULE_NAME//"@getNew(): "//num2str(mc_MaxNumDomainCheckToStop)//&
     339             :                                         " proposals were drawn out of the objective function's Domain. As per the value set for the&
     340         687 :                                         & simulation specification variable 'maxNumDomainCheckToStop', "//mc_methodName//" will abort now."
     341             : 
     342             : #if defined UNIFORM
     343         120 :         if (ndim>1_IK) then
     344         108 :             mc_negLogVolUnitBall = -getLogVolUnitBall(ndim)
     345             :         else
     346          12 :             mc_negLogVolUnitBall = 0._RK
     347             :         end if
     348             : #endif
     349             : 
     350             :         ! read/write the first entry of the restart file
     351             : 
     352         687 :         if (mc_Image%isLeader) then
     353             :             block
     354         269 :                 real(RK) :: meanAccRateSinceStart
     355         269 :                 if (isFreshRun) then
     356         237 :                     call writeRestartFile(meanAccRateSinceStart=1._RK)
     357         237 :                     call writeRestartFile()
     358             :                 else
     359          32 :                     call readRestartFile(meanAccRateSinceStart)
     360          32 :                     call readRestartFile()
     361             :                 end if
     362             :             end block
     363             :         end if
     364             : 
     365         687 :     end function constructProposalSymmetric
     366             : 
     367             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     368             : 
     369             :     !> \brief
     370             :     !> This procedure is a static method of the [ParaDXXXProposalNormal_type](@ref paradxxxproposalnormal_type)
     371             :     !> or [ParaDXXXProposalUniform_type](@ref paradxxxproposaluniform_type) classes.\n
     372             :     !> Return a newly sampled state to the kernel routine for acceptance check.
     373             :     !>
     374             :     !>
     375             :     !> @param[in]   nd          :   The number of dimensions of the objective function.
     376             :     !> @param[in]   counterDRS  :   The Delayed Rejection Stage counter.
     377             :     !> @param[in]   StateOld    :   The old sampled state.
     378             :     !>
     379             :     !> \return
     380             :     !> `StateNew` : The newly sampled state.
     381      416105 :     function getNew ( nd            &
     382             :                     , counterDRS    &
     383      416105 :                     , StateOld      &
     384     1566938 :                     ) result (StateNew)
     385             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     386             :         !DEC$ ATTRIBUTES DLLEXPORT :: getNew
     387             : #endif
     388         687 :         use Statistics_mod, only: GET_RANDOM_PROPOSAL
     389             :         use Constants_mod, only: IK, RK
     390             :         use Err_mod, only: warn, abort
     391             :         use ParaMonteLogFunc_mod, only: getLogFunc_proc
     392             : 
     393             :         implicit none
     394             : 
     395             :         character(*), parameter                         :: PROCEDURE_NAME = MODULE_NAME // "@getNew()"
     396             : 
     397             :         integer(IK), intent(in)                         :: nd
     398             :         integer(IK), intent(in)                         :: counterDRS
     399             :         real(RK)   , intent(in)                         :: StateOld(nd)
     400             :         real(RK)                                        :: StateNew(nd)
     401             :         integer(IK)                                     :: domainCheckCounter
     402             : 
     403      416105 :         domainCheckCounter = 0_IK
     404             : 
     405       65279 :         loopBoundaryCheck: do ! Check for the support Region consistency:
     406             : #if defined UNIFORM || defined NORMAL
     407             :             StateNew(1:nd) = GET_RANDOM_PROPOSAL( nd                                    & ! LCOV_EXCL_LINE
     408             :                                                 , StateOld                              & ! LCOV_EXCL_LINE
     409             :                                                 ! ATTN: The colon index in place of 1:nd below avoids the temporary array creation
     410             :                                                 , comv_CholDiagLower(:,1:nd,counterDRS) & ! LCOV_EXCL_LINE
     411             :                                                 , comv_CholDiagLower(1:nd,0,counterDRS) & ! LCOV_EXCL_LINE
     412      481384 :                                                 )
     413             : #endif
     414     2120606 :             if ( any(StateNew(1:nd)<=mc_DomainLowerLimitVec) .or. any(StateNew(1:nd)>=mc_DomainUpperLimitVec) ) then
     415       65297 :                 domainCheckCounter = domainCheckCounter + 1
     416       65330 :                 if (domainCheckCounter==mc_MaxNumDomainCheckToWarn) then
     417          33 :                     call warn( prefix = mc_methodBrand, outputUnit = mc_logFileUnit, msg = mc_MaxNumDomainCheckToWarnMsg )
     418             :                 end if
     419       65297 :                 if (domainCheckCounter==mc_MaxNumDomainCheckToStop) then
     420          18 :                     ProposalErr%occurred = .true.
     421          18 :                     ProposalErr%msg = mc_MaxNumDomainCheckToStopMsg
     422             : #if !defined CODECOV_ENABLED && ((!defined MATLAB_ENABLED && !defined PYTHON_ENABLED && !defined R_ENABLED) || defined CAF_ENABLED && defined MPI_ENABLED )
     423             :                     call abort( Err = ProposalErr, prefix = mc_methodBrand, newline = "\n", outputUnit = mc_logFileUnit )
     424             : #endif
     425          18 :                     return
     426             :                 end if
     427       65279 :                 cycle loopBoundaryCheck
     428             :             end if
     429      416087 :             exit loopBoundaryCheck
     430             :         end do loopBoundaryCheck
     431             : 
     432      832210 :     end function getNew
     433             : 
     434             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     435             : 
     436             : #if defined PARADISE
     437             :     !> \brief
     438             :     !> This procedure is a static method of the [ParaDXXXProposalNormal_type](@ref paradxxxproposalnormal_type)
     439             :     !> or [ParaDXXXProposalUniform_type](@ref paradxxxproposaluniform_type) classes.\n
     440             :     !> Return the natural logarithm of the probability of acceptance of the new state given the old state.
     441             :     !>
     442             :     !> @param[in]   nd          :   The number of dimensions of the objective function.
     443             :     !> @param[in]   counterDRS  :   The Delayed Rejection Stage counter.
     444             :     !> @param[in]   StateOld    :   The old sampled state.
     445             :     !> @param[in]   StateNew    :   The new sampled state.
     446             :     !>
     447             :     !> \return
     448             :     !> `logProb` : The log probability of obtaining obtaining the new sample given the old sample.
     449             :     ! LCOV_EXCL_START
     450             :     pure function getLogProb( nd                &
     451             :                             , counterDRS        &
     452             :                             , StateOld          &
     453             :                             , StateNew          &
     454             :                             ) result(logProb)
     455             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     456             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProb
     457             : #endif
     458             : #if defined NORMAL
     459             :         use Statistics_mod, only: getLogProbMVN
     460             : #elif defined UNIFORM
     461             :         use Statistics_mod, only: isInsideEllipsoid
     462             : #endif
     463             :         use Constants_mod, only: IK, RK, NEGINF_RK
     464             :         implicit none
     465             :         integer(IK), intent(in)             :: nd
     466             :         integer(IK), intent(in)             :: counterDRS
     467             :         real(RK)   , intent(in)             :: StateOld(nd)
     468             :         real(RK)   , intent(in)             :: StateNew(nd)
     469             :         real(RK)                            :: logProb
     470             : #if defined UNIFORM
     471             :             if (isInsideEllipsoid(nd,StateNew-StateOld,mv_InvCovMat(1:mc_ndim,1:mc_ndim,counterDRS))) then
     472             :                 logProb = mc_negLogVolUnitBall + mv_logSqrtDetInvCovMat(counterDRS)
     473             :             else
     474             :                 logProb = NEGINF_RK
     475             :             end if
     476             : #elif defined NORMAL
     477             :             logProb = getLogProbMVN ( nd = nd &
     478             :                                     , MeanVec = StateOld &
     479             :                                     , InvCovMat = mv_InvCovMat(1:nd,1:nd,counterDRS) &
     480             :                                     , logSqrtDetInvCovMat = mv_logSqrtDetInvCovMat(counterDRS) &
     481             :                                     , Point = StateNew &
     482             :                                     )
     483             : #endif
     484             :     end function getLogProb
     485             :     ! LCOV_EXCL_STOP
     486             : #endif
     487             : 
     488             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     489             : 
     490             :     !> \brief
     491             :     !> This procedure is a static method of the [ParaDXXXProposalNormal_type](@ref paradxxxproposalnormal_type)
     492             :     !> or [ParaDXXXProposalUniform_type](@ref paradxxxproposaluniform_type) classes.\n
     493             :     !> Perform adaptation of the proposal distribution of the MCMC sampler.
     494             :     !>
     495             :     !> @param[in]       nd                      :   The number of dimensions of the objective function.
     496             :     !> @param[in]       chainSize               :   The size of the input chain in compact format.
     497             :     !> @param[in]       Chain                   :   The two dimensional `(nd, chainSize)` array of accepted states.
     498             :     !> @param[in]       ChainWeight             :   The one dimensional integer array of the weights of the sampled states in the chain.
     499             :     !> @param[in]       isFreshRun              :   The logical flag indicating whether the simulation is in regular (non-restart) mode.
     500             :     !> @param[in]       samplerUpdateIsGreedy   :   The logical flag indicating whether the adaptation is in greedy mode..
     501             :     !> @param[inout]    meanAccRateSinceStart   :   The mean acceptance rate since the start of the sampling.
     502             :     !> @param[out]      samplerUpdateSucceeded  :   The output logical flag indicating whether the sampling succeeded.
     503             :     !> @param[out]      adaptationMeasure       :   The output real number in the range `[0,1]` indicating the amount of adaptation,
     504             :     !>                                              with zero indicating no adaptation and one indicating extreme adaptation to the extent
     505             :     !>                                              that the new adapted proposal distribution is completely different from the previous proposal.
     506             :     !> \warning
     507             :     !> This routine must be exclusively called by the leader images.
     508             :     !>
     509             :     !> \remark
     510             :     !> For information on the meaning of `adaptationMeasure`, see the paper by Shahmoradi and Bagheri, 2020, whose PDF is available at:
     511             :     !> [https://www.cdslab.org/paramonte/notes/overview/preface/#the-paradram-sampler](https://www.cdslab.org/paramonte/notes/overview/preface/#the-paradram-sampler)
     512      122107 :     subroutine doAdaptation ( nd                        &
     513             :                             , chainSize                 &
     514      122107 :                             , Chain                     &
     515      122107 :                             , ChainWeight               &
     516             :                             , isFreshRun                &
     517             :                             , samplerUpdateIsGreedy     &
     518             :                             , meanAccRateSinceStart     &
     519             :                             , samplerUpdateSucceeded    &
     520             :                             , adaptationMeasure         &
     521             :                             )
     522             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     523             :         !DEC$ ATTRIBUTES DLLEXPORT :: doAdaptation
     524             : #endif
     525             : 
     526             :         use Statistics_mod, only: getSamCovUpperMeanTrans, getWeiSamCovUppMeanTrans, mergeMeanCovUpper ! LCOV_EXCL_LINE
     527             :         use Matrix_mod, only: getCholeskyFactor, getLogSqrtDetPosDefMat
     528             :         use Constants_mod, only: RK, IK ! , EPS_RK
     529             :         use String_mod, only: num2str
     530             :         use Err_mod, only: abort, warn
     531             :         implicit none
     532             : 
     533             :         character(*), parameter                         :: PROCEDURE_NAME = MODULE_NAME // "@doAdaptation()"
     534             : 
     535             :         integer(IK), intent(in)                         :: nd
     536             :         integer(IK), intent(in)                         :: chainSize
     537             :         real(RK)   , intent(in)                         :: Chain(nd,chainSize)
     538             :         integer(IK), intent(in)                         :: ChainWeight(chainSize)
     539             :         logical    , intent(in)                         :: isFreshRun
     540             :         logical    , intent(in)                         :: samplerUpdateIsGreedy
     541             :         real(RK)   , intent(inout)                      :: meanAccRateSinceStart ! is intent(out) in restart mode, intent(in) in fresh mode.
     542             :         logical    , intent(out)                        :: samplerUpdateSucceeded
     543             :         real(RK)   , intent(out)                        :: adaptationMeasure
     544             : 
     545             :         integer(IK)                                     :: i, j
     546      352663 :         real(RK)                                        :: MeanNew(nd)
     547      474770 :         real(RK)                                        :: MeanCurrent(nd)
     548      922224 :         real(RK)                                        :: CovMatUpperOld(nd,nd)
     549      922224 :         real(RK)                                        :: CovMatUpperCurrent(nd,nd)
     550      122107 :         real(RK)                                        :: logSqrtDetNew
     551      122107 :         real(RK)                                        :: logSqrtDetSum
     552      122107 :         real(RK)                                        :: adaptiveScaleFactor
     553             :         logical                                         :: adaptationMeasureComputationNeeded ! only used to avoid redundant affinity computation, if no update occurs
     554             :         logical                                         :: singularityOccurred, scalingNeeded
     555             :         integer(IK)                                     :: sampleSizeOld, sampleSizeCurrent
     556             : 
     557      122107 :         scalingNeeded = .false.
     558      122107 :         sampleSizeOld = mv_sampleSizeOld_save ! this is kept only for restoration of mv_sampleSizeOld_save, if needed.
     559             : 
     560             :         ! read/write meanAccRateSinceStart from/to restart file
     561             : 
     562      122107 :         if (mc_Image%isLeader) then
     563      122107 :             if (isFreshRun) then
     564       74420 :                 call writeRestartFile(meanAccRateSinceStart)
     565             :             else
     566       47687 :                 call readRestartFile(meanAccRateSinceStart)
     567             :             end if
     568             :         end if
     569             : 
     570             :         ! First if there are less than nd+1 points for new covariance computation, then just scale the covariance and return
     571             : 
     572      122107 :         blockSufficientSampleSizeCheck: if (chainSize>nd) then
     573             : 
     574             :             ! get the upper covariance matrix and Mean of the new sample
     575             : 
     576       30429 :             if (samplerUpdateIsGreedy) then
     577          43 :                 sampleSizeCurrent = chainSize
     578          43 :                 call getSamCovUpperMeanTrans(sampleSizeCurrent,nd,Chain,CovMatUpperCurrent,MeanCurrent)
     579             :             else
     580      138542 :                 sampleSizeCurrent = sum(ChainWeight)
     581       30386 :                 call getWeiSamCovUppMeanTrans(chainSize,sampleSizeCurrent,nd,Chain,ChainWeight,CovMatUpperCurrent,MeanCurrent)
     582             :             end if
     583             : 
     584             :             !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     585             : 
     586             :             ! combine old and new covariance matrices if both exist
     587             : 
     588       30429 :             blockMergeCovMat: if (mv_sampleSizeOld_save==1_IK) then
     589             : 
     590             :                 ! There is no prior old Covariance matrix to combine with the new one from the new chain
     591             : 
     592         665 :                 mv_MeanOld_save(1:nd) = MeanCurrent
     593         257 :                 mv_sampleSizeOld_save = sampleSizeCurrent
     594             : 
     595             :                 ! copy and then scale the new covariance matrix by the default scale factor, which will be then used to get the Cholesky Factor
     596             : 
     597         665 :                 do j = 1, nd
     598        1224 :                     do i = 1, j
     599         559 :                         CovMatUpperOld(i,j) = comv_CholDiagLower(i,j,0)   ! This will be used to recover the old covariance in case of update failure, and to compute the adaptation measure
     600         967 :                         comv_CholDiagLower(i,j,0) = CovMatUpperCurrent(i,j) * mc_defaultScaleFactorSq
     601             :                     end do
     602             :                 end do
     603             : 
     604             :             else blockMergeCovMat
     605             : 
     606             :                 ! first scale the new covariance matrix by the default scale factor, which will be then used to get the Cholesky Factor
     607             : 
     608       79016 :                 do j = 1, nd
     609      146532 :                     do i = 1, j
     610       67516 :                         CovMatUpperOld(i,j) = comv_CholDiagLower(i,j,0)   ! This will be used to recover the old covariance in case of update failure, and to compute the adaptation measure
     611      116360 :                         CovMatUpperCurrent(i,j) = CovMatUpperCurrent(i,j) * mc_defaultScaleFactorSq
     612             :                     end do
     613             :                 end do
     614             : 
     615             :                 ! now combine it with the old covariance matrix.
     616             :                 ! Do not set the full boundaries' range `(1:nd)` for `comv_CholDiagLower` in the following subroutine call.
     617             :                 ! Setting the boundaries forces the compiler to generate a temporary array.
     618             : 
     619             :                 call mergeMeanCovUpper  ( nd            = nd                            & ! LCOV_EXCL_LINE
     620             :                                         , npA           = mv_sampleSizeOld_save         & ! LCOV_EXCL_LINE
     621             :                                         , MeanVecA      = mv_MeanOld_save               & ! LCOV_EXCL_LINE
     622             :                                         , CovMatUpperA  = CovMatUpperOld                & ! LCOV_EXCL_LINE
     623             :                                         , npB           = sampleSizeCurrent             & ! LCOV_EXCL_LINE
     624             :                                         , MeanVecB      = MeanCurrent                   & ! LCOV_EXCL_LINE
     625             :                                         , CovMatUpperB  = CovMatUpperCurrent            & ! LCOV_EXCL_LINE
     626             :                                         , MeanVecAB     = MeanNew                       & ! LCOV_EXCL_LINE
     627             :                                         , CovMatUpperAB = comv_CholDiagLower(:,1:nd,0)  & ! LCOV_EXCL_LINE
     628       30172 :                                         )
     629       79016 :                 mv_MeanOld_save(1:nd) = MeanNew
     630             : 
     631             :             end if blockMergeCovMat
     632             : 
     633             :             !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     634             : 
     635             :             ! now get the Cholesky factorization
     636             : 
     637             :             ! WARNING: Do not set the full boundaries' range `(1:nd)` for the first index of `comv_CholDiagLower` in the following subroutine call.
     638             :             ! WARNING: Setting the boundaries forces the compiler to generate a temporary array.
     639             : 
     640       30429 :             call getCholeskyFactor( nd, comv_CholDiagLower(:,1:nd,0), comv_CholDiagLower(1:nd,0,0) )
     641             : 
     642       30429 :             blockPosDefCheck: if (comv_CholDiagLower(1,0,0)>0._RK) then
     643             : 
     644             :                 !singularityOccurred = .false.
     645       30429 :                 samplerUpdateSucceeded = .true.
     646       30429 :                 adaptationMeasureComputationNeeded = .true.
     647       30429 :                 mv_sampleSizeOld_save = mv_sampleSizeOld_save + sampleSizeCurrent
     648             : 
     649             :             ! LCOV_EXCL_START
     650             :             else blockPosDefCheck
     651             : 
     652             :                 adaptationMeasure = 0._RK
     653             :                 !singularityOccurred = .true.
     654             :                 samplerUpdateSucceeded = .false.
     655             :                 adaptationMeasureComputationNeeded = .false.
     656             :                 mv_sampleSizeOld_save = sampleSizeOld
     657             : 
     658             :                 ! it may be a good idea to add a warning message printed out here for the singularity occurrence
     659             : 
     660             :                 call warn   ( prefix = mc_methodBrand       & ! LCOV_EXCL_LINE
     661             :                             , outputUnit = mc_logFileUnit   & ! LCOV_EXCL_LINE
     662             :                             , marginTop = 0_IK              & ! LCOV_EXCL_LINE
     663             :                             , marginBot = 0_IK              & ! LCOV_EXCL_LINE
     664             :                             , msg = "Singularity occurred while updating the proposal distribution's covariance matrix." & ! LCOV_EXCL_LINE
     665             :                             )
     666             : 
     667             :                 ! recover the old upper covariance matrix
     668             : 
     669             :                 do j = 1, nd
     670             :                     do i = 1, j
     671             :                         comv_CholDiagLower(i,j,0) = CovMatUpperOld(i,j)
     672             :                     end do
     673             :                 end do
     674             : 
     675             :                 ! ensure the old Cholesky factorization can be recovered
     676             : 
     677             :                 ! WARNING: Do not set the full boundaries' range `(1:nd)` for the first index of `comv_CholDiagLower` in the following subroutine call.
     678             :                 ! WARNING: Setting the boundaries forces the compiler to generate a temporary array.
     679             : 
     680             :                 call getCholeskyFactor( nd, comv_CholDiagLower(:,1:nd,0), comv_CholDiagLower(1:nd,0,0) ) ! avoid temporary array creation by using : indexer
     681             :                 if (comv_CholDiagLower(1,0,0)<0._RK) then
     682             :                     write(mc_logFileUnit,"(A)")
     683             :                     write(mc_logFileUnit,"(A)") "Singular covariance matrix detected:"
     684             :                     write(mc_logFileUnit,"(A)")
     685             :                     do j = 1, nd
     686             :                         write(mc_logFileUnit,"(*(E25.15))") comv_CholDiagLower(1:j,j,0)
     687             :                     end do
     688             :                     write(mc_logFileUnit,"(A)")
     689             :                     ProposalErr%occurred = .true.
     690             :                     ProposalErr%msg = PROCEDURE_NAME // &
     691             :                                     ": Error occurred while attempting to compute the Cholesky factorization of the &
     692             :                                     &covariance matrix of the proposal distribution of " // mc_methodName // "'s sampler. &
     693             :                                     &This is highly unusual, and can be indicative of some major underlying problems.\n&
     694             :                                     &It may also be due to a runtime computational glitch, in particular, for high-dimensional simulations. &
     695             :                                     &In such case, consider increasing the value of the input variable adaptiveUpdatePeriod.\n&
     696             :                                     &It may also be that your input objective function has been incorrectly implemented.\n&
     697             :                                     &For example, ensure that you are passing a correct value of ndim to the ParaMonte sampler routine,\n&
     698             :                                     &the same value that is expected as input to your objective function's implementation.\n&
     699             :                                     &Otherwise, restarting the simulation might resolve the error."
     700             :                     call abort( Err = ProposalErr, prefix = mc_methodBrand, newline = "\n", outputUnit = mc_logFileUnit )
     701             :                     return
     702             :                 end if
     703             : 
     704             :             end if blockPosDefCheck
     705             :             ! LCOV_EXCL_STOP
     706             : 
     707             :             !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     708             : 
     709             :             !! perform global adaptive scaling is requested
     710             :             !if (mc_scalingRequested) then
     711             :             !    if (meanAccRateSinceStart<mc_targetAcceptanceRate) then
     712             :             !        mv_adaptiveScaleFactorSq_save = mc_maxScaleFactorSq**((mc_targetAcceptanceRate-meanAccRateSinceStart)/mc_targetAcceptanceRate)
     713             :             !    else
     714             :             !        mv_adaptiveScaleFactorSq_save = mc_maxScaleFactorSq**((mc_targetAcceptanceRate-meanAccRateSinceStart)/(1._RK-mc_targetAcceptanceRate))
     715             :             !    end if
     716             :             !    mv_adaptiveScaleFactorSq_save = mv_adaptiveScaleFactorSq_save * (meanAccRateSinceStart/mc_targetAcceptanceRate)**mc_ndimInverse
     717             :             !end if
     718       30429 :             if (mc_scalingRequested) scalingNeeded = .true.
     719             : 
     720             :         else blockSufficientSampleSizeCheck
     721             : 
     722             :             ! singularity has occurred. If the first covariance merging has not occurred yet, set the scaling factor appropriately to shrink the covariance matrix.
     723             : 
     724       91678 :             samplerUpdateSucceeded = .false.
     725       91678 :             if (mv_sampleSizeOld_save==1_IK .or. mc_scalingRequested) then
     726       82834 :                 scalingNeeded = .true.
     727       82834 :                 adaptationMeasureComputationNeeded = .true.
     728             :                 ! save the old covariance matrix for the computation of the adaptation measure
     729      248446 :                 do j = 1, nd
     730      496836 :                     do i = 1,j
     731      414002 :                         CovMatUpperOld(i,j) = comv_CholDiagLower(i,j,0)
     732             :                     end do
     733             :                 end do
     734             :             else
     735        8844 :                 adaptationMeasure = 0._RK
     736        8844 :                 adaptationMeasureComputationNeeded = .false.
     737             :             end if
     738             : 
     739             : 
     740             :         end if blockSufficientSampleSizeCheck
     741             : 
     742             :         ! adjust the scale of the covariance matrix and the Cholesky factor, if needed
     743             : 
     744             : !scalingNeeded = .true.
     745      122107 :         if (scalingNeeded) then
     746       91951 :             if ( meanAccRateSinceStart < mc_TargetAcceptanceRateLimit(1) .or. meanAccRateSinceStart > mc_TargetAcceptanceRateLimit(2) ) then
     747        9851 :                 mv_adaptiveScaleFactorSq_save = (meanAccRateSinceStart/mc_targetAcceptanceRate)**mc_ndimInverse
     748             : !block
     749             : !    use Statistics_mod, only: getRandUniform
     750             : !    integer, save :: counter = 0_IK
     751             : !    counter = counter - 1
     752             : !    mv_adaptiveScaleFactorSq_save = mv_adaptiveScaleFactorSq_save * exp(-counter*getRandUniform(-1.e0_RK,1.e0_RK)/1.e4_RK)
     753             : !    !use Statistics_mod, only: getRandInt
     754             : !    !mv_adaptiveScaleFactorSq_save = mv_adaptiveScaleFactorSq_save * exp(real(getRandInt(-1_IK,1_IK),kind=RK))
     755             : !    !write(*,*) counter, mv_adaptiveScaleFactorSq_save
     756             : !end block
     757        9851 :                 adaptiveScaleFactor = sqrt(mv_adaptiveScaleFactorSq_save)
     758       29534 :                 do j = 1, nd
     759             :                     ! update the Cholesky diagonal elements
     760       19683 :                     comv_CholDiagLower(j,0,0) = comv_CholDiagLower(j,0,0) * adaptiveScaleFactor
     761             :                     ! update covariance matrix
     762       49198 :                     do i = 1,j
     763       49198 :                         comv_CholDiagLower(i,j,0) = comv_CholDiagLower(i,j,0) * mv_adaptiveScaleFactorSq_save
     764             :                     end do
     765             :                     ! update the Cholesky factorization
     766       39366 :                     do i = j+1, nd
     767       29515 :                         comv_CholDiagLower(i,j,0) = comv_CholDiagLower(i,j,0) * adaptiveScaleFactor
     768             :                     end do
     769             :                 end do
     770             :             end if
     771             :         end if
     772             : 
     773             :         ! compute the adaptivity only if any updates has occurred
     774             : 
     775             : 
     776      122107 :         blockAdaptationMeasureComputation: if (adaptationMeasureComputationNeeded) then
     777             : 
     778      328127 :             logSqrtDetNew = sum(log( comv_CholDiagLower(1:nd,0,0) ))
     779             : 
     780             :             ! use the universal upper bound
     781             : 
     782      328127 :             do j = 1, nd
     783      644592 :                 do i = 1,j
     784      531329 :                     CovMatUpperCurrent(i,j) = 0.5_RK * ( comv_CholDiagLower(i,j,0) + CovMatUpperOld(i,j) ) ! dummy
     785             :                 end do
     786             :             end do
     787      113263 :             call getLogSqrtDetPosDefMat(nd,CovMatUpperCurrent,logSqrtDetSum,singularityOccurred)
     788      113263 :             if (singularityOccurred) then
     789             :                 ! LCOV_EXCL_START
     790             :                 write(mc_logFileUnit,"(A)")
     791             :                 write(mc_logFileUnit,"(A)") "Singular covariance matrix detected while computing the Adaptation measure:"
     792             :                 write(mc_logFileUnit,"(A)")
     793             :                 do j = 1, nd
     794             :                     write(mc_logFileUnit,"(*(E25.15))") CovMatUpperCurrent(1:j,j)
     795             :                 end do
     796             :                 ProposalErr%occurred = .true.
     797             :                 ProposalErr%msg = PROCEDURE_NAME // &
     798             :                                 ": Error occurred while computing the Cholesky factorization of &
     799             :                                 &a matrix needed for the computation of the Adaptation measure. &
     800             :                                 &Such error is highly unusual, and requires an in depth investigation of the case.\n&
     801             :                                 &It may also be due to a runtime computational glitch, in particular, for high-dimensional simulations. &
     802             :                                 &In such case, consider increasing the value of the input variable adaptiveUpdatePeriod.\n&
     803             :                                 &It may also be that your input objective function has been incorrectly implemented.\n&
     804             :                                 &For example, ensure that you are passing a correct value of ndim to the ParaMonte sampler routine,\n&
     805             :                                 &the same value that is expected as input to your objective function's implementation.\n&
     806             :                                 &Otherwise, restarting the simulation might resolve the error."
     807             :                 call abort( Err = ProposalErr, prefix = mc_methodBrand, newline = "\n", outputUnit = mc_logFileUnit )
     808             :                 return
     809             :                 ! LCOV_EXCL_STOP
     810             :             end if
     811             : 
     812             :             !adaptationMeasure = 1._RK - exp( 0.5_RK*(mv_logSqrtDetOld_save+logSqrtDetNew) - logSqrtDetSum )
     813      113263 :             adaptationMeasure = 1._RK - exp( 0.5*(mv_logSqrtDetOld_save + logSqrtDetNew) - logSqrtDetSum ) ! totalVariationUpperBound
     814      113263 :             if (adaptationMeasure>0._RK) then
     815       40431 :                 adaptationMeasure = sqrt(adaptationMeasure) ! totalVariationUpperBound
     816             :             ! LCOV_EXCL_START
     817             :             elseif (adaptationMeasure<0._RK) then
     818             :                 call warn   ( prefix = mc_methodBrand &
     819             :                             , outputUnit = mc_logFileUnit &
     820             :                             , msg = mc_negativeTotalVariationMsg//num2str(adaptationMeasure) )
     821             :                 adaptationMeasure = 0._RK
     822             :             end if
     823             :             ! LCOV_EXCL_STOP
     824      113263 :             mv_logSqrtDetOld_save = logSqrtDetNew
     825             : 
     826             : !block
     827             : !integer, save :: counter = 0
     828             : !counter = counter + 1
     829             : !!if (counter==1) then
     830             : !if (adaptationMeasure>1._RK) then
     831             : !write(*,*)
     832             : !write(*,*) mv_logSqrtDetOld_save
     833             : !write(*,*) logSqrtDetNew
     834             : !write(*,*) logSqrtDetSum
     835             : !write(*,*) mv_logSqrtDetOld_save + logSqrtDetNew - 2_IK * logSqrtDetSum
     836             : !write(*,*) exp( mv_logSqrtDetOld_save + logSqrtDetNew - 2_IK * logSqrtDetSum )
     837             : !write(*,*)
     838             : !end if
     839             : !end block
     840             : 
     841             :             ! update the higher-stage delayed-rejection Cholesky Lower matrices
     842             : 
     843      113263 :             if (mc_delayedRejectionRequested) call updateDelRejCholDiagLower()
     844      113263 :             call getInvCovMat()
     845             : 
     846             :         end if blockAdaptationMeasureComputation
     847             : 
     848             :         ! read/write restart file
     849             : 
     850      122107 :         if (mc_Image%isLeader) then
     851      122107 :             if (isFreshRun) then
     852       74420 :                 call writeRestartFile()
     853             :             else
     854       47687 :                 call readRestartFile()
     855             :             end if
     856             :         end if
     857             : 
     858      122107 :     end subroutine doAdaptation
     859             : 
     860             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     861             : 
     862             :     ! LCOV_EXCL_START
     863             :     ! ATTN: This routine needs further correction for the delayed rejection method
     864             : #if false
     865             :     subroutine doAutoTune   ( adaptationMeasure &
     866             :                             , AutoTuneScaleSq   &
     867             :                             )
     868             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     869             :         !DEC$ ATTRIBUTES DLLEXPORT :: doAutoTune
     870             : #endif
     871             :         use Matrix_mod, only: getLogSqrtDetPosDefMat
     872             :         use Constants_mod, only: RK, IK
     873             :         use Err_mod, only: abort
     874             :         implicit none
     875             :         character(*), parameter                         :: PROCEDURE_NAME = MODULE_NAME // "@doAutoTune()"
     876             : 
     877             :         real(RK)   , intent(in)                         :: AutoTuneScaleSq(1)
     878             :         real(RK)   , intent(inout)                      :: adaptationMeasure
     879             :         real(RK)                                        :: logSqrtDetSum, logSqrtDetOld, logSqrtDetNew
     880             :         real(RK)                                        :: CovMatUpperOld(1,1), CovMatUpperCurrent(1,1)
     881             :         logical                                         :: singularityOccurred
     882             : 
     883             :         CovMatUpperOld = comv_CholDiagLower(1:mc_ndim,1:mc_ndim,0)
     884             :         logSqrtDetOld = sum(log( comv_CholDiagLower(1:mc_ndim,0,0) ))
     885             : 
     886             :         if (AutoTuneScaleSq(1)==0._RK) then
     887             :             comv_CholDiagLower(1,1,0) = 0.25_RK*comv_CholDiagLower(1,1,0)
     888             :             comv_CholDiagLower(1,0,0) = sqrt(comv_CholDiagLower(1,1,0))
     889             :         else
     890             :             comv_CholDiagLower(1,1,0) = AutoTuneScaleSq(1)
     891             :             comv_CholDiagLower(1,0,0) = sqrt(AutoTuneScaleSq(1))
     892             :         end if
     893             : 
     894             :         ! compute the adaptivity
     895             : 
     896             :         logSqrtDetNew = sum(log( comv_CholDiagLower(1:mc_ndim,0,0) ))
     897             :         CovMatUpperCurrent = 0.5_RK * ( comv_CholDiagLower(1:mc_ndim,1:mc_ndim,0) + CovMatUpperOld )
     898             :         call getLogSqrtDetPosDefMat(1_IK,CovMatUpperCurrent,logSqrtDetSum,singularityOccurred)
     899             :         if (singularityOccurred) then
     900             :             ProposalErr%occurred = .true.
     901             :             ProposalErr%msg = PROCEDURE_NAME // &
     902             :                             ": Error occurred while computing the Cholesky factorization of &
     903             :                             &a matrix needed for the computation of the proposal distribution's adaptation measure. &
     904             :                             &Such error is highly unusual, and requires an in depth investigation of the case. &
     905             :                             &It may also be that your input objective function has been incorrectly implemented.\n&
     906             :                             &For example, ensure that you are passing a correct value of ndim to the ParaMonte sampler routine,\n&
     907             :                             &the same value that is expected as input to your objective function's implementation.\n&
     908             :                             &Otherwise, restarting the simulation might resolve the error."
     909             :             call abort( Err = ProposalErr, prefix = mc_methodBrand, newline = "\n", outputUnit = mc_logFileUnit )
     910             :             return
     911             :         end if
     912             :         adaptationMeasure = 1._RK - exp( 0.5_RK*(logSqrtDetOld+logSqrtDetNew) - logSqrtDetSum )
     913             : 
     914             :     end subroutine doAutoTune
     915             :     ! LCOV_EXCL_STOP
     916             : #endif
     917             : 
     918             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     919             : 
     920             :     ! Note: based on some benchmarks with ndim = 1, the new design with merging cholesky diag and lower is faster than the original
     921             :     ! Note: double-communication, implementation. Here are some timings on 4 images:
     922             :     ! Note: new single communication:
     923             :     ! Note: image 2: avgTime =  6.960734060198531E-006
     924             :     ! Note: image 3: avgTime =  7.658279491640721E-006
     925             :     ! Note: image 4: avgTime =  9.261191328273417E-006
     926             :     ! Note: avg(avgTime): 7.960068293370891e-06
     927             :     ! Note: old double communication:
     928             :     ! Note: image 4: avgTime =  1.733505615104837E-005
     929             :     ! Note: image 3: avgTime =  1.442608268140151E-005
     930             :     ! Note: image 2: avgTime =  1.420345299036357E-005
     931             :     ! Note: avg(avgTime): 1.532153060760448e-05
     932             :     ! Note: avg(speedup): 1.924798889020109
     933             :     ! Note: One would expect this speed up to diminish as ndim goes to infinity,
     934             :     ! Note: since data transfer will dominate the communication overhead.
     935             :     !> \brief
     936             :     !> Broadcast adaptation to all images.
     937             :     !> \warning
     938             :     !> When CAF parallelism is used, this routine must be exclusively called by the rooter images.
     939             :     !> When MPI parallelism is used, this routine must be called by all images.
     940             : #if defined CAF_ENABLED || MPI_ENABLED
     941       99555 :     subroutine bcastAdaptation()
     942             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     943             :         !DEC$ ATTRIBUTES DLLEXPORT :: bcastAdaptation
     944             : #endif
     945             : #if defined CAF_ENABLED
     946             :         implicit none
     947             :         comv_CholDiagLower(1:mc_ndim,0:mc_ndim,0) = comv_CholDiagLower(1:mc_ndim,0:mc_ndim,0)[1]
     948             :         if (mc_delayedRejectionRequested) call updateDelRejCholDiagLower()  ! update the higher-stage delayed-rejection Cholesky Lower matrices
     949             : #elif defined MPI_ENABLED
     950             :         use mpi ! LCOV_EXCL_LINE
     951             :         implicit none
     952             :         integer :: ierrMPI
     953             :         call mpi_bcast  ( comv_CholDiagLower    & ! LCOV_EXCL_LINE ! buffer: XXX: first element need not be shared. This may need a fix in future.
     954             :                         , mc_ndimSqPlusNdim     & ! LCOV_EXCL_LINE ! count
     955             :                         , mpi_double_precision  & ! LCOV_EXCL_LINE ! datatype
     956             :                         , 0                     & ! LCOV_EXCL_LINE ! root: broadcasting rank
     957             :                         , mpi_comm_world        & ! LCOV_EXCL_LINE ! comm
     958             :                         , ierrMPI               & ! LCOV_EXCL_LINE ! ierr
     959       99555 :                         )
     960             :         ! It is essential for the following to be exclusively done by the rooter images. The leaders have had their updates in `doAdaptation()`.
     961       99555 :         if (mc_Image%isRooter .and. mc_delayedRejectionRequested) call updateDelRejCholDiagLower()
     962             : #endif
     963       99555 :         call getInvCovMat()
     964       99555 :     end subroutine bcastAdaptation
     965             : #endif
     966             : 
     967             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     968             : 
     969             :     !> \brief
     970             :     !> Return the inverse covariance matrix of the current covariance of the proposal distribution.
     971             :     !>
     972             :     !> \warning
     973             :     !> This procedure is `impure` and sets the values of multiple global variables upon exit.
     974             :     !>
     975             :     !> \remark
     976             :     !> This procedure is required for the proposal probabilities
     977             :     !>
     978             :     !> \todo
     979             :     !> The specification of the array bounds causes the compiler to generate a temporary copy of the array, despite being unnecessary.
     980             :     !> Right now, this is resolved by replacing the array bounds with `:`. A better solution is to add the `contiguous` attribute to
     981             :     !> the corresponding argument of [getInvMatFromCholFac](ref matrix_mod::getinvmatfromcholfac) to guarantee it to the compiler.
     982             :     !> More than improving performance, this would turn off the pesky compiler warnings about temporary array creation.
     983      213505 :     subroutine getInvCovMat()
     984             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     985             :         !DEC$ ATTRIBUTES DLLEXPORT :: getInvCovMat
     986             : #endif
     987             :         use Matrix_mod, only: getInvMatFromCholFac ! LCOV_EXCL_LINE
     988             :         implicit none
     989             :         integer(IK) :: istage
     990             :         ! update the inverse covariance matrix of the proposal from the computed Cholesky factor
     991      213505 :         do concurrent(istage=0:mc_DelayedRejectionCount)
     992             :             ! WARNING: Do not set the full boundaries' range `(1:mc_ndim)` for the first index of `comv_CholDiagLower` in the following subroutine call.
     993             :             ! WARNING: Setting the boundaries forces the compiler to generate a temporary array.
     994             :             mv_InvCovMat(1:mc_ndim,1:mc_ndim,istage) = getInvMatFromCholFac ( nd = mc_ndim & ! LCOV_EXCL_LINE
     995             :                                                                             , CholeskyLower = comv_CholDiagLower(:,1:mc_ndim,istage) & ! LCOV_EXCL_LINE
     996             :                                                                             , Diagonal = comv_CholDiagLower(1:mc_ndim,0,istage) & ! LCOV_EXCL_LINE
     997     2286752 :                                                                             )
     998     1233758 :             mv_logSqrtDetInvCovMat(istage) = -sum(log( comv_CholDiagLower(1:mc_ndim,0,istage) ))
     999             :         end do
    1000             : !if (mc_ndim==1 .and. abs(log(sqrt(mv_InvCovMat(1,1,0)))-mv_logSqrtDetInvCovMat(0))>1.e-13_RK) then
    1001             : !write(*,"(*(g0,:,' '))") "log(sqrt(mv_InvCovMat(1,1,0))) /= mv_logSqrtDetInvCovMat(0)"
    1002             : !write(*,"(*(g0,:,' '))") log(sqrt(mv_InvCovMat(1,1,0))), mv_logSqrtDetInvCovMat(0)
    1003             : !write(*,"(*(g0,:,' '))") abs(-log(sqrt(mv_InvCovMat(1,1,0)))-mv_logSqrtDetInvCovMat(0))
    1004             : !error stop
    1005             : !endif
    1006      213505 :     end subroutine getInvCovMat
    1007             : 
    1008             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1009             : 
    1010             :     !> \brief
    1011             :     !> Update the Cholesky factorization of the higher-stage delayed rejection covariance matrices.
    1012             :     !>
    1013             :     !> \warning
    1014             :     !> This procedure is impure and sets the values of multiple global variables upon exit.
    1015             :     !>
    1016             :     !> \todo
    1017             :     !> The performance of this update could be improved by only updating the higher-stage covariance, only when needed.
    1018             :     !> However, the gain will be likely minimal, especially in low-dimensions.
    1019       30197 :     subroutine updateDelRejCholDiagLower()
    1020             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1021             :         !DEC$ ATTRIBUTES DLLEXPORT :: updateDelRejCholDiagLower
    1022             : #endif
    1023             :         implicit none
    1024             :         integer(IK) :: j, istage
    1025             :         ! update the Cholesky factor of the delayed-rejection-stage proposal distributions
    1026      147623 :         do istage = 1, mc_DelayedRejectionCount
    1027      346413 :             comv_CholDiagLower(1:mc_ndim,0,istage) = comv_CholDiagLower(1:mc_ndim,0,istage-1) * mc_DelayedRejectionScaleFactorVec(istage)
    1028      376610 :             do j = 1, mc_ndim
    1029      457974 :                 comv_CholDiagLower(j+1:mc_ndim,j,istage) = comv_CholDiagLower(j+1:mc_ndim,j,istage-1) * mc_DelayedRejectionScaleFactorVec(istage)
    1030             :             end do
    1031             :         end do
    1032             :         ! There is no need to check for positive-definiteness of the comv_CholDiagLower, it is already checked on the first image.
    1033      243702 :     end subroutine updateDelRejCholDiagLower
    1034             : 
    1035             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1036             : 
    1037             :     !> \brief
    1038             :     !> This procedure is a static method of the [ParaDXXXProposalNormal_type](@ref paradxxxproposalnormal_type)
    1039             :     !> or [ParaDXXXProposalUniform_type](@ref paradxxxproposaluniform_type) classes.\n
    1040             :     !> This procedure is called by the sampler kernel routines.\n
    1041             :     !> Write the restart information to the output file.
    1042             :     !>
    1043             :     !> @param[in]   meanAccRateSinceStart : The current mean acceptance rate of the sampling (**optional**).
    1044      149314 :     subroutine writeRestartFile(meanAccRateSinceStart)
    1045             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1046             :         !DEC$ ATTRIBUTES DLLEXPORT :: writeRestartFile
    1047             : #endif
    1048             :         implicit none
    1049             :         real(RK), intent(in), optional  :: meanAccRateSinceStart
    1050             :         integer(IK)                     :: i, j
    1051      149314 :         if (present(meanAccRateSinceStart)) then
    1052       74657 :             if (mc_isBinaryRestartFileFormat) then
    1053       73302 :                 write(mc_restartFileUnit) meanAccRateSinceStart
    1054             :             else
    1055        1355 :                 write(mc_restartFileUnit,mc_restartFileFormat) "meanAcceptanceRateSinceStart", meanAccRateSinceStart
    1056             :             end if
    1057       74657 :         elseif (mc_isAsciiRestartFileFormat) then
    1058        1355 :             write( mc_restartFileUnit, mc_restartFileFormat ) "sampleSize" & ! sampleSizeOld
    1059             :                                                             , mv_sampleSizeOld_save &
    1060        1355 :                                                             , "logSqrtDeterminant" & ! logSqrtDetOld
    1061             :                                                             , mv_logSqrtDetOld_save &
    1062        1355 :                                                             , "adaptiveScaleFactorSquared" & ! adaptiveScaleFactorSq
    1063             :                                                             , mv_adaptiveScaleFactorSq_save * mc_defaultScaleFactorSq &
    1064        1355 :                                                             , "meanVec" & ! MeanOld(1:ndim)
    1065             :                                                             , mv_MeanOld_save(1:mc_ndim) & ! LCOV_EXCL_LINE
    1066        1355 :                                                             , "covMat" & ! CholDiagLower(1:ndim,0:ndim,0)
    1067        8696 :                                                             , ((comv_CholDiagLower(i,j,0),i=1,j),j=1,mc_ndim)
    1068             :                                                            !, (comv_CholDiagLower(1:mc_ndim,0:mc_ndim,0)
    1069             :         end if
    1070      149314 :         flush(mc_restartFileUnit)
    1071      179511 :     end subroutine writeRestartFile
    1072             : 
    1073             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1074             : 
    1075             :     !> \brief
    1076             :     !> This procedure is a static method of the [ParaDXXXProposalNormal_type](@ref paradxxxproposalnormal_type)
    1077             :     !> or [ParaDXXXProposalUniform_type](@ref paradxxxproposaluniform_type) classes.\n
    1078             :     !> This procedure is called by the sampler kernel routines.\n
    1079             :     !> Read the restart information from the restart file.
    1080             :     !>
    1081             :     !> @param[out]  meanAccRateSinceStart : The current mean acceptance rate of the sampling (**optional**).
    1082       95438 :     subroutine readRestartFile(meanAccRateSinceStart)
    1083             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1084             :         !DEC$ ATTRIBUTES DLLEXPORT :: readRestartFile
    1085             : #endif
    1086             :         implicit none
    1087             :         real(RK), intent(out), optional :: meanAccRateSinceStart
    1088             :         integer(IK)                     :: i
    1089       95438 :         if (present(meanAccRateSinceStart)) then
    1090       47719 :             if (mc_isBinaryRestartFileFormat) then
    1091       47173 :                 read(mc_restartFileUnit) meanAccRateSinceStart
    1092             :             else
    1093         546 :                 read(mc_restartFileUnit,*)
    1094         546 :                 read(mc_restartFileUnit,*) meanAccRateSinceStart
    1095             :             end if
    1096       47719 :         elseif (mc_isAsciiRestartFileFormat) then
    1097        7644 :             do i = 1, 8 + mc_ndim * (mc_ndim+3) / 2
    1098             :                 !read( mc_restartFileUnit, mc_restartFileFormat )
    1099        7644 :                 read( mc_restartFileUnit, * )
    1100             :             end do
    1101             :         end if
    1102      244752 :     end subroutine readRestartFile
    1103             : 
    1104             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1105             : 
    1106             : 
    1107             : #undef ParaDXXX
    1108             : #undef GET_RANDOM_PROPOSAL
    1109             : 

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