Error message for non-sparse Jacobian pattern is not obvious

[DanielVandH]

If I provided a non-sparse sparsity pattern for the Jacobian, the warning is not so obvious that it's about not having a sparsity pattern:

using FiniteDiff
function f(dx, x)
    dx .= [x[1]^2 + x[2]^2]
end
J = zeros(1, 2)
θ = [5.0, 3.0]
y0 = [0.0]
cache = FiniteDiff.JacobianCache(copy(θ), copy(y0), copy(y0), Val(:forward); sparsity = [1 1])
FiniteDiff.finite_difference_jacobian!(J, f, θ, cache)

ERROR: MethodError: no method matching length(::Nothing)
Closest candidates are:
  length(::Union{Base.KeySet, Base.ValueIterator}) at abstractdict.jl:58
  length(::Union{LinearAlgebra.Adjoint{T, <:Union{StaticArraysCore.StaticArray{Tuple{var"#s2"}, T, 1} where var"#s2", StaticArraysCore.StaticArray{Tuple{var"#s3", var"#s4"}, T, 2} where {var"#s3", var"#s4"}}}, LinearAlgebra.Diagonal{T, <:StaticArraysCore.StaticArray{Tuple{var"#s13"}, T, 1} where var"#s13"}, LinearAlgebra.Hermitian{T, <:StaticArraysCore.StaticArray{Tuple{var"#s10", var"#s11"}, T, 2} where {var"#s10", var"#s11"}}, LinearAlgebra.LowerTriangular{T, <:StaticArraysCore.StaticArray{Tuple{var"#s18", var"#s19"}, T, 2} where {var"#s18", var"#s19"}}, LinearAlgebra.Symmetric{T, <:StaticArraysCore.StaticArray{Tuple{var"#s7", var"#s8"}, T, 2} where {var"#s7", var"#s8"}}, LinearAlgebra.Transpose{T, <:Union{StaticArraysCore.StaticArray{Tuple{var"#s2"}, T, 1} where var"#s2", StaticArraysCore.StaticArray{Tuple{var"#s3", var"#s4"}, T, 2} where {var"#s3", var"#s4"}}}, LinearAlgebra.UnitLowerTriangular{T, <:StaticArraysCore.StaticArray{Tuple{var"#s24", var"#s25"}, T, 2} where {var"#s24", var"#s25"}}, LinearAlgebra.UnitUpperTriangular{T, <:StaticArraysCore.StaticArray{Tuple{var"#s21", var"#s22"}, T, 2} where {var"#s21", var"#s22"}}, LinearAlgebra.UpperTriangular{T, <:StaticArraysCore.StaticArray{Tuple{var"#s15", var"#s16"}, T, 2} where {var"#s15", var"#s16"}}, StaticArraysCore.StaticArray{<:Tuple, T}, StaticArraysCore.StaticArray{Tuple{var"#s25"}, T, 1} where var"#s25", StaticArraysCore.StaticArray{Tuple{var"#s1", var"#s3"}, T, 2} where {var"#s1", var"#s3"}} where T) at C:\Users\licer\.julia\packages\StaticArrays\5bAMi\src\abstractarray.jl:1    
  length(::Union{LinearAlgebra.Adjoint{T, S}, LinearAlgebra.Transpose{T, S}} where {T, S}) at C:\Users\licer\AppData\Local\Programs\Julia-1.8.0-rc1\share\julia\stdlib\v1.8\LinearAlgebra\src\adjtrans.jl:172
  ...
Stacktrace:
 [1] _colorediteration!
   @ C:\Users\licer\.julia\packages\FiniteDiff\dhnni\src\iteration_utils.jl:2 [inlined]
 [2] finite_difference_jacobian!(J::Matrix{Float64}, f::typeof(f), x::Vector{Float64}, cache::FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Matrix{Int64}, Val{:forward}(), Float64}, f_in::Nothing; relstep::Float64, absstep::Float64, colorvec::UnitRange{Int64}, sparsity::Matrix{Int64}, dir::Bool)
   @ FiniteDiff C:\Users\licer\.julia\packages\FiniteDiff\dhnni\src\jacobians.jl:393
 [3] finite_difference_jacobian!(J::Matrix{Float64}, f::Function, x::Vector{Float64}, cache::FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Matrix{Int64}, Val{:forward}(), Float64}, f_in::Nothing) 
