EconForge / interpolation.py

BSD 2-Clause "Simplified" License
123 stars 35 forks source link

Not possible to use `eval_linear` with interpolation option ("xto.*") within jitted function #52

Closed gboehl closed 5 years ago

gboehl commented 5 years ago

Hi and thanks for this very excellent work. This is an extremely useful project!

Big picture: I am numbafying a complete policy function iteration, which I suppose is a rather standard use case for this package. Therefore I am interpolating within a jitted function (via eval_linear). The return vector has the same dimensionality as the input vector, which is multivariate. While eval_linear with univariate output can be successfully called within a jitted function independently of whether interpolation options are used (such as xto.Linear), this does not work for multivariate output.

Here a simple example based on your quick guide:

import numpy as np
from numba import njit
from interpolation.splines import eval_linear, UCGrid, nodes
from interpolation.splines import extrap_options as xto

f = lambda x,y: np.sin(x**3+y**2+0.00001)/np.sqrt(x**2+y**2+0.00001)
g = lambda x,y: np.sin(x**3+y**2+0.00001)/np.sqrt(x**2+y**2+0.00001)

grid = UCGrid((-1.0, 1.0, 10), (-1.0, 1.0, 10))
gp = nodes(grid)   # 100x2 matrix

mvalues = np.concatenate([
   f(gp[:,0], gp[:,1]).reshape((10,10))[:,:,None],
   g(gp[:,0], gp[:,1]).reshape((10,10))[:,:,None]
],axis=2) # 10x10x2 array

points = np.random.random((1000,2))

@njit
def fun():
    eval_linear(grid, mvalues, points)

@njit
def no_fun():
    eval_linear(grid, mvalues, points, xto.LINEAR)

fun()       # works happily
no_fun()    # does not :'(

Environment is arch linux with: Python 3.7.4 numba 0.45 numpy 1.17 interpolation.py master branch

I tried to dig into your code but I admit that I got lost in codegen.py. I am happy about hints.

albop commented 5 years ago

@gboehl Thanks for your interest. The xto.LINEAR, xto.CONSTANT and other options are objects which are just designed so as to hold different types so multiple dispatch can occur. Initial idea was to define them as different jitted classes which are then automatically of a different type (see https://github.com/EconForge/interpolation.py/blob/master/interpolation/splines/option_types.py). But, as I just learnt, jitclasses don't seem to be unboxed correctly in njitted context, https://github.com/numba/numba/issues/2501 . I have pushed a temporary workaround where xto.LINEAR is just a random tuple value. The implementation is ugly but it seems to work. Could you try master a few minutes from now ?

gboehl commented 5 years ago

Thanks for the quick fix. It seems that you closed your PR instead of merging it. Anyways, I checked out albob/fix_xto and now using the xto.LINEAR option does not throw errors. Results seem reasonable.

Many thanks!

albop commented 5 years ago

Merged.

Le dim. 11 août 2019 à 02:11, Gregor Boehl notifications@github.com a écrit :

Thanks for the quick fix. It seems that you closed your PR instead of merging it. Anyways, I checked out albob/fix_xto and now using the xto.LINEAR option does not throw errors. Results seem reasonable.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/EconForge/interpolation.py/issues/52?email_source=notifications&email_token=AACDSKOJAL6UX6G52S4EM4DQD5KKPA5CNFSM4IKXPVCKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD4AXKTY#issuecomment-520189263, or mute the thread https://github.com/notifications/unsubscribe-auth/AACDSKODMJLGOSKXU5G3VQ3QD5KKPANCNFSM4IKXPVCA .

gboehl commented 5 years ago

Sorry, this is not a fix as it seems to be messing with the shape of the output of the interpolation. I'll have closer look next week.

albop commented 5 years ago

Strange, I would have assume output shape would be independent from the interpolation method. I don't have any obvious anomaly on my laptop. If you want to look I would suggest to start with splines/eval_splines.py . This is were the time dispatch is implemented. It is a bit heavy because of the need to overload function. Dispatch rules determine the source code of the actual spline interpolation logic. You can just look at the output of lines 102 and 127 to see exactly what is being generated (maybe we should have a debug option somewhere instead of commented print statements...)

gboehl commented 5 years ago

I can't reproduce this on my workstation with the latest versions of the related packages, so I assume that the error was somehow related to my python environment. Thanks a lot for the quick reaction!

(closing)