jakeret / hope

HOPE: A Python Just-In-Time compiler for astrophysical computations
GNU General Public License v3.0
382 stars 27 forks source link

Issue with numpy.interp #43

Closed herbelj closed 8 years ago

herbelj commented 8 years ago

I found an issue when using the function numpy.interp within a Python function compiled using HOPE. The issue is that interpolation produces different results depending on whether HOPE is used or not. The following example illustrates this.

    import numpy as np
    import hope
    import pickle

    def interpolation(x, xp, fp):
        intp = np.interp(x, xp, fp)
        return intp

    @hope.jit
    def interpolation_hope(x, xp, fp):
        intp = np.interp(x, xp, fp)
        return intp

    with open('np.interp_github_issue.pkl', 'r') as f:
        x, xp, fp = pickle.load(f)
    x = x.astype(np.float32)
    xp = xp.astype(np.float32)
    fp = fp.astype(np.float32)

    intp = interpolation(x, xp, fp)
    intp_hope = interpolation_hope(x, xp, fp)
    np.allclose(intp, intp_hope)

np.interp_github_issue.pkl.zip

cosmo-ethz commented 8 years ago

Thanks for reporting this. Seems like HOPE is not handling the boundaries correctly

import numpy as np
import hope

def interpolation(x, xp, fp):
    intp = np.interp(x, xp, fp)
    return intp
interpolation_hope = hope.jit(interpolation)

xp = np.arange(5,15)
fp = -xp
x = np.arange(0, 20)

intp = interpolation(x, xp, fp)
intp_hope = interpolation_hope(x, xp, fp)
np.allclose(intp, intp_hope)
cosmo-ethz commented 8 years ago

@herbelj I modified the bounds checking in the HOPE interp implementation. May I ask you to checkout the interp_bounds branch and quickly check if this resolves your problem? Thanks

herbelj commented 8 years ago

@cosmo-ethz The problem is fixed now, thanks! :)