Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
scipy vs numba interp1d.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "scipy vs numba interp1d.ipynb",
"version": "0.3.2",
"views": {},
"default_view": {},
"provenance": [],
"collapsed_sections": []
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
}
},
"cells": [
{
"metadata": {
"id": "f97hCEt8qUWG",
"colab_type": "code",
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
},
"base_uri": "https://localhost:8080/",
"height": 68
},
"outputId": "99ab9033-4800-4d50-9319-482b81e787a6",
"executionInfo": {
"status": "ok",
"timestamp": 1528420189439,
"user_tz": 420,
"elapsed": 2477,
"user": {
"displayName": "Stephan Hoyer",
"photoUrl": "//lh4.googleusercontent.com/-bwQVXpRw0z8/AAAAAAAAAAI/AAAAAAAAACQ/obT9z9YnNnc/s50-c-k-no/photo.jpg",
"userId": "100105766565685654482"
}
}
},
"cell_type": "code",
"source": [
"! pip install numba"
],
"execution_count": 29,
"outputs": [
{
"output_type": "stream",
"text": [
"Requirement already satisfied: numba in /usr/local/lib/python3.6/dist-packages (0.38.1)\r\n",
"Requirement already satisfied: llvmlite>=0.23.0dev0 in /usr/local/lib/python3.6/dist-packages (from numba) (0.23.2)\r\n",
"Requirement already satisfied: numpy in /usr/local/lib/python3.6/dist-packages (from numba) (1.14.3)\r\n"
],
"name": "stdout"
}
]
},
{
"metadata": {
"id": "8JxQ0wO8qrZ7",
"colab_type": "code",
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
}
}
},
"cell_type": "code",
"source": [
"from numba import guvectorize\n",
"import numpy as np\n",
"import scipy.interpolate\n",
"\n",
"\n",
"def _interp1d(xnew, xvals, yvals, ynew):\n",
" i = 0\n",
" N = len(xvals)\n",
" if xnew[0] < xvals[0]:\n",
" x_a = 0.0\n",
" y_a = 0.0\n",
" x_b = xvals[0]\n",
" y_b = yvals[0]\n",
" else:\n",
" while xnew[0] >= xvals[i] and i < N:\n",
" i += 1\n",
" if xnew[0] == xvals[i]:\n",
" ynew[0] = yvals[i]\n",
" return\n",
" if i == N:\n",
" i = N-1\n",
" x_a = xvals[i-1]\n",
" y_a = yvals[i-1]\n",
" x_b = xvals[i]\n",
" y_b = yvals[i]\n",
" slope = (xnew[0] - x_a)/(x_b - x_a)\n",
" ynew[0] = slope * (y_b-y_a) + y_a\n",
" return\n",
"\n",
" \n",
"interp1d_numba = guvectorize(\n",
" ['float64[:], float64[:], float64[:], float64[:]'],\n",
" \"(),(n),(n) -> ()\", nopython=True)(_interp1d)\n",
"\n",
"\n",
"def _interp1d_np(xnew, xvals, yvals):\n",
" xnew = np.asarray(xnew)[..., np.newaxis]\n",
" out = np.empty_like(xnew)\n",
" _interp1d(xnew, xvals, yvals, out)\n",
" return out.squeeze(axis=-1)\n",
"\n",
"interp1d_numpy = np.vectorize(_interp1d_np, signature='(),(n),(n)->()')"
],
"execution_count": 0,
"outputs": []
},
{
"metadata": {
"id": "nqHSTRSHqxR5",
"colab_type": "code",
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
}
}
},
"cell_type": "code",
"source": [
"xnew = np.arange(900) / 900\n",
"xvals = yvals = np.arange(1000.0) / 1000"
],
"execution_count": 0,
"outputs": []
},
{
"metadata": {
"id": "FtlNs5Qorxzg",
"colab_type": "code",
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
}
}
},
"cell_type": "code",
"source": [
"np.testing.assert_allclose(interp1d_numpy(xnew, xvals, yvals), scipy.interpolate.interp1d(xvals, yvals)(xnew))\n",
"np.testing.assert_allclose(interp1d_numba(xnew, xvals, yvals), scipy.interpolate.interp1d(xvals, yvals)(xnew))"
],
"execution_count": 0,
"outputs": []
},
{
"metadata": {
"id": "K1HvMlcQq8sh",
"colab_type": "code",
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
},
"base_uri": "https://localhost:8080/",
"height": 34
},
"outputId": "f5f42e58-899b-46f6-b011-e24cf7ca75d2",
"executionInfo": {
"status": "ok",
"timestamp": 1528420813197,
"user_tz": 420,
"elapsed": 6096,
"user": {
"displayName": "Stephan Hoyer",
"photoUrl": "//lh4.googleusercontent.com/-bwQVXpRw0z8/AAAAAAAAAAI/AAAAAAAAACQ/obT9z9YnNnc/s50-c-k-no/photo.jpg",
"userId": "100105766565685654482"
}
}
},
"cell_type": "code",
"source": [
"%timeit interp1d_numpy(xnew, xvals, yvals)"
],
"execution_count": 67,
"outputs": [
{
"output_type": "stream",
"text": [
"10 loops, best of 3: 133 ms per loop\n"
],
"name": "stdout"
}
]
},
{
"metadata": {
"id": "PjXt7_Ay87Y2",
"colab_type": "code",
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
},
"base_uri": "https://localhost:8080/",
"height": 51
},
"outputId": "1525d6fa-9e04-4d7a-c37e-3708f2fa8d3a",
"executionInfo": {
"status": "ok",
"timestamp": 1528420816086,
"user_tz": 420,
"elapsed": 2095,
"user": {
"displayName": "Stephan Hoyer",
"photoUrl": "//lh4.googleusercontent.com/-bwQVXpRw0z8/AAAAAAAAAAI/AAAAAAAAACQ/obT9z9YnNnc/s50-c-k-no/photo.jpg",
"userId": "100105766565685654482"
}
}
},
"cell_type": "code",
"source": [
"%timeit interp1d_numba(xnew, xvals, yvals)"
],
"execution_count": 68,
"outputs": [
{
"output_type": "stream",
"text": [
"The slowest run took 13.55 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
"1000 loops, best of 3: 384 µs per loop\n"
],
"name": "stdout"
}
]
},
{
"metadata": {
"id": "pNZz374orr6f",
"colab_type": "code",
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
},
"base_uri": "https://localhost:8080/",
"height": 51
},
"outputId": "8bf28ed5-2219-4cd4-f7e8-044e61da766f",
"executionInfo": {
"status": "ok",
"timestamp": 1528420823759,
"user_tz": 420,
"elapsed": 4294,
"user": {
"displayName": "Stephan Hoyer",
"photoUrl": "//lh4.googleusercontent.com/-bwQVXpRw0z8/AAAAAAAAAAI/AAAAAAAAACQ/obT9z9YnNnc/s50-c-k-no/photo.jpg",
"userId": "100105766565685654482"
}
}
},
"cell_type": "code",
"source": [
"%timeit scipy.interpolate.interp1d(xvals, yvals)(xnew)"
],
"execution_count": 69,
"outputs": [
{
"output_type": "stream",
"text": [
"The slowest run took 51.79 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
"10000 loops, best of 3: 89.1 µs per loop\n"
],
"name": "stdout"
}
]
},
{
"metadata": {
"id": "pvdv6DYA-kuw",
"colab_type": "code",
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
}
}
},
"cell_type": "code",
"source": [
""
],
"execution_count": 0,
"outputs": []
}
]
}
@teoliphant

This comment has been minimized.

Copy link

teoliphant commented Jun 8, 2018

This is very helpful. Definitely binary-search would be much better for large N. In my case, N was about 10 with most interpolations being done in the first 2-3 and so binary-search didn't really help. But, in general, of course it's much better.

In my use-case, also scipy.interpolate.interp1d doesn't do what I needed for the situation where yvals.shape = (len(xnew), len(xvals)) with both xvals and xnew being 1-d. I needed the output to be 1-d with length the same as xnew. But, the output of scipy.interpolate.interp1d would have shape (yvals.shape[0], len(xnew)) --- i.e. an cross-product instead of a "zip" operation.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.