JuliaMath / Interpolations.jl

Fast, continuous interpolation of discrete datasets in Julia
http://juliamath.github.io/Interpolations.jl/
Other
524 stars 110 forks source link

update coefficient matrix #464

Closed koehlerson closed 2 years ago

koehlerson commented 2 years ago

I have a large interpolation in 4 dimensions. I want to iterate over each grid point and update the value of the 4 dimensional coefficient matrix, such that I don't have to construct a new interpolation. Is there a more recommended way than doing this:

itp = LinearInterpolation((sample, sample, sample, sample), [f([x1 x2;y1 y2]) for x1 in sample, x2 in sample, y1 in sample, y2 in sample],extrapolation_bc = Line());

itp.itp.itp.coefs[1,1,1,1] = 5.0

Does this in general the thing I would like to do? Or is there some reordering going on under the hood? If I check for the first grid point, it yields the exact result, but I'm unsure if that's just by chance

itp(0.1,0.,0.,0.1) = 5.0 # first grid point is at (0.1, 0.0, 0.0, 0.1)
mkitti commented 2 years ago

I think you might be able to only do this for the Constant and Linear case. No prefiltering is done.

https://github.com/JuliaMath/Interpolations.jl/blob/d11ddd16aa36ba9b4397be0624f74ef55da40b3a/src/b-splines/prefiltering.jl#L33

I would not say this is supported usage, but it may work in your particular case.

koehlerson commented 2 years ago

I reopened this, since I maybe need to go quadratic with my interpolation. Is it possible to update the coefficient matrix, i.e. recompute the original, non-filtered coef matrix indices? Would be happy to contribute, but I need some guidance

koehlerson commented 2 years ago

Maybe to coin the question differently: Is there a way to turn off the prefiltering system? Looking at https://ieeexplore.ieee.org/document/875199 which is linked in https://github.com/JuliaMath/Interpolations.jl/blob/master/doc/Interpolations.jl.ipynb it seems that this step tries to overcome something similar (or equivalent?) to https://en.wikipedia.org/wiki/Runge%27s_phenomenon . I'd be happy to have the trade off introducing an error while being able to reuse the objects and just update their coef matrices.

Maybe someone has a totally different approach for me. My core problem is, that my coefs are changing (in some places) and reconstructing the interpolation everytime seems like a lot of overhead.

koehlerson commented 2 years ago

Sorry for spamming so much, but I need a conceptual sketch quite urgent and I think I found a workaround myself:

julia> itp = Interpolations.BSplineInterpolation(Float64,A,BSpline(Quadratic(Line(OnGrid()))),axes(A))
julia> A == itp.coefs
true

scaling, evaluating, updating the coefs etc. works for me this way. For other readers, note: this does not preserve interpolation properties, it's just a really quick and dirty sketch for me to see if this helps to solve my problem. It's still an open problem how to update the coefficients of a (true) bspline interpolation efficiently. Since prefiltering the system solves an system of equations, will be quite hard to reconstruct what happened while solving.

mkitti commented 2 years ago

If you are trying to address Runge phenomenon, then switching to Chybebshev interpolation might be the better strategy. See J.P. Boyd, Chebfun, Approxfun.jl.

koehlerson commented 2 years ago

Ah no no! I misunderstood why we do prefiltering. I'm used to lagrange interpolations where you have by construction of the basis function the interpolation property fulfilled, which seems not the case for bsplines, thus the prefiltering. I'm fine with any kind of quadratic/cubic interpolation if I can update the coefficients