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 file contains implementations of procedures [pm_distanceEuclid](@ref pm_distanceEuclid).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, September 1, 2017, 12:00 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
28 : #if getDisEuclid_ENABLED && XYZ_ENABLED
29 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 :
31 : real(TKC) :: absx, absy, absz, maximum
32 2 : absx = abs(x)
33 2 : absy = abs(y)
34 2 : absz = abs(z)
35 2 : maximum = max(absx, absy, absz)
36 2 : if(maximum == 0._TKC .or. maximum > huge(maximum)) then
37 : ! maximum can be zero for max(0,nan,0) adding all three entries together will make sure nan will not disappear.
38 0 : distance = absx + absy + absz
39 : else
40 2 : distance = maximum * sqrt((absx / maximum)**2 + (absy / maximum)**2 + (absz / maximum)**2)
41 : end if
42 :
43 : !%%%%%%%%%%%%%%%%%%%
44 : #elif getDisEuclid_ENABLED
45 : !%%%%%%%%%%%%%%%%%%%
46 :
47 : #if D1_D1_ENABLED || D1_D2_ENABLED || D2_D1_ENABLED || D2_D2_ENABLED
48 : #define REF ref,
49 : #elif D1_XX_ENABLED || D2_XX_ENABLED
50 : #define REF
51 : #else
52 : #error "Unrecognized interface."
53 : #endif
54 28593 : if (present(method)) then
55 : select type (method)
56 : type is (euclid_type)
57 9529 : call setDisEuclid(distance, point, REF method)
58 : type is (euclidu_type)
59 9529 : call setDisEuclid(distance, point, REF method)
60 : type is (euclidsq_type)
61 9529 : call setDisEuclid(distance, point, REF method)
62 : class default
63 : error stop MODULE_NAME//SK_"@getDisEuclid(): Unsupported input value for `method`." ! LCOV_EXCL_LINE
64 : end select
65 6000 : return
66 : end if
67 6 : call setDisEuclid(distance, point, REF euclid)
68 : #undef REF
69 :
70 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
71 : #elif setDisEuclid_ENABLED && (D1_XX_ENABLED || D1_D1_ENABLED)
72 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
73 :
74 : integer(IK) :: idim
75 : #if D1_XX_ENABLED
76 : #define GET_DIFF(PNT,REF)PNT
77 : #elif D1_D1_ENABLED
78 : #define GET_DIFF(PNT,REF)(PNT - REF)
79 13579266 : CHECK_ASSERTION(__LINE__, size(point, 1, IK) == size(ref, 1, IK), SK_"@setDisEuclid(): The condition `size(point) == size(ref)` must hold. size(point), size(ref) = "//getStr([size(point, 1, IK), size(ref, 1, IK)]))
80 : #else
81 : #error "Unrecognized interface."
82 : #endif
83 : #if MED_ENABLED
84 : block
85 : real(TKC) :: invmax, maximum
86 : ! First find the maximum.
87 757585 : maximum = -huge(maximum)
88 9976976 : do idim = 1_IK, size(point, 1, IK)
89 8444552 : invmax = abs(GET_DIFF(point(idim),ref(idim))) ! placeholder.
90 9976976 : if (maximum < invmax) maximum = invmax
91 : end do
92 1532424 : if(maximum == 0._TKC .or. huge(maximum) < maximum) then
93 248366 : distance = sum(GET_DIFF(point,ref)) ! Ensure propagation of potential nan in the vector.
94 : else
95 1494278 : distance = 0._TKC
96 1494278 : invmax = 1._TKC / maximum
97 9728610 : do idim = 1_IK, size(point, 1, IK)
98 9728610 : distance = distance + (GET_DIFF(point(idim),ref(idim)) * invmax)**2
99 : end do
100 1494278 : distance = maximum * sqrt(distance)
101 : end if
102 : end block
103 : #elif MEU_ENABLED || MEQ_ENABLED
104 3017333 : distance = 0._TKC
105 19950818 : do idim = 1_IK, size(point, 1, IK)
106 19950818 : distance = distance + GET_DIFF(point(idim),ref(idim))**2
107 : end do
108 : #if MEU_ENABLED
109 1531923 : distance = sqrt(distance)
110 : #endif
111 : #else
112 : #error "Unrecognized interface."
113 : #endif
114 : #undef GET_DIFF
115 :
116 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
117 : #elif setDisEuclid_ENABLED && D1_D2_ENABLED
118 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119 :
120 : integer(IK) :: iref, ndim
121 285225 : ndim = size(point, 1, IK)
122 1711350 : CHECK_ASSERTION(__LINE__, size(ref, 2, IK) == size(distance, 1, IK), SK_"@getDisEuclid(): The condition `size(ref, 2) == size(distance)` must hold. shape(ref), size(distance) = "//getStr([shape(ref, IK), size(distance, 1, IK)]))
123 3291138 : do iref = 1_IK, size(ref, 2, IK)
124 3291138 : call setDisEuclid(distance(iref), point, ref(1:ndim, iref), method)
125 : end do
126 :
127 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
128 : #elif setDisEuclid_ENABLED && (D2_XX_ENABLED || D2_D1_ENABLED)
129 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
130 :
131 : integer(IK) :: ipnt, ndim
132 12002 : ndim = size(point, 1, IK)
133 72012 : CHECK_ASSERTION(__LINE__, size(point, 2, IK) == size(distance, 1, IK), SK_"@getDisEuclid(): The condition `size(point, 2) == size(distance)` must hold. shape(point), size(distance) = "//getStr([shape(point, IK), size(distance, 1, IK)]))
134 137246 : do ipnt = 1_IK, size(point, 2, IK)
135 : #if D2_XX_ENABLED
136 69997 : call setDisEuclid(distance(ipnt), point(1:ndim, ipnt), method)
137 : #elif D2_D1_ENABLED
138 67249 : call setDisEuclid(distance(ipnt), point(1:ndim, ipnt), ref, method)
139 : #else
140 : #error "Unrecognized interface."
141 : #endif
142 : end do
143 :
144 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145 : #elif setDisEuclid_ENABLED && D2_D2_ENABLED
146 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147 :
148 : integer(IK) :: ipnt, ndim, nref
149 16682 : nref = size(ref, 2, IK)
150 16682 : ndim = size(point, 1, IK)
151 250230 : CHECK_ASSERTION(__LINE__, all(shape(distance, IK) == [nref, size(point, 2, IK)]), SK_"@getDisEuclid(): The condition `all(shape(distance) == [size(ref, 2), size(point, 2)])` must hold. shape(distance), shape(ref), shape(point) = "//getStr([shape(distance, IK), shape(ref, IK), shape(point, IK)]))
152 192451 : do ipnt = 1_IK, size(point, 2, IK)
153 192451 : call setDisEuclid(distance(1 : nref, ipnt), point(1 : ndim, ipnt), ref, method)
154 : end do
155 :
156 : !%%%%%%%%%%%%%%%%%%%%%%
157 : #elif getDisMatEuclid_ENABLED
158 : !%%%%%%%%%%%%%%%%%%%%%%
159 :
160 : #if FUL_ENABLED
161 : type(rdpack_type), parameter :: pack = rdpack_type()
162 : type(uppLowDia_type), parameter :: subset = uppLowDia_type()
163 : #endif
164 9660 : if (present(method)) then
165 : select type (method)
166 : type is (euclidsq_type)
167 3215 : call setDisMatEuclid(distance, pack, subset, point, method)
168 : type is (euclidu_type)
169 3205 : call setDisMatEuclid(distance, pack, subset, point, method)
170 : type is (euclid_type)
171 3210 : call setDisMatEuclid(distance, pack, subset, point, method)
172 : class default
173 : error stop MODULE_NAME//SK_"@getDisMatEuclid(): Unsupported input value for `method`." ! LCOV_EXCL_LINE
174 : end select
175 : return
176 : end if
177 30 : call setDisMatEuclid(distance, pack, subset, point, euclid)
178 :
179 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
180 : #elif setDisMatEuclid_ENABLED && RDP_ENABLED && ULD_ENABLED
181 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
182 :
183 : ! Construct the full distance matrix.
184 : integer(IK) :: ndim, npnt
185 : integer(IK) :: ipnt, jpnt
186 10845 : ndim = size(point, 1, IK)
187 10845 : npnt = size(point, 2, IK)
188 119295 : CHECK_ASSERTION(__LINE__, all(shape(distance, IK) == [npnt, npnt]), \
189 : SK_"@setDisMatEuclid(): The condition `all(shape(distance) == shape(point))` must hold. shape(distance), shape(point) = "//getStr([shape(distance, IK), shape(point, IK)]))
190 125407 : do jpnt = 1_IK, npnt
191 114562 : distance(jpnt, jpnt) = 0._TKC
192 852036 : do ipnt = jpnt + 1, npnt
193 : ! construct the lower triangle element and transpose to upper triangle.
194 726629 : call setDisEuclid(distance(ipnt, jpnt), point(1 : ndim, ipnt), point(1 : ndim, jpnt), method)
195 841191 : distance(jpnt, ipnt) = distance(ipnt, jpnt)
196 : end do
197 : end do
198 :
199 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
200 : #elif setDisMatEuclid_ENABLED && RDP_ENABLED && ULX_ENABLED
201 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
202 :
203 : ! Construct the full distance matrix excluding diagonal.
204 : integer(IK), parameter :: doff = 1_IK ! diagonal offset.
205 : integer(IK) :: ndim, npnt
206 : integer(IK) :: ipnt, jpnt
207 10845 : ndim = size(point, 1, IK)
208 10845 : npnt = size(point, 2, IK)
209 119295 : CHECK_ASSERTION(__LINE__, all(shape(distance, IK) == [npnt - 1_IK, npnt]), \
210 : SK_"@setDisMatEuclid(): The condition `all(shape(distance) == [size(point, 2) - 1_IK, size(point, 2)])` must hold. shape(distance), shape(point) = "//\
211 : getStr([shape(distance, IK), shape(point, IK)]))
212 125407 : do jpnt = 1_IK, npnt
213 852036 : do ipnt = jpnt + 1, npnt
214 : ! construct the lower triangle element and transpose to upper triangle.
215 726629 : call setDisEuclid(distance(ipnt - doff, jpnt), point(1 : ndim, ipnt), point(1 : ndim, jpnt), method)
216 841191 : distance(jpnt, ipnt) = distance(ipnt - doff, jpnt)
217 : end do
218 : end do
219 : #else
220 : !%%%%%%%%%%%%%%%%%%%%%%%%
221 : #error "Unrecognized interface."
222 : !%%%%%%%%%%%%%%%%%%%%%%%%
223 : #endif
224 : ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
225 : !#elif getDisEuclid_ENABLED && (PNT_ENABLED || REF_ENABLED)
226 : ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
227 : !
228 : ! integer(IK) :: i
229 : ! real(TKC) :: invmax
230 : ! real(TKC) :: maximum
231 : !#if REF_ENABLED
232 : ! real(TKC) :: diff(size(point, 1, IK))
233 : ! CHECK_ASSERTION(__LINE__, size(point, 1, IK) == size(ref, 1, IK), SK_"@getDisEuclid(): The condition `size(point) == size(ref)` must hold. size(point), size(ref) = "//getStr([size(point, 1, IK) == size(ref, 1, IK)]))
234 : !#else
235 : !#define DIFF point
236 : !#endif
237 : ! maximum = -huge(maximum)
238 : ! do i = 1_IK, size(point, 1, IK)
239 : !#if PNT_ENABLED
240 : ! invmax = abs(point(i))
241 : ! if (maximum < invmax) maximum = invmax
242 : !#elif REF_ENABLED
243 : ! diff(i) = abs(point(i) - ref(i))
244 : ! if (maximum < diff(i)) maximum = diff(i)
245 : !#else
246 : !#error "Unrecognized interface."
247 : !#endif
248 : ! end do
249 : ! if(maximum == 0._TKC .or. maximum > huge(maximum)) then
250 : ! distance = sum(point) ! Ensure propagation of potential nan in the vector.
251 : ! else
252 : ! distance = 0._TKC
253 : ! invmax = 1._TKC / maximum
254 : ! do i = 1_IK, size(point, 1, IK)
255 : ! distance = distance + (DIFF(i) * invmax)**2
256 : ! end do
257 : ! distance = maximum * sqrt(distance)
258 : ! end if
259 : !
260 : ! !%%%%%%%%%%%%%%%%%%%%%%%
261 : !#elif getDisMatEuclid_ENABLED
262 : ! !%%%%%%%%%%%%%%%%%%%%%%%
263 : !
264 : !!#if LFP_ENABLED
265 : !!#define OPTIONAL_ARG
266 : !!#elif Sym_ENABLED
267 : !!#define OPTIONAL_ARG diag,
268 : !!#elif REF_ENABLED
269 : !!#define OPTIONAL_ARG ref,
270 : !!#else
271 : !!#error "Unrecognized interface."
272 : !!#endif
273 : ! if (present(method)) then
274 : ! select type (method)
275 : ! type is (euclidsq_type)
276 : ! call setDisMatEuclid(distance, point, OPTIONAL_ARG euclidsq)
277 : ! return
278 : ! type is (euclidu_type)
279 : ! call setDisMatEuclid(distance, point, OPTIONAL_ARG euclidu)
280 : ! return
281 : ! call setDisMatEuclid(distance, point, OPTIONAL_ARG euclid)
282 : ! end select
283 : ! end if
284 : ! ! type is (euclid_type)
285 : ! call setDisMatEuclid(distance, point, OPTIONAL_ARG euclidu)
286 : !#endif
287 : !
288 : ! !%%%%%%%%%%%%%%%%%%%%%%%
289 : !#elif setDisMatEuclid_ENABLED
290 : ! !%%%%%%%%%%%%%%%%%%%%%%%
291 : !
292 : ! integer(IK) :: idim
293 : !
294 : ! !%%%%%%%%%%%%%%%%%%%%%%%
295 : !#elif setDisMatEuclid_ENABLED
296 : ! !%%%%%%%%%%%%%%%%%%%%%%%
297 : !
298 : ! integer(IK) :: ndim, npnt
299 : ! integer(IK) :: ipnt, jref
300 : !#if MUF_ENABLED || MEQ_ENABLED
301 : ! integer(IK) :: idim
302 : !#endif
303 : !#if LFP_ENABLED
304 : ! integer(IK) :: indx
305 : ! indx = 0_IK
306 : !#endif
307 : ! ndim = size(point, 1, IK)
308 : ! npnt = size(point, 2, IK)
309 : !#if DEF_ENABLED || ULD_ENABLED
310 : !#define REFERENCE point
311 : !#define LBND jref + 1_IK
312 : !#define NREF npnt
313 : !#elif REF_ENABLED
314 : !#define NREF size(ref, 2, IK)
315 : !#define REFERENCE ref
316 : !#define LBND 1_IK
317 : ! CHECK_ASSERTION(__LINE__, ndim == size(ref, 1, IK), SK_"@setDisMatEuclid(): The condition `size(point, 1) == size(ref, 1)` must hold. size(point, 1), size(ref, 1) = "//getStr([ndim, size(ref, 1, IK)]))
318 : !#else
319 : !#error "Unrecognized interface."
320 : !#endif
321 : ! ! \warning
322 : ! ! This loops is well-reasoned. Do not change it blindly.
323 : ! do jref = 1_IK, NREF
324 : !#if Sym_ENABLED
325 : ! distance(jref, jref) = diag
326 : !#endif
327 : ! do ipnt = LBND, npnt
328 : !#if LFP_ENABLED
329 : ! indx = indx + 1_IK
330 : !#define INDEX indx
331 : !#else
332 : !#define INDEX ipnt, jref
333 : !#endif
334 : !#if MUF_ENABLED || MEQ_ENABLED
335 : ! distance(INDEX) = 0._TKC
336 : ! do idim = 1_IK, ndim
337 : ! distance(INDEX) = distance(INDEX) + (point(idim, ipnt) - REFERENCE(idim, jref))**2
338 : ! end do
339 : !#if MUF_ENABLED
340 : ! distance(INDEX) = sqrt(distance(INDEX))
341 : !#endif
342 : !#elif S_ENABLED
343 : ! distance(INDEX) = getDisEuclid(point(1:ndim, ipnt), REFERENCE(1:ndim, jref))
344 : !#else
345 : !#error "Unrecognized interface."
346 : !#endif
347 : !#if Sym_ENABLED
348 : ! distance(jref, ipnt) = distance(ipnt, jref)
349 : !#endif
350 : ! end do
351 : ! end do
352 : !
353 : !#else
354 : ! !%%%%%%%%%%%%%%%%%%%%%%%%
355 : !#error "Unrecognized interface."
356 : ! !%%%%%%%%%%%%%%%%%%%%%%%%
357 : !#endif
358 : !#undef OPTIONAL_ARG
359 : !#undef REFERENCE
360 : !#undef INDEX
361 : !#undef NREF
362 : !#undef LBND
|
