Line data    Source code

       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       3             : !!!!                                                                                                                            !!!!
       4             : !!!!    ParaMonte: Parallel Monte Carlo and Machine Learning Library.                                                           !!!!
       5             : !!!!                                                                                                                            !!!!
       6             : !!!!    Copyright (C) 2012-present, The Computational Data Science Lab                                                          !!!!
       7             : !!!!                                                                                                                            !!!!
       8             : !!!!    This file is part of the ParaMonte library.                                                                             !!!!
       9             : !!!!                                                                                                                            !!!!
      10             : !!!!    LICENSE                                                                                                                 !!!!
      11             : !!!!                                                                                                                            !!!!
      12             : !!!!       https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md                                                          !!!!
      13             : !!!!                                                                                                                            !!!!
      14             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      15             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      16             : 
      17             : !>  \brief
      18             : !>  This include file contains implementations of the procedures in module [pm_distNormShell](@ref pm_distNormShell).
      19             : !>
      20             : !>  \finmain
      21             : !>
      22             : !>  \author
      23             : !>  \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
      24             : 
      25             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      26             : 
      27             : #define CHECK_POSITIVE_RADIUS \
      28             : CHECK_ASSERTION(__LINE__, all([0._RKC < radius]), \
      29             : SK_"@getNormShellLogUDF(): The condition `all([0._RKC < radius])` must hold. width = "//getStr(radius)) ! fpp
      30             : 
      31             : #define CHECK_POSITIVE_WIDTH \
      32             : CHECK_ASSERTION(__LINE__, all([0._RKC < width]), \
      33             : SK_"@getNormShellLogUDF(): The condition `all([0._RKC < width])` must hold. width = "//getStr(width)); ! fpp
      34             : 
      35             : #define CHECK_LENGTH_RADIUS \
      36             : CHECK_ASSERTION(__LINE__, size(radius, 1, IK) == size(center, 2, IK), \
      37             : SK_"@getNormShellLogUDF(): The condition `size(radius, 1) == size(center, 2)` must hold. size(radius, 1), size(center, 2) = "//\
      38             : getStr([size(radius, 1, IK), size(center, 2, IK)])) ! fpp
      39             : 
      40             : #define CHECK_LENGTH_WIDTH \
      41             : CHECK_ASSERTION(__LINE__, size(width, 1, IK) == size(center, 2, IK), \
      42             : SK_"@getNormShellLogUDF(): The condition `size(width, 1) == size(center, 2)` must hold. size(width, 1), size(center, 2) = "//\
      43             : getStr([size(width, 1, IK), size(center, 2, IK)])) ! fpp
      44             : 
      45             :         !%%%%%%%%%%%%%%%%%%%%%%%%%
      46             : #if     getNormShellLogUDF_ENABLED
      47             :         !%%%%%%%%%%%%%%%%%%%%%%%%%
      48             : 
      49             :         real(RKC), parameter :: ZERO = 0._RKC, ONE = 1._RKC !, TWO = 2._RKC !, PI = acos(-1._RKC)
      50             : #if     CI_ENABLED
      51     1473006 :         CHECK_ASSERTION(__LINE__, size(X, 1, IK) == size(center, 1, IK), \
      52             :         SK_"@getNormShellLogUDF(): The condition `size(X, 1) == size(center, 1)` must hold. size(X, 1), size(center, 1) = "//\
      53             :         getStr([size(X, 1, IK), size(center, 1, IK)])) ! fpp
      54             : 
      55     3928012 :         CHECK_ASSERTION(__LINE__, all(size(X, 1, IK) == [size(invCov, 1, IK), size(invCov, 2, IK)]), \
      56             :         SK_"@getNormShellLogUDF(): The condition `all(size(X, 1) == [size(invCov, 1), size(invCov, 2)])` must hold. size(X, 1), shape(invCov) = "//\
      57             :         getStr([size(X, 1), shape(invCov)])) ! fpp
      58             : 
      59     5401014 :         CHECK_ASSERTION(__LINE__, size(invCov, rank(invCov)) == size(center, rank(center)), \
      60             :         SK_"@getNormShellLogUDF(): The condition `size(invCov, rank(invCov)) == size(center, rank(invCov))` must hold. shape(invCov), shape(center) = "//\
      61             :         getStr([shape(invCov), shape(center)])) ! fpp
      62             : #elif   !DD_ENABLED
      63             : #error  "Unrecognized interface."
      64             : #endif
      65             :         ! Compute the density function(s).
      66             : #if     D1_ENABLED && (One_ENABLED || (Mix_ENABLED && CI_ENABLED))
      67             : #if     DD_ENABLED
      68             : #define MAHAL_SQ dot_product(X, X)
      69             : #elif   CI_ENABLED
      70             : #define MAHAL_SQ getMahalSq(X, invCov, center)
      71             : #else
      72             : #error  "Unrecognized interface."
      73             : #endif
      74      491003 :         if (present(width)) then
      75     1475002 :             CHECK_POSITIVE_WIDTH
      76      491002 :             if (present(radius)) then
      77     1475002 :                 CHECK_POSITIVE_RADIUS
      78     1475002 :                 logUDF = -.5_RKC * ((sqrt(MAHAL_SQ) - radius) / width)**2
      79             :             else
      80           0 :                 logUDF = -.5_RKC * ((sqrt(MAHAL_SQ) - ONE) / width)**2
      81             :             end if
      82             :         else
      83           1 :             if (present(radius)) then
      84           0 :                 CHECK_POSITIVE_RADIUS
      85           0 :                 logUDF = -.5_RKC *  (sqrt(MAHAL_SQ) - radius)**2
      86             :             else
      87           2 :                 logUDF = -.5_RKC *  (sqrt(MAHAL_SQ) - ONE)**2
      88             :             end if
      89             :         end if
      90             : #elif   Mix_ENABLED && D1_ENABLED && DD_ENABLED
      91           0 :         CHECK_ASSERTION(__LINE__, size(width, 1, IK) == size(radius, 1, IK), \
      92             :         SK_"@getNormShellLogUDF(): The condition `size(width) == size(radius)` must hold. size(width), size(radius) = "//\
      93             :         getStr([size(width, 1, IK), size(radius, 1, IK)])) ! fpp
      94           0 :         CHECK_POSITIVE_RADIUS
      95           0 :         CHECK_POSITIVE_WIDTH
      96           0 :         logUDF = -.5_RKC * ((sqrt(sum(X**2)) - radius) / width)**2
      97             : #else
      98             : #error  "Unrecognized interface."
      99             : #endif
     100             : 
     101             : #else
     102             :         !%%%%%%%%%%%%%%%%%%%%%%%%
     103             : #error  "Unrecognized interface."
     104             :         !%%%%%%%%%%%%%%%%%%%%%%%%
     105             : #endif
     106             : 
     107             : #undef  CHECK_POSITIVE_RADIUS
     108             : #undef  CHECK_POSITIVE_WIDTH
     109             : #undef  CHECK_LENGTH_RADIUS
     110             : #undef  CHECK_LENGTH_WIDTH
     111             : #undef  MAHAL_SQ