Sparse matrix/vectors with fixed sparsity pattern.

[SobhanMP]

This pull request adds

• an abstract type like AbstractMatrixCSC for sparse vectors (AbstractSparseVectorC)
• an abstract type for AbstractMatrixCSC whose colptr and rowval cannot be changed
• an abstract type for AbstractSparseVectorC whose colptr and rowval cannot be changed
• broadcast is the union of pattern in the presence of fixed elements
• through the magic of constprop return _is_fixed(A, Bs...) ? fixed(r) : r is typesafe.
• the tests don't check for every thing but check that the operations are correct and do not allocate.
• the fixed function can be used as a shortcut for creating sparse things with fixed sparsity.

once all green this should be ready to merge. this commit is NOT breaking.

apart from the bike-shed possible, i have a few concerns

• should the nonzeros field be non-resizable? right now, since rowvals is non-resizable, it is technically non-resizable (in correct arrays).
• i broadened the definitions in the higher order fns file. Mainly i changed the definition from SparseMatrixCSC to AbstractSparseMatrixCSC. is this correct?
• when this get merged, i think it would be nice to have an operator on non zero values, something nicer than

x = sprandn(100, 100, 0.1)
y = typeof(x)(sizeof(x)..., copy(getcolptr(x)), copy(rowvals(x)), nonzeros(x) .+ 2)

• should we share the colptr between the instances of FixedSparseCSC

i too feel that it's unfortunately long, 🙃. i thought it would be easier than this.

Ref: #190

[SobhanMP]

after some hacking, i've come with some cases that needs discussing.
let $x$ be some CSC matrix with fixed pattern, what would x .+ 1 be? should it have no zeros or should it just add 1 to the nonzero entries? if y is another matrix like x, what should x .+ y look like? should it error if the structure of one does not fit in the other?

[rayegun]

X .+ Y should be a union of the patterns I believe. X .+ 1 is tricky, and I'll defer to @ChrisRackauckas here since he asked for this initially.

[ChrisRackauckas]

X .+ Y should be a union of the patterns I believe.

Yes.

X .+ 1 is tricky

That's just dense: there's no other representation. So make it a dense sparse matrix, but of course, almost certainly the user would never do that.

[SobhanMP]

for reasons i don't fully understand this seems to work.
broadcast now takes the union of the structure of the FixedSparseCSCs

I propose adding

fill!(x::FixedSparseCSC, v) = fill!(nonzeros(x), v)

as fill!(::FixedSparseCSC, v::Any) only works if iszero(v) or x is already a dense sparse matrix(sparse matrix with all entries). In either case, this is equivalent and only changes the behavior from error to something weird in some cases.

furthermore, i think it would be nice to have a shortcut for

FixedSparseCSC(size(x)..., copy(getcolptr(x)), copy(rowvals(x)), nonzeros(x) .+ 1)

i.e. add 1 to all non structural zeros. Since

x .+ 1

can't be that (as it would not work with the array api), maybe there could be another operator? maybe a macro could do the job?

[codecov-commenter]

Codecov Report

Merging #201 (bbdc86f) into main (7147438) will decrease coverage by 0.44%.
The diff coverage is 86.06%.

@@            Coverage Diff             @@
##             main     #201      +/-   ##
==========================================
- Coverage   92.13%   91.68%   -0.45%     
==========================================
  Files          11       12       +1     
  Lines        7210     7302      +92     
==========================================
+ Hits         6643     6695      +52     
- Misses        567      607      +40     

Impacted Files	Coverage Δ
src/solvers/cholmod.jl	90.66% <ø> (ø)
src/sparsematrix.jl	94.27% <67.39%> (-1.05%)	⬇️
src/solvers/umfpack.jl	88.10% <71.42%> (ø)
src/readonly.jl	73.68% <73.68%> (ø)
src/sparsevector.jl	95.11% <88.88%> (-0.49%)	⬇️
src/SparseArrays.jl	100.00% <100.00%> (ø)
src/abstractsparse.jl	47.36% <100.00%> (+2.92%)	⬆️
src/higherorderfns.jl	96.84% <100.00%> (+0.14%)	⬆️
src/solvers/spqr.jl	88.53% <100.00%> (+0.14%)	⬆️
... and 1 more

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[SobhanMP]

@ChrisRackauckas could you test it a bit? the instruction to use a dev version of this package is in the README, and you can use fixed to convert to fixed patterns.

[ViralBShah]

There's a lot going on in this PR - I'll need time to understand.

[rayegun]

Is AbstractSparseVectorC any different from AbstractSparseVector? I'm not sure it makes sense to differentiate between storage orders for a 1-d array.

I think storage order should be a trait anyway but that's a bigger issue.

[SobhanMP]

i assumed one can implement a AbstractSparseVector with a dictionary or start from 0 instead of 1. Hence AbstractSparseVectorC is the subset that work like SparseVector, with one indexed nonzeros/nonzeroinds. The point is that AbstractSparseArray{Tv, Ti, 1} does not say anything about the structure of the array.

[SobhanMP]

case in point, look at #219

[ViralBShah]

I am wondering if this code can live in a separate package? Or rather what do we need to do in SparseArrays.jl so that it can be a separate package?

[SobhanMP]

i don't think it can be a separate package w/o much redundancy (and piracy) in higherorderfns. having fixed arrays changes the meaning of broadcast/map/map!. also it should be much easier to check for breakage if they are part of this package.

[rayegun]

A few comments on the actual code inline.

You mentioned that this can't be split out because that would be piracy. I don't think there would be any need for piracy (these changes involve a new type anyway).

I would argue this should look like:

1. Make higherorderfns use AbstractSparseMatrixCSC rather than SparseMatrixCSC in the SparseOrStructuredMatrix union. Most other instances of SparseMatrixCSC in signatures can be turned into AbstractSparseMatrixCSC.
2. Add trait to determine whether new indices can be inserted. _is_fixed is fine as is, _can_insert might be more clear.
3. _noshapecheck_map can be overloaded in a new package for AbstractFixed..., perhaps as simply as _noshapecheck_map(f, A::AbstractFixed, Bs...) = fixed(invoke(_noshapecheck_map, AbstractSparseMatrixCSC, ..., A, Bs...).

That avoids any piracy (as far as I can tell), and lets us avoid significant additions to the SparseArrays API before Julia 1.10 when we can really look at major changes (SparseArrays will still be an stdlib for 1.9 unfortunately).

[rayegun]

Is this ready @SobhanMP ? Everything including concrete types can be merged once green. Please be sure to document that this is not public right now.

[SobhanMP]

i'd like to have a bit more time to check stuff, will merge by monday.