CliMA / AtmosphericProfilesLibrary.jl

A library of atmospheric profiles
Apache License 2.0
2 stars 0 forks source link

Switch from Dierckx.jl to Interpolations.jl #43

Closed Sbozzolo closed 2 days ago

Sbozzolo commented 5 days ago

Dierckx does not preserve type https://github.com/kbarbary/Dierckx.jl/issues/48 and it is not GPU compatible. It also has some loading costs that can be avoided.

szy21 commented 5 days ago

Yes. I tried to fix it here some time ago: https://github.com/CliMA/AtmosphericProfilesLibrary.jl/pull/39. But it didn't work.

charleskawczynski commented 5 days ago

This is an upstream issue: https://github.com/kbarbary/Dierckx.jl/issues/48

Sbozzolo commented 5 days ago

Is it fine to switch to linear interpolation?

If yes, we can move to Interpolations.jl, which is GPU compatible, and stop using Dierckx.

charleskawczynski commented 5 days ago
julia> x = Float32[1,2];

julia> z = Float32[1,2];

julia> Dierckx.Spline1D(z, x; k = 1)(Float32(2))
2.0
charleskawczynski commented 5 days ago

Perhaps we could switch to Interpolations.jl. Although: https://github.com/JuliaMath/Interpolations.jl/issues/598

charleskawczynski commented 5 days ago

Also, it has quite a few dependencies: https://github.com/JuliaMath/Interpolations.jl/blob/master/Project.toml#L6-L17

charleskawczynski commented 5 days ago

I had actually thought about writing our own, e.g., SimpleInterpolations.jl. Make it GPU compatible and type-stable.

charleskawczynski commented 5 days ago

Honestly, we could roll our own very simply interpolation package. Ah right: https://github.com/CliMA/SimpleInterpolations.jl. I just never got it off the ground

Sbozzolo commented 5 days ago

We are already using Interpolations.jl because it is the only package that is GPU compatible and supports 3D interpolations on irregular grids (needed to read in aerosol data for example).

It seems to be type-preserving:

B = Float32[1, 2, 10, 20, 30]
5-element Vector{Float32}:
  1.0
  2.0
 10.0
 20.0
 30.0

julia> A = Float32[1, 2, 10, 20, 30]
5-element Vector{Float32}:
  1.0
  2.0
 10.0
 20.0
 30.0

julia> itp = interpolate((A, ), B, Gridded(Linear()))
5-element interpolate((::Vector{Float32},), ::Vector{Float32}, Gridded(Linear())) with element type Float32:
  1.0
  2.0
 10.0
 20.0
 30.0

julia> itp(2.0)
2.0

julia> itp(2.f0)
2.0f0
charleskawczynski commented 5 days ago

Sure, switching to Interpolations is fine by me. Alternatively, we could wrap types to recover type stability in this package, which might be a quicker solution, but more of a bandaid.

charleskawczynski commented 5 days ago

Should we rename the issue to "Switch from Dierckx.jl to Interpolations.jl"?