Fix implementation of mul! for AbstractMatrix and AbstractVector

[odow]

Closes #251

The previous implementations were incorrect because they did not account for non-Base.OneTo axes.

[odow]

This now results in:

julia> using JuMP

julia> model = Model()
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.

julia> @variable(model, x[0:1, 1:2], container = DenseAxisArray)
2-dimensional DenseAxisArray{VariableRef,2,...} with index sets:
    Dimension 1, 0:1
    Dimension 2, Base.OneTo(2)
And data, a 2×2 Matrix{VariableRef}:
 x[0,1]  x[0,2]
 x[1,1]  x[1,2]

julia> @variable(model, y[[:a, :b]], container = DenseAxisArray)
1-dimensional DenseAxisArray{VariableRef,1,...} with index sets:
    Dimension 1, [:a, :b]
And data, a 2-element Vector{VariableRef}:
 y[a]
 y[b]

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> A * x
2×2 Matrix{AffExpr}:
 x[0,1] + 2 x[1,1]    x[0,2] + 2 x[1,2]
 3 x[0,1] + 4 x[1,1]  3 x[0,2] + 4 x[1,2]

julia> A * y
2-element Vector{AffExpr}:
 y[a] + 2 y[b]
 3 y[a] + 4 y[b]

julia> x * A
2×2 Matrix{AffExpr}:
 x[0,1] + 3 x[0,2]  2 x[0,1] + 4 x[0,2]
 x[1,1] + 3 x[1,2]  2 x[1,1] + 4 x[1,2]

julia> y' * A
1×2 adjoint(::Vector{AffExpr}) with eltype AffExpr:
 y[a] + 3 y[b]  2 y[a] + 4 y[b]

julia> x * y
ERROR: MethodError: no method matching similar(::JuMP.Containers.DenseAxisArray{VariableRef, 1, Tuple{Vector{Symbol}}, Tuple{JuMP.Containers._AxisLookup{Dict{Symbol, Int64}}}}, ::Type{QuadExpr}, ::Tuple{UnitRange{Int64}})

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Comparison is base (69a675d) 90.93% compared to head (3f35060) 90.86%.
Report is 1 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #252      +/-   ##
==========================================
- Coverage   90.93%   90.86%   -0.08%     
==========================================
  Files          23       23              
  Lines        2163     2157       -6     
==========================================
- Hits         1967     1960       -7     
- Misses        196      197       +1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[blegat]

Is dropping @inbounds intentional ? I guess it's because you assume we are manipulating expensive types so the cost of add_mul will greatly outweigh the cost of getindex ?
We should double check that there is no regression, e.g., with the benchmarks of the README

[odow]

Is dropping @inbounds intentional?

Yes. I'm just not in favor of using @inbounds anywhere. The small performance gains do not outweigh the potential for bugs.

[odow]

We should double check that there is no regression, e.g., with the benchmarks of the README

No. The current implementation was wrong. Let's focus on correctness over performance.

[odow]

Bump @blegat

[odow]

Bump bump @blegat

[blegat]

We should double check that the README example and examples of the paper still work

[odow]

README example still works. It is:

julia> using LinearAlgebra

julia> trial = @benchmark LinearAlgebra.mul!($c2, $A2, $b2)
BenchmarkTools.Trial: 378 samples with 1 evaluation.
 Range (min … max):   5.647 ms … 290.737 ms  ┊ GC (min … max):  0.00% … 72.52%
 Time  (median):      6.052 ms               ┊ GC (median):     0.00%
 Time  (mean ± σ):   13.577 ms ±  35.051 ms  ┊ GC (mean ± σ):  31.47% ± 11.47%

  █                                                             
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▄▁▁▁▁▁▁▁▁▁▁▁▁▄▄▁▁▁▁▄▁▁▁▄▁▅▄▁▁▅▄▅ ▅
  5.65 ms       Histogram: log(frequency) by time       173 ms <

 Memory estimate: 3.67 MiB, allocs estimate: 198144.

julia> trial2 = @benchmark MA.add_mul!!($c2, $A2, $b2)
BenchmarkTools.Trial: 2637 samples with 1 evaluation.
 Range (min … max):  1.292 ms … 35.666 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.628 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.888 ms ±  1.027 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▃                                                          
  ██▇▅▅▄▄▄▃▃▃▃▃▃▃▃▃▃▃▄▄▄▅▅▅▄▄▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂ ▃
  1.29 ms        Histogram: frequency by time        3.91 ms <

 Memory estimate: 48 bytes, allocs estimate: 3.

julia> buffer = MA.buffer_for(MA.add_mul, typeof(c2), typeof(A2), typeof(b2))
0

julia> trial3 = @benchmark MA.buffered_operate!!($buffer, MA.add_mul, $c2, $A2, $b2)
BenchmarkTools.Trial: 3354 samples with 1 evaluation.
 Range (min … max):  1.292 ms …   3.408 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.372 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.487 ms ± 270.838 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▅▇▅▅▄▄▄▃▃▃▂▁▁▁▁▁▁▁▁▁▁                                      ▁
  ██████████████████████████████▇▇▇▇▇▆▆▆▇▅▇▇▆▅▇▆▅▃▅▇▃▇▆▆▅▄▆▆▅ █
  1.29 ms      Histogram: log(frequency) by time      2.49 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.