concretize methods

[vpuri3]

fix #196

@ChrisRackauckas what is the expected heuristic here? Currently this only checks if we can convert to AbstractMatrix. Not if that's efficient or preferable.

Does isconcrete(L) == true implies we can getindex/ setindex!? Do we want separate traits can_setindex, can_getindex?

EDIT

We can simplify the workflow to just two functions

• isconvertible indicating if L can be converted (cheaply) to an AbstractMatrix. This is true for lazy combinations of ScalarOp, MatrixOp, but not for TensorProdOp, InvertedOp. For FunctionOp, the behaviour can be set with kwarg isconvertible.
• concretize(L) returns a concrete type for L. It is a no-op when L is already concrete.

workflow in downstream packages:

if need_concrete_mat
    if !isconvertible(L)
        throw(ArgumentError("method needs concrete matrix, but $L cannot be converted to a concrete type." *
                            "Please initialize $L as an `AbstractMatrix`, or `Number`"))
    end
    L = concretize(L)
end

[codecov]

Codecov Report

Merging #205 (03d6d39) into master (303622a) will decrease coverage by 0.89%.
The diff coverage is 12.90%.

@@            Coverage Diff             @@
##           master     #205      +/-   ##
==========================================
- Coverage   71.12%   70.24%   -0.89%     
==========================================
  Files          10       10              
  Lines        1562     1583      +21     
==========================================
+ Hits         1111     1112       +1     
- Misses        451      471      +20     

Impacted Files	Coverage Δ
src/SciMLOperators.jl	100.00% <ø> (ø)
src/basic.jl	71.73% <0.00%> (-0.18%)	⬇️
src/batch.jl	63.01% <0.00%> (-4.64%)	⬇️
src/func.jl	83.94% <0.00%> (-0.62%)	⬇️
src/scalar.jl	54.24% <0.00%> (-0.36%)	⬇️
src/tensor.jl	78.62% <0.00%> (-0.28%)	⬇️
src/matrix.jl	67.52% <12.50%> (-2.53%)	⬇️
src/interface.jl	51.45% <23.07%> (-1.07%)	⬇️

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

[vpuri3]

I think a good set of traits would be

• can_convert indicating if a convert method is the same storage cost as the operator itself. true for lazy combinations of MatrixOp, ScalarOp but not tensorprod.
• isconcrete indicating if getindex/ setindex is possible (only true for MatrixOp)

both will be false for FunctionOp

[vpuri3]

@ChrisRackauckas does this match what you need for ode.jl?

[ChrisRackauckas]

Looks like it.

[vpuri3]

@ChrisRackauckas take a look?

[vpuri3]

@ChrisRackauckas can you take a look? The interface is described in the PR description

[vpuri3]

Question: are convert methods expected to return copies or the original object? i.e. Should convert(L::MatrixOperator) return L.A or copy(L.A)?

EDIT

Answer: it returns an aliased array

julia> x = ones(4,4);

julia> y = convert(AbstractArray, x);

julia> y[:] .= 0;

julia> x
4×4 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  0.0

[vpuri3]

@ChrisRackauckas ping

[ChrisRackauckas]

Sorry I missed this one so much, I'm not sure how. What the reason for the preference of concretize(L) instead of just convert(AbstractMatrix, L)? The convert method would be more general and would cover AbstractArray types.

[vpuri3]

SciMLOperators are like matrices, and ScalarOperators are like numbers. So concretize(L) is convenience for converting the former to AbstractMatrix, and the latter to Number. I don't see how it would be possible to convert to AbstractArray of any other shape.

[ChrisRackauckas]

I see what you mean. May be hard to use in practice unless we make `concretize(L::Union{Number,AbstractArray}) = L

[ChrisRackauckas]

Suggested change

	AbstractMatrix,
	Factorization,
	# SciMLOperators
	AbstractSciMLOperator,
	}
	) = convert(AbstractMatrix, L)
	AbstractArray,
	Factorization,
	# SciMLOperators
	AbstractSciMLOperator,
	}
	) = convert(AbstractArray, L)

[ChrisRackauckas]

Should that be extended to array?

[vpuri3]

No because, an ND array cannot be used as an operator, i.e. you can't mul! or * it to another array.

julia> ones(4,4,4) * ones(4)
ERROR: MethodError: no method matching *(::Array{Float64, 3}, ::Vector{Float64})

So we strictly want to return AbstractMatrix types

[vpuri3]

And, in the package, I've defined convert(::Type{AbstractMatrix}, L::MyOperator) methods for each operator type like

Base.convert(::Type{AbstractMatrix}, ii::IdentityOperator) = Diagonal(ones(Bool, ii.len))

So calling convert with AbstractArray would lead to errors:

julia> convert(AbstractArray, DiagonalOperator(ones(4)))                                                                                                                                                 [0/110]
ERROR: MethodError: Cannot `convert` an object of type 
  MatrixOperator{Float64, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}} to an object of type  
  AbstractArray

[ChrisRackauckas]

Well this works for now.

[gaurav-arya]

add to docs?

edit: #211

[ChrisRackauckas]

Thanks, missed that.