Add functionality to use arbitrary singular values in ConcatDiskArrays (#256)

* Add test case for zero fill concatdiskarray

* Enable arbitrary single values in concatDiskArray

* Update src/cat.jl

Co-authored-by: Rafael Schouten <[email protected]>

* Update src/cat.jl

Co-authored-by: Rafael Schouten <[email protected]>

* Apply review

Co-authored-by: Rafael Schouten <[email protected]>

* Wrap fill values in MissingTile for ConcatDiskArray

This introduces a MissingTile type to wrap the values for ConcatDiskArray.
This allows to use vectors as elements in the fillvalue.

---------

Co-authored-by: Felix Cremer <[email protected]>
Co-authored-by: Rafael Schouten <[email protected]>

Author: 3 people

src/cat.jl

@@ -1,4 +1,3 @@
-
 """
     ConcatDiskArray <: AbstractDiskArray
     
@@ -31,16 +30,17 @@ function ConcatDiskArray(arrays::AbstractArray{Union{<:AbstractArray,Missing}})
     M = ndims(et)
     _ConcatDiskArray(arrays, T, Val(N), Val(M))
 end
+
 function infer_eltypes(arrays)
     foldl(arrays, init=(-1, Union{})) do (M, T), a
-        if ismissing(a)
-            (M, promote_type(Missing, T))
-        else
+        if a isa AbstractArray
             M == -1 || ndims(a) == M || throw(ArgumentError("All arrays to concatenate must have equal ndims"))
-            (ndims(a), promote_type(eltype(a), T))
+            M = ndims(a)
+        end
+            (M, promote_type(eltype(a), T))
         end
     end
-end
+
 function ConcatDiskArray(arrays::AbstractArray{<:AbstractArray})
     N = ndims(arrays)
     T = eltype(eltype(arrays))
@@ -78,6 +78,11 @@ function ConcatDiskArray(arrays1::AbstractArray, T, ::Val{D},::Val{ID}) where {D
     return ConcatDiskArray{T,D,typeof(arrays1),typeof(chunks),typeof(hc),ID}(arrays1, startinds, sizes, chunks, hc,Val(ID))
 end
 
+struct MissingTile{F}
+    fillvalue::F
+end
+Base.eltype(::Type{MissingTile{F}}) where F = F
+
 function extenddims(a::Tuple{Vararg{Any,N}}, b::Tuple{Vararg{Any,M}}, fillval) where {N,M} 
     length(a) > length(b) && error("Wrong")
     extenddims((a..., fillval), b, fillval)
@@ -90,7 +95,7 @@ function arraysize_and_startinds(arrays1)
     sizes = map(i -> zeros(Int, i), size(arrays1))
     for i in CartesianIndices(arrays1)
         ai = arrays1[i]
-        ismissing(ai) && continue
+        ai isa MissingTile && continue
         sizecur = extenddims(size(ai), size(arrays1), 1)
         foreach(sizecur, i.I, sizes) do si, ind, sizeall
             if sizeall[ind] == 0
@@ -123,10 +128,10 @@ function readblock!(a::ConcatDiskArray, aout, inds::AbstractUnitRange...)
     # Find affected blocks and indices in blocks
     _concat_diskarray_block_io(a, inds...) do outer_range, array_range, I
         vout = view(aout, outer_range...)
-        if ismissing(I)
-            vout .= missing
-        else
+        if I isa CartesianIndex
             readblock!(a.parents[I], vout, array_range...)
+        else 
+            vout .= (I.fillvalue,)
         end
     end
 end
@@ -170,8 +175,8 @@ function _concat_diskarray_block_io(f, a::ConcatDiskArray, inds...)
         #Shorten array range to shape of actual array
         array_range = ntuple(j -> array_range[j], ID)
         outer_range = fix_outerrangeshape(outer_range, array_range)
-        if ismissing(myar)
-            f(outer_range, array_range, missing)
+        if myar isa MissingTile
+            f(outer_range, array_range, myar)
         else
             f(outer_range, array_range, cI)
         end
@@ -189,13 +194,13 @@ function concat_chunksize(parents)
     newchunks = map(s -> Vector{Union{RegularChunks,IrregularChunks}}(undef, s), size(parents))
     for i in CartesianIndices(parents)
         array = parents[i]
-        ismissing(array) && continue
+        array isa MissingTile && continue
         chunks = eachchunk(array)
         foreach(chunks.chunks, i.I, newchunks) do c, ind, newc
             if !isassigned(newc, ind)
                 newc[ind] = c
             elseif c != newc[ind]
-                throw(ArgumentError("Chunk sizes don't forma grid"))
+                throw(ArgumentError("Chunk sizes don't form a grid"))
             end
         end
     end

src/subarray.jl

@@ -44,7 +44,7 @@ function eachchunk_view(::Chunked, vv)
 end
 eachchunk_view(::Unchunked, a) = estimate_chunksize(a)
 
-# Implementaion macro
+# Implementation macro
 
 macro implement_subarray(t)
     t = esc(t)

test/runtests.jl

@@ -492,22 +492,52 @@ end
         @test slic == Float64[1, 2, 3, 4, 1, 2, 3, 4]
     end
 
+    @testset "Concat DiskArray with fill zero tiles" begin
+        a = zeros(Int, 3, 4)
+        b = ones(Int, 2, 4)
+        c = fill(2, 3, 5)
+        d = fill(0, 2, 5)
+        aconc = DiskArrays.ConcatDiskArray(reshape([a, b, c, DiskArrays.MissingTile(0)], 2, 2))
+        abase = [a c; b d]
+        @test all(isequal.(aconc[:, :], abase))
+        @test all(isequal.(aconc[3:4, 4:6], abase[3:4, 4:6]))
+        ch = DiskArrays.eachchunk(aconc)
+        @test ch.chunks[1] == [1:3, 4:5]
+        @test ch.chunks[2] == [1:4, 5:9]
+        @test eltype(aconc) == Int
+
+        a = ones(100, 50)
+        b = [rem(i.I[3], 5) == 0 ? DiskArrays.MissingTile(0) : a for i in CartesianIndices((1, 1, 100))]
+        b[1] = DiskArrays.MissingTile(0)
+        a_conc = DiskArrays.ConcatDiskArray(b)
+        ch = eachchunk(a_conc)
+        @test ch.chunks[1] == [1:100]
+        @test ch.chunks[2] == [1:50]
+        @test ch.chunks[3] === DiskArrays.RegularChunks(1, 0, 100)
+
+        @test all(isequal.(a_conc[2, 2, 1:5], [0, 1.0, 1.0, 1.0, 0]))
+        @test all(isequal.(a_conc[end, end, 95:100], [0, 1.0, 1.0, 1.0, 1.0, 0]))
+
+    end
+
+
     @testset "Concat DiskArray with missing tiles" begin
         a = zeros(Int, 3, 4)
         b = ones(Int, 2, 4)
         c = fill(2, 3, 5)
         d = fill(missing, 2, 5)
-        aconc = DiskArrays.ConcatDiskArray(reshape([a, b, c, missing], 2, 2))
+        aconc = DiskArrays.ConcatDiskArray(reshape([a, b, c, DiskArrays.MissingTile(missing)], 2, 2))
         abase = [a c; b d]
         @test all(isequal.(aconc[:, :], abase))
         @test all(isequal.(aconc[3:4, 4:6], abase[3:4, 4:6]))
         ch = DiskArrays.eachchunk(aconc)
         @test ch.chunks[1] == [1:3, 4:5]
         @test ch.chunks[2] == [1:4, 5:9]
+        @test eltype(aconc) == Union{Int, Missing}
 
         a = ones(100, 50)
-        b = [rem(i.I[3], 5) == 0 ? missing : a for i in CartesianIndices((1, 1, 100))]
-        b[1] = missing
+        b = [rem(i.I[3], 5) == 0 ? DiskArrays.MissingTile(missing) : a for i in CartesianIndices((1, 1, 100))]
+        b[1] = DiskArrays.MissingTile(missing)
         a_conc = DiskArrays.ConcatDiskArray(b)
         ch = eachchunk(a_conc)
         @test ch.chunks[1] == [1:100]
@@ -518,6 +548,35 @@ end
         @test all(isequal.(a_conc[end, end, 95:100], [missing, 1.0, 1.0, 1.0, 1.0, missing]))
 
     end
+
+    @testset "Concat DiskArray with fill zero vector tiles" begin
+        a = fill([1,1], 3, 4)
+        b = fill([1,2], 2, 4)
+        c = fill([2,1], 3, 5)
+        d = fill([2,2], 2, 5)
+        aconc = DiskArrays.ConcatDiskArray(reshape([a, b, c, DiskArrays.MissingTile([2,2])], 2, 2))
+        abase = [a c; b d]
+        @test all(isequal.(aconc[:, :], abase))
+        @test all(isequal.(aconc[3:4, 4:6], abase[3:4, 4:6]))
+        ch = DiskArrays.eachchunk(aconc)
+        @test ch.chunks[1] == [1:3, 4:5]
+        @test ch.chunks[2] == [1:4, 5:9]
+        @test eltype(aconc) == Vector{Int}
+
+        a = fill([1,1], 100, 50)
+        b = [rem(i.I[3], 5) == 0 ? DiskArrays.MissingTile([0,0]) : a for i in CartesianIndices((1, 1, 100))]
+        b[1] = DiskArrays.MissingTile([0,0])
+        a_conc = DiskArrays.ConcatDiskArray(b)
+        ch = eachchunk(a_conc)
+        @test ch.chunks[1] == [1:100]
+        @test ch.chunks[2] == [1:50]
+        @test ch.chunks[3] === DiskArrays.RegularChunks(1, 0, 100)
+
+        @test all(isequal.(a_conc[2, 2, 1:5], [[0,0], [1,1],[1,1] , [1,1], [0,0]]))
+        @test all(isequal.(a_conc[end, end, 95:100], [[0,0], [1,1], [1,1], [1,1],[1,1], [0,0]]))
+
+    end
+
 end
 
 @testset "Broadcast with length 1 and 0 final dim" begin
@@ -929,8 +988,10 @@ struct TestArray{T,N} <: AbstractArray{T,N} end
     DiskArrays.@implement_array_methods TestArray
     DiskArrays.@implement_permutedims TestArray
     DiskArrays.@implement_subarray TestArray
-    DiskArrays.@implement_diskarray TestArray
     @test DiskArrays.isdisk(TestArray) == true
+    DiskArrays.@implement_diskarray TestArray2
+    @test DiskArrays.isdisk(TestArray2) == true
+
 end
 
 # issue #123