`getdata2d` relies on matplotlib

[henry2004y]

https://github.com/henry2004y/Batsrus.jl/blob/a540abcb258cf878522074e071fd57f69f4b8f09/src/plot/utility.jl#L11

Currently getdata2d depends on matplotlib to do interpolations. In the future we really need to figure out how to do it native Julia. Related to #14.

[henry2004y]

This is also related to a newly wanted feature of extracting variables at a fixed spatial location from slices, presumably in an unstructured mesh.

[henry2004y]

I found one solution using NaturalNeighbours.jl. It is about 10 times slower than Matplotlibs's triangular interpolation on a small dataset, but at least we have something in native Julia now.

[DanielVandH]

I found one solution using NaturalNeighbours.jl. It is about 10 times slower than Matplotlibs's triangular interpolation on a small dataset, but at least we have something in native Julia now.

You'll probably get a faster speed if you replace https://github.com/henry2004y/Batsrus.jl/blob/2ae4d98052efc27f5c8951ae9038b2d4ddf55632/src/plot/utility.jl#L130 to instead be itp(x, y), where x and y are the vectors you're evaluating at, similar to what is done in the README https://github.com/DanielVandH/NaturalNeighbours.jl. You'll be able to hook into the multithreading that way.

[henry2004y]

I found one solution using NaturalNeighbours.jl. It is about 10 times slower than Matplotlibs's triangular interpolation on a small dataset, but at least we have something in native Julia now.

You'll probably get a faster speed if you replace

https://github.com/henry2004y/Batsrus.jl/blob/2ae4d98052efc27f5c8951ae9038b2d4ddf55632/src/plot/utility.jl#L130

to instead be itp(x, y), where x and y are the vectors you're evaluating at, similar to what is done in the README https://github.com/DanielVandH/NaturalNeighbours.jl. You'll be able to hook into the multithreading that way.

Thank you for your tips! I tried and it is indeed a bit faster, with 4 thread giving a speedup of ~ 2. However, the current way I did it requires extra array allocations:

function interpolate2d_generalized_coords2(X::T, Y::T, W::T,
   plotrange::Vector{<:AbstractFloat}, plotinterval::Real) where
   {T<:AbstractVector}
   xi = range(plotrange[1], stop=plotrange[2], step=plotinterval)
   yi = range(plotrange[3], stop=plotrange[4], step=plotinterval)
   itp = interpolate(X, Y, W)
   #Wi = [itp(x, y; method=Triangle()) for y in yi, x in xi]::Matrix{eltype(W)}
   Xi, Yi = Batsrus.meshgrid(xi, yi)
   x_ = vec(Xi)
   y_ = vec(Yi)
   Wi = itp(x_, y_; method=Triangle())

   xi, yi, Wi
end

which is something I may wish to avoid (say using LazyGrids.jl, but I haven't tried yet). However, even with multithreading it is still a bit slower than the Matplotlib triangulation interpolation implementation. Do you know the reason @DanielVandH ?

Here is a minimal working example for you to test:

using NaturalNeighbours
using StableRNGs
using PyPlot

function meshgrid(x, y)
   X = [x for _ in y, x in x]
   Y = [y for y in y, _ in x]

   X, Y
end

## The data 
rng = StableRNG(123)
f = (x, y) -> sin(x * y) - cos(x - y) * exp(-(x - y)^2)
x = vec([(i - 1) / 9 for i in (1, 3, 4, 5, 8, 9, 10), j in (1, 2, 3, 5, 6, 7, 9, 10)])
y = vec([(j - 1) / 9 for i in (1, 3, 4, 5, 8, 9, 10), j in (1, 2, 3, 5, 6, 7, 9, 10)])
z = f.(x, y)

## The interpolant and grid 
itp = interpolate(x, y, z; derivatives=true)
xg = LinRange(0, 1, 100)
yg = LinRange(0, 1, 100)
_x = vec([x for x in xg, _ in yg])
_y = vec([y for _ in xg, y in yg])

## Evaluate some interpolants 
triangle_vals = itp(_x, _y; method=Triangle())
println("NaturalNeighbours.jl:")
@time itp(_x, _y; method=Triangle());

triang = matplotlib.tri.Triangulation(x, y)
# Perform linear interpolation on the triangle mesh
interpolator = matplotlib.tri.LinearTriInterpolator(triang, z)
Xi, Yi = meshgrid(xg, yg)
Wi = interpolator(Xi, Yi)
println("Matplotlib:")
@time Wi = interpolator(Xi, Yi);
println("")

What I get is

NaturalNeighbours.jl:
  0.010370 seconds (37.41 k allocations: 1.056 MiB)
Matplotlib:
  0.000652 seconds (40 allocations: 80.031 KiB)

I will post an issue in your repo~

[DanielVandH]

which is something I may wish to avoid (say using LazyGrids.jl, but I haven't tried yet)

Yeah the need for vec is annoying. I could probably fix this at some point so that matrices can also be provided if I found the time.