ParaMonte Fortran 2.0.0
Parallel Monte Carlo and Machine Learning Library
See the latest version documentation.
pm_clusDensity.F90
Go to the documentation of this file.
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
50module pm_clusDensity
51
52#if CHECK_ENABLED
53 use pm_err, only: getFine
54 use pm_val2str, only: getStr
55 use pm_err, only: setAsserted
56#define CHECK_ASSERTION(LINE,ASSERTION,MSG) \
57call setAsserted(ASSERTION,getFine(__FILE__,LINE)//MODULE_NAME//MSG);
58#else
59#define CHECK_ASSERTION(LINE,ASSERTION,MSG) continue;
60#endif
61
62 use pm_err, only: err_type
63 use pm_kind, only: SK, IK, LK, RK
64 use iso_fortran_env, only: output_unit
65
66 implicit none
67
68#define MAXDEN_ENABLED 1
69#define MINVOL_ENABLED 0
70 character(*, SK), parameter :: MODULE_NAME = "@pm_clusDensity"
71 real(RK) , parameter :: POSINF = +huge(0._RK) / 10
72 real(RK) , parameter :: NEGINF = -huge(0._RK) / 10
73
74 type :: KvolTry_type
75 integer(IK) , allocatable :: EffectiveSize(:)
76 integer(IK) , allocatable :: Membership(:)
77 integer(IK) , allocatable :: CumSumSize(:)
78 integer(IK) , allocatable :: Size(:)
79 real(RK) , allocatable :: Center(:,:)
80 real(RK) , allocatable :: potential(:)
81 real(RK) , allocatable :: MahalSq(:,:)
82 real(RK) , allocatable :: InvCovMat(:,:,:)
83 real(RK) , allocatable :: ChoLowCovUpp(:,:,:)
84 real(RK) , allocatable :: ScaleFactorSq(:)
85 real(RK) , allocatable :: ScaleFactor(:)
86 real(RK) , allocatable :: LogVolRatio(:)
87 real(RK) , allocatable :: LogVolNormed(:)
88 real(RK) , allocatable :: LogLikeFitness(:)
89 end type
90
91 type, extends(KvolTry_type) :: ParTry_type
92 integer(IK) :: ic
93 integer(IK) :: nc
94 integer(IK) :: ne
95 integer(IK) :: np
96 integer(IK) :: nps
97 integer(IK) :: npe
98 integer(IK) :: itp ! parent partition level
99 integer(IK) :: its ! subpartition level
100 logical(LK) :: deallocated = .true._LK
101 real(RK) , allocatable :: SumLogLikeFitness(:)
102 logical(LK) , allocatable :: SubclusteringWarranted(:)
103 end type ParTry_type
104
105 type KmeansTry_type
106 real(RK) , allocatable :: potential(:)
107 real(RK) , allocatable :: Center(:,:)
108 integer(IK) , allocatable :: Membership(:)
109 integer(IK) , allocatable :: Size(:)
110 type(err_type) :: err
111 end type
112
115 type :: partition_type
116 integer(IK) :: nsim
117 integer(IK) :: nd,np
118 integer(IK) :: nc,nt
119 integer(IK) :: nemax
120 integer(IK) :: minSize
121 integer(IK) :: kmeansNumFailMax
122 integer(IK) :: kmeansNumRecursionMax
123 integer(IK) :: kvolumeNumRecursionMax
124 logical(LK) :: scaleOptimizationEnabled
125 logical(LK) :: shapeOptimizationEnabled
126 logical(LK) :: shapeAdaptationEnabled
127 logical(LK) :: biasCorrectionEnabled
128 logical(LK) :: mahalSqWeightEnabled
129 logical(LK) :: optimizationEnabled
130 real(RK) :: mahalSqWeightExponent
131 real(RK) :: parentLogVolNormed
132 real(RK) :: pointLogVolNormed
133#if MINVOL_ENABLED
134 real(RK) :: inclusionFraction
135#endif
136 !real(RK) :: scaleFactor !< See the interface of [partition_typer()](@ref constructpartition).
137 !real(RK) :: scaleFactorSq !< See the interface of [partition_typer()](@ref constructpartition).
138 real(RK) :: logExpansion
139 real(RK) :: logShrinkage
140 real(RK) :: kmeansRelTol
141 real(RK) :: nplog
142 integer(IK) :: ne
143 integer(IK) :: numRecursiveCall
144 integer(IK) :: convergenceFailureCount
145 logical(LK) :: stanEnabled
146 type(ParTry_type), allocatable :: Try(:)
147 integer(IK) , allocatable :: Size(:)
148 integer(IK) , allocatable :: Final(:,:)
149 integer(IK) , allocatable :: CumSumSize(:)
150 integer(IK) , allocatable :: Membership(:)
151 integer(IK) , allocatable :: PointIndex(:)
152 integer(IK) , allocatable :: EffectiveSize(:)
153 real(RK) , allocatable :: LogLikeFitness(:)
154 real(RK) , allocatable :: LogVolNormed(:)
155 real(RK) , allocatable :: ChoLowCovUpp(:,:,:)
156 real(RK) , allocatable :: InvCovMat(:,:,:)
157 real(RK) , allocatable :: Center(:,:)
158 real(RK) , allocatable :: VolNormedCDF(:)
159 real(RK) , allocatable :: BiasCorrectionScaleFactorSq(:)
160 type(err_type) :: err
161 contains
162 procedure, pass :: write => writePartition
163 procedure, pass :: run => runPartition
164 end type
165
166 interface partition_type
167 module procedure :: partition_typer
168 end interface
169
170!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
171
172contains
173
174!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
175
220 function partition_typer( Point & ! LCOV_EXCL_LINE
221 , nc,nt & ! LCOV_EXCL_LINE
222 , nsim & ! LCOV_EXCL_LINE
223 , nemax & ! LCOV_EXCL_LINE
224 , minSize & ! LCOV_EXCL_LINE
225 , stanEnabled & ! LCOV_EXCL_LINE
226 , trimEnabled & ! LCOV_EXCL_LINE
227 , kmeansRelTol & ! LCOV_EXCL_LINE
228 , logExpansion & ! LCOV_EXCL_LINE
229 , logShrinkage & ! LCOV_EXCL_LINE
230#if MINVOL_ENABLED
231 , inclusionFraction & ! LCOV_EXCL_LINE
232#endif
233 , parentLogVolNormed & ! LCOV_EXCL_LINE
234 , mahalSqWeightExponent & ! LCOV_EXCL_LINE
235 , kmeansNumFailMax & ! LCOV_EXCL_LINE
236 , kmeansNumRecursionMax & ! LCOV_EXCL_LINE
237 , kvolumeNumRecursionMax & ! LCOV_EXCL_LINE
238 , scaleOptimizationEnabled & ! LCOV_EXCL_LINE
239 , shapeOptimizationEnabled & ! LCOV_EXCL_LINE
240 , shapeAdaptationEnabled & ! LCOV_EXCL_LINE
241 , biasCorrectionEnabled & ! LCOV_EXCL_LINE
242 ) result(Partition)
243#if INTEL_COMPILER_ENABLED && DLL_ENABLED && (WINDOWS_ENABLED || DARWIN_ENABLED)
244 !DEC$ ATTRIBUTES DLLEXPORT :: partition_typer
245#endif
246 real(RK) , intent(inout), contiguous :: Point(:,:)
247 integer(IK) , intent(in) , optional :: nc,nt
248 integer(IK) , intent(in) , optional :: nsim
249 integer(IK) , intent(in) , optional :: nemax
250 integer(IK) , intent(in) , optional :: minSize
251 logical(LK) , intent(in) , optional :: stanEnabled
252 logical(LK) , intent(in) , optional :: trimEnabled
253 real(RK) , intent(in) , optional :: kmeansRelTol
254 real(RK) , intent(in) , optional :: logExpansion
255 real(RK) , intent(in) , optional :: logShrinkage
256#if MINVOL_ENABLED
257 real(RK) , intent(in) , optional :: inclusionFraction
258#endif
259 real(RK) , intent(in) , optional :: parentLogVolNormed
260 real(RK) , intent(in) , optional :: mahalSqWeightExponent
261 integer(IK) , intent(in) , optional :: kmeansNumFailMax
262 integer(IK) , intent(in) , optional :: kmeansNumRecursionMax
263 integer(IK) , intent(in) , optional :: kvolumeNumRecursionMax
264 logical(LK) , intent(in) , optional :: scaleOptimizationEnabled
265 logical(LK) , intent(in) , optional :: shapeOptimizationEnabled
266 logical(LK) , intent(in) , optional :: shapeAdaptationEnabled
267 logical(LK) , intent(in) , optional :: biasCorrectionEnabled
268 type(partition_type) :: Partition
269
270 Partition%nd = size(Point(:,1))
271 Partition%np = size(Point(1,:))
272
273 Partition%nc = 2; if (present(nc)) Partition%nc = max(Partition%nc, nc)
274 Partition%nt = 1; if (present(nt)) Partition%nt = max(Partition%nt, nt)
275 Partition%nsim = 0; if (present(nsim)) Partition%nsim = max(Partition%nsim, nsim)
276 Partition%minSize = 0; if (present(minSize)) Partition%minSize = minSize ! Partition%nd + 1
277 Partition%nemax = Partition%np / max(1, Partition%minSize); if (present(nemax)) Partition%nemax = max(1, nemax)
278 if (Partition%nemax > Partition%np) Partition%nemax = Partition%np
279 Partition%kvolumeNumRecursionMax = 3; if (present(kvolumeNumRecursionMax)) Partition%kvolumeNumRecursionMax = kvolumeNumRecursionMax
280 Partition%kmeansNumRecursionMax = 300; if (present(kmeansNumRecursionMax)) Partition%kmeansNumRecursionMax = kmeansNumRecursionMax
281 Partition%kmeansNumFailMax = 10; if (present(kmeansNumFailMax)) Partition%kmeansNumFailMax = kmeansNumFailMax
282 Partition%mahalSqWeightExponent = 0._RK; if (present(mahalSqWeightExponent)) Partition%mahalSqWeightExponent = mahalSqWeightExponent
283 Partition%mahalSqWeightEnabled = abs(Partition%mahalSqWeightExponent) > 1.e-4_RK
284#if MINVOL_ENABLED
285 Partition%inclusionFraction = 0._RK; if (present(inclusionFraction)) then; if (abs(inclusionFraction)>1.e-5_RK) Partition%inclusionFraction = inclusionFraction; endif
286#endif
287 !Partition%scaleFactor = 1._RK; if (present(expansionFac)) Partition%scaleFactor = expansionFac**(1._RK/real(Partition%nd,RK))
288 !Partition%scaleFactorSq = Partition%scaleFactor**2_IK
289 Partition%logShrinkage = 0._RK; if (present(logExpansion)) Partition%logExpansion = logExpansion
290 Partition%logShrinkage = 0._RK; if (present(logShrinkage)) Partition%logShrinkage = logShrinkage
291 Partition%kmeansRelTol = 1.e-4_RK; if (present(kmeansRelTol)) Partition%kmeansRelTol = kmeansRelTol
292 Partition%stanEnabled = .true.; if (present(stanEnabled)) Partition%stanEnabled = stanEnabled
293 Partition%nplog = log(real(Partition%np,RK))
294
295 Partition%scaleOptimizationEnabled = .true.; if (present(scaleOptimizationEnabled)) Partition%scaleOptimizationEnabled = scaleOptimizationEnabled
296 Partition%shapeOptimizationEnabled = .true.; if (present(shapeOptimizationEnabled)) Partition%shapeOptimizationEnabled = shapeOptimizationEnabled
297 Partition%shapeAdaptationEnabled = .true.; if (present(shapeAdaptationEnabled)) Partition%shapeAdaptationEnabled = shapeAdaptationEnabled
298 Partition%biasCorrectionEnabled = .true.; if (present(biasCorrectionEnabled)) Partition%biasCorrectionEnabled = biasCorrectionEnabled
299 Partition%optimizationEnabled = Partition%shapeAdaptationEnabled .or. Partition%scaleOptimizationEnabled .or. Partition%shapeOptimizationEnabled
300 !if (Partition%optimizationEnabled) Partition%logShrinkage = NEGINF
301
302 if (Partition%mahalSqWeightExponent == 0._RK .and. .not. present(parentLogVolNormed)) then ! \todo: could we instead approximate `parentLogVolNormed`?
303 Partition%err%msg = "If `mahalSqWeightExponent` is given as input argument, then `parentLogVolNormed` must also be given as input."
304 Partition%err%occurred = .true._LK
305 error stop
306 end if
307
308 allocate( Partition%Try(0:Partition%nemax) & ! LCOV_EXCL_LINE
309 , Partition%Center(Partition%nd, Partition%nemax) & ! LCOV_EXCL_LINE
310 , Partition%InvCovMat(Partition%nd, Partition%nd + 1, Partition%nemax) & ! LCOV_EXCL_LINE
311 , Partition%ChoLowCovUpp(Partition%nd, Partition%nd + 1, Partition%nemax) & ! LCOV_EXCL_LINE
312 , Partition%LogVolNormed(1 : Partition%nemax) & ! LCOV_EXCL_LINE
313 , Partition%LogLikeFitness(1 : Partition%nemax) & ! LCOV_EXCL_LINE
314 , Partition%EffectiveSize(1 : Partition%nemax) & ! LCOV_EXCL_LINE
315 , Partition%VolNormedCDF(1 : Partition%nemax) & ! LCOV_EXCL_LINE
316 , Partition%CumSumSize(0:Partition%nemax) & ! LCOV_EXCL_LINE
317 , Partition%Final (2,Partition%nemax) & ! LCOV_EXCL_LINE
318 , Partition%size (1 : Partition%nemax) & ! LCOV_EXCL_LINE
319 , Partition%Membership(Partition%np) & ! LCOV_EXCL_LINE
320 , Partition%PointIndex(Partition%np) & ! LCOV_EXCL_LINE
321 , Partition%BiasCorrectionScaleFactorSq(Partition%nd+1 : Partition%np) & ! LCOV_EXCL_LINE
322 )
323
324 Partition%BiasCorrectionScaleFactorSq(Partition%nd+1 : Partition%np) = getBiasCorrectionScaleFactorSq(nd = Partition%nd, npmin = Partition%nd + 1, npmax = Partition%np)
325
326 Partition%CumSumSize(0) = 0 ! this must never change.
327 call Partition%run(Point, parentLogVolNormed)
328
329 if (present(trimEnabled)) then
330 if (trimEnabled) then
331 Partition%Try = Partition%Try(0 : Partition%ne)
332 Partition%Center = Partition%Center(1 : Partition%nd, 1 : Partition%ne)
333 Partition%InvCovMat = Partition%InvCovMat(1 : Partition%nd, 1 : Partition%nd , 1 : Partition%ne)
334 Partition%ChoLowCovUpp = Partition%ChoLowCovUpp(1 : Partition%nd, 1 : Partition%nd + 1 , 1 : Partition%ne)
335 Partition%LogVolNormed = Partition%LogVolNormed(1 : Partition%ne)
336 Partition%LogLikeFitness = Partition%LogLikeFitness(1 : Partition%ne)
337 Partition%EffectiveSize = Partition%EffectiveSize(1 : Partition%ne)
338 Partition%CumSumSize = Partition%CumSumSize(0 : Partition%ne)
339 Partition%Final = Partition%Final (1 : 2, 1 : Partition%ne)
340 Partition%size = Partition%size (1 : Partition%ne)
341 end if
342 end if
343
344 end function partition_typer
345
346!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
347
360 subroutine runPartition(Partition, Point, parentLogVolNormed)
361#if INTEL_COMPILER_ENABLED && DLL_ENABLED && (WINDOWS_ENABLED || DARWIN_ENABLED)
362 !DEC$ ATTRIBUTES DLLEXPORT :: runPartition
363#endif
364 use pm_matrixChol, only: setMatChol, lowDia, transHerm
365 use pm_mathCumPropExp, only: setCumPropExp, sequence
366 use pm_mathLogSumExp, only: getLogSumExp
367 use pm_clusKmeans, only: setKmeans, rngf
368 use pm_matrixInv, only: setMatInv, choLow
369 use pm_sampleCov, only: setCov, uppDia
370 use pm_arraySort, only: setSorted
371 use pm_val2str, only: getStr
372
373 implicit none
374
375 class(partition_type) , intent(inout) :: Partition
376 real(RK) , intent(inout) , contiguous :: Point(:,:) ! Partition%nd, Partition%np - must be after Partition declaration.
377 real(RK) , intent(in) , optional :: parentLogVolNormed
378
379 character(*, SK), parameter :: PROCEDURE_NAME = MODULE_NAME//SK_"@runPartition()"
380
381 real(RK) :: ndHalf
382 real(RK) :: unifrnd
383 real(RK) :: ndInverse
384 real(RK) :: expanFactor
385 real(RK) :: scaleFactor
386 real(RK) :: scaleFactorSq
387 real(RK) :: ndHalfInverse
388 real(RK) :: sumPotentialNew
389 real(RK) :: sumPotentialOld
390 real(RK) :: logLikeFitnessGain
391 real(RK) :: scaleFactorSqInverse
392 real(RK) :: PointStd(Partition%nd)
393 real(RK) :: PointStdInv(Partition%nd)
394 real(RK) :: MinDistanceSq(Partition%np)
395 real(RK) :: CumSumMinDistanceSq(Partition%np)
396 real(RK) :: NormedPoint(Partition%nd, Partition%np)
397
398 integer(IK) :: itp
399 integer(IK) :: itmax
400 integer(IK) :: nemax
401 integer(IK) :: ndPlusOne
402 integer(IK) :: minSizePartition
403 integer(IK) :: biasedEffectiveSize
404 integer(IK) :: PointIndex(Partition%np)
405 integer(IK) :: KmeansMemberCounter(Partition%nc)
406 integer(IK) :: i,j, ip, ic,ie,nd,kt,it,ktmin
407 integer(IK) :: ipstart,ipend
408 integer(IK) :: niter
409
410 type(KmeansTry_type) :: KmeansTry(Partition%nt)
411
412 ! kvolume variables
413
414 integer(IK) :: kvnew,kvopt,kvdum
415 integer(IK) :: kvolRecursionCounter
416 integer(IK) :: KvolPointIndex(Partition%np)
417 logical(LK) :: reassignmentOccurred
418 real(RK) :: MahalSqWeight(Partition%nc)
419 real(RK) :: logVolNormedExpected
420 real(RK) :: LogSumVolNormed(2)
421 real(RK) :: maxLogVolNormed
422 type(KvolTry_type) :: KvolTry(2)
423
424 ! optimization variables
425
426 real(RK) :: MahalSq(Partition%np) ! Partition%np, Partition%ne
427 real(RK) :: logVolRatio
428 real(RK) :: logExpanFactor
429 real(RK) :: logScaleFactorSq
430 real(RK) :: logLikeFitnessRef
431 real(RK) :: logLikeOccurrence
432 real(RK) :: PointNormed(Partition%nd)
433 real(RK), allocatable :: SortedValue(:) ! Partition%nemax
434 integer , allocatable :: SortedIndex(:) ! Partition%nemax
435 integer(IK) :: MahalSqSortedIndex(Partition%np)
436 integer(IK) :: optimizationCounter
437 integer(IK) :: npc,nps,npe,je,jes!, counter
438 logical(LK) :: optimizationWarranted
439 logical(LK) :: shapeAdaptationEnabled
440 logical(LK) :: subclusteringWarranted
441 logical(LK) :: scaleOptimizationSucceeded
442 logical(LK) :: shapeOptimizationSucceeded
443
444#if MINVOL_ENABLED
445 Partition%mahalSqWeightEnabled = mahalSqWeightExponent /= 0._RK .and. parentLogVolNormed > NEGINF ! \todo: why is the latter needed? what have you done to humanity, alzheimer...
446#endif
447 nd = Partition%nd
448 ndPlusOne = nd + 1_IK
449 ndHalf = 0.5_RK * nd
450 ndInverse = 1._RK / nd
451 ndHalfInverse = 1._RK / ndHalf
452
453 minSizePartition = max(1,Partition%minSize)
454
455 Partition%err%occurred = .false._LK
456 Partition%numRecursiveCall = 0_IK
457 Partition%convergenceFailureCount = 0_IK
458 Partition%Membership(1 : Partition%np) = 0_IK ! \todo: is this needed?
459 if (present(parentLogVolNormed)) Partition%parentLogVolNormed = parentLogVolNormed + Partition%logExpansion
460 Partition%pointLogVolNormed = Partition%parentLogVolNormed - Partition%nplog
461 !partitionScalingEnabled = abs(Partition%scaleFactor - 1._RK) < 1.e-6_RK
462
463 do concurrent(ip = 1_IK:Partition%np)
464 Partition%PointIndex(ip) = ip ! \todo: this component must be removed
465 end do
466
467 if (Partition%nt > 1_IK) then
468 do kt = 1_IK, Partition%nt
469 allocate( KmeansTry(kt)%size(Partition%nc) & ! LCOV_EXCL_LINE
470 , KmeansTry(kt)%potential(Partition%nc) & ! LCOV_EXCL_LINE
471 , KmeansTry(kt)%Center(nd, Partition%nc) & ! LCOV_EXCL_LINE
472 , KmeansTry(kt)%Membership(Partition%np) & ! LCOV_EXCL_LINE
473 )
474 end do
475 end if
476
477 if (Partition%kvolumeNumRecursionMax > 0_IK) then ! .or. Partition%biasCorrectionEnabled) then ! \todo The size must be corrected in case of `Partition%biasCorrectionEnabled`.
478 do kt = 1_IK, 2_IK ! size(KvolTry)
479 allocate( KvolTry(kt)%EffectiveSize(Partition%nc) & ! LCOV_EXCL_LINE
480 , KvolTry(kt)%Membership(Partition%np) & ! LCOV_EXCL_LINE
481 , KvolTry(kt)%CumSumSize(0 : Partition%nc) & ! LCOV_EXCL_LINE
482 , KvolTry(kt)%size(Partition%nc) & ! LCOV_EXCL_LINE
483 , KvolTry(kt)%Center(nd, Partition%nc) & ! LCOV_EXCL_LINE
484 , KvolTry(kt)%MahalSq(Partition%np,Partition%nc) & ! LCOV_EXCL_LINE
485 , KvolTry(kt)%InvCovMat(nd, nd, Partition%nc) & ! LCOV_EXCL_LINE
486 , KvolTry(kt)%ChoLowCovUpp(nd, nd + 1, Partition%nc) & ! LCOV_EXCL_LINE
487 , KvolTry(kt)%ScaleFactorSq(Partition%nc) & ! LCOV_EXCL_LINE
488 , KvolTry(kt)%ScaleFactor(Partition%nc) & ! LCOV_EXCL_LINE
489 , KvolTry(kt)%LogVolRatio(Partition%nc) & ! LCOV_EXCL_LINE
490 , KvolTry(kt)%LogVolNormed(Partition%nc) & ! LCOV_EXCL_LINE
491 , KvolTry(kt)%LogLikeFitness(Partition%nc) & ! LCOV_EXCL_LINE
492 )
493 end do
494 end if
495
496 it = 0_IK
497 itmax = 1_IK
498 Partition%ne = 0_IK
499 Partition%Try(0_IK)%ic = 0_IK
500 Partition%Try(0_IK)%nc = 1_IK
501 Partition%Try(0_IK)%ne = 0_IK
502 Partition%Try(0_IK)%itp = 0_IK
503 nemax = Partition%nemax
504
505 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
506
507 loopRecursivePartition: do
508
509 Partition%numRecursiveCall = Partition%numRecursiveCall + 1_IK
510
511 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
512 blockForwardBackwardStream: if (Partition%Try(it)%ic == 0_IK) then !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
513 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
514
515 nemax = nemax - Partition%Try(it)%nc + 1_IK
516 if (Partition%Try(it)%deallocated) then
517 allocate( Partition%Try(it)%size(Partition%nc) & ! LCOV_EXCL_LINE
518 , Partition%Try(it)%CumSumSize(0_IK:Partition%nc) & ! LCOV_EXCL_LINE
519 , Partition%Try(it)%Center(nd, Partition%nc) & ! LCOV_EXCL_LINE
520 , Partition%Try(it)%LogVolRatio(Partition%nc) & ! LCOV_EXCL_LINE
521 , Partition%Try(it)%ScaleFactor(Partition%nc) & ! LCOV_EXCL_LINE
522 , Partition%Try(it)%LogVolNormed(Partition%nc) & ! LCOV_EXCL_LINE
523 , Partition%Try(it)%ScaleFactorSq(Partition%nc) & ! LCOV_EXCL_LINE
524 , Partition%Try(it)%LogLikeFitness(Partition%nc) & ! LCOV_EXCL_LINE
525 , Partition%Try(it)%InvCovMat(nd, nd, Partition%nc) & ! LCOV_EXCL_LINE
526 , Partition%Try(it)%SumLogLikeFitness(Partition%nc) & ! LCOV_EXCL_LINE
527 , Partition%Try(it)%ChoLowCovUpp(nd, nd + 1, Partition%nc) & ! LCOV_EXCL_LINE
528 , Partition%Try(it)%MahalSq(Partition%np,Partition%nc) & ! LCOV_EXCL_LINE
529 , Partition%Try(it)%SubclusteringWarranted(Partition%nc) & ! LCOV_EXCL_LINE
530 , Partition%Try(it)%EffectiveSize(Partition%nc) & ! LCOV_EXCL_LINE
531 , Partition%Try(it)%Membership(Partition%np) & ! LCOV_EXCL_LINE
532 )
533 Partition%Try(it)%deallocated = .false._LK
534 end if
535 Partition%Try(it)%SumLogLikeFitness(1_IK:Partition%Try(it)%nc) = 0._RK ! POSINF
536
537 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
538 ! Perform Kmeans clustering or set the single cluster properties.
539 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
540
541 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
542 blockKmeansOrSingularPartition: if (Partition%Try(it)%nc > 1_IK) then !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
543 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
544
545 Partition%Try(it)%nps = Partition%Try(Partition%Try(it)%itp)%CumSumSize(Partition%Try(Partition%Try(it)%itp)%ic-1) + 1_IK
546 Partition%Try(it)%npe = Partition%Try(Partition%Try(it)%itp)%CumSumSize(Partition%Try(Partition%Try(it)%itp)%ic)
547 Partition%Try(it)%np = Partition%Try(it)%npe - Partition%Try(it)%nps + 1_IK
548
549 blockStandardization: if (Partition%stanEnabled) then
550#define STAN_ENABLED 1
551 if (Partition%kvolumeNumRecursionMax > 0_IK) then ! Store the partitions in the temporary `KvolTry` objects.
552#define KVOL_ENABLED 1
553#include "pm_clusDensity.inc.kmeans.inc.F90"
554#include "pm_clusDensity.inc.prop.inc.F90"
555#undef KVOL_ENABLED
556 else
557#include "pm_clusDensity.inc.kmeans.inc.F90"
558#include "pm_clusDensity.inc.prop.inc.F90"
559 end if
560#undef STAN_ENABLED
561 else blockStandardization
562 if (Partition%kvolumeNumRecursionMax > 0_IK) then ! Store the partitions in the temporary `KvolTry` objects.
563#define KVOL_ENABLED 1
564#include "pm_clusDensity.inc.kmeans.inc.F90"
565#include "pm_clusDensity.inc.prop.inc.F90"
566#undef KVOL_ENABLED
567 else
568#include "pm_clusDensity.inc.kmeans.inc.F90"
569#include "pm_clusDensity.inc.prop.inc.F90"
570 end if
571 end if blockStandardization
572
573 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
574 ! If optimization is enabled, optimize clusters via kvolume.
575 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
576
577 blockRecursiveKvolumeOptimization: if (it > 0_IK .and. Partition%kvolumeNumRecursionMax > 0_IK) then
578
579 kvnew = 2_IK
580 kvolRecursionCounter = 0_IK
581 loopRecursivePartitionOptimization: do
582
583 kvolRecursionCounter = kvolRecursionCounter + 1_IK
584
585 KvolTry(kvnew)%size(1 : Partition%Try(it)%nc) = KvolTry(kvopt)%size(1 : Partition%Try(it)%nc)
586 do concurrent(ic = 1 : Partition%Try(it)%nc)
587 KvolTry(kvnew)%Center(1 : nd, ic) = KvolTry(kvopt)%Center(1 : nd, ic) * KvolTry(kvopt)%size(ic)
588 end do
589
590 ! Determine the cluster to which the point is the closest and switch the point membership if needed.
591
592 blockReassignPoint: if (Partition%mahalSqWeightEnabled) then
593#define MAHALSQ_WEIGHT_ENABLED 1
594#include "pm_clusDensity.inc.kvolume.inc.F90"
595#undef MAHALSQ_WEIGHT_ENABLED
596 else blockReassignPoint
597#include "pm_clusDensity.inc.kvolume.inc.F90"
598 end if blockReassignPoint
599
600 ! Restart the process if anything has changed.
601
602 !if (any(KvolTry(kvnew)%size(1 : Partition%Try(it)%nc) == Partition%Try(it)%np)) write(*,*) "Zero Size Cluster Identified"
603 !write(*,*) "KVolume Size, np:", &
604 !KvolTry(kvnew)%size(1 : Partition%Try(it)%nc), Partition%Try(it)%np, &
605 !KvolTry(kvnew)%LogVolNormed(1 : Partition%Try(it)%nc), &
606 !KvolTry(kvnew)%LogLikeFitness(1 : Partition%Try(it)%nc), sum(KvolTry(kvnew)%LogLikeFitness(1 : Partition%Try(it)%nc))
607
608 blockReclusteringNeeded: if (reassignmentOccurred) then
609
610 ! Check for the existence of more than one viable cluster.
611
612 do ic = 1, Partition%Try(it)%nc
613 if (KvolTry(kvnew)%size(ic) == Partition%Try(it)%np) exit loopRecursivePartitionOptimization
614 end do
615 !if (any(KvolTry(kvnew)%size(1 : Partition%Try(it)%nc) == Partition%Try(it)%np)) exit loopRecursivePartitionOptimization !then
616 ! it = Partition%Try(it)%itp
617 ! Partition%err%occurred = .true._LK
618 ! cycle loopRecursivePartition
619 !end if
620
621 ! Convert sum(Point) to center of clusters.
622
623 do concurrent(ic = 1 : Partition%Try(it)%nc)
624 if (KvolTry(kvnew)%size(ic) > 0_IK) KvolTry(kvnew)%Center(1 : nd, ic) = KvolTry(kvnew)%Center(1 : nd, ic) / KvolTry(kvnew)%size(ic)
625 end do
626
627 ! Reorder Point based on the identified clusters and recompute the cluster properties.
628
629#define KVOL_ENABLED 2
630#include "pm_clusDensity.inc.prop.inc.F90"
631#undef KVOL_ENABLED
632
633 ! Check for the stability of the optimization. ATTN: Any exit after property calculations must reshuffle points.
634
635 if (LogSumVolNormed(kvnew) > LogSumVolNormed(kvopt)) then ! reshuffle Point to match that of kvopt partitioning.
636 Point(1 : nd, Partition%Try(it)%nps : Partition%Try(it)%npe) = Point(1 : nd,KvolPointIndex(Partition%Try(it)%nps : Partition%Try(it)%npe))
637 Partition%PointIndex(Partition%Try(it)%nps : Partition%Try(it)%npe) = Partition%PointIndex(KvolPointIndex(Partition%Try(it)%nps : Partition%Try(it)%npe))
638 exit loopRecursivePartitionOptimization
639 end if
640
641 !do ic = 1, Partition%Try(it)%nc
642 ! if (KvolTry(kvopt)%LogVolNormed(ic) > KvolTry(kvnew)%LogVolNormed(ic)) then
643 ! Point(1 : nd, Partition%Try(it)%nps : Partition%Try(it)%npe) = Point(1 : nd,KvolPointIndex(Partition%Try(it)%nps : Partition%Try(it)%npe))
644 ! Partition%PointIndex(Partition%Try(it)%nps : Partition%Try(it)%npe) = Partition%PointIndex(KvolPointIndex(Partition%Try(it)%nps : Partition%Try(it)%npe))
645 ! exit loopRecursivePartitionOptimization
646 ! end if
647 !end do
648
649 !KvolSumLogLikeFitness(kvolRecursionCounter) = real(sum(KvolTry(kvnew)%LogLikeFitness(1 : Partition%Try(it)%nc)),RK1)
650 !do i = 1, kvolRecursionCounter - 1
651 ! if (KvolSumLogLikeFitness(kvolRecursionCounter) == KvolSumLogLikeFitness(i)) then ! reshuffle Point to the last arrangement and exit
652 ! Point(1 : nd, Partition%Try(it)%nps : Partition%Try(it)%npe) = Point(1 : nd,KvolPointIndex(Partition%Try(it)%nps : Partition%Try(it)%npe))
653 ! Partition%PointIndex(Partition%Try(it)%nps : Partition%Try(it)%npe) = Partition%PointIndex(KvolPointIndex(Partition%Try(it)%nps : Partition%Try(it)%npe))
654 ! exit loopRecursivePartitionOptimization
655 ! end if
656 !end do
657
658 !logLikeFitnessGain = sum(KvolTry(kvnew)%LogLikeFitness) - sum(KvolTry(kvopt)%LogLikeFitness)
659 !if (logLikeFitnessGain < 0._RK) then
660 ! call random_number(unifrnd)
661 ! if (unifrnd > 0._RK) then; unifrnd = log(unifrnd); else; unifrnd = -huge(unifrnd); end if ! LCOV_EXCL_LINE
662 ! logLikeFitnessGain = logLikeFitnessGain - unifrnd
663 !end if
664 !blockOptimalPartitionSelection: if (.true. .or. logLikeFitnessGain > 0._RK) then
665 !blockOptimalPartitionSelection: if (sum(KvolTry(kvnew)%LogLikeFitness) > sum(KvolTry(kvopt)%LogLikeFitness)) then
666 !blockOptimalPartitionSelection: if (sum(KvolTry(kvnew)%LogVolNormed) < sum(KvolTry(kvopt)%LogVolNormed)) then
667
668 kvdum = kvopt
669 kvopt = kvnew
670 kvnew = kvdum
671
672 if (kvolRecursionCounter < Partition%kvolumeNumRecursionMax) cycle loopRecursivePartitionOptimization
673
674 !if (kvolRecursionCounter < Partition%kvolumeNumRecursionMax) then
675 ! ! Replace the larger volume partitions with smaller volume partitions, then cycle.
676 ! do concurrent(ic = 1 : Partition%Try(it)%nc)
677 ! if (KvolTry(kvopt)%LogVolNormed(ic) > KvolTry(kvnew)%LogVolNormed(ic)) then
678 ! KvolTry(kvopt)%EffectiveSize (1 : Partition%Try(it)%nc) = KvolTry(kvnew)%EffectiveSize (1 : Partition%Try(it)%nc)
679 ! KvolTry(kvopt)%Membership (Partition%Try(it)%nps : Partition%Try(it)%npe) = KvolTry(kvnew)%Membership (Partition%Try(it)%nps : Partition%Try(it)%npe)
680 ! KvolTry(kvopt)%CumSumSize (0:Partition%Try(it)%nc) = KvolTry(kvnew)%CumSumSize (0:Partition%Try(it)%nc)
681 ! KvolTry(kvopt)%size (1 : Partition%Try(it)%nc) = KvolTry(kvnew)%size (1 : Partition%Try(it)%nc)
682 ! KvolTry(kvopt)%Center (1 : nd, 1 : Partition%Try(it)%nc) = KvolTry(kvnew)%Center (1 : nd, 1 : Partition%Try(it)%nc)
683 ! KvolTry(kvopt)%ChoDia (1 : nd, 1 : Partition%Try(it)%nc) = KvolTry(kvnew)%ChoDia (1 : nd, 1 : Partition%Try(it)%nc)
684 ! KvolTry(kvopt)%MahalSq (1 : Partition%np, 1 : Partition%Try(it)%nc) = KvolTry(kvnew)%MahalSq (1 : Partition%np, 1 : Partition%Try(it)%nc)
685 ! KvolTry(kvopt)%InvCovMat (1 : nd,1 : nd, 1 : Partition%Try(it)%nc) = KvolTry(kvnew)%InvCovMat (1 : nd,1 : nd, 1 : Partition%Try(it)%nc)
686 ! KvolTry(kvopt)%ChoLowCovUpp (1 : nd,1 : nd + 1, 1 : Partition%Try(it)%nc) = KvolTry(kvnew)%ChoLowCovUpp (1 : nd,1 : nd + 1, 1 : Partition%Try(it)%nc)
687 ! KvolTry(kvopt)%ScaleFactorSq (1 : Partition%Try(it)%nc) = KvolTry(kvnew)%ScaleFactorSq (1 : Partition%Try(it)%nc)
688 ! KvolTry(kvopt)%ScaleFactor (1 : Partition%Try(it)%nc) = KvolTry(kvnew)%ScaleFactor (1 : Partition%Try(it)%nc)
689 ! KvolTry(kvopt)%LogVolRatio (1 : Partition%Try(it)%nc) = KvolTry(kvnew)%LogVolRatio (1 : Partition%Try(it)%nc)
690 ! KvolTry(kvopt)%LogVolNormed (1 : Partition%Try(it)%nc) = KvolTry(kvnew)%LogVolNormed (1 : Partition%Try(it)%nc)
691 ! KvolTry(kvopt)%LogLikeFitness (1 : Partition%Try(it)%nc) = KvolTry(kvnew)%LogLikeFitness (1 : Partition%Try(it)%nc)
692 ! end if
693 ! end do
694 ! cycle loopRecursivePartitionOptimization
695 !end if
696
697 !else blockOptimalPartitionSelection ! reshuffle Point to the last arrangement
698 ! Point(1 : nd, Partition%Try(it)%nps : Partition%Try(it)%npe) = Point(1 : nd,KvolPointIndex(Partition%Try(it)%nps : Partition%Try(it)%npe))
699 ! Partition%PointIndex(Partition%Try(it)%nps : Partition%Try(it)%npe) = Partition%PointIndex(KvolPointIndex(Partition%Try(it)%nps : Partition%Try(it)%npe))
700 !end if blockOptimalPartitionSelection
701
702 end if blockReclusteringNeeded
703
704 exit loopRecursivePartitionOptimization
705
706 end do loopRecursivePartitionOptimization
707
708 !write(*,*) "kvolRecursionCounter", kvolRecursionCounter
709 !if (kvolRecursionCounter == Partition%kvolumeNumRecursionMax) read(*,*)
710
711 ! assign KvolTry results to the current partition try.
712
713 Partition%Try(it)%EffectiveSize(1 : Partition%Try(it)%nc) = KvolTry(kvopt)%EffectiveSize(1 : Partition%Try(it)%nc)
714 Partition%Try(it)%Membership(Partition%Try(it)%nps : Partition%Try(it)%npe) = KvolTry(kvopt)%Membership(Partition%Try(it)%nps : Partition%Try(it)%npe)
715 Partition%Try(it)%CumSumSize(0 : Partition%Try(it)%nc) = KvolTry(kvopt)%CumSumSize(0 : Partition%Try(it)%nc)
716 Partition%Try(it)%size (1 : Partition%Try(it)%nc) = KvolTry(kvopt)%size (1 : Partition%Try(it)%nc)
717 Partition%Try(it)%Center(1 : nd, 1 : Partition%Try(it)%nc) = KvolTry(kvopt)%Center(1 : nd, 1 : Partition%Try(it)%nc)
718 Partition%Try(it)%MahalSq(1 : Partition%np, 1 : Partition%Try(it)%nc) = KvolTry(kvopt)%MahalSq(1 : Partition%np, 1 : Partition%Try(it)%nc)
719 Partition%Try(it)%InvCovMat(1 : nd, 1 : nd, 1 : Partition%Try(it)%nc) = KvolTry(kvopt)%InvCovMat(1 : nd, 1 : nd, 1 : Partition%Try(it)%nc)
720 Partition%Try(it)%ChoLowCovUpp(1 : nd, 1 : nd + 1, 1 : Partition%Try(it)%nc) = KvolTry(kvopt)%ChoLowCovUpp(1 : nd, 1 : nd + 1, 1 : Partition%Try(it)%nc)
721 Partition%Try(it)%ScaleFactorSq(1 : Partition%Try(it)%nc) = KvolTry(kvopt)%ScaleFactorSq(1 : Partition%Try(it)%nc)
722 Partition%Try(it)%ScaleFactor(1 : Partition%Try(it)%nc) = KvolTry(kvopt)%ScaleFactor(1 : Partition%Try(it)%nc)
723 Partition%Try(it)%LogVolRatio(1 : Partition%Try(it)%nc) = KvolTry(kvopt)%LogVolRatio(1 : Partition%Try(it)%nc)
724 Partition%Try(it)%LogVolNormed(1 : Partition%Try(it)%nc) = KvolTry(kvopt)%LogVolNormed(1 : Partition%Try(it)%nc)
725 Partition%Try(it)%LogLikeFitness(1 : Partition%Try(it)%nc) = KvolTry(kvopt)%LogLikeFitness(1 : Partition%Try(it)%nc)
726
727 end if blockRecursiveKvolumeOptimization
728
729 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
730 else blockKmeansOrSingularPartition !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
731 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
732
733 CHECK_ASSERTION(__LINE__, it == 0, SK_": The condition `it == 0` must hold. \todo: this must be fixed. This condition happens when Partition%nemax ~ Partition%nc. it = "//getStr(it))
734 Partition%Try(it)%npe = Partition%np
735 Partition%Try(it)%nps = Partition%CumSumSize(Partition%ne) + 1_IK
736 Partition%Try(it)%np = Partition%Try(it)%npe - Partition%Try(it)%nps + 1_IK
737 Partition%Try(it)%size(1) = Partition%Try(it)%np
738 Partition%Try(it)%Membership(Partition%Try(it)%nps : Partition%Try(it)%npe) = 1_IK ! \todo: is this necessary? YES, in pm_clusDensity.inc.prop.inc.F90
739 Partition%Try(it)%Center(1 : nd,1) = 0._RK
740 do ip = Partition%Try(it)%nps, Partition%Try(it)%npe
741 Partition%Try(it)%Center(1 : nd,1) = Partition%Try(it)%Center(1 : nd,1) + Point(1 : nd, ip)
742 end do
743 Partition%Try(it)%Center(1 : nd,1) = Partition%Try(it)%Center(1 : nd,1) / Partition%Try(it)%np
744
745 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
746 ! Compute cluster properties.
747 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
748
749#include "pm_clusDensity.inc.prop.inc.F90"
750
751 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
752 end if blockKmeansOrSingularPartition !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
753 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
754
755#if CHECK_ENABLED
756 block
757 integer :: sizeDum
758 real(RK) :: CenterDum(nd)
759 do ic = 1, Partition%Try(it)%nc
760 sizeDum = 0
761 CenterDum = 0._RK
762 do ip = Partition%Try(it)%CumSumSize(ic-1) + 1, Partition%Try(it)%CumSumSize(ic)
763 CenterDum(1 : nd) = CenterDum(1 : nd) + Point(1 : nd, ip)
764 sizeDum = sizeDum + 1_IK
765 end do
766 if (sizeDum > 0) then
767 CenterDum(1 : nd) = CenterDum(1 : nd) / sizeDum
768 if (any(abs(CenterDum(1 : nd) - Partition%Try(it)%Center(1 : nd, ic)) > 1.e-6_RK)) then
769 if (sizeDum /= Partition%Try(it)%size(ic)) then; write(*,*) "sizeDum, Kmeans%size = ", sizeDum, Partition%Try(it)%size; error stop; endif
770 write(*,*) "nd,minSize, ic,nc,it = ", nd, Partition%minSize, ic, Partition%Try(it)%nc, it
771 write(*,*) "Partition%Try(it)%Center(:) = ", Partition%Try(it)%Center(1 : nd, ic)
772 write(*,*) "CenterDum(:) = ", CenterDum(1 : nd)
773 write(*,*) "Partition%Try(it)%size(:) = ", Partition%Try(it)%size
774 write(*,*) "Partition%err%stat = ", Partition%err%stat
775 write(*,*) "Partition%err%occurred = ", Partition%err%occurred
776 error stop
777 endif
778 endif
779 end do
780 end block
781#endif
782
783 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
784 elseif (Partition%err%occurred) then blockForwardBackwardStream ! failed subcluster !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
785 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
786
787 !write(*,*) "failed subclustering"
788 nemax = nemax + Partition%Try(Partition%Try(it)%its)%nc - 1_IK
789 Partition%ne = Partition%Try(it)%ne + 1
790 Partition%err%occurred = .false._LK
791 itmax = itmax - 1
792 CHECK_ASSERTION(__LINE__, Partition%Try(it)%size(Partition%Try(it)%ic) > 0, SK_": The condition `Partition%Try(it)%size(Partition%Try(it)%ic) > 0` must hold. it, ic, Partition%Try(it)%ic, Partition%Try(it)%size(Partition%Try(it)%ic) = "//getStr([it, ic, Partition%Try(it)%ic, Partition%Try(it)%size(Partition%Try(it)%ic)]))
793#include "pm_clusDensity.inc.terminal.inc.F90"
794
795 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
796 end if blockForwardBackwardStream !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
797 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
798
799 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
800 ! Compute the likelihood and initiate the next search.
801 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
802
803 loopOverSubclusters: do
804
805 blockForwardBackwardDecision: if (Partition%Try(it)%ic < Partition%Try(it)%nc) then
806
807 Partition%Try(it)%ic = Partition%Try(it)%ic + 1_IK
808
809 blockComputeLikelihoodForNonEmptyCluster: if (Partition%Try(it)%size(Partition%Try(it)%ic) > 0_IK) then
810
811 ! Compute the LogVolume and occurrence likelihood of the cluster based on the input estimated volume.
812
813 !Partition%Try(it)%LogVolRatio(Partition%Try(it)%ic) = Partition%Try(it)%LogVolNormed(Partition%Try(it)%ic) - Partition%pointLogVolNormed - log(real(Partition%Try(it)%EffectiveSize(Partition%Try(it)%ic),RK))
814 !Partition%Try(it)%LogLikeFitness(Partition%Try(it)%ic) = getLogLikeFitness(Partition%Try(it)%EffectiveSize(Partition%Try(it)%ic), Partition%Try(it)%LogVolRatio(Partition%Try(it)%ic))
815
816 Partition%Try(it)%SubclusteringWarranted(Partition%Try(it)%ic) = Partition%Try(it)%LogVolRatio(Partition%Try(it)%ic) > 0._RK .and. Partition%Try(it)%size(Partition%Try(it)%ic) > nd .and. nemax > 1
817 if (Partition%Try(it)%SubclusteringWarranted(Partition%Try(it)%ic)) then
818 call random_number(unifrnd)
819 if (unifrnd > 0._RK) then; unifrnd = log(unifrnd); else; unifrnd = -huge(unifrnd); end if ! LCOV_EXCL_LINE
820#if CHECK_ENABLED
821 if (Partition%Try(it)%LogLikeFitness(Partition%Try(it)%ic) < NEGINF .or. Partition%Try(it)%LogLikeFitness(Partition%Try(it)%ic) > 0._RK) then
822 write(*,*) .or."Internal error occurred: Partition%Try(it)%LogLikeFitness(Partition%Try(it)%ic) < NEGINF Partition%Try(it)%LogLikeFitness(Partition%Try(it)%ic) > 0._RK"
823 write(*,*) "Partition%Try(it)%LogLikeFitness(Partition%Try(it)%ic):", Partition%Try(it)%LogLikeFitness(Partition%Try(it)%ic)
824 write(*,*) "Partition%Try(it)%EffectiveSize(Partition%Try(it)%ic):", Partition%Try(it)%EffectiveSize(Partition%Try(it)%ic)
825 write(*,*) "Partition%Try(it)%LogVolRatio(Partition%Try(it)%ic):", Partition%Try(it)%LogVolRatio(Partition%Try(it)%ic)
826 error stop
827 end if
828#endif
829 Partition%Try(it)%SubclusteringWarranted(Partition%Try(it)%ic) = Partition%Try(it)%LogLikeFitness(Partition%Try(it)%ic) < unifrnd
830 end if
831
832 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
833 ! Perform bias correction before moving forward.
834 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
835
836 if (Partition%Try(it)%size(Partition%Try(it)%ic) > nd) then
837 associate(ic => Partition%Try(it)%ic)
838#include "pm_clusDensity.inc.biasCorrection.inc.F90"
839 end associate
840 end if
841
842 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
843 ! Search for subclusters or do not.
844 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
845
846 if (Partition%Try(it)%SubclusteringWarranted(Partition%Try(it)%ic)) then ! forward
847 itp = it
848 it = itmax
849 itmax = itmax + 1_IK
850 Partition%Try(it)%ic = 0_IK
851 Partition%Try(it)%itp = itp
852 Partition%Try(itp)%its = it
853 Partition%Try(itp)%ne = Partition%ne
854 Partition%Try(it)%nc = min(Partition%nc, nemax)
855 cycle loopRecursivePartition
856 else
857 Partition%ne = Partition%ne + 1_IK
858#include "pm_clusDensity.inc.terminal.inc.F90"
859 cycle loopOverSubclusters
860 end if
861
862 else blockComputeLikelihoodForNonEmptyCluster ! empty cluster
863
864 cycle loopOverSubclusters
865
866 end if blockComputeLikelihoodForNonEmptyCluster
867
868 else blockForwardBackwardDecision
869
870 if (it > 0_IK) then ! backward
871 it = Partition%Try(it)%itp
872 Partition%Try(it)%SumLogLikeFitness(Partition%Try(it)%ic) = sum(Partition%Try(Partition%Try(it)%its)%SumLogLikeFitness(1 : Partition%Try(Partition%Try(it)%its)%nc))
873 if (Partition%Try(it)%LogLikeFitness(Partition%Try(it)%ic) > Partition%Try(it)%SumLogLikeFitness(Partition%Try(it)%ic) + Partition%logShrinkage) then
874 !if (Partition%Try(it)%LogLikeFitness(Partition%Try(it)%ic) + 0.5_RK*nd**2*log(real(Partition%Try(it)%size(Partition%Try(it)%ic),RK)) > Partition%Try(it)%SumLogLikeFitness(Partition%Try(it)%ic) + Partition%logShrinkage) then
875 Partition%Try(Partition%Try(it)%its)%SumLogLikeFitness(1 : Partition%Try(Partition%Try(it)%its)%nc) = 0._RK ! POSINF ! This and the other is likely unnecessary. Remove!
876 nemax = nemax + Partition%Try(Partition%Try(it)%its)%nc - 1_IK
877 Partition%ne = Partition%Try(it)%ne + 1_IK
878 itmax = itmax - 1_IK
879#if CHECK_ENABLED
880 if (Partition%Try(it)%size(Partition%Try(it)%ic) <= 0_IK) error stop
881#endif
882#include "pm_clusDensity.inc.terminal.inc.F90"
883 end if
884 cycle loopRecursivePartition
885 else
886 exit loopRecursivePartition
887 end if
888
889 end if blockForwardBackwardDecision
890
891 end do loopOverSubclusters
892
893 end do loopRecursivePartition
894
895 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
896 ! Optimize the partitions.
897 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
898
899#include "pm_clusDensity.inc.optimize.inc.F90"
900
901 if(Partition%err%occurred) error stop
902
903 end subroutine runPartition
904
905!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
906
907 PURE function getLogLikeFitness(count, logVolRatio) result(logLikeFitness)
908#if INTEL_COMPILER_ENABLED && DLL_ENABLED && (WINDOWS_ENABLED || DARWIN_ENABLED)
909 !DEC$ ATTRIBUTES DLLEXPORT :: getLogLikeFitness
910#endif
911 integer(IK) , intent(in) :: count
912 real(RK) , intent(in) :: logVolRatio
913 real(RK) :: logLikeFitness
914 logLikeFitness = count * (1._RK + logVolRatio - exp(logVolRatio))
915#if CHECK_ENABLED
916 if (logLikeFitness < NEGINF .or. logLikeFitness > 0._RK) then
917 write(*,*) .or."Internal error occurred: logLikeFitness < NEGINF logLikeFitness > 0._RK"
918 write(*,*) "count, logVolRatio, logLikeFitness:", count, logVolRatio, logLikeFitness
919 error stop
920 end if
921#endif
922 end function getLogLikeFitness
923
924!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
925
928 function getBiasCorrectionScaleFactorSq(nd,npmin,npmax) result(BiasCorrectionScaleFactorSq)
929#if INTEL_COMPILER_ENABLED && DLL_ENABLED && (WINDOWS_ENABLED || DARWIN_ENABLED)
930 !DEC$ ATTRIBUTES DLLEXPORT :: getBiasCorrectionScaleFactorSq
931#endif
932 integer(IK) , intent(in) :: nd
933 integer(IK) , intent(in) :: npmin
934 integer(IK) , intent(in) :: npmax
935 real(RK) :: BiasCorrectionScaleFactorSq(npmin:npmax)
936 integer(IK) :: ip, ndim !, npminDef, icase
937 real(RK) :: LogCount(npmin:npmax)
938
939 type :: Case_type
940 real(RK), allocatable :: LogNumPointInSpheroid(:)
941 real(RK), allocatable :: LogMeanLogScaleFactor(:)
942 end type Case_type
943#include "pm_clusDensity.inc.scaleFactor.inc.F90"
944
945 ndim = nd
946 if (npmin < nd + 1_IK) then
947 write(output_unit,"(*(g0,:,' '))") "Internal error occurred. npmin < nd + 1_IK :", npmin, nd + 1_IK
948 error stop
949 end if
950 if (ndim < lbound(Case, dim = 1)) then
951 write(output_unit,"(*(g0,:,' '))") "Internal error occurred. ndim < lbound(Case, dim = 1) :", ndim, lbound(Case, dim = 1)
952 error stop
953 end if
954 if (ndim > ubound(Case, dim = 1)) ndim = ubound(Case, dim = 1) ! assume the largest available dimension.
955
956 LogCount(npmin:npmax) = [( log(real(ip,RK)), ip = npmin, npmax )]
957#if INTERP_ENABLED
958 block
959 use InterpLinear_mod, only: InterpLinear_type
960 type(InterpLinear_type) :: Interp
961 Interp = InterpLinear_type(SortedX = Case(ndim)%LogNumPointInSpheroid, SortedY = Case(ndim)%LogMeanLogScaleFactor)
962 BiasCorrectionScaleFactorSq(npmin:npmax) = exp(2_IK*exp(Interp%predict(QueryX = LogCount)))
963 end block
964#else
965 BiasCorrectionScaleFactorSq = 0
966#endif
967#if CHECK_ENABLED
968 ! WARNING: scaleFactorSq must be larger than 1 in all circumstances. This is crucial later on.
969 if (any(BiasCorrectionScaleFactorSq < 1._RK)) error stop
970#else
971 where (BiasCorrectionScaleFactorSq < 1._RK); BiasCorrectionScaleFactorSq = 1._RK; end where
972#endif
973
974 end function getBiasCorrectionScaleFactorSq
975
976!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
977
978 subroutine writePartition ( Partition & ! LCOV_EXCL_LINE
979 , fileUnit & ! LCOV_EXCL_LINE
980 , Point & ! LCOV_EXCL_LINE
981 )
982#if INTEL_COMPILER_ENABLED && DLL_ENABLED && (WINDOWS_ENABLED || DARWIN_ENABLED)
983 !DEC$ ATTRIBUTES DLLEXPORT :: writePartition
984#endif
985 class(partition_type) , intent(in) :: Partition
986 integer(IK) , intent(in) :: fileUnit
987 real(RK) , intent(in) :: Point(Partition%nd, Partition%np)
988 character(*,SK) , parameter :: fileFormat = "(*(g0,:,','))"
989 integer(IK) :: j, ic
990
991 write(fileUnit,"(*(g0,:,'"//new_line("a")//"'))") "nd", Partition%nd, "np", Partition%np, "nc", Partition%ne
992
993 write(fileUnit,"(A)") "numRecursiveCall"
994 write(fileUnit, fileFormat) Partition%numRecursiveCall
995
996 write(fileUnit,"(A)") "pointLogVolNormed"
997 write(fileUnit, fileFormat) Partition%pointLogVolNormed
998
999 write(fileUnit,"(A)") "parentLogVolNormed"
1000 write(fileUnit, fileFormat) Partition%parentLogVolNormed
1001
1002 write(fileUnit,"(A)") "Size"
1003 write(fileUnit, fileFormat) Partition%size(1 : Partition%ne)
1004
1005 write(fileUnit,"(A)") "Center"
1006 write(fileUnit, fileFormat) Partition%Center(1 : Partition%nd, 1 : Partition%ne)
1007
1008 write(fileUnit,"(A)") "LogVolNormed"
1009 write(fileUnit, fileFormat) Partition%LogVolNormed(1 : Partition%ne)
1010
1011 write(fileUnit,"(A)") "ChoLow"
1012 write(fileUnit, fileFormat) ((Partition%ChoLowCovUpp(j : Partition%nd, j, ic), j = 1, Partition%nd), ic = 1, Partition%ne)
1013
1014 write(fileUnit,"(A)") "Point"
1015 write(fileUnit, fileFormat) Point
1016
1017 write(fileUnit,"(A)") "LogLikeFitness"
1018 write(fileUnit, fileFormat) Partition%LogLikeFitness(1 : Partition%ne)
1019
1020 write(fileUnit,"(A)") "Membership"
1021 write(fileUnit, fileFormat) Partition%Membership(1 : Partition%np)
1022
1023 end subroutine writePartition
1024
1025!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1026
1027end module pm_clusDensity
pm_arraySort::setSorted
Sort the input scalar string or contiguous vector in ascending order, or return the sorted indices of...
Definition: pm_arraySort.F90:5540
pm_clusKmeans::setKmeans
Compute and return an iteratively-refined set of cluster centers given the input sample using the k-m...
Definition: pm_clusKmeans.F90:1661
pm_err::getFine
Generate and return a decorated string resulting from the concatenation of getFile(FILE) and getLine(...
Definition: pm_err.F90:382
pm_err::setAsserted
Verify the input assertion holds and if it does not, print the (optional) input message on stdout and...
Definition: pm_err.F90:735
pm_mathCumPropExp::setCumPropExp
Return the cumulative sum of the proportions of the exponential of the input array,...
Definition: pm_mathCumPropExp.F90:569
pm_mathLogSumExp::getLogSumExp
Generate and return the natural logarithm of the sum of the exponential of the input array robustly (...
Definition: pm_mathLogSumExp.F90:147
pm_matrixChol::setMatChol
Compute and return the lower/upper-triangle of the Cholesky factorization of the input Symmetric/Her...
Definition: pm_matrixChol.F90:827
pm_matrixInv::setMatInv
Generate and return the full inverse of a general or triangular matrix or a subset of the inverse of ...
Definition: pm_matrixInv.F90:1786
pm_sampleCov::setCov
Return the covariance matrix corresponding to the input (potentially weighted) correlation matrix or ...
Definition: pm_sampleCov.F90:1694
pm_val2str::getStr
Generate and return the conversion of the input value to an output Fortran string,...
Definition: pm_val2str.F90:167
pm_arraySort
This module contains procedures and generic interfaces for various sorting tasks.
Definition: pm_arraySort.F90:89
pm_clusDensity
This module contains procedures and routines for finding the minimum-volume bounding ellipsoid partit...
Definition: pm_clusDensity.F90:50
pm_clusDensity::POSINF
real(RK), parameter POSINF
Definition: pm_clusDensity.F90:71
pm_clusDensity::partition_typer
function partition_typer(Point, nc, nt, nsim, nemax, minSize, stanEnabled, trimEnabled, kmeansRelTol, logExpansion, logShrinkage if MINVOL_ENABLED
This function can serve as the constructor for partition_type while also performing the partitioning ...
Definition: pm_clusDensity.F90:231
pm_clusDensity::runPartition
subroutine runPartition(Partition, Point, parentLogVolNormed)
This procedure is a method of the class partition_type. Perform recursive clustering of the input Poi...
Definition: pm_clusDensity.F90:361
pm_clusDensity::NEGINF
real(RK), parameter NEGINF
Definition: pm_clusDensity.F90:72
pm_clusDensity::getLogLikeFitness
PURE real(RK) function getLogLikeFitness(count, logVolRatio)
Definition: pm_clusDensity.F90:908
pm_clusDensity::writePartition
subroutine writePartition(Partition, fileUnit, Point)
Definition: pm_clusDensity.F90:982
pm_clusDensity::MODULE_NAME
character(*, SK), parameter MODULE_NAME
Definition: pm_clusDensity.F90:70
pm_clusDensity::getBiasCorrectionScaleFactorSq
real(RK) function, dimension(npmin:npmax) getBiasCorrectionScaleFactorSq(nd, npmin, npmax)
Definition: pm_clusDensity.F90:929
pm_clusKmeans
This module contains procedures and routines for the computing the Kmeans clustering of a given set o...
Definition: pm_clusKmeans.F90:113
pm_err
This module contains classes and procedures for reporting and handling errors.
Definition: pm_err.F90:52
pm_kind
This module defines the relevant Fortran kind type-parameters frequently used in the ParaMonte librar...
Definition: pm_kind.F90:268
pm_kind::RK
integer, parameter RK
The default real kind in the ParaMonte library: real64 in Fortran, c_double in C-Fortran Interoperati...
Definition: pm_kind.F90:543
pm_kind::LK
integer, parameter LK
The default logical kind in the ParaMonte library: kind(.true.) in Fortran, kind(....
Definition: pm_kind.F90:541
pm_kind::IK
integer, parameter IK
The default integer kind in the ParaMonte library: int32 in Fortran, c_int32_t in C-Fortran Interoper...
Definition: pm_kind.F90:540
pm_kind::SK
integer, parameter SK
The default character kind in the ParaMonte library: kind("a") in Fortran, c_char in C-Fortran Intero...
Definition: pm_kind.F90:539
pm_mathCumPropExp
This module contains the procedures and interfaces for computing the cumulative sum of the exponentia...
Definition: pm_mathCumPropExp.F90:47
pm_mathLogSumExp
This module contains the procedures and interfaces for computing the natural logarithm of the sum of ...
Definition: pm_mathLogSumExp.F90:39
pm_matrixChol
This module contains procedures and generic interfaces for computing the Cholesky factorization of po...
Definition: pm_matrixChol.F90:133
pm_matrixInv
This module contains abstract and concrete derived types and procedures related to the inversion of s...
Definition: pm_matrixInv.F90:212
pm_sampleCov
This module contains classes and procedures for computing the properties related to the covariance ma...
Definition: pm_sampleCov.F90:170
pm_val2str
This module contains the generic procedures for converting values of different types and kinds to For...
Definition: pm_val2str.F90:58
pm_clusDensity::KmeansTry_type
Definition: pm_clusDensity.F90:105
pm_clusDensity::KvolTry_type
Definition: pm_clusDensity.F90:74
pm_clusDensity::ParTry_type
Definition: pm_clusDensity.F90:91
pm_clusDensity::partition_type
The partition_type class. Partitions an input array Point(nd,np) with nd attributes and np observatio...
Definition: pm_clusDensity.F90:115
pm_err::err_type
This is the derived type for generating objects to gracefully and verbosely handle runtime unexpected...
Definition: pm_err.F90:157