leonfoks
commented
3 years ago First stab. Also, @jvdp1 should we add the proposal label?
FWIW. I would be up for doing this, I have yet to contribute directly to stdlib because I have not had the mental capacity to overcome my learning curve with fypp, but this is an easy enough example to get my toes wet and stuck in. Once we have this, I can work on an implementation of Shewchuks adaptive floating point predicates (exact computation of the point inside/on/outside circle problem) and then KdTrees.
Description Hoare's O(n) quick select algorithm to find the kth sorted entry in an unsorted array, without performing the full sort. We can develop two types of algorithm based on this. quick_select and arg_quick_select (or ord_quick_select, but former matches numpy argmin, argpartition, argsort etc.).
Benefits Can be used to quickly find the median, or perhaps the value of the kth percentile in an array. The quick select is O(n) and in-place so is faster than first sorting the array and then choosing the kth value. This will be very useful for KdTrees where at each iteration, the dimension with the highest variance is split at its median value to create nicely balanced trees. A fast median (quick_select) here is key, and also a fast arg_median (and hence arg_quick_select).
API
I'll put what I have in Coretran, this will change to match the stdlib format.
Also, I don't know how to use fypp, or the syntax, but ill try to use suitable placeholders for possible changeable declarations.
The implementations will go in a submodule. By
interface select
!!Use an in-place quick select on an array of numbers
module subroutine quick_select_<array_type>_<k_type>(array, k, out, left, right)
!! Interfaced with select()
<array_type>, intent(inout) :: array(:)
!! Array to choose kth smallest from
<k_type>, intent(in) :: k
!! kth smallest element, where 1 <= k <= size(array)
<array_type> :: out
!! Return the value of the kth element
<k_type>, intent(in), optional :: left
!! If array has already been partially sorted, use this as the starting left index
<k_type>, intent(in), optional :: right
!! If array has already been partially sorted, use this as the starting right index
end subroutine
end interface
interface arg_select
!!Use an in-place quick select on an array of numbers using an index into the array.
module subroutine arg_quick_select_(array, indx, k, out, left, right)
!! Interfaced with argSelect()
<array_type>, intent(in) :: array(:)
!! Array to choose kth smallest from
<k_type>, intent(inout) :: indx(size(array))
!! Index into array.
<k_type>, intent(in) :: k
!! kth smallest element, where 1 <= k <= size(array)
<array_type> :: out
!! Return the value of the kth element
<k_type>, intent(in), optional :: left
!! If array has already been partially sorted, use this as the starting left index
<k_type>, intent(in), optional :: right
!! If array has already been partially sorted, use this as the starting right index
end subroutine
end interfaceExample
program select
use m_select, only: select
real(r64), allocatable :: d1D(:)
real(64) :: value
integer(i64), allocatable :: i1D(:)
integer(i64) :: ivalue
allocate(d1D(10000))
d1D = Assign values
k = (1 + size(d1D)) / 2
value = select(d1D, k)
write(*,'(a)') 'kth element? '//str(value)
allocate(i1D(10000))
i1D = Assign integers
ivalue = select(i1D, k)
write(*,'(a)') 'kth element? '//str(ivalue)
end programCurrent implementations I have an implementation in Coretran, and during a discussion that started here @gareth-nx brought up the Matlab version with a permissive license, @gareth-nx did a quick translation and ill attach that below.
Timing I timed the matlab version and my coretran version using the same arrays and compile flags. Timings are below and are very close with coretran being surprisingly faster for sorted arrays but maybe this is just timing fluctuations.
It might be easiest at this point to simply translate the matlab code, and add the arg_select variant which should be easy enough.
Left column is array length, right is time in seconds. Timing was done around a loop of 10 quick selections and the average taken. The array used had samples from a random normal.
==========================
Testing : Selection
==========================
Double precision quick select
16777216 1.694E-01
8388608 9.460E-02
4194304 5.000E-02
2097152 2.530E-02
1048576 1.030E-02
524288 3.500E-03
Already sorted double precision quick select
16777216 1.030E-02
8388608 4.500E-03
4194304 2.800E-03
2097152 1.000E-03
1048576 6.000E-04
524288 2.000E-04
==========================
Testing : Matlab Selection
==========================
Double precision quick select
16777216 1.579E-01
8388608 8.530E-02
4194304 4.260E-02
2097152 2.130E-02
1048576 9.300E-03
524288 4.100E-03
Already sorted double precision quick select
16777216 1.310E-02
8388608 6.200E-03
4194304 3.600E-03
2097152 1.500E-03
1048576 1.200E-03
524288 4.000E-04Matlab translated version
module qselect_mod
!
! This code is a translation of a matlab implementation "qselect" by Manolis Lourakis
! https://www.mathworks.com/matlabcentral/fileexchange/68947-qselect
! Below is the licence of qselect
!
!
! Copyright (c) 2018, Manolis Lourakis
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! * Redistributions of source code must retain the above copyright notice, this
! list of conditions and the following disclaimer.
!
! * Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution
! * Neither the name of Foundation for Research and Technology - Hellas nor the names of its
! contributors may be used to endorse or promote products derived from this
! software without specific prior written permission.
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
use iso_c_binding, only: C_DOUBLE, C_INT
implicit none
integer, parameter :: dp = C_DOUBLE, ip = C_INT
contains
function qselect(a, k) result(kth_largest)
! qselect - select the k-th smallest out of n numbers
! Implements Hoare s Quickselect algorithm (https://en.wikipedia.org/wiki/Quickselect)
! with the median-of-3 pivot strategy. Operates in-place, avoiding sorting and recursion
!
! kth = qselect(a, k)
!
! a - array of n elements
! k - specifies the desired k-th smallest element
! the k-th *largest* element can be found by passing length(a)+1-k
!
! Returns
!
! kth - the sought element
! Manolis Lourakis 2007-18
! Institute of Computer Science, Foundation for Research & Technology - Hellas
! Heraklion, Crete, Greece
! Sep 2018 - Original version. (v. 1.0)
!
real(dp), intent(inout) :: a(:)
integer(ip), intent(in) :: k
real(dp) :: kth_largest
integer(ip) :: l, r, s, i, j, k_local
real(dp) :: temp, pivot
l = 1_ip
r = size(a)
if(k < 1 .or. k > r) stop "in qselect, k must be between 1 and size(a)";
k_local = k
do while(.true.)
s = (l+r)/2_ip ! Deliberate integer division
if(a(s) < a(r)) then
temp = a(s); a(s) = a(r); a(r) = temp
end if
if(a(s) < a(l)) then
temp = a(s); a(s) = a(l); a(l) = temp
end if
if(a(r) < a(l)) then
temp = a(r); a(r) = a(l); a(l) = temp;
end if
pivot = a(r) ! Median
i = l
do j = l, r-1
if(a(j) <= pivot) then
temp = a(i); a(i) = a(j); a(j) = temp
i = i+1_ip
end if
end do
temp = a(r); a(r) = a(i); a(i) = temp
s = i-l+1
if(k_local < s) then
r = i-1_ip
else if(k_local > s) then
l=i+1_ip; k_local=k_local-s;
else
kth_largest = a(i)
return
end if
end do
end function
end module
gareth-nx
jvdp1
Beliavsky
ghbrown
In #426, selection algorithms were mentioned. See e.g., these comment, comment, and comment.
Description
To be completed
Prior Art
To be completed