The ParaMonte Documentation Website
Current view: top level - kernel - Statistics_mod.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: MPI Parallel Kernel - Code Coverage Report Lines: 941 1003 93.8 %
Date: 2021-01-08 13:07:16 Functions: 83 83 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       3             : !!!!
       4             : !!!!   MIT License
       5             : !!!!
       6             : !!!!   ParaMonte: plain powerful parallel Monte Carlo library.
       7             : !!!!
       8             : !!!!   Copyright (C) 2012-present, The Computational Data Science Lab
       9             : !!!!
      10             : !!!!   This file is part of the ParaMonte library.
      11             : !!!!
      12             : !!!!   Permission is hereby granted, free of charge, to any person obtaining a
      13             : !!!!   copy of this software and associated documentation files (the "Software"),
      14             : !!!!   to deal in the Software without restriction, including without limitation
      15             : !!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
      16             : !!!!   and/or sell copies of the Software, and to permit persons to whom the
      17             : !!!!   Software is furnished to do so, subject to the following conditions:
      18             : !!!!
      19             : !!!!   The above copyright notice and this permission notice shall be
      20             : !!!!   included in all copies or substantial portions of the Software.
      21             : !!!!
      22             : !!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
      23             : !!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
      24             : !!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
      25             : !!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
      26             : !!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
      27             : !!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
      28             : !!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
      29             : !!!!
      30             : !!!!   ACKNOWLEDGMENT
      31             : !!!!
      32             : !!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
      33             : !!!!   As per the ParaMonte library license agreement terms, if you use any parts of
      34             : !!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
      35             : !!!!   work (education/research/industry/development/...) by citing the ParaMonte
      36             : !!!!   library as described on this page:
      37             : !!!!
      38             : !!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
      39             : !!!!
      40             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      41             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      42             : 
      43             : !> \brief This module contains the classes and procedures for various statistical computations.
      44             : !> \author Amir Shahmoradi
      45             : 
      46             : module Statistics_mod
      47             : 
      48             :     use Constants_mod, only: RK, IK
      49             : 
      50             :     implicit none
      51             : 
      52             :     !logical, save :: paradramPrintEnabled = .false.
      53             :     !logical, save :: paradisePrintEnabled = .false.
      54             : 
      55             :     character(len=*), parameter :: MODULE_NAME = "@Statistics_mod"
      56             : 
      57             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      58             : 
      59             :     interface getSNormCDF
      60             :         module procedure :: getSNormCDF_SPR, getSNormCDF_DPR
      61             :     end interface getSNormCDF
      62             : 
      63             :     interface getBetaCDF
      64             :         module procedure :: getBetaCDF_SPR, getBetaCDF_DPR
      65             :     end interface getBetaCDF
      66             : 
      67             :     interface getBetaContinuedFraction
      68             :         module procedure :: getBetaContinuedFraction_SPR, getBetaContinuedFraction_DPR
      69             :     end interface getBetaContinuedFraction
      70             : 
      71             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      72             : 
      73             :     interface getMean
      74             :         module procedure :: getMean_2D
      75             :     end interface getMean
      76             : 
      77             :     interface flatten
      78             :         module procedure :: flatten_2D
      79             :     end interface flatten
      80             : 
      81             :     interface getNormData
      82             :         module procedure :: getNormData_2D, normalizeWeightedData_2d
      83             :     end interface getNormData
      84             : 
      85             :     interface getVariance
      86             :         module procedure :: getVariance_1D, getVariance_2D
      87             :     end interface getVariance
      88             : 
      89             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      90             : 
      91             :     interface getLogProbLogn
      92             :         module procedure :: getLogProbLognSP, getLogProbLognMP
      93             :     end interface getLogProbLogn
      94             : 
      95             :     interface getLogProbLognorm
      96             :         module procedure :: getLogProbLognSP, getLogProbLognMP
      97             :     end interface getLogProbLognorm
      98             : 
      99             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     100             : 
     101             :     interface getLogProbNormSP
     102             :         module procedure :: getLogProbNormSP_RK, getLogProbNormSP_CK
     103             :     end interface getLogProbNormSP
     104             : 
     105             :     interface getLogProbNormMP
     106             :         module procedure :: getLogProbNormMP_RK, getLogProbNormMP_CK
     107             :     end interface getLogProbNormMP
     108             : 
     109             :     interface getLogProbMVNSP
     110             :         module procedure :: getLogProbMVNSP_RK, getLogProbMVNSP_CK
     111             :     end interface getLogProbMVNSP
     112             : 
     113             :     interface getLogProbMVNMP
     114             :         module procedure :: getLogProbMVNMP_RK, getLogProbMVNMP_CK
     115             :     end interface getLogProbMVNMP
     116             : 
     117             :     interface getLogProbNorm
     118             :         module procedure :: getLogProbNormSP_RK, getLogProbNormMP_RK
     119             :         module procedure :: getLogProbNormSP_CK, getLogProbNormMP_CK
     120             :     end interface getLogProbNorm
     121             : 
     122             :     interface getLogProbMVN
     123             :         module procedure :: getLogProbMVNSP_RK, getLogProbMVNMP_RK
     124             :         module procedure :: getLogProbMVNSP_CK, getLogProbMVNMP_CK
     125             :     end interface getLogProbMVN
     126             : 
     127             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     128             : 
     129             :     interface getLogProbMixNorm
     130             :         module procedure :: getLogProbMixNormSP_RK, getLogProbMixNormSP_CK
     131             :         module procedure :: getLogProbMixNormMP_RK, getLogProbMixNormMP_CK
     132             :     end interface getLogProbMixNorm
     133             : 
     134             :     interface getLogProbMixMVN
     135             :         module procedure :: getLogProbMixMVNSP_RK, getLogProbMixMVNSP_CK
     136             :         module procedure :: getLogProbMixMVNMP_RK, getLogProbMixMVNMP_CK
     137             :     end interface getLogProbMixMVN
     138             : 
     139             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     140             : 
     141             :     interface getMahalSqSP
     142             :         module procedure :: getMahalSqSP_RK, getMahalSqSP_CK
     143             :     end interface getMahalSqSP
     144             : 
     145             :     interface getMahalSqMP
     146             :         module procedure :: getMahalSqMP_RK, getMahalSqMP_CK
     147             :     end interface getMahalSqMP
     148             : 
     149             :     interface getMahalSq
     150             :         module procedure :: getMahalSqSP_RK, getMahalSqMP_RK
     151             :         module procedure :: getMahalSqSP_CK, getMahalSqMP_CK
     152             :     end interface getMahalSq
     153             : 
     154             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     155             : 
     156             : contains
     157             : 
     158             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     159             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     160             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     161             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     162             : 
     163             : 
     164             :     !> \brief
     165             :     !> Return the square of Mahalanobis distance for a single point. The output is a scalar variable.
     166             :     !>
     167             :     !> \param[in]   nd          :   The number of dimensions of the input `Point`.
     168             :     !> \param[in]   MeanVec     :   The mean vector of the sample.
     169             :     !> \param[in]   InvCovMat   :   The inverse covariance matrix of the sample.
     170             :     !> \param[in]   Point       :   The `Point` whose distance from the sample is to computed.
     171             :     !>
     172             :     !> \return
     173             :     !> `mahalSq` : The Mahalanobis distance squared of the point from
     174             :     !> the sample characterized by the input `MeanVec` and `InvCovMat`.
     175          12 :     pure function getMahalSqSP_RK(nd,MeanVec,InvCovMat,Point) result(mahalSq)
     176             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     177             :         !DEC$ ATTRIBUTES DLLEXPORT :: getMahalSqSP_RK
     178             : #endif
     179             :         use Constants_mod, only: IK, RK
     180             :         implicit none
     181             :         integer(IK), intent(in) :: nd
     182             :         real(RK)   , intent(in) :: MeanVec(nd)
     183             :         real(RK)   , intent(in) :: InvCovMat(nd,nd)        ! Inverse of the covariance matrix
     184             :         real(RK)   , intent(in) :: Point(nd)               ! input data points
     185             :         real(RK)                :: mahalSq
     186          48 :         real(RK)                :: NormedPoint(nd)
     187          48 :         NormedPoint = Point - MeanVec
     188          48 :         mahalSq = dot_product( NormedPoint , matmul(InvCovMat,NormedPoint) )
     189          12 :     end function getMahalSqSP_RK
     190             : 
     191             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     192             : 
     193             :     !> \brief
     194             :     !> Return the square of Mahalanobis distances for an row-wise array of points.
     195             :     !>
     196             :     !> \param[in]   nd          :   The number of dimensions of the input `Point` array.
     197             :     !> \param[in]   np          :   The number of points in the input input `Point` array.
     198             :     !> \param[in]   MeanVec     :   The mean vector of length `nd` of the sample.
     199             :     !> \param[in]   InvCovMat   :   The inverse covariance matrix `(nd,nd)` of the sample.
     200             :     !> \param[in]   Point       :   The `(nd,np)` array of points whose distances squared from the sample are to computed.
     201             :     !>
     202             :     !> \return
     203             :     !> `MahalSq` : A vector of length `np` containing the squares of the Mahalanobis distances
     204             :     !> of the input points from the sample characterized by the input `MeanVec` and `InvCovMat`.
     205             :     !>
     206             :     !> \warning
     207             :     !> For the sake of preserving the purity and computational efficiency of the function,
     208             :     !> if the computation fails at any stage, the first element of output will be returned negative.
     209          60 :     pure function getMahalSqMP_RK(nd,np,MeanVec,InvCovMat,Point) result(MahalSq)
     210             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     211             :         !DEC$ ATTRIBUTES DLLEXPORT :: getMahalSqMP_RK
     212             : #endif
     213          12 :         use Constants_mod, only: IK, RK
     214             :         implicit none
     215             :         integer(IK), intent(in) :: nd,np
     216             :         real(RK), intent(in)    :: MeanVec(nd)
     217             :         real(RK), intent(in)    :: InvCovMat(nd,nd) ! Inverse of the covariance matrix
     218             :         real(RK), intent(in)    :: Point(nd,np)     ! input data points
     219             :         real(RK)                :: MahalSq(np)      ! function return
     220          48 :         real(RK)                :: NormedPoint(nd)
     221             :         integer(IK)             :: ip
     222          36 :         do ip = 1,np
     223          96 :             NormedPoint = Point(1:nd,ip) - MeanVec
     224          96 :             MahalSq(ip) = dot_product( NormedPoint , matmul(InvCovMat,NormedPoint) )
     225          36 :             if (MahalSq(ip)<0._RK) then
     226             :             ! LCOV_EXCL_START
     227             :                 MahalSq(1) = -1._RK
     228             :                 return
     229             :             end if
     230             :             ! LCOV_EXCL_STOP
     231             :         end do
     232          24 :     end function getMahalSqMP_RK
     233             : 
     234             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     235             : 
     236             :     !> \brief
     237             :     !> Return the square of Mahalanobis distance for a single complex point. The output is a scalar variable.
     238             :     !>
     239             :     !> \param[in]   nd          :   The number of dimensions of the input `Point`.
     240             :     !> \param[in]   MeanVec     :   The mean vector of the sample.
     241             :     !> \param[in]   InvCovMat   :   The inverse covariance matrix of the sample.
     242             :     !> \param[in]   Point       :   The `Point` whose distance from the sample is to computed.
     243             :     !>
     244             :     !> \return
     245             :     !> `mahalSq` : The Mahalanobis distance squared of the point from
     246             :     !> the sample characterized by the input `MeanVec` and `InvCovMat`.
     247          12 :     pure function getMahalSqSP_CK(nd,MeanVec,InvCovMat,Point) result(mahalSq)
     248             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     249             :         !DEC$ ATTRIBUTES DLLEXPORT :: getMahalSqSP_CK
     250             : #endif
     251          12 :         use Constants_mod, only: IK, RK, CK
     252             :         implicit none
     253             :         integer(IK), intent(in)  :: nd
     254             :         complex(CK), intent(in)  :: MeanVec(nd)
     255             :         complex(CK), intent(in)  :: InvCovMat(nd,nd) ! Inverse of the covariance matrix
     256             :         complex(CK), intent(in)  :: Point(nd)        ! input data points
     257             :         complex(CK)              :: mahalSq             ! function return
     258          84 :         mahalSq = sum( (Point-MeanVec) * matmul(InvCovMat,Point-MeanVec) )
     259          12 :     end function getMahalSqSP_CK
     260             : 
     261             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     262             : 
     263             :     !> \brief
     264             :     !> Return the square of Mahalanobis distances for an row-wise array of complex-valued points.
     265             :     !>
     266             :     !> \param[in]   nd          :   The number of dimensions of the input `Point` array.
     267             :     !> \param[in]   np          :   The number of points in the input input `Point` array.
     268             :     !> \param[in]   MeanVec     :   The mean vector of length `nd` of the sample.
     269             :     !> \param[in]   InvCovMat   :   The inverse covariance matrix `(nd,nd)` of the sample.
     270             :     !> \param[in]   Point       :   The `(nd,np)` array of points whose distances squared from the sample are to computed.
     271             :     !>
     272             :     !> \return
     273             :     !> `MahalSq` : A vector of length `np` containing the squares of the Mahalanobis distances
     274             :     !> of the input points from the sample characterized by the input `MeanVec` and `InvCovMat`.
     275             :     !>
     276             :     !> \warning
     277             :     !> For the sake of preserving the purity and computational efficiency of the function,
     278             :     !> if the computation fails at any stage, the first element of output will be returned negative.
     279          60 :     pure function getMahalSqMP_CK(nd,np,MeanVec,InvCovMat,Point) result(MahalSq)
     280             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     281             :         !DEC$ ATTRIBUTES DLLEXPORT :: getMahalSqMP_CK
     282             : #endif
     283          12 :         use Constants_mod, only: IK, RK, CK
     284             :         implicit none
     285             :         integer(IK), intent(in)  :: nd,np
     286             :         complex(CK), intent(in)  :: MeanVec(nd)
     287             :         complex(CK), intent(in)  :: InvCovMat(nd,nd) ! Inverse of the covariance matrix
     288             :         complex(CK), intent(in)  :: Point(nd,np)     ! input data points
     289             :         complex(CK)              :: MahalSq(np)         ! function return
     290             :         integer(IK)              :: ip
     291          36 :         do ip = 1,np
     292          48 :             MahalSq(ip) = sum( (Point(1:nd,ip)-MeanVec) * &
     293         192 :             matmul(InvCovMat,Point(1:nd,ip)-MeanVec) )
     294          36 :             if (real(MahalSq(ip))<0._RK) then
     295             :             ! LCOV_EXCL_START
     296             :                 MahalSq(1) = (-1._RK, -1._RK)
     297             :                 return
     298             :             end if
     299             :             ! LCOV_EXCL_STOP
     300             :         end do
     301          24 :     end function getMahalSqMP_CK
     302             : 
     303             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     304             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     305             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     306             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     307             : 
     308           9 :     pure function getLogProbNormSP_RK(mean,inverseVariance,logSqrtInverseVariance,point) result(logProbNorm)
     309             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     310             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbNormSP_RK
     311             : #endif
     312          12 :         use Constants_mod, only: RK, LOGINVSQRT2PI
     313             :         implicit none
     314             :         real(RK), intent(in) :: mean,inverseVariance,logSqrtInverseVariance,point
     315             :         real(RK)             :: logProbNorm
     316           9 :         logProbNorm = LOGINVSQRT2PI + logSqrtInverseVariance - 0.5_RK * inverseVariance * (point-mean)**2
     317          18 :     end function getLogProbNormSP_RK
     318             : 
     319             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     320             : 
     321          45 :     pure function getLogProbNormMP_RK(np,mean,inverseVariance,logSqrtInverseVariance,Point) result(logProbNorm)
     322             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     323             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbNormMP_RK
     324             : #endif
     325           9 :         use Constants_mod, only: LOGINVSQRT2PI
     326             :         implicit none
     327             :         integer(IK), intent(in) :: np
     328             :         real(RK)   , intent(in) :: mean,inverseVariance,logSqrtInverseVariance,Point(np)
     329             :         real(RK)                :: logProbNorm(np)
     330          27 :         logProbNorm = LOGINVSQRT2PI + logSqrtInverseVariance - 0.5_RK * inverseVariance * (Point-mean)**2
     331           9 :     end function getLogProbNormMP_RK
     332             : 
     333             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     334             : 
     335             :     ! NOTE: if MahalSq computation fails, output probability will be returned as NullVal%RK from module Constant_mod.
     336           9 :     pure function getLogProbMVNSP_RK(nd,MeanVec,InvCovMat,logSqrtDetInvCovMat,Point) result(logProbNorm)
     337             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     338             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMVNSP_RK
     339             : #endif
     340           9 :         use Constants_mod, only: LOGINVSQRT2PI, NullVal
     341             :         implicit none
     342             :         integer(IK), intent(in) :: nd
     343             :         real(RK)   , intent(in) :: MeanVec(nd)
     344             :         real(RK)   , intent(in) :: InvCovMat(nd,nd)
     345             :         real(RK)   , intent(in) :: logSqrtDetInvCovMat
     346             :         real(RK)   , intent(in) :: Point(nd)
     347           9 :         real(RK)                :: logProbNorm, dummy
     348           9 :         dummy = getMahalSqSP_RK(nd,MeanVec,InvCovMat,Point)
     349           9 :         if (dummy<0._RK) then
     350           0 :             logProbNorm = NullVal%RK
     351             :         else
     352           9 :             logProbNorm = nd*LOGINVSQRT2PI + logSqrtDetInvCovMat - 0.5_RK * dummy
     353             :         end if
     354           9 :     end function getLogProbMVNSP_RK
     355             : 
     356             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     357             : 
     358             :     ! NOTE: if MahalSq computation fails, output probability will be returned as NullVal%RK from module Constant_mod.
     359          45 :     pure function getLogProbMVNMP_RK(nd,np,MeanVec,InvCovMat,logSqrtDetInvCovMat,Point) result(logProbNorm)
     360             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     361             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMVNMP_RK
     362             : #endif
     363           9 :         use Constants_mod, only: LOGINVSQRT2PI, NullVal
     364             :         implicit none
     365             :         integer(IK), intent(in) :: nd,np
     366             :         real(RK)   , intent(in) :: MeanVec(nd)
     367             :         real(RK)   , intent(in) :: InvCovMat(nd,nd)
     368             :         real(RK)   , intent(in) :: logSqrtDetInvCovMat
     369             :         real(RK)   , intent(in) :: Point(nd,np)
     370          36 :         real(RK)                :: logProbNorm(np), Dummy(np)
     371           9 :         Dummy = getMahalSqMP_RK(nd,np,MeanVec,InvCovMat,Point)
     372           9 :         if (Dummy(1)<0._RK) then
     373           0 :             logProbNorm = NullVal%RK
     374             :         else
     375          27 :             logProbNorm = nd*LOGINVSQRT2PI + logSqrtDetInvCovMat - 0.5_RK * Dummy
     376             :         end if
     377           9 :     end function getLogProbMVNMP_RK
     378             : 
     379             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     380             : 
     381           9 :     function getLogProbNormSP_CK(mean,inverseVariance,logSqrtInverseVariance,point) result(logProbNorm)
     382             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     383             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbNormSP_CK
     384             : #endif
     385           9 :         use Constants_mod, only: RK, CK, LOGINVSQRT2PI
     386             :         implicit none
     387             :         complex(CK), intent(in) :: mean,inverseVariance,logSqrtInverseVariance,point
     388             :         complex(CK)             :: logProbNorm
     389           9 :         logProbNorm = LOGINVSQRT2PI + logSqrtInverseVariance - 0.5_RK * inverseVariance * (point-mean)**2
     390          18 :     end function getLogProbNormSP_CK
     391             : 
     392             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     393             : 
     394          45 :     function getLogProbNormMP_CK(np,mean,inverseVariance,logSqrtInverseVariance,Point) result(logProbNorm)
     395             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     396             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbNormMP_CK
     397             : #endif
     398           9 :         use Constants_mod, only: IK, RK, CK, LOGINVSQRT2PI
     399             :         implicit none
     400             :         integer(IK), intent(in) :: np
     401             :         complex(CK)   , intent(in) :: mean,inverseVariance,logSqrtInverseVariance,Point(np)
     402             :         complex(CK)                :: logProbNorm(np)
     403          27 :         logProbNorm = LOGINVSQRT2PI + logSqrtInverseVariance - 0.5_RK * inverseVariance * (Point-mean)**2
     404           9 :     end function getLogProbNormMP_CK
     405             : 
     406             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     407             : 
     408           9 :     function getLogProbMVNSP_CK(nd,MeanVec,InvCovMat,logSqrtDetInvCovMat,Point) result(logProbMVN)
     409             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     410             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMVNSP_CK
     411             : #endif
     412           9 :         use Constants_mod, only: IK, RK, CK, LOGINVSQRT2PI, NullVal
     413             :         implicit none
     414             :         integer(IK), intent(in) :: nd
     415             :         complex(CK), intent(in) :: MeanVec(nd)
     416             :         complex(CK), intent(in) :: InvCovMat(nd,nd)
     417             :         complex(CK), intent(in) :: logSqrtDetInvCovMat
     418             :         complex(CK), intent(in) :: Point(nd)
     419             :         complex(CK)             :: logProbMVN, dummy
     420           9 :         dummy = getMahalSqSP(nd,MeanVec,InvCovMat,Point)
     421           9 :         if (real(dummy)<0._RK) then
     422           0 :             logProbMVN = NullVal%RK
     423             :         else
     424           9 :             logProbMVN  = nd*LOGINVSQRT2PI + logSqrtDetInvCovMat - 0.5_RK * dummy
     425             :         end if
     426           9 :     end function getLogProbMVNSP_CK
     427             : 
     428             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     429             : 
     430          45 :     function getLogProbMVNMP_CK(nd,np,MeanVec,InvCovMat,logSqrtDetInvCovMat,Point) result(logProbMVN)
     431             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     432             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMVNMP_CK
     433             : #endif
     434           9 :         use Constants_mod, only: IK, RK, CK, LOGINVSQRT2PI, NullVal
     435             :         implicit none
     436             :         integer(IK), intent(in) :: nd,np
     437             :         complex(CK), intent(in) :: MeanVec(nd)
     438             :         complex(CK), intent(in) :: InvCovMat(nd,nd)
     439             :         complex(CK), intent(in) :: logSqrtDetInvCovMat
     440             :         complex(CK), intent(in) :: Point(nd,np)
     441          36 :         complex(CK)             :: logProbMVN(np), Dummy(np)
     442           9 :         Dummy = getMahalSqMP(nd,np,MeanVec,InvCovMat,Point)
     443           9 :         if (real(Dummy(1))<0._RK) then
     444           0 :             logProbMVN = NullVal%RK
     445             :         else
     446          27 :             logProbMVN  = nd*LOGINVSQRT2PI + logSqrtDetInvCovMat - 0.5_RK * Dummy
     447             :         end if
     448           9 :     end function getLogProbMVNMP_CK
     449             : 
     450             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     451             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     452             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     453             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     454             : 
     455             :     ! SDSP stands for Single-Dimensional Gaussian mixture with Single Point input
     456             :     ! For a proper probability normalization, the sum of the amplitudes must equal one.
     457           3 :     function getLogProbMixNormSP_RK(nmode,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,point) result(logProbMixNorm)
     458             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     459             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixNormSP_RK
     460             : #endif
     461           9 :         use Constants_mod, only: IK, RK, LOGTINY_RK
     462             :         implicit none
     463             :         integer(IK), intent(in) :: nmode
     464             :         real(RK)   , intent(in) :: LogAmplitude(nmode),MeanVec(nmode)
     465             :         real(RK)   , intent(in) :: InvCovMat(nmode),LogSqrtDetInvCovMat(nmode)
     466             :         real(RK)   , intent(in) :: point
     467             :         real(RK)                :: logProbMixNorm
     468           9 :         real(RK)                :: normFac,LogProb(nmode)
     469             :         integer(IK)             :: imode
     470           9 :         do imode = 1, nmode
     471           9 :             LogProb(imode)  = LogAmplitude(imode) + getLogProbNormSP_RK(MeanVec(imode),InvCovMat(imode),LogSqrtDetInvCovMat(imode),point)
     472             :         end do
     473          12 :         normFac = maxval(LogProb)
     474           9 :         LogProb = LogProb - normFac
     475          27 :         where(LogProb<LOGTINY_RK)
     476           3 :             LogProb = 0._RK
     477             :         elsewhere
     478           3 :             LogProb = exp(LogProb)
     479             :         end where
     480           9 :         logProbMixNorm = normFac + log(sum(LogProb))
     481           3 :     end function getLogProbMixNormSP_RK
     482             : 
     483             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     484             : 
     485          15 :   function getLogProbMixNormMP_RK(nmode,np,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,Point) result(logProbMixNorm)
     486             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     487             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixNormMP_RK
     488             : #endif
     489           3 :         use Constants_mod, only: IK, RK, LOGTINY_RK
     490             :         implicit none
     491             :         integer(IK), intent(in) :: nmode,np
     492             :         real(RK)   , intent(in) :: LogAmplitude(nmode),MeanVec(nmode)
     493             :         real(RK)   , intent(in) :: InvCovMat(nmode),LogSqrtDetInvCovMat(nmode)
     494             :         real(RK)   , intent(in) :: Point(np)
     495             :         real(RK)                :: logProbMixNorm(np)
     496          30 :         real(RK)                :: NormFac(np),LogProb(nmode,np)
     497             :         integer(IK)             :: imode, ip
     498           9 :         do imode = 1, nmode
     499          21 :             LogProb(imode,1:np) = LogAmplitude(imode) + getLogProbNormMP_RK(np,MeanVec(imode),InvCovMat(imode),LogSqrtDetInvCovMat(imode),Point)
     500             :         end do
     501           3 :         NormFac = maxval(LogProb,dim=1)
     502           9 :         do ip = 1,np
     503          18 :             LogProb(1:nmode,ip) = LogProb(1:nmode,ip) - NormFac(ip)
     504          18 :             do imode = 1,nmode
     505          18 :                 if ( LogProb(imode,ip) < LOGTINY_RK ) then
     506           0 :                     LogProb(imode,ip) = 0._RK
     507             :                 else
     508          12 :                     LogProb(imode,ip) = exp( LogProb(imode,ip) )
     509             :                 end if
     510             :             end do
     511          21 :             logProbMixNorm(ip) = NormFac(ip) + log(sum(LogProb(1:nmode,ip)))
     512             :         end do
     513           3 :     end function getLogProbMixNormMP_RK
     514             : 
     515             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     516             : 
     517           3 :     function getLogProbMixMVNSP_RK(nmode,nd,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,Point) result(logProbMixMVN)
     518             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     519             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixMVNSP_RK
     520             : #endif
     521           3 :         use Constants_mod, only: IK, RK, LOGTINY_RK
     522             :         implicit none
     523             :         integer(IK), intent(in) :: nmode,nd
     524             :         real(RK)   , intent(in) :: LogAmplitude(nmode), MeanVec(nd,nmode)
     525             :         real(RK)   , intent(in) :: InvCovMat(nd,nd,nmode), LogSqrtDetInvCovMat(nmode)
     526             :         real(RK)   , intent(in) :: Point(nd)
     527             :         real(RK)                :: logProbMixMVN
     528           9 :         real(RK)                :: normFac,LogProb(nmode)
     529             :         integer(IK)             :: imode
     530           9 :         do imode = 1, nmode
     531           9 :             LogProb(imode) = LogAmplitude(imode) + getLogProbMVNSP_RK(nd,MeanVec(1:nd,imode),InvCovMat(1:nd,1:nd,imode),LogSqrtDetInvCovMat(imode),Point)
     532             :         end do
     533          12 :         normFac = maxval(LogProb)
     534           9 :         LogProb = LogProb - normFac
     535          27 :         where(LogProb<LOGTINY_RK)
     536           3 :             LogProb = 0._RK
     537             :         elsewhere
     538           3 :             LogProb = exp(LogProb)
     539             :         end where
     540           9 :         logProbMixMVN = normFac + log(sum(LogProb))
     541           3 :     end function getLogProbMixMVNSP_RK
     542             : 
     543             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     544             : 
     545          15 :   function getLogProbMixMVNMP_RK(nmode,nd,np,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,Point) result(logProbMixMVN)
     546             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     547             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixMVNMP_RK
     548             : #endif
     549           3 :         use Constants_mod, only: IK, RK, LOGTINY_RK
     550             :         implicit none
     551             :         integer(IK), intent(in) :: nmode,nd,np
     552             :         real(RK)   , intent(in) :: LogAmplitude(nmode),MeanVec(nd,nmode)
     553             :         real(RK)   , intent(in) :: InvCovMat(nd,nd,nmode), LogSqrtDetInvCovMat(nmode)
     554             :         real(RK)   , intent(in) :: Point(nd,np)
     555             :         real(RK)                :: logProbMixMVN(np)
     556          30 :         real(RK)                :: NormFac(np),LogProb(nmode,np)
     557             :         integer(IK)             :: imode, ip
     558           9 :         do imode = 1, nmode
     559          12 :             LogProb(imode,1:np) = LogAmplitude(imode) + &
     560          33 :             getLogProbMVNMP_RK(nd,np,MeanVec(1:nd,imode),InvCovMat(1:nd,1:nd,imode),LogSqrtDetInvCovMat(imode),Point)
     561             :         end do
     562           3 :         NormFac = maxval(LogProb,dim=1)
     563           9 :         do ip = 1,np
     564          18 :             LogProb(1:nmode,ip) = LogProb(1:nmode,ip) - NormFac(ip)
     565          18 :             do imode = 1,nmode
     566          18 :                 if ( LogProb(imode,ip)<LOGTINY_RK ) then
     567           6 :                     LogProb(imode,ip) = 0._RK
     568             :                 else
     569           6 :                     LogProb(imode,ip) = exp( LogProb(imode,ip) )
     570             :                 end if
     571             :             end do
     572          21 :         logProbMixMVN(ip) = NormFac(ip) + log(sum(LogProb(1:nmode,ip)))
     573             :         end do
     574           3 :     end function getLogProbMixMVNMP_RK
     575             : 
     576             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     577             : 
     578             :     ! SDSP stands for 1-dimensional Gaussian mixture with scalar input point
     579           3 :     function getLogProbMixNormSP_CK(nmode,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,point) result(logProbMixNorm)
     580             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     581             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixNormSP_CK
     582             : #endif
     583           3 :         use Constants_mod, only: IK, RK, CK, LOGTINY_RK
     584             :         implicit none
     585             :         integer(IK), intent(in) :: nmode
     586             :         complex(CK), intent(in) :: LogAmplitude(nmode),MeanVec(nmode)
     587             :         complex(CK), intent(in) :: InvCovMat(nmode),LogSqrtDetInvCovMat(nmode)
     588             :         complex(CK), intent(in) :: point
     589             :         complex(CK)             :: logProbMixNorm
     590           9 :         complex(CK)             :: normFac,LogProb(nmode)
     591             :         integer(IK)             :: imode
     592           9 :         do imode = 1, nmode
     593           9 :             LogProb(imode) = LogAmplitude(imode) + getLogProbNorm(MeanVec(imode),InvCovMat(imode),LogSqrtDetInvCovMat(imode),point)
     594             :         end do
     595          12 :         normFac = maxval(real(LogProb))
     596           9 :         LogProb = LogProb - normFac
     597          27 :         where(real(LogProb)<LOGTINY_RK)
     598           3 :             LogProb = 0._RK
     599             :         elsewhere
     600           3 :             LogProb = exp(LogProb)
     601             :         end where
     602           9 :         logProbMixNorm = normFac + log(sum(LogProb))
     603           3 :     end function getLogProbMixNormSP_CK
     604             : 
     605             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     606             : 
     607          15 :     function getLogProbMixNormMP_CK(nmode,np,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,Point) result(logProbMixNorm)
     608             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     609             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixNormMP_CK
     610             : #endif
     611           3 :         use Constants_mod, only: IK, RK, CK, LOGTINY_RK
     612             :         implicit none
     613             :         integer(IK), intent(in) :: nmode,np
     614             :         complex(CK), intent(in) :: LogAmplitude(nmode),MeanVec(nmode)
     615             :         complex(CK), intent(in) :: InvCovMat(nmode),LogSqrtDetInvCovMat(nmode)
     616             :         complex(CK), intent(in) :: Point(np)
     617             :         complex(CK)             :: logProbMixNorm(np)
     618          30 :         complex(CK)             :: normFac(np),LogProb(nmode,np)
     619             :         integer(IK)             :: imode, ip
     620           9 :         do imode = 1, nmode
     621          21 :             LogProb(imode,1:np) = LogAmplitude(imode) + getLogProbNorm(np,MeanVec(imode),InvCovMat(imode),LogSqrtDetInvCovMat(imode),Point)
     622             :         end do
     623          27 :         normFac = maxval(real(LogProb),dim=1)
     624           9 :         do ip = 1,np
     625          18 :             LogProb(1:nmode,ip) = LogProb(1:nmode,ip) - normFac(ip)
     626          18 :             do imode = 1,nmode
     627          18 :                 if ( real(LogProb(imode,ip)) < LOGTINY_RK ) then
     628           0 :                     LogProb(imode,ip) = 0._RK
     629             :                 else
     630          12 :                     LogProb(imode,ip) = exp( LogProb(imode,ip) )
     631             :                 end if
     632             :             end do
     633          21 :             logProbMixNorm(ip) = normFac(ip) + log(sum(LogProb(1:nmode,ip)))
     634             :         end do
     635           3 :     end function getLogProbMixNormMP_CK
     636             : 
     637             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     638             : 
     639           3 :   function getLogProbMixMVNSP_CK(nmode,nd,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,Point) result(logProbMixMVN)
     640             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     641             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixMVNSP_CK
     642             : #endif
     643           3 :         use Constants_mod, only: IK, RK, CK, LOGTINY_RK
     644             :         implicit none
     645             :         integer(IK), intent(in) :: nmode,nd
     646             :         complex(CK), intent(in) :: LogAmplitude(nmode), MeanVec(nd,nmode)
     647             :         complex(CK), intent(in) :: InvCovMat(nd,nd,nmode), LogSqrtDetInvCovMat(nmode)
     648             :         complex(CK), intent(in) :: Point(nd)
     649             :         complex(CK)             :: logProbMixMVN
     650           9 :         complex(CK)             :: normFac,LogProb(nmode)
     651             :         integer(IK)             :: imode
     652           9 :         do imode = 1, nmode
     653           9 :             LogProb(imode) = LogAmplitude(imode) + getLogProbMVN(nd,MeanVec(1:nd,imode),InvCovMat(1:nd,1:nd,imode),LogSqrtDetInvCovMat(imode),Point)
     654             :         end do
     655          12 :         normFac = maxval(real(LogProb))
     656           9 :         LogProb = LogProb - normFac
     657          27 :         where(real(LogProb)<LOGTINY_RK)
     658           3 :             LogProb = 0._RK
     659             :         elsewhere
     660           3 :             LogProb = exp(LogProb)
     661             :         end where
     662           9 :         logProbMixMVN = normFac + log(sum(LogProb))
     663           3 :     end function getLogProbMixMVNSP_CK
     664             : 
     665             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     666             : 
     667          15 :     function getLogProbMixMVNMP_CK(nmode,nd,np,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,Point) result(logProbMixMVN)
     668             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     669             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixMVNMP_CK
     670             : #endif
     671           3 :         use Constants_mod, only: IK, RK, CK, LOGTINY_RK
     672             :         implicit none
     673             :         integer(IK), intent(in) :: nmode,nd,np
     674             :         complex(CK), intent(in) :: LogAmplitude(nmode),MeanVec(nd,nmode)
     675             :         complex(CK), intent(in) :: InvCovMat(nd,nd,nmode), LogSqrtDetInvCovMat(nmode)
     676             :         complex(CK), intent(in) :: Point(nd,np)
     677             :         complex(CK)             :: logProbMixMVN(np)
     678          30 :         complex(CK)             :: normFac(np),LogProb(nmode,np)
     679             :         integer(IK)             :: imode, ip
     680           9 :         do imode = 1, nmode
     681          21 :             LogProb(imode,1:np) = LogAmplitude(imode) + getLogProbMVN(nd,np,MeanVec(1:nd,imode),InvCovMat(1:nd,1:nd,imode),LogSqrtDetInvCovMat(imode),Point)
     682             :         end do
     683          27 :         normFac = maxval(real(LogProb),dim=1)
     684           9 :         do ip = 1,np
     685          18 :             LogProb(1:nmode,ip) = LogProb(1:nmode,ip) - normFac(ip)
     686          18 :             do imode = 1,nmode
     687          18 :                 if ( real(LogProb(imode,ip))<LOGTINY_RK ) then
     688           6 :                     LogProb(imode,ip) = 0._RK
     689             :                 else
     690           6 :                     LogProb(imode,ip) = exp( LogProb(imode,ip) )
     691             :                 end if
     692             :             end do
     693          21 :             logProbMixMVN(ip) = normFac(ip) + log(sum(LogProb(1:nmode,ip)))
     694             :         end do
     695           3 :     end function getLogProbMixMVNMP_CK
     696             : 
     697             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     698             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     699             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     700             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     701             : 
     702             :     !include "Statistics_mod@MahalSq_RK.inc.f90"
     703             :     !include "Statistics_mod@MahalSq_CK.inc.f90"
     704             :     !include "Statistics_mod@LogProbGaus_RK.inc.f90"
     705             :     !include "Statistics_mod@LogProbGaus_CK.inc.f90"
     706             :     !include "Statistics_mod@LogProbGausMix_RK.inc.f90"
     707             :     !include "Statistics_mod@LogProbGausMix_CK.inc.f90"
     708             : 
     709             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     710             : 
     711             :     !> \brief
     712             :     !> Flatten the input `Point` array such that each element of the output 
     713             :     !> `FlattenedPoint` array has the same unity weight as elements in the array.
     714             :     !>
     715             :     !> \param[in]       nd      :   The number of dimensions of the input sample.
     716             :     !> \param[in]       np      :   The number of points in the sample.
     717             :     !> \param[in]       Point   :   The array of shape `(nd,np)` containing the sample.
     718             :     !> \param[in]       Weight  :   The vector of length `np` containing the weights of points in the sample.
     719             :     !>                              The values of elements of Weight are allowed to be negative, in which case,
     720             :     !>                              the corresponding elements will be excluded from the output `FlattenedPoint`.
     721             :     !>
     722             :     !> \warning
     723             :     !> Note the shape of the input argument `Point(nd,np)`.
     724             :     !>
     725             :     !> \return
     726             :     !> `FlattenedPoint` : The flattened array whose elements all have the same weight.
     727          11 :     pure function flatten_2D(nd,np,Point,Weight) result(FlattenedPoint)
     728             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     729             :         !DEC$ ATTRIBUTES DLLEXPORT :: flatten_2D
     730             : #endif
     731             :         ! the implementation for one-dimension is very concise and nice: mean = sum(Weight*Point) / sum(Weight)
     732             :         implicit none
     733             :         integer(IK), intent(in) :: np,nd            ! np: number of observations, nd: number of variables for each observation
     734             :         real(RK)   , intent(in) :: Point(nd,np)     ! Point is the data matrix
     735             :         integer(IK), intent(in) :: Weight(np)       ! sample weight
     736             :         integer(IK)             :: ip, iweight, sumWeight, counter
     737             :         real(RK), allocatable   :: FlattenedPoint(:,:)
     738          11 :         sumWeight = 0_IK
     739         103 :         do ip = 1, np
     740         103 :             if (Weight(ip)>0_IK) sumWeight = sumWeight + Weight(ip)
     741             :         end do
     742          11 :         allocate(FlattenedPoint(nd,sumWeight))
     743          11 :         counter = 0_IK
     744         103 :         do ip = 1, np
     745         192 :             do iweight = 1, Weight(ip)
     746          89 :                 counter = counter + 1_IK
     747         439 :                 FlattenedPoint(1:nd,counter) = Point(1:nd,ip)
     748             :             end do
     749             :         end do
     750           3 :     end function flatten_2D
     751             : 
     752             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     753             : 
     754             :     !> \brief
     755             :     !> Return the mean of a sample of multidimensional points.
     756             :     !>
     757             :     !> \param[in]       nd      :   The number of dimensions of the input sample.
     758             :     !> \param[in]       np      :   The number of points in the sample.
     759             :     !> \param[in]       Point   :   The array of shape `(nd,np)` containing the sample.
     760             :     !> \param[in]       Weight  :   The vector of length `np` containing the weights of points in the sample (**optional**, default = vector of ones).
     761             :     !>
     762             :     !> \warning
     763             :     !> Note the shape of the input argument `Point(nd,np)`.
     764             :     !>
     765             :     !> \return
     766             :     !> `Mean` : The output mean vector of length `nd`.
     767          71 :     pure function getMean_2D(nd,np,Point,Weight) result(Mean)
     768             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     769             :         !DEC$ ATTRIBUTES DLLEXPORT :: getMean_2D
     770             : #endif
     771             :         ! the implementation for one-dimension is very concise and nice: mean = sum(Weight*Point) / sum(Weight)
     772             :         implicit none
     773             :         integer(IK), intent(in)             :: np,nd            ! np: number of observations, nd: number of variables for each observation
     774             :         real(RK)   , intent(in)             :: Point(nd,np)     ! Point is the data matrix
     775             :         integer(IK), intent(in), optional   :: Weight(np)       ! sample weight
     776             :         real(RK)                            :: Mean(nd)         ! output mean vector
     777             :         integer(IK)                         :: ip, sumWeight
     778          45 :         Mean = 0._RK
     779          13 :         if (present(Weight)) then
     780          10 :             sumWeight = 0_IK
     781          43 :             do ip = 1,np
     782          33 :                 sumWeight = sumWeight + Weight(ip)
     783         124 :                 Mean = Mean + Weight(ip) * Point(1:nd,ip)
     784             :             end do
     785          33 :             Mean = Mean / sumWeight
     786             :         else
     787          18 :             do ip = 1,np
     788          63 :                 Mean = Mean + Point(1:nd,ip)
     789             :             end do
     790          12 :             Mean = Mean / real(np,kind=RK)
     791             :         end if
     792          24 :     end function getMean_2D
     793             : 
     794             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     795             : 
     796             :     !> \brief
     797             :     !> Return the normalized data with respect to the input mean vector of length `nd`.
     798             :     !>
     799             :     !> \param[in]       nd      :   The number of dimensions of the input sample.
     800             :     !> \param[in]       np      :   The number of points in the sample.
     801             :     !> \param[in]       Mean    :   The mean vector of length `nd`.
     802             :     !> \param[in]       Point   :   The array of shape `(nd,np)` containing the sample.
     803             :     !>
     804             :     !> \return
     805             :     !> `NormData` : The output normalized points array of shape `(np,nd)`.
     806             :     !>
     807             :     !> \remark
     808             :     !> Note the difference in the shape of the input `Point` vs. the output `NormData`.
     809          63 :     pure function getNormData_2D(nd,np,Mean,Point) result(NormData)
     810             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     811             :         !DEC$ ATTRIBUTES DLLEXPORT :: getNormData_2D
     812             : #endif
     813             :         implicit none
     814             :         integer(IK), intent(in)  :: np,nd           ! np is the number of observations, nd is the number of parameters for each observation
     815             :         real(RK)   , intent(in)  :: Mean(nd)        ! Mean vector
     816             :         real(RK)   , intent(in)  :: Point(nd,np)    ! Point is the matrix of the data, CovMat contains the elements of the sample covariance matrix
     817             :         real(RK)                 :: NormData(np,nd)
     818             :         integer(IK)              :: i
     819          18 :         do i = 1,np
     820          63 :             NormData(i,1:nd) = Point(1:nd,i) - Mean
     821             :         end do
     822          16 :     end function getNormData_2D
     823             : 
     824             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     825             : 
     826             :     !> \brief
     827             :     !> Return the normalized 2D data of size `(nd,np)` with respect to the mean of the data along the second dimension of length `np`.
     828             :     !>
     829             :     !> \param[in]   nd          :   The length of the input `Data` matrix along the first dimension.
     830             :     !> \param[in]   np          :   The length of the input `Data` matrix along the second dimension.
     831             :     !> \param[in]   Data        :   The input data series data vector.
     832             :     !> \param[in]   Weight      :   The vector of weights of the input data points (**optional**, default = array of ones).
     833             :     !> \param[in]   tenabled    :   A logical value that, if `.true.` will cause the output `NormData` to have transposed shape
     834             :     !>                              of the input `Point(nd,np)` matrix, that is `(np,nd)` (**optional**, default = `.false.`).
     835             :     !>
     836             :     !> \return
     837             :     !> `NormData` : The integrated autocorrelation (IAC) via the BatchMeans method.
     838             :     !>
     839             :     !> \remark
     840             :     !> Note that np must be large enough to get a meaningful answer.
     841          57 :     pure function normalizeWeightedData_2D(nd, np, Data, Weight, tenabled) result(NormData)
     842             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     843             :         !DEC$ ATTRIBUTES DLLEXPORT :: normalizeWeightedData_2D
     844             : #endif
     845           3 :         use Constants_mod, only: IK, RK
     846             :         implicit none
     847             :         integer(IK) , intent(in)            :: nd, np
     848             :         real(RK)    , intent(in)            :: Data(nd,np)
     849             :         integer(IK) , intent(in), optional  :: Weight(np)
     850             :         logical     , intent(in), optional  :: tenabled
     851             :         real(RK)    , allocatable           :: NormData(:,:)
     852             :         logical                             :: tenabledDefault
     853         153 :         real(RK)                            :: Mean(nd)
     854          51 :         real(RK)                            :: sumWeight
     855             :         integer(IK)                         :: id, ip
     856             : 
     857          51 :         tenabledDefault = .false.
     858          51 :         if (present(tenabled)) tenabledDefault = tenabled
     859             : 
     860         102 :         Mean = 0._RK
     861          51 :         if (present(Weight)) then
     862           6 :             sumWeight = 0._RK
     863       30261 :             do ip = 1, np
     864       30255 :                 sumWeight = sumWeight + Weight(ip)
     865       60516 :                 Mean = Mean + Weight(ip) * Data(1:nd,ip)
     866             :             end do
     867             :         else
     868          45 :             sumWeight = np
     869      449370 :             do ip = 1, np
     870      898695 :                 Mean = Mean + Data(1:nd,ip)
     871             :             end do
     872             :         end if
     873         102 :         Mean = Mean / sumWeight
     874             : 
     875          51 :         if (tenabledDefault) then
     876           6 :             allocate(NormData(np,nd))
     877           6 :             do concurrent(id = 1:nd)
     878       30267 :                 NormData(1:np,id) = Data(id,1:np) - Mean(id)
     879             :             end do
     880             :         else
     881          45 :             allocate(NormData(nd,np))
     882          45 :             do concurrent(id = 1:nd)
     883      449415 :                 NormData(id,1:np) = Data(id,1:np) - Mean(id)
     884             :             end do
     885             :         end if
     886             : 
     887         102 :     end function normalizeWeightedData_2D
     888             : 
     889             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     890             : 
     891             :     !> \brief
     892             :     !> Return the variance of the input vector of values of length `np`.
     893             :     !>
     894             :     !> \param[in]       np          :   The number of points in the sample.
     895             :     !> \param[in]       mean        :   The mean of the vector.
     896             :     !> \param[in]       Point       :   The array of shape `np` containing the sample.
     897             :     !> \param[in]       Weight      :   The vector of weights of the points in the sample (**optional**).
     898             :     !> \param[in]       sumWeight   :   The sum of the vector of weights (**optional**, if `Weight` is also missing).
     899             :     !>
     900             :     !> \return
     901             :     !> `variance` : The variance of the input sample.
     902             :     !>
     903             :     !> \warning
     904             :     !> If `Weight` is present, then `sumWeight` must be also present.
     905             :     !> Why? if mean is already given, it implies that sumWeight is also computed a priori.
     906             :     !>
     907             :     !> \remark
     908             :     !> One also use the concise Fortran syntax to achieve the same goal as this function:
     909             :     !> ```
     910             :     !> mean = sum(Weight*Point) / sum(Weight)
     911             :     !> variance = sum( (Weight*(Point-mean))**2 ) / (sum(Weight)-1)
     912             :     !> ```
     913             :     !> But the above concise version will be slightly slower as it involves three loops instead of two.
     914           9 :     pure function getVariance_1D(np,mean,Point,Weight,sumWeight) result(variance)
     915             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     916             :         !DEC$ ATTRIBUTES DLLEXPORT :: getVariance_1D
     917             : #endif
     918             :         implicit none
     919             :         integer(IK), intent(in)             :: np                       ! np is the number of observations (points) whose variance is to be computed
     920             :         real(RK)   , intent(in)             :: mean                     ! the mean value of the vector Point
     921             :         real(RK)   , intent(in)             :: Point(np)                ! Point is the vector of data
     922             :         integer(IK), intent(in), optional   :: Weight(np), sumWeight    ! sample weight
     923             :         real(RK)                            :: variance                 ! output mean vector
     924             :         integer(IK)                         :: ip
     925           6 :         variance = 0._RK
     926           6 :         if (present(Weight)) then
     927          48 :             do ip = 1,np
     928          48 :                 variance = variance + Weight(ip) * ( Point(ip) - mean )**2
     929             :             end do
     930           3 :             variance = variance / real(sumWeight-1_IK,kind=RK)
     931             :         else
     932          48 :             do ip = 1,np
     933          48 :                 variance = variance + ( Point(ip) - mean )**2
     934             :             end do
     935           3 :             variance = variance / real(np-1_IK,kind=RK)
     936             :         end if
     937          57 :     end function getVariance_1D
     938             : 
     939             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     940             : 
     941             :     !> \brief
     942             :     !> Return the vector of variance of a set of `np` points of `nd` dimensions.
     943             :     !>
     944             :     !> \param[in]       nd          :   The number of dimensions of the input sample.
     945             :     !> \param[in]       np          :   The number of points in the sample.
     946             :     !> \param[in]       Mean        :   The mean vector of the sample.
     947             :     !> \param[in]       Point       :   The array of shape `(nd,np)` containing the sample.
     948             :     !> \param[in]       Weight      :   The vector of weights of the points in the sample of shape `(nd,np)` (**optional**).
     949             :     !>
     950             :     !> \return
     951             :     !> `Variance` : The vector of length `nd` of the variances of the input sample.
     952          36 :     pure function getVariance_2D(nd,np,Mean,Point,Weight) result(Variance)
     953             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     954             :         !DEC$ ATTRIBUTES DLLEXPORT :: getVariance_2D
     955             : #endif
     956             :         ! returns the variance of each row in Point weighted by the corresponding Weight if provided
     957             :         ! pass the Mean vector by calling getMean() or getMean_2D()
     958             :         implicit none
     959             :         integer(IK), intent(in)             :: nd, np           ! np is the number of observations (points) whose variance is to be computed
     960             :         real(RK)   , intent(in)             :: Mean(nd)         ! the Mean value of the vector Point
     961             :         real(RK)   , intent(in)             :: Point(nd,np)     ! Point is the vector of data
     962             :         integer(IK), intent(in), optional   :: Weight(np)       ! sample weight
     963             :         real(RK)                            :: Variance(nd)     ! output Mean vector
     964             :         integer(IK)                         :: ip, sumWeight
     965          24 :         Variance = 0._RK
     966           6 :         if (present(Weight)) then
     967           3 :             sumWeight = 0_IK
     968          18 :             do ip = 1,np
     969          15 :                 sumWeight = sumWeight + Weight(ip)
     970          63 :                 Variance = Variance + Weight(ip) * ( Point(1:nd,ip) - Mean )**2
     971             :             end do
     972          12 :             Variance = Variance / real(sumWeight-1_IK,kind=RK)
     973             :         else
     974          18 :             do ip = 1,np
     975          63 :                 Variance = Variance + ( Point(1:nd,ip) - Mean )**2
     976             :             end do
     977          12 :             Variance = Variance / real(np-1_IK,kind=RK)
     978             :         end if
     979          12 :     end function getVariance_2D
     980             : 
     981             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     982             : 
     983             :     !> \brief
     984             :     !> Return the lower triangle Cholesky Factor of the covariance matrix of a set of points in the lower part of `CholeskyLower`.
     985             :     ! The upper part of `CholeskyLower`, including the diagonal elements of it, will contain the covariance matrix of the sample.
     986             :     ! The output argument `Diagonal`, contains the diagonal terms of Cholesky Factor.
     987             :     !>
     988             :     !> \param[in]       nd              :   The number of dimensions of the input sample.
     989             :     !> \param[in]       np              :   The number of points in the sample.
     990             :     !> \param[in]       Mean            :   The mean vector of the sample.
     991             :     !> \param[in]       Point           :   The array of shape `(nd,np)` containing the sample.
     992             :     !> \param[out]      CholeskyLower   :   The output matrix of shape `(nd,nd)` whose lower triangle contains elements of the Cholesky factor.
     993             :     !>                                      The upper triangle of the matrix contains the covariance matrix of the sample.
     994             :     !> \param[out]      Diagonal        :   The diagonal elements of the Cholesky factor.
     995             :     !>
     996             :     !> \todo
     997             :     !> The efficiency of this code can be further improved.
     998           3 :     subroutine getSamCholFac(nd,np,Mean,Point,CholeskyLower,Diagonal)
     999             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1000             :         !DEC$ ATTRIBUTES DLLEXPORT :: getSamCholFac
    1001             : #endif
    1002           6 :         use Matrix_mod, only: getCholeskyFactor
    1003             :         implicit none
    1004             :         integer(IK), intent(in)  :: nd,np                  ! np is the number of observations, nd is the number of parameters for each observation
    1005             :         real(RK)   , intent(in)  :: Mean(nd)               ! Mean vector
    1006             :         real(RK)   , intent(in)  :: Point(nd,np)           ! Point is the matrix of the data, CovMat contains the elements of the sample covariance matrix
    1007             :         real(RK)   , intent(out) :: CholeskyLower(nd,nd)   ! Lower Cholesky Factor of the covariance matrix
    1008             :         real(RK)   , intent(out) :: Diagonal(nd)           ! Diagonal elements of the Cholesky factorization
    1009          57 :         real(RK)                 :: NormData(np,nd), npMinusOneInverse
    1010             :         integer(IK)              :: i,j
    1011             : 
    1012          18 :         do i = 1,np
    1013          63 :             NormData(i,1:nd) = Point(1:nd,i) - Mean
    1014             :         end do
    1015             : 
    1016             :         ! Only upper half of CovMat is needed
    1017           3 :         npMinusOneInverse = 1._RK / real(np-1,kind=RK)
    1018          12 :         do j = 1,nd
    1019          30 :             do i = 1,j
    1020             :                 ! Get the covariance matrix elements: only the upper half of CovMat is needed
    1021         117 :                 CholeskyLower(i,j) = dot_product( NormData(1:np,i) , NormData(1:np,j) ) * npMinusOneInverse
    1022             :             end do
    1023             :         end do
    1024             : 
    1025             :         ! XXX: The efficiency can be improved by merging it with the above loops
    1026           3 :         call getCholeskyFactor(nd,CholeskyLower,Diagonal)
    1027             : 
    1028           3 :   end subroutine getSamCholFac
    1029             : 
    1030             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1031             : 
    1032             :     !> \brief
    1033             :     !> Return the sample mean, covariance matrix, the Mahalanobis distances squared of the points with respect to the sample,
    1034             :     !> and the square root of the determinant of the inverse covariance matrix of the sample.
    1035             :     !>
    1036             :     !> \param[in]       nd                  :   The number of dimensions of the input sample.
    1037             :     !> \param[in]       np                  :   The number of points in the sample.
    1038             :     !> \param[in]       Point               :   The array of shape `(np,nd)` containing the sample.
    1039             :     !> \param[out]      CovMat              :   The output matrix of shape `(nd,nd)` representing the covariance matrix of the input data.
    1040             :     !> \param[out]      Mean                :   The output mean vector of the sample.
    1041             :     !> \param[out]      MahalSq             :   The output Mahalanobis distances squared of the points with respect to the sample (**optional**).
    1042             :     !> \param[out]      InvCovMat           :   The output inverse covariance matrix of the input data (**optional**).
    1043             :     !> \param[out]      sqrtDetInvCovMat    :   The output square root of the determinant of the inverse covariance matrix of the sample (**optional**).
    1044             :     !>
    1045             :     !> \warning
    1046             :     !> If `sqrtDetInvCovMat` is present, then `MahalSq` and `InvCovMat` must be also present.
    1047             :     !>
    1048             :     !> \remark
    1049             :     !> Note the shape of the input `Point(np,nd)`.
    1050             :     !>
    1051             :     !> \remark
    1052             :     !> See also, [getSamCovMeanTrans](@ref getsamcovmeantrans).
    1053             :     !>
    1054             :     !> \remark
    1055             :     !> For more information, see Geisser & Cornfield (1963) "Posterior distributions for multivariate normal parameters".
    1056             :     !> Also, Box and Tiao (1973), "Bayesian Inference in Statistical Analysis" Page 421.
    1057             :     !>
    1058             :     !> \author
    1059             :     !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
    1060           9 :     subroutine getSamCovMean(np,nd,Point,CovMat,Mean,MahalSq,InvCovMat,sqrtDetInvCovMat)
    1061             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1062             :         !DEC$ ATTRIBUTES DLLEXPORT :: getSamCovMean
    1063             : #endif
    1064           3 :         use Matrix_mod, only: getInvPosDefMatSqrtDet
    1065             :         implicit none
    1066             :         integer(IK), intent(in)             :: np,nd                  ! np is the number of observations, nd is the number of parameters for each observation
    1067             :         real(RK)   , intent(in)             :: Point(np,nd)           ! Point is the matrix of the data, CovMat contains the elements of the sample covariance matrix
    1068             :         real(RK)   , intent(out)            :: CovMat(nd,nd)          ! Covariance matrix of the input data
    1069             :         real(RK)   , intent(out)            :: Mean(nd)               ! Mean vector
    1070             :         real(RK)   , intent(out), optional  :: MahalSq(np)            ! Vector of Mahalanobis Distances Squared, with respect to the mean position of the sample
    1071             :         real(RK)   , intent(out), optional  :: InvCovMat(nd,nd)       ! Inverse Covariance matrix of the input data
    1072             :         real(RK)   , intent(out), optional  :: sqrtDetInvCovMat       ! sqrt determinant of the inverse covariance matrix
    1073          57 :         real(RK)                            :: NormData(np,nd)
    1074          15 :         real(RK)                            :: DummyVec(nd)
    1075             :         integer(IK)                         :: i,j
    1076             : 
    1077          12 :         do j = 1,nd
    1078          54 :             Mean(j) = sum(Point(1:np,j)) / real(np,kind=RK)
    1079          57 :             NormData(1:np,j) = Point(1:np,j) - Mean(j)
    1080             :         end do
    1081          12 :         do i = 1,nd
    1082          39 :             do j = 1,nd
    1083         171 :                 CovMat(i,j) = dot_product(NormData(1:np,i),NormData(1:np,j))/real(np-1,kind=RK)
    1084             :             end do
    1085             :         end do
    1086             : 
    1087           3 :         if (present(sqrtDetInvCovMat)) then   ! Calculate InCovMat, MahalSq, and sqrt Determinant of InCovMat
    1088          12 :             do j = 1,nd
    1089          30 :                 do i = 1,j
    1090          27 :                 InvCovMat(i,j) = CovMat(i,j)    ! Only the upper half of CovMat is needed
    1091             :                 end do
    1092             :             end do
    1093           3 :             call getInvPosDefMatSqrtDet(nd,InvCovMat,sqrtDetInvCovMat)
    1094          18 :             do i = 1,np
    1095          60 :                 do j = 1,nd
    1096         195 :                 DummyVec(j) = dot_product(InvCovMat(1:nd,j),NormData(i,1:nd))
    1097             :                 end do
    1098          63 :                 MahalSq(i) = dot_product(NormData(i,1:nd),DummyVec)
    1099             :             end do
    1100             :         end if
    1101             : 
    1102           3 :     end subroutine getSamCovMean
    1103             : 
    1104             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1105             : 
    1106             :     !> \brief
    1107             :     !> Return the sample mean, covariance matrix, the Mahalanobis distances squared of the points with respect to the sample,
    1108             :     !> and the square root of the determinant of the inverse covariance matrix of the sample.
    1109             :     !>
    1110             :     !> \param[in]       nd                  :   The number of dimensions of the input sample.
    1111             :     !> \param[in]       np                  :   The number of points in the sample.
    1112             :     !> \param[in]       Point               :   The array of shape `(nd,np)` containing the sample.
    1113             :     !> \param[out]      CovMat              :   The output matrix of shape `(nd,nd)` representing the covariance matrix of the input data.
    1114             :     !> \param[out]      Mean                :   The output mean vector of the sample.
    1115             :     !> \param[out]      MahalSq             :   The output Mahalanobis distances squared of the points with respect to the sample (**optional**).
    1116             :     !> \param[out]      InvCovMat           :   The output inverse covariance matrix of the input data (**optional**).
    1117             :     !> \param[out]      sqrtDetInvCovMat    :   The output square root of the determinant of the inverse covariance matrix of the sample (**optional**).
    1118             :     !>
    1119             :     !> \warning
    1120             :     !> If `sqrtDetInvCovMat` is present, then `MahalSq` and `InvCovMat` must be also present.
    1121             :     !>
    1122             :     !> \remark
    1123             :     !> Note the shape of the input `Point(nd,np)`.
    1124             :     !>
    1125             :     !> \remark
    1126             :     !> This subroutine has the same functionality as [getSamCovMean](@ref getsamcovmean), with the only difference that input data is transposed here on input.
    1127             :     !> Based on the preliminary benchmarks with Intel 2017 ifort, `getSamCovMean()` is slightly faster than `getSamCovMeanTrans()`.
    1128             :     !>
    1129             :     !> \remark
    1130             :     !> For more information, see Geisser & Cornfield (1963) "Posterior distributions for multivariate normal parameters".
    1131             :     !> Also, Box and Tiao (1973), "Bayesian Inference in Statistical Analysis" Page 421.
    1132             :     !>
    1133             :     !> \author
    1134             :     !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
    1135          18 :     subroutine getSamCovMeanTrans(np,nd,Point,CovMat,Mean,MahalSq,InvCovMat,sqrtDetInvCovMat)
    1136             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1137             :         !DEC$ ATTRIBUTES DLLEXPORT :: getSamCovMeanTrans
    1138             : #endif
    1139             : 
    1140           3 :         use Matrix_mod, only: getInvPosDefMatSqrtDet
    1141             :         implicit none
    1142             :         integer(IK), intent(in)            :: np,nd                  ! np is the number of observations, nd is the number of parameters for each observation
    1143             :         real(RK)   , intent(in)            :: Point(nd,np)           ! Point is the matrix of the data, CovMat contains the elements of the sample covariance matrix
    1144             :         real(RK)   , intent(out)           :: CovMat(nd,nd)          ! Covariance matrix of the input data
    1145             :         real(RK)   , intent(out)           :: Mean(nd)               ! Mean vector
    1146             :         real(RK)   , intent(out), optional :: MahalSq(np)            ! Vector of Mahalanobis Distances Squared, with respect to the mean position of the sample
    1147             :         real(RK)   , intent(out), optional :: InvCovMat(nd,nd)       ! Inverse Covariance matrix of the input data
    1148             :         real(RK)   , intent(out), optional :: sqrtDetInvCovMat       ! sqrt determinant of the inverse covariance matrix
    1149          60 :         real(RK)   , dimension(nd)         :: DummyVec
    1150         432 :         real(RK)   , dimension(nd,np)      :: NormData
    1151             :         integer(IK)                        :: i,j
    1152             : 
    1153          48 :         Mean = 0._RK
    1154         117 :         do i = 1,np
    1155         432 :             do j = 1,nd
    1156         420 :                 Mean(j) = Mean(j) + Point(j,i)
    1157             :             end do
    1158             :         end do
    1159          48 :         Mean = Mean / real(np,kind=RK)
    1160             : 
    1161         117 :         do i = 1,np
    1162         432 :             NormData(1:nd,i) = Point(1:nd,i) - Mean
    1163             :         end do
    1164             : 
    1165          48 :         do i = 1,nd
    1166         156 :             do j = 1,nd
    1167        1089 :                 CovMat(i,j) = dot_product(NormData(i,1:np),NormData(j,1:np)) / real(np-1,kind=RK)
    1168             :             end do
    1169             :         end do
    1170             : 
    1171          12 :         if (present(sqrtDetInvCovMat)) then   ! Calculate InCovMat, MahalSq, and sqrt Determinant of InCovMat
    1172          12 :             do j = 1,nd
    1173          30 :                 do i = 1,j
    1174          27 :                 InvCovMat(i,j) = CovMat(i,j)    ! Only the upper half of CovMat is needed
    1175             :                 end do
    1176             :             end do
    1177           3 :             call getInvPosDefMatSqrtDet(nd,InvCovMat,sqrtDetInvCovMat)
    1178          18 :             do i = 1,np
    1179          60 :                 do j = 1,nd
    1180         195 :                 DummyVec(j) = dot_product(InvCovMat(1:nd,j),NormData(1:nd,i))
    1181             :                 end do
    1182          63 :                 MahalSq(i) = dot_product(NormData(1:nd,i),DummyVec)
    1183             :                 !MahalSq = dot_product(NormData(1:nd,i),DummyVec)
    1184             :                 !if (maxMahal<MahalSq) maxMahal = MahalSq
    1185             :             end do
    1186             :             !maxMahal = maxval(MahalSq)
    1187             :         end if
    1188             : 
    1189          12 :     end subroutine getSamCovMeanTrans
    1190             : 
    1191             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1192             : 
    1193             :     !> \brief
    1194             :     !> Return the sample mean and the upper triangle of the covariance matrix of the input sample.
    1195             :     !>
    1196             :     !> \param[in]       nd                  :   The number of dimensions of the input sample.
    1197             :     !> \param[in]       np                  :   The number of points in the sample.
    1198             :     !> \param[in]       Point               :   The array of shape `(nd,np)` containing the sample.
    1199             :     !> \param[out]      CovMatUpper         :   The output matrix of shape `(nd,nd)` whose upper triangle represents the covariance matrix of the input data.
    1200             :     !> \param[out]      Mean                :   The output mean vector of the sample.
    1201             :     !>
    1202             :     !> \remark
    1203             :     !> Note the shape of the input `Point(nd,np)`.
    1204             :     !>
    1205             :     !> \remark
    1206             :     !> This subroutine has the same functionality as [getSamCovMeanTrans](@ref getsamcovmeantrans), with the only difference
    1207             :     !> only the upper triangle of the covariance matrix is returned. Also, optional arguments are not available.
    1208             :     !> This subroutine is specifically optimized for use in the ParaMonte samplers.
    1209             :     !>
    1210             :     !> \remark
    1211             :     !> For more information, see Geisser & Cornfield (1963) "Posterior distributions for multivariate normal parameters".
    1212             :     !> Also, Box and Tiao (1973), "Bayesian Inference in Statistical Analysis" Page 421.
    1213             :     !>
    1214             :     !> \author
    1215             :     !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
    1216          64 :     subroutine getSamCovUpperMeanTrans(np,nd,Point,CovMatUpper,Mean)
    1217             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1218             :         !DEC$ ATTRIBUTES DLLEXPORT :: getSamCovUpperMeanTrans
    1219             : #endif
    1220          12 :         use Matrix_mod, only: getInvPosDefMatSqrtDet
    1221             :         implicit none
    1222             :         integer(IK), intent(in)            :: np,nd                  ! np is the number of observations, nd is the number of parameters for each observation
    1223             :         real(RK)   , intent(in)            :: Point(nd,np)           ! Point is the matrix of the data, CovMatUpper contains the elements of the sample covariance matrix
    1224             :         real(RK)   , intent(out)           :: CovMatUpper(nd,nd)     ! Covariance matrix of the input data
    1225             :         real(RK)   , intent(out)           :: Mean(nd)               ! Mean vector
    1226        1180 :         real(RK)                           :: npMinusOneInvReal, NormData(nd,np)
    1227             :         integer(IK)                        :: i,j
    1228             : 
    1229         194 :         Mean = 0._RK
    1230         391 :         do i = 1,np
    1231        1180 :             do j = 1,nd
    1232        1116 :                 Mean(j) = Mean(j) + Point(j,i)
    1233             :             end do
    1234             :         end do
    1235         194 :         Mean = Mean / real(np,kind=RK)
    1236             : 
    1237         391 :         do i = 1,np
    1238        1180 :             NormData(1:nd,i) = Point(1:nd,i) - Mean
    1239             :         end do
    1240             : 
    1241          64 :         npMinusOneInvReal = 1._RK / real(np-1,kind=RK)
    1242         194 :         do j = 1,nd
    1243         411 :             do i = 1,j
    1244        1793 :                 CovMatUpper(i,j) = dot_product(NormData(i,1:np),NormData(j,1:np)) * npMinusOneInvReal
    1245             :             end do
    1246             :         end do
    1247             : 
    1248          64 :   end subroutine getSamCovUpperMeanTrans
    1249             : 
    1250             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1251             : 
    1252             :     !> \brief
    1253             :     !> Return the mean and the upper triangle of the covariance matrix of the input *weighted* sample.
    1254             :     !>
    1255             :     !> \param[in]       nd                  :   The number of dimensions of the input sample.
    1256             :     !> \param[in]       sumWeight           :   The sum of all sample weights.
    1257             :     !> \param[in]       np                  :   The number of points in the sample.
    1258             :     !> \param[in]       Point               :   The array of shape `(nd,np)` containing the sample.
    1259             :     !> \param[in]       Weight              :   The integer array of length `np`, representing the weights of individual points in the sample.
    1260             :     !> \param[out]      CovMatUpper         :   The output matrix of shape `(nd,nd)` whose upper triangle represents the covariance matrix of the input data.
    1261             :     !> \param[out]      Mean                :   The output mean vector of the sample.
    1262             :     !>
    1263             :     !> \warning
    1264             :     !> Note the shape of the input argument `Point(nd,np)`.
    1265             :     !>
    1266             :     !> \remark
    1267             :     !> This subroutine has the same functionality as [getSamCovUpperMeanTrans](@ref getsamcovuppermeantrans), with the only difference
    1268             :     !> only the upper triangle of the covariance matrix is returned. Also, optional arguments are not available.
    1269             :     !> This subroutine is specifically optimized for use in the ParaMonte samplers.
    1270             :     !>
    1271             :     !> \remark
    1272             :     !> For more information, see Geisser & Cornfield (1963) "Posterior distributions for multivariate normal parameters".
    1273             :     !> Also, Box and Tiao (1973), "Bayesian Inference in Statistical Analysis" Page 421.
    1274             :     !>
    1275             :     !> \author
    1276             :     !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
    1277       30906 :     subroutine getWeiSamCovUppMeanTrans(np,sumWeight,nd,Point,Weight,CovMatUpper,Mean)
    1278             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1279             :         !DEC$ ATTRIBUTES DLLEXPORT :: getWeiSamCovUppMeanTrans
    1280             : #endif
    1281             : 
    1282          64 :         use Matrix_mod, only: getInvPosDefMatSqrtDet
    1283             :         implicit none
    1284             :         integer(IK), intent(in)             :: np,nd,sumWeight        ! np is the number of observations, nd is the number of parameters for each observation
    1285             :         integer(IK), intent(in)             :: Weight(np)             ! weights of the points
    1286             :         real(RK)   , intent(in)             :: Point(nd,np)           ! Point is the matrix of the data, CovMatUpper contains the elements of the sample covariance matrix
    1287             :         real(RK)   , intent(out)            :: CovMatUpper(nd,nd)     ! Covariance matrix of the input data
    1288             :         real(RK)   , intent(out)            :: Mean(nd)               ! Mean vector
    1289       30906 :         real(RK)                            :: sumWeightMinusOneInvReal
    1290      696554 :         real(RK)                            :: NormData(nd,np)
    1291             :         integer(IK)                         :: i,j,ip
    1292             : 
    1293       80924 :         Mean = 0._RK
    1294      278138 :         do i = 1,np
    1295      696554 :             do j = 1,nd
    1296      665648 :                 Mean(j) = Mean(j) + Weight(i)*Point(j,i)
    1297             :             end do
    1298             :         end do
    1299       80924 :         Mean = Mean / real(sumWeight,kind=RK)
    1300             : 
    1301      278138 :         do i = 1,np
    1302      696554 :             NormData(1:nd,i) = Point(1:nd,i) - Mean
    1303             :         end do
    1304             : 
    1305       30906 :         sumWeightMinusOneInvReal = 1._RK / real(sumWeight-1,kind=RK)
    1306       80924 :         do j = 1,nd
    1307      150057 :             do i = 1,j
    1308       69133 :                 CovMatUpper(i,j) = 0
    1309      658748 :                 do ip = 1,np
    1310      658748 :                     CovMatUpper(i,j) = CovMatUpper(i,j) + Weight(ip)*NormData(i,ip)*NormData(j,ip)
    1311             :                 end do
    1312      119151 :                 CovMatUpper(i,j) = CovMatUpper(i,j) * sumWeightMinusOneInvReal
    1313             :             end do
    1314             :         end do
    1315             : 
    1316       30906 :     end subroutine getWeiSamCovUppMeanTrans
    1317             : 
    1318             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1319             : 
    1320             :     !> \brief
    1321             :     !> Given two input sample means and covariance matrices, return the combination of them as a single mean and covariance matrix.
    1322             :     !>
    1323             :     !> \param[in]       nd          :   The number of dimensions of the input sample.
    1324             :     !> \param[in]       npA         :   The number of points in sample `A`.
    1325             :     !> \param[in]       MeanVecA    :   The mean vector of sample `A`.
    1326             :     !> \param[in]       CovMatA     :   The covariance matrix of sample `A`.
    1327             :     !> \param[in]       npB         :   The number of points in sample `B`.
    1328             :     !> \param[in]       MeanVecB    :   The mean vector of sample `B`.
    1329             :     !> \param[in]       CovMatB     :   The covariance matrix of sample `B`.
    1330             :     !> \param[out]      MeanVecAB   :   The output mean vector of the combined sample.
    1331             :     !> \param[out]      CovMatAB    :   The output covariance matrix of the combined sample.
    1332             :     !>
    1333             :     !> \author
    1334             :     !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
    1335             :     !>
    1336             :     !> \remark
    1337             :     !> An exact implementation of this algorithm which needs only the upper triangles of the input matrices and
    1338             :     !> yields only the upper triangle of the covariance matrix is given in [mergeMeanCovUpper](@ref mergemeancovupper).
    1339             :     !> The alternative implementation is much more efficient, by a factor of 6-7 with all compiler optimization flags on.
    1340             :     !>
    1341             :     ! This subroutine uses a recursion equation similar to http://stats.stackexchange.com/questions/97642/how-to-combine-sample-means-and-sample-variances
    1342             :     ! See also: O’Neill, (2014), "Some Useful Moment Results in Sampling Problems".
    1343             :     ! See also: https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Covariance
    1344             :     ! See also: https://stats.stackexchange.com/questions/43159/how-to-calculate-pooled-variance-of-two-groups-given-known-group-variances-mean
    1345           3 :     subroutine mergeMeanCov(nd,npA,MeanVecA,CovMatA,npB,MeanVecB,CovMatB,MeanVecAB,CovMatAB)
    1346             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1347             :         !DEC$ ATTRIBUTES DLLEXPORT :: mergeMeanCov
    1348             : #endif
    1349             :         implicit none
    1350             :         integer(IK), intent(in)       :: nd
    1351             :         integer(IK), intent(in)       :: npA,npB
    1352             :         real(RK)   , intent(in)       :: MeanVecA(nd),CovMatA(nd,nd)
    1353             :         real(RK)   , intent(in)       :: MeanVecB(nd),CovMatB(nd,nd)
    1354             :         real(RK)   , intent(out)      :: MeanVecAB(nd),CovMatAB(nd,nd)
    1355          42 :         real(RK)   , dimension(nd,1)  :: MeanMatA,MeanMatB,MeanMat
    1356          42 :         real(RK)                      :: DistanceSq(nd,nd)
    1357           3 :         real(RK)                      :: npABinverse
    1358             :         integer(IK)                   :: npAB
    1359             : 
    1360           3 :         npAB = npA + npB
    1361           3 :         npABinverse = 1._RK / real(npAB, kind=RK)
    1362          12 :         MeanMatA(1:nd,1) = MeanVecA
    1363          12 :         MeanMatB(1:nd,1) = MeanVecB
    1364             : 
    1365             :         ! Compute the new Mean
    1366             : 
    1367          12 :         MeanVecAB = ( npA * MeanVecA + npB * MeanVecB ) * npABinverse
    1368          12 :         MeanMat(1:nd,1) = MeanVecAB
    1369             : 
    1370             :         ! Compute the new Covariance matrix
    1371             : 
    1372          69 :         DistanceSq = matmul( (MeanMatA-MeanMatB), transpose((MeanMatA-MeanMatB)) ) * npA * npB * npABinverse
    1373          39 :         CovMatAB = ( (npA-1) * CovMatA + (npB-1) * CovMatB + DistanceSq ) / real(npAB-1, kind=RK)
    1374             : 
    1375       30906 :     end subroutine mergeMeanCov
    1376             : 
    1377             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1378             : 
    1379             : !    !> \brief
    1380             : !    !> Given two input sample means and covariance matrices, return the combination of them as a single mean and covariance matrix.
    1381             : !    !>
    1382             : !    !> \param[in]       nd              :   The number of dimensions of the input sample.
    1383             : !    !> \param[in]       npA             :   The number of points in sample `A`.
    1384             : !    !> \param[in]       MeanVecA        :   The mean vector of sample `A`.
    1385             : !    !> \param[in]       CovMatUpperA    :   The covariance matrix of sample `A`.
    1386             : !    !> \param[in]       npB             :   The number of points in sample `B`.
    1387             : !    !> \param[in]       MeanVecB        :   The mean vector of sample `B`.
    1388             : !    !> \param[in]       CovMatUpperB    :   The covariance matrix of sample `B`.
    1389             : !    !> \param[out]      MeanVecAB       :   The output mean vector of the combined sample.
    1390             : !    !> \param[out]      CovMatUpperAB   :   The output covariance matrix of the combined sample.
    1391             : !    !>
    1392             : !    !> \todo
    1393             : !    !> The efficiency of this algorithm might still be improved by converting the upper triangle covariance matrix to a packed vector.
    1394             : !    !>
    1395             : !    !> \remark
    1396             : !    !> This subroutine is the same as [mergeMeanCov](@ref mergeMeanCov), with the **important difference** that only the
    1397             : !    !> upper triangles and diagonals of the input covariance matrices need to be given by the user: `CovMatUpperA`, `CovMatUpperB`
    1398             : !    !> This alternative implementation is 6-7 times faster, with all compiler optimization flags on.
    1399             : !    !>
    1400             : !    !> \author
    1401             : !    !> Amir Shahmoradi, Nov 24, 2020, 4:19 AM, Dallas, TX
    1402             : !    subroutine mergeMeanCovUpperSlow(nd,npA,MeanVecA,CovMatUpperA,npB,MeanVecB,CovMatUpperB,MeanVecAB,CovMatUpperAB)
    1403             : !#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1404             : !        !DEC$ ATTRIBUTES DLLEXPORT :: mergeMeanCovUpperSlow
    1405             : !#endif
    1406             : !        implicit none
    1407             : !        integer(IK), intent(in)             :: nd
    1408             : !        integer(IK), intent(in)             :: npA,npB
    1409             : !        real(RK)   , intent(in)             :: MeanVecA(nd),CovMatUpperA(nd,nd)
    1410             : !        real(RK)   , intent(in)             :: MeanVecB(nd),CovMatUpperB(nd,nd)
    1411             : !        real(RK)   , intent(out)            :: MeanVecAB(nd),CovMatUpperAB(nd,nd)
    1412             : !        real(RK)                            :: npABinverse, npAnpB2npAB
    1413             : !        real(RK)                            :: npA2npAB, npB2npAB
    1414             : !        real(RK)                            :: MeanVecDiffAB(nd)
    1415             : !        integer(IK)                         :: npAB, i, j
    1416             : !
    1417             : !        npAB = npA + npB
    1418             : !        npABinverse = 1._RK / real(npAB, kind=RK)
    1419             : !        npAnpB2npAB = npA * npB * npABinverse
    1420             : !        npA2npAB = npA * npABinverse
    1421             : !        npB2npAB = npB * npABinverse
    1422             : !
    1423             : !        ! Compute the new Mean and Covariance matrix
    1424             : !
    1425             : !        do j = 1, nd
    1426             : !            MeanVecDiffAB(j) = MeanVecA(j) - MeanVecB(j)
    1427             : !            !MeanVecAB(j) = ( npA * MeanVecA(j) + npB * MeanVecB(j) ) * npABinverse
    1428             : !            MeanVecAB(j) = npA2npAB * MeanVecA(j) + npB2npAB * MeanVecB(j)
    1429             : !            do i = 1, j
    1430             : !                CovMatUpperAB(i,j) = ( (npA-1) * CovMatUpperA(i,j) + (npB-1) * CovMatUpperB(i,j) + MeanVecDiffAB(i) * MeanVecDiffAB(j) * npAnpB2npAB ) / real(npAB-1, kind=RK)
    1431             : !            end do
    1432             : !        end do
    1433             : !
    1434             : !
    1435             : !    end subroutine mergeMeanCovUpperSlow
    1436             : 
    1437             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1438             : 
    1439             :     !> \brief
    1440             :     !> Given two input sample means and covariance matrices, return the combination of them as a single mean and covariance matrix.
    1441             :     !>
    1442             :     !> \param[in]       nd              :   The number of dimensions of the input sample.
    1443             :     !> \param[in]       npA             :   The number of points in sample `A`.
    1444             :     !> \param[in]       MeanVecA        :   The mean vector of sample `A`.
    1445             :     !> \param[in]       CovMatUpperA    :   The covariance matrix of sample `A`.
    1446             :     !> \param[in]       npB             :   The number of points in sample `B`.
    1447             :     !> \param[in]       MeanVecB        :   The mean vector of sample `B`.
    1448             :     !> \param[in]       CovMatUpperB    :   The covariance matrix of sample `B`.
    1449             :     !> \param[out]      MeanVecAB       :   The output mean vector of the combined sample.
    1450             :     !> \param[out]      CovMatUpperAB   :   The output covariance matrix of the combined sample.
    1451             :     !>
    1452             :     !> \remark
    1453             :     !> This subroutine is the same as [mergeMeanCov](@ref mergemeancov), with the **important difference** that only the
    1454             :     !> upper triangles and diagonals of the input covariance matrices need to be given by the user: `CovMatUpperA`, `CovMatUpperB`
    1455             :     !> This alternative implementation is 6-7 times faster, with all compiler optimization flags on.
    1456             :     !> In addition, all computational coefficients are predefined in this implementation,
    1457             :     !> resulting in an extra 10%-15% efficiency gain.
    1458             :     !>
    1459             :     !> \todo
    1460             :     !> The efficiency of this algorithm might still be improved by converting the upper triangle covariance matrix to a packed vector.
    1461             :     !>
    1462             :     !> \author
    1463             :     !> Amir Shahmoradi, Nov 24, 2020, 4:19 AM, Dallas, TX
    1464       30175 :     subroutine mergeMeanCovUpper(nd,npA,MeanVecA,CovMatUpperA,npB,MeanVecB,CovMatUpperB,MeanVecAB,CovMatUpperAB)
    1465             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1466             :         !DEC$ ATTRIBUTES DLLEXPORT :: mergeMeanCovUpper
    1467             : #endif
    1468             :         implicit none
    1469             :         integer(IK), intent(in)             :: nd
    1470             :         integer(IK), intent(in)             :: npA,npB
    1471             :         real(RK)   , intent(in)             :: MeanVecA(nd),CovMatUpperA(nd,nd)
    1472             :         real(RK)   , intent(in)             :: MeanVecB(nd),CovMatUpperB(nd,nd)
    1473             :         real(RK)   , intent(out)            :: MeanVecAB(nd),CovMatUpperAB(nd,nd)
    1474       79028 :         real(RK)                            :: MeanVecDiffAB(nd)
    1475       30175 :         real(RK)                            :: npAnpB2npAB2npABMinusOne
    1476       30175 :         real(RK)                            :: npAMinusOne2npABMinusOne
    1477       30175 :         real(RK)                            :: npBMinusOne2npABMinusOne
    1478       30175 :         real(RK)                            :: npABMinusOneInverse
    1479       30175 :         real(RK)                            :: npABinverse
    1480       30175 :         real(RK)                            :: npA2npAB
    1481       30175 :         real(RK)                            :: npB2npAB
    1482             :         integer(IK)                         :: npAB, i, j
    1483             : 
    1484       30175 :         npAB = npA + npB
    1485       30175 :         npABinverse = 1._RK / real(npAB, kind=RK)
    1486       30175 :         npABMinusOneInverse = 1._RK / real(npAB-1, kind=RK)
    1487       30175 :         npAnpB2npAB2npABMinusOne = npA * npB * npABinverse * npABMinusOneInverse
    1488       30175 :         npA2npAB = npA * npABinverse
    1489       30175 :         npB2npAB = npB * npABinverse
    1490       30175 :         npAMinusOne2npABMinusOne = (npA - 1_IK) * npABMinusOneInverse
    1491       30175 :         npBMinusOne2npABMinusOne = (npB - 1_IK) * npABMinusOneInverse
    1492             : 
    1493             :         ! Compute the new Mean and Covariance matrix
    1494             : 
    1495       79028 :         do j = 1, nd
    1496       48853 :             MeanVecDiffAB(j) = MeanVecA(j) - MeanVecB(j)
    1497       48853 :             MeanVecAB(j) = npA2npAB * MeanVecA(j) + npB2npAB * MeanVecB(j)
    1498      146562 :             do i = 1, j
    1499       67534 :                 CovMatUpperAB(i,j)  = npAMinusOne2npABMinusOne * CovMatUpperA(i,j) &
    1500       67534 :                                     + npBMinusOne2npABMinusOne * CovMatUpperB(i,j) &
    1501      116387 :                                     + npAnpB2npAB2npABMinusOne * MeanVecDiffAB(i) * MeanVecDiffAB(j)
    1502             :             end do
    1503             :         end do
    1504             : 
    1505             : 
    1506           3 :     end subroutine mergeMeanCovUpper
    1507             : 
    1508             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1509             : 
    1510             :     !> \brief
    1511             :     !> Given two input sample means and covariance matrices, return the combination of them as a single mean and covariance matrix.
    1512             :     !>
    1513             :     !> \param[in]       nd              :   The number of dimensions of the input sample.
    1514             :     !> \param[in]       npA             :   The number of points in sample `A`.
    1515             :     !> \param[in]       MeanVecA        :   The mean vector of sample `A`.
    1516             :     !> \param[in]       CovMatUpperA    :   The covariance matrix of sample `A`.
    1517             :     !> \param[in]       npB             :   The number of points in sample `B`.
    1518             :     !> \param[inout]    MeanVecB        :   The mean vector of sample `B`.
    1519             :     !> \param[inout]    CovMatUpperB    :   The covariance matrix of sample `B`.
    1520             :     !>
    1521             :     !> \remark
    1522             :     !> This subroutine is the same as [mergeMeanCovUpper](@ref mergemeancovupper), with the **important difference** that
    1523             :     !> the resulting output mean and covariance matrices are written to the input arguments `MeanVecB`, `CovMatUpperB`.
    1524             :     !> This alternative implementation results in another extra 15%-20% efficiency gain. This result is based on
    1525             :     !> the benchmarks with Intel Fortran compiler 19.4 with all compiler optimization flags on.
    1526             :     !>
    1527             :     !> \todo
    1528             :     !> The efficiency of this algorithm might still be improved by converting the upper triangle covariance matrix to a packed vector.
    1529             :     !>
    1530             :     !> \author
    1531             :     !> Amir Shahmoradi, Nov 25, 2020, 1:00 AM, Dallas, TX
    1532           6 :     subroutine mergeMeanCovUpperDense(nd,npA,MeanVecA,CovMatUpperA,npB,MeanVecB,CovMatUpperB)
    1533             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1534             :         !DEC$ ATTRIBUTES DLLEXPORT :: mergeMeanCovUpperDense
    1535             : #endif
    1536             :         implicit none
    1537             :         integer(IK), intent(in)             :: nd
    1538             :         integer(IK), intent(in)             :: npA,npB
    1539             :         real(RK)   , intent(in)             :: MeanVecA(nd),CovMatUpperA(nd,nd)
    1540             :         real(RK)   , intent(inout)          :: MeanVecB(nd),CovMatUpperB(nd,nd)
    1541          24 :         real(RK)                            :: MeanVecDiffAB(nd)
    1542           6 :         real(RK)                            :: npAnpB2npAB2npABMinusOne
    1543           6 :         real(RK)                            :: npAMinusOne2npABMinusOne
    1544           6 :         real(RK)                            :: npBMinusOne2npABMinusOne
    1545           6 :         real(RK)                            :: npABMinusOneInverse
    1546           6 :         real(RK)                            :: npABinverse
    1547           6 :         real(RK)                            :: npA2npAB
    1548           6 :         real(RK)                            :: npB2npAB
    1549             :         integer(IK)                         :: npAB, i, j
    1550             : 
    1551           6 :         npAB = npA + npB
    1552           6 :         npABinverse = 1._RK / real(npAB, kind=RK)
    1553           6 :         npABMinusOneInverse = 1._RK / real(npAB-1, kind=RK)
    1554           6 :         npAnpB2npAB2npABMinusOne = npA * npB * npABinverse * npABMinusOneInverse
    1555           6 :         npA2npAB = npA * npABinverse
    1556           6 :         npB2npAB = npB * npABinverse
    1557           6 :         npAMinusOne2npABMinusOne = (npA - 1_IK) * npABMinusOneInverse
    1558           6 :         npBMinusOne2npABMinusOne = (npB - 1_IK) * npABMinusOneInverse
    1559             : 
    1560             :         ! Compute the new Mean and Covariance matrix
    1561             : 
    1562          24 :         do j = 1, nd
    1563          18 :             MeanVecDiffAB(j) = MeanVecA(j) - MeanVecB(j)
    1564          18 :             MeanVecB(j) = npA2npAB * MeanVecA(j) + npB2npAB * MeanVecB(j)
    1565          60 :             do i = 1, j
    1566          36 :                 CovMatUpperB(i,j)   = npAMinusOne2npABMinusOne * CovMatUpperA(i,j) &
    1567          36 :                                     + npBMinusOne2npABMinusOne * CovMatUpperB(i,j) &
    1568          54 :                                     + npAnpB2npAB2npABMinusOne * MeanVecDiffAB(i) * MeanVecDiffAB(j)
    1569             :             end do
    1570             :         end do
    1571             : 
    1572       30175 :     end subroutine mergeMeanCovUpperDense
    1573             : 
    1574             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1575             : 
    1576             :     !> \brief
    1577             :     !> Return a random standard Gaussian deviate with zero mean and unit variance.
    1578             :     !>
    1579             :     !> \return
    1580             :     !> `randGaus` : The random standard Gaussian deviate with zero mean and unit variance.
    1581             :     !>
    1582             :     !> \remark
    1583             :     !> See also, Numerical Recipes in Fortran, by Press et al. (1990)
    1584             :     !>
    1585             :     !> \author
    1586             :     !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
    1587      891872 :     function getRandGaus() result(randGaus)
    1588             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1589             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandGaus
    1590             : #endif
    1591             : 
    1592             :         implicit none
    1593             :         integer(IK), save :: iset=0
    1594             :         real(RK)   , save :: gset
    1595      891872 :         real(RK)          :: fac,rsq,vec(2)
    1596             :         real(RK)          :: randGaus
    1597             : 
    1598      891872 :         if (iset == 0_IK) then
    1599      122007 :             do
    1600             :                 !block
    1601             :                 !integer :: i, n
    1602             :                 !real(RK) :: unifrnd(30)
    1603             :                 !!integer, dimension(:), allocatable :: seed
    1604             :                 !if (paradramPrintEnabled .or. paradisePrintEnabled) then
    1605             :                 !    !do i = 1, 22
    1606             :                 !    call random_number(unifrnd)
    1607             :                 !    write(*,"(*(g0,:,'"//new_line("a")//"'))") unifrnd
    1608             :                 !    !end do
    1609             :                 !    !call random_seed(size = n); allocate(seed(n))
    1610             :                 !    !call random_seed(get = seed)
    1611             :                 !    !write(*,"(*(g0,:,' '))") seed
    1612             :                 !    !write(*,"(*(g0,:,' '))") StateOld
    1613             :                 !    !write(*,"(*(g0,:,' '))") StateNew
    1614             :                 !    !write(*,"(*(g0,:,' '))") CholeskyLower
    1615             :                 !    !write(*,"(*(g0,:,' '))") domainCheckCounter
    1616             :                 !    paradisePrintEnabled = .false.
    1617             :                 !    paradramPrintEnabled = .false.
    1618             :                 !end if
    1619             :                 !end block
    1620      567944 :                 call random_number(vec)
    1621     1703830 :                 vec = 2._RK*vec - 1._RK
    1622      567944 :                 rsq = vec(1)**2 + vec(2)**2
    1623      567944 :                 if ( rsq > 0._RK .and. rsq < 1._RK ) exit
    1624             :             end do
    1625      445937 :             fac = sqrt(-2._RK*log(rsq)/rsq)
    1626      445937 :             gset = vec(1)*fac
    1627      445937 :             randGaus = vec(2)*fac
    1628      445937 :             iset = 1_IK
    1629             :         else
    1630      445935 :             randGaus = gset
    1631      445935 :             iset = 0_IK
    1632             :         endif
    1633             : 
    1634      891878 :     end function getRandGaus
    1635             : 
    1636             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1637             : 
    1638             :     !> \brief
    1639             :     !> Return a random Gaussian deviate with the given mean and standard deviation.
    1640             :     !>
    1641             :     !> \param[in]       mean    :   The mean of the Gaussian distribution.
    1642             :     !> \param[in]       std     :   The standard deviation of the Gaussian distribution. It must be a positive real number.
    1643             :     !>
    1644             :     !> \return
    1645             :     !> `randNorm` : A normally distributed deviate with the given mean and standard deviation.
    1646             :     !>
    1647             :     !> \author
    1648             :     !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
    1649           3 :     function getRandNorm(mean,std) result(randNorm)
    1650             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1651             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandNorm
    1652             : #endif
    1653             :         implicit none
    1654             :         real(RK), intent(in)    :: mean, std
    1655             :         real(RK)                :: randNorm
    1656           6 :         randNorm = mean + std * getRandGaus()
    1657      891875 :     end function getRandNorm
    1658             : 
    1659             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1660             : 
    1661             :     !> \brief
    1662             :     !> Return a log-normally distributed deviate with the given mean and standard deviation.
    1663             :     !>
    1664             :     !> \param[in]       mean    :   The mean of the Lognormal distribution.
    1665             :     !> \param[in]       std     :   The standard deviation of the Lognormal distribution. It must be a positive real number.
    1666             :     !>
    1667             :     !> \return
    1668             :     !> `randLogn` : A Lognormally distributed deviate with the given mean and standard deviation.
    1669             :     !>
    1670             :     !> \author
    1671             :     !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
    1672           3 :     function getRandLogn(mean,std) result(randLogn)
    1673             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1674             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandLogn
    1675             : #endif
    1676             :         implicit none
    1677             :         real(RK), intent(in)    :: mean, std
    1678             :         real(RK)                :: randLogn
    1679           6 :         randLogn = exp( mean + std*getRandGaus() )
    1680           6 :     end function getRandLogn
    1681             : 
    1682             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1683             : 
    1684             :     ! This subroutine is legacy and slow. use getRandMVN() in this same module.
    1685             :     ! Given the mean vector MeanVec and the covariance matrix CovMat, this subroutine generates a random vector x (of length nd>=2)
    1686             :     ! from an nd-dimensional multivariate normal distribution.
    1687             :     ! First a vector of nd standard normal random deviates is generated,
    1688             :     ! then this vector is multiplied by the Cholesky decomposition of the covariance matrix,
    1689             :     ! then the vector MeanVec is added to the resulting vector, and is stored in the output vector x.
    1690             :     ! ATTENTION: Only the upper half of the covariance matrix (plus the diagonal terms) need to be given in the input.
    1691             :     ! in the ouput, the upper half and diagonal part will still be the covariance matrix, while the lower half will be
    1692             :     ! the Cholesky decomposition elements (excluding its diagonal terms that are provided only in the vector Diagonal).
    1693             :     ! USES choldc.f90, getRandGaus.f90
    1694             :     !> Amir Shahmoradi, March 22, 2012, 2:21 PM, IFS, UTEXAS
    1695           3 :     subroutine getMVNDev(nd,MeanVec,CovMat,X)
    1696             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1697             :         !DEC$ ATTRIBUTES DLLEXPORT :: getMVNDev
    1698             : #endif
    1699             : 
    1700             :         use iso_fortran_env, only: output_unit
    1701           3 :         use Matrix_mod, only: getCholeskyFactor
    1702             : 
    1703             :         implicit none
    1704             :         integer(IK), intent(in)  :: nd
    1705             :         real(RK)   , intent(in)  :: MeanVec(nd), CovMat(nd,nd)
    1706             :         real(RK)   , intent(out) :: X(nd)
    1707          33 :         real(RK)                 :: CholeskyLower(nd,nd), Diagonal(nd), DummyVec(nd)
    1708             :         integer(IK)              :: i
    1709             : 
    1710          21 :         CholeskyLower = CovMat
    1711           3 :         call getCholeskyFactor(nd,CholeskyLower,Diagonal)
    1712           3 :         if (Diagonal(1)<0._RK) then
    1713             :         ! LCOV_EXCL_START
    1714             :             write(output_unit,"(A)") "getCholeskyFactor() failed in getMVNDev()"
    1715             :             error stop
    1716             :         end if
    1717             :         ! LCOV_EXCL_STOP
    1718           9 :         do i=1,nd
    1719           6 :             DummyVec(i) = getRandGaus()
    1720           9 :             x(i) = DummyVec(i) * Diagonal(i)
    1721             :         end do
    1722           6 :         do i=2,nd
    1723           9 :             x(i) = x(i) + dot_product(CholeskyLower(i,1:i-1),DummyVec(1:i-1))
    1724             :         end do
    1725           9 :         x = x + MeanVec
    1726             : 
    1727           3 :   end subroutine getMVNDev
    1728             : 
    1729             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1730             : 
    1731             :     ! This subroutine is legacy and slow. use getRandMVU() in this same module.
    1732             :     ! Given the mean vector MeanVec and the covariance matrix CovMat, this subroutine generates a random vector X (of length nd>=2)
    1733             :     ! from an nd-dimensional multivariate ellipsoidal uniform distribution, such that getMVUDev() is randomly distributed inside the nd-dimensional ellipsoid.
    1734             :     ! ATTENTION: Only the upper half of the covariance matrix (plus the diagonal terms) need to be given in the input.
    1735             :     ! in the ouput, the upper half and diagonal part will still be the covariance matrix, while the lower half will be
    1736             :     ! the Cholesky decomposition elements (excluding its diagonal terms that are provided only in the vector Diagonal).
    1737             :     ! USES getCholeskyFactor.f90, getRandGaus.f90
    1738             :     !> Amir Shahmoradi, April 25, 2016, 2:21 PM, IFS, UTEXAS
    1739           3 :     subroutine getMVUDev(nd,MeanVec,CovMat,X)
    1740             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1741             :         !DEC$ ATTRIBUTES DLLEXPORT :: getMVUDev
    1742             : #endif
    1743             : 
    1744           3 :         use Matrix_mod, only: getCholeskyFactor
    1745             : 
    1746             :         implicit none
    1747             :         integer(IK), intent(in)  :: nd
    1748             :         real(RK)   , intent(in)  :: MeanVec(nd)
    1749             :         real(RK)   , intent(in)  :: CovMat(nd,nd)
    1750             :         real(RK)   , intent(out) :: X(nd)
    1751          33 :         real(RK)                 :: Diagonal(nd), DummyVec(nd), CholeskyLower(nd,nd), dummy
    1752             :         integer(IK)              :: i
    1753             : 
    1754          21 :         CholeskyLower = CovMat
    1755           3 :         call getCholeskyFactor(nd,CholeskyLower,Diagonal)
    1756           3 :         if (Diagonal(1)<0._RK) then
    1757             :         ! LCOV_EXCL_START
    1758             :             error stop
    1759             :             !call abortProgram( output_unit , 1 , 1 , 'Statitistics@getMVUDev()@getCholeskyFactor() failed.' )
    1760             :         end if
    1761             :         ! LCOV_EXCL_STOP
    1762           9 :         do i=1,nd
    1763           9 :             DummyVec(i) = getRandGaus()
    1764             :         end do
    1765           3 :         call random_number(dummy)
    1766           9 :         dummy = (dummy**(1._RK/real(nd,kind=RK)))/norm2(DummyVec)  ! Now DummyVec is a uniformly drawn point from inside of nd-D sphere.
    1767           9 :         DummyVec = dummy*DummyVec
    1768             : 
    1769             :         ! Now transform this point to a point inside the ellipsoid.
    1770           9 :         do i=1,nd
    1771           9 :             X(i) = DummyVec(i)*Diagonal(i)
    1772             :         end do
    1773             : 
    1774           6 :         do i=2,nd
    1775           9 :             X(i) = X(i) + dot_product(CholeskyLower(i,1:i-1),DummyVec(1:i-1))
    1776             :         end do
    1777             : 
    1778           9 :         X = X + MeanVec
    1779             : 
    1780           3 :     end subroutine getMVUDev
    1781             : 
    1782             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1783             : 
    1784             :     !> \brief
    1785             :     !> Return a MultiVariate Normal (MVN) random vector with the given mean and
    1786             :     !> covariance matrix represented by the the input Cholesky factorization.
    1787             :     !>
    1788             :     !> \param[in]   nd              :   The number of dimensions of the MVN distribution.
    1789             :     !> \param[in]   MeanVec         :   The mean vector of the MVN distribution.
    1790             :     !> \param[in]   CholeskyLower   :   The Cholesky lower triangle of the covariance matrix of the MVN distribution.
    1791             :     !> \param[in]   Diagonal        :   The Diagonal elements of the Cholesky lower triangle of the covariance matrix of the MVN distribution.
    1792             :     !>
    1793             :     !> \return
    1794             :     !> `RandMVN` : The randomly generated MVN vector.
    1795             :     !>
    1796             :     !> \author
    1797             :     !> Amir Shahmoradi, April 23, 2017, 12:36 AM, ICES, UTEXAS
    1798     1587650 :     function getRandMVN(nd,MeanVec,CholeskyLower,Diagonal) result(RandMVN)
    1799             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1800             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandMVN
    1801             : #endif
    1802             :         implicit none
    1803             :         integer(IK), intent(in) :: nd
    1804             :         real(RK)   , intent(in) :: MeanVec(nd)
    1805             :         real(RK)   , intent(in) :: CholeskyLower(nd,nd), Diagonal(nd)   ! Cholesky lower triangle and its diagonal terms, calculated from the input CovMat.
    1806      337279 :         real(RK)                :: RandMVN(nd), dummy
    1807             :         integer(IK)             :: j,i
    1808      913094 :         RandMVN = MeanVec
    1809      913094 :         do j = 1,nd
    1810      575815 :             dummy = getRandGaus()
    1811      575815 :             RandMVN(j) = RandMVN(j) + Diagonal(j) * dummy
    1812     1151630 :             do i = j+1,nd
    1813      814351 :                 RandMVN(i) = RandMVN(i) + CholeskyLower(i,j) * dummy
    1814             :             end do
    1815             :         end do
    1816      337282 :     end function getRandMVN
    1817             : 
    1818             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1819             : 
    1820             : !    ! Given an input Mean vector of length nd, Covariance Matrix of dimension (nd,nd), as well as a vector of integers representing
    1821             : !    ! the variables (corresponding to CovMat columns) that are given
    1822             : !    ! This subroutine gives out a conditional Multivariate Normal Random deviate.
    1823             : !    ! random p-tivariate normal deviate, given that the first pg variables x1 are given (i.e. fixed).
    1824             : !    ! For a review of Multivariate Normal distribution: Applied Multivariate Statistical Analysis, Johnson, Wichern, 1998, 4th ed.
    1825             : !    !> Amir Shahmoradi, Oct 20, 2009, 9:12 PM, MTU
    1826             : !    function getCondRandMVN(nd,MeanVec,CovMat,nIndIndx,IndIndx) result(CondRandMVN)
    1827             : !        use Matrix_mod, only: getRegresCoef
    1828             : !        implicit none
    1829             : !        integer(IK), intent(in) :: nd, nIndIndx, IndIndx(nIndIndx)
    1830             : !        real(RK)   , intent(in) :: MeanVec(nd), CovMat(nd,nd)
    1831             : !        real(RK)                :: CondRandMVN(nd),
    1832             : !        integer(IK)             :: j, i
    1833             : !    end function getCondRandMVN
    1834             : 
    1835             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1836             : 
    1837             :     !> \brief
    1838             :     !> Given the Cholesky Lower triangle and diagonals of a given covariance matrix, this function return one point uniformly
    1839             :     !> randomly drawn from inside of an `nd`-ellipsoid, whose `nd` elements `RandMVU(i), i=1:nd` are guaranteed
    1840             :     !> to be in the range:
    1841             :     !> ```
    1842             :     !> MeanVec(i) - sqrt(CovMat(i,i)) < RandMVU(i) < MeanVec(i) + sqrt(CovMat(i,i))
    1843             :     !> ```
    1844             :     !>
    1845             :     !> \param[in]   nd              :   The number of dimensions of the MVU distribution.
    1846             :     !> \param[in]   MeanVec         :   The mean vector of the MVU distribution.
    1847             :     !> \param[in]   CholeskyLower   :   The Cholesky lower triangle of the covariance matrix of the MVU distribution.
    1848             :     !> \param[in]   Diagonal        :   The Diagonal elements of the Cholesky lower triangle of the covariance matrix of the MVU distribution.
    1849             :     !>
    1850             :     !> \return
    1851             :     !> `RandMVU` : The randomly generated MVU vector.
    1852             :     !>
    1853             :     !> \author
    1854             :     !> Amir Shahmoradi, April 23, 2017, 1:36 AM, ICES, UTEXAS
    1855      716850 :     function getRandMVU(nd,MeanVec,CholeskyLower,Diagonal) result(RandMVU)
    1856             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1857             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandMVU
    1858             : #endif
    1859             :         implicit none
    1860             :         integer(IK), intent(in) :: nd
    1861             :         real(RK)   , intent(in) :: MeanVec(nd)
    1862             :         real(RK)   , intent(in) :: CholeskyLower(nd,nd) ! Cholesky lower triangle, calculated from the input CovMat.
    1863             :         real(RK)   , intent(in) :: Diagonal(nd)         ! Cholesky diagonal terms, calculated from the input CovMat.
    1864      572736 :         real(RK)                :: RandMVU(nd), dummy, DummyVec(nd), sumSqDummyVec
    1865             :         integer(IK)             :: i,j
    1866      144114 :         sumSqDummyVec = 0._RK
    1867      428622 :         do j=1,nd
    1868      284508 :             DummyVec(j) = getRandGaus()
    1869      428622 :             sumSqDummyVec = sumSqDummyVec + DummyVec(j)**2
    1870             :         end do
    1871      144114 :         call random_number(dummy)
    1872      144114 :         dummy = dummy**(1._RK/nd) / sqrt(sumSqDummyVec)
    1873      428622 :         DummyVec = DummyVec * dummy  ! a uniform random point from inside of nd-sphere
    1874      428622 :         RandMVU = MeanVec
    1875      428622 :         do j = 1,nd
    1876      284508 :             RandMVU(j) = RandMVU(j) + Diagonal(j) * DummyVec(j)
    1877      569016 :             do i = j+1,nd
    1878      424902 :                 RandMVU(i) = RandMVU(i) + CholeskyLower(i,j) * DummyVec(j)
    1879             :             end do
    1880             :         end do
    1881      481393 :     end function getRandMVU
    1882             : 
    1883             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1884             : 
    1885             :     !> \brief
    1886             :     !> Return `.true.` if the input `NormedPoint` (normalized with respect to the center of the target ellipsoid) is
    1887             :     !> within or on the boundary of the ellipsoid whose boundary is described by the representative matrix
    1888             :     !> \f$ \Sigma \f$ (`RepMat`), such that,
    1889             :     !> \f{equation}{
    1890             :     !> X^T ~ \Sigma^{-1} ~ X = 1 ~,
    1891             :     !> \f}
    1892             :     !> for all \f$X\f$ on the boundary.
    1893             :     !>
    1894             :     !> \param[in]   nd          :   The number of dimensions of the ellipsoid (the size of `NormedPoint`).
    1895             :     !> \param[in]   NormedPoint :   The input point, normalized with respect to the center of the ellipsoid,
    1896             :     !>                                  whose location with respect to the ellipsoid boundary is to be determined.
    1897             :     !> \param[in]   InvRepMat   :   The inverse of the representative matrix of the target ellipsoid.
    1898             :     !>
    1899             :     !> \return
    1900             :     !> `isInsideEllipsoid` : The logical value indicating whether the input point is inside or on the boundary of the target ellipsoid.
    1901             :     !>
    1902             :     !> \remark
    1903             :     !> Note that the input matrix is the inverse of `RepMat`: `InvRepMat`.
    1904             :     !>
    1905             :     !> \author
    1906             :     !> Amir Shahmoradi, April 23, 2017, 1:36 AM, ICES, UTEXAS
    1907          15 :     pure function isInsideEllipsoid(nd,NormedPoint,InvRepMat)
    1908             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1909             :         !DEC$ ATTRIBUTES DLLEXPORT :: isInsideEllipsoid
    1910             : #endif
    1911      144114 :         use Math_mod, only: getLogVolEllipsoid
    1912             :         implicit none
    1913             :         integer(IK), intent(in) :: nd
    1914             :         real(RK)   , intent(in) :: NormedPoint(nd)
    1915             :         real(RK)   , intent(in) :: InvRepMat(nd,nd)
    1916             :         logical                 :: isInsideEllipsoid
    1917          45 :         isInsideEllipsoid = dot_product(NormedPoint,matmul(InvRepMat,NormedPoint)) <= 1._RK
    1918          15 :     end function isInsideEllipsoid
    1919             : 
    1920             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1921             : 
    1922             :     !> \brief
    1923             :     !> Return the natural logarithm of the probability density function value of a point uniformly distributed within an ellipsoid,
    1924             :     !> whose logarithm of the square root of the determinant of its representative covariance matrix is given by `logSqrtDetCovMat`.
    1925             :     !>
    1926             :     !> \param[in]   nd                  :   The number of dimensions of the MVU distribution.
    1927             :     !> \param[in]   logSqrtDetCovMat    :   The logarithm of the square root of the determinant of
    1928             :     !>                                      the inverse of the representative covariance matrix of the ellipsoid.
    1929             :     !>
    1930             :     !> \return
    1931             :     !> `logProbMVU` :   The natural logarithm of the probability density function value
    1932             :     !>                  of a point uniformly distributed within the target ellipsoid.
    1933             :     !>
    1934             :     !> \author
    1935             :     !> Amir Shahmoradi, April 23, 2017, 1:36 AM, ICES, UTEXAS
    1936           3 :     pure function getLogProbMVU(nd,logSqrtDetCovMat) result(logProbMVU)
    1937             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1938             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMVU
    1939             : #endif
    1940          15 :         use Math_mod, only: getLogVolEllipsoid
    1941             :         implicit none
    1942             :         integer(IK), intent(in) :: nd
    1943             :         real(RK)   , intent(in) :: logSqrtDetCovMat
    1944             :         real(RK)                :: logProbMVU
    1945           3 :         logProbMVU = -getLogVolEllipsoid(nd,logSqrtDetCovMat)
    1946           6 :     end function getLogProbMVU
    1947             : 
    1948             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1949             : 
    1950             :     !> \brief
    1951             :     !> Return a random point on the target ellipsoid by projecting a random point uniformly distributed within the ellipsoid on its boundary.
    1952             :     !>
    1953             :     !> \param[in]   nd              :   The number of dimensions of the ellipsoid.
    1954             :     !> \param[in]   MeanVec         :   The mean vector (center) of the ellipsoid.
    1955             :     !> \param[in]   CholeskyLower   :   The Cholesky lower triangle of the representative covariance matrix of the ellipsoid.
    1956             :     !> \param[in]   Diagonal        :   The Diagonal elements of the Cholesky lower triangle of the representative covariance matrix of the ellipsoid.
    1957             :     !>
    1958             :     !> \return
    1959             :     !> `RandPointOnEllipsoid` : A random point on the target ellipsoid by projecting a random
    1960             :     !>                          point uniformly distributed within the ellipsoid on its boundary.
    1961             :     !>
    1962             :     !> \remark
    1963             :     !> This is algorithm is similar to [getRandMVU](@ref getrandmvu), with the only difference that
    1964             :     !> points are drawn randomly from the surface of the ellipsoid instead of inside of its interior.
    1965             :     !>
    1966             :     !> \remark
    1967             :     !> Note that the distribution of points on the surface of the ellipsoid is **NOT** uniform.
    1968             :     !> Regions of high curvature will have more points randomly sampled from them.
    1969             :     !> Generating uniform random points on arbitrary-dimension ellipsoids is not a task with trivial solution!
    1970             :     !>
    1971             :     !> \todo
    1972             :     !> The performance of this algorithm can be further improved.
    1973             :     !>
    1974             :     !> \author
    1975             :     !> Amir Shahmoradi, April 23, 2017, 1:36 AM, ICES, UTEXAS
    1976          15 :     function getRandPointOnEllipsoid(nd,MeanVec,CholeskyLower,Diagonal) result(RandPointOnEllipsoid)
    1977             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    1978             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandPointOnEllipsoid
    1979             : #endif
    1980             :         implicit none
    1981             :         integer(IK), intent(in) :: nd
    1982             :         real(RK)   , intent(in) :: MeanVec(nd)
    1983             :         real(RK)   , intent(in) :: CholeskyLower(nd,nd) ! Cholesky lower triangle, calculated from the MVN CovMat.
    1984             :         real(RK)   , intent(in) :: Diagonal(nd)         ! Cholesky diagonal terms, calculated from the MVN CovMat.
    1985          12 :         real(RK)                :: RandPointOnEllipsoid(nd), DummyVec(nd), sumSqDummyVec
    1986             :         integer(IK)             :: i,j
    1987           3 :         sumSqDummyVec = 0._RK
    1988           9 :         do j=1,nd
    1989           6 :             DummyVec(j) = getRandGaus()
    1990           9 :             sumSqDummyVec = sumSqDummyVec + DummyVec(j)**2
    1991             :         end do
    1992           9 :         DummyVec = DummyVec / sqrt(sumSqDummyVec)    ! DummyVec is a random point on the surface of nd-sphere.
    1993           9 :         RandPointOnEllipsoid = 0._RK
    1994           9 :         do j = 1,nd
    1995           6 :             RandPointOnEllipsoid(j) = RandPointOnEllipsoid(j) + Diagonal(j) * DummyVec(j)
    1996          12 :             do i = j+1,nd
    1997           9 :                 RandPointOnEllipsoid(i) = RandPointOnEllipsoid(i) + CholeskyLower(i,j) * DummyVec(j)
    1998             :             end do
    1999             :         end do
    2000           9 :         RandPointOnEllipsoid = RandPointOnEllipsoid + MeanVec
    2001           6 :     end function getRandPointOnEllipsoid
    2002             : 
    2003             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2004             : 
    2005             :     !> \brief
    2006             :     !> Return the natural logarithm of the Lognormal probability density function.
    2007             :     !>
    2008             :     !> \param[in]   logMean                 :   The mean of the Lognormal distribution.
    2009             :     !> \param[in]   inverseVariance         :   The inverse variance of the Lognormal distribution.
    2010             :     !> \param[in]   logSqrtInverseVariance  :   The natural logarithm of the square root of the inverse variance of the Lognormal distribution.
    2011             :     !> \param[in]   logPoint                :   The natural logarithm of the point at which the Lognormal PDF must be computed.
    2012             :     !>
    2013             :     !> \return
    2014             :     !> `logProbLogn` : The natural logarithm of the Lognormal probability density function.
    2015             :     !>
    2016             :     !> \author
    2017             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2018        3894 :     pure function getLogProbLognSP(logMean,inverseVariance,logSqrtInverseVariance,logPoint) result(logProbLogn)
    2019             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2020             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbLognSP
    2021             : #endif
    2022           3 :         use Constants_mod, only: LOGINVSQRT2PI
    2023             :         implicit none
    2024             :         real(RK), intent(in) :: logMean,inverseVariance,logSqrtInverseVariance,logPoint
    2025             :         real(RK)             :: logProbLogn
    2026        3894 :         logProbLogn = LOGINVSQRT2PI + logSqrtInverseVariance - logPoint - 0.5_RK * inverseVariance * (logPoint-logMean)**2
    2027        7788 :     end function getLogProbLognSP
    2028             : 
    2029             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2030             : 
    2031             :     !> \brief
    2032             :     !> Return the natural logarithm of the Lognormal probability density function.
    2033             :     !>
    2034             :     !> \param[in]   np                      :   The size of the input vector of points represented by `LogPoint`.
    2035             :     !> \param[in]   logMean                 :   The mean of the Lognormal distribution.
    2036             :     !> \param[in]   inverseVariance         :   The inverse variance of the Lognormal distribution.
    2037             :     !> \param[in]   logSqrtInverseVariance  :   The natural logarithm of the square root of the inverse variance of the Lognormal distribution.
    2038             :     !> \param[in]   LogPoint                :   The natural logarithm of the vector of points at which the Lognormal PDF must be computed.
    2039             :     !>
    2040             :     !> \return
    2041             :     !> `logProbLogn` : The natural logarithm of the Lognormal probability density function.
    2042             :     !>
    2043             :     !> \author
    2044             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2045          36 :     pure function getLogProbLognMP(np,logMean,inverseVariance,logSqrtInverseVariance,LogPoint) result(logProbLogn)
    2046             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2047             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbLognMP
    2048             : #endif
    2049        3894 :         use Constants_mod, only: LOGINVSQRT2PI
    2050             :         implicit none
    2051             :         integer(IK), intent(in) :: np
    2052             :         real(RK)   , intent(in) :: logMean,inverseVariance,logSqrtInverseVariance,LogPoint(np)
    2053             :         real(RK)                :: logProbLogn(np)
    2054          24 :         logProbLogn = LOGINVSQRT2PI + logSqrtInverseVariance - LogPoint - 0.5_RK * inverseVariance * (LogPoint-logMean)**2
    2055           6 :     end function getLogProbLognMP
    2056             : 
    2057             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2058             : 
    2059             :     !> \brief
    2060             :     !> Return a single-precision uniformly-distributed random real-valued number in the range `[0,1]`.
    2061             :     !>
    2062             :     !> \param[inout]    idum    :   The input integer random seed with the `save` attribute.
    2063             :     !>
    2064             :     !> \return
    2065             :     !> `randRealLecuyer`        :   A single-precision uniformly-distributed random real-valued number in the range `[0,1]`.
    2066             :     !>
    2067             :     !> \warning
    2068             :     !> Do not change the value of `idum` between calls.
    2069             :     !>
    2070             :     !> \remark
    2071             :     !> This routine is guaranteed to random numbers with priodicity larger than `~2*10**18` random numbers.
    2072             :     !> For more info see Press et al. (1990) Numerical Recipes.
    2073             :     !>
    2074             :     !> \remark
    2075             :     !> This routine is solely kept for backward compatibility reasons.
    2076             :     !> The Fortran intrinsic subroutine `random_number()` is recommended to be used against this function.
    2077             :     !>
    2078             :     !> \author
    2079             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2080         600 :     function getRandRealLecuyer(idum) result(randRealLecuyer)
    2081             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2082             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandRealLecuyer
    2083             : #endif
    2084             :         implicit none
    2085             :         integer(IK), intent(inout) :: idum
    2086             :         integer(IK), parameter     :: im1=2147483563, im2=2147483399, imm1=im1-1, ia1=40014, ia2=40692
    2087             :         integer(IK), parameter     :: iq1=53668, iq2=52774, ir1=12211, ir2=3791, ntab=32, ndiv=1+imm1/ntab
    2088             :         real(RK)   , parameter     :: am=1._RK/real(im1,kind=RK), eps=1.2e-7_RK, rnmx=1._RK-eps
    2089             :         real(RK)                   :: randRealLecuyer
    2090             :         integer(IK)                :: idum2,j,k,iv(ntab),iy
    2091             :         save                       :: iv, iy, idum2
    2092             :         data idum2/123456789/, iv/ntab*0/, iy/0/
    2093         600 :         if (idum <= 0) then
    2094           0 :             idum = max(-idum,1)
    2095           0 :             idum2 = idum
    2096           0 :             do j = ntab+8,1,-1
    2097           0 :                 k = idum/iq1
    2098           0 :                 idum = ia1*(idum-k*iq1)-k*ir1
    2099           0 :                 if (idum < 0) idum = idum+im1
    2100           0 :                 if (j <= ntab) iv(j) = idum
    2101             :             end do
    2102           0 :             iy = iv(1)
    2103             :         endif
    2104         600 :         k = idum/iq1
    2105         600 :         idum = ia1*(idum-k*iq1)-k*ir1
    2106         600 :         if (idum < 0) idum=idum+im1
    2107         600 :         k = idum2/iq2
    2108         600 :         idum2 = ia2*(idum2-k*iq2)-k*ir2
    2109         600 :         if (idum2 < 0) idum2=idum2+im2
    2110         600 :         j = 1+iy/ndiv
    2111         600 :         iy = iv(j)-idum2
    2112         600 :         iv(j) = idum
    2113         600 :         if(iy < 1)iy = iy+imm1
    2114         600 :         randRealLecuyer = min(am*iy,rnmx)
    2115         606 :     end function getRandRealLecuyer
    2116             : 
    2117             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2118             : 
    2119             :     !> \brief
    2120             :     !> Return an integer uniformly-distributed random integer-valued number in the range `[lowerBound , upperBound]`.
    2121             :     !>
    2122             :     !> \param[in]       lowerBound  :   The inclusive integer lower bound of the integer uniform distribution.
    2123             :     !> \param[in]       upperBound  :   The inclusive integer upper bound of the integer uniform distribution.
    2124             :     !> \param[inout]    idum        :   The input integer random seed with the `save` attribute.
    2125             :     !>
    2126             :     !> \return
    2127             :     !> `randRealLecuyer`    : A uniformly-distributed random integer-valued number in the range `[lowerBound , upperBound]`.
    2128             :     !>
    2129             :     !> \warning
    2130             :     !> Do not change the value of `idum` between calls.
    2131             :     !>
    2132             :     !> \remark
    2133             :     !> This routine is guaranteed to random numbers with priodicity larger than `~2*10**18` random numbers.
    2134             :     !> For more info see Press et al. (1990) Numerical Recipes.
    2135             :     !>
    2136             :     !> \remark
    2137             :     !> This routine is solely kept for backward compatibility reasons.
    2138             :     !> The [getRandInt](@ref getrandint) is recommended to be used against this routine.
    2139             :     !>
    2140             :     !> \author
    2141             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2142         300 :     function getRandIntLecuyer(lowerBound,upperBound,idum)
    2143             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2144             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandIntLecuyer
    2145             : #endif
    2146             :         implicit none
    2147             :         integer(IK), intent(in)    :: lowerBound,upperBound
    2148             :         integer(IK), intent(inout) :: idum
    2149             :         integer(IK)                :: getRandIntLecuyer
    2150         600 :         getRandIntLecuyer = lowerBound + nint( getRandRealLecuyer(idum)*real(upperBound-lowerBound,kind=RK) )
    2151         900 :     end function getRandIntLecuyer
    2152             : 
    2153             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2154             : 
    2155             :     !> \brief
    2156             :     !> Return an integer uniformly-distributed random integer-valued number in the range `[lowerBound , upperBound]`.
    2157             :     !>
    2158             :     !> \param[in]       lowerBound  :   The lower bound of the integer uniform distribution.
    2159             :     !> \param[in]       upperBound  :   The upper bound of the integer uniform distribution.
    2160             :     !>
    2161             :     !> \return
    2162             :     !> `randInt` : A uniformly-distributed random integer-valued number in the range `[lowerBound , upperBound]`.
    2163             :     !>
    2164             :     !> \author
    2165             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2166       60113 :     function getRandInt(lowerBound,upperBound) result(randInt)
    2167             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2168             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandInt
    2169             : #endif
    2170             :         implicit none
    2171             :         integer(IK), intent(in) :: lowerBound,upperBound
    2172       60113 :         real(RK)                :: dummy
    2173             :         integer(IK)             :: randInt
    2174       60113 :         call random_number(dummy)
    2175       60113 :         randInt = lowerBound + nint( dummy*real(upperBound-lowerBound,kind=RK) )
    2176       60413 :     end function getRandInt
    2177             : 
    2178             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2179             : 
    2180             :     !> \brief
    2181             :     !> Return an integer uniformly-distributed random integer-valued number in the range `[lowerBound , upperBound]` using
    2182             :     !> the built-in random number generator of Fortran.
    2183             :     !>
    2184             :     !> \param[in]       lowerBound  :   The lower bound of the integer uniform distribution.
    2185             :     !> \param[in]       upperBound  :   The upper bound of the integer uniform distribution.
    2186             :     !>
    2187             :     !> \return
    2188             :     !> `randUniform`    : A uniformly-distributed random real-valued number in the range `[lowerBound , upperBound]`.
    2189             :     !>
    2190             :     !> \author
    2191             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2192         300 :     function getRandUniform(lowerBound,upperBound) result(randUniform)
    2193             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2194             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandUniform
    2195             : #endif
    2196             :         implicit none
    2197             :         real(RK), intent(in)    :: lowerBound, upperBound
    2198             :         real(RK)                :: randUniform
    2199         300 :         call random_number(randUniform)
    2200         300 :         randUniform = lowerBound + randUniform * (upperBound - lowerBound)
    2201       60413 :     end function getRandUniform
    2202             : 
    2203             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2204             : 
    2205             :     !> \brief
    2206             :     !> Return a Gamma-distributed random number, following the prescription in the GSL library.
    2207             :     !>
    2208             :     !> \param[in]       alpha   :   The shape parameter of the Gamma distribution.
    2209             :     !>
    2210             :     !> \return
    2211             :     !> `randGamma`    : A Gamma-distributed random real-valued number in the range `[0,+Infinity]`.
    2212             :     !>
    2213             :     !> \warning
    2214             :     !> The condition `alpha > 0` must hold, otherwise the negative value `-1` will be returned to indicate the occurrence of error.
    2215             :     !>
    2216             :     !> \author
    2217             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2218        1503 :     function getRandGamma(alpha) result(randGamma)
    2219             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2220             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandGamma
    2221             : #endif
    2222             :         implicit none
    2223             :         real(RK), intent(in) :: alpha
    2224             :         real(RK)             :: randGamma
    2225        1503 :         real(RK)             :: c,u,v,z
    2226        1503 :         if (alpha<=0._RK) then  ! illegal value of alpha
    2227           3 :             randGamma = -1._RK
    2228           3 :             return
    2229             :         else
    2230        1500 :             randGamma = alpha
    2231        1500 :             if (randGamma<1._RK) randGamma = randGamma + 1._RK
    2232        1500 :             randGamma = randGamma - 0.3333333333333333_RK
    2233        1500 :             c = 3._RK*sqrt(randGamma)
    2234        1500 :             c = 1._RK / c
    2235          22 :             do
    2236           0 :                 do
    2237        1522 :                     z = getRandGaus()
    2238        1522 :                     v = 1._RK + c*z
    2239        1522 :                     if (v<=0._RK) cycle
    2240        1522 :                     exit
    2241             :                 end do
    2242        1522 :                 v = v**3
    2243        1522 :                 call random_number(u)
    2244        1522 :                 if ( log(u) >= 0.5_RK * z**2 + randGamma * ( 1._RK - v + log(v) ) ) cycle
    2245        1500 :                 randGamma = randGamma * v
    2246        1500 :                 exit
    2247             :             end do
    2248        1500 :             if (alpha<1._RK) then
    2249           0 :                 call random_number(u)
    2250           0 :                 randGamma = randGamma * u**(1._RK/alpha)
    2251             :             end if
    2252             :         end if
    2253        1803 :     end function getRandGamma
    2254             : 
    2255             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2256             : 
    2257             :     !> \brief
    2258             :     !> Return a Gamma-distributed random number, whose shape parameter `alpha` is an integer.
    2259             :     !>
    2260             :     !> \param[in]       alpha   :   The shape integer parameter of the Gamma distribution.
    2261             :     !>
    2262             :     !> \return
    2263             :     !> `randGamma`    : A Gamma-distributed random real-valued number in the range `[0,+Infinity]`.
    2264             :     !>
    2265             :     !> \warning
    2266             :     !> The condition `alpha > 1` must hold, otherwise the negative value `-1` will be returned to indicate the occurrence of error.
    2267             :     !>
    2268             :     !> \author
    2269             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2270         303 :     function getRandGammaIntShape(alpha) result(randGamma)
    2271             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2272             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandGammaIntShape
    2273             : #endif
    2274             :         implicit none
    2275             :         integer(IK), intent(in) :: alpha
    2276             :         real(RK)                :: randGamma
    2277         303 :         real(RK)                :: am,e,h,s,x,y,Vector(2),Array(5)
    2278         303 :         if (alpha < 1) then  ! illegal value of alpha
    2279           3 :             randGamma = -1._RK
    2280           3 :             return
    2281         300 :         elseif (alpha < 6) then
    2282         300 :             call random_number(Array(1:alpha))
    2283         900 :             x = -log(product(Array(1:alpha)))
    2284             :         else    ! use rejection sampling
    2285           0 :             do
    2286           0 :                 call random_number(Vector)
    2287           0 :                 Vector(2) = 2._RK*Vector(2)-1._RK
    2288           0 :                 if (dot_product(Vector,Vector) > 1._RK) cycle
    2289           0 :                 y = Vector(2) / Vector(1)
    2290           0 :                 am = alpha - 1
    2291           0 :                 s = sqrt(2._RK*am + 1._RK)
    2292           0 :                 x = s*y + am
    2293           0 :                 if (x <= 0.0) cycle
    2294           0 :                 e = (1._RK+y**2) * exp(am*log(x/am)-s*y)
    2295           0 :                 call random_number(h)
    2296           0 :                 if (h <= e) exit
    2297             :             end do
    2298             :         end if
    2299         300 :         randGamma = x
    2300        1806 :     end function getRandGammaIntShape
    2301             : 
    2302             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2303             : 
    2304             :     !> \brief
    2305             :     !> Return a random Beta-distributed variable.
    2306             :     !>
    2307             :     !> \param[in]       alpha   :   The first shape parameter of the Beta distribution.
    2308             :     !> \param[in]       beta    :   The second shape parameter of the Beta distribution.
    2309             :     !>
    2310             :     !> \return
    2311             :     !> `randBeta`   : A Beta-distributed random real-valued number in the range `[0,1]`.
    2312             :     !>
    2313             :     !> \warning
    2314             :     !> The conditions `alpha > 0` and `beta > 0` must hold, otherwise the negative
    2315             :     !> value `-1` will be returned to indicate the occurrence of error.
    2316             :     !>
    2317             :     !> \author
    2318             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2319         609 :     function getRandBeta(alpha,beta) result(randBeta)
    2320             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2321             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandBeta
    2322             : #endif
    2323             :         implicit none
    2324             :         real(RK), intent(in) :: alpha,beta
    2325             :         real(RK)             :: randBeta
    2326         609 :         real(RK)             :: x
    2327         609 :         if ( alpha>0._RK .and. beta>0._RK ) then
    2328         600 :             x = getRandGamma(alpha)
    2329         600 :             randBeta = x / ( x + getRandGamma(beta) )
    2330             :         else ! illegal value of alpha or beta
    2331           9 :             randBeta = -1._RK
    2332             :         end if
    2333         912 :     end function getRandBeta
    2334             : 
    2335             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2336             : 
    2337             :     !> \brief
    2338             :     !> Return a random Exponential-distributed value whose inverse mean is given as input.
    2339             :     !>
    2340             :     !> \param[in]   invMean :   The inverse of the mean of the Exponential distribution.
    2341             :     !>
    2342             :     !> \return
    2343             :     !> `randExp` : An Exponential-distributed random real-valued number in the range `[0,+Infinity]` with mean `1 / invMean`.
    2344             :     !>
    2345             :     !> \warning
    2346             :     !> It is the user's onus to ensure `invMean > 0` on input. This condition will NOT be checked by this routine.
    2347             :     !>
    2348             :     !> \author
    2349             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2350         300 :     function getRandExpWithInvMean(invMean) result(randExp)
    2351             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2352             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandExpWithInvMean
    2353             : #endif
    2354             :         implicit none
    2355             :         real(RK), intent(in)    :: invMean
    2356             :         real(RK)                :: randExp
    2357         300 :         call random_number(randExp)
    2358         300 :         randExp = -log(randExp) * invMean
    2359         909 :     end function getRandExpWithInvMean
    2360             : 
    2361             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2362             : 
    2363             :     !> \brief
    2364             :     !> Return a random Exponential-distributed value whose mean is \f$\lambda = 1\f$.
    2365             :     !>
    2366             :     !> \return
    2367             :     !> `randExp` : A random Exponential-distributed value whose mean \f$\lambda = 1\f$.
    2368             :     !>
    2369             :     !> \author
    2370             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2371         300 :     function getRandExp() result(randExp)
    2372             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2373             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandExp
    2374             : #endif
    2375             :         implicit none
    2376             :         real(RK) :: randExp
    2377         300 :         call random_number(randExp)
    2378         300 :         randExp = -log(randExp)
    2379         600 :     end function getRandExp
    2380             : 
    2381             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2382             : 
    2383             :     !> \brief
    2384             :     !> Return a random correlation matrix.
    2385             :     !>
    2386             :     !> \param[in]   nd  :   The rank of the correlation matrix.
    2387             :     !> \param[in]   eta :   The parameter that roughly represents the shape parameters of the Beta distribution.
    2388             :     !>                      The larger the value of `eta` is, the more homogeneous the correlation matrix will look.
    2389             :     !>                      In other words, set this parameter to some small number to generate strong random correlations
    2390             :     !>                      in the output random correlation matrix.
    2391             :     !>
    2392             :     !> \return
    2393             :     !> `RandCorMat` : A random correlation matrix.
    2394             :     !>
    2395             :     !> \warning
    2396             :     !> The conditions `nd > 1` and `eta > 0.0` must hold, otherwise the first element of
    2397             :     !> output, `getRandCorMat(1,1)`, will be set to `-1` to indicate the occurrence of an error.
    2398             :     !>
    2399             :     !> \author
    2400             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2401        2781 :     function getRandCorMat(nd,eta) result(RandCorMat)
    2402             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2403             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandCorMat
    2404             : #endif
    2405         300 :         use Matrix_mod, only: getCholeskyFactor
    2406             :         implicit none
    2407             :         integer(IK), intent(in) :: nd
    2408             :         real(RK)   , intent(in) :: eta
    2409         309 :         real(RK)                :: RandCorMat(nd,nd), dummy
    2410        1545 :         real(RK)                :: beta,sumSqDummyVec,DummyVec(nd-1),W(nd-1),Diagonal(nd-1)
    2411             :         integer(IK)             :: m,j,i
    2412             : 
    2413         309 :         if (nd<2_IK .or. eta<=0._RK) then  ! illegal value for eta.
    2414           9 :             RandCorMat(1,1) = -1._RK
    2415           9 :             return
    2416             :         end if
    2417             : 
    2418         900 :         do m = 1,nd
    2419         900 :             RandCorMat(m,m) = 1._RK
    2420             :         end do
    2421         300 :         beta = eta + 0.5_RK*(nd-2._RK)
    2422         300 :         dummy = getRandBeta(beta,beta)
    2423         300 :         if (dummy<=0._RK .or. dummy>=1._RK) then
    2424             :         ! LCOV_EXCL_START
    2425             :             error stop
    2426             :             !call abortProgram( output_unit , 1 , 1 , 'Statitistics@getRandCorMat() failed. Random Beta variable out of bound: ' // num2str(dummy) )
    2427             :         end if
    2428             :         ! LCOV_EXCL_STOP
    2429         300 :         RandCorMat(1,2) = 2._RK * dummy - 1._RK ! for the moment, only the upper half of RandCorMat is needed, the lower half will contain cholesky lower triangle.
    2430             : 
    2431         300 :         do m = 2,nd-1
    2432           0 :             beta = beta - 0.5_RK
    2433           0 :             sumSqDummyVec = 0._RK
    2434           0 :             do j=1,m
    2435           0 :                 DummyVec(j) = getRandGaus()
    2436           0 :                 sumSqDummyVec = sumSqDummyVec + DummyVec(j)**2
    2437             :             end do
    2438           0 :             DummyVec(1:m) = DummyVec(1:m) / sqrt(sumSqDummyVec)   ! DummyVec is now a uniform random point from inside of m-sphere
    2439           0 :             dummy = getRandBeta(0.5e0_RK*m,beta)
    2440           0 :             W(1:m) = sqrt(dummy) * DummyVec(1:m)
    2441           0 :             call getCholeskyFactor(m,RandCorMat(1:m,1:m),Diagonal(1:m))
    2442           0 :             if (Diagonal(1)<0._RK) then
    2443           0 :                 error stop
    2444             :                 !call abortProgram( output_unit , 1 , 1 , 'Statitistics@getRandCorMat()@getCholeskyFactor() failed.' )
    2445             :             end if
    2446           0 :             DummyVec(1:m) = 0._RK
    2447           0 :             do j = 1,m
    2448           0 :                 DummyVec(j) = DummyVec(j) + Diagonal(j) * W(j)
    2449           0 :                 do i = j+1,m
    2450           0 :                     DummyVec(i) = DummyVec(i) + RandCorMat(i,j) * DummyVec(j)
    2451             :                 end do
    2452             :             end do
    2453         300 :             RandCorMat(1:m,m+1) = DummyVec(1:m)
    2454             :         end do
    2455         600 :         do i=1,nd-1
    2456         900 :             RandCorMat(i+1:nd,i) = RandCorMat(i,i+1:nd)
    2457             :         end do
    2458         618 :   end function getRandCorMat
    2459             : 
    2460             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2461             : 
    2462             : !  function getRandCorMat(nd,eta)    ! based on the idea of LKJ (2007). But there is something wrong with this routine
    2463             : !    use Matrix_mod, only: getCholeskyFactor
    2464             : !    implicit none
    2465             : !    !integer, intent(in) :: nd,eta
    2466             : !    integer, intent(in) :: nd
    2467             : !    real(RK), intent(in) :: eta
    2468             : !    integer :: m,mNew,j,i
    2469             : !    real(RK) :: getRandCorMat(nd,nd), dummy, failureCounter
    2470             : !    real(RK) :: beta,sumSqDummyVec,DummyVec(nd-1),W(nd-1),Diagonal(nd-1)
    2471             : !
    2472             : !    if (nd<2 .or. eta<=0._RK) then  ! illegal value for eta. set getRandCorMat=0, return
    2473             : !      getRandCorMat = -1._RK
    2474             : !      return
    2475             : !    end if
    2476             : !
    2477             : !    do m = 1,nd
    2478             : !      getRandCorMat(m,m) = 1._RK
    2479             : !    end do
    2480             : !    beta = eta + 0.5_RK*(nd-2._RK)
    2481             : !
    2482             : !    do
    2483             : !      dummy = getRandBeta(beta,beta)
    2484             : !      if (dummy>0._RK .and. dummy<1._RK) exit
    2485             : !      write(*,*) "**Warning** random Beta variable out of bound.", dummy
    2486             : !      write(*,*) "Something is wrong with getRandBeta()."
    2487             : !      cycle
    2488             : !    end do
    2489             : !    getRandCorMat(1,2) = 2._RK * dummy - 1._RK    ! for the moment, only the upper half of getRandCorMat is needed, the lower half will contain cholesky lower triangle.
    2490             : !
    2491             : !    m = 2
    2492             : !    call getCholeskyFactor(m,getRandCorMat(1:m,1:m),Diagonal(1:m))
    2493             : !
    2494             : !    failureCounter = 0
    2495             : !    onionLayer: do
    2496             : !
    2497             : !      beta = beta - 0.5_RK
    2498             : !
    2499             : !      sumSqDummyVec = 0._RK
    2500             : !      do j=1,m
    2501             : !        DummyVec(j) = getRandGaus()
    2502             : !        sumSqDummyVec = sumSqDummyVec + DummyVec(j)**2
    2503             : !      end do
    2504             : !      DummyVec(1:m) = DummyVec(1:m) / sqrt(sumSqDummyVec)   ! DummyVec is now a uniform random point from inside of m-sphere
    2505             : !
    2506             : !      mNew = m + 1
    2507             : !      posDefCheck: do
    2508             : !
    2509             : !        do
    2510             : !          dummy = getRandBeta(0.5_RK*m,beta)
    2511             : !          if (dummy>0._RK .and. dummy<1._RK) exit
    2512             : !          write(*,*) "**Warning** random Beta variable out of bound.", dummy
    2513             : !          write(*,*) "Something is wrong with getRandBeta()."
    2514             : !          read(*,*)
    2515             : !          cycle
    2516             : !        end do
    2517             : !        W(1:m) = sqrt(dummy) * DummyVec(1:m)
    2518             : !
    2519             : !        getRandCorMat(1:m,mNew) = 0._RK
    2520             : !        do j = 1,m
    2521             : !          getRandCorMat(j,mNew) = getRandCorMat(j,mNew) + Diagonal(j) * W(j)
    2522             : !          do i = j+1,m
    2523             : !            getRandCorMat(i,mNew) = getRandCorMat(i,mNew) + getRandCorMat(i,j) * getRandCorMat(j,mNew)
    2524             : !          end do
    2525             : !        end do
    2526             : !
    2527             : !
    2528             : !        call getCholeskyFactor(mNew,getRandCorMat(1:mNew,1:mNew),Diagonal(1:mNew))  ! Now check if the new matrix is positive-definite, then proceed with the next layer
    2529             : !        if (Diagonal(1)<0._RK) then
    2530             : !          failureCounter = failureCounter + 1
    2531             : !          cycle posDefCheck
    2532             : !          !write(*,*) "Cholesky factorization failed in getRandCorMat()."
    2533             : !          !write(*,*) m
    2534             : !          !write(*,*) getRandCorMat(1:m,1:m)
    2535             : !          !stop
    2536             : !        end if
    2537             : !        exit posDefCheck
    2538             : !
    2539             : !      end do posDefCheck
    2540             : !
    2541             : !      if (mNew==nd) exit onionLayer
    2542             : !      m = mNew
    2543             : !
    2544             : !    end do onionLayer
    2545             : !
    2546             : !    if (failureCounter>0) write(*,*) 'failureRatio: ', dble(failureCounter)/dble(nd-2)
    2547             : !    do i=1,nd-1
    2548             : !      getRandCorMat(i+1:nd,i) = getRandCorMat(i,i+1:nd)
    2549             : !    end do
    2550             : !
    2551             : !  end function getRandCorMat
    2552             : 
    2553             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2554             : 
    2555             :     !> \brief
    2556             :     !> Returns a random correlation matrix using Monte Carlo rejection method.
    2557             :     !>
    2558             :     !> \param[in]   nd      :   The rank of the correlation matrix.
    2559             :     !> \param[in]   minRho  :   The minimum correlation coefficient to be expected in the output random correlation matrix.
    2560             :     !> \param[in]   maxRho  :   The maximum correlation coefficient to be expected in the output random correlation matrix.
    2561             :     !>
    2562             :     !> \return
    2563             :     !> `RandCorMat` : A random correlation matrix. Only the upper half of
    2564             :     !> `RandCorMat` is the correlation matrix, lower half is NOT set on output.
    2565             :     !>
    2566             :     !> \warning
    2567             :     !> The conditions `nd >= 1` and `maxRho < minRho` must hold, otherwise, `RandCorMat(1,1) = -1._RK` will be returned.
    2568             :     !>
    2569             :     !> \remark
    2570             :     !> This subroutine is very slow for high matrix dimensions ( `nd >~ 10` ).
    2571             :     !>
    2572             :     !> \author
    2573             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2574        2754 :     function getRandCorMatRejection(nd,minRho,maxRho) result(RandCorMat)
    2575             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2576             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRandCorMatRejection
    2577             : #endif
    2578             : 
    2579         309 :         use Matrix_mod, only: isPosDef
    2580             :         implicit none
    2581             :         integer(IK), intent(in) :: nd
    2582             :         real(RK)   , intent(in) :: minRho,maxRho
    2583         918 :         real(RK)                :: RandCorMat(nd,nd), RhoVec(nd*(nd-1))
    2584             :         integer(IK)             :: i,j,irho
    2585         306 :         if (maxRho<minRho .or. nd<1_IK) then
    2586           6 :             RandCorMat(1,1) = -1._RK
    2587           6 :             return
    2588             :         end if
    2589         300 :         if (nd==1_IK) then
    2590           0 :           RandCorMat = 1._RK
    2591             :         else
    2592           0 :             rejection: do
    2593         300 :                 call random_number(RhoVec)
    2594         900 :                 RhoVec = minRho + RhoVec*(maxRho-minRho)
    2595         300 :                 irho = 0
    2596         900 :                 do j=1,nd
    2597         600 :                     RandCorMat(j,j) = 1._RK
    2598        1200 :                     do i=1,j-1
    2599         300 :                         irho = irho + 1
    2600         900 :                         RandCorMat(i,j) = RhoVec(irho)
    2601             :                     end do
    2602             :                 end do
    2603         300 :                 if (isPosDef(nd,RandCorMat)) exit rejection
    2604           0 :                 cycle rejection
    2605             :             end do rejection
    2606             :         end if
    2607         600 :         do j=1,nd-1
    2608         900 :             RandCorMat(j+1:nd,j) = RandCorMat(j,j+1:nd)
    2609             :         end do
    2610         612 :   end function getRandCorMatRejection
    2611             : 
    2612             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2613             : 
    2614             :     !> \brief
    2615             :     !> Convert the upper-triangle covariance matrix to the upper-triangle correlation matrix.
    2616             :     !>
    2617             :     !> \param[in]   nd          :   The rank of the covariance matrix.
    2618             :     !> \param[in]   CovMatUpper :   The upper-triangle covariance matrix. The lower-triangle will not be used.
    2619             :     !>
    2620             :     !> \return
    2621             :     !> `CorMatUpper` : An upper-triangle correlation matrix. The lower-triangle will NOT be set.
    2622             :     !>
    2623             :     !> \author
    2624             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2625        3840 :     pure function getCorMatUpperFromCovMatUpper(nd,CovMatUpper) result(CorMatUpper)
    2626             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2627             :         !DEC$ ATTRIBUTES DLLEXPORT :: getCorMatUpperFromCovMatUpper
    2628             : #endif
    2629             :         implicit none
    2630             :         integer(IK)  , intent(in) :: nd
    2631             :         real(RK)     , intent(in) :: CovMatUpper(nd,nd)
    2632             :         real(RK)                  :: CorMatUpper(nd,nd)
    2633        1350 :         real(RK)                  :: InverseStdVec(nd)
    2634             :         integer(IK)               :: i,j
    2635        1350 :         do j = 1, nd
    2636         830 :             InverseStdVec(j) = 1._RK / sqrt(CovMatUpper(j,j))
    2637        2490 :             do i = 1, j
    2638        1970 :                 CorMatUpper(i,j) = CovMatUpper(i,j) * InverseStdVec(j) * InverseStdVec(i)
    2639             :             end do
    2640             :         end do
    2641         826 :     end function getCorMatUpperFromCovMatUpper
    2642             : 
    2643             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2644             : 
    2645             :     !> \brief
    2646             :     !> Convert the upper-triangle correlation matrix to the upper-triangle covariance matrix.
    2647             :     !>
    2648             :     !> \param[in]   nd          :   The rank of the correlation matrix.
    2649             :     !> \param[in]   StdVec      :   The input standard deviation vector of length `nd`.
    2650             :     !> \param[in]   CorMatUpper :   The upper-triangle correlation matrix. The lower-triangle will not be used.
    2651             :     !>
    2652             :     !> \return
    2653             :     !> `CovMatUpper` : An upper-triangle covariance matrix. The lower-triangle will NOT be set.
    2654             :     !>
    2655             :     !> \author
    2656             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2657          27 :     pure function getCovMatUpperFromCorMatUpper(nd,StdVec,CorMatUpper) result(CovMatUpper)
    2658             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2659             :         !DEC$ ATTRIBUTES DLLEXPORT :: getCovMatUpperFromCorMatUpper
    2660             : #endif
    2661             :         implicit none
    2662             :         integer(IK)  , intent(in) :: nd
    2663             :         real(RK)     , intent(in) :: StdVec(nd), CorMatUpper(nd,nd)
    2664             :         real(RK)                  :: CovMatUpper(nd,nd)
    2665             :         integer(IK)               :: i,j
    2666           9 :         do j=1,nd
    2667          18 :             do i=1,j
    2668          15 :                 CovMatUpper(i,j) = CorMatUpper(i,j) * StdVec(j) * StdVec(i)
    2669             :             end do
    2670             :         end do
    2671         523 :     end function getCovMatUpperFromCorMatUpper
    2672             : 
    2673             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2674             : 
    2675             :     !> \brief
    2676             :     !> Convert the lower-triangle correlation matrix to the upper-triangle covariance matrix.
    2677             :     !>
    2678             :     !> \param[in]   nd          :   The rank of the correlation matrix.
    2679             :     !> \param[in]   StdVec      :   The input standard deviation vector of length `nd`.
    2680             :     !> \param[in]   CorMatLower :   The lower-triangle correlation matrix. The upper-triangle will not be used.
    2681             :     !>
    2682             :     !> \return
    2683             :     !> `CovMatUpper` : An upper-triangle covariance matrix. The lower-triangle will NOT be set.
    2684             :     !>
    2685             :     !> \author
    2686             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2687          27 :     pure function getCovMatUpperFromCorMatLower(nd,StdVec,CorMatLower) result(CovMatUpper)
    2688             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2689             :         !DEC$ ATTRIBUTES DLLEXPORT :: getCovMatUpperFromCorMatLower
    2690             : #endif
    2691             :         implicit none
    2692             :         integer(IK)  , intent(in) :: nd
    2693             :         real(RK)     , intent(in) :: StdVec(nd), CorMatLower(nd,nd)
    2694             :         real(RK)                  :: CovMatUpper(nd,nd)
    2695             :         integer(IK)               :: i,j
    2696           9 :         do j=1,nd
    2697           6 :             CovMatUpper(j,j) = StdVec(j)**2
    2698          12 :             do i=1,j-1
    2699           9 :                 CovMatUpper(i,j) = CorMatLower(j,i) * StdVec(j) * StdVec(i)
    2700             :             end do
    2701             :         end do
    2702           6 :     end function getCovMatUpperFromCorMatLower
    2703             : 
    2704             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2705             : 
    2706             :     !> \brief
    2707             :     !> Convert the upper-triangle correlation matrix to the lower-triangle covariance matrix.
    2708             :     !>
    2709             :     !> \param[in]   nd          :   The rank of the correlation matrix.
    2710             :     !> \param[in]   StdVec      :   The input standard deviation vector of length `nd`.
    2711             :     !> \param[in]   CorMatUpper :   The upper-triangle correlation matrix. The lower-triangle will not be used.
    2712             :     !>
    2713             :     !> \return
    2714             :     !> `CovMatLower` : An lower-triangle covariance matrix. The upper-triangle will NOT be set.
    2715             :     !>
    2716             :     !> \author
    2717             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2718          27 :     pure function getCovMatLowerFromCorMatUpper(nd,StdVec,CorMatUpper) result(CovMatLower)
    2719             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2720             :         !DEC$ ATTRIBUTES DLLEXPORT :: getCovMatLowerFromCorMatUpper
    2721             : #endif
    2722             :         implicit none
    2723             :         integer(IK)  , intent(in) :: nd
    2724             :         real(RK)     , intent(in) :: StdVec(nd), CorMatUpper(nd,nd)
    2725             :         real(RK)                  :: CovMatLower(nd,nd)
    2726             :         integer(IK)               :: i,j
    2727           9 :         do j=1,nd
    2728           6 :             CovMatLower(j,j) = StdVec(j)**2
    2729          12 :             do i=1,j-1
    2730           9 :                 CovMatLower(j,i) = CorMatUpper(i,j) * StdVec(j) * StdVec(i)
    2731             :             end do
    2732             :         end do
    2733           6 :     end function getCovMatLowerFromCorMatUpper
    2734             : 
    2735             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2736             : 
    2737             :     !> \brief
    2738             :     !> Convert the lower-triangle correlation matrix to the lower-triangle covariance matrix.
    2739             :     !>
    2740             :     !> \param[in]   nd          :   The rank of the correlation matrix.
    2741             :     !> \param[in]   StdVec      :   The input standard deviation vector of length `nd`.
    2742             :     !> \param[in]   CorMatLower :   The lower-triangle correlation matrix. The upper-triangle will not be used.
    2743             :     !>
    2744             :     !> \return
    2745             :     !> `CovMatLower` : An lower-triangle covariance matrix. The upper-triangle will NOT be set.
    2746             :     !>
    2747             :     !> \author
    2748             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2749          27 :     pure function getCovMatLowerFromCorMatLower(nd,StdVec,CorMatLower) result(CovMatLower)
    2750             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2751             :         !DEC$ ATTRIBUTES DLLEXPORT :: getCovMatLowerFromCorMatLower
    2752             : #endif
    2753             :         implicit none
    2754             :         integer(IK)  , intent(in) :: nd
    2755             :         real(RK)     , intent(in) :: StdVec(nd), CorMatLower(nd,nd)
    2756             :         real(RK)                  :: CovMatLower(nd,nd)
    2757             :         integer(IK)               :: i,j
    2758           9 :         do j=1,nd
    2759           6 :             CovMatLower(j,j) = StdVec(j)**2
    2760          12 :             do i=1,j-1
    2761           9 :                 CovMatLower(j,i) = CorMatLower(j,i) * StdVec(j) * StdVec(i)
    2762             :             end do
    2763             :         end do
    2764           6 :     end function getCovMatLowerFromCorMatLower
    2765             : 
    2766             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2767             : 
    2768             :     !> \brief
    2769             :     !> Convert the input correlation matrix to the output covariance matrix.
    2770             :     !>
    2771             :     !> \param[in]   nd          :   The rank of the correlation matrix.
    2772             :     !> \param[in]   StdVec      :   The input standard deviation vector of length `nd`.
    2773             :     !> \param[in]   CorMatUpper :   The upper-triangle correlation matrix. The lower-triangle will not be used.
    2774             :     !>
    2775             :     !> \return
    2776             :     !> `CovMatFull` : The full covariance matrix.
    2777             :     !>
    2778             :     !> \author
    2779             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2780        7098 :     pure function getCovMatFromCorMatUpper(nd,StdVec,CorMatUpper) result(CovMatFull)
    2781             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2782             :         !DEC$ ATTRIBUTES DLLEXPORT :: getCovMatFromCorMatUpper
    2783             : #endif
    2784             :         implicit none
    2785             :         integer(IK)  , intent(in) :: nd
    2786             :         real(RK)     , intent(in) :: StdVec(nd), CorMatUpper(nd,nd)   ! only upper half needed
    2787             :         real(RK)                  :: CovMatFull(nd,nd)
    2788             :         integer(IK)               :: i,j
    2789        2562 :         do j=1,nd
    2790        1512 :             CovMatFull(j,j) = StdVec(j)**2
    2791        3024 :             do i=1,j-1
    2792         462 :                 CovMatFull(i,j) = CorMatUpper(i,j) * StdVec(j) * StdVec(i)
    2793        1974 :                 CovMatFull(j,i) = CovMatFull(i,j)
    2794             :             end do
    2795             :         end do
    2796        1053 :     end function getCovMatFromCorMatUpper
    2797             : 
    2798             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2799             : 
    2800             : !    !> \brief
    2801             : !    !> Return the Geometric distribution PDF values for a range of trials, starting at index `1`.
    2802             : !    !> If the probability of success on each trial is `successProb`, then the probability that
    2803             : !    !> the `k`th trial (out of `k` trials) is the first success is `GeoLogPDF(k)`.
    2804             : !    !>
    2805             : !    !> \param[in]   successProb     :   The probability of success.
    2806             : !    !> \param[in]   logPdfPrecision :   The precision value below which the PDF is practically considered to be zero (**optional**).
    2807             : !    !> \param[in]   minSeqLen       :   The minimum length of the range of `k` values for which the PDF will be computed (**optional**).
    2808             : !    !> \param[in]   seqLen          :   The length of the range of `k` values for which the PDF will be computed (**optional**).
    2809             : !    !>                                  If provided, it will overwrite the the output sequence length as inferred from
    2810             : !    !>                                  the combination of `minSeqLen` and `logPdfPrecision`.
    2811             : !    !>
    2812             : !    !> \return
    2813             : !    !> `GeoLogPDF`  :   An allocatable representing the geometric PDF over a range of `k` values, whose length is
    2814             : !    !>                  `seqLen`, or if not provided, is determined from the values of `logPdfPrecision` and `minSeqLen`.
    2815             : !    !>
    2816             : !    !> \author
    2817             : !    !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2818             : !    function getGeoLogPDF_old(successProb,logPdfPrecision,minSeqLen,seqLen) result(GeoLogPDF)
    2819             : !#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2820             : !        !DEC$ ATTRIBUTES DLLEXPORT :: getGeoLogPDF_old
    2821             : !#endif
    2822             : !        use Constants_mod, only: IK, RK
    2823             : !        implicit none
    2824             : !        real(RK)    , intent(in)            :: successProb
    2825             : !        real(RK)    , intent(in), optional  :: logPdfPrecision
    2826             : !        integer(IK) , intent(in), optional  :: minSeqLen
    2827             : !        integer(IK) , intent(in), optional  :: seqLen
    2828             : !        real(RK)    , allocatable           :: GeoLogPDF(:)
    2829             : !        real(RK)    , parameter             :: LOG_PDF_PRECISION = log(0.001_RK)
    2830             : !        real(RK)                            :: logProbFailure
    2831             : !        integer(IK)                         :: lenGeoLogPDF, i
    2832             : !        logProbFailure = log(1._RK - successProb)
    2833             : !        if (present(seqLen)) then
    2834             : !            lenGeoLogPDF = seqLen
    2835             : !        else
    2836             : !            if (present(logPdfPrecision)) then
    2837             : !                lenGeoLogPDF = ceiling(  logPdfPrecision / logProbFailure)
    2838             : !            else
    2839             : !                lenGeoLogPDF = ceiling(LOG_PDF_PRECISION / logProbFailure)
    2840             : !            end if
    2841             : !            if (present(minSeqLen)) lenGeoLogPDF = max(minSeqLen,lenGeoLogPDF)
    2842             : !        end if
    2843             : !        allocate(GeoLogPDF(lenGeoLogPDF))
    2844             : !        GeoLogPDF(1) = log(successProb)
    2845             : !        do i = 2, lenGeoLogPDF
    2846             : !            GeoLogPDF(i) = GeoLogPDF(i-1) + logProbFailure
    2847             : !        end do
    2848             : !    end function getGeoLogPDF_old
    2849             : 
    2850             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2851             : 
    2852             :     !> \brief
    2853             :     !> Return the Geometric distribution PDF values for the input number of trials,
    2854             :     !> the trials at which first success happens, and the success probability.
    2855             :     !>
    2856             :     !> \param[in]   numTrial    :   The number of trials.
    2857             :     !> \param[in]   SuccessStep :   The vector of trials of length `numTrial` at which the first success happens.
    2858             :     !> \param[in]   successProb :   The probability of success.
    2859             :     !>
    2860             :     !> \return
    2861             :     !> `LogProbGeo` :   An output vector of length `numTrial` representing the geometric PDF values.
    2862             :     !>
    2863             :     !> \author
    2864             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2865          39 :     pure function getLogProbGeo(numTrial, SuccessStep, successProb) result(LogProbGeo)
    2866             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2867             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbGeo
    2868             : #endif
    2869        1050 :         use Constants_mod, only: IK, RK
    2870             :         implicit none
    2871             :         integer(IK) , intent(in)    :: numTrial
    2872             :         integer(IK) , intent(in)    :: SuccessStep(numTrial)
    2873             :         real(RK)    , intent(in)    :: successProb
    2874             :         real(RK)                    :: LogProbGeo(numTrial)
    2875           3 :         real(RK)                    :: logProbSuccess, logProbFailure
    2876           3 :         logProbSuccess = log(successProb)
    2877           3 :         logProbFailure = log(1._RK - successProb)
    2878          33 :         LogProbGeo = logProbSuccess + (SuccessStep-1_IK) * logProbFailure
    2879           3 :     end function getLogProbGeo
    2880             : 
    2881             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2882             : 
    2883             :     !> \brief
    2884             :     !> Compute the natural logarithm of the Geometric distribution PDF of a limited range of Bernoulli trials,
    2885             :     !> starting at index `1` up to `maxNumTrial`. In other words, upon reaching the trial `maxNumTrial`,
    2886             :     !> the Bernoulli trials count restart from index `1`. This Cyclic Geometric distribution is
    2887             :     !> particularly useful in the parallelization studies of Monte Carlo simulation.
    2888             :     !>
    2889             :     !> \param[in]   successProb :   The probability of success.
    2890             :     !> \param[in]   maxNumTrial :   The maximum number of trails possible.
    2891             :     !>                              After `maxNumTrial` tries, the Geometric distribution restarts from index `1`.
    2892             :     !> \param[in]   numTrial    :   The length of the array `SuccessStep`. Note that `numTrial < maxNumTrial` must hold.
    2893             :     !> \param[in]   SuccessStep :   A vector of length `(1:numTrial)` of integers that represent
    2894             :     !>                              the steps at which the Bernoulli successes occur.
    2895             :     !>
    2896             :     !> \return
    2897             :     !> `LogProbGeoCyclic`       :   A real-valued vector of length `(1:numTrial)` whose values represent the probabilities
    2898             :     !>                              of having Bernoulli successes at the corresponding SuccessStep values.
    2899             :     !>
    2900             :     !> \warning
    2901             :     !> Any value of SuccessStep must be an integer numbers between `1` and `maxNumTrial`.
    2902             :     !> The onus is on the user to ensure this condition holds.
    2903             :     !>
    2904             :     !> \author
    2905             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2906      575452 :     pure function getLogProbGeoCyclic(successProb,maxNumTrial,numTrial,SuccessStep) result(LogProbGeoCyclic)
    2907             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2908             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbGeoCyclic
    2909             : #endif
    2910           3 :         use Constants_mod, only: IK, RK, NEGLOGINF_RK
    2911             :         implicit none
    2912             :         real(RK)    , intent(in)    :: successProb
    2913             :         integer(IK) , intent(in)    :: maxNumTrial
    2914             :         integer(IK) , intent(in)    :: numTrial
    2915             :         integer(IK) , intent(in)    :: SuccessStep(numTrial)
    2916             :         real(RK)                    :: LogProbGeoCyclic(numTrial)
    2917       73608 :         real(RK)                    :: failureProb, logProbSuccess, logProbFailure, logDenominator, exponentiation
    2918       73608 :         if (successProb>0._RK .and. successProb<1._RK) then ! tolerate log(0)
    2919       73608 :             failureProb = 1._RK - successProb
    2920       73608 :             logProbSuccess = log(successProb)
    2921       73608 :             logProbFailure = log(failureProb)
    2922       73608 :             exponentiation = maxNumTrial * logProbFailure
    2923       73608 :             if (exponentiation<NEGLOGINF_RK) then ! tolerate underflow
    2924           0 :                 logDenominator = 0._RK
    2925             :             else
    2926       73608 :                 exponentiation = exp(exponentiation)
    2927       73608 :                 if (exponentiation<1._RK) then ! tolerate log(0)
    2928       73608 :                     logDenominator = log(1._RK-exponentiation)
    2929             :                 else
    2930           0 :                     logDenominator = NEGLOGINF_RK
    2931             :                 end if
    2932             :             end if
    2933      428236 :             LogProbGeoCyclic = logProbSuccess + (SuccessStep-1) * logProbFailure - logDenominator
    2934           0 :         elseif (successProb==0._RK) then
    2935           0 :             LogProbGeoCyclic = -log(real(maxNumTrial,kind=RK))
    2936           0 :         elseif (successProb==1._RK) then
    2937           0 :             LogProbGeoCyclic(1) = 0._RK
    2938           0 :             LogProbGeoCyclic(2:numTrial) = NEGLOGINF_RK
    2939             :         else
    2940           0 :             LogProbGeoCyclic = NEGLOGINF_RK
    2941             :         end if
    2942       73608 :     end function getLogProbGeoCyclic
    2943             : 
    2944             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2945             : 
    2946             :     !> \brief
    2947             :     !> Return the standard normal distribution PDF value.
    2948             :     !>
    2949             :     !> \param[in]   z   :   The input value at which the PDF will be computed.
    2950             :     !>
    2951             :     !> \return
    2952             :     !> `snormPDF`   :   The standard normal distribution PDF value.
    2953             :     !>
    2954             :     !> \author
    2955             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2956           3 :     function getSNormPDF(z) result(snormPDF)
    2957             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2958             :         !DEC$ ATTRIBUTES DLLEXPORT :: getSNormPDF
    2959             : #endif
    2960       73608 :         use Constants_mod, only: INVSQRT2PI
    2961             :         implicit none
    2962             :         real(RK), intent(in) :: z
    2963             :         real(RK)             :: snormPDF
    2964           3 :         snormPDF = INVSQRT2PI * exp( -0.5_RK*z**2 )
    2965           6 :     end function getSNormPDF
    2966             : 
    2967             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2968             : 
    2969             :     !> \brief
    2970             :     !> Return the non-standard normal distribution PDF value.
    2971             :     !>
    2972             :     !> \param[in]   avg :   The mean of the Normal distribution.
    2973             :     !> \param[in]   std :   The standard deviation of the Normal distribution.
    2974             :     !> \param[in]   var :   The variance of the Normal distribution.
    2975             :     !> \param[in]   x   :   The point at which the PDF will be computed.
    2976             :     !>
    2977             :     !> \return
    2978             :     !> `normPDF`        :   The normal distribution PDF value at the given input point.
    2979             :     !>
    2980             :     !> \author
    2981             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    2982           3 :     function getNormPDF(avg,std,var,x) result(normPDF)
    2983             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    2984             :         !DEC$ ATTRIBUTES DLLEXPORT :: getNormPDF
    2985             : #endif
    2986           3 :         use Constants_mod, only: INVSQRT2PI
    2987             :         implicit none
    2988             :         real(RK), intent(in) :: avg,std,var,x
    2989             :         real(RK)             :: normPDF
    2990           3 :         normPDF = INVSQRT2PI * exp( -(x-avg)**2/(2._RK*var) ) / std
    2991           6 :     end function getNormPDF
    2992             : 
    2993             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2994             : 
    2995             :     !> \brief
    2996             :     !> Return the non-standard normal distribution Cumulative Probability Density function (CDF) value.
    2997             :     !>
    2998             :     !> \param[in]   avg :   The mean of the Normal distribution.
    2999             :     !> \param[in]   std :   The standard deviation of the Normal distribution.
    3000             :     !> \param[in]   x   :   The point at which the PDF will be computed.
    3001             :     !>
    3002             :     !> \return
    3003             :     !> `normCDF`        :   The normal distribution CDF value at the given input point.
    3004             :     !>
    3005             :     !> \author
    3006             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    3007           3 :     pure function getNormCDF(avg,std,x) result(normCDF)
    3008             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3009             :         !DEC$ ATTRIBUTES DLLEXPORT :: getNormCDF
    3010             : #endif
    3011           3 :         use Constants_mod, only: RK, SQRT2
    3012             :         implicit none
    3013             :         real(RK), intent(in) :: avg,std,x
    3014             :         real(RK)             :: normCDF
    3015           3 :         normCDF = 0.5_RK * ( 1._RK + erf( real((x-avg)/(SQRT2*std),kind=RK) ) )
    3016           6 :     end function getNormCDF
    3017             : 
    3018             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3019             : 
    3020             :     !> \brief
    3021             :     !> Return the standard normal distribution Cumulative Probability Density function (CDF) value.
    3022             :     !>
    3023             :     !> \param[in]   x       :   The point at which the PDF will be computed.
    3024             :     !>
    3025             :     !> \return
    3026             :     !> `normCDF`   :   The normal distribution CDF value at the given input point.
    3027             :     !>
    3028             :     !> \remark
    3029             :     !> This procedure performs all calculations in `real32` real kind. If 64-bit accuracy matters more than performance,
    3030             :     !> then use the [getSNormCDF_DPR](@ref getsnormcdf_dpr) for a more-accurate double-precision but slower results.
    3031             :     !>
    3032             :     !> \author
    3033             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    3034           3 :     pure function getSNormCDF_SPR(x) result(normCDF)
    3035             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3036             :         !DEC$ ATTRIBUTES DLLEXPORT :: getSNormCDF_SPR
    3037             : #endif
    3038             :         use iso_fortran_env, only: RK => real32
    3039             :         implicit none
    3040             :         real(RK)    , intent(in)    :: x
    3041             :         real(RK)                    :: normCDF
    3042             :         real(RK), parameter         :: INVSQRT2 = 1._RK / sqrt(2._RK)
    3043           3 :         normCDF = 0.5_RK * ( 1._RK + erf( real(x*INVSQRT2,kind=RK) ) )
    3044           3 :     end function getSNormCDF_SPR
    3045             : 
    3046             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3047             : 
    3048             :     !> \brief
    3049             :     !> Return the standard normal distribution Cumulative Probability Density function (CDF) value.
    3050             :     !>
    3051             :     !> \param[in]   x   :   The point at which the PDF will be computed.
    3052             :     !>
    3053             :     !> \return
    3054             :     !> `normCDF`   :   The normal distribution CDF value at the given input point.
    3055             :     !>
    3056             :     !> \remark
    3057             :     !> This procedure performs all calculations in `DPR` (`real64`) real kind. If performance matters more than 64-bit accuracy,
    3058             :     !> then use the [getSNormCDF_SPR](@ref getsnormcdf_spr) for a faster, but less-accurate single-precision results.
    3059             :     !>
    3060             :     !> \author
    3061             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    3062         153 :     pure function getSNormCDF_DPR(x) result(normCDF)
    3063             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3064             :         !DEC$ ATTRIBUTES DLLEXPORT :: getSNormCDF_DPR
    3065             : #endif
    3066             :         use iso_fortran_env, only: RK => real64
    3067             :         implicit none
    3068             :         real(RK), intent(in)    :: x
    3069             :         real(RK)                :: normCDF
    3070             :         real(RK), parameter     :: INVSQRT2 = 1._RK / sqrt(2._RK)
    3071         153 :         normCDF = 0.5_RK * ( 1._RK + erf( real(x*INVSQRT2,kind=RK) ) )
    3072           3 :     end function getSNormCDF_DPR
    3073             : 
    3074             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3075             : 
    3076             :     !> \brief
    3077             :     !> Return the Beta distribution Cumulative Probability Density function (CDF) value.
    3078             :     !>
    3079             :     !> \param[in]   alpha   :   The first shape parameter of the Beta distribution.
    3080             :     !> \param[in]   beta    :   The second shape parameter of the Beta distribution.
    3081             :     !> \param[in]   x       :   The point at which the CDF will be computed.
    3082             :     !>
    3083             :     !> \return
    3084             :     !> `betaCDF`   :   The Beta distribution CDF value at the given input point.
    3085             :     !>
    3086             :     !> \warning
    3087             :     !> If `x` is not in the range `[0,1]`, a negative value for `betaCDF` will be returned to indicate the occurrence of error.
    3088             :     !>
    3089             :     !> \warning
    3090             :     !> The onus is on the user to ensure that the input (`alpha`, `beta`) shape parameters are positive.
    3091             :     !>
    3092             :     !> \remark
    3093             :     !> This procedure performs all calculations in `real32` real kind. If 64-bit accuracy matters more than performance,
    3094             :     !> then use the [getBetaCDF_DPR](@ref getbetacdf_dpr) for a more-accurate double-precision but slower results.
    3095             :     !>
    3096             :     !> \todo
    3097             :     !> The efficiency of this code can be improved by making `x` a vector on input.
    3098             :     !>
    3099             :     !> \author
    3100             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    3101           9 :     function getBetaCDF_SPR(alpha,beta,x) result(betaCDF)
    3102             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3103             :         !DEC$ ATTRIBUTES DLLEXPORT :: getBetaCDF_SPR
    3104             : #endif
    3105             :         use iso_fortran_env, only: RK => real32
    3106             :         implicit none
    3107             :         real(RK), intent(in)    :: alpha, beta, x
    3108           9 :         real(RK)                :: bt
    3109             :         real(RK)                :: betaCDF
    3110           9 :         if (x < 0._RK .or. x > 1._RK) then
    3111           6 :             betaCDF = -1._RK
    3112           6 :             return
    3113             :         end if
    3114           3 :         if (x == 0._RK .or. x == 1._RK) then
    3115           0 :             bt = 0._RK
    3116             :         else
    3117             :             bt = exp( log_gamma(real(alpha+beta,kind=RK)) &
    3118             :                     - log_gamma(real(alpha,kind=RK)) - log_gamma(real(beta,kind=RK)) &
    3119           3 :                     + alpha*log(x) + beta*log(1._RK-x) )
    3120             :         end if
    3121           3 :         if ( x < (alpha+1._RK) / (alpha+beta+2._RK) ) then
    3122           0 :             betaCDF = bt * getBetaContinuedFraction(alpha,beta,x) / alpha
    3123             :         else
    3124           3 :             betaCDF = 1._RK - bt * getBetaContinuedFraction(beta,alpha,1._RK-x) / beta
    3125             :         end if
    3126         162 :     end function getBetaCDF_SPR
    3127             : 
    3128             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3129             : 
    3130             :     !> \brief
    3131             :     !> Return the Beta distribution Cumulative Probability Density function (CDF) value.
    3132             :     !>
    3133             :     !> \param[in]   alpha   :   The first shape parameter of the Beta distribution.
    3134             :     !> \param[in]   beta    :   The second shape parameter of the Beta distribution.
    3135             :     !> \param[in]   x       :   The point at which the CDF will be computed.
    3136             :     !>
    3137             :     !> \return
    3138             :     !> `betaCDF`   :   The Beta distribution CDF value at the given input point.
    3139             :     !>
    3140             :     !> \warning
    3141             :     !> If `x` is not in the range `[0,1]`, a negative value for `betaCDF` will be returned to indicate the occurrence of error.
    3142             :     !>
    3143             :     !> \warning
    3144             :     !> The onus is on the user to ensure that the input (`alpha`, `beta`) shape parameters are positive.
    3145             :     !>
    3146             :     !> \remark
    3147             :     !> This procedure performs all calculations in `DPR` (`real64`) real kind. If performance matters more than 64-bit accuracy,
    3148             :     !> then use the [getBetaCDF_SPR](@ref getbetacdf_spr) for a faster, but less-accurate single-precision results.
    3149             :     !>
    3150             :     !> \todo
    3151             :     !> The efficiency of this code can be improved by making `x` a vector on input.
    3152             :     !>
    3153             :     !> \author
    3154             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    3155          15 :     function getBetaCDF_DPR(alpha,beta,x) result(betaCDF)
    3156             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3157             :         !DEC$ ATTRIBUTES DLLEXPORT :: getBetaCDF_DPR
    3158             : #endif
    3159             :         use iso_fortran_env, only: RK => real64
    3160             :         implicit none
    3161             :         real(RK), intent(in)    :: alpha, beta, x
    3162          15 :         real(RK)                :: bt
    3163             :         real(RK)                :: betaCDF
    3164          15 :         if (x < 0._RK .or. x > 1._RK) then
    3165           6 :             betaCDF = -1._RK
    3166           6 :             return
    3167             :         end if
    3168           9 :         if (x == 0._RK .or. x == 1._RK) then
    3169           0 :             bt = 0._RK
    3170             :         else
    3171             :             bt = exp( log_gamma(real(alpha+beta,kind=RK)) &
    3172             :                     - log_gamma(real(alpha,kind=RK)) - log_gamma(real(beta,kind=RK)) &
    3173           9 :                     + alpha*log(x) + beta*log(1._RK-x) )
    3174             :         end if
    3175           9 :         if ( x < (alpha+1._RK) / (alpha+beta+2._RK) ) then
    3176           6 :             betaCDF = bt * getBetaContinuedFraction(alpha,beta,x) / alpha
    3177             :         else
    3178           3 :             betaCDF = 1._RK - bt * getBetaContinuedFraction(beta,alpha,1._RK-x) / beta
    3179             :         end if
    3180          24 :     end function getBetaCDF_DPR
    3181             : 
    3182             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3183             : 
    3184             :     !> \brief
    3185             :     !> Return the Beta Continued Fraction (BCF).
    3186             :     !>
    3187             :     !> \param[in]   alpha   :   The first shape parameter of the Beta distribution.
    3188             :     !> \param[in]   beta    :   The second shape parameter of the Beta distribution.
    3189             :     !> \param[in]   x       :   The point at which the BCF will be computed.
    3190             :     !>
    3191             :     !> \return
    3192             :     !> `betaCDF`   :   The BCF value at the given input point.
    3193             :     !>
    3194             :     !> \author
    3195             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    3196           3 :     function getBetaContinuedFraction_SPR(alpha,beta,x) result(betaContinuedFraction)
    3197             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3198             :         !DEC$ ATTRIBUTES DLLEXPORT :: getBetaContinuedFraction_SPR
    3199             : #endif
    3200             :         use iso_fortran_env, only: RK => real32
    3201             :         implicit none
    3202             :         real(RK)   , intent(in) :: alpha,beta,x
    3203             :         real(RK)   , parameter  :: eps = epsilon(x), fpmin = tiny(x)/eps
    3204             :         integer(IK), parameter  :: maxit = 100_IK
    3205           3 :         real(RK)                :: aa,c,d,del,qab,qam,qap
    3206             :         real(RK)                :: betaContinuedFraction
    3207             :         integer(IK)             :: m,m2
    3208           3 :         qab = alpha+beta
    3209           3 :         qap = alpha+1._RK
    3210           3 :         qam = alpha-1._RK
    3211           3 :         c = 1._RK
    3212           3 :         d = 1._RK-qab*x/qap
    3213           0 :         if (abs(d) < fpmin) d = fpmin
    3214           3 :         d = 1._RK/d
    3215           3 :         betaContinuedFraction = d
    3216           6 :         do m = 1,maxit
    3217           6 :             m2 = 2*m
    3218           6 :             aa = m*(beta-m)*x/((qam+m2)*(alpha+m2))
    3219           6 :             d = 1._RK+aa*d
    3220           6 :             if (abs(d) < fpmin) d = fpmin
    3221           6 :             c = 1._RK+aa/c
    3222           6 :             if (abs(c) < fpmin) c = fpmin
    3223           6 :             d = 1._RK/d
    3224           6 :             betaContinuedFraction = betaContinuedFraction*d*c
    3225           6 :             aa = -(alpha+m)*(qab+m)*x/((alpha+m2)*(qap+m2))
    3226           6 :             d = 1._RK+aa*d
    3227           6 :             if (abs(d) < fpmin) d = fpmin
    3228           6 :             c = 1._RK+aa/c
    3229           6 :             if (abs(c) < fpmin) c = fpmin
    3230           6 :             d = 1._RK/d
    3231           6 :             del = d*c
    3232           6 :             betaContinuedFraction = betaContinuedFraction*del
    3233           6 :             if (abs(del-1._RK) <=  eps) exit
    3234             :         end do
    3235           3 :         if (m > maxit) then
    3236             :         ! LCOV_EXCL_START
    3237             :             error stop
    3238             :             !call abortProgram( output_unit , 1 , 1 , &
    3239             :             !'Statitistics@getBetaContinuedFraction_SPR() failed: alpha or beta too big, or maxit too small.' )
    3240             :         end if
    3241             :         ! LCOV_EXCL_STOP
    3242          18 :     end function getBetaContinuedFraction_SPR
    3243             : 
    3244             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3245             : 
    3246             :     !> \brief
    3247             :     !> Return the Beta Continued Fraction (BCF).
    3248             :     !>
    3249             :     !> \param[in]   alpha   :   The first shape parameter of the Beta distribution.
    3250             :     !> \param[in]   beta    :   The second shape parameter of the Beta distribution.
    3251             :     !> \param[in]   x       :   The point at which the BCF will be computed.
    3252             :     !>
    3253             :     !> \return
    3254             :     !> `betaCDF`   :   The BCF value at the given input point.
    3255             :     !>
    3256             :     !> \author
    3257             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    3258           9 :     function getBetaContinuedFraction_DPR(alpha,beta,x) result(betaContinuedFraction)
    3259             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3260             :         !DEC$ ATTRIBUTES DLLEXPORT :: getBetaContinuedFraction_DPR
    3261             : #endif
    3262             :         use iso_fortran_env, only: RK => real64
    3263             :         implicit none
    3264             :         real(RK)   , intent(in) :: alpha,beta,x
    3265             :         real(RK)   , parameter  :: eps = epsilon(x), fpmin = tiny(x)/eps
    3266             :         integer(IK), parameter  :: maxit = 100_IK
    3267           9 :         real(RK)                :: aa,c,d,del,qab,qam,qap
    3268             :         real(RK)                :: betaContinuedFraction
    3269             :         integer(IK)             :: m,m2
    3270           9 :         qab = alpha+beta
    3271           9 :         qap = alpha+1._RK
    3272           9 :         qam = alpha-1._RK
    3273           9 :         c = 1._RK
    3274           9 :         d = 1._RK-qab*x/qap
    3275           0 :         if (abs(d) < fpmin) d = fpmin
    3276           9 :         d = 1._RK/d
    3277           9 :         betaContinuedFraction = d
    3278          90 :         do m = 1,maxit
    3279          90 :             m2 = 2*m
    3280          90 :             aa = m*(beta-m)*x/((qam+m2)*(alpha+m2))
    3281          90 :             d = 1._RK+aa*d
    3282          90 :             if (abs(d) < fpmin) d = fpmin
    3283          90 :             c = 1._RK+aa/c
    3284          90 :             if (abs(c) < fpmin) c = fpmin
    3285          90 :             d = 1._RK/d
    3286          90 :             betaContinuedFraction = betaContinuedFraction*d*c
    3287          90 :             aa = -(alpha+m)*(qab+m)*x/((alpha+m2)*(qap+m2))
    3288          90 :             d = 1._RK+aa*d
    3289          90 :             if (abs(d) < fpmin) d = fpmin
    3290          90 :             c = 1._RK+aa/c
    3291          90 :             if (abs(c) < fpmin) c = fpmin
    3292          90 :             d = 1._RK/d
    3293          90 :             del = d*c
    3294          90 :             betaContinuedFraction = betaContinuedFraction*del
    3295          90 :             if (abs(del-1._RK) <=  eps) exit
    3296             :         end do
    3297           9 :         if (m > maxit) then
    3298             :         ! LCOV_EXCL_START
    3299             :             error stop
    3300             :             !call abortProgram( output_unit , 1 , 1 , &
    3301             :             !'Statitistics@getBetaContinuedFraction_DPR() failed: alpha or beta too big, or maxit too small.' )
    3302             :         end if
    3303             :         ! LCOV_EXCL_STOP
    3304          12 :     end function getBetaContinuedFraction_DPR
    3305             : 
    3306             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3307             : 
    3308             :     !> \brief
    3309             :     !> Return the one-sample Kolmogorov–Smirnov (KS) test results for the null hypothesis that the data
    3310             :     !> in vector `Point` comes from a distribution whose CDF is specified by the input `getCDF()` function.
    3311             :     !>
    3312             :     !> \param[in]       np      :   The number of points in the input vector `Point`.
    3313             :     !> \param[inout]    Point   :   The sample. On return, this array will be sorted in Ascending order.
    3314             :     !> \param[in]       getCDF  :   The function returning the Cumulative Distribution Function (CDF) of the sample.
    3315             :     !> \param[out]      statKS  :   The KS test statistic.
    3316             :     !> \param[out]      probKS  :   The `p`-value of the test, returned as a scalar value in the range `[0,1]`.
    3317             :     !>                              The output `probKS` is the probability of observing a test statistic as extreme as,
    3318             :     !>                              or more extreme than, the observed value under the null hypothesis.
    3319             :     !>                              Small values of `probKS` cast doubt on the validity of the null hypothesis.
    3320             :     !> \param[out]      Err     :   An object of class [Err_type](@ref err_mod::err_type) containing information
    3321             :     !>                              about the occurrence of any error during the KS test computation.
    3322             :     !>
    3323             :     !> \author
    3324             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    3325           3 :     subroutine doKS1(np,Point,getCDF,statKS,probKS,Err)
    3326             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3327             :         !DEC$ ATTRIBUTES DLLEXPORT :: doKS1
    3328             : #endif
    3329           9 :         use Sort_mod, only : sortAscending
    3330             :         use Err_mod, only: Err_type
    3331             : 
    3332             :         implicit none
    3333             : 
    3334             :         integer(IK)     , intent(in)    :: np
    3335             :         real(RK)        , intent(out)   :: statKS,probKS
    3336             :         real(RK)        , intent(inout) :: Point(np)
    3337             :         type(Err_type)  , intent(out)   :: Err
    3338             : 
    3339             :         character(*)    , parameter     :: PROCEDURE_NAME = MODULE_NAME//"@doKS1"
    3340           3 :         real(RK)                        :: npSqrt
    3341           3 :         real(RK)                        :: cdf,cdfObserved,dt,frac
    3342             :         integer(IK)                     :: j
    3343             : 
    3344             :         interface
    3345             :             function getCDF(x)
    3346             :                 use Constants_mod, only: RK
    3347             :                 real(RK), intent(in) :: x
    3348             :                 real(RK)             :: getCDF
    3349             :             end function getCDF
    3350             :         end interface
    3351             : 
    3352           3 :         call sortAscending(np,Point,Err)
    3353           3 :         if (Err%occurred) then
    3354             :         ! LCOV_EXCL_START
    3355             :             Err%msg = PROCEDURE_NAME//Err%msg
    3356             :             return
    3357             :         end if
    3358             :         ! LCOV_EXCL_STOP
    3359             : 
    3360           3 :         statKS = 0._RK
    3361           3 :         cdfObserved = 0._RK
    3362           3 :         npSqrt = np
    3363         153 :         do j = 1,np
    3364         150 :             frac = j/npSqrt
    3365         150 :             cdf = getCDF(Point(j))
    3366         150 :             dt = max( abs(cdfObserved-cdf) , abs(frac-cdf) )
    3367         150 :             if( dt > statKS ) statKS = dt
    3368         153 :             cdfObserved = frac
    3369             :         end do
    3370           3 :         npSqrt = sqrt(npSqrt)
    3371           3 :         probKS = getProbKS( (npSqrt+0.12_RK+0.11_RK/npSqrt)*statKS )
    3372             : 
    3373           3 :     end subroutine doKS1
    3374             : 
    3375             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3376             : 
    3377             :     !> \brief
    3378             :     !> Return the one-sample Kolmogorov–Smirnov (KS) test results for the assumption that the
    3379             :     !> points originate from a uniform distribution in [0,1]. So, all Point must be in [0,1] on input.
    3380             :     !>
    3381             :     !> \param[in]       np      :   The number of points in the input vector `Point`.
    3382             :     !> \param[inout]    Point   :   The sample. On return, this array will be sorted in Ascending order.
    3383             :     !> \param[out]      statKS  :   The KS test statistic.
    3384             :     !> \param[out]      probKS  :   The `p`-value of the test, returned as a scalar value in the range `[0,1]`.
    3385             :     !>                              The output `probKS` is the probability of observing a test statistic as extreme as,
    3386             :     !>                              or more extreme than, the observed value under the null hypothesis.
    3387             :     !>                              Small values of `probKS` cast doubt on the validity of the null hypothesis.
    3388             :     !> \param[out]      Err     :   An object of class [Err_type](@ref err_mod::err_type) containing information
    3389             :     !>                              about the occurrence of any error during the KS test computation.
    3390             :     !>
    3391             :     !> \author
    3392             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    3393           3 :     pure subroutine doUniformKS1(np,Point,statKS,probKS,Err)
    3394             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3395             :         !DEC$ ATTRIBUTES DLLEXPORT :: doUniformKS1
    3396             : #endif
    3397           3 :         use Sort_mod, only : sortAscending
    3398             :         use Err_mod, only: Err_type
    3399             : 
    3400             :         implicit none
    3401             : 
    3402             :         integer(IK)     , intent(in)    :: np
    3403             :         real(RK)        , intent(out)   :: statKS,probKS
    3404             :         real(RK)        , intent(inout) :: Point(np)
    3405             :         type(Err_type)  , intent(out)   :: Err
    3406             : 
    3407             :         character(*)    , parameter     :: PROCEDURE_NAME = MODULE_NAME//"@doUniformKS1"
    3408           3 :         real(RK)                        :: npSqrt
    3409           3 :         real(RK)                        :: cdf,cdfObserved,dt,frac
    3410             :         integer(IK)                     :: j
    3411             : 
    3412           3 :         call sortAscending(np,Point,Err)
    3413           3 :         if (Err%occurred) then
    3414             :             ! LCOV_EXCL_START
    3415             :             Err%msg = PROCEDURE_NAME//Err%msg
    3416             :             return
    3417             :             ! LCOV_EXCL_STOP
    3418             :         end if
    3419             : 
    3420           3 :         statKS = 0._RK
    3421           3 :         cdfObserved = 0._RK
    3422           3 :         npSqrt = np
    3423         153 :         do j = 1,np
    3424         150 :             frac = j/npSqrt
    3425         150 :             cdf = Point(j)
    3426         150 :             dt = max( abs(cdfObserved-cdf) , abs(frac-cdf) )
    3427         150 :             if( dt > statKS ) statKS = dt
    3428         153 :             cdfObserved = frac
    3429             :         end do
    3430           3 :         npSqrt = sqrt(npSqrt)
    3431           3 :         probKS = getProbKS( (npSqrt+0.12_RK+0.11_RK/npSqrt)*statKS )
    3432           3 :     end subroutine doUniformKS1
    3433             : 
    3434             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3435             : 
    3436             :     !> \brief
    3437             :     !> Return the two-sample Kolmogorov–Smirnov (KS) test results under the assumption that the
    3438             :     !> points originate from the same distribution. It is assumed that the two input arrays are sorted ascending on input.
    3439             :     !>
    3440             :     !> \param[in]       np1             :   The number of points in the input vector `SortedPoint1`.
    3441             :     !> \param[in]       np2             :   The number of points in the input vector `SortedPoint2`.
    3442             :     !> \param[inout]    SortedPoint1    :   The first input sorted sample. On input, it must be sorted in ascending-order.
    3443             :     !> \param[inout]    SortedPoint2    :   The second input sorted sample. On input, it must be sorted in ascending-order.
    3444             :     !> \param[out]      statKS          :   The KS test statistic.
    3445             :     !> \param[out]      probKS          :   The `p`-value of the test, returned as a scalar value in the range `[0,1]`.
    3446             :     !>                                      The output `probKS` is the probability of observing a test statistic as extreme as,
    3447             :     !>                                      or more extreme than, the observed value under the null hypothesis.
    3448             :     !>                                      Small values of `probKS` cast doubt on the validity of the null hypothesis.
    3449             :     !>
    3450             :     !> \author
    3451             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    3452         315 :     subroutine doSortedKS2(np1,np2,SortedPoint1,SortedPoint2,statKS,probKS)
    3453             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3454             :         !DEC$ ATTRIBUTES DLLEXPORT :: doSortedKS2
    3455             : #endif
    3456             :         integer(IK) , intent(in)    :: np1, np2
    3457             :         real(RK)    , intent(in)    :: SortedPoint1(np1), SortedPoint2(np2)
    3458             :         real(RK)    , intent(out)   :: statKS, probKS
    3459         315 :         real(RK)                    :: dummy1,dummy2,dt,np1_RK,np2_RK,npEffective,cdf1,cdf2
    3460             :         integer(IK)                 :: j1,j2
    3461         315 :         np1_RK = np1
    3462         315 :         np2_RK = np2
    3463         315 :         j1 = 1_IK
    3464         315 :         j2 = 1_IK
    3465         315 :         cdf1 = 0._RK
    3466         315 :         cdf2 = 0._RK
    3467         315 :         statKS = 0._RK
    3468      965805 :         do
    3469      966120 :             if ( j1<=np1 .and. j2<=np2 )then
    3470      965805 :                 dummy1 = SortedPoint1(j1)
    3471      965805 :                 dummy2 = SortedPoint2(j2)
    3472      965805 :                 if( dummy1 <= dummy2 ) then
    3473      483969 :                   cdf1 = j1 / np1_RK
    3474      483969 :                   j1 = j1 + 1
    3475             :                 endif
    3476      965805 :                 if( dummy2 <= dummy1 ) then
    3477      483978 :                   cdf2 = j2 / np2_RK
    3478      483978 :                   j2 = j2 + 1_IK
    3479             :                 endif
    3480      965805 :                 dt = abs(cdf2-cdf1)
    3481      965805 :                 if (dt>statKS) statKS = dt
    3482      965805 :                 cycle
    3483             :             end if
    3484         315 :             exit
    3485             :         end do
    3486         315 :         npEffective = sqrt( np1_RK * np2_RK / ( np1_RK + np2_RK ) )
    3487         315 :         probKS = getProbKS( ( npEffective + 0.12_RK + 0.11_RK / npEffective ) * statKS )
    3488           3 :     end subroutine doSortedKS2
    3489             : 
    3490             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3491             : 
    3492             :     !> \brief
    3493             :     !> Return the Kolmogorov–Smirnov (KS) probability.
    3494             :     !>
    3495             :     !> \author
    3496             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    3497         321 :     pure function getProbKS(lambda) result(probKS)
    3498             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3499             :         !DEC$ ATTRIBUTES DLLEXPORT :: getProbKS
    3500             : #endif
    3501             :         implicit none
    3502             :         real(RK)   , intent(in) :: lambda
    3503             :         real(RK)   , parameter  :: EPS1 = 0.001_RK, EPS2 = 1.e-8_RK
    3504             :         integer(IK), parameter  :: NITER = 100
    3505             :         integer(IK)             :: j
    3506         321 :         real(RK)                :: a2,fac,term,termbf
    3507             :         real(RK)                :: probKS
    3508         321 :         a2 = -2._RK*lambda**2
    3509         321 :         fac = 2._RK
    3510         321 :         probKS = 0._RK
    3511         321 :         termbf = 0._RK
    3512         938 :         do j = 1, NITER
    3513         938 :             term = fac*exp(a2*j**2)
    3514         938 :             probKS = probKS+term
    3515         938 :             if (abs(term) <= EPS1*termbf .or. abs(term) <= EPS2*probKS) return
    3516         617 :             fac = -fac
    3517         617 :             termbf = abs(term)
    3518             :         end do
    3519           0 :         probKS = 1._RK
    3520         636 :     end function getProbKS
    3521             : 
    3522             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3523             : 
    3524             :     !> \brief
    3525             :     !> Returns the uniform CDF on support [0,1). This is rather redundant, aint it? but sometimes, needed.
    3526             :     !>
    3527             :     !> \param[in] x : The point at which the CDF must be computed.
    3528             :     !>
    3529             :     !> \author
    3530             :     !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
    3531           3 :     pure function getUniformCDF(x)
    3532             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3533             :         !DEC$ ATTRIBUTES DLLEXPORT :: getUniformCDF
    3534             : #endif
    3535             :         implicit none
    3536             :         real(RK), intent(in) :: x
    3537             :         real(RK)             :: getUniformCDF
    3538           3 :         getUniformCDF = x
    3539         321 :     end function getUniformCDF
    3540             : 
    3541             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3542             : 
    3543             :     !> \brief
    3544             :     !> Return the 1-D histogram (Density plot) of the input vector `X`.
    3545             :     !> The number of bins in the `X` range (`nxbin`) is determined by the user.
    3546             :     !> The range of `X`, `[xmin, xmax]`, should be also given by the user.
    3547             :     !> The program returns two arrays of `Xbin` and `Density(x)` as output.
    3548             :     !>
    3549             :     !> \param[in]       method          :   The method by which the hist count should be returned:
    3550             :     !>                                      + `"pdf"`   : Return the probability density function histogram.
    3551             :     !>                                      + `"count"` : Return the count histogram.
    3552             :     !> \param[in]       xmin            :   The minimum of the histogram binning.
    3553             :     !> \param[in]       xmax            :   The maximum of the histogram binning.
    3554             :     !> \param[in]       nxbin           :   The number of histogram bins.
    3555             :     !> \param[in]       np              :   The length of input vector `X`.
    3556             :     !> \param[in]       X               :   The vector of length `nxbin` of values to be binned.
    3557             :     !> \param[out]      Xbin            :   The vector of length `nxbin` of values representing the bin left corners.
    3558             :     !> \param[out]      Density         :   The vector of length `nxbin` of values representing the densities in each bin.
    3559             :     !> \param[out]      errorOccurred   :   The logical output flag indicating whether error has occurred.
    3560             :     !>
    3561             :     !> \author
    3562             :     !> Amir Shahmoradi, Sep 1, 2017, 12:30 AM, ICES, UT Austin
    3563           9 :     subroutine getHist1D(method, xmin, xmax, nxbin, np, X, Xbin, Density, errorOccurred)
    3564             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3565             :         !DEC$ ATTRIBUTES DLLEXPORT :: getHist1D
    3566             : #endif
    3567             :         implicit none
    3568             :         character(*), intent(in)  :: method
    3569             :         integer(IK) , intent(in)  :: np,nxbin
    3570             :         real(RK)    , intent(in)  :: xmin,xmax
    3571             :         real(RK)    , intent(in)  :: X(np)
    3572             :         real(RK)    , intent(out) :: Xbin(nxbin), Density(nxbin)
    3573             :         logical     , intent(out) :: errorOccurred
    3574           9 :         real(RK)                  :: xbinsize
    3575             :         integer(IK)               :: i, ip, thisXbin, npEffective
    3576             : 
    3577           9 :         xbinsize = (xmax-xmin) / real(nxbin,kind=RK)
    3578         189 :         Xbin = [ (xmin + real(i-1,kind=RK)*xbinsize,i=1,nxbin) ]
    3579             : 
    3580          99 :         Density = 0._RK
    3581           9 :         npEffective = 0_IK
    3582         459 :         do ip = 1, np
    3583         459 :             if (X(ip)>=xmin .and. X(ip)<xmax) then
    3584         450 :                 npEffective = npEffective + 1_IK
    3585         450 :                 thisXbin = getBin(X(ip),xmin,nxbin,xbinsize)
    3586         450 :                 Density(thisXbin) = Density(thisXbin) + 1._RK
    3587             :             end if
    3588             :         end do
    3589             : 
    3590          99 :         Xbin = Xbin + 0.5_RK * xbinsize
    3591             : 
    3592           9 :         if(method=="count") then
    3593           3 :             errorOccurred = .false.
    3594           3 :             return
    3595           6 :         elseif(method=="pdf") then
    3596          33 :             Density = Density / real(npEffective,kind=RK)
    3597           3 :             errorOccurred = .false.
    3598             :         else
    3599           3 :             errorOccurred = .true.
    3600             :         end if
    3601             : 
    3602          12 :     end subroutine getHist1D
    3603             : 
    3604             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3605             : 
    3606             :     !> \brief
    3607             :     !> Return the 2-D histogram (Density plot) of a set of data points with (X,Y) coordinates.
    3608             :     !> The number of bins in the `X` and `Y` directions (`[nxbin, nybin]`) are determined by the user.
    3609             :     !> The range of `X` and `Y` (`xmin`,`xmax`,`ymin`,`ymax`) should also be given by the user.
    3610             :     !> The program returns three arrays of `Xbin`, `Ybin`, and `Density(y,x)` as the output.
    3611             :     !>
    3612             :     !> \param[in]       histType        :   The method by which the normalization of the histogram counts should be done:
    3613             :     !>                                      + `"count"`     : Return the count histogram.
    3614             :     !>                                      + `"pdf"`       : Return the probability density function (PDF) histogram.
    3615             :     !>                                      + `"pdf(y|x)"`  : Return the conditional PDF of `y` given `x` histogram.
    3616             :     !>                                      + `"pdf(x|y)"`  : Return the conditional PDF of `x` given `y` histogram.
    3617             :     !> \param[in]       xmin            :   The minimum of the histogram binning along the x-axis.
    3618             :     !> \param[in]       xmax            :   The maximum of the histogram binning along the x-axis.
    3619             :     !> \param[in]       ymin            :   The minimum of the histogram binning along the y-axis.
    3620             :     !> \param[in]       ymax            :   The maximum of the histogram binning along the y-axis.
    3621             :     !> \param[in]       nxbin           :   The number of histogram bins along the x-axis.
    3622             :     !> \param[in]       nybin           :   The number of histogram bins along the y-axis.
    3623             :     !> \param[in]       np              :   The length of input vector `X`.
    3624             :     !> \param[in]       X               :   The vector of length `nxbin` of values to be binned.
    3625             :     !> \param[in]       Y               :   The vector of length `nybin` of values to be binned.
    3626             :     !> \param[out]      Xbin            :   The vector of length `nxbin` of values representing the bin centers.
    3627             :     !> \param[out]      Ybin            :   The vector of length `nybin` of values representing the bin centers.
    3628             :     !> \param[out]      Density         :   The array of shape `(nybin,nxbin)` of values representing the densities per bin.
    3629             :     !> \param[out]      errorOccurred   :   A logical output value indicating whether an error has occurred.
    3630             :     !>
    3631             :     !> \author
    3632             :     !> Amir Shahmoradi, Sep 1, 2017, 12:00 AM, ICES, UT Austin
    3633          15 :     subroutine getHist2D(histType,xmin,xmax,ymin,ymax,nxbin,nybin,np,X,Y,Xbin,Ybin,Density,errorOccurred)
    3634             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3635             :         !DEC$ ATTRIBUTES DLLEXPORT :: getHist2D
    3636             : #endif
    3637             :         !use, intrinsic :: iso_fortran_env, only: output_unit
    3638           9 :         use String_mod, only: getLowerCase
    3639             :         implicit none
    3640             :         character(*), intent(in)  :: histType
    3641             :         integer(IK) , intent(in)  :: np,nxbin,nybin
    3642             :         real(RK)    , intent(in)  :: xmin,xmax,ymin,ymax
    3643             :         real(RK)    , intent(in)  :: X(np),Y(np)
    3644             :         real(RK)    , intent(out) :: Xbin(nxbin),Ybin(nybin),Density(nybin,nxbin)
    3645             :         logical     , intent(out) :: errorOccurred
    3646          15 :         character(:), allocatable :: method
    3647          15 :         real(RK)                  :: xbinsize,ybinsize
    3648             :         integer(IK)               :: i, ip, thisXbin, thisYbin, npEffective
    3649             : 
    3650          15 :         errorOccurred = .false.
    3651             : 
    3652          15 :         xbinsize = (xmax-xmin) / real(nxbin,kind=RK)
    3653          15 :         ybinsize = (ymax-ymin) / real(nybin,kind=RK)
    3654         255 :         Xbin = [ (xmin+real(i-1,kind=RK)*xbinsize,i=1,nxbin) ]
    3655         225 :         Ybin = [ (ymin+real(i-1,kind=RK)*ybinsize,i=1,nybin) ]
    3656             : 
    3657         975 :         Density = 0._RK
    3658          15 :         npEffective = 0_IK
    3659         765 :         do ip = 1,np
    3660         765 :             if (X(ip)>=xmin .and. X(ip)<xmax .and. Y(ip)>=ymin .and. Y(ip)<ymax) then
    3661         750 :                 npEffective = npEffective + 1_IK
    3662         750 :                 thisXbin = getBin(X(ip),xmin,nxbin,xbinsize)
    3663         750 :                 thisYbin = getBin(Y(ip),ymin,nybin,ybinsize)
    3664         750 :                 Density(thisYbin,thisXbin) = Density(thisYbin,thisXbin) + 1._RK
    3665             :             end if
    3666             :         end do
    3667             : 
    3668         135 :         Xbin = Xbin + 0.5_RK * xbinsize
    3669         120 :         Ybin = Ybin + 0.5_RK * ybinsize
    3670          15 :         method = getLowerCase(trim(adjustl(histType)))
    3671          15 :         if(method=="pdf") then
    3672         195 :             Density = Density / real(npEffective,kind=RK)
    3673          12 :         elseif(method=="pdf(y|x)") then
    3674          27 :             do i = 1,nxbin
    3675         363 :                 Density(1:nybin,i) = Density(1:nybin,i) / sum(Density(1:nybin,i))
    3676             :             end do
    3677           9 :         elseif(method=="pdf(x|y)") then
    3678          24 :             do i = 1,nybin
    3679         360 :                 Density(i,1:nxbin) = Density(i,1:nxbin) / sum(Density(i,1:nxbin))
    3680             :             end do
    3681           6 :         elseif(method=="count") then
    3682           3 :             return
    3683             :         else
    3684           3 :             errorOccurred = .true.
    3685             :         end if
    3686             : 
    3687          15 :     end subroutine getHist2D
    3688             : 
    3689             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3690             : 
    3691             :     !> \brief
    3692             :     !> Given the range of the variable `x`, `xmin:xmin+binsize*nbin`, and the number of bins, `nbin`, with which
    3693             :     !> this range is divided, find which bin the input value `x` falls among `[1:nbin]` bins.
    3694             :     !> The output `ibin` is the number that identifies the bin.
    3695             :     !>
    3696             :     !> \param[in]       x               :   The input value whose bin ID is to be found.
    3697             :     !> \param[in]       lowerBound      :   The lower limit on the value of `x`.
    3698             :     !> \param[in]       nbin            :   The number of bins to be considered starting from `lowerBound`.
    3699             :     !> \param[in]       binsize         :   The size of the bins. It must be exactly `(xmax - xmin) / nbin`.
    3700             :     !>
    3701             :     !> \return
    3702             :     !> `ibin` : The ID of the bin to which the input value `x` belongs.
    3703             :     !>
    3704             :     !> \warning
    3705             :     !> If `x <= xmin` or `x xmin + nbin * binsize`, then `ibin = -1` will be returned to indicate error.
    3706             :     !>
    3707             :     !> \remark
    3708             :     !> If `bmin < x <= bmax` then `x` belongs to this bin.
    3709             :     !>
    3710             :     !> \todo
    3711             :     !> The performance and interface of this routine can be significantly improved.
    3712             :     !> It is more sensible to pass a contiguous array of bin edges as input instead of `lowerBound` and `binsize`.
    3713             :     !> If the input point is not within any bins, an index of zero should be returned.
    3714             :     !>
    3715             :     !> \author
    3716             :     !> Version 3.0, Sep 1, 2017, 11:12 AM, Amir Shahmoradi, ICES, The University of Texas at Austin.
    3717        1950 :     pure function getBin(x,lowerBound,nbin,binsize) result(ibin)
    3718             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3719             :         !DEC$ ATTRIBUTES DLLEXPORT :: getBin
    3720             : #endif
    3721             : 
    3722             :         implicit none
    3723             :         integer(IK), intent(in) :: nbin
    3724             :         real(RK)   , intent(in) :: x,lowerBound,binsize
    3725        1950 :         real(RK)                :: xmin,xmid
    3726             :         integer(IK)             :: ibin,minbin,midbin,maxbin
    3727             : 
    3728        1950 :         if (x<lowerBound .or. x>=lowerBound+nbin*binsize) then
    3729             :             ! LCOV_EXCL_START
    3730             :             ibin = -1_IK
    3731             :             return
    3732             :             ! LCOV_EXCL_STOP
    3733             :         end if
    3734             : 
    3735        1950 :         minbin = 1
    3736        1950 :         maxbin = nbin
    3737        1950 :         xmin = lowerBound
    3738        5547 :         loopFindBin: do
    3739        7497 :             midbin = (minbin+maxbin) / 2
    3740        7497 :             xmid = xmin + midbin*binsize
    3741        7497 :             if (x<xmid) then
    3742        3021 :                 if (minbin==midbin) then
    3743          54 :                     ibin = minbin
    3744          54 :                     exit loopFindBin
    3745             :                 end if
    3746        2967 :                 maxbin = midbin
    3747        2967 :                 cycle loopFindBin
    3748             :             else
    3749        4476 :                 if (minbin==midbin) then
    3750        1896 :                     ibin = maxbin
    3751        1896 :                     exit loopFindBin
    3752             :                 end if
    3753        2580 :                 minbin = midbin
    3754        2580 :                 cycle loopFindBin
    3755             :             end if
    3756             :         end do loopFindBin
    3757             : 
    3758        1965 :     end function getBin
    3759             : 
    3760             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3761             : 
    3762             :     !> \brief
    3763             :     !> Return the quantiles of an input sample of points, given the input quantile probabilities.
    3764             :     !>
    3765             :     !> \param[in]       np                          :   The number of points in the input sample.
    3766             :     !> \param[in]       nq                          :   The number of output quantiles.
    3767             :     !> \param[in]       SortedQuantileProbability   :   A sorted ascending-order vector of probabilities at which the quantiles will be returned.
    3768             :     !> \param[in]       Point                       :   The vector of length `np` representing the input sample.
    3769             :     !> \param[in]       Weight                      :   The vector of length `np` representing the weights of the points in the input sample.
    3770             :     !> \param[in]       sumWeight                   :   The sum of the vector weights of the points: `sum(Weight)`.
    3771             :     !>
    3772             :     !> \return
    3773             :     !> `Quantile` : The output vector of length `nq`, representing the quantiles corresponding to the input `SortedQuantileProbability` probabilities.
    3774             :     !>
    3775             :     !> \author
    3776             :     !> Amir Shahmoradi, Sep 1, 2017, 12:00 AM, ICES, UT Austin
    3777       10969 :     pure function getQuantile(np,nq,SortedQuantileProbability,Point,Weight,sumWeight) result(Quantile)
    3778             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
    3779             :         !DEC$ ATTRIBUTES DLLEXPORT :: getQuantile
    3780             : #endif
    3781        1950 :         use Constants_mod, only: IK, RK, NEGINF_RK
    3782             :         use Sort_mod, only: indexArray
    3783             :         use Err_mod, only: Err_type
    3784             :         implicit none
    3785             :         integer(IK) , intent(in)            :: np, nq
    3786             :         real(RK)    , intent(in)            :: SortedQuantileProbability(nq), Point(np)
    3787             :         integer(IK) , intent(in), optional  :: Weight(np), sumWeight
    3788             :         real(RK)                            :: Quantile(nq)
    3789        1688 :         integer(IK)                         :: ip, iq, iw, weightCounter, Indx(np), SortedQuantileDensity(nq)
    3790         844 :         type(Err_type)                      :: Err
    3791         844 :         call indexArray(np,Point,Indx,Err)
    3792         844 :         if (Err%occurred) then
    3793             :             ! LCOV_EXCL_START
    3794             :             Quantile = NEGINF_RK
    3795             :             return
    3796             :             ! LCOV_EXCL_STOP
    3797             :         end if
    3798         844 :         iq = 1_IK
    3799         844 :         if (present(sumWeight)) then
    3800         841 :             weightCounter = 0_IK
    3801        8410 :             SortedQuantileDensity = nint( SortedQuantileProbability * sumWeight )
    3802      235888 :             loopWeighted: do ip = 1, np
    3803      937957 :                 do iw = 1, Weight(Indx(ip))
    3804      702859 :                     weightCounter = weightCounter + 1_IK
    3805      937906 :                     if (weightCounter>=SortedQuantileDensity(iq)) then
    3806        7468 :                         Quantile(iq) = Point(Indx(ip))
    3807        7468 :                         iq = iq + 1_IK
    3808        7468 :                         if (iq>nq) exit loopWeighted
    3809             :                     end if
    3810             :                 end do
    3811             :             end do loopWeighted
    3812             :         else
    3813          30 :             SortedQuantileDensity = nint( SortedQuantileProbability * np )
    3814         150 :             loopNonWeighted: do ip = 1, np
    3815         150 :                 if (ip>=SortedQuantileDensity(iq)) then
    3816          27 :                     Quantile(iq) = Point(Indx(ip))
    3817          27 :                     iq = iq + 1_IK
    3818          27 :                     if (iq>nq) exit loopNonWeighted
    3819             :                 end if
    3820             :             end do loopNonWeighted
    3821             :         end if
    3822         945 :         if (iq<=nq) Quantile(iq:nq) = Point(Indx(np))
    3823        2532 :     end function getQuantile
    3824             : 
    3825             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3826             : 
    3827             : end module Statistics_mod ! LCOV_EXCL_LINE

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