add option for writable off-diagonal elements

[bjarthur]

@mkitti this is the first time i've used value types. can you please review?

the idea is that LinearAlgebra.Symmetric does not allow the off diagonal elements to be set, and so neither should the default PackedArrays.SymmetricPacked. but some might want this so add it as an option.

[mkitti]

I've been thinking of the reason why the type in Base is read-only. Thinking about it, it is problematic to use indexing to change the matrix. If you iterated over the indices, you might end up changing things twice.

Perhaps a solution then is to create a view of the matrix that is just the upper or lower triangle. Now it is perfectly valid to iterate over this triangle and mutate it.

[bjarthur]

not being able to set off diagonal elements is frequently raised as an issue. see e.g. JuliaLang/julia#43818 from just 5 days ago. your logic to not permit this is the reason it was decided against. see JuliaLang/julia#33071 (comment)

what i'm proposing here is to allow the user flexibility in choosing to permit it or not. there are use cases where it is desirable i imagine. and not being able to do so just seems strange to me.

if you wanted to iterated over just the triangle, you can simply reach into the struct and iterate over the field 'tri'. if this were commonly desired, creating an interface to do that would be best.

most importantly, what do you think of using a value type to parameterize whether the off diagonal elements are writeable? the alternative is to create two types.

[mkitti]

Yes, but tri is just a vector. What I want is a type that knows it is a triangle and that has two-dimensional indices reflecting where in the triangle we are. Trying to access two-dimensional indices on the other triangle of the matrix is a BoundsError.

Simply adding a :RW interface does not seem to actually address the issue. What is missing is the triangle abstraction I just described. A symmetric matrix itself does not satisfy the issue because it is trying to look like a matrix with full 2D indices.

I think the read-write interface should be a TriangleView. It could wrap around a symmetric matrix to provide native read-write access to the underlying vector. Perhaps it could also provide a way of mapping the other triangle to a native triangle.

[bjarthur]

you know about https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.UpperTriangular ?

[mkitti]

Yes, but that describes an upper triangular matrix. I just want a triangle. It should not represent a matrix.

julia> ut = UpperTriangular(zeros(5,5))
5×5 UpperTriangular{Float64, Matrix{Float64}}:
 0.0  0.0  0.0  0.0  0.0
  ⋅   0.0  0.0  0.0  0.0
  ⋅    ⋅   0.0  0.0  0.0
  ⋅    ⋅    ⋅   0.0  0.0
  ⋅    ⋅    ⋅    ⋅   0.0

julia> CartesianIndices(ut)
5×5 CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}:
 CartesianIndex(1, 1)  CartesianIndex(1, 2)  CartesianIndex(1, 3)  CartesianIndex(1, 4)  CartesianIndex(1, 5)
 CartesianIndex(2, 1)  CartesianIndex(2, 2)  CartesianIndex(2, 3)  CartesianIndex(2, 4)  CartesianIndex(2, 5)
 CartesianIndex(3, 1)  CartesianIndex(3, 2)  CartesianIndex(3, 3)  CartesianIndex(3, 4)  CartesianIndex(3, 5)
 CartesianIndex(4, 1)  CartesianIndex(4, 2)  CartesianIndex(4, 3)  CartesianIndex(4, 4)  CartesianIndex(4, 5)
 CartesianIndex(5, 1)  CartesianIndex(5, 2)  CartesianIndex(5, 3)  CartesianIndex(5, 4)  CartesianIndex(5, 5)

julia> ut isa AbstractMatrix
true

[mkitti]

Is there a lazy iterator for this sequence?

julia> f(n,m) = [CartesianIndex(i,j) for i in 1:n, j in 1:m if i <= j]
f (generic function with 1 method)

julia> f(1,1)
1-element Vector{CartesianIndex{2}}:
 CartesianIndex(1, 1)

julia> f(2,2)
3-element Vector{CartesianIndex{2}}:
 CartesianIndex(1, 1)
 CartesianIndex(1, 2)
 CartesianIndex(2, 2)

julia> f(3,3)
6-element Vector{CartesianIndex{2}}:
 CartesianIndex(1, 1)
 CartesianIndex(1, 2)
 CartesianIndex(2, 2)
 CartesianIndex(1, 3)
 CartesianIndex(2, 3)
 CartesianIndex(3, 3)

julia> f(4,4)
10-element Vector{CartesianIndex{2}}:
 CartesianIndex(1, 1)
 CartesianIndex(1, 2)
 CartesianIndex(2, 2)
 CartesianIndex(1, 3)
 CartesianIndex(2, 3)
 CartesianIndex(3, 3)
 CartesianIndex(1, 4)
 CartesianIndex(2, 4)
 CartesianIndex(3, 4)
 CartesianIndex(4, 4)

julia> f(5,5)
15-element Vector{CartesianIndex{2}}:
 CartesianIndex(1, 1)
 CartesianIndex(1, 2)
 CartesianIndex(2, 2)
 CartesianIndex(1, 3)
 CartesianIndex(2, 3)
 CartesianIndex(3, 3)
 CartesianIndex(1, 4)
 CartesianIndex(2, 4)
 CartesianIndex(3, 4)
 CartesianIndex(4, 4)
 CartesianIndex(1, 5)
 CartesianIndex(2, 5)
 CartesianIndex(3, 5)
 CartesianIndex(4, 5)
 CartesianIndex(5, 5)

julia> f(3,2)
3-element Vector{CartesianIndex{2}}:
 CartesianIndex(1, 1)
 CartesianIndex(1, 2)
 CartesianIndex(2, 2)

julia> f(2,3)
5-element Vector{CartesianIndex{2}}:
 CartesianIndex(1, 1)
 CartesianIndex(1, 2)
 CartesianIndex(2, 2)
 CartesianIndex(1, 3)
 CartesianIndex(2, 3)

[mkitti]

@bjarthur Do you have any thoughts about changing iteration when you turn on RW mode? If RW is true, then my thought is that eachindex should only return indices for either the lower or upper triangle.

[mkitti]

For a lazy iterator, here's what have so far:

julia> struct TriangleIndices{N}
           dims::Dims{N}
       end

julia> Base.iterate(TI::TriangleIndices{N}) where N = (t = CartesianIndex(ntuple(i->1, N)); (t,t))

julia> function Base.iterate(TI::TriangleIndices{N}, state) where N
           out = nothing
           for d in 1:N
               n = state[d] + 1
               if n <= TI.dims[d] && (d == N) || n <= state[N]
                   out = ntuple(i-> i  < d ? 1 :
                                    i == d ? n :
                                    state[i],
                                N)
                   break
               end
               if d == N
                   return nothing
               end
           end
           out = CartesianIndex(out)
           return (out, out)
       end

julia> ti = TriangleIndices(4,4)
TriangleIndices{2}((4, 4))

julia> for i in ti
           println(i)
       end
CartesianIndex(1, 1)
CartesianIndex(1, 2)
CartesianIndex(2, 2)
CartesianIndex(1, 3)
CartesianIndex(2, 3)
CartesianIndex(3, 3)
CartesianIndex(1, 4)
CartesianIndex(2, 4)
CartesianIndex(3, 4)
CartesianIndex(4, 4)