Sampling a Multivariate Normal (MVN) distribution | ParaMonte: Parallel Monte Carlo Library

The mathematical objective function

Suppose we want to sample random vectors

\[\mathbf{x}=[x_1,x_2,x_3,x_4] ~~,\]

from a 4-dimensional Multivariate Normal (MVN) Probability Density Function (PDF),

\[\mathrm{PDF}_\mathrm{MVN} \equiv \pi ( \mathbf{x} ~|~ \mu, \Sigma) = (2\pi)^{-\frac{k}{2}} \sqrt{|\Sigma|} \exp\bigg( -\frac{1}{2} (\mathbf{x}-\mu)^\mathrm{T}\Sigma^{-1}(\mathbf{x}-\mu) \bigg) ~~,\]

that has a mean vector,

\[\mu = [0,0,0,0]\]

with a covariance matrix,

\[\Sigma = \begin{bmatrix} 1.0 & 0.5 & 0.5 & 0.5 \\ 0.5 & 1.0 & 0.5 & 0.5 \\ 0.5 & 0.5 & 1.0 & 0.5 \\ 0.5 & 0.5 & 0.5 & 1.0 \end{bmatrix}\]

This will serve as our mathematical objective function. Clearly, all of the four dimensions of this MVN are correlated with each other with a Pearson’s correlation coefficient of $0.5$.

Working with the natural logarithm of the objective function

In general, since probabilities are often tiny numbers, we’d want to always work with their natural logarithms in the world of computers to avoid unwanted overflow or underflow of numbers, which can crash our simulations. To highlight the importance of taking the natural logarithms of all input objective functions to the ParaMonte routines, we always name them like getLogFunc().

Important: All ParaMonte routines assume that the user-provided objective functions return the natural logarithms of the actual values of the objective functions.

Implementing the objective function in Fortran

Here is a minimalistic implementation of the 4-D multivariate normal distribution objective function in Fortran, which can also be downloaded as logfunc.f90.

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%  
!%  Description:
!%      +   Return the natural logarithm of an ndim-dimensional Multivariate Normal (MVN) 
!%          probability density function (PDF) with the Mean and Covariance Matrix as defined below.
!%          Reference: https://en.wikipedia.org/wiki/Multivariate_normal_distribution
!%  Input:
!%      +   ndim:       The number of dimensions of the domain of the objective function.
!%      +   point:      The input 64-bit real-valued vector of length ndim, 
!%                      at which the natural logarithm of objective function is computed.
!%  Output:
!%      +   logFunc:    A 64-bit real scalar number representing the natural logarithm of the objective function.
!%  Author:
!%      +   Computational Data Science Lab, Monday 9:03 AM, May 16 2016, ICES, UT Austin
!%  Visit:
!%      +   https://www.cdslab.org/paramonte
!%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

module LogFunc_mod

    use paramonte, only: IK, RK

    implicit none

    ! The number of dimensions of the domain of the objective function

    integer(IK), parameter  :: NDIM = 4_IK

    ! The mean vector of the MVN

    real(RK), parameter     :: MEAN(NDIM) = [0._RK, 0._RK, 0._RK, 0._RK]

    ! The coefficient of the Multivariate Normal Distribution: log(1/sqrt(2*Pi)^ndim)

    real(RK), parameter     :: LOG_INVERSE_SQRT_TWO_PI = log(1._RK/sqrt(2._RK*acos(-1._RK)))

    ! The covariance matrix of the MVN

    real(RK), parameter     :: COVMAT(NDIM,NDIM) =  reshape([ 1.0_RK,0.5_RK,0.5_RK,0.5_RK &
                                                            , 0.5_RK,1.0_RK,0.5_RK,0.5_RK &
                                                            , 0.5_RK,0.5_RK,1.0_RK,0.5_RK &
                                                            , 0.5_RK,0.5_RK,0.5_RK,1.0_RK ], shape=shape(COVMAT))

    ! The inverse covariance matrix of the MVN

    real(RK), parameter     :: INVCOVMAT(NDIM,NDIM) = reshape(  [ +1.6_RK, -0.4_RK, -0.4_RK, -0.4_RK &
                                                                , -0.4_RK, +1.6_RK, -0.4_RK, -0.4_RK &
                                                                , -0.4_RK, -0.4_RK, +1.6_RK, -0.4_RK &
                                                                , -0.4_RK, -0.4_RK, -0.4_RK, +1.6_RK ], shape=shape(INVCOVMAT))

    ! The logarithm of square root of the determinant of the inverse covariance matrix of the Multivariate Normal Distribution 

    real(RK), parameter     :: LOG_SQRT_DET_INV_COV = 0.581575404902840_RK

    ! MVN_COEF

    real(RK), parameter     :: MVN_COEF = NDIM * LOG_INVERSE_SQRT_TWO_PI + LOG_SQRT_DET_INV_COV

contains

    function getLogFunc(ndim,Point) result(logFunc)
        ! Return the negative natural logarithm of MVN distribution evaluated at the input vector point.
        implicit none
        integer(IK), intent(in) :: ndim
        real(RK), intent(in)    :: Point(ndim)
        real(RK)                :: NormedPoint(ndim)
        real(RK)                :: logFunc
        NormedPoint = Point - MEAN
        logFunc = MVN_COEF - 0.5_RK * ( dot_product(NormedPoint,matmul(INVCOVMAT,NormedPoint)) )
    end function getLogFunc

