ALG11-Regresion_Lineal_COMPLETO.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8-  Machine Learning\n",
    "\n",
    "El nombre de este apartado es un poco pomposo pero es importante entender que el **álgebra** juega un papel fundamental en el Machine Learning y para ello vamos a ver algunos conceptos que nos van a ser útiles en un futuro.\n",
    "\n",
    "Hay muchas formas de dividir los problemas de Machine Learning, algunas de ellas:<br>\n",
    "\n",
    "**<center>Supervised Learning | Unsupervised Learning</center>**\n",
    "**<center>Probabilistic Methods | Non-probabilistic Methods</center>**\n",
    "**<center>Modeling Approach | Optimization Techniques</center>**\n",
    "**<center>...</center>**\n",
    "\n",
    "Dentro de los problemas de **Supervised Learning** podemos encontrar (nos vamos a centrar solo en este caso porque vamos a explicar algunas herramientas que nos proporciona la álgebra para resolver ciertos problemas):\n",
    "\n",
    "<img src=\"Images/reg_clas.png\" style=\"width: 75%\"/>\n",
    "\n",
    "- **Regresión**: Dado un conjunto de datos, predecir una variable continua\n",
    "- **Clasificación**: Dado un conjunto de datos, predecir un valor discreto (label/class...)\n",
    "\n",
    "## 8.1 - Regresión Lineal\n",
    "\n",
    "La **regresión lineal** es el caso más simple de **regresión**. Imaginaros que nos hemos ido a Islandia a ver géiseres (Geyser) y hemos medido durante un día cada vez que hace erupción y cuanto tiempo dura la erupción. Con estos datos nos gustaría predecir cuando va a ser la siguiente erupción.\n",
    "\n",
    "Si observamos los datos, vemos como en el **eje X** tenemos el tiempo que el géiser estuvo en erupción, y luego en el **eje y**, tenemos el tiempo que tenemos que esperar a la próxima erupción. Vemos como claramente existe una **relación entre ambas variables, ¿no?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"Images/erupcion.png\" style=\"width: 340px;\" align=\"left\"/><img src=\"Images/erupcion_graph.png\" style=\"width: 290px;\" align=\"left\"/><img src=\"Images/reg_erupcion_graph.png\" style=\"width: 265px;\" align=\"left\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Así que mirando estos datos, **podríamos querer llegar a algún tipo de función**, una simplificación\n",
    "donde simplemente tienes que tomar la entrada, **que es el tiempo que el géiser ha estado en erupción** y entonces la salida será **una predicción de cuánto tiempo creemos que vamos a tener que esperar a la próxima erupción.** Y mirando esto, la elección más natural es modelar esto con una función lineal. \n",
    "\n",
    "Este es el ajuste de regresión más simple, donde tienen una entrada que creemos que está relacionada de acuerdo a alguna función lineal con la salida. \n",
    "\n",
    "#### Función lineal\n",
    "\n",
    "El modelo para este problema será una función lineal donde tendremos que la variable $y$ es el **tiempo de espera** y la variable $x$ es la **duración de la erupción**.\n",
    "\n",
    "<center>$y = w_0 + x \\cdot w_1$</center><br>\n",
    "<center>$tiempo_{espera} = w_0 + duracion_{erupcion} \\cdot w_1$</center>\n",
    "\n",
    "$w_1$: *Pendiente / Slope*<br>\n",
    "$w_0$: *Término independiente / Bias / Offset / Shift...* \n",
    "\n",
    "Donde $w_0$ y $w_1$ deben ser determinados a partir de los datos de entrada.\n",
    "\n",
    "### ¿Cómo trasladamos esto a problemas con un mayor número de dimensiones?\n",
    "\n",
    "Imagina que tenemos dos entradas y queremos predecir una salida basada en esas dos entradas. En ese caso, podríamos seguir utilizando una función lineal donde:\n",
    "\n",
    "<center>$output = w_0 + input_{1} \\cdot w_1 + input_{2} \\cdot w_2$</center>\n",
    "\n",
    "Geométricamente podría ser representado como:\n",
    "<img src=\"Images/map_geometric.png\" style=\"width: 340px;\"/>\n",
    "\n",
    "### Caso genérico\n",
    "\n",
    "Por tanto, el objetivo de una **regresión** es:\n",
    "\n",
    "Encontrar una función $f: \\mathbb{R}^d \\rightarrow \\mathbb{R}$ tal que $y \\approx f(x;w)$ para el conjunto de datos $(x, y)$. Una regresión es denominada **lineal** si la predicción $f$ es una función lineal con parámetros no conocidos $w$.\n",
    "\n",
    "De una forma más genérica podemos escribir la **regresión** como:\n",
    "\n",
    "<center>$y_{i} \\approx  f(x_{i}; w) = w_0 + \\sum_{j=1}^d x_{ij}w_{j}$</center>\n",
    "\n",
    "*Nota:* Consideramos que tenemos más observaciones que dimensiones $d<n$\n",
    "\n",
    "## 8.2 - Mínimos Cuadrados (Least Squares)\n",
    "\n",
    "**Genial!!** Ya sabemos que es una regresión y como se clasifican algunos problemas de ML en el caso del géiser, ¿cómo podemos determinar los parámetros $w$ que nos devuelve la recta que mejor se ajusta a los datos?\n",
    "\n",
    "<img src=\"Images/interesting2.gif\" style=\"width: 340px;\"/>\n",
    "\n",
    "Si recapitulamos un poco... Tenemos un conjunto de datos $(x_1; y_1) ... (x_n, y_n)$ y queremos usar estos datos para determinar $w$, tal que $y_i \\approx f(x_i, w)$ pero antes debemos tener una **\"objective function\"** para determinar que valores de $w$ consideramos por buenos. Para ello, vamos a usar **Mínimos cuadrados**\n",
    "\n",
    "### Mínimos Cuadrados\n",
    "\n",
    "Los **mínimos cuadrados** nos ayuda a elegir los valores de $w$ que minimizan la suma de los errores al cuadrado.<br><br>\n",
    "\n",
    "<center>$\\DeclareMathOperator*{\\argmin}{argmin} w_{LS} = \\argmin _w \\sum_{i=1}^n (y_{i} - f(x_{i}; w))^2$</center><br>\n",
    "\n",
    "<center>$\\DeclareMathOperator*{\\argmin}{argmin} argmin$ = Minimizar</center>\n",
    "\n",
    "Basicamente, el objetivo de los **mínimos cuadrados $L$** es minimizar la suma de los errores al cuadrado. Si observamos el siguiente plano, la **distancia vertical es el error.**\n",
    "\n",
    "<img src=\"Images/least_squares.png\" style=\"width: 340px;\"/>\n",
    "\n",
    "#### Ejemplo\n",
    "\n",
    "Aquí hay un ejemplo muy sencillo donde tenemos como entrada un vector de 2 dimensiones con medidas de **educación** y **antigüedad** y como salida tenemos el **salario**. Basicamente queremos predecir los ingresos de una persona en base a su educación y antigüedad.\n",
    "\n",
    "<img src=\"Images/ex_least.png\" style=\"width: 340px;\" align=\"left\"/><br><br><br>\n",
    "<b>Input</b>: (Educación, Antigüedad) $\\in \\mathbb{R}^2$<br><br>\n",
    "<b>Output</b>: Salario $\\in \\mathbb{R}$<br><br>\n",
    "<b>Model</b>: $salario \\approx w_0 + (educacion) \\cdot w_1 + (antigüedad) \\cdot w_2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.3 - Mínimos Cuadrados - Matrices\n",
    "\n",
    "Existen distintas formas para la ecuación: \n",
    "\n",
    "<center>$\\DeclareMathOperator*{\\argmin}{argmin} w_{LS} = \\argmin _w \\sum_{i=1}^n (y_{i} - f(x_{i}; w))^2$</center><br>\n",
    "\n",
    "Si desarrollamos dicha fórmula podemos llegar a la ecuación de **mínimos cuadrados** en forma matricial [aquí viene el desarrollo, por si le queréis echar un ojo](http://pillowlab.princeton.edu/teaching/statneuro2018/slides/notes03b_LeastSquaresRegression.pdf):\n",
    "\\\n",
    "<center>$X^{T}X W'=X^{T}y$</center><br>\n",
    "\n",
    "O lo que es lo mismo:<br>\n",
    "<center>$W'=(X^{T}X)^{-1}X^{T}y$</center>\n",
    "    \n",
    "Donde la matrix $X$, viene definida como:\n",
    "\n",
    "<center>$\n",
    "  X =\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   1 & x_{11} & ... & x_{1d} \\\\\n",
    "   1 & x_{21} & ... & x_{2d} \\\\\n",
    "   ... & ... & ... & ... \\\\\n",
    "   1 & x_{nd} & ... & x_{nd} \\\\\n",
    "  \\end{array} } \\right]\n",
    "$</center><br>\n",
    "\n",
    "Se le añade un $1$ a cada fila de la matriz por convención y por la componente $w_0$ de la ecuación - Indicado por $\\mathbb{R}^{d+1}$\n",
    "\n",
    "### 8.3.1 - ¿Por qué n > d?\n",
    "\n",
    "$n > d$ indica que el **número de observaciones** sea mayor que el **el número de dimensiones**, vamos a poner un ejemplo para entenderlo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/.pyenv/versions/3.8.3/lib/python3.8/site-packages/pandas/compat/__init__.py:124: UserWarning: Could not import the lzma module. Your installed Python is incomplete. Attempting to use lzma compression will result in a RuntimeError.\n",
      "  warnings.warn(msg)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Name</th>\n",
       "      <th>Media bachiller - Grade</th>\n",
       "      <th>Media selectividad - Test</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Teresa</td>\n",
       "      <td>8.7</td>\n",
       "      <td>7.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Alex</td>\n",
       "      <td>6.7</td>\n",
       "      <td>8.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Sara</td>\n",
       "      <td>9.0</td>\n",
       "      <td>9.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Pepe</td>\n",
       "      <td>7.0</td>\n",
       "      <td>5.1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     Name  Media bachiller - Grade  Media selectividad - Test\n",
       "0  Teresa                      8.7                        7.0\n",
       "1    Alex                      6.7                        8.0\n",
       "2    Sara                      9.0                        9.1\n",
       "3    Pepe                      7.0                        5.1"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "data = [['Teresa', 8.7, 7], ['Alex', 6.7, 8], ['Sara', 9, 9.1], ['Pepe', 7, 5.1]]\n",
    "df = pd.DataFrame(data, columns=['Name', 'Media bachiller - Grade', 'Media selectividad - Test'])\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El número `n` es el número de observaciones y el número `d` son las dimensiones. Es decir, si en nuestro problema tenemos **3 personas** a estudiar y vamos a analizar la influencia de **media de bachiller** y de su **media en selectividad**, tendremos que:\n",
    "\n",
    "- Número observaciones `n`: 3 personas (Teresa, Alex, Sara)\n",
    "- Número de dimensiones `d`: 2 variables (Media bachiller, Media selectividad)\n",
    "\n",
    "<img src=\"./images/nd_x.jpg\" width=\"400\">\n",
    "\n",
    "En el caso de que $n < d+1$ no podemos aplicar **mínimos cuadrados**, ya que $(X^TX)^{-1}$ no existe y habría un número infinito de soluciones.\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "    <b>Calcula $(X^TX)^{-1}$</b><br><br>\n",
    "    \n",
    "<center>$\n",
    "  X =\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   1 & 8.7 & 7 \\\\\n",
    "   1 & 6.7 & 8 \\\\\n",
    "   1 & 9 & 9.1 \\\\\n",
    "   1 & 7 & 5.1 \\\\\n",
    "  \\end{array} } \\right]\n",
    "$</center><br>\n",
    "</div> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.  8.7 7. ]\n",
      " [1.  6.7 8. ]\n",
      " [1.  9.  9.1]\n",
      " [1.  7.  5.1]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "X = np.array([[1, 8.7, 7],\n",
    "              [1, 6.7, 8],\n",
    "              [1, 9, 9.1],\n",
    "              [1, 7, 5.1]])\n",
    "\n",
    "print(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1. , 1. , 1. , 1. ],\n",
       "       [8.7, 6.7, 9. , 7. ],\n",
       "       [7. , 8. , 9.1, 5.1]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_t = np.transpose(X)\n",
    "X_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[15.79028203, -1.73113364, -0.26724424],\n",
       "       [-1.73113364,  0.31926267, -0.10617512],\n",
       "       [-0.26724424, -0.10617512,  0.15078341]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prod = np.dot(X_t, X)\n",
    "inv_prod = np.linalg.inv(prod)\n",
    "inv_prod"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>Calcula $(X^TX)^{-1}$</b><br><br>\n",
    "    \n",
    "<center>$\n",
    "  X =\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   1 & 8.7 & 7 \\\\\n",
    "   1 & 6.7 & 8 \\\\\n",
    "  \\end{array} } \\right]\n",
    "$</center><br>\n",
    "</div> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 6.38948197e+15, -2.81474977e+14, -5.62949953e+14],\n",
       "       [-2.81474977e+14,  1.23997787e+13,  2.47995574e+13],\n",
       "       [-5.62949953e+14,  2.47995574e+13,  4.95991148e+13]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.array([[1, 8.7, 7],\n",
    "              [1, 6.7, 8]])\n",
    "\n",
    "X_t = np.transpose(X)\n",
    "prod = np.dot(X_t, X)\n",
    "inv_prod = np.linalg.inv(prod)\n",
    "inv_prod"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3.2 - La filosofía detrás de los mínimos cuadrados\n",
    "\n",
    "¿Os acordais de cuando hablamos de los sistemas de ecuaciones? \n",
    "\n",
    "<img src=\"./Images/secuaciones.png\" width=50%>\n",
    "\n",
    "Cuando un sistema era incompatible significaba que tenía infinitas soluciones, pero hay una forma de encontrar una solución aproximada y aquí es donde entra en juego: **Los Mínimos Cuadrados**.\n",
    "\n",
    "**Supongamos que tratamos de resolver un sistema de la forma Ax = b** y que no existe solución para este sistema. En este caso, puede ser interesante estudiar para qué vector o vectores $x'$, $Ax' ≈ b$, esto es, tratamos de buscar una solución aproximada del sistema.\n",
    "\n",
    "Una forma de interpretar esta aproximación es la siguiente:\n",
    "Si Ax' ≈ b, entonces podemos intentar encontrar x' tal que kAx' − bk sea lo menor posible.\n",
    "Es decir, tratamos de resolver el problema minx ‖Ax − b‖ \n",
    "\n",
    "Si el sistema tiene solución, entonces el mínimo del problema anterior es cero, y por tanto Ax = b; pero si no hay solución, tiene sentido buscar el vector x que haga menor esa norma, y por tanto que más se parece a una solución.\n",
    "\n",
    "Al final, al igual que hemos hecho antes, desarrollamos dicha expresión, se llega a la conclusión de que se trata de resolver el sistema:\n",
    "\n",
    "$$A^TAx'=A^Tb$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.4 - \"Ampliando\" la Regresión Lineal\n",
    "\n",
    "Lo bueno de la **regresión lineal** y de los **mínimos cuadrados** es que es muy versátil. Seguro que pensáis que calcular una función lineal es un problema super simple, teniendo en cuenta que solo estamos dibujando una línea a través de nuestro hiperplano, pero eso no es así. Las funciones lineales son la base de **TODO**, incluso de las redes neuronales.\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "    <b>Os he dejado un extra sobre los perceptrones: EXTRA5_Perceptron</b>\n",
    "</div> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.4.1 - \"Ampliando\" las regresiones lineales\n",
    "\n",
    "Bueno, volvamos a las **regresiones lineales**😃. Como hemos dicho, las **regresiones lineales** y los **mínimos cuadrados** son herramientas muy potentes. Recuerda que la **regresión lineal** se llama regresión lineal porque la salida es lineal:<br><br>\n",
    "\n",
    "<center>$y_{i} \\approx  f(x_{i}; w) = w_0 + \\sum_{j=1}^d x_{ij}w_{j}$</center>\n",
    "\n",
    "Así que, tenemos una función lineal de los pesos $w$, no necesariamente de las entradas $x$. En el siguiente ejemplo, la regresión lineal va a ser algo que **NO** queramos usar, porque no va a modelar correctamente nuestros datos, ¿qué hacemos entonces?:\n",
    "\n",
    "<img src=\"./images/pol_regresion.png\" width=50%>\n",
    "\n",
    "Si miramos la ecuación original de la regresión lineal, tenemos la salida $y$, que es aproximadamente igual a la entrada $x$\n",
    "multiplicado por un peso $w_1$, más un desplazamiento $w_0$. Básicamente estamos diciendo que los datos tienen que ser modelados por una línea, pero no es lo ideal en este caso.\n",
    "\n",
    "Pero podríamos decir ahora que nuestra salida $y$ sea de la forma:\n",
    "\n",
    "<img src=\"./images/pol_regresion2.png\" width=50%>\n",
    "\n",
    "Esto es lo que se denomina `regresión polinomial` usando una función de 3º orden. Lo mejor de todo es que se resuelve de la misma forma que la `regresión lineal`, usando: \n",
    "\n",
    "<center>$X^{T}X W'=X^{T}y$</center><br>\n",
    "\n",
    "O lo que es lo mismo:<br>\n",
    "<center>$W'=(X^{T}X)^{-1}X^{T}y$</center>\n",
    "    \n",
    "La única diferencia es que nuestra matriz $X$ con un polinomio de grado $p$, será de la forma:\n",
    "\n",
    "<div class=\"alert alert-block alert-danger\">\n",
    "    <b>Nota:</b> $d=1$ indica que solo tiene una variable, por ejemplo: edad / dinero; ...\n",
    "</div><br><br>\n",
    "\n",
    "<center>$\n",
    "  X =\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   1 & x_1 & x_1^2 & ... & x_n^p \\\\\n",
    "   1 & x_2 & x_2^2 & ... & x_n^p \\\\\n",
    "   1 & ... & ... & ... & ... \\\\\n",
    "   1 & x_n & x_n^2 & ... & x_n^p \\\\\n",
    "  \\end{array} } \\right]\n",
    "$</center>\n",
    "\n",
    "Donde `p` será el orden de nuestro polinomio.\n",
    "\n",
    "#### Influencia del grado del polinomio en nuestra regresión\n",
    "\n",
    "<img src=\"./images/gif-pol-reg.gif\" width=50%>\n",
    "\n",
    "### - Regresión polinomial en $\\mathbb{R}^2$ (2 dimensiones)\n",
    "\n",
    "Como hemos visto la regresión polinomial se resuelve de la misma forma que una regresión linea, la única diferencia es que tendremos más componentes en nuestra matriz $X$. El `número de columnas` vendrá dado por:\n",
    "\n",
    "<center>$(orden \\,\\, polinomio)\\, x \\,(dimensiones) + 1$</center><br>\n",
    "\n",
    "En ese caso para un **polinomio de orden 2**:<br>\n",
    "<center>$\n",
    "  y_i = w_0 + w_1x_{i1} + w_2x_{i2} + w_3x_{i1}^2 + w_4x_{i2}^2\n",
    "$</center>\n",
    "\n",
    "Para un **polinomio de orden 3**:<br>\n",
    "\n",
    "<center>$\n",
    "  y_i = w_0 + w_1x_{i1} + w_2x_{i2} + w_3x_{i1}^2 + w_4x_{i2}^2 + w_5x_{i1}^3 + w_6x_{i2}^3\n",
    "$</center>\n",
    "\n",
    "<img src=\"./images/pol_3order.png\" width=50%>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.4.2 - Resumen\n",
    "\n",
    "#### ¿Qué es la regresión?\n",
    "\n",
    "La `regresión` es un modelo matemático que busca determinar la relación entre una variable dependiente $y$ con respecto a otras variables independientes $x_n$.\n",
    "\n",
    "En la regresión, dado un conjunto de datos, nos va a predecir una variable continua. \n",
    "\n",
    "<img src=\"./images/regresion.jpg\" width=50%>\n",
    "\n",
    "#### ¿Cuál es el objetivo de la regresión?\n",
    "\n",
    "El objetivo de una **regresión** es encontrar una función $f: \\mathbb{R}^d \\rightarrow \\mathbb{R}$ tal que $y \\approx f(x;w)$ para el conjunto de datos $(x, y)$. Una regresión es denominada **lineal** si la predicción $f$ es una función lineal con parámetros no conocidos $w$.\n",
    "\n",
    "De una forma más genérica podemos escribir la **regresión** como:\n",
    "\n",
    "<center>$y_{i} \\approx  f(x_{i}; w) = w_0 + \\sum_{j=1}^d x_{ij}w_{j}$</center>\n",
    "\n",
    "*Nota:* Consideramos que tenemos más observaciones que dimensiones $n > d$<br><br>\n",
    "<center>$\n",
    "  X =\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   1 & x_{11} & ... & x_{1d} \\\\\n",
    "   1 & x_{21} & ... & x_{2d} \\\\\n",
    "   ... & ... & ... & ... \\\\\n",
    "   1 & x_{nd} & ... & x_{nd} \\\\\n",
    "  \\end{array} } \\right]\n",
    "$</center>\n",
    "\n",
    "#### ¿Cómo la resolvemos?\n",
    "\n",
    "Usando **mínimos cuadrados** cuyo objetivo es encontrar los valores de $w$ que minimizan la suma de los errores al cuadrado.\n",
    "Existen distintas formas para la ecuación: \n",
    "\n",
    "<center>$\\DeclareMathOperator*{\\argmin}{argmin} w_{LS} = \\argmin _w \\sum_{i=1}^n (y_{i} - f(x_{i}; w))^2$</center><br>\n",
    "\n",
    "Si desarrollamos dicha fórmula podemos llegar a la ecuación de **mínimos cuadrados** en forma matricial [aquí viene el desarrollo, por si le queréis echar un ojo](http://pillowlab.princeton.edu/teaching/statneuro2018/slides/notes03b_LeastSquaresRegression.pdf):\n",
    "\n",
    "<center>$X^{T}X W'=X^{T}y$</center><br>\n",
    "\n",
    "O lo que es lo mismo:<br>\n",
    "<center>$W'=(X^{T}X)^{-1}X^{T}y$</center>\n",
    "    \n",
    "    \n",
    "#### Regresion polinomial\n",
    "\n",
    "Para casos en los que no podamos modelar nuestro problema con una **función lineal** podemos aplicar un polinomio para ello. En este caso se llamaría **regresión polinomial** y se resuelve de la misma forma:\n",
    "\n",
    "<center>$X^{T}X W'=X^{T}y$</center><br>\n",
    "<center>$W'=(X^{T}X)^{-1}X^{T}y$</center>\n",
    "\n",
    "#### Matriz X\n",
    "\n",
    "La única diferencia entre la regresión lineal y la regresión polinomial es como escribimos la matriz $X$, aquí algunos ejemplos:\n",
    "\n",
    "- Problema con 2 dimensiones\n",
    "\n",
    "<center>$\n",
    "  X =\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   1 & x_{1d} & x_{1d} & ... & x_{nd} \\\\\n",
    "   1 & x_{2d} & x_{2d} & ... & x_{nd} \\\\\n",
    "   1 & ... & ... & ... & ... \\\\\n",
    "   1 & x_{nd} & x_{nd} & ... & x_{nd} \\\\\n",
    "  \\end{array} } \\right] = \\left[ {\\begin{array}{cc}\n",
    "   1 & x_{11} & x_{12}\\\\\n",
    "   1 & x_{21} & x_{22}\\\\\n",
    "   1 & ... & ...\\\\\n",
    "   1 & x_{n1} & x_{n2}\\\\\n",
    "  \\end{array} } \\right] \n",
    "$</center>\n",
    "\n",
    "- Problema con 3 dimensiones y 4 observaciones\n",
    "\n",
    "<center>$\n",
    "  X =\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   1 & x_{1d} & x_{1d} & x_{1d} &... & x_{nd} \\\\\n",
    "   1 & x_{2d} & x_{2d} & x_{2d} &... & x_{nd} \\\\\n",
    "   1 & x_{3d} & x_{3d} & x_{3d} &... & x_{nd} \\\\\n",
    "   1 & x_{4d} & x_{4d} & x_{4d} &... & x_{nd} \\\\\n",
    "  \\end{array} } \\right] = \\left[ {\\begin{array}{cc}\n",
    "   1 & x_{11} & x_{12} & x_{13} \\\\\n",
    "   1 & x_{21} & x_{22} & x_{23} \\\\\n",
    "   1 & x_{31} & x_{32} & x_{33} \\\\\n",
    "   1 & x_{41} & x_{42} & x_{43} \\\\\n",
    "  \\end{array} } \\right] \n",
    "$</center>\n",
    "\n",
    "- Problema con 4 observaciones y 2 dimensiones, modelado con un polinomio de 2º grado\n",
    "\n",
    "<center>$\n",
    "  X =\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   1 & x_{1d} & x_{1d} & x_{1d}^2 &... & x_{nd}^p \\\\\n",
    "   1 & x_{2d} & x_{2d} & x_{2d}^2 &... & x_{nd}^p \\\\\n",
    "   1 & x_{3d} & x_{3d} & x_{3d}^2 &... & x_{nd}^p \\\\\n",
    "   1 & x_{4d} & x_{4d} & x_{4d}^2 &... & x_{nd}^p \\\\\n",
    "  \\end{array} } \\right] = \\left[ {\\begin{array}{cc}\n",
    "   1 & x_{11} & x_{12} & x_{11}^2 & x_{12}^2 \\\\\n",
    "   1 & x_{21} & x_{22} & x_{21}^2 & x_{22}^2 \\\\\n",
    "   1 & x_{31} & x_{32} & x_{31}^2 & x_{32}^2 \\\\\n",
    "   1 & x_{41} & x_{42} & x_{41}^2 & x_{42}^2 \\\\\n",
    "  \\end{array} } \\right] \n",
    "$</center>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.5 Gradient Descent\n",
    "\n",
    "Hemos visto como **Least Squares/Mínimos Cuadrados** nos permite encontrar la recta que mejor aproxima nuestros datos. Los métodos que nos permiten encontrar la solución mínima o máxima se llaman **métodos de optimización**. Existe numerosos tipos de métodos de optimización, pero nos vamos a centrar en uno de los más usados hoy en día en el mundo del Machine Learning - **Gradient Descent**.\n",
    "\n",
    "**Gradient Descent** es un método de optimización iterativo (cada iteración que ejecutamos estamos un poco más cerca del mínimo de la función - solución). \n",
    " \n",
    "\n",
    "<center>$y_{i} \\approx  f(x_{i}; w) = w_0 + \\sum_{j=1}^d x_{ij}w_{j}$</center>\n",
    "\n",
    "Podemos convertirlo a matrices:\n",
    "\n",
    "<center>$y_{i} \\approx  f(x_{i}; w) = XW$</center>\n",
    "\n",
    "Y queremos encontrar:\n",
    "\n",
    "* **Los argumentos $w$ que minimiza la función** - $\\underset{w}{\\operatorname{argmin}}f(w)$\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "    <b>Nota</b> $\\min f(x)$ y $\\arg \\min f(x)$ son cosas distintas. Un ejemplo:\n",
    "</div> \n",
    "\n",
    "Considera la función $f(x) = 100 + (x - 3)^2$\n",
    "\n",
    "$\\min f(x) = f(\\arg \\min f(x))$\n",
    "\n",
    "Entonces:\n",
    "\n",
    "$\\frac{df(x)}{x} = 2(x - 3)$<br>\n",
    "$\\arg \\min 2(x - 3) = 0$<br>\n",
    "$\\arg \\min f(x) = 3$\n",
    "\n",
    "Por otro lado:\n",
    "\n",
    "$\\min f(x) = f(\\arg \\min f(x)) = 100 + (3 - 3)^2$<br>\n",
    "$\\min f(x) = 100$\n",
    "\n",
    "Podemos graficarlo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x13153d7c0>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD4CAYAAADsKpHdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAm9klEQVR4nO3dd3yV9d3G8c83m0AGkEAChCF7rwhWQBxVUangQHGiaN1tHa2z2j4dT63W+rgFFamKsw5wixOrrKCsYICwE0YCgRAICQn5PX+cgwYaWebkPuN6v155cc7vvu9zrgN6rtzbnHOIiIhEeR1ARESCgwpBREQAFYKIiPipEEREBFAhiIiIX4zXAY5UWlqaa9++vdcxRERCyrx58zY759LrmhayhdC+fXtycnK8jiEiElLMbM2PTdMmIxERAVQIIiLip0IQERFAhSAiIn4qBBERAVQIIiLip0IQEREgAguhfHc1f/8gj3Ul5V5HEREJKhFXCKW7qnju69XcPXUxuheEiMgPIq4QMlMacfMpXfl8aTHvL97odRwRkaARcYUAMO5n7ejZKpn/eTuXsooqr+OIiASFiCyEmOgo/npWb4rKKnngo2VexxERCQoRWQgA/bJSueSYdjw3czWLCkq9jiMi4rmILQSA357aleZN4rnjzYVU76nxOo6IiKciuhCSE2L54y96srhwO5O/Xu11HBERT0V0IQCc3juDE7qm88/pyyjctsvrOCIinon4QjAz/jSqF87BH3RugohEsIAVgplNMrMiM1u83/ivzCzPzHLN7L5a43eYWb6ZLTWzUwOVqy5ZzRK5+eQufPxdER/o3AQRiVCBXEOYDIyoPWBmJwCjgL7OuZ7AP/zjPYCxQE//Mo+bWXQAs/2Xy4e0p2erZP4wLZfSXTo3QUQiT8AKwTk3AyjZb/ha4F7nXKV/niL/+CjgZedcpXNuFZAPDApUtrrEREdx79l92Lyjkr9/kNeQby0iEhQaeh9CF2CYmc02sy/M7Gj/eGtgXa35CvxjDap3mxTGD+nAi7PXMmfV/l0mIhLeGroQYoBmwDHA74BXzcwOdWEzu8rMcswsp7i4OCABbz6lC22aNuKONxZSWb0nIO8hIhKMGroQCoA3nM8coAZIAwqBrFrztfGP7cM5N9E5l+2cy05PTw9IwMS4GP4yuhcrinfy2GcrAvIeIiLBqKEL4S3gBAAz6wLEAZuBacBYM4s3sw5AZ2BOA2f73vFdWzC6Xyue+DyfpRvLvIohItKgAnnY6UvATKCrmRWY2RXAJOAo/6GoLwPj/GsLucCrwBLgA+B655yn22vuHtmDpIRYbnt9IXtqdG6CiIQ/C9UTsbKzs11OTk5A32Pq/EJ+8/J87h7ZgyuGdgjoe4mINAQzm+ecy65rWsSfqXwgZ/ZtxQld0/nHh0t1y00RCXsqhAMwM/5yVm+iDO58c5EuayEiYU2FcBCtUxtx+2nd+HL5Zl6bV+B1HBGRgFEhHIKLBrdjUIdm/OWdJRRtr/A6johIQKgQDkFUlPH3c/pQWV3DXW/piqgiEp5UCIeoQ1pjbj65C9OXbOKdhRu8jiMiUu9UCIfhiqEd6NsmhT9My2XLjkqv44iI1CsVwmGIiY7ivnP7UlZRxR+m5XodR0SkXqkQDlPXjCR+c1Jn3lm4gQ8Wa9ORiIQPFcIRuHp4R3q2Sub3by1m687dXscREakXKoQjEBsdxT/G9GVbeRV/fFubjkQkPKgQjlD3zGRuOLETU+ev58Nc3YdZREKfCuEnuP6ETvTITOauNxdRok1HIhLiVAg/QWx0FA+c15fSXTrqSERCnwrhJ+qemcyvT+zM2wvW8/4iHXUkIqFLhVAPrjm+I71bp/D7txazWSesiUiIUiHUg72bjsoqqvn9m7rWkYiEJhVCPenSMolbTunCB7kbeWt+oddxREQOmwqhHl057CgGtmvKPVNz2VC6y+s4IiKHRYVQj6KjjAfG9KV6j+O213WHNREJLSqEetY+rTF3nN6NGcuKmTJ7rddxREQOmQohAC45ph3DOqfx13e/Y9XmnV7HERE5JAErBDObZGZFZra41tgfzazQzOb7f06vNe0OM8s3s6VmdmqgcjUEM+P+c/sSG23c/Op8qvfUeB1JROSgArmGMBkYUcf4g865fv6f9wDMrAcwFujpX+ZxM4sOYLaAy0hJ4M+je/Ht2m1MmLHS6zgiIgcVsEJwzs0ASg5x9lHAy865SufcKiAfGBSobA1lVL/WjOyTyYPTl7G4sNTrOCIiB+TFPoQbzGyhf5NSU/9Ya2BdrXkK/GP7MLOrzCzHzHKKi4sbIutP9pfRvUhrEs+Nr8ynomqP13FERH5UQxfCE0BHoB+wAXjgcBZ2zk10zmU757LT09MDEK/+pSbG8Y8xfckv2sG97+d5HUdE5Ec1aCE45zY55/Y452qAp/hhs1AhkFVr1jb+sbAwtHMalw9pz+SvVzNjWWis2YhI5GnQQjCzzFpPzwL2HoE0DRhrZvFm1gHoDMxpyGyBdtuIbnRu0YTfvrZAt90UkaAUyMNOXwJmAl3NrMDMrgDuM7NFZrYQOAG4CcA5lwu8CiwBPgCud86F1Qb3hNho/m9sP7aW7+b2NxbqLGYRCToWql9M2dnZLicnx+sYh23ijBX873t5/P2c3px/dFuv44hIhDGzec657Lqm6UzlBnbl0KM4tmNz/jhtCSuLd3gdR0TkeyqEBhYVZTxwXl/iYqK48ZX5VOksZhEJEioED2SmNOLes3uzsKCUf05f5nUcERFAheCZ03pnMvboLJ78YgVf52/2Oo6IiArBS/f8ogcd0hpz06vzdSiqiHhOheChxLgYHh7bn5Kdu7n1dR2KKiLeUiF4rFfrFG4b0Y3pSzbxgm6oIyIeUiEEgfFDOjC8Szp/fmcJeRu3ex1HRCKUCiEI7D0UNaVRLDe8+C3lu6u9jiQiEUiFECTSmsTz4Hn9WFG8gz+9vcTrOCISgVQIQWRo5zSuHd6Rl+eu4+0F672OIyIRRoUQZG46uQsD2qZyxxuLWLNlp9dxRCSCqBCCTGx0FA9f0J/oKOOGF7+lsjqsLvoqIkFMhRCE2jRN5P5z+7CosFR3WRORBqNCCFKn9Mzg8iHtefar1XyYu9HrOCISAVQIQez207rRp00Kv3ttAetKyr2OIyJhToUQxOJjonnswgE44PoXv9H+BBEJKBVCkMtqlsj95/ZlYUEpf3tP+xNEJHBUCCFgRK8Mxg/pwOSvV/Puwg1exxGRMKVCCBG3n9aNflmp3Pb6Qt16U0QCQoUQIuJionjsogHERhvXTfmGXbu1P0FE6pcKIYS0Tm3Eg+f3Y+mmMn7/1mLdP0FE6lXACsHMJplZkZktrmPaLWbmzCzN/9zM7GEzyzezhWY2IFC5Qt3xXVvwqxM78/o3Bbwyd53XcUQkjARyDWEyMGL/QTPLAk4Bat8N5jSgs//nKuCJAOYKeb85qTNDO6Vxz7RcFhWUeh1HRMJEwArBOTcDKKlj0oPArUDt7R2jgOeczywg1cwyA5Ut1EVHGQ+N7Uda4ziunTKPbeW6H7OI/HQNug/BzEYBhc65BftNag3U3v5R4B/bf/mrzCzHzHKKi4sDmDT4NW8Sz+MXD6RoeyU3vjKfmhrtTxCRn6bBCsHMEoE7gXuO9DWccxOdc9nOuez09PT6Cxei+mWlcs8vevD50mIe+mS513FEJMQ15BpCR6ADsMDMVgNtgG/MLAMoBLJqzdvGPyYHcdHgtpwzoA0PfbKcT77b5HUcEQlhDVYIzrlFzrkWzrn2zrn2+DYLDXDObQSmAZf6jzY6Bih1zumU3ENgZvz1rF70bJXMja/MZ/Vm3VRHRI5MIA87fQmYCXQ1swIzu+IAs78HrATygaeA6wKVKxwlxEbz5MUDiYkyrn5+Hjsrq72OJCIhyEL15Kbs7GyXk5PjdYyg8uXyYsZNmsNpvTN59IL+mJnXkUQkyJjZPOdcdl3TdKZyGBnWOZ1bR3Tj3YUbmDBjpddxRCTEqBDCzNXHHcXIPpn8/YM8Pl9a5HUcEQkhKoQwY2bcd24furZM4tcvfaudzCJyyFQIYSgxLoanLs0mKsq46vkcdmgns4gcAhVCmMpqlshjFw5gRfFObtKZzCJyCFQIYWxIpzTuPqM705ds4sGPl3kdR0SCXIzXASSwxh3bnu82lPHIp/l0zUhiZJ9WXkcSkSClNYQwZ2b8aXRPsts15bevLWBxoS6XLSJ1UyFEgPiYaJ64eCDNG8dz5b9yKNpe4XUkEQlCKoQIkZ4Uz1OXZlO6q4pfPj+Piirdk1lE9qVCiCA9WiXz4Pn9WLBuG7e9vlD3ZBaRfagQIsyIXhn87tSuTJ2/nkc+zfc6jogEER1lFIGuO74jK4p28M/py+iQ1phf9NWRRyKiNYSIZGb87ZzeHN3ed+TRt2u3eh1JRIKACiFCxcdEM+GSbFomJ/DL5+ZRsLXc60gi4jEVQgRr1jiOSZdlU1m9h/GT57K9osrrSCLiIRVChOvUIokJFw9kZfFOrnvhG6r21HgdSUQ8okIQju2Uxt/O7s1/8jfz+zcX63BUkQilo4wEgDHZWawtKeeRT/Np2zyR60/o5HUkEWlgKgT53s0nd6Fg6y7u/3AprVITOKt/G68jiUgDUiHI98yMv5/Th42lFdz674W0TE7g2I5pXscSkQaifQiyj7iYKJ68ZCDtmzfm6ufnsXRjmdeRRKSBBKwQzGySmRWZ2eJaY382s4VmNt/MPjKzVv5xM7OHzSzfP31AoHLJwaU0imXy+EEkxkUzbtIc1m/b5XUkEWkAgVxDmAyM2G/sfudcH+dcP+Ad4B7/+GlAZ//PVcATAcwlh6B1aiMmXz6InZXVXPbsHErLdY6CSLgLWCE452YAJfuNba/1tDGw9/jGUcBzzmcWkGpmmYHKJoeme2YyEy4ZyKrNO/nl8zm6ZLZImGvwfQhm9lczWwdcxA9rCK2BdbVmK/CP7b/sVWaWY2Y5xcXFgQ8rHNspjX+M6cucVSXc+PJ89tToHAWRcNXgheCcu8s5lwVMAW44zGUnOueynXPZ6enpgQko/2VUv9bcPbIHH+Ru5O6pOnFNJFx5eZTRFOAc/+NCIKvWtDb+MQkSVwztwLXHd+TF2Wt58OPlXscRkQBo0EIws861no4C8vyPpwGX+o82OgYodc5taMhscnC3ntqVMQPb8PAny3lu5mqv44hIPQvYiWlm9hJwPJBmZgXAH4DTzawrUAOsAa7xz/4ecDqQD5QDlwcqlxw5M+NvZ/dma3kVf5iWS0qjWEb1+69dPSISoixUtwdnZ2e7nJwcr2NEpIqqPYybNId5a7by1KXZnNCthdeRROQQmdk851x2XdN0prIctoTYaJ4el023zCSunTKPuatLDr6QiAQ9FYIckaSEWP51+SBapTZi/LNzWVxY6nUkEfmJVAhyxJo3ieeFKwaT3CiWSyfNIb9I1z0SCWUqBPlJWqU2YsqVg4mOMi56ejbrSnRvZpFQpUKQn6x9WmNeuGIwFVU1XPj0LDaU6mJ4IqFIhSD1omtGEs9fMYhtO6u46KnZFJVVeB1JRA6TCkHqTZ82qTx7+dFs3F7BxU/PpmTnbq8jichhUCFIvcpu34ynx2WzZks5Fz89m23lKgWRUKFCkHp3bMc0Jl6aTX7RDi55Zg6lu3QvBZFQoEKQgBjeJZ0nLxlA3sbtXDppDtsrVAoiwU6FIAFzYreWPHbhAHILSxk3aQ5lKgWRoKZCkIA6pWcGj17Yn0UFpVyqUhAJaioECbgRvTJVCiIhQIUgDWL/UtA+BZHgo0KQBjOiVyaPXTSAxYWlXPz0bErLVQoiwUSFIA3q1J4ZPHnxQPI2lHHBU7N08ppIEFEhSIM7qXtLnhqXzYriHVwwcZYucyESJFQI4onhXdJ59rKjWVtSztgJuiCeSDBQIYhnju2UxvNXDKK4rJIxT85k7RZdOlvESyoE8VR2+2ZM+eVgdlRWM2bC1yzfpJvsiHhFhSCe69MmlZevOoYaB+dNmMnCgm1eRxKJSCoECQrdMpJ57eqf0Tg+hgufms2slVu8jiQScQJWCGY2ycyKzGxxrbH7zSzPzBaa2Ztmllpr2h1mlm9mS83s1EDlkuDVPq0x/77mWDJSEhg3aQ7Tl2zyOpJIRAnkGsJkYMR+Y9OBXs65PsAy4A4AM+sBjAV6+pd53MyiA5hNglRGSgKvXv0zumUmc80L83gtZ53XkUQiRsAKwTk3AyjZb+wj51y1/+ksoI3/8SjgZedcpXNuFZAPDApUNgluzRrH8eKVgzm2Y3N+9++FTPhihdeRRCKCl/sQxgPv+x+3Bmr/KljgH9uHmV1lZjlmllNcXNwAEcUrjeNjeGbc0Yzsk8nf3s/jz+8soabGeR1LJKzFePGmZnYXUA1MOZzlnHMTgYkA2dnZ+nYIc3ExUTw8tj/pSfE8859VFJVV8o8xfYiP0dZEkUBo8EIws8uAkcBJzrm9X+qFQFat2dr4xyTCRUUZ94zsQUZyAn97P4/NZZVMuHQgyQmxXkcTCTsNusnIzEYAtwJnOudqn5Y6DRhrZvFm1gHoDMxpyGwSvMyMq4d35J/n9WXu6hLGPDGT9dt0qQuR+hbIw05fAmYCXc2swMyuAB4FkoDpZjbfzJ4EcM7lAq8CS4APgOudc3sClU1C09kD2vCv8YNYv20XZz3+FUvWb/c6kkhYsR+22oSW7Oxsl5OT43UM8UDexu1c/uxctu+q4rGLBnB81xZeRxIJGWY2zzmXXdc0naksIadbRjJvXjeEds0bM37yXJ6fudrrSCJhQYUgISkjJYHXrvkZJ3Zrwd1Tc/nT20vYo8NSRX4SFYKErMbxMUy4JJvxQzow6atVXPmvuZTpXs0iR0yFICEtOsq45xc9+MvoXsxYvplznvha91UQOUIqBAkLFx/TjufGD2LT9kpGP/4Vs3W1VJHDpkKQsDGkUxpvXT+E1MRYLnp6Ni/MWuN1JJGQokKQsNIhrTFvXT+EYZ3T+P1bi7nzzUXsrq7xOpZISFAhSNhJTojl6XFHc83wjrw4ey0XPjWLou0VXscSCXoqBAlL0VHG7ad14+EL+pO7fjsjH/kP89aUHHxBkQimQpCwdmbfVrxx3bEkxEYzduIsnp+5mlA9O18k0FQIEva6Zybz9g1DGdopjbun5nLTK/Mp31198AVFIowKQSJCSmIsz4w7mltO7sLUBesZ/dhXrCje4XUskaCiQpCIERVl/Oqkzjw3fhCbd+zmzEf+w9T5uu2GyF4qBIk4wzqn8+6vh9I9M5nfvDyfO99cREWVrrYuokKQiJSZ0oiXrjrm+0NTRz/2FflFZV7HEvGUCkEiVmx0FLef1o1nLzuaorJKfvHIV7w6d52OQpKIpUKQiHdCtxa8/5th9MtK5dbXF/Krl76ldJeumiqRR4UgArRMTuCFKwfzu1O78sHijZz+0Je6QJ5EHBWCiF90lHH9CZ3497XHEhttXPDULO7/ME/XQpKIoUIQ2U+/rFTe/fUwxgzM4rHPVnDW41+xbJN2OEv4UyGI1KFxfAx/P7cPEy8ZyMbSCkY+8h8mzlih23RKWFMhiBzAKT0z+PCm4zi+Szr/+14e502YyUqd4SxhKmCFYGaTzKzIzBbXGhtjZrlmVmNm2fvNf4eZ5ZvZUjM7NVC5RA5XWpN4JlwykAfP78vyTWWc9tCXPP3lSq0tSNgJ5BrCZGDEfmOLgbOBGbUHzawHMBbo6V/mcTOLDmA2kcNiZpzVvw3Tbx7O0E5p/OXd7xjz5Ncs174FCSMBKwTn3AygZL+x75xzS+uYfRTwsnOu0jm3CsgHBgUqm8iRapmcwNPjsvm/8/uxavNOznj4PzzyyXIdiSRhIVj2IbQG1tV6XuAf24eZXWVmOWaWU1xc3GDhRGozM0b3b830m4dzSs+WPDB9GSMf+ZKc1boBj4S2YCmEQ+Kcm+icy3bOZaenp3sdRyJcWpN4Hr1wAJMuy2Zn5R7OfXImd765iNJyneUsoSlYCqEQyKr1vI1/TCTonditJR/ddBxXDu3AK3PXceIDn/Najq6JJKEnWAphGjDWzOLNrAPQGZjjcSaRQ9Y4Pobfj+zBtBuG0K55Ir/790LOmzCTJeu3ex1N5JAF8rDTl4CZQFczKzCzK8zsLDMrAH4GvGtmHwI453KBV4ElwAfA9c45XaBeQk7PVin8+5pjue+cPuQX7WDkI19y91uL2Va+2+toIgdlobpam52d7XJycryOIfKjSsur+Of0pTw/aw3JjWK5+eQuXDioLTHRwbJiLpHIzOY557Lrmqb/MkUCJCUxlv8Z1Yt3fz2M7hnJ3DM1lxEPfcnnS4u8jiZSJxWCSIB1z0zmxV8OZuIlA6neU8Nlz87lkmdmk7u+1OtoIvtQIYg0ADPjlJ4ZfHTTcO4e2YNFhaWMfOQ/3PzKfNaVlHsdTwTQPgQRT5TuquKJz1cw6atVOOe4aHA7rj+hE+lJ8V5HkzAX+H0IZmdidvthLvNrzL7DbIr/+WjM7jnIMv/A7MQjzikSJFIaxXL7ad344nfHc+7ANjw/aw3D7/+M+z7IY+tOHZEk3vBuDcEsD/g5zhX4n38NnIlzmw+wTDvgKZw7RWsIEk5WFu/gwY+X887C9TSOi+GyY9tz5bAOpCbGeR1NwsyRryGYtccsD7PJmC3DbApmP8fsK8yWYzbIP99lmD3qfzwZs4cx+xqzlZidW8frPgkcBbyP2U2YdQEqvy8Ds6mYXep/fPX3axHOrQGaY5ZxBH8PIkHrqPQmPHJBfz688TiGd0nn0c/yGXLvp9z7fh6bd1R6HU8iRMwhzNMJGAOMB+YCFwJDgTOBO4HRdSyT6Z+nG76zkP+9z1TnrsFsBHACzm3G7HLgm1pzXAV8hdkq4BbgmFrTvgGGHEJukZDTpWUSj100gF9t3M6jn+YzYcYKJn+9irFHt+XKYR1o0zTR64gSxg6lEFbh3CIAzHKBT3DOYbYIaP8jy7yFczXAEsxaHsJ7ZAI/XL7UuU3+/QmfAWfhXO3LSBYBrQ7hNUVCVreMZB69cAA3Fe/g8c9W8MKsNTw/aw1n9m3F1cOPoltGstcRJQwdSiHUXl+tqfW85gDL117GDuE9dgEp+431Brbw31/+Cf75RcJex/QmPHBeX24+pQvPfLmKl+eu5c1vCxnWOY0rhnZgeJd0zA7lfzGRgwuW8xC+w7dpyse3b+I0oD/wW3wXvNurC747r4lEjNapjbjnFz34+vYTuXVEV5ZtKuOyZ+fy839+wXMzV7OjstrriBIGgqUQZgD9MTPM4oGngPE4tx7fPoRJ/mmx+IpDhxdJREpNjOO64zvx5a0n8uD5fWkSH8M9U3M55n8/4Q9TF7N0o27pGe62le/muw2BuYpu8JyYZvYQ8DbOfXyAec4CBuDc3TrsVMRn/rpt/Ovr1by7cAO799SQ3a4pFwxqy+m9M2kUp1uThwPnHLNXlfDynLW8t3gjHdOb8N6vhx7R5sIDHXYaTIXQEhiMc9MOMM8YYDrObVMhiOyrZOduXp9XwItz1rJq806S4mMY2bcV52W3oV9WqvY1hKDCbbt485sCXv+m0PdvmhDDWf1bM/botvRodWQHFoRGIRwmFYJI3WpqHHNWl/BqzjreW7SBiqoajkprzOj+rTmrf2uymunQ1WC2vaKKDxZvZOr8Qr5esQXnYHCHZozJzuKMeljrUyGIRKjtFVV8sGgjb3xbwKyVvqO3+2WlMrJPJmf0ySQzpZHHCQWgfHc1n+YV8e7CDXySV8Tu6hraNktkdP/WnDugDW2b11+JqxBEhMJtu5g2fz3vLFxPrv/WngPapjKiVwYjembW65eOHFzprio+X1rEh7kb+TSviIqqGtKaxHNG7wxG9W9N/wBt5lMhiMg+Vhbv4N2FG/hwyUYWF/rKoWvLJE7q3oKTurekX1Yq0VHa51Df1mzZyad5RXyaV8TMFVuornGkNYlnRK+WnNG7FYM6NAv437sKQUR+1LqScj7M3cgn3xUxZ3UJe2ocqYmxDO2UxvAu6QztnKZNS0doZ2U1s1Zu4cvlm5mxvJiVxTsBOCq9MSd3b8kpPTPon5VKVAOWrwpBRA5JaXkVXywvZsayYr5YVkxxme+iAx3SGnNsx+YMPqo5g9o3IyMlweOkwal8dzXfrt3GzBVbmLlyCwvWbaO6xpEQG8WgDs05oWs6J3ZrQbvmjT3LqEIQkcPmnOO7DWV8vWIzM1dsYfaqku/PiG7bLJGB7ZoyoG0q/ds2pWtGErHRwXKea8NwzrGuZBffrtvKt2u3MW/NVpZs2M6eGkd0lNG7dQrHHNWcYZ3TGNiuKQmxwXFOiApBRH6y6j01fLehjNmrtjB3dQnfrN32/RpEXEwU3TOS6Nk6hZ6tkumWkUSXlkkkJcR6nLp+VO2pYdXmneRtLCN3fSlL1m9ncWEpW8urAEiIjaJfVirZ7ZoxsH1Tsts1DdrP7kkhmNkkYCRQ5Jzr5R9rBryC7yqpq4HznHNbzbcr/SHgdKAcuMw5901dr7uXCkHEW845Crft4tu121hUWMqiglIWry+lrOKH6yq1Tm1ExxZN6JTehA7pjWnXLJF2zRNpldoo6NYonHMUl1Wybms5qzeXs2rzTlZu3kF+0Q5Wbd5J1R7fd2VstNGlZRK9WqXQNyuVvlkpdGkZOmtIXhXCccAO4LlahXAfUOKcu9d8t9xs6py7zcxOB36FrxAGAw855wYf6PVVCCLBxzlHwdZdLN1YRt7G7Swv8n2hrizeya6qPd/PF2XQMjmBVqmNyEhJoEVSPOlJ8aQ3iadZ4zhSE+NomhhLk4QYkuJjSYiNOqJDMHdX11C+u5rtu6rZtms3pbuqKNm5m807drN5RyWbtlewsbSCjdsrKNy6i8rqmu+XjY4y2jZL5Ki0xnTJSKJLyyZ0aZlE5xZJxMWExpd/XTzbZGRm7YF3ahXCUuB459wGM8sEPnfOdTWzCf7HL+0/34+9tgpBJHTU1DiKyipZs2Una0rKKSgpp3BbBeu37WLj9gqKyyoPeMXW6CgjISaK+NhoEmKiiI42YqKi2HtwjnNQ4xxVexyV1TXsrt7Drqo93/9WX5eYKKNlcgItk+PJTGlEq9QEspolktU0kaxmibRtlhjSX/w/5kCFcCj3Q6hPLWt9yW8E9t48pzWwrtZ8Bf6xfQrBzK7Cdzc12rZtG9ikIlJvoqKMjJQEMlISGHxU8zrn2VlZzZYdu9lavpuS8t2UlldRVlnNjopqdlRWUVFVQ0XVHiqra9hT46iucdTUODDfTVeizIiNjiIuJor4mCgaxUXTOC6axLgYkhvFkuL/adY4jrQmcaQ0itX1nfbT0IXwPeecM7PDWj1xzk0EJoJvDSEgwUTEE43jY2gcH3NoZ0xfeSXcfDP06HFoL15ZCWecAZs3wx13wPnnw7nnwn33wVFH1b3M7t3w85/Dp59CjGdflQ2qoT/lJjPLrLXJqMg/Xghk1ZqvjX9MRMLBnj0QXY+HXT799OHN/+23vj/nz/f9mZvry/RjZQAQFwcnnQSvvAIXXXREMUNNQ28gmwaM8z8eB0ytNX6p+RwDlB5o/4GIBJHRo2HgQOjZEyZO/GG8SRO45Rbo2xdmzoRnnoEuXWDQIPjlL+GGG/77tf74Rxg3DoYNg3bt4I034NZboXdvGDECqnyHeXL88bB3H2KTJnDXXb73OeYY2LRp39csKoKLL4a5c6FfP1ixAqZMgVGjfNPXrIHOnX1rDzU1vvf+6KMfPtuUKfX2VxXsAlYIZvYSMBPoamYFZnYFcC9wspktB37ufw7wHrASyMd3t7TrApVLROrZpEkwb57vC/rhh2HLFt/4zp0weDAsWOD7TfzPf4ZZs+CrryAv78dfb8UK32aaadN8X+QnnACLFkGjRvDuu/89/86dviJYsACOOw6eemrf6S1a+NYohg3zrSF07OjLMHCgb3q7dnDbbXDttfDAA77NUKec4pvWq5evSCJEwDYZOecu+JFJJ9UxrwOuD1QWEQmghx+GN9/0PV63DpYvh+bNfZuIzjnHNz5nDgwfDs2a+Z6PGQPLltX9eqedBrGxvrWCPXt8awbge7569X/PHxcHI0f6Hg8cCNOnHzzzhg2Qnv7D8yuvhNdegyef/GGzEvg+Q1wclJVBUtLBXzfERcaeEhEJjM8/h48/9m0SSkz0bcqpqPBNS0g4sv0G8fG+P6OifMWw90igqCioruPQ1NrzREfXPc/+GjX6ISdAeTkUFPge79ix75d/ZaXvs0SA8DvIVkQaTmkpNG3qK4O8PN8mobocfTR88QVs3er7wn799YbNub/u3SE//4fnt93m23H8pz/59m/stWULpKX5SicCqBBE5MiNGOH7gu/eHW6/3bctvy6tW8Odd/p2KA8ZAu3bQ0pKg0bdxxln+NZuwFdUc+f+UApxcfDss75pn33mmzdC6OJ2ItIwduzwHRFUXQ1nnQXjx/v+9MKuXb6d1V99deDNWmefDffe6zs6KkyE5dVOzawYWON1jiOQBmz2OkQD02eODAf8zBOgzXGQHA/2BWwfD+u8/PY5G5IXQkU+7K5rejzYldDsMdhygJcJxX/nds659LomhGwhhCozy/mxdg5X+syRQZ859GkfgoiIACoEERHxUyE0vIkHnyXs6DNHBn3mEKd9CCIiAmgNQURE/FQIIiICqBA8ZWa3mJkzszSvswSamd1vZnlmttDM3jSzVK8zBYKZjTCzpWaW779veFgzsywz+8zMlphZrpn9xutMDcXMos3sWzN7x+ss9UWF4BEzywJOAdZ6naWBTAd6Oef6AMuAOzzOU+/MLBp4DDgN6AFcYGaHeEuvkFUN3OKc6wEcA1wfAZ95r98A33kdoj6pELzzIHArEBF79Z1zHznn9l6Gcha+u+KFm0FAvnNupXNuN/AyMMrjTAHlnNvgnPvG/7gM3xdka29TBZ6ZtQHOAA7z1m3BTYXgATMbBRQ65xZ4ncUj44H3vQ4RAK2BdbWeFxABX457mVl7oD8w2+MoDeH/8P1CV+Nxjnql+yEEiJl9DGTUMeku4E58m4vCyoE+s3Nuqn+eu/BtZoic+xJGADNrArwO3Oic2+51nkAys5FAkXNunpkd73GceqVCCBDn3M/rGjez3kAHYIH5burRBvjGzAY55zY2YMR692OfeS8zuwwYCZzkwvMEmEIgq9bzNv6xsGZmsfjKYIpz7g2v8zSAIcCZZnY6kAAkm9kLzrmLPc71k+nENI+Z2Wog2zkXaldMPCxmNgL4JzDcOVfsdZ5AMLMYfDvMT8JXBHOBC51zuZ4GCyDz/VbzL6DEOXejx3EanH8N4bfOuZEeR6kX2ocgDeVRIAmYbmbzzexJrwPVN/9O8xuAD/HtXH01nMvAbwhwCXCi/991vv83ZwlBWkMQERFAawgiIuKnQhAREUCFICIifioEEREBVAgiIuKnQhAREUCFICIifv8P59IUnrxjVYsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "fig, axs = plt.subplots() \n",
    "\n",
    "x = np.linspace(-5, 5, 100)\n",
    "\n",
    "# calculate the y value for each element of the x vector\n",
    "y = 100 + (x - 3)**2 \n",
    "\n",
    "axs.annotate('arg min f(x)',\n",
    "            xy=(0.85, 0), xycoords='axes fraction',\n",
    "            xytext=(-0, 0), textcoords='offset pixels',\n",
    "            horizontalalignment='right',\n",
    "            verticalalignment='bottom',\n",
    "            color=\"r\")\n",
    "\n",
    "axs.annotate('min f(x)',\n",
    "            xy=(0, 0.1), xycoords='axes fraction',\n",
    "            xytext=(-0, 0), textcoords='offset pixels',\n",
    "            horizontalalignment='right',\n",
    "            verticalalignment='bottom',\n",
    "            color=\"r\")\n",
    "\n",
    "axs.plot(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como queremos encontrar el valor de $w$ que hace mínimo $f(x)$, es decir $\\arg \\min f(x)$, los pasos a seguir son:\n",
    "\n",
    "* Empezar en algún punto aleatorio de la curva\n",
    "* Iterativamente, movernos en la dirección de mayor pendiente $\\nabla f(x)$ con pasos con tamaño $\\eta$\n",
    "\n",
    "\\begin{equation}\n",
    "    x := x - \\eta \\cdot \\nabla f(x)\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "    x_t = x_{t-1} - \\eta \\cdot \\nabla f(x)\n",
    "\\end{equation}\n",
    "\n",
    "El parámetro $\\eta$ se denomina **learning rate** y es el parámetro que se encarga de medir el \"step size\" en cada iteración mientras nos movemos hacia el mínimo de la función. \n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "    <b>¿Por qué el valor negativo en la ecuación?</b>\n",
    "</div> \n",
    "\n",
    "\n",
    "**DEMO GRADIENT DESCENT:** https://developers.google.com/machine-learning/crash-course/fitter/graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<>:30: SyntaxWarning: \"is\" with a literal. Did you mean \"==\"?\n",
      "<>:30: SyntaxWarning: \"is\" with a literal. Did you mean \"==\"?\n",
      "/var/folders/zb/r_57byp97tx0mcljd3fd0_640000gn/T/ipykernel_4152/3173065574.py:30: SyntaxWarning: \"is\" with a literal. Did you mean \"==\"?\n",
      "  if plot_method is 'scatter':\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "class GD(): \n",
    "    def __init__(self, eta=0.1, x0=0, max_iter=50, diff_to_stop=0.01, delta=0.01): \n",
    "        self.x0 = x0 #randomly initialize any value\n",
    "        self.delta = delta #used for gradient calculations\n",
    "        self.eta = eta #learning rate \n",
    "        self.diff_to_stop = diff_to_stop #stop the algorithms if steps are smaller than this value\n",
    "        self.max_iter = max_iter #when to stop\n",
    "\n",
    "    def run(self, f, plot_method=None):       \n",
    "        plt.figure(figsize=(8, 4))\n",
    "        xx = np.linspace(-10, 10, 100)\n",
    "        yy = f(xx)\n",
    "        plt.plot(xx, yy)\n",
    "        plt.grid()\n",
    "        \n",
    "        x_new = self.x0\n",
    "        iter_no = 0\n",
    "        while True:\n",
    "            iter_no += 1\n",
    "\n",
    "            #main algorithm\n",
    "            x = x_new #note: x_new is the position of the previous move\n",
    "            grad = (f(x+self.delta) - f(x))/self.delta #evaluate the gradient at f(x)\n",
    "            x_new = x - self.eta*grad #move in the direction of gradient\n",
    "            #end of main algorithm\n",
    "\n",
    "            if plot_method is 'scatter':\n",
    "                plt.scatter(x_new, f(x_new), color='r')\n",
    "            else:\n",
    "                plt.arrow(x, f(x), x_new-x, f(x_new)-f(x), head_width=0.3, head_length=2, color='r')\n",
    "\n",
    "            step_size = np.abs(x_new - x)\n",
    "            if step_size <= self.diff_to_stop:\n",
    "                print('Successfully converged with a step size of {} after {} iterations!'.format(step_size, iter_no))\n",
    "                plt.scatter(x_new, f(x_new), s=200, marker='*')\n",
    "                break\n",
    "            elif iter_no == self.max_iter:\n",
    "                print('Max iterations completed! Convergence cannot be guaranteed. Step size={}'.format(step_size))\n",
    "                plt.scatter(x_new, f(x_new), s=200, marker='*')\n",
    "                break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>Vamos a analizar que sucede cuando modificamos los parámetros</b>\n",
    "</div> \n",
    "\n",
    "* Incrementa el número de iteraciones a 50 Increase the maximum number of iterations in order to 50.\n",
    "* Cambia el learning rate a diferentes valores $\\eta \\in \\{0.05, 0.1, 0.8, 1.1\\}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Successfully converged with a step size of 0.009450635923843244 after 24 iterations!\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAD4CAYAAAA5OEWQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA2xUlEQVR4nO3deVwV1f/H8dfhssoqi4iAAuKC4gruaZK7lpqpWVlqmZX5zSwt20tbbNP2zLKyb6aWZS6pqYVbpSnuGwquuIDmiuxwfn9w+33JQBEuzL2Xz/PxuA/unZk78z7MHT7M3JkzSmuNEEIIIayDg9EBhBBCCPE/UpiFEEIIKyKFWQghhLAiUpiFEEIIKyKFWQghhLAijkYHAPD399dhYWEWm9/ly5dxd3e32PyMJG2xTvbSFntpB0hbrJG9tAMs35aEhIQzWuuA4sZZRWEOCwtj8+bNFpvf6tWr6dy5s8XmZyRpi3Wyl7bYSztA2mKN7KUdYPm2KKWOlDRODmULIYQQVkQKsxBCCGFFpDALIYQQVkQKsxBCCGFFpDALIYQQVuSahVkp9blSKk0ptavIMF+l1Eql1AHzz+rm4Uop9Z5SKkkptUMp1bIiwwshhBD2pjR7zF8CPa8YNhH4RWtdD/jF/BqgF1DP/BgFfGyZmEIIIUTVcM3CrLVeC5y9YnA/YJb5+Sygf5HhX+lCGwAfpVSQhbKWSlJaOnP2ZlNQILezFEIIUX6/J53h58O55OQVVMryVGnux6yUCgOWaK2jza/Pa619zM8VcE5r7aOUWgJM0VqvN4/7BXhSa/2v3kOUUqMo3KsmMDAwZu7cuRZp0B8n8vhkRzYPNnOhbZBV9J9SLunp6Xh4eBgdwyKkLdbHXtoB0hZrZA/t0Frz4h9ZpGfn8/qN7jg6KIvMNy4uLkFrHVvcuHJXLq21Vkpd9+6p1noGMAMgNjZWW6pHlU4Fmp8OLmN5ionHB3fCyWTb57dJzznWyV7aYi/tAGmLNbKHdizdeZIjF7cwsokLXW+Kq5RllrVqpf59iNr8M808/DgQWmS6EPOwSuPgoLitvjOH/8pgfkJKZS5aCCGEHcnLL+DtFYlE1vCgfa3KOwJb1sK8CBhmfj4MWFhk+D3ms7PbAhe01ifLmfG6NQ8w0bK2D++uOkBWbn5lL14IIYQdWLD1OMmnLzO+e30clGUOYZdGaS6XmgP8ATRQSqUope4DpgDdlFIHgK7m1wBLgYNAEvApMLpCUl+DUooJPRpy6mIWX28osZ9wIYQQoljZefm8s+oATUO86dG4ZqUu+5r75lrrO0oY1aWYaTXwcHlDWUK7un50rOfPh/FJ3N4qFE9XJ6MjCSGEsBFzNh7l+PlMptzWBFWJe8tg5z1/TejRgHMZucxcf8joKEIIIWzE5ew8PohPol2EHzdE+lf68u26MDcN8aFn45p8tu4Qf6VnGx1HCCGEDfjy98OcSc9hfI8Glb63DHZemAEe716fjJw8PoxPNjqKEEIIK3fucg7TVyfTNaoGMXWqG5LB7gtzvUBPBsaE8PWGI6ScyzA6jhBCCCv20eokLufkMaFHQ8My2H1hBni0a31QMHXlfqOjCCGEsFLHz2cy648jDGgZQoOanoblqBKFuZaPG8Pbh7Fg63H2nbpodBwhhBBW6B3zztu4bvUNzVElCjPA6M518XBx5M3liUZHEUIIYWUOpF7i+y0p3NO2DsE+boZmqTKF2aeaMw91rssv+9LYdPjKm2UJIYSoyt74ORF3Z0cejos0OkrVKcwAI9qHE+jlwpRl+yjNXbWEEELYv4QjZ1m5J5UHboyguruz0XGqVmF2czYxtkt9Eo6cY+WeVKPjCCGEMJjWmteXJRLg6cK9N4QbHQeoYoUZYHBsCHUD3Hl9+T7y8ivnptdCCCGs08o9qfx5+Cxju9SjmnPl3UHqaqpcYXY0OfBkz4Ykn77MvM3HjI4jhBDCIHn5BUxZvo+IAHeGtAq99hsqSZUrzADdGgXSKqw601Ye4HJ2ntFxhBBCGGDe5mMcPH2ZiT0b4miynnJoPUkqkVKKp3tHcSY9m0/XHTQ6jhBCiEp2OTuPaSsP0DrMl26NAo2O8w9VsjADtKhdnT5Ngpix9iBpl7KMjiOEEKISfbruIGfSs3mqd0NDblRxNVW2MEPhbSFz8gp4d9UBo6MIIYSoJGmXspix9iB9mgTRorYxN6q4mipdmMP83Rnatg5zNx0jKe2S0XGEEEJUgndWHSAnr4AJPRoYHaVYVbowA/znpkjcnExMWbbP6ChCCCEq2IHUS8z98yh3talNmL+70XGKVeULs5+HC6Pj6rJqbxq/J58xOo4QQogK9OrSvbi7ODK2q7E3qriaKl+YAe7tEE6wjxuv/LSXggLpqlMIIezRugOniU88zZi4SHytoOvNkkhhBlydTDzRswG7T1zkh63HjY4jhBDCwvILNK/8tJdQXzeGtQ8zOs5VSWE269usFs1CfXjr50QycqTTESGEsCfzE46x79QlnuzZEFcnk9FxrkoKs5lSiuf6RHHqYhafrj1kdBwhhBAWcjk7j7dW7KdlbR/6NAkyOs41SWEuIjbMl17RNflkbTJpF6XTESGEsAefrD3I6UvZPNOnkdV1JlIcKcxXmNirIbn5Bby1ItHoKEIIIcrp5IVMZqxNpk/TIGLqWF9nIsWRwnyFOn7ujOgQzncJKew6fsHoOEIIIcrh9WX7KNAwsWdDo6OUmhTmYoy5KRLfas5MWrIHreXyKSGEsEVbjp7jx20nuL9jOKG+1YyOU2pSmIvh5erEY93r8+ehsyzbdcroOEIIIa6T1ppJi/cQ4OnCQ50jjY5zXaQwl+D22FAa1vTk1aV7ycrNNzqOEEKI67Bw2wm2HTvPEz0a4OHiaHSc6yKFuQSOJgeev7kRKecymbleLp8SQghbkZGTx5Rl+2gS7M1tLUOMjnPdylWYlVLjlFK7lVK7lFJzlFKuSqlwpdRGpVSSUmqeUsp6+z27hvaR/nRrFMhH8Uly+ZQQQtiIGWsPcupiFs/f0ggHB+u/POpKZS7MSqlg4BEgVmsdDZiAIcDrwDStdSRwDrjPEkGN8nTvKHLyC3jjZ7l8SgghrN2J85lMX1N4eVSrMF+j45RJeQ9lOwJuSilHoBpwErgJmG8ePwvoX85lGCrc3517O4QzPyGFbcfOGx1HCCHEVby6dC9aw1O9bOfyqCup8lwOpJQaC7wCZAIrgLHABvPeMkqpUGCZeY/6yveOAkYBBAYGxsydO7fMOa6Unp6Oh4eHxeaXmaeZuC4TP1fFs21dcajEnmMs3RYjSVusj720A6Qt1qiy25F4Np/X/syif6QT/SMt+y2qpdsSFxeXoLWOLXak1rpMD6A68CsQADgBPwJDgaQi04QCu641r5iYGG1J8fHxFp2f1lrP33xM13lyif5u8zGLz/tqKqItRpG2WB97aYfW0hZrVJntyM3L1z2mrdHtX/tFZ2TnWXz+lm4LsFmXUBPLcyi7K3BIa31aa50L/AB0AHzMh7YBQgC7uI/irS2CaR7qw5Rl+7iUlWt0HCGEEEXM2VR496hn+kTh5mzdd4+6lvIU5qNAW6VUNVXYK3gXYA8QDww0TzMMWFi+iNbBwUHxUt/GnEnP5oNfk4yOI4QQwux8Rg5vr0ikbUThjYhsXZkLs9Z6I4UneW0BdprnNQN4EnhMKZUE+AEzLZDTKjQL9WFwbAif/3aI5NPpRscRQggBTF25n4uZubzYt7FN3D3qWsp1VrbW+gWtdUOtdbTW+m6tdbbW+qDWurXWOlJrPUhrnW2psNZgQo+GuDqamLRY+tEWQgij7T15ka83HOHutnVoWNPL6DgWIT1/XacATxce7VafNftPs2JPqtFxhBCiytJa8/zCXXi7OTGuW32j41iMFOYyGNauDg0CPZm0eA+ZOdKPthBCGGHB1uNsOnyOib0a4lPNZjuZ/BcpzGXgaHJgUr/GHD+fycer5UQwIYSobBezcnl16T6ah/owKCbU6DgWJYW5jNpE+NGveS2mrz3Ikb8uGx1HCCGqlHdWHuCvy9lM6tfYJvvDvhopzOXwdO8onBwULy7aLSeCCSFEJdl36iKz/jjMHa1r0zTEx+g4FieFuRwCvVx5tGt94hNPs2pvmtFxhBDC7mmtef7H3Xi6OjKhewOj41QIKczlNLxDGPVqePDS4t1yIpgQQlSwH7cd58/DZ3miR0Oqu9vPCV9FSWEuJyeTA5P7R5NyLpMP4g8YHUcIIezWhYxcXvlpL81CfRjSyr5O+CpKCrMFtI3wY0CLYGasPUhS2iWj4wghhF16c8U+zl7O4ZX+0XZ3wldRUpgt5Ok+Ubg5mXjuRzkRTAghLG37sfPM3niUe9qFER3sbXScCiWF2UL8PVx4omdD/jj4Fwu3nTA6jhBC2I38As2zP+4iwMOFx7vbTw9fJZHCbEF3tq5Ns1AfXv5pDxcy5daQQghhCV9vOMLO4xd47uZGeLo6GR2nwklhtiAHB8Ur/aM5ezmHt35ONDqOEELYvLRLWbz1cyId6/lzc9Mgo+NUCinMFhYd7M097cL4euMRth49Z3QcIYSwaZMW7yE7v4CX7OSWjqUhhbkCPN69PoGerjz1w05y8wuMjiOEEDYpPjGNJTtOMiYukogAD6PjVBopzBXA09WJF/s2Zt+pS3y+/pDRcYQQwuZk5OTx3I+7qBvgzgM3Rhgdp1JJYa4gPaNr0q1RINNW7efY2Qyj4wghhE15d9UBUs5l8uqtTXBxNBkdp1JJYa5AL/VtjEkpnv1xl1zbLIQQpbTnxEU+W3+I22NDaRPhZ3ScSieFuQLV8nHj8e4NWLP/NEt2nDQ6jhBCWL38As1TC3bi4+bEU70bGh3HEFKYK9iw9mE0CfbmpcV7uJAh1zYLIcTVfL3hCNuPnee5mxvhU80+b1JxLVKYK5jJQfHagCacy8jh1aV7jY4jhBBW6/j5TN5Yvo+O9fzp17yW0XEMI4W5EkQHe3N/xwjmbT7G70lnjI4jhBBWR2vNswt2ooFXb21SZa5ZLo4U5kryaNd6hPu7M/GHnXLfZiGEuMKi7SeITzzN+O4NCPWtZnQcQ0lhriSuTiZeG9CEo2czeGfVfqPjCCGE1Th7OYeXFu+heagPw9qHGR3HcFKYK1HbCD/uaB3Kp+sOsjPlgtFxhBDCKkxavJtLWbm8fltTTHZ8n+XSksJcySb2iiq8ReT3O6S7TiFElRe/bAM/bjvBQzfWpUFNT6PjWAUpzJXM282Jyf2j2XvyIp+sSTY6jhBCGObS5WyeXryPyDNHefimSKPjWA0pzAbo0bgmNzcN4r1fktifesnoOEIIYYhX73mBVA9f3nzy1irX7ebVSGE2yEt9G+Pp6siE77aTJ4e0hRBVzG/v/5c59ToyMtyFFvVqGh3HqkhhNoifhwsv9WvM9pQLfCZ3oBJCVCGXU07y5O4cIrLP89jIrkbHsTrlKsxKKR+l1Hyl1D6l1F6lVDullK9SaqVS6oD5Z3VLhbU3fZoE0bNxTaau3E9SWrrRcYQQolK8PmISx71r8MbYXrg6ySHsK5V3j/ldYLnWuiHQDNgLTAR+0VrXA34xvxbFUEoxqX9jqjmbeGL+dvIL5A5UQgj7tvGe//BVzM0MbxFIbJiv0XGsUpkLs1LKG+gEzATQWudorc8D/YBZ5slmAf3LF9G+1fB05YVbGrHl6Hlmrj9odBwhhKgwl+PXMsGtGbUdc5kwoKXRcayWKut9gpVSzYEZwB4K95YTgLHAca21j3kaBZz7+/UV7x8FjAIIDAyMmTt3bplyFCc9PR0PDw+Lza+iaa15b2s2O8/k81J7N4I9/vf/kq215WqkLdbHXtoB0hZr9I92aM1Xa44Tn+XDxNZuNPC1rUPYll4ncXFxCVrr2GJHaq3L9ABigTygjfn1u8Bk4PwV05271rxiYmK0JcXHx1t0fpUh7WKWbv7Sz/qW99fp3Lz8/x9ui20pibTF+thLO7SWtlijou1YF9Zc13lyiZ60cKdxgcrB0usE2KxLqInl+Y45BUjRWm80v54PtARSlVJBAOafaeVYRpUR4OnCy/2bsCPlAh+vlo5HhBD24+KUN3mi1yNEeJqY0CvK6DhWr8yFWWt9CjimlGpgHtSFwsPai4Bh5mHDgIXlSliF9GkaVNjxyK8H2HPiotFxhBCi/A4f5uW1KZzyDODtu9vIWdilUN6zsv8DzFZK7QCaA68CU4BuSqkDQFfza1FKk/tF4+3mzGPfbiMnTzoeEULYtl+7DOLbpt15MC6SFrXl6tnScCzPm7XW2yj8rvlKXcoz36qsurszUwY0YeRXm3ln1X5auxqdSAghyiZ9zwFe6vUIDQKqMbZrPaPj2Azp+csKdW0UyO2xoUxfk8yBc/lGxxFCiOumFy3my1MenPeszrQ7YqQv7OsghdlKPXdLI4KruzFjRzbp2XlGxxFCiNK7dIkFT09jc44Xj/eMolEtL6MT2RQpzFbKw8WRqYObcyZT8/KSPUbHEUKIUksJqcsL3R6kfnUH7u8YYXQcmyOF2Yq1CvOld7gTczcdY+WeVKPjCCHENRWMG8fjfR6jwN2D+5u4YHJQRkeyOVKYrdyt9ZyICvJi4vc7OJOebXQcIYQo2c6dzPztCBtrN+GF/k0IqCYlpizkt2blHB0U7w5pzqXsPCZ8t/3v3tSEEMK65Oezu2s/3ux0D90bBTIoJsToRDZLCrMNqB/oydO9GhKfeJqv/jhidBwhhPiXzAZRPNL3CXx83JlyW1MKb5UgykIKs40Y1j6MuAYBvLJ0L4mnLhkdRwgh/ufLL3m5bjeS/UKZensLfN2djU5k06Qw2wilFG8OaoaXqyOPzNlKVq5c3yyEsAKnT7PilenMbtGbUZ0iuKGev9GJbJ4UZhvi7+HCm4OakZh6iSnL9hkdRwghSI1oyJO9HqFxLS/Gd29w7TeIa5LCbGPiGtRgePswvvz9MPGJcuMuIYRxCu64g/G9HyXTy4d3h7TA2VFKiiXIb9EGTezVkIY1PRn/7XbSLmYZHUcIURX99hufHMplXXhLnr8lmsgaHkYnshtSmG2Qq5OJD+5sQUZOPo/O20Z+gVxCJYSoRDk5JAweyVud7qFPkyDuaB1qdCK7IoXZRkXW8OSlvo35Pfkvpq9JNjqOEKIKueBXg0f6TiDI151XBzSRS6MsTAqzDRsUG0K/5rWYunI/mw+fNTqOEKIK0G++xcSuo0n1CeT9O1vi7eZkdCS7I4XZhimleLl/NCHV3Rg7dxvnM3KMjiSEsGdHjvDN7F9Y1qADE3o2pEXt6kYnsktSmG2cp6sT79/RgrRLWUyYv0O67BRCVAyt2dMqjkld7qdT/QC5a1QFksJsB5qG+PBUryhW7kll5vpDRscRQtihS1178HD/iXhX92Dq4GY4yF2jKowUZjsxokMYPRvXZMqyfSQcOWd0HCGEHdFLfuIpjxYc9a3F+3fG4O/hYnQkuyaF2U4opXhjUFNq+bgx5pstnL0s3zcLISwgPZ2vn/+IJVGdeLxnQ9pE+BmdyO5JYbYjXq5OfHRXS/5Kz+Gxb7dRINc3CyHKaWdkcybfdD9xDQJ4sFNdo+NUCVKY7Ux0sDfP3dKI1Ymn+ViubxZClMOF8RMZ3f8p/L3dmDq4uXyvXEmkMNuhoW1q07dZLd5ekchvSWeMjiOEsEEFu3bx2GEXTvnU4P17WlNdbuVYaaQw2yGlFK8NaEJkDQ/+M2crJ85nGh1JCGFLCgr4aNRkfolszXP9mhBTR65XrkxSmO2Uu4sj04fGkJNXwEOzt5CdJ/dvFkKUztqOfXm741D6Nwnk7rZ1jI5T5UhhtmMRAR68NagZ24+dZ9LiPUbHEULYgJSZs3kk9i4aeDjw6qDm0g+2AaQw27me0TV54MYIZm88yvyEFKPjCCGsWFbqaUavTSPf2YWPH+pMNWdHoyNVSVKYq4AJ3RvQLsKPpxfsZEfKeaPjCCGskNaaZ+99hR1B9Xl7eHvC/d2NjlRlSWGuAhxNDnxwZwsCPFx44L8JnEnPNjqSEMLKfPXQZOY36cYjrQLp3rim0XGqNCnMVYSfhwuf3B3DuYwcRs/eQm5+gdGRhBBW4o/Fa5nk1ZyuLpd59NYYo+NUeeUuzEopk1Jqq1Jqifl1uFJqo1IqSSk1TyklF79Ziehgb16/rSl/HjrLy0vkZDAhBBw/fZGHV6ZQ58Ippj11q3QiYgUsscc8Fthb5PXrwDStdSRwDrjPAssQFtKveTD3dwxn1h9H+HbTMaPjCCGMsHw5ODiQ5eLGAxO+INfkyKc/vIKntwf89JPR6aq8chVmpVQI0Af4zPxaATcB882TzAL6l2cZwvKe7NmQGyL9efbHXSQcOWt0HCFEZUtLQ7u7M6HHf9gdGMG0xW9TNz0Nhg6FTp2MTlflKa3LfqMDpdR84DXAExgPDAc2mPeWUUqFAsu01tHFvHcUMAogMDAwZu7cuWXOcaX09HQ8PDwsNj8jVVRb0nM0kzdkkpmneb6dG/5uFX+6gawX62Mv7QBpy3W5cIFFO87xw+UABrqf5mbvixAeDp6eFl2MrJOSxcXFJWitY4sdqbUu0wO4GfjI/LwzsATwB5KKTBMK7LrWvGJiYrQlxcfHW3R+RqrIthxIvaSjX1iue0xbo9OzcitsOX+T9WJ97KUdWktbrseyERN0nSeX6LF9J+iCUaO0Tk+vkOXIOikZsFmXUBPLs5vUAeirlDoMzKXwEPa7gI9S6u+r0kOA4+VYhqhAkTU8eP+OFuxPvcS4eXKbSCGqgt0ffME4nzY0/+sQU8bdjPrkE3CXa5atSZkLs9b6Ka11iNY6DBgC/Kq1vguIBwaaJxsGLCx3SlFhOjeowbN9GrFiTypvr0w0Oo4QogKdXryc+/eAj6mAGVNH4tpZvk+2RhXxxeKTwGNKqSTAD5hZAcsQFjSiQxhDWoXyYXyydNsphJ3K3LKNkd/t4ay7D58+2p0anq5GRxIlsEhHqFrr1cBq8/ODQGtLzFdUDqUUk/tHc/RsBk/9sINgHzfa1fUzOpYQwkIKUlJ4bPK37Kjfjun3tCI62NvoSOIqpOcvAYCTyYGPh8ZQx8+dB79OIPl0utGRhBCWcOkSrw99jmUNOvBMn0b0kO42rZ4UZvH/vN2c+GJ4KxwdFCO+2MRf0qe2ELYtL485HQfxSduBDG0Vwn0dI4xOJEpBCrP4h1Dfanw6LJbUi1mM+m8CWbn5RkcSQpSF1qyr14pnu4/mxnAfXuzfRO6tbCOkMIt/aVm7OlMHNyfhyDnGzdtGvlxGJYTN2d0ghof6P009/2p8MKw1jib5c28rZE2JYvVpGsSzfaJYtusUk5fs+bvDGCGEDUjp1Z/h3cfh5ePBlw/cgKerk9GRxHWwyFnZwj6N7BjByQtZzFx/iFo+rozqVNfoSEKIazg/ZhzDQnqS7V2d2Q/eQE1vuSzK1khhFlf1TO8oTl3M4tWl+wj0cqVf82CjIwkhSpD13geMPB/EsZBg/ntfO+oHWrbva1E5pDCLq3JwULw9qBmnL2Uz/rvt+Hu40CHS3+hYQogr5C1cxNj4kyTUb8cHd8bQJkL6IrBV8h2zuCZXJxOf3h1L3QAPRn21mR0p542OJIQoQm/axDPTV/Fzg/Y837cxfZoGGR1JlIMUZlEq3tWcmHVva6q7OzP8i00kpUkHJEJYhSNHeGP8h8xr1oP/3BTJiA7hRicS5SSFWZRaoJcrX9/XBgcF98zcyMkLmUZHEqJqu3CBzwY9ysftBnFn61Ae61bf6ETCAqQwi+sS5u/OlyNacykrj7tn/snZyzlGRxKiasrNZf4NA3j5ppH0bhzIZOlAxG5IYRbXLTrYm0+HxXLsbAbDPv+Ti1m5RkcSomrRmmXRnXmi11g6hPkw7Y4WmBykKNsLKcyiTNpG+DF9aAz7Tl3k3i82kZGTZ3QkIaqM+KgOPNJ3Ai1qeTJjRBtcHE1GRxIWJIVZlFlcwxq8O6QFW46eY9RX0q+2EJXh956382Cf8TTwr8bnD9yAu4tc9WpvpDCLcundJIg3BzZjfdIZxnyzhdz8AqMjCWG3tjwwnpGNBlHb24WvxnTG20262rRHUphFud0WE8Lk/tGs2pvG2LlbpTgLUQF2vPoew9xaEVDNkdmPdsHX3dnoSKKCyDEQYRF3t61DTl4Bk5fsQaltvHt7c7mbjRAWsmvWfIamBeLjYuKb8T2o4SX9X9szKczCYu67IZyCAs0rS/fioBTTBjeT4ixEOe1atpa7tuThadLMmXgLwT5uRkcSFUwKs7Co+ztFUKA1ry3bh4OCqYOby2UcQpTRnj93M3T5cTzyc5n7wkBCqlczOpKoBFKYhcU9cGNd8rXmjeWJALw9SPachSiVggLYtg0yM9m1ZDV3rzhBtbxs5gyPIbSabENVhaxpUSFGd47kiZ4NWLjtBI/ICWFClM6MGdC+PQd3HuLOValUy81mzrxnqd2mGaxYYXQ6UUlkj1lUmNGdI3E2OfDyT3vJydvC7aHa6EhCWLedO0nwD+eNcyHUyEzjmwUvEWLKhe++g1tuMTqdqCRSmEWFGtkxAhdHB55buJu00yY6dszH1Ul6KRKiOBt2HuXewZPxcchj3oIXCbrvbnjuOagm3y1XJVKYRYW7u10YTiYHnvphJyO+2MSnw2LxkN6KhPiH+Cen8FCr4YReSOPhes4EbVwH4XILx6pIvmMWlWJI69rc39SFPw+f5a5PN3BO7kolxP9bfPsY7i+Iom7eJeb0C8cnqq4U5SpMCrOoNO1rOfLJ0Bj2nrrE4E/+4NSFLKMjCWEsrfmm7a08EtaDlt6KOa/fhf8tPY1OJQwmhVlUqq6NApk1ojUnzmcycPrvHD5z2ehIQhgjJ4fpbQfxdOeRdA5yZdYTffBylb6vhRRmYYB2df2YM6otl7PzGDj9d3amXDA6khCVquDMGV7uOZopcSO4pZ43n4zpgpuznBQpCklhFoZoGuLDdw+2x8XRxJAZf7B2/2mjIwlRKbJ372XsvW/wWetbGd4qmHdGdMDZUf4Ui/8p86dBKRWqlIpXSu1RSu1WSo01D/dVSq1USh0w/6xuubjCnkTW8OCH0e2p7efOvV9u4octKUZHEqJCXVy+iuEvL2BxoxuZ2KM+LwxoJl3Win8pz79pecDjWutGQFvgYaVUI2Ai8IvWuh7wi/m1EMUK9HJl3gNtaR3uy2Pfbuej1UloLR2RCPtz6sNPuX3eXjaFRjN1cDMejKuHUlKUxb+VuTBrrU9qrbeYn18C9gLBQD9glnmyWUD/cmYUds7L1YkvRrSib7NavLE8kYnf75QuPIVd2TX2GfrtceFoQCifj2zHgJYhRkcSVkxZYu9EKRUGrAWigaNaax/zcAWc+/v1Fe8ZBYwCCAwMjJk7d265c/wtPT0dDw8Pi83PSFWpLQVasyApl8XJuUT5OjCmhSvuTta5R2Ev68Ve2gHW25atW48wPdUXdycY19qdUM9r7w9Za1uul720Ayzflri4uAStdWyxI7XW5XoAHkACMMD8+vwV489dax4xMTHakuLj4y06PyNVxbZ8t/mYjnz6J33TW/H6yJnLFRuqjOxlvdhLO7S2vrYU5OfrmbH9dPiEhfrmSYt16oXMUr/X2tpSVvbSDq0t3xZgsy6hJpbrVECllBPwPTBba/2DeXCqUirIPD4ISCvPMkTVMzAmhK/ubcOZ9Bz6f/QbGw7+ZXQkIa5LTnoGz/R6hEld7qdbbXfmPdmTGl6uRscSNqI8Z2UrYCawV2s9tcioRcAw8/NhwMKyxxNVVbu6fiwY3R6fak4M/WwjX284YnQkIUrlzKEUhj7wPt+06MWDsYF8NPomqjlL3/Ci9Mqzx9wBuBu4SSm1zfzoDUwBuimlDgBdza+FuG4RAR78+HAHOtbz59kfd/H0gp3k5MlJYcJ67V67hX5vrGR7zXq8O6AREwfGyuVQ4rqV+d84rfV6oKRPXJeyzleIorxcnfhsWCve/DmR6WuSSUpN54O7WlDDUw4LCuuy6MufeHJHFj5KMf/ROJqEShcOomykuxlh9UwOiom9GvLukObsOH6em99bz+bDZ42OJaq6o0ehVStyopvy4i1jeWQfNE5LZuGOWTRZv9zodMKGSWEWNqNf82AWjO5ANWcTQ2Zs4PP1h6QzEmGcmTNJPXCEO5rfzZeNu3Pvph+ZM/cZapw5CWFhRqcTNkwKs7ApUUFeLBxzA3ENazBpyR7+M2cr6dl5RscSVY3W/P79L/QZOpW9NcJ5f9k0nk9eidO382DvXmjXzuiEwoZJYRY2x9vNiU+GxvBEzwYs3XmSm99bx67jcocqUTnycnKZ2mc0d/V5Eq/syyz86RVumTAckpNhwACQbjZFOUlhFjbJwUExunMkc0e1Iyu3gAEf/c4Xv8mhbVGxTv22mTuHv817TW9mQOI6FrdU1Nv1JwwdCia5baOwDLm4Tti01uG+LBvbkfHfbeelxXv4PfkvXr+tKb7uzkZHE3bmlwefZrxzNNk1I5naP4oBrXpKMRYVQvaYhc2r7u7MZ8Niee7mRqxOTKPHO2tZnSgdzgnLyDhxiqd7jOE+nw7U9HRh8VM9GdA2QoqyqDCyxyzsglKK+24Ip12EH4/O28rwLzYxrF0dJvaKws1Z/oCKstn23ueM25XH4eY9eKBFAI/dFoOLo3yeRMWSPWZhVxrV8mLRmBu4t0M4s/44ws3vr2PbsfNGxxI2Jiczm3e6jOC2FH+yPb35ZlR7nrq9tRRlUSmkMAu74+pk4vlbGvH1fW3IyMlnwEe/8drSvWTl5hsdTdiAXT+toe/Dn/JOq4H0DXVl2csDaFfXz+hYogqRwizs1g31/Pl5XCdubxXKJ2sP0vu9dSQckR7DRPGy8/J5a+TL9FtzgbPuPnx6Z3Om/ac73m5ORkcTVYwUZmHXvFydeG1AU/57X2uycwsYOP0PXli4i4tZuUZHE1bkz437uPnBT/jAvwX9vbJY+frtdGsabHQsUUVJYRZVQsd6Afw8rhN3t63DVxuO0PXtNSzdeVKue66qMjIgIoJzbl482esRBi9IJsPJjS/mv8jbi9/Gu5rsJQvjSGEWVYaHiyOT+kWzYHQH/D1cGD17C/d+uYmjf2UYHU1UsoLnn2e+Vz263D+d+U268sDG71n51Vjiwn1g6tRrvl+IiiSFWVQ5zUN9WDSmA8/2iWLjobN0nbaGt35O5LL0uV0lbJ8+m9tSgxjf8xHqnDvJkm8m8FQ4VNuWACtXQocORkcUVZwUZlElOZocGNkxgl8f70zv6Jp8EJ9El7fXsHDbcTm8badOJx1lQu9H6XfYh2Pegby56iO+r3OeqK3r4fPPoV49oyMKAUhhFlVcTW9X3hnSgu8fakeApwtj525jwMe/8+chOXvbXmRczuSDIU8Q9+FGfmzcmVFHfie+cQaD1s3H4eXJEBBgdEQh/kF6/hICiKnjy48Pd2B+wjGmrtzP4E/+oGtUIBN7NSCyhqfR8UQZ5OUXMP/Fj5l61pO0sBvp5pbBxNFdqRvQ3+hoQlyVFGYhzEwOittb1aZvs2A+/+0QH69Opvu0tQyKCWXMTZGE+lYzOqIohYICzbKvlzNt7WGS/MNo6ZTGh6Pa0CrC3+hoQpSKFGYhruDmbOLhuEiGtArlg/gkZm84yvdbUhgUW1igg33cjI4oisrIgJkzKQiowbKJb/Jumhv7aoRTVymmt/Kgx4DeKLlHsrAhUpiFKIGfhwsv3NKYUZ0i+Cg+mXmbjjE/4RgDY0Jp7lJgdDwBkJVFXvceLDvvyBv9xnCMRkSYjvHuoje5+eBGTDHTQN1odEohrosUZiGuIcjbjcn9o3moc10+jE/iu80pzM0vYO2FLTx0Y12ig72NjlglZaVnMP/WB5jReBhHqwdRk2ymLX+XvplHMY19BO5aCN6yboTtkcIsRCnV8nHjlVubMLZLPV6cs5a1iaf5acdJOkT6MaJ9OHENa2BykEOmFS3tzEW+mfguX7uGcyZmCM1OJPL0H7Nxfmw4N816E5o1MzqiEOUihVmI61TDy5XBDZyZMqwDszccZdbvhxn51WZCfd0Y1i6MQbGhcuODCrBtWzJfvv5ffgppTq5/Szonb2ZU4k7aDR+AmrGM1Rs2SFEWdkEKsxBl5OXqxEOd6zKyYzgrdqcy6/fDvPzTXt5akUjvJkHcHhtK63BfOfGoLDIz4cUXuXDgEAtVDeZ5N2B3jQg8ajbmrkN/cM/Tw4nwiYKICKOTCmFxUpiFKCcnkwN9mgbRp2kQu45f4Js/j7J42wl+2HKcCH93BsaG0K95sNWezV1QoHGwokPwedt38Md9j/F9UHOW1R9CtpMLjU8lMXnFR9yasgWPe+6CllFGxxSiwkhhFsKCooO9efXWJjzbJ4qlO0/x7aZjvLE8kTeWJ9IqrDp9mwfTO7omfh4uRkcFIL9A027KL/z8aCd8qjkblkNrzdZtySx67TOW1IjmTNdxeGalM3j3r9ye/BvRXdvCe09B+/bgIB0WCvsmhVmIClDN2ZGBMSEMjAnh6F8ZLNp+nIXbTvDcj7t4cdFuWoVVp3ujmnRrFGhoxyUXs/JIvZjDyj2pDIoNrbwFr19P7py5/JnvwYrzJlb6RnLCqwbOtdvQJXkT/Y5spnP7hri+cBe0eU+KsahSpDALUcFq+1VjzE31eDgukn2nLrFkxwlW7E5l0pI9TFqyh6ggL26sH0Cnev7EhFXHxdFUadnOZeQAMD8hpVIK84k1G1j/6sescw1kTXhHLnp64OKeTadDW3hs3dd0378Br2lvwkOvgHw3L6ooKcxCVBKlFFFBXkQFeTGhR0MOnbnMyj2nWLU3jc/WHWT6mmTcnEy0DvelTYQvrcJ8aRLsjatTxRTqS1m5ZOTkAya2Hj3PpaxcPF0tdDZ5QQHEx3Ny3gI2p5vYnOnEb351SfKvDS0GE5B+lu5JG+mWsp2OjYKo1icO+o8DPz/LLF8IG1YhhVkp1RN4FzABn2mtp1TEcoSwZeH+7ozqVJdRneqSnp3HhuS/WJ90hnUHTrNm+WkAnE0ONAnxplmID9HBXkQHexPh746j6foO7V7IyOX4+cx/DPs9+Qx/75M6OSrmbTpG+7r/7E862McN72rXKNYFBZCQwIXdiexOvcyuzXvZle1EQnAUx317gS9Uy8kk5vhehuxYwQ2HttLgxljU6y9Bo0ayZyzEFSxemJVSJuBDoBuQAmxSSi3SWu+x9LKEsBceLo50bRRI10aBAPyVnk3CkXMkHDnHpsNn+ebPI2TlFnYD6urkQGQND+oGFD4iAtyp4+tOLR9XfN2di70866PVSXyy9iAujg44FSnqD9QvvPf05ex83ll1gHdWHQAgN7+A7LwCHugUwVO9o0BrOHuW7IuXSE3YTcr6TRxMPsHBfGeSPQJI8gvluHcgUB3qhhB08TQtTiRy36aFtErdT1SoL46DBsL0z6F69Qr+bQph2ypij7k1kKS1PgiglJoL9AOkMAtRSn4eLnRvXJPujWsCkJebx6G0S+xKOc+ulPMcSEtnc/JpFm478Y/3uShNsCmXgIJsfPMyqZ59Gd/Mi/hmXKCrcwBr/OqRbvrfZp988X9FPD077/+fO+bnEXtyP6c/+pURX3hy1s2bE17+nPbwLZzAtRU0BrecLOqeTSHm+F6Gbl1K49OHaOzjhF/HtvBIX4h8CGrXrsDflBD2R2mtLTtDpQYCPbXWI82v7wbaaK3HXDHdKGAUQGBgYMzcuXMtliE9PR0PDw+Lzc9I0hbrVOq27NoF2dkVliNbK07lO3Mm34mz+U78VeDImXwnLmpH0gtMXCowka5NaK7/cLEzBXg65OPhkI+nyqO6KQ8/h1z8zD8DTTlUdyzAwcsTfH3BywscjTttpUp+vqycvbQDLN+WuLi4BK11bHHjDNuKtNYzgBkAsbGxunPnzhab9+rVq7Hk/IwkbbFOpW7L/v0wbhzk5V172gpSgCLd2Y0sJ2eyHF34q5oXD976DKme/gyNzOPrpP/9GQi6kMa33zxJQMYFXPMKz9gmMBDatIFu3aBDh8LvhfPyCouwi3Vcjw1V9PNl5eylHVC5bamIwnwcKHrdRYh5mBBVz6hRhY+/HTsGubnw8svg7Ay//AJJSRUawQGNV04GXjkZANRIP8v5aoV3XQq8ojOyc76B1EhNweXwwcLCGxFR/DXEVlSQhbA3FVGYNwH1lFLhFBbkIcCdFbAcIWxPqPl/1s8/v/a0SUng5gY1a4LJfMnUpEmFe7AtWkD9+oW3Ndy8Gd55Bzp2hBtugLp1C9/39/TdukHLlv9fTNfvScVp3lays/NxUAonk8LkoMjKLcDR0cT6Q+foEtXQ8m0XQpSKxQuz1jpPKTUG+JnCy6U+11rvtvRyhLB7kZH/Hvb88/8e1qoVzJ5d/DyKmf6HrcdJz87H1ckBX3cT257vyhvL9zFv8zHSs/P4YetxukQFljO8EKKsKqSfO631Uq11fa11Xa31KxWxDCHE9cvJKyB+XxoeLo58cncsQd6uuLs48lK/aKYPjcHDxZFf96aRm19gdFQhqizpgFaIKqRAa26LCSF+fGdurB/wj3GdG9QgfnxnbosJIb/AsldrCCFKT7rkFKIKcXUy8XL/6BLHB3i6XHW8EKLiyR6zEEIIYUWkMAshhBBWRAqzEEIIYUWkMAshhBBWRAqzEEIIYUUsfhOLMoVQ6jRwxIKz9AfOWHB+RpK2WCd7aYu9tAOkLdbIXtoBlm9LHa11QHEjrKIwW5pSanNJd+2wNdIW62QvbbGXdoC0xRrZSzugctsih7KFEEIIKyKFWQghhLAi9lqYZxgdwIKkLdbJXtpiL+0AaYs1spd2QCW2xS6/YxZCCCFslb3uMQshhBA2SQqzEEIIYUVstjArpQYppXYrpQqUUrFXjHtKKZWklEpUSvUo4f3hSqmN5unmKaWcKyf51ZmzbDM/DiultpUw3WGl1E7zdJsrOWapKKVeVEodL9Ke3iVM19O8rpKUUhMrO+e1KKXeVErtU0rtUEotUEr5lDCd1a6Ta/2OlVIu5s9eknm7CDMg5jUppUKVUvFKqT3m7X9sMdN0VkpdKPK5e96IrKVxrc+MKvSeeb3sUEq1NCLn1SilGhT5XW9TSl1USj16xTRWu06UUp8rpdKUUruKDPNVSq1USh0w/6xewnuHmac5oJQaZrFQWmubfABRQANgNRBbZHgjYDvgAoQDyYCpmPd/CwwxP58OPGR0m4rJ+DbwfAnjDgP+Rme8Rv4XgfHXmMZkXkcRgLN53TUyOvsVGbsDjubnrwOv29I6Kc3vGBgNTDc/HwLMMzp3CW0JAlqan3sC+4tpS2dgidFZS9meq35mgN7AMkABbYGNRme+RntMwCkKO8+wiXUCdAJaAruKDHsDmGh+PrG4bR7wBQ6af1Y3P69uiUw2u8estd6rtU4sZlQ/YK7WOltrfQhIAloXnUAppYCbgPnmQbOA/hUY97qZMw4G5hidpYK1BpK01ge11jnAXArXodXQWq/QWueZX24AQozMUwal+R33o3A7gMLtoov5M2hVtNYntdZbzM8vAXuBYGNTVah+wFe60AbARykVZHSoq+gCJGutLdmTY4XSWq8Fzl4xuOj2UFJ96AGs1Fqf1VqfA1YCPS2RyWYL81UEA8eKvE7h3xuuH3C+yB/b4qYxWkcgVWt9oITxGlihlEpQSo2qxFzXa4z5ENznJRwOKs36sib3UrgHUxxrXSel+R3//zTm7eIChduJ1TIfbm8BbCxmdDul1Hal1DKlVOPKTXZdrvWZsbXtYwgl70zYyjoBCNRanzQ/PwUEFjNNha0bR0vMpKIopVYBNYsZ9YzWemFl57GUUrbrDq6+t3yD1vq4UqoGsFIptc/8n1+lulpbgI+ByRT+8ZlM4aH5eysvXemVZp0opZ4B8oDZJczGKtZJVaCU8gC+Bx7VWl+8YvQWCg+lppvPa/gRqFfJEUvLbj4z5vN0+gJPFTPaltbJP2ittVKqUq8rturCrLXuWoa3HQdCi7wOMQ8r6i8KDwk5mvcOipumwlyrXUopR2AAEHOVeRw3/0xTSi2g8HBlpW/QpV1HSqlPgSXFjCrN+qpwpVgnw4GbgS7a/AVTMfOwinVSjNL8jv+eJsX8+fOmcDuxOkopJwqL8myt9Q9Xji9aqLXWS5VSHyml/LXWVnczhVJ8Zqxi+yilXsAWrXXqlSNsaZ2YpSqlgrTWJ81fHaQVM81xCr87/1sIhec8lZs9HspeBAwxn2UaTuF/ZX8WncD8hzUeGGgeNAywpj3wrsA+rXVKcSOVUu5KKc+/n1N4ctKu4qY10hXfhd1K8Rk3AfVU4VnyzhQeCltUGflKSynVE3gC6Ku1zihhGmteJ6X5HS+icDuAwu3i15L+ATGS+XvvmcBerfXUEqap+ff340qp1hT+nbO6fzJK+ZlZBNxjPju7LXChyCFWa1PiUT5bWSdFFN0eSqoPPwPdlVLVzV/TdTcPKz+jz4gr64PCP/QpQDaQCvxcZNwzFJ6Fmgj0KjJ8KVDL/DyCwoKdBHwHuBjdpiI5vwQevGJYLWBpkezbzY/dFB5uNTx3Me34L7AT2EHhBz3oyraYX/em8OzaZGtsi/kzcgzYZn78ffayzayT4n7HwCQK/9kAcDVvB0nm7SLC6MwltOMGCr8a2VFkffQGHvx7mwHGmNfBdgpP1mtvdO4S2lLsZ+aKtijgQ/N620mRK1Cs6QG4U1hovYsMs4l1QuE/EyeBXHNNuY/C8yt+AQ4AqwBf87SxwGdF3nuveZtJAkZYKpN0ySmEEEJYEXs8lC2EEELYLCnMQgghhBWRwiyEEEJYESnMQgghhBWRwiyEEEJYESnMQgghhBWRwiyEEEJYkf8DYQl6GEN1Hm8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def f(x, p=1, q=0, r=0):\n",
    "    return p*(x-q)**2 + r\n",
    "\n",
    "GD(eta=0.1, x0=8, max_iter=100, diff_to_stop=0.01, delta=0.01).run(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "¿Qué sucede en este caso? \n",
    "$$f(x) = -\\exp(-(x-2)^2) - 0.5 \\exp(-(x+2)^2)$$\n",
    "\n",
    "Sugerencia: max_iter en $100$, x0 en $-3$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Successfully converged with a step size of 2.9556418361664782e-06 after 19 iterations!\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAD4CAYAAAAuGtVZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0JUlEQVR4nO3deXxc5X33/c9Po12yrdWyLVmSbbyAnbDY2GDA2IEESJNCEpKQuAktIc5G29x9eBpSt7nbpH6aNm1umju0DSVNIXHq0GzQhITVjiGAjU3AK14l2xLetFmrtc31/DEjIxvJkjzLOWf0fb9eemnmnEvnXNecmfnpWo855xARERF/S/M6AyIiIjIyBWwREZEAUMAWEREJAAVsERGRAFDAFhERCYB0rzMwnJKSElddXR3XY3Z0dJCXlxfXY3ohVcoBKotfpUpZUqUcoLL4VbzLsnXr1gbnXOlQ+3wbsKurq9myZUtcj7lhwwaWL18e12N6IVXKASqLX6VKWVKlHKCy+FW8y2Jmh4bbpyZxEZFxLBzWWhxBoYAtIjJONXf0cPXXn6Wls8frrMgoKGCLiIxTT+8+zvHWbp7eddzrrMgoKGCLiIxTP95ad9Zv8TcFbBGRcajtdC+vHW4B4HeHW2g73etthmREcRklbmY3A/8MhICHnHNfP2d/FvAIsBBoBD7qnKuNx7lFROT8TnX2Ut/Sdda2F//hO2RYNT1ZuWR0tvOj++5n6Z9/5qw05QU5TMrNSGZW5TxiDthmFgIeAN4N1AGvmNnjzrldg5J9Cmh2zl1kZncAfw98NNZzi4jIyP5lw36+s/EgWelpZITS4HQXWBUdWbkAdGTlcn93Fff/83OQnUNvf5juvjCfWTaTL7/3Yo9zLwPiUcNeDOx3zh0EMLN1wK3A4IB9K/DX0cc/Br5tZuaSdG/PcNhx/zN7qT3Uw6s9e5JxyoQKZDnMzn4a3XSotodt/ftIMzAzQmlGyIyMkJER/XLJSk9jYnYGE3PSmZidwaTcDEryskhLs6HPJRIQXT391Ld00tzZS3t3H53d/XR099Hd109f2NHX7+gLO8LOceBg5LMSdo6Bb85hv0DP+WrNDKVx9cxiNtc20t0XBsuArLNrzu1Z0cU/uvsIpcHVM4vJDKXxzafi/10TyO+wYUzpCSftXBZrzDSz24GbnXN3R59/AljinLtnUJod0TR10ecHomkazjnWKmAVQFlZ2cJ169bFlLcBYef41JOdRN7eqfAlH6xyJOK/sswQlOWmUZZrTM1L45LiEHMK0wilGTQ1QX099PRAZiaUl0NRUQJycbb29nby8/MTfp5kSJWy+KkcHb2O7Q39bD/Zz9GOMA1dYVoTOJtquG8IN4pPpCX8+yVY32Hn8+mLHUur4vceW7FixVbn3KKh9vlqpTPn3IPAgwCLFi1y8Vw9pmZF6qyukwrlcNFawobfbGDZsutxQH84sq0vHKav39HbH6anP8zp3jDt3X20dvXSerqXpo4eDjV2UtPQQU1DB7+r7eR/DvZSmJvBjVnt3PTI/Vz/xotkhPsjJ8vNhQcfhJUrE1qmVLguA1KlLF6Xo7uvnx9uOsyvdxxjy6Fm+sOOwtwMFpQXsaQwh4rCXMoLcijOzyQvK528zHTyskJkpYfICEVanDJCaZjB8xs3smL58jOtU2YXFvDaTvfy3i8+wpGCKW/bN73lGE/c/0kmZCe239rr6xJPySxLPAJ2PTB90POK6Lah0tSZWTowicjgMxmnzAwzSDMjPRSZrJARGtgbGvbvhtLR3cfGvSd5cucxfr35FP/9/i8xb0kNa556gIX1b0BnJ6xenfCAHVhr18JnPgMdHZHnaWnwve95m6cUsLmmift+uo2DJzu4eOpEPnf9LFbMm8xl0wsiLUFjlJ5mF/R358oIpXFi0pBLVXNyUimZ6Zo85FfxuDKvALPNbIaZZQJ3AI+fk+Zx4M7o49uB55LVfy2pLy8rnVveMZX777icrd9ayQM//ztas/P40B/8I1++6QucysqDw4e9zqb/rF0LJSXwB3/wVrAGCIfh5EmYMCGSRsbkVFcvX/7pdj7ynZfo6Qvz8F2L+dWfXse9N81lYVVhXIJuLF7Y10BGdiYAOb2nyejrJbv3NADp2Vm8sK/hfH8uHoo5YDvn+oB7gCeB3cCjzrmdZvZVM/v9aLLvAsVmth/4M+C+WM8rMpTMimn83p7f8vRDn+fuzT/j0Xe+hxs+/W9sXPwer7PmL2vXwqpV0Bhp6GrOnsATc6/hyMTJb6Vpb4+kUdAetYMn23n3N3/Do1uO8JllM3nqfy3j+jlD12a98tPf1dPe3U92RhofvX4ur615Hx9dNpfsjDTau/v46e/ObSAVv4hLH7Zz7gngiXO2fWXQ49PAh+NxLpHzWrMGVq0ir7OTv1z/XT6w8zn+n/ffy6oV9/DDw81cUVnodQ69t3Yt3Hknpy3EM/Ou5eeXrGDDzIX0hSJfB4sPb2d+1ySuyMxloroTRq3tdC+rvr+VvrDjsS9cw4LySV5n6W16+sKsf+ME+VnpPLDyijP/TPzNrQtYMW8y9/zwdzy3+wS9/eHI9C/xFV0RSS0rV0YGmFVVgRnzc8L84LoCyoryuOs/X2H/iXavc+itaM36xfL5LP7CI9xz631sn3IRd215jHU/vI97Nz5CQ14h32ufypX3fJ+fXbJc3QmjEA47/uzR16lp6ODbH7/cl8EaIjNmPrSwgvX3Ln9bzX/53Mmsv3c5H1pYQb/u4OVLCtiSelauhNraSF9sbS0ld36cR+5aTHqacef/Xc/xeZdGBlZVV4+/5t7Vq9k2cRqf/uBfMqWtkR/+11/w4r/+EX+x4XtcdWQH97z0KM8+9Fm+UlDLpcf28aVb/pRX514Zea3G62s2Ct96bh9P7zrOX/7exSydVeJ1doaVnRHib29bQOmErCH3l07I4m9vW0B2xtgGfkpyKGDLuFBVnMd/TmmkpaObO5d+mlOZuXDo0Ljro93fHuYPP/w3FHa18v1H/4qlh7cRcmcv/GDAzIzTfOena5jS3shnb7iH441tkcU4xuFrNpKndh7j/mf28aErKvjDpdVeZ0dSmAK2jBsL/vY+/u1nazhQXMGXbvmTyMaBPtpUtnYtVFfz5qTJfPIjXyXNhfnBj/6Ksvamt9KEQvCDH0R+qqoAKCwr4t+f/RbtmTms+sBqToeic3PHw2s2SkeaOvmzR1/n0opJrPnAggueGy0yGgrYMn4cPsx1ta/xJ79dx6/nXsPmivlntqesaJ91y/FGPvGRr9KWlcvDj36F6pajb6XJzYWHH450JQx0JyxcCLW1zD2wnW/+8pu8Pm0uq2/6wltrZKXyazYGD248SE9fmAdWXqFmZEk4BWwZPyorAbj7lZ8zpa2BNSs+RRg7sz0lrV4NnZ3cf83HqS2cxkM//irzT9REatRmkdr0+VaBq6zk5r0v8Se//SE/eceNPPrOd5/ZPt41d/Tw31uPcOtl06gozPU6OzIOKGDL+LFmDeTmktPXzb0bH+H1aXP4xaU3RLanqsOHOTJxMmsvv4WPbHuKJXU7I9vD4TOD8s47ZSv6mn3xhf/i0jf38K9LPkw4Ny+1X7NRWrvpEKd7w9x93UyvsyLjhAK2jB+Dpnx9YNdvuKT5CP/w/j+m+6N3eJ2zxKms5J+v/RjmHH/y4rqzto9K9DVLq6rkj179H2qLpvH8P3133M/L7u7r5+GXDrFsTilzp0zwOjsyTihgy/gS7aMN9fex+v+9nbreEI+8eMjrXCXM/r/6O346/wbufPUXTG2LLt+fmzu2GnL0NbvltWcozsvk+zmqUT722pucbOvm09fN8DorMo4oYMu4dc1FJSyfW8r/fW4fzR0JvM+hh/4paw456cbn6l4eXZ/1eWSlh7hj8XSefeMER5o6E5DbYHDO8d3na5g3ZQLXXuTfOdeSehSwZVz78i0X097dx7fX7/c6K3G3ra6FX+04xt0r5lC0d+fo+qxH8PElVRiwdtP4HSW+cV8De4638enrZmoalySVAraMa3OnTODWy8r50StH6Orp9zo7cfWPT+2lMDeDu+PYbFtekMONF5fx6JYjnO5NrddrtB56/iBlE7N4/6XTvM6KjDMK2DLuffTK6bR39/HrnUdHThwQmw42snHvST6//CImZGfE9difvLqapo4entieOq/XaL1xrJXn9zVw59Jq3Tdakk7vOBn3FlcXUVmUy39vqfM6K3GzdtNhCnIz+MTVVXE/9jUXFTOzNI9HXkrdwXrD+eW2o6QZfOxKzUOX5FPAlnEvLc340BUVvHigkbrmgA+mWruW9ovm8tSWGn7v9efIfnTdyH8zRmbGJ66q4rUjLWyvOxX34/vZxr0nubyykMK8TK+zIuOQArYI8KGF5ZjBT7bWe52VCxddhvSp7HJOZ2TxgU2PJ+xGHR9aWEFuZogfbh4/g8+aOnrYVn+KZbNLR04skgAK2CJARWEuS2cV8+NXjxAO6r2Ao8uQ/vyS5VS0HGNh/e6E3ahjYnYG180u4fl9J+N+bL96YX8DzsGyOZrKJd5QwBaJun1hBUeauthc2zRyYj86fJgTeQW8UH0Zt+3agA3anghXzyymrrlr3MzJ3rj3JAW5GbyzosDrrMg4pYAtEnXz/KlMyEoP7uCzykp+MW8Z4bQQt+3acNb2RFgaXTTkpQONCTm+nzjneH7fSa65qIRQmuZeizcUsEWicjJDvO/SqTyx/Sjt3X1eZ2ds1q6F9nZ+Pn85C47t56LG6D8dY12GdAxmT86nJD+TFw80JOT4frLneBvHW7u5Xv3X4iEFbJFBbl9YQVdvP09sC9Ac4+hgswMum21T57xVuy4uvuBlSEfDzLhqZjEvHmjEuYD2+4/Sxr2Rvvrr1H8tHlLAFhnkispCZpbm8eOtAWoWjw42e+yS5aSF+3n/7o2R7fn5Cb+r1tJZJZxo6+ZgQ0dCz+O1jXsbmFOWz9RJOV5nRcYxBWyRQcyM971zGq8cagrODUEOH8YBP79kOUsPb6OsvenM9kRbOqsYgBdTuB+7q6efzbVNms4lnlPAFjnHirmlOAcbgzJlqbKSV6fN43DhVG7buf6s7YlWVZzL1EnZvJTC/dgv1zTS0xdm2RwFbPGWArbIOd5ZUUBhbgYb9gQkYK9Zwy/esYKs3m5u2vtSZFsCB5sNZmZcPauYlw82BXf++gg27j1JVnoai2cUeZ0VGecUsEXOEUozrp9Tym/2ngxGEFq5ko1LbmZJwwEm9J6O6Z7XF2LprBKaOnrYc7wtKedLto17T7J4RhHZGSGvsyLjnAK2yBBWzJt8ZilKv6tv6eJAd4hld38oLve8HqurU7gfu76liwMnO7hezeHiAwrYIkO4bnYpZrBhzwmvszKi56NTjrwKKuUFOVQV56bkAiovRMcxqP9a/EABW2QIRXmZXDa9gPUB6MfeuO8kUyZmc9HkfM/ysHRWMZsONtLXH/YsD4nw2pFTTMrJYLaHr63IAAVskWEsnzOZbXUtNLZ3e52VYfX1h3lhXwPL5pRg5t2SmVfNLKatu4+db7Z6lodE2PXmKeZPm+jpaysyQAFbZBgr5vl/ete2+lO0nu7jOo/nCA/0Y790MHWaxfv6w7xxrI1Lpk70OisiQIwB28yKzOxpM9sX/V04RJrLzOwlM9tpZtvM7KOxnFMkWRZMm0RJfqavp3dt3HsSM7j2Im+XzJw8IdIk/3IKBeyahg66+8JcMk0BW/wh1hr2fcCzzrnZwLPR5+fqBD7pnJsP3Azcb2YFMZ5XJOHS0oxl0eld/T6d3rVx78nIvPG8TK+zwjsrJrErhZrEB5r350+b5HFORCJiDdi3Ag9HHz8M3HZuAufcXufcvujjN4ETgIZcSiAsnzuZls5eXq9r8Torb3Oqq5fXjrSwbLY/bkhx8ZSJnGjrpikoS7qOYNfRVjLT05hZmud1VkQAsFjusmNmLc65guhjA5oHng+TfjGRwD7fOfe24aRmtgpYBVBWVrZw3bp1F5y3obS3t5OfH/zRnqlSDvB/Wdp7HH/8XCfvn5XBB2efvxab7LK8cqyPB17r5i+WZDOnML6LelxIWXY09PGPW7r50pXZXFzsj0VGYrkm33ili85e+N9L/XHDD79/VsZCZRneihUrtjrnFg250zl33h/gGWDHED+3Ai3npG0+z3GmAnuAq0Y6p3OOhQsXunhbv3593I/phVQph3PBKMsH/+W37n3fen7EdMkuy30/ed3N/8qvXU9ff9yPfSFlOd7a5aq+9Av33ecPxj0/F+pCr0k4HHaX/c2T7ks/fj2+GYpBED4ro6WyDA/Y4oaJi+kjRXvn3I3D7TOz42Y21Tl31MymEmnuHirdROCXwGrn3MsjnVPET66fU8r/eWYvLZ09FOR631cMkX+0N+5tYOmsYjJC/pjsUZqfRVFeJnuOBX+J0mOtp2nu7NWAM/GVWD/pjwN3Rh/fCTx2bgIzywR+BjzinPtxjOcTSbolM4pwDjbXNHmdlTMONnRQ39LlqxW4zIx5UybwRgqsKb6zfmDAmQK2+EesAfvrwLvNbB9wY/Q5ZrbIzB6KpvkIsAz4QzN7LfpzWYznFUmaS6cXkJmexiYfBeyN0eVI/XaP5rlTJrD3WFswbppyHruOtmIGc6coYIt/jNgkfj7OuUbghiG2bwHujj7+AfCDWM4j4qXsjBCXTy/wVQ37t/sbqSrOpbI41+usnOXiKRPp6u3ncFMn1SXBHV29681WqovzyM+K6StSJK780fkl4nNLZhaz881TtJ7u9TorhMOOV2qbuGpGsddZeZu5UyYA8MaxYM/H3nn0lPqvxXcUsEVG4aoZRYQdbK1t9jor7DnexqmuXhbPKPI6K28zp2wCZvBGgAeenerq5UhTl5YkFd9RwBYZhcsrC8kIGS/XeL/05kDTvB8Ddk5miOriPN44GtyAvftopHVANWzxGwVskVHIyQxxaUUBmw5634+9uaaJaZOyqSj0x4Ie55o3ZQJ7AjxSfNebGiEu/qSALTJKS2YWsb3+FB3dfZ7lwTnHppomFs8o8u0tH+dOmUBtYwedPd69TrHYdbSVkvwsJk/I9jorImdRwBYZpSUziukPO7Ye8q4fu6ahg4b2bhb7cMDZgHlTJuIc7Dve7nVWLsjON1vVHC6+pIAtMkpXVBUSSjM2ediP7ef+6wHzoiPFg7jiWU9fmP0n2tQcLr6kgC0ySvlZ6Swon+RpP/bmmiaK8zKZ5eM7SFUW5ZKTEWJ3AKd27TvRRm+/0whx8SUFbJExuGpGEa/XtdDV0+/J+f3efw2R+4jPmTIhkDXsnRpwJj6mgC0yBktmFtHb7/jd4eT3Y9c1d1Lf0sUSHzeHD5hXNoE3jrUN3KkvMHYfbSUnI0RVsX9bMGT8UsAWGYNFm58hLRzm5T/6IlRXw9q1STv3K7UD/df+HXA2YN7UCTR19HCyvdvrrIxJbUMHM0ryCKX5twVDxi8FbJHRWruWiZ9bxSUnDrK5Yj4cOgSrViUtaG+uaWJidvqZ5T/97MwSpQFbQKW2sZPqEn+tzy4yQAFbZLRWr4bOTpYc2cGr5fPoDqVDZ2dkexJsqmniyuqiQNT+5kXvchWkfuy+/jBHmjqpVnO4+JQCtshoHT4MwJLD2+lJz2TblDlnbU+kk23dHDzZ4evpXIMV5WUyeUJWoEaK1zV30Rd2gb7LmKQ2BWyR0aqsBODKul0AbKpccNb2RHqr/zoYARtg3tSJgaph1zZ2ADBDAVt8SgFbZLTWrIHcXApPtzH3ZC2bpi+A3NzI9gTbXNNETkaIBeWTEn6ueJlZkkdtQ0dgRorXNkQCdpXP7jEuMkABW2S0Vq6EBx+EqioWH9nJ1or59H3nwcj2BHv5YCNXVBWQEQrOR7aqOJeOnn4aO3q8zsqo1DZ2kpcZojQ/y+usiAwpOJ9+ET9YuRJqa1n8jb+kMyOLndf/XsJP2dLZw57jbVwVgOlcgw0M3joUbWr2u5qGDqpL8ny9KI2MbwrYIhdgye82ALD59rsSPh97c00TzsGSmcEK2JXRpuVDjZ0e52R0DjV2aIS4+JoCtshYrV3L5C98mhlN9WyqWPDWfOymxKwx/vLBJrLS07h0enD6rwEqCnMwC0bA7u0Pc6S5S3OwxdcUsEXGKjofe/GRHWyePp8wFpmPXV+fkNNtqmnkispCstJDCTl+omSlh5g2KYfDTf4P2HXNXfSHnWrY4msK2CJjFZ13vfjITlqz89lTWhXZ3hP/wVWnunrZdbSVJTODM51rsKri3DPTpfxMU7okCBSwRcYqOu968ZEdAGyeHp2PnZkZ91O9Eu2/vipg/dcDqopzORyAJvG3pnQpYIt/KWCLjFV0PnZF6wmmtZ5g8/T5kfnY5eVxP9WmmkYy09O4bHpB3I+dDJVFeTR29NB2utfrrJxXbUMH+VnplOTH/58ukXhRwBYZq+h8bKuqYsmRnWyqvhT3nQehKP7N1ptqmrhsegHZGcHqvx5QHZCR4jXRm35oSpf4mQK2yIUYmI/9V39CQ/ZEDt50a9xP0Xq6lx31pwLbHA5vTe3y+8CzQ40dag4X31PAFonBwNrem2viP6Vra20zYQdXBWj98HNVnVk8xb8Bu7c/TF1zFzMUsMXnFLBFYjCzJI+S/Mz4B+y1a3n5S2vI6O/l8puuTto9t+MtPyud4rxMDjf5d6T4kabOyJQujRAXn1PAFomBmbF4RlF8A/batbBqFS9PquKyN/eSc3B/ZGGWgAbtquJcahv8W8N+a0qXFk0Rf1PAFonR4uoi6lu6ONkZjs8BV6+mvc+xY8pFLDmyPbKtszOyYEsAVRXn+boPe+CfCfVhi98pYIvE6NrZpQBsb+iPzwEPH2Zr+cX0p4VYEp3rPbA9iCqLcnnzVBfdfXF6feKstrGDCdGmexE/iylgm1mRmT1tZvuivwvPk3aimdWZ2bdjOaeI38wqzaOyKJfXT8YpIFVW8vL0d5De38fC+t1nbQ+i6pJcnIss/+lHukuXBEWsNez7gGedc7OBZ6PPh/M1YGOM5xPxHTNjxdxSdjf2c7o3DkF7zRperr6Udx7bR25vd2Rbbm5kwZYAqizy9202DzV2UlWs/mvxv1gD9q3Aw9HHDwO3DZXIzBYCZcBTMZ5PxJdWzJtMTxhePtgY87Eabr2d16fO5trmGjCDqip48MHI3O8AqvLx4ik9fWHqmju1hrgEgjnnLvyPzVqccwXRxwY0DzwflCYNeA74A+BGYJFz7p5hjrcKWAVQVla2cN26dRect6G0t7eTn58f12N6IVXKAalTlp5+xxee7WBZRQafuCQrpmNtrOvlP3b08DdLs6ma6M0KZ/G8Ls45PvdMJ9dVpLPy4them7EaqRxH28N8+YUuPv2OTK4pz0hizsYuVT4roLKcz4oVK7Y65xYNtS99pD82s2eAKUPsOmvIqnPOmdlQ0f/zwBPOubqR+oiccw8CDwIsWrTILV++fKTsjcmGDRuI9zG9kCrlgNQqyyWbH2ffrqNcf89dWGVlpAn7AmrFP3h4C+UFrXzy/Ss861eN93WZue15+nOyWb78yrgdczRGKsezu4/DC1u4+dpFLKwadgiOL6TSZ0VluTAjBmzn3I3D7TOz42Y21Tl31MymAieGSHY1cJ2ZfR7IBzLNrN05d77+bpFgWbuWS8PdPJJfxsHCacw6dCgydxrGFLQ7e/p4ft9JPra4MqUGQVUV57LneJvX2Xib2mgzfbX6sCUAYu3Dfhy4M/r4TuCxcxM451Y65yqdc9XAvcAjCtaSclav5h0Z7QCsnxWtRV7A3OmNexvo7gvznkvK4p1DT1UW51LX1EV/+MK74BKhtqGDCdnpFGlKlwRArAH768C7zWwfkf7prwOY2SIzeyjWzIkExuHDlIZ6md1wiPUzF521fSye3nWcSTkZXBng9cOHUl2cR09/mGOtp73OylkON0VGiKdSa4akrpgCtnOu0Tl3g3NutnPuRudcU3T7Fufc3UOk/8/hBpyJBFp0jvS7Dmxh8/T5tGfmnLV9NPr6wzz7xnHeNW8yGaHUWtOoqig6UrzBX1O76po7qShQc7gEQ2p9K4h4Zc0aSEtj+YEt9IYy+G3VpWOeO73lUDMtnb0p1xwOb91m85CPlih1zlHf0kV5YY7XWREZFQVskXhYuRKqqlgUamdCdwfrL10x5rnTT+08TmZ6GsvmlCYwo96YOimHjJD5ai52Y0cPp3vDVChgS0CMOEpcREapqIiMmoNct3Yr60uLcB+/gdH2jDrneHr3Ma69qIS8rNT7WIbSjOlFub66zWZ9dKnU8gIFbAkG1bBF4mz53Mkcb+1me/2pUf/NG8faONLUxbtTsDl8QFWRv26zObC2eUWh+rAlGBSwReLspkumkJ+VzoMbD476b57edRwzuOHiyQnMmbcGbrMZy+qK8VTfEvnnQX3YEhQK2CJxNik3g09cXcUvtx9l/4n2EdM75/j1jmNcUVnI5AnZScihNyoKc2jv7qO1q8/rrACRJvEJ2elMyvH3kqQiAxSwRRLgU9fOICs9jX/ZsH/EtE/vOs6uo63cdnl5EnLmnfLtWwComzkPqqth7VpP81PX3KX+awkUBWyRBCjJz2Llkioee+1NDp9nZPTp3n6+9stdzCnL52NXTk9iDpNs7VrKv/E1AOonlMLA0q0eBu36li71X0ugKGCLJMiqZTMJpdl5a9kPPX+QI01d/PX755OeYoulnGX1asqPR1Z9q58U7ae/gKVb48U5R11zl6Z0SaCk8DeEiLfKJmbz0UXT+cmrddS3dL1t/5stXTyw/gC3LJjC0otKPMhhEh0+TFFXK9m9p6mfOPms7V5o7eqjvbtPAVsCRQFbJIE+u3wWzsF3fnPgbfv+7ldvEHaOv3jvxR7kLMkqKzGgvPXkWzXs6HYvHGmOjhBXH7YEiAK2SAKVF+TwoSsqWPfKEZ7YfpT6li6cc2w62Mj/vP4mn71+FtOLxkE/6po1kJtL+akT1E+MruQ2xqVb42mgxUN92BIkqbekkojPfGHFRTy56xifX/sqAIW5GZgZ0yZl89nrZ3mcuySJLtFa/rMd7Jw4C6qqIsF6DEu3xtPAoimagy1BooAtkmCVxbm8/OUb2H20lR1vtrKz/hR7j7fxxzfMJicz5HX2kmflSiqm7afxyT107T3gadnrm7vIyQhRmKs52BIcCtgiSZCdEeLyykIuryz0OiueGugzrm/p4qLJ+Z7lo665k4rCHN0HWwJFfdgikjQDTdBDjZpPJt1WU4JIAVtEkuZMDbvZ+4CtKV0SNArYIpI0ZROzSU+zMzfe8EJ7dx8tnb2UF2iEuASLAraIJE0ozZgyKdvTGnb9mdtqqoYtwaKALSJJVV6Q42kfdl2zbqspwaSALSJJVV6Y420Nu0U1bAkmBWwRSaqKghyOtZ6mtz/syfnrmrvITE+jJC/Lk/OLXCgFbBFJqvLCHMIOjp067cn565u7qCjIIS1Nc7AlWBSwRSSpBkZne9WPXdfcqf5rCSQFbBFJqjOLp3jUj6052BJUCtgiklRTJ2UD3tSwT/f209Deo9tqSiApYItIUmVnhCidkOVJDbuuWbfVlOBSwBaRpPNqLvbAOdWHLUGkgC0iSVde6E3APrNoiprEJYAUsEUk6SqiNexw2CX1vPXNXaSnGWUTs5N6XpF4iClgm1mRmT1tZvuiv4e82a+ZVZrZU2a228x2mVl1LOcVkWArL8yhpy9MQ0d3Us9b19zF1IJsQpqDLQEUaw37PuBZ59xs4Nno86E8AnzDOXcxsBg4EeN5RSTAvLrNZn1Ll5rDJbBiDdi3Ag9HHz8M3HZuAjO7BEh3zj0N4Jxrd855d289EfHcmbnYSe7Hrm/u0ghxCaxYA3aZc+5o9PExoGyINHOAFjP7qZn9zsy+YWahGM8rIgHmRQ27py/M8bbTqmFLYJlz5x/0YWbPAFOG2LUaeNg5VzAobbNz7qx+bDO7HfgucDlwGPgR8IRz7rtDnGsVsAqgrKxs4bp168ZUmJG0t7eTn58f12N6IVXKASqLXyWjLJ9/poOrp6XziUsSdxOOweU40Rnmzzd28akFmVxXkZGwcyaK3l/+FO+yrFixYqtzbtGQO51zF/wD7AGmRh9PBfYMkeYq4DeDnn8CeGCkYy9cuNDF2/r16+N+TC+kSjmcU1n8Khllufn+je6u721O6DkGl+O3+066qi/9wv1238mEnjNR9P7yp3iXBdjihomLsTaJPw7cGX18J/DYEGleAQrMrDT6/F3ArhjPKyIBl+zFU+q0aIoEXKwB++vAu81sH3Bj9DlmtsjMHgJwzvUD9wLPmtl2wIB/j/G8IhJwFYU51DV3DbS8JVx9cxdmMHWSArYEU3osf+ycawRuGGL7FuDuQc+fBt4Zy7lEJLWUF+TQ3t1Ha1cfk3IT36dc39LF5AlZZKZrvSgJJr1zRcQTA7e4rGtJzixPTemSoFPAFhFPDATPuiRN7dKiKRJ0Ctgi4okzi6ckIWD3hx1HT3VpwJkEmgK2iHiiMDeD3MxQUmrYJ9pO09vvVMOWQFPAFhFPmFl0pHji+7AHavGqYUuQKWCLiGeSNRd74BwVqmFLgClgi4hnKgpzk9IkXqcatqQABWwR8UxFYQ6nunppO92b0PPUt3RF+8xjWnpCxFMK2CLimWTdZrO+WSPEJfgUsEXEM2fmYjclOGBrDrakAAVsEfHMmdXOEjhS3DkXqWEXaJUzCTYFbBHxTHFeJtkZaQltEm/u7KWrt//MPwciQaWALSKeMTPKC3ISOlJ8oPauPmwJOgVsEfFUoqd2nVk0RX3YEnAK2CLiqfLCxC6ecmbRFNWwJeAUsEXEUxWFOTR19NDR3ZeQ49c1d5GXGWJSTuLvuS2SSArYIuKpgaldiapl17dE5mCbWUKOL5IsCtgi4qmBvuVE3WYzMqVLzeESfArYIuKp6Qmeiz1QwxYJOgVsEfFUSX4Wmelp1CWgSbyrz3Gqq1eLpkhKUMAWEU+lpSVuLnZjlwM0B1tSgwK2iHiuojAxAbuhK3zm+CJBp4AtIp6rKMxJyKCzhmgNu0KDziQFKGCLiOfKC3JoaO/mdG9/XI/beNqRGUqjJD8rrscV8YICtoh47sxtNuNcy27sCjOtIJu0NM3BluBTwBYRzw30Mcd78ZSGLqcBZ5IyFLBFxHPlCZqL3dAVpkJTuiRFKGCLiOcmT8gmI2RxbRJvPd1Law9Ul+TF7ZgiXlLAFhHPhdKMqZPiO1K8tqEDgBklqmFLalDAFhFfiMzFjl+TeM2ZgJ0ft2OKeEkBW0R8Id6Lp9Q2RIJ/VbFq2JIaYgrYZlZkZk+b2b7o78Jh0v2Dme00s91m9i3Tfe5E5BzlBbmcaOumuy8+c7FrGzsoyjayM0JxOZ6I12KtYd8HPOucmw08G31+FjNbClwDvBNYAFwJXB/jeUUkxVScGSken1p2TUMHZbmqG0jqiDVg3wo8HH38MHDbEGkckA1kAllABnA8xvOKSIqZWRoZzV1zsiMux6tt7KAsT71+kjrMOXfhf2zW4pwriD42oHng+Tnp/hG4GzDg28651cMcbxWwCqCsrGzhunXrLjhvQ2lvbyc/P/gDUFKlHKCy+JUXZenodXzh2U4+MjeD987IjOlY7T2Oe57r5APVjlvn6Zr4jcoyvBUrVmx1zi0aal/6SH9sZs8AU4bYdVbQdc45M3tb9Dezi4CLgYropqfN7Drn3PPnpnXOPQg8CLBo0SK3fPnykbI3Jhs2bCDex/RCqpQDVBa/8qos/3vTMzChlOXLL43pOK8ebobnXqSyKFvXxIdUlgszYsB2zt043D4zO25mU51zR81sKnBiiGQfAF52zrVH/+ZXwNXA2wK2iIxvs0rzOBCHJvGBOdhluWoSl9QR67v5ceDO6OM7gceGSHMYuN7M0s0sg8iAs90xnldEUtCsyfkcPNke83FqGzpIM5isQWeSQmIN2F8H3m1m+4Abo88xs0Vm9lA0zY+BA8B24HXgdefc/8R4XhFJQbNK82nu7KWpoyem49Q0dlJemEO67tIlKWTEJvHzcc41AjcMsX0LkUFmOOf6gc/Ech4RGR8GRoofONlOUV7RBR+ntqGD6uI8IL53/xLxkjp4RMQ3LiqNjLY9cOLCm8Wdc9Q2dDBTN/2QFKOALSK+Ma0gh6z0NA42XPjAs8aOHtq6+3SXLkk5Ctgi4huhNGNGSV5MNeyBm34oYEuqUcAWEV+ZVZrPgRhGip+5S1exArakFgVsEfGVWaV5HG7qvOCbgNQ2dJCeZmfWJhdJFQrYIuIrsybnE3ZwqPHC7o1d29jB9KJc0kP6epPUone0iPjKzJLISPELXUClpqGTGeq/lhSkgC0ivvLWXOyxjxR3znGocWAOtkhqUcAWEV/Jy0pn6qTsCxopfqKtm86efmaU5CYgZyLeUsAWEd+50JHiB09qSpekLgVsEfGdmaV5HDzZgXNvu2PvedU2RgO2msQlBSlgi4jvzCrNp627j5Nt3WP6u9qGDjLT05hWoCldknoUsEXEd2ZF1xTfP8Zm8ZqGDqqKcgnpLl2SghSwRcR3Zk2+sJHitY0d6r+WlKWALSK+M2ViNrmZoTHNxe7u66e2sVN36ZKUpYAtIr5jZswszRtTDXtb3Sl6+sJcUVWYwJyJeEcBW0R8aVZp/pjmYm+uaQJgcXVRorIk4ikFbBHxpVml+dS3dNHVM7qbgGyqaWJu2QQK8zITnDMRbyhgi4gvvbVE6ci17L7+MFtrm1g8Q7VrSV0K2CLiS5dNLwDg5YONI6bd+WYrHT39LJmpgC2pSwFbRHypojCXOWX5PPfGiRHTqv9axgMFbBHxrRXzJrO5pom2073nTbepppEZJXlMnpidpJyJJJ8Ctoj41rvmTqYv7HhhX8OwacJhx+aaJtWuJeUpYIuIby2sKmRidvp5m8X3HG+j9XSf+q8l5Slgi4hvpYfSWDanlPV7ThIOD33nrk3RQWkaIS6pTgFbRHztXfMm09DezY43Tw25f3NtE+UFOVQU5iY5ZyLJpYAtIr52/ZxSzBiyWdy5SP/1EtWuZRxQwBYRXyvOz+LSigLWDxGwD5zsoKG9R83hMi4oYIuI771r3mRerzvFybbus7afmX+tgC3jgAK2iPjeu+ZNBmDDnrNr2ZtqGinJz2KGbqkp44ACtoj43vxpE5k8IYv1gwL2qc5eNh1sYsnMIszMw9yJJEdMAdvMPmxmO80sbGaLzpPuZjPbY2b7zey+WM4pIuOPmbFi7mSe39vAb/c38MV1v+PK/+8ZjrWe5pYFU7zOnkhSpMf49zuADwLfGS6BmYWAB4B3A3XAK2b2uHNuV4znFpFxZMW8yfxoyxFWPrSJCdnp3HHldD6yaDoLyid5nTWRpIgpYDvndgMjNUctBvY75w5G064DbgUUsEVk1JbPLeWua2awoHwityyYSk5myOssiSSVOTf06kFjOojZBuBe59yWIfbdDtzsnLs7+vwTwBLn3D1DpF0FrAIoKytbuG7dupjzNlh7ezv5+flxPaYXUqUcoLL4VaqUJVXKASqLX8W7LCtWrNjqnBuyi3nEGraZPQMM1Um02jn3WKyZG8w59yDwIMCiRYvc8uXL43l4NmzYQLyP6YVUKQeoLH6VKmVJlXKAyuJXySzLiAHbOXdjjOeoB6YPel4R3SYiIiKjlIxpXa8As81shpllAncAjyfhvCIiIikj1mldHzCzOuBq4Jdm9mR0+zQzewLAOdcH3AM8CewGHnXO7Ywt2yIiIuNLrKPEfwb8bIjtbwLvHfT8CeCJWM4lIiIynmmlMxERkQBQwBYREQkABWwREZEAiMvCKYlgZieBQ3E+bAnQEOdjeiFVygEqi1+lSllSpRygsvhVvMtS5ZwrHWqHbwN2IpjZluFWkAmSVCkHqCx+lSplSZVygMriV8ksi5rERUREAkABW0REJADGW8B+0OsMxEmqlANUFr9KlbKkSjlAZfGrpJVlXPVhi4iIBNV4q2GLiIgEkgK2iIhIAKRcwDazD5vZTjMLm9mic/Z92cz2m9keM7tpmL+fYWaboul+FL3DmKei+Xgt+lNrZq8Nk67WzLZH021JcjZHxcz+2szqB5XnvcOkuzl6nfab2X3JzudomNk3zOwNM9tmZj8zs4Jh0vnyuoz0GptZVvS9tz/6maj2IJsjMrPpZrbezHZFP/t/OkSa5WZ2atD77ite5HU0Rnq/WMS3otdlm5ld4UU+R2Jmcwe93q+ZWauZffGcNL69Lmb2H2Z2wsx2DNpWZGZPm9m+6O/CYf72zmiafWZ2Z9wy5ZxLqR/gYmAusAFYNGj7JcDrQBYwAzgAhIb4+0eBO6KP/w34nNdlOid//wR8ZZh9tUCJ13kcIf9/Ddw7QppQ9PrMBDKj1+0Sr/M+RD7fA6RHH/898PdBuS6jeY2BzwP/Fn18B/Ajr/M9TFmmAldEH08A9g5RluXAL7zO6yjLc973C5EbK/0KMOAqYJPXeR5FmULAMSKLggTiugDLgCuAHYO2/QNwX/TxfUN95oEi4GD0d2H0cWE88pRyNWzn3G7n3J4hdt0KrHPOdTvnaoD9wOLBCczMgHcBP45uehi4LYHZHZNo/j4C/JfXeUmwxcB+59xB51wPsI7I9fMV59xTLnL7WICXgQov8zNGo3mNbyXyGYDIZ+KG6HvQV5xzR51zr0YftxG5jW+5t7lKqFuBR1zEy0CBmU31OlMjuAE44JyL9+qVCeOc2wg0nbN58GdiuPhwE/C0c67JOdcMPA3cHI88pVzAPo9y4Mig53W8/UNdDLQM+hIeKo2XrgOOO+f2DbPfAU+Z2VYzW5XEfI3VPdGmvP8YpklpNNfKb+4iUusZih+vy2he4zNpop+JU0Q+I74Vbba/HNg0xO6rzex1M/uVmc1Pbs7GZKT3SxA/H3cwfEUjKNcFoMw5dzT6+BhQNkSahF2fmO6H7RUzewaYMsSu1c65x5Kdn3gYZZk+xvlr19c65+rNbDLwtJm9Ef0vManOVxbgX4GvEflS+hqRJv67kpe7sRnNdTGz1UAfsHaYw/jiuqQ6M8sHfgJ80TnXes7uV4k0x7ZHx038HJid5CyOVkq9X6LjgH4f+PIQu4N0Xc7inHNmltR50YEM2M65Gy/gz+qB6YOeV0S3DdZIpHkpPVqjGCpNQoxUJjNLBz4ILDzPMeqjv0+Y2c+INHsm/YM+2utjZv8O/GKIXaO5Vkkxiuvyh8D7gBtctANriGP44rqcYzSv8UCauuj7bxKRz4jvmFkGkWC91jn303P3Dw7gzrknzOxfzKzEOee7G1CM4v3im8/HKN0CvOqcO37ujiBdl6jjZjbVOXc02g1xYog09UT65gdUEBlTFbPx1CT+OHBHdOTrDCL/xW0enCD6hbseuD266U7ALzX2G4E3nHN1Q+00szwzmzDwmMiAqB1DpfXSOX1tH2DoPL4CzLbIiP1MIs1pjycjf2NhZjcDfw78vnOuc5g0fr0uo3mNHyfyGYDIZ+K54f4p8VK0X/27wG7n3DeHSTNloP/dzBYT+e7z3T8fo3y/PA58Mjpa/Crg1KBmWj8atmUwKNdlkMGfieHiw5PAe8ysMNrl957otth5PRIv3j9EgkAd0A0cB54ctG81kZGxe4BbBm1/ApgWfTyTSCDfD/w3kOV1maL5+k/gs+dsmwY8MSjfr0d/dhJpsvU830OU4/vAdmAbkTf/1HPLEn3+XiKjfQ/4uCz7ifRVvRb9GRhRHYjrMtRrDHyVyD8gANnRz8D+6Gdiptd5HqYc1xLpYtk26Fq8F/jswGcGuCf6+r9OZIDgUq/zPUxZhny/nFMWAx6IXrftDJoN47cfII9IAJ40aFsgrguRfzKOAr3RmPIpImM4ngX2Ac8ARdG0i4CHBv3tXdHPzX7gj+KVJy1NKiIiEgDjqUlcREQksBSwRUREAkABW0REJAAUsEVERAJAAVtERCQAFLBFREQCQAFbREQkAP5/i5s0W4T6z6MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def f_multimodal(x):\n",
    "    return -np.exp(-(x-2)**2) - 0.5*np.exp(-(x+2)**2)\n",
    "\n",
    "GD(eta=6, x0=-3, max_iter=100, diff_to_stop=0.01, delta=0.01).run(f_multimodal, plot_method='scatter')\n",
    "\n",
    "# Local minimum - Alguna solución?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradient Descent para Linear Regression\n",
    "\n",
    "Aplicar Gradient Descent para un modelo de regression lineal no tiene ninguna complejidad. Nuestro modelo viene dado por la siguiente ecuación:\n",
    "\n",
    "$$f(x)=w_0+w^T x.$$\n",
    "\n",
    "Debemos optimizar esa función para encontrar los valores de $w$ que la hacen mínima, ¿cómo medimos \"qué la hace mínima y qué no\"? Usando **Residual Sum Squares (RSS)**, que basicamente mide cuanto se aleja el modelo de los valores reales.\n",
    "\n",
    "<img src=\"Images/least_squares.png\" style=\"width: 340px;\"/>\n",
    "\n",
    "**¿Pero esto no es lo mismo que Least Squares?** Si y no... **Least Squares es un método para construir un modelo y RSS es una medida**. Se basan en lo mismo, en medir la distancia entre los valores predichos y los valores verdaderos. Entonces, aplicando RSS a nuestra regresión:\n",
    "\n",
    "$$ RSS(w) = \\sum_{n=1}^{N}[y_n-f(x_n)]^2 =  \\sum_{n=1}^{N}[y_n- (w_0 + \\sum_{d=1}^{D}w_dx_{nd}) ]^2 .$$\n",
    "\n",
    "$$ RSS(w) = \\frac{1}{2}\\sum_{n=1}^{N}[y_n-f(x_n)]^2$$\n",
    "\n",
    "Y lo que queremos es minimizar esta distancia, para que el modelo se acerque lo máximo posible a los valores verdaderos.\n",
    "\n",
    "$$\\nabla RSS(w) = X^T(Xw^t-y)$$\n",
    "\n",
    "En resumen, el gradient descendiente para una regresión lineal, se basa en resolver esta ecuación de forma iterativa:\n",
    "\n",
    "$$w^{t+1} = w^t - \\eta * \\nabla RSS(w)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}