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 the implementation of procedures in [pm_distUnifSphere](@ref pm_distUnifSphere).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, April 23, 2017, 1:36 AM, Institute for Computational Engineering and Sciences (ICES), University of Texas at Austin
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : ! Set the uniform RNG
28 : #if RNGD_ENABLED
29 : #define RNG
30 : #elif RNGF_ENABLED || RNGX_ENABLED
31 : #define RNG rng,
32 : #elif getUSR_ENABLED || setUSR_ENABLED
33 : #error "Unrecognized interface."
34 : #endif
35 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
36 : #if getUnifSphereLogPDF_ENABLED
37 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
38 :
39 : #if D0_ENABLED
40 4 : CHECK_ASSERTION(__LINE__, 0_IK < ndim, SK_"@getUnifSphereLogPDF(): The condition `0 < ndim` must hold. ndim = "//getStr(ndim))
41 4 : logPDF = (1 - ndim) * logRadius - getLogVolUnitBall(real(ndim, TKC)) - log(real(ndim, TKC))
42 : #else
43 : #error "Unrecognized interface."
44 : #endif
45 :
46 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
47 : #elif D1_ENABLED && setUSR_ENABLED
48 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
49 :
50 : integer(IK) :: idim, ndim
51 : real(TKC) :: normfac
52 : #if AC_ENABLED
53 : #define UNIFSPHERERAND unifSphereRand
54 1008 : real(TKC) :: unifSphereRand(size(rand, 1, IK))
55 : #elif DC_ENABLED
56 : #define UNIFSPHERERAND rand
57 : #else
58 : #error "Unrecognized interface."
59 : #endif
60 : ! Perform runtime bound checks.
61 : #if AM_ENABLED
62 1512 : CHECK_ASSERTION(__LINE__, size(rand, 1, IK) == size(mean, 1, IK), SK_"@setUnifSphereRand(): The condition `size(rand) == size(mean)` must hold. size(rand), size(mean) = "//getStr([size(rand, 1, IK), size(mean, 1, IK)]))
63 : #elif AC_ENABLED
64 2016 : CHECK_ASSERTION(__LINE__, all(size(rand, 1, IK) == shape(chol, IK)), SK_"@setUnifSphereRand(): The condition `all(size(rand) == shape(chol))` must hold. size(rand), shape(chol) = "//getStr([size(rand, 1, IK), shape(chol, IK)]))
65 : #elif !(DM_ENABLED || DC_ENABLED)
66 : #error "Unrecognized interface."
67 : #endif
68 151 : ndim = size(rand, 1, IK)
69 : do
70 0 : normfac = 0._TKC
71 2721 : do idim = 1_IK, ndim
72 1814 : call setNormRand(RNG UNIFSPHERERAND(idim))
73 2721 : normfac = normfac + UNIFSPHERERAND(idim)**2
74 : end do
75 : ! Ensure the vector is not origin. Highly unlikely but possible.
76 907 : if (0._TKC < normfac) exit
77 : end do
78 403 : normfac = 1._TKC / sqrt(normfac)
79 : #if DM_ENABLED && !AC_ENABLED
80 453 : UNIFSPHERERAND = UNIFSPHERERAND * normfac ! a uniform random point on the surface of nd-sphere.
81 : #elif AM_ENABLED && !AC_ENABLED
82 756 : UNIFSPHERERAND = UNIFSPHERERAND * normfac + mean ! a uniform random point on the surface of nd-sphere with given mean.
83 : #elif AC_ENABLED
84 : ! Define the indexing rules.
85 : #if XLD_ENABLED
86 : #define GET_INDEX(I,J)I,J
87 : #elif UXD_ENABLED
88 : #define GET_INDEX(I,J)J,I
89 : #elif D1_ENABLED && setUSR_ENABLED
90 : #error "Unrecognized interface."
91 : #endif
92 : ! Define the default mean.
93 : #if DM_ENABLED
94 : #define MEAN_PLUS(I)
95 : #elif AM_ENABLED
96 : #define MEAN_PLUS(I) mean(I) +
97 : #else
98 : #error "Unrecognized interface."
99 : #endif
100 : ! Separate the first to allow the possibility of adding an optional `mean`.
101 1512 : rand(1 : ndim) = MEAN_PLUS(1 : ndim) chol(GET_INDEX(1 : ndim, 1)) * unifSphereRand(1) * normfac
102 1008 : do idim = 2_IK, ndim
103 1512 : rand(idim : ndim) = rand(idim : ndim) + chol(GET_INDEX(idim : ndim, idim)) * unifSphereRand(idim) * normfac
104 : end do
105 : #else
106 : #error "Unrecognized interface."
107 : #endif
108 : #undef UNIFSPHERERAND
109 : #undef MEAN_PLUS
110 : #undef GET_INDEX
111 : #undef MEAN
112 :
113 : !%%%%%%%%%%%%%
114 : #elif getUSR_ENABLED
115 : !%%%%%%%%%%%%%
116 :
117 : #if AM_ENABLED && DC_ENABLED
118 101 : call setUnifSphereRand(RNG rand, mean)
119 : #elif DM_ENABLED && AC_ENABLED
120 101 : call setUnifSphereRand(RNG rand, chol, subset)
121 : #elif AM_ENABLED && AC_ENABLED
122 101 : call setUnifSphereRand(RNG rand, mean, chol, subset)
123 : #else
124 : #error "Unrecognized interface."
125 : #endif
126 :
127 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
128 : #elif D2_ENABLED && setUSR_ENABLED
129 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
130 :
131 : integer(IK) :: ipnt, ndim
132 0 : ndim = size(rand, 1, IK)
133 0 : do ipnt = 1_IK, size(rand, 2, IK)
134 : #if DM_ENABLED && DC_ENABLED
135 0 : call setUnifSphereRand(RNG rand(1:ndim, ipnt))
136 : #elif AM_ENABLED && DC_ENABLED
137 0 : call setUnifSphereRand(RNG rand(1:ndim, ipnt), mean)
138 : #elif DM_ENABLED && AC_ENABLED
139 0 : call setUnifSphereRand(RNG rand(1:ndim, ipnt), chol, subset)
140 : #elif AM_ENABLED && AC_ENABLED
141 0 : call setUnifSphereRand(RNG rand(1:ndim, ipnt), mean, chol, subset)
142 : #else
143 : #error "Unrecognized interface."
144 : #endif
145 : end do
146 :
147 : #else
148 : !%%%%%%%%%%%%%%%%%%%%%%%%
149 : #error "Unrecognized interface."
150 : !%%%%%%%%%%%%%%%%%%%%%%%%
151 : #endif
152 : #undef RNG
|
