add StepRange support for CartesianIndices

[johnnychen94]

This PR intended to only add StepRange support for CartesianIndices, but more things are added here as I was working on it.

Three major changes are involved here:

StepRange CartesianIndices

The biggest change here is to make it possible to pass a step to CartesianIndices so that more flexible loop behavior can be done

R = CartesianIndices(A)
Δ = CartesianIndex(ntuple(_->2, ndims(A)))

for I in first(R):Δ:last(R)
    ...
end

Another new thing introduced by this PR is that step(::CartesianIndices) now returns a CartesianIndex, this enables us to treat CartesianIndices as a multidimensional range type.

Reduce CartesianIndices overhead

By reworking the iterating part, the overhead issue #9080 is accidentally fixed.

using Test
using BenchmarkTools

function sumcart_iter(A)
    s = zero(eltype(A))
    for I in CartesianIndices(A)
        @inbounds s += A[I]
    end
    s
end
function sumcart_manual(A::AbstractMatrix)
    s = zero(eltype(A))
    ax1, ax2 = axes(A)
    for j in ax2, i in ax1
        @inbounds s += A[i, j]
    end
    s
end

julia> A = rand(4, 4);

julia> @btime sumcart_iter($A)
  6.757 ns (0 allocations: 0 bytes) # this PR
  15.315 ns (0 allocations: 0 bytes) # master

julia> @btime sumcart_manual($A) # baseline
  8.013 ns (0 allocations: 0 bytes)

julia> A = rand(256,20);

julia> @btime sumcart_iter($A);
  4.543 μs (0 allocations: 0 bytes) # this PR
  5.017 μs (0 allocations: 0 bytes) # master

julia> @btime sumcart_manual($A);
  4.544 μs (0 allocations: 0 bytes) # baseline

So 🎉 looping over CartesianIndices is marginally faster than a manual loop.

There's still some overhead for StepRange:

function sumcart_iter(A)
    s = zero(eltype(A))
    R = CartesianIndices(A)
    Δ = CartesianIndex(ntuple(_->2, ndims(R)))
    for I in first(R):Δ:last(R)
        @inbounds s += A[I]
    end
    s
end

function sumcart_manual(A::AbstractMatrix)
    s = zero(eltype(A))
    ax1, ax2 = axes(A)
    for j in first(ax2):2:last(ax2), i in first(ax1):2:last(ax1)
        @inbounds s += A[i, j]
    end
    s
end

julia> A = rand(8, 8);

julia> @btime sumcart_iter($A)
  36.418 ns (0 allocations: 0 bytes)

julia> @btime sumcart_manual($A)
  32.492 ns (0 allocations: 0 bytes)

but that's more difficult to deal with because (IIUC) step can be either positive or negative. I'm not very familiar with this kind of low-level optimization...

tests and bug fixes

A lot of tests are spread in other files, I tried to collect most of the test patterns into test/cartesian.jl so it is easier to first check if include("test/cartesian.jl") works before triggering the whole heavy make test workflow.

A bug is found and fixed (see also #37928) :

# master
julia> CartesianIndices((2, 1:2))
1×2 CartesianIndices{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}}:
 CartesianIndex(2, 1)  CartesianIndex(2, 2)

# this PR
julia> CartesianIndices((2, 1:2))
2×2 CartesianIndices{2, Tuple{Base.OneTo{Int64}, UnitRange{Int64}}}:
 CartesianIndex(1, 1)  CartesianIndex(1, 2)
 CartesianIndex(2, 1)  CartesianIndex(2, 2)

# master/this PR
julia> CartesianIndices((2, 2))
2×2 CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}:
 CartesianIndex(1, 1)  CartesianIndex(1, 2)
 CartesianIndex(2, 1)  CartesianIndex(2, 2)

New tests are added to test common operations and properties. Old tests are still kept so as to not break anything unexpectedly.

---

There is still something I want to tackle, but I'm a bit running out of my spare time to keep working on this, which I'll leave them as future work:

• StepRange support to LinearIndices (add StepRange support for CartesianIndices #37829 (comment))
• tweak the StepRange CartesianIndices overhead

cc: @timholy @mbauman

closes #37577
closes #9080

[timholy]

I think this is a good change. But the biggest worry about this is that a lot of code written for CartesianIndices has formerly assumed that the step size would always be 1, and that will no longer be true. Sort of like the transition to indexing starting at something other than 1, this could trigger a bigger shakeup across the ecosystem. It's probably the right thing to do, but be aware that this change may have a "tail" that could last for years. (It's probably a lot less disruptive than the 1-indexing change, though.)

