In the following code

``````import Base: ^

struct cccbbb
v::Matrix{Float64}
end

^(cb::cccbbb, p::Int) = cb.v^p

n = 5
r = LinRange(0.0, 10.0, 5)
cb = cccbbb(rand(n,n))

tmp = cb.^2
@show typeof(cb.^2), typeof(cb.^2 ./r), typeof(tmp./r)

``````

the result is

``````(typeof(cb .^ 2), typeof(cb .^ 2 ./ r), typeof(tmp ./ r)) = (Matrix{Float64}, Vector{Matrix{Float64}}, Matrix{Float64})
``````

The type of (cb .^ 2 ./ r) is not as expected. Is this a bug?

Yeah, that’s a little confusing, but I think the result is correct. Here’s a slightly simplified version of your example:

``````import Base: ^

struct A
m::Matrix{Float64}
end

^(a::A, p::Int) = a.m ^ p

m = [1 2
3 4]

a = A(m)
``````
``````julia> m ^ 2
2×2 Array{Int64,2}:
7  10
15  22

julia> m .^ 2
2×2 Array{Int64,2}:
1   4
9  16

julia> m .^ 2 ./ [1, 2]
2×2 Array{Float64,2}:
1.0  4.0
4.5  8.0

julia> a .^ 2 ./ [1, 2]
2-element Array{Array{Float64,2},1}:
[7.0 10.0; 15.0 22.0]
[3.5 5.0; 7.5 11.0]
``````

The `Base.broadcastable(a::A) = Ref(a)` line tells Julia to treat an object of type `A` as a scalar for broadcasting. So

``````a .^ 2 ./ [1, 2]
``````

is effectively equivalent to

``````[(a^2)/1, (a^2)/2]
``````

And they way you’ve defined `^` on your new type means that it is doing a matrix power (`m * m`), not an element by element power.

Can you elaborate on what output you expected and/or what output you are trying to achieve?

3 Likes

``````julia> b = a.^2
2×2 Matrix{Float64}:
7.0  10.0
15.0  22.0
``````

and

``````julia> b ./ [ 1, 2]
2×2 Matrix{Float64}:
7.0  10.0
7.5  11.0
``````

The result is the Matrix{Float64} type as expected. But the

``````a .^ 2 ./ [1, 2]
``````

type is

``````2-element Array{Array{Float64,2},1}
``````

which is not expected.

Can the `Base.broadcastable(a::A) = Ref(a)` just affect the operator `^`?

What we’re trying to tell you is that `b ./ [1, 2]` and `a .^ 2 ./ [1, 2]` are two distinct computations. An important feature of the dot syntax is that vectorized functions are fused into a single loop. In the second case, you’re thinking Julia would first broadcast exponentiation to obtain the matrix, but as @CameronBieganek noted in his reply, the operations are fused into `[(a^2)/1, (a^2/2)]`. This is the expected behavior, since element-wise exponentiation is different from exponentiation after all (take, matrices, for example). The behavior here is more confusing than usual because you’re wrapping a matrix, but also tell Julia to treat your struct as a scalar.

If you want the same result with the first case, try `a ^ 2 ./ [1,2]` instead.

3 Likes