pm_sampling::getErrSampling Interface Reference

Generate and return .true. if the procedure fails to fully accomplish the task of Monte Carlo sampling of the specified input mathematical objective function, otherwise, return .false..
More...

Detailed Description

Generate and return .true. if the procedure fails to fully accomplish the task of Monte Carlo sampling of the specified input mathematical objective function, otherwise, return .false..

This generic interface is the entry point to all ParaMonte Monte Carlo samplers of mathematical density functions.
Although the procedures of this generic interface return a single scalar object of type err_type, the procedures generate massive amounts of information about each simulation which are stored in appropriate external hard drive files.

See,

1. this generic documentation page for more information on the generated output files for samplings performed using the ParaDRAM sampler.

Parameters

[in]	sampler	: The input scalar object that can be of, type paradram_type indicating the use of the ParaDRAM sampler, type paradise_type indicating the use of the ParaDISE sampler, type paranest_type indicating the use of the ParaNest sampler,
[in]	getLogFunc	: The input user-specified function that takes an input contiguous vector of shape state(1:ndim) of type real of kind any supported by the processor (e.g., RK, RK32, RK64, or RK128). On input, state represents a state (point) from within the domain of the user-specified target density function whose function value must be returned. On output, the user-specified procedure getLogFunc() must return the function value corresponding to the input state(1:ndim). The following illustrate the generic interface of getLogFunc(state), function getLogFunc(state) result(logFunc) real(RKG), intent(in), contiguous :: state(:) ! (1 : ndim) real(RKG) :: logFunc end function where the condition all(shape(state) == [ndim]) holds where ndim is the number of dimensions of the target density function. If you are an end user, this is where you stop reading about this input argument. If you are a developer, continue reading about the alternative internal interfaces of this argument for dynamic (MATLAB/Python/R) languages. Setting the preprocessing flag CFI_ENABLED to a non-zero value at the time of library compilation activates the C-Fortran Interoperability (CFI) interface of this generic interface. This requires getLogFunc() to have the following interfaces: When the library is built for OpenMP parallelism (i.e., OMP_ENABLED=1 holds) for MATLAB/Python/R programming languages, the input getLogFunc() must have the following interface: function getLogFunc(logFuncState, ndim, njob, avgTimePerFunCall, avgCommPerFunCall) result(mold) integer(IK), value :: ndim, njob real(RKG), intent(inout) :: logFuncState(ndim + 1, njob) real(RKD), intent(inout) :: avgTimePerFunCall, avgCommPerFunCall real(RKG) :: mold end function where, the variable ndim is the number of dimensions of the target density function and the variable njob is the number of states for which the function must be evaluated concurrently. the constant RKD refers to the real kind double precision supported by the compiler. the output mold merely dictates the real kind of the interface and its value is ignored. on input, the subset logFuncState(2:ndim+1, 1:njob) represents the set of njob states (points) from within the domain of the user-specified target density function whose function values must be returned. On output, the user-specified procedure getLogFunc() must return njob function values in elements logFuncState(1, 1:njob) corresponding to the input states logFuncState(2:ndim+1, 1:njob). the two extra intent(inout) arguments avgTimePerFunCall, avgCommPerFunCall represent: the average time (in seconds) it takes to compute the user-specified function for a single input state. the average time (in seconds) it takes to create images/processes/threads for parallel computation of the user-specified function and communicate information among them, collectively called parallelization overhead. These two extra arguments are essential for accurate post-processing the performance of the parallel computations in the simulation. If there is no parallelization within he user-specified getLogFunc, then the parallelization overhead can be set to zero. This procedure interface is entirely different from the density function interfaces in the previous versions of the ParaMonte library. The new interface has been primarily developed to facilitate the concurrent or parallel evaluation of target density function in higher-level languages for a set of njob input states. This interface is particularly vital in higher-level programming languages such as MATLAB, Python, and R where ParaMonte thread-level parallelism is impossible due to the global interpreter lock (GIL). This interface is currently activated if and only if the preprocessors OMP_ENABLED=1 and either MATLAB_ENABLED=1 or PYTHON_ENABLED=1 or R_ENABLED=1 are set. When the library is built for any build configuration other than the above with preprocessor CFI_ENABLED=1, particularly, for the C and C++ programming languages, the input getLogFunc() must have the following interface: function getLogFunc(state, ndim) result(logFunc) integer(IK), value :: ndim real(RKG), intent(in) :: state(ndim) real(RKG) :: logFunc end function where, the variable ndim is the number of dimensions of the target density function. the variable state is an input point from within the domain of the density function. the constant RKG refers to the real kind used for computations. the output logFunc is the computed value of the density function at the specified input state.
[in]	ndim	: The input scalar integer of default kind IK representing the number of dimensions of the domain of the objective function that is to be sampled.

Returns

err : The output scalar of type err_type whose component occur is set to the .true. if and only if the procedure fails to accomplish the sampling task, in which case, the msg component of err is set to a description of the error origin.
In such a case, the description of the possible cause(s) of the error will be written to the output report file that is automatically generated for all simulations.
The primary reason for simulation failures is a wrong syntactic or semantic input value(s) for the simulation settings specified within the components of the input argument sampler or within the external output file specified by the inputFile component of the input sampler.
On output, the msg component of err is an empty string if the sampling completes successfully.

Possible calling interfaces ⛓

use pm_sampling, only: getErrSampling, err_type, IK
type(err_type) :: err
integer(IK) :: ndim
 
err = getErrSampling(sampler, getLogFunc, ndim)

Warning

The condition 0 < ndim must hold for the corresponding input arguments.
These conditions are verified only if the library is built with the preprocessor macro CHECK_ENABLED=1.

Remarks

The procedures under discussion are impure.

Example usage ⛓

program example
 
    use pm_kind, only: RKG => RK
    use pm_kind, only: SK, IK, LK
    use pm_io, only: display_type
    use pm_sampling, only: getErrSampling, paradram_type
    use pm_matrixInit, only: getMatInit, uppLowDia
    use pm_arrayFill, only: getFilled
    use pm_sysPath, only: glob
    use pm_err, only: err_type
 
    implicit none
 
    type(err_type) :: err
    type(display_type) :: disp
    integer(IK), parameter :: NDIM = 4
    real(RKG), parameter :: MEAN(*) = 4 * [real(RKG) :: -3, -1, 1, 3]
 
    disp = display_type(file = SK_"main.out.F90")
 
    call disp%skip()
    call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
    call disp%show("! Simple minimally-informed simulation.")
    call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
    call disp%skip()
 
    block
        call disp%skip()
        call disp%show("NDIM")
        call disp%show( NDIM )
        call disp%show("err = getErrSampling(paradram_type(parallelismMpiFinalizeEnabled = .false._LK), getLogFunc, NDIM)")
                        err = getErrSampling(paradram_type(parallelismMpiFinalizeEnabled = .false._LK), getLogFunc, NDIM)
        call disp%show("if (err%occurred) error stop err%msg")
                        if (err%occurred) error stop err%msg
        call disp%skip()
    end block
 
    call disp%skip()
    call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
    call disp%show("! Optionally-informed simulation.")
    call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
    call disp%skip()
 
    block
        type(paradram_type) :: sampler
        call disp%skip()
        call disp%show("NDIM")
        call disp%show( NDIM )
       ! GNU gfortran yields a segfault with the following default constructor.
       !call disp%show("sampler = paradram_type(parallelismMpiFinalizeEnabled = .false._LK, outputFileName = 'mvn', outputChainSize = 100000, outputStatus = 'retry')")
       !                sampler = paradram_type(parallelismMpiFinalizeEnabled = .false._LK, outputFileName = 'mvn', outputChainSize = 100000, outputStatus = 'retry')
        call disp%show("sampler = paradram_type()")
                        sampler = paradram_type()
        call disp%show("sampler%parallelismMpiFinalizeEnabled = .false._LK; sampler%outputFileName = 'mvn'; sampler%outputChainSize = 100000; sampler%outputStatus = 'retry'")
                        sampler%parallelismMpiFinalizeEnabled = .false._LK; sampler%outputFileName = 'mvn'; sampler%outputChainSize = 100000; sampler%outputStatus = 'retry'
        call disp%show("err = getErrSampling(sampler, getLogFunc, NDIM)")
                        err = getErrSampling(sampler, getLogFunc, NDIM)
        call disp%show("if (err%occurred) error stop err%msg")
                        if (err%occurred) error stop err%msg
        call disp%show("glob(pattern = sampler%outputFileName)")
        call disp%show( glob(pattern = sampler%outputFileName) )
        call disp%skip()
    end block
 
contains
 
    recursive function getLogFunc(state) result(logFunc)
        use pm_distMultiNorm, only: getMultiNormLogPDF
        real(RKG), intent(in), contiguous :: state(:)
        real(RKG) :: logFunc
        integer :: i
        logFunc = getMultiNormLogPDF(state, mean = real(MEAN, RKG))
