{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Tutorial A11 – Solutions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__1__ Write a function that takes three vectors $\\mathbf{a}$, $\\mathbf{b}$, $\\mathbf{c}$ in 2D (NumPy arrays) and checks whether they form a triangle, i.e. if $a + b = \\pm c$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def is_triangle(a, b, c):\n", " \"\"\"Check if a, b, c form a triangle\n", " \n", " Args:\n", " a, b, c (numpy.ndarray): Input vectors\n", " \n", " Returns:\n", " True or False \n", " \"\"\"\n", " \n", " if np.allclose(a + b, np.absolute(c)):\n", " return True\n", " return False" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "is_triangle(\n", " np.array([0, 1]), # a\n", " np.array([1, 0]), # b\n", " np.array([1, 1]) # c\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__2__ Consider the vector $\\mathbf{v} = \\begin{pmatrix}0& 0& 1\\end{pmatrix}^\\mathrm{T}$ and the rotation matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{equation*}\n", "\\mathbf{R} = \\begin{pmatrix} 1 & 0 & 0\\\\\n", "0 & 0 & -1\\\\\n", "0 & 1 & 0\n", "\\end{pmatrix}\n", "\\end{equation*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " - Rotate $\\mathbf{v}$ by $\\mathbf{R}$ to obtain $\\mathbf{v}^\\prime$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0 -1 0]\n" ] } ], "source": [ "v = np.array([0, 0, 1])\n", "R = np.array([[1, 0, 0],\n", " [0, 0, -1],\n", " [0, 1, 0]])\n", "v_ = np.dot(R, v)\n", "print(v_)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " - Compute the scalar product of $\\mathbf{v}$ and $\\mathbf{v}^\\prime$" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0\n" ] } ], "source": [ "print(np.dot(v, v_))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " - By which angle did we rotate along which axis?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*90° around x-axis*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{equation*}\n", "\\mathbf{R} = \\begin{pmatrix} 1 & 0 & 0\\\\\n", "0 & \\cos \\theta & -\\sin \\theta\\\\\n", "0 & \\sin \\theta & \\cos \\theta\n", "\\end{pmatrix}\n", "\\end{equation*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__3__ Consider the second Laplace spherical harmonic\n", "\n", "$$Y_1^{-1}(\\theta, \\phi) = \\frac{1}{2}\\sqrt{\\frac{3}{2\\pi}}\\sin\\theta\\mathrm{e}^{-i\\phi}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " - Compute $Y$ on a grid of 10 points per dimension ($0 \\leq \\theta \\leq \\pi$, $0 \\leq \\phi < 2\\pi$)." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "$\\theta$ [0. 0.34906585 0.6981317 1.04719755 1.3962634 1.74532925\n", " 2.0943951 2.44346095 2.7925268 3.14159265]\n", "$\\phi$ [0. 0.6981317 1.3962634 2.0943951 2.7925268 3.4906585\n", " 4.1887902 4.88692191 5.58505361 6.28318531]\n" ] } ], "source": [ "theta = np.linspace(0, np.pi, 10)\n", "phi = np.linspace(0, 2*np.pi, 10)\n", "print(r\"$\\theta$\", theta)\n", "print(r\"$\\phi$\", phi)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 0.00000000e+00+0.00000000e+00j 1.18165959e-01+0.00000000e+00j\n", " 2.22079358e-01+0.00000000e+00j 2.99206710e-01+0.00000000e+00j\n", " 3.40245317e-01+0.00000000e+00j 3.40245317e-01+0.00000000e+00j\n", " 2.99206710e-01+0.00000000e+00j 2.22079358e-01+0.00000000e+00j\n", " 1.18165959e-01+0.00000000e+00j 4.23108304e-17+0.00000000e+00j]\n", " [ 0.00000000e+00+0.00000000e+00j 9.05203759e-02-7.59556140e-02j\n", " 1.70122659e-01-1.42749860e-01j 2.29205638e-01-1.92326366e-01j\n", " 2.60643034e-01-2.18705474e-01j 2.60643034e-01-2.18705474e-01j\n", " 2.29205638e-01-1.92326366e-01j 1.70122659e-01-1.42749860e-01j\n", " 9.05203759e-02-7.59556140e-02j 3.24119765e-17-2.71968776e-17j]\n", " [ 0.00000000e+00+0.00000000e+00j 2.05193034e-02-1.16370752e-01j\n", " 3.85636759e-02-2.18705474e-01j 5.19567000e-02-2.94661088e-01j\n", " 5.90829793e-02-3.35076226e-01j 5.90829793e-02-3.35076226e-01j\n", " 5.19567000e-02-2.94661088e-01j 3.85636759e-02-2.18705474e-01j\n", " 2.05193034e-02-1.16370752e-01j 7.34719860e-18-4.16680338e-17j]\n", " [ 0.00000000e+00-0.00000000e+00j -5.90829793e-02-1.02334722e-01j\n", " -1.11039679e-01-1.92326366e-01j -1.49603355e-01-2.59120612e-01j\n", " -1.70122659e-01-2.94661088e-01j -1.70122659e-01-2.94661088e-01j\n", " -1.49603355e-01-2.59120612e-01j -1.11039679e-01-1.92326366e-01j\n", " -5.90829793e-02-1.02334722e-01j -2.11554152e-17-3.66422540e-17j]\n", " [ 0.00000000e+00-0.00000000e+00j -1.11039679e-01-4.04151381e-02j\n", " -2.08686334e-01-7.59556140e-02j -2.81162338e-01-1.02334722e-01j\n", " -3.19726014e-01-1.16370752e-01j -3.19726014e-01-1.16370752e-01j\n", " -2.81162338e-01-1.02334722e-01j -2.08686334e-01-7.59556140e-02j\n", " -1.11039679e-01-4.04151381e-02j -3.97591751e-17-1.44711563e-17j]\n", " [-0.00000000e+00+0.00000000e+00j -1.11039679e-01+4.04151381e-02j\n", " -2.08686334e-01+7.59556140e-02j -2.81162338e-01+1.02334722e-01j\n", " -3.19726014e-01+1.16370752e-01j -3.19726014e-01+1.16370752e-01j\n", " -2.81162338e-01+1.02334722e-01j -2.08686334e-01+7.59556140e-02j\n", " -1.11039679e-01+4.04151381e-02j -3.97591751e-17+1.44711563e-17j]\n", " [-0.00000000e+00+0.00000000e+00j -5.90829793e-02+1.02334722e-01j\n", " -1.11039679e-01+1.92326366e-01j -1.49603355e-01+2.59120612e-01j\n", " -1.70122659e-01+2.94661088e-01j -1.70122659e-01+2.94661088e-01j\n", " -1.49603355e-01+2.59120612e-01j -1.11039679e-01+1.92326366e-01j\n", " -5.90829793e-02+1.02334722e-01j -2.11554152e-17+3.66422540e-17j]\n", " [ 0.00000000e+00+0.00000000e+00j 2.05193034e-02+1.16370752e-01j\n", " 3.85636759e-02+2.18705474e-01j 5.19567000e-02+2.94661088e-01j\n", " 5.90829793e-02+3.35076226e-01j 5.90829793e-02+3.35076226e-01j\n", " 5.19567000e-02+2.94661088e-01j 3.85636759e-02+2.18705474e-01j\n", " 2.05193034e-02+1.16370752e-01j 7.34719860e-18+4.16680338e-17j]\n", " [ 0.00000000e+00+0.00000000e+00j 9.05203759e-02+7.59556140e-02j\n", " 1.70122659e-01+1.42749860e-01j 2.29205638e-01+1.92326366e-01j\n", " 2.60643034e-01+2.18705474e-01j 2.60643034e-01+2.18705474e-01j\n", " 2.29205638e-01+1.92326366e-01j 1.70122659e-01+1.42749860e-01j\n", " 9.05203759e-02+7.59556140e-02j 3.24119765e-17+2.71968776e-17j]\n", " [ 0.00000000e+00+0.00000000e+00j 1.18165959e-01+2.89423126e-17j\n", " 2.22079358e-01+5.43937551e-17j 2.99206710e-01+7.32845080e-17j\n", " 3.40245317e-01+8.33360677e-17j 3.40245317e-01+8.33360677e-17j\n", " 2.99206710e-01+7.32845080e-17j 2.22079358e-01+5.43937551e-17j\n", " 1.18165959e-01+2.89423126e-17j 4.23108304e-17+1.03631646e-32j]]\n", "(10, 10)\n" ] } ], "source": [ "def Y(theta, phi):\n", " return 1/2 * np.sqrt(3 / (2*np.pi)) * np.sin(theta) * np.exp(-1j * phi)\n", "\n", "T, P = np.meshgrid(theta, phi)\n", "\n", "y = Y(T, P) \n", "print(y)\n", "print(y.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__4__ Given 1000 randomly places points in 2D, compute all the pairwise distances according to the formula below.\n", "\n", "$d_{ij} = \\sqrt{(x_i-x_j)^2+(y_i-y_j)^2}$ with $i, j \\in (1,2,...,1000)$\n", "\n", "Use three different Methods for this and compare the times:\n", "\n", " - __a)__ A naive double for loop\n", "\n", "How many distances are larger than 0.5?" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "import time\n", "import numpy as np\n", "\n", "samples = 1000 # number of points\n", "\n", "p = (\n", " np.random.random(samples), # x\n", " np.random.random(samples), # y\n", ")\n", "x = p[0]\n", "y = p[1]" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Double for loop needs 1.968029 seconds\n", "Number of elements larger than 0.5: 523098\n" ] } ], "source": [ "t0 = time.time()\n", "\n", "D = np.zeros((samples, samples))\n", "# Array container to store the distances\n", "\n", "for i in range(samples):\n", " for j in range(samples):\n", " if i >= j:\n", " # We make use of the symmetry\n", " # and only calculate whats necessary\n", " continue\n", "\n", " D[i, j] = D[j, i] = np.sqrt((x[i] - x[j])**2 + (y[i] - y[j])**2)\n", " \n", "print(f\"Double for loop needs {time.time() - t0:.6f} seconds\")\n", "print(\"Number of elements larger than 0.5: \", (D > 0.5).sum())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " - __Advanced b)__ A vectorized version using meshgrids" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Meshgrid needs 0.025411 seconds\n", "Number of elements larger than 0.5: 523098\n" ] } ], "source": [ "t0 = time.time()\n", "\n", "X, Y = np.meshgrid(x, y)\n", "# Get all combinations of x and y\n", "\n", "D = np.sqrt((X.T - X)**2 + (Y.T - Y)**2)\n", "\n", "print(f\"Meshgrid needs {time.time() - t0:.6f} seconds\")\n", "print(\"Number of elements larger than 0.5: \", (D > 0.5).sum())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " - __Advanced c)__ Outer subtraction" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Substract.outer needs 0.020332 seconds\n", "Number of elements larger than 0.5: 523098\n" ] } ], "source": [ "t0 = time.time()\n", "\n", "D = np.sqrt(\n", " np.subtract.outer(x, x)**2 + np.subtract.outer(y, y)**2\n", " )\n", "\n", "print(f\"Substract.outer needs {time.time() - t0:.6f} seconds\")\n", "print(\"Number of elements larger than 0.5: \", (D > 0.5).sum())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Advanced d)__: The `cdist` function from the `scipy` module (`scipy.spatial.distance.cdist`)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "from scipy.spatial.distance import cdist\n", "\n", "P = np.array(p).T\n", "# We need our points in a different format, though" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Substract.outer needs 0.009225 seconds\n", "Number of elements larger than 0.5: 523098\n" ] } ], "source": [ "t0 = time.time()\n", "\n", "D = cdist(P, P)\n", "\n", "print(f\"Scipy's cdist needs {time.time() - t0:.6f} seconds\")\n", "print(\"Number of elements larger than 0.5: \", (D > 0.5).sum())" ] } ], "metadata": { "anaconda-cloud": {}, "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.8" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }