constant timing for convolution methods (fft and direct)

fftconv-conv-timings-2d.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import asarray, array, ndarray\n",
    "from scipy.signal import fftconvolve, correlate\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def _inputs_swap_needed(mode, shape1, shape2):\n",
    "    if mode == 'valid':\n",
    "        ok1, ok2 = True, True\n",
    "\n",
    "        for d1, d2 in zip(shape1, shape2):\n",
    "            if not d1 >= d2:\n",
    "                ok1 = False\n",
    "            if not d2 >= d1:\n",
    "                ok2 = False\n",
    "\n",
    "        if not (ok1 or ok2):\n",
    "            raise ValueError(\"For 'valid' mode, one must be at least \"\n",
    "                             \"as large as the other in every dimension\")\n",
    "\n",
    "        return not ok1\n",
    "\n",
    "    return False\n",
    "\n",
    "def _numeric_arrays(arrays, kinds='buifc'):\n",
    "    if type(arrays) == ndarray:\n",
    "        return arrays.dtype.kind in kinds\n",
    "    for array_ in arrays:\n",
    "        if array_.dtype.kind not in kinds:\n",
    "            return False\n",
    "    return True\n",
    "\n",
    "\n",
    "def _prod(iterable):\n",
    "    product = 1\n",
    "    for x in iterable:\n",
    "        product *= x\n",
    "    return product\n",
    "\n",
    "\n",
    "def _fftconv_faster(x, h, mode):\n",
    "    if mode == 'full':\n",
    "        small_dims = set([min(h.shape), min(x.shape)])\n",
    "        large_dims = set([max(h.shape), max(x.shape)])\n",
    "        if h.ndim == x.ndim == 2 \\\n",
    "        and (min(h.shape) < 50 or min(x.shape) < 50)\\\n",
    "        and max(large_dims) > 40:\n",
    "            # (40, 2*n), (2, n)\n",
    "            return True\n",
    "        if h.size <= 90 or x.size <= 90:  # orange lines\n",
    "            return False\n",
    "        if x.size*h.size < 600**2:  # green line\n",
    "            return False\n",
    "        out_shape = [n + k - 1 for n, k in zip(x.shape, h.shape)]\n",
    "    elif mode == 'same':\n",
    "        if h.size <= 100 or x.size <= 6:  # orange lines\n",
    "            return False\n",
    "        if x.size*h.size < 100**2:  # green line\n",
    "            return False\n",
    "        out_shape = x.shape\n",
    "    elif mode == 'valid':\n",
    "        shape_diff = _prod([n - k for n, k in zip(x.shape, h.shape)])\n",
    "        # this decision was calculated with a single dimension\n",
    "        shape_diff = np.power(np.abs(shape_diff), 1 / x.ndim)\n",
    "        if (x.ndim == 1 and shape_diff < 300) or \\\n",
    "           (x.ndim < 6 and shape_diff <= 2) or (shape_diff == 0):\n",
    "            return False\n",
    "        out_shape = [n - k + 1 for n, k in zip(x.shape, h.shape)]\n",
    "        if h.size <= 90 or x.size <= 90:  # orange lines\n",
    "            return False\n",
    "        if x.size <= 1600 and h.size <= 1600:  # cyan box\n",
    "            return False\n",
    "        if x.size == h.size:  # purple diagonal\n",
    "            return False\n",
    "        if x.size*h.size < 800**2:  # green line\n",
    "            return False\n",
    "    else:\n",
    "        raise ValueError('mode is invalid')\n",
    "\n",
    "    # see whether the Fourier transform convolution method or the direct\n",
    "    # convolution method is faster (discussed in scikit-image PR #1792)\n",
    "    big_O_constant = 5e-6 if x.ndim > 1 else 4.81e-4\n",
    "    direct_time = big_O_constant * (x.size * h.size * _prod(out_shape))\n",
    "    fft_time = np.sum(n * np.log(n) for n in (x.shape + h.shape +\n",
    "                                               tuple(out_shape)))\n",
    "    return fft_time < direct_time\n",
    "\n",
    "\n",
    "def _reverse_and_conj(x):\n",
    "    reverse = [slice(None, None, -1)] * x.ndim\n",
    "    return x[reverse].conj()\n",
    "\n",
    "\n",
    "def _np_conv_ok(volume, kernel, mode):\n",
    "    np_conv_ok = volume.ndim == kernel.ndim == 1\n",
    "    return np_conv_ok and (volume.size >= kernel.size or mode != 'same')\n",
    "\n",
    "\n",
    "def _choose_conv_method(volume, kernel, mode):\n",
    "    # fftconvolve doesn't support complex256\n",
    "    if hasattr(np, \"complex256\"):\n",
    "        if volume.dtype == 'complex256' or kernel.dtype == 'complex256':\n",
    "            return 'direct'\n",
    "\n",
    "    # for integer input,\n",
    "    # catch when more precision required than float provides (representing a\n",
    "    # integer as float can lose precision in fftconvolve if larger than 2**52)\n",
    "    if any([_numeric_arrays([x], kinds='ui') for x in [volume, kernel]]):\n",
    "        max_value = int(np.abs(volume).max()) * int(np.abs(kernel).max())\n",
    "        max_value *= int(min(volume.size, kernel.size))\n",
    "        if max_value > 2**np.finfo('float').nmant - 1:\n",
    "            return 'direct'\n",
    "\n",
    "    if _numeric_arrays([volume, kernel]) and _fftconv_faster(volume,\n",
    "                                                             kernel, mode):\n",
    "        return 'fft'\n",
    "    return 'direct'\n",
    "\n",
    "\n",
    "def convolve(in1, in2, mode='full', method='auto'):\n",
    "    volume = asarray(in1)\n",
    "    kernel = asarray(in2)\n",
    "\n",
    "    if volume.ndim == kernel.ndim == 0:\n",
    "        return volume * kernel\n",
    "\n",
    "    if _inputs_swap_needed(mode, volume.shape, kernel.shape):\n",
    "        # Convolution is commutative; order doesn't have any effect on output\n",
    "        volume, kernel = kernel, volume\n",
    "\n",
    "    if method == 'fft' and volume.dtype.kind not in 'buifc' \\\n",
    "                       and kernel.dtype.kind not in 'buifc':\n",
    "        raise ValueError(\"fftconvolve only supports numeric dtypes\")\n",
    "\n",
    "    if method == 'auto':\n",
    "        method = _choose_conv_method(volume, kernel, mode)\n",
    "\n",
    "    if method == 'fft':\n",
    "        out = fftconvolve(volume, kernel, mode=mode)\n",
    "        if volume.dtype.kind in 'ui':\n",
    "            out = np.around(out)\n",
    "        return out.astype(volume.dtype)\n",
    "\n",
    "    # fastpath to faster numpy convolve for 1d inputs\n",
    "    if _np_conv_ok(volume, kernel, mode):\n",
    "        return np.convolve(volume, kernel, mode)\n",
    "\n",
    "    return correlate(volume, _reverse_and_conj(kernel), mode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "randn = np.random.randn\n",
    "import timeit\n",
    "def time(x, h, mode='valid', method='auto', repeat=20, number=2):\n",
    "    t = timeit.repeat(\"convolve(x, h, mode='{}', method='{}')\".format(mode, method), \n",
    "                      \"import numpy as np\\n\" +\n",
    "                      \"from __main__ import convolve, x, h\", repeat=repeat, number=number)\n",
    "    return np.array(t).min() / number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "valid\n",
      "2\n",
      "3\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "8\n",
      "9\n",
      "11\n",
      "14\n",
      "17\n",
      "21\n",
      "25\n",
      "31\n",
      "37\n",
      "46\n",
      "55\n",
      "67\n",
      "82\n",
      "100\n"
     ]
    }
   ],
   "source": [
    "import pickle\n",
    "import pandas as pd\n",
    "\n",
    "ndims = 2\n",
    "N = 20  # we run N^2 test points\n",
    "Ns = {'valid': 2, 'same':2, 'full': 2}\n",
    "\n",
    "for mode in ['valid']:#, 'full', 'valid']:#, 'same', 'full']:\n",
    "    print(mode)\n",
    "    to_dump = {}\n",
    "    max_log_n = Ns[mode]\n",
    "    for n in np.logspace(0.4, max_log_n, num=N, dtype=int):\n",
    "        print(n)\n",
    "        m = n\n",
    "        for a in np.logspace(0.4, max_log_n, num=N, dtype=int):\n",
    "            b = a\n",
    "            x = randn(n, a)\n",
    "            h = randn(a, n)\n",
    "            if mode == 'valid':\n",
    "                x = randn(n, n//2)\n",
    "                h = randn(a, a//2)\n",
    "            t_fft = time(x, h, mode=mode, method='fft')\n",
    "            t_direct = time(x, h, mode=mode, method='direct')\n",
    "            to_dump[(n, m, a, b)] = np.log10(t_fft / t_direct)\n",
    "\n",
    "    filename = 'better-timings-2d/mode={}_log=True_N={}_max_log_n={}.pkl'.format(mode, N, max_log_n, ndims)\n",
    "    pickle.dump(to_dump, open(filename, 'wb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pickle\n",
    "import pandas as pd\n",
    "mode = 'same'\n",
    "N = 20\n",
    "max_log_n = 2\n",
    "ndims = 2\n",
    "filename = 'better-timings-2d/mode={}_log=True_N={}_max_log_n={}.pkl'.format(mode, N, max_log_n, ndims)\n",
    "d = pickle.load(open(filename, 'rb'))\n",
    "# if ndims == 2:\n",
    "#     d = {(n, n, m, m): d[(n, m)] for (n, m) in d}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal constant = [ 34519.21021589]\n",
      "0.00261603004747 2.66528811381 0.495785210527\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "def log_timend(sizes, O):\n",
    "    x_shape, h_shape = sizes[:, :2].copy(), sizes[:, 2:].copy()\n",
    "    out_shape = {'full': [n + k - 1 for n, k in zip(x_shape, h_shape)],\n",
    "                 'same': [n if n[0] > k[0] else k\n",
    "                              for n, k in zip(x_shape, h_shape)],\n",
    "                 'valid': [n-k+1 if n[0] > k[0] else k-n+1\n",
    "                              for n, k in zip(x_shape, h_shape)]}\n",
    "    out_shape = np.array(out_shape[mode])\n",
    "    shapes = np.hstack((x_shape, h_shape, out_shape)).T\n",
    "    \n",
    "    direct_time = np.product(shapes, axis=0)\n",
    "    fft_time = np.sum(shapes * np.log(shapes), axis=0)\n",
    "    return np.log10(O * fft_time / direct_time)\n",
    "\n",
    "x = np.array([[*sizes] for sizes in d])\n",
    "y = np.array([d[sizes] for sizes in d])\n",
    "t_ratio = np.exp(y)\n",
    "\n",
    "kwargs = {'sigma': 1 / weights} if False else {}\n",
    "# weights = np.abs(1 / (1 - t_ratio + 1e-3))\n",
    "# weights = np.power(weights, 1/5)\n",
    "# weights = np.ones_like(weights) if 'sigma' not in kwargs.keys() else weights\n",
    "\n",
    "y_hat = log_timend(x, 1)\n",
    "error = np.abs(y - y_hat)\n",
    "popt, pcov = curve_fit(log_timend, x, y, p0=[2000],#, 0, 0],#, [1, 0, 0, 0, 0, 0]],\n",
    "                       **kwargs, maxfev=int(1e4))\n",
    "\n",
    "print(\"Optimal constant = {}\".format(popt))\n",
    "\n",
    "y_hat = log_timend(x, *popt)\n",
    "error = np.abs(y - y_hat)\n",
    "print(error.min(), error.max(), np.median(error))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAECdJREFUeJzt3W2MXGd5xvH/5ThYhDSWW7AXNRAXqryAVJKqCkUgsVVI\nwVSt86EK0DcHhFRVIJAqVTioxf5UQb8gqn5CvMitoCWggk0JiYmcUUXbhEASESWOC60IL8KDKIQq\nTRWR5u6HHaNl8e6ceduZefz/SUc+c/acPfc+tq59fM85Z1JVSJLasGPeBUiSpsdQl6SGGOqS1BBD\nXZIaYqhLUkMMdUlqyNBQT3JlkgeS3D/480dJ3pFkT5KTSc4kuTPJ7u0oWJK0uYxynXqSHcC3gZcD\nbwf+q6r+Ksm7gD1VdXg2ZUqSuhi1/fIa4D+q6lvAQeDYYPsx4KZpFiZJGt2oof4G4OOD9X1V1Qeo\nqrPA3mkWJkkaXedQT3Ix8DvAJwebNvZtfN6AJM3ZzhH2PQB8paq+P3jdT7KvqvpJVoDvne+gJIa9\nJI2hqjLqMaO0X94E/P261yeAWwbrh4DjWxS2tMuRI0fmXsOFWv8y127981+Wvf5xdQr1JJew9ibp\nP67b/D7gxiRngBuA945dhSRpKjq1X6rqSeB5G7b9gLWglyQtCO8oHWJ1dXXeJUxkmetf5trB+udt\n2esf10g3H411gqRmfQ5Jak0SasZvlEqSFpyhLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtS\nQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXE\nUJekhnQK9SS7k3wyyekkDyd5eZI9SU4mOZPkziS7Z12sJGlrXWfqHwBur6prgJcBjwKHgbuq6irg\nFHDrLApcWdlPkrGXlZX9syhLkhZSqmrrHZLLgAeq6sUbtj8KvLqq+klWgF5VXX2e42vYOYacHxj/\neAiTnF+S5iEJVZVRj+syU/8l4PtJPprk/iQfTHIJsK+q+gBVdRbYO+rJJUnTtbPjPr8KvK2qvpzk\n/ay1XjZOfzedDh89evQn66urq6yuro5cqCS1rNfr0ev1Jv4+Xdov+4B/q6oXDV6/irVQfzGwuq79\ncveg577xeNsvkjSimbVfBi2WbyW5crDpBuBh4ARwy2DbIeD4qCeXJE3X0Jk6QJKXAR8CLgb+E3gz\ncBFwG/AC4DHg5qp6/DzHOlOXpBGNO1PvFOqTMNQlaXSzvPpFkrQkDHVJaoihLkkNMdQlqSGGuiQ1\nxFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMM\ndUlqiKEuSQ0x1CWpITu34yTvec+R7TiNJF3wUlWzPUFScHTMo78GfAyYpMYw659RkqYtCVWVkY/b\nnlAf9xx3AAcw1CVdaMYN9U7tlyTfAH4EPAP8uKquT7IH+ARwBfAN4Oaq+tGoBUiSpqfrG6XPAKtV\ndV1VXT/Ydhi4q6quAk4Bt86iwMntIsnYy8rK/nn/AJLUWddQz3n2PQgcG6wfA26aVlHT9RRr7Zvx\nln7/sTnULEnj6RrqBXwhyX1J3jrYtq+q+gBVdRbYO4sCJUnddb2k8ZVV9d0kzwNOJjnDz757ucW7\nkUfXra8OFknSOb1ej16vN/H3GfnqlyRHgCeAt7LWZ+8nWQHurqprzrP/3K9+8eoZSctm3KtfhrZf\nklyS5NLB+nOA3wQeAk4Atwx2OwQcH/XkkqTp6tJ+2Qd8em3GzU7gY1V1MsmXgduSvAV4DLh5hnVK\nkjq4IG4+sv0iadnMrP0iSVoehrokNcRQn7GVlf3e0Spp22zLo3cvZGt3pI7fk+/3R26pSbqAOVOX\npIYY6pLUEENdkhpiqEtSQ3yjdKi157FL0jIw1Ic69zz2cfkLQdL2sf0iSQ0x1CWpIYa6JDXEUJek\nhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUkM6hnmRHkvuTnBi83pPk\nZJIzSe5Msnt2ZUqSuhhlpv5O4JF1rw8Dd1XVVcAp4NZpFiZJGl2nUE9yOfB64EPrNh8Ejg3WjwE3\nTbc0SdKous7U3w/8GT/9aRH7qqoPUFVngb1Trk2SNKKhn3yU5LeAflU9mGR1i123+Higo+vWVweL\nJOmcXq9Hr9eb+PukauuPakvyl8AfAE8DzwZ+Dvg08GvAalX1k6wAd1fVNec5vsb/OLg7gANM/nFy\ny338sL8jSe1JQlWN/HmYQ9svVfXuqnphVb0IeCNwqqr+EPgscMtgt0PA8VFPLkmarkmuU38vcGOS\nM8ANg9eaul0kGXtZWdk/7x9A0jYa2n6Z+AS2X+Z+vO0bafnMrP0iSVoehrokNcRQl6SGGOqS1BBD\nXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQl\nqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWrI0FBPsivJvUkeSPJQkiOD7XuSnExyJsmdSXbP\nvlxJ0laGhnpVPQX8RlVdB1wLHEhyPXAYuKuqrgJOAbfOtFJJ0lCd2i9V9eRgdRewEyjgIHBssP0Y\ncNPUq5MkjaRTqCfZkeQB4Czwhaq6D9hXVX2AqjoL7J1dmZKkLnZ22amqngGuS3IZ8OkkL2Vttv5T\nu23+HY6uW18dLJKkc3q9Hr1eb+Lvk6otsvh8ByR/ATwJvBVYrap+khXg7qq65jz715Z5v6U7gAOM\nfzxALvjjR/07ljR/SaiqjHpcl6tfnnvuypYkzwZuBE4DJ4BbBrsdAo6PenJJ0nR1ab88HziWZAdr\nvwQ+UVW3J7kHuC3JW4DHgJtnWKckqYOR2y8jn8D2y9yPt/0iLZ+ZtV8kScvDUJekhhjqktQQQ12S\nGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakh\nhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0ZGupJLk9yKsnDSR5K8o7B9j1JTiY5k+TO\nJLtnX64kaStdZupPA39aVS8FXgG8LcnVwGHgrqq6CjgF3Dq7MiVJXQwN9ao6W1UPDtafAE4DlwMH\ngWOD3Y4BN82qSElSNyP11JPsB64F7gH2VVUf1oIf2Dvt4iRJo9nZdccklwKfAt5ZVU8kqQ27bHy9\nztF166uDRZJ0Tq/Xo9frTfx9UrVFFp/bKdkJ/BPw+ar6wGDbaWC1qvpJVoC7q+qa8xxbW+b9lu4A\nDjD+8QC54I/v8ncsabEkoaoy6nFd2y8fAR45F+gDJ4BbBuuHgOOjnlySNF1DZ+pJXgn8M/AQa1PG\nAt4NfAm4DXgB8Bhwc1U9fp7jnak7U5c0onFn6kN76lX1L8BFm3z5NaOeUJI0O95RKkkNMdQlqSGG\nuiQ1xFCXpIYY6pLUEEO9ebtIMvaysrJ/3j+ApBF0fkyAltVTTHKde78/8mWykubImbokNcRQl6SG\nGOqS1BBDXZIaYqhriPGvnvHKGWn7efWLhhj/6hmvnJG2nzN1SWqIoS5JDTHUJakhhrokNcRQl6SG\nGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0ZGupJPpykn+Sr67btSXIyyZkkdybZPdsyJUlddJmpfxR4\n7YZth4G7quoq4BRw67QLkySNbmioV9UXgR9u2HwQODZYPwbcNOW6JEljGLenvreq+gBVdRbYO72S\nJEnjmtajd4c8m/XouvXVwSJJOqfX69Hr9Sb+Pqka/qzsJFcAn62qXxm8Pg2sVlU/yQpwd1Vds8mx\nNf6n2d8BHGD84wHi8XM7PnT59yXpZyWhqkb+UIKu7ZcMlnNOALcM1g8Bx0c9sSRp+rpc0vhx4F+B\nK5N8M8mbgfcCNyY5A9wweC1JmrOhPfWq+r1NvvSaKdciSZqQd5RKUkMMdUlqiKEuSQ0x1CWpIYa6\nJDXEUJekhhjqktQQQ12SGmKoq1krK/tJMvaysrJ/3j+CNLJpPaVRWjj9/mNM8jCzfn/kZylJc+dM\nXZIaYqhrYU3aPpEuRLZftLAmbZ/89NOipQuDM3VJaoihLs3IpO2jiy56jlfvaGS2X6QZmbR99Mwz\nk30UoVfvXJicqUtSQwx1zdCuJb96Zdnrn69J2k+2jsZn+0Uz9BTLffXKstc/X5O0n2wdjc+ZuiQ1\nxFCXmjVZ+8gWyHKy/SI1a7L2kS2Q5TTRTD3J65I8muTfk7xrWkVJksYzdqgn2QH8DfBa4KXAm5Jc\nPa3CJGkSvV5v3iXMxSQz9euBr1XVY1X1Y+AfgIPTKUvShW3y9wMmCfVlfhb/JD31XwS+te71t1kL\nekma0HzfD1jmZ/Fvyxull13222Md9/TT3+PJJ6dcjCQ1LFXj/TZK8uvA0ap63eD1YaCq6n0b9pvk\n7g1JumBV1chT/klC/SLgDHAD8F3gS8Cbqur0WN9QkjSxsdsvVfV/Sd4OnGTtDdcPG+iSNF9jz9Ql\nSYtnao8J6HIjUpK/TvK1JA8muXZa557UsNqTvDrJ40nuHyx/Po86N5Pkw0n6Sb66xT4LOfYwvP5F\nHv8klyc5leThJA8leccm+y3k+Hepf8HHf1eSe5M8MKj/yCb7Ldz4d6l9rLGvqokX1n45fB24ArgY\neBC4esM+B4DPDdZfDtwzjXNvU+2vBk7Mu9YtfoZXAdcCX93k6ws59iPUv7DjD6wA1w7WL2Xtfaal\n+Lc/Qv0LO/6D+i4Z/HkRcA9w/RKN/7DaRx77ac3Uu9yIdBD4W4CquhfYnWTflM4/ia43US3sgzCq\n6ovAD7fYZVHHHuhUPyzo+FfV2ap6cLD+BHCatXs41lvY8e9YPyzo+ANU1bkLn3ex9j7hxp7yIo//\nsNphxLGfVqif70akjf8wNu7znfPsMw9dagd4xeC/bp9L8pLtKW1qFnXsR7Hw459kP2v/47h3w5eW\nYvy3qB8WePyT7EjyAHAW+EJV3bdhl4Ud/w61w4hj71Mau/kK8MKqejLJAeAzwJVzrulCsvDjn+RS\n4FPAOwcz3qUypP6FHv+qega4LsllwGeSvKSqHpl3XV10qH3ksZ/WTP07wAvXvb58sG3jPi8Yss88\nDK29qp4499+kqvo8cHGSn9++Eie2qGPfyaKPf5KdrAXi31XV8fPsstDjP6z+RR//c6rqv4G7gddt\n+NJCjz9sXvs4Yz+tUL8P+OUkVyR5FvBG4MSGfU4AfwQ/uRv18arqT+n8kxha+/r+W5LrWbsU9Afb\nW+ZQYfPe26KO/Xqb1r8E4/8R4JGq+sAmX1/08d+y/kUe/yTPTbJ7sP5s4Ebg0Q27LeT4d6l9nLGf\nSvulNrkRKckfr325PlhVtyd5fZKvA/8DvHka555Ul9qB303yJ8CPgf8F3jC/in9Wko8Dq8AvJPkm\ncAR4Fgs+9ucMq58FHv8krwR+H3ho0Bst4N2sXU218OPfpX4WePyB5wPHsvYo8B3AJwbjvfDZQ4fa\nGWPsvflIkhriZ5RKUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12SGvL/isoDKwtRfjsA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11b396080>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "error.sort()\n",
    "plt.figure()\n",
    "plt.hist(error, bins=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}

This notebook re-implements a lot of SciPy's convolve and correlate. That can probably be ripped out now.

This notebook also measures the average time, not the minimum time.