#if     OMP_ENABLED
        ! We are done. But kill some time for an illustration of parallel sampling.
        do i = 1, 500
            logFunc = logFunc + getMultiNormLogPDF(state, mean = 4 * [real(RKG) :: -3, -1, 1, 3]) * merge(1, -1, mod(i, 2) == 0)
        end do
#endif
    end function
 
    !subroutine setLogFunc(logFuncState)
    !    use pm_distMultiNorm, only: getMultiNormLogPDF
    !    real(RKG), intent(inout), contiguous :: logFuncState(0:,:)
    !    integer(IK) :: ithread
    !    do ithread = 1, size(logFuncState, 2)
    !        logFuncState(0, ithread) = getMultiNormLogPDF(logFuncState(1:, ithread), mean = real(MEAN, RKG))
    !    end do
    !end subroutine
 
end program example

Example Unix compile command via Intel ifort compiler ⛓

#!/usr/bin/env sh
rm main.exe
ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
./main.exe

Example Windows Batch compile command via Intel ifort compiler ⛓

del main.exe
set PATH=..\..\..\lib;%PATH%
ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
main.exe

Example Unix / MinGW compile command via GNU gfortran compiler ⛓

#!/usr/bin/env sh
rm main.exe
gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
./main.exe

Example output ⛓

 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Simple minimally-informed simulation.
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
NDIM
+4
err = getErrSampling(paradram_type(parallelismMpiFinalizeEnabled = .false._LK), getLogFunc, NDIM)
if (err%occurred) error stop err%msg
 
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Optionally-informed simulation.
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
NDIM
+4
sampler = paradram_type()
sampler%parallelismMpiFinalizeEnabled = .false._LK; sampler%outputFileName = 'mvn'; sampler%outputChainSize = 100000; sampler%outputStatus = 'retry'
err = getErrSampling(sampler, getLogFunc, NDIM)
if (err%occurred) error stop err%msg
glob(pattern = sampler%outputFileName)
 
 
 

Example usage ⛓

module logfunc
    use pm_kind, only: IK, RKG => RKS ! All real kinds are supported.
    implicit none
    integer(IK) , parameter :: NDIM = 4_IK
    real(RKG)   , parameter :: MEAN(NDIM) = [-6, -2, 2, 6] ! The mean vector of the MVN
    real(RKG)   , parameter :: LOG_INVERSE_SQRT_TWO_PI = log(1./sqrt(2.*acos(-1.))) ! MVN distribution coefficient: log(1/sqrt(2*Pi)^ndim).
    real(RKG)   , parameter :: COVMAT(NDIM, NDIM) = reshape([ +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 ], shape(COVMAT)) ! The covariance matrix of the MVN
    real(RKG)   , parameter :: INVCOV(NDIM, NDIM) = reshape([ +1.6, -0.4, -0.4, -0.4 &
                                                            , -0.4, +1.6, -0.4, -0.4 &
                                                            , -0.4, -0.4, +1.6, -0.4 &
                                                            , -0.4, -0.4, -0.4, +1.6 ], shape(INVCOV)) ! The inverse covariance matrix of the MVN
    real(RKG)   , parameter :: MVN_COEF = NDIM * LOG_INVERSE_SQRT_TWO_PI + 0.581575404902840 ! log(sqrt(det(INVCOV)))
contains
    function getLogFunc(state) result(logFunc)
        ! Return the negative natural logarithm of MVN distribution 
        ! evaluated at the input vector `state` of length `ndim`.
        real(RKG), intent(in), contiguous :: state(:)
        real(RKG) :: stateNormed(size(state))
        real(RKG) :: logFunc
        stateNormed = state - MEAN
        logFunc = MVN_COEF - 0.5_RKG * (dot_product(stateNormed, matmul(INVCOV, stateNormed)))
    end function
end module logfunc
 
program example
    use logfunc, only: RKG, NDIM, getLogFunc
    use pm_sampling, only: getErrSampling, paradram_type, err_type
    type(err_type) :: err
    err = getErrSampling(paradram_type(inputFile = 'input.nml'), getLogFunc, NDIM)
    if (err%occurred) error stop "sampler failed: "//err%msg
end

The contents of example specifications namelist input file input.nml ⛓

!   DESCRIPTION:
!
!       The external input file for sampling the 4-dimensional Multivariate Normal distribution function as implemented in the accompanying source files.
!       This file is common between all supported programming language environments.
!
!   NOTE:
!
!       All simulation specifications (including this whole file) are optional and can be nearly safely commented out.
!       However, if domain boundaries are finite, set them explicitly.
!
!   USAGE:
!
!       --  Comments must begin with an exclamation mark `!`.
!       --  Comments can appear anywhere on an empty line or, after a variable assignment
!           (but not in the middle of a variable assignment whether in single or multiple lines).
!       --  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/generic/latest/usage/sampling/paradram/input/
!
!      To see detailed descriptions of each of variables, visit:
!
!           https://www.cdslab.org/paramonte/generic/latest/usage/sampling/paradram/specifications/
!
&paradram
 
    ! Base specifications.
 
    description                         = "
This\n
    is a\n
        multi-line\n
            description.\\n"                                    ! strings must be enclosed with "" or '' and can be continued on multiple lines.
                                                                ! No comments within strings are allowed.
    domain                              = "cube"
    domainAxisName                      = "variable1"
                                          "variable2"           ! values can appear in multiple lines.
    domainBallAvg                       = 0 0 0 0               ! values can be separated with blanks or commas.
    domainBallCor                       = 1 0 0 0
                                          0 1 0 0
                                          0 0 1 0
                                          0 0 0 1
    domainBallCov                       = 1 0 0 0
                                          0 1 0 0
                                          0 0 1 0
                                          0 0 0 1
    domainBallStd                       = 1 1 1 1
    domainCubeLimitLower                = 4*-1.e10              ! repetition pattern rules apply here. 4 dimensions => 4-element vector of values.
    domainCubeLimitUpper(1)             = +1.e10                ! Elements of vectors can be set individually.
    domainCubeLimitUpper(2:4)           = 3*+1.e10              ! Elements of vectors can be set individually.
    domainErrCount                      = 100
    domainErrCountMax                   = 1000
    inputFileHasPriority                = FALSE                 ! This is relevant only to simulations within the Fortran programming language.
    outputChainFileFormat               = "compact"
   !outputColumnWidth                   = 25                    ! This is an example of a variable that is commented out.
                                                                ! Therefore, its value will not be read by the sampler routine.
                                                                ! To pass it to the routine, simply remove the `!` mark at the beginning of the line.
    outputFileName                      = "./out/mvn"           ! A forward-slash character at the end of the string value would indicate the specified path
                                                                ! is to be interpreted as the name of the folder to contain the simulation output files.
                                                                ! The base name for the simulation output files will be generated from the current date and time.
                                                                ! Otherwise, the specified base name at the end of the string will be used in naming the simulation output files.
    outputPrecision                     = 17
    outputReportPeriod                  = 1000
    outputRestartFileFormat             = "ascii"
    outputSampleSize                    = -1
    outputSeparator                     = ","
    outputSplashMode                    = "normal"              ! or quiet or silent.
    outputStatus                        = "retry"               ! or extend.
    parallelism                         = "multi chain"         ! "singleChain" would also work. Similarly, "multichain", "multi chain", or "multiChain".
    parallelismMpiFinalizeEnabled       = false                 ! TRUE, true, .true., .t., and t would be also all valid logical values representing truth.
    parallelismNumThread                = 3                     ! number of threads to use in shared-memory parallelism.
   !randomSeed                          = 2136275,
   !targetAcceptanceRate                = 0.23e0
 
    ! MCMC specifications.
 
    outputChainSize                     = 10000
    outputSampleRefinementCount         = 10
    outputSampleRefinementMethod        = "BatchMeans"
    proposal                            = "normal"              ! or "uniform" as you wish.
    proposalCor(:, 1)                   = 1 0 0 0               ! first matrix column.
    proposalCor(:, 2:4)                 = 0 1 0 0
                                          0 0 1 0
                                          0 0 0 1               ! other matrix columns.
    proposalCov                         = 1 0 0 0
                                          0 1 0 0
                                          0 0 1 0
                                          0 0 0 1               ! or specify all matrix elements in one statement.
    proposalScale                       = "2*0.5*Gelman"        ! The asterisk here means multiplication since it is enclosed within quotation marks.
   !proposalStart                       = 4*1.e0                ! four values of 1.e0 are specified here by the repetition pattern symbol *
    proposalStartDomainCubeLimitLower   = 4*-10.e0              ! repetition pattern rules apply again here. 4 dimensions => 4-element vector of values.
    proposalStartDomainCubeLimitUpper   = 4*+10.e0              ! repetition pattern rules apply again here. 4 dimensions => 4-element vector of values.
    proposalStartRandomized             = false
    proposalStd                         = 4*1.0                 ! repetition pattern rules apply again here. 4 dimensions => 4-element vector of values.
 
    ! DRAM specifications.
 
    proposalAdaptationBurnin            = 1.
    proposalAdaptationCount             = 10000000
    proposalAdaptationCountGreedy       = 0
    proposalAdaptationPeriod            = 35
    proposalDelayedRejectionCount       = 5
    proposalDelayedRejectionScale       = 4*1., 2.              ! The first four elements are 1, followed by 2.
 
