{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "4ec40081b048ce2f34f3f4fedbb0be10",
     "grade": false,
     "grade_id": "cell-98f724ece1aacb67",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "# CDS: Numerical Methods Assignments\n",
    "\n",
    "- See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.\n",
    "\n",
    "- Solutions must be submitted via the Jupyter Hub.\n",
    "\n",
    "- Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\".\n",
    "\n",
    "## Submission\n",
    "\n",
    "1. Name all team members in the the cell below\n",
    "2. make sure everything runs as expected\n",
    "3. **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart)\n",
    "4. **run all cells** (in the menubar, select Cell$\\rightarrow$Run All)\n",
    "5. Check all outputs (Out[\\*]) for errors and **resolve them if necessary**\n",
    "6. submit your solutions  **in time (before the deadline)**"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "team_members = \"Koen Vendrig, Kees van Kempen\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "2d40509e944f45e5c23cc3e732547aad",
     "grade": false,
     "grade_id": "cell-2a0ab25436430f05",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "## Discrete and Fast Fourier Transforms (DFT and FFT)\n",
    "\n",
    "In the following we will implement a DFT algorithm and, based on that, a FFT algorithm. Our aim is to experience the drastic improvement of computational time in the FFT case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "90044d2a3233721361765db06afba03b",
     "grade": true,
     "grade_id": "cell-abee6fbaf30772f2",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "bbe05b6dd0bc5b479cf66199114d7e4d",
     "grade": false,
     "grade_id": "cell-a1c3327dc1cad0db",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 1\n",
    "\n",
    "Implement a Python function $\\text{DFT(yk)}$ which returns the Fourier transform defined by\n",
    "\n",
    "\\begin{equation}\n",
    "\\beta_j = \\sum^{N-1}_{k=0} f(x_k) e^{-ij x_k}\n",
    "\\end{equation}\n",
    "\n",
    "with $x_k = \\frac{2\\pi k}{N}$ and $j = 0, 1, ..., N-1$. The $\\text{yk}$ should represent the array corresponding to $y_k = f(x_k)$. Please note that this definition is slightly different to the one we introduced in the lecture. Here we follow the notation of Numpy and Scipy.\n",
    "\n",
    "Hint: try to write the sum as a matrix-vector product and use $\\text{numpy.dot()}$ to evaluate it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "2fded00d23ce12e99ba32c09a6370b3c",
     "grade": true,
     "grade_id": "cell-5f6638846212c9d1",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def DFT(yk):\n",
    "    N = len(yk)\n",
    "    x = 2*np.pi*np.arange(N)/N\n",
    "    beta = yk.T * np.exp(-np.dot(np.arange(N), x)*1j)\n",
    "    return beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.-8.57252759e-16j, -2.-1.71450552e-15j, -7.-6.00076932e-15j,\n",
       "       -6.-5.14351656e-15j])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DFT(np.array([1,2,7,6]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "15164465870c44b9b4e6328f56eaed20",
     "grade": false,
     "grade_id": "cell-74e9ce917ff9d690",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 2 \n",
    "\n",
    "Make sure your function $\\text{DFT(yk)}$ and Numpy’s FFT function $\\text{numpy.fft.fft(yk)}$ return\n",
    "the same data by plotting $|\\beta_j|$ vs. $j$ for\n",
    "\n",
    "\\begin{equation}\n",
    "    y_k = f(x_k) = e^{20i x_k} + e^{40 i x_k}\n",
    "\\end{equation}\n",
    "and\n",
    "\\begin{equation}\n",
    "    y_k = f(x_k) = e^{i 5 x_k^2}\n",
    "\\end{equation}\n",
    "\n",
    "using $N = 128$ for both routines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "14cef25098dec15c058916f4fc58c133",
     "grade": true,
     "grade_id": "cell-7cc28776346ee714",
     "locked": false,
     "points": 1,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "1674c10abab84f4d7c15bc7e5fea53b2",
     "grade": false,
     "grade_id": "cell-e04e5bc7ed412f64",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 3\n",
    "\n",
    "Analyze the evaluation-time scaling of your $\\text{DFT(yk)}$ function with the help of the timeit\n",
    "module. Base your code on the following example:\n",
    "\n",
    "```python\n",
    "import timeit\n",
    "\n",
    "tOut = timeit.repeat(stmt=lambda: DFT(yk), number=10, repeat=5)\n",
    "tMean = np.mean(tOut)\n",
    "```\n",
    "This example evaluates $\\text{DFT(yk)}$ 5 × 10 times and stores the resulting 5 evaluation times in tOut. Afterwards we calculate the mean value of these 5 repetitions. \n",
    "Use this example to calculate and plot the evaluation time of your $\\text{DFT(yk)}$ function for $N = 2^2, 2^3, ..., 2^M$. Depending on your implementation you might be able to go up to $M = 10$. Be careful and increase M just step by step!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "0e92c1c2d548a1c9b1476dd295415c2e",
     "grade": true,
     "grade_id": "cell-0ab81532ab86e322",
     "locked": false,
     "points": 4,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "3063d5b408e133b1930c56fa45fed54e",
     "grade": false,
     "grade_id": "cell-73ca4de972164356",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 4\n",
    "\n",
    "A very simple FFT algorithm can be derived by the following separation of the sum from\n",
    "above:\n",
    "\n",
    "\\begin{align}\n",
    "     \\beta_j = \\sum^{N-1}_{k=0} f(x_k) e^{-ij \\frac{2\\pi k}{N}} \n",
    "         &= \\sum^{N/2 - 1}_{k=0} f(x_{2k}) e^{-ij \\frac{2\\pi 2k}{N}} \n",
    "          + \\sum^{N/2 - 1}_{k=0} f(x_{2k+1}) e^{-ij \\frac{2\\pi (2k+1)}{N}}\\\\ \n",
    "         &= \\sum^{N/2 - 1}_{k=0} f(x_{2k}) e^{-ij \\frac{2\\pi k}{N/2}}\n",
    "          + \\sum^{N/2 - 1}_{k=0} f(x_{2k+1}) e^{-ij \\frac{2\\pi k}{N/2}} e^{-ij \\frac{2\\pi}{N}}\\\\\n",
    "         &= \\beta^{\\text{even}}_j + \\beta^{\\text{odd}}_j e^{-ij \\frac{2\\pi}{N}}\n",
    "\\end{align}\n",
    "\n",
    "where $\\beta^{\\text{even}}_j$ is the Fourier transform based on only even $k$ (or $x_k$) and $\\beta^{\\text{odd}}_j$ the Fourier transform based on only odd $k$. In case $N = 2^M$ this even/odd separation can be done again and again in a recursive way. \n",
    "\n",
    "Use the template below to implement a $\\text{FFT(yk)}$ function, making use of your $\\text{DFT(yk)}$ function from above. Make sure that you get the same results as before by comparing the results from $\\text{DFT(yk)}$\n",
    "and $\\text{FFT(yk)}$ for both functions defined in task 2.\n",
    "\n",
    "```python\n",
    "def FFT(yk):\n",
    "    \"\"\"Don't forget to write a docstring ...\n",
    "    \"\"\"\n",
    "    N = # ... get the length of yk\n",
    "    \n",
    "    assert # ... check if N is a power of 2. Hint: use the % (modulo) operator\n",
    "    \n",
    "    if(N <= 2):\n",
    "        return # ... call DFT with all yk points\n",
    "    \n",
    "    else:\n",
    "        betaEven = # ... call FFT but using just even yk points\n",
    "        betaOdd = # ... call FFT but using just odd yk points\n",
    "        \n",
    "        expTerms = np.exp(-1j * 2.0 * np.pi * np.arange(N) / N)\n",
    "        \n",
    "        # Remember : beta_j is periodic in j !\n",
    "        betaEvenFull = np.concatenate([betaEven, betaEven])\n",
    "        betaOddFull = np.concatenate([betaOdd, betaOdd])\n",
    "        \n",
    "        return betaEvenFull + expTerms * betaOddFull\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "fdad5c990a077e2bc2c58d4c4da46973",
     "grade": true,
     "grade_id": "cell-ce8233802d8ccb83",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def FFT(yk):\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "33b415b478341202e57888621063fc35",
     "grade": true,
     "grade_id": "cell-a27cf0cb147b31ae",
     "locked": false,
     "points": 1,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# Do your plotting here ...\n",
    "\n",
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "c1d2dd0c2543e228d48d84e94720a3a8",
     "grade": false,
     "grade_id": "cell-c61537808f7cce68",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 5\n",
    "\n",
    "Analyze the evaluation-time scaling of your $\\text{FFT(yk)}$ function with the help of the timeit module and compare it to the scaling of the $\\text{DFT(yk)}$ function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "bb4f1eee977aad469245f60fca596938",
     "grade": true,
     "grade_id": "cell-aaf90559928426bf",
     "locked": false,
     "points": 1,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}