https://www.cdslab.org/paramonte/fortran/2
Current view: top level - main - pm_clustering@routines.inc.F90 (source / functions) Hit Total Coverage
Test: ParaMonte 2.0.0 :: Serial Fortran - Code Coverage Report Lines: 169 176 96.0 %
Date: 2024-04-08 03:18:57 Functions: 0 0 -
Legend: Lines: hit not hit

          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_clustering](@ref pm_clustering).
      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             : #define INTRINSIC_ENABLED 0
      28             : 
      29             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      30             : #if     setMember_ENABLED && D0_D1_ENABLED
      31             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      32             : 
      33             :         real(RKC) :: temp
      34             :         integer(IK) :: icls
      35           1 :         CHECK_ASSERTION(__LINE__, 0_IK < size(center, 1, IK), SK_"@setMember(): The condition `0 < size(center, rank(center))` must hold. shape(center) = "//getStr(shape(center, IK)))
      36             :         ! First cluster.
      37             :         disq = 0._RKC
      38           1 :         disq = (sample - center(1))**2
      39           1 :         membership = 1
      40             :         ! All the rest of clusters.
      41           3 :         do icls = 2, size(center, 1, IK)
      42           2 :             temp = (sample - center(icls))**2
      43           3 :             if (temp < disq) then
      44           0 :                 disq = temp
      45           0 :                 membership = icls
      46             :             end if
      47             :         end do
      48             : 
      49             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      50             : #elif   setMember_ENABLED && D1_D1_ENABLED
      51             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      52             : 
      53             :         real(RKC) :: temp
      54             :         integer(IK) :: icls, isam, nsam
      55           1 :         nsam = size(sample, 1, IK)
      56           1 :         CHECK_ASSERTION(__LINE__, 0_IK <  size(center, 1, IK), SK_"@setMember(): The condition `0 < size(center, 2)` must hold. size(center, 1) = "//getStr(size(center, 1, IK)))
      57           4 :         CHECK_ASSERTION(__LINE__, nsam == size(disq, 1, IK), SK_"@setMember(): The condition `size(sample, rank(sample)) == size(disq)` must hold. shape(sample), size(disq) = "//getStr([shape(sample, IK), size(disq, 1, IK)]))
      58           4 :         CHECK_ASSERTION(__LINE__, nsam == size(membership, 1, IK), SK_"@setMember(): The condition `size(sample, rank(sample)) == size(membership)` must hold. shape(sample), size(membership) = "//getStr([shape(sample, IK), size(membership, 1, IK)]))
      59             :         ! First cluster.
      60           5 :         do isam = 1, nsam
      61           4 :             disq(isam) = (sample(isam) - center(1))**2
      62           5 :             membership(isam) = 1
      63             :         end do
      64             :         ! All the rest of clusters.
      65           3 :         do icls = 2, size(center, 1, IK)
      66          11 :             do isam = 1, nsam
      67           8 :                 temp = (sample(isam) - center(icls))**2
      68          10 :                 if (temp < disq(isam)) then
      69           2 :                     disq(isam) = temp
      70           2 :                     membership(isam) = icls
      71             :                 end if
      72             :             end do
      73             :         end do
      74             : 
      75             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      76             : #elif   setMember_ENABLED && D1_D2_ENABLED
      77             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      78             : 
      79             :         real(RKC) :: temp
      80             :         integer(IK) :: icls, idim, ndim
      81           1 :         ndim = size(sample, 1, IK)
      82           1 :         CHECK_ASSERTION(__LINE__, 0_IK <  size(center, rank(center), IK), SK_"@setMember(): The condition `0 < size(center, rank(center))` must hold. shape(center) = "//getStr(shape(center, IK)))
      83           3 :         CHECK_ASSERTION(__LINE__, ndim == size(center, 1, IK), SK_"@setMember(): The condition `size(sample, 1) == size(center, 1)` must hold. size(sample, 1), size(center, 1) = "//getStr([size(sample, 1, IK), size(center, 1, IK)]))
      84             :         ! First cluster.
      85             :         !disq = sum((sample - center(1 : ndim, 1))**2)
      86           1 :         disq = 0._RKC
      87           4 :         do idim = 1, ndim
      88           4 :             disq = disq + (sample(idim) - center(idim, 1))**2
      89             :         end do
      90           1 :         membership = 1
      91             :         ! All the rest of clusters.
      92           3 :         do icls = 2, size(center, 2, IK)
      93             :            !temp = sum((sample  - center(1 : ndim, icls))**2)
      94           0 :             temp = 0._RKC
      95           8 :             do idim = 1, ndim
      96           8 :                 temp = temp + (sample(idim) - center(idim, icls))**2
      97             :             end do
      98           3 :             if (temp < disq) then
      99           1 :                 disq = temp
     100           1 :                 membership = icls
     101             :             end if
     102             :         end do
     103             : 
     104             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     105             : #elif   setMember_ENABLED && D2_D2_ENABLED
     106             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     107             : 
     108             :         real(RKC) :: temp
     109             :         integer(IK) :: icls, idim, ndim, isam, nsam
     110          38 :         ndim = size(sample, 1, IK)
     111          38 :         nsam = size(sample, 2, IK)
     112          38 :         CHECK_ASSERTION(__LINE__, 0_IK <  size(center, 2, IK), SK_"@setMember(): The condition `0 < size(center, rank(center))` must hold. shape(center) = "//getStr(shape(center, IK)))
     113         228 :         CHECK_ASSERTION(__LINE__, nsam == size(disq, 1, IK), SK_"@setMember(): The condition `size(sample, rank(sample)) == size(disq)` must hold. shape(sample), size(disq) = "//getStr([shape(sample, IK), size(disq, 1, IK)]))
     114         114 :         CHECK_ASSERTION(__LINE__, ndim == size(center, 1, IK), SK_"@setMember(): The condition `size(sample, 1) == size(center, 1)` must hold. size(sample, 1), size(center, 1) = "//getStr([size(sample, 1, IK), size(center, 1, IK)]))
     115         228 :         CHECK_ASSERTION(__LINE__, nsam == size(membership, 1, IK), SK_"@setMember(): The condition `size(sample, rank(sample)) == size(membership, 1)` must hold. shape(sample), size(membership, 1) = "//getStr([shape(sample, IK), size(membership, 1, IK)]))
     116             :         ! First cluster.
     117             :         icls = 1
     118      130086 :         do isam = 1, nsam
     119             :             !disq(isam) = sum((center(1 : ndim, icls) - sample(1 : ndim, isam))**2)
     120      130048 :             disq(isam) = 0._RKC
     121      390180 :             do idim = 1, ndim
     122      390180 :                 disq(isam) = disq(isam) + (sample(idim, isam) - center(idim, icls))**2
     123             :             end do
     124      130086 :             membership(isam) = icls
     125             :         end do
     126             :         ! All the rest of clusters.
     127         165 :         do icls = 2, size(center, 2, IK)
     128      520283 :             do isam = 1, nsam
     129             :                 !temp = sum((center(1 : ndim, icls) - sample(1 : ndim, isam))**2)
     130           0 :                 temp = 0._RKC
     131     1560452 :                 do idim = 1, ndim
     132     1560452 :                     temp = temp + (sample(idim, isam) - center(idim, icls))**2
     133             :                 end do
     134      520245 :                 if (temp < disq(isam)) then
     135      184792 :                     disq(isam) = temp
     136      184792 :                     membership(isam) = icls
     137             :                 end if
     138             :             end do
     139             :         end do
     140             : 
     141             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     142             : #elif   setCenter_ENABLED && D1_D1_ENABLED
     143             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     144             : 
     145             :         integer(IK) :: icls, ncls, isam, nsam
     146             : 
     147           1 :         nsam = ubound(sample, rank(sample), IK)
     148           1 :         ncls = ubound(center, rank(center), IK)
     149             : 
     150           5 :         CHECK_ASSERTION(__LINE__, all(0._RKC <= disq), SK_"@setMember(): The condition `all(0 <= disq)` must hold. disq = "//getStr(disq))
     151           5 :         CHECK_ASSERTION(__LINE__, nsam == ubound(disq, 1, IK), SK_"@setMember(): The condition `size(sample, rank(sample)) == size(disq)` must hold. shape(sample), size(disq) = "//getStr([shape(sample, IK), ubound(disq, 1, IK)]))
     152             :        !CHECK_ASSERTION(__LINE__, nsam >= ubound(center, rank(center), IK), SK_"@setMember(): The condition `size(sample, rank(sample)) >= size(center, rank(center))` must hold. shape(sample), shape(center) = "//getStr([shape(sample, IK), shape(center, IK)]))
     153           5 :         CHECK_ASSERTION(__LINE__, nsam == ubound(membership, 1, IK), SK_"@setMember(): The condition `size(sample, rank(sample)) == size(membership)` must hold. shape(sample), size(membership) = "//getStr([shape(sample, IK), ubound(membership, 1, IK)]))
     154           5 :         CHECK_ASSERTION(__LINE__, ncls == ubound(potential, 1, IK), SK_"@setMember(): The condition `size(center, rank(center)) == size(potential)` must hold. shape(center), size(potential) = "//getStr([shape(center, IK), ubound(potential, 1, IK)]))
     155           5 :         CHECK_ASSERTION(__LINE__, ncls == ubound(size, 1, IK), SK_"@setMember(): The condition `size(center, rank(center)) == size(size)` must hold. shape(center), size(size) = "//getStr([shape(center, IK), ubound(size, 1, IK)]))
     156           1 :         CHECK_ASSERTION(__LINE__, 0_IK <  ncls, SK_"@setMember(): The condition `0 < size(center, rank(center))` must hold. shape(center) = "//getStr(shape(center, IK)))
     157          15 :         CHECK_ASSERTION(__LINE__, all(0_IK < membership .and. membership <= ncls), SK_"@setMember(): The condition `all(0 < membership .and. membership <= size(center, rank(center)))` must hold. size(center, rank(center)), membership = "//getStr([ncls, membership]))
     158             : 
     159             :         ! Initialize.
     160             : 
     161             :         do concurrent(icls = 1 : ncls)
     162           2 :             potential(icls) = 0._RKC
     163           2 :             center(icls) = 0._RKC
     164           3 :             size(icls) = 0_IK
     165             :         end do
     166             : 
     167             :         ! compute new centers.
     168             : 
     169           5 :         do isam = 1, nsam
     170           4 :             center(membership(isam)) = center(membership(isam)) + sample(isam)
     171           4 :             potential(membership(isam)) = potential(membership(isam)) + disq(isam)
     172           5 :             size(membership(isam)) = size(membership(isam)) + 1
     173             :         end do
     174             : 
     175             :         ! normalize new centers.
     176             : 
     177             :         do concurrent(icls = 1 : ncls)
     178           3 :             if (0_IK < size(icls)) center(icls) = center(icls) / real(size(icls), RKC)
     179             :         end do
     180             : 
     181             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     182             : #elif   setCenter_ENABLED && D2_D2_ENABLED
     183             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     184             : 
     185             :         real(RKC) :: sizeinv
     186             :         integer(IK) :: icls, ncls, ndim, isam, nsam
     187             : 
     188          36 :         ndim = ubound(sample, 1, IK)
     189          36 :         nsam = ubound(sample, 2, IK)
     190          36 :         ncls = ubound(center, 2, IK)
     191             : 
     192      125080 :         CHECK_ASSERTION(__LINE__, all(0._RKC <= disq), SK_"@setMember(): The condition `all(0 <= disq)` must hold. disq = "//getStr(disq))
     193         252 :         CHECK_ASSERTION(__LINE__, nsam == ubound(disq, 1, IK), SK_"@setMember(): The condition `ubound(sample, rank(sample)) == ubound(disq, 1)` must hold. shape(sample), ubound(disq, 1) = "//getStr([shape(sample, IK), ubound(disq, 1, IK)]))
     194         144 :         CHECK_ASSERTION(__LINE__, ndim == ubound(center, 1, IK), SK_"@setMember(): The condition `ubound(sample, 1) == ubound(center, 1)` must hold. ubound(sample, 1), ubound(center, 1) = "//getStr([ubound(sample, 1, IK), ubound(center, 1, IK)]))
     195             :        !CHECK_ASSERTION(__LINE__, nsam >= ubound(center, rank(center), IK), SK_"@setMember(): The condition `ubound(sample, rank(sample)) >= ubound(center, rank(center))` must hold. shape(sample), shape(center) = "//getStr([shape(sample, IK), shape(center, IK)]))
     196         252 :         CHECK_ASSERTION(__LINE__, nsam == ubound(membership, 1, IK), SK_"@setMember(): The condition `ubound(sample, rank(sample)) == ubound(membership, 1)` must hold. shape(sample), ubound(membership, 1) = "//getStr([shape(sample, IK), ubound(membership, 1, IK)]))
     197         252 :         CHECK_ASSERTION(__LINE__, ncls == ubound(potential, 1, IK), SK_"@setMember(): The condition `ubound(center, rank(center)) == ubound(potential, 1)` must hold. shape(center), ubound(potential, 1) = "//getStr([shape(center, IK), ubound(potential, 1, IK)]))
     198         252 :         CHECK_ASSERTION(__LINE__, ncls == ubound(size, 1, IK), SK_"@setMember(): The condition `ubound(center, rank(center)) == ubound(size, 1)` must hold. shape(center), ubound(size, 1) = "//getStr([shape(center, IK), ubound(size, 1, IK)]))
     199          36 :         CHECK_ASSERTION(__LINE__, 0_IK <  ncls, SK_"@setMember(): The condition `0 < size(center, rank(center))` must hold. shape(center) = "//getStr(shape(center, IK)))
     200      375240 :         CHECK_ASSERTION(__LINE__, all(0 < membership .and. membership <= ncls), SK_"@setMember(): The condition `all(0 < membership .and. membership <= size(center, rank(center)))` must hold. size(center, rank(center)), membership = "//getStr([ubound(center, 2, IK), membership]))
     201             : 
     202             :         ! Initialize.
     203             : 
     204             :         do concurrent(icls = 1 : ncls)
     205         489 :             center(1 : ndim, icls) = 0._RKC
     206         157 :             potential(icls) = 0._RKC
     207         193 :             size(icls) = 0_IK
     208             :         end do
     209             : 
     210             :         ! compute new centers.
     211             : 
     212      125080 :         do isam = 1, nsam
     213      375164 :             center(1 : ndim, membership(isam)) = center(1 : ndim, membership(isam)) + sample(1 : ndim, isam)
     214      125044 :             potential(membership(isam)) = potential(membership(isam)) + disq(isam)
     215      125080 :             size(membership(isam)) = size(membership(isam)) + 1
     216             :         end do
     217             : 
     218             :         ! normalize new centers.
     219             : 
     220             :         do concurrent(icls = 1 : ncls)
     221         193 :             if (0_IK < size(icls)) then
     222         157 :                 sizeinv = 1._RKC / real(size(icls), RKC)
     223         489 :                 center(1 : ndim, icls) = center(1 : ndim, icls) * sizeinv
     224             :             end if
     225             :         end do
     226             : 
     227             :         !%%%%%%%%%%%%%%%%%%
     228             : #elif   setKmeansPP_ENABLED
     229             :         !%%%%%%%%%%%%%%%%%%
     230             : 
     231             :         real(RKC)   :: temp
     232          39 :         real(RKC)   :: cumsumdisq(0 : ubound(sample, 2, IK)) ! Sum of distance squared of points up to the given index of the vector.
     233             :         integer(IK) :: idim, icls, isam, ndim, nsam, cid
     234             : #if     Optional_ENABLED
     235             :         integer(IK) :: ncls
     236          12 :         ncls = ubound(center, 2, IK)
     237          60 :         CHECK_ASSERTION(__LINE__, ubound(sample, 1, IK) == ubound(center, 1, IK), SK_"@setKmeansPP(): The condition `ubound(sample, 1) == ubound(center, 1)` must hold. ubound(sample, 1), ubound(center, 1) = "//getStr([ubound(sample, 1, IK), ubound(center, 1, IK)]))
     238          48 :         CHECK_ASSERTION(__LINE__, ncls == ubound(potential, 1, IK), SK_"@setKmeansPP(): The condition `ncls == ubound(potential, 1)` must hold. ncls, ubound(potential, 1) = "//getStr([ncls, ubound(potential, 1, IK)]))
     239          48 :         CHECK_ASSERTION(__LINE__, ncls == ubound(size, 1, IK), SK_"@setKmeansPP(): The condition `ncls == ubound(size, 1)` must hold. ncls, ubound(size, 1) = "//getStr([ncls, ubound(size, 1, IK)]))
     240             : #elif   !Default_ENABLED
     241             : #error  "Unrecognized interface."
     242             : #endif
     243          10 :         ndim = ubound(sample, 1, IK)
     244          22 :         nsam = ubound(sample, 2, IK)
     245          22 :         CHECK_ASSERTION(__LINE__, 0_IK  < ncls, SK_"@setKmeansPP(): The condition `0 < ncls` must hold. ncls = "//getStr(ncls))
     246          22 :         CHECK_ASSERTION(__LINE__, 0_IK  < ubound(sample, 2, IK), SK_"@setKmeansPP(): The condition `0 < ubound(sample, rank(sample))` must hold. shape(sample) = "//getStr(shape(sample, IK)))
     247         154 :         CHECK_ASSERTION(__LINE__, nsam == ubound(disq, 1, IK), SK_"@setMember(): The condition `ubound(sample, rank(sample)) == ubound(disq, 1)` must hold. shape(sample), size(disq) = "//getStr([shape(sample, IK), ubound(disq, 1, IK)]))
     248          88 :         CHECK_ASSERTION(__LINE__, nsam == ubound(membership, 1, IK), SK_"@setKmeansPP(): The condition `ubound(sample, rank(sample)) == ubound(membership, 1)` must hold. ubound(sample, rank(sample)), ubound(membership, 1) = "//getStr([nsam, ubound(membership, 1, IK)]))
     249          66 :         CHECK_ASSERTION(__LINE__, ncls <= nsam, SK_"@setKmeansPP(): The condition `ncls < ubound(sample, rank(sample))` must hold. ncls, ubound(sample, rank(sample)) = "//getStr([ncls, nsam]))
     250             : 
     251          22 :         cumsumdisq(0) = 0._RKC ! This element must always remain zero.
     252             : 
     253             :         ! Define the first cluster center.
     254             : 
     255             :         icls = 1
     256          22 :         call setUnifRand(rng, cid, 1_IK, nsam) ! This random assignment is unnecessary if the `sample` is in random order.
     257             : #if     Optional_ENABLED
     258          52 :         center(1 : ndim, icls) = sample(1 : ndim, cid)
     259             : #endif
     260       10095 :         do isam = 1, nsam
     261             :             ! Find the distance-squared of each sample to the #icls cluster cid.
     262             : #if         INTRINSIC_ENABLED
     263             :             disq(isam) = sum((sample(1 : ndim, isam) - sample(1 : ndim, cid))**2)
     264             : #else
     265       30298 :             disq(isam) = 0._RKC; do idim = 1, ndim; disq(isam) = disq(isam) + (sample(idim, isam) - sample(idim, cid))**2; end do
     266             : #endif
     267       10073 :             cumsumdisq(isam) = cumsumdisq(isam - 1) + disq(isam)
     268       10095 :             membership(isam) = icls
     269             :         end do
     270             : 
     271             :         ! Find the second cluster center.
     272             : 
     273          22 :         call setUnifRand(rng, temp)
     274          22 :         temp = temp * cumsumdisq(nsam)
     275             : #if     INTRINSIC_ENABLED
     276             :         cid = minloc(cumsumdisq, 1, temp < cumsumdisq) - 1 ! `-1` compensates for the 0 starting index of `cumsumdisq`.
     277             :         !if (cid < 1) cid = nsam + 1
     278             : #else
     279        7062 :         do cid = 1, nsam; if (temp < cumsumdisq(cid)) exit; end do
     280             : #endif
     281             :         !print *, "ndim, nsam, ncls, cid, temp", ndim, nsam, ncls, cid, temp, cumsumdisq(1:nsam)
     282             : 
     283             :         ! Take care of the unusual case of only one cluster with a single sample.
     284             : 
     285          22 :         if (1_IK == ncls) then
     286             : #if         Optional_ENABLED
     287          13 :             potential(1) = sum(disq)
     288           5 :             size(1) = nsam
     289             : #endif
     290             :             return
     291             :         end if
     292             : 
     293             :         ! Find the rest of cluster centers.
     294             : 
     295          40 :         do icls = 2, ncls - 1
     296             : #if         Optional_ENABLED
     297          53 :             center(1 : ndim, icls) = sample(1 : ndim, cid)
     298             : #endif
     299             :             !if (nsam < cid) exit ! only if nsam < ncls.
     300       30144 :             do isam = 1, nsam
     301             : #if             INTRINSIC_ENABLED
     302             :                 temp = sum((sample(1 : ndim, isam) - sample(1 : ndim, cid))**2)
     303             : #else
     304       90519 :                 temp = 0._RKC; do idim = 1, ndim; temp = temp + (sample(idim, isam) - sample(idim, cid))**2; end do
     305             : #endif
     306       30119 :                 if (temp < disq(isam)) then
     307       10608 :                     disq(isam) = temp
     308       10608 :                     membership(isam) = icls
     309             :                 end if
     310       30144 :                 cumsumdisq(isam) = cumsumdisq(isam - 1) + disq(isam)
     311             :             end do
     312          25 :             call setUnifRand(rng, temp)
     313          25 :             temp = temp * cumsumdisq(nsam)
     314             : #if         INTRINSIC_ENABLED
     315             :             cid = minloc(cumsumdisq, 1, temp < cumsumdisq) - 1
     316             :             !if (cid < 1) cid = nsam + 1
     317             : #else
     318       14958 :             do cid = 1, nsam; if (temp < cumsumdisq(cid)) exit; end do
     319             : #endif
     320             :         end do
     321             : 
     322             :         ! Initialize the output.
     323             : 
     324             : #if     Optional_ENABLED
     325             :         do concurrent(icls = 1 : ncls)
     326             :             !center(1 : ndim, icls) = 0._RKC
     327          27 :             potential(icls) = 0._RKC
     328          34 :             size(icls) = 0_IK
     329             :         end do
     330          29 :         center(1 : ndim, ncls) = sample(1 : ndim, cid)
     331             : #endif
     332             : 
     333             :         ! Final re-computation of the cluster centers and sizes and sample memberships.
     334             : 
     335       10077 :         do isam = 1, nsam
     336             : #if         INTRINSIC_ENABLED
     337             :             temp = sum((sample(1 : ndim, isam) - sample(1 : ndim, cid))**2)
     338             : #else
     339       30253 :             temp = 0._RKC; do idim = 1, ndim; temp = temp + (sample(idim, isam) - sample(idim, cid))**2; end do
     340             : #endif
     341       10070 :             if (temp < disq(isam)) then
     342        1475 :                 disq(isam) = temp
     343        1475 :                 membership(isam) = ncls
     344             :             end if
     345             : #if         Optional_ENABLED
     346             :            !center(1 : ndim, membership(isam)) = center(1 : ndim, membership(isam)) + sample(1 : ndim, isam)
     347       10025 :             potential(membership(isam)) = potential(membership(isam)) + disq(isam)
     348       10032 :             size(membership(isam)) = size(membership(isam)) + 1
     349             : #endif
     350             :         end do
     351             : 
     352             :         ! Normalize sample. This requires the sample size to be larger than the number of clusters to have `all(size > 0)` hold.
     353             : 
     354             :         !do concurrent(icls = 1 : ncls)
     355             :         !    temp = 1._RKC / real(size(icls), RKC)
     356             :         !    center(1 : ndim, icls) = center(1 : ndim, icls) * temp
     357             :         !end do
     358             : 
     359             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     360             : #elif   setKmeans_ENABLED && Init_ENABLED
     361             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     362             : 
     363          11 :         integer(IK) :: retin, sizemin, retinxam, membersnew(ubound(membership, 1, IK))
     364             : 
     365        5051 :         CHECK_ASSERTION(__LINE__, all(0._RKC <= disq), SK_"@setKmeans(): The condition `all(0 <= disq)` must hold. disq = "//getStr(disq))
     366          88 :         CHECK_ASSERTION(__LINE__, ubound(sample, 2, IK) == ubound(disq, 1, IK), SK_"@setKmeans(): The condition `size(sample, rank(sample)) == size(disq)` must hold. shape(sample), size(disq) = "//getStr([shape(sample, IK), ubound(disq, 1, IK)]))
     367          55 :         CHECK_ASSERTION(__LINE__, ubound(sample, 1, IK) == ubound(center, 1, IK), SK_"@setKmeans(): The condition `size(sample, 1) == size(center, 1)` must hold. size(sample, 1), size(center, 1) = "//getStr([ubound(sample, 1, IK), ubound(center, 1, IK)]))
     368             :        !CHECK_ASSERTION(__LINE__, ubound(sample, 2, IK) >= ubound(center, rank(center), IK), SK_"@setKmeans(): The condition `size(sample, rank(sample)) >= size(center, rank(center))` must hold. shape(sample), shape(center) = "//getStr([shape(sample, IK), shape(center, IK)]))
     369          77 :         CHECK_ASSERTION(__LINE__, ubound(sample, 2, IK) == ubound(membership, 1, IK), SK_"@setKmeans(): The condition `size(sample, rank(sample)) == size(membership)` must hold. shape(sample), size(membership) = "//getStr([shape(sample, IK), ubound(membership, 1, IK)]))
     370          88 :         CHECK_ASSERTION(__LINE__, ubound(center, 2, IK) == ubound(potential, 1, IK), SK_"@setKmeans(): The condition `size(center, rank(center)) == size(potential)` must hold. shape(center), size(potential) = "//getStr([shape(center, IK), ubound(potential, 1, IK)]))
     371          77 :         CHECK_ASSERTION(__LINE__, ubound(center, 2, IK) == ubound(size, 1, IK), SK_"@setKmeans(): The condition `size(center, rank(center)) == size(size)` must hold. shape(center), size(size) = "//getStr([shape(center, IK), ubound(size, 1, IK)]))
     372          11 :         CHECK_ASSERTION(__LINE__, ubound(center, 2, IK) >  0_IK, SK_"@setKmeans(): The condition `0 < size(center, rank(center))` must hold. shape(center) = "//getStr(shape(center, IK)))
     373       15153 :         CHECK_ASSERTION(__LINE__, all(0 < membership .and. membership <= ubound(center, 2, IK)), SK_"@setKmeans(): The condition `all(0 < membership .and. membership <= size(center, rank(center)))` must hold. size(center, rank(center)), membership = "//getStr([ubound(center, 2, IK), membership]))
     374             : 
     375          11 :         if (present(maxniter)) then
     376           0 :             CHECK_ASSERTION(__LINE__, 0_IK <= maxniter, SK_"@setKmeans(): The condition `0 <= maxniter` must hold. maxniter = "//getStr(maxniter))
     377             :             retinxam = maxniter
     378             :         else
     379             :             retinxam = 300_IK
     380             :         end if
     381          11 :         if (present(minsize)) then
     382           0 :             CHECK_ASSERTION(__LINE__, 0_IK <= minsize, SK_"@setKmeans(): The condition `0 <= minsize` must hold. minsize = "//getStr(minsize))
     383             :             sizemin = minsize
     384             :         else
     385             :             sizemin = 1_IK
     386             :         end if
     387             : 
     388             :         retin = 1
     389          11 :         do
     390             : 
     391             :             ! compute the new centers, based on the updated memberships.
     392             : 
     393          22 :             call setCenter(membership, disq, sample, center, size, potential)
     394         112 :             failed = any(size < sizemin)
     395          22 :             if (failed) exit
     396             : 
     397             :             ! compute the new memberships, based on the input initial centers.
     398             : 
     399          22 :             call setMember(membersnew, disq, sample, center)
     400        5909 :             if (all(membersnew == membership)) exit
     401          12 :             failed = retinxam < retin
     402          12 :             if (failed) exit
     403             : 
     404          12 :             retin = retin + 1
     405             : 
     406             :             ! compute the new centers, based on the updated memberships.
     407             : 
     408          12 :             call setCenter(membersnew, disq, sample, center, size, potential)
     409          72 :             failed = any(size < sizemin)
     410          12 :             if (failed) exit
     411             : 
     412             :             ! compute the new memberships, based on the input initial centers.
     413             : 
     414          12 :             call setMember(membership, disq, sample, center)
     415       10282 :             if (all(membersnew == membership)) exit
     416          11 :             failed = retinxam < retin
     417          11 :             if (failed) exit
     418             : 
     419          22 :             retin = retin + 1
     420             : 
     421             :         end do
     422          11 :         if (present(niter)) niter = retin
     423             : 
     424             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     425             : #elif   setKmeans_ENABLED && (RNGF_ENABLED || RNGX_ENABLED)
     426             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     427             : 
     428             :         integer(IK) :: nfail_def, itry, ncls
     429             : 
     430          10 :         if (present(nfail)) then
     431           0 :             CHECK_ASSERTION(__LINE__, 0_IK <= nfail, SK_"@setKmeans(): The condition `0 <= nfail` must hold. nfail = "//getStr(nfail))
     432             :             nfail_def = nfail
     433             :         else
     434             :             nfail_def = 1
     435             :         end if
     436             : 
     437          20 :         ncls = ubound(center, 2, IK)
     438          10 :         do itry = 1, nfail_def
     439          10 :             call setKmeansPP(rng, membership, disq, sample, ncls)
     440          10 :             call setKmeans(membership, disq, sample, center, size, potential, failed, niter, maxniter, minsize)
     441          10 :             if (.not. failed) exit
     442             :         end do
     443             : 
     444             : #else
     445             :         !%%%%%%%%%%%%%%%%%%%%%%%%
     446             : #error  "Unrecognized interface."
     447             :         !%%%%%%%%%%%%%%%%%%%%%%%%
     448             : #endif
     449             : 
     450             : #undef  INTRINSIC_ENABLED

ParaMonte: Parallel Monte Carlo and Machine Learning Library 
The Computational Data Science Lab
© Copyright 2012 - 2024