/

Example Unix compile command via Intel ifort compiler ⛓

#!/usr/bin/env sh
rm main.exe
ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
./main.exe

Example Windows Batch compile command via Intel ifort compiler ⛓

del main.exe
set PATH=..\..\..\lib;%PATH%
ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
main.exe

Example Unix / MinGW compile command via GNU gfortran compiler ⛓

#!/usr/bin/env sh
rm main.exe
gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
./main.exe

Postprocessing of the example output ⛓

#!/usr/bin/env python
 
import matplotlib.pyplot as plt
#import paramonte as pm
import pandas as pd
import numpy as np
import glob
import sys
import os
funcname = os.path.basename(os.getcwd())
 
linewidth = 2
fontsize = 17
files = glob.glob("./out/*_chain.txt")
 
for file in files:
 
    df = pd.read_csv(file, delimiter = ",")
    sindex = list(df.columns).index("sampleLogFunc") + 1 # start column index.
 
    # traceplot
 
    #print(df.values)
    fig = plt.figure(figsize = (8, 6))
    ax = plt.subplot(1,1,1)
    ax.plot ( range(len(df.values[:, 0]))
            , df.values[:, sindex:]
            , zorder = 1000
            )
    plt.minorticks_on()
    ax.set_xscale("log")
    ax.set_xlabel("MCMC Count", fontsize = 17)
    ax.set_ylabel("MCMC State", fontsize = 17)
    ax.tick_params(axis = "x", which = "minor")
    ax.tick_params(axis = "y", which = "minor")
    plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
    #ax.legend(df.columns.values, fontsize = fontsize)
    plt.tight_layout()
    plt.savefig(funcname + ".traceplot.png")
 
    # scatterplot
 
    if len(df.values[1, sindex:]) > 1:
        #print(df.values)
        fig = plt.figure(figsize = (8, 6))
        ax = plt.subplot(1,1,1)
        ax.scatter  ( df.values[:, sindex]
                    , df.values[:, sindex + 1]
                    , zorder = 1000
                    , alpha = 0.5
                    , s = 1
                    )
        plt.minorticks_on()
        #ax.set_xscale("log")
        ax.set_xlabel(df.columns.values[sindex], fontsize = 17)
        ax.set_ylabel(df.columns.values[sindex + 1], fontsize = 17)
        ax.tick_params(axis = "x", which = "minor")
        ax.tick_params(axis = "y", which = "minor")
        plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
        #ax.legend(df.columns.values, fontsize = fontsize)
        plt.tight_layout()
        plt.savefig(funcname + ".scatterplot.png")
 
    # adaptation measure.
 
    #print(df.values)
    if "proposalAdaptation" in df.columns.values:
        if any(df["proposalAdaptation"].values != 0):
            fig = plt.figure(figsize = (8, 6))
            ax = plt.subplot(1,1,1)
            ax.scatter  ( range(len(df.values[:, 0]))
                        , df["proposalAdaptation"].values
                        , zorder = 1000
                        , alpha = 0.5
                        , c = "red"
                        , s = 1
                        )
            plt.minorticks_on()
            #ax.set_xscale("log")
            ax.set_yscale("log")
            ax.set_xlabel("MCMC Count", fontsize = 17)
            ax.set_ylabel("proposalAdaptation", fontsize = 17)
            ax.tick_params(axis = "x", which = "minor")
            ax.tick_params(axis = "y", which = "minor")
            plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
            #ax.legend(df.columns.values, fontsize = fontsize)
            plt.tight_layout()
            plt.savefig(funcname + ".proposalAdaptation.png")
 
    #sim = pm.ParaDRAM()
    #sample = sim.readSample("./out/", renabled = True)
    #sample[0].plot.grid()

Visualization of the example output ⛓

Example usage ⛓

module logfunc
    use pm_kind, only: IK, RKG => RKD ! all processor kinds are supported.
    integer(IK), parameter :: NDIM = 2
contains
    recursive function getLogFunc(state) result(logFunc)
        ! Return the negative natural logarithm of Himmelblau function evaluated at the input vector state.
        real(RKG), intent(in), contiguous :: state(:)
        real(RKG) :: logFunc
        if (size(state) /= 2) error stop "The input state vector size must be 2."
        logFunc = -log((state(1)**2 + state(2) - 11)**2 + (state(1) + state(2)**2 - 7)**2 + 0.1_RKG)
    end function
end module logfunc
 
program example
    use logfunc, only: NDIM, getLogFunc
    use pm_sampling, only: getErrSampling, paradram_type, err_type
    type(paradram_type) :: sampler
    type(err_type) :: err
    ! \bug
    ! gfortran bug requires setting sampler components explicitly.
    sampler = paradram_type()
    sampler%inputFile = "input.nml"
    sampler%outputFileName = "./out/himmelblau"
    err = getErrSampling(sampler, getLogFunc, NDIM)
    if (err%occurred) error stop "sampler failed: "//err%msg
end program example

The contents of example specifications namelist input file input.nml ⛓

!   DESCRIPTION:
!
!       The external input file for sampling the 2-dimensional Himmelblau function as implemented in the accompanying source files.
!       This file is common between all supported programming language environments.
!
!   NOTE:
!
!       All simulation specifications (including this whole file) are optional and can be nearly safely commented out.
!       However, if domain boundaries are finite, set them explicitly.
!
!   USAGE:
!
!       --  Comments must begin with an exclamation mark `!`.
!       --  Comments can appear anywhere on an empty line or, after a variable assignment
!           (but not in the middle of a variable assignment whether in single or multiple lines).
!       --  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/generic/latest/usage/sampling/paradram/input/
!
!      To see detailed descriptions of each of variables, visit:
!
!           https://www.cdslab.org/paramonte/generic/latest/usage/sampling/paradram/specifications/
!
&paradram
 
    ! Base specifications.
 
    description                         = "
This\n
    is a\n
        multi-line\n
            description.\\n"                                    ! strings must be enclosed with "" or '' and can be continued on multiple lines.
                                                                ! No comments within strings are allowed.
    domain                              = "cube"
    domainAxisName                      = "variable1"
                                          "variable2"           ! values can appear in multiple lines.
    domainBallAvg                       = 0 0                   ! values can be separated with blanks or commas.
    domainBallCor                       = 1 0
                                            0 1
    domainBallCov                       = 1, 0
                                            0, 1
    domainBallStd                       = 1, 1
    domainCubeLimitLower                = 2*-1.e30              ! repetition pattern rules apply here. 2 dimensions => 2-element vector of values.
    domainCubeLimitUpper(1)             = +1.e30                ! Elements of vectors can be set individually.
    domainCubeLimitUpper(2)             = +1.e30                ! Elements of vectors can be set individually.
    domainErrCount                      = 100
    domainErrCountMax                   = 1000
    inputFileHasPriority                = FALSE                 ! This is relevant only to simulations within the Fortran programming language.
    outputChainFileFormat               = "compact"
   !outputColumnWidth                   = 25                    ! This is an example of a variable that is commented out.
                                                                ! Therefore, its value will not be read by the sampler routine.
                                                                ! To pass it to the routine, simply remove the `!` mark at the beginning of the line.
    outputFileName                      = "./out/himmelblau"    ! A forward-slash character at the end of the string value would indicate the specified path
                                                                ! is to be interpreted as the name of the folder to contain the simulation output files.
                                                                ! The base name for the simulation output files will be generated from the current date and time.
                                                                ! Otherwise, the specified base name at the end of the string will be used in naming the simulation output files.
    outputPrecision                     = 17
    outputReportPeriod                  = 1000
    outputRestartFileFormat             = "ascii"
    outputSampleSize                    = -1
    outputSeparator                     = ","
    outputSplashMode                    = "normal"              ! or quiet or silent.
    outputStatus                        = "extend"              ! or retry.
    parallelism                         = "multi chain"         ! "singleChain" would also work. Similarly, "multichain", "multi chain", or "multiChain".
    parallelismMpiFinalizeEnabled       = false                 ! TRUE, true, .true., .t., and t would be also all valid logical values representing truth.
    parallelismNumThread                = 3                     ! number of threads to use in shared-memory parallelism.
   !randomSeed                          = 2136275,
   !targetAcceptanceRate                = 0.23e0
 
    ! MCMC specifications.
 
    outputChainSize                     = 10000
    outputSampleRefinementCount         = 10
    outputSampleRefinementMethod        = "BatchMeans"
    proposal                            = "normal"              ! or "uniform" as you wish.
    proposalCor(:, 1)                   = 1 0                   ! first matrix column.
    proposalCor(:, 2)                   = 0 1                   ! second matrix column.
    proposalCov                         = 1 0 0 1               ! or specify all matrix elements in one statement.
    proposalScale                       = "2*0.5*Gelman"        ! The asterisk here means multiplication since it is enclosed within quotation marks.
    proposalStart                       = 2*1.e0                ! four values of 1.e0 are specified here by the repetition pattern symbol *
    proposalStartDomainCubeLimitLower   = 2*-100.e0             ! repetition pattern rules apply again here. 2 dimensions => 2-element vector of values.
    proposalStartDomainCubeLimitUpper   = 2*+100.e0             ! repetition pattern rules apply again here. 2 dimensions => 2-element vector of values.
    proposalStartRandomized             = false
    proposalStd                         = 2*1.0                 ! repetition pattern rules apply again here. 2 dimensions => 2-element vector of values.
 
    ! DRAM specifications.
 
    proposalAdaptationBurnin            = 1.
    proposalAdaptationCount             = 10000000
    proposalAdaptationCountGreedy       = 0
    proposalAdaptationPeriod            = 35
    proposalDelayedRejectionCount       = 5
    proposalDelayedRejectionScale       = 4*1., 2.              ! The first four elements are 1, followed by 2.
 