end module LogFunc_mod

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Calling the ParaMonte samplers

Suppose we want to randomly sample the above objective function via the Delayed-Rejection Adaptive Metropolis-Hastings Markov Chain Monte Carlo sampler of ParaMonte (the ParaDRAM sampler).

Depending on whether you have already built the ParaMonte library from scratch on your system or you intend to use the prebuilt ParaMonte libraries, as well as the compiler you are using to build your application, there are two interfaces to the ParaDRAM sampler. These two interfaces are separated via the compiler preprocessor flag IS_COMPATIBLE_COMPILER in the main.f90 code below,

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%  
!%  Description:
!%      +   Run the Monte Carlo sampler of the ParaMonte library given the input log-target density function `getLogFunc()`.
!%  Output:
!%      +   The simulation output files.
!%  Author:
!%      +   Computational Data Science Lab, Monday 9:03 AM, May 16 2016, ICES, UT Austin
!%  Visit:
!%      +   https://www.cdslab.org/paramonte
!%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

!   Keep in mind that Fortran is case-insensitive, except for character values and string values. 
!   So, feel free to use upper-case or lower-case for the Fortran syntax and entity names.
!   The ParaMonte library uses camelCase convention for variable naming.

#if defined IS_COMPATIBLE_COMPILER

    program main

    ! This is the Object-Oriented-Programming (OOP) style interface to the ParaMonte routines.
    ! This is more flexible but less portable, as its compilation requires the same compiler 
    ! brand and version with which the ParaMonte library has been built.

    use paramonte, only: ParaDRAM, IK
    use LogFunc_mod, only: getLogFunc, NDIM

    implicit none
    type(ParaDRAM) :: pmpd

    call pmpd%runSampler( ndim = NDIM &
                        , getLogFunc = getLogFunc &
                        , inputFile = "./paramonte.in" &    ! this is optional argument
                        ! You can also specify simulation specifications as input arguments, like 
                        ! the following. This is possible only from the OOP interface to ParaDRAM.
                        , greedyAdaptationCount = 0_IK &    ! this is optional argument
                        , description = "an example run" &  ! this is optional argument
                        ! More optional arguments can appear here. 
                        ! See the ParaDRAM routine's list of input arguments.
                        )

    end program main

#else

    ! This is the default simple procedural interface to the ParaMonte routines.
    ! The first two arguments to the sampler routines are mandatory.
    ! The inputFile argument is optional.

    program main

    use paramonte, only: runParaDRAM
    use LogFunc_mod, only: getLogFunc, NDIM

    implicit none

    call runParaDRAM( NDIM &
                    , getLogFunc &
                    , "./paramonte.in" & ! this is optional argument
                    )

    end program main

#endif

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Therefore, if you intend to use the Object-Oriented Programming (OOP) calling interface, you will have to define and pass the preprocessor flag IS_COMPATIBLE_COMPILER at compile-time to your compiler (to be discussed later). If this flag is not defined at compile-time, the compiler will automatically use the latter procedural-style calling method enclosed between #else and #endif Fortran preprocessor directives in the main code above.

Calling the ParaMonte samplers: The signature module file

Here is the contents of the Fortran paramonte.90 signature module file for the ParaMonte library routines,

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!
!!!!   MIT License
!!!!
!!!!   ParaMonte: plain powerful parallel Monte Carlo library.
!!!!
!!!!   Copyright (C) 2012-present, The Computational Data Science Lab
!!!!
!!!!   This file is part of the ParaMonte library.
!!!!
!!!!   Permission is hereby granted, free of charge, to any person obtaining a 
!!!!   copy of this software and associated documentation files (the "Software"), 
!!!!   to deal in the Software without restriction, including without limitation 
!!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense, 
!!!!   and/or sell copies of the Software, and to permit persons to whom the 
!!!!   Software is furnished to do so, subject to the following conditions:
!!!!
!!!!   The above copyright notice and this permission notice shall be 
!!!!   included in all copies or substantial portions of the Software.
!!!!
!!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 
!!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 
!!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 
!!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, 
!!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR 
!!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE 
!!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
!!!!
!!!!   ACKNOWLEDGMENT
!!!!
!!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
!!!!   As per the ParaMonte library license agreement terms, if you use any parts of 
!!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your 
!!!!   work (education/research/industry/development/...) by citing the ParaMonte 
!!!!   library as described on this page:
!!!!
!!!!       https://github.com/cdslaborg/paramonte/blob/3548c097f2a25dfc0613061800656d27d0e0ddbe/ACKNOWLEDGMENT.md
!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  
!  Description:
!       - The Fortran module source file containing the prototypes of the ParaMonte library routines to be called from Fortran.
!  Prototypes:
!       - runParaDRAM   : The procedural interface to the ParaDRAM sampler routine. This is the basic default interface.
!       - ParaDRAM_type : The ParaDRAM sampler class (derived type).
!       - ParaDRAM      : An alias for ParaDRAM_type().
!  Author:
!       - Computational Data Science Lab, Monday 9:03 AM, May 16 2016, ICES, UT Austin
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module paramonte

    ! ParaMonte default integer, real, and complex kinds are defined by IK, RK, and CK below.
    ! Use IK and RK from this module for integer and real kind specifications in your Fortran codes
    ! to ensure the highest level of consistency with the integer and real kinds of the ParaMonte library.

