# Intersection of ProductIterator

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}}`?

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

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)
``````