/

Example Unix compile command via Intel ifort compiler ⛓

#!/usr/bin/env sh
rm main.exe
ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
./main.exe

Example Windows Batch compile command via Intel ifort compiler ⛓

del main.exe
set PATH=..\..\..\lib;%PATH%
ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
main.exe

Example Unix / MinGW compile command via GNU gfortran compiler ⛓

#!/usr/bin/env sh
rm main.exe
gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
./main.exe

Postprocessing of the example output ⛓

#!/usr/bin/env python
 
import matplotlib.pyplot as plt
#import paramonte as pm
import pandas as pd
import numpy as np
import glob
import sys
import os
funcname = os.path.basename(os.getcwd())
 
linewidth = 2
fontsize = 17
files = glob.glob("./out/*_chain.txt")
 
for file in files:
 
    df = pd.read_csv(file, delimiter = ",")
    sindex = list(df.columns).index("sampleLogFunc") + 1 # start column index.
 
    # traceplot
 
    #print(df.values)
    fig = plt.figure(figsize = (8, 6))
    ax = plt.subplot(1,1,1)
    ax.plot ( range(len(df.values[:, 0]))
            , df.values[:, sindex:]
            , zorder = 1000
            )
    plt.minorticks_on()
    ax.set_xscale("log")
    ax.set_xlabel("MCMC Count", fontsize = 17)
    ax.set_ylabel("MCMC State", fontsize = 17)
    ax.tick_params(axis = "x", which = "minor")
    ax.tick_params(axis = "y", which = "minor")
    plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
    #ax.legend(df.columns.values, fontsize = fontsize)
    plt.tight_layout()
    plt.savefig(funcname + ".traceplot.png")
 
    # scatterplot
 
    if len(df.values[1, sindex:]) > 1:
        #print(df.values)
        fig = plt.figure(figsize = (8, 6))
        ax = plt.subplot(1,1,1)
        ax.scatter  ( df.values[:, sindex]
                    , df.values[:, sindex + 1]
                    , zorder = 1000
                    , alpha = 0.5
                    , s = 1
                    )
        plt.minorticks_on()
        #ax.set_xscale("log")
        ax.set_xlabel(df.columns.values[sindex], fontsize = 17)
        ax.set_ylabel(df.columns.values[sindex + 1], fontsize = 17)
        ax.tick_params(axis = "x", which = "minor")
        ax.tick_params(axis = "y", which = "minor")
        plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
        #ax.legend(df.columns.values, fontsize = fontsize)
        plt.tight_layout()
        plt.savefig(funcname + ".scatterplot.png")
 
    # adaptation measure.
 
    #print(df.values)
    if "proposalAdaptation" in df.columns.values:
        if any(df["proposalAdaptation"].values != 0):
            fig = plt.figure(figsize = (8, 6))
            ax = plt.subplot(1,1,1)
            ax.scatter  ( range(len(df.values[:, 0]))
                        , df["proposalAdaptation"].values
                        , zorder = 1000
                        , alpha = 0.5
                        , c = "red"
                        , s = 1
                        )
            plt.minorticks_on()
            #ax.set_xscale("log")
            ax.set_yscale("log")
            ax.set_xlabel("MCMC Count", fontsize = 17)
            ax.set_ylabel("proposalAdaptation", fontsize = 17)
            ax.tick_params(axis = "x", which = "minor")
            ax.tick_params(axis = "y", which = "minor")
            plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
            #ax.legend(df.columns.values, fontsize = fontsize)
            plt.tight_layout()
            plt.savefig(funcname + ".proposalAdaptation.png")
 
    #sim = pm.ParaDRAM()
    #sample = sim.readSample("./out/", renabled = True)
    #sample[0].plot.grid()

Visualization of the example output ⛓

Example usage ⛓

module auxil_mod
 
    use pm_kind, only: RKG => RKD ! All real kinds are supported.
 
    implicit none
 
    real(RKG), parameter :: EPS = sqrt(epsilon(0._RKG))
    real(RKG), parameter :: LARGE = sqrt(huge(0._RKG))
 
end module auxil_mod
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
module data_mod
 
    use auxil_mod, only: RKG
    use pm_kind, only: IK
    implicit none
 
    private
    public :: obs_type, posuniv_type, stat_type, data_type
 
    type, abstract :: obs_type
    end type
 
    type, extends(obs_type) :: univ_type
        real(RKG) :: val
    end type
 
    type, extends(univ_type) :: posuniv_type
        real(RKG) :: log
    end type
 
    interface posuniv_type
        module procedure :: posuniv_typer
    end interface
 
    type stat_type
        real(RKG) :: lim(2)
    end type
 
    type :: data_type
        integer(IK) :: nobs
        type(stat_type) :: stat
        class(obs_type), allocatable :: obs(:)
    end type
 
    interface data_type
        module procedure :: data_typer
    end interface
 
contains
 
    function posuniv_typer(val) result(self)
        real(RKG), intent(in), contiguous :: val(:)
        type(posuniv_type), allocatable :: self(:)
        integer(IK) :: iobs
        allocate(self(size(val)))
        do concurrent(iobs = 1 : size(val))
            self(iobs)%log = log(val(iobs))
            self(iobs)%val = val(iobs)
        end do
    end function
 
    function data_typer(obs, stat) result(self)
        class(obs_type), intent(in), contiguous :: obs(:)
        type(stat_type), intent(in) :: stat
        type(data_type) :: self
        self%stat = stat
        self%nobs = size(obs, 1, IK)
        allocate(self%obs, source = obs)
    end function
 
end module data_mod
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
module domain_mod
 
    !use pm_container, only: css_type
    use auxil_mod, only: RKG
    use pm_kind, only: SK
    implicit none
 
    type :: domain_type
        real(RKG), allocatable :: lower(:)
        real(RKG), allocatable :: upper(:)
        !type(css_type), allocatable :: names(:)
        character(63, SK), allocatable :: names(:)
    end type
 
end module domain_mod
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
module model_mod
 
    use pm_kind, only: SK, IK
    use data_mod, only: obs_type
    use domain_mod, only: domain_type
    use auxil_mod, only: RKG
    implicit none
 
    type, abstract :: model_type
        type(domain_type) :: domain
        !real(RKG), allocatable :: param(:)
        integer(IK) :: npar
    contains
        procedure(getParam_proc), deferred :: getParam
        procedure(setParam_proc), deferred :: setParam
        procedure(getLogPDF_proc), deferred :: getLogPDF
    end type
 
    type :: con_type
        class(model_type), allocatable :: model
    end type
 
    interface con_type
        module procedure :: con_typer
    end interface
 
    abstract interface
        function getParam_proc(self) result(param)
            import :: RKG, model_type
            class(model_type), intent(in) :: self
            real(RKG) :: param(self%npar)
        end function
        subroutine setParam_proc(self, param)
            import :: RKG, model_type
            class(model_type), intent(inout) :: self
            real(RKG), intent(in), contiguous :: param(:)
        end subroutine
        impure elemental function getLogPDF_proc(self, obs) result(logPDF)
            import :: RKG, model_type, obs_type
            class(model_type), intent(in) :: self
            class(obs_type), intent(in) :: obs
            real(RKG) :: logPDF
        end function
    end interface
 
