copy(::MatrixElem) can return a different type

[rfourquet]

It was decided recently (#638) that similar(mat, ring) doesn't necessarily return a matrix of the the same type as mat. With the current implementation of copy, which calls out to similar, it also means that copy(mat) may be of a different type than mat. This also means that deepcopy can fail, because Base.deepcopy puts a typeassert on the return type (assumed to be of the same type as the input). For example:

julia> using Nemo

Welcome to Nemo version 0.14.3

Nemo comes with absolutely no warranty whatsoever


julia> m = one(Generic.MatSpace{fmpz}(FlintZZ, 2, 2, false))
[1  0]
[0  1]

julia> typeof(m)
AbstractAlgebra.Generic.MatSpaceElem{fmpz}

julia> copy(m)
[1  0]
[0  1]

julia> typeof(ans)
fmpz_mat

julia> deepcopy(m)
ERROR: TypeError: in typeassert, expected AbstractAlgebra.Generic.MatSpaceElem{fmpz}, got fmpz_mat
Stacktrace:
 [1] deepcopy(::AbstractAlgebra.Generic.MatSpaceElem{fmpz}) at ./deepcopy.jl:30
 [2] top-level scope at REPL[17]:1

I find copy(mat) not returning the same type as mat slightly annoying, but this is consistent with some other matrix constructors which will favor fmpz_mat when the elem_type is fmpz. What do you think?
On the other hand, deepcopy should be fixed to not fail.

[thofma]

Yes, copy and deepcopy should not change the type. This is a bug.

Whatever similar does, it should not affect copy or deepcopy.

[wbhart]

Agreed.

[thofma]

Do we have to fix it also in AbstractAlgebra?

[rfourquet]

When Nemo is not loaded, I'm not aware of instances of similar where it returns a different type (so this bug doesn't happen). But depending on the chosen fix, it may have to be applied in AbstractAlgebra.

It seems we need a function which returns a matrix of the same type as the input matrix, with the same dimensions (both for copy and deepcopy). But similar has a slightly different behavior now, so we may need a new function... or change back similar defined in Nemo to do just that.

[rfourquet]

Besides copy (now fixed), a bunch of functions won't return the same type as the input:

julia> m = one(Generic.MatSpace{fmpq}(FlintQQ, 2, 2, false));

julia> typeof(m)
AbstractAlgebra.Generic.MatSpaceElem{fmpq}

julia> typeof(-m)
fmpq_mat

julia> for mat in [2m, m+m, m*m, m^2, sub(m, 1:2, 1:2), m[1:2,:], divexact(m, 1),
                   kronecker_product(m, m), transpose(m), gram(m),
                   perm"(1, 2)" * m, solve_triu(m, m), solve(m, m),
                   pseudo_inv(m)[1], inv(m), nullspace(m)[2],
                   left_kernel(m)[2], kernel(m)[2], hcat(m, m), [m; m], map(identity, m)]
          @assert fmpq_mat == typeof(mat)
end

julia> typeof(m^1) == typeof(m^2)
false

julia> powers(m, 3)
ERROR: MethodError: Cannot `convert` an object of type fmpq_mat to an object of type AbstractAlgebra.Generic.MatSpaceElem{fmpq}
Closest candidates are:
[...]

Bug? not bug?

[wbhart]

Definitely a bug. The first case alone is clearly bad. If you start with something in a given matrix space, -m should not give you something in another matrix space.

We definitely need to fix this.

It's not at all clear why it is doing that. I'm pretty sure Nemo didn't used to do that, did it?

[wbhart]

It seems that the -m bug must be in similar. It looks like similar is not actually creating a similar matrix.

Let me look at the definition of similar and see if this can be fixed.

[thofma]

It started when we added the option to change the base ring in similar.

[wbhart]

The issue is the definition of similar in Nemo:

function similar(::MatElem, R::FlintIntegerRing, r::Int, c::Int)
   z = fmpz_mat(r, c)
   z.base_ring = R
   return z
end

You just need to change the first argument to have type fmpz_mat.

One clearly doesn't get a "similar" matrix in any other case.

[wbhart]

Or you can change the general definition of similar(M) in AbstractAlgebra to not call the version which "changes" the base ring.

[thofma]

But we still want similar(::fmpq_mat, ::FlintIntegerRing, ::Int, ::Int) to return fmpz_mat.

