scipy vs numba interp1d.ipynb

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": []
    }
  ]
}

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.