fbin.f : Fast equally-populated bins using the quickselect algorithm
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
module fbin | |
use iso_c_binding | |
implicit none | |
contains | |
subroutine eqpop_bin(X, Nb, Ntrl, Xout) | |
! ARG | |
integer, intent(in) :: Nb, Ntrl | |
real(c_double), intent(in) :: X(Ntrl) | |
integer(c_int8_t), intent(out) :: Xout(Ntrl) | |
! LOC | |
integer :: ks(Nb-1) | |
real(c_double) :: arr(Ntrl) | |
integer :: idx(Ntrl), i, j, left, right, startidx | |
! copy for in place manipulation | |
arr = X | |
! calculate kth elements required | |
do i=1,(Nb-1) | |
ks(i) = nint( (Ntrl*i) / real(Nb) ) | |
end do | |
! fill index | |
do i=1,Ntrl | |
idx(i) = i | |
end do | |
left = 1 | |
right = Ntrl | |
call selectkth(ks(1), arr, idx, left, right) | |
do i=2,(Nb-1) | |
left = ks(i-1) + 1; | |
call selectkth(ks(i), arr, idx, left, right) | |
end do | |
! now assign bin labels | |
startidx = 1 | |
do i=1,(Nb-1) | |
do j=startidx,ks(i) | |
Xout(idx(j)) = i-1 | |
end do | |
startidx = ks(i) + 1 | |
end do | |
do i=startidx,Ntrl | |
Xout(idx(i)) = Nb - 1 | |
end do | |
end subroutine eqpop_bin | |
subroutine selectkth(k, arr, idx, left_in, right_in) | |
! ARG | |
real(c_double), intent(inout) :: arr(:) | |
integer, intent(inout) :: idx(:) | |
integer, intent(in) :: k, left_in, right_in | |
! LOC | |
integer :: r, w, mididx, tmp_i, left, right | |
real(c_double) :: piv, tmp_r | |
left = left_in | |
right = right_in | |
! median of three approach to choose pivot | |
do while (left < right) | |
r = left | |
w = right | |
mididx = floor( (r+w) / 2.0 ) | |
! sort with three conditional swaps | |
if( mididx > r ) then | |
if( arr(mididx) < arr(r) ) then | |
tmp_r = arr(mididx); arr(mididx) = arr(r); arr(r) = tmp_r | |
tmp_i = idx(mididx); idx(mididx) = idx(r); idx(r) = tmp_i | |
endif | |
if( arr(w) < arr(mididx) ) then | |
tmp_r = arr(mididx); arr(mididx) = arr(w); arr(w) = tmp_r | |
tmp_i = idx(mididx); idx(mididx) = idx(w); idx(w) = tmp_i | |
endif | |
if( arr(mididx) < arr(r) ) then | |
tmp_r = arr(mididx); arr(mididx) = arr(r); arr(r) = tmp_r | |
tmp_i = idx(mididx); idx(mididx) = idx(r); idx(r) = tmp_i | |
endif | |
endif | |
piv = arr(mididx) | |
do while (r < w) | |
if( arr(r) >= piv ) then | |
tmp_r = arr(w); arr(w) = arr(r); arr(r) = tmp_r | |
tmp_i = idx(w); idx(w) = idx(r); idx(r) = tmp_i | |
w = w - 1 | |
else | |
r = r + 1 | |
endif | |
end do | |
if( arr(r) > piv ) then | |
r = r - 1 | |
endif | |
if( k <= r ) then | |
right = r; | |
else | |
left = r+1; | |
endif | |
end do | |
end subroutine selectkth | |
subroutine eqpop_bin_sorted(X, Nb, Ntrl, Xout) | |
! ARG | |
integer, intent(in) :: Nb, Ntrl | |
real(c_double), intent(in) :: X(Ntrl) | |
integer(c_int16_t), intent(out) :: Xout(Ntrl) | |
! LOC | |
integer :: ks(Nb-1) | |
integer :: i, j, startidx | |
! calculate kth elements required | |
do i=1,(Nb-1) | |
ks(i) = nint( (Ntrl*i) / real(Nb) ) | |
end do | |
! assign bin labels | |
startidx = 1 | |
do i=1,(Nb-1) | |
do j=startidx,ks(i) | |
Xout(j) = i-1 | |
end do | |
startidx = ks(i) + 1 | |
end do | |
do i=startidx,Ntrl | |
Xout(i) = Nb - 1 | |
end do | |
end subroutine eqpop_bin_sorted | |
end module fbin |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment