选择

stdlib_selection 模块

选择概述

假设您希望查找大小为 N 的数组中第 k 小的条目的值,或该值的索引。虽然可以通过使用 sortsort_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))) 为真。用户可以选择指定 leftright 索引以限制对第 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)))。用户可以选择指定 leftright 索引以限制对第 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 元素;它会运行,但与其他条目相比,对 arrayNaN 条目的顺序没有一致的解释。

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 元素;它会运行,但与其他条目相比,对 arrayNaN 条目的顺序没有一致的解释。

虽然 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 的比较

以下程序比较了 selectarg_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