Ported from Unpack.jl

.gitignore

@@ -1,3 +1,5 @@
 *.jl.cov
 *.jl.*.cov
 *.jl.mem
+*.ipynb
+src/scratch/.ipynb_checkpoints/*

LICENSE.md

@@ -1,6 +1,6 @@
 The Destruct.jl package is licensed under the MIT "Expat" License:
 
-> Copyright (c) 2018: SAM-VI.
+> Copyright (c) 2018: Samuel Palato.
 >
 >
 > Permission is hereby granted, free of charge, to any person obtaining a copy

README.md

@@ -1 +1,53 @@
-# Destruct
+# Destruct.jl
+Destructuring arrays of tuples in Julia (v0.6).
+
+## Overview
+Using julia's 'dot-call' syntax on functions with multiple return arguments
+results in an array of tuples. Sometimes, you want the tuple of arrays instead,
+preserving array shape.
+This can be achieved using `destruct`, which converts an array of tuple to a
+tuple of arrays.
+
+Works with any tuples (ie: with elements of different types).
+
+## Example
+```julia
+julia> using Destruct; using BenchmarkTools
+julia> f(a, b) = a+1im*b, a*b, convert(Int, round(a-b)); # some transform returing multiple values
+julia> v = f.(rand(3,1), rand(1,4));
+julia> typeof(v)
+Array{Tuple{Complex{Float64},Float64,Int64},2}
+julia> x, y, z = destruct(v);
+julia> z
+3×4 Array{Int64,2}:
+ 0  0  0  0
+ 1  0  1  1
+ 1  0  1  1
+julia> v = f.(rand(500,1,1), rand(1,500,500));
+julia> @btime destruct($v); # using BenchmarkTools
+  1.396 s (7 allocations: 3.73 GiB)
+```
+Getting this out of the way:
+```julia
+julia> x, y, z = f.(rand(100,1,1), rand(1,100,100)) |> destruct;
+```
+## Performance
+A common way to unpack Arrays of tuples uses the broadcast dot-call:
+```julia
+unpack_broadcast(w::Array{<:Tuple}) = Tuple((v->v[i]).(w) for i=1:length(w[1]))
+```
+However, this approach suffers from two problems: it doesn't access the elements
+in the order they are stored in memory and has huge memory consumption for
+Tuples with varying types (`Tuples` instead of `NTuples`).
+
+This "broadcast unpack" takes between 1.5x and 2x longer than `destruct`
+supplied here for arrays of `NTuples`. The performance gain is much larger 
+for tuples of heterogenous types: in the 10x to 75x range, using 1/10th
+of the memory.
+
+See timing scripts: `timing.jl` and `comparative_timing.jl`.
+
+## How does it work?
+The `destruct` function uses macros from `Base.Cartesian` to allocate
+destination arrays and iterate over all the things. The alternative
+implementations using broadcast dot-call is available as `Destruct.unpack_broadcast`.

src/Destruct.jl

@@ -1,5 +1,55 @@
 module Destruct
+export destruct
 
-# package code goes here
+using Base.Cartesian
 
-end # module
+unpack_broadcast(w::Array{<:Tuple}) = Tuple((v->v[i]).(w) for i=1:length(w[1]))
+
+"""
+    destruct(v::Array{<:Tuple,N})
+
+Destructure an array of tuples to a tuple of arrays. Works for tuples with
+elements of varying types.
+## Examples
+```julia-repl
+julia> f(a, b) = a+b, a*b, a-b;
+julia> v = f.(rand(3,1), rand(1,4));
+julia> x, y, z = destruct(v);
+julia> x
+3×4 Array{Float64,2}:
+ 0.301013  0.888299  1.03866  1.0867
+ 0.853248  1.44053   1.5909   1.63894
+ 0.687546  1.27483   1.4252   1.47324
+julia> v = f.(rand(100,1,1), rand(1,100,100));
+julia> @btime destruct(v);
+  7.138 ms (7 allocations: 22.89 MiB)
+julia> x, y, z = f.(rand(100,1,1), rand(1,100,100)) |> destruct;
+```
+"""
+# TODO: get rit of nloops, use eachindex instead.
+@generated function destruct(v::Array{T,N}) where {T <: Tuple, N}
+    TT = T.types
+    M = length(TT)
+    quote
+        @nexprs $M i -> out_i = similar(v,$TT[i])
+        @inbounds @nloops $N j v begin
+            @nexprs $M i -> ((@nref $N out_i j) = (@nref $N v j)[i])
+        end
+        return @ntuple $M out
+    end
+end
+
+@generated function destruct_cartesian(v::Array{T,N}) where {T <: Tuple, N}
+    TT = T.types
+    M = length(TT)
+    quote
+        @nexprs $M i -> out_i = similar(v,$TT[i])
+        @inbounds @nloops $N j v begin
+            @nexprs $M i -> ((@nref $N out_i j) = (@nref $N v j)[i])
+        end
+        return @ntuple $M out
+    end
+end
+
+
+end # module Destruct

src/scratch/comparative_timing.jl

@@ -0,0 +1,27 @@
+using Destruct
+using BenchmarkTools
+
+unpack_bc(w::Array{<:Tuple}) = Tuple((v->v[i]).(w) for i=1:length(w[1]))
+
+f(a, b) = a+b, a*b, a-b
+
+println("Homogeneous types")
+packed = f.(rand(2,2,1), rand(1,1,2))
+println("$(typeof(packed))")
+for sz in [10 50 100]
+    packed = f.(rand(sz,sz,1), rand(1,1,sz))
+    r = (@belapsed unpack_bc($packed))/(@belapsed destruct($packed))
+    println("$sz^3 : $r")
+end
+println("")
+
+g(a, b) = a+b+1im*(a-b), a*b, convert(Int, round(a-b))
+
+println("Heterogeneous types")
+packed = g.(rand(2,2,1), rand(1,1,2))
+println("$(typeof(packed))")
+for sz in [10 50 100]
+    packed = g.(rand(sz,sz,1), rand(1,1,sz))
+    r = (@belapsed unpack_bc($packed))/(@belapsed destruct($packed))
+    println("$sz^3 : $r")
+end

src/scratch/compare_tup_ntup.jl

@@ -0,0 +1,32 @@
+using Base.Cartesian
+using BenchmarkTools
+
+@generated function unpack_n(v::Array{T, N}) where {M, U, T <: NTuple{M, U}, N}
+    quote
+        @nexprs $M i -> out_i = similar(v,U)
+        @inbounds @nloops $N j v begin
+            @nexprs $M i -> ((@nref $N out_i j) = (@nref $N v j)[i])
+        end
+        return @ntuple $M out
+    end
+end
+
+@generated function unpack_t(v::Array{T,N}) where {T <: Tuple, N}
+    TT = T.types
+    M = length(TT)
+    quote
+        @nexprs $M i -> out_i = similar(v,$TT[i])
+        @inbounds @nloops $N j v begin
+            @nexprs $M i -> ((@nref $N out_i j) = (@nref $N v j)[i])
+        end
+        return @ntuple $M out
+    end
+end
+
+f(a, b) = a+b, a*b, a-b
+
+for sz in [10 50 100 200]
+    packed = f.(rand(sz,sz,1), rand(1,1,sz))
+    r = (@belapsed unpack_t($packed))/(@belapsed unpack_n($packed))
+    println("$sz^3 : $r")
+end

src/scratch/quick_timing.jl

@@ -0,0 +1,13 @@
+using Unpack
+using BenchmarkTools
+
+function unpack_bc(w::Array{<:Tuple})
+    Tuple((v->v[i]).(w) for i=1:length(w[1]))
+end
+
+g(a, b) = a+1im*b, a*b, convert(Int, round(a-b))
+packed = g.(rand(100,100,1), rand(1,1,100))
+println("Broadcast")
+@btime unpack_bc($packed)
+println("Cartesian")
+@btime unpack($packed)

src/scratch/smoke.jl

@@ -0,0 +1,16 @@
+using Unpack
+
+transform(x, y) = x+y, x*y, x-y
+
+x = 1.0:10000.0
+y = 1.0:10000.0
+x_grid = reshape(x, (length(x), 1))
+y_grid = reshape(y, (1, length(y)))
+#v = transform.(x, y)
+#typeof(v)  # Array{Tuple{Int64,Int64,Int64},1}
+
+#w = transform.(
+#    reshape(x, (length(x), 1)),
+#    reshape(y, (1, length(y)))
+#)
+#typeof(w)  # Array{Tuple{Int64,Int64,Int64},2}

src/scratch/timing.jl

@@ -0,0 +1,46 @@
+using Unpack
+using BenchmarkTools
+
+function unpack_bc(w::Array{<:Tuple})
+    Tuple((v->v[i]).(w) for i=1:length(w[1]))
+end
+
+# some transform
+f(a, b) = a+b, a*b, a-b
+# some transform with different return types
+g(a, b) = a+b+1im*(a-b), a*b, convert(Int, round(a-b))
+
+shape = (10,10)
+a_1 = rand(shape)
+a_2 = rand(shape)
+a_3 = rand(shape)
+unpacked = (a_1, a_2, a_3)
+packed = collect(zip(a_1, a_2, a_3))
+
+@assert unpack_bc(packed) == unpacked
+@assert unpack(packed) == unpacked
+
+println("Array of NTuple")
+println("---------------")
+for sz=[10 100 200]
+    packed = f.(rand(sz,sz,1), rand(1,1,sz))
+    shape = size(packed)
+    println("shape: $shape")
+    println("broadcast")
+    @btime unpack_bc($packed)
+    println("cartesian")
+    @btime unpack($packed)
+    println("")
+end
+println("Array of Tuple")
+println("--------------")
+for sz=[10 100 200]
+    packed = g.(rand(sz,sz,1), rand(1,1,sz))
+    shape = size(packed)
+    println("shape: $shape")
+    println("broadcast")
+    @btime unpack_bc($packed)
+    println("cartesian")
+    @btime unpack($packed)
+    println("")
+end

test/runtests.jl

@@ -1,5 +1,30 @@
 using Destruct
 using Base.Test
+import Base.rand
 
-# write your own tests here
-@test 1 == 2
+rand(rng::AbstractRNG, T::Type{String}) = randstring(rng)
+types = [Int32, Float64, Complex128, Bool]
+
+@testset "NTuple" begin
+for T=types, N=[1,2,3,4,5]
+    sz = fill(5, N)
+    a = rand(T, sz...); b = rand(T, sz...); z = collect(zip(a,b))
+    @assert size(z) == size(a) # make sure collect preserved array shape
+    @test (a,b) == destruct(z)
+end end
+
+@testset "Tuple" begin
+for T1=types, T2=types, N=[1,2,3]
+    sz = fill(3,N)
+    a = rand(T1, sz...); b = rand(T2, sz...); z = collect(zip(a,b))
+    @assert size(z) == size(a) # make sure collect preserved array shape
+    @test (a,b) == destruct(z)
+end end
+
+@testset "Tuple, Non square" begin
+for T1=types, T2=types, N=[2,3]
+    sz = 2:N+1
+    a = rand(T1, sz...); b = rand(T2, sz...); z = collect(zip(a,b))
+    @assert size(z) == size(a) # make sure collect preserved array shape
+    @test (a,b) == destruct(z)
+end end