Specialised `issubset` for `CartesianIndices`

[andreasKroepelin]

It seems like there is currently no specialised implementation of issubset for CartesianIndices, while there is one for AbstractUnitRange. This means that what could be a constant-time task (check if lower/upper bounds of potential subset are greater/less than lower/upper bounds of potential superset) instead checks every value individually and even allocates.

I think something like this could solve this issue quite easily:

function issubset(l::CartesianIndices{N, <: NTuple{N, AbstractUnitRange{Int}}}, r::CartesianIndices{N, <: NTuple{N, AbstractUnitRange{Int}}}) where N
    all(issubset.(l.indices, r.indices))
end

[Seelengrab]

Are you sure the current method allocates? It doesn't for me and is quite fast (though I agree that it could check the bounds for a speedup in very high dimensions/large spans).

julia> a                                                                                                          
3×5 CartesianIndices{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}}:                                               
 CartesianIndex(1, 10)  CartesianIndex(1, 11)  CartesianIndex(1, 12)  CartesianIndex(1, 13)  CartesianIndex(1, 14)
 CartesianIndex(2, 10)  CartesianIndex(2, 11)  CartesianIndex(2, 12)  CartesianIndex(2, 13)  CartesianIndex(2, 14)
 CartesianIndex(3, 10)  CartesianIndex(3, 11)  CartesianIndex(3, 12)  CartesianIndex(3, 13)  CartesianIndex(3, 14)
                                                                                                                  
julia> b                                                                                                          
2×4 CartesianIndices{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}}:                                               
 CartesianIndex(2, 10)  CartesianIndex(2, 11)  CartesianIndex(2, 12)  CartesianIndex(2, 13)                       
 CartesianIndex(3, 10)  CartesianIndex(3, 11)  CartesianIndex(3, 12)  CartesianIndex(3, 13)                       
                                                                                                                  
julia> @btime issubset($b,$a)                                                                                     
  22.166 ns (0 allocations: 0 bytes)                                                                              
true                                                                                                              

julia> versioninfo()                           
Julia Version 1.7.0-DEV.1072                   
Commit 364660c5fd (2021-05-07 06:23 UTC)       
Platform Info:                                 
  OS: Linux (x86_64-linux-gnu)                 
  CPU: Intel(R) Core(TM) i7-6600U CPU @ 2.60GHz
  WORD_SIZE: 64                                
  LIBM: libopenlibm                            
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)       
Environment:                                   
  JULIA_PKG_SERVER =                           
  JULIA_NUM_THREADS = 4                        

[andreasKroepelin]

It does so for "large" CartesianIndices. The code apparently uses an heuristic that defines "large" as > 70 elements.

You can try that with

a = CartesianIndices((1:1,))
b = CartesianIndices((1:71,))
@btime issubset($a, $b)

[Seelengrab]

Ah I see! I missed that threshold.

I couldn't get your fix to allocate (I was worried because of the broadcast) even with very large CartesianIndices (200 dimensions) and the old version almost immediately complained due to OutOfMemory, so looks very good to me! Mind opening a PR?

Check for benchmarks

julia> using BenchmarkTools                                                  
                                                                             
julia> @eval Base function issubset(l::CartesianIndices, r::CartesianIndices)
           all(issubset.(l.indices, r.indices))                              
       end                                                                   
issubset (generic function with 5 methods)                                   
                                                                             
julia> a = CartesianIndices((1:1,));                                         
                                                                             
julia> b = CartesianIndices((Iterators.repeated(1:71, 100)...,));            
                                                                             
julia> @benchmark issubset($a, $b)                                           
BechmarkTools.Trial: 10000 samples with 878 evaluations.                     
 Range (min … max):  131.891 ns …   3.866 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     142.938 ns               ┊ GC (median):    0.00%        
 Time  (mean ± σ):   166.098 ns ± 111.798 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%
                                                                             
  █▆▅▅▄▄▄▃▃▃▃▂▂▂▂▂▁▁                                            ▂            
  ██████████████████████▇█▇▇▆▇▇▆▅▆▆▆▆▅▄▄▅▆▅▁▃▄▄▆▄▅▅▃▃▄▄▃▄▁▄▄▄▄▄ █            
  132 ns        Histogram: log(frequency) by time        498 ns <            
                                                                             
 Memory estimate: 0 bytes, allocs estimate: 0.                               
                                                                             