One thing we can probably do to ease the transition is define

const UnitCartesianIndices{N,R<:NTuple{N,AbstractUnitRange{Int}}} = CartesianIndices{N,R}

and encourage people to use it when they want to enforce a step size of 1. (Not sure whether to export this, though, probably not.)

[codecov-commenter]

Codecov Report

Merging #37829 into master will decrease coverage by 0.06%.
The diff coverage is 91.66%.

@@            Coverage Diff             @@
##           master   #37829      +/-   ##
==========================================
- Coverage   87.64%   87.58%   -0.07%     
==========================================
  Files         368      368              
  Lines       73053    73038      -15     
==========================================
- Hits        64030    63968      -62     
- Misses       9023     9070      +47     

Impacted Files	Coverage Δ
base/Base.jl	46.66% <0.00%> (ø)
stdlib/Distributed/src/cluster.jl	69.28% <ø> (ø)
stdlib/Test/src/Test.jl	86.13% <ø> (-5.11%)	⬇️
base/timing.jl	74.00% <75.00%> (+1.36%)	⬆️
base/binaryplatforms.jl	83.38% <100.00%> (+0.05%)	⬆️
base/errorshow.jl	88.26% <100.00%> (ø)
base/initdefs.jl	57.81% <100.00%> (ø)
base/methodshow.jl	82.86% <100.00%> (-0.07%)	⬇️
base/process.jl	91.27% <100.00%> (ø)
base/sort.jl	98.31% <100.00%> (+0.48%)	⬆️
... and 19 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 8a9666a...133f07f. Read the comment docs.

[johnnychen94]

Where should the new test goes, test/cartesian.jl?

[timholy]

A very impressive PR! Solving #9080 alone is amazing, but doing so while cleaning up implementations, adding more flexibility, testing more thoroughly is a tour de force.

Subject to fixing the test errors this seems basically ready from what I can tell, pending how you want to resolve the LinearIndices situation. To be clear, I think that could be a separate PR, if you want to go with the errorring version here first and get it out in the wild.

[timholy]

Actually, thinking about this more, maybe we can handle this? After all, the definition of LinearIndices is

struct LinearIndices{N,R<:NTuple{N,AbstractUnitRange{Int}}} <: AbstractArray{Int,N}
    indices::R
end

and so some of the same tricks you're using here could be used there too. You might have to reset the integer at the end of each "column," perhaps by maintaining a hybrid linear/cartesian iterator:

julia> inds = LinearIndices((1:3, 1:4))
3×4 LinearIndices{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}}:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> inds[1:2:3,:]   # can't just increase by a constant
2×4 Matrix{Int64}:
 1  4  7  10
 3  6  9  12

Only, of course, when it's a StepRange.

Keep in mind that converting Cartesian to Linear is pretty fast (multiply and add); it's the opposite transition, Linear to Cartesian, that requires div and hence is dirt-slow. (That's why Cartesian is the fallback.)

[johnnychen94]

Does step(R::CartesianIndices) = CartesianIndex(step.(R.indices)) be crazy to you? With this PR done, I feel it might be pretty useful to think CartesianIndices as a multidimensional range now.

[timholy]

Let's find out if generalizing this type causes package problems: @nanosoldier runtests(ALL, vs = ":master")

[nanosoldier]

Your package evaluation job has completed - possible new issues were detected. A full report can be found here. cc @maleadt

[timholy]

Interesting failure in https://github.com/JuliaCI/NanosoldierReports/blob/4086959a886ed94fbfd32417489e909b652a9952/pkgeval/by_hash/3386b42_vs_3559b72/logs/Images/1.6.0-DEV-84648ad804.log

[johnnychen94]

I think that's due to https://github.com/JuliaArrays/TiledIteration.jl/blob/master/src/TiledIteration.jl#L102 and is fixable.

[timholy]

Good. That seems like the only real one I've seen. I think we can merge this.

[timholy]

Thank you @johnnychen94! 🎆

[johnnychen94]

My very first newsworthy PR 🎉