05_Gradientenabstiegsverfahren.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> Gradientenabstiegsverfahren</h1>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Motivation</h2>\n",
    "\n",
    "<font size=\"3\" face=\"Verdana\">\n",
    "<p style=\"text-align:justify\">Beim Gradientenabstiegsverfahren nähert man sich, ausgehend von geschätzten Startparametern, schrittweise dem Optimum an. Dabei geht man immer in die Richtung des steilsten Abstiegs also des negativen Gradienten. Im Gegensatz zum Gauß-Newton Verfahren muss nicht die Pseudo-Inverse der Jakobi-Matrix berechnet werden und ist daher effizienter.</p>\n",
    "</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Verfahren</h2>\n",
    "\n",
    "<font size=\"3\" face=\"Verdana\">\n",
    "<p style=\"text-align:justify\">Bevor das Gradientenabstiegsverfahren angewendet werden kann, muss eine Funktion gefunden werden, die das Verhalten der Kurve bestimmt:</p>\n",
    "    \n",
    "\\begin{align}\n",
    "f(p)\n",
    "\\end{align}\n",
    "\n",
    "<p style=\"text-align:justify\">Daraus lässt sich dann eine Fehlerfunktion bestimmen, die minimiert werden muss um die optimalen Parameter zu finden:</p>\n",
    "\n",
    "\\begin{align}\n",
    "r = y - f(p) = F(p) = [F_1(p_1,...,p_n),...,F_m(p_1,...,p_n)]^T \n",
    "\\end{align}\n",
    "\n",
    "<p style=\"text-align:justify\">$F_1$ bis $F_m$ unterscheiden sich dabei nur in den eingesetzten Messpunkten.</p>\n",
    "\n",
    "<p style=\"text-align:justify\">Außerdem müssen noch die $n$ Startparameter geschätzt werden:</p>\n",
    "\n",
    "\\begin{align}\n",
    "p^{(0)} = [p_1^{(0)},...,p_n^{(0)}]\n",
    "\\end{align}\n",
    "\n",
    "<h4>1. Berechnen der Fehlerfunktion mit den Startparametern: </h4>\n",
    "\n",
    "<p style=\"text-align:justify\">Im ersten Schritt werden die Werte der Fehlerfunktion mithilfe der geschätzten Startparameter und den $m$ Messpunkten berechnet:</p>\n",
    "\n",
    "\\begin{align}\n",
    "r = \n",
    "\\begin{bmatrix} \n",
    "       r_{1}   \\\\\n",
    "       \\dots   \\\\\n",
    "       r_{m}   \\\\\n",
    "       \\end{bmatrix} = \n",
    "       \\begin{bmatrix} \n",
    "       F_{1}(p_1^{(0)},...,p_n^{(0)})   \\\\\n",
    "       \\dots   \\\\\n",
    "       F_{m}(p_1^{(0)},...,p_n^{(0)})   \\\\\n",
    "       \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "<h4>2. Berechnen der Jakobi-Matrix für die Fehlerfunktion:</h4>\n",
    "\n",
    "<p style=\"text-align:justify\">Dannach wird die Jakobi-Matrix der Fehlerfunktion berechnet.\n",
    "Dazu werden zuerst die partiellen Ableitungen der Fehlerfunktion nach den Parametern $p_1$ bis $p_n$ gebildet und dann durch Einsetzen der Startparameter und Messpunkte die Einträge der Jakobi-Matrix berechnet.</p>\n",
    "\n",
    "\\begin{align}\n",
    "       J =\n",
    "       \\begin{bmatrix} \n",
    "       \\frac{\\partial r_{1}}{\\partial p_{1}} & \\ldots  & \\frac{\\partial r_{1}}{\\partial p_{n}}  \\\\\n",
    "       \\vdots & \\ddots & \\vdots  \\\\\n",
    "       \\frac{\\partial r_{m}}{\\partial p_{1}} & \\ldots & \\frac{\\partial r_{n}}{\\partial p_{n}}  \\\\\n",
    "       \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "<h4>3. Aktualisieren der Parameter</h4>\n",
    "\n",
    "<p style=\"text-align:justify\">Nun kann mithilfe der Jakobi-Matrix und den Werten der Fehlerfunktion der Kostengradient berechnet werden:</p>\n",
    "    \n",
    "\\begin{align}\n",
    "C_{gradient} = J^{T(s)} * r\n",
    "\\end{align}\n",
    "\n",
    "<p style=\"text-align:justify\">und dann die Parameter aktualisiert werden:</p>\n",
    "\n",
    "\\begin{align}\n",
    "p^{(s+1)} = p^{(s)} - \\gamma * C_{gradient}\n",
    "\\end{align}\n",
    "\n",
    "<p style=\"text-align:justify\">$\\gamma$ gibt hier die Schrittweite an. Es kann konstant gewählt sein, oder in jedem Schritt neu berechnet werden.\n",
    "Wird $\\gamma$ sehr klein gewählt konvergiert das Verfahren nur sehr langsam. Wird es dagegen sehr groß gewählt, kann der Schritt das Optimum überspringen und im schlimmsten Fall divergieren (Siehe Bild 1). Daher sollte $\\gamma$ sorgfältig gewählt werden.</p>\n",
    "\n",
    "<figure>\n",
    "    <img src=\"Stepsize.png\" style=\"width:40%\">\n",
    "    <figcaption style=\"text-align:center\"> Bild 2 </figcaption>\n",
    "</figure>\n",
    "\n",
    "<p style=\"text-align:justify\">Das $\\gamma$ kann zum Beispiel mit der Barzilai-Borwein Methode in jedem Schritt folgendermaßen angepasst werden:</p>\n",
    "\n",
    "\\begin{align}\n",
    "\\gamma_n = \\frac{|(p_n-p_{n-1})^T * (C_{gradient, n}-C_{gradient, n-1})|}{\\|(C_{gradient, n}-C_{gradient, n-1})\\|^2}\n",
    "\\end{align}\n",
    "\n",
    "<p style=\"text-align:justify\">Diese Möglichkeit ist auch im Code weiter unten implementiert.</p>\n",
    "\n",
    "<br />\n",
    "<p style=\"text-align:justify\">Mit den neuen Parametern $p_1,...,p_n$ werden nun die Schritte 1 bis 3 so lange wiederholt bis die Fehlerfunktion klein genug ist, also $\\|r\\| < \\varepsilon$ </p>\n",
    "</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Beispiel </h2>\n",
    "\n",
    "<font size=\"3\" face=\"Verdana\">\n",
    "<p style=\"text-align:justify\">Gegeben sind vier Messpunkte (Tabelle 1), deren Lage duch einen Kreis beschrieben werden kann (siehe Bild 2).</p>\n",
    "\n",
    "<div class=\"container\">\n",
    "    \n",
    "   <div style=\"float:left;width:30%\">\n",
    "        <table style=\"width:30%\">\n",
    "          <caption align=\"bottom\">Tabelle 1: Messpunkte</caption>\n",
    "          <tr>\n",
    "            <th>x1</th>\n",
    "            <th>x2</th>\n",
    "          </tr>\n",
    "          <tr>\n",
    "            <td>3</td>\n",
    "            <td>4</td>\n",
    "          </tr>\n",
    "          <tr>\n",
    "            <td>3.982</td>\n",
    "            <td>2.279</td>\n",
    "          </tr>\n",
    "          <tr>\n",
    "            <td>5</td>\n",
    "            <td>6</td>\n",
    "          </tr>\n",
    "          <tr>\n",
    "            <td>6.275</td>\n",
    "            <td>5.541</td>\n",
    "          </tr>\n",
    "        </table>\n",
    "   </div>\n",
    "    \n",
    "   <div style=\"float:left;width:30%\">\n",
    "        <figure>\n",
    "            <img src=\"Kreis.png\" style=\"width:50%\">\n",
    "            <figcaption style=\"text-align:center\"> Bild 2 </figcaption>\n",
    "        </figure>\n",
    "   </div>\n",
    "    \n",
    "</div>\n",
    "\n",
    "<p style=\"text-align:justify\">Es sollen die Werte für z1, z2 und R (siehe Bild 3) optimiert werden.</p>\n",
    "\n",
    "<figure>\n",
    "    <img src=\"KreisVariablen.png\" style=\"width:20%\">\n",
    "    <figcaption style=\"text-align:center\"> Bild 3 </figcaption>\n",
    "</figure>\n",
    "\n",
    "<p style=\"text-align:justify\">Ein Kreis kann beschrieben werden durch: </p>\n",
    "\n",
    "\\begin{align}\n",
    "f(x_{1}, x_{2}) = (z_{1} - x_{1})^{2} + (z_{2} - x_{2})^{2} - R^{2}\n",
    "\\end{align}\n",
    "\n",
    "<p style=\"text-align:justify\">oder in Bezug auf die zu optimierenden Parameter durch:</p>\n",
    "\n",
    "\\begin{align}\n",
    "f(z_{1}, z_{2}, R) = (z_{1} - x_{1})^{2} + (z_{2} - x_{2})^{2} - R^{2}\n",
    "\\end{align}\n",
    "\n",
    "<p style=\"text-align:justify\">Weiteres Umformen ergibt die Fehlerfunktion:</p>\n",
    "\n",
    "\\begin{align}\n",
    "R^{2} = (z_{1} - x_{1})^{2} + (z_{2} - x_{2})^{2}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "R = \\sqrt{(z_{1} - x_{1})^{2} + (z_{2} - x_{2})^{2}}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "R = \\|z-x\\|\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "F(z_{1}, z_{2}, R) = R - \\|z-x\\|\n",
    "\\end{align}\n",
    "\n",
    "<p style=\"text-align:justify\">Als Startparameter werden die folgenden Werte festgelegt:</p>\n",
    "\n",
    "\\begin{align}\n",
    "p^{0} = [5.5, 4.3, 1.9]^{T}\n",
    "\\end{align}\n",
    "\n",
    "<h4>1. Berechnen der Fehlerfunktion mit den Startparametern:</h4>\n",
    "\n",
    "<p style=\"text-align:justify\">Für jeden der vier Messpunkte wird mithilfe der Startparameter der Wert der Fehlerfunktion berechnet und als Vektor dargestellt:</p>\n",
    "\n",
    "\\begin{align}\n",
    "r_{i} = (1.9)^{2} - \\sqrt{(5.5-x_{i1})^{2} + (4.3-x_{i2})^{2}}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "       F(p^{(0)}) =\n",
    "       \\begin{bmatrix} \n",
    "       r_{1}   \\\\\n",
    "       r_{2}   \\\\\n",
    "       r_{3}   \\\\\n",
    "       r_{4}   \\\\\n",
    "       \\end{bmatrix} =\n",
    "       \\begin{bmatrix} \n",
    "       1.9 - \\sqrt{(5.5-3)^{2} + (4.3-4)^{2}}   \\\\\n",
    "       1.9 - \\sqrt{(5.5-3.982)^{2} + (4.3-2.279)^{2}}   \\\\\n",
    "       1.9 - \\sqrt{(5.5-5)^{2} + (4.3-6)^{2}}   \\\\\n",
    "       1.9 - \\sqrt{(5.5-6.275)^{2} + (4.3-5.541)^{2}}   \\\\\n",
    "       \\end{bmatrix} =\n",
    "       \\begin{bmatrix} \n",
    "       -0.61794   \\\\\n",
    "       -0.62760   \\\\\n",
    "       0.12800   \\\\\n",
    "       0.43688   \\\\\n",
    "       \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "<h4>2. Berechnen der Jacobi-Matrix für die Fehlerfunktion:</h4>\n",
    "\n",
    "<p style=\"text-align:justify\">Mithilfe der Ableitungsregeln:</p>\n",
    "\n",
    "\\begin{align}\n",
    "f(x) = \\sqrt{x} \\longrightarrow f'(x) = \\frac{1}{2 \\sqrt{x}}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "f(x) = g(h(x)) \\longrightarrow f'(x) = g'(h(x)) * h*(x)\n",
    "\\end{align}\n",
    "\n",
    "<p style=\"text-align:justify\">kann die Fehlerfunktion jeweils nach einem der drei Parameter abgeleitet werden:</p>\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial r_{i}}{\\partial z_{1}} = - \\frac{1}{2 \\sqrt{(z_{1}-x_{1})^{2} + (z_{2}-x_{2})^{2}}} * 2(z_{1}-x_{1}) = - \\frac{(z_{1}-x_{1})}{\\|z-x\\|}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial r_{i}}{\\partial z_{2}} = - \\frac{(z_{2}-x_{2})}{\\|z-x\\|}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial r_{i}}{\\partial r} = 1\n",
    "\\end{align}\n",
    "\n",
    "<p style=\"text-align:justify\">Daraus können dann die Einträge der Jacobi-Matrix berechnet werden:</p>\n",
    "\n",
    "\\begin{align}\n",
    "       J =\n",
    "       \\begin{bmatrix} \n",
    "       \\frac{\\partial r_{1}}{\\partial p_{1}} & \\ldots  & \\frac{\\partial r_{1}}{\\partial p_{n}}  \\\\\n",
    "       \\vdots & \\ddots & \\vdots  \\\\\n",
    "       \\frac{\\partial r_{m}}{\\partial p_{1}} & \\ldots & \\frac{\\partial r_{n}}{\\partial p_{n}}  \\\\\n",
    "       \\end{bmatrix} =\n",
    "       \\begin{bmatrix} \n",
    "       \\frac{\\partial r_{1}}{\\partial z_{1}} & \\frac{\\partial r_{1}}{\\partial z_{2}} & \\frac{\\partial r_{1}}{\\partial r}   \\\\\n",
    "       \\frac{\\partial r_{2}}{\\partial z_{1}} & \\frac{\\partial r_{2}}{\\partial z_{2}} & \\frac{\\partial r_{2}}{\\partial r}   \\\\\n",
    "       \\frac{\\partial r_{3}}{\\partial z_{1}} & \\frac{\\partial r_{3}}{\\partial z_{2}} & \\frac{\\partial r_{3}}{\\partial r}   \\\\\n",
    "       \\frac{\\partial r_{4}}{\\partial z_{1}} & \\frac{\\partial r_{4}}{\\partial z_{2}} & \\frac{\\partial r_{4}}{\\partial r}   \\\\\n",
    "       \\end{bmatrix} = \n",
    "       \\begin{bmatrix} \n",
    "       -\\frac{z_1-x_{11}}{\\|z-x_1\\|} & -\\frac{z_2-x_{12}}{\\|z-x_2\\|} & 1   \\\\\n",
    "       -\\frac{z_1-x_{21}}{\\|z-x_2\\|} & -\\frac{z_2-x_{22}}{\\|z-x_2\\|} & 1   \\\\\n",
    "       -\\frac{z_1-x_{31}}{\\|z-x_3\\|} & -\\frac{z_2-x_{32}}{\\|z-x_3\\|} & 1   \\\\\n",
    "       -\\frac{z_1-x_{41}}{\\|z-x_4\\|} & -\\frac{z_2-x_{42}}{\\|z-x_4\\|} & 1   \\\\\n",
    "       \\end{bmatrix} = \n",
    "       \\begin{bmatrix} \n",
    "       -\\frac{5.5-3}{\\sqrt{(5.5-3)^2+(4.3-4)^2}} & -\\frac{4.3-4}{\\sqrt{(5.5-3)^2+(4.3-4)^2}} & 1   \\\\\n",
    "       -\\frac{5.5-3.982}{\\sqrt{(5.5-3.982)^2+(4.3-2.279)^2}} & -\\frac{4.3-2.279}{\\sqrt{(5.5-3)^2+(4.3-2.279)^2}} & 1   \\\\\n",
    "       -\\frac{5.5-5}{\\sqrt{(5.5-5)^2+(4.3-6)^2}} & -\\frac{4.3-6}{\\sqrt{(5.5-5)^2+(4.3-6)^2}} & 1   \\\\\n",
    "       -\\frac{5.5-6.275}{\\sqrt{(5.5-6.275)^2+(4.3-5.541)^2}} & -\\frac{4.3-5.541}{\\sqrt{(5.5-6.275)^2+(4.3-5.541)^2}} & 1   \\\\\n",
    "       \\end{bmatrix} =\n",
    "       \\begin{bmatrix} \n",
    "       -0.99288 & -0.11915 & 1   \\\\\n",
    "       -0.60057 & -0.79957 & 1   \\\\\n",
    "       -0.28217 & 0.95937 & 1   \\\\\n",
    "       0.52969 & 0.84819 & 1   \\\\\n",
    "       \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "<h4>3. Aktualisieren der Parameter:</h4>\n",
    "\n",
    "<p style=\"text-align:justify\">Der Schrittweite $\\gamma$ wird in diesem Beispiel der Einfachkeit halber ein fester Wert von $0.1$ zugeordnet</p>\n",
    "\n",
    "<p style=\"text-align:justify\">Mit den Werten der Fehlerfunktion, der Jacobi-Matrix und der Schrittweite $\\gamma$ können nun die Parameter aktualisiert werden:</p>\n",
    "\n",
    "\\begin{align}\n",
    "C_{gradient} = J^{T(s)} * r\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "       C_{gradient} =\n",
    "       \\begin{bmatrix} \n",
    "       -0.99288 & -0.60057 & -0.28217 & 0.52969   \\\\\n",
    "       -0.11915 & -0.79957 & 0.95937 & 0.84819   \\\\\n",
    "       1 & 1 & 1 & 1   \\\\\n",
    "       \\end{bmatrix} \n",
    "       \\begin{bmatrix} \n",
    "       -0.61794    \\\\\n",
    "       -0.62760    \\\\\n",
    "       0.12800    \\\\\n",
    "       0.43688    \\\\\n",
    "       \\end{bmatrix} =\n",
    "       \\begin{bmatrix} \n",
    "       1.18575    \\\\\n",
    "       1.06879    \\\\\n",
    "       -0.68066    \\\\\n",
    "       \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "p^{(s+1)} = p^{(s)} - \\gamma * C_{gradient}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "       p^{(s+1)} =\n",
    "       \\begin{bmatrix} \n",
    "       5.5   \\\\\n",
    "       4.3   \\\\\n",
    "       1.9   \\\\\n",
    "       \\end{bmatrix} - 0.1\n",
    "       \\begin{bmatrix} \n",
    "       1.18575    \\\\\n",
    "       1.06879    \\\\\n",
    "       -0.68066    \\\\\n",
    "       \\end{bmatrix} =\n",
    "       \\begin{bmatrix} \n",
    "       5.381427    \\\\\n",
    "       4.1926    \\\\\n",
    "       1.968066    \\\\\n",
    "       \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "<p style=\"text-align:justify\">Nun wird mit den aktualisierten Parametern wieder bei Schritt 1 begonnen und so lange wiederholt bis die Fehlerfunktion ausreichend klein ist: $\\|r\\| < \\varepsilon$</p>\n",
    "\n",
    "</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Übungen</h2>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font size=\"3\" face=\"Verdana\">\n",
    "<h4>1. Code-Beispiel</h4>\n",
    "<p style=\"text-align:justify\">Im folgenden befindet sich eine implementierung des Gradientenabstiegsverfahrens.</p>\n",
    "</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# Zielfunktion\n",
    "def function(p, x):\n",
    "    return p[2] - math.sqrt((p[0]-x[0])**2 + (p[1]-x[1])**2)\n",
    "    \n",
    "    \n",
    "# Partielle Ableitungen der Zielfunktion\n",
    "def f_par_deriv_z1(p, x):\n",
    "    return -(p[0]-x[0]) / math.sqrt((p[0]-x[0])**2 + (p[1]-x[1])**2)\n",
    "\n",
    "def f_par_deriv_z2(p, x):\n",
    "    return -(p[1]-x[1]) / math.sqrt((p[0]-x[0])**2 + (p[1]-x[1])**2)\n",
    "\n",
    "def f_par_deriv_R(p,x):\n",
    "    return 1\n",
    "\n",
    "nabla_function = [f_par_deriv_z1, f_par_deriv_z2, f_par_deriv_R]\n",
    "\n",
    "\n",
    "# Berechnen der Jakobi Matrix\n",
    "def jacobian(p, points, nabla_f):\n",
    "    rows = np.shape(points)[0] \n",
    "    cols = np.shape(p)[0]\n",
    "    J = np.zeros((rows,cols))\n",
    "    for row in range(rows):\n",
    "        for col in range(cols):\n",
    "            J[row][col] = nabla_f[col](p, points[row])\n",
    "    return J\n",
    "\n",
    "\n",
    "# Berechnen des Residual\n",
    "def residual(p, points):\n",
    "    rows = np.shape(points)[0]\n",
    "    # Residual Vektor initialisieren\n",
    "    r = np.zeros(rows)\n",
    "    for row in range(rows):\n",
    "        r[row] = function(p, points[row])\n",
    "    return r\n",
    "\n",
    "\n",
    "# Gradientenabstiegsverfahren mit festem gamma\n",
    "def Gradientenabstieg_festesGamma(nabla_f, p0, x, gamma):\n",
    "    # Residual und Parameter initialisieren\n",
    "    r = residual(p0, x)\n",
    "    p = p0\n",
    "    \n",
    "    # Abbruchbedingungen\n",
    "    MAX_ITER = 1000\n",
    "    last_r = r\n",
    "    MIN_DELTA = 0.000001\n",
    "    EPSILON = 0.000001\n",
    "    \n",
    "    i=0\n",
    "    while i < MAX_ITER:\n",
    "        i = i+1\n",
    "        \n",
    "        # aktuellen Kreis darstellen\n",
    "        circle = plt.Circle((p[0], p[1]), radius=p[2], fc='none', ec=[max(0.9-0.05*i, 0), max(0.9-0.05*i, 0), max(0.9-0.05*i, 0)])\n",
    "        plt.gca().add_patch(circle)\n",
    "        \n",
    "        # neue Jakobi Matrix berechnen\n",
    "        J = jacobian(p, x, nabla_f)\n",
    "        Jt = np.matrix.transpose(J)\n",
    "        \n",
    "        # Parameter aktualisieren\n",
    "        p = p - gamma*(np.dot(Jt, r))\n",
    "        \n",
    "        # Residual aktualisieren\n",
    "        r = residual(p, x)\n",
    "        \n",
    "        #Abbruchbedingungen prüfen\n",
    "        if np.linalg.norm(r) < EPSILON or np.linalg.norm(r - last_r) < MIN_DELTA:\n",
    "            break\n",
    "            \n",
    "        last_r = r\n",
    "        \n",
    "    # Parameter und Anzahl der Iterationen zurückgeben\n",
    "    ret = [p, i]\n",
    "    return ret\n",
    "    \n",
    "    \n",
    "# Gradientenabstiegsverfahren mit angepasstem gamma mithilfe der Barzilai-Borwein Methode\n",
    "def Gradientenabstieg(nabla_f, p0, x, gamma):\n",
    "    # Residual und Parameter initialisieren\n",
    "    r = residual(p0, x)\n",
    "    p = p0\n",
    "    J = jacobian(p, x, nabla_f)\n",
    "    Jt = np.matrix.transpose(J)\n",
    "    c_gradient = np.dot(Jt, r)\n",
    "    \n",
    "    # Abbruchbedingungen\n",
    "    MAX_ITER = 1000\n",
    "    last_r = r\n",
    "    MIN_DELTA = 0.000001\n",
    "    EPSILON = 0.000001\n",
    "    \n",
    "    i=0\n",
    "    while i < MAX_ITER:\n",
    "        i = i+1\n",
    "        \n",
    "        # aktuellen Kreis darstellen\n",
    "        circle = plt.Circle((p[0], p[1]), radius=p[2], fc='none', ec=[max(0.9-0.05*i, 0), max(0.9-0.05*i, 0), max(0.9-0.05*i, 0)])\n",
    "        plt.gca().add_patch(circle)\n",
    "\n",
    "        # Parameter aktualisieren\n",
    "        last_p = p\n",
    "        p = p - gamma*c_gradient\n",
    "        print(gamma)\n",
    "        # Residual aktualisieren\n",
    "        r = residual(p, x)\n",
    "        \n",
    "        # Jakobi Matrix aktualisieren\n",
    "        J = jacobian(p, x, nabla_f)\n",
    "        Jt = np.matrix.transpose(J)\n",
    "        last_c_gradient = c_gradient\n",
    "        c_gradient = np.dot(Jt, r)\n",
    "        \n",
    "        # Gamma aktualisieren\n",
    "        delta_p = p - last_p\n",
    "        delta_c_gradient = c_gradient - last_c_gradient\n",
    "        norm_delta_c_gradient = np.linalg.norm(delta_c_gradient)\n",
    "        if norm_delta_c_gradient != 0:\n",
    "            gamma = np.absolute(np.dot(delta_p, delta_c_gradient)) / norm_delta_c_gradient**2\n",
    "            \n",
    "        #Abbruchbedingungen prüfen\n",
    "        if np.linalg.norm(r) < EPSILON or np.linalg.norm(r - last_r) < MIN_DELTA:\n",
    "            break\n",
    "            \n",
    "        last_r = r\n",
    "    \n",
    "    # Parameter und Anzahl der Iterationen zurückgeben\n",
    "    ret = [p, i]\n",
    "    return ret\n",
    "            \n",
    "\n",
    "# Funktion um verrauschte Messpunkte zu generieren\n",
    "def getPoints(x, std, R, z_eins, z_zwei):\n",
    "    for i in range(np.size(x,0)):\n",
    "        x[i][1] = (- math.sqrt(R**2 - (z_eins - x[i][0])**2) + z_zwei) + np.random.normal(0, std) \n",
    "    return x\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font size=\"3\" face=\"Verdana\">\n",
    "<p style=\"text-align:justify\">Welche Auswirkungen können Sie erkennen, wenn Sie die Startparameter verändern?</p>\n",
    "<p style=\"text-align:justify\">Wie wirkt sich eine Veränderung der Messpunkte und deren Qualität aus?</p>\n",
    "<p style=\"text-align:justify\">Was bewirkt eine Änderung der Schrittweite Gamma für das Verfahren mit einem festen Gamma?</p>\n",
    "<p style=\"text-align:justify\">Wie verändert sich das Ergebnis, wenn das Verfahren ohne festes Gamma verwendet wird?</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAD3CAYAAADrGWTVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9eXxcd33v/f6d2TWrRvtqSd7tbBAnBEhyKdCnabmQXqAFLvQBHtbC7cZz4dVLywVCt4ct7SUsDXtCSylQaAglhIYQtsSJbeJFtmxLsqyRRsvs+z6/54/xOWi3bMnSyP69X695SZpzzszv2Gc+8z3fVUgpUSgUCsXWRNvsBSgUCoXi8lEirlAoFFsYJeIKhUKxhVEirlAoFFsYJeIKhUKxhTFv5Js1NzfLvr6+jXxLhUKh2PIcPnw4LKVsWWrbhop4X18fhw4d2si3VCgUii2PEOL8ctuUO0WhUCi2MErEFQqFYgujRFyhUCi2MErEFQqFYgujRFyhUCi2MErEFQqFYgujRFyhUCi2MBuaJ65QbBZSSkqlEoVCgWq1SrVaRUppPIQQCCHQNA0hBCaTCZvNhtlsRgix2ctXKJZFibhiyyOlJJPJkMlkSKfTpFIpUqkU2WyWXC5HoVCgWCwyt3e+Lt7APJHWhVwXb03TsNls2O12GhoacDqdOBwOGhoacDgcxt9K6BWbhRJxxZZCSkkqlWJiYoLp6Wmi0SipVIpSqUSpVDKs7Gq1SrlcplKpGI+5wj1X0HUrXBdtk8lk/NQfuoVuNpuNh9VqxWKxYLFYcLvduFwuvF4vzc3N+Hw+TCbTZv0zKa4hlIgr6ppcLse5c+cMwU4kEmSzWcrlMqVSiXw+P0+8deGuVquGeOvbL4Yu0haLxRBu/TFXzE0mkyHgVqvVEHWbzWa4YCwWC06nE7fbjc/no6OjA5fLhaapMJRifVEirqgrpJRMT09z6tQpxsfHSaVS8wS7UChQKpUoFosUCgXDx60Ltu7rLpfLlMtlY/tcC3yhFa7/1C1x3bo2m82YTCbjeV3ALRYLNpvN2E8X8qWEXXe9nD59Grvdjsfjobm5mebmZjwejxJ1xZpRIq7YdMrlMqOjo5w9e5bJyUnS6TT5fJ5MJrNIsPW/S6US2WzWsLKllIaYF4tFisXiPKt8Jeb6s3Wh1i1rXcznirwuzna73RB1XdgXCrrD4cDhcBjbk8kkMzMz2O127HY7TU1NtLW14fP5sFgsV/qfWnEVokRcsSkUi0VOnTrF2bNnmZ6eplAokMlkyGaz5PN5UqmUIcbFYpFcLmf4vYvFItls1hBq3Y2ii7wu3PqjUqkY77tUQFMX6Lm+77mCPVecdYvcZrPhdDqxWq2G9T5X2K1Wq7GP/mXgcDhIJpOGwDudTmKxGOPj49jtdvx+vyHodrtdBUsVq0KJuGJDmZmZ4emnn+bcuXPk83nS6bQh4JlMhnw+TzabJZvNGlklcwU7l8sZf88VbB0hxDx3yeWgi6emafNcMUIIw12iaRp2ux2n04nFYlkk7BaLBbvdjsvlmifoVqsVu91uZLzoVrrb7SabzRIOh3E6nTidTlpaWmhublYWumJFlIgrrjiVSoWTJ09y5MgRwuEw8XjcsLZzuRy5XM4Q8FwuRzKZNH6fK9ilUumi77VWAZ/7GnMteP15/c5A0zSSySQAVqt1kbDrYu31eg1r3uPxGILe0NBgCLvD4SCfz2O1WnE4HLjdbsxms5GF09TUREtLCy6XS1nnikUoEVdcMZLJJAcPHuT06dOk02kSiYTxM5vNkkqlKBQKZLNZ0um0IeB6bvdqMko2i7lrW0rYHQ4HNpsNl8tFQ0MDVqsVj8djZLA0NDQYD5fLZVj0brebXC6H2Wwml8tht9tJJpOEQiGcTidtbW34/X6VvqgwUCKuWHcmJyf5xS9+wcTEBMlkkmw2SywWI5vNEo/HyWQyhpDPtbwLhcJmL31NzBV2/YsoHo8bwq1b2VarFafTicvlwuFw4PF4jOKhTCZjiHuhUMBkMuFyucjlckSjUTKZzDzr3OFwbOIZK+qBNYu4ECIHHLzw54NSyi+u9TUVW5NwOMyPfvQjxsfHSSaThlinUinDCk8mkyQSCcPiLpfLm73sK4Yu6rpvP5lMEolEDF+4z+fD4XDg9XqNalBd2F0uFy6XC5vNZsQM3G431WoVs9mMlJJoNIrH46GzsxO73b7JZ6vYLNbDEp+UUr5oHV5HsUVJJpM89thjnD17llgsRjweN4RbF/JkMkksFjMCmdcaup+9XC4b7QESiQQNDQ34fD4aGhqMqk/9d13Yi8UiVquVcrlMKpXC6/UaYl6tVkkkEvj9fjo6OrBarZt8poqNZj1EvF0I8QQQAd4jpRybu1EI8Xbg7QC9vb3r8HaKeiGXy/HEE09w4sQJQ7zj8TiRSIR0Om2IeTweJ51O17WPeyPRBV33pScSCSMjRbfK9Z8ej4dsNovT6aRUKmGz2ahWqySTSbxeL6VSCYfDYVjmLS0ttLe3YzYrT+m1wnr8T/dJKcNCiN8Cvgi8ZO5GKeX9wP0ABw4cWHvqgGLTKRaLPPXUUxw+fJhIJEIymSQajRKJREgkEsRiMaOnSS6X2+zl1j16A6+51rnX68XtdpNOp0mn00YKosvlolAoGMKdSqXw+XzkcjmcTidQc2u1t7fT2tqqKkKvAdYs4lLK8IWfPxRCfHrtS1LUM2fPnuWRRx4hFAqRTqcJh8NGT5NwOEwoFCKRSNSFr9sCWKld5PpDAuU5jwJQWe4FNhDdOtcDvIlEAo/HQyqVwuPx4PP5DFdKOp3G4/EYwi2lxGKxIKUkm83S1NREKBRiZmaG7u5u/H6/Sk28ilmTiAshXEBOSlkRQtwAhNdnWYp6I5/P8/3vf5/BwUGSyaQh3rr7RM//LhaLG7Ke3guPTqBjzs+5v7uoibQu2LpYzxV1KxAHpoDggp/676PA7Aack46UkkqlQiwWM9wmuoDrrhU9RbOxsZFsNktjYyOlUgmfz0cikSCXy+H3+5mYmCASidDf36+Khq5S1mqJ7wP+UQiRombkvGPtS1LUG2fPnuXhhx9mdnaWUChENBolGo0SDocN10k+n79i798H3LzgUaAmrnPFdpD5QvzIU09x6/Oet+Jr/8O99/I373nPoi+AfdT8gp3ATiAHHAYOXfh5GJhZx3NcDl3M5+bZz7XCS6WSkZmSyWTQNI1sNovf7yeXy9Hd3U2hUODEiRP09PTQ1NSkrPKrDLEeFW6r5cCBA/LQoUMb9n6KtTHX+o7FYszOzjIzM0MoFCIWixGJRMhms+v6nlbgRRceumDn+bVw6o+pDbxuf/6zn/EHd95prOfAhZ9zhf0/qeXZXsnQrRDCqPxsbGw0OiK6XC78fj9erxeXy0VjY6OR9eJ2u3E4HLS2tpLNZnE4HAwMDCirfIshhDgspTyw5DYl4oqlWGh9z87OEolEmJ6eZnZ2lnQ6vW7v1QT8DvBy4DepWdQ/4tdW70YK9mqZK+y3AncB7cDDwEPU1r++X2/zcTgcNDY2GgMoPB4Pra2t2O122tvbsVqttLW1YTab8fv9aJpGb28vlUqFcrmsrPItxkoirvKQFPMolUo8/PDDHD9+fJ71rf8Mh8Prkiq4i5povwK4CXiMmvi9G5iVkheu+R2uLLffcQfnFny59AvBy6mdwwPAz6id08PU3DvrSS6XM/rL6N0fK5UKLpeLSqWC1+ulWCwavnKv18vk5CQul4u2tjbDV66s8q2PEnGFQSKR4J/+6Z8IBAKEw2FmZmaIRCLMzMwwPT1NJpNZ0+u7gT8A3knN+n4I+DvgoMNBJJvlv635DDaXuaLuFYK7qH1J/S1wEvgM8G1gvUK/Usp57QuKxaKRblgqlUin02iaRiaTMYqFNE1jbGyMvXv3EolEGBwcZPv27bjd7nValWKjUe4UBQCBQIB/+Zd/YWpqilAoxPT0tCHekUhkTdb39cAfAq+h5jv+LPCe736Xl9999/osvs55wS230HboEO8CbqBWTPGPwPg6vocQAp/Ph9/vp7GxkZaWlnmj4SwWC319fUgpaW1tRUrJDTfcQC6XIx6P09XVRVtb2zquSLGeKHeKYkUOHTrEI488wvT0NNPT00xNTa3Z+rYCrwTeBfRTE6391Pzbv79uK98a/PKZZ4zfdwvBO4EjwC+oWeePUkvtWgtLWeV+v598Po8QAq/Xy+joKF1dXQSDQVpbWzl27Bg7duwwnstms/T19Sk/+RZDifg1TLVa5T/+4z84ePAgMzMzBINBgsEgk5OTl+379gLvodZn4QRwLxB+/vP56S9/yUfWd/lbktMX7nydQvBa4G+A+4D/Q+2Lbi2ulmq1SqFQYGpqyigayufzxoSjfD6Ppmn4fD5mZ2dpbGxkdHSU5uZm+vv7OX/+PKdOnWLXrl2qbH8Lof6nrlFyuRzf+MY3GBoaMgR8cnKSQCBg9MS+FOzAHwHvpebrfhEwJCUvXddVXz1kLoj5RCDAq3t7+QtqX34fBL7G2lIVpZRG4ZU+sq5ardLc3EypVELTNKNnue4nT6fT3HjjjQQCAQYHB9m5cycNDQ1rP1HFFUc1VrgGiUQi3H///Rw/fpzJyUnOnz/P2NgY586du2QBNwFvAc4AzwPuBN4iJUN1mBZYj3T39PCUlLxcSt4AvA04Si1zZy1Uq1UymQzBYJCpqSmmpqaMNNGxsTEKhQKBQIBisUgoFKJYLPL000/T3d2N0+lkaGiIaDS69hNUXHGUiF9jzMzM8PnPf57R0VGCwSDnzp3j3LlzjI2NXXLhziupuUxeD7waeJWUnFLifdn8XErkE0/wv4C/ppaiuJZUS32c3MzMDJOTkwSDQSKRCLOzswQCAfL5PMPDw2iaxvT0NCaTiYMHD9Lc3ExbWxtjY2NMT0+vz8kprhhKxK8hgsEgX/ziFzl//jyBQIDh4WFGR0cZHx+/pJ4nL6BWnfiXwJ8AnadOcVCJ97pwx5138j0p+czb3sb91FwrD1HLq79cKpUKoVCIyclJpqenjd9nZmbIZrOcOXMGq9XK5OQkZrOZp59+GqfTSW9vr+FmU9QvSsSvEQKBAF/+8pcZHx9nYmKC0dFRzp07x9TU1KoDmA5qgcpvAn8P/Nv7388PpWT3nj1XcOXXJp+9/34ekJI7Wlp4nFomy3u4/A+s3m98YmLCqLqdmJggGo2STqc5ffr0PCE/ePAgmqYxMDDAzMwMgUBg/U5Osa6owOY1QCAQ4Ktf/SqBQICJiQnOnz/P+Pg4qVRq1a9xO/Al4Glqed8RZXlvCIHZWv/E7ULwReBVwJupxSAuFSklyWSSUqlEqVSiUqkgpURKSbVa5fTp0+zZs4eJiQl6enp45plnuPXWW9m+fTujo6NIKdVglzpEWeJXOcFgkK9+9auGBT42NsbY2NiqBVy3vr9BLfPk9VIqAd8ERqTkpx/8IP/M2qxyKSW5XG5eK4XJyUljCtPQ0BANDQ1MTExgtVp55kKO+/bt24lEIsoir0OUiF/FzMzM8JWvfMUQcD2AudoCntupZUq0ULO+v6vEe1P53x/6EPdJyfOoZa/8jMvzlc8NeOpFXbqQp9NpTp48idPpZHJyEqvVytNPPw3Ajh07DH+6on5QIn6VEolE+NKXvsT58+eZnJy8pAwUM/BxlPVdryy0yv/oMl5DSkmpVDKscV3IY7EYuVyOwcFBPB4Pk5OTNDQ0cPDgQcxmM9u3bzfqChT1gRLxq5B8Pm8IuJ5GeP78+VXNu/QDPwT2oqzvemauVf7/UItXXOqc+7lCPte1EolEKBaLHD16FK/XSyAQwOv1cujQIVwuF/39/UaWi2LzUSJ+lVGtVnnwwQcZGxszLPDVCvh+aoHLZ4Av3n23sr63ACMX2va6gceBS21hJaWkXC4bQj47O0swGGRmpja36MSJE1itViYmJjCbzRw9ehSfz0dPTw+BQGBd+8orLg8l4lcZ3/ve9xgcHCQYDBpBzNW4UF4B/Jha2ff7qlW+/d3vXumlKtaJjJRkvvQlHqH2JfzcSzx+rkWuDwCZmppiYmICKSWBQIBsNkskEiGVSjE8PExLSwtNTU0MDw9v2FxVxdIoEb+KeOaZZ/jFL35BMBgkEAisOoj5fmpNmP4r8JVSSXWx24K88c1v5h4p+RPgEWptfy8FPdg5OztrDMAOhUKcP3+eTCZDqVQiHA5TLBaNSs6enh7sdjtnzpxZl0EhistDifhVwvj4OA899JBRYTc6OnrRW10H8HVqVvjzgJ9eaIqk2Lp8R0peQm0QxV8Bl/J1LKWkUCgYPVaSySTJZJJgMMjExAQej4dgMIjD4eDYsWNkMhl27txJtVpldHT0ypyQ4qIoEb8KSCaTPPjgg0xMTBiBzFgstuIxPmrukxLwX4DTyaQxNV2xtTkmJbdQ67vydS6toq9arZJOp4088mg0SiaTIR6PMzg4SHNzM+fPn8fr9fLMM89QqVTYuXMnqVRKpR5uEkrEtzjlctkop9cDmRdrWtRETcB/CYQ+/nEmQiE1nusqIywlr2powEGtTcKlZK7orWynpqaMCU+JRIJyucyzzz6Lx+NhfHwct9vN4cOHsdls9Pf3MzMzc1HjQbH+rIuICyF2CSFKQojb1+P1FKvn3/7t3xgeHiYYDBo54SuN3GsFfgL8AGh94AFe89rX0tTUtEGrVWwkkUyG9/b1UQa+S63n+2qpVqtEo1EjY0Wf+iSlZHh4mEqlYvw9OjqKz+ejvb2dsbEx8vn8FTojxVKslyX+AeCJdXotxSoZGhri6aefNgKZgUCAUqm07P66gP8r8HunT/PC22+no6NDBTKvYk6fO8d9L3whUeB7XLqQRyIRwuEwkUjEKAjKZDKYTCbC4TAmk4nR0VFSqRSdnZ04nU5GRkZWNCQU68uaRVwIcSswDUysfTmK1VIoFPjWt75ljFQbGxtbMRe8idqQ4n8BPlguo5lMtLW1oWnKo3a185Of/5yHXvUqQsC3uDTXSqlUIhQKEYlEiMfjRKNRQqEQQ0NDdHV1MTIyQkdHB88++yzVapWBgQGKxSJTU1NX6GwUC1mPT/BfAn+33EYhxNuFEIeEEIdUhdf68e1vf5vx8XGmpqYYHR0lkUgsu6+P2jDeh4HXDw0xOzuL2+1W47euIb7xrW/x9895DnlqX+SrDXbqDbNCoRDRaJRkMkk6naZcLnPs2DFMJhNjY2OGRW42m+nt7WV6elq5VTaINYm4EOJlwCEpZWS5faSU90spD0gpD7S0tKzl7RQXGBoa4siRI8bt7UpfjnbgP6j5um58+GE6OjvJ5/M0Nzdv1HIVdcLBI0d4b1cXFuArl3CclJJEIjGvqnN2dpZcLoemaYvcKk1NTbhcLuVW2SDWaonfBLxICPEI8JvAx4UQ29a+LMVyzHWjTE9PG1V1y/EFYBTo+Od/5rfuuovZ2VnlRrmGGZ2Y4L9bLGynVuS1WvSMlVgsRiKRIBqNEg6HlVulDljTJ1lK+ddSyhdLKe8CfgT8Tynl+fVZmmIp5rpRzp07t6If/H3UWpW+FXjt616n3CgKAJLFIq8E/pBaoddqKZfLhm98ObeK2WxWbpUNZt3MMSnlm6SUP1+v11MsZqEbJRJZ1ovFy4A/Bv4bkJOSVCql3CgKg6CUvBL4PLXGZ6uhWq2SSqWMQiC910ruQqVvKBTC4XAot8oGo+6ptwiVSoVvfetbxrDbldwoe6m1Jn0VMHFh9JZyoygW8rSUvAf4d2otiFfDQreKnoKou1VOnTpFf38/J0+eBDDcKhcrQFNcPuoTvUX4yU9+wvj4ONPT0yu2lm2k9qF8L/DUBZGPxWI0NDQoN4piEV+Tkm9Tq+pcbcbKXLdKJpOhWCySz+cZGxujUCgQi8XIZDKEw2HMZjOdnZ1MT0+rJllXCCXiW4BisciPf/xjo9fzctkoJmrTeB4CvnpBwMvlMrFYTFVlKpbl2GtfSx745Cr3n+tW0VvWRiIRotEonZ2dnDlzhm3btjE0NISUktbWVkwmExMTqpTkSqBEfAvwgx/8wCh7DgQCy7pRPgJUge/edpvxXDQaxePxYLVe6twXxbXC177+dV5HLb3sDas8RkpJLBYjlUpRKpUoFovGoGWoDeiuVqtMT08jhKC7u5twOEy5XL5Sp3HNokS8zslkMjz55JPGHMTl3CjPA94M/AHwsyefBGrVdslkEr9/tR5PxbVKQkpeB3wC6FjlMZVKxSjHn56eNtwrLS0tnDt3jp6eHk6fPk21WsXv92Oz2RgfH7+CZ3FtokS8znnooYeYmpoypq0shQ34MrVslNk5VnokEsHn86ke4YpV8Ssp+Sxw/yr3r1arJBIJkskk+XyebDZLIpEwRrqNjIwYo90Aenp6iMfjFAqFK3YO1yJKxOuYaDTKkSNHjOniy43B+jAwCPzTnA9HoVAgk8nQ2Ni4MYtVXBX807ZtdFO7o1sNejVnOBwmHA6TSqUoFAp4PB4mJibo6OjgzJkzVCoVPB4PDodDWePrjBLxOuY73/mOYYXrg2sX8jzgjcC7AMscv3c4HMbv92MymTZkrYqrgzNjY7wJ+Dirc6tUq1WjMVY6nTaqOY8fP47L5WJ4eBifz8fY2BgAvb29pFKpVc19VawOJeJ1ytTUFIODg0YGwFIBoeXcKLlcjkKhgNfr3bD1Kq4enr1Et4q8UEwWi8WIx+NGhabNZiMQCNDT08PIyAjlchmn04nL5VLW+DqiRLxO+f73v29Y4culFH4YOAH864JslXg8js/nU4U9isvmUtwqUkqSySThcJhEIkEsFiMYDDI4OIjL5eLs2bM0NjYSDAaBmjWezWZXbBmhWD3qU16H5PN5Tp06ZVjhlUpl0T4HqLlR3r3g+UqlQiaTUVa4Yk3Mdau0rWJ/KSWZTIZkMkk8HqdcLmOxWHA6nUamytjYGFJK7HY7DodDNcdaJ5SI1yGPP/640ZsiGo0uuc8nqHWhm11ghScSCVwul/KFK9bMs1LyAPC/V7Gv7lKZnZ0llUqRy+UIh8OcP3+eSqVCLBajXC4Tj8cBaGtrIx6PqyrOdUCJeJ0hpeTJJ58kHA4TCoWWtMJ/G2gGhm65ZdGxiUQCn8+3MYtVXPX8LfB7wPZV7KsPkEgkEsbQ5EwmQ3NzMydPnmTbtm1GgLOxsRGTycTs7OwVXP21gRLxOmNoaIjp6WkikciSk8MFtQ/W+4FfPv30vG3ZbBZN07DbL2WSokKxPBEpuRf4q1Xsq1vj6XSaSqWCyWQin8+TTqeJxWI4nU5mZ2cpFosIIfD7/SsONFGsDiXidcajjz5KOBxmZmZmyaHHrwOywOu+/vVF2/SApkKxnvwDcCfwnFXsq2ka8XicmZkZJicnjfGBfr+fEydO0NbWZmSmtLe3UyqVSKfTV3L5Vz1KxOuIRCLB6OgooVCIcDi8aLuFWn+UPwde89rXzttWKpXI5XK43e4NWavi2iEjJR+hdgd4McrlstHZsFQqYTKZsFgsWCwWI91wfHwcKSUWiwWXy2VkrSguDyXidcSjjz5qBDSXKoZ4BzAEPLFEA6x4PI7X61VphYorws927WI78Bur2LdUKhEKhYy88VAoxPDwMCaTiZmZGSwWi+EL7+joMCYEKS4P9YmvE6rVKkeOHCESiSw5sccF/AVLz0WUUpJOp/F4PFd6mYprlBOnT/OXwN+tYt9qtWq0fYjFYiSTSYQQuN1uhoaG2LZtG4FAAAC3243FYlm2IllxcZSI1wlnz541CnuSyeSi7e8Gfkwt7WshxWIRKaVqN6u4ouz5wAcwAy9fxb6pVGpe9abeICsajdLU1EQ4HDYyr3w+37KptIqLo0S8TnjqqaeIRCLMzs4uCmhq1HqjfHyZY9PpNC6XCyHElV6m4hrmQ/fcwyeA/7GKfTVNI5PJEIlEmJqaIhgMMj4+jt1uJxAI4PF4jLhPU1MTpVJpyUC+4uIoEa8TTp8+TSwWW9KV8jJgEjiyzDCITCaD0+m8sgtUKID/9Hq5Edh5kf0qlQrJZJJSqYSUkkKhgMPhwOFwMDIyQltbm+FCcTgcmM1mZY1fJkrE64CZmRlCodCyvZbfBXxmmWPL5TLFYlHNz1RsCDPxOF8C3rmKfavVKpFIhEwmgxCCcrlMpVJhdnaW1tZWZmZmkFIihMDlcikRv0yUiNcBTz31FPF4fMm0wu3Ac6kNsl0K3QpXrhTFRvGPwP8NOFaxb6lUIhqNEo1Gjb74UEunNZvNJBIJoOZSyeVyqgz/MliTiAshPEKIXwohfiKEeFoI8ZL1Wti1xNGjR410rIW8k1q72fwyrpR0Oq1cKYoNZUxKngRee5H99GZsc8vxy+UybrebM2fOzHOpuN1uNE0zRF2xetZqiaeBO6WUL6L2f7qaDCTFHLLZLNPT08RisUWuFDu1ToX/uMyx1WqVXC6nRFyx4XyGxR00l6JYLBKLxSiVSmiahsViQUpJMBikvb3dEHFN02hoaFgyJqRYmTWJuJSyKqXUs/Q9wLGF+wgh3i6EOCSEOKT6JCzmmWeeIR6PE4lEFjW7eg3wNDC6jBWey+Ww2WyqY6Fiw7nxve+lEbjlIvsJIUin04TDYSNvPBKJUCwWMZvNRuoh1JpipVIp5DLXu2Jp1uwTF0J0CSF+DjwKfGfhdinl/VLKA1LKAy0tLWt9u6uOw4cPL+tKeQfw2RWOzeVyOByr8UwqFOvL3330o3yO2jW6EpVKhVKpZLSj1e86rVYrY2NjtLa2GtWbeg98NSzi0liziEspJ6WUtwO3AvetfUnXFoFAwGjZOZdOYDdwtG35lvyFQgGbzXZlF6hQLMM3gFdwcRHJZDLkcjnMZjOapuF2uxFCEAwGaWxsNAwYi8WC2WxesnunYnnWGticqyBJILW25VxbJBIJY8jswt4R/xV4BDhz7tyyx+fzedV2VrFpnJeSSeD5F9nPZDIRDoeJxWLk83kqlYqRfujz+YxgphACu92uuhpeImu1xK8TQvxUCPE48BDwp+uwpmuG06dPk06nl3SlvILaP+hyIl0ul5FSYjabr+wiFYoVeIjatboS5XKZarVq3HEmEgljxqbNZiOTyRhGjNPpVO6US2Stgc3DUso7pZS/IaW8XR3daAEAACAASURBVEr52Hot7Frg1KlTpNPpRa6UBmr9mx+BZfO/dStc5YcrNpPViLhesalpGqlUipmZGSKRCDabjWg0itvtNvoFOZ1OqtWqKsG/BFSxzyYSCARIpVKLUgt/EzgIHB4eXvZY5UpR1ANvue8+3Kxchi+EIJfLUalUsFqtVCoVLBYLALOzs3i9XsOl4nA4MJlMS7ZiViyNEvFNJBQKGaOs5qK7UlYa8KCCmop64A/f/W6+x8qdDaWUaJpGNBqlUqkYvm9N05ienp4n4haLBU3TlF/8ElAivkkkEgmSyeSitrMataDm92DFfijKElfUC6txqWiaRi6XM/oDVatVKpWKEdzU40IquHnpKBHfJPSg5sKL9XnADPDgz362rKWtgpqKeuIXVis3Af4V9tF93IlEwrju9dRDm81GNps1gpsNDQ0quHkJKBHfJJYLat4J/Ce1IonlRLpcLmOxWFRQU1EXJAoFngReeJH95gY3p6eniUajNDQ0LApuOhwOFdy8BJSIbxITExNLBjVvBg4BVqt1WZEul8uq1F5RVxymdu0uhx7c1A2QarWKyWTCZDIRi8Vwu93GXak+WFkFN1eHEvFNIh6PLxnUvJnaB2KloGW5XFauFEVdcTERl1JiMpmIx+OGK9BisVCpVMhms9jtdsOgMZlMaJpGsVjckLVvdZSIbxLZbHaRpdEINAPPe8MbVpyXuZKrRaHYDC4m4lALbmazWSOwaTKZKJVKZLNZbDabMY/TbDYjhFAivkqUiG8CevGDftHq3Az8Cvjrv/1bZYkrthRf//nPsQAdK+xTLpeN2ZvlctmY9qNb4nNFXEqpfOKrRIn4JpDL5SgUCkv6ww9TCwCtZIkrEVfUGy944QtX5VIRQpDNZikWi0YPlUwmM8+dommacqdcAkrEN4FwOEyhUFhkaegiXiwWLyriKrCpqDdWI+KVSsXwg+szNzOZzDx3CtSCm8oSXx1KxDeBUChEoVBYZGnoIl4qlYyy5KVQlriiHjnExTNUisUiJpOJTCZDPp83Bn1bLBYKhYIxEMJsNisRXyVKxDeB2dlZisXivIvUA7QAL3nnO43bzqXQrRkl4op64wjwnIvso49pKxaLpNNp8vk8NpuNQqGA2Ww2DBu9x4qa8nNxlIhvArolPje9sAcIAJ/+7GdXFHEdVeijqDde+sY30gas5OgrlUpUq1VjKHIul8Nuty9KM9Q0zcheUayMEvFNQB+KPNfK6ASCF36/mCWuBFxRj3zhK18hDCw3i0q/3nWXihCCarWK2WxelGYohFAulVWiRHwTmJsrq9MBTF34XQm1YqsyxcpphnPbzFqtViNAr8d55t6dapqmRHwVKBHfBMrl8qKLs4NfW+IrCbgSeEU9E2RlEdc0jXw+T7VaNa5j3SIXQhjWuhBi3t+K5VEivgno+bFz6eTXlri6eBVblSlq1/JKzL32dbHWRXzhNvU5uDhKxDeBarW6SMQXWuLq4lVsRS5miUspjWt/7h2lfoc597pXn4PVoUR8E1iLJa4ubEU9s1pLfKFY6xkr+vNLibpiaZSIbwJSykUX59zA5sLbSoViq3CxwCYw79qf6wNX8Z7LQ4n4JrBUyXwLMHvhd6vVumzfCP1iVxaKoh6ZBVovss9cK1ufv6lb4bqIzxV3xcqsScSFEM8RQvxCCPFTIcSPhRAD67Wwqxk9mDMXK6C3w7qYiGuapix1RV1SAJZvGFFD07R5Yj03E0V/fmGffcXyrNUSnwLuklLeCXwc+PDal3T1oxc6zMUM3Pjc5wK1gRArdXAzmUzGPEKFop4os7KI64I91/891zDRPxd6q1pliV+cNTXgkFJOz/mzSO3/UHER9Kb3OoJaqfJb3/52pJRYrdZFbWoXHq8sFUU9UmJlURFCGJN75gYxdRHXtJpdqV/fSsQvzrr4xIUQTuCvgY8tse3tQohDQohDoVBoPd5uy+NwOOb5xU3Uvv2sThtdH+3i3mfuJZVNLXu82WxWlriiLrmYJT6XUqlEpVKhVCrhcDgoFovYbDaklMb1rYu6YnnW/C8khLAA3wD+Vkp5cuF2KeX9UsoDUsoDLS0ta327q4KmpqZ5rWb1EGWlWsFcNfOlY19iJjnDD87+YMnjlYgr6hUBrBSt0e80gXktaBsaGoyOhpVKxWh+tVJLZkWNtQY2NeBrwHellN9dnyVd/bS0tGC3242/K9RuQX/0o0dpNDcym5/FbXbz+9/8fUaiI4uOVz5xRb1iYWWfarVaxW63U6lUsNvtNDQ0kM1m51ni+tATJeKrY62W+CuBlwFvEEL8RAjxqXVY01VPZ2fnPBGH2oX/00cfp9HcSJUqqXIKp8nJvU/du+h4ZYkr6hUzNb/4clQqFSP7yul00tDQANQsdJPJhMlkMixx1Td/daxJxKWU35JSuqSUL7rw+KP1WtjVTFtbG3a7fV7QpgxEZ0M0mhsBCBfCNFoaefDYg4uOV4FNRb1iZnlLXM82MZvNWK1WbDabYY0Xi0XDsNEHKi9MAFAsjYoabAKNjY3Y7fZ5t4o5wAn02foAGMuMMeAcIF1MLzpeWeKKeqWB2rW8FEIIbDYbAG63G5PJRENDA06n0/CHA8a1rVwpq0OJ+CZgt9ux2W3zbhWngHZgwFGrlzqTPsMu9y5cVtei43URVwU/inpjbvuIpbBarZTLZZxOpyHqDoeDQqFgWOJ6jYQS8dWhRHyTcDgc8yba642D3MKNQ3NwOnWa3e7d/MENf7DoWE3TLppLrlBsBnMbuS1En+Jjt9uxWq2YzWbMZvO8zBSAfD6PpmlKxFeJihpsEk6nc15wU2/hGSvG6Lf1M5wepq+hj5fe+NIlj9eHyzocjo1ZsEKxCua2VF4KKSUuV+3uUs8B190pLpeLSqVi3GUqEV8dyhLfJLwe7zwR1y3xfCnPDscONJOG1Walzbb0xEK73W7MI1Qo6oWVLHGTyUS1WjXEWhdst9ttpBkWCgUjzVCJ+OpQIr5J9Pf3G+lVMKeZfh5e0vESjr3zGB1NHaRSS1duKhFX1CPLWeJ6Vor+KBaLaJpGKpWiqamJRCKB1+s1RHyue0WxMkrEN4mbbrqJxsZG42/dEi8ny/Tb+tnu347b7SadXpydArUAUalUUsFNRV2xnCWuT/Tx+/0UCgXcbjcdHR3zLG7dMLHZbGSz2XlGjmJ5lIhvEnv37qWpqcnwC+rN9MfHx8lmsxSLRVwu17KWuApuKuqRlbJT9GpNp9OJ1+vF6/XS3NxMMpnE6/UihCCfz8+z2hUXR4n4JuFwOPD7/YZfPEjNihkZGcFqtRKNRnG5XGQymWULe/TgpkJRDzQ1NOAAokts03PCpZT4/X6sVitWq5WWlhYSiQQ+n8/wkZdKJZxO50Yvf8uiRHwT6ejowO12AxCglif+5OOPU6lUmJqawmw243a7icfjSx6v/OKKeqI7l2N0mW16KX1LS4sxAKJardLc3Ew8Hp/nD8/lcsqVcgkoEd9EduzYgd/vB2rN2IeAG4HJyUkmJyeBWsfDcDi85PF68yA1qk1RD9wMHF5mm54X7nQ6jQZXqVSK5uZmI6ipZ6gof/iloUR8E7npppsMEYfaB+Bman7xQqFAMpmkubmZSCSypFBbrVaEEMqloqgLDrC8iEspaW5uplAo4PF4lgxqptNpnE6nEvFLRIn4JrIwuKmL+OjoKBaLhcnJSRoaGjCZTMtmqbhcrmW3KRQbyXKWuBDCCGo2NDTg9/vxeDzzgpr6gAiolduroObqUSK+iTgcDqMZFvxaxI8fP044HDZcKs3Nzcu6VJxOJ5lMZoNWrFAszXOvv57rgF8tsU3TNJxOJ5VKBZ/Ph81mw2Kx0NLSQjwex+fzkclkcDqdyh9+GSgR32QGBgbweDwAHAN2ApMjI4yOjhKLxSgUCjQ1NRGJRJY83uFwUCqVKJVW6uKsUFxZKidOMA4sZU5UKhU0TaO5uRlN08jn8xQKBbq7u5mdnaW5uZl0Om1kYykRvzSUiG8yL33pS+ns7ATmBzfPnTuHw+EgGAzi8XjI5/NLZqIIIZQ1rth0VgpqWiwWvF4vFosFh8NBR0cH2WwWn89HLpfD4/EYfYD0TBXF6lEivsnceuutdHd3G3/rLpWTJ0+SyWSYnJxE07QVrXHlF1dsNhfzh7vdbqNKs7m5me7ubkKhEK2treRyOaNviqZpi6ZeKVZGifgmYzab2blzp2F9HAJuBc6cOcPw8DBTU1OUSiWam5uZnZ1d8jX0Vp6qBF+xWdxK7dpdCpfLhcvlwuv1Gtfotm3bmJmZoa2tzfCH61a4muZzaSgRrwNe+MIXGi6VHwF3AflslvHxcRoaGjh//jxNTU3kcrkl3SYmk8lI0VIoNpp2IdgBHFxim8lkwmaz4fV60TTNKF7r6OggGo3S1NREJpPB5XIZlZuKS0OJeB3w0pe+lN7eXgDOAbPULJvBwUGi0Shnz55FCEFHR4eRsbIQr9dLIpHYsDUrFDovA37I0gOSK5UKra2tmM1menp66OrqorW1lWQyid/vJ5/P43A4qFarFAoFo9e4YvUoEa8D3G43fX19xoCH7wEvB06dOsXQ0BD5fJ5wOExHRwezs7NLztd0uVwUi0VV+KPYcF5O7ZpdiKZpeDwezGYzbW1tWCwWLBbLPFeK7kLRqzaVK+XSUSJeJ9x88810dHQA8BDwCiAajTI8PIwQgrNnz2K32/H5fMzMzCw6XgiB1+tdts+KQnElcAjBbwA/WGa7nloItarNZDJJb28vs7Oz+Hw+qtXqPH+44tJRIl4n3HXXXUaWytNAMzBArfBnaGiIyclJ8vk8XV1dBIPBJcvwvV4vqVRKBTgVG8ZLgCNAbIltmqbR1taGyWSio6ODbdu24XA4KJfLRvM2PdiZTqeViF8maxZxIcQPhRAhIcRfrseCrlV6enoYGBjAYrEggYep3aaePXvWCHCOjIwY1stS/m+LxUJDQwPJZHKjl6+4RnkFtTvHhej97nVXit1ux263s3PnTgKBAO3t7YZwJ5NJnE4nJpNpo5d/VbAelvhbgPeuw+tc8/zO7/yOEeD8HrUPSLlc5sSJEwSDQYaHh5FSGtb4Uvh8PuLxuOpsqLji/Pn73resP1wIgd1up62tDZvNhsPhIJlMMjAwwMzMDF6v1xDuUChEU1PTRi//qmHNIi6lnFiPhSjg7rvvZt++fQD8J3AL4KPmUtF948FgkLa2NqLR6JJBTIfDgZRS9RlXXHF+/LGPEQNGltne0tKCx+Oht7eXbdu20dfXZxT46BWb+XyeXC43b1Sh4tK44j5xIcTbhRCHhBCHQqHQlX67LY3VauXWW2+lqamJLPAfwOuBeDzO6dOnyefzDA4OYjabaW9v5/z584teQwhhWOMKxZXkTcA3lnhebwXR3NxMY2OjYVTs3buX8+fP09bWZlRmhkIho6eK4vK44v9yUsr7pZQHpJQHWlparvTbbXne+MY3snfvXgA+A7zrwvPHjx/n5MmTFAoFAoGAEeHPZrOLXsPj8ZDNZlW6oeKK4RGC1wKfX2KbyWSiubkZl8tFW1sbAwMDOBwOY3ZmsViksbGRarVKJBIxslcUl4f6+qszenp6uO6667DZbPwUkMB/odYQ6/Tp06RSKY4dO4bZbKa7u5uxsbFFr2EymfD7/cu2r1Uo1sobgMdYPBRZ94V7PB66u7ux2+3YbDb279/P+Pg47e3tCCFwu93EYjFcLhc2m20TzuDqYT2yUz5PLbD5JiHEd9e+JMWrXvUqtm/fDvzaGpdScuzYMU6dOoWUktHRUbq7u4nH46RSqUWvoc8szOVyG7t4xVXPN//1X3kXtWtzISaTCYfDQXNzMzabjaamJtLpNJ2dnUSjUcNKB5idnUXdna+d9Qhsvk1KuV9KuUNK+bvrsahrnRe/+MVGgPNB4DeBDmrphqdOnSIcDnP8+HGg1khodHTxeFq982EoFFKZKop15VOveQ0m4CcLnterLZuammhtbaWrq4v29nZ27dpFMBikpaUFq9VqzIatVCpGL33F5aPcKXWIpmnceeeddHR0kAL+BXgrUCqVOHbsGGfOnMFqtXLmzBk6OjrI5XLEYovLLTweD9VqVfUaV6wry1nhUOuo2d7ebrSTTafT7N69m/HxccNCF0IYVrgqs187SsTrlDe/+c3ccMMNAHwWeDtgBkZGRhgcHGRiYoJTp05RKpXo7+9ndHR0kcUthDBGuylrXLEetAvBbwEPLLFN0zRaW1vx+Xx0d3ezfft2urq6iMfjRitah8NBsVgkkUiogOY6oUS8TnG5XLziFa+gvb2d49S6G76Cmm/86NGjDA8PY7fbOXnyJK2trVSrVZZK4XQ6nWiatqTfXKG4VN4G/CuwsCZYCEFDQwNdXV00NTXhcDjIZrPcdNNNRqWxXtAzNTVFc3OzGoa8TigRr2Pe+ta3cssttwBwL/CXgADGx8c5fvw4Z86cYWRkhHg8zs6dOxkeHl40a1MIQUtLC+FwWPVUUawJrxD8D+AfltimaRrt7e20tbUZLSR27NhBKBTC6/Xi9/ux2Wzk83ni8Tjt7e0bvfyrFiXidYzVauV1r3sdvb29fIdav+bfv7Dt6NGjjI6OomkaTz31FB6Ph5aWFoaHhxe9jsPhwOFwLDveTaFYDf8vtW6FpxY8rw976O7uxuv14na7yWazRlqh1+s1XCeTk5O0tbUpK3wdUSJe57zmNa/htttuA+DPgb8CLMD09DTHjh3j2LFjlMtlBgcHGRgYIJFILJkf3tLSQjKZVCmHisuiTQjeDXxwiW26Fd7R0UFXVxfd3d3s27ePiYkJvF4vnZ2dmM1mMpkMmUyG1tbWjV7+VY0S8TpH0zTe9KY3sWPHDh6n1qfirRe2HT16lFOnTjE6Osrp06dJJpPs2bOHM2fOLHKrmM1mWltbmZ6eVm4VxSXzAeCrwPiC5zVNw+v10tXVZbhNcrkcO3bsIBgM0tHRgdvtBmpWeEdHhyqxX2fUv+YW4K677uIFL3gBQgj+FzXfeAOQSCR49tlnGRkZWZVbxe12Y7PZlFtFcUkMXCix/5sFzwshjF7hnZ2dtLW10dbWxk033cTw8DCNjY10dXUhhCCZTFIsFlVGyhVAifgWQAjBO9/5Tq677jp+BfwU+NML24aHhzl69CjHjh2jUqlc1K2izzdUbhXFavkItWDmwq9+IQR+v5/29nZ8Ph/t7e0Ui0W6urqIRqPs2rULs9mMlJLJyUk6OztVXvgVQIn4FuH5z38+d9xxB1arlb8E/gzwA9Vq1XCrjIyMKLeKYl25SQheTC07aiFms5nOzk46Ojrw+/00NDRw6623cuLECTo6OozJ9Xohmmo3e2VQIr6F+PCHP8ztt9/OCPBNan5KmO9WEULw1FNP4XK5aGlp4fTp04sKfZRbRbEaHnzgAf4/4K+BpWp+u7q66Onpobu7m507d9LS0kKlUqFarbJv3z6EEJTLZSYmJuju7lZW+BVCifgWorm5mXe84x309fXxv4HXALdd2Ka7VY4fP065XObgwYMMDAxQLBaX7Dve1tZGMplUJfmKZXnsjW+kFfjHBc/rbpTe3l46OjpobW2lXC5z0003MTY2xv79+7FYLAAEAgF8Pp8R3FSsP0rEtxi/93u/x2//9m+TsFj4I+ArgJ1fu1UGBwc5ceIE4XCYoaEh9u/fz9TU1KJqTpPJRGdnJ9PT0xSLxU04E0U90yUEH6U2+KG8YJvVamXbtm10d3fT1tZGU1MTt912GydPnjSeg9owk3Q6TVdX1wav/tpCifgWQwjBPffcwx133MG3gaPAPRe2JRIJfvWrXzE0NMS5c+cYGhoiFAqxf/9+zpw5QzqdnvdaDoeDpqYmJicnqVQqG30qijrlwQce4PPAp4FjS2zv7Oykp6eH3t5edu3aRUtLC7lcDk3T2L9/P1CbDTs+Pk5fX58agHyFUSK+BZnrVnk3tQb9ultlZGSEZ599ltOnT5NIJDh48CDVapUdO3Zw4sSJRVa3z+fD6XQyNTWlmmQpgJobpYOlUwqXcqPs2bOHqakpbr75ZiMHXLlRNg4l4luU5dwqACdPnjQyViqVCj/96U/x+Xy0trYyODi4KCulpaUFKaWaBKRYlRulp6fHcKMcOHCA4eFhdu7caQi2cqNsLErEtygruVXK5TJHjx7l6NGjnDhxAoCf//znbNu2DbPZvKgQSAhBZ2cn6XSaZHJhfzrFtcJKbhT9Gunp6aGnp4ddu3bR3NxMOp3Gbrcbk6iUG2XjUSK+hWlubuaP//iP2bFjh+FWef6FbdlsliNHjnD8+HGOHj1q/L1nzx7i8TiBQGDea+mBzlAopAqBrlEeX8aNArVpPdu2bTOm9VSrVTo7O4nFYjz3uc9FCIGUUrlRNgHVSmyLc/fddzM0NMRHP/pR3hGN8g3gVmAaiEQiHDp0CIfDgc1mo1wuYzabuf766zl69Ciaps275bXZbLS1tREMBunp6cFqtW7WaSk2mOcKwQ+BF7PYjeJyudi5cyd9fX10d3fjcrm4/vrrmZ6e5jnPeQ4NDQ0AhEIhstkse/bs2ejlX9MoEb8KeN/73sf4+Dhf+MIXuL9Y5DvAi4ACMDExwa9+9StsNptRbGEymbjxxht59tln0TSNjo4O47VcLheVSoWJiQl6enqMfF/F1UubEDwNvBM4sWCbzWZj165d9Pf309vbS2dnJ3v37iUSibB9+3YjnTCZTDI1NcWePXuUG2WDUSJ+FSCE4GMf+xgzMzP81be/zfXA54A3X9h+5swZbDYbmqahaRpSykVCrn8YAbxeL9Vq1ai0U0J+9WITgh8DXwb+bcE2s9lMf38/27dvN352dnZSKpVobGw0/OCFQoFz587R39+PzWbb6FO45lEifpXQ0NDAJz/5ScLhMG9+4gl+Tq2/yr3UCoFOnDhhiLjuv5RScv3113P8+HGklPOmrTQ2NiKlVEJ+FfPQv/87nwFm+HVQXEcIQXd3Nzt27KC/v58dO3bgdrvxer3k83luvPFGhBBUKhWGh4fp6OhQk+s3iTWLuBDiTdTm+Ergj6SUR9b6morLo7e3l3vuuYd3v/vd3H3iBE8Bg8CjQKVSMYRcF3H9Q3jDDTdw/PhxI1il4/f7AZSQX6X85+/+Lm8FXkDtw6sjhKCtrY0dO3YwMDBAf38/ZrOZvr4+wuEwd9xxByaTCSkl586dM/r0KDaHNYm4EKIR+GNqtSZdwIPA7euwLsVlcuedd/Knf/qnfPCDH+T3Jyf5NnAHcBYolUocP34cIQTVatUo7qlUKtx4440cPXqUarVKd3e38Xp+vx8hBIFAgO7ubhXs3IKMREf4xJOf4GvHvka6mMZldfG896d4gFo208LuOV6vlx07dtDX12cEuPfv308oFOKWW27Bbq9VJASDQSqVCgMDA6q51SayVkv8ecDPpJRF4JwQwiWEsEkpC+uwNsVl8pa3vIWRkRHuu+8+/iKV4iHghUAUKBaLHDt2zBBxvfBHv0U+ceIEuVyO7du3G9V3jY2NCCGYmJigq6tL+T23ED84+wNe/c1XU6qUKFVrbYk735/ia9TmtS5sjeZwONi5cyf9/f309PTg8/m47rrriMVi7Nu3z2gnG41GiUaj7NmzR03q2WTWKuJ+IDbn78SF56b0J4QQb6fmbqG3t3eNb6dYLR/5yEfIZrN84QtfYHsmw6PAS6j9BxWLRcN9Uq1WjZ7jqVSK5z//+Zw7d47jx4+zb98+w4Xi8/nQNI2JiQna2tpwuVybdm6K1TESHeHV33w12VLWeG7gQ/Aj4H3UhovMxW63s2/fPnbt2sXAwADt7e3s37+fbDZLd3e3cYeWSqUIBALs3LlTudjqgLWKeBTwzfnbe+E5Aynl/cD9AAcOHFDNOTYIk8nEJz7xCcrlMu//9KdxAI8AvwmkqblWTpw4QaVSoVgsUi6X2b17N4899hh33HEH8XicI0eOcN111+F0OgHweDxYLBaCwSDFYtGw0BX1ySee/ASlyq+HgvR+CB6jNmz7gQX72mw29u3bx+7duxkYGKCzs5MbbriBYrFIe3s7O3fuBCCTyTA6Okp/f7+RH67YXNZ6H3QQuF0IYRFC9AJp5UqpH0wmE3//93/P2972Nv4MeBb4PrX5nFArkR4cHOSZZ54xRrzF43Eef/xxTCYT27Zt49lnn53XU8XhcNDb20sqlVLTgeqcrx37GqVqCSEFnR+qCfi9XLCo5mC2mA0B7+/vp7u7m5tuuolyuUxzczO7d+9GCEE2m2V4eJht27apTJQ6Yk2WuJQyJoT4DPAEtQD3n6zLqhTrhtls5lOf+hRSSt71hS/weWoW+cuAFLX0w7Nnz5LJZMhkMuRyOQYGBjh8+DDbt29n//79nDx5kkwmQ29vL0IILBYLPT09zMzMEAgE6OzsVLfVdUi6mMZZddJ0T4bHqA13+D8L9rFYLOzcs9MQ8N7eXp7znOdQKBTwer3s3bsXIQS5XI6zZ8/S29trjF1T1AdrTjGUUn4J+NI6rEVxhbDZbHzmM5/BbDbzts99jvuo+UXvAuIX9gkGg+TyObK5LJlMhmw2i6ZpJBIJbr75Zs6cOUMmkzGG32qaRnt7O7FYjEAgQEdHBw6HY/NOUrGIfXIfuXsGeQz4JPCpBdstVgt79u5h1+5dbN++nd7eXm655RbS6TR+v98YsZbNZjl79iw9PT1qTmYdoop9rhEsFgv33XcfFouFd3/qU3wS+DHwfwG6syQWjXH40GEy2ZpFXiqV2LVrF4899hgHDhwgnU5z6NAhdu/ebfjD/X4/VquVYDBIU1MTXq9X+cnrgAcffJDKhwf5CbVp9Z9fsN1qs7Jn7x5279nN7h276enp4bbbbiMWi9Ha2sqePXsQQpDJZBgeHqa3t1cJeJ2iRPwawmQyce+992K2mHnPJ+/lw9SCGr8LHL+wTy6X49jRY2SzWfKFvJGZ8MtfxRHqxQAAE/hJREFU/pJt27bR19fH0NAQTU1NDAwMYDabcblcWCwWpqenSafTtLW1KffKJhGLxXj/+9/P6Oc+xxPA/6RWvDGXBmcDu/fuZs/uPewY2EFXdxe33XYbkUiErq4udu7ciRCCdDrNyMgIfX19eL3eTTgbxWpQIn6NYTKZyP9GntebX88D336KoZERHqOWA/rdC/tUyhVOnzpNPpsnmUgSiUTo7e1FSsnU1BTPfe5zyWaz86xym81Gb28vsViM8fFxZZVvAg899BB33303f0atF8orgV8s2Mff5Gfv3r3s3r2b3r5e2trauP0FtzM5OcmuXbvo6+sDMP4f+/v7VRCzzlEifg3yteNfI9WQ4sX//cWMPtrMbx88yHeA66mln+l5oOfPnyeRSBCJREgkEmQyGbq7uykWi8ta5X6/H6fTqazyDSQWi/GBD3yAL3z603wZuJFaCXVgwX6d3Z3s21vLQunc1kl3Rzc333QzwWCQm2++maamJuOLOhKJsHPnzv+/vXuPrfOu7zj+/p37ec7dzzk+9jnHPo5vcSSS9ZIRWKHaEBODQdEmpGlSO6AasE3bgO0fpk1qp45usAEDVoKq0oJaolSdmjZtEgalrCvNPUYhtJWT2Ln6fj/3+29/PD5PbOfSQmOf2P69pKMT24/t3+PYn/N9fs/vooYRrgEqxDegTMnYMPll68v0/EEPmyMf5d0vvsizGEH+KaA+PWRubo6TJ08yMzPD7OzsDavy3t5empqaVFW+ivbu3cvXv/51Bl55hZ9hBPf7uPL/B8YIpa6uLvr6+uju7qanp4f29naam5uZmprirrvuQtM0qtUq58+fp1wu09fXp1581wgV4huQ1+ElXUoDcEacYfzOce7R7+X3fvAU38W4BP84cHHh+Eqlwrlz58wQr1fl8XicfD5Pa2srHR0dnD59GrfbTWdnJz6f76qqPBKJqCn7N8nw8DAPP/wwzz33HLGREY5iDCH88rLjXC4XfX199Pb20t3dTTKZZOvWrVSrVaSU3HXXXdhsNorFIoODg3g8HjZt2qSm0q8hKsQ3oHu33ctj/Y+Za2mkRIpdm3Zx399+kk9/6wd8ATgMfAZjclDd8qp8amqKRCJBJpNhZGSEjo4OvF4vp06dIhAImLP62tvbmZub4/Lly3g8HnRdV1Xeb2h2dpZHH32UAwcO8Morr/A5jNEni+9p1AWDQfr6+ti8eTM9PT1Eo1F27NjBxMQEiUSC3t5ehBCk02mGhoZobW0lEomoK6Y1RtRXslsN27dvl8ePH1+176dc2+DMINu+u23Jmhp1Hy9+nEv/cwlffz+PA68CX+DKePK6eki3tbURj8fNTXSdTic9PT0EAgHGx8dpbm4mmUzidDqpVqvMzs4yNzeH3+9H13W1C8zblMvleOKJJ3jppZc4dOgQzvFxvgf4Mbq/3lx0rBCCeDzO5s2bzQo8HA5zxx13MDIywrve9S5zyeHJyUlGRkbUDcxbnBDihJRy+zU/pkJ8Y7rW6nYAdoudPksfW9/cyvM/3MVXMLpWPgfsX/Y1rFYroVCIZDJpbt3V2tpKMpnEbrezefNm3G43k5OTZsjb7XYqlQrT09Ok02lCoRChUEhdvl9HqVTi6aefZt++fRw7doyhoSE+i3ED+mvAfwDVRcdrmkZvb6+5EmEymaSnpwePx0O5XOb222/H7/dTqVS4dOkSuVyO7u5u1c11i1MhrlzT4Mwg3zj8DZ785ZPmOtP3bbuPL77ni2wKbmLnzp08/vjj+Pv7+R7Gqndf5Oqq3GazEYlEzPWnW1tbaW5upqOjA7vdTnd3N1arlVQqRTQaJRaL4fF4KJVKTE9Pk8vlaGpqIhAIqDBfUCwW2bt3Ly+++CLHjh1jYGCARK1mVt+fBt5YdHy9+q7fvOzq6kLXdbZu3crMzIwZ5vVZuBcuXCAYDBKPx9XV0BqgQlz5jQ0ODvLggw+y56mnbliVgzErtKWlhWQyaYZ5KBQiFovhcDhIJBIEg0Gy2Swej4dYLEY4HDbDPJ/P4/f7CQaDG3bzieHhYZ5//nn6+/s5ceIEr7/+OuVy+W1V3/Vx3sur79tuu41AIGBW35lMho6ODnw+X0POUfn1qRBX3pFarXZVVX4c+Efg7DWOdzqdxGIxEokE0WiUcDhMNBolFAqh6zo+n494PA4Y3QWtra3EYjGzSpyfn8fhcBAMBvF6vev+Rlu1WuXQoUP85Cc/4fTp05w5c4bXX3+dQqHADuCrgIsbV989PT10dnai6zrbtm1jenpaVd/riApx5aYYHBzkgQceYM8Pf8jngb8D/htjk93RaxzvcDjQdZ3W1lZaWlqIRqM0Nzfj9/uJx+PY7XaSySRut5tsNksoFCIajZqb8c7Pz1MulwkEAgQCAWy29TWYampqiv3793Ps2DEGBwc5e/YsFy5coFQq0Qc8DGwHHgR+wNLq2+fz0d3dbW5knEwm6e7uxuv1qup7HVIhrtw0tVqNXbt2sXv3bg7t28eXgPsx1qj+CsbOQctZrVZ8Ph/RaNQMdF3XaW5uJhQK4fP5cDqdRKNRHA4H1WqVYDBIOBzG5/ORz+dJp9O4XC68Xi8ej2dNDlGUUjIyMsLRo0c5deqUGdyDg4OMj48DkAD+GfgoRgX+CFBY9DU0TaOzs5Ouri5z1mz9/kM9qOtrn8zNzXHp0iVVfa8DKsSVm65QKPDYY4+xZ88eTr/8Mg8A9wD/DvwXS4NnMbfbTSQSoaWlhebmZnRdx+/3EwqF8Hq9hMNhrFYr0WgUj8dDtVrF4/GYE4eq1Sq5XA673W4GutPpvGW7XCqVCgMDAxw8cpAL5y4wMT7B2NgYQ0NDnD9/nmzW2Ka4CfgHjC6T72L8HBe/IDqdTjo6OpZU3rqu093dTS6Xo7W1lZ6eHtxuN6lUiuHhYaSUtLW1qep7HVAhrqyY+fl5du7cyQsvvMDMwYN8GfhtjLWrHwemr/N5NpsNn8+Hruvour6k8g6Hw2bY2+12s1qXUmKz2QiFQrjdbqxWK7VaDYvFgsfjwePxmO9vlFKpxOjoKAMDAwwMDHDx4kXGx8cZHx9nbGyMsbExpqenqVaNzpEkxo3izwDPYHRNjS36ejabjfb2drq7u+ns7DTDu7e3l2KxSFNTE319fXi9XnK5HMPDwxSLRWKxmNo+bx1RIa6suLGxMb797W9z4MABxC9+wV9jLHH7AvAdjCVvr8dqtaJpGk1NTYRCIcLhMKFQiEAgYAa20+kkEAjgdDqx2+0EAgFcLhcWiwWn02kuhyuEwO12mw+Xy4XT6VyRYK8H9uDgIJcuXWJycpK5uTlSqRSTk5OMjo4yNjbGxMQE8/Pz5lZ2AvgQ8FfAezH2u3wEGFr2M4nH42Z4t7e3E4lE6O3tpVar4fF42LJlC8FgkEKhwMjICJlMhtbWVnRdV0M11xkV4sqqGRoaYufOnRw+fJg3fv5zPgX8JZDCCPNdQP4Gn2+xWHC73fj9fjPUNU3D5/Ph8XjMf9eD3e/343K5sFqtOJ1OHA4HFosFm82Gw+HAbrdjs9lwOp1omoamaTidTlwul/kxq9WKEAIhBFJKarUa5XKZfD5PNptlbm6O+fl5ZuZnmE3NUsgVyKfzpFIp0uk0qVSKTCbD3Nwc6XSaiYkJxsfHyeWWzojVMbpL/gJjrP0jwO5lPw9N04jH4+Zs2EQiQTgcNifkWK1WtmzZQjgcplKpMDIywuzsrHnTWPV7r08qxJVVNzY2xu7du3n11VfZ8+yz/D5G5fk+jE0KdmEMU7zRb58Qwgx1l8tFIBDA7XYTDAaXBLvT6TQDvB7M9eCuv68e1BaLBWERCITZ1VCr1ZBSIqWkUqlQKpUol8uUSiUKhQLlcplyuUyxWFwS7Ol0mtnZWVLpFJl0hnw+b3aT1NmAu4E/w7hn8DzGi9mxZeeq67o5xjsej5NIJNB1nUQiQa1WIxQK0dHRQTQaJZfLmVV/OBympaVl3Y3cUZZSIa40TLlc5kc/+hH79+/n6NGjTPX38xmMDQtCGN0tL2DsxH6jCr3OYrGYXSYulwu/34/dbjera5fLZVbXmqaZb9e7FyQSi7BcCXBZA3nluR7WhUKBYrFIuVymUCiQL+QpFopkc1nS6TSFfOGqwK4LYuxfes/C82mM/u4ngJlFx9ntduLxOMmOJLFEjI42Y6RJW1ubOTywra2NZDKJpmnMzs4yMTFBpVIhEomohcQ2EBXiyi1hYGCAZ555hoMHD3Lq1Cmcly/zMYywuxP4GbAXeBGY+DW/dj2khRBm5W21WnE4HDidTiwWCxJpdpsIIajJGrImlwR4qVSiUq5QrVYpl8tmlf5WNi2cxz0YY7v/d+Fc9rH0RqXVZjWH/LW1tRGLxQi2BAk3hdm2eRtCCDRNMyvySqXC5OQk09PTaJpGJBJRa7NvQCrElVtKJpPhwIEDHD16lNOnTzMwMMDEwIBZvX4IuIDR3XJi4fFLoNi4Ji/hBW7HeOG5E2M0TgDjimIvV19VOJ1OwuEw8UTcHFoZDAfRQzolrYTm0+gKdbGpzRg66Pf7SafTTE5Oksvl0HVdrcW+wakQV25Z8/PzvPbaaxw7dswM9F+eOMFtGAG5feG5FxjgSqifAM5hDGFcqd9gK9AM9HAlsO/EmJBzallbfrWsHV6vl2g0SiKRIBKN0Bxpxt9krAtT1IroPh2f28fx1HFOZk7y6CceRbfpZl+7pmnmKB010kRZsRAXQvwRxuzgTVJK11sdr0JcuZFSqUR/fz8HDx7kzTff5OLFi4yOjjI6Okp6aoptXAnSO4B2jBX9xoERjKn/i5/HMar3ysKjjDG8zwbYF57dQOuiR2zRcxjjReIcSwP7DZZOgbfZbHi9XpqamohEIuZzc3OzMfY96GPOOkfMHyNHjiPzRzgyfwSscHfL3dy/5X5cGP3763WJAeWdWckQ14Es8CspZfdbHa9CXHm7arUaQ0NDnDx5kvPnzy+ZMFMP9rm5OSyVCi1cCd/FQRwFHFwJbTtQY2moFzFCf/kLwCjGi8DyW5eLAzscDhNqChnL6AYDBINBAj5jbLvf7zf748PhMBavhQOXD9A/2k+vr5e7o3fjc/rQQzrJaBKv16sqbuW6Vrw7RQhxVoW4spKklMzMzHD27FmGhoYYHh5mYmKCyclJstms+Uin02QyGTKZDLlcjkKhQKVSATAn2yz+na8POwTMseb18eSapuH1enG73WgezZhA5HLj9Xvx+r0EvUGEJrBpNkKeEGFvmOZIM8Fg0FzfpVqtks/nyeVy5ogZj8eD3+/H7XarG5TK23KjEF/xazYhxGcxtgCkvb19pb+dsk4JIcwp+jt27ACMUJ6amuLixYtMTEyQSqXIZrNmaNbHdOdyOarVKtValWq1aoS5xBwvbrEaY8etFisulwu3y20OWawPU7RYLWAFl8NFnjxZmSVLlpzIkSfP9r7tdDZ1UiwWze9XfyGIRqNomqaGAyor4i0rcSGEFWMD9OX2SSkfWjhGVeLKLUFKSblcNivzXC5HLpcjm82SyqaolCsIBMVqkVK1ZG5NJ4Wkyd2E5tDMoYn1SUP1cebFSpFMMYPT4iRbyTJdnMblcNGhd6B7dXO8ugps5WZ7R5W4lLIKvOemt0pRVoAQAofDgcPhIBQKXfXxwZlBvnX4W/z43I9xW9zEPDE+0PEBPtz1YSKeiDlzc/F48vrDbreb4a76r5VbxTvqThFCvB94AIgJIV4CviOlfPamtExRVkBXUxff/Mg3G90MRblp3lGISylfBT54k9qiKIqi/JrUNaGiKMoapkJcURRlDVMhriiKsoapEFcURVnDVIgriqKsYSrEFUVR1jAV4oqiKGvYqq4nLoSYxFjv/1YWBqYa3YgG2Ijnrc55Y1gP55yUUkau9YFVDfG1QAhx/HprFKxnG/G81TlvDOv9nFV3iqIoyhqmQlxRFGUNUyF+tUcb3YAG2Yjnrc55Y1jX56z6xBVFUdYwVYkriqKsYSrEFUVR1jAV4tchhPhPIcThhceXGt2elSKE+JQQ4qAQ4jUhxB2Nbs9KE0LcvnCu/yeEeFkI0dnoNq0WIUSvEKIshHhfo9uyWoQQdwohfiyE+JkQ4quNbs9KUH3i1yGE6JFSnhFCWDD2GL1XSjnY6HbdTEKIEPBTjO334sCTUsp1/QcuhGgBslLKtBDiI8CfSinva3S7VoMQ4kmgFXhQSvnzRrdnpQkhHMA+4I+llOlGt2elqEr8OqSUZxaea0B14bHe7ABelVKWpJTnAK8QwtnoRq0kKeXYoj/oElBpZHtWixDi3cAYcLnRbVlF7wUywK6Fq673N7pBK+Edbc+2EQgh7gMGpZTnG92WFdAEzC56e37hfaONac7qEUJ4gC8Dn250W1bJP2Gc69ca3ZBVFAN+C7gN8AE/FUJskeus+2FDh7gQworRVbLcPinlQ0KIDwKfBD62ui1bNTNAcNHbgYX3rWtCCDvwNPCvUso3Gt2elSaE+EPguJRyWgjR6OasphngoJQyBaSEEFNABJhobLNurg0d4lLKKkZ/8FWEEDuAh4APSynzq9qw1XME+JeFUGsFMlLKYoPbtKIW7nE8BTwnpXyu0e1ZJbcBvyuE+B1gK9AnhPgTKeWtvhjdO3UEeEgIYQPcQDMw3dgm3XzqxuZ1CCF+tfDP+upnfy+lPNGo9qwUIcT9wJ8DEvi8lPJ4g5u0ooQQnwC+D9TP85SU8m8a16LVJYT4PvDYRrixCWZ36OcAO/BvUso9DW7STadCXFEUZQ1To1MURVHWMBXiiqIoa5gKcUVRlDVMhbiiKMoapkJcURRlDVMhriiKsoapEFcURVnD/h90iWIq8b0EuwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameter:\n",
      "z1:  1.9122505537257344\n",
      "z2:  1.9721518658963153\n",
      "R:  2.9282300052564363\n",
      "Anzahl der Interationen: 223\n"
     ]
    }
   ],
   "source": [
    "########################################################################\n",
    "# Startparameter (z1, z2, R)\n",
    "p0 = [2.3794, 2.2532, 3.0370]\n",
    "\n",
    "# Schrittweite\n",
    "gamma = 0.2\n",
    "\n",
    "# Messpunkte: x Werte gegeben, y Werte werden mit get Points berechnet und mit einer Gaußschen Normalverteilung verrauscht\n",
    "#(0 = kein Rauschen)\n",
    "# getPoionts(Messpunkte, Standardabweichung, Radius, z1, z2)\n",
    "x = getPoints([[-1, 0], [0, 0], [2, 0], [2.5, 0], [3, 0], [4, 0]], 0.1, 3, 2, 2)\n",
    "    \n",
    "# Gradientenabstiegsverfahren anwenden\n",
    "#festes Gamma:\n",
    "res = Gradientenabstieg_festesGamma(nabla_function, p0, x, gamma)\n",
    "#dynamisches Gamma:\n",
    "#res = Gradientenabstieg(nabla_function, p0, x, gamma)\n",
    "\n",
    "########################################################################\n",
    "\n",
    "# Messpunkte anzeigen\n",
    "for i in range(np.shape(x)[0]):\n",
    "    plt.scatter(x[i][0], x[i][1], s=100, facecolors='g', edgecolors='g')\n",
    "\n",
    "# den letzten Kreis in rot darstellen\n",
    "circle = plt.Circle((res[0][0], res[0][1]), radius=res[0][2], fc='none', ec='r')\n",
    "plt.gca().add_patch(circle)\n",
    "\n",
    "# Schaubild anzeigen\n",
    "plt.tick_params(axis='both', which='major', labelsize=9)\n",
    "plt.ylim(-10, 10)\n",
    "plt.xlim(-10, 10)\n",
    "plt.axis('equal')\n",
    "plt.show()\n",
    "\n",
    "# Resultate anzeigen\n",
    "print(\"Parameter:\")\n",
    "print(\"z1: \", res[0][0])\n",
    "print(\"z2: \", res[0][1])\n",
    "print(\"R: \", res[0][2])\n",
    "print(\"Anzahl der Interationen:\", res[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}