ParaMonte Fortran 2.0.0
Parallel Monte Carlo and Machine Learning Library
See the latest version documentation.
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,
  1. type paradram_type indicating the use of the ParaDRAM sampler,
  2. type paradise_type indicating the use of the ParaDISE sampler,
  3. 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:
  1. 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,
    1. 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.
    2. the constant RKD refers to the real kind double precision supported by the compiler.
    3. the output mold merely dictates the real kind of the interface and its value is ignored.
    4. 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).
    5. the two extra intent(inout) arguments avgTimePerFunCall, avgCommPerFunCall represent:
      1. the average time (in seconds) it takes to compute the user-specified function for a single input state.
      2. 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.
  2. 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,
    1. the variable ndim is the number of dimensions of the target density function.
    2. the variable state is an input point from within the domain of the density function.
    3. the constant RKG refers to the real kind used for computations.
    4. 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)
pm_sampling::getErrSampling
Generate and return .true. if the procedure fails to fully accomplish the task of Monte Carlo samplin...
Definition: pm_sampling.F90:674
pm_sampling
This module contains procedures and generic interfaces for the ParaMonte library sampler routines.
Definition: pm_sampling.F90:33
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 ⛓

1program example
2
3 use pm_kind, only: RKG => RK
4 use pm_kind, only: SK, IK, LK
5 use pm_io, only: display_type
6 use pm_sampling, only: getErrSampling, paradram_type
7 use pm_matrixInit, only: getMatInit, uppLowDia
8 use pm_arrayFill, only: getFilled
9 use pm_sysPath, only: glob
10 use pm_err, only: err_type
11
12 implicit none
13
14 type(err_type) :: err
15 type(display_type) :: disp
16 integer(IK), parameter :: NDIM = 4
17 real(RKG), parameter :: MEAN(*) = 4 * [real(RKG) :: -3, -1, 1, 3]
18
19 disp = display_type(file = SK_"main.out.F90")
20
21 call disp%skip()
22 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
23 call disp%show("! Simple minimally-informed simulation.")
24 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
25 call disp%skip()
26
27 block
28 call disp%skip()
29 call disp%show("NDIM")
30 call disp%show( NDIM )
31 call disp%show("err = getErrSampling(paradram_type(parallelismMpiFinalizeEnabled = .false._LK), getLogFunc, NDIM)")
32 err = getErrSampling(paradram_type(parallelismMpiFinalizeEnabled = .false._LK), getLogFunc, NDIM)
33 call disp%show("if (err%occurred) error stop err%msg")
34 if (err%occurred) error stop err%msg
35 call disp%skip()
36 end block
37
38 call disp%skip()
39 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
40 call disp%show("! Optionally-informed simulation.")
41 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
42 call disp%skip()
43
44 block
45 type(paradram_type) :: sampler
46 call disp%skip()
47 call disp%show("NDIM")
48 call disp%show( NDIM )
49 ! GNU gfortran yields a segfault with the following default constructor.
50 !call disp%show("sampler = paradram_type(parallelismMpiFinalizeEnabled = .false._LK, outputFileName = 'mvn', outputChainSize = 100000, outputStatus = 'retry')")
51 ! sampler = paradram_type(parallelismMpiFinalizeEnabled = .false._LK, outputFileName = 'mvn', outputChainSize = 100000, outputStatus = 'retry')
52 call disp%show("sampler = paradram_type()")
53 sampler = paradram_type()
54 call disp%show("sampler%parallelismMpiFinalizeEnabled = .false._LK; sampler%outputFileName = 'mvn'; sampler%outputChainSize = 100000; sampler%outputStatus = 'retry'")
55 sampler%parallelismMpiFinalizeEnabled = .false._LK; sampler%outputFileName = 'mvn'; sampler%outputChainSize = 100000; sampler%outputStatus = 'retry'
56 call disp%show("err = getErrSampling(sampler, getLogFunc, NDIM)")
57 err = getErrSampling(sampler, getLogFunc, NDIM)
58 call disp%show("if (err%occurred) error stop err%msg")
59 if (err%occurred) error stop err%msg
60 call disp%show("glob(pattern = sampler%outputFileName)")
61 call disp%show( glob(pattern = sampler%outputFileName) )
62 call disp%skip()
63 end block
64
65contains
66
67 recursive function getLogFunc(state) result(logFunc)
68 use pm_distMultiNorm, only: getMultiNormLogPDF
69 real(RKG), intent(in), contiguous :: state(:)
70 real(RKG) :: logFunc
71 integer :: i
72 logFunc = getMultiNormLogPDF(state, mean = real(MEAN, RKG))
73#if OMP_ENABLED
74 ! We are done. But kill some time for an illustration of parallel sampling.
75 do i = 1, 500
76 logFunc = logFunc + getMultiNormLogPDF(state, mean = 4 * [real(RKG) :: -3, -1, 1, 3]) * merge(1, -1, mod(i, 2) == 0)
77 end do
78#endif
79 end function
80
81 !subroutine setLogFunc(logFuncState)
82 ! use pm_distMultiNorm, only: getMultiNormLogPDF
83 ! real(RKG), intent(inout), contiguous :: logFuncState(0:,:)
84 ! integer(IK) :: ithread
85 ! do ithread = 1, size(logFuncState, 2)
86 ! logFuncState(0, ithread) = getMultiNormLogPDF(logFuncState(1:, ithread), mean = real(MEAN, RKG))
87 ! end do
88 !end subroutine
89
90end program example
pm_arrayFill::getFilled
Generate and return an array of the specified rank and shape of arbitrary intrinsic type and kind wit...
Definition: pm_arrayFill.F90:132
pm_distMultiNorm::getMultiNormLogPDF
Generate and return the natural logarithm of the Probability Density Function (PDF) of the MultiVaria...
Definition: pm_distMultiNorm.F90:449
pm_err::show
Generate and return an object of type stop_type with the user-specified input attributes.
Definition: pm_err.F90:1618
pm_io::skip
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11508
pm_matrixInit::getMatInit
Generate and return a matrix of shape (shape(1), shape(2)) with the upper/lower triangle and the diag...
Definition: pm_matrixInit.F90:289
pm_sysPath::glob
Generate and return the result of globing the specified input pattern.
Definition: pm_sysPath.F90:954
pm_arrayFill
This module contains procedures and generic interfaces for convenient allocation and filling of array...
Definition: pm_arrayFill.F90:42
pm_distMultiNorm
This module contains classes and procedures for computing various statistical quantities related to t...
Definition: pm_distMultiNorm.F90:60
pm_err
This module contains classes and procedures for reporting and handling errors.
Definition: pm_err.F90:52
pm_io
This module contains classes and procedures for input/output (IO) or generic display operations on st...
Definition: pm_io.F90:252
pm_io::disp
type(display_type) disp
This is a scalar module variable an object of type display_type for general display.
Definition: pm_io.F90:11393
pm_kind
This module defines the relevant Fortran kind type-parameters frequently used in the ParaMonte librar...
Definition: pm_kind.F90:268
pm_kind::RK
integer, parameter RK
The default real kind in the ParaMonte library: real64 in Fortran, c_double in C-Fortran Interoperati...
Definition: pm_kind.F90:543
pm_kind::LK
integer, parameter LK
The default logical kind in the ParaMonte library: kind(.true.) in Fortran, kind(....
Definition: pm_kind.F90:541
pm_kind::IK
integer, parameter IK
The default integer kind in the ParaMonte library: int32 in Fortran, c_int32_t in C-Fortran Interoper...
Definition: pm_kind.F90:540
pm_kind::SK
integer, parameter SK
The default character kind in the ParaMonte library: kind("a") in Fortran, c_char in C-Fortran Intero...
Definition: pm_kind.F90:539
pm_matrixInit
This module contains procedures and generic interfaces relevant to generating and initializing matric...
Definition: pm_matrixInit.F90:116
pm_sysPath
This module contains classes and procedures for manipulating system file/folder paths.
Definition: pm_sysPath.F90:274
pm_err::err_type
This is the derived type for generating objects to gracefully and verbosely handle runtime unexpected...
Definition: pm_err.F90:157
pm_io::display_type
Generate and return an object of type display_type.
Definition: pm_io.F90:10282
pm_sampling::paradram_type
This is a derived type for constructing objects containing the optional simulation properties of the ...
Definition: pm_sampling.F90:250

Example Unix compile command via Intel ifort compiler ⛓
1#!/usr/bin/env sh
2rm main.exe
3ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example Windows Batch compile command via Intel ifort compiler ⛓
1del main.exe
2set PATH=..\..\..\lib;%PATH%
3ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
4main.exe

Example Unix / MinGW compile command via GNU gfortran compiler ⛓
1#!/usr/bin/env sh
2rm main.exe
3gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example output ⛓
1
2!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3! Simple minimally-informed simulation.
4!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5
6
7NDIM
8+4
9err = getErrSampling(paradram_type(parallelismMpiFinalizeEnabled = .false._LK), getLogFunc, NDIM)
10if (err%occurred) error stop err%msg
11
12
13!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14! Optionally-informed simulation.
15!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
16
17
18NDIM
19+4
20sampler = paradram_type()
21sampler%parallelismMpiFinalizeEnabled = .false._LK; sampler%outputFileName = 'mvn'; sampler%outputChainSize = 100000; sampler%outputStatus = 'retry'
22err = getErrSampling(sampler, getLogFunc, NDIM)
23if (err%occurred) error stop err%msg
24glob(pattern = sampler%outputFileName)
25
26
27


Example usage ⛓

1module logfunc
2 use pm_kind, only: IK, RKG => RKS ! All real kinds are supported.
3 implicit none
4 integer(IK) , parameter :: NDIM = 4_IK
5 real(RKG) , parameter :: MEAN(NDIM) = [-6, -2, 2, 6] ! The mean vector of the MVN
6 real(RKG) , parameter :: LOG_INVERSE_SQRT_TWO_PI = log(1./sqrt(2.*acos(-1.))) ! MVN distribution coefficient: log(1/sqrt(2*Pi)^ndim).
7 real(RKG) , parameter :: COVMAT(NDIM, NDIM) = reshape([ +1.0, +0.5, +0.5, +0.5 &
8 , +0.5, +1.0, +0.5, +0.5 &
9 , +0.5, +0.5, +1.0, +0.5 &
10 , +0.5, +0.5, +0.5, +1.0 ], shape(COVMAT)) ! The covariance matrix of the MVN
11 real(RKG) , parameter :: INVCOV(NDIM, NDIM) = reshape([ +1.6, -0.4, -0.4, -0.4 &
12 , -0.4, +1.6, -0.4, -0.4 &
13 , -0.4, -0.4, +1.6, -0.4 &
14 , -0.4, -0.4, -0.4, +1.6 ], shape(INVCOV)) ! The inverse covariance matrix of the MVN
15 real(RKG) , parameter :: MVN_COEF = NDIM * LOG_INVERSE_SQRT_TWO_PI + 0.581575404902840 ! log(sqrt(det(INVCOV)))
16contains
17 function getLogFunc(state) result(logFunc)
18 ! Return the negative natural logarithm of MVN distribution
19 ! evaluated at the input vector `state` of length `ndim`.
20 real(RKG), intent(in), contiguous :: state(:)
21 real(RKG) :: stateNormed(size(state))
22 real(RKG) :: logFunc
23 stateNormed = state - MEAN
24 logFunc = MVN_COEF - 0.5_RKG * (dot_product(stateNormed, matmul(INVCOV, stateNormed)))
25 end function
26end module logfunc
27
28program example
29 use logfunc, only: RKG, NDIM, getLogFunc
30 use pm_sampling, only: getErrSampling, paradram_type, err_type
31 type(err_type) :: err
32 err = getErrSampling(paradram_type(inputFile = 'input.nml'), getLogFunc, NDIM)
33 if (err%occurred) error stop "sampler failed: "//err%msg
34end
pm_kind::RKS
integer, parameter RKS
The single-precision real kind in Fortran mode. On most platforms, this is an 32-bit real kind.
Definition: pm_kind.F90:567

The contents of example specifications namelist input file input.nml ⛓
1! DESCRIPTION:
2!
3! The external input file for sampling the 4-dimensional Multivariate Normal distribution function as implemented in the accompanying source files.
4! This file is common between all supported programming language environments.
5!
6! NOTE:
7!
8! All simulation specifications (including this whole file) are optional and can be nearly safely commented out.
9! However, if domain boundaries are finite, set them explicitly.
10!
11! USAGE:
12!
13! -- Comments must begin with an exclamation mark `!`.
14! -- Comments can appear anywhere on an empty line or, after a variable assignment
15! (but not in the middle of a variable assignment whether in single or multiple lines).
16! -- All variable assignments are optional and can be commented out. In such cases, appropriate default values will be assigned.
17! -- Use ParaDRAM namelist (group) name to group a set of ParaDRAM simulation specification variables.
18! -- The order of the input variables in the namelist groups is irrelevant and unimportant.
19! -- Variables can be defined multiple times, but only the last definition will be considered as input.
20! -- All variable names are case insensitive. However, for clarity, this software follows the camelCase code-writing practice.
21! -- String values must be enclosed with either single or double quotation marks.
22! -- 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.
23! -- All vectors and arrays in the input file begin with index 1. This is following the convention of
24! the majority of science-oriented programming languages: Fortran, Julia, Mathematica, MATLAB, and R.
25!
26! For comprehensive guidelines on the input file organization and rules, visit:
27!
28! https://www.cdslab.org/paramonte/generic/latest/usage/sampling/paradram/input/
29!
30! To see detailed descriptions of each of variables, visit:
31!
32! https://www.cdslab.org/paramonte/generic/latest/usage/sampling/paradram/specifications/
33!
34&paradram
35
36 ! Base specifications.
37
38 description = "
39This\n
40 is a\n
41 multi-line\n
42 description.\\n" ! strings must be enclosed with "" or '' and can be continued on multiple lines.
43 ! No comments within strings are allowed.
44 domain = "cube"
45 domainAxisName = "variable1"
46 "variable2" ! values can appear in multiple lines.
47 domainBallAvg = 0 0 0 0 ! values can be separated with blanks or commas.
48 domainBallCor = 1 0 0 0
49 0 1 0 0
50 0 0 1 0
51 0 0 0 1
52 domainBallCov = 1 0 0 0
53 0 1 0 0
54 0 0 1 0
55 0 0 0 1
56 domainBallStd = 1 1 1 1
57 domainCubeLimitLower = 4*-1.e10 ! repetition pattern rules apply here. 4 dimensions => 4-element vector of values.
58 domainCubeLimitUpper(1) = +1.e10 ! Elements of vectors can be set individually.
59 domainCubeLimitUpper(2:4) = 3*+1.e10 ! Elements of vectors can be set individually.
60 domainErrCount = 100
61 domainErrCountMax = 1000
62 inputFileHasPriority = FALSE ! This is relevant only to simulations within the Fortran programming language.
63 outputChainFileFormat = "compact"
64 !outputColumnWidth = 25 ! This is an example of a variable that is commented out.
65 ! Therefore, its value will not be read by the sampler routine.
66 ! To pass it to the routine, simply remove the `!` mark at the beginning of the line.
67 outputFileName = "./out/mvn" ! A forward-slash character at the end of the string value would indicate the specified path
68 ! is to be interpreted as the name of the folder to contain the simulation output files.
69 ! The base name for the simulation output files will be generated from the current date and time.
70 ! Otherwise, the specified base name at the end of the string will be used in naming the simulation output files.
71 outputPrecision = 17
72 outputReportPeriod = 1000
73 outputRestartFileFormat = "ascii"
74 outputSampleSize = -1
75 outputSeparator = ","
76 outputSplashMode = "normal" ! or quiet or silent.
77 outputStatus = "retry" ! or extend.
78 parallelism = "multi chain" ! "singleChain" would also work. Similarly, "multichain", "multi chain", or "multiChain".
79 parallelismMpiFinalizeEnabled = false ! TRUE, true, .true., .t., and t would be also all valid logical values representing truth.
80 parallelismNumThread = 3 ! number of threads to use in shared-memory parallelism.
81 !randomSeed = 2136275,
82 !targetAcceptanceRate = 0.23e0
83
84 ! MCMC specifications.
85
86 outputChainSize = 10000
87 outputSampleRefinementCount = 10
88 outputSampleRefinementMethod = "BatchMeans"
89 proposal = "normal" ! or "uniform" as you wish.
90 proposalCor(:, 1) = 1 0 0 0 ! first matrix column.
91 proposalCor(:, 2:4) = 0 1 0 0
92 0 0 1 0
93 0 0 0 1 ! other matrix columns.
94 proposalCov = 1 0 0 0
95 0 1 0 0
96 0 0 1 0
97 0 0 0 1 ! or specify all matrix elements in one statement.
98 proposalScale = "2*0.5*Gelman" ! The asterisk here means multiplication since it is enclosed within quotation marks.
99 !proposalStart = 4*1.e0 ! four values of 1.e0 are specified here by the repetition pattern symbol *
100 proposalStartDomainCubeLimitLower = 4*-10.e0 ! repetition pattern rules apply again here. 4 dimensions => 4-element vector of values.
101 proposalStartDomainCubeLimitUpper = 4*+10.e0 ! repetition pattern rules apply again here. 4 dimensions => 4-element vector of values.
102 proposalStartRandomized = false
103 proposalStd = 4*1.0 ! repetition pattern rules apply again here. 4 dimensions => 4-element vector of values.
104
105 ! DRAM specifications.
106
107 proposalAdaptationBurnin = 1.
108 proposalAdaptationCount = 10000000
109 proposalAdaptationCountGreedy = 0
110 proposalAdaptationPeriod = 35
111 proposalDelayedRejectionCount = 5
112 proposalDelayedRejectionScale = 4*1., 2. ! The first four elements are 1, followed by 2.
113
114/

Example Unix compile command via Intel ifort compiler ⛓
1#!/usr/bin/env sh
2rm main.exe
3ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example Windows Batch compile command via Intel ifort compiler ⛓
1del main.exe
2set PATH=..\..\..\lib;%PATH%
3ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
4main.exe

Example Unix / MinGW compile command via GNU gfortran compiler ⛓
1#!/usr/bin/env sh
2rm main.exe
3gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Postprocessing of the example output ⛓
1#!/usr/bin/env python
2
3import matplotlib.pyplot as plt
4#import paramonte as pm
5import pandas as pd
6import numpy as np
7import glob
8import sys
9import os
10funcname = os.path.basename(os.getcwd())
11
12linewidth = 2
13fontsize = 17
14files = glob.glob("./out/*_chain.txt")
15
16for file in files:
17
18 df = pd.read_csv(file, delimiter = ",")
19 sindex = list(df.columns).index("sampleLogFunc") + 1 # start column index.
20
21 # traceplot
22
23 #print(df.values)
24 fig = plt.figure(figsize = (8, 6))
25 ax = plt.subplot(1,1,1)
26 ax.plot ( range(len(df.values[:, 0]))
27 , df.values[:, sindex:]
28 , zorder = 1000
29 )
30 plt.minorticks_on()
31 ax.set_xscale("log")
32 ax.set_xlabel("MCMC Count", fontsize = 17)
33 ax.set_ylabel("MCMC State", fontsize = 17)
34 ax.tick_params(axis = "x", which = "minor")
35 ax.tick_params(axis = "y", which = "minor")
36 plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
37 #ax.legend(df.columns.values, fontsize = fontsize)
38 plt.tight_layout()
39 plt.savefig(funcname + ".traceplot.png")
40
41 # scatterplot
42
43 if len(df.values[1, sindex:]) > 1:
44 #print(df.values)
45 fig = plt.figure(figsize = (8, 6))
46 ax = plt.subplot(1,1,1)
47 ax.scatter ( df.values[:, sindex]
48 , df.values[:, sindex + 1]
49 , zorder = 1000
50 , alpha = 0.5
51 , s = 1
52 )
53 plt.minorticks_on()
54 #ax.set_xscale("log")
55 ax.set_xlabel(df.columns.values[sindex], fontsize = 17)
56 ax.set_ylabel(df.columns.values[sindex + 1], fontsize = 17)
57 ax.tick_params(axis = "x", which = "minor")
58 ax.tick_params(axis = "y", which = "minor")
59 plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
60 #ax.legend(df.columns.values, fontsize = fontsize)
61 plt.tight_layout()
62 plt.savefig(funcname + ".scatterplot.png")
63
64 # adaptation measure.
65
66 #print(df.values)
67 if "proposalAdaptation" in df.columns.values:
68 if any(df["proposalAdaptation"].values != 0):
69 fig = plt.figure(figsize = (8, 6))
70 ax = plt.subplot(1,1,1)
71 ax.scatter ( range(len(df.values[:, 0]))
72 , df["proposalAdaptation"].values
73 , zorder = 1000
74 , alpha = 0.5
75 , c = "red"
76 , s = 1
77 )
78 plt.minorticks_on()
79 #ax.set_xscale("log")
80 ax.set_yscale("log")
81 ax.set_xlabel("MCMC Count", fontsize = 17)
82 ax.set_ylabel("proposalAdaptation", fontsize = 17)
83 ax.tick_params(axis = "x", which = "minor")
84 ax.tick_params(axis = "y", which = "minor")
85 plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
86 #ax.legend(df.columns.values, fontsize = fontsize)
87 plt.tight_layout()
88 plt.savefig(funcname + ".proposalAdaptation.png")
89
90 #sim = pm.ParaDRAM()
91 #sample = sim.readSample("./out/", renabled = True)
92 #sample[0].plot.grid()

Visualization of the example output ⛓


Example usage ⛓

1module logfunc
2 use pm_kind, only: IK, RKG => RKD ! all processor kinds are supported.
3 integer(IK), parameter :: NDIM = 2
4contains
5 recursive function getLogFunc(state) result(logFunc)
6 ! Return the negative natural logarithm of Himmelblau function evaluated at the input vector state.
7 real(RKG), intent(in), contiguous :: state(:)
8 real(RKG) :: logFunc
9 if (size(state) /= 2) error stop "The input state vector size must be 2."
10 logFunc = -log((state(1)**2 + state(2) - 11)**2 + (state(1) + state(2)**2 - 7)**2 + 0.1_RKG)
11 end function
12end module logfunc
13
14program example
15 use logfunc, only: NDIM, getLogFunc
16 use pm_sampling, only: getErrSampling, paradram_type, err_type
17 type(paradram_type) :: sampler
18 type(err_type) :: err
19 ! \bug
20 ! gfortran bug requires setting sampler components explicitly.
21 sampler = paradram_type()
22 sampler%inputFile = "input.nml"
23 sampler%outputFileName = "./out/himmelblau"
24 err = getErrSampling(sampler, getLogFunc, NDIM)
25 if (err%occurred) error stop "sampler failed: "//err%msg
26end program example
pm_kind::RKD
integer, parameter RKD
The double precision real kind in Fortran mode. On most platforms, this is an 64-bit real kind.
Definition: pm_kind.F90:568

The contents of example specifications namelist input file input.nml ⛓
1! DESCRIPTION:
2!
3! The external input file for sampling the 2-dimensional Himmelblau function as implemented in the accompanying source files.
4! This file is common between all supported programming language environments.
5!
6! NOTE:
7!
8! All simulation specifications (including this whole file) are optional and can be nearly safely commented out.
9! However, if domain boundaries are finite, set them explicitly.
10!
11! USAGE:
12!
13! -- Comments must begin with an exclamation mark `!`.
14! -- Comments can appear anywhere on an empty line or, after a variable assignment
15! (but not in the middle of a variable assignment whether in single or multiple lines).
16! -- All variable assignments are optional and can be commented out. In such cases, appropriate default values will be assigned.
17! -- Use ParaDRAM namelist (group) name to group a set of ParaDRAM simulation specification variables.
18! -- The order of the input variables in the namelist groups is irrelevant and unimportant.
19! -- Variables can be defined multiple times, but only the last definition will be considered as input.
20! -- All variable names are case insensitive. However, for clarity, this software follows the camelCase code-writing practice.
21! -- String values must be enclosed with either single or double quotation marks.
22! -- 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.
23! -- All vectors and arrays in the input file begin with index 1. This is following the convention of
24! the majority of science-oriented programming languages: Fortran, Julia, Mathematica, MATLAB, and R.
25!
26! For comprehensive guidelines on the input file organization and rules, visit:
27!
28! https://www.cdslab.org/paramonte/generic/latest/usage/sampling/paradram/input/
29!
30! To see detailed descriptions of each of variables, visit:
31!
32! https://www.cdslab.org/paramonte/generic/latest/usage/sampling/paradram/specifications/
33!
34&paradram
35
36 ! Base specifications.
37
38 description = "
39This\n
40 is a\n
41 multi-line\n
42 description.\\n" ! strings must be enclosed with "" or '' and can be continued on multiple lines.
43 ! No comments within strings are allowed.
44 domain = "cube"
45 domainAxisName = "variable1"
46 "variable2" ! values can appear in multiple lines.
47 domainBallAvg = 0 0 ! values can be separated with blanks or commas.
48 domainBallCor = 1 0
49 0 1
50 domainBallCov = 1, 0
51 0, 1
52 domainBallStd = 1, 1
53 domainCubeLimitLower = 2*-1.e30 ! repetition pattern rules apply here. 2 dimensions => 2-element vector of values.
54 domainCubeLimitUpper(1) = +1.e30 ! Elements of vectors can be set individually.
55 domainCubeLimitUpper(2) = +1.e30 ! Elements of vectors can be set individually.
56 domainErrCount = 100
57 domainErrCountMax = 1000
58 inputFileHasPriority = FALSE ! This is relevant only to simulations within the Fortran programming language.
59 outputChainFileFormat = "compact"
60 !outputColumnWidth = 25 ! This is an example of a variable that is commented out.
61 ! Therefore, its value will not be read by the sampler routine.
62 ! To pass it to the routine, simply remove the `!` mark at the beginning of the line.
63 outputFileName = "./out/himmelblau" ! A forward-slash character at the end of the string value would indicate the specified path
64 ! is to be interpreted as the name of the folder to contain the simulation output files.
65 ! The base name for the simulation output files will be generated from the current date and time.
66 ! Otherwise, the specified base name at the end of the string will be used in naming the simulation output files.
67 outputPrecision = 17
68 outputReportPeriod = 1000
69 outputRestartFileFormat = "ascii"
70 outputSampleSize = -1
71 outputSeparator = ","
72 outputSplashMode = "normal" ! or quiet or silent.
73 outputStatus = "extend" ! or retry.
74 parallelism = "multi chain" ! "singleChain" would also work. Similarly, "multichain", "multi chain", or "multiChain".
75 parallelismMpiFinalizeEnabled = false ! TRUE, true, .true., .t., and t would be also all valid logical values representing truth.
76 parallelismNumThread = 3 ! number of threads to use in shared-memory parallelism.
77 !randomSeed = 2136275,
78 !targetAcceptanceRate = 0.23e0
79
80 ! MCMC specifications.
81
82 outputChainSize = 10000
83 outputSampleRefinementCount = 10
84 outputSampleRefinementMethod = "BatchMeans"
85 proposal = "normal" ! or "uniform" as you wish.
86 proposalCor(:, 1) = 1 0 ! first matrix column.
87 proposalCor(:, 2) = 0 1 ! second matrix column.
88 proposalCov = 1 0 0 1 ! or specify all matrix elements in one statement.
89 proposalScale = "2*0.5*Gelman" ! The asterisk here means multiplication since it is enclosed within quotation marks.
90 proposalStart = 2*1.e0 ! four values of 1.e0 are specified here by the repetition pattern symbol *
91 proposalStartDomainCubeLimitLower = 2*-100.e0 ! repetition pattern rules apply again here. 2 dimensions => 2-element vector of values.
92 proposalStartDomainCubeLimitUpper = 2*+100.e0 ! repetition pattern rules apply again here. 2 dimensions => 2-element vector of values.
93 proposalStartRandomized = false
94 proposalStd = 2*1.0 ! repetition pattern rules apply again here. 2 dimensions => 2-element vector of values.
95
96 ! DRAM specifications.
97
98 proposalAdaptationBurnin = 1.
99 proposalAdaptationCount = 10000000
100 proposalAdaptationCountGreedy = 0
101 proposalAdaptationPeriod = 35
102 proposalDelayedRejectionCount = 5
103 proposalDelayedRejectionScale = 4*1., 2. ! The first four elements are 1, followed by 2.
104
105/

Example Unix compile command via Intel ifort compiler ⛓
1#!/usr/bin/env sh
2rm main.exe
3ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example Windows Batch compile command via Intel ifort compiler ⛓
1del main.exe
2set PATH=..\..\..\lib;%PATH%
3ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
4main.exe

Example Unix / MinGW compile command via GNU gfortran compiler ⛓
1#!/usr/bin/env sh
2rm main.exe
3gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Postprocessing of the example output ⛓
1#!/usr/bin/env python
2
3import matplotlib.pyplot as plt
4#import paramonte as pm
5import pandas as pd
6import numpy as np
7import glob
8import sys
9import os
10funcname = os.path.basename(os.getcwd())
11
12linewidth = 2
13fontsize = 17
14files = glob.glob("./out/*_chain.txt")
15
16for file in files:
17
18 df = pd.read_csv(file, delimiter = ",")
19 sindex = list(df.columns).index("sampleLogFunc") + 1 # start column index.
20
21 # traceplot
22
23 #print(df.values)
24 fig = plt.figure(figsize = (8, 6))
25 ax = plt.subplot(1,1,1)
26 ax.plot ( range(len(df.values[:, 0]))
27 , df.values[:, sindex:]
28 , zorder = 1000
29 )
30 plt.minorticks_on()
31 ax.set_xscale("log")
32 ax.set_xlabel("MCMC Count", fontsize = 17)
33 ax.set_ylabel("MCMC State", fontsize = 17)
34 ax.tick_params(axis = "x", which = "minor")
35 ax.tick_params(axis = "y", which = "minor")
36 plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
37 #ax.legend(df.columns.values, fontsize = fontsize)
38 plt.tight_layout()
39 plt.savefig(funcname + ".traceplot.png")
40
41 # scatterplot
42
43 if len(df.values[1, sindex:]) > 1:
44 #print(df.values)
45 fig = plt.figure(figsize = (8, 6))
46 ax = plt.subplot(1,1,1)
47 ax.scatter ( df.values[:, sindex]
48 , df.values[:, sindex + 1]
49 , zorder = 1000
50 , alpha = 0.5
51 , s = 1
52 )
53 plt.minorticks_on()
54 #ax.set_xscale("log")
55 ax.set_xlabel(df.columns.values[sindex], fontsize = 17)
56 ax.set_ylabel(df.columns.values[sindex + 1], fontsize = 17)
57 ax.tick_params(axis = "x", which = "minor")
58 ax.tick_params(axis = "y", which = "minor")
59 plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
60 #ax.legend(df.columns.values, fontsize = fontsize)
61 plt.tight_layout()
62 plt.savefig(funcname + ".scatterplot.png")
63
64 # adaptation measure.
65
66 #print(df.values)
67 if "proposalAdaptation" in df.columns.values:
68 if any(df["proposalAdaptation"].values != 0):
69 fig = plt.figure(figsize = (8, 6))
70 ax = plt.subplot(1,1,1)
71 ax.scatter ( range(len(df.values[:, 0]))
72 , df["proposalAdaptation"].values
73 , zorder = 1000
74 , alpha = 0.5
75 , c = "red"
76 , s = 1
77 )
78 plt.minorticks_on()
79 #ax.set_xscale("log")
80 ax.set_yscale("log")
81 ax.set_xlabel("MCMC Count", fontsize = 17)
82 ax.set_ylabel("proposalAdaptation", fontsize = 17)
83 ax.tick_params(axis = "x", which = "minor")
84 ax.tick_params(axis = "y", which = "minor")
85 plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
86 #ax.legend(df.columns.values, fontsize = fontsize)
87 plt.tight_layout()
88 plt.savefig(funcname + ".proposalAdaptation.png")
89
90 #sim = pm.ParaDRAM()
91 #sample = sim.readSample("./out/", renabled = True)
92 #sample[0].plot.grid()

Visualization of the example output ⛓


Example usage ⛓

1module auxil_mod
2
3 use pm_kind, only: RKG => RKD ! All real kinds are supported.
4
5 implicit none
6
7 real(RKG), parameter :: EPS = sqrt(epsilon(0._RKG))
8 real(RKG), parameter :: LARGE = sqrt(huge(0._RKG))
9
10end module auxil_mod
11
12!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14module data_mod
15
16 use auxil_mod, only: RKG
17 use pm_kind, only: IK
18 implicit none
19
20 private
21 public :: obs_type, posuniv_type, stat_type, data_type
22
23 type, abstract :: obs_type
24 end type
25
26 type, extends(obs_type) :: univ_type
27 real(RKG) :: val
28 end type
29
30 type, extends(univ_type) :: posuniv_type
31 real(RKG) :: log
32 end type
33
34 interface posuniv_type
35 module procedure :: posuniv_typer
36 end interface
37
38 type stat_type
39 real(RKG) :: lim(2)
40 end type
41
42 type :: data_type
43 integer(IK) :: nobs
44 type(stat_type) :: stat
45 class(obs_type), allocatable :: obs(:)
46 end type
47
48 interface data_type
49 module procedure :: data_typer
50 end interface
51
52contains
53
54 function posuniv_typer(val) result(self)
55 real(RKG), intent(in), contiguous :: val(:)
56 type(posuniv_type), allocatable :: self(:)
57 integer(IK) :: iobs
58 allocate(self(size(val)))
59 do concurrent(iobs = 1 : size(val))
60 self(iobs)%log = log(val(iobs))
61 self(iobs)%val = val(iobs)
62 end do
63 end function
64
65 function data_typer(obs, stat) result(self)
66 class(obs_type), intent(in), contiguous :: obs(:)
67 type(stat_type), intent(in) :: stat
68 type(data_type) :: self
69 self%stat = stat
70 self%nobs = size(obs, 1, IK)
71 allocate(self%obs, source = obs)
72 end function
73
74end module data_mod
75
76!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
77
78module domain_mod
79
80 !use pm_container, only: css_type
81 use auxil_mod, only: RKG
82 use pm_kind, only: SK
83 implicit none
84
85 type :: domain_type
86 real(RKG), allocatable :: lower(:)
87 real(RKG), allocatable :: upper(:)
88 !type(css_type), allocatable :: names(:)
89 character(63, SK), allocatable :: names(:)
90 end type
91
92end module domain_mod
93
94!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
95
96module model_mod
97
98 use pm_kind, only: SK, IK
99 use data_mod, only: obs_type
100 use domain_mod, only: domain_type
101 use auxil_mod, only: RKG
102 implicit none
103
104 type, abstract :: model_type
105 type(domain_type) :: domain
106 !real(RKG), allocatable :: param(:)
107 integer(IK) :: npar
108 contains
109 procedure(getParam_proc), deferred :: getParam
110 procedure(setParam_proc), deferred :: setParam
111 procedure(getLogPDF_proc), deferred :: getLogPDF
112 end type
113
114 type :: con_type
115 class(model_type), allocatable :: model
116 end type
117
118 interface con_type
119 module procedure :: con_typer
120 end interface
121
122 abstract interface
123 function getParam_proc(self) result(param)
124 import :: RKG, model_type
125 class(model_type), intent(in) :: self
126 real(RKG) :: param(self%npar)
127 end function
128 subroutine setParam_proc(self, param)
129 import :: RKG, model_type
130 class(model_type), intent(inout) :: self
131 real(RKG), intent(in), contiguous :: param(:)
132 end subroutine
133 impure elemental function getLogPDF_proc(self, obs) result(logPDF)
134 import :: RKG, model_type, obs_type
135 class(model_type), intent(in) :: self
136 class(obs_type), intent(in) :: obs
137 real(RKG) :: logPDF
138 end function
139 end interface
140
141contains
142
143 pure elemental function con_typer(model) result(self)
144 class(model_type), intent(in) :: model
145 type(con_type) :: self
146 self%model = model
147 end function
148
149end module model_mod
150
151!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
152
153module lognorm_mod
154
155 use model_mod, only: domain_type, model_type, IK
156 use auxil_mod, only: RKG, LARGE
157 implicit none
158
159 type, extends(model_type) :: lognorm_type
160 real(RKG), private :: avg, invstd, loginvstd
161 contains
162 procedure :: getParam, setParam, getLogPDF
163 end type
164
165 interface lognorm_type
166 module procedure :: lognorm_typer
167 end interface
168
169contains
170
171 function lognorm_typer(param, domain) result(self)
172 real(RKG), intent(in), contiguous, optional :: param(:)
173 type(domain_type), intent(in), optional :: domain
174 type(lognorm_type) :: self
175 self%npar = 2
176 if (present(domain)) then
177 self%domain = domain
178 else
179 self%domain = domain_type([-LARGE, -LARGE], [+LARGE, +LARGE], [character(63) :: "lognorm_avg", "lognorm_logstd"])
180 end if
181 if (present(param)) call self%setParam(param)
182 end function
183
184 subroutine setParam(self, param)
185 use pm_kind, only: SK, IK
186 class(lognorm_type), intent(inout) :: self
187 real(RKG), intent(in), contiguous :: param(:)
188 if (size(param) /= self%npar) error stop "`lognorm_type%setParam()` takes a vector of two parameters representing `avg` and `logstd`."
189 self%invstd = exp(-param(2))
190 self%loginvstd = -param(2)
191 self%avg = param(1)
192 end subroutine
193
194 function getParam(self) result(param)
195 class(lognorm_type), intent(in) :: self
196 real(RKG) :: param(self%npar)
197 param = [self%avg, -self%loginvstd]
198 end function
199
200 impure elemental function getLogPDF(self, obs) result(logPDF)
201 use pm_distNorm, only: setNormLogPDF
202 use data_mod, only: obs_type, posuniv_type
203 class(lognorm_type), intent(in) :: self
204 class(obs_type), intent(in):: obs
205 real(RKG) :: logPDF
206 select type (obs)
207 type is (posuniv_type)
208 call setNormLogPDF(logPDF, obs%log, mu = self%avg, invSigma = self%invstd, logInvSigma = self%loginvstd)
209 class default
210 error stop "Unrecognized data."
211 end select
212 end function
213
214end module lognorm_mod
215
216!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
217
218module flatPoweto_mod
219
220 use auxil_mod, only: RKG, LARGE
221 use model_mod, only: domain_type, model_type
222 implicit none
223
224 type, extends(model_type) :: flatPoweto_type
225 real(RKG), private :: logbreak, alpha, break, alphap1, logPDFNF, limx(2), loglimx(2)
226 contains
227 procedure, pass :: getParam, setParam, getLogPDF
228 end type
229
230 interface flatPoweto_type
231 module procedure :: flatPoweto_typer
232 end interface
233
234contains
235
236 function flatPoweto_typer(limx, param, domain) result(self)
237 real(RKG), intent(in), contiguous, optional :: param(:)
238 type(domain_type), intent(in), optional :: domain
239 real(RKG), intent(in), contiguous :: limx(:)
240 character(:), allocatable :: msg
241 type(flatPoweto_type) :: self
242 self%npar = 2
243 self%limx = limx
244 self%loglimx = log(limx)
245 if (present(domain)) then
246 self%domain = domain
247 else
248 self%domain = domain_type([-LARGE, -LARGE], [+LARGE, +LARGE], [character(63) :: "flatPoweto_logbreak", "flatPoweto_alpha"])
249 end if
250 if (present(param)) call self%setParam(param)
251 end function
252
253 function getParam(self) result(param)
254 class(flatPoweto_type), intent(in) :: self
255 real(RKG) :: param(self%npar)
256 param = [self%logbreak, self%alpha]
257 end function
258
259 subroutine setParam(self, param)
260 use pm_kind, only: SK, IK
261 use pm_val2str, only: getStr
262 use pm_mathMinMax, only: getMinMax
263 use pm_mathLogAddExp, only: getLogAddExp
264 use pm_mathLogSubExp, only: getLogSubExp
265 real(RKG), intent(in), contiguous :: param(:)
266 class(flatPoweto_type), intent(inout) :: self
267 real(RKG) :: logModelInt1, logModelInt2, small, large
268 if (size(param) /= self%npar) error stop "`lognorm_type%setParam()` takes a vector of two parameters representing `avg` and `logstd`."
269 self%alpha = param(2)
270 self%logbreak = param(1)
271 self%alphap1 = self%alpha + 1
272 self%break = exp(self%logbreak)
273 logModelInt1 = log(self%break - self%limx(1))
274 if (self%alphap1 < 0._RKG) then
275 logModelInt2 = -log(-self%alphap1)
276 large = self%alphap1 * self%logbreak
277 small = self%alphap1 * self%loglimx(2)
278 elseif (0._RKG < self%alphap1) then
279 logModelInt2 = -log(self%alphap1)
280 small = self%alphap1 * self%logbreak
281 large = self%alphap1 * self%loglimx(2)
282 else
283 logModelInt2 = 0._RKG
284 small = self%logbreak
285 large = self%loglimx(2)
286 end if
287 logModelInt2 = logModelInt2 - self%alpha * self%logbreak + getLogSubExp(smaller = small, larger = large)
288 self%logPDFNF = -getLogAddExp(getMinMax(logModelInt1, logModelInt2))
289 end subroutine
290
291 pure elemental function getLogPDF(self, obs) result(logPDF)
292 use data_mod, only: obs_type, posuniv_type
293 class(flatPoweto_type), intent(in) :: self
294 class(obs_type), intent(in) :: obs
295 real(RKG) :: logPDF
296 select type (obs)
297 type is (posuniv_type)
298 if (obs%log < self%logbreak) then
299 logPDF = obs%log
300 else
301 logPDF = self%alphap1 * obs%log - self%alpha * self%logbreak
302 end if
303 logPDF = logPDF + self%logPDFNF
304 class default
305 error stop "Unrecognized data."
306 end select
307 end function
308
309end module flatPoweto_mod
310
311!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
312
313module flatPowetoTapered_mod
314
315 use auxil_mod, only: RKG, LARGE
316 use model_mod, only: domain_type, model_type
317 implicit none
318
319 type, extends(model_type) :: flatPowetoTapered_type
320 real(RKG), private :: logbreak, alpha, beta, break, alphap1, logPDFNF, limx(2), loglimx(2)
321 contains
322 procedure, private :: getLogUDF
323 procedure :: getParam, setParam, getLogPDF
324 end type
325
326 interface flatPowetoTapered_type
327 module procedure :: flatPowetoTapered_typer
328 end interface
329
330contains
331
332 function flatPowetoTapered_typer(limx, param, domain) result(self)
333 real(RKG), intent(in), contiguous, optional :: param(:)
334 type(domain_type), intent(in), optional :: domain
335 real(RKG), intent(in), contiguous :: limx(:)
336 type(flatPowetoTapered_type) :: self
337 self%npar = 3
338 self%limx = limx
339 self%loglimx = log(limx)
340 if (present(domain)) then
341 self%domain = domain
342 else
343 self%domain = domain_type([-LARGE, -LARGE, -LARGE], [+LARGE, +LARGE, +LARGE], [character(63) :: "flatPowetoTapered_logbreak", "flatPowetoTapered_alpha", "flatPowetoTapered_beta"])
344 end if
345 if (present(param)) call self%setParam(param)
346 end function
347
348 function getParam(self) result(param)
349 class(flatPowetoTapered_type), intent(in) :: self
350 real(RKG) :: param(self%npar)
351 param = [self%logbreak, self%alpha, self%beta]
352 end function
353
354 subroutine setParam(self, param)
355 use pm_val2str, only: getStr
356 use pm_kind, only: SK, IK, LK
357 use pm_mathMinMax, only: getMinMax
358 use pm_mathLogAddExp, only: getLogAddExp
359 use pm_quadPack, only: isFailedQuad, getQuadErr, weps
360 class(flatPowetoTapered_type), intent(inout) :: self
361 real(RKG), intent(in), contiguous :: param(:)
362 real(RKG) :: logModelInt1, logModelInt2, logNF
363 character(255, SK) :: msg
364 logical(LK) :: failed
365 if (size(param) /= self%npar) error stop "`flatPowetoFlatPowetoTapered_type` takes a vector of three parameters"&
366 //"representing `logbreak`, `alpha`, `beta`. size(param), npar = "//getStr([size(param, 1, IK), self%npar])
367 self%logbreak = param(1)
368 self%alpha = param(2)
369 self%beta = param(3)
370 self%alphap1 = self%alpha + 1
371 self%break = exp(self%logbreak)
372 logModelInt1 = log(self%break - self%limx(1))
373 logNF = max(self%logbreak, self%alphap1 * self%loglimx(2) - self%alpha * self%logbreak - self%beta * (self%limx(2) - self%break))
374 failed = isFailedQuad(getDensity, lb = self%logbreak, ub = self%loglimx(2), help = weps, integral = logModelInt2, msg = msg)
375 if (failed) then
376 !error stop "@getLogModelInt(): "//trim(msg)
377 self%logPDFNF = sqrt(huge(self%logPDFNF))
378 else
379 logModelInt2 = log(logModelInt2) + logNF
380 self%logPDFNF = -getLogAddExp(getMinMax(logModelInt1, logModelInt2))
381 end if
382 contains
383 function getDensity(logx) result(density)
384 real(RKG), intent(in) :: logx
385 real(RKG) :: density
386 density = exp(self%getLogUDF(logx, exp(logx)) - logNF)
387 end function
388 end subroutine
389
390 pure elemental function getLogPDF(self, obs) result(logPDF)
391 use data_mod, only: obs_type, posuniv_type
392 class(flatPowetoTapered_type), intent(in) :: self
393 class(obs_type), intent(in):: obs
394 real(RKG) :: logPDF
395 select type (obs)
396 type is (posuniv_type)
397 if (obs%log < self%logbreak) then
398 logPDF = obs%log
399 else
400 logPDF = self%getLogUDF(obs%log, obs%val)
401 end if
402 logPDF = logPDF + self%logPDFNF
403 class default
404 error stop "Unrecognized data."
405 end select
406 end function
407
408 pure elemental function getLogUDF(self, logx, x) result(logUDF)
409 class(flatPowetoTapered_type), intent(in) :: self
410 real(RKG), intent(in) :: logx, x
411 real(RKG) :: logUDF
412 logUDF = self%alphap1 * logx - self%alpha * self%logbreak - self%beta * (x - self%break)
413 end function
414
415end module flatPowetoTapered_mod
416
417!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
418
419module mixture_mod
420
421 use pm_val2str, only: getStr
422 use pm_arrayResize, only: setResized
423 use pm_mathLogAddExp, only: getLogAddExp
424 use pm_mathLogSumExp, only: getLogSumExp
425 use model_mod, only: domain_type, model_type, con_type
426 use auxil_mod, only: RKG, EPS, LARGE
427 use data_mod, only: obs_type
428 use pm_kind, only: SK, IK
429 implicit none
430
431 type, extends(model_type) :: mixture_type
432 type(con_type), allocatable :: comp(:)
433 real(RKG), allocatable :: logfrac(:)
434 integer(IK), private :: ncomp
435 contains
436 procedure :: getParam, setParam, getLogPDF
437 end type
438
439 interface mixture_type
440 module procedure :: mixture_typer
441 end interface
442
443contains
444
445 function mixture_typer(comp, param) result(self)
446 type(con_type), intent(in), contiguous :: comp(:)
447 real(RKG), intent(in), contiguous, optional :: param(:)
448 character(:, SK), allocatable :: prefix
449 integer(IK) :: icomp, lpar, upar, ipar
450 type(mixture_type) :: self
451 self%ncomp = size(comp, 1, IK)
452 allocate(self%comp, source = comp)
453 call setResized(self%logfrac, self%ncomp)
454 self%npar = sum([(self%comp(icomp)%model%npar, icomp = 1, self%ncomp)]) + self%ncomp - 1
455 call setResized(self%domain%lower, self%npar)
456 call setResized(self%domain%upper, self%npar)
457 call setResized(self%domain%names, self%npar)
458 lpar = 0
459 do icomp = 1, self%ncomp - 1
460 upar = lpar + self%comp(icomp)%model%npar
461 prefix = SK_"comp"//getStr(icomp)//SK_"_"
462 self%domain%lower(lpar + 1 : upar + 1) = [self%comp(icomp)%model%domain%lower, -LARGE]
463 self%domain%upper(lpar + 1 : upar + 1) = [self%comp(icomp)%model%domain%upper, +LARGE]
464 self%domain%names(lpar + 1 : upar + 1) = [character(63, SK) :: prefix//self%comp(icomp)%model%domain%names, prefix//SK_"fisherz(frac)"]
465 lpar = upar + 1
466 end do
467 prefix = SK_"comp"//getStr(icomp)//SK_"_"
468 self%domain%lower(lpar + 1 :) = self%comp(icomp)%model%domain%lower
469 self%domain%upper(lpar + 1 :) = self%comp(icomp)%model%domain%upper
470 self%domain%names(lpar + 1 :) = prefix//self%comp(icomp)%model%domain%names
471 if (present(param)) call self%setParam(param)
472 end function
473
474 subroutine setParam(self, param)
475 use pm_mathFisher, only: getFisherInv
476 class(mixture_type), intent(inout) :: self
477 real(RKG), intent(in), contiguous :: param(:)
478 integer(IK) :: icomp, lpar, upar
479 real(RKG) :: sumfrac, mixfrac
480 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])
481 lpar = 0
482 sumfrac = 0
483 do icomp = 1, self%ncomp - 1
484 upar = lpar + self%comp(icomp)%model%npar
485 call self%comp(icomp)%model%setParam(param(lpar + 1 : upar))
486 mixfrac = getFisherInv(param(upar + 1), 0._RKG, 1._RKG)
487 self%logfrac(icomp) = log(mixfrac)
488 sumfrac = sumfrac + mixfrac
489 lpar = upar + 1
490 end do
491 self%logfrac(icomp) = log(1 - sumfrac)
492 call self%comp(icomp)%model%setParam(param(lpar + 1 :))
493 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)
494 end subroutine
495
496 function getParam(self) result(param)
497 use pm_mathFisher, only: getFisher
498 class(mixture_type), intent(in) :: self
499 integer(IK) :: icomp, lpar, upar
500 real(RKG) :: param(self%npar)
501 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()]
502 end function
503
504 impure elemental function getLogPDF(self, obs) result(logPDF)
505 class(mixture_type), intent(in) :: self
506 real(RKG) :: logPDFS(self%ncomp), maxLogPDF
507 class(obs_type), intent(in) :: obs
508 real(RKG) :: logPDF
509 integer(IK) :: icomp
510 maxLogPDF = -huge(0._RKG)
511 do icomp = 1, self%ncomp
512 logPDFS(icomp) = self%logfrac(icomp) + self%comp(icomp)%model%getLogPDF(obs)
513 maxLogPDF = max(maxLogPDF, logPDFS(icomp))
514 end do
515 logPDF = getLogSumExp(logPDFS, maxLogPDF)
516 end function
517
518end module mixture_mod
519
520!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
521
522module fit_mod
523
524 use auxil_mod, only: RKG
525 use pm_sampling, only: sampler_type
526 use pm_sampling, only: paradram_type
527 use model_mod, only: model_type
528 use data_mod, only: data_type
529 implicit none
530
531 private
532 public :: fit_type, paradram_type
533
534 type :: best_type
535 real(RKG) :: loglike
536 real(RKG), allocatable :: param(:) ! best fit model parameters.
537 end type
538
539 type :: fit_type
540 real(RKG) :: bic
541 type(best_type) :: best
542 type(data_type) :: data
543 class(model_type), allocatable :: model
544 class(sampler_type), allocatable :: sampler
545 contains
546 procedure, pass :: run
547 !procedure, pass :: setHist
548 end type
549
550 interface fit_type
551 module procedure :: fit_typer
552 end interface
553
554contains
555
556 function fit_typer(data, model, sampler) result(self)
557 class(sampler_type), intent(in), optional :: sampler
558 class(model_type), intent(in) :: model
559 type(data_type), intent(in) :: data
560 type(fit_type) :: self
561 if (present(sampler)) then
562 self%sampler = sampler
563 else
564 self%sampler = paradram_type()
565 end if
566 self%model = model
567 self%data = data
568 end function
569
570 subroutine run(self)
571
572 use pm_kind, only: SK, IK
573 use pm_err, only: err_type
574 use pm_sampling, only: getErrSampling, paradram_type
575 class(fit_type), intent(inout) :: self
576 type(err_type) :: err
577
578 if (.not. allocated(self%sampler%outputStatus)) self%sampler%outputStatus = "retry"
579 if (.not. allocated(self%sampler%domainAxisName)) self%sampler%domainAxisName = self%model%domain%names
580 if (.not. allocated(self%sampler%domainCubeLimitLower)) self%sampler%domainCubeLimitLower = self%model%domain%lower
581 if (.not. allocated(self%sampler%domainCubeLimitUpper)) self%sampler%domainCubeLimitUpper = self%model%domain%upper
582
583 select type (sampler => self%sampler)
584 type is (paradram_type)
585 if (.not. allocated(sampler%outputChainSize)) sampler%outputChainSize = 30000
586 if (.not. allocated(sampler%proposalScale)) sampler%proposalScale = "gelman"
587 if (.not. allocated(sampler%proposalStart)) sampler%proposalStart = self%model%getParam()
588 err = getErrSampling(sampler, getLogLike, self%model%npar)
589 end select
590 if (err%occurred) error stop "sampler failed: "//err%msg
591
592 ! Find the mean best fit parameters and write post-prediction data to output file for visualization.
593
594 block
595
596 use pm_io, only: getErrTableRead
597 use pm_sysPath, only: glob, css_type
598 use pm_parallelism, only: getImageID
599
600 integer(IK) :: stat, ibest!, offset = 1
601 type(css_type), allocatable :: path(:)
602 real(RKG), allocatable :: logx(:), logPDF(:)
603 real(RKG), allocatable :: table(:,:), state(:)
604
605 if (getImageID() == 1) write(*, "(A)") "Searching for files: "//self%sampler%outputFileName//SK_"*"
606 path = glob(self%sampler%outputFileName//SK_"*")
607 if (size(path) == 0) error stop "There is no sample file in the output folder."
608
609 stat = getErrTableRead(path(size(path))%val, table, roff = 1_IK)
610 if (stat /= 0) error stop "Failed to read output sample."
611 ibest = maxloc(table(:, 1), dim = 1_IK)
612 self%best%loglike = table(ibest, 1)
613 self%best%param = table(ibest, 2:)
614 self%bic = self%model%npar * log(real(self%data%nobs, RKG)) - 2 * self%best%loglike
615
616 end block
617
618 contains
619
620 function getLogLike(param) result(logLike)
621 real(RKG), intent(in), contiguous :: param(:)
622 real(RKG) :: logLike
623 call self%model%setParam(param)
624 loglike = sum(self%model%getLogPDF(self%data%obs))
625 end function
626
627 end subroutine
628
629end module fit_mod
630
631!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
632
633program example
634
635 use pm_kind, only: SK, IK
636 use fit_mod, only: fit_type
637 use model_mod, only: con_type
638 use auxil_mod, only: RKG, LARGE
639 use lognorm_mod, only: lognorm_type
640 use mixture_mod, only: mixture_type
641 use pm_sampling, only: paradram_type
642 use pm_arraySpace, only: getLogSpace
643 use flatPoweto_mod, only: flatPoweto_type
644 use data_mod, only: posuniv_type, data_type, stat_type
645 use flatPowetoTapered_mod, only: flatPowetoTapered_type
646 use pm_io, only: display_type, css_type, getErrTableWrite
647
648 integer(IK) :: imodel, stat
649 type(paradram_type) :: sampler
650 type(data_type) :: data, pred
651 type(display_type) :: disp
652 type(fit_type) :: fit
653
654 disp = display_type(file = SK_"main.out.F90")
655
656 ! read data.
657
658 block
659 use pm_io, only: getErrTableRead
660 real(RKG), allocatable :: table(:,:)
661 character(255, SK) :: iomsg
662 stat = getErrTableRead("data.csv", table, roff = 1_IK, iomsg = iomsg)
663 if (stat /= 0) error stop "Failed to read data: "//trim(iomsg)
664 associate(val => table(:, 4))
665 data = data_type(obs = posuniv_type(val), stat = stat_type([minval(val), maxval(val)]))
666 end associate
667 end block
668
669 ! perform fitting.
670
671 do imodel = 1, 4
672
673 sampler = paradram_type()
674 sampler%parallelismMpiFinalizeEnabled = .false.
675
676 if (imodel == 1) then
677 sampler%outputFileName = "./mixLogNormLogNorm"
678 sampler%proposalStart = [-.2, 0.4, .4, 3.6, -0.2]
679 sampler%domainCubeLimitLower = [real(RKG) :: -LARGE, -LARGE, -LARGE, log(3.), -LARGE]
680 sampler%domainCubeLimitUpper = [real(RKG) :: log(3.), +LARGE, +LARGE, +LARGE, +LARGE]
681 !sampler%proposalScale = "0.1 * gelman"
682 sampler%proposalCov = reshape ( &
683 [ 1.6E-2, 5.5E-3, 4.1E-3, 3.9E-3,-2.5E-3 &
684 , 5.5E-3, 3.1E-3, 1.7E-3, 1.4E-3,-9.5E-4 &
685 , 4.1E-3, 1.7E-3, 1.9E-3, 1.2E-3,-8.8E-4 &
686 , 3.9E-3, 1.4E-3, 1.2E-3, 2.0E-3,-8.9E-4 &
687 , -2.5E-3,-9.5E-4,-8.8E-4,-8.9E-4, 1.0E-3 &
688 ], shape = [size(sampler%proposalStart), size(sampler%proposalStart)])
689 fit = fit_type(data, mixture_type([con_type(lognorm_type()), con_type(lognorm_type())]), sampler)
690 elseif (imodel == 2) then
691 sampler%proposalScale = "0.1 * gelman"
692 sampler%outputFileName = "./mixFlatPowetoFlatPowetoTapered"
693 sampler%proposalStart = [log(.37), -1.4, .3, log(21.3), -.5, 0.0156]
694 sampler%proposalCov = reshape ( &
695 [ 1.E-2, -3.E-3, 1.E-2, -7.E-3, -1.E-4, -6.E-4 &
696 , -3.E-3, 3.E-3, -4.E-3, 6.E-3, 1.E-4, 1.E-3 &
697 , 1.E-2, -4.E-3, 1.E-1, -7.E-2, -8.E-4, 3.E-4 &
698 , -7.E-3, 6.E-3, -7.E-2, 7.E-2, 8.E-4, 3.E-3 &
699 , -1.E-4, 1.E-4, -8.E-4, 8.E-4, 1.E-5, 6.E-5 &
700 , -6.E-4, 1.E-3, 3.E-4, 3.E-3, 6.E-5, 2.E-3 &
701 ], shape = [size(sampler%proposalStart), size(sampler%proposalStart)])
702 fit = fit_type(data, mixture_type([con_type(flatPoweto_type(data%stat%lim)), con_type(flatPowetoTapered_type(data%stat%lim))]), sampler)
703 elseif (imodel == 3) then
704 sampler%proposalScale = "0.1 * gelman"
705 sampler%outputFileName = "./mixLogNormFlatPowetoTapered"
706 sampler%proposalStart = [-.2, 0.4, .3, log(21.3), -.5, 0.0156]
707 fit = fit_type(data, mixture_type([con_type(lognorm_type()), con_type(flatPowetoTapered_type(data%stat%lim))]), sampler)
708 elseif (imodel == 4) then
709 sampler%proposalScale = "0.1 * gelman"
710 sampler%outputFileName = "./mixFlatPowetoLogNorm"
711 sampler%proposalStart = [log(.37), -1.4, .3, 3.6, -0.2]
712 fit = fit_type(data, mixture_type([con_type(flatPoweto_type(data%stat%lim)), con_type(lognorm_type())]), sampler)
713 else
714 sampler%outputFileName = "./logNorm"
715 fit = fit_type(data, lognorm_type([real(RKG) :: 1, 1]), sampler)
716 end if
717
718 call fit%run()
719
720 call disp%show(css_type(fit%sampler%domainAxisName))
721 call disp%show(fit%best%param)
722 call disp%skip()
723 call disp%show("[fit%best%loglike, fit%bic]")
724 call disp%show( [fit%best%loglike, fit%bic] )
725 call disp%skip()
726
727 ! generate synthetic data for histogram reconstruction.
728
729 associate(val => getLogSpace(logx1 = log(data%stat%lim(1)), logx2 = log(data%stat%lim(2)), count = 1000_IK))
730 pred = data_type(obs = posuniv_type(val), stat = stat_type([minval(val), maxval(val)]))
731 call fit%model%setParam(fit%best%param)
732 stat = getErrTableWrite(fit%sampler%outputFileName//SK_".hist", reshape([log(val), fit%model%getLogPDF(pred%obs)], [size(val), 2]), header = SK_"logx,logPDF")
733 if (stat /= 0) error stop "Failed to write the histogram visualization data."
734 end associate
735
736 end do
737
738end
pm_arrayResize::setResized
Allocate or resize (shrink or expand) an input allocatable scalar string or array of rank 1....
Definition: pm_arrayResize.F90:249
pm_arraySpace::getLogSpace
Generate count evenly-logarithmically-spaced points over the interval [base**logx1,...
Definition: pm_arraySpace.F90:538
pm_distNorm::setNormLogPDF
Generate the natural logarithm of probability density function (PDF) of the univariate Normal distrib...
Definition: pm_distNorm.F90:437
pm_io::getErrTableRead
Generate and return the iostat code resulting from reading the contents of the specified file or unit...
Definition: pm_io.F90:2390
pm_io::getErrTableWrite
Generate and return the iostat code resulting from writing the input table of rank 1 or 2 to the spec...
Definition: pm_io.F90:5940
pm_io::show
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11726
pm_mathFisher::getFisherInv
Generate and return the inverse Fisher transformation of the input Fisher z value.
Definition: pm_mathFisher.F90:335
pm_mathFisher::getFisher
Generate and return the Fisher transformation of the input Fisher z value.
Definition: pm_mathFisher.F90:153
pm_mathLogAddExp::getLogAddExp
Generate and return the logarithm of the sum of the exponential of two input logarithmic values robus...
Definition: pm_mathLogAddExp.F90:166
pm_mathLogSubExp::getLogSubExp
Generate and return the natural logarithm of the subtraction of the exponential of the smaller input ...
Definition: pm_mathLogSubExp.F90:166
pm_mathLogSumExp::getLogSumExp
Generate and return the natural logarithm of the sum of the exponential of the input array robustly (...
Definition: pm_mathLogSumExp.F90:147
pm_mathMinMax::getMinMax
Generate and return an array of size two containing the minimum and maximum of the two input values i...
Definition: pm_mathMinMax.F90:117
pm_quadPack::getQuadErr
Compute the 1D integral of the input scalar (potentially singular) integrand getFunc on a finite or s...
Definition: pm_quadPack.F90:5643
pm_quadPack::isFailedQuad
Compute the 1D integral of the input scalar (potentially singular) integrand getFunc on a finite or s...
Definition: pm_quadPack.F90:4989
pm_val2str::getStr
Generate and return the conversion of the input value to an output Fortran string,...
Definition: pm_val2str.F90:167
pm_arrayResize
This module contains procedures and generic interfaces for resizing allocatable arrays of various typ...
Definition: pm_arrayResize.F90:81
pm_arraySpace
This module contains procedures and generic interfaces for generating arrays with linear or logarithm...
Definition: pm_arraySpace.F90:33
pm_distNorm
This module contains classes and procedures for computing various statistical quantities related to t...
Definition: pm_distNorm.F90:179
pm_mathFisher
This module contains procedures and generic interfaces for evaluating the Fisher transformation and i...
Definition: pm_mathFisher.F90:83
pm_mathLogAddExp
This module contains procedures and generic interfaces for adding two real or complex values without ...
Definition: pm_mathLogAddExp.F90:38
pm_mathLogSubExp
This module contains procedures and generic interfaces for subtracting two real or complex values wit...
Definition: pm_mathLogSubExp.F90:37
pm_mathLogSumExp
This module contains the procedures and interfaces for computing the natural logarithm of the sum of ...
Definition: pm_mathLogSumExp.F90:39
pm_mathMinMax
This module contains procedures and generic interfaces for finding the minimum and maximum of two inp...
Definition: pm_mathMinMax.F90:39
pm_parallelism
This module contains procedures and generic interfaces for facilitating parallel computations or comp...
Definition: pm_parallelism.F90:35
pm_parallelism::getImageID
integer(IK) function getImageID()
Generate and return the ID of the current Coarray image / MPI process / OpenMP thread,...
Definition: pm_parallelism.F90:524
pm_quadPack
This module contains classes and procedures for non-adaptive and adaptive global numerical quadrature...
Definition: pm_quadPack.F90:141
pm_quadPack::weps
type(weps_type), parameter weps
The scalar constant object of type weps_type that indicates the use of Epsilon extrapolation method o...
Definition: pm_quadPack.F90:4551
pm_val2str
This module contains the generic procedures for converting values of different types and kinds to For...
Definition: pm_val2str.F90:58
pm_kind::model_type
This is the abstract derived type for creating objects of class model_type that contain the character...
Definition: pm_kind.F90:1504
pm_sampling::sampler_type
This is a derived type for constructing objects containing the optional simulation properties of the ...
Definition: pm_sampling.F90:80

Example Unix compile command via Intel ifort compiler ⛓
1#!/usr/bin/env sh
2rm main.exe
3ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example Windows Batch compile command via Intel ifort compiler ⛓
1del main.exe
2set PATH=..\..\..\lib;%PATH%
3ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
4main.exe

Example Unix / MinGW compile command via GNU gfortran compiler ⛓
1#!/usr/bin/env sh
2rm main.exe
3gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Postprocessing of the example output ⛓
1#!/usr/bin/env python
2
3import matplotlib.pyplot as plt
4#import paramonte as pm
5import pandas as pd
6import numpy as np
7import glob
8import sys
9import os
10
11plt.rcParams['text.usetex'] = True
12# Uncomment above line if LaTeX is installed on the system.
13# On Ubuntu, this requires: sudo apt install dvipng texlive-latex-extra texlive-fonts-recommended cm-super
14examname = os.path.basename(os.getcwd())
15
16
17
18linewidth = 2
19fontsize = 17
20files = glob.glob("./*_chain.txt")
21
22for file in files:
23
24 basename = file.split("_")[0]
25 df = pd.read_csv(file, delimiter = ",")
26 if "_chain.txt" in file:
27 sindex = 7 # start column index.
28 elif "_sample.txt" in file:
29 sindex = 1 # start column index.
30 else:
31 sys.exit("Unrecognized simulation output file: " + file)
32
33 # traceplot
34
35 #print(df.values)
36 fig = plt.figure(figsize = (8, 6))
37 ax = plt.subplot(1,1,1)
38 ax.plot ( range(len(df.values[:, 0]))
39 , df.values[:, sindex:]
40 , zorder = 1000
41 )
42 plt.minorticks_on()
43 ax.set_xscale("log")
44 ax.set_xlabel("MCMC Count", fontsize = 17)
45 ax.set_ylabel("MCMC State", fontsize = 17)
46 ax.tick_params(axis = "x", which = "major", labelsize = fontsize)
47 ax.tick_params(axis = "y", which = "major", labelsize = fontsize)
48 plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
49 #ax.legend(df.columns.values, fontsize = fontsize)
50 #plt.title(basename)
51 plt.tight_layout()
52 plt.savefig(basename + ".traceplot.png")
53
54 # scatterplot
55
56 if len(df.values[1, sindex:]) == 2:
57 #print(df.values)
58 fig = plt.figure(figsize = (8, 6))
59 ax = plt.subplot(1,1,1)
60 ax.scatter ( df.values[:, sindex]
61 , df.values[:, sindex + 1]
62 , zorder = 1000
63 , s = 1
64 )
65 plt.minorticks_on()
66 ax.set_xscale("log")
67 ax.set_xlabel(df.columns.values[sindex], fontsize = 17)
68 ax.set_ylabel(df.columns.values[sindex + 1], fontsize = 17)
69 ax.tick_params(axis = "x", which = "major", labelsize = fontsize)
70 ax.tick_params(axis = "y", which = "major", labelsize = fontsize)
71 plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
72 #ax.legend(df.columns.values, fontsize = fontsize)
73 #plt.title(basename)
74 plt.tight_layout()
75 plt.savefig(basename + ".scatterplot.png")
76
77 # adaptation measure.
78
79 #print(df.values)
80 if "proposalAdaptation" in df.columns.values:
81 if any(df["proposalAdaptation"].values != 0):
82 fig = plt.figure(figsize = (8, 6))
83 ax = plt.subplot(1,1,1)
84 ax.scatter ( range(len(df.values[:, 0]))
85 , df["proposalAdaptation"].values
86 , zorder = 1000
87 , c = "red"
88 , s = 1
89 )
90 plt.minorticks_on()
91 #ax.set_xscale("log")
92 ax.set_yscale("log")
93 ax.set_xlabel("MCMC Count", fontsize = 17)
94 ax.tick_params(axis = "x", which = "major", labelsize = fontsize)
95 ax.tick_params(axis = "y", which = "major", labelsize = fontsize)
96 ax.set_ylabel("proposalAdaptation", fontsize = 17)
97 plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
98 #ax.legend(df.columns.values, fontsize = fontsize)
99 #plt.title(basename)
100 plt.tight_layout()
101 plt.savefig(basename + ".proposalAdaptation.png")
102
103 #sim = pm.ParaDRAM()
104 #sample = sim.readSample("./out/", renabled = True)
105 #sample[0].plot.grid()
106
107
110
111# from Functions import mergeBins, findPlateauLocation
112
113def mergeBins(counts, bins, threshold = 5):
114
115 if len(bins)%2 == 0:
116 halfway = int(np.ceil(len(bins)/2))
117 else:
118 halfway = int(np.floor(len(bins)/2))
119
120 threshold = threshold
121 leftCounts = []
122 leftBins = [bins[0]]
123 rightCounts = []
124 rightBins = [bins[-1]]
125 leftSumCounts = 0
126 rightSumCounts = 0
127 leftBinstart = 0
128 rightBinstart = 0
129 for i in range(0 , halfway): #goes to the middle of the distribution
130 # print(counts[i],sumCounts)
131
132 leftSumCounts += counts[i]
133
134
135
136 if leftSumCounts >= threshold:
137 # print("added value to list")
138 leftCounts += [leftSumCounts]
139 leftSumCounts = 0
140 leftBins += [bins[i+1]]
141
142 #if counts[-i-1]==3:print(rightSumCounts,-i-1)
143
144
145 rightSumCounts += counts[-i-1]
146 if rightSumCounts >= threshold:
147 rightCounts += [rightSumCounts]
148 rightSumCounts = 0
149 rightBins += [bins[-i-2]]
150
151 if len(bins)%2==0:
152 if len(rightCounts) != 0:
153 rightCounts.pop()
154 rightBins.pop()
155 leftBins.pop()
156 else:
157 rightCounts += [counts[halfway]]
158 #rightBins += [bins[halfway]]
159
160 newCounts = leftCounts + rightCounts[-1::-1]
161 newBins = leftBins + rightBins[-1::-1]
162 #print(len(newCounts),len(newBins))
163 origionaldNdT = np.array(counts) / np.array([bins[i + 1] - bins[i] for i in range(0, len(bins) - 1) ])
164 newdNdT = np.array(newCounts) / np.array([newBins[i + 1] - newBins[i] for i in range(0, len(newBins) - 1) ])
165
166 return np.array(newCounts), np.array(newBins), np.array(origionaldNdT), np.array(newdNdT)
167
168# Read the duration data.
169
170data = pd.read_csv("data.csv", delimiter = ",")
171durs = data["T90"].values
172
173nbins = 50
174binning = 10**np.linspace(np.log10(min(durs)), np.log10(max(durs)), nbins)
175countsAndBins = np.histogram(durs, bins = binning)
176counts = countsAndBins[0]
177bins = countsAndBins[1]
178bins = bins.tolist()
179newCounts, newBins, origionaldNdT, dNdT = mergeBins(counts, binning)
180
181# get the middle of the bin edges for location of uncertainties.
182#midBins = [newBins[i] + (newBins[i + 1] - newBins[i]) / 2 for i in range(0, len(newBins) - 1)]
183# use algorithm to find the longest plateau within one std of right most point
184#plateau, binEdge, plateauErrors = findPlateauLocation(dNdT, dNdTerrors, newBins)
185# find the best value for the plateau
186#bestParams = opt.minimize(chisqfunc, [plateau[-1]], args = (plateau, plateauErrors))
187#x = np.logspace(np.log10(min(binEdge)),np.log10(max(binEdge)),len(midBins))
188#plt.plot(x,(bestParams["x"]).tolist()*len(x),color = 'k',lw = 2)
189#x = np.logspace(np.log10(min(binEdge)), np.log10(max(binEdge)), len(midBins))
190
191histFiles = glob.glob("./*.hist")
192
193for histFile in histFiles:
194
195 fig = plt.figure()
196 plt.stairs ( dNdT
197 , newBins
198 , linewidth = 2
199 , baseline = None
200 , color = "green"
201 , label = 'BATSE'
202 )
203 plt.xscale('log')
204 plt.yscale('log')
205
206 # Read the fitting data to plot.
207
208 model = pd.read_csv(histFile, delimiter = ",")
209 x = np.exp(model["logx"].values)
210 y = len(durs) * np.exp(model["logPDF"].values) / x
211 plt.grid(which = "both")
212 plt.plot(x, y, c = "magenta", label = histFile)
213 plt.xlabel(r"$T_{90} ~ [s]$", fontsize = fontsize)
214 plt.ylabel(r"$dN ~ / ~ dT_{90} ~ [1 / s]$", fontsize = fontsize)
215 #plt.legend(fontsize = 13)
216 plt.xticks(size = fontsize)
217 plt.yticks(size = fontsize)
218 plt.ylim(7*10**-3, 2*10**3)
219 plt.tight_layout()
220 plt.savefig(histFile + ".png")
221 #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: