Fixed#938: Speed up njit-decorated function for sliding dot product [WIP]

docs/OTFFT.ipynb

@@ -0,0 +1,178 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "5b88fa48",
+   "metadata": {},
+   "source": [
+    "This notebook implements the 6-step / 8-step FFT (fft) algorithm as provided in [OTFFT](http://wwwa.pikara.ne.jp/okojisan/otfft-en/stockham2.html). Accordingly, the inverse FFT (ifft) algorithm will be implemented as well. When the input of FFT, `T` with length `n`, consists of real-valued elements only, we can take advantage of real FFT (rfft) as explained in [2.6.2](https://www.researchgate.net/profile/Christos-Bechlioulis/publication/341270520_FFT_algorithms_are_not_mine_However_I_am_going_to_convince_you_soon_regarding_the_visit_of_RMS_to_our_university_Believe_it_or_not_this_is_me_This_is_us_Univeristy_of_Patras_you_have_chosen_a_quite_wr/links/5fa53ce7299bf10f7328c33b/FFT-algorithms-are-not-mine-However-I-am-going-to-convince-you-soon-regarding-the-visit-of-RMS-to-our-university-Believe-it-or-not-this-is-me-This-is-us-Univeristy-of-Patras-you-have-chosen-a-quite.pdf). "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "c9a886c2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import math\n",
+    "import time\n",
+    "\n",
+    "import numba\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import scipy\n",
+    "\n",
+    "from numba import njit, prange\n",
+    "import numpy.testing as npt\n",
+    "\n",
+    "from stumpy import core"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "80765bca",
+   "metadata": {},
+   "source": [
+    "Let's start with `rfft`. First, we write the test!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "7bb4242c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def test_rfft(n_powers_list):\n",
+    "    seed = 0\n",
+    "    np.random.seed(seed)\n",
+    "    for p in n_powers_list:\n",
+    "        n = 2 ** p\n",
+    "        T = np.random.rand(n)\n",
+    "        \n",
+    "        ref = scipy.fft.rfft(T)\n",
+    "        comp = rfft(T)\n",
+    "        \n",
+    "        npt.assert_almost_equal(ref, comp)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5a1c553e",
+   "metadata": {},
+   "source": [
+    "We now implement `rfft` function according to the steps provided in [2.6.2](https://www.researchgate.net/profile/Christos-Bechlioulis/publication/341270520_FFT_algorithms_are_not_mine_However_I_am_going_to_convince_you_soon_regarding_the_visit_of_RMS_to_our_university_Believe_it_or_not_this_is_me_This_is_us_Univeristy_of_Patras_you_have_chosen_a_quite_wr/links/5fa53ce7299bf10f7328c33b/FFT-algorithms-are-not-mine-However-I-am-going-to-convince-you-soon-regarding-the-visit-of-RMS-to-our-university-Believe-it-or-not-this-is-me-This-is-us-Univeristy-of-Patras-you-have-chosen-a-quite.pdf)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "d0033890",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def _rfft(T):\n",
+    "    \"\"\"\n",
+    "    For the input `T` with length `n=len(T)`, this function returns its\n",
+    "    real fast fourier transform (rfft) with length of `(n // 2) + 1`.\n",
+    "    \n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    T : numpy.ndarray\n",
+    "        A time series of interest, with real-valued numbers\n",
+    "    \n",
+    "    Returns\n",
+    "    -------\n",
+    "    out : numpy.ndarray\n",
+    "        the real fast fourier transform (rfft) of input `T`\n",
+    "    \"\"\"\n",
+    "    n = len(T)\n",
+    "    half_n = int(n // 2)\n",
+    "    \n",
+    "    x = T[::2] + 1j * T[1::2]\n",
+    "    x[:] = scipy.fft.fft(x)  # we will implement our fft shortly!\n",
+    "    \n",
+    "    out = np.empty(half_n + 1, dtype=np.complex_)\n",
+    "    out[0] = x[0].real + x[0].imag\n",
+    "    out[half_n] = x[0].real - x[0].imag\n",
+    "    out[n // 4] = x[n // 4].conjugate()\n",
+    "    \n",
+    "    theta0 = 2 * math.pi / n\n",
+    "    for k in range(1, n // 4):\n",
+    "        theta = theta0 * k\n",
+    "        a =  x[half_n - k].conjugate()\n",
+    "        b = 0.5 * (x[k] - a) * (1.0 + complex(math.sin(theta), math.cos(theta)))\n",
+    "        out[k] = x[k] - b\n",
+    "        out[half_n - k] = (a + b).conjugate()\n",
+    "    \n",
+    "    return out"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "d5db0eb9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def rfft(T):\n",
+    "    \"\"\"\n",
+    "    For the input `T` with length `n=len(T)`, this function returns its\n",
+    "    real fast fourier transform (rfft) with length of `(n // 2) + 1`.\n",
+    "    \n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    T : numpy.ndarray\n",
+    "        A time series of interest, with real-valued numbers\n",
+    "    \n",
+    "    Returns\n",
+    "    -------\n",
+    "    out : numpy.ndarray\n",
+    "        the real fast fourier transform (rfft) of input `T`\n",
+    "    \"\"\"\n",
+    "    return _rfft(T)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "ded21f89",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n_powers_list = np.arange(2, 11)\n",
+    "test_rfft(n_powers_list)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "335ddf79",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.10.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}