Add new rewrite that does not make strong assumptions on result type #170
Conversation
Codecov ReportBase: 83.66% // Head: 84.75% // Increases project coverage by
Additional details and impacted files@@ Coverage Diff @@
## master #170 +/- ##
==========================================
+ Coverage 83.66% 84.75% +1.08%
==========================================
Files 20 21 +1
Lines 1898 2033 +135
==========================================
+ Hits 1588 1723 +135
Misses 310 310
Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here. ☔ View full report at Codecov. |
That's a bit surprising but as this was never benchmarked properly, we might be doing something inefficient somewhere indeed. |
I need to check with JuMP types. It's probably much more efficient if the cost of creating a new object is high. |
This comment was marked as outdated.
This comment was marked as outdated.
|
This PR also closes #44 julia> MA2.@rewrite(1 .+ sum(1 for i in 1:0) * 1^2)
1
julia> MA.@rewrite(1 .+ sum(1 for i in 1:0) * 1^2)
ERROR: ArgumentError: reducing over an empty collection is not allowedjulia> @macroexpand MA2.@rewrite(1 .+ sum(1 for i in 1:0) * 1^2)
quote
#= /Users/oscar/.julia/dev/MutableArithmetics/src/new_rewrite.jl:33 =#
let
#= /Users/oscar/.julia/dev/MutableArithmetics/src/new_rewrite.jl:34 =#
begin
#= /Users/oscar/.julia/dev/MutableArithmetics/src/new_rewrite.jl:30 =#
var"#9192###2317" = MutableArithmetics.Zero()
for i = 1:0
#= /Users/oscar/.julia/dev/MutableArithmetics/src/new_rewrite.jl:233 =#
var"#9192###2317" = (MutableArithmetics.operate!!)(MutableArithmetics.MutableArithmetics2.:+, var"#9192###2317", 1)
end
var"#9193###2318" = 1 ^ 2
var"#9194###2316" = (MutableArithmetics.operate!!)(MutableArithmetics.add_mul, MutableArithmetics.Zero(), var"#9192###2317", var"#9193###2318")
var"#9195###2319" = 1 .+ var"#9194###2316"
end
#= /Users/oscar/.julia/dev/MutableArithmetics/src/new_rewrite.jl:35 =#
var"#9195###2319"
end
end |
|
The @rewrite(sum(i * x[i] for i in 1:n) * 2)
# becomes
y = MA.Zero()
for i in 1:n
y = MA.operate!!(MA.add_mul, y, i, x[i])
end
y = MA.operate!!(MA.add_mul, MA.Zero(), y, 2)So the last operation allocates a new expression. Ideally, we'd do something like y = MA.operate!!(*, y, 2)but this isn't what Edit: I've done this and updated the benchmarks |
odow
force-pushed
the
od/rewrite-2
branch
4 times, most recently
from
42f7017 to
4e01cde
Compare
|
I think it's safe to say that this is now even slightly more efficient than the existing stats = DataFrames.combine(
DataFrames.groupby(df, [:type, :method]),
:result => Statistics.mean => :mean,
:result => Statistics.median => :median,
:result => (x -> exp(Statistics.mean(log.(x)))) => :geomean,
)
julia> DataFrames.sort!(stats, [:method, :type])
15×5 DataFrame
Row │ type method mean median geomean
│ String String Float64 Float64 Float64
─────┼─────────────────────────────────────────────────────────
1 │ BigFloat IM 0.644393 0.543036 0.574129
2 │ BigInt IM 0.805905 0.617568 0.662454
3 │ Float64 IM 0.801565 0.901345 0.698718
4 │ Int64 IM 0.781128 0.955712 0.688505
5 │ JuMP.VariableRef IM 65.746 1.19753 5.23297
6 │ BigFloat MA 1.0 1.0 1.0
7 │ BigInt MA 1.0 1.0 1.0
8 │ Float64 MA 1.0 1.0 1.0
9 │ Int64 MA 1.0 1.0 1.0
10 │ JuMP.VariableRef MA 1.0 1.0 1.0
11 │ BigFloat MA2 0.798462 0.96 0.759529
12 │ BigInt MA2 0.80857 0.977416 0.76787
13 │ Float64 MA2 0.965268 0.997299 0.959214
14 │ Int64 MA2 1.0056 0.989867 0.998947
15 │ JuMP.VariableRef MA2 0.970589 0.987654 0.968822The main decision is whether to add it as a new feature, or whether to replace the existing |
|
So one major thing this PR doesn't do is rewrite broadcasts. I need to:
|
|
@blegat can you find the benchmarks that demonstrated the need to rewrite broadcasts?
|
|
I don't see that |
Yes. I fixed this by rewriting |
The changes break JuMP's tests by changing printing and the error types that are thrown, so I'd be in favor of keeping it as a new feature. |
If that is the case and given that MA2 seems strictly an improvement of MA(?), why not do a full replace of MA and tag version v2.0, then JuMP can update to the new version in the next release. Unless there is still some uncertainty around the design or possible drawbacks of MA2? |
|
Releasing v2.0 seems a bit drastic. There are 748 dependent packages, so releasing a breaking change has cost: https://juliahub.com/ui/Packages/MutableArithmetics/EoEec/1.0.5. Perhaps we could have
Main drawbacks are a lack of testing. This is the problem with Julia's interface design; it's hard to tell which methods we are missing without running examples. We know that the current design works for people; we can't be 100% confident that the new design also works for the same use-cases. |
Wow, I see your point. I had no idea this package was so widely used. |
odow
changed the title
Fix rewrite with views Update Support broadcasting Fix handling of * Fix generators Add tests for repeated sums Fix formatting Updates Place new rewrite behind opt-in kwarg More coverage Update docstrings
|
Okay, I think I'm happy for this to be merged. It's strictly a feature addition hidden behind an opt-in flag. |
|
Hmm. Found a few more edge cases that we need to fix before merging. julia> using JuMP
julia> using LinearAlgebra
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[1:3])
3-element Vector{VariableRef}:
x[1]
x[2]
x[3]
julia> A = LowerTriangular(rand(3, 3))
3×3 LowerTriangular{Float64, Matrix{Float64}}:
0.507834 ⋅ ⋅
0.178216 0.774911 ⋅
0.498208 0.438567 0.165303
julia> const MA = JuMP._MA
MutableArithmetics
julia> MA.@rewrite(A * x, move_factors_into_sums = true)
3-element Vector{AffExpr}:
0.5078344563174482 x[1]
0.1782159773716443 x[1] + 0.7749111979573089 x[2]
0.49820804791789297 x[1] + 0.4385670932233323 x[2] + 0.16530329346630457 x[3]
julia> MA.@rewrite(A * x, move_factors_into_sums = false)
ERROR: MethodError: Cannot `convert` an object of type QuadExpr to an object of type AffExpr
Closest candidates are:
convert(::Type{GenericAffExpr{T, V}}, ::GenericAffExpr{T, V}) where {T, V} at /Users/oscar/.julia/packages/JuMP/Z1pVn/src/aff_expr.jl:538
convert(::Type{GenericAffExpr{T, V}}, ::GenericAffExpr{S, V}) where {S, T, V} at /Users/oscar/.julia/packages/JuMP/Z1pVn/src/aff_expr.jl:545
convert(::Type{GenericAffExpr{T, V}}, ::Union{Number, UniformScaling}) where {T, V} at /Users/oscar/.julia/packages/JuMP/Z1pVn/src/aff_expr.jl:534
...
Stacktrace:
[1] setindex!
@ ./array.jl:845 [inlined]
[2] lmul!(A::LowerTriangular{AffExpr, Matrix{AffExpr}}, B::Vector{AffExpr})
@ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/triangular.jl:951
[3] *(A::LowerTriangular{Float64, Matrix{Float64}}, B::Vector{VariableRef})
@ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/triangular.jl:1657
[4] top-level scope
@ ~/.julia/dev/MutableArithmetics/src/rewrite.jl:317 |
|
By not doing fancy rewrites, this now fixes (or at least, improves) the examples from #66 as well: julia> using JuMP
julia> using LinearAlgebra
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[1:3])
3-element Vector{VariableRef}:
x[1]
x[2]
x[3]
julia> A = LowerTriangular(rand(3, 3))
3×3 LowerTriangular{Float64, Matrix{Float64}}:
0.840885 ⋅ ⋅
0.944976 0.178466 ⋅
0.440437 0.252618 0.637874
julia> c = [1, 2, 3]
3-element Vector{Int64}:
1
2
3
julia> const MA = JuMP._MA
MutableArithmetics
julia> MA.@rewrite(A * x, move_factors_into_sums = true)
3-element Vector{AffExpr}:
0.8408850607032472 x[1]
0.944975752029531 x[1] + 0.17846559026575726 x[2]
0.44043661564033143 x[1] + 0.25261817336168857 x[2] + 0.6378741155843315 x[3]
julia> MA.@rewrite(A * x, move_factors_into_sums = false)
3-element Vector{AffExpr}:
0.8408850607032472 x[1]
0.944975752029531 x[1] + 0.17846559026575726 x[2]
0.44043661564033143 x[1] + 0.25261817336168857 x[2] + 0.6378741155843315 x[3]
julia> MA.@rewrite((A * x) .* c, move_factors_into_sums = false)
3-element Vector{AffExpr}:
0.8408850607032472 x[1]
1.889951504059062 x[1] + 0.3569311805315145 x[2]
1.3213098469209943 x[1] + 0.7578545200850657 x[2] + 1.9136223467529945 x[3]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[1:2])
2-element Vector{VariableRef}:
x[1]
x[2]
julia> @variable(model, y[1:2, 1:2])
2×2 Matrix{VariableRef}:
y[1,1] y[1,2]
y[2,1] y[2,2]
julia> A = ones(2, 2)
2×2 Matrix{Float64}:
1.0 1.0
1.0 1.0
julia> MA.@rewrite(y - Diagonal(1 .- x) * A, move_factors_into_sums = true)
ERROR: MethodError: Cannot `convert` an object of type QuadExpr to an object of type AffExpr
Closest candidates are:
convert(::Type{GenericAffExpr{T, V}}, ::GenericAffExpr{T, V}) where {T, V} at /Users/oscar/.julia/packages/JuMP/Z1pVn/src/aff_expr.jl:538
convert(::Type{GenericAffExpr{T, V}}, ::GenericAffExpr{S, V}) where {S, T, V} at /Users/oscar/.julia/packages/JuMP/Z1pVn/src/aff_expr.jl:545
convert(::Type{GenericAffExpr{T, V}}, ::Union{Number, UniformScaling}) where {T, V} at /Users/oscar/.julia/packages/JuMP/Z1pVn/src/aff_expr.jl:534
...
Stacktrace:
[1] setindex!
@ ./array.jl:845 [inlined]
[2] setindex!
@ ./multidimensional.jl:645 [inlined]
[3] macro expansion
@ ./broadcast.jl:984 [inlined]
[4] macro expansion
@ ./simdloop.jl:77 [inlined]
[5] copyto!
@ ./broadcast.jl:983 [inlined]
[6] copyto!
@ ./broadcast.jl:936 [inlined]
[7] materialize!
@ ./broadcast.jl:894 [inlined]
[8] materialize!
@ ./broadcast.jl:891 [inlined]
[9] lmul!(D::Diagonal{AffExpr, Vector{AffExpr}}, B::Matrix{AffExpr})
@ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/diagonal.jl:212
[10] *(D::Diagonal{AffExpr, Vector{AffExpr}}, A::Matrix{Float64})
@ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/diagonal.jl:201
[11] sub_mul
@ ~/.julia/dev/MutableArithmetics/src/MutableArithmetics.jl:39 [inlined]
[12] operate
@ ~/.julia/dev/MutableArithmetics/src/interface.jl:199 [inlined]
[13] operate_fallback!!
@ ~/.julia/dev/MutableArithmetics/src/interface.jl:571 [inlined]
[14] operate!!(::typeof(MutableArithmetics.sub_mul), ::Matrix{VariableRef}, ::Diagonal{AffExpr, Vector{AffExpr}}, ::Matrix{Float64})
@ MutableArithmetics ~/.julia/dev/MutableArithmetics/src/rewrite.jl:89
[15] top-level scope
@ ~/.julia/dev/MutableArithmetics/src/rewrite.jl:322
julia> MA.@rewrite(y - Diagonal(1 .- x) * A, move_factors_into_sums = false)
2×2 Matrix{AffExpr}:
y[1,1] + x[1] - 1 y[1,2] + x[1] - 1
y[2,1] + x[2] - 1 y[2,2] + x[2] - 1 |
|
Since this is opt-in, perhaps we should merge this, and then continue to tweak in follow-up PRs. |
|
I agree, let's merge so that we can see changes to this more easily |
Closes #44
Closes #169
This branch is interesting to try out.
First, the current
@rewriteproduces different results to what we would obtain using the non-mutating API. I'm surprised this hasn't come up before.Second, this PR has better performance than the current
@rewrite, at least for theBigXXXtypes. The differences can be quite large (2x). Noticeably, just evaluating things is faster, likely because of type stability?Benchmarks
Code
Results
Values in
()are the ratio relative to "MA" column.