#if defined IS_COMPATIBLE_COMPILER

    ! Use the identifiers here to access the Object-Oriented interface to the ParaMonte routines.
    ! To enable the object-oriented interface, you must pass the IS_COMPATIBLE_COMPILER 
    ! preprocessor flag to the compiler at the time of compiling this module.

    use ParaDRAM_mod, only: IK, RK, CK
    use ParaDRAM_mod, only: getLogFunc_proc
    use ParaDRAM_mod, only: getLogFunc_interface => getLogFunc_proc
    use ParaDRAM_mod, only: ParaDRAM => ParaDRAM_type
    use ParaDRAM_mod, only: ParaDRAM_type

#else

    ! This is the procedural interface to the ParaDRAM sampler routine of the ParaMonte library.
    ! By default, if the preprocessor flag IS_COMPATIBLE_COMPILER is not passed to the compiler,
    ! the following procedural interfaces will be used, depending on the choice of Fortran compiler.

    use, intrinsic :: iso_fortran_env, only: IK => int32, RK => real64
    implicit none
    public :: runParaDRAM, getLogFunc_proc

    ! The Fortran objective function interface (getLogFunc). Here, `proc` stands for the procedure interface.

    abstract interface
        function getLogFunc_proc(ndim,Point) result(logFunc)
            import :: IK, RK
            integer(IK) , intent(in)    :: ndim
            real(RK)    , intent(in)    :: Point(ndim)
            real(RK)                    :: logFunc
        end function getLogFunc_proc
    end interface

#if defined __WIN64__ && __GFORTRAN__

    ! NOTE: __WIN64__ is not automatically predefined by GNU Fortran compiler.
    ! NOTE: The onus is on the user to define __WIN64__ at the time of compiling this code.

    ! This is the procedural interface to the ParaDRAM sampler routine of the ParaMonte library, 
    ! to be used by the GNU Fortran compiler in combination with the ParaMonte library prebuilt 
    ! via the Intel compiler suite on Windows. This particular separate interface exists here 
    ! because of the fundamental symbol exporting convention differences of the GNU and Intel 
    ! Fortran compilers, with Intel making all symbols uppercase, while GNU making all 
    ! symbols lowercase. 

    abstract interface
        ! C-style Fortran interface for the the objective function
        ! This is to be used only to bind the ParaMonte library compiled by the 
        ! Intel Compilers with Fortran applications compiled with GNU compilers on Windows.
        function getLogFuncIntelGNU_proc(ndim,Point) result(logFunc) bind(C)
            import :: IK, RK
            integer(IK), intent(in)         :: ndim
            real(RK), intent(in)            :: Point(ndim)
            real(RK)                        :: logFunc
        end function getLogFuncIntelGNU_proc
    end interface

    interface
        subroutine runParaDRAMIntelGNU  ( ndim                  &
                                        , getLogFuncIntelGNU    &
                                        , inputFileVec          &
                                        , inputFileLen          &
                                        ) bind(C, name = "runParaDRAMIntelGNU")
            import :: IK, getLogFuncIntelGNU_proc
            implicit none
            integer(IK) , intent(in)                :: ndim
            procedure(getLogFuncIntelGNU_proc)      :: getLogFuncIntelGNU
            character(1), dimension(*), intent(in)  :: inputFileVec
            integer(IK) , intent(in)                :: inputFileLen
        end subroutine runParaDRAMIntelGNU
    end interface

contains

    subroutine runParaDRAM  ( ndim          &
                            , getLogFunc    &
                            , inputFile     &
                            )
        implicit none
        integer(IK) , intent(in)            :: ndim
        procedure(getLogFunc_proc)          :: getLogFunc
        character(*), intent(in), optional  :: inputFile
        character(1)                        :: inputFileVec(10000) ! (kind=c_char)
        integer(IK)                         :: inputFileLen
        integer(IK)                         :: i
        if (present(inputFile)) then
            inputFileLen = len(inputFile)
            do i = 1, inputFileLen
                inputFileVec(i) = inputFile(i:i)
            end do
        else
            inputFileLen = 0
        end if
        call runParaDRAMIntelGNU( ndim = ndim & ! int(ndim, kind=CIK) &
                                , getLogFuncIntelGNU = getLogFuncIntelGNU &
                                , inputFileVec = inputFileVec(1:inputFileLen) &
                                , inputFileLen = inputFileLen &
                                )
    contains
        function getLogFuncIntelGNU(ndim,Point) result(logFunc) bind(C)
            ! The C-style Fortran objective function wrapper.
            implicit none
            integer(IK) , intent(in)    :: ndim
            real(RK)    , intent(in)    :: Point(ndim)
            real(RK)                    :: logFunc
            logFunc = getLogFunc(ndim,Point)
        end function getLogFuncIntelGNU
    end subroutine runParaDRAM

#else

    interface
        subroutine runParaDRAM  ( ndim          &
                                , getLogFunc    &
                                , inputFile     &
                                )
            import :: IK, getLogFunc_proc
            implicit none
            integer(IK) , intent(in)            :: ndim
            procedure(getLogFunc_proc)          :: getLogFunc
            character(*), intent(in), optional  :: inputFile
        end subroutine runParaDRAM
    end interface

#endif

#endif

end module paramonte