And similar for similar(::Generic.MatSpaceElem{nf_elem}, ::FlintIntegerRing, ::Int, ::Int).

[thofma]

There is no version which does not change the base ring, since it is a keyword argument. So it will always call similar(M, base_ring = FlintZZ) which will always produce fmpz_mat.

[thofma]

I see that you did not follow all the discussions about similar that we had in the last weeks :)

[wbhart]

But we still want similar(::fmpq_mat, ::FlintIntegerRing, ::Int, ::Int) to return fmpz_mat.

And similar for similar(::Generic.MatSpaceElem{nf_elem}, ::FlintIntegerRing, ::Int, ::Int).

What!? The first one I can live with. The second is surely a bug! I'll come and discuss with you.

[rfourquet]

since it is a keyword argument

Minor remark, it's not a keyword arg, just a normal arg with default value.

What!? The first one I can live with.

This is what I had implemented initially for this behavior: function similar(::NemoMatrices, R::FlintIntegerRing, r::Int, c::Int) ... end, i.e. similar would return a Nemo matrix when the input is a Nemo matrix. But Thommy wished the current version, in order to get this property:

similar(similar(M, R), base_ring(M)) == similar(M)

Which looks nice, but unfortunately leads to this questionable behavior of similar at the root of the bugs in this issue. See #638 for context.

[wbhart]

@thofma and I had a (very) long discussion about what behaviour we ultimately want.

In his case, he wants change_base_ring(ZZ, M) in Nemo, where M is a Generic.MatSpaceElem{nf_elem} to return an fmpz_mat. Reason: both fmpz_mat and Generic.MatSpaceElem{nf_elem} are matrices that Nemo would create.

In my case, I want change_base_ring(Nemo.ZZ, p) in Singular, where p is a Singular spoly{fmpq} to return an spoly{fmpz}. Reason: both polynomials are what Singular would create over a Nemo coefficient ring.

So we don't disagree on what we want, it's just not implemented that way at present (and probably wasn't implemented correctly before the new similar functionality, either).

The following simple 7 steps will lead to balance being restored in the force:

1. Remove change base ring functionality from similar, i.e. similar will only create a similar kind of matrix over the same base ring. This will fix the -m issue.

2. Create abstract types NemoMatrixElem{T} <: AbstractAlgebra.MatrixElem{T}, NemoPolyElem{T} <: AbstractAlgebra.PolyElem{T}, NemoMPolyElem{T} <: AbstractAlgebra.MPolyElem{T} and make all Nemo matrix and polynomial element types belong to one of these, instead of directly to the AbstractAlgebra abstract types.

2a. Create an abstract type NemoRing <: AbstractAlgebra.Ring and make all Nemo ring element types belong to this instead of directly to the AbstractAlgebra abstract type.

3. In AbstractAlgebra, implement generic fallback change_base_ring functions, which always creates generic Mat/Poly/MPoly's. It will have signature, e.g: change_base_ring(::Ring, ::AbstractAlgebra.MatElem).

4. For Nemo rings with native matrix, poly or mpoly types, e.g. fmpz, implement change_base_ring(::fmpz, ::NemoMat). This should always create a Nemo matrix (fmpz_mat, in the example).

5. For Nemo rings without native matrix, poly or mpoly types, e.g. nf_elem, implement change_base_ring(::NemoRing, Generic.MatSpaceElem{nf_elem}) which will always construct a Nemo matrix. This will ensure @thofma gets the functionality he wants.

6. map_entries/coeffs(f, m) should dispatch to map_entries/coeffs(R, f, M) where R = parent(zero(base_ring(M))).

7. Implement a similar strategy for map_entries/coeffs as for change_base_ring, based on NemoMat, NemoRing. Note that we never want to explicitly create MatrixSpace's in the code, whereas we don't mind creating parent objects for polynomials, so the matrix versions might be handled slightly differently (by calling matrix(...) to construct the matrix, presumably).

I don't discuss here how we handle the Singular functionality, which allows for creating Singular polynomials over Nemo coefficient rings and vice versa. But it is clear that the strategy above is consistent across all our systems and will work there too.

[wbhart]

I've opened a separate ticket for this. Nemocas/AbstractAlgebra.jl#490

[rfourquet]

Re-opening, for tracking progress on fixing the examples above.

[rfourquet]

The examples above are now fixed, with tests for some of them.