PDAF_pf_resampling.F90 Source File


Source Code

! Copyright (c) 2019-2024 Lars Nerger
!
! This file is part of PDAF.
!
! PDAF is free software: you can redistribute it and/or modify
! it under the terms of the GNU Lesser General Public License
! as published by the Free Software Foundation, either version
! 3 of the License, or (at your option) any later version.
!
! PDAF 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 Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public
! License along with PDAF.  If not, see <http://www.gnu.org/licenses/>.
!
!$Id$
!BOP
!
! !ROUTINE: PDAF_pf_resampling --- Get particle indices for resampling
!
! !INTERFACE:
SUBROUTINE PDAF_pf_resampling(method, Nin, Nout, weights, IDs, screen)

! !DESCRIPTION:
! Determine particle indices for resampling. Implemented
! are three sampling schemes:
! (1) probabilistic resampling
! (2) stochastic unversal resampling
! (3) residual resampling
! The schemes follow the algorithms in 
! Vetra-Carvalho et al., State-of-the-art stochastic data
! assimilation methods for high-dimensional non-Gaussian problems.
! Tellus A, 70:1, 1445364, doi:10.1080/16000870.2018.1445364
!
! !  This is a core routine of PDAF and
!    should not be changed by the user   !
!
! !REVISION HISTORY:
! 2019-05 - Lars Nerger initial code
! Later revisions - see svn log
!
! !USES:
! Include definitions for real type of different precision
! (Defines BLAS/LAPACK routines and MPI_REALTYPE)
#include "typedefs.h"

  USE PDAF_timer, &
       ONLY: PDAF_timeit
  USE PDAF_memcounting, &
       ONLY: PDAF_memcount
  USE PDAF_mod_filtermpi, &
       ONLY: mype

  IMPLICIT NONE

! !ARGUMENTS:
  INTEGER, INTENT(in) :: method       ! Choose resampling method
                                      ! (1) probabilistic resampling
  INTEGER, INTENT(in) :: Nin          ! number of particles
  INTEGER, INTENT(in) :: Nout         ! number of particles to be resampled
  REAL, INTENT(in)    :: weights(Nin) ! Weights
  INTEGER, INTENT(out) :: IDs(Nout)   ! Indices of resampled ensmeble states
  INTEGER, INTENT(in) :: screen       ! Verbosity flag

! !CALLING SEQUENCE:
!EOP
       
! *** local variables ***
  INTEGER :: i, j                    ! Loop counters
  INTEGER :: c, Nr                   ! Counter for resampling
  INTEGER, SAVE :: first = 1         ! flag for init of random number seed
  INTEGER, SAVE :: iseed(4)          ! seed array for random number routine
  REAL :: rndval                     ! Random value
  REAL, ALLOCATABLE :: w_acc(:)      ! accumulated weights
  INTEGER, ALLOCATABLE :: w_i(:)     ! Integer weights
  REAL, ALLOCATABLE :: w_r(:)        ! residual weights


! **********************
! *** INITIALIZATION ***
! **********************

  IF (mype == 0 .AND. screen > 0) THEN
     WRITE (*, '(a, 5x, a)') &
          'PDAF', 'Resample particles'
     IF (method == 1) THEN
        WRITE (*, '(a, 5x, a)') &
             'PDAF', '--- Probabilistic resampling'
     ELSE IF (method == 2) THEN
        WRITE (*, '(a, 5x, a)') &
             'PDAF', '--- Stachastic universal resampling'
     ELSE IF (method == 3) THEN
        WRITE (*, '(a, 5x, a)') &
             'PDAF', '--- Residual resampling'
     END IF
  END IF

  IF (first==1) THEN
     ! Set random seed
     iseed(1)=15
     iseed(2)=25
     iseed(3)=47
     iseed(4)=17
     first = 0
  END IF

  ! Initialize accumulated weights
  ALLOCATE(w_acc(Nin))
  w_acc = 0.0

  ! Initialize rasampling IDs
  IDs = 0

   
  Rtype: IF (method == 1) THEN

! ********************************
! *** Probabilistic resampling ***
! ********************************

     ! Get accumulated weights
     w_acc(1) = weights(1)

     DO i = 2, Nin
        W_acc(i) = w_acc(i-1) + weights(i)
     END DO

     c = 1
     DO j = 1, Nout

        ! Init random number
        CALL larnvTYPE(1, iseed, 1, rndval)

        ! Determine index
        checkacc: DO i = 1, Nin
           IF (rndval > w_acc(i)) THEN
              c = c + 1
           ELSE
              EXIT checkacc
           END IF
        END DO checkacc

        IDs(j) = c
        c = 1

     END DO


  ELSE IF (method == 2) THEN

! ***************************************
! *** Stochastic universal resampling ***
! ***************************************

     ! Get accumulated weights
     w_acc(1) = weights(1)

     DO i = 2, Nin
        w_acc(i) = w_acc(i-1) + weights(i)
     END DO

     ! Init random number
     CALL larnvTYPE(1, iseed, 1, rndval)
     rndval = rndval / Nin


     c = 1

     DO j = 1, Nout

        ! Determine index
        checkaccB: DO i = 1, Nin
           IF (rndval > w_acc(i)) THEN
              c = c + 1
           ELSE
              EXIT checkaccB
           END if
        END DO checkaccB

        IDs(j) = c
        
        rndval = rndval + 1.0/REAL(Nin)

        c = 1

     END DO

  ELSE IF (method == 3) THEN

! ***************************
! *** Residual resampling ***
! ***************************

     ALLOCATE(w_i(Nin))
     ALLOCATE(w_r(Nin))

     DO i = 1, Nin
        w_i(i) = FLOOR(weights(i) * Nout)
     END DO

     ! First round setting resamping indices
     Nr = Nout
     c = 1
     IDs = 0
     DO j = 1, Nout
        IF (w_i(j) > 0) THEN
           IDs(c:c+w_i(j)-1) = j
           c = c + w_i(j)
           Nr = Nr - w_i(j)
        END IF
     END DO

     ! Now perform probabilistic resampling for remaining Nr indices

     ! Determine residual weights 

     IF (Nr > 0) THEN
        DO i = 1, Nin
           w_r(i) = (weights(i)*Nout - REAL(w_i(i))) / REAL (Nr)
        END DO
     ELSE
        w_r(:) = 0
     END IF

     ! Get accumulated weights
     w_acc(1) = w_r(1)

     DO i = 2, Nin
        w_acc(i) = w_acc(i-1) + w_r(i)
     END DO

     c = 1
     DO j = Nout-Nr+1, Nout

        ! Init random number
        CALL larnvTYPE(1, iseed, 1, rndval)
        rndval = rndval * REAL(Nout) / REAL(Nin)

        ! Determine index
        checkaccC: DO i = 1, Nin
           IF (rndval > w_acc(i)) THEN
              c = c + 1
           ELSE
              EXIT checkaccC
           END IF
        END DO checkaccC

        IDs(j) = c
        c = 1

     END DO

     DEALLOCATE(w_i, w_r)

    END IF Rtype


! ********************
! *** Finishing up ***
! ********************

  DEALLOCATE(w_acc)

END SUBROUTINE PDAF_pf_resampling