What is the idiomatic way to manipulate Matrix types and Matrix{Elem}?

[quachpas]

Hello! Thanks for this awesome package. I am manipulating matrices in and out, and I've hit a syntactic bump where it feels awkward to combine AbstractArray operations with Nemo matrix types, e.g.,

In [251]: A = QQ[0 1; 2 0]
Out[251]: [0   1]
[2   0]

In [252]: @which rank(A[1, :])
Out[252]: rank(x::Union{Number, AbstractVector})
     @ LinearAlgebra ~/.julia/juliaup/julia-1.11.0+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1003

In [253]: @which rank(A[2, :])

Out[253]: rank(x::Union{Number, AbstractVector})
     @ LinearAlgebra ~/.julia/juliaup/julia-1.11.0+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1003

In [254]: B = hcat(A[1, :], A[2, :])
Out[254]: 2×2 Matrix{QQFieldElem}:
 0  2
 1  0

In [255]: @info typeof(B)
[ Info: Matrix{QQFieldElem}

In [256]: rank(B)
ERROR: MethodError: no method matching AbstractFloat(::QQFieldElem)
The type `AbstractFloat` exists, but no method is defined for this combination of argument types when trying to construct it.

In [257]: rank(matrix(QQ, B))
Out[257]: 2

I end up using matrix(R, arr) a lot, not sure if this is the correct way to do things?

[joschmitt]

In principle, there is nothing wrong with using matrix(R, arr). Converting between Nemo matrices and Julia Matrix/Vector is usually inefficient though.

I am not sure what you are aiming to do in the example above. First of all, B ends up being the transpose of A, so you might want to do B = transpose(A) in this case. Further, accessing rows of A via A[1, :] gives you a Vector and not a matrix; this is consistent with what Julia does for Matrix. If you want to have the row as a (Nemo) matrix, you can do A[1:1, :] (so ranges in both indices).

[quachpas]

Thanks for the answer!

The examples above do not lead anywhere. The full code is given below, as a translation from MATLAB to Julia from [1].
Thanks for the reminder about accessing matrices rows/columns. Writing things this way, I managed to preserve most operations as a RMatrix. Although I was surprised by a few things:

Nemo

In [564]: A = QQ[0; 1]
Out[564]: [0]
[1]
In [565]: [A zeros(QQ,2)] 
ERROR: DimensionMismatch: mismatch in dimension 1 (expected 1 got 2)
In [566]: [A zeros(QQ, 2, 1)]
ERROR: DimensionMismatch: mismatch in dimension 1 (expected 1 got 2)
In [567]: [A zero_matrix(QQ, 2, 1)]
Out[567]: [0   0]
[1   0]
In [568]: hcat(A, I(2))
ERROR: DimensionMismatch: mismatch in dimension 1 (expected 1 got 2)
In [569]: vcat(A, I(2))
Out[569]: 3×2 Matrix{Any}:
      [0; 1]       [0; 1]
  true        false
 false         true

With Matrix(A), errors above are the opposite.

Compared to Julia

In [570]: B = [0; 1;;]
Out[570]: 2×1 Matrix{Int64}:
 0
 1
In [571]: [B zeros(2)]
Out[571]: 2×2 Matrix{Float64}:
 0.0  0.0
 1.0  0.0
In [572]: [B zeros(2, 1)]
Out[572]: 2×2 Matrix{Float64}:
 0.0  0.0
 1.0  0.0
In [573]: hcat(B, I(2))
Out[573]: 2×3 Matrix{Int64}:
 0  1  0
 1  0  1
In [574]: vcat(B, I(2))
ERROR: DimensionMismatch: number of columns of each array must match (got (1, 2))

I am unsure whether this is a bug, or a desired separation between RMatrix types and julia Matrix. It is certainly surprising to have to replace zeros with zero_matrix, I(d) with identity_matrix(R, d), etc. I feel like RMatrix should not be used with Julia Matrices if possible. It is treated as a completely separate object, hence why vcat(A, I(2)) would work.

My working code for the curious

function build_dictionary(Ā)
    Ā = matrix(QQ, Ā)
    @debug "Constructing initial dictionary without objective function"
    A0 = Ā[:, 1:(end - 1)]
    b = Ā[:, end:end]

    N = Int[0]
    (n, d) = size(A0)

    A_N = QQMatrix[]
    for i in 1:n
        if rank(A0[1:i, :]) > rank(A0[1:(i - 1), :])
            push!(A_N, A0[i:i, :])
            push!(N, i)
            @debug "A_N" A_N
            if rank(vcat(A_N...)) == d
                break
            end
        end
    end
    A_N = vcat(A_N...)

    B = setdiff(1:n, N[2:end])
    B = B[1:end]

    Ainv = inv(A_N)
    A1 = [b[B, :] zero_matrix(QQ, n - d, d)]
    A2 = Ainv * [b[N[2:end], :] -identity_matrix(QQ, d)]
    A3 = A0[B, :] * A2
    A = A1 - A3

    @debug "Initial dictionary constructed"
    return B, N, A
end

[1] https://github.com/Github-DongZelin/Examples-and-Command-of-the-algorithm/blob/main/functions_standard/Symbolic%20basis-enumeration%20function.mlx

[thofma]

Yes, what you observe is correct. The Nemo matrices and julia Array's do not play nice together, which is on purpose. There are some design decisions, which make Matrix useless for our applications, hence the separation.

In principle, we could make these block matrices work with Nemo matrices mixed with julia arrays, but you observed, these are two different worlds and so mixing might not be a good idea.

[quachpas]

I understand! I feel like it is a missed opportunity not to keep the same interface for both Nemo matrices and Julia arrays, but that would be an issue for AbstractAlgebra.jl probably.
I am closing this issue as completed then! Thank you both for your answers

[thofma]

Can you expand on the "missed interface"? Do you mean how zeros and Id work?

[fingolfin]

I feel like it is a missed opportunity not to keep the same interface for both Nemo matrices and Julia arrays, that would be an issue for AbstractAlgebra.jl probably.

You'll talk to the same people there. And you won't get a fundamentally different answer: if they could use the same interface, there'd be no reason for having our own kind of matrices.

But e.g. the Julia matrix infrastructure is heavily geared around the assumption that all you need to create a matrix is a type (for the elements) and the dimensions. But for our work, a type for entries is almost always insufficient: you need a parent ring for the elements (aka in our terminology as "base ring" of the matrix) instead.

There are other issues (see e.g. the FAQ entry at https://docs.oscar-system.org/dev/General/faq/)

That said, we generally try to reduce friction by imitating certain APIs where it is sensibly possible, and we are always open about hearing about specific use cases, and looking into supporting them. But it is always a tradeoff:

As @thofma already implied, we certainly could support e.g. [A zeros(QQ,2)] or hcat(A, I(2)) etc., although there is an inherent inefficiency to this. One can debate whether it is better to ease onboarding for people who are already familiar with Julia this way, as the cost of letting them learn and use "bad" patterns that are inherently inefficient.

But in the end it is also a question of "who does the work", i.e. someone has to implement this :-)

[quachpas]

You'll talk to the same people there.

I hadn't realised this!

Can you expand on the "missed interface"? Do you mean how zeros and Id work?

Sure! I meant that RMatrix standard functions have different names from Julia's matrix interface, i.e. zero_matrix vs zeros.

The current interface is grossly the following

Code	Type
zeros(T, n, m)	n x m Matrix{T}
zeros(R, n, m)	n x m Matrix{RElem}
zeros(RElem, n, m)	n x m Matrix{RElem}
zero_matrix(R,n ,m)	n x m RMatrix

This is valid for other function other than zeros. ATM I think there is a confusion between R and RElem in constructing Matrix{RElem}. For ease of use, zeros(R, r...) is defined, but it should be technically be disallowed, as R is not the Matrix element type in the output array, but RElem is. If zeros(R, n, m) would give a RMatrix instead, it would make "more sense". Switching things around would enable removing additional named functions such as zero_matrix and keep a "unified interface".

Code	Type	Function
zeros(T, n, m)	n x m Matrix{T}	Base.array.jl
zeros(RElem, n, m)	n x m Matrix{RElem}	Base.array.jl
zeros(R, n, m)	n x m Rmatrix	AbstractAlgebra.Matrix.jl

There are other issues (see e.g. the FAQ entry at https://docs.oscar-system.org/dev/General/faq/)

Thanks for the link, I'll take a look.

But in the end it is also a question of "who does the work", i.e. someone has to implement this :-)

The bane of all FOSS, actually finding time to implement!! If you're willing, I can open another issue for the above suggestion, or re-open this one. However, as you said, the tradeoff lies in the user onboarding and making a huge breaking change.