Make sure it supports SparseConnectivityTracer

[ChrisRackauckas]

We will want to make sure SCT works on these functions like it does for DataInterpolations.jl

Pinging @adrhill so he's in the know.

[adrhill]

Since we have some ongoing discussions in adrhill/SparseConnectivityTracer.jl#234 related to DataInterpolations.jl's very dynamic behavior, it would be good if NDInterpolations restricted itself to somewhat predictable input and output types.

A negative example from DI that was given in the issue is:

t = [1.0, 2.0, 3.0]
u = ["one", (2,2), fill(3, 3, 3)]
interp = ConstantInterpolation(u, t)

interp(1.0)              # "one"
interp(2.0)              # (2, 2)
interp(3.0)              # fill(3, 3, 3)
interp([1.0, 2.0, 3.0])  # Any["one", (2, 2), fill(3, 3, 3)]

While local sparsity patterns probably work on this code, it is simply impossible for SCT to ever provide a global sparsity pattern in this example.

For us to be able to implement alternative evaluation via operator-overloading, we need to (at least) be able to infer the output type and size without running the interpolation.

[SouthEndMusic]

I have implemented 2 types of evaluation for an NDInterpolation:

• Single point evaluation, this is analogous to how interpolation objects from DataInterpolations are callable:

NDInterpolations.jl/src/NDInterpolations.jl

Lines 53 to 69 in 718dbe1

	# In place single input evaluation
	function (interp::NDInterpolation{N_in})(
	out::Union{Number, AbstractArray{<:Number}},
	t::Tuple{Vararg{Number, N_in}};
	derivative_orders::NTuple{N_in, <:Integer} = ntuple(_ -> 0, N_in)
	) where {N_in}
	validate_derivative_orders(derivative_orders)
	idx = get_idx(interp.interp_dims, t)
	@assert size(out)==size(interp.u)[(N_in + 1):end] "The size of out must match the size of the last N_out dimensions of u."
	_interpolate!(out, interp, t, idx, derivative_orders)
	end
	# Out of place single input evaluation
	function (interp::NDInterpolation)(t::Tuple{Vararg{Number}}; kwargs...)
	out = make_out(interp, t)
	interp(out, t; kwargs...)
	end

This already works with ForwardDiff, so it should also work with SparseConnectivityTracer.

• Multi point evaluation, which uses pre-allocated input stored in the interpolation dimensions:

NDInterpolations.jl/src/interpolation_parallel.jl

Lines 105 to 132 in 718dbe1

	@kernel function eval_kernel(
	out,
	@Const(A),
	derivative_orders,
	eval_grid
	)
	N_in = length(A.interp_dims)
	N_out = ndims(A.u) - N_in
	k = @index(Global, NTuple)
	if eval_grid
	t_eval = ntuple(i -> A.interp_dims[i].t_eval[k[i]], N_in)
	idx_eval = ntuple(i -> A.interp_dims[i].idx_eval[k[i]], N_in)
	else
	t_eval = ntuple(i -> A.interp_dims[i].t_eval[only(k)], N_in)
	idx_eval = ntuple(i -> A.interp_dims[i].idx_eval[only(k)], N_in)
	end
	if iszero(N_out)
	out[k...] = _interpolate!(
	make_out(A, t_eval), A, t_eval, idx_eval, derivative_orders)
	else
	_interpolate!(
	view(out, k..., ..),
	A, t_eval, idx_eval, derivative_orders)
	end
	end

This probably poses a problem for differentiating the multi point output with respect to the pre-allocated input for packages that need special number types like ForwardDiff and SparseConnectivityTracer, unless the interpolation object is constructed within the function to be differentiated.

Regarding the discussion around DataInterpolations: within an NDInterpolation object, u is an N dimensional AbstractArray, and there are N_in <= N interpolation dimensions. Hence the output is an N_out = N - N_in dimensional array, which is a scalar for N_out = 0. N_in and N_out are type parameters of the NDInterpolation type. Moreover, for N_out > 0, the output of the interpolation has the same size as the last N_out dimensions of u:

NDInterpolations.jl/src/interpolation_utils.jl

Lines 14 to 16 in 718dbe1

	function get_output_size(interp::NDInterpolation{N_in}) where {N_in}
	size(interp.u)[(N_in + 1):end]
	end

[SouthEndMusic]

@adrhill how can we make sure that problems are avoided? What kinds of tests should be added? What does the restriction you mention look like? For single point out of place evaluation the output is allocated here:

NDInterpolations.jl/src/interpolation_utils.jl

Lines 90 to 103 in 0cb8e4c

	function make_out(
	interp::NDInterpolation{N_in, 0},
	t::NTuple{N_in, >:Number}
	) where {N_in}
	zero(promote_type(eltype(interp.u), map(typeof, t)...))
	end
	function make_out(
	interp::NDInterpolation{N_in},
	t::NTuple{N_in, >:Number}
	) where {N_in}
	similar(
	interp.u, promote_type(eltype(interp.u), map(eltype, t)...), get_output_size(interp))
	end

with

NDInterpolations.jl/src/interpolation_utils.jl

Lines 79 to 81 in 0cb8e4c

	function get_output_size(interp::NDInterpolation{N_in}) where {N_in}
	size(interp.u)[(N_in + 1):end]
	end

To be safe it might be better to add map(interp_dim -> eltype(interp_dim.t), interp.interp_dims) to promote_type.

Should it be validated that the element type of u (and t) is a concrete scalar type? I guess the above already assumes that. The second type parameter of an NDInterpolation is the output dimensionality, yielding a scalar output if the output dimensionality is $0$.

Restricting the element type of u to a scalar still allows for vector output, just by making ndims(u) larger than the number of input dimensions.

[SouthEndMusic]

Although, I guess for scalar computations more care has to be taken that the type of the zero that gets 'allocated' is preserved, avoiding e.g.

x = 0
x += 1.0
typeof(x) == Float64

[adrhill]

@adrhill how can we make sure that problems are avoided?

SCT implements an "alternative evaluation" of DI's (or NDI's) methods, in some cases completely side-stepping the original code. Since DI/NDI's input-output behavior is very dynamic (see "This example is pretty cursed" in adrhill/SparseConnectivityTracer.jl#234 (comment)), an ideal SCT extension would have to implement a ton of methods overloads.

Instead, I would suggest starting out by supporting the most relevant subset of input types using strict typing, i.e. only supporting interpolations on matrices and vectors of <:Number.

What kinds of tests should be added?

You want to make sure that for different scalar/array/matrix/tuple inputs of tracers t, the output interp(t) returns a conservative sparsity pattern of the correct size. If the sparsity pattern has too many zeros, the sparse AD computation of Jacobians/Hessian will be incorrect.