Adopt MutableArithmetics.jl ?

[fingolfin]

Perhaps we should reconsider adopting https://github.com/jump-dev/MutableArithmetics.jl which has gotten some traction by now. To be clear: this does not invalidate all work done one the AA/Nemo/... mutable API, but rather it can complement it, and help us to fix some issues with its design and implementation.

For some background, this was originally brought up in kalmarek/GroupsCore.jl#21 by @blegat on @kalmarek's GroupsCore repository. I've been hesitant to support this, and @thofma and @fieker seemed rather sceptical about it.

But since then I've had a chance to watch the JuliaCon 2021 presentation on MutableArithmetics, and looked at their code and documentation. And at the same time, I tried to use some of our own mutable arithmetic, and found many gaps in what we provide, and some bugs (I think -- PRs will be submitted once I had a chance to write some test for and clean up what I added to Nemo)

Overall I think there are a few advantages to be had:

• they carefully thought about dealing with aliasing issues, which we currently do not seem to do well in all cases (I'll have some upcoming PRs)
• they already provide and documents APIs like MA.mutable_copy and MA.copy_if_mutable, which we really should also offer for our types -- see also fix some deepcopy_internals: they need to call #955 were it is made clear that many of our deepcopy methods violate the Julia specification
• as far as I understand, their approach makes it easy (or at least easier) to find places where one "forgot" to implement mutable arithmetics, resulting in less than optimal code
• their @rewrite macro for automatically rewriting an arithmetic expression to use mutating APIs seems pretty sweet

In addition, it seems (see here) like adding support for this won't be that much work, we just have to provide a few "adapter" functions and then "our" mutable API and "theirs" could be used in parallel.

The main work would be to add MA.mutable_copy methods and use them, but that's something I believe we have to do anyway (if not MA.mutable_copy then via our own analogue).

The main drawback in my eyes is that this adds an external dependency on something out of our control. However, the JuMP people seem reasonable enough, so I'd be willing to risk it.

Anyway, this would still be a major undertaking; but it can be done in separate steps

• for now add our own analogues of mutable_copy and copy_if_mutable and start adding and using them; I'd try to make sure they match their MA counterparts exactly. Then if we later start using MutableArithmetics, it'll be a simple switch; and even if we don't use MA, then for people familiar with MA it'll be easier to use our code
• have a simple prototype which explores this idea that we can provide both APIs simultaneously with minimal effort; say for starts, just deal with fmpz in Nemo and see how that works out. If it does, one can use it more in Nemo (and Singular, and elsewhere, and eventually in AA). If not, we can give it up.

[thofma]

I don't doubt that this package is the best thing since sliced bread. What I am not sure about is whether it helps at all with our matrix pains. Indeed, it also gets everything wrong if matrix entries are allowed to alias each other. The problem with matrices is not just missing pretty syntax.

julia> A = rand(-10:10, n, n); A2 = big.(A);

julia> b = rand(-10:10, n); b2 = big.(b);

julia> c = rand(-10:10, n); c2 = big.(c);

julia> c2 = [BigInt(0), BigInt(0)]; MA.add_mul!!(c2, A2, b2); c2 == A2 * b2
true

julia> c2 .= 0; MA.add_mul!!(c2, A2, b2); c2 == A2 * b2
false

Our pains with the matrices are of a different nature. I am happy to discuss it (again). It boils down to what A[i, j] = b should do.

1. Disallow aliasing between matrix entries

Since this is just too easy to do wrong as a user (M[1, 1] = a and M[2, 1] = a), it is unreasonable to expect users to think of this. To mitigate this, we could

• copy on write always, i.e., turn all M[i, j] = a) internally into M[i, j] = deepcopy(a) (mutable_copy(a)).
• check if the new entry aliases another one and copy if necessary (it is not clear, whether "aliasing" can actually be detected or not (is it a recursive concept?).

To circumvent the performance problems, one could introduce

setindex!(A, i, j, a, copy = false)

or

@mcopy A[i, j] = a

and use this everywhere in library code. But this means that user code dealing with matrix entries will always be slow, since no user will ever use this complicated syntax.

2. Allow aliasing of matrix entries

This means no copy on write (ever) and the user can set entries in whichever way he wants without performance penalty. The disadvantage is, that mutable_operation!(M[i, j], x, y) is illegal in general. Implementing matrix methods efficiently gets quite tedious and is prone to errors (see the bugs you found).

3. Conclusion

Over the last years I had many discussions with Bill, Claus and Dan. At one point, Bill tried out the version 1. and the setindex! with the copy keyword. I think we decided that it was too cumbersome to use (we did not think of the macro syntax back then).

[blegat]

I agree that having aliases as entries of arrays can be a confusing mistake:

julia> c = big.([0, 0])
2-element Vector{BigInt}:
 0
 0

julia> c[1] === c[2]
false

julia> c .= 0
2-element Vector{BigInt}:
 0
 0

julia> c[1] === c[2]
true

julia> c = zeros(BigInt, 2)
2-element Vector{BigInt}:
 0
 0

julia> c[1] === c[2]
true

In JuMP, the way we deal with this is that the user calling the MutableArithmetics (MA) API is expected to be an advanced user that knows about these aliasing issues.
JuMP macros uses the rewriting mechanism of MA (MA.@rewrite) to rewrite the expressions of its macro to use the MA API but it never assumes that anything the user gives is mutable and calls MA.copy_if_mutable to the first term of the expression if we will call mutable operation on it and mutable_copy is recursive so there is no alias in the end.
The linear and sparse algebra implemented in MA assumes that there is no aliasing in the output.
Unless the user explicitly calls the MA API wrongly as in your example above, it will be the case.
MA.mul(A2, b2) in your example above will create an output vector with no aliasing and then do the efficient mutable-aware matrix multiplication.
MA.mul!(c2, A2, b2) will first call MA.zero!(c2) which will create non-alias zeros on c2 so it will work as well.