03_Scipy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# scipy\n",
    "\n",
    "Stellt verschiedene Untermodule für verschiedene Anwendungszwecke bereit. Eine Auswahl:\n",
    "\n",
    "- `scipy.special` für spezielle Funktionen wie\n",
    "\n",
    "- `scipy.linalg` für Kernfunktionen der linearen Algebra (LGS lösen, Determinante usw.)\n",
    "\n",
    "- `scipy.integrate` für numerische Integration und Lösen von ODEs.\n",
    "\n",
    "- `scipy.optimize` u.a. für Nullstellen- und Minimierungsprobleme.\n",
    "\n",
    "- `scipy.interpolate` zum Interpolieren von Funktionen (Splines etc.)\n",
    "\n",
    "Bemerkung: Der Großteil der Funktionen sind lediglich *Wrapper* für in C oder Fortran geschriebene Routinen. Der Rest ist durch vektorisierte numpy-Funktionen implementiert."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Lineare Algebra mit `scipy.linalg`\n",
    "\n",
    "Eigenwerte, Matrixinverse, Lösen linearer Gleichungssyteme und vieles mehr..\n",
    "\n",
    "**Bemerkung:** Einige Standardfunktionen bzw. -methoden der Linearen Algebra sind (aus Gründen der Rückwärtskompatibilität) sowohl in scipy als auch in numpy vorhanden. **Grundsätzlich sollten jene aus scipy verwendet werden!**\n",
    "\n",
    "- [Offizielle Dokumentation](https://docs.scipy.org/doc/scipy/reference/linalg.html)\n",
    "- [Offizielles Tutorial](https://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[30. 11. 12. 13. 14.]\n",
      " [15. 36. 17. 18. 19.]\n",
      " [20. 21. 42. 23. 24.]\n",
      " [25. 26. 27. 48. 29.]\n",
      " [30. 31. 32. 33. 54.]] \n",
      "\n",
      " 18799999.999999996\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import eigvals, inv, solve, det\n",
    "\n",
    "A = 20*np.eye(5) + np.arange(10,35).reshape(5,5)\n",
    "print(A, \"\\n\\n\", det(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.04893617 -0.00255319 -0.00404255 -0.00553191 -0.00702128]\n",
      " [-0.00425532  0.04478723 -0.00617021 -0.00712766 -0.00808511]\n",
      " [-0.00744681 -0.00787234  0.04170213 -0.0087234  -0.00914894]\n",
      " [-0.0106383  -0.01053191 -0.01042553  0.03968085 -0.01021277]\n",
      " [-0.01382979 -0.01319149 -0.01255319 -0.01191489  0.0387234 ]]\n"
     ]
    }
   ],
   "source": [
    "print(inv(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.00000000e+00  1.24032729e-16 -1.11022302e-16  1.47451495e-17\n",
      "   1.90819582e-17]\n",
      " [-6.24500451e-17  1.00000000e+00  0.00000000e+00  4.51028104e-17\n",
      "  -6.93889390e-18]\n",
      " [-9.36750677e-17 -1.30104261e-16  1.00000000e+00 -3.64291930e-17\n",
      "   3.12250226e-17]\n",
      " [-3.81639165e-17  6.41847686e-17 -5.55111512e-17  1.00000000e+00\n",
      "  -2.42861287e-17]\n",
      " [ 4.16333634e-17 -9.02056208e-17  2.22044605e-16  9.02056208e-17\n",
      "   1.00000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "print(inv(A) @ A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prüfen, ob $A^{-1}A = E$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.all(inv(A) @ A == np.eye(5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Da wir mit Gleitkommazahlen reichnen, macht ein solcher Vergleich natürlich wenig Sinn, besser geeignet ist die Funktion `np.allclose(A,B)`, welche $A$ und $B$ innerhalb einer gewissen Toleranz elementweise auf Gleichheit überprüft"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(inv(A) @ A, np.eye(5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[132.22761571+0.j  17.77238429+0.j  20.        +0.j  20.        +0.j\n",
      "  20.        +0.j]\n"
     ]
    }
   ],
   "source": [
    "lambdas = eigvals(A)\n",
    "print(lambdas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.23617021 -0.01968085 -0.00319149  0.16329787  0.22978723]\n"
     ]
    }
   ],
   "source": [
    "b = np.array([-2, 3 , 4, 8, 10])\n",
    "sol = solve(A,b)\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Optimierung mit `scipy.optimize`\n",
    "\n",
    "Stellt u.A. Funktionen für\n",
    "\n",
    "- lineare und nichtlineare (restringierte) Optimierungsprobleme \n",
    "\n",
    "- Nullstellenprobleme\n",
    "\n",
    "- Kleinste-Quadrate- und Funktionsfittingprobleme\n",
    "\n",
    "bereit.\n",
    "\n",
    "- [Offizielle Dokumentation](https://docs.scipy.org/doc/scipy/reference/optimize.html)\n",
    "- [Offizielles Tutorial](https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Als kleines Beispiel wollen wir die Rosenbrockfunktion mit $n=200$ minimieren und wählen als Startpunkt $x_0 = (\\frac{11}{10}, \\frac{11}{10},\\ldots,\\frac{11}{10})^{\\top}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Möglichkeit 1\n",
    "def rosenbrock1(x):\n",
    "    x = np.asarray(x)\n",
    "    res = 0.0\n",
    "    for i in range(x.size-1):\n",
    "        res += 100*(x[i+1]-x[i]**2)**2 + (1-x[i])**2\n",
    "    return res\n",
    "\n",
    "# Möglichkeit 2 (vektorisiert)\n",
    "def rosenbrock2(x):\n",
    "    x = np.asarray(x)\n",
    "    return np.sum(100*(x[1:]-x[:-1]**2)**2 + (1-x[:-1])**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "27.2 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -r 1 -n 1\n",
    "res1 = minimize(rosenbrock1, x0=1.1*np.ones(200))\n",
    "print(res1.message)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "1.38 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -r 1 -n 1 \n",
    "res2 = minimize(rosenbrock2, x0=1.1*np.ones(200))\n",
    "print(res2.message)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hier sieht man bereits, dass es sich lohnt, Funktionen so effizient wie möglich zu implementieren."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Die Funktion `minimize` gibt ein `OptimizeResult`-Objekt zurück, in welchem alle Informationen in Attributen gespeichert sind, siehe [hier](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.OptimizeResult.html#scipy.optimize.OptimizeResult) für eine Übersicht. Den Wert eines Objektattributes enthält man stets durch `objekt.attribut`.\n"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}