Just use BroadcastArray from LazyArrays?

[dlfivefifty]

Note you can accomplish the same thing via, eg,

f = (k,j) -> …
BroadcastArray(f, 1:n, (1:m)’)

This supports arbitrary axes including offset.

[timholy]

I knew something like this had to exist. That said, it looks like it's not a subtype of AbstractArray. That's going to be a compatibility problem for many uses.

[dlfivefifty]

Of course it’s a subtype of AbstractArray

[dlfivefifty]

Note also it supports parameters so avoids closure issues:

BroadcastVector((k,x) -> cos(k*x), 1:5, 0.1)

[timholy]

Oh, I misunderstood the meaning of the README (it was referring to Broadcasted, not BroadcastArray).

I'm noticing a performance issue, though. Modifying

GetindexArrays.jl/test/runtests.jl

Lines 74 to 107 in 45cae43

	function special_getindex(A::AbstractArray{T,N}, tforms, i::Vararg{Int,N}) where {T,N}
	# Get the spatial slice and transformation for this time point
	tax, ti = timeaxis(A), i[timedim(A)]
	Aspatial = view(A, tax(ti))
	tform = tforms[ti]
	# Compute the transformed center position
	si = ImageAxes.filter_space_axes(AxisArrays.axes(A), i)
	tformi = tform(SVector(si))
	ci = CartesianIndex(Tuple(round.(Int, tformi)))
	# Average over a spatial 3x3x... region, considering just those points that are in-bounds
	sinds = CartesianIndices(axes(Aspatial))
	s, n = float(zero(eltype(A))), 0
	@inbounds for j in max(first(sinds), ci-oneunit(ci)):min(last(sinds), ci+oneunit(ci))
	s += Aspatial[j] # cumulative sum
	n += 1 # number of in-bounds voxels
	end
	# square-root transform and return the result
	# Both `s` and `n` will be zero if there were no in-bounds points, in which case we return NaN (which is good)
	return sqrt(s/n)
	end
	# Because tforms is captured, construct everything inside a function
	function make_all(P, tforms)
	# forced specialization
	capture_getindex(i::Vararg{Int,N}) where {N} = special_getindex(P, tforms, i...)
	A = @inferred(GetindexArray(capture_getindex, axes(P)))
	return A
	end
	tforms = [LinearMap(RotMatrix(θ)) for θ = 0:π/16:2π]
	P = zeros(8, 13, length(tforms))
	P[3:6, 4:9, :] .= 1
	P = AxisArray(OffsetArray(P, -3:4, -6:6, 1:length(tforms)), :y, :x, :time)
	A = make_all(P, tforms)
	@test @inferred(A[0, 0, lastindex(tforms)]) === 1.0

to use BroadcastArray, which requires something like

using ImageFiltering   # for kernelfactors (adjoint only covers up to 2d, this is 3d so higher-dimensional axis reshaping is needed)
function make_all(P, tforms)
    # forced specialization
    capture_getindex(i::Vararg{Int,N}) where {N} = special_getindex(P, tforms, i...)
    A = @inferred(BroadcastArray(capture_getindex, kernelfactors(axes(P))...))
    return A
end

Then I get

julia> using BenchmarkTools

julia> @btime $A[i,j,k] setup=(i=rand(axes(A,1)); j=rand(axes(A,2)); k=rand(axes(A,3)))
  178.280 ns (14 allocations: 480 bytes)
0.0

whereas with this package it's

julia> @btime $A[i,j,k] setup=(i=rand(axes(A,1)); j=rand(axes(A,2)); k=rand(axes(A,3)))
  113.538 ns (11 allocations: 352 bytes)

Update: SimpleTraits has a performance issue; if I stop using it in ImageAxes and just handle the traits manually, then I get
about 48ms for BroadcastArray and 26ms for this package. So it's almost a factor of 2 difference. That's enough for me to care about, though if BroadcastArray can be fixed then I'd be happy to delete this repo.

Thoughts?

[dlfivefifty]

I didn't try your explicit example but note that getindex at a single index is going to be hit by overhead by checking for special broadcasting due to 1-length. Note it's reducing to getindex(::Base.Broadcasting.Broadcast, k::Integer, j::Integer) so any slowness is due to slowness in broadcasting. It might be possible to special case here where at compile time we can recognise the dimensions which are 1.

Note LazyArrays.jl does much better with slices than with single getindex calls. If you take a slice they are around the same time:

julia> n = 1000; A = GetindexArray((k,j) -> cos(k+j), (Base.OneTo(n), Base.OneTo(n)))

julia> @btime A[1:500,1:500];
  2.661 ms (2 allocations: 1.91 MiB)

julia> n = 1000; A = BroadcastArray((k,j) -> cos(k+j), Base.OneTo(n), Base.OneTo(n)')

julia> @btime A[1:500,1:500];
  2.652 ms (2 allocations: 1.91 MiB)

And if you want parameters it's a 4x performance improvement usng BroadcastArray:

julia> let x = 0.1
       n = 1000; A = GetindexArray((k,j) -> cos(x * (k+j)), (Base.OneTo(n), Base.OneTo(n)));
       end

julia> @btime A[1:500,1:500];
  14.769 ms (619807 allocations: 11.36 MiB)

julia> x = 0.1; n = 1000; A = BroadcastArray((k,j,x) -> cos(x*(k+j)), Base.OneTo(n), Base.OneTo(n)',x);

julia> @btime A[1:500,1:500];
  3.040 ms (2 allocations: 1.91 MiB)

[timholy]

Yep, though #1 does the same parameters trick.

This package was inspired by a demo where I'm making the point that scalar-indexing performance matters, so I don't want to sweep it under the rug. And it's possible to have no-overhead scalar indexing here.

[timholy]

Hah, we can of course do the same thing with MappedArrays, and there it seems there's no performance penalty. Construction is with mappedarray(f, CartesianIndices(axs)). Maybe we really don't need this package.

[dlfivefifty]

Your Extras seems like it's trying to recreate broadcasting, just without the 1-length special case...

Ideally in Julia v2.0 the following changes:

1. Scalars have special axes, say called Scalar (or OneToOne)
2. transpose/adjoints have axes(a', 1) isa Scalar
3. Broadcasting which propagates along other axes is only supported when the axes is either past the axes or is Scalar.

This would mean special cases can be done at compile-time with no performance penalty. It would also make clearer that a'*b returns a Number.

[timholy]

I totally agree that changes to broadcasting have wide-ranging benefits when there's an unambiguous improvement, and if some of those can't be done in Julia 1.x because they would be breaking, then we should plan on fixing what we can for Julia 2.0.

But for this particular effort, I'm actually not focused on broadcasting; it's more about supporting things like spatial deformations of the underlying grid. For example, image registration of a timeseries 2d or 3d movie, where the original array might be many TB. JuliaImages & JuliaArrays support a ton of lazy operations in composable fashion, but layering has its disadvantages sometimes, and so fusing them together into a single closure does have advantages.