How make .+ work with A::Vector{Float64} and B::OffsetArray(::Vector{Float64}, 0:1)

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]); 

    return result

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 *)

Hey elvispy! Welcome to the julia discourse.

For your use-case you should just use:

julia> A .+ (B.parent)
2-element Vector{Float64}:

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:

julia> A
2-element Vector{Float64}:

julia> A .+ B
ERROR: DimensionMismatch: arrays could not be broadcast to a common size; got a dimension with lengths 2 and 2
 [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}:

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:

julia> A .+ C
2-element OffsetArray(::Vector{Float64}, 1:2) with eltype Float64 with indices 1:2:

This discussion seems to be related: Check that axes start with 1 for AbstractRange operations by timholy · Pull Request #30950 · JuliaLang/julia · GitHub

1 Like

Thank you! This solution is really practical to me.

On another note, could you tell me why the following code works

julia> collectPl.(0.1, lmax = 1).parent .* collectPl.(0.1, lmax = 1).parent
2-element Vector{Float64}:

but not

julia> LegendrePolys(x::Float64; lmax::Integer) = (collectPl(x, lmax = lmax)).parent
julia> @. LegendrePolys(0.1, lmax = 1) * LegendrePolys(0.1, lmax=1)
MethodError: no method matching *(::Vector{Float64}, ::Vector{Float64})

or use parent function instead of field access:
A .+ parent(B)


Weird! I would expect that to work! Let’s see if someone else can answer this one!

@kristoffer.carlsson shared this two discussions on Slack about why such a thing does not work:

Broadcasted rand.(::Int64), similar to rand(::Int64), but errs with .< · Issue #32338 · JuliaLang/julia · GitHub and
Cannot broadcast when rand(n) is involved · Issue #32744 · JuliaLang/julia · GitHub

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).

1 Like

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.

1 Like

And noting that the following will avoid problem with @. macro:
LegendrePolys(0.1, lmax=1) .* LegendrePolys(0.1, lmax=1)