contains
 
    pure elemental function con_typer(model) result(self)
        class(model_type), intent(in) :: model
        type(con_type) :: self
        self%model = model
    end function
 
end module model_mod
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
module lognorm_mod
 
    use model_mod, only: domain_type, model_type, IK
    use auxil_mod, only: RKG, LARGE
    implicit none
 
    type, extends(model_type) :: lognorm_type
        real(RKG), private :: avg, invstd, loginvstd
    contains
        procedure :: getParam, setParam, getLogPDF
    end type
 
    interface lognorm_type
        module procedure :: lognorm_typer
    end interface
 
contains
 
    function lognorm_typer(param, domain) result(self)
        real(RKG), intent(in), contiguous, optional :: param(:)
        type(domain_type), intent(in), optional :: domain
        type(lognorm_type) :: self
        self%npar = 2
        if (present(domain)) then
            self%domain = domain
        else
            self%domain = domain_type([-LARGE, -LARGE], [+LARGE, +LARGE], [character(63) :: "lognorm_avg", "lognorm_logstd"])
        end if
        if (present(param)) call self%setParam(param)
    end function
 
    subroutine setParam(self, param)
        use pm_kind, only: SK, IK
        class(lognorm_type), intent(inout) :: self
        real(RKG), intent(in), contiguous :: param(:)
        if (size(param) /= self%npar) error stop "`lognorm_type%setParam()` takes a vector of two parameters representing `avg` and `logstd`."
        self%invstd = exp(-param(2))
        self%loginvstd = -param(2)
        self%avg = param(1)
    end subroutine
 
    function getParam(self) result(param)
        class(lognorm_type), intent(in) :: self
        real(RKG) :: param(self%npar)
        param = [self%avg, -self%loginvstd]
    end function
 
    impure elemental function getLogPDF(self, obs) result(logPDF)
        use pm_distNorm, only: setNormLogPDF
        use data_mod, only: obs_type, posuniv_type
        class(lognorm_type), intent(in) :: self
        class(obs_type), intent(in):: obs
        real(RKG) :: logPDF
        select type (obs)
        type is (posuniv_type)
            call setNormLogPDF(logPDF, obs%log, mu = self%avg, invSigma = self%invstd, logInvSigma = self%loginvstd)
        class default
            error stop "Unrecognized data."
        end select
    end function
 
end module lognorm_mod
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
module flatPoweto_mod
 
    use auxil_mod, only: RKG, LARGE
    use model_mod, only: domain_type, model_type
    implicit none
 
    type, extends(model_type) :: flatPoweto_type
        real(RKG), private :: logbreak, alpha, break, alphap1, logPDFNF, limx(2), loglimx(2)
    contains
        procedure, pass :: getParam, setParam, getLogPDF
    end type
 
    interface flatPoweto_type
        module procedure :: flatPoweto_typer
    end interface
 
contains
 
    function flatPoweto_typer(limx, param, domain) result(self)
        real(RKG), intent(in), contiguous, optional :: param(:)
        type(domain_type), intent(in), optional :: domain
        real(RKG), intent(in), contiguous :: limx(:)
        character(:), allocatable :: msg
        type(flatPoweto_type) :: self
        self%npar = 2
        self%limx = limx
        self%loglimx = log(limx)
        if (present(domain)) then
            self%domain = domain
        else
            self%domain = domain_type([-LARGE, -LARGE], [+LARGE, +LARGE], [character(63) :: "flatPoweto_logbreak", "flatPoweto_alpha"])
        end if
        if (present(param)) call self%setParam(param)
    end function
 
    function getParam(self) result(param)
        class(flatPoweto_type), intent(in) :: self
        real(RKG) :: param(self%npar)
        param = [self%logbreak, self%alpha]
    end function
 
    subroutine setParam(self, param)
        use pm_kind, only: SK, IK
        use pm_val2str, only: getStr
        use pm_mathMinMax, only: getMinMax
        use pm_mathLogAddExp, only: getLogAddExp
        use pm_mathLogSubExp, only: getLogSubExp
        real(RKG), intent(in), contiguous :: param(:)
        class(flatPoweto_type), intent(inout) :: self
        real(RKG) :: logModelInt1, logModelInt2, small, large
        if (size(param) /= self%npar) error stop "`lognorm_type%setParam()` takes a vector of two parameters representing `avg` and `logstd`."
        self%alpha = param(2)
        self%logbreak = param(1)
        self%alphap1 = self%alpha + 1
        self%break = exp(self%logbreak)
        logModelInt1 = log(self%break - self%limx(1))
        if (self%alphap1 < 0._RKG) then
            logModelInt2 = -log(-self%alphap1)
            large = self%alphap1 * self%logbreak
            small = self%alphap1 * self%loglimx(2)
        elseif (0._RKG < self%alphap1) then
            logModelInt2 = -log(self%alphap1)
            small = self%alphap1 * self%logbreak
            large = self%alphap1 * self%loglimx(2)
        else
            logModelInt2 = 0._RKG
            small = self%logbreak
            large = self%loglimx(2)
        end if
        logModelInt2 = logModelInt2 - self%alpha * self%logbreak + getLogSubExp(smaller = small, larger = large)
        self%logPDFNF = -getLogAddExp(getMinMax(logModelInt1, logModelInt2))
    end subroutine
 
    pure elemental function getLogPDF(self, obs) result(logPDF)
        use data_mod, only: obs_type, posuniv_type
        class(flatPoweto_type), intent(in) :: self
        class(obs_type), intent(in) :: obs
        real(RKG) :: logPDF
        select type (obs)
        type is (posuniv_type)
            if (obs%log < self%logbreak) then
                logPDF = obs%log
            else
                logPDF = self%alphap1 * obs%log - self%alpha * self%logbreak
            end if
            logPDF = logPDF + self%logPDFNF
        class default
            error stop "Unrecognized data."
        end select
    end function
 
end module flatPoweto_mod
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
module flatPowetoTapered_mod
 
    use auxil_mod, only: RKG, LARGE
    use model_mod, only: domain_type, model_type
    implicit none
 
    type, extends(model_type) :: flatPowetoTapered_type
        real(RKG), private  :: logbreak, alpha, beta, break, alphap1, logPDFNF, limx(2), loglimx(2)
    contains
        procedure, private :: getLogUDF
        procedure :: getParam, setParam, getLogPDF
    end type
 
    interface flatPowetoTapered_type
        module procedure :: flatPowetoTapered_typer
    end interface
 
contains
 
    function flatPowetoTapered_typer(limx, param, domain) result(self)
        real(RKG), intent(in), contiguous, optional :: param(:)
        type(domain_type), intent(in), optional :: domain
        real(RKG), intent(in), contiguous :: limx(:)
        type(flatPowetoTapered_type) :: self
        self%npar = 3
        self%limx = limx
        self%loglimx = log(limx)
        if (present(domain)) then
            self%domain = domain
        else
            self%domain = domain_type([-LARGE, -LARGE, -LARGE], [+LARGE, +LARGE, +LARGE], [character(63) :: "flatPowetoTapered_logbreak", "flatPowetoTapered_alpha", "flatPowetoTapered_beta"])
        end if
        if (present(param)) call self%setParam(param)
    end function
 
    function getParam(self) result(param)
        class(flatPowetoTapered_type), intent(in) :: self
        real(RKG) :: param(self%npar)
        param = [self%logbreak, self%alpha, self%beta]
    end function
 
    subroutine setParam(self, param)
        use pm_val2str, only: getStr
        use pm_kind, only: SK, IK, LK
        use pm_mathMinMax, only: getMinMax
        use pm_mathLogAddExp, only: getLogAddExp
        use pm_quadPack, only: isFailedQuad, getQuadErr, weps
        class(flatPowetoTapered_type), intent(inout) :: self
        real(RKG), intent(in), contiguous :: param(:)
        real(RKG) :: logModelInt1, logModelInt2, logNF
        character(255, SK) :: msg
        logical(LK) :: failed
        if (size(param) /= self%npar) error stop "`flatPowetoFlatPowetoTapered_type` takes a vector of three parameters"&
        //"representing `logbreak`, `alpha`, `beta`. size(param), npar = "//getStr([size(param, 1, IK), self%npar])
        self%logbreak = param(1)
        self%alpha = param(2)
        self%beta = param(3)
        self%alphap1 = self%alpha + 1
        self%break = exp(self%logbreak)
        logModelInt1 = log(self%break - self%limx(1))
        logNF = max(self%logbreak, self%alphap1 * self%loglimx(2) - self%alpha * self%logbreak - self%beta * (self%limx(2) - self%break))
        failed = isFailedQuad(getDensity, lb = self%logbreak, ub = self%loglimx(2), help = weps, integral = logModelInt2, msg = msg)
        if (failed) error stop "@getLogModelInt(): "//trim(msg)
        logModelInt2 = log(logModelInt2) + logNF
        self%logPDFNF = -getLogAddExp(getMinMax(logModelInt1, logModelInt2))
    contains
        function getDensity(logx) result(density)
            real(RKG), intent(in) :: logx
            real(RKG) :: density
            density = exp(self%getLogUDF(logx, exp(logx)) - logNF)
        end function
    end subroutine
 
    pure elemental function getLogPDF(self, obs) result(logPDF)
        use data_mod, only: obs_type, posuniv_type
        class(flatPowetoTapered_type), intent(in) :: self
        class(obs_type), intent(in):: obs
        real(RKG) :: logPDF
        select type (obs)
        type is (posuniv_type)
            if (obs%log < self%logbreak) then
                logPDF = obs%log
            else
                logPDF = self%getLogUDF(obs%log, obs%val)
            end if
            logPDF = logPDF + self%logPDFNF
        class default
            error stop "Unrecognized data."
        end select
    end function
 
    pure elemental function getLogUDF(self, logx, x) result(logUDF)
        class(flatPowetoTapered_type), intent(in) :: self
        real(RKG), intent(in) :: logx, x
        real(RKG) :: logUDF
        logUDF = self%alphap1 * logx - self%alpha * self%logbreak - self%beta * (x - self%break)
    end function
 