Again, note that there are two interfaces to the ParaMonte library routines, separated via the compiler preprocessor flag IS_COMPATIBLE_COMPILER in the signature module file.

In the following sections, we will describe which interface is suitable for what circumstances.

Calling the ParaMonte samplers: The modern object-oriented style

The Object-Oriented-Programing (OOP) style interface is the parts of the codes in the above that are enclosed by the IS_COMPATIBLE_COMPILER preprocessor flag. The OOP interface can be used only if,

you have built the ParaMonte library from scratch on your system using a Fortran compiler installed on your system (ideally, the Intel ifort Fortran compiler), or,
you intend to use the prebuilt ParaMonte libraries which can be downloaded from the ParaMonte release page on GitHub, but you promise that you are using a compatible compiler to compile your example. Typically this means using the same compiler/library suite and version that has been used to build ParaMonte.

If one of the above two scenarios describes your situation, then you can safely use the OOP-style interface (enclosed between the two preprocessor directives #if defined IS_COMPATIBLE_COMPILER and #else) in the above code.

Tip: As you may have noticed in the above main code, one advantage of the object-oriented style is that upon finishing the simulation, all simulation properties and results are stored in a single object (e.g., pmpd in the above program) which belongs to the ParaDRAM_type class. This information can be then readily used and processed for other purposes within your simulation at runtime.

Using this modern Fortran style, you can also provide optional arguments at the time of calling the ParaMonte routines. For example, the variables chainSize, adaptiveUpdateCount, …, specified in the OOP-style interface of the runSampler() method in the above code would overwrite any values for the corresponding variables specified taken from the input file: ./paramonte.in (unless you set the input specification variable inputFileHasPriority to true from within the input file).

In addition to being able to assign simulation specifications from within an input file, the current major release the ParaMonte library also enables the assignment of all simulation specifications from the Fortran Object-Oriented interface to the ParaMonte sampler routines (except the input specification variable inputFileHasPriority which is always available only from within the input file).

Note: You can view the full object-oriented interface of the runSampler() method of the ParaDRAM_type sampler class in ParaDRAM_mod.inc.f90.

Calling the ParaMonte samplers: The old, less flexible, but more portable style

If you intend to use the prebuilt ParaMonte libraries to build your application (or for some reason, you are building your example with a compiler other than the one you used for building the ParaMonte library), then the modern object-oriented style of calling ParaMonte routines described in the previous section may not work for you. This is mainly due to the potential mismatch between the existing symbols in the ParaMonte library and the symbols that your different compiler expects to see in the ParaMonte library.

In such a case, you can use the simpler procedural Fortran interface to the ParaMonte library routines as shown in the above main code, enclosed between #else and #endif Fortran preprocessor directives. You can view the runParaDRAM() sampler routine’s interface in the signature module file paramonte.f90 shown in the previous sections in the above.

Note: Although, the procedural interface provides more portability, there is, in general, no guarantee that this procedural interface compiled by a certain compiler would match an example main code compiled with a different compiler. However, as of today, the procedural interface to the ParaMonte library routines is interoperable between the Intel and GNU Fortran compilers. This means that you can link your main procedural-interface codes compiled via the GNU Fortran compiler against the ParaMonte libraries prebuilt via the Intel Fortran compiler, and vice versa. While this is true at the moment, the situation may change in the future releases of the two compilers (although it seems unlikely).

Calling the ParaMonte samplers: The input file

We will store the simulation specifications in a separate external input file whose path "./paramonte.in" is given by the variable inputFile in the main code above. See this page for a detailed description of the structure and convention rules of the input files. Also, see this page for detailed descriptions of the simulation specifications of the ParaDRAM sampler.

Tip: Note that the entire contents of the input file is optional. If not provided, all simulation variables will be assigned appropriate default values. Nevertheless, in some cases the user may need to provide further information about the simulation via the input simulation variables to help the ParaMonte routines yield better performance or more accurate results.

Tip: You can also pass the entire contents of the input file as a string value to the input variable inputFile of the ParaMonte routines, instead of passing the-path-to-this-input-file.

Here is the contents of this paramonte.in input file,

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%  
!%  Description:
!%      +   Run the Monte Carlo sampler of the ParaMonte library given the input log-target density function `getLogFunc()`.
!%  Output:
!%      +   The simulation output files.
!%  Author:
!%      +   Computational Data Science Lab, Monday 9:03 AM, May 16 2016, ICES, UT Austin
!%  Visit:
!%      +   https://www.cdslab.org/paramonte
!%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
!   USAGE:
!
!       --  Comments must begin with an exclamation mark (!).
!       --  Comments can appear anywhere on an empty line or, after a value assignment.
!       --  All variable assignments are optional and can be commented out. In such cases, appropriate default values will be assigned.
!       --  Use ParaDRAM namelist (group) name to group a set of ParaDRAM simulation specification variables.
!       --  The order of the input variables in the namelist groups is irrelevant and unimportant.
!       --  Variables can be defined multiple times, but only the last definition will be considered as input.
!       --  All variable names are case insensitive. However, for clarity, this software follows the camelCase code-writing practice.
!       --  String values must be enclosed with either single or double quotation marks.
!       --  Logical values are case-insensitive and can be either .true., true, or t for a TRUE value, and .false., false, or f for a FALSE value.
!       --  All vectors and arrays in the input file begin with index 1. This is following the convention of 
!           the majority of science-oriented programming languages: Fortran, Julia, Mathematica, MATLAB, and R.
!
!      For comprehensive guidelines on the input file organization and rules, visit: 
!   
!           https://www.cdslab.org/paramonte/notes/usage/paradram/input/
!   
!      To see detailed descriptions of each of variables, visit:
!   
!           https://www.cdslab.org/paramonte/notes/usage/paradram/specifications/
!

&ParaDRAM

    ! Base specifications

    description                         = "
This\n
    is a\n
        multi-line\n
            description."                                   ! strings must be enclosed with "" or '' and can be continued on multiple lines.
                                                            ! No comments within strings are allowed.
   !outputColumnWidth                   = 25                ! this is an example of a variable that is commented out and 
                                                            ! therefore, its value won't be read by the sampler routine.
                                                            ! To pass it to the routine, simply remove the ! mark at 
                                                            ! the beginning of the line.
    outputRealPrecision                 = 17
   !outputDelimiter                     = ","
    outputFileName                      = "./out/"          ! the last forward-slash character indicates that this 
                                                            ! is the folder where the output files will have to stored.
                                                            ! However, since no output filename prefix has been specified,
                                                            ! the output filenames will be assigned a randomly-generated prefix.
   !sampleSize                          = 111
    randomSeed                          = 2136275,
    chainFileFormat                     = "compact"
    variableNameList                    = "variable1"       ! Notice the missing fourth variable name here. 
                                        , "variable2"       ! Any missing variable name will be automatically assigned an appropriate name. 
                                        , "variable3"
    domainLowerLimitVec                 = 4*-1.e300         ! repetition pattern rules apply again here. 4 dimensions => 4-element vector of values.
    domainUpperLimitVec                 = 4*1.e300          ! repetition pattern rules apply again here. 4 dimensions => 4-element vector of values.
    parallelizationModel                = "single chain"    ! "singleChain" would also work. Similarly, "multichain", "multi chain", or "multiChain".
   !targetAcceptanceRate                = 0.23e0
    progressReportPeriod                = 1000
    maxNumDomainCheckToWarn             = 100
    maxNumDomainCheckToStop             = 1000
    restartFileFormat                   = "binary"
    overwriteRequested                  = true              ! FALSE, false, .false., .f., and f would be also all valid logical values representing False
    silentModeRequested                 = false             ! FALSE, false, .false., .f., and f would be also all valid logical values representing False
   !mpiFinalizeRequested                = true              ! TRUE, true, .true., .t., and t would be also all valid logical values representing True

    ! MCMC specifications

    chainSize                           = 30000
    startPointVec                       = 4*1.e0            ! four values of 1.e0 are specified here by the repetition pattern symbol *
   !sampleRefinementCount               = 10
    sampleRefinementMethod              = "BatchMeans"
    randomStartPointDomainLowerLimitVec = 4*-100.e0         ! repetition pattern rules apply again here. 4 dimensions => 4-element vector of values.
    randomStartPointDomainUpperLimitVec = 4*100.0           ! repetition pattern rules apply again here. 4 dimensions => 4-element vector of values.
    randomStartPointRequested           = false

    ! DRAM specifications

    scaleFactor                         = "2*0.5*Gelman"    ! The asterisk here means multiplication since it is enclosed within quotation marks.
    proposalModel                       = "normal"          ! or "uniform" as you wish.
    adaptiveUpdateCount                 = 10000000
    adaptiveUpdatePeriod                = 35
    greedyAdaptationCount               = 0
   !delayedRejectionCount               = 5
   !delayedRejectionScaleFactorVec      = 5*1.
   !burninAdaptationMeasure             = 1.
   !proposalStartStdVec                 = 4*1.0             ! repetition pattern rules apply again here. 4 dimensions => 4-element vector of values.
   !proposalStartCorMat                 =   1,0,0,0,        ! 2-dimensional correlation-matrix definition, although it is commented out and won't be read.
   !                                        0,1,0,0,
   !                                        0,0,1,0,
   !                                        0,0,0,1,
   !proposalStartCovMat                 =   100,0,0,0,
   !                                        0,100,0,0,
   !                                        0,0,100,0,
   !                                        0,0,0,100,

/

Compiling and linking to generate the executable

Now, follow the instructions on this page on how to build your Fortran example for serial and parallel simulations on the Operating System (OS) of you choice: Windows, Linux, or macOS / Darwin.

Running the ParaMonte simulations

Once you have built your application, running the simulation is as simple as calling the generated executable on the command-prompt for serial applications, or invoking the MPI/Coarray launcher on the command-prompt for MPI/Coarray parallel simulations. Follow the instructions and tips on this page on how to run your Fortran example on either Windows, Linux, or macOS / Darwin Operating System (OS), on either a single processor or in parallel, on multiple processors or clusters of nodes of processors.

Tip: If you linking your applications (on Linux or macOS) against a shared build of the ParaMonte library, you may encounter runtime errors similar in nature to the following,

main.exe: error while loading shared libraries: libparamonte_c_linux_x64_gnu_release_shared_heap.so: cannot open shared object file: No such file or directory

Such an error message imply that the shared runtime library mentioned in the error message is missing from the known search paths to application. To resolve such runtime errors, you will have to add the path to the required runtime shared library to your Bash environmental variable LD_LIBRARY_PATH. For example, the following will prepend the current directory to the list of paths for the runtime library search.

export LD_LIBRARY_PATH=.:$LD_LIBRARY_PATH

Running the ParaMonte simulations on a single processor

If the executable main.exe has been built for serial simulations (using one processor), simply call the executable on the command line. The calling method is identical on all three Operating Systems: Windows, Linux, and Darwin (macOS). Here is an example call on Windows OS,

D:\example>main.exe

************************************************************************************************************************************
************************************************************************************************************************************
****                                                                                                                            ****
****                                                                                                                            ****
****                                                         ParaMonte                                                          ****
****                                                       Version 1.0.0                                                        ****
****                                                        May 23 2018                                                         ****
****                                                                                                                            ****
****                                                   Department of Physics                                                    ****
****                                              Computational & Data Science Lab                                              ****
****                                          Data Science Program, College of Science                                          ****
****                                            The University of Texas at Arlington                                            ****
****                                                                                                                            ****
****                                                  originally developed at                                                   ****
****                                                                                                                            ****
****                                                 Multiscale Modeling Group                                                  ****
****                                          Center for Computational Oncology (CCO)                                           ****
****                                 Oden Institute for Computational Engineering and Sciences                                  ****
****                               Department of Aerospace Engineering and Engineering Mechanics                                ****
****                                     Department of Neurology, Dell-Seton Medical School                                     ****
****                                            Department of Biomedical Engineering                                            ****
****                                             The University of Texas at Austin                                              ****
****                                                                                                                            ****
****                                   For questions and further information, please contact:                                   ****
****                                                                                                                            ****
****                                                      Amir Shahmoradi                                                       ****
****                                                                                                                            ****
****                                                   shahmoradi@utexas.edu                                                    ****
****                                                  amir.shahmoradi@uta.edu                                                   ****
****                                                   ashahmoradi@gmail.com                                                    ****
****                                                                                                                            ****
****                                                       cdslab.org/pm                                                        ****
****                                                                                                                            ****
****                                             https://www.cdslab.org/paramonte/                                              ****
****                                                                                                                            ****
****                                                                                                                            ****
************************************************************************************************************************************
************************************************************************************************************************************



************************************************************************************************************************************
****                                                                                                                            ****
****                                              Setting up ParaDRAM environment                                               ****
****                                                                                                                            ****
************************************************************************************************************************************


        ParaDRAM - NOTE: Variable outputFileName detected among the input variables to ParaDRAM:
        ParaDRAM - NOTE: ./out/
        ParaDRAM - NOTE:
        ParaDRAM - NOTE: Absolute path to the current working directory:
        ParaDRAM - NOTE: D:\example
        ParaDRAM - NOTE:
        ParaDRAM - NOTE: Generating the requested directory for ParaDRAM output files:
        ParaDRAM - NOTE: .\out\
        ParaDRAM - NOTE:
        ParaDRAM - NOTE: No user-input filename prefix for ParaDRAM output files detected.
        ParaDRAM - NOTE: Generating appropriate filenames for ParaDRAM output files from the current date and time...
        ParaDRAM - NOTE:
        ParaDRAM - NOTE: ParaDRAM output files will be prefixed with:
        ParaDRAM - NOTE: .\out\ParaDRAM_run_20200322_002335_512



        ParaDRAM - NOTE: Searching for previous runs of ParaDRAM...



        ParaDRAM - NOTE: No pre-existing ParaDRAM run detected.
        ParaDRAM - NOTE: Starting a fresh ParaDRAM run...



        ParaDRAM - NOTE: Generating the output report file:
        ParaDRAM - NOTE: .\out\ParaDRAM_run_20200322_002335_512_process_1_report.txt



        ParaDRAM - NOTE: Please see the output report and progress files for further realtime simulation details.



                         Accepted/Total Func. Call   Dynamic/Overall Acc. Rate   Elapsed/Remained Time [s]
                         =========================   =========================   =========================
                                    30000 / 130490              0.179 / 0.2291              2.2880 / 0.000



        ParaDRAM - NOTE: Computing the Markov chain's statistical properties...



        ParaDRAM - NOTE: Computing the final decorrelated sample size...



        ParaDRAM - NOTE: Generating the output sample file:
        ParaDRAM - NOTE: .\out\ParaDRAM_run_20200322_002335_512_process_1_sample.txt



        ParaDRAM - NOTE: Computing the output sample's statistical properties...





        ParaDRAM - NOTE: Mission Accomplished.

Running the ParaMonte simulations in parallel on multiple processors

If the executable main.exe has been built for MPI-enabled parallel simulations (on multiple processors), you will need to invoke the MPI launcher to start the simulation.

On Windows
For parallel simulations on a single node (of multiple processors), for example on a single computer, the command typically looks like the following,
```
mpiexec -localonly -n 3 main.exe
```
where the flag -localonly ensures the simulation is run only on a single node. This option avoids the invocation of the Hydra service which would require prior registration. The flag -n 3 assigns three MPI tasks to three physical processors for the simulation. To run your simulation on a cluster of nodes, follow the guidelines and instructions provided by the Intel MPI library development team.
On Linux / Darwin (Mac)
The MPI launcher command for parallel simulations on Linux and Darwin are identical and look very similar to the Windows command (except for the flag -localonly which is not needed),
```
mpiexec -n 3 main.exe
```
where the flag -n 3 assigns three MPI tasks to three physical processors for the simulation.

Warning: Do not specify more MPI tasks than the available number of physical processors on your system. If there are fewer physical processors available on the machine than the requested MPI tasks, some processors will be assigned more than one MPI tasks. This will significantly degrade the performance of your simulation and will potentially eliminate any benefits gained from running your simulation in parallel. For example, if you are running a parallel simulation on a quad-core laptop (with four physical processors), then the maximum number of MPI tasks should not go beyond 4 under any circumstances.

For example, the following command will launch the simulation on 3 processors on a Windows OS,

D:\example>mpiexec -localonly -n 3 main.exe

************************************************************************************************************************************
************************************************************************************************************************************
****                                                                                                                            ****
****                                                                                                                            ****
****                                                         ParaMonte                                                          ****
****                                                       Version 1.0.0                                                        ****
****                                                        May 23 2018                                                         ****
****                                                                                                                            ****
****                                                   Department of Physics                                                    ****
****                                              Computational & Data Science Lab                                              ****
****                                          Data Science Program, College of Science                                          ****
****                                            The University of Texas at Arlington                                            ****
****                                                                                                                            ****
****                                                  originally developed at                                                   ****
****                                                                                                                            ****
****                                                 Multiscale Modeling Group                                                  ****
****                                          Center for Computational Oncology (CCO)                                           ****
****                                 Oden Institute for Computational Engineering and Sciences                                  ****
****                               Department of Aerospace Engineering and Engineering Mechanics                                ****
****                                     Department of Neurology, Dell-Seton Medical School                                     ****
****                                            Department of Biomedical Engineering                                            ****
****                                             The University of Texas at Austin                                              ****
****                                                                                                                            ****
****                                   For questions and further information, please contact:                                   ****
****                                                                                                                            ****
****                                                      Amir Shahmoradi                                                       ****
****                                                                                                                            ****
****                                                   shahmoradi@utexas.edu                                                    ****
****                                                  amir.shahmoradi@uta.edu                                                   ****
****                                                   ashahmoradi@gmail.com                                                    ****
****                                                                                                                            ****
****                                                       cdslab.org/pm                                                        ****
****                                                                                                                            ****
****                                             https://www.cdslab.org/paramonte/                                              ****
****                                                                                                                            ****
****                                                                                                                            ****
************************************************************************************************************************************
************************************************************************************************************************************



************************************************************************************************************************************
****                                                                                                                            ****
****                                              Setting up ParaDRAM environment                                               ****
****                                                                                                                            ****
************************************************************************************************************************************


        ParaDRAM - NOTE: Variable outputFileName detected among the input variables to ParaDRAM:
        ParaDRAM - NOTE: ./out/
        ParaDRAM - NOTE:
        ParaDRAM - NOTE: Absolute path to the current working directory:
        ParaDRAM - NOTE: D:\example
        ParaDRAM - NOTE:
        ParaDRAM - NOTE: Generating the requested directory for ParaDRAM output files:
        ParaDRAM - NOTE: .\out\
        ParaDRAM - NOTE:
        ParaDRAM - NOTE: No user-input filename prefix for ParaDRAM output files detected.
        ParaDRAM - NOTE: Generating appropriate filenames for ParaDRAM output files from the current date and time...
        ParaDRAM - NOTE:
        ParaDRAM - NOTE: ParaDRAM output files will be prefixed with:
        ParaDRAM - NOTE: .\out\ParaDRAM_run_20200321_193047_392



        ParaDRAM - NOTE: Searching for previous runs of ParaDRAM...



        ParaDRAM - NOTE: No pre-existing ParaDRAM run detected.
        ParaDRAM - NOTE: Starting a fresh ParaDRAM run...



        ParaDRAM - NOTE: Generating the output report file:
        ParaDRAM - NOTE: .\out\ParaDRAM_run_20200321_193047_392_process_1_report.txt



        ParaDRAM - NOTE: Please see the output report and progress files for further realtime simulation details.



                         Accepted/Total Func. Call   Dynamic/Overall Acc. Rate   Elapsed/Remained Time [s]
                         =========================   =========================   =========================
                                    30000 / 130260              0.200 / 0.2298              1.3140 / 0.000



        ParaDRAM - NOTE: The time cost of calling the user-provided objective function must be at least 52.5143 times more
        ParaDRAM - NOTE: (that is, ~24.6651 microseconds) to see any performance benefits from singlechain parallelization
        ParaDRAM - NOTE: model for this simulation. Consider simulating in the serial mode in the future (if used within the
        ParaDRAM - NOTE: same computing environment and with the same configuration as used here).



        ParaDRAM - NOTE: Computing the Markov chain's statistical properties...



        ParaDRAM - NOTE: Computing the final decorrelated sample size...



        ParaDRAM - NOTE: Generating the output sample file:
        ParaDRAM - NOTE: .\out\ParaDRAM_run_20200321_193047_392_process_1_sample.txt



        ParaDRAM - NOTE: Computing the output sample's statistical properties...





        ParaDRAM - NOTE: Mission Accomplished.

You may have noticed the following note in the above simulation output,

        ParaDRAM - NOTE: The time cost of calling the user-provided objective function must be at least 52.5143 times more
        ParaDRAM - NOTE: (that is, ~24.6651 microseconds) to see any performance benefits from singlechain parallelization
        ParaDRAM - NOTE: model for this simulation. Consider simulating in the serial mode in the future (if used within the
        ParaDRAM - NOTE: same computing environment and with the same configuration as used here).

This indicates that the simulation in parallel does not lead to any efficiency gain. Given the computer specifications on which the above simulation is performed (a decently powerful computer), this may not be surprising. Based on the CPU power and inter-processor communication costs, our example simulation in the above seems to be a good candidate for serial simulations instead.

Tip: The optimal number of processors for parallel simulations can vary widely from one MPI-implementation, processor architecture, and simulation problem to another. For large-scale simulations, it is therefore advisable that you run a small simulation for a limited number of samples that can be specified by the appropriate input simulation variable in the input file (e.g., chainSize=2000 for MCMC simulations performed via ParaDRAM sampler). In general, the number of sampled points should be small enough to finish the test-simulation in a few minutes. Then, check the contents of the output file for this small simulation that has the suffix _report.txt and look up the part corresponding to the following simulation information near the bottom of the file,

Average time cost of inter-process communications per used (accepted or rejected or delayed-rejection) function call, in seconds:
2.446639524998980E-005
 
 
Average pure time cost of each function call, in seconds:
4.696832225884184E-007
 
 
Number of processes (images):
3
 
 
Estimated maximum speedup gained via singlechain parallelization model compared to serial mode:
1.904244758247332E-002
 
 
Predicted optimal number of physical computing processes for singlechain parallelization model, given the current MCMC sampling efficiency:
1
 
 
Predicted optimal maximum speedup gained via singlechain parallelization model, given the current MCMC sampling efficiency:
1.00000000000000
 
        ParaDRAM - NOTE: The time cost of calling the user-provided objective function must be at least 52.5143 times more 
        ParaDRAM - NOTE: (that is, ~24.6651 microseconds) to see any performance benefits from singlechain parallelization 
        ParaDRAM - NOTE: model for this simulation. Consider simulating in the serial mode in the future (if used within the 
        ParaDRAM - NOTE: same computing environment and with the same configuration as used here).
 
 
Predicted absolute optimal number of physical computing processes for singlechain parallelization model, under any MCMC sampling efficiency:
1
 
 
Predicted absolute optimal maximum speedup gained via singlechain parallelization model, under any MCMC sampling efficiency:
1.00000000000000
 
        ParaDRAM - NOTE: This simulation will likely NOT benefit from any additional computing processor beyond the predicted 
        ParaDRAM - NOTE: absolute optimal number, 1, in the above. This is true for any value of MCMC sampling efficiency.

At least two numbers in the information provided in the above box merit your attention before you go ahead to set up your full-production simulation:
1. The predicted optimal number of physical computing processes gives the simulation's best estimate for the optimal number of processes that could give you the highest speed-up, given the current sampling efficiency of your sampler. This is important because if your sampler's efficiency is high, then much of the work will be performed by only the very few first processors and all the work done by the rest of processors is wasted. Therefore, if you believe that your sampler's efficiency in the production mode will remain at the same level as in the test run, you may want to set the number of MPI tasks equal to the predicted optimal number by your test simulation.
2. The predicted absolute optimal number of physical computing processes gives the simulation's best estimate for the absolute-maximum number of processes that could give you the highest speed-up, for any conceivable sampling efficiency of your sampler. You may want to make sure that whatever number of processors you request for your production simulation, it would not go beyond this number as it will gradually become detrimental to your simulation's overall performance. Beyond this predicted absolute processor count, the cost of inter-processor communications gradually becomes comparable and eventually more than any time-savings that could be earned by simulating in parallel.
Keep in mind that these predictions for the optimal number of processors are only estimates and can vary widely depending on the runtime setup of your simulation, or the topology of the inter-processor network and communications.

Final notes

Mission accomplished. We have now used the adaptive Markov chain sampler routine of the ParaMonte library to generate random numbers from our objective function of interest. This simulation will generate five files in a folder named ./out/ in the current directory, each of which contains some specific information about the simulation. The purpose of each file is indicated by the suffix at the end of its name: _chain, _sample, _report, _progress, _restart. If multichain parallelism is requested via the input simulation specification variable parallelizationModel = 'multiChain', then there will be five simulation output files generated for each processor separately and independently of each other. This is as if we have simulated several independent but simultaneous ParaDRAM simulations of our objective function. More details about the contents of each of the ParaDRAM simulation output files can be found on this page.

Tip: Self-contained examples with Multivariate Normal Distribution as the objective function are available on the the release page of the ParaMonte library on GitHub. These examples come with build-scripts for both Windows and Linux/macOS that can seamlessly build them for you by simply calling the scripts. Depending on the programming language of your choice and your system's architecture, these self-contained examples minimally require a C/C++/Fortran compiler (e.g., Intel Parallel Studio ≥ 2019.1.1 or Microsoft C/C++ compiler ≥ 2019 on Windows, or GNU C/C++/Fortran compilers ≥ 8.4 on Unix-like Operating Systems). If you follow the installation instructions of the ParaMonte library, then all of these build requirements will be automatically installed on your system.

Post-processing the ParaMonte simulation results

The post-processing of the output files generated by the ParaMonte routines can be nicely and quickly done via either the ParaMonte Python library or the ParaMonte MATLAB library.

If you have any questions about the topics discussed on this page, feel free to ask in the comment section below, or raise an issue on the GitHub page of the library, or reach out to the ParaMonte library authors.