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:652
pm_sampling
This module contains procedures and generic interfaces for the ParaMonte library sampler routines.
Definition: pm_sampling.F90:30
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:244

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) error stop "@getLogModelInt(): "//trim(msg)
376 logModelInt2 = log(logModelInt2) + logNF
377 self%logPDFNF = -getLogAddExp(getMinMax(logModelInt1, logModelInt2))
378 contains
379 function getDensity(logx) result(density)
380 real(RKG), intent(in) :: logx
381 real(RKG) :: density
382 density = exp(self%getLogUDF(logx, exp(logx)) - logNF)
383 end function
384 end subroutine
385
386 pure elemental function getLogPDF(self, obs) result(logPDF)
387 use data_mod, only: obs_type, posuniv_type
388 class(flatPowetoTapered_type), intent(in) :: self
389 class(obs_type), intent(in):: obs
390 real(RKG) :: logPDF
391 select type (obs)
392 type is (posuniv_type)
393 if (obs%log < self%logbreak) then
394 logPDF = obs%log
395 else
396 logPDF = self%getLogUDF(obs%log, obs%val)
397 end if
398 logPDF = logPDF + self%logPDFNF
399 class default
400 error stop "Unrecognized data."
401 end select
402 end function
403
404 pure elemental function getLogUDF(self, logx, x) result(logUDF)
405 class(flatPowetoTapered_type), intent(in) :: self
406 real(RKG), intent(in) :: logx, x
407 real(RKG) :: logUDF
408 logUDF = self%alphap1 * logx - self%alpha * self%logbreak - self%beta * (x - self%break)
409 end function
410
411end module flatPowetoTapered_mod
412
413!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
414
415module mixture_mod
416
417 use pm_val2str, only: getStr
418 use pm_arrayResize, only: setResized
419 use pm_mathLogAddExp, only: getLogAddExp
420 use pm_mathLogSumExp, only: getLogSumExp
421 use model_mod, only: domain_type, model_type, con_type
422 use auxil_mod, only: RKG, EPS, LARGE
423 use data_mod, only: obs_type
424 use pm_kind, only: SK, IK
425 implicit none
426
427 type, extends(model_type) :: mixture_type
428 type(con_type), allocatable :: comp(:)
429 real(RKG), allocatable :: logfrac(:)
430 integer(IK), private :: ncomp
431 contains
432 procedure :: getParam, setParam, getLogPDF
433 end type
434
435 interface mixture_type
436 module procedure :: mixture_typer
437 end interface
438
439contains
440
441 function mixture_typer(comp, param) result(self)
442 type(con_type), intent(in), contiguous :: comp(:)
443 real(RKG), intent(in), contiguous, optional :: param(:)
444 character(:, SK), allocatable :: prefix
445 integer(IK) :: icomp, lpar, upar, ipar
446 type(mixture_type) :: self
447 self%ncomp = size(comp, 1, IK)
448 allocate(self%comp, source = comp)
449 call setResized(self%logfrac, self%ncomp)
450 self%npar = sum([(self%comp(icomp)%model%npar, icomp = 1, self%ncomp)]) + self%ncomp - 1
451 call setResized(self%domain%lower, self%npar)
452 call setResized(self%domain%upper, self%npar)
453 call setResized(self%domain%names, self%npar)
454 lpar = 0
455 do icomp = 1, self%ncomp - 1
456 upar = lpar + self%comp(icomp)%model%npar
457 prefix = SK_"comp"//getStr(icomp)//SK_"_"
458 self%domain%lower(lpar + 1 : upar + 1) = [self%comp(icomp)%model%domain%lower, -LARGE]
459 self%domain%upper(lpar + 1 : upar + 1) = [self%comp(icomp)%model%domain%upper, +LARGE]
460 self%domain%names(lpar + 1 : upar + 1) = [character(63, SK) :: prefix//self%comp(icomp)%model%domain%names, prefix//SK_"fisherz(frac)"]
461 lpar = upar + 1
462 end do
463 prefix = SK_"comp"//getStr(icomp)//SK_"_"
464 self%domain%lower(lpar + 1 :) = self%comp(icomp)%model%domain%lower
465 self%domain%upper(lpar + 1 :) = self%comp(icomp)%model%domain%upper
466 self%domain%names(lpar + 1 :) = prefix//self%comp(icomp)%model%domain%names
467 if (present(param)) call self%setParam(param)
468 end function
469
470 subroutine setParam(self, param)
471 use pm_mathFisher, only: getFisherInv
472 class(mixture_type), intent(inout) :: self
473 real(RKG), intent(in), contiguous :: param(:)
474 integer(IK) :: icomp, lpar, upar
475 real(RKG) :: sumfrac, mixfrac
476 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])
477 lpar = 0
478 sumfrac = 0
479 do icomp = 1, self%ncomp - 1
480 upar = lpar + self%comp(icomp)%model%npar
481 call self%comp(icomp)%model%setParam(param(lpar + 1 : upar))
482 mixfrac = getFisherInv(param(upar + 1), 0._RKG, 1._RKG)
483 self%logfrac(icomp) = log(mixfrac)
484 sumfrac = sumfrac + mixfrac
485 lpar = upar + 1
486 end do
487 self%logfrac(icomp) = log(1 - sumfrac)
488 call self%comp(icomp)%model%setParam(param(lpar + 1 :))
489 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)
490 end subroutine
491
492 function getParam(self) result(param)
493 use pm_mathFisher, only: getFisher
494 class(mixture_type), intent(in) :: self
495 integer(IK) :: icomp, lpar, upar
496 real(RKG) :: param(self%npar)
497 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()]
498 end function
499
500 impure elemental function getLogPDF(self, obs) result(logPDF)
501 class(mixture_type), intent(in) :: self
502 real(RKG) :: logPDFS(self%ncomp), maxLogPDF
503 class(obs_type), intent(in) :: obs
504 real(RKG) :: logPDF
505 integer(IK) :: icomp
506 maxLogPDF = -huge(0._RKG)
507 do icomp = 1, self%ncomp
508 logPDFS(icomp) = self%logfrac(icomp) + self%comp(icomp)%model%getLogPDF(obs)
509 maxLogPDF = max(maxLogPDF, logPDFS(icomp))
510 end do
511 logPDF = getLogSumExp(logPDFS, maxLogPDF)
512 end function
513
514end module mixture_mod
515
516!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
517
518module fit_mod
519
520 use auxil_mod, only: RKG
521 use pm_sampling, only: sampler_type
522 use pm_sampling, only: paradram_type
523 use model_mod, only: model_type
524 use data_mod, only: data_type
525 implicit none
526
527 private
528 public :: fit_type, paradram_type
529
530 type :: best_type
531 real(RKG) :: loglike
532 real(RKG), allocatable :: param(:) ! best fit model parameters.
533 end type
534
535 type :: fit_type
536 real(RKG) :: bic
537 type(best_type) :: best
538 type(data_type) :: data
539 class(model_type), allocatable :: model
540 class(sampler_type), allocatable :: sampler
541 contains
542 procedure, pass :: run
543 !procedure, pass :: setHist
544 end type
545
546 interface fit_type
547 module procedure :: fit_typer
548 end interface
549
550contains
551
552 function fit_typer(data, model, sampler) result(self)
553 class(sampler_type), intent(in), optional :: sampler
554 class(model_type), intent(in) :: model
555 type(data_type), intent(in) :: data
556 type(fit_type) :: self
557 if (present(sampler)) then
558 self%sampler = sampler
559 else
560 self%sampler = paradram_type()
561 end if
562 self%model = model
563 self%data = data
564 end function
565
566 subroutine run(self)
567
568 use pm_kind, only: SK, IK
569 use pm_err, only: err_type
570 use pm_sampling, only: getErrSampling, paradram_type
571 class(fit_type), intent(inout) :: self
572 type(err_type) :: err
573
574 if (.not. allocated(self%sampler%outputStatus)) self%sampler%outputStatus = "retry"
575 if (.not. allocated(self%sampler%domainAxisName)) self%sampler%domainAxisName = self%model%domain%names
576 if (.not. allocated(self%sampler%domainCubeLimitLower)) self%sampler%domainCubeLimitLower = self%model%domain%lower
577 if (.not. allocated(self%sampler%domainCubeLimitUpper)) self%sampler%domainCubeLimitUpper = self%model%domain%upper
578
579 select type (sampler => self%sampler)
580 type is (paradram_type)
581 if (.not. allocated(sampler%outputChainSize)) sampler%outputChainSize = 30000
582 if (.not. allocated(sampler%proposalScale)) sampler%proposalScale = "gelman"
583 if (.not. allocated(sampler%proposalStart)) sampler%proposalStart = self%model%getParam()
584 err = getErrSampling(sampler, getLogLike, self%model%npar)
585 end select
586 if (err%occurred) error stop "sampler failed: "//err%msg
587
588 ! Find the mean best fit parameters and write post-prediction data to output file for visualization.
589
590 block
591
592 use pm_io, only: getErrTableRead
593 use pm_sysPath, only: glob, css_type
594
595 integer(IK) :: stat, ibest!, offset = 1
596 type(css_type), allocatable :: path(:)
597 real(RKG), allocatable :: logx(:), logPDF(:)
598 real(RKG), allocatable :: table(:,:), state(:)
599
600 write(*, "(A)") "Searching for files: "//self%sampler%outputFileName//SK_"*"
601 path = glob(self%sampler%outputFileName//SK_"*")
602 if (size(path) == 0) error stop "There is no sample file in the output folder."
603
604 stat = getErrTableRead(path(size(path))%val, table, roff = 1_IK)
605 if (stat /= 0) error stop "Failed to read output sample."
606 ibest = maxloc(table(:, 1), dim = 1_IK)
607 self%best%loglike = table(ibest, 1)
608 self%best%param = table(ibest, 2:)
609 self%bic = self%model%npar * log(real(self%data%nobs, RKG)) - 2 * self%best%loglike
610
611 end block
612
613 contains
614
615 function getLogLike(param) result(logLike)
616 real(RKG), intent(in), contiguous :: param(:)
617 real(RKG) :: logLike
618 call self%model%setParam(param)
619 loglike = sum(self%model%getLogPDF(self%data%obs))
620 end function
621
622 end subroutine
623
624end module fit_mod
625
626!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
627
628program example
629
630 use pm_kind, only: SK, IK
631 use fit_mod, only: fit_type
632 use model_mod, only: con_type
633 use auxil_mod, only: RKG, LARGE
634 use lognorm_mod, only: lognorm_type
635 use mixture_mod, only: mixture_type
636 use pm_sampling, only: paradram_type
637 use pm_arraySpace, only: getLogSpace
638 use flatPoweto_mod, only: flatPoweto_type
639 use data_mod, only: posuniv_type, data_type, stat_type
640 use flatPowetoTapered_mod, only: flatPowetoTapered_type
641 use pm_io, only: display_type, css_type, getErrTableWrite
642
643 integer(IK) :: imodel, stat
644 type(paradram_type) :: sampler
645 type(data_type) :: data, pred
646 type(display_type) :: disp
647 type(fit_type) :: fit
648
649 disp = display_type(file = SK_"main.out.F90")
650
651 ! read data.
652
653 block
654 use pm_io, only: getErrTableRead
655 real(RKG), allocatable :: table(:,:)
656 character(255, SK) :: iomsg
657 stat = getErrTableRead("data.csv", table, roff = 1_IK, iomsg = iomsg)
658 if (stat /= 0) error stop "Failed to read data: "//trim(iomsg)
659 associate(val => table(:, 4))
660 data = data_type(obs = posuniv_type(val), stat = stat_type([minval(val), maxval(val)]))
661 end associate
662 end block
663
664 ! perform fitting.
665
666 do imodel = 1, 4
667
668 sampler = paradram_type()
669
670 if (imodel == 1) then
671 sampler%outputFileName = "./mixLogNormLogNorm"
672 sampler%proposalStart = [-.2, 0.4, .4, 3.6, -0.2]
673 sampler%domainCubeLimitLower = [real(RKG) :: -LARGE, -LARGE, -LARGE, log(3.), -LARGE]
674 sampler%domainCubeLimitUpper = [real(RKG) :: log(3.), +LARGE, +LARGE, +LARGE, +LARGE]
675 !sampler%proposalScale = "0.1 * gelman"
676 sampler%proposalCov = reshape ( &
677 [ 1.6E-2, 5.5E-3, 4.1E-3, 3.9E-3,-2.5E-3 &
678 , 5.5E-3, 3.1E-3, 1.7E-3, 1.4E-3,-9.5E-4 &
679 , 4.1E-3, 1.7E-3, 1.9E-3, 1.2E-3,-8.8E-4 &
680 , 3.9E-3, 1.4E-3, 1.2E-3, 2.0E-3,-8.9E-4 &
681 , -2.5E-3,-9.5E-4,-8.8E-4,-8.9E-4, 1.0E-3 &
682 ], shape = [size(sampler%proposalStart), size(sampler%proposalStart)])
683 fit = fit_type(data, mixture_type([con_type(lognorm_type()), con_type(lognorm_type())]), sampler)
684 elseif (imodel == 2) then
685 sampler%proposalScale = "0.1 * gelman"
686 sampler%outputFileName = "./mixFlatPowetoFlatPowetoTapered"
687 sampler%proposalStart = [log(.37), -1.4, .3, log(21.3), -.5, 0.0156]
688 sampler%proposalCov = reshape ( &
689 [ 1.E-2, -3.E-3, 1.E-2, -7.E-3, -1.E-4, -6.E-4 &
690 , -3.E-3, 3.E-3, -4.E-3, 6.E-3, 1.E-4, 1.E-3 &
691 , 1.E-2, -4.E-3, 1.E-1, -7.E-2, -8.E-4, 3.E-4 &
692 , -7.E-3, 6.E-3, -7.E-2, 7.E-2, 8.E-4, 3.E-3 &
693 , -1.E-4, 1.E-4, -8.E-4, 8.E-4, 1.E-5, 6.E-5 &
694 , -6.E-4, 1.E-3, 3.E-4, 3.E-3, 6.E-5, 2.E-3 &
695 ], shape = [size(sampler%proposalStart), size(sampler%proposalStart)])
696 fit = fit_type(data, mixture_type([con_type(flatPoweto_type(data%stat%lim)), con_type(flatPowetoTapered_type(data%stat%lim))]), sampler)
697 elseif (imodel == 3) then
698 sampler%proposalScale = "0.1 * gelman"
699 sampler%outputFileName = "./mixLogNormFlatPowetoTapered"
700 sampler%proposalStart = [-.2, 0.4, .3, log(21.3), -.5, 0.0156]
701 fit = fit_type(data, mixture_type([con_type(lognorm_type()), con_type(flatPowetoTapered_type(data%stat%lim))]), sampler)
702 elseif (imodel == 4) then
703 sampler%proposalScale = "0.1 * gelman"
704 sampler%outputFileName = "./mixFlatPowetoLogNorm"
705 sampler%proposalStart = [log(.37), -1.4, .3, 3.6, -0.2]
706 fit = fit_type(data, mixture_type([con_type(flatPoweto_type(data%stat%lim)), con_type(lognorm_type())]), sampler)
707 else
708 sampler%outputFileName = "./logNorm"
709 fit = fit_type(data, lognorm_type([real(RKG) :: 1, 1]), sampler)
710 end if
711
712 call fit%run()
713
714 call disp%show(css_type(fit%sampler%domainAxisName))
715 call disp%show(fit%best%param)
716 call disp%skip()
717 call disp%show("[fit%best%loglike, fit%bic]")
718 call disp%show( [fit%best%loglike, fit%bic] )
719 call disp%skip()
720
721 ! generate synthetic data for histogram reconstruction.
722
723 associate(val => getLogSpace(logx1 = log(data%stat%lim(1)), logx2 = log(data%stat%lim(2)), count = 1000_IK))
724 pred = data_type(obs = posuniv_type(val), stat = stat_type([minval(val), maxval(val)]))
725 call fit%model%setParam(fit%best%param)
726 stat = getErrTableWrite(fit%sampler%outputFileName//SK_".hist", reshape([log(val), fit%model%getLogPDF(pred%obs)], [size(val), 2]), header = SK_"logx,logPDF")
727 if (stat /= 0) error stop "Failed to write the histogram visualization data."
728 end associate
729
730 end do
731
732end
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_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:76

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


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 at Austin

Definition at line 652 of file pm_sampling.F90.

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