Matrix multiplication on specified indices

[mcabbott]

It would be nice to be able to specify which indices to multiply on, instead of trusting their position (and giving an error if they don't match). Like sum(A, dims=:k), this should let you write code which only works on NamedDimsArrays, but doesn't care about the order of the underlying Array.

This PR is a first attempt to add this, to see what you think. For now it's written *(:k, A, B) but perhaps *(A,B; dims=:k) or even Contract{(:k,)}(A,B) would be better. The docstring example uses #62 to create infix *ⱼ easily. Ideally something similar would work for \ or / too, #6.

Haven't thought about speed, constant propagation etc. Nor tried to make it work for ndims(A)==3 etc, which possibly should be handled by @requireing OMEinsum.jl, or something.

[oxinabox]

we definately do want some kind of einsum/contract API.
I think it would be OK to just add a dependency on an einsum package,
assuming we can identify one that is fast and general enough.

[codecov]

Codecov Report

Merging #63 into master will increase coverage by 0.05%.
The diff coverage is 82.6%.

@@            Coverage Diff             @@
##           master      #63      +/-   ##
==========================================
+ Coverage   86.42%   86.47%   +0.05%     
==========================================
  Files           8        8              
  Lines         221      244      +23     
==========================================
+ Hits          191      211      +20     
- Misses         30       33       +3

Impacted Files	Coverage Δ
src/functions_math.jl	92.18% <82.6%> (-5.38%)	⬇️
src/functions_dims.jl	100% <0%> (+6.66%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 41d55be...ca56785. Read the comment docs.

[mcabbott]

OK here's a slightly more refined version, still only up to matrices with one contraction.

Is this the right sort of thing to test whether constant propagation is working?

julia> @code_warntype ((B,AB) -> *(:b, B, AB))(B, AB)
Variables
  #self#::Core.Compiler.Const(var"##9#10"(), false)
  B::NamedDimsArray{(:b,),Float64,1,Array{Float64,1}}
  AB::NamedDimsArray{(:a, :b),Float64,2,Array{Float64,2}}

Body::NamedDimsArray{(:a,),Float64,1,Array{Float64,1}}
1 ─ %1 = (:b * B * AB)::NamedDimsArray{(:a,),Float64,1,Array{Float64,1}}
└──      return %1

It seems equally confident if given :z (which is an error). For a case like @code_warntype (AB -> *(:b, AB, AB))(AB) where the :b matters it gives Union{NamedDimsArray{(:b, :b),...}, NamedDimsArray{(:a, :a),... which I guess is a bad sign.

[mcabbott]

More complicated contractions are now handled by OMEinsum, which has a pretty straightforward non-macro interface. No tests, but testing would be welcomed.

[mcabbott]

Latest draft treats correctly arrays of arrays, and deletes the attempt to use OMEinsum.jl for more general contractions. I now think that overloading * like this should only do things * can do, i.e. should insert transpose as required, but not allow contraction of higer-dim arrays, nor contraction of more than two arrays. Other things can be some other contract function.

Perhaps it would be nicer to write *(A, B; dims=:j), but making that work will involve re-writing how * acts on named arrays (without dims), instead of just building on top. And I can’t figure out how to make keywords behave well with constant propagation.

[oxinabox]

At some point we should come back to this idea,
but I think it will want some careful thinking about how to do it.

[mcabbott]

Yes. I have a version I'm happier with at https://github.com/mcabbott/NamedPlus.jl/blob/master/src/mul.jl but should use it a bit first.