Generic abstract array construction without instances

[ChrisRackauckas]

We do not have a good way of fully specifying abstract arrays without giving the full instance. For example, in DiffEq I would like the user to be able to pass whatever kind of Jacobian matrix type suits their needs. This would allow them to specify a banded matrix from BandedMatrices.jl, or a sparse matrix, etc. However, this is going to be part of the high level problem type which is something that I assume holds little memory and thus is simple to copy around (since it is in every other case). I would like to not hold instances of giant sparse matrices here.

The problem is that there is no way to re-construct that sparse matrix generically only from type information. If I know the instance, I can do:

similar(A)

to get a new cache of it, or I can do

similar(A,T)

to get a new cache with a different element type (which can be necessary for handling units). But, if I just have the type, I don't know how to re-create their sprasity pattern, and thus I am stuck with the problem type holding instances of these large objects in a very memory inefficient manner.

Now, the way this is supposed to work is that someone in this case uses the sparse constructor along with some array which specifies the non-zero elements. And for dense arrays one can pass a tuple for size. And for BandedMatrix ...? You can see that handling this using the given methods is not generic: I would have to setup each cache construction to handle every specific case the user can throw at me.

But it doesn't need to be like this. If instead we had an idea of "the abstract information which is required to re-build this type", then we could specify the way to construct any abstract type without providing instances. The simplest case would be dense arrays, where

struct DenseArrayConstructorInformation{T,N} <: AbstractArrayConstructorInformation
  size::NTuple{T,N}
end

and thus (I am not wedded to these names at all)

construct(DenseArrayConstructorInformation((2,3),Float64)

is sufficient. But for sparse arrays, this could be different:

struct SparseArrayConstructorInformation{T,N,A} <: AbstractArrayConstructorInformation
  size::NTuple{T,N}
  nonzeros::A
end

where nonzeros can be any sufficiently nice collection (so this could be specified with a lazy array if this pattern is simple enough, making it take almost zero memory) and thus construct can be used on this. But then users can extend this. For example, the BandedMatrices.jl package could define the constructor information

struct BandedMatrioxConstructorInformation{T,N,A} <: AbstractArrayConstructorInformation
  size::NTuple{T,N}
  upper_size::Int
  lower_size::Int
end

And then this can be used with construct. Each AbstractArray can have a much more memory effect "layout" information object from which it can be constructed, and thus this makes it easier to pass around ideas about AbstractArrays without passing and re-creating instances themselves.

[Sacha0]

Linking #18161. Best!

[dlfivefifty]

I'm not convinced a special type is necessary here: I'm pretty happy with

BandedMatrix(uninitialized, (n,m), (l,u))

which could be generalized with a simple routine arrayparameters:

arrayparameters(A::Array) = size(A) # or maybe (size(A),)
arrayparameters(A::BandedMatrix) = (size(A), bandwidths(A))

then general code would work via

function foo(A:: AbstractArray)
   # make empty array of type typeof(A) with correct parameters:
   ret = typeof(A)(uninitialized, arrayparameters(A)...)
   # ... do stuff
end

[ChrisRackauckas]

which could be generalized with a simple routine arrayparameters

Yes, arrayparameters is pretty similar to what I am suggesting, except it requires an instance of the type and splatting in order to work. I don't see why you'd want to build the array in order to get the generic information to build another one like it (since that means you'd need to be passing around arrays instead of a generalized idea of "size" into cache construction routines), and it won't work well for a heterogeneously sized array of arrays. At some point, isn't it easier to just make something a type instead of having an API deal with tuples of tuples of ... and document every combination? So arrayparameters is the same in spirit but can be a lot worse performance-wise and IMO is much harder to document or extend.

[dlfivefifty]

I don't think the splat is an issue: just wrap the parameters in a tuple and the splat is fine:

arrayparameters(A::ArrayOfArrays{<:AbstractArray}) = (arrayparameters.(A),)
function ArrayOfArrays{T}(::Uninitialized, pars) where T <: AbstractArray 
   ret = Array{T}(uninitialized, size(pars))
   for (k,j) in eachindex(ret)
        ret[k,j] = T(uninitialized, pars[k,j]...)
   end
   ArrayOfArrays(ret)
end

[ChrisRackauckas]

Here's another way to describe it @dlfivefifty. Let's say I allow the user to tell me what type of matrix the Jacobian will be, and I will build a few caches for that. If I want the user to be able to pass me construction information like (size(A), bandwidths(A)), then how do I know it's a banded matrix? Is the proper API to have the user pass a datatype and a tuple of tuple of numbers, and document what a tuple of tuples to pass along with which datatype to build a banded matrix, and do the same for a Toeplitz matrix, and ... ? How is that not just begging to be a type?

The other way of course is to require that the user just passes A, the already constructed matrix. But this means that your input a lot heavier since now it can hold at most a dense matrix, when in reality all you needed to know here was some generalized idea of size. It also doesn't help you with the first construction.

Instead of phrasing it as whether a special type is necessary, I think the better question is: why should we have to special case for every possible sparsity structure instead of just formalizing an API for sparsity structures?

[dlfivefifty]

Is there an example not of the form foo(MyMatrixType, pars)? That is, you always have a type and a list of parameters...of course you can wrap this in a type but I think it's overkill as then every matrix needs it's own type.

[ChrisRackauckas]

Well everything can always be a type and a list of parameters otherwise you couldn't construct it in the first place, but that gets painful. For example, what about a non-binary tree abstract array, like MultiScaleArrays.jl? That could easily be a huge type that's hard to construct, so I don't see why we should be constructing a false tree to be able to arrayparameters to then construct a new one with appropriate units. You can make that API be that you're passing tuple of tuple of ... of tuples of sizes into typeof, but here (it's not sparse structure but "contiguousness structure") but Julia begins to have issues when you start having an API with huge tuples or function signatures. What about with RaggedArrays.jl if we wanted to allow OffsetArrays.jl?

While you can keep on passing on tuples of unknown structure information, one thing that begins to happen is you cannot query this information until the type is actually created. For example, if we have a function where we get the sparsity structure from the user

function f(x,y;pars=nothing,pars_type = nothing)
  # do some stuff
  J = pars_type(pars)
end

# User code
pars=(...) # all of the stuff for a SparseMatrixCSC
f(x,y,pars=(...),pars_type=SparseMatrixCSC) 

One thing to note is that you cannot get the size of the matrix they want to construct until after you construct it, unless you special case for every pars_type. You're always going to naturally use the two together, but now on you to pass around two things for every sparsity/continguousness specification, and any place you want to use it you need do define g(pars_type,pars) and dispatch on pars_type (since they need to handle pars differently, unless all you're going to do is construct of course).

[dlfivefifty]

OK, MultiScaleArrays.jl is a good example (where you can use a type to hide the types of the parameters, unlike Tuple which exposes the types and requires recompilation).

But you could just still fit in the above framework:

arrayparameters(A::MultiScaleArray) = (MultiScaleArrayInformation(A),)

thus avoiding the need for lots of extra types for Base arrays.

[ChrisRackauckas]

But that still requires having the constructed array in order to call arrayparameters, which is exactly the thing I would like to try and avoid.

[dlfivefifty]

Where do you get the Array parameters then? In the banded matrix case, I can only imagine two scenarios: 1. Generic code where you have a matrix but don’t know it’s type 2. You know you want a banded matrix so you specify the bandwidths. I cannot imagine a scenario where you have generic code but need to will into existence an array with specific parameters, because where do those parameters come from?

…

Sent from my iPhone

On 8 Mar 2018, at 15:16, Christopher Rackauckas ***@***.***> wrote: But that still requires having the constructed array in order to call arrayparameters, which is exactly the thing I would like to try and avoid. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[ChrisRackauckas]

The Jacobian or the bbprime matrix in stiff ODE/SDE solvers are usually the only size n^2 matrices (and n^2 * m for m Wiener process for bbprime) in the entire program. Less if sparse of course, but these are by far the largest arrays. When units are involved its tricky, so it would be tough to have that as part of the API (and then you can't resize! most arrays too), so there are really two options. Either take in a "prototype" array which you ignore to change the element types when necessary (which is what's currently done), or have a way to construct the matrix they want from minimal information. The latter is how SUNDIALS does it, but they have to special case for each of the matrices they support. The former is how we currently do it because we don't have generic matrix constructors in Julia. This means that we have two copies of the largest arrays when we really only need one. arrayparameters doesn't fix the problem since it's still asking for a "prototype" array in order to do the construction.

[dlfivefifty]

Right, so you are given the necessary information as an argument to the function. I don't see why there needs to be extra types in Base to do this, as you can always take the type and the parameters as two arguments (or a 2-tuple):

function construct(::Type{AT}, pars) where AT<:AbstractArray =
   AT(uninitialized, pars...)

will do exactly why you want, you just need to give it two arguments instead of one.

[ChrisRackauckas]

But going back to one of the first questions, how does that let me query details of the future matrix, like its size, directly from pars? Do I need to write size(::Type{AT},pars) dispatches? That's type-piracy so it's even against recommendations to do this.

[dlfivefifty]

It's nothing that can't be worked around: mysize(pars) = first(pars) will work for 90% cases.

The issue is that adding a special type in Base just adds more code and more complexity...and it seems unlikely to happen in v1.0 given the tight timetable.

[ChrisRackauckas]

I don't get why we should settle for a method that only works on some arrays (and works without a defined interface).

And I never asked for 1.0. This issue becomes pretty apparent when writing code for a generic matrix input (and only generic, so most Base arrays are fine), so I assumed I'd have to wait until there are more people running into this problem for it to pick up steam. FWIW the proposal for Jacobian types in DiffEq right now uses the prototype approach in order to get around this for now (SciML/DifferentialEquations.jl#220)

[dlfivefifty]

I've changed my mind now that BandedMatrices.jl supports general backends, and how hard it is to create general matrices. I think @ChrisRackauckas suggestion is a good one.