julia> b = CartesianIndices((Iterators.repeated(1:71, 200)...,));            
                                                                             
julia> @benchmark issubset($a, $b)                                           
BechmarkTools.Trial: 10000 samples with 346 evaluations.                     
 Range (min … max):  261.272 ns …  16.332 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     278.902 ns               ┊ GC (median):    0.00%        
 Time  (mean ± σ):   333.423 ns ± 322.689 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%
                                                                             
  ██▄▅▄▄▄▃▃▃▄▃▃▂▂▂▂▁ ▁ ▁                                        ▂            
  █████████████████████████▇▇▆▆▆▆▅▆▆▄▅▅▅▅▅▄▆▃▅▄▃▄▄▃▄▄▄▃▄▅▄▁▃▄▄▄ █            
  261 ns        Histogram: log(frequency) by time        959 ns <            
                                                                             
 Memory estimate: 0 bytes, allocs estimate: 0.                               
                                                                             
julia> @benchmark issubset($b, $b)                                           
BechmarkTools.Trial: 10000 samples with 210 evaluations.                     
 Range (min … max):  363.810 ns …  44.885 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     388.571 ns               ┊ GC (median):    0.00%        
 Time  (mean ± σ):   481.775 ns ± 601.198 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%
                                                                             
  █▇▅▄▄▄▃▃▃▃▃▂▂▃▂▁▁▁▁                                           ▂            
  ████████████████████████▇▇▇▆▆▆▆▅▆▅▄▄▅▄▅▄▅▄▃▁▅▄▅▅▄▄▄▄▁▃▃▄▃▅▄▁▄ █            
  364 ns        Histogram: log(frequency) by time       1.76 μs <            
                                                                             
 Memory estimate: 0 bytes, allocs estimate: 0.                               

[johnnychen94]

A good proposal to me. Since the start:step:range feature for CartesianIndices in Julia 1.6, CartesianIndices becomes more like an N-dimensional range type.

Optionally we could short-circult the check:

all(t->issubset(t[1], t[2]), zip(l.indices, r.indices)) # short-circult

[Seelengrab]

Nice idea, that gives another factor ~2, even when it can't be short circuited!

Benchmark

julia> @eval Base function issubset(l::CartesianIndices, r::CartesianIndices)
           all(issubset.(l.indices, r.indices))                              
       end                                                                   
issubset (generic function with 5 methods)                                   
                                                                             
julia> @benchmark issubset($b, $b)                                           
BechmarkTools.Trial: 10000 samples with 202 evaluations.                     
 Range (min … max):  395.545 ns …  20.680 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     422.277 ns               ┊ GC (median):    0.00%        
 Time  (mean ± σ):   517.444 ns ± 445.408 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%
                                                                             
  █▇▅▄▄▃▃▃▃▂▂▂▂▂▂▁▁▁                                            ▂            
  ███████████████████████▇▇▇▆▆▇▆▆▅▅▅▅▅▄▅▄▄▅▄▄▄▃▄▃▃▃▁▄▁▃▃▃▅▄▃▅▄▁ █            
  396 ns        Histogram: log(frequency) by time       1.95 μs <            
                                                                             
 Memory estimate: 0 bytes, allocs estimate: 0.        

julia> @eval Base function issubset(l::CartesianIndices, r::CartesianIndices)
           all(t->issubset(t[1], t[2]), zip(l.indices, r.indices))           
       end                                                                   
issubset (generic function with 5 methods)                                   
                                                                             
julia> @benchmark issubset($b, $b)                                           
BechmarkTools.Trial: 10000 samples with 506 evaluations.                     
 Range (min … max):  219.368 ns …   8.683 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     233.399 ns               ┊ GC (median):    0.00%        
 Time  (mean ± σ):   280.988 ns ± 224.409 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%
                                                                             
  █▇▅▄▄▄▃▃▃▃▂▂▂▂▂▂▁▁▁ ▁                                         ▂            
  ███████████████████████▇▇▇▇▆▆▆▆▆▅▅▆▅▅▅▅▅▅▄▅▅▁▄▅▄▅▄▄▁▃▄▃▄▄▄▃▄▁ █            
  219 ns        Histogram: log(frequency) by time        908 ns <            
                                                                             
 Memory estimate: 0 bytes, allocs estimate: 0.                          

[devmotion]

It seems this was fixed by #56282