end module flatPowetoTapered_mod
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
module mixture_mod
 
    use pm_val2str, only: getStr
    use pm_arrayResize, only: setResized
    use pm_mathLogAddExp, only: getLogAddExp
    use pm_mathLogSumExp, only: getLogSumExp
    use model_mod, only: domain_type, model_type, con_type
    use auxil_mod, only: RKG, EPS, LARGE
    use data_mod, only: obs_type
    use pm_kind, only: SK, IK
    implicit none
 
    type, extends(model_type) :: mixture_type
        type(con_type), allocatable :: comp(:)
        real(RKG), allocatable :: logfrac(:)
        integer(IK), private :: ncomp
    contains
        procedure :: getParam, setParam, getLogPDF
    end type
 
    interface mixture_type
        module procedure :: mixture_typer
    end interface
 
contains
 
    function mixture_typer(comp, param) result(self)
        type(con_type), intent(in), contiguous :: comp(:)
        real(RKG), intent(in), contiguous, optional :: param(:)
        character(:, SK), allocatable :: prefix
        integer(IK) :: icomp, lpar, upar, ipar
        type(mixture_type) :: self
        self%ncomp = size(comp, 1, IK)
        allocate(self%comp, source = comp)
        call setResized(self%logfrac, self%ncomp)
        self%npar = sum([(self%comp(icomp)%model%npar, icomp = 1, self%ncomp)]) + self%ncomp - 1
        call setResized(self%domain%lower, self%npar)
        call setResized(self%domain%upper, self%npar)
        call setResized(self%domain%names, self%npar)
        lpar = 0
        do icomp = 1, self%ncomp - 1
            upar = lpar + self%comp(icomp)%model%npar
            prefix = SK_"comp"//getStr(icomp)//SK_"_"
            self%domain%lower(lpar + 1 : upar + 1) = [self%comp(icomp)%model%domain%lower, -LARGE]
            self%domain%upper(lpar + 1 : upar + 1) = [self%comp(icomp)%model%domain%upper, +LARGE]
            self%domain%names(lpar + 1 : upar + 1) = [character(63, SK) :: prefix//self%comp(icomp)%model%domain%names, prefix//SK_"fisherz(frac)"]
            lpar = upar + 1
        end do
        prefix = SK_"comp"//getStr(icomp)//SK_"_"
        self%domain%lower(lpar + 1 :) = self%comp(icomp)%model%domain%lower
        self%domain%upper(lpar + 1 :) = self%comp(icomp)%model%domain%upper
        self%domain%names(lpar + 1 :) = prefix//self%comp(icomp)%model%domain%names
        if (present(param)) call self%setParam(param)
    end function
 
    subroutine setParam(self, param)
        use pm_mathFisher, only: getFisherInv
        class(mixture_type), intent(inout) :: self
        real(RKG), intent(in), contiguous :: param(:)
        integer(IK) :: icomp, lpar, upar
        real(RKG) :: sumfrac, mixfrac
        if (size(param) /= self%npar) error stop "@mixture_type%setParam(): The condition `size(param) == self%npar` must hold. size(param), self%npar = "//getStr([size(param), self%npar])
        lpar = 0
        sumfrac = 0
        do icomp = 1, self%ncomp - 1
            upar = lpar + self%comp(icomp)%model%npar
            call self%comp(icomp)%model%setParam(param(lpar + 1 : upar))
            mixfrac = getFisherInv(param(upar + 1), 0._RKG, 1._RKG)
            self%logfrac(icomp) = log(mixfrac)
            sumfrac = sumfrac + mixfrac
            lpar = upar + 1
        end do
        self%logfrac(icomp) = log(1 - sumfrac)
        call self%comp(icomp)%model%setParam(param(lpar + 1 :))
        if (.not. abs(1 - sum(exp(self%logfrac))) < EPS) error stop "The condition `sum(exp(self%logfrac)) == 1` must hold. self%logfrac = "//getStr(self%logfrac)
    end subroutine
 
    function getParam(self) result(param)
        use pm_mathFisher, only: getFisher
        class(mixture_type), intent(in) :: self
        integer(IK) :: icomp, lpar, upar
        real(RKG) :: param(self%npar)
        param = [(self%comp(icomp)%model%getParam(), getFisher(exp(self%logfrac(icomp)), 0._RKG, 1._RKG), icomp = 1, self%ncomp - 1), self%comp(self%ncomp)%model%getParam()]
    end function
 
    impure elemental function getLogPDF(self, obs) result(logPDF)
        class(mixture_type), intent(in) :: self
        real(RKG) :: logPDFS(self%ncomp), maxLogPDF
        class(obs_type), intent(in) :: obs
        real(RKG) :: logPDF
        integer(IK) :: icomp
        maxLogPDF = -huge(0._RKG)
        do icomp = 1, self%ncomp
            logPDFS(icomp) = self%logfrac(icomp) + self%comp(icomp)%model%getLogPDF(obs)
            maxLogPDF = max(maxLogPDF, logPDFS(icomp))
        end do
        logPDF = getLogSumExp(logPDFS, maxLogPDF)
    end function
 
end module mixture_mod
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
module fit_mod
 
    use auxil_mod, only: RKG
    use pm_sampling, only: sampler_type
    use pm_sampling, only: paradram_type
    use model_mod, only: model_type
    use data_mod, only: data_type
    implicit none
 
    private
    public :: fit_type, paradram_type
 
    type :: best_type
        real(RKG) :: loglike
        real(RKG), allocatable :: param(:) ! best fit model parameters.
    end type
 
    type :: fit_type
        real(RKG) :: bic
        type(best_type) :: best
        type(data_type) :: data
        class(model_type), allocatable :: model
        class(sampler_type), allocatable :: sampler
    contains
        procedure, pass :: run
        !procedure, pass :: setHist
    end type
 
    interface fit_type
        module procedure :: fit_typer
    end interface
 
contains
 
    function fit_typer(data, model, sampler) result(self)
        class(sampler_type), intent(in), optional :: sampler
        class(model_type), intent(in) :: model
        type(data_type), intent(in) :: data
        type(fit_type) :: self
        if (present(sampler)) then
            self%sampler = sampler
        else
            self%sampler = paradram_type()
        end if
        self%model = model
        self%data = data
    end function
 
    subroutine run(self)
 
        use pm_kind, only: SK, IK
        use pm_err, only: err_type
        use pm_sampling, only: getErrSampling, paradram_type
        class(fit_type), intent(inout) :: self
        type(err_type) :: err
 
        if (.not. allocated(self%sampler%outputStatus)) self%sampler%outputStatus = "retry"
        if (.not. allocated(self%sampler%domainAxisName)) self%sampler%domainAxisName = self%model%domain%names
        if (.not. allocated(self%sampler%domainCubeLimitLower)) self%sampler%domainCubeLimitLower = self%model%domain%lower
        if (.not. allocated(self%sampler%domainCubeLimitUpper)) self%sampler%domainCubeLimitUpper = self%model%domain%upper
 
        select type (sampler => self%sampler)
        type is (paradram_type)
            if (.not. allocated(sampler%outputChainSize)) sampler%outputChainSize = 30000
            if (.not. allocated(sampler%proposalScale)) sampler%proposalScale = "gelman"
            if (.not. allocated(sampler%proposalStart)) sampler%proposalStart = self%model%getParam()
            err = getErrSampling(sampler, getLogLike, self%model%npar)
        end select
        if (err%occurred) error stop "sampler failed: "//err%msg
 
        ! Find the mean best fit parameters and write post-prediction data to output file for visualization.
 
        block
 
            use pm_io, only: getErrTableRead
            use pm_sysPath, only: glob, css_type
 
            integer(IK) :: stat, ibest!, offset = 1
            type(css_type), allocatable :: path(:)
            real(RKG), allocatable :: logx(:), logPDF(:)
            real(RKG), allocatable :: table(:,:), state(:)
 
            write(*, "(A)") "Searching for files: "//self%sampler%outputFileName//SK_"*"
            path = glob(self%sampler%outputFileName//SK_"*")
            if (size(path) == 0) error stop "There is no sample file in the output folder."
 
            stat = getErrTableRead(path(size(path))%val, table, roff = 1_IK)
            if (stat /= 0) error stop "Failed to read output sample."
            ibest = maxloc(table(:, 1), dim = 1_IK)
            self%best%loglike = table(ibest, 1)
            self%best%param = table(ibest, 2:)
            self%bic = self%model%npar * log(real(self%data%nobs, RKG)) - 2 * self%best%loglike
 
        end block
 
    contains
 
        function getLogLike(param) result(logLike)
            real(RKG), intent(in), contiguous :: param(:)
            real(RKG) :: logLike
            call self%model%setParam(param)
            loglike = sum(self%model%getLogPDF(self%data%obs))
        end function
 
    end subroutine
 
end module fit_mod
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
program example
 
    use pm_kind, only: SK, IK
    use fit_mod, only: fit_type
    use model_mod, only: con_type
    use auxil_mod, only: RKG, LARGE
    use lognorm_mod, only: lognorm_type
    use mixture_mod, only: mixture_type
    use pm_sampling, only: paradram_type
    use pm_arraySpace, only: getLogSpace
    use flatPoweto_mod, only: flatPoweto_type
    use data_mod, only: posuniv_type, data_type, stat_type
    use flatPowetoTapered_mod, only: flatPowetoTapered_type
    use pm_io, only: display_type, css_type, getErrTableWrite
 
    integer(IK) :: imodel, stat
    type(paradram_type) :: sampler
    type(data_type) :: data, pred
    type(display_type) :: disp
    type(fit_type) :: fit
 
    disp = display_type(file = SK_"main.out.F90")
 
    ! read data.
 
    block
        use pm_io, only: getErrTableRead
        real(RKG), allocatable :: table(:,:)
        character(255, SK) :: iomsg
        stat = getErrTableRead("data.csv", table, roff = 1_IK, iomsg = iomsg)
        if (stat /= 0) error stop "Failed to read data: "//trim(iomsg)
        associate(val => table(:, 4))
            data = data_type(obs = posuniv_type(val), stat = stat_type([minval(val), maxval(val)]))
        end associate
    end block
 
    ! perform fitting.
 
    do imodel = 1, 4
 
        sampler = paradram_type()
 
        if (imodel == 1) then
            sampler%outputFileName = "./mixLogNormLogNorm"
            sampler%proposalStart = [-.2, 0.4, .4, 3.6, -0.2]
            sampler%domainCubeLimitLower = [real(RKG) :: -LARGE, -LARGE, -LARGE, log(3.), -LARGE]
            sampler%domainCubeLimitUpper = [real(RKG) :: log(3.), +LARGE, +LARGE, +LARGE, +LARGE]
            !sampler%proposalScale = "0.1 * gelman"
            sampler%proposalCov = reshape   ( &
                                            [  1.6E-2, 5.5E-3, 4.1E-3, 3.9E-3,-2.5E-3 &
                                            ,  5.5E-3, 3.1E-3, 1.7E-3, 1.4E-3,-9.5E-4 &
                                            ,  4.1E-3, 1.7E-3, 1.9E-3, 1.2E-3,-8.8E-4 &
                                            ,  3.9E-3, 1.4E-3, 1.2E-3, 2.0E-3,-8.9E-4 &
                                            , -2.5E-3,-9.5E-4,-8.8E-4,-8.9E-4, 1.0E-3 &
                                            ], shape = [size(sampler%proposalStart), size(sampler%proposalStart)])
            fit = fit_type(data, mixture_type([con_type(lognorm_type()), con_type(lognorm_type())]), sampler)
        elseif (imodel == 2) then
            sampler%proposalScale = "0.1 * gelman"
            sampler%outputFileName = "./mixFlatPowetoFlatPowetoTapered"
            sampler%proposalStart = [log(.37), -1.4, .3, log(21.3), -.5, 0.0156]
            sampler%proposalCov = reshape   ( &
                                            [  1.E-2, -3.E-3,  1.E-2, -7.E-3, -1.E-4, -6.E-4 &
                                            , -3.E-3,  3.E-3, -4.E-3,  6.E-3,  1.E-4,  1.E-3 &
                                            ,  1.E-2, -4.E-3,  1.E-1, -7.E-2, -8.E-4,  3.E-4 &
                                            , -7.E-3,  6.E-3, -7.E-2,  7.E-2,  8.E-4,  3.E-3 &
                                            , -1.E-4,  1.E-4, -8.E-4,  8.E-4,  1.E-5,  6.E-5 &
                                            , -6.E-4,  1.E-3,  3.E-4,  3.E-3,  6.E-5,  2.E-3 &
                                            ], shape = [size(sampler%proposalStart), size(sampler%proposalStart)])
            fit = fit_type(data, mixture_type([con_type(flatPoweto_type(data%stat%lim)), con_type(flatPowetoTapered_type(data%stat%lim))]), sampler)
        elseif (imodel == 3) then
            sampler%proposalScale = "0.1 * gelman"
            sampler%outputFileName = "./mixLogNormFlatPowetoTapered"
            sampler%proposalStart = [-.2, 0.4, .3, log(21.3), -.5, 0.0156]
            fit = fit_type(data, mixture_type([con_type(lognorm_type()), con_type(flatPowetoTapered_type(data%stat%lim))]), sampler)
        elseif (imodel == 4) then
            sampler%proposalScale = "0.1 * gelman"
            sampler%outputFileName = "./mixFlatPowetoLogNorm"
            sampler%proposalStart = [log(.37), -1.4, .3, 3.6, -0.2]
            fit = fit_type(data, mixture_type([con_type(flatPoweto_type(data%stat%lim)), con_type(lognorm_type())]), sampler)
        else
            sampler%outputFileName = "./logNorm"
            fit = fit_type(data, lognorm_type([real(RKG) :: 1, 1]), sampler)
        end if
 
        call fit%run()
 
        call disp%show(css_type(fit%sampler%domainAxisName))
        call disp%show(fit%best%param)
        call disp%skip()
        call disp%show("[fit%best%loglike, fit%bic]")
        call disp%show( [fit%best%loglike, fit%bic] )
        call disp%skip()
 
        ! generate synthetic data for histogram reconstruction.
 
        associate(val => getLogSpace(logx1 = log(data%stat%lim(1)), logx2 = log(data%stat%lim(2)), count = 1000_IK))
            pred = data_type(obs = posuniv_type(val), stat = stat_type([minval(val), maxval(val)]))
            call fit%model%setParam(fit%best%param)
            stat = getErrTableWrite(fit%sampler%outputFileName//SK_".hist", reshape([log(val), fit%model%getLogPDF(pred%obs)], [size(val), 2]), header = SK_"logx,logPDF")
            if (stat /= 0) error stop "Failed to write the histogram visualization data."
        end associate
 
    end do
 
end

Example Unix compile command via Intel ifort compiler ⛓

#!/usr/bin/env sh
rm main.exe
ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
./main.exe

Example Windows Batch compile command via Intel ifort compiler ⛓

del main.exe
set PATH=..\..\..\lib;%PATH%
ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
main.exe

Example Unix / MinGW compile command via GNU gfortran compiler ⛓

#!/usr/bin/env sh
rm main.exe
gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
./main.exe

Postprocessing of the example output ⛓

#!/usr/bin/env python
 
import matplotlib.pyplot as plt
#import paramonte as pm
import pandas as pd
import numpy as np
import glob
import sys
import os
 
plt.rcParams['text.usetex'] = True
# Uncomment above line if LaTeX is installed on the system.
# On Ubuntu, this requires: sudo apt install dvipng texlive-latex-extra texlive-fonts-recommended cm-super
examname = os.path.basename(os.getcwd())
 
 
 
linewidth = 2
fontsize = 17
files = glob.glob("./*_chain.txt")
 
for file in files:
 
    basename = file.split("_")[0]
    df = pd.read_csv(file, delimiter = ",")
    if "_chain.txt" in file:
        sindex = 7 # start column index.
    elif "_sample.txt" in file:
        sindex = 1 # start column index.
    else:
        sys.exit("Unrecognized simulation output file: " + file)
 
    # traceplot
 
    #print(df.values)
    fig = plt.figure(figsize = (8, 6))
    ax = plt.subplot(1,1,1)
    ax.plot ( range(len(df.values[:, 0]))
            , df.values[:, sindex:]
            , zorder = 1000
            )
    plt.minorticks_on()
    ax.set_xscale("log")
    ax.set_xlabel("MCMC Count", fontsize = 17)
    ax.set_ylabel("MCMC State", fontsize = 17)
    ax.tick_params(axis = "x", which = "major", labelsize = fontsize)
    ax.tick_params(axis = "y", which = "major", labelsize = fontsize)
    plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
    #ax.legend(df.columns.values, fontsize = fontsize)
    #plt.title(basename)
    plt.tight_layout()
    plt.savefig(basename + ".traceplot.png")
 
    # scatterplot
 
    if len(df.values[1, sindex:]) == 2:
        #print(df.values)
        fig = plt.figure(figsize = (8, 6))
        ax = plt.subplot(1,1,1)
        ax.scatter  ( df.values[:, sindex]
                    , df.values[:, sindex + 1]
                    , zorder = 1000
                    , s = 1
                    )
        plt.minorticks_on()
        ax.set_xscale("log")
        ax.set_xlabel(df.columns.values[sindex], fontsize = 17)
        ax.set_ylabel(df.columns.values[sindex + 1], fontsize = 17)
        ax.tick_params(axis = "x", which = "major", labelsize = fontsize)
        ax.tick_params(axis = "y", which = "major", labelsize = fontsize)
        plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
        #ax.legend(df.columns.values, fontsize = fontsize)
        #plt.title(basename)
        plt.tight_layout()
        plt.savefig(basename + ".scatterplot.png")
 
    # adaptation measure.
 
    #print(df.values)
    if "proposalAdaptation" in df.columns.values:
        if any(df["proposalAdaptation"].values != 0):
            fig = plt.figure(figsize = (8, 6))
            ax = plt.subplot(1,1,1)
            ax.scatter  ( range(len(df.values[:, 0]))
                        , df["proposalAdaptation"].values
                        , zorder = 1000
                        , c = "red"
                        , s = 1
                        )
            plt.minorticks_on()
            #ax.set_xscale("log")
            ax.set_yscale("log")
            ax.set_xlabel("MCMC Count", fontsize = 17)
            ax.tick_params(axis = "x", which = "major", labelsize = fontsize)
            ax.tick_params(axis = "y", which = "major", labelsize = fontsize)
            ax.set_ylabel("proposalAdaptation", fontsize = 17)
            plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
            #ax.legend(df.columns.values, fontsize = fontsize)
            #plt.title(basename)
            plt.tight_layout()
            plt.savefig(basename + ".proposalAdaptation.png")
 
    #sim = pm.ParaDRAM()
    #sample = sim.readSample("./out/", renabled = True)
    #sample[0].plot.grid()
 
 
 
# from Functions import mergeBins, findPlateauLocation
 
def mergeBins(counts, bins, threshold = 5):
 
    if len(bins)%2 == 0:
        halfway = int(np.ceil(len(bins)/2))
    else:
        halfway = int(np.floor(len(bins)/2))
 
    threshold = threshold
    leftCounts = []
    leftBins = [bins[0]]
    rightCounts = []
    rightBins = [bins[-1]]
    leftSumCounts = 0
    rightSumCounts = 0
    leftBinstart = 0
    rightBinstart = 0
    for i in range(0 , halfway): #goes to the middle of the distribution
       # print(counts[i],sumCounts)
 
        leftSumCounts += counts[i]
 
 
 
        if leftSumCounts >= threshold:
           # print("added value to list")
            leftCounts += [leftSumCounts]
            leftSumCounts = 0
            leftBins += [bins[i+1]]
 
        #if counts[-i-1]==3:print(rightSumCounts,-i-1)
 
 
        rightSumCounts += counts[-i-1]
        if rightSumCounts >= threshold:
            rightCounts += [rightSumCounts]
            rightSumCounts = 0
            rightBins += [bins[-i-2]]
 
    if len(bins)%2==0:
        if len(rightCounts) != 0:
            rightCounts.pop()
            rightBins.pop()
        leftBins.pop()
    else:
        rightCounts += [counts[halfway]]
        #rightBins += [bins[halfway]]
 
    newCounts = leftCounts + rightCounts[-1::-1]
    newBins = leftBins + rightBins[-1::-1]
    #print(len(newCounts),len(newBins))
    origionaldNdT = np.array(counts) / np.array([bins[i + 1] - bins[i] for i in range(0, len(bins) - 1) ])
    newdNdT = np.array(newCounts) / np.array([newBins[i + 1] - newBins[i] for i in range(0, len(newBins) - 1) ])
 
    return np.array(newCounts), np.array(newBins), np.array(origionaldNdT), np.array(newdNdT)
 
# Read the duration data.
 
data = pd.read_csv("data.csv", delimiter = ",")
durs = data["T90"].values
 
nbins = 50
binning = 10**np.linspace(np.log10(min(durs)), np.log10(max(durs)), nbins)
countsAndBins = np.histogram(durs, bins = binning)
counts = countsAndBins[0]
bins = countsAndBins[1]
bins = bins.tolist()
newCounts, newBins, origionaldNdT, dNdT = mergeBins(counts, binning)
 
# get the middle of the bin edges for location of uncertainties.
#midBins = [newBins[i] + (newBins[i + 1] - newBins[i]) / 2 for i in range(0, len(newBins) - 1)]
# use algorithm to find the longest plateau within one std of right most point
#plateau, binEdge, plateauErrors = findPlateauLocation(dNdT, dNdTerrors, newBins)
# find the best value for the plateau
#bestParams = opt.minimize(chisqfunc, [plateau[-1]], args = (plateau, plateauErrors))
#x = np.logspace(np.log10(min(binEdge)),np.log10(max(binEdge)),len(midBins))
#plt.plot(x,(bestParams["x"]).tolist()*len(x),color = 'k',lw = 2)
#x = np.logspace(np.log10(min(binEdge)), np.log10(max(binEdge)), len(midBins))
 
histFiles = glob.glob("./*.hist")
 
for histFile in histFiles:
 
    fig = plt.figure()
    plt.stairs  ( dNdT
                , newBins
                , linewidth = 2
                , baseline = None
                , color = "green"
                , label = 'BATSE'
                )
    plt.xscale('log')
    plt.yscale('log')
 
    # Read the fitting data to plot.
 
    model = pd.read_csv(histFile, delimiter = ",")
    x = np.exp(model["logx"].values)
    y = len(durs) * np.exp(model["logPDF"].values) / x
    plt.grid(which = "both")
    plt.plot(x, y, c = "magenta", label = histFile)
    plt.xlabel(r"$T_{90} ~ [s]$", fontsize = fontsize)
    plt.ylabel(r"$dN ~ / ~ dT_{90} ~ [1 / s]$", fontsize = fontsize)
    #plt.legend(fontsize = 13)
    plt.xticks(size = fontsize)
    plt.yticks(size = fontsize)
    plt.ylim(7*10**-3, 2*10**3)
    plt.tight_layout()
    plt.savefig(histFile + ".png")
    #plt.show()

Visualization of the example output ⛓

Test:

test_pm_sampling

Bug:

Status: Unresolved
Source: Intel LLVM Fortran Compiler ifx version 2024.0.2 20231213
Description: Intel LLVM Fortran Compiler ifx cannot handle compilation of the pm_sampling implementation include file.

/home/amir/git/paramonte/src/fortran/main/pm_sampling@routines.inc.F90(454): error #5623: **Internal compiler error: internal abort** Please report this error along with the circumstances in which it occurred in a Software Problem Report.  Note: File and line given may not be explicit cause of this error.
                   if (spec%parallelism%is%forkJoin .and. scalings%current%scaling(spec%image%count) < 1._RKG) then
------------------------------------------------------------------------------------------------------^
compilation aborted for /home/amir/git/paramonte/src/fortran/main/pm_sampling@routines.F90 (code 3)

Remedy (as of ParaMonte Library version 2.0.0): The new interface allocatable output arguments entirely, thus obviating the need to handle this gracefully.

Final Remarks ⛓

---

If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

1. If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
   
2. If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.
   

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Copyright

Computational Data Science Lab

Author:

Amir Shahmoradi, September 1, 2012, 12:00 AM, National Institute for Fusion Studies, The University of Texas Austin

Definition at line 674 of file pm_sampling.F90.

---

The documentation for this interface was generated from the following file:

• src/fortran/main/pm_sampling.F90