stdlib_selection 模块假设您希望查找大小为 N 的数组中第 k 小的条目的值,或该值的索引。虽然可以通过使用 sort 或 sort_index 来自 stdlib_sorting 对整个数组进行排序,然后查找第 k 个条目来完成,但这需要 O(N x LOG(N)) 的时间。但是,单个条目的选择可以在 O(N) 时间内完成,对于大型数组来说,这要快得多。例如,这对于快速查找数组的中位数或其他百分位数很有用。
因此,Fortran 标准库提供了一个模块 stdlib_selection,它实现了选择算法。
模块 stdlib_selection 定义了两个泛型子例程
select 用于查找数组的第 k 小的条目。输入数组也会就地修改,并在返回时将被部分排序,使得 all(array(1:k) <= array(k))) 和 all(array(k) <= array((k+1):size(array))) 为真。用户可以选择指定 left 和 right 索引以限制对第 k 小值的搜索。如果您之前已调用 select 来查找更小或更大的秩(这将导致 array 的部分排序,从而暗示某些约束条件),这将很有用。
arg_select 用于查找数组的第 k 小的条目的索引。在这种情况下,输入数组不会被修改,但用户必须提供一个与 array 大小相同的输入索引数组,其索引是 1:size(array) 的排列,该数组将被修改。返回后,索引数组将被修改,使得 all(array(index(1:k)) <= array(index(k))) 和 all(array(k) <= array(k+1:size(array)))。用户可以选择指定 left 和 right 索引以限制对第 k 小值的搜索。如果您之前已调用 arg_select 来查找更小或更大的秩(这将导致 index 的部分排序,从而暗示某些约束条件),这将很有用。
select - 在输入数组中查找第 k 小的值实验性
返回 array(:) 的第 k 小的值,并部分排序 array(:),使得 all(array(1:k) <= array(k)) 和 all(array(k) <= array((k+1):size(array)))
call select ( array, k, kth_smallest [, left, right ] )
泛型子例程。
array:应为以下任意类型的秩一数组:integer(int8)、integer(int16)、integer(int32)、integer(int64)、real(sp)、real(dp)、real(xdp)、real(qp)。它是一个 intent(inout) 参数。
k:应为以下任意类型的标量:integer(int8)、integer(int16)、integer(int32)、integer(int64)。它是一个 intent(in) 参数。我们搜索 array(:) 的第 k 小的条目。
kth_smallest:应为与 array 类型相同的标量。它是一个 intent(out) 参数。返回时,它包含 array(:) 的第 k 小的条目。
left(可选):应为与 k 类型相同的标量。它是一个 intent(in) 参数。如果指定,则我们假设第 k 小的值肯定包含在 array(left:size(array)) 中。如果 left 不存在,则默认为 1。如果进行了多次 select 调用,这通常很有用,因为 array 的部分排序意味着我们需要搜索的位置存在约束。
right(可选):应为与 k 类型相同的标量。它是一个 intent(in) 参数。如果指定,则我们假设第 k 小的值肯定包含在 array(1:right) 中。如果 right 不存在,则默认为 size(array)。如果进行了多次 select 调用,这通常很有用,因为 array 的部分排序意味着我们需要搜索的位置存在约束。
单个值的选取应具有 O(size(array)) 的运行时,因此它在渐近意义上比完全排序 array 快。本文档末尾的测试程序表明情况确实如此。
代码不支持 array 中的 NaN 元素;它会运行,但与其他条目相比,对 array 的 NaN 条目的顺序没有一致的解释。
select 源自 Leon Foks 编写的 Coretran 库中的代码,https://github.com/leonfoks/coretran。Leon Foks 已许可将此处的代码在 stdlib 的 MIT 许可证下发布。
program example_select
use stdlib_selection, only: select
implicit none
real, allocatable :: array(:)
real :: kth_smallest
integer :: k
array = [3., 2., 7., 4., 5., 1., 4., -1.]
k = 2
call select(array, k, kth_smallest)
print *, kth_smallest ! print 1.0
k = 7
! Due to the previous call to select, we know for sure this is in an
! index >= 2
call select(array, k, kth_smallest, left=2)
print *, kth_smallest ! print 5.0
k = 6
! Due to the previous two calls to select, we know for sure this is in
! an index >= 2 and <= 7
call select(array, k, kth_smallest, left=2, right=7)
print *, kth_smallest ! print 4.0
end program example_select
arg_select - 在输入数组中查找第 k 小的值的索引实验性
返回 array(:) 的第 k 小的值的索引,并部分排序索引数组 indx(:),使得 all(array(indx(1:k)) <= array(indx(k))) 和 all(array(indx(k)) <= array(indx((k+1):size(array))))
call arg_select ( array, indx, k, kth_smallest [, left, right ] )
泛型子例程。
array:应为以下任意类型的秩一数组:integer(int8)、integer(int16)、integer(int32)、integer(int64)、real(sp)、real(dp)、real(xdp)、real(qp)。它是一个 intent(in) 参数。在输入时,它是我们在其中搜索第 k 小的条目的数组。
indx:应为与 array 大小相同的秩一数组,包含任意顺序的 1:size(array) 中的所有整数。它是以下任意类型:integer(int8)、integer(int16)、integer(int32)、integer(int64)。它是一个 intent(inout) 参数。返回时,其元素将定义 array(:) 的部分排序,使得:all( array(indx(1:k-1)) <= array(indx(k)) ) 和 all(array(indx(k)) <= array(indx(k+1:size(array))))。
k:应为与 indx 类型相同的标量。它是一个 intent(in) 参数。我们搜索 array(:) 的第 k 小的条目。
kth_smallest:与 indx 类型相同的标量。它是一个 intent(out) 参数,返回时,它包含 array(:) 的第 k 小的条目的索引。
left(可选):应为与 k 类型相同的标量。它是一个 intent(in) 参数。如果指定,则我们假设第 k 小的值肯定包含在 array(indx(left:size(array))) 中。如果 left 不存在,则默认为 1。如果进行了多次 arg_select 调用,这通常很有用,因为 indx 的部分排序意味着我们需要搜索的位置存在约束。
right(可选):应为与 k 类型相同的标量。它是一个 intent(in) 参数。如果指定,则我们假设第 k 小的值肯定包含在 array(indx(1:right)) 中。如果 right 不存在,则默认为 size(array)。如果进行了多次 arg_select 调用,这通常很有用,因为 indx 的重新排序意味着我们需要搜索的位置存在约束。
arg_select 不会修改 array,这与 select 不同。
indx 的部分排序不稳定,即映射到数组相等值的索引可能会重新排序。
代码不支持 array 中的 NaN 元素;它会运行,但与其他条目相比,对 array 的 NaN 条目的顺序没有一致的解释。
虽然 indx 包含 1:size(array) 的整数排列是必要的,但代码不会检查这一点。例如,如果 size(array) == 4,则我们可以有 indx = [4, 2, 1, 3] 或 indx = [1, 2, 3, 4],但不能有 indx = [2, 1, 2, 4]。避免此类错误是用户的责任。
单个值的选取应具有 O(size(array)) 的运行时,因此它在渐近意义上比完全排序 array 快。这些文档末尾的测试程序证实了这一点。
arg_select 是使用 Leon Foks 编写的 Coretran 库中的代码派生的,https://github.com/leonfoks/coretran。Leon Foks 已许可将此处的代码在 stdlib 的 MIT 许可证下发布。
program example_arg_select
use stdlib_selection, only: arg_select
implicit none
real, allocatable :: array(:)
integer, allocatable :: indx(:)
integer :: kth_smallest
integer :: k
array = [3., 2., 7., 4., 5., 1., 4., -1.]
indx = [(k, k=1, size(array))]
k = 2
call arg_select(array, indx, k, kth_smallest)
print *, array(kth_smallest) ! print 1.0
k = 7
! Due to the previous call to arg_select, we know for sure this is in an
! index >= 2
call arg_select(array, indx, k, kth_smallest, left=2)
print *, array(kth_smallest) ! print 5.0
k = 6
! Due to the previous two calls to arg_select, we know for sure this is in
! an index >= 2 and <= 7
call arg_select(array, indx, k, kth_smallest, left=2, right=7)
print *, array(kth_smallest) ! print 4.0
end program example_arg_select
sort 的比较以下程序比较了 select 和 arg_select 计算数组中位数的时间,以及使用 stdlib 中的 sort。理论上,我们应该看到选择例程的速度提高,这类似于 LOG(size(array))。
program selection_vs_sort
use stdlib_kinds, only: dp, sp, int64
use stdlib_selection, only: select, arg_select
use stdlib_sorting, only: sort
implicit none
call compare_select_sort_for_median(1)
call compare_select_sort_for_median(11)
call compare_select_sort_for_median(101)
call compare_select_sort_for_median(1001)
call compare_select_sort_for_median(10001)
call compare_select_sort_for_median(100001)
contains
subroutine compare_select_sort_for_median(N)
integer, intent(in) :: N
integer :: i, k, result_arg_select, indx(N), indx_local(N)
real :: random_vals(N), local_random_vals(N)
integer, parameter :: test_reps = 100
integer(int64) :: t0, t1
real :: result_sort, result_select
integer(int64) :: time_sort, time_select, time_arg_select
logical :: select_test_passed, arg_select_test_passed
! Ensure N is odd
if (mod(N, 2) /= 1) stop
time_sort = 0
time_select = 0
time_arg_select = 0
select_test_passed = .true.
arg_select_test_passed = .true.
indx = (/(i, i=1, N)/)
k = (N + 1)/2 ! Deliberate integer division
do i = 1, test_reps
call random_number(random_vals)
! Compute the median with sorting
local_random_vals = random_vals
call system_clock(t0)
call sort(local_random_vals)
result_sort = local_random_vals(k)
call system_clock(t1)
time_sort = time_sort + (t1 - t0)
! Compute the median with selection, assuming N is odd
local_random_vals = random_vals
call system_clock(t0)
call select(local_random_vals, k, result_select)
call system_clock(t1)
time_select = time_select + (t1 - t0)
! Compute the median with arg_select, assuming N is odd
local_random_vals = random_vals
indx_local = indx
call system_clock(t0)
call arg_select(local_random_vals, indx_local, k, result_arg_select)
call system_clock(t1)
time_arg_select = time_arg_select + (t1 - t0)
if (result_select /= result_sort) select_test_passed = .FALSE.
if (local_random_vals(result_arg_select) /= result_sort) arg_select_test_passed = .FALSE.
end do
print *, "select ; N=", N, '; ', merge('PASS', 'FAIL', select_test_passed), &
'; Relative-speedup-vs-sort:', (1.0*time_sort)/(1.0*time_select)
print *, "arg_select; N=", N, '; ', merge('PASS', 'FAIL', arg_select_test_passed), &
'; Relative-speedup-vs-sort:', (1.0*time_sort)/(1.0*time_arg_select)
end subroutine
end program
当 array 很大时,结果似乎与预期一致;程序打印
select ; N= 1 ; PASS; Relative-speedup-vs-sort: 1.90928173
arg_select; N= 1 ; PASS; Relative-speedup-vs-sort: 1.76875830
select ; N= 11 ; PASS; Relative-speedup-vs-sort: 1.14835048
arg_select; N= 11 ; PASS; Relative-speedup-vs-sort: 1.00794709
select ; N= 101 ; PASS; Relative-speedup-vs-sort: 2.31012774
arg_select; N= 101 ; PASS; Relative-speedup-vs-sort: 1.92877376
select ; N= 1001 ; PASS; Relative-speedup-vs-sort: 4.24190664
arg_select; N= 1001 ; PASS; Relative-speedup-vs-sort: 3.54580402
select ; N= 10001 ; PASS; Relative-speedup-vs-sort: 5.61573362
arg_select; N= 10001 ; PASS; Relative-speedup-vs-sort: 4.79348087
select ; N= 100001 ; PASS; Relative-speedup-vs-sort: 7.28823519
arg_select; N= 100001 ; PASS; Relative-speedup-vs-sort: 6.03007460