(repeats 2 times)
   @ FiniteDiff C:\Users\licer\.julia\packages\FiniteDiff\dhnni\src\jacobians.jl:324
 [4] top-level scope
   @ c:\Users\licer\using FiniteDiff.jl:9

using SparseArrays
cache = FiniteDiff.JacobianCache(copy(θ), copy(y0), copy(y0), Val(:forward); sparsity = sparse([1 1]))
FiniteDiff.finite_difference_jacobian!(J, f, θ, cache) # works!

It could be nice to throw an error in this case, e.g. defining something like (code below isn't tested)

const INVALID_PROTOTYPE_ERROR = """
                                A sparsity pattern was provided that was a matrix but was not sparse. Consider calling `SparseArrays.sparse` on your pattern and re-defining the cache.
                                """
struct InvalidPrototypeError <: Exception end 

function Base.showerror(io::IO, ::InvalidPrototypeError)
    print(io, INVALID_PROTOTYPE_ERROR)
end

mutable struct JacobianCache{CacheType1,CacheType2,CacheType3,CacheType4,ColorType,SparsityType,fdtype,returntype}
    x1  :: CacheType1
    x2  :: CacheType2
    fx  :: CacheType3
    fx1 :: CacheType4
    colorvec :: ColorType
    sparsity :: SparsityType
    JacobianCache(x1, x2, fx, fx1, colorvec, sparsity::AbstractMatrix) = !(sparsity isa AbstractSparseMatrix) ? throw(INVALID_PROTOTYPE_ERROR) : new(x1, x2, fx, fx1, colorvec, sparsity)
end

(I think that's how new would be used? Haven't used it really.)

Worth a pull request? Maybe the error message isn't really about the sparsity pattern at all since jacobians.jl actually has a case for J not being an AbstractSparseMatrix, in which case the proposal above is not needed but this is still worth the issue I guess. I'm just assuming it should be sparse since the README claims "sparsity should be a sparse or structured matrix (Tridiagonal, Banded, etc. according to the ArrayInterfaceCore.jl specs)".

[DanielVandH]

The issue, if you do want to allow for non-sparse patterns with no problem, probably comes from this part of jacobians.jl

# J is a sparse matrix, so decompress on the fly
                @. vfx1 = (vfx1 - vfx) / epsilon
                if ArrayInterfaceCore.fast_scalar_indexing(x1)
                    if sparseCSC_common_sparsity
                        _colorediteration!(J,vfx1,colorvec,color_i,n)
                    else
                        _colorediteration!(J,sparsity,rows_index,cols_index,vfx1,colorvec,color_i,n)
                    end
                else
                    #=
                    J.nzval[rows_index] .+= (colorvec[cols_index] .== color_i) .* vfx1[rows_index]
                    or
                    J[rows_index, cols_index] .+= (colorvec[cols_index] .== color_i) .* vfx1[rows_index]
                    += means requires a zero'd out start
                    =#
                    if J isa AbstractSparseMatrix
                        @. void_setindex!((J.nzval,), getindex((J.nzval,), rows_index) + (getindex((_color,), cols_index) == color_i) * getindex((vfx1,), rows_index), rows_index)
                    else
                        @. void_setindex!((J,), getindex((J,), rows_index, cols_index) + (getindex((_color,), cols_index) == color_i) * getindex((vfx1,), rows_index), rows_index, cols_index)
                    end
                end

which says that J is sparse, and then also checks that J is an abstract sparse matrix anyway - maybe some of the ideas were confused at the time? The problematic part of this is _coloriteration, which calls this method:

@inline function _colorediteration!(J,sparsity,rows_index,cols_index,vfx,colorvec,color_i,ncols)
    @inbounds for i in 1:length(cols_index)
        if colorvec[cols_index[i]] == color_i
            J[rows_index[i],cols_index[i]] = vfx[rows_index[i]]
        end
    end
end

But rows_index and cols_index are both nothing in this case (the sparsity argument isn't even used actually - nor is it in the other _colorediteration! for sparse matrices?). It seems like there just isn't intended to be any support for non-sparse matrices, despite in some places in the code looking like that was maybe intended initially.

[ChrisRackauckas]

Yeah it was kind of intended, but then we realized everything like Diagonal, BlockBandedMatrix, etc. all went the sparse route. Your fix is good.