EconForge / interpolation.py

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

Errors with multidimensional outputs #104

Open Mv77 opened 1 year ago

Mv77 commented 1 year ago

Hi @albop,

I have been working on enhancing our HARK wrapper of interpolation to work with functions of multiple outputs, $f:\mathbb{R}^n\rightarrow \mathbb{R}^m, m>1$. I have run into a strange behavior where calling eval_linear with parameters representing a function of multiple outputs can return interpolations of a single output. Say, I am approximating [g(x), h(x)] and it might only return g(x).

I created the following script that reproduces the error

from interpolation.splines import CGrid, eval_linear
from interpolation.splines import extrap_options as xto
import numpy as np

# Create grids
x = np.linspace(0, 10, 11)
y = np.linspace(0, 10, 11)
grid = CGrid(x, y)

# Make two linear functions
f1 = lambda x, y: 2 * x + y
f2 = lambda x, y: 3 * x + 2 * y

# Values for the two interpolated functions
x_tiled, y_tiled = np.meshgrid(x, y, indexing='ij')

z1 = f1(x_tiled, y_tiled)
z2 = f2(x_tiled, y_tiled)
z1z2 = np.stack([z1,z2], axis=-1)

# Evaluate always on the same points
npoints = 300
eval_points = np.column_stack([np.linspace(0, 10, npoints)]*2)

for i in range(1000):

    print(i)

    # Evaluate 1d interpolator
    if True:
        z1_eval = eval_linear(
            grid,
            z1,
            eval_points,
            xto.LINEAR,
        )
        assert(z1_eval.shape == (npoints,))

    # Evaluate 2d interpolator
    z1z2_eval = eval_linear(
        grid,
        z1z2,
        eval_points,
        xto.LINEAR,
    )
    assert(z1z2_eval.shape == (npoints,2))

    print('ok')

The script interpolates a 2-d function a thousand times and raises an error if the output is not of the expected shape.

Here is the output I get

...
ok
163
ok
164
ok
165
Traceback (most recent call last):
  File "/home/mvg/Dropbox/bug_example.py", line 49, in <module>
    assert(z1z2_eval.shape == (npoints,2))
AssertionError

There are various strange things about the behavior of this script:

This makes me think that what is going on is some bug in the overloading of functions with different inputs, but that only explains part of the strange behaviors above.

Do you know what might be going on? I'd be happy to help fix/test this issue.

albop commented 1 year ago

Hmm, it doesn't look immediately obvious. I need to think a little bit. Will get back to it later this week.

albop commented 1 year ago

@Mv77, I have, I believe found the source of the error and fixed it in master. Do you want to check ?

albop commented 1 year ago

Actually, I haven't. The bug is still there...

Mv77 commented 1 year ago

Hey @albop , thanks for looking into it.

I ran the test anyway and yes, it looks like the error is still there.