Einsum.jl icon indicating copy to clipboard operation
Einsum.jl copied to clipboard

Does not work with Unitful.jl?

Open singularitti opened this issue 1 year ago • 0 comments

For example, I have the following sum rule:

@einsum T′[i, j] := T[k, l] * Q[i, k] * Q[j, l]

and the macro-expanded code for above is:

julia> @macroexpand @einsum T′[i, j] := T[k, l] * Q[i, k] * Q[j, l]
quote
    #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:207 =#
    local var"##T#342" = promote_type(eltype(T), eltype(Q), eltype(Q))
    #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:208 =#
    T′ = Array{var"##T#342"}(undef, size(Q, 1), size(Q, 1))
    #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:209 =#
    if size(T, 2) == size(Q, 2)
        nothing
    else
        Base.throw(Base.AssertionError("size(T, 2) == size(Q, 2)"))
    end
    if size(T, 1) == size(Q, 2)
        nothing
    else
        Base.throw(Base.AssertionError("size(T, 1) == size(Q, 2)"))
    end
    #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:212 =#
    let i, j, k, l
        #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:213 =#
        begin
            $(Expr(:inbounds, true))
            local var"#170#val" = for j = 1:size(Q, 1)
                        #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:278 =#
                        for i = 1:size(Q, 1)
                            #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:278 =#
                            begin
                                #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:176 =#
                                local var"##s#343" = zero(var"##T#342")
                                #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:177 =#
                                for l = 1:size(T, 2)
                                    #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:278 =#
                                    for k = 1:size(T, 1)
                                        #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:278 =#
                                        var"##s#343" += T[k, l] * Q[i, k] * Q[j, l]
                                        #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:279 =#
                                    end
                                    #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:279 =#
                                end
                                #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:178 =#
                                T′[i, j] = var"##s#343"
                            end
                            #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:279 =#
                        end
                        #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:279 =#
                    end
            $(Expr(:inbounds, :pop))
            var"#170#val"
        end
    end
    #= ~/.julia/packages/Einsum/AVMOj/src/Einsum.jl:216 =#
    T′
end

However, in this line:

local var"##T#342" = promote_type(eltype(T), eltype(Q), eltype(Q))

It does not consider that T, Q may have units defined in Unitful.jl. For example, say T is a matrix filled with quantities like 1u"m", and Q is filled with quantities like 2.0u"kg". Then it will return type Quantity{Float64}:

julia> S = promote_type(typeof(1u"m"), typeof(2.0u"kg"), typeof(2.0u"kg"))
Quantity{Float64}

The final eltype of the results will have dimension m * kg * kg,

As I mentioned in https://github.com/PainterQubits/Unitful.jl/issues/616, S will be used in

local var"##s#343" = zero(var"##T#342")

as the initial summation value. However, zero for Quantity{Float64} is not defined. So it will throw an ArgumentError.

singularitti avatar Feb 16 '23 10:02 singularitti