sisl/gridinterpolations.jl

Multidimensional grid interpolation in arbitrary dimensions

GridInterpolations

This package performs multivariate interpolation on a rectilinear grid. At the moment, it provides implementations of multilinear and simplex interpolation. As of benchmarks in December 2016, multilinear interpolation performs fastest and with the most accuracy.

The following image visualizes grid-based interpolation in two dimensions, with shape of interpolater for (−0.3,0.8) inscribed. The small dots reflect the interpolation's estimate for sin(x)+2cos(y)+sin(5xy), which is the underlying reward function approximated by the large dot lattice.

For a description of multilinear and simplex interpolation see: Scott Davies, Multidimensional Triangulation and Interpolation for Reinforcement Learning, Advances in Neural Information Processing Systems, Cambridge, MA: MIT Press, 1997. pdf

There are some related packages, such as Grid.jl and Interpolations.jl.

Installation

Start Julia and run the following command:

Pkg.add("GridInterpolations")

Usage

To use the GridInterpolations module, begin your code with

using GridInterpolations

Interpolation

Create two-dimensional interpolation grids, a data array, and a point of interest:

grid = RectangleGrid([0., 0.5, 1.],[0., 0.5, 1.])  	# rectangular grid
sGrid = SimplexGrid([0., 0.5, 1.],[0., 0.5, 1.])	# simplex grid
gridData = [8., 1., 6., 3., 5., 7., 4., 9., 2.]   	# vector of value data at each cut
x = [0.25, 0.75]  									# point at which to perform interpolation

Perform interpolation on the rectangular grid:

julia> interpolate(grid,gridData,x)
5.25

Or interpolate on the simplex grid:

julia> interpolate(sGrid,gridData,x)
6.0

Compute interpolants for the grids:

julia> sGrid = SimplexGrid([0., 0.5, 1.],[0., 0.5, 1.])
[[0.0,0.5,1.0],[0.0,0.5,1.0]]

julia> interpolants(sGrid, x)
([4,5,8],[0.5,0.0,0.5])

Convert an index to a Grid coordinate:

julia> ind2x(grid, 3)
2-element Array{Float64,1}:
 1.0
 0.0

Number of vertices in the grid:

julia> length(grid)
9

Number of dimensions:

julia> dimensions(grid)
2

Multi-dimensional indexing using Cartesian coordinates:

julia> [grid[c] for c in CartesianIndices((3,3))]
3×3 Array{Array{Float64,1},2}:
 [0.0, 0.0]  [0.0, 0.5]  [0.0, 1.0]
 [0.5, 0.0]  [0.5, 0.5]  [0.5, 1.0]
 [1.0, 0.0]  [1.0, 0.5]  [1.0, 1.0]

or multi-dimensional indices

julia> grid[2,2]
2-element Array{Float64,1}:
 0.5
 0.5

Sequential iteration over grid points:

julia> for x in grid
           # do stuff
       end

Credits

Contributors to this package include Maxim Egorov, Eric Mueller, and Mykel Kochenderfer.

gridinterpolations.jl's Issues

NearestGrid

It would be nice to have something like this:

type NearestGrid <: AbstractGrid

It just does nearest neighbor interpolation. It should be pretty simple to implement.

iteration

It would be slick to be able to iterate through the grid points like this

for p in grid
    ...
end

though it is not too difficult now with ind2x

Error tagging new release

The tag name "v0.1" is not of the appropriate SemVer form (vX.Y.Z).
cc: @MaximeBouton

Change tests to use Interpolations.jl instead of Grid.jl

Grid.jl is deprecated after 0.5. See JuliaLang/METADATA.jl#5749

Grid should be removed from REQUIRE and runtests.jl should be modified to use Interpolations.jl.

Tag new version

I think this package has enough changes (especially those from March 18 and 19) to warrant tagging a new version, especially since the types are parameterized by the dimension now.

@Shushman , you probably have been using this the most. What do you think? are there any changes coming up?

Compat with StaticArrays

For some reason the general registry requires StaticArrays = "0.0.0 - 0.12" for this package, hence it cannot be used together with packages requiring StaticArrays >1.0

complex conjugation when interpolating complex valued function

While I can interpolate complex valued functions, e.g.,
grid = RectangleGrid([0., 1., 2.])
gridData = [1.0im, 2.0im, 3.0im]

the interpolated function returns complex conjugated values: E.g.,
interpolate(grid,gridData, [0.0])
gives 0.0-1.0im instead of 1.0im

Improve runtime performance of simplex interpolation

Simplex interpolation runs surprisingly slowly in benchmarks from December 2016. See images below. Simplex is blue; multilinear interpolation is red. Simplex interpolation is always slower than multilinear interpolation for up to 1 million points in a grid (the largest number of points tested).

I profiled the code to look into this and I marked the lines in the code that execute a large number of times. sortperm! is the biggest offender (currently line 276 of GridInterpolations.jl). Switching to the native Julia sortperm has minimal effect, which suggests the call needs to be removed entirely to stand a chance of faster performance with simplex than multilinear interpolation. To my fast-and-dirty read of the code, this may in fact be possible. It seems like sortperm might only be used to present the data in a conveniently sorted way to the user, a piece that I suspect we could drop from the API without further ado since multilinear interpolation doesn't uphold that requirement.

There is a branch named simplex_optimization, which I haven't looked at. There is a chance this issue is being addressed there.

"Where multilinear interpolation is almost instantaneous, simplex interpolation is between twice and ten times slower for all tests. It is possible that above 1 million points – when the lattice stops fitting in RAM – that simplex interpolation may perform better, although at that point, a coarser lattice may be more appropriate."

Use AbstractArray where possible

Since StaticArrays is pretty useful, we should use AbstractArray instead of Array in method argument type annotations so that this package is compatible with them. It should cause no performance degradation.

Cannot interpolate complex numbers when passing a 2d array

When I pass a 2d array with complex numbers to interpolate, the code fails. I tries to convert the complex numbers to Float64.

Grid with Dates

I tried to create a RectangleGrid(t,p), where t is a Date, p is a float with no luck, I got error "ERROR: MethodError: Cannot convert an object of type Date to an object of type Float64". Is there way to resolve it?

I tried to redefine RectangleGrid struct, but I also failed with it at the moment.

ind2x!(x, grid ind)

To be consistent with the julia style guide https://docs.julialang.org/en/v1/manual/style-guide/index.html#Write-functions-with-argument-ordering-similar-to-Julia-Base-1, we should have

ind2x!(x, grid, ind)

rather than

ind2x!(grid, ind, x)

This should not be a terrible change since we can do

@deprecate ind2x!(grid::AbstractGrid, ind, x) ind2x!(x, grid, ind)

There may be other functions that should be changed as well.

handling of NaNs in interpolants()

What is the expected / desired behavior for the last line in this snippet?

using GridInterpolations
r = GridInterpolations.RectangleGrid(1:4, 1:4)
interpolants(r, [2.5, 3.5])  # business as usual
interpolants(r, [2.5, NaN])  # expected behavior?

The output is ([-2, -1, 2, 3], [NaN, NaN, NaN, NaN]).
While NaN weights seem to make sense, the indices seem kind of dangerous (imagine someone is just checking for the support, not the actual weights).

What would be a sane behavior? Maybe ([0], [NaN]) or ([NaN],[NaN])? Or just checking for NaNs at the beginning of the call? Throwing an error?

@MaximeBouton

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Generate rectangle of dimension N

This might be already doable but I cannot figure it out.
I am trying to generate a N-dimensional rectangular grid where N is a parameter.

So if my grid is

grid_1d = [0., 0.5, 1.]
N = 2
grid = RectangleGrid(grid_1d, grid_1d)

Of course this is a little annoying for more dimensions

N=6
grid = RectangleGrid(grid_1d, grid_1d, grid_1d, grid_1d, grid_1d, grid_1d)

I looked into the constructor but was not able to figure out an elegant way of programmatically implement my rectangle.

Any help would be appreciated!

---

Thank you for the great and useful package.

Base.show() doesn't differentiate between simplex and rectangular grids

Hi all,

I've been using GridInterpolations in one of my projects, and I have been saving string representations of the grids that I've been using so that I can remember what I did later. As the code stands, Base.show() simply returns a vector of vectors representing the grid, but it doesn't indicate whether the grid uses simplex or rectangular interpolation. Can we change Base.show() so that it returns an unambiguous description of the grid, or, even more ideally a string that can be parsed and used directly to construct an identical grid object?

i.e. Base.show(mySimplexGrid) returns "SimplexGrid([1.0,2.0,3.0],[1.0,2.0,3.0])"

If this sounds like a good idea, I'll submit a pull request for it this afternoon.

Thanks,

Zach

Scipy Performance Comparison

This issue discusses performance comparison for multidimensional grid interpolation. GridInterpolations.jl is compared to interpolation done using scipy in Python.

2D interpolation of a 3x3 grid for 10,000 points:
GridInterpolations.jl with Simplex: ~0.065s
Scipy Interpolation: ~0.5s

3D interpolation of a 3x3x3 grid for 1,000,000 points:
GridInterpolations.jl with Simplex: ~3.8s
Scipy Interpolation: ~54s

GridInterpolations.jl is ~7.5X faster for 2D interpolation.
GridInterpolations.jl is ~14X faster for 3D interpolation.

Some details about the scipy implementation:

• Using scipy.spatial.Delaunay for triangulation (not timing this part)
• Using scipy.interpolate.LinearNDInterpolator for the interpolation (timing only this part)
• LinearNDInterpolation is calling the qhull lib

Interface for gridData

Wouldn't it be better if the gridData is provided as a mufti-dimensional array rather than a vector? Consider the example given in the README:

gridData = [8. 1. 6. 3. 5. 7. 4. 9. 2.]             # value data at each cut

It is not clear whether the data is provided in row-major or col-major form (and this gets more complicated for higher dimensions). An alternative is to have an interface where

gridData = [8. 1. 6.; 3. 5. 7.; 4. 9. 2.]

so that it is clear that gridData[0, 0.5] = 3. etc. Internally, this could be as simple vectorizing the array and then proceedings as one is doing now.