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