kmeans_clustering.F90

Go to the documentation of this file.

!! Copyright (C) 2024. A Buccheri.
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
 
#include "global.h"
 
module kmeans_clustering_oct_m
  use, intrinsic :: iso_fortran_env
  use debug_oct_m
  use global_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use mpi_distribute_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use space_oct_m
  use sort_oct_m
  use quickrnd_oct_m
  implicit none
  private
 
  public :: assign_points_to_centroids_finite_bc, &
    update_centroids, &
    weighted_kmeans, &
    sample_initial_centroids
 
contains
 
  ! TODO(Alex) Issue #1004. Implement `assign_points_to_centroids` for periodic boundary conditions
 
  subroutine assign_points_to_centroids_finite_bc(mesh, centroids, ip_to_ic)
    class(mesh_t), intent(in)  :: mesh
    real(real64),  intent(in)  :: centroids(:, :)
    integer,       intent(out) :: ip_to_ic(:)
 
    integer                   :: ip, ic, icen
    integer                   :: n_centroids
    real(real64)              :: min_dist, dist
    real(real64), allocatable :: point(:)
    real(real64), parameter   :: tol = 1.0e-13_real64 
    ! required to distinguish degenerate points, else `if (dist < min_dist)` can vary on
    ! different hardware due to numerical noise. One could equally choose tol to be some
    ! percentage of the grid spacing.
 
    push_sub(assign_points_to_centroids_finite_bc)
 
    ! Grid to centroid index map should have size of the grid
    assert(size(ip_to_ic) == mesh%np)
 
    n_centroids = size(centroids, 2)
    safe_allocate(point(1:size(centroids, 1)))
 
    !$omp parallel do default(shared) private(point, icen, min_dist, dist)
    do ip = 1, mesh%np
      ip_to_ic(ip) = 0
      ! Compute which centroid, grid point `ip` is closest to
      point = mesh%x(:, ip)
      icen = 1
      min_dist = sum((centroids(:, 1) - point(:))**2)
      do ic = 2, n_centroids
        dist = sum((centroids(:, ic) - point(:))**2)
        if (dist < min_dist - tol) then
          min_dist = dist
          icen = ic
        endif
      enddo
      ip_to_ic(ip) = icen
    enddo
    !$omp end parallel do
 
    safe_deallocate_a(point)
 
    pop_sub(assign_points_to_centroids_finite_bc)
 
  end subroutine assign_points_to_centroids_finite_bc
 
 
  subroutine update_centroids(mesh, weight, ip_to_ic, centroids)
    class(mesh_t), intent(in)    :: mesh
    real(real64),  intent(in)    :: weight(:)
    integer,       intent(in)    :: ip_to_ic(:)
    real(real64), contiguous,  intent(inout) :: centroids(:, :)
    !                                                          Out: Updated centroid positions
 
    integer                   :: n_centroids, ic, ip
    real(real64)              :: one_over_denom
    real(real64), allocatable :: denominator(:)
 
    push_sub(update_centroids)
 
    ! The indexing of weight and grid must be consistent => must belong to the same spatial distribution
    ! This does not explicit assert this, but if the grid sizes differ, this is a clear indicator of a problem
    assert(mesh%np == size(weight))
 
    n_centroids = size(centroids, 2)
    safe_allocate(denominator(1:n_centroids))
 
    do ic = 1, n_centroids
      centroids(:, ic) = 0._real64
      denominator(ic) = 0._real64
    enddo
 
    !$omp parallel do private(ic) reduction(+ : centroids, denominator)
    do ip = 1, mesh%np
      ic = ip_to_ic(ip)
      ! Initially accumulate the numerator in `centroids`
      centroids(:, ic) = centroids(:, ic) + (mesh%x(:, ip) * weight(ip))
      denominator(ic) = denominator(ic) + weight(ip)
    enddo
    !$omp end parallel do
 
    ! Gather contributions to numerator and denominator of all centroids, from all domains of the mesh/grid
    call mesh%allreduce(centroids)
    call mesh%allreduce(denominator)
 
    ! If division by zero occurs here it implies that the sum(weight) = 0 for all grid points
    ! in cluster ic. This can occur if the initial centroid is poorly chosen at a point with no
    !! associated weight (such as the vacuum of a crystal cell)
    !$omp parallel do private(one_over_denom) reduction(* : centroids)
    do ic = 1, n_centroids
      one_over_denom = 1._real64 / denominator(ic)
      centroids(:, ic) = centroids(:, ic) * one_over_denom
    enddo
    !$omp end parallel do
 
    safe_deallocate_a(denominator)
    pop_sub(update_centroids)
 
  end subroutine update_centroids
 
 
  subroutine compute_grid_difference(points, updated_points, tol, points_differ)
    real(real64), intent(in)    :: points(:, :)
    real(real64), intent(in)    :: updated_points(:, :)
    real(real64), intent(in)    :: tol
    logical,      intent(out)   :: points_differ(:)
 
    integer                     :: ip, n_dim
    real(real64), allocatable   :: diff(:)
 
    push_sub(compute_grid_difference)
 
    n_dim = size(points, 1)
    allocate(diff(n_dim))
 
    !$omp parallel do default(shared) private(diff)
    do ip = 1, size(points, 2)
      diff(:) = abs(updated_points(:, ip) - points(:, ip))
      points_differ(ip) = any(diff > tol)
    enddo
    !$omp end parallel do
 
    if(debug%info) then
      call report_differences_in_grids(points, updated_points, tol, points_differ)
    endif
 
    pop_sub(compute_grid_difference)
 
  end subroutine compute_grid_difference
 
 
  subroutine report_differences_in_grids(points, updated_points, tol, points_differ)
    real(real64), intent(in)  :: points(:, :)
    real(real64), intent(in)  :: updated_points(:, :)
    real(real64), intent(in)  :: tol
    logical,      intent(in)  :: points_differ(:)
 
    integer,      allocatable :: indices(:)
    integer                   :: i, j, n_unconverged, ndim
    character(len=50)         :: f_string
    real(real64), allocatable :: diff(:)
 
    push_sub(report_differences_in_grids)
 
    indices = pack([(i, i=1,size(points_differ))], points_differ)
    n_unconverged = size(indices)
    ndim = size(points, 1)
    allocate(diff(ndim))
 
    write(f_string, '(A, I1, A, I1, A, I1, A)') '(', &
    & ndim, '(F16.10, X), ', &
    & ndim, '(F16.10, X), ', &
    & ndim, '(F16.10, X), F16.10)'
 
    write(message(1), '(a)') "# Current Point  ,  Prior Point  ,  |ri - r_{i-1}|  ,  tol"
    call messages_info(1)
    do j = 1, n_unconverged
      i = indices(j)
      diff(:) = abs(updated_points(:, i) - points(:, i))
      write(message(1), f_string) updated_points(:, i), points(:, i), diff, tol
      call messages_info(1)
    enddo
    write(message(1), *) "Summary:", n_unconverged, "of out", size(points, 2), "are not converged"
    call messages_info(1)
 
    pop_sub(report_differences_in_grids)
 
  end subroutine report_differences_in_grids
 
 
  subroutine weighted_kmeans(space, mesh, weight, centroids, n_iter, centroid_tol, discretize, inertia)
    class(space_t),  intent(in)    :: space
    class(mesh_t),   intent(in)    :: mesh
    real(real64),    intent(in)    :: weight(:)
    real(real64), contiguous,   intent(inout) :: centroids(:, :)
    !                                                           Out: Final centroids
    integer,      optional, intent(in ) :: n_iter
    real(real64), optional, intent(in ) :: centroid_tol
    logical,      optional, intent(in ) :: discretize
    real(real64), optional, intent(out) :: inertia
 
    logical                   :: discretize_centroids
    integer                   :: n_iterations, n_centroid, i
    real(real64)              :: tol
    integer,      allocatable :: ip_to_ic(:)
    real(real64), allocatable :: prior_centroids(:, :)
    logical,      allocatable :: points_differ(:)
 
    push_sub(weighted_kmeans)
 
    n_iterations = optional_default(n_iter, 200)
    tol = optional_default(centroid_tol, 1.e-4_real64)
    discretize_centroids = optional_default(discretize, .true.)
 
    ! Should use a positive number of iterations
    assert(n_iterations >= 1)
    ! Number of weights inconsistent with number of grid points
    assert(size(weight) == mesh%np)
    ! Spatial dimensions of centroids array is inconsistent
    assert(size(centroids, 1) == space%dim)
    ! Assignment of points to centroids only implemented for finite BCs
    assert(.not. space%is_periodic())
 
    ! Work arrays
    n_centroid = size(centroids, 2)
    safe_allocate_source(prior_centroids(space%dim, size(centroids, 2)), centroids)
    safe_allocate(ip_to_ic(1:mesh%np))
    safe_allocate(points_differ(1:n_centroid))
 
    write(message(1), '(a)') 'Debug: Performing weighted Kmeans clustering '
    call messages_info(1, debug_only=.true.)
 
    do i = 1, n_iterations
      write(message(1), '(a, I3)') 'Debug: Iteration ', i
      call messages_info(1, debug_only=.true.)
      ! TODO(Alex) Issue #1004. Implement `assign_points_to_centroids` for periodic boundary conditions
      call assign_points_to_centroids_finite_bc(mesh, centroids, ip_to_ic)
 
      call update_centroids(mesh, weight, ip_to_ic, centroids)
      call compute_grid_difference(prior_centroids, centroids, tol, points_differ)
 
      if (any(points_differ)) then
        prior_centroids = centroids
      else
        write(message(1), '(a)') 'Debug: All centroid points converged'
        call messages_info(1, debug_only=.true.)
        ! Break loop
        exit
      endif
 
    enddo
 
    if (discretize_centroids) then
      call mesh_discretize_values_to_mesh(mesh, centroids)
    endif
 
    if (present(inertia)) then
      call compute_centroid_inertia(mesh, centroids, weight, ip_to_ic, inertia)
    endif
 
    safe_deallocate_a(prior_centroids)
    safe_deallocate_a(ip_to_ic)
    safe_deallocate_a(points_differ)
 
    pop_sub(weighted_kmeans)
 
  end subroutine weighted_kmeans
 
 
  subroutine sample_initial_centroids(mesh, centroids, seed_value)
    class(mesh_t),             intent(in )           :: mesh
    real(real64), contiguous,  intent(out)           :: centroids(:, :)
    integer(int64),            intent(inout), optional :: seed_value
    ! This will get mutated by the Fisher Yates shuffle
 
    integer(int32)              :: n_centroids
    integer(int64), allocatable :: centroid_idx(:)
    integer(int64)              :: ipg
    integer(int32)              :: ic
 
    push_sub(sample_initial_centroids)
 
    n_centroids = size(centroids, 2)
    safe_allocate(centroid_idx(1:n_centroids))
 
    !  Choose n_centroids indices from [1, np_global]
    call fisher_yates_shuffle(n_centroids, mesh%np_global, seed_value, centroid_idx)
 
    ! Convert ip_global to (x,y,z)
    do ic = 1, n_centroids
      ipg = centroid_idx(ic)
      centroids(:, ic) = mesh_x_global(mesh, ipg)
    enddo
 
    safe_deallocate_a(centroid_idx)
 
    pop_sub(sample_initial_centroids)
 
  end subroutine sample_initial_centroids
 
 
  subroutine compute_centroid_inertia(mesh, centroids, weight, ip_to_ic, inertia)
    class(mesh_t),   intent(in)    :: mesh
    real(real64),    intent(in)    :: centroids(:, :)
    real(real64),    intent(in)    :: weight(:)
    integer,         intent(in)    :: ip_to_ic(:)
    real(real64),    intent(out)   :: inertia
    integer :: ip, ic
 
    inertia = 0.0_real64
 
    !$omp parallel do private(ip) reduction(+ : inertia)
    do ip = 1, mesh%np
      ic = ip_to_ic(ip)
      inertia = inertia + weight(ip) * sum((centroids(:, ic) - mesh%x(:, ip))**2)
    enddo
    !$omp end parallel do
 
    call mesh%allreduce(inertia)
 
  end subroutine compute_centroid_inertia
 
end module kmeans_clustering_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: