The LegendrePolynomials.jl package has a function collectPl which returns an OffsetArray Object.
It appears that, for some reason, the dot operations are not working properly with this type of data:
julia> A = [1.0, 2.0];
julia> B = collectPl(0.1, lmax = 1) # Just two arrays of the same size
julia> A .+ B
ERROR: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 2 and 2")
I was able to dispatch another method for the broadcast operator when the second input has the desired data type:
function broadcast(f::Any, amplitudes::Vector{Float64}, A::T) where T <: OffsetArrays.OffsetArray
result = fill(0.0, length(amplitudes));
for (ii, jj) in zip(eachindex(amplitudes), eachindex(A))
result[ii] = f(amplitudes[ii], A[jj]);
end
return result
end
But the same error appears, probably because .+ is calling the usual broadcast function.
If, however, I try to dispatch another method for .+, Julia complaints that this is not a valid function name.
julia> .+(amplitudes::Vector{Floag64}, A::T) where T <: OffsetArrays.OffsetArray = broadcast(+, amplitudes, A)
syntax: invalid function name ".+" around ...
How can I implement in not that many lines of code the two dot operations that I’m interested in? (+ and *)
julia> A .+ (B.parent)
2-element Vector{Float64}:
2.0
2.1
Reasoning
Let me speculate.
I would think that the reason this is not working is that one of the arrays is an OffsetArray with indices 0 and 1. So if you interpret the broadcasted addition .+ as an “element-wise addition” you would like to add A[i] + B[i] for all the indices in A in B. But in this case A has indices 0 and 1 and B has indices 1 and 2. So even though they have the same number of elements they don’t “match”. I agree that the error message is not very useful in this case, but making it work by force in general could be a foot gun.
If you are sure that the operation you want to do is the element-wise addition disregarding the fact that both vectors have different indices you could acces the underlying array of B, which you can do with B.parent:
julia> using LegendrePolynomials
julia> A = [1.0, 2.0];
julia> B = collectPl(0.1, lmax = 1) # Just two arrays of the same size
2-element OffsetArray(::Vector{Float64}, 0:1) with eltype Float64 with indices 0:1:
1.0
0.1
julia> A
2-element Vector{Float64}:
1.0
2.0
julia> A .+ B
ERROR: DimensionMismatch: arrays could not be broadcast to a common size; got a dimension with lengths 2 and 2
Stacktrace:
[1] _bcs1
@ ./broadcast.jl:516 [inlined]
[2] _bcs
@ ./broadcast.jl:510 [inlined]
[3] broadcast_shape
@ ./broadcast.jl:504 [inlined]
[4] combine_axes
@ ./broadcast.jl:499 [inlined]
[5] instantiate
@ ./broadcast.jl:281 [inlined]
[6] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(+), Tuple{Vector{Float64}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}}})
@ Base.Broadcast ./broadcast.jl:860
[7] top-level scope
@ REPL[8]:1
julia> A .+ B.parent
2-element Vector{Float64}:
2.0
2.1
The problem then is not that broadcasting does not work for OffsetArrays. For example if we re-wrap the contents of B in an OffsetArray with indices 1:2 (the same indices as in A) things work as you would expect.
# Code continues from previous snippet
julia> using OffsetArrays
julia> C = OffsetArray(B.parent, 1:2)
2-element OffsetArray(::Vector{Float64}, 1:2) with eltype Float64 with indices 1:2:
1.0
0.1
julia> A .+ C
2-element OffsetArray(::Vector{Float64}, 1:2) with eltype Float64 with indices 1:2:
2.0
2.1
Basically @. LegendrePolys(0.1, lmax = 1) * LegendrePolys(0.1, lmax=1) transforms to LegendrePolys.(0.1, lmax = 1) .* LegendrePolys.(0.1, lmax=1) (see the dots also between LegendrePolys and their arguments).
Indeed, accessing internal fields is not generally recommended. There is also OffsetArrays.no_offset_view, not sure which is better suited between this and the parent function.