Only calculate the lower triangle (#91)

* Only the lower triangle needs to be calculated

* Make benchmarks prettier

* Wrap matrices in `Symmetric`

* Add fun stuff to distancematrix

Author: asinghvi17

benchmarks/parallel_vs_serial.jl

@@ -4,21 +4,71 @@ This file benchmarks the parallel and serial computation of a recurrence matrix
 using DynamicalSystemsBase, RecurrenceAnalysis
 using Statistics, BenchmarkTools, Base.Threads
 
+function pretty_time(b::BenchmarkTools.Trial)
+    med = median(b)
+
+    return string(
+        BenchmarkTools.prettytime(BenchmarkTools.time(med)),
+        " (",
+        BenchmarkTools.prettypercent(BenchmarkTools.gcratio(med)),
+        " GC)"
+    )
+
+end
+
 Ns = round.(Int, 10.0 .^ (2:0.5:4.5))
 ro = Systems.roessler()
 
-println("I am using $(nthreads()) threads. I am reporting results as:")
-println("time of (parallel/serial) for 3D and 1D trajectories.")
+println("I am using $(nthreads()) threads.")
 for N in Ns
+    # set up datasets
     tr = trajectory(ro, N*0.1; dt = 0.1, Tr = 10.0)
-    println("For N = $(length(tr))...")
+    printstyled("For N = $(length(tr))..."; color = :blue, bold = true)
+    println()
     x = tr[:, 1]
+
+    # calculate only the lower triangle
     b_serial   = @benchmark RecurrenceMatrix($(tr), 5.0; metric = Euclidean(), parallel = false)
     b_parallel = @benchmark RecurrenceMatrix($(tr), 5.0; metric = Euclidean(), parallel = true)
     b_s_1 = @benchmark RecurrenceMatrix($(x), 2.0; parallel = false)
     b_p_1 = @benchmark RecurrenceMatrix($(x), 2.0; parallel = true)
+
+    # calculate triangular metrics
     threeD = median(b_serial.times)/median(b_parallel.times)
     oneD = median(b_s_1.times)/median(b_p_1.times)
     threeD, oneD = round.((threeD, oneD); sigdigits = 3)
-    println("3D=$(threeD), 1D=$(oneD)")
+    printstyled("Lower triangle only"; color = :green)
+    println("\nSpeedups: \n    3D=$threeD\n    1D=$oneD")
+    println(
+        """
+        Raw timings:
+            3D serial:   $(pretty_time(b_serial))
+            3D parallel: $(pretty_time(b_parallel))
+            1D serial:   $(pretty_time(b_s_1))
+            1D parallel: $(pretty_time(b_p_1))
+        """
+    )
+
+    # calculate the full matrix
+    b_serial_full   = @benchmark CrossRecurrenceMatrix($(tr), $(tr), 5.0; metric = Euclidean(), parallel = false)
+    b_parallel_full = @benchmark CrossRecurrenceMatrix($(tr), $(tr), 5.0; metric = Euclidean(), parallel = true)
+    b_s_1_full = @benchmark CrossRecurrenceMatrix($(x), $(x), 2.0; parallel = false)
+    b_p_1_full = @benchmark CrossRecurrenceMatrix($(x), $(x), 2.0; parallel = true)
+
+    # calculate full metrics
+    threeD = median(b_serial_full.times)/median(b_parallel_full.times)
+    oneD = median(b_s_1_full.times)/median(b_p_1_full.times)
+    threeD, oneD = round.((threeD, oneD); sigdigits = 3)
+    printstyled("Full matrix"; color = :green)
+    println("\nSpeedups: \n    3D=$threeD\n    1D=$oneD")
+    println(
+        """
+        Raw timings:
+            3D serial:   $(pretty_time(b_serial_full))
+            3D parallel: $(pretty_time(b_parallel_full))
+            1D serial:   $(pretty_time(b_s_1_full))
+            1D parallel: $(pretty_time(b_p_1_full))
+        """
+    )
+    println()
 end

src/matrices.jl

@@ -32,36 +32,39 @@ The list of strings available to define the metric are:
 * `"manhattan"`, `"cityblock"`, `"taxicab"` or `"min"` for the Manhattan or L1 norm
   (`Cityblock()` in `Distances`).
 """
-distancematrix(x, metric::Union{Metric,String}=DEFAULT_METRIC) = distancematrix(x, x, metric)
+distancematrix(x, metric::Union{Metric,String}=DEFAULT_METRIC, parallel = size(x)[1] > 500) = _distancematrix(x, getmetric(metric), Val(parallel))
 
 # For 1-dimensional arrays (vectors), the distance does not depend on the metric
-distancematrix(x::Vector, y::Vector, metric=DEFAULT_METRIC) = abs.(x .- y')
+distancematrix(x::Vector, y::Vector, metric=DEFAULT_METRIC, parallel = size(x)[1] > 500) = abs.(x .- y')
 
 # If the metric is supplied as a string, get the corresponding Metric from Distances
-distancematrix(x, y, metric::String) = distancematrix(x, y, getmetric(metric))
+distancematrix(x, y, metric::String, parallel = size(x)[1] > 500) = distancematrix(x, y, getmetric(metric), parallel)
 
 const MAXDIM = 9
-function distancematrix(x::Tx, y::Ty, metric::Metric=DEFAULT_METRIC) where
+function distancematrix(x::Tx, y::Ty, metric::Metric=DEFAULT_METRIC, parallel = size(x)[1] > 500) where
          {Tx<:Union{AbstractMatrix, Dataset}} where {Ty<:Union{AbstractMatrix, Dataset}}
     sx, sy = size(x), size(y)
-    if sx[2] != sy[2]
-        error("the dimensions of `x` and `y` data points must be the equal")
-    end
+    @assert sx[2] == sy[2] """
+        The dimensions of the data points in `x` and `y` must be equal!
+        Found dim(x)=$(sx[2]), dim(y)=$(sy[2]).
+        """
+
     if sx[2] ≥ MAXDIM && typeof(metric) == Euclidean # Blas optimization
-        return _distancematrix(Matrix(x), Matrix(y), metric)
+        return _distancematrix(Matrix(x), Matrix(y), metric, Val(parallel))
     elseif Tx <: Matrix && Ty <: Matrix && metric == Chebyshev()
-        return _distancematrix(x, y, metric)
+        return _distancematrix(x, y, metric, Val(parallel))
     else
-        return _distancematrix(Dataset(x), Dataset(y), metric)
+        return _distancematrix(Dataset(x), Dataset(y), metric, Val(parallel))
     end
 end
 
 # Core function for Matrices (wrapper of `pairwise` from the Distances package)
-_distancematrix(x::AbstractMatrix, y::AbstractMatrix, metric::Metric) =
-pairwise(metric, x', y', dims=2)
+_distancematrix(x::AbstractMatrix, y::AbstractMatrix, metric::Metric, ::Val{false}) = pairwise(metric, x', y', dims=2)
+
+# First we define the serial versions of the functions.
 # Core function for Datasets
 function _distancematrix(x::Dataset{S,Tx}, y::Dataset{S,Ty},
-    metric::Metric) where {S, Tx, Ty}
+    metric::Metric, ::Val{false}) where {S, Tx, Ty}
 
     x = x.data
     y = y.data
@@ -74,6 +77,127 @@ function _distancematrix(x::Dataset{S,Tx}, y::Dataset{S,Ty},
     return d
 end
 
+# Now, we define the parallel versions.
+
+function _distancematrix(x::Dataset{S,Tx}, y::Dataset{S,Ty},
+    metric::Metric, ::Val{true}) where {S, Tx, Ty}
+
+    x = x.data
+    y = y.data
+    d = zeros(promote_type(Tx,Ty), length(x), length(y))
+    Threads.@threads for j in 1:length(y)
+        for i in 1:length(x)
+            @inbounds d[i,j] = evaluate(metric, x[i], y[j])
+        end
+    end
+    return d
+end
+
+function _distancematrix(x::Matrix{Tx}, y::Matrix{Ty},
+    metric::Metric, ::Val{true}) where {Tx, Ty}
+
+    x = x.data
+    y = y.data
+    d = zeros(promote_type(Tx,Ty), length(x), length(y))
+    Threads.@threads for j in 1:length(y)
+        for i in 1:length(x)
+            @inbounds d[i,j] = evaluate(metric, x[i, :], y[j, :])
+        end
+    end
+    return d
+end
+
+
+# Now, we define methods for the single trajectory.
+# We can speed this up by only calculating the lower triangle,
+# which halves the number of computations needed.
+
+# Again, we'll define the serial version first:
+function _distancematrix(x::Vector{T}, metric::Metric, ::Val{false}) where T
+    d = zeros(T, length(x), length(x))
+
+    for j in 1:length(x)
+        for i in 1:j
+            @inbounds d[i, j] = abs(x[i] - x[j])
+        end
+    end
+
+    return Symmetric(d, :U)
+
+end
+
+function _distancematrix(x::Dataset{S, T}, metric::Metric, ::Val{false}) where T where S
+    d = zeros(T, length(x), length(x))
+
+    for j in 1:length(x)
+        for i in 1:j
+            @inbounds d[i, j] = evaluate(metric, x[i], x[j])
+        end
+    end
+
+    return Symmetric(d, :U)
+
+end
+
+# Now, we define the parallel version.  There's a twist, though.
+
+# The single trajectory case is a little tricky.  If you try it naively,
+# using the method we used for the serial computation above, the points are
+# partioned out unequally between threads.  This means that performance actually
+# **decreases** compared to the full parallel version.  To mitigate this, we need
+# to partition the computation equally among all threads.
+# The function I've written below does essentially this - given a length,
+# it calculates the appropriate partitioning, and returns a list of indices,
+# by utilizing the fact that the area is proportional to the square of the height.
+# It partitions the "triangle" which needs to be computed into "trapezoids",
+# which all have an equal area.
+function partition_indices(len)
+    indices = Vector{UnitRange{Int}}(undef, Threads.nthreads())
+    length = len
+    offset = 1
+
+    # Here, we have to iterate in reverse, to get an equal area every time
+    for n in Threads.nthreads():-1:1
+        partition = round(Int, length / sqrt(n)) # length varies as square root of area
+        indices[n] = offset:(partition .+ offset .- 1)
+        length -= partition
+        offset += partition
+    end
+
+    return indices
+end
+
+# And now for the actual implementation:
+function _distancematrix(x::Vector{T}, metric::Metric, ::Val{true}) where T
+    d = zeros(T, length(x), length(x))
+
+    Threads.@threads for k in partition_indices(length(x))
+        for j in k
+            for i in 1:j
+                @inbounds d[i, j] = abs(x[i] - x[j])
+            end
+        end
+    end
+
+    return Symmetric(d, :U)
+
+end
+
+function _distancematrix(x::Dataset{S, T}, metric::Metric, ::Val{true}) where T where S
+    d = zeros(T, length(x), length(x))
+
+    Threads.@threads for k in partition_indices(length(x))
+        for j in k
+            for i in 1:j
+                @inbounds d[i, j] = evaluate(metric, x[i], x[j])
+            end
+        end
+    end
+
+    return Symmetric(d, :U)
+
+end
+
 ################################################################################
 # AbstractRecurrenceMatrix type definitions and documentation strings
 ################################################################################
@@ -308,12 +432,9 @@ Return a sparse matrix which encodes recurrence points.
 
 Note that `parallel` may be either `Val(true)` or `Val(false)`.
 """
-recurrence_matrix(x, metric::Metric, ε::Real, parallel::Val) =
-recurrence_matrix(x, x, metric, ε, parallel)
-
-# Convert Matrices to Datasets (better performance in all cases)
 function recurrence_matrix(x::AbstractMatrix, y::AbstractMatrix,
                            metric::Metric, ε, parallel::Val)
+    # Convert Matrices to Datasets (better performance in all cases)
     return recurrence_matrix(Dataset(x), Dataset(y), metric, ε, parallel)
 end
 
@@ -357,6 +478,42 @@ function recurrence_matrix(x::AbstractVector, y::AbstractVector, metric, ε, par
     return sparse(rowvals, colvals, nzvals, length(x), length(y))
 end
 
+function recurrence_matrix(x::AbstractVector, metric::Metric, ε::Real, parallel::Val{false})
+    rowvals = Vector{Int}()
+    colvals = Vector{Int}()
+    for j in 1:length(x)
+        nzcol = 0
+        for i in 1:j
+            if @inbounds abs(x[i] - x[j]) ≤ ε
+                push!(rowvals, i)
+                nzcol += 1
+            end
+        end
+        append!(colvals, fill(j, (nzcol,)))
+    end
+    nzvals = fill(true, (length(rowvals),))
+    return Symmetric(sparse(rowvals, colvals, nzvals, length(x), length(x)), :U)
+end
+
+function recurrence_matrix(xx::Dataset, metric::Metric, ε::Real, parallel::Val{false})
+    x = xx.data
+    rowvals = Vector{Int}()
+    colvals = Vector{Int}()
+    for j in 1:length(x)
+        nzcol = 0
+        for i in 1:j
+            @inbounds if evaluate(metric, x[i], x[j]) ≤ ε
+                push!(rowvals, i)
+                nzcol += 1
+            end
+        end
+        append!(colvals, fill(j, (nzcol,)))
+    end
+    nzvals = fill(true, (length(rowvals),))
+    return Symmetric(sparse(rowvals, colvals, nzvals, length(x), length(x)), :U)
+end
+
+
 # Now, we define the parallel versions of these functions.
 
 # Core function
@@ -412,3 +569,61 @@ function recurrence_matrix(x::AbstractVector, y::AbstractVector, metric, ε, par
     nzvals = fill(true, (length(finalrows),))
     return sparse(finalrows, finalcols, nzvals, length(x), length(y))
 end
+
+
+
+function recurrence_matrix(xx::AbstractVector, metric::Metric, ε::Real, parallel::Val{true})
+    x = xx
+    # We create an `Array` of `Array`s, for each thread to have its
+    # own array to push to.  This avoids race conditions with
+    # multiple threads pushing to the same `Array` (`Array`s are not atomic).
+    rowvals = [Vector{Int}() for _ in 1:Threads.nthreads()]
+    colvals = [Vector{Int}() for _ in 1:Threads.nthreads()]
+
+    # This is the same logic as the serial function, but parallelized.
+    Threads.@threads for k in partition_indices(length(x))
+        threadn = Threads.threadid()
+        for j in k
+            nzcol = 0
+            for i in 1:j
+                @inbounds if abs(x[i] - x[j]) ≤ ε
+                    push!(rowvals[threadn], i) # push to the thread-specific row array
+                    nzcol += 1
+                end
+            end
+            append!(colvals[threadn], fill(j, (nzcol,)))
+        end
+    end
+    finalrows = vcat(rowvals...) # merge into one array
+    finalcols = vcat(colvals...) # merge into one array
+    nzvals = fill(true, (length(finalrows),))
+    return Symmetric(sparse(finalrows, finalcols, nzvals, length(x), length(x)), :U)
+end
+
+function recurrence_matrix(xx::Dataset, metric::Metric, ε::Real, parallel::Val{true})
+    x = xx.data
+    # We create an `Array` of `Array`s, for each thread to have its
+    # own array to push to.  This avoids race conditions with
+    # multiple threads pushing to the same `Array` (`Array`s are not atomic).
+    rowvals = [Vector{Int}() for _ in 1:Threads.nthreads()]
+    colvals = [Vector{Int}() for _ in 1:Threads.nthreads()]
+
+    # This is the same logic as the serial function, but parallelized.
+    Threads.@threads for k in partition_indices(length(x))
+        threadn = Threads.threadid()
+        for j in k
+            nzcol = 0
+            for i in 1:j
+                @inbounds if evaluate(metric, x[i], x[j]) ≤ ε
+                    push!(rowvals[threadn], i) # push to the thread-specific row array
+                    nzcol += 1
+                end
+            end
+            append!(colvals[threadn], fill(j, (nzcol,)))
+        end
+    end
+    finalrows = vcat(rowvals...) # merge into one array
+    finalcols = vcat(colvals...) # merge into one array
+    nzvals = fill(true, (length(finalrows),))
+    return Symmetric(sparse(finalrows, finalcols, nzvals, length(x), length(x)), :U)
+end