Skip to content

Commit

Permalink
comment changes, formatted document
Browse files Browse the repository at this point in the history
  • Loading branch information
PGS62 committed Jan 27, 2021
1 parent 8a59392 commit fa98fe7
Showing 1 changed file with 19 additions and 20 deletions.
39 changes: 19 additions & 20 deletions src/rankcorr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,8 @@ function corkendall!(x::RealVector, y::RealVector, permx=sortperm(x))
x[:] = x[permx]
y[:] = y[permx]

npairs = float(n) * (n - 1)/ 2
#ntiesx, ntiesy, ndoubleties are floats to avoid overflows on 32bit
npairs = float(n) * (n - 1) / 2
# ntiesx, ntiesy, ndoubleties are floats to avoid overflows on 32bit
ntiesx, ntiesy, ndoubleties, k, nswaps = 0.0, 0.0, 0.0, 0, 0

@inbounds for i 2:n
Expand All @@ -31,18 +31,18 @@ function corkendall!(x::RealVector, y::RealVector, permx=sortperm(x))
# sorted first on x, then (where x values are tied) on y. Hence
# double ties can be counted by calling countties.
sort!(view(y, (i - k - 1):(i - 1)))
ntiesx += float(k) * (k + 1)/2
ntiesx += float(k) * (k + 1) / 2
ndoubleties += countties(y, i - k - 1, i - 1)
k = 0
end
end
if k > 0
sort!(view(y, ((n - k):n)))
ntiesx += float(k) * (k + 1)/ 2
ntiesx += float(k) * (k + 1) / 2
ndoubleties += countties(y, n - k, n)
end

nswaps = mergesort!(y, 1, n)
nswaps = msort!(y, 1, n)
ntiesy = countties(y, 1, n)

(npairs + ndoubleties - ntiesx - ntiesy - 2 * nswaps) /
Expand All @@ -55,20 +55,20 @@ end
Assumes `x` is sorted. Returns the number of ties within `x[lo:hi]`.
"""
function countties(x::AbstractVector, lo::Integer, hi::Integer)
#avoid ovorflows on 32 bit by using floats
# avoid overflows on 32 bit by using floats
thistiecount, result = 0.0, 0.0
(lo < 1 || hi > length(x)) && error("Bounds error")
@inbounds for i (lo + 1):hi
if x[i] == x[i - 1]
thistiecount += 1.0
elseif thistiecount > 0
result += thistiecount * (thistiecount + 1)/2
result += thistiecount * (thistiecount + 1) / 2
thistiecount = 0.0
end
end

if thistiecount > 0
result += thistiecount * (thistiecount + 1)/2
result += thistiecount * (thistiecount + 1) / 2
end
result
end
Expand Down Expand Up @@ -113,29 +113,29 @@ end

# Auxilliary functions for Kendall's rank correlation

# Tests appeared to show that a value of 64 is optimal,
# Tests appear to show that a value of 64 is optimal,
# but note that the equivalent constant in base/sort.jl is 20.
const SMALL_THRESHOLD = 64

# Copy was from https://github.com/JuliaLang/julia/commit/28330a2fef4d9d149ba0fd3ffa06347b50067647 dated 20 Sep 2020
"""
mergesort!(v::AbstractVector, lo::Integer, hi::Integer, t=similar(v, 0))
msort!(v::AbstractVector, lo::Integer, hi::Integer, t=similar(v, 0))
Mutates `v` by sorting elements `x[lo:hi]` using the mergesort algorithm.
Mutates `v` by sorting elements `x[lo:hi]` using the merge sort algorithm.
This method is a copy-paste-edit of sort! in base/sort.jl (the method specialised on MergeSortAlg),
but amended to return the bubblesort distance.
"""
function mergesort!(v::AbstractVector, lo::Integer, hi::Integer, t=similar(v, 0))
#avoid overflow errors (if length(v)> 2^16) on 32 bit by using float
function msort!(v::AbstractVector, lo::Integer, hi::Integer, t=similar(v, 0))
# avoid overflow errors (if length(v)> 2^16) on 32 bit by using float
nswaps = 0.0
@inbounds if lo < hi
hi - lo <= SMALL_THRESHOLD && return insertionsort!(v, lo, hi)
hi - lo <= SMALL_THRESHOLD && return isort!(v, lo, hi)

m = midpoint(lo, hi)
(length(t) < m - lo + 1) && resize!(t, m - lo + 1)

nswaps = mergesort!(v, lo, m, t)
nswaps += mergesort!(v, m + 1, hi, t)
nswaps = msort!(v, lo, m, t)
nswaps += msort!(v, m + 1, hi, t)

i, j = 1, lo
while j <= m
Expand Down Expand Up @@ -165,20 +165,19 @@ function mergesort!(v::AbstractVector, lo::Integer, hi::Integer, t=similar(v, 0)
return nswaps
end

# This implementation of `midpoint` is performance-optimized but safe only if `lo <= hi`.
# This function is also copied from base/sort.jl
midpoint(lo::T, hi::T) where T <: Integer = lo + ((hi - lo) >>> 0x01)
midpoint(lo::Integer, hi::Integer) = midpoint(promote(lo, hi)...)

# Copy was from https://github.com/JuliaLang/julia/commit/28330a2fef4d9d149ba0fd3ffa06347b50067647 dated 20 Sep 2020
"""
insertionsort!(v::AbstractVector, lo::Integer, hi::Integer)
isort!(v::AbstractVector, lo::Integer, hi::Integer)
Mutates `v` by sorting elements `x[lo:hi]` using the insertionsort algorithm.
Mutates `v` by sorting elements `x[lo:hi]` using the insertion sort algorithm.
This method is a copy-paste-edit of sort! in base/sort.jl (the method specialised on InsertionSortAlg),
amended to return the bubblesort distance.
"""
function insertionsort!(v::AbstractVector, lo::Integer, hi::Integer)
function isort!(v::AbstractVector, lo::Integer, hi::Integer)
if lo == hi return 0.0 end
nswaps = 0.0
@inbounds for i = lo + 1:hi
Expand Down

0 comments on commit fa98fe7

Please sign in to comment.