Intersection of ProductIterator

Internals & Design
#1

I noticed that intersect on Iterators.ProductIterator{<:NTuple{N,UnitRange}} generates all the values and stores them in a Vector. I was expecting it to return an Iterators.ProductIterator{<:NTuple{N,UnitRange}} since intersect on UnitRange returns a UnitRange.

For instance, the following function

foo(a::Iterators.ProductIterator{<:NTuple{N,UnitRange}}, b::Iterators.ProductIterator{<:NTuple{N,UnitRange}}) where {N} =
    Iterators.product(map(intersect, a.iterators, b.iterators)...)

satisfies intersect(a, b) == vec( collect( foo(a, b) ) ).

I must admit that my usage of Iterators.ProductIterator is quite limited so perhaps I am overlooking something here?
If not, why not having intersect of Iterators.ProductIterator{<:NTuple{N,UnitRange}} return an Iterators.ProductIterator{<:NTuple{N,UnitRange}}?

#2

Maybe nobody wanted this specific thing before? It would not be hard to write the method, and could apply to intersect(::ProductIterator, ::ProductIterator) not just ranges I guess. However, in the case of unit ranges, this is also covered by CartesianIndices:

julia> a = Iterators.product(1:10, 1:10);

julia> b = Iterators.product(2:3, 9:99)
Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}, UnitRange{Int64}}}((2:3, 9:99))

julia> intersect(a,b)
4-element Vector{Tuple{Int64, Int64}}:
 (2, 9)
 (3, 9)
 (2, 10)
 (3, 10)

julia> @which intersect(aa, bb)
intersect(itr, itrs...) in Base at array.jl:2472

julia> ca = CartesianIndices((1:10, 1:10));

julia> cb = CartesianIndices((2:3, 9:99));

julia> intersect(ca, cb)
2×2 CartesianIndices{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}}:
 CartesianIndex(2, 9)  CartesianIndex(2, 10)
 CartesianIndex(3, 9)  CartesianIndex(3, 10)
1 Like
#3

Probably, though I am not sure what kind of Iterators.ProductIterator people use in general. For some general a::Iterators.ProductIterator and b::Iterators.ProductIterator, it may not be easy to represent their intersection as a Iterators.ProductIterator?

At least for Iterators.ProductIterator{<:NTuple{N,AbstractRange}} where {N} it should be safe that the intersection can be expressed as a Iterators.ProductIterator{<:NTuple{N,AbstractRange}} where {N}.

Indeed, the behaviour of CartesianIndices could also be changed accordingly.

FWIW, here’s a small benchmark

julia> foo(a::Iterators.ProductIterator{<:NTuple{N,AbstractRange}}, b::Iterators.ProductIterator{<:NTuple{N,AbstractRange}}) where {N} =
           Iterators.product(map(intersect, a.iterators, b.iterators)...)
foo (generic function with 1 method)

julia> a = Iterators.product(1:1_000, 1:1_000);

julia> b = Iterators.product(200:1_000, 800:1_000);

julia> ca = CartesianIndices((1:1_000, 1:1_000));

julia> cb = CartesianIndices((200:1_000, 800:1_000));

julia> intersect(a, b) == vec( collect(foo(a, b)) )
true

julia> @btime intersect($a, $b);
  245.417 ms (46 allocations: 60.25 MiB)

julia> @btime foo(Ref($a)[], Ref($b)[]);
  3.590 ns (0 allocations: 0 bytes)

julia> @btime intersect($ca, $cb);
  260.656 ms (46 allocations: 60.25 MiB)