# Developer documentation

Interpolations provides flexibility without compromising on performance by exploiting metaprogramming to generate streamlined code. However, for people new to metaprogramming this can can be a barrier. Fortunately, with a few tips a lot of the mystique goes away.

## Looking under the hood

First let's create an interpolation object:

```
julia> using Interpolations
julia> A = rand(5)
5-element Array{Float64,1}:
0.74838
0.995383
0.978916
0.134746
0.430876
julia> yitp = interpolate(A, BSpline(Linear()), OnGrid())
5-element Interpolations.BSplineInterpolation{Float64,1,Float64,Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid}:
0.74838
0.995383
0.978916
0.134746
0.430876
```

We can use this object to learn a lot about how Interpolations works.
For example, the key functionality provided by `yitp`

is `getindex`

, i.e., `itp[3.2]`

.
Where is this implemented?

```
julia> @which yitp[3.2]
getindex{T,N}(itp::Interpolations.BSplineInterpolation{T,N,TCoefs,IT<:Interpolations.BSpline{D<:Interpolations.Degree{N}},GT<:Interpolations.GridType},xs::Real) at /home/tim/.julia/v0.4/Interpolations/src/b-splines/indexing.jl:42
```

Your specific output (and especially the line number) may differ, but the point is that you've now found out where this is implemented. If you take a look at that function definition, you might see something like this:

```
@generated function getindex{T,N}(itp::BSplineInterpolation{T,N}, xs::Real)
if N > 1
error("Linear indexing is not supported for interpolation objects")
end
getindex_impl(itp)
end
```

This is a generated function, and you'll need to familiarize yourself with how these work.
The "interesting" part of the function is the call to `getindex_impl`

; we can see the code that gets generated like this:

```
julia> Interpolations.getindex_impl(typeof(yitp))
quote # /home/tim/.julia/v0.4/Interpolations/src/b-splines/indexing.jl, line 7:
@nexprs 1 (d->begin # /home/tim/.julia/v0.4/Interpolations/src/b-splines/indexing.jl, line 7:
x_d = xs[d]
end) # line 11:
begin # /home/tim/.julia/v0.4/Interpolations/src/b-splines/linear.jl, line 5:
@nexprs 1 (d->begin # /home/tim/.julia/v0.4/Interpolations/src/b-splines/linear.jl, line 5:
begin # /home/tim/.julia/v0.4/Interpolations/src/b-splines/linear.jl, line 6:
ix_d = clamp(floor(Int,real(x_d)),1,size(itp,d) - 1) # line 7:
ixp_d = ix_d + 1 # line 8:
fx_d = x_d - ix_d
end
end)
end # line 14:
@nexprs 1 (d->begin # /home/tim/.julia/v0.4/Interpolations/src/b-splines/linear.jl, line 14:
begin # /home/tim/.julia/v0.4/Interpolations/src/b-splines/linear.jl, line 20:
c_d = 1 - fx_d # line 21:
cp_d = fx_d
end
end) # line 17:
@inbounds ret = c_1 * itp.coefs[ix_1] + cp_1 * itp.coefs[ixp_1] # line 18:
ret
end
```

You can see that this code makes use of Base.Cartesian, which you may also need to study.
However, the impact of these macros can be gleaned through `macroexpand`

:

```
julia> using Base.Cartesian
julia> macroexpand(Interpolations.getindex_impl(typeof(yitp)))
quote # /home/tim/.julia/v0.4/Interpolations/src/b-splines/indexing.jl, line 7:
begin
x_1 = xs[1]
end # line 11:
begin # /home/tim/.julia/v0.4/Interpolations/src/b-splines/linear.jl, line 5:
begin
begin # /home/tim/.julia/v0.4/Interpolations/src/b-splines/linear.jl, line 6:
ix_1 = clamp(floor(Int,real(x_1)),1,size(itp,1) - 1) # line 7:
ixp_1 = ix_1 + 1 # line 8:
fx_1 = x_1 - ix_1
end
end
end # line 14:
begin
begin # /home/tim/.julia/v0.4/Interpolations/src/b-splines/linear.jl, line 20:
c_1 = 1 - fx_1 # line 21:
cp_1 = fx_1
end
end # line 17:
begin
$(Expr(:boundscheck, false))
begin
ret = c_1 * itp.coefs[ix_1] + cp_1 * itp.coefs[ixp_1]
$(Expr(:boundscheck, :(Base.pop)))
end
end # line 18:
ret
end
```

This is probably starting to look like something you can read. Briefly, what's happening is:

`floor(Int,x_1)`

gets clamped to the range`1:size(itp,1)-1`

and assigned to`ix_1`

; this is the lower-bound integer grid point for the*first dimension*(this is a one-dimensional problem, but in two or higher dimensions you'd have`ix_2`

, etc.)`ixp_1`

is defined as`ix_1+1`

; this is the upper-bound integer grid point. In Interpolations,`m`

and`p`

often mean "minus" and "plus", meaning the lower or upper grid point.- The fractional part is stored in
`fx_1`

- Position-coefficients
`c_1`

and`cp_1`

associated with the lower and upper grid point are computed from`fx_1`

- The interpolation is performed using the position-coefficients, grid points, and data-coefficients and stored in
`ret`

, which is returned.

As useful exercises:

- Try creating a 2-dimensional linear interpolation object and examine the created code
- Create a
`Quadratic`

interpolation object and do the same

Once you've gotten this far, you probably understand quite a lot about how Interpolations works.
At this point, your best bet is to start looking into the helper functions used by `getindex_impl`

;
once you learn how to define these, you should be able to extend Interpolations to support new algorithms.