Introduce types, e.g. DivDiffData, to get around complicated broadcast machinery

[DanielVandH]

@dlfivefifty Is this roughly what you had in mind for trying to get around some of the broadcast machinery? I'm simply caching the broadcast results in a CachedBroadcastVector with cache_filldata! implemented as

function LazyArrays.cache_filldata!(K::CachedBroadcastVector{T, F, A, B}, inds) where {T, F, A, B}
    @show inds, K.datasize
    f, a, b = K.f, K.A, K.B
    resizedata!(a, maximum(inds)) 
    resizedata!(b, maximum(inds))
    @inbounds for k in inds
        K.data[k] = f(a[k], b[k])
    end
end

It's currently only implemented in divdiff to get an idea of what it looks like:

function divdiff(Q::SemiclassicalJacobi, P::SemiclassicalJacobi)
    A,B,_ = recurrencecoefficients(P)
    α,β,_ = recurrencecoefficients(Q)

    d = AccumulateAbstractVector(*, CBV(/, A, Vcat(1,α)))
    v1 = AccumulateAbstractVector(+, CBV(/, B, A))
    v2 = MulAddAccumulate(CBV(/, Vcat(0,0,α[2:∞]), α), CBV(/, Vcat(0,CBV(/, β,  α)),  α))
    v3 = AccumulateAbstractVector(*, Vcat(A[1]A[2], CBV(/, A[3:∞], α)))
    op1 = CBV(*, 1:∞, d)
    op2 = CBV(/, B[2:end], A[2:end])
    op3 = CBV(+, v1, op2)
    op4 = CBV(*, 1:∞, op3)
    op5 = CBV(/, β, α)
    op6 = CBV(*, α, v2)
    op7 = CBV(+, op5, op6)
    op8 = CBV(*, 2:∞, op7)
    op9 = CBV(-, op4, op8)
    op10 = CBV(*, op9, v3) 
    return _BandedMatrix(Vcat(op1', op10'), ∞, 2,-1)'
    # return _BandedMatrix(Vcat(((1:∞) .* d)', (((1:∞) .* (v1 .+ B[2:end]./A[2:end]) .- (2:∞) .* (α .* v2 .+ β ./ α)) .* v3)'), ∞, 2,-1)'
end

which is a bit annoying as it requires a bunch of individual CBV calls. I did initially start with a CachedFusedBroadcastVector but it's somewhat unreadable as I need to manually build the Tuple of functions and vectors. I also tried specific structs for each type of operation, e.g. RangeTimesAccumulate, but getting this type of struct correct seems nicer.

There's an issue with this implementation though that I'm not sure I understand. It keeps trying to indefinitely resize the vector in some cases. For example, putting a @show inds, K.datasize inside LazyArrays.cache_filldata!(K::CachedBroadcastVector{T, F, A, B}, inds) where {T, F, A, B} and differentiating a weighted basis,

julia> diff(Weighted(SemiclassicalJacobi(2.0, 2, 2, 2)))
...
[truncated]
...
(inds, K.datasize) = (454:454, (453,))
(inds, K.datasize) = (455:455, (454,))
(inds, K.datasize) = ERROR: InterruptException:
Stacktrace:
  [1] try_yieldto(undo::typeof(Base.ensure_rescheduled))
    @ Base .\task.jl:1128
  [2] wait()
    @ Base .\task.jl:1200
  [3] uv_write(s::Base.TTY, p::Ptr{UInt8}, n::UInt64)
    @ Base .\stream.jl:1081
  [4] unsafe_write(s::Base.TTY, p::Ptr{UInt8}, n::UInt64)
    @ Base .\stream.jl:1154
  [5] write
    @ .\strings\io.jl:246 [inlined]
  [6] print
    @ .\strings\io.jl:248 [inlined]
  [7] print(::Base.TTY, ::String, ::String, ::String)
    @ Base .\strings\io.jl:46
  [8] println(x1::String, x2::String)
    @ Base .\coreio.jl:6
  [9] cache_filldata!(K::SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{…}, inds::UnitRange{…})
    @ SemiclassicalOrthogonalPolynomials c:\Users\djv23\.julia\dev\SemiclassicalOrthogonalPolynomials.jl\src\derivatives.jl:52
 [10] _vec_resizedata!(B::SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{…}, n::Int64)
    @ LazyArrays C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\src\cache.jl:246
 [11] resizedata!
    @ C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\src\cache.jl:254 [inlined]
 [12] resizedata!
    @ C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\src\cache.jl:222 [inlined]
 [13] getindex
    @ C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\src\cache.jl:79 [inlined]
 [14] getindex
    @ C:\Users\djv23\.julia\juliaup\julia-1.12.1+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\adjtrans.jl:342 [inlined]
 [15] copyto_unaliased!(deststyle::IndexLinear, dest::LinearAlgebra.Adjoint{…}, srcstyle::IndexLinear, src::LinearAlgebra.Adjoint{…})
    @ Base .\abstractarray.jl:1090
 [16] copyto!(dest::LinearAlgebra.Adjoint{…}, src::LinearAlgebra.Adjoint{…})
    @ Base .\abstractarray.jl:1070
 [17] copyto!_layout(::ArrayLayouts.DualLayout{…}, ::LazyArrays.LazyLayout, dest::LinearAlgebra.Adjoint{…}, src::LinearAlgebra.Adjoint{…})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\ArrayLayouts.jl:266
 [18] copyto!_layout(dlay::ArrayLayouts.DualLayout{…}, ::ArrayLayouts.DualLayout{…}, dest::LinearAlgebra.Adjoint{…}, src::LinearAlgebra.Adjoint{…})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\memorylayout.jl:310
 [19] copyto!_layout(dest::LinearAlgebra.Adjoint{…}, src::LinearAlgebra.Adjoint{…})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\ArrayLayouts.jl:272
 [20] copyto!(dest::LinearAlgebra.Adjoint{…}, src::LinearAlgebra.Adjoint{…})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\ArrayLayouts.jl:285
 [21] copymutable
    @ .\abstractarray.jl:1201 [inlined]
 [22] copy
    @ .\abstractarray.jl:1144 [inlined]
 [23] map
    @ .\tuple.jl:359 [inlined]
 [24] copy(f::ApplyArray{Float64, 2, typeof(vcat), Tuple{LinearAlgebra.Adjoint{…}, LinearAlgebra.Adjoint{…}}})
    @ LazyArrays C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\src\lazyconcat.jl:97
 [25] copy(B::BandedMatrices.BandedMatrix{Float64, ApplyArray{…}, InfiniteArrays.OneToInf{…}})
    @ BandedMatrices C:\Users\djv23\.julia\packages\BandedMatrices\Vz6Ql\src\banded\BandedMatrix.jl:167
 [26] copy(A::LinearAlgebra.Adjoint{Float64, BandedMatrices.BandedMatrix{Float64, ApplyArray{…}, InfiniteArrays.OneToInf{…}}})
    @ InfiniteArraysBandedMatricesExt C:\Users\djv23\.julia\packages\InfiniteArrays\LdDTO\ext\InfiniteArraysBandedMatricesExt.jl:416
 [27] _copy_oftype(A::LinearAlgebra.Adjoint{Float64, BandedMatrices.BandedMatrix{…}}, ::Type{Float64})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\ArrayLayouts.jl:70
 [28] copy(M::Mul{…})
    @ LazyArraysBandedMatricesExt C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\ext\LazyArraysBandedMatricesExt.jl:571
 [29] materialize
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\mul.jl:137 [inlined]
 [30] mul
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\mul.jl:138 [inlined]
 [31] *(A::LinearAlgebra.Diagonal{…}, B::LinearAlgebra.Adjoint{…})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\mul.jl:303
 [32] _copy_ldiv_mul
    @ C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\src\linalg\inv.jl:115 [inlined]
 [33] copy
    @ C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\src\linalg\inv.jl:119 [inlined]
 [34] copy
    @ C:\Users\djv23\.julia\packages\ContinuumArrays\TDola\src\bases\bases.jl:368 [inlined]
 [35] copy(L::Ldiv{…})
    @ SemiclassicalOrthogonalPolynomials c:\Users\djv23\.julia\dev\SemiclassicalOrthogonalPolynomials.jl\src\SemiclassicalOrthogonalPolynomials.jl:416
 [36] materialize
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\ldiv.jl:22 [inlined]
 [37] ldiv
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\ldiv.jl:98 [inlined]
 [38] \(A::SemiclassicalJacobi{Float64}, B::QuasiArrays.ApplyQuasiMatrix{Float64, typeof(*), Tuple{…}})
    @ QuasiArrays C:\Users\djv23\.julia\packages\QuasiArrays\wlEWh\src\matmul.jl:34
 [39] divdiff(wP::Weighted{Float64, SemiclassicalJacobi{Float64}}, wQ::Weighted{Float64, SemiclassicalJacobi{Float64}})
    @ SemiclassicalOrthogonalPolynomials c:\Users\djv23\.julia\dev\SemiclassicalOrthogonalPolynomials.jl\src\derivatives.jl:110
 [40] diff(wQ::Weighted{Float64, SemiclassicalJacobi{Float64}}; dims::Int64)
    @ SemiclassicalOrthogonalPolynomials c:\Users\djv23\.julia\dev\SemiclassicalOrthogonalPolynomials.jl\src\derivatives.jl:115
 [41] diff(wQ::Weighted{Float64, SemiclassicalJacobi{Float64}})
    @ SemiclassicalOrthogonalPolynomials c:\Users\djv23\.julia\dev\SemiclassicalOrthogonalPolynomials.jl\src\derivatives.jl:113
 [42] top-level scope
    @ REPL[19]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> diff(Weighted(SemiclassicalJacobi(2.0, 2, 2, 2)))

Any idea what's going on with this implementation that it keeps wanting to resize? I thought _vec_resizedata!(B::AbstractVector, n) should get around this automatically. Note that the indexing does work sometimes though, e.g.

julia> diff(SemiclassicalJacobi(2.0, 2, 2, 2)).args[2]
ℵ₀×ℵ₀ adjoint(::BandedMatrices.BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(vcat), Tuple{LinearAlgebra.Adjoint{Float64, SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(*), InfiniteArrays.InfUnitRange{Int64}, Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}, LinearAlgebra.Adjoint{Float64, SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(*), SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(-), SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(*), InfiniteArrays.InfUnitRange{Int64}, SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(+), Cumsum{Float64, 1, AbstractVector{Float64}}, SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(/), BroadcastVector{Float64, typeof(-), Tuple{BroadcastVector{Float64, typeof(/), Tuple{SubArray{Float64, 1, InfiniteLinearAlgebra.TridiagonalConjugationBand{Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}}, false}, SubArray{Float64, 1, InfiniteLinearAlgebra.TridiagonalConjugationBand{Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}}, false}}}}}, BroadcastVector{Float64, typeof(inv), Tuple{SubArray{Float64, 1, InfiniteLinearAlgebra.TridiagonalConjugationBand{Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}}, false}}}}}}, SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(*), InfiniteArrays.InfUnitRange{Int64}, SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(+), SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(/), BroadcastVector{Float64, typeof(-), Tuple{BroadcastVector{Float64, typeof(/), Tuple{InfiniteLinearAlgebra.TridiagonalConjugationBand{Float64}, InfiniteLinearAlgebra.TridiagonalConjugationBand{Float64}}}}}, BroadcastVector{Float64, typeof(inv), Tuple{InfiniteLinearAlgebra.TridiagonalConjugationBand{Float64}}}}, SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(*), BroadcastVector{Float64, typeof(inv), Tuple{InfiniteLinearAlgebra.TridiagonalConjugationBand{Float64}}}, MulAddAccumulate{Float64}}}}}, Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}}}, InfiniteArrays.OneToInf{Int64}}) with eltype Float64 with indices OneToInf()×OneToInf():
  ⋅   5.40256  -1.09972    ⋅         ⋅         ⋅         ⋅         ⋅         ⋅       …  
  ⋅    ⋅        8.0885   -1.72961    ⋅         ⋅         ⋅         ⋅         ⋅
  ⋅    ⋅         ⋅       10.4884   -2.26406    ⋅         ⋅         ⋅         ⋅
  ⋅    ⋅         ⋅         ⋅       12.7791   -2.74673    ⋅         ⋅         ⋅
  ⋅    ⋅         ⋅         ⋅         ⋅       15.0146   -3.19782    ⋅         ⋅
  ⋅    ⋅         ⋅         ⋅         ⋅         ⋅       17.2173   -3.62826    ⋅       …
  ⋅    ⋅         ⋅         ⋅         ⋅         ⋅         ⋅       19.3988   -4.04445
  ⋅    ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅       21.5657
 ⋮                                            ⋮                                      ⋱

works and uses CachedBroadcastVector.

[dlfivefifty]

What's the difference between a CachedBroadcastVector and a CachedVector{<:Any,<:BroadcastVector}?

We definitely do not want to cache every single operation

[DanielVandH]

I think the issue was with type inference for nested BroadcastVectors (or maybe just for caches of them) gives an eltype of Union{}.

Maybe if instead of using a cache for all the individual broadcasts I just build up all the BroadcastVectors and cache those

[dlfivefifty]

No the issue has nothing to do with type inference. It' that if you broadcast cached vectors then call to getindex have no way to tell the underlying vectors to call resizedata!. And we don't have the ability to reduce a getindex to a broadcast of the underlying cached data.

The proposal was to do a single type to wrap and cache the entire vector, not to introduce more generic lazy arrays that are composed.

[DanielVandH]

I think this is closer to what you want? I still need to figure out why e.g. (taken directly from the tests)

t = 2
P = SemiclassicalJacobi(t, -0.5, -0.5, -0.5)
Q = SemiclassicalJacobi(t,  0.5,  0.5,  0.5, P)
x = axes(P,1)
D = Derivative(x)
D = Q \ (D*P)

tries to indefinitely resize the cache though.