The julia DSL for this is closer to what you want: "M_ij = A_ik * B_kj" reads

    @einsum M[i,j] := A[i,k] * B[k,j]

Or faster with @tensor from https://github.com/Jutho/TensorOperations.jl & the same notation. I'm not aware of any attempt to distinguish co-vectors; you can have g = Diagonal([-1,1,1,1]) and insert A[i,k]g[k,k′]B[k′,j] and this is handled efficiently.