Writing a new subtype of AbstractArray

New to Julia
1

I’ve just defined a new type, which represents a lower-triangular matrix and is a subtype of AbstractArray. I extended all the methods required by that interface. When I create an object of this type in the REPL, Julia tries to print a representation off it. This results in a bounds error, since the default display function tries to print a full square matrix, and mine is only lower triangular.

I tried to fix this problem by redefining Base.show. That function works when called directly, but it doesn’t solve the problem. Julia seems not to call “show” when displaying a subtype of AbstractArray. Any suggestions?

My code is below:

using StaticArrays

"""
Lower triangular square matrix stored in column-major order.

    x        offset of col 1 = 0
    x x      offset of col 2 = 4
    x x x    offset of col 3 = 7
    x x x x  offset of col 4 = 9

    A[i, j] = data[i-j+1 + offset[j]]
    A[4, 4] = data[1 + 9] = data[10]
"""
mutable struct LTriMatrix{N,M,T} <: AbstractMatrix{T}
    dim :: Int
    data :: MVector{M,T}
    offset :: SVector{N,Int}

    function LTriMatrix(dim, x)
        n = Int(dim)
        m = (n*(n+1)) ÷ 2
        mat = new{n, m, typeof(x)}()
        mat.dim = n

        mat.data = MVector{m, typeof(x)}(undef)
    
        offset = MVector{n,Int}(undef)
        offset[1] = 0
        for i in 2:n
            offset[i] = offset[i-1] + n - i + 2
        end
        mat.offset = SVector(offset)
        mat
    end
end

Base.size(A::LTriMatrix{N,M,T}) where {N,M,T} = (A.dim, A.dim)

Base.length(A::LTriMatrix{N,M,T}) where {N,M,T} = length(A.data)

Base.IndexStyle(::LTriMatrix{N,M,T}) where {N,M,T} = IndexCartesian()

function Base.getindex(A::LTriMatrix{N,M,T}, i::Int, j::Int) where {N,M,T}
    if(i < j)
        error("column index exceeds row index")
    end
    A.data[i-j+1 + A.offset[j]]
end

# Fails if value cannot be converted to type T.
function Base.setindex!(A::LTriMatrix{N,M,T}, value, i::Int,
                        j::Int) where {N,M,T}
    if(i < j)
        error("column index exceeds row index")
    end
    A.data[i-j+1 + A.offset[j]] = value
end

"""
Write a representation of an LTriMatrix to io. This works when called
explicitly, but is not called by default when Julia tries to display
an object of this type on the terminal.
"""
function Base.show(io::IO, A::LTriMatrix{N,M,T}) where {N,M,T}
    for i in 1:A.dim
        for j in 1:i
            print(io, j>1 ? " " : "", A[i,j])
        end
        println(io)
    end
end
2

To be an AbstractArray, the type must be rectangular in its shape. Effectively all methods that work on AbstractArrays demand this property. You can still build a LTriMatrix, but just don’t make it a subtype of AbstractMatrix. Defining it as a subtype like this gains you nothing but bugs (seriously!).

For your root question, AbstractArray defines its show methods with a particular MIME type: show(io::IO, ::MIME{Symbol("text/plain")}, X::AbstractArray), but again, this is just the first symptom of a much larger problem lurking behind it.

1 Like
3

Thanks for the help.

4

Usually “lower-triangular” is understood to mean that there are zeros above the diagonal. If that’s what you want, and you simply return 0 instead of error in this case, it should work fine. But you can also simply use the existing LowerTriangular from LinearAlgebra.

4 Likes
5

Based on my reading of the docs, the version of LowerTriangular that LinearAlgebra provides is a view of a full rectangular matrix. I wanted to avoid storing the whole matrix.

I agree with your other point. It would be better to return 0 than to throw an error. Then I could subtype AbstractArray.

Thanks.

1 Like