cannot @cast on views of other @cast results

[prittjam]

This does not work,

using TensorCast, StaticArrays
world_basis_to_P(R,c; K=SMatrix{3,3}(one(eltype(R))I)) = K*[R' -R'*c]
R = rand(3,3,100)
X₁ = rand(3,100)
@cast P[i,j,k] := world_basis_to_P(R[:,:,k],X₁[:,k])[i,j]
HH = view(P,:,[1,2,4],:)
@cast invH[i,j,k] := inv(HH[:,:,k])[i,j]

but this does

using TensorCast, StaticArrays
world_basis_to_P(R,c; K=SMatrix{3,3}(one(eltype(R))I)) = K*[R' -R'*c]
R = rand(3,3,100)
X₁ = rand(3,100)
@cast P[i,j,k] := world_basis_to_P(R[:,:,k],X₁[:,k])[i,j]
HH = P,[:,[1,2,4],:]
@cast invH[i,j,k] := inv(HH[:,:,k])[i,j]

[mcabbott]

The problem is that inv doesn't accept sufficiently exotic AbstractMatrix types:

julia> HH = view(P,:,[1,2,4],:)
3×3×100 view(lazystack(::Vector{Matrix{Float64}}), :, [1, 2, 4], :) ...

julia> inv(HH[:,:,1])  # inv(::Matrix)
3×3 Matrix{Float64}:
  6.7336   -18.7875   21.6366
 -4.88027    1.76205   0.148484
  2.17953  -14.5834   16.9896

julia> inv(@view HH[:,:,1])  # same error as @cast
ERROR: MethodError: no method matching factorize(::SubArray{Float64, 2, LazyStack.Stacked{...

julia> @cast invH[i,j,k] := inv(HH[:,:,k])[i,j] lazy=false
3×3×100 Array{Float64, 3}:

@cast is quite lazy by default, which saves allocations, and scalar broadcasting tends to work fine with this. But after a few operations you can get a very complicated weird type, and this may be slow to deal with, or as here, fail completely. lazy=false should opt out of all of this.

[prittjam]

Thanks for the response.

What do you mean by "slow to deal with?" I assumed that it would not affect the runtime to resolve the type (after the first run).

[mcabbott]

Yes complicated types alone aren't a problem. But complicated array types are, partly because accessing their elements involves more steps (like translating indices for views, and following different pointers for slices) and partly because they prevent the use of optimised routines which expect a simple Array:

julia> @btime sum(x)  setup=(x=randn(30,30));
  min 247.120 ns, mean 249.191 ns (0 allocations)

julia> @btime x*x  setup=(x=randn(30,30));
  min 1.483 μs, mean 3.319 μs (1 allocation, 7.19 KiB)

julia> using LazyStack

julia> @btime sum(x)  setup=(x=view(lazystack([randn(30) for _ in 1:30]), 30:-1:1, :));
  min 4.411 μs, mean 4.446 μs (0 allocations)  # complicated to access each element

julia> @btime x*x  setup=(x=lazystack([randn(30) for _ in 1:30]));  # won't use BLAS
  min 13.375 μs, mean 16.052 μs (3 allocations, 27.69 KiB)

I don't think there's a great general rule for when such costs outweigh the costs of making an array. This package's defaults are quite lazy, picturing things like this:

julia> @btime abs2.(Compat.stack(x))  setup=(x=[randn(30) for _ in 1:30]);
  min 900.615 ns, mean 3.502 μs (2 allocations, 14.38 KiB)

julia> @btime abs2.(lazystack(x))  setup=(x=[randn(30) for _ in 1:30]);
  min 549.690 ns, mean 1.710 μs (1 allocation, 7.19 KiB)