fill! behavior for structurally-constrained storage types?

[Sacha0]

Ref. discussion between @martinholters and @tkelman starting at #16740 (comment). How should fill! behave for storage types with structural constraints such as Bidiagonal and SparseMatrixCSC? Brief illustration:

julia> fill!(Diagonal(rand(4)), 1) # fills within structural constraints
4×4 Diagonal{Float64}:
 1.0   ⋅    ⋅    ⋅
  ⋅   1.0   ⋅    ⋅
  ⋅    ⋅   1.0   ⋅
  ⋅    ⋅    ⋅   1.0

julia> fill!(Bidiagonal(rand(4), rand(3), true), 1) # fails noting structural constraints
ERROR: ArgumentError: Array A of type Bidiagonal{Float64} and size (4,4) can
    not be filled with x=1, since some of its entries are constrained.
 in fill!(::Bidiagonal{Float64}, ::Int64) at ./linalg/bidiag.jl:553

julia> spfoo = sprand(4, 4, 0.5)
4×4 sparse matrix with 7 Float64 nonzero entries:
    [2, 1]  =  0.403982
    [3, 1]  =  0.521604
    [4, 1]  =  0.857491
    [1, 2]  =  0.318858
    [3, 3]  =  0.893084
    [2, 4]  =  0.646681
    [4, 4]  =  0.601919

julia> fill!(spfoo, 1) # fills without respecting structural constraints
4×4 sparse matrix with 16 Float64 nonzero entries:
    [1, 1]  =  1.0
    [2, 1]  =  1.0
    [3, 1]  =  1.0
    [4, 1]  =  1.0
    [1, 2]  =  1.0
    [2, 2]  =  1.0
    [3, 2]  =  1.0
    [4, 2]  =  1.0
    [1, 3]  =  1.0
    [2, 3]  =  1.0
    [3, 3]  =  1.0
    [4, 3]  =  1.0
    [1, 4]  =  1.0
    [2, 4]  =  1.0
    [3, 4]  =  1.0
    [4, 4]  =  1.0

In other words, does fill! mean 'fill every addressable entry in the provided AbstractArray', or 'fill every stored entry in the provided AbstractArray'? Should fill! bifurcate into, e.g., fillall! and fillstored!? Best!

[martinholters]

Bifurcation sounds good, guess I'd prefer fill and fillnz.

[tkelman]

There is a fillslots! helper function but it's not exported yet I believe. I think fillstored! would be a somewhat better name.

[Sacha0]

fillnz! seems ambiguous: What should fillnz!(Diagonal([0, 0, 0]), 1) do? Yield Diagonal([0, 0, 0]) or Diagonal([1, 1, 1])? Best!

[JaredCrean2]

I would reverse the convention: fill! = fillstored! a separate fillall!. I think that leads to easier dense-sparse generalizations (at least for the class of problems I usually deal with). Densifying a sparse matrix should be opt-in.

[tkelman]

fill!(A, 1) not meaning the same thing as A[:] = 1 would be really problematic and error-prone from a generic programming perspective. Getting storage-dependent behavior is the non-generic case and should not be the standard API.

[Sacha0]

I would reverse the convention: fill! = fillstored! a separate fillall!. I think that leads to easier dense-sparse generalizations (at least for the class of problems I usually deal with). Densifying a sparse matrix should be opt-in.

On the one hand I like this argument. Being explicit (fillstored!/fillall!) and using fillstored! for existing uses of fill! strikes me as worth typing the extra few characters though.

On the other hand I also like @tkelman's argument. Which makes the explicit names seem still more attractive, requiring the author to make clear her intention.

[JeffBezanson]

A think a good way to address this is to use representations of the relevant indices. For example, instead of having Xstored! for many functions X, do something like A[stored(A)] = 1. There has recently been some related discussion here: JuliaData/IndexedTables.jl#23
Also in https://github.com/timholy/ArrayIteration.jl @timholy has been experimenting with these kinds of abstractions.

[JaredCrean2]

For one-liners like fill! the stored(A) approach looks good, but for more involved operations the question of how to define the Xstored, Xall, and X methods is still relevant.

@tkelman yeah, breaking that identity would be bad, although I'm not sure using colons on sparse matrices is a good idea at all. I guess the next best thing would be fill! = fillall! for sparse matrices and make fillstored! fall back to fill! for dense matrices.

As a general principle, I think the sparse-specific functions should have fallbacks for the dense case to enable generic programming.

[tkelman]

I'm not sure using colons on sparse matrices is a good idea at all.

It isn't, but unless it's going to be a method error I think we should keep the existing "linear algebraic" concepts consistent between dense and sparse, as most of the AbstractArray interface tends to be designed around that. New abstractions for storage-dependent behavior should probably be the special cases, and need careful thought to migrate code to use them correctly - somewhat like (otherwise dense) unconventionally-indexed arrays.

[Sacha0]

Per triage, the long-term direction should be fill!(stored(A), val). For the short term, perhaps fillslots! -> fillstored! and remains unexported from LinAlg. Best!

[StefanKarpinski]

+1 fillstored! and eventually fill!(stored(A), val).

[ChrisRackauckas]

Bumping on this. I just hit it, and from the looks of it it seems other scientific computing groups are hitting it too. It's very non-intuitive for fill! to densify.