522 lines
37 KiB
Plaintext
522 lines
37 KiB
Plaintext
{
|
||
"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": 1,
|
||
"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\n",
|
||
"from matplotlib import pyplot as plt\n",
|
||
"import timeit"
|
||
]
|
||
},
|
||
{
|
||
"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": 2,
|
||
"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",
|
||
" xk = 2*np.pi*np.arange(N)/N\n",
|
||
" beta = np.dot(yk, np.exp(np.outer(-np.arange(N), xk*1j)))\n",
|
||
" return beta"
|
||
]
|
||
},
|
||
{
|
||
"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": 3,
|
||
"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": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA6UAAAGICAYAAAC5jEaHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABAAklEQVR4nO3de5SlZXXv+998VwGtiMilIYbWNAYGcovdpuPG42WAGgMGheyNCBElXsLhDN0KbqMQo8Zztjl6dBjjIGKIuIGEKAbdigx1S1S2aEBshKEIEi5CbIPSoIgXbr3eef54L+tdq1ZVra56q9/1zPX9jNGjq1atqnoeVlOz5vM8cz7m7gIAAAAAoAtZ1wMAAAAAAMwuklIAAAAAQGdISgEAAAAAnSEpBQAAAAB0hqQUAAAAANAZklIAAAAAQGdISgEAAAAAnSEpBQAAAAB0hqQU2IHM7BlmdrWZfc3MPm5mO5WPv9fMrjKzf2g89nozO7DbERfM7GQz2zry2FSPGQCAaVLFSDPbamZXln/WLvL8F5rZW83s76o4C0RFUgrsWD+U9Dx3f66kOyUdZ2ZPk7Sfuz9H0vclnSBJ7n6Ou9/a2UhLZtaT9FIVY68em+oxAwAwbdz9HEmPSvrf7n5k+WfrIs//kru/V9KvJO28o8YJdIGkFNiB3P1ud3+wfPcRSbmk/0PSl8rHvijpWZJkZleWf+9sZtea2Roz+yMzO3cHD/tkSf9cjrWy4JinYLwAAEydKq5LelZ50uivrLBg3DSz0yV9yd1/1cmggR2EpBTogJn9lqQXSvqcpD0kPVB+6OeS9jSzvSXdI0nu/oik8yV9SNIrJP3XHTjOnqQTJV0y8qEFx9zleAEAmEaNuH63pAMkPVfSPpL+80Jxs0xIXyDpqWa2RxfjBnaUua4HAMwaM3u8pH+Q9Cfu/qiZ3S/p8eWHd5f0U0m/I+m7jU/7V0kfVnFkdtt2fK+/lCR3/8tlDvcUSZ9099zMmo8vNeZljRcAgKB+R9J33f1hSQ9Lkpl9WtIRkj6lMXHT3T8i6SPdDBfYsdgpRXhm9jgz65vZExuPHWZmd5vZbqv0Pf/UzG4ys5+b2RfMbJ/y8TlJn5D0Lne/pXz6v6pYCZWkP5D0DTUSvHJ19YOS3q4iSZSZ7Wpmf74DxnyIpFea2RclHWhmH1pqzDtqvAAALGaa4r8GMbL5fZ8j6bZxcROYNSSlCM/df6miGc/TGw+/R9JfufsvFvtcM7vczO5f4M/lC3zOn0s6XdJLJK2V9CNJ/7388MmS/pOkt5f1ly9z9xsk/cTMrpJ0qIoV08MlfcfMdpF0kaQzJL1P0svMbFdJvyfp383sbWa2fnv/m0w6Znd/q7u/0N2PlnSru79BkhYZ8807YrwAACxlyuL/4ZK+I+nZZnZdGT/3U1EeMy9uspiLWWPu3vUYgFVnZhdIut3d/x8ze66kCyQ91d0fMbO/lnSRu1/fwvfZR9IPJG10938rH3umpA+7+8aVfv3G93mLiiM/Jy52PHaS47s7YsyTjhcAgDalGv/N7EhJ6yT9lqSL3f3OlY4RmGbUlGJWfEuD46b/n6S3l40FJOlgFTt8bXi+pF0kXduowTRJKw54I/aQ9FZJZ2mwClt8s2IF99nlu2vKx84o3/+6ux/bwZgXHC8AAKso1fj/DLGYixlCUopZ8S1JbzGz/6IiUfunxsce4+4PmdlZknaV9A4vjxCY2RdU1HyMc5W7HzPy2J6S/qe7v7Td4Q9z97PLN+cleM2kc8JGR6s+5sXGCwDAKko1/rOYi5nC8V3MhLI28wEVrdj/L3f/Qvn4PpI+JmmrpCvc/Z8W/ioTfZ9nqrjm5YXu/u2y0+5Rki7zDv5nm/D47lSNGQCAtsxq/AdSQ6MjzISyBft3Jd1ZBaTS4SoaINyx0oBUfp+rJf3fkj5lZr+UdJOko6c5IKU4ZgAAJkH8B9LATilmgpntLOk2FbUZ1zQeP0PSFkmvknS2u3+nmxECAIC2Ef+BNLBTilnxTknfaAak0uGSvi3pTyT9jZntuaMHBgAAVg3xH0gAO6UIzcyeLumrKu4G+yN3v7fjIQEAgFVG/AfSQlIKAAAAAOgMx3cBAAAAAJ0hKQUAAAAAdIakFAAAAADQmbmuByBJe++9t69fv77rYQAAgrjuuuvudfe1XY8jZcRmAECbFovNU5GUrl+/Xps3b+56GACAIMzsrq7HkDpiMwCgTYvFZo7vAgAAAAA6Q1IKAAAAAOgMSSkAAAAAoDNTUVMKACl69NFHtWXLFj300ENdD2VmrVmzRuvWrdNOO+3U9VAAAFOA2Ny95cRmklIAWKYtW7Zot9120/r162VmXQ9n5ri77rvvPm3ZskX7779/18MBAEwBYnO3lhubOb4LAMv00EMPaa+99iLodcTMtNdee7EaDgCoEZu7tdzYTFIKACtA0OsW//0BAKOIDd1azn9/ju8CQMJ6vZ4OP/xwPfroo5qbm9MrX/lKnXnmmcqyTFdeeaWOO+64+vjM3nvvrYMOOkjf+MY39Mgjj+gHP/iBDjroIEnSX/zFX+iEE07ocioAAIRAbN5+JKUAkLDHPOYxuuGGGyRJ99xzj/74j/9YDzzwgN71rndJkp7znOfo8ssvn/d5d955p4499tj6cwEAQDuIzduP47sAEMQ+++yj8847T+ecc47cvevhAAAw84jNk2GnFABa8K7PfU83/ccDrX7NQ37z8Xrniw/drs95ylOeon6/r3vuuUeSdNVVV2nDhg2SpJe+9KV629ve1uoYAQCYVsTmdJCUJug/fvB97fkbT9Kax+za9VBW1U+23K5dH7+nHvf4PboeCpCshY4IAZFtue1G7fvkA7XTzrt0PRQAmIfYPB9JaWLyfl+7XXCkbjjodTrij9/e9XBW1SPnH6s7fuN5eub/+bddDwVY0vaumq6WO+64Q71eT/vss49uvvnmrocD7HC/fOBnWvsPR+qGje/S7x3/uq6HA6BDxOZ0UFOamDzvazd7UP7g/V0PZdU9zn+p3sM/73oYQDK2bt2q008/Xa9//etph4+Z9fCDv9Iu9qj6v/5Z10MBAGLzhNgpTUy/v6140fJ+10NZdZn6MzFPYCUefPBBbdiwoW47/4pXvEJvetObuh4W0Bnv96s3uh0IgJlFbN5+JKWJ8Tyv3uh2IDuAyWWiSxmwmH5/4V+8jzzySB155JFjP7Z+/XrdeOONqzQqoDu5z06cBDCdiM3bj+O7ien3txVvzMAKcM/zmZgnAKA9eRkn60VcAMDUIylNTF4GWZuBe45MLmOlGwCwHfL6RBGLmgCQCpLSxOQzVCvTU05SCgDYLp6XJ4roSQAAySApTYxXQXYGkjVTPhPzBAC0Z7BTGv9EEQBEQVKamKqm1GZgBZidUgDA9pqlE0UAEAVJaWK8TtLirwD3zIvdUgAAJjVDJ4oAIAqS0sTMygrwYJ78UgGk5Pvf/742bNigjRs36vbbb9eHPvQhHXzwwXr5y1+uz3zmM7rpppvGft7DDz+sF7zgBdqwYYMuueQSXXXVVTr00EO1YcMGXX311fr85z+/g2eCVOXUlALAkBRiM0lpYvIyyEY/1jor8wSi+cxnPqMTTjhB119/vX77t39bH/7wh3XFFVfo4osvXjTwXX/99ZKkG264QS972ct08cUX6+yzz9YNN9ygW265haQUE8vz8iQRNaUAICmN2ExSmph6BzH4CnBdOxt8RxhYqTvvvFMHH3yw/vRP/1SHHnqoXvjCF+rBBx/UkUceqc2bN0uS7r33Xq1fv16SdMEFF+j444/X7//+72v9+vU655xz9IEPfEAbN27UEUccoZ/+9KeSisu93/jGN2rDhg067LDDdO211yrPcx144IHaunWrpKKhzAEHHFC///nPf14f/OAHde655+qoo47S6aefrjvuuEPHHHOM3v3ud+uyyy7Tn/3Zn2nDhg26/fbb6zncc889OuWUU/Stb31LGzZs0N/93d/pk5/8pN7+9rfr5JNP1jve8Q5dcskl9UotsJiq+y7xA0BXiM3bb27FXwE7VHUZuAWvKZ2VeSKQL5wl/fi77X7N3zhcOuY9Sz7t1ltv1cc//nH9/d//vU488UR96lOfWvT5N954o66//no99NBDOuCAA/Te975X119/vc4880xddNFFOuOMMyRJv/71r3XDDTfoa1/7ml796lfrxhtv1CmnnKKLL75YZ5xxhv7lX/5FT3va07R27VpJ0ote9CKdfvrpetzjHqc3v/nNkqQvfvGL+upXv6q9995bt956q4499lidcMIJQ+PZZ5999NGPflTvf//7dfnll0uSrr766vq5F1xwgTZv3qxzzjlne/8LYgZ53X2XkzbAzCM2JxOb2SlNTF0rE3wFmJ1SYHL777+/NmzYIEn63d/9Xd15552LPv+oo47SbrvtprVr12r33XfXi1/8YknS4YcfPvS5J598siTpuc99rh544AHdf//9evWrX62LLrpIkvSxj31Mr3rVq1qfD7ASeX824iSA6UZs3j5L7pSa2cckHSvpHnc/rHzsfZJeLOkRSbdLepW7319+7GxJr5HUl/QGd/9fqzP02eQzUmuZs1OK1Eywarpadtlll/rtXq+nBx98UHNzc/X/Rw899NCCz8+yrH4/yzJt27at/piZDX2emelJT3qS9t13X33lK1/Rtddeq4svvrj1+QArUXWpjx4nAUyA2JyMSXZKL5B09MhjV0g6zN1/R9K/STpbkszsEEknSTq0/JwPm1mvtdFC3p+RpLSeJyvdwHKsX79e1113nSTp0ksvXdbXqGpEvv71r2v33XfX7rvvLkl67Wtfq1NOOUUvfelL1etN/iN+t9120y9+8YvtHsdyPw+zybkSBsCUIjYvbMmk1N2/JumnI499yd2rlP0aSevKt4+T9Al3f9jdfyDpNknPaG20aBzfjR1sq+NXWfB5AqvlzW9+s84991xt3LhR995777K+xpo1a7Rx40adfvrpOv/88+vHX/KSl+iXv/xlfTzota99bd24YTEnnXSS3ve+99Ut6T/ykY/oIx/5yJKfd9RRR+mmm26i0REmMqgpZVETwHQhNi/MfIKW6Wa2XtLl1fHdkY99TtIl7v6PZnaOpGvc/R/Lj50v6QvuvuhSwKZNm3yS/2iQbv/uNfrtT/2Bvr3rc/X0P/tc18NZNff9ZIv2OvdQfX+nQ/TUt13d9XCAsW6++WYdfPDBXQ9jVRx55JF6//vfr02bNs372ObNm3XmmWfqqquu6mBk8417HczsOnefP3hMLNXYfOM3PqfDrjhF1+75Yj3jDf/Y9XAA7GDE5jRj84q675rZ2yRtk7TdB5fN7DRJp0nSk5/85JUMY6bUre4VewfROb4LTKX3vOc9Ovfcc5OsV8HSQsTmaqc0+NVpAFCJEJuX3X3XzP5ERQOkl/tgu/VHkp7UeNq68rF53P08d9/k7puqlsVY2qy0us+dRkdAl6688sqxK7FnnXWW7rrrLj372c/uYFRYbRFic90QkPgBIJjIsXlZSamZHS3pLZJe4u6/bnzoMkknmdkuZra/pAMlXbvyYaKS11elBE9K65pSVroBAJOjphQA0jPJlTAfl3SkpL3NbIukd6rotruLpCvKtsTXuPvp7v49M/ukpJtUHOt9nTtRoU2DVvex/7NyJQxS4e7z2rNjx5mkLwJmS9WHMfriLYCFEZu7tZzYvGRS6u4nj3n4/DGPVc9/t6R3b/dIMJFZOZZU1c7SfRfTbM2aNbrvvvu01157Efw64O667777tGbNmq6HginieREfSUqB2URs7tZyY/OKGh1hx6uOJc3OTim/VGB6rVu3Tlu2bNHWrVu7HsrMWrNmjdatW7f0EzEzuKcUmG3E5u4tJzaTlCYmr7vvxt4pzavuuySlmGI77bST9t9//66HAaDJq+7txA9gFhGb07Ts7rvoyIzslFat/Dm+CwDYHnmfkzYAkBqS0sT4jCRr+YzcxwoAaFfV6IjjuwCQDpLSxPiMdKXNy0YVWfB5AgBaVp8oIikFgFSQlCZm0Oo+9vFduu8CAJaj7lIfPE4CQCQkpYmpW90H30F0uu8CAJbDZ+NEEQBEQlKamEFNaewV4Lxf7pSSlAIAtgc7pQCQHJLS1FSt7oOvADsr3QCAZSB+AEB6SEoTU7W6j15rWe0I99gpBQBsB2pKASA9JKWpqXdKYydr1JQCAJal2il1dkoBIBUkpYmpV4CDJ2vVPaXUlAIAtkvVe4H4AQDJIClNjc/G8d3qnrmMlW4AwHaoa0qjx0kACISkNDGzslPqrHQDAJZjRuIkAERCUpqYagU4C95VsKopJSkFAGwXdkoBIDkkpamp7ymNHWzdqSkFACwDJ20AIDkkpYmZla60nhc7wdF3hAEA7aKmFADSQ1KaGp+NFeCqpnTOYs8TANCyGYmTABAJSWlqqhXg6DuIjUvPq91hAACWVHZtDx8nASAQktLEVDuIveArwHl/ML9+f1uHIwEAJKXuvdBf4okAgGlBUpoan42a0uZOac5OKQBgQjYrJ4oAIBCS0tTMSFfBakdYknJ2SgEAE3JqSgEgOSSlqanuKfXgK8CNrol5zhEsAMCE6L4LAMkhKU2Mz+BOab9PUgoAmIzNSJwEgEhIShNjXt3fGTvY+tBOaey5AgBaRE0pACSHpDQxM1Mr09gpdWpKAQCTqspcosdJAAiEpDQ1dbCNvQLcvJuUmlIAwMRISgEgOSSlialqZeYseLAduhKGpBQAMBmblRNFABAISWlimrWWHrnWclbmCQBo14ycKAKASEhKU+PNrrRxay2Hu+/GnScAoGX11WksaAJAKkhKUzMrXWmHdko5vgsAmIxRUwoAySEpTYw1ErQ88g7i0Dz5xQIAMCGSUgBIDklpcgY1MqEbALFTCgBYBqOmFACSQ1KamqFay7jJWjMRzfPAO8IAgFbRfRcA0kNSmhibkZpS88EKN913AQCTK+IHO6UAkA6S0tQ0uu964JpS555SAMAy1Dul5ixqAkAiSEqTM4s1pfxSAQCYzKycKAKASEhKEzPUfTdwUjozXYYBAO1qJKXccw0AaVgyKTWzj5nZPWZ2Y+OxPc3sCjO7tfx7j/JxM7MPmdltZvYdM3v6ag5+Js3IDqI358kF6ACACQ3vlMZdvAWASCbZKb1A0tEjj50l6cvufqCkL5fvS9Ixkg4s/5wm6dx2homKzcoKcLN2ll8qAAATskbXXU7aAEAalkxK3f1rkn468vBxki4s375Q0vGNxy/ywjWSnmBmT2xprJAkzcj9nax0AwCWgZpSAEjPcmtK93X3u8u3fyxp3/Lt/ST9sPG8LeVj85jZaWa22cw2b926dZnDmD1DwbYfN9g2a0oV+D5WAJgmEWLz8Iki4gcApGDFjY7c3aXtvwzM3c9z903uvmnt2rUrHcbMMJ+RndJml2FqSgFgh4gRmxsxI3ScBIA4lpuU/qQ6llv+fU/5+I8kPanxvHXlY2iJDd3fGbhWJm/ex8ovFQCAyWSNOBm69wIABLLcpPQySaeWb58q6bONx19ZduE9QtLPG8d80QYf7CBG7r47tCPsJKUAgMmYN+/zjhsnASCSuaWeYGYfl3SkpL3NbIukd0p6j6RPmtlrJN0l6cTy6Z+X9CJJt0n6taRXrcKYZ9rwTmngZI3uuwCAZbBZaQgIAIEsmZS6+8kLfOj5Y57rkl630kFhMbOxUzo78wQAtIl7SgEgPStudIQdq1krE/n+tWb3XVa6AQCTGu5SHzdOAkAkJKWpGaq1DLyDOCvzBAC0avj47nZfDgAA6ABJaWJmpVZmqNERK90AgAllYqcUAFJDUpqYbEYaHTUbOonuuwCACQ3XlJKUAkAKSEoT02x1r9D3d9LoCACw/YZPFBE/ACAFJKWJaQbbPHCtpXElDABgGbKhe0qJHwCQApLSxJjnyt0kSR54p3RonoGTbwBAu0z9QfxgpxQAkkBSmhjzXNvUkyR55FrLxjzFSjcAYEImr+MHjY4AIA0kpYkxNZLSwMnaUPLNSjcAYELZ0OIt8QMAUkBSmphMubbZnKTYyZo15xl5RxgA0Kqh+BF48RYAIiEpTczwDmLcYGtDx3fjJt8AgHZljRNFNDoCgDSQlCamOL5b7SDGTdaKpJSVbgDA9snc6/gR++o0AIiDpDQxmbv61U5p6AYOeT1PBU6+AQDtytSv40fkq9MAIBKS0sSY+urXtZa+xLPTlXlf22yn4h1qSgEAEzL5IE6yUwoASSApTUwz2CqPu1Nq7uob3XcBANsnU95YvCUpBYAUkJQmJvPBsdbQNaWNXyrYKQUATCpTo8yFngQAkASS0sQ0k7XIwdY8V14npXGTbwBAu7JmmQsnbQAgCSSliSmOJZW1loGDrXmuftU9MXDyDQBoV1HmUsTJyIu3ABAJSWliMnfl2QzslCofzDNwQycAQLt63owfcRdvASASktLEZOort/g1pZnyep4WOPkGALTL5IM4GfrqNACIg6Q0MUWwrY7vxk3WiprS8vgVjY4AABPqqRk/OGkDACkgKU1M1jjWGrkBkM3IPAEA7RqKH4GvTgOASEhKE5PJ6660kWtKM3d5efyKpBQAMKlip5SaUgBICUlpYjL15TOwg2jqy5Wp7xb6mDIAoF098zpORl68BYBISEoTY3L5DOyUFvPM1FcWOvkGALQn7xdxsYqTka9OA4BISEoTU7S6LxsdBU7WMs8ly+TKJBodAQAmkJeLtVWcjLx4CwCRkJQmxuSN47txg60pr3dKje6JAIAJ9MsrYJx7SgEgKSSliekpl2agpjQrk1KXhU6+AQDt8eq4bt19l/gBACkgKU2MKZfXx5ICJ6XukvXUN2pKAQCTGeyUxi9zAYBISEoTM7xTGncFOFNfLpPLZKx0AwAmkI/slFJTCgBpIClNTM9c3itXgAMHW5NLWU+5epKoKQUALK3uvttjpxQAUkJSmpAq2Ko+lhQ3WatqSnNqSgEAE6p3Rum+CwBJISlNSF4H217xd+BkLVNRU5ork7HSDQCYQFVTql78hoAAEAlJaUKqYGvWU98tdLDN1JfKmtLIx5QBAO2proCx+kQR8QMAUkBSmpCq265nPfWVhU7WivtYi5pSo6YUADCBusyFmlIASApJaUIGO6UmVxa6prTnuWSZcqOmFAAwmarMxXrxr04DgEhIShMyaHVf7JRa4GTNqCkFAGyn+TulceMkAESyoqTUzM40s++Z2Y1m9nEzW2Nm+5vZN83sNjO7xMx2bmuws64OtpYVtZaBk7Wecnm5I0xSCgCYRLUzms3A1WkAEMmyk1Iz20/SGyRtcvfDJPUknSTpvZL+2t0PkPQzSa9pY6BotLa3nvqWhV4BNuXFTqlloZNvAEB78rzqvhv/6jQAiGSlx3fnJD3GzOYkPVbS3ZKeJ+nS8uMXSjp+hd8DpbqmNCt2SiPvIPaU1zvCkecJAGiP1zWl1ZUwcRdvASCSZSel7v4jSe+X9O8qktGfS7pO0v3uXi5Vaouk/VY6SBSqVvdFrWUv9A5iz1zKenL1il1TAACW4GWZi2Xxr04DgEhWcnx3D0nHSdpf0m9K2lXS0dvx+aeZ2WYz27x169blDmOmNGtK88A1pUPztLjzBIBpk3pszqt7Si1THvzqNACIZCXHd18g6QfuvtXdH5X0aUnPkvSE8jivJK2T9KNxn+zu57n7JnfftHbt2hUMY3bUre6zqittzGCb582GTjQ6AoAdJfXYXC9qZnNFUkpNKQAkYSVJ6b9LOsLMHmtmJun5km6S9FVJJ5TPOVXSZ1c2RFQGwTZ2993Bfaw9ucVNvgEA7aprSs2Uy4gfAJCIldSUflNFQ6NvS/pu+bXOk/RWSW8ys9sk7SXp/BbGCQ1a3VtZUxp1B7Gap1c7wmKlGwCwtGajo2KnNGacBIBo5pZ+ysLc/Z2S3jny8B2SnrGSr4vxqlb3lsWutRzslFb3lLLSDQBYmg/VlBrddwEgESu9EgY7UH1Pafia0jLZzsrju+yUAgAmkNdXpxX3XEc9UQQA0ZCUJqRudV82AIrawKHZfZeaUgDAxKoylyzj+C4AJISkNCF1q/t6BThmslbvCFuvPL4bM/kGALRrtEs9SSkApIGkNCF5fSl40X036rGkuqY0y+RmyhRzngCAdtW9F6xXxsmYi7cAEA1JaUKGdxB7sqDJWtWoQtZTbnG7DAMAWlYd3+2xUwoAKSEpTYg3jyUF7r7brCmVsrDJNwCgXVWczCwrGwISPwAgBSSlCfFGTakHDrbNmiCneyIAYELVPdfqFb0Xoi7eAkA0JKUJqVvdB+9KW++UZsU8qSkFAEzCvYiTWUZNKQCkhKQ0JSO1MlHv76xWus1i7wgDANrleREXi+67vbBXpwFANCSlCamPtZb3lEZdAa67J7JTCgDYDt6Mk8ZOKQCkgqQ0IfUOYjZXHN8Nu1NaHd/tSYHnCQBoWZmEZr258kQRi5oAkAKS0oT4yA5i1BVg7zdXujNlQecJAGhX3q8WbzPKPwAgISSlCam672Z1992YO4h5s8uw9dgpBQBMptopzei+CwApISlNiA91pbWwtZZV913LMklx5wkAaNfgPm92SgEgJSSlCalb3ffmlFsvbLCtjinLip3SLOg8AQAtq08UzYUucwGAaEhKE1K3urdMCtzAYdDQqVc2dIo5TwBAuwY7peU9pZR/AEASSEoTMhRsLe6xpLxfNXQquu9yfBcAMIm690IvU64eO6UAkAiS0pT4cFIaNVmrfqmouu+y0g0AmMjQ4i07pQCQCpLShFSt7rNeL3QDh3pHuFfslPaCJt8AgHZV5R+93lzRvZ2dUgBIAklpShqt7kPvlDZqSmU9akoBAJPx0e677JQCQApIShPSbHWvwMda87L7rlns5BsA0LLmfd6Br04DgGhISlMy0uo+i3osqdop7fWkrKeMlW4AwASqxdss64W+Og0AoiEpTchw991e2J3S+pcKyySx0g0AmFDVKK83p8hXpwFANCSlCWm2uo+crFU1per15Fkv7DwBAC2rd0qz0FenAUA0JKUpGd0pDRps3Yua0ox7SgEA28Or7rv0JACAlJCUJmS41X3cYOt5cSy56r6bBT2mDABoV7OmNPLVaQAQDUlpSnw2uu/WtbOWySzTnPFLBQBgadUVMFnwxVsAiIakNCVDre4DB9vqPtbenDzrFQ/lQecKAGiN+6CmNPLiLQBEQ1KakOaxpMi1lnm/7J5Y/lIhSf3+ti6HBABIgY+UuUS9Og0AgiEpTUmz1X3gpLTeKc16Mit2SnN2SgEAS7B8UOYS+eo0AIiGpDQlQ63u4zYA8qFfKkySlLNTCgBYgnuubV79ahP36jQAiIakNCWNVveyTFnUroJ17exc0YFXUp5zBAsAsATvK1exmBn56jQAiIakNCGzUlPqjftYBzWlJKUAgCV4rrz81SZ0Q0AACIakNCHNVvehk9Jqp7SXSRk1pQCAyVjer5NSuu8CQDpIShPSbHXvWdyaUo3ZKXVqSgEAS3J2SgEgQSSlKWm0ug+9U5oP5klNKQBgYnlfedkgL3KcBIBoSEoT0mx1b4G771ZXwjTvKSUpBQAsxRo1pSSlAJCOFSWlZvYEM7vUzL5vZjeb2TPNbE8zu8LMbi3/3qOtwc66Zqt7t0yZeb2rGErdfbcnlfeUhpwnAKBd3pfXx3cDL94CQDAr3Sn9G0lfdPenSnqapJslnSXpy+5+oKQvl++jDY1W95EbADW7DFtWdd+lphQAsBRXv7lTypUwAJCEZSelZra7pOdKOl+S3P0Rd79f0nGSLiyfdqGk41c2RNSGjiUVyWnIZK38JcJ6c3Xy7RzfBQAswfK+XNSUAkBqVrJTur+krZL+h5ldb2YfNbNdJe3r7neXz/mxpH1XOkgUhlvdB24AlA+6DFtVU9rnFwsAwBKoKQWAJK0kKZ2T9HRJ57r7Rkm/0shRXXd3aXxBh5mdZmabzWzz1q1bVzCMWTJodT9I1uLulPZ6vbr7LjulALD6Uo/N5vmgpjTy1WkAEMxKktItkra4+zfL9y9VkaT+xMyeKEnl3/eM+2R3P8/dN7n7prVr165gGDOk2ep+RmpKB913AybfADBl0o/NuXJjpxQAUrPspNTdfyzph2Z2UPnQ8yXdJOkySaeWj50q6bMrGiFqo63uJanfj7eDaF6sbGeNe0rpvgsAWEoRJ4vF29BXpwFAMHMr/Pz/KuliM9tZ0h2SXqUi0f2kmb1G0l2STlzh90Cl0eq+2ilVwGOt7oOa0tC1swCAVg0d321cnVZ1cgcATKcVJaXufoOkTWM+9PyVfF0sZNDq3izwVSl1Temcsl75ywU7pQCAJZj3B8d3G2UuPZJSAJhq/JROyFCr+8A1pVbuilpzpzRi8g0AaJf74ERR5KvTACAYktKUjKkpjdiV1j3XNi93hKuaUi5ABwAswZplLpR/AEAySEoT0qyVsSxwsPX+oFFFVvwdMfkGALTN6+O7oa9OA4BgSEqTks9GsG3sCJsVZc8hk28AQKsyn40yFwCIhqQ0Ic1W96qvSonX7t7y/iAprZpTBLz6BgDQMs+Vl8d2I1+dBgDRkJQmpHl8V5F3SuWNpLRc6aamFACwBFM+b6c04tVpABANSWlCmq3uBw2AAgbbvK/cql8qyoZOrHQDAJZgnsvLndLQV6cBQDAkpSlptrrPqp3SeMmaNWpKs6yoKQ2ZfAMAWpX5/J1SakoBYPqRlCak2eq+agAU8qqU5jyzuFffAADaNmgIGPnqNACIhqQ0KY1W9+VVKRF3SiVXf6Sm1FnpBgAsIfNcmoWr0wAgGJLShDRb3Vt1rDVgsLW8MU9WugEAEzLuKQWAJJGUpqTR6r4OthGTtWZNaS/wMWUAQKvM+/Lq+G7gq9MAIBqS0oQ0W91br9pBjJesNa++qWtKWekGACzB5DNydRoAxEJSmpBmq3tZVSsTMdjmjdrZaqeUlW4AwOIy79fJaOir0wAgGJLShDRb3Wf1peAxd0rzep5la/+QyTcAoE1FTWkZHwNfnQYA0ZCUJqXZ6j5uV8Hm8d3BPaXxkm8AQLuK7rtVozziBwCkgqQ0IUOt7kPXlPYHyXfdqCJe8g0AaJdpUOYS++o0AIiFpDQhw63uq2Qt4LFWHzSqyHpxjykDANplyiWLf3UaAERDUpqQZqv7KlmL2OrevN84vss9pQCAyWSNhoChr04DgGBIShPSbHVfBduQO6XNHWFqSgEAEyqO78YvcwGAaEhKEzK+1X28YJt5v9FluPwnyko3AGAJmbxOSmNfnQYAsZCUJqTZ6t4iNwDyvJ5n1purHwMAYDHNhoCRr04DgGhIShPSbHUfuQGQqXEfq1FTCgCYTPP4buSr0wAgGpLShAy1ug/cwMGajSqq5JudUgDAEjLlgzIXakoBIBkkpQlptrqvd0o9XlKaeWOnNPIxZQBAq0wuz6rF28BXpwFAMCSlCWm2us8C379myuvuu4Oa0njzBAC0q6e8PrY7uDqN+AEA046kNCFDtTKBu+82j+9mHN8FAEzINOi9MGgISPwAgGlHUpqQrJGU1slaP94KsA0d3+VKGADAZDLljeO7NMoDgFSQlCYkcx8cSyqTtZA7pY2GTr3y+K67dzkkAEACijg5XP4RMU4CQDQkpQnJ1K93EC1wTWnWuBLGyuTbAs4TANCurFFTWsUP79PoCACmHUlpQkxe15JWO6UWcAW4WVM62CklKQUALK4ocxlZvOWkDQBMPZLShAzXlMbdKW02dKquhKHREQBgKc2d0iwrklNxJQwATD2S0oRkatSUWtya0uGaIJJSAMBkijhZLWpSUwoAqSApTUimvupW91WyFnKntC9v/NPsu4WcJwCgXXOW12Uu9dVpxA8AmHokpQkxed3qPvKxVpMP7mOV1FcWcp4AgPbU95GOnrThnlIAmHokpQnpeT6/1X3AFeCsMU9Jxa4pjY4AAIvoV1126+O73FMKAKkgKU2INWpKq1b3EXcQm42OpGKn1OieCABYRF7uiFodJ6kpBYBUkJQmpNdodd+rGwDFWwHORpJSl4WcJwCgPXm5U+ojV6fRkwAApt+Kk1Iz65nZ9WZ2efn+/mb2TTO7zcwuMbOdVz5MSMUO4qDVfdxamaL7bq9+v2/UlAIAFpeXyaeVi7dVmQvxAwCmXxs7pW+UdHPj/fdK+mt3P0DSzyS9poXvARU7paMNHDzgDmKmfrE7WnKZjJVuAMAi+v0yTmQjV6cRPwBg6q0oKTWzdZL+UNJHy/dN0vMkXVo+5UJJx6/ke2CgZz4ItsG779Yt/SXl6kmiphQAsLB8pPuucc81ACRjpTulH5T0FknVT/y9JN3v7mULPG2RtN+4TzSz08xss5lt3rp16wqHEV9erQCXwbZXHUuKeHx3pKY0p6YUAHaIlGOz1913hxdv2SkFgOm37KTUzI6VdI+7X7ecz3f389x9k7tvWrt27XKHMTOqWpl6BbismYm4ApxpuKY0VyYLOE8AmDYpx+a6pjQbvjotYpwEgGjmVvC5z5L0EjN7kaQ1kh4v6W8kPcHM5srd0nWSfrTyYaLf36Y5NVvdZ+p7zB3ETH1ppKaU7okAgMUMFm9Hr04jfgDAtFv2Tqm7n+3u69x9vaSTJH3F3V8u6auSTiifdqqkz654lJCXx3Q9G95BjLgCbPKRefaKOlMAABbgIzWlPWpKASAZq3FP6VslvcnMblNRY3r+KnyPmdMva2XqY7sqk9KANaU9H3QZlqTcYu4IAwDaU8fJbLSmNF6cBIBoVnJ8t+buV0q6snz7DknPaOPrYqDuKji0g2iygMmaUVMKANhOdUOjbPjqNBY1AWD6rcZOKVbBaPddKe7x3Z5yuTVrSklKAQCLy/tFnLCR7rv0JACA6UdSmggfaeAgVVelxEvWTPnwPC1m8g0AaI+PdN+tr05zehIAwLQjKU3EoFamWWuZhTy+29NwTanL2CkFACwqz8t7Sssd0sHVafHiJABEQ1KaCK+SMovffbdnPlQ76+oVu6cAACzAyzIXq+7zrq9OI34AwLQjKU3ErNSUjp2n8UsFAGBxeRknbPTqNGpKAWDqkZQmIq9rZZo7iPG67w4uP28e36XREQBgcdWi5lCZizJqSgEgASSliah3EOcF21jJ2uA+1kbyHbR2FgDQnoUaAhI/AGD6kZQmorr824Lf31nN00eOX5lY6QYALMzHnCiKuHgLABGRlCai6io42n03WrAd7JSO3lPKSjcAYGE+tqbU6L4LAAkgKU1EfSwpfE1pmWRnI8d32SkFACwirxc1R69Oi7V4CwARkZQmYrTVvSTl6oVr4DCu+y41pQCAJVVlLj2O7wJAakhKEzGu1b1bvJ3ScY0qiuO7sZJvAEC76i71wa9OA4CISEoTsVCre1OsYFvXlGbNnVJTFmyeAIB21Q0Bs7nBYwHLXAAgIpLSRCy8gxgrWasaVQy19LdeuHkCANrl4xoCslMKAEkgKU3E2Fb3AbvvjqspVcAdYQBAu6pFzWz0SrFgcRIAIiIpTcS4VvcRd0rzMcm30z0RALCEqiGggl+dBgARkZQmYlyr++KqlFjBNh/zS4VbRk0pAGBR9U5pj5pSAEgNSWkqxrS6L4JtrGStblQRvHYWANAuH9t9N97VaQAQEUlpIsa3uo/XACgf06iCnVIAwFLG9V6IeHUaAEREUpqIsa3uzcId3627DDd+qZBlMrHSDQBYhI9pCEijPABIAklpIsa1uveAV6VUjSpGa2czVroBAIuoFm+zHuUfAJAaktJEjGt1H7GmNB/XZdh67JQCABZVnbTJjO67AJAaktJEjGt1H7n7bnNHWDJqSgEAi6qvTmOnFACSQ1KaiPGt7jNlwYJtXVNqIzulweYJAGjZQvdcs6gJAFOPpDQR41rdF8daYwXb8d0T6b4LAFhcXVOajd5TSvwAgGlHUpqIhVvdxwq2PqamlO67AICluM8v/4h4dRoARERSmooxre4jNgDK+2WX4dHuu+yUAgAWM677bsCr0wAgIpLSRIxvdW/xrkrJ5zeqEEkpAGAJdffd0cVbdkoBYOqRlCZiXKv7kDulY2pnSUoBAEvy8Yu3JKUAMP1IShMxrtW9LN5VKdWOsDUbVVhPWbDkGwDQMh/XEJDuuwCQApLSVNTHkoaTtWjB1vOyprSZfGe9cFffAADaVS1q9hpXp+VG/ACAFJCUJmKwgzh8KXi0YFvfx9o8vqt4O8IAgJaNaQgosVMKACkgKU1E1eo+6w3XWkYLtt4vGzc1a4KyHkkpAGBx42pKLaOmFAASQFKainE7pRZ4pzQb7b5LTSkAYBF1mctoTSnxAwCmHUlpIqruu81amYjdd6t5Dh2/MnZKAQCLcy/i4VCcVBbv6jQACIikNBV1rUzz+G68WksfcyWMWaY5izVPAEC7LJ8fJ50rxQAgCSSlqRhzrNUj7iDWtbONle5yzlWzJwAARlW9F5o7peL4LgAkYdlJqZk9ycy+amY3mdn3zOyN5eN7mtkVZnZr+fce7Q13dnldKzNaaxkrUcv7Ve3scEMnSer3t3UxJABACsYu3nJ8FwBSsJKd0m2S/pu7HyLpCEmvM7NDJJ0l6cvufqCkL5fvY6XKYGvRV4B9fvJtVryds1MKAFjI2O678XovAEBEy05K3f1ud/92+fYvJN0saT9Jx0m6sHzahZKOX+EYIY3tKijL1Au2U+pja4JMkpSzUwoAWEjeV99t5MF4vRcAIKJWakrNbL2kjZK+KWlfd7+7/NCPJe3bxveYeeUKcG90BTjYlTCD41eDHeGqE2+ecwQLALAAz9Uf+bUmZJwEgIBWnJSa2eMkfUrSGe7+QPNjXvRnH3tuxsxOM7PNZrZ569atKx1GeAvVlMbdKR2epyT1+ySlALCako7N3pfPS0rj9V4AgIhWlJSa2U4qEtKL3f3T5cM/MbMnlh9/oqR7xn2uu5/n7pvcfdPatWtXMoyZYOX9a1nwmlKva4Ia/zQzakoBYEdIOTab+7yd0ohxEgAiWkn3XZN0vqSb3f0DjQ9dJunU8u1TJX12+cNDxX1+TalbL9xOqRbZKXVqSgEAC/G+XMM1peyUAkAa5pZ+yoKeJekVkr5rZjeUj/25pPdI+qSZvUbSXZJOXNEIUahrShsvWRZvBbi6i7Q5T2pKAQBL8lx9m79TSlIKANNv2Umpu39d0mibu8rzl/t1MZ6N6Uobsaa0uhJm3D2lJKUAgIVYPn+nlKQUANLQSvddrD73XNt8dAW4p8y83l0MYczl5yrvKQ01TwBAy1y5esOPWE9ZsBNFABARSWkqvK98zAqwFKsB0Lguw9WuaZ+aUgDAQhaIkxlXwgDA1CMpTYXnykdfrizgsdbylwcbqp2tdkoDzRMA0CobFyc5vgsASSApTYTl/bHBVgq2g5jP7zJs1Y5wn18sAAALoKYUAJJFUpoMn5eUWsRay7rLcPP4LjulAIDF2bia0oyaUgBIAUlpKvK+chtZAQ5YazmupnRQOxtnngCAlvmYOMlOKQAkgaQ0EeNrZar7O+MEXPNiRTsbc09pqB1hAECrxsVJo/suACSBpDQV3pcvUFPqkXZKfX5N6SD55vguAGA883xenHTL4l2dBgABkZQmw9UfXQHO4u2UDmpKBzulWa9MviPNEwDQLs+V22iX+oBxEgACIilNhC3QVVCKVVNq5W6ojdspDTRPAEC7ip3S0ThZvB8pTgJARCSlqRh7T2nZDCjQCrB7rm0+fkfYuQAdALAAUy4f6b5L+QcApIGkNBHjamWq+zv7kbrSel/5yEq3ZcX7XAkDAFiQ5/O67w7uuQ4UJwEgIJLSZMyvlQnZlXZs98SivpSVbgDAQsYt3lJTCgBpIClNRNHqfngF2OsV4DjJmuX9+UlpVV8aaJ4AgHaZ9+u4OHiw6r1A/ACAaUZSmoixx3frWstIwdbHJKXlSjc1pQCABdiY+DHovRApTgJAPCSliTDvzz++G7FWJu/PqwlSVt3Hyi8VAIDxbMx93hawSz0ARERSmgr3+bUyvapWxjsY0OqwMTWlWVbUlMbaEQYAtMnk84/vUlMKAEkgKU3EYivAHqz77vxjytU8SUoBAOMtVlNK/ACA6UZSmgxfpPtupGDr6i9UO8tKNwBgATbmRFHdkyBUnASAeEhKE5F5Xz7Sfbe+FDxQraXl8+dprHQDAJaQKZdzTykAJImkNBWeKy+T0EpWd98NtIM4rqa0V9WUBponAKBVNiZOqj5pE6f3AgBERFKaCFM+f6c0C3hP6dirb6ruu6x0AwDGM+Wa92sNO6UAkASS0kSY5/KRFeCYNaX5mNrZaqeUlW4AwHhFnJyF+7wBIB6S0kRkPn+nNGKtZXElzPA8s6x8P1KXYQBAq4qa0tErYeKdKAKAiEhKk7HYDmKcWstxx3ezgPMEALRrbPmHET8AIAUkpYnIfH6tjJU7iJFWgM3785JvhTymDABo07id0ohxEgAiIilNhI29p7RaAQ4UbMfcM5f1ylpa7ikFACzA5HVjo/qxKk6yqAkAU42kNBHm/TErwFVX2jjJmnl/zPHdeLWzAIB2ZePiZNV9l/gBAFONpDQRpjE7iOWxVkXaKV10RzhO8g0AaJfJ53ep71WLmsQPAJhmJKWJyLw/71hS3VUw0Apw5v15XYarnVIFmicAoF2Zcmn0Pu8ySc3p3g4AU42kNBFFTenwCnDVlTbUTqnn8+fZm6s/BgDAOOPu865PFLFTCgBTjaQ0EUX33ZF7SiPWlGr+faxZwPtYAQDtGntPab1TSvwAgGlGUpoI08IrwJFqLcetdFvVfTfQPAEA7RrbfZeaUgBIAklpIky5ZCM7iL3q/s44tTKZj9kp5Z5SAMASMuXybGRR0+LFSQCIiKQ0Edm4HcQ6WfMuhrRK8nndd6kpBQAsZVyjo8HibaQ4CQDxkJQmwsbUygyuSomzAlzUzo6/jzVUQycAQKsy5XUNacXqngRx4iQARERSmohMPi8pzbJyRThQrYyNuae0R00pAGAJ4+KkBey9AAARrVpSamZHm9ktZnabmZ21Wt9nVmQ+fwW4OtYaqdbSvD/vmHKvnie/VAAAxss8l7L4cRIAIlqVpNSKzgJ/K+kYSYdIOtnMDlmN7zUrxl2VUjVwiLSDaPL58+T4LgBgCeNqSuv4waImAEy11dopfYak29z9Dnd/RNInJB23St9rJmQatwJcXQoeJ1nLFtkpjTRPAEC7xnbfLXsvcE8pAEy3uVX6uvtJ+mHj/S2S/tMqfS9J0k/v+ZG2nvefV/NbdOop/sCYWpni/f2++7e65aZ/6mJYrXvSti26ZZd9hx6rroR5yl3/rFv++5UdjArAjrLtqHfo0Gf9YdfDQEvu/Y+7dN/HTtwh3+sge3TePaVV74W11/yVbvnWOTtkHAAQTf8F79IhRxy9qt9jtZLSJZnZaZJOk6QnP/nJLXy9TA/P7brirzOtbp7bqMc+7Y+GHtv7N35L39r9aO3y8NaORtW+2+YOU37YCUOPZb2ertnnRD32gds7GhWAHWWnuc7CEtR+bFa242Lzd+Z+T3tseMnQY098ymG6brfnaadHf75DxgAAEe3cW/3YbO7t391lZs+U9Jfu/gfl+2dLkrv/v+Oev2nTJt+8eXPr4wAAzCYzu87dN3U9jpQRmwEAbVosNq9WTem3JB1oZvub2c6STpJ02Sp9LwAAAABAolZlL9bdt5nZ6yX9L0k9SR9z9++txvcCAAAAAKRr1Q4Iu/vnJX1+tb4+AAAAACB9q3V8FwAAAACAJZGUAgAAAAA6Q1IKAAAAAOgMSSkAAAAAoDMkpQAAAACAzpCUAgAAAAA6Q1IKAAAAAOgMSSkAAAAAoDMkpQAAAACAzpCUAgAAAAA6Y+7e9RhkZlsl3dXSl9tb0r0tfa1pxjxjYZ6xMM/u/Za7r+16ECkjNi8L84yFecbCPLu3YGyeiqS0TWa22d03dT2O1cY8Y2GesTBPYNis/FthnrEwz1iY53Tj+C4AAAAAoDMkpQAAAACAzkRMSs/regA7CPOMhXnGwjyBYbPyb4V5xsI8Y2GeUyxcTSkAAAAAIB0Rd0oBAAAAAIkIk5Sa2dFmdouZ3WZmZ3U9nraY2ZPM7KtmdpOZfc/M3lg+vqeZXWFmt5Z/79H1WNtgZj0zu97MLi/f39/Mvlm+rpeY2c5dj3GlzOwJZnapmX3fzG42s2dGfD3N7Mzy3+yNZvZxM1sT5fU0s4+Z2T1mdmPjsbGvoRU+VM75O2b29O5Gvn0WmOf7yn+73zGz/2lmT2h87OxynreY2R90MmhMFWJz+j/LJWJzpNeT2Exs7mTQEwiRlJpZT9LfSjpG0iGSTjazQ7odVWu2Sfpv7n6IpCMkva6c21mSvuzuB0r6cvl+BG+UdHPj/fdK+mt3P0DSzyS9ppNRtetvJH3R3Z8q6Wkq5hvq9TSz/SS9QdImdz9MUk/SSYrzel4g6eiRxxZ6DY+RdGD55zRJ5+6gMbbhAs2f5xWSDnP335H0b5LOlqTy59JJkg4tP+fD5c9mzChic/o/yxuIzQFeT2IzsVlTHJtDJKWSniHpNne/w90fkfQJScd1PKZWuPvd7v7t8u1fqPghuZ+K+V1YPu1CScd3MsAWmdk6SX8o6aPl+ybpeZIuLZ+S/DzNbHdJz5V0viS5+yPufr8Cvp6S5iQ9xszmJD1W0t0K8nq6+9ck/XTk4YVew+MkXeSFayQ9wcyeuEMGukLj5unuX3L3beW710haV759nKRPuPvD7v4DSbep+NmM2UVsTvRnXBOxOdbrKWLz8Y3Hic1TJEpSup+kHzbe31I+FoqZrZe0UdI3Je3r7neXH/qxpH27GleLPijpLZLy8v29JN3f+J8swuu6v6Stkv5HeRTqo2a2q4K9nu7+I0nvl/TvKgLezyVdp3ivZ9NCr2Hkn0+vlvSF8u3I88TyzMS/CWJziNeV2Bzr9WwiNicyzyhJaXhm9jhJn5J0hrs/0PyYFy2Uk26jbGbHSrrH3a/reiyrbE7S0yWd6+4bJf1KI8eBgryee6hYndtf0m9K2lXzj5qEFeE1XIqZvU3FEcaLux4L0BVicxjE5hkQ4TVcSsqxOUpS+iNJT2q8v658LAQz20lF0LvY3T9dPvyT6phB+fc9XY2vJc+S9BIzu1PFEa/nqajveEJ5xESK8bpukbTF3b9Zvn+pikAY7fV8gaQfuPtWd39U0qdVvMbRXs+mhV7DcD+fzOxPJB0r6eU+uFcs3DyxYqH/TRCbQ/0sJzbHej2biM2JzDNKUvotSQeW3cN2VlHQe1nHY2pFWbtxvqSb3f0DjQ9dJunU8u1TJX12R4+tTe5+truvc/f1Kl6/r7j7yyV9VdIJ5dMizPPHkn5oZgeVDz1f0k0K9nqqOBp0hJk9tvw3XM0z1Os5YqHX8DJJryw7/R0h6eeNo0TJMbOjVRzle4m7/7rxocsknWRmu5jZ/iqaR1zbxRgxNYjNif+MIzbHej1FbCY2T3NsdvcQfyS9SEW3qdslva3r8bQ4r2erOGrwHUk3lH9epKKm48uSbpX0L5L27HqsLc75SEmXl28/RcX/PLdJ+mdJu3Q9vhbmt0HS5vI1/YykPSK+npLeJen7km6U9A+Sdonyekr6uIp6nEdVrLC/ZqHXUJKp6EB6u6Tvquh62PkcVjDP21TUp1Q/jz7SeP7bynneIumYrsfPn+7/EJvT/1nemDOxOcDrSWwmNnc9/oX+WDlYAAAAAAB2uCjHdwEAAAAACSIpBQAAAAB0hqQUAAAAANAZklIAAAAAQGdISgEAAAAAnSEpBRJjZv/a9RgAAMAAsRlYGa6EAQAAAAB0hp1SIDFm9suuxwAAAAaIzcDKkJQCAAAAADpDUgoAAAAA6AxJKQAAAACgMySlAAAAAIDOkJQC6aFlNgAA04XYDKwASSmQEDPbS9JPux4HAAAoEJuBlSMpBRJhZr8p6WpJ7+96LAAAgNgMtMXcOW0AAAAAAOgGO6UAAAAAgM6QlAIAAAAAOkNSCgAAAADoDEkpAAAAAKAzJKUAAAAAgM6QlAIAAAAAOvP/AxVAbdTJLxFzAAAAAElFTkSuQmCC\n",
|
||
"text/plain": [
|
||
"<Figure size 1152x432 with 2 Axes>"
|
||
]
|
||
},
|
||
"metadata": {
|
||
"needs_background": "light"
|
||
},
|
||
"output_type": "display_data"
|
||
}
|
||
],
|
||
"source": [
|
||
"N = 128\n",
|
||
"\n",
|
||
"xk = 2*np.pi*np.arange(N)/N\n",
|
||
"yk0 = np.exp(20j*xk) + np.exp(40j*xk)\n",
|
||
"yk1 = np.exp(5j*xk*2)\n",
|
||
"\n",
|
||
"fig, ax = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(16,6))\n",
|
||
"\n",
|
||
"ax[0].set_xlabel(\"j\")\n",
|
||
"ax[1].set_xlabel(\"j\")\n",
|
||
"\n",
|
||
"ax[0].set_title(\"$y_k = e^{20ix_k} + e^{40ix_k}$\")\n",
|
||
"ax[0].plot(np.abs(DFT(yk0)), label=\"DFT\")\n",
|
||
"ax[0].plot(np.abs(np.fft.fft(yk0)), label=\"numpy.fft.fft\")\n",
|
||
"ax[0].legend(loc=\"upper right\")\n",
|
||
"\n",
|
||
"ax[1].set_title(\"$y_k = e^{i5x^2_k}$\")\n",
|
||
"ax[1].plot(np.abs(DFT(yk1)), label=\"DFT\")\n",
|
||
"ax[1].plot(np.abs(np.fft.fft(yk1)), label=\"numpy.fft.fft\")\n",
|
||
"ax[1].legend(loc=\"upper right\")\n",
|
||
"\n",
|
||
"# TODO: So the graphs overlap completely. Is this good enough?\n",
|
||
"# To make it more clear, we could mirror one of the graphs (multiply by -1),\n",
|
||
"# like what is often done in spectroscopy, or we could add the difference.\n",
|
||
"\n",
|
||
"fig.show()"
|
||
]
|
||
},
|
||
{
|
||
"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": 4,
|
||
"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": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"M = 2 gives\n",
|
||
"tOut = [0.0010804971680045128, 0.0005500782281160355, 0.0005198698490858078, 0.0005232971161603928, 0.0006225360557436943]\n",
|
||
"\n",
|
||
"M = 3 gives\n",
|
||
"tOut = [0.0005145492032170296, 0.00047505367547273636, 0.000535009428858757, 0.0003488045185804367, 0.00028009340167045593]\n",
|
||
"\n",
|
||
"M = 4 gives\n",
|
||
"tOut = [0.0002853143960237503, 0.00024075806140899658, 0.00021690316498279572, 0.00021761376410722733, 0.00023322459310293198]\n",
|
||
"\n",
|
||
"M = 5 gives\n",
|
||
"tOut = [0.000910954549908638, 0.0005088187754154205, 0.0005268435925245285, 0.0005177054554224014, 0.0005334652960300446]\n",
|
||
"\n",
|
||
"M = 6 gives\n",
|
||
"tOut = [0.002634277567267418, 0.0023930883035063744, 0.0024030981585383415, 0.002396835945546627, 0.0023911641910672188]\n",
|
||
"\n",
|
||
"M = 7 gives\n",
|
||
"tOut = [0.0092502785846591, 0.008834738284349442, 0.00884521659463644, 0.008864735253155231, 0.008844294585287571]\n",
|
||
"\n",
|
||
"M = 8 gives\n",
|
||
"tOut = [0.034397597424685955, 0.03395726904273033, 0.03396071586757898, 0.03401851560920477, 0.03401240427047014]\n",
|
||
"\n",
|
||
"M = 9 gives\n",
|
||
"tOut = [0.14702877961099148, 0.148259031586349, 0.14664556831121445, 0.14639647398144007, 0.14646969363093376]\n",
|
||
"\n",
|
||
"M = 10 gives\n",
|
||
"tOut = [0.5762987565249205, 0.5738314474001527, 0.5741532389074564, 0.5737893972545862, 0.574077476747334]\n",
|
||
"\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"for M in range(2, 10+1):\n",
|
||
" N = 2**M\n",
|
||
" xk = 2*np.pi*np.arange(N)/N\n",
|
||
" # Using the first equation for yk from the previous exercise.\n",
|
||
" yk = np.exp(20j*xk) + np.exp(40j*xk)\n",
|
||
" tOut = timeit.repeat(stmt=lambda: DFT(yk), number=10, repeat=5)\n",
|
||
" tMean = np.mean(tOut)\n",
|
||
" print(\"M =\", M, \"gives\")\n",
|
||
" print(\"tOut =\", tOut)\n",
|
||
" print()"
|
||
]
|
||
},
|
||
{
|
||
"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
|
||
}
|