A question on extrapolation with Interpolations.jl

Specific Domains Numerics
1

Say that I make an array that I want to do extrapolation on:

using Interpolations
a = [1,2,3,4]
a = interpolate!(a, BSpline(Linear()), OnGrid())
a = extrapolate(a, 3)

Then a[0] = 3 and likewise for any index outside of 1:4.

But what if I wanted indices less than 1 to return 3 and indices greater than 4 to return 2? If I try

a = extrapolate(a, [3,2])

then I get an error.

The reason that I want to do this is because I have Dirichlet boundary conditions that are different on each side. Currently, the closest thing is Flat() boundary conditions, but to my mind that more closely represents Neumann boundary conditions with a zero gradient.

2

I recently faced a similar problem in my project. Eventually I just implemented a new Array type that overrides getindex like:

getindex(A::DirichletArray, I...) =
    checkbounds(Bool, A.a, I...) ? A.a[I...] : A.oob_val

This kind of approach could easily be extended for your case by checking indices appropriately.

You can find details on the Array interface in the manual.

4

Slightly a hack, but setting a knot infinitesimally to the left and right of the ranges (in the Float64 sense), might do the job with Flat(). Sample code:

left, right = 3,2
a = [left,1,2,3,4,right]
knots = ([prevfloat(1.0);1.0:1.0:4.0;nextfloat(4.0)],)
a_itp = interpolate(knots, a, Gridded(Linear()))
a_ext = extrapolate(a_itp,Flat())
@assert a_ext[0.0]==3.0
@assert a_ext[4.5]==2.0
@assert a_ext[4]==4
5

Thanks, that definitely seems like a possibility :smiley:.

Although, I quite like some of the other features that Interpolations.jl has like using the index 1.5 for getting the value half way between indexes 1 and 2. Having said that, I have only been using linear interpolation so far. So your solution would actually work if I don’t try anything more complicated than what I already have.

6

Thanks, that is actually a really cool idea!

There are two weird edge cases that I have noticed so far:

@assert length(a_ext) == 6  # I would expect this to be 4.
@assert a_ext[end] == 2.0   # I would expect this to be 4.0.

I just looked through my code and neither of these would cause problems. But if I change things in the future, then I may expect the above things to behave differently which might introduce subtle bugs.

Just to be clear, I am not saying that your code is buggy, I am just saying that there are certain cases where my expectations differ from the behavior that you defined.

7

You can totally combine Interpolations.jl with my approach. Just wrap the interpolate() object in the DirichletArray and add your boundary logic in the getindex method. It will simply pass on the Float64 indices and do the interpolation.

1 Like
8

Oh cool thanks! I didn’t think of that… Here is what I put together based on your recommendation.

using Interpolations

# I will add the methods to make this a subtype of AbstractArray later...
type DirichletArray # <: AbstractArray
    a::AbstractArray
    oob_val::AbstractArray
end

function Base.getindex(A::DirichletArray, I...)
    checkbounds(Bool, A.a, I...) && return A.a[I...]
    I... < 1 ? A.oob_val[1] : A.oob_val[2]
end

a = [1,2,3,4]
a = interpolate!(a, BSpline(Linear()), OnGrid())

A = DirichletArray(a, [3,2])

A[1]    # 1
A[3.5]  # 3.5
A[0]    # 3
A[5.1]  # 2

This array behaves just the way that I want it to!

9

Yeah, that’s basically what I meant. Btw if you want type stability (and therefore good performance) I would recommend something like this:

type DirichletArray{S <: AbstractArray{T, N}} <: AbstractArray{T, N}
    a::S
    oob_val::Tuple{T, T}
end
10

That code gives the error

ERROR: UndefVarError: T not defined
11

Oh sorry, I didn’t test it. The important thing is not having unstable fields. By using introducing S you get proper type inference for functions operating on the array.

immutable DirichletArray{S <: AbstractArray, T <: Number} <: AbstractArray
    a::S
    oob_val::Tuple{T, T}
end
1 Like
12

Thanks for your help! :grinning:

13

You need to include the AbstractArray’s two parameters when you subtype it. This is now an error on 0.6, and on 0.5 you’ll end up with a bunch of broken fallbacks. See https://docs.julialang.org/en/stable/manual/interfaces/#abstract-arrays for more details.

That means that it’s a bit more complicated:

immutable DirichletArray{T, N, S <: AbstractArray} <: AbstractArray{T, N}
    a::S
    oob_val::Tuple{T, T}
end
# And now you need an outer constructor that computes those parameters for you:
DirichletArray{T,N}(a::AbstractArray{T,N}, oob_val::Tuple{T,T}) = DirichletArray{T,N,typeof(a)}(a, oob_val)
1 Like
14

Thanks @mbauman! With the previous example I had to do things like

Base.print_matrix(io::IO, mat::DirichletArray, args...) =
                   Base.print_matrix(io, mat.a, args...)

To stop Julia from complaining about method errors, but now things seem quite hunky-dory! Looking at the code in https://github.com/JuliaMath/Interpolations.jl/tree/master/src/extrapolation, I see that this is not how the Interpolations.jl folk do things. However, I am happy with this current solution because it behaves the way that I want it to for my particular problem.

If anyone reading this is interested in a MWE, then this is what I have put together so far:

using Interpolations

immutable DirichletArray{T, N, S <: AbstractArray} <: AbstractArray{T, N}
    a::S
    oob_val::Tuple{T, T}
end
# And now you need an outer constructor that computes those parameters for you:
DirichletArray{T,N}(a::AbstractArray{T,N}, oob_val::Tuple{T,T}) = DirichletArray{T,N,typeof(a)}(a, oob_val)

Base.size(A::DirichletArray) = size(A.a)
Base.linearindexing{T<:DirichletArray}(::Type{T}) = Base.LinearFast()

function Base.getindex(A::DirichletArray, I...)
    checkbounds(Bool, A.a, I...) && return A.a[I...]
    I... < 1 ? A.oob_val[1] : A.oob_val[2]
end

a = [1.0,2.0,3.0,4.0]
a = interpolate!(a, BSpline(Linear()), OnGrid())

A = DirichletArray(a, (3.0,2.0))

A[1]    # 1.0
A[3.5]  # 3.5
A[0]    # 3.0
A[5.1]  # 2.0

Thanks again for all of the help everyone. :grinning:

1 Like
15

Thanks for figuring it out. This has actually been bothering me for a while. Maybe this could be added to the documentation as a “how to build wrapper arrays?” tutorial. Even a blog post would probably help a lot of people trying to figure this out.

16

Fortunately 0.6 will make that particular snag much more obvious.

But you’re right — the documentation here can always be improved. I know the tables on the interfaces page are intimidating, but the prose is meant to be quite accessible and tutorial-like. I chose the examples I did because I thought they were the best at gradually introducing the functionality… but the wrapper array is indeed a very common case and a frequent question. Often those who have recently learned how things work are the best at documenting issues like this since you know where you looked and what would have been helpful… wanna give it a shot?

17

Sure, I will draft something up. One question though: Why does the naive version not work? I.e. this version:

type DirichletArray{S <: AbstractArray{T, N}} <: AbstractArray{T, N}
    a::S
    oob_val::Tuple{T, T}
end

I really don’t get why the compiler cannot infer T during construction. What am I missing?