# A question on extrapolation with Interpolations.jl

**URL:** <https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669>\
**Category:** Numerics\
**Created:** [May 12, 2017, 7:55am UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669 "2017-05-12T07:55:33Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![JackDevine](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jackdevine/32/1048_2.png) [@JackDevine](https://discourse.julialang.org/u/JackDevine)\
**Post date:** [May 12, 2017, 7:55am UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/1 "2017-05-12T07:55:34Z")

</div>

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

```julia
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

```julia
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.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [May 12, 2017, 9:01am UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/2 "2017-05-12T09:01:26Z")

</div>

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

```julia
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](https://docs.julialang.org/en/release-0.5/manual/interfaces/#abstract-arrays).

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [May 12, 2017, 4:40pm UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/4 "2017-05-12T16:40:41Z")

</div>

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:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![JackDevine](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jackdevine/32/1048_2.png) [@JackDevine](https://discourse.julialang.org/u/JackDevine)\
**Post date:** [May 12, 2017, 11:08pm UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/5 "2017-05-12T23:08:56Z")

</div>

Thanks, that definitely seems like a possibility 😃.

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.

---

<div class="post-metadata">

**Author:** ![JackDevine](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jackdevine/32/1048_2.png) [@JackDevine](https://discourse.julialang.org/u/JackDevine)\
**Post date:** [May 12, 2017, 11:22pm UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/6 "2017-05-12T23:22:32Z")

</div>

Thanks, that is actually a really cool idea!

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

```julia
@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.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [May 13, 2017, 9:06am UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/7 "2017-05-13T09:06:19Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![JackDevine](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jackdevine/32/1048_2.png) [@JackDevine](https://discourse.julialang.org/u/JackDevine)\
**Post date:** [May 13, 2017, 9:41am UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/8 "2017-05-13T09:41:10Z")

</div>

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

```julia
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!

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [May 13, 2017, 10:19am UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/9 "2017-05-13T10:19:21Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![JackDevine](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jackdevine/32/1048_2.png) [@JackDevine](https://discourse.julialang.org/u/JackDevine)\
**Post date:** [May 13, 2017, 10:52am UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/10 "2017-05-13T10:52:13Z")

</div>

That code gives the error

```julia
ERROR: UndefVarError: T not defined

```

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [May 13, 2017, 11:19am UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/11 "2017-05-13T11:19:49Z")

</div>

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.

```julia
immutable DirichletArray{S <: AbstractArray, T <: Number} <: AbstractArray
    a::S
    oob_val::Tuple{T, T}
end

```

---

<div class="post-metadata">

**Author:** ![JackDevine](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jackdevine/32/1048_2.png) [@JackDevine](https://discourse.julialang.org/u/JackDevine)\
**Post date:** [May 13, 2017, 11:29am UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/12 "2017-05-13T11:29:56Z")

</div>

Thanks for your help! 😀

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [May 13, 2017, 4:45pm UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/13 "2017-05-13T16:45:47Z")

</div>

> [@josuagrw](#):
>
> immutable DirichletArray{S \<: AbstractArray, T \<: Number} \<: AbstractArray

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](https://docs.julialang.org/en/stable/manual/interfaces/#abstract-arrays) for more details.

That means that it’s a bit more complicated:

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![JackDevine](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jackdevine/32/1048_2.png) [@JackDevine](https://discourse.julialang.org/u/JackDevine)\
**Post date:** [May 14, 2017, 12:54am UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/14 "2017-05-14T00:54:35Z")

</div>

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

```julia
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](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:

```julia
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. 😀

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [May 14, 2017, 9:32am UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/15 "2017-05-14T09:32:41Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [May 14, 2017, 9:15pm UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/16 "2017-05-14T21:15:21Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [May 15, 2017, 9:20am UTC](https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/17 "2017-05-15T09:20:37Z")

</div>

> [@mbauman](#):
>
> 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?

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

```julia
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?