TesisGradoSuribresta/Simulaciones/Detector_Orden.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generacion de muestras"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def generar_muestra_coseno(N, num_cosenos, varianza_ruido):\n",
    "    \"\"\"\n",
    "    Genera un vector de muestras como suma de cosenos con parámetros aleatorios y ruido gaussiano.\n",
    "    \n",
    "    Parámetros:\n",
    "    -----------\n",
    "    N : int\n",
    "        Cantidad de muestras a generar.\n",
    "    num_cosenos : int\n",
    "        Número de cosenos a sumar.\n",
    "    varianza_ruido : float\n",
    "        Varianza del ruido gaussiano agregado.\n",
    "    \n",
    "    Retorna:\n",
    "    --------\n",
    "    x : ndarray de shape (N,)\n",
    "        Vector de muestras resultante.\n",
    "    parametros : list of dict\n",
    "        Lista con las constantes de cada coseno (amplitud, frecuencia y fase).\n",
    "    varianza_ruido : float\n",
    "        Varianza del ruido utilizado.\n",
    "    \"\"\"\n",
    "    t = np.arange(N) / N   # tiempo normalizado entre 0 y 1\n",
    "    \n",
    "    # Parámetros aleatorios de los cosenos\n",
    "    amplitudes = np.random.uniform(0.5, 2.0, num_cosenos)      # amplitudes aleatorias\n",
    "    frecuencias = np.random.uniform(1, 20, num_cosenos)        # frecuencias aleatorias (en Hz normalizado)\n",
    "    fases = np.random.uniform(0, 2*np.pi, num_cosenos)         # fases aleatorias\n",
    "    \n",
    "    # Generar señal base como suma de cosenos\n",
    "    x_clean = np.zeros(N)\n",
    "    for A, f, phi in zip(amplitudes, frecuencias, fases):\n",
    "        x_clean += A * np.cos(2 * np.pi * f * t + phi)\n",
    "    \n",
    "    # Agregar ruido gaussiano\n",
    "    ruido = np.random.normal(0, np.sqrt(varianza_ruido), N)\n",
    "    x = x_clean + ruido\n",
    "    \n",
    "    # Empaquetar parámetros\n",
    "    parametros = []\n",
    "    for A, f, phi in zip(amplitudes, frecuencias, fases):\n",
    "        parametros.append({\n",
    "            \"amplitud\": A,\n",
    "            \"frecuencia\": f,\n",
    "            \"fase\": phi\n",
    "        })\n",
    "    \n",
    "    return x, parametros, varianza_ruido\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimador de "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "def modelo_cosenos(t, *params):\n",
    "    \"\"\"\n",
    "    Modelo de suma de cosenos.\n",
    "    \n",
    "    params = [A1, f1, phi1, A2, f2, phi2, ..., An, fn, phin]\n",
    "    \"\"\"\n",
    "    num_cosenos = len(params) // 3\n",
    "    y = np.zeros_like(t)\n",
    "    for i in range(num_cosenos):\n",
    "        A = params[3*i]\n",
    "        f = params[3*i + 1]\n",
    "        phi = params[3*i + 2]\n",
    "        y += A * np.cos(2 * np.pi * f * t + phi)\n",
    "    return y\n",
    "\n",
    "def estimar_parametros(y, num_cosenos, max_iter=10000):\n",
    "    \"\"\"\n",
    "    Estima los parámetros de una señal que es suma de cosenos mediante mínimos cuadrados no lineales.\n",
    "    \n",
    "    Parámetros:\n",
    "    -----------\n",
    "    y : ndarray\n",
    "        Vector de muestras observadas.\n",
    "    num_cosenos : int\n",
    "        Número de cosenos que componen la señal.\n",
    "    max_iter : int\n",
    "        Número máximo de iteraciones del optimizador.\n",
    "    \n",
    "    Retorna:\n",
    "    --------\n",
    "    parametros : list of dict\n",
    "        Estimación de los parámetros (amplitud, frecuencia, fase).\n",
    "    \"\"\"\n",
    "    N = len(y)\n",
    "    t = np.arange(N) / N\n",
    "    \n",
    "    # Valores iniciales aproximados (necesarios para que converja)\n",
    "    A0 = np.ones(num_cosenos) * (np.std(y) / num_cosenos)\n",
    "    f0 = np.linspace(1, 10, num_cosenos)   # inicializar frecuencias en un rango\n",
    "    phi0 = np.zeros(num_cosenos)\n",
    "    \n",
    "    p0 = []\n",
    "    for A, f, phi in zip(A0, f0, phi0):\n",
    "        p0 += [A, f, phi]\n",
    "    p0 = np.array(p0)\n",
    "    \n",
    "    # Ajuste no lineal\n",
    "    popt, _ = curve_fit(modelo_cosenos, t, y, p0=p0, maxfev=max_iter)\n",
    "    \n",
    "    # Empaquetar parámetros\n",
    "    parametros = []\n",
    "    for i in range(num_cosenos):\n",
    "        A = popt[3*i]\n",
    "        f = popt[3*i+1]\n",
    "        phi = popt[3*i+2]\n",
    "        parametros.append({\n",
    "            \"amplitud\": A,\n",
    "            \"frecuencia\": f,\n",
    "            \"fase\": phi\n",
    "        })\n",
    "    \n",
    "    return parametros\n",
    "\n",
    "def reconstruir_senal(N, parametros):\n",
    "    \"\"\"\n",
    "    Reconstruye la señal a partir de parámetros de cosenos.\n",
    "    \n",
    "    parametros : list of dict con {\"amplitud\", \"frecuencia\", \"fase\"}\n",
    "    \"\"\"\n",
    "    t = np.arange(N) / N\n",
    "    y = np.zeros(N)\n",
    "    for p in parametros:\n",
    "        y += p[\"amplitud\"] * np.cos(2 * np.pi * p[\"frecuencia\"] * t + p[\"fase\"])\n",
    "    return t, y\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def estimador_de_orden(y, orden_max, C, max_iter=10000, epsilon=1e-12, devolver_parametros=False):\n",
    "    \"\"\"\n",
    "    Estima el orden óptimo (número de cosenos) usando un criterio:\n",
    "        Criterio(k) = N * log(var_residuo_k + epsilon) + C * k\n",
    "\n",
    "    Parámetros\n",
    "    ----------\n",
    "    y : ndarray (N,)\n",
    "        Señal observada.\n",
    "    orden_max : int\n",
    "        Orden máximo a evaluar (número máximo de cosenos).\n",
    "    C : float\n",
    "        Constante que multiplica al término de penalización por el orden k.\n",
    "    max_iter : int\n",
    "        Máx. iteraciones para el optimizador en estimar_parametros.\n",
    "    epsilon : float\n",
    "        Pequeño valor para estabilizar el logaritmo.\n",
    "    devolver_parametros : bool\n",
    "        Si True, devuelve también los parámetros del orden óptimo.\n",
    "\n",
    "    Retorna\n",
    "    -------\n",
    "    orden_optimo : int\n",
    "        Orden estimado.\n",
    "    varianzas : list[float]\n",
    "        Varianzas del residuo para k = 1..orden_max.\n",
    "    criterios : list[float]\n",
    "        Valores del criterio para k = 1..orden_max.\n",
    "    (opcional) parametros_optimos : list[dict]\n",
    "        Parámetros {amplitud, frecuencia, fase} del orden óptimo (si devolver_parametros=True).\n",
    "    \"\"\"\n",
    "    N = len(y)\n",
    "    varianzas = []\n",
    "    criterios = []\n",
    "    mejores_parametros = None\n",
    "    parametros_k_opt = None\n",
    "\n",
    "    for k in range(1, orden_max + 1):\n",
    "        try:\n",
    "            # Estimar parámetros para orden k\n",
    "            params_k = estimar_parametros(y, num_cosenos=k, max_iter=max_iter)\n",
    "\n",
    "            # Reconstruir y calcular residuo\n",
    "            _, y_hat = reconstruir_senal(N, params_k)\n",
    "            residuo = y - y_hat\n",
    "            var_est = float(np.mean(residuo**2))\n",
    "\n",
    "            # Guardar métricas\n",
    "            varianzas.append(var_est)\n",
    "            criterio_k = N * np.log(var_est + epsilon) + C * k\n",
    "            criterios.append(criterio_k)\n",
    "\n",
    "            # Guardar los parámetros si van siendo los mejores\n",
    "            if devolver_parametros:\n",
    "                if (parametros_k_opt is None) or (criterio_k < criterios[parametros_k_opt - 1]):\n",
    "                    parametros_k_opt = k\n",
    "                    mejores_parametros = params_k\n",
    "\n",
    "        except Exception:\n",
    "            # Si falla el ajuste para este orden, lo penalizamos fuertemente\n",
    "            varianzas.append(np.inf)\n",
    "            criterios.append(np.inf)\n",
    "\n",
    "    # Selección del orden óptimo\n",
    "    orden_optimo = int(np.argmin(criterios) + 1)\n",
    "\n",
    "    if devolver_parametros:\n",
    "        # Si no se setearon (por fallos), intentar reevaluar los del orden óptimo\n",
    "        if mejores_parametros is None or parametros_k_opt != orden_optimo:\n",
    "            try:\n",
    "                mejores_parametros = estimar_parametros(y, num_cosenos=orden_optimo, max_iter=max_iter)\n",
    "            except Exception:\n",
    "                mejores_parametros = None\n",
    "        return orden_optimo, varianzas, criterios, mejores_parametros\n",
    "\n",
    "    return orden_optimo, varianzas, criterios\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Orden real: 3\n",
      "Orden estimado: 4\n",
      "Varianzas: [3.7188377032816926, 3.7099064438728937, 3.6910485044075623, 2.4082977259715626, 2.4043273936036753, inf, 2.389950171990722, inf]\n",
      "Criterios: [np.float64(662.2055870548606), np.float64(666.503329508359), np.float64(669.455282613183), np.float64(461.46008007889515), np.float64(466.1350972119629), inf, np.float64(474.1362585935285), inf]\n",
      "Parámetros óptimos: [{'amplitud': np.float64(-0.0480616545387829), 'frecuencia': np.float64(-0.01254727613082784), 'fase': np.float64(5.6630225145901925)}, {'amplitud': np.float64(1.6294124254561178), 'frecuencia': np.float64(3.242768470130581), 'fase': np.float64(1.6525282377898725)}, {'amplitud': np.float64(-0.08429321831662702), 'frecuencia': np.float64(6.9129993235227385), 'fase': np.float64(0.5002733545284419)}, {'amplitud': np.float64(-0.11051057544789585), 'frecuencia': np.float64(9.999389491130865), 'fase': np.float64(0.3362456153446631)}]\n"
     ]
    }
   ],
   "source": [
    "# Señal sintética\n",
    "x, params_true, _ = generar_muestra_coseno(N=500, num_cosenos=3, varianza_ruido=0.05)\n",
    "\n",
    "# Estimo orden con constante C=5.0 (ajustá C a tu gusto)\n",
    "orden_opt, vars_k, crits_k, params_opt = estimador_de_orden(\n",
    "    x, orden_max=8, C=5.50, devolver_parametros=True\n",
    ")\n",
    "\n",
    "print(\"Orden real:\", len(params_true))\n",
    "print(\"Orden estimado:\", orden_opt)\n",
    "print(\"Varianzas:\", vars_k)\n",
    "print(\"Criterios:\", crits_k)\n",
    "print(\"Parámetros óptimos:\", params_opt)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Orden real: 3\n",
      "Tasa de acierto: 0.13333333333333333\n",
      "Órdenes estimados: [5, 4, 1, 6, 5, 4, 4, 4, 7, 2, 1, 1, 3, 4, 5, 3, 1, 5, 3, 5, 5, 4, 6, 2, 4, 6, 4, 4, 4, 3]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAHHCAYAAACle7JuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABPh0lEQVR4nO3deXxM9/4/8NckkskuCQmCLJIgIhSxL6FCGqGWWqpcBKVEI7Q00VpSS2hVbbe2fivaaylaS6mIXe2JtUiJPVXEmkgi63x+f+SXqTFJJDHj5IzX8/E4D2fOnDnnfc6MyWs+53POUQghBIiIiIhkyEjqAoiIiIjKikGGiIiIZItBhoiIiGSLQYaIiIhki0GGiIiIZItBhoiIiGSLQYaIiIhki0GGiIiIZItBhgxSVlYWZs2ahZ07d0pdik48evQIkZGROH78uNSlEBGVKwwyMjdt2jQoFIrXsq727dujffv26sf79++HQqHAxo0bX8v6n6dQKDBt2rQinx8/fjxWr16N5s2bv5Z6hgwZAldXV70sWwiBQYMGYf/+/WjUqNErL8/V1RVDhgx59cLeAC9+5ssDvn+G4WXfYQVe53e8XDHIlCPR0dFQKBTqwczMDE5OTggICMDChQvx9OlTnaznn3/+wbRp03DmzBmdLK+8Wb9+PTZv3owdO3bA1tZW6nJe2VdffYUbN25g06ZNMDU1lbocg3Px4kVMmzYNN27ckLqUN9qzZ88wbNgw1K9fHxUrVoSVlRUaNmyIBQsWICcnR2v+J0+eYMSIEXBwcIClpSU6dOiAU6dOSVA5Sa2C1AWQti+//BJubm7IycnB3bt3sX//foSFhWHevHnYunUrGjRooJ73iy++QHh4eKmW/88//yAyMhKurq546623Svy62NjYUq1Hn549e4YKFbQ/vkII/P3339ixYwecnZ0lqEy3MjMzkZubi99//90gQll5dPHiRURGRqJ9+/ZarWrl6TNv6J49e4YLFy6gS5cucHV1hZGREY4cOYJx48bh+PHjWLNmjXpelUqFoKAgnD17FhMmTEDlypXx3XffoX379jh58iQ8PT0l3JKSKeo7jEqPe7EcCgwMhK+vr/pxREQE9u7di65du+Ldd99FQkICzM3NAQAVKlTQ+3+GjIwMWFhYlKvWADMzs0KnKxQKjB8//jVXoz9mZmb4/PPPS/Wa9PR0WFpa6qmiN0t5+swbOnt7exw7dkxj2kcffYSKFSti8eLFmDdvHqpWrQoA2LhxI44cOYINGzagd+/eAIC+ffuidu3amDp1qkbokUJR/wdVKhWys7NhZmZW5HcYlR4PLcnE22+/jcmTJ+PmzZv43//+p55e2PHTXbt2oU2bNrC1tYWVlRXq1KmDSZMmAcjv19K0aVMAQHBwsPowVnR0NID8PgH169fHyZMn0a5dO1hYWKhfW1R/gby8PEyaNAlVq1aFpaUl3n33XSQlJWnMU9Rx/cKWmZmZiWnTpqF27dowMzNDtWrV0KtXL1y9elU9T2HHl0+fPo3AwEDY2NjAysoKHTt21PpiLDh8d/jwYYwfP17dLN2zZ0/cv39fq77CbN68GfXr14eZmRnq16+PTZs2FTqfSqXC/Pnz4e3tDTMzM1SpUgUjR47E48ePS7SevXv3om3btrC0tIStrS26d++OhIQEjXkK3v+LFy/igw8+gJ2dHdq0aQMgv3VqxowZqFGjBiwsLNChQwdcuHCh0HU9efIEYWFhqFmzJpRKJTw8PDBnzhyoVCr1PDdu3IBCocDcuXOxfPlyuLu7Q6lUomnTpoiLi9Na5l9//YXevXvD3t4eZmZm8PX1xdatWzXmycnJQWRkJDw9PWFmZoZKlSqhTZs22LVr10v3T0lqBoB169ahSZMmsLa2ho2NDXx8fLBgwQIA+Z+HPn36AAA6dOig/v+wf/9+AEX3C1u/fj0iIyNRvXp1WFtbo3fv3khJSUFWVhbCwsLg6OgIKysrBAcHIysrS6OelStX4u2334ajoyOUSiXq1auHJUuWaG1fad6/a9euoU+fPrC3t4eFhQVatGiB7du3a823aNEieHt7w8LCAnZ2dvD19S3RH/2srCxMnToVHh4eUCqVqFmzJiZOnKi1bQqFAmPGjFH/H1EqlfD29kZMTMxL11GUglayJ0+eqKdt3LgRVapUQa9evdTTHBwc0LdvX2zZskWrrsLs2LEDfn5+6s9F06ZNtfbFhg0b0KRJE5ibm6Ny5coYOHAgbt++rTHPkCFDYGVlhatXr6JLly6wtrbGgAEDAPy7P1avXg1vb28olUr1vijsO+zQoUNo2rQpzMzM4O7ujmXLlhVae0k/Q/Hx8QgICEDlypVhbm4ONzc3DB069KX7Rm7YIiMj//nPfzBp0iTExsbiww8/LHSeCxcuoGvXrmjQoAG+/PJLKJVKXLlyBYcPHwYAeHl54csvv8SUKVMwYsQItG3bFgDQqlUr9TIePnyIwMBAvP/++xg4cCCqVKlSbF0zZ86EQqHAZ599huTkZMyfPx/+/v44c+aMuuWopPLy8tC1a1fs2bMH77//PsaOHYunT59i165dOH/+PNzd3Yvc7rZt28LGxgYTJ06EiYkJli1bhvbt2+PAgQNanX4//vhj2NnZYerUqbhx4wbmz5+PMWPG4Oeffy62vtjYWLz33nuoV68eoqKi8PDhQwQHB6NGjRpa844cORLR0dEIDg5GaGgorl+/jsWLF+P06dM4fPgwTExMilzP7t27ERgYiFq1amHatGl49uwZFi1ahNatW+PUqVNah0D69OkDT09PzJo1C0IIAMCUKVMwY8YMdOnSBV26dMGpU6fQuXNnZGdna7w2IyMDfn5+uH37NkaOHAlnZ2ccOXIEERERuHPnDubPn68x/5o1a/D06VOMHDkSCoUCX331FXr16oVr166pt+nChQto3bo1qlevjvDwcFhaWmL9+vXo0aMHfvnlF/Ts2RNAfhCLiorC8OHD0axZM6SmpiI+Ph6nTp1Cp06ditw/Ja15165d6N+/Pzp27Ig5c+YAABISEnD48GGMHTsW7dq1Q2hoKBYuXIhJkybBy8sLANT/FiUqKgrm5uYIDw/HlStXsGjRIpiYmMDIyAiPHz/GtGnTcOzYMURHR8PNzQ1TpkxRv3bJkiXw9vbGu+++iwoVKuC3337D6NGjoVKpEBISop6vpO/fvXv30KpVK2RkZCA0NBSVKlXCqlWr8O6772Ljxo3qfb1ixQqEhoaid+/eGDt2LDIzM3Hu3DkcP34cH3zwQZHbqlKp8O677+LQoUMYMWIEvLy88Oeff+Lbb7/F5cuXsXnzZo35Dx06hF9//RWjR4+GtbU1Fi5ciPfeew+3bt1CpUqVit2vAJCdnY3U1FQ8e/YM8fHxmDt3LlxcXODh4aGe5/Tp02jcuDGMjDR/izdr1gzLly/H5cuX4ePjU+Q6oqOjMXToUHh7eyMiIgK2trY4ffo0YmJi1Pui4P9u06ZNERUVhXv37mHBggU4fPgwTp8+rXGoNzc3FwEBAWjTpg3mzp0LCwsL9XN79+7F+vXrMWbMGFSuXLnIkwL+/PNPdO7cGQ4ODpg2bRpyc3MxderUQr9/S/IZSk5OVi8vPDwctra2uHHjBn799deXvgeyI6jcWLlypQAg4uLiipynYsWKolGjRurHU6dOFc+/jd9++60AIO7fv1/kMuLi4gQAsXLlSq3n/Pz8BACxdOnSQp/z8/NTP963b58AIKpXry5SU1PV09evXy8AiAULFqinubi4iMGDB790mT/88IMAIObNm6c1r0qlUo8DEFOnTlU/7tGjhzA1NRVXr15VT/vnn3+EtbW1aNeunXpawT729/fXWN64ceOEsbGxePLkidZ6n/fWW2+JatWqacwXGxsrAAgXFxf1tD/++EMAEKtXr9Z4fUxMTKHTC1uPo6OjePjwoXra2bNnhZGRkRg0aJB6WsH7379/f43XJycnC1NTUxEUFKSxnZMmTRIANN6L6dOnC0tLS3H58mWNZYSHhwtjY2Nx69YtIYQQ169fFwBEpUqVxKNHj9TzbdmyRQAQv/32m3pax44dhY+Pj8jMzFRPU6lUolWrVsLT01M9rWHDhiIoKKjYfVGYktY8duxYYWNjI3Jzc4tc1oYNGwQAsW/fPq3nivrM169fX2RnZ6un9+/fXygUChEYGKjx+pYtW2p8LoQQIiMjQ2s9AQEBolatWurHpXn/wsLCBADxxx9/qKc9ffpUuLm5CVdXV5GXlyeEEKJ79+7C29u7yP1QlJ9++kkYGRlpLF8IIZYuXSoAiMOHD6unARCmpqbiypUr6mlnz54VAMSiRYtKtL61a9cKAOrB19dXnDt3TmMeS0tLMXToUK3Xbt++XQAQMTExRS7/yZMnwtraWjRv3lw8e/ZM47mCfZ2dnS0cHR1F/fr1NebZtm2bACCmTJminjZ48GABQISHh2utC4AwMjISFy5cKPS5F7/DzMzMxM2bN9XTLl68KIyNjTW+44Uo2Wdo06ZNL/17Yih4aElmrKysij17qeBXwpYtW7Sa2EtKqVQiODi4xPMPGjQI1tbW6se9e/dGtWrV8Pvvv5d63b/88gsqV66Mjz/+WOu5ok5BzMvLQ2xsLHr06IFatWqpp1erVg0ffPABDh06hNTUVI3XjBgxQmN5bdu2RV5eHm7evFlkbXfu3MGZM2cwePBgVKxYUT29U6dOqFevnsa8GzZsQMWKFdGpUyc8ePBAPTRp0gRWVlbYt2/fS9czZMgQ2Nvbq6c3aNAAnTp1KnS/fvTRRxqPd+/ejezsbHz88cca2xkWFqb12g0bNqBt27aws7PTqNXf3x95eXk4ePCgxvz9+vWDnZ2d+nFBq961a9cA5F/zZu/evejbty+ePn2qXt7Dhw8REBCAxMREdfO8ra0tLly4gMTExCL3R2FKWrOtrS3S09NLdKiqNAYNGqTRota8eXMIIbSa7Zs3b46kpCTk5uaqpz3fSpmSkoIHDx7Az88P165dQ0pKCoDSvX+///47mjVrpj6kCOR/T4wYMQI3btzAxYsXAeTvi7///rvQw4DF2bBhA7y8vFC3bl2Nff32228DgNZn2d/fX6PltEGDBrCxsVF/Pl6mQ4cO2LVrFzZs2ICPPvoIJiYmSE9P15jn2bNnUCqVWq8t6Hfy7NmzIpe/a9cuPH36FOHh4Vr9VAr2dXx8PJKTkzF69GiNeYKCglC3bt1CD9uNGjWq0PX5+flpfT+8KC8vDzt37kSPHj00TlLw8vJCQECA1vwl+QwV/C3Ytm1boWd9GRIGGZlJS0vTCA0v6tevH1q3bo3hw4ejSpUqeP/997F+/fpShZrq1auXqpPji2cIKBQKeHh4lOl01qtXr6JOnTql6sB8//59ZGRkoE6dOlrPeXl5QaVSafXZefGMpoI/zMX1XykIOYWdEfHiuhMTE5GSkgJHR0c4ODhoDGlpaUhOTn7peorangcPHmh9sbu5uZWoVgcHB40QUlBrTEyMVp3+/v4AoFXry/bdlStXIITA5MmTtZY5depUjWV++eWXePLkCWrXrg0fHx9MmDAB586dK3LflLbm0aNHo3bt2ggMDESNGjUwdOjQV+qvUdQ+KAi2NWvW1JquUqnUf1wA4PDhw/D391f3fXJwcFD3QyuYrzTv382bN4v8rDy/rM8++wxWVlZo1qwZPD09ERISoj7kXJzExERcuHBBa1/Xrl0bwMs/H0D+Z6SkfcOqVKkCf39/9O7dG0uWLEHXrl3RqVMn3L17Vz2Publ5of1gMjMz1c8XpaCvXf369Yucp7j/g3Xr1tX6wVOhQoVCDy8D2v83C3P//n08e/asRN8tQMk+Q35+fnjvvfcQGRmJypUro3v37li5cmWJ+g/JDfvIyMjff/+NlJQUjWPFLzI3N8fBgwexb98+bN++HTExMfj555/x9ttvIzY2FsbGxi9dT2n7tZREca0pJalJ14pap/j//UtelUqlgqOjI1avXl3o8w4ODjpZT4FXec9UKhU6deqEiRMnFvp8wR+sAi/bdwWh+dNPPy301yQA9We4Xbt2uHr1KrZs2YLY2Fh8//33+Pbbb7F06VIMHz78lWt2dHTEmTNnsHPnTuzYsQM7duzAypUrMWjQIKxatarI5b9MUfvgZfvm6tWr6NixI+rWrYt58+ahZs2aMDU1xe+//45vv/22zK2oJeHl5YVLly5h27ZtiImJwS+//ILvvvsOU6ZMQWRkZJGvU6lU8PHxwbx58wp9/sXwpuv/W71798bnn3+OLVu2YOTIkQDyW1vv3LmjNW/BNCcnpzKtq6yUSqVWf50Cuv4+LelnqOBipceOHcNvv/2GnTt3YujQofjmm29w7NgxWFlZ6bQuKTHIyMhPP/0EAEX+cShgZGSEjh07omPHjpg3bx5mzZqFzz//HPv27YO/v7/OrxL54mEBIQSuXLmicb0bOzs7jbMOCty8eVPjcJC7uzuOHz+OnJycYjvDPs/BwQEWFha4dOmS1nN//fUXjIyMtL5sy8LFxQWA9vYC0Fq3u7s7du/ejdatW5f6i6xgPUVtT+XKlV96evXztT6/f+/fv6/1y9jd3R1paWnq1oxXVbA+ExOTEi3T3t4ewcHBCA4ORlpaGtq1a4dp06YVG2RKU7OpqSm6deuGbt26QaVSYfTo0Vi2bBkmT54MDw+P13rV1N9++w1ZWVnYunWrRsvFi4dnSvP+ubi4FPlZeX5ZAGBpaYl+/fqhX79+yM7ORq9evTBz5kxEREQUeTqwu7s7zp49i44dO0pyhdmCw0TPt2q99dZb+OOPP6BSqTQCxPHjx2FhYaEVvp9XcNjr/PnzRf4ofP7/YMEhtAKXLl3S2Ke64ODgAHNz8xJ9t5T0M1SgRYsWaNGiBWbOnIk1a9ZgwIABWLduXbH/v+SGh5ZkYu/evZg+fTrc3NzUp/YV5tGjR1rTCi56V9CkWPBHsLBgURY//vijRr+djRs34s6dOwgMDFRPc3d3x7FjxzTOuNi2bZvWIZ/33nsPDx48wOLFi7XWU9QvOmNjY3Tu3BlbtmzROJx17949rFmzBm3atIGNjU1ZN0+tWrVqeOutt7Bq1SqNL9Vdu3ap+yEU6Nu3L/Ly8jB9+nSt5eTm5ha7759fz/PznT9/HrGxsejSpctLa/X394eJiQkWLVqksd9ePAOpoNajR48Wel+qJ0+eaPTvKAlHR0e0b98ey5YtK/RX8/OnuT98+FDjOSsrK3h4eLy0+bukNb+4fCMjI3XA1tf/h+IUtFY8/56kpKRg5cqVGvOV5v3r0qULTpw4gaNHj6qnpaenY/ny5XB1dVX3z3hxX5iamqJevXoQQhTbh6Jv3764ffs2VqxYofXcs2fPtA5zltWDBw8K/T/+/fffA4DGtbV69+6Ne/fuaZyB8+DBA2zYsAHdunUrtP9Mgc6dO8Pa2hpRUVHqQ1EFCtbv6+sLR0dHLF26VOOzuGPHDiQkJCAoKKhsG1kEY2NjBAQEYPPmzbh165Z6ekJCgtZnvKSfocePH2vtzxf/FhgKtsiUQzt27MBff/2F3Nxc3Lt3D3v37sWuXbvg4uKCrVu3FnshpS+//BIHDx5EUFAQXFxckJycjO+++w41atRQdwZ0d3eHra0tli5dCmtra1haWqJ58+YlOpZbGHt7e7Rp0wbBwcG4d+8e5s+fDw8PD41TxIcPH46NGzfinXfeQd++fXH16lX873//0zqdetCgQfjxxx8xfvx4nDhxAm3btkV6ejp2796N0aNHo3v37oXWMGPGDPX1c0aPHo0KFSpg2bJlyMrKwldffVWm7SpMVFQUgoKC0KZNGwwdOhSPHj1SX5sjLS1NPZ+fnx9GjhyJqKgonDlzBp07d4aJiQkSExOxYcMGLFiwQH0hr8J8/fXXCAwMRMuWLTFs2DD16dcVK1Ys0f1ZHBwc8OmnnyIqKgpdu3ZFly5dcPr0aezYsQOVK1fWmHfChAnYunUrunbtiiFDhqBJkyZIT0/Hn3/+iY0bN+LGjRtar3mZ//73v2jTpg18fHzw4YcfolatWrh37x6OHj2Kv//+G2fPngUA1KtXD+3bt0eTJk1gb2+P+Ph4bNy4EWPGjCl2+SWtefjw4Xj06BHefvtt1KhRAzdv3sSiRYvw1ltvqfuQvPXWWzA2NsacOXOQkpICpVKpvkaHrnXu3FndQjRy5EikpaVhxYoVcHR01Ah9pXn/wsPDsXbtWgQGBiI0NBT29vZYtWoVrl+/jl9++UXdYtG5c2dUrVoVrVu3RpUqVZCQkIDFixcjKCio2H53//nPf7B+/Xp89NFH2LdvH1q3bo28vDz89ddfWL9+PXbu3KkRMsrqf//7H5YuXarutP/06VPs3LkTu3btQrdu3TRaRnr37o0WLVogODgYFy9eVF/ZNy8vr9jDZABgY2ODb7/9FsOHD0fTpk3V1186e/YsMjIysGrVKpiYmGDOnDkIDg6Gn58f+vfvrz792tXVFePGjXvl7X1RZGQkYmJi0LZtW4wePRq5ubnq75bn+42V9DO0atUqfPfdd+jZsyfc3d3x9OlTrFixAjY2NiX6MSQrr/9EKSpKwanBBYOpqamoWrWq6NSpk1iwYIHGKc4FXjz9es+ePaJ79+7CyclJmJqaCicnJ9G/f3+t01S3bNki6tWrJypUqKBxKrafn1+Rp2gWdSrq2rVrRUREhHB0dBTm5uYiKChI4xTCAt98842oXr26UCqVonXr1iI+Pl5rmULkn1r4+eefCzc3N2FiYiKqVq0qevfurXFqNV44dVEIIU6dOiUCAgKElZWVsLCwEB06dBBHjhwpdB+/eEpiwbYUdgrui3755Rfh5eUllEqlqFevnvj111/F4MGDtU6zFUKI5cuXiyZNmghzc3NhbW0tfHx8xMSJE8U///zz0vXs3r1btG7dWpibmwsbGxvRrVs3cfHiRY15Ct7/wk63z8vLE5GRkaJatWrC3NxctG/fXpw/f77QU+GfPn0qIiIihIeHhzA1NRWVK1cWrVq1EnPnzlWfZlxw+vXXX3+tta7C3o+rV6+KQYMGiapVqwoTExNRvXp10bVrV7Fx40b1PDNmzBDNmjUTtra2wtzcXNStW1fMnDlT49TmopSk5o0bN4rOnTsLR0dHYWpqKpydncXIkSPFnTt3NJa1YsUKUatWLfWprgWfg6I+8xs2bNB4fVGfq8Len61bt4oGDRoIMzMz4erqKubMmaO+7MD169fV85Xm/bt69aro3bu3sLW1FWZmZqJZs2Zi27ZtGvMsW7ZMtGvXTlSqVEkolUrh7u4uJkyYIFJSUl66r7Ozs8WcOXOEt7e3UCqVws7OTjRp0kRERkZqvB6ACAkJ0Xp9UZdfeF5cXJzo06ePcHZ2FkqlUlhaWorGjRuLefPmiZycHK35Hz16JIYNGyYqVaokLCwshJ+fX6lONd66dato1aqV+v9Xs2bNxNq1azXm+fnnn0WjRo2EUqkU9vb2YsCAAeLvv//WmGfw4MHC0tKy0HUUtT8Knnvx/8yBAwdEkyZNhKmpqahVq5ZYunSp1nd8Qe0v+wydOnVK9O/fX70/HR0dRdeuXUV8fHyJ95FcKITQUe9GIiIioteMfWSIiIhIthhkiIiISLYYZIiIiEi2GGSIiIhIthhkiIiISLYYZIiIiEi2DP6CeCqVCv/88w+sra0lubw2ERERlZ4QAk+fPoWTk1OR97IC3oAg888//+jkPjtERET0+iUlJRV5d3HgDQgyBZfeTkpK0sn9dohKLD0dKLgL7z//AC+50SMREf0rNTUVNWvWLPYWGsAbEGQKDifZ2NgwyNDr9f9v7gYAsLFhkCEiKoOXdQthZ18iIiKSLQYZIiIiki2DP7REJJkKFYDBg/8dJyIineO3K5G+KJVAdLTUVRDJlkqlQnZ2ttRlkJ6YmJjA+Pm+hGXEIENEROVOdnY2rl+/DpVKJXUppEe2traoWrXqK13njUGGSF+EADIy8sctLABekJGoRIQQuHPnDoyNjVGzZs1iL4ZG8iSEQEZGBpKTkwEA1apVK/OyGGSI9CUjA7Cyyh9PS+Pp10QllJubi4yMDDg5OcHCwkLqckhPzM3NAQDJyclwdHQs82EmxlwiIipX8vLyAACmpqYSV0L6VhBUc3JyyrwMBhkiIiqXeH88w6eL95hBhoiIiGRL0iBz8OBBdOvWDU5OTlAoFNi8ebPG80IITJkyBdWqVYO5uTn8/f2RmJgoTbFERESvWXR0NGxtbaUuo8z2798PhUKBJ0+e6G0dkgaZ9PR0NGzYEP/9738Lff6rr77CwoULsXTpUhw/fhyWlpYICAhAZmbma66UiIjo5ZKSkjB06FA4OTnB1NQULi4uGDt2LB4+fCh1aQZL0rOWAgMDERgYWOhzQgjMnz8fX3zxBbp37w4A+PHHH1GlShVs3rwZ77///usslYiIqFjXrl1Dy5YtUbt2baxduxZubm64cOECJkyYgB07duDYsWOwt7cv9LXZ2dnlqnNzeaunOOW2j8z169dx9+5d+Pv7q6dVrFgRzZs3x9GjRyWsjKiEjI2B3r3zBx1cvZKIyreQkBCYmpoiNjYWfn5+cHZ2RmBgIHbv3o3bt2/j888/V8/r6uqK6dOnY9CgQbCxscGIESMA5B9KcnZ2hoWFBXr27FloS86WLVvQuHFjmJmZoVatWoiMjERubq76eYVCge+//x49e/aEhYUFPD09sXXr1mJrL6qeQ4cOoW3btjA3N0fNmjURGhqK9PR09et++ukn+Pr6wtraGlWrVsUHH3ygvjbM61Jug8zdu3cBAFWqVNGYXqVKFfVzhcnKykJqaqrGQCQJMzNgw4b8wcxM6mqI5C89vejhxS4Hxc377FnJ5i2FR48eYefOnRg9erT6+igFqlatigEDBuDnn3+GEEI9fe7cuWjYsCFOnz6NyZMn4/jx4xg2bBjGjBmDM2fOoEOHDpgxY4bGsv744w8MGjQIY8eOxcWLF7Fs2TJER0dj5syZGvNFRkaib9++OHfuHLp06YIBAwbg0aNHxW7Di/VcvXoV77zzDt577z2cO3cOP//8Mw4dOoQxY8aoX5OTk4Pp06fj7Nmz2Lx5M27cuIEhQ4aUat+9KoO7IF5UVBQiIyOlLoOo3HMN3y51CS91Y3aQ1CVQeVJwgcnCdOkCbH/uM+3o+O+VtV/k5wfs3//vY1dX4MED7fmeCx0vk5iYCCEEvLy8Cn3ey8sLjx8/xv379+Ho6AgAePvtt/HJJ5+o55k8eTLeeecdTJw4EQBQu3ZtHDlyBDExMep5IiMjER4ejsH//4a0tWrVwvTp0zFx4kRMnTpVPd+QIUPQv39/AMCsWbOwcOFCnDhxAu+8806R2/BiPcOHD8eAAQMQFhYGAPD09MTChQvh5+eHJUuWwMzMDEOHDlXPX6tWLSxcuBBNmzZFWloarIp7v3So3LbIVK1aFQBw7949jen37t1TP1eYiIgIpKSkqIekpCS91klERFRAlCL8+Pr6ajxOSEhA8+bNNaa1bNlS4/HZs2fx5ZdfwsrKSj18+OGHuHPnDjKeC24NGjRQj1taWsLGxualh3xerOfs2bOIjo7WWFdAQABUKhWuX78OADh58iS6desGZ2dnWFtbw8/PDwBw69atEu6FV1duW2Tc3NxQtWpV7NmzB2+99RYAIDU1FcePH8eoUaOKfJ1SqYRSqXxNVRIVIz2dtygg0qW0tKKfe7EfWnF/tF+8d9ONG2UuqYCHhwcUCgUSEhLQs2dPrecTEhJgZ2cHBwcH9TTLMnwnpKWlITIyEr169dJ6zuy5Q9gmJiYazykUipfegPPFetLS0jBy5EiEhoZqzevs7Iz09HQEBAQgICAAq1evhoODA27duoWAgIDXetdySYNMWloarly5on58/fp1nDlzBvb29nB2dkZYWBhmzJgBT09PuLm5YfLkyXByckKPHj2kK5qIiKRRmj/8+pq3CJUqVUKnTp3w3XffYdy4cRr9ZO7evYvVq1dj0KBBxV7J1svLC8ePH9eYduzYMY3HjRs3xqVLl+Dh4fHKNb9M48aNcfHixSLX9eeff+Lhw4eYPXs2atasCQCIj4/Xe10vkvTQUnx8PBo1aoRGjRoBAMaPH49GjRphypQpAICJEyfi448/xogRI9TH3GJiYjRSJxERUXmwePFiZGVlISAgAAcPHkRSUhJiYmLQqVMnVK9eXatD7otCQ0MRExODuXPnIjExEYsXL9boHwMAU6ZMwY8//ojIyEhcuHABCQkJWLduHb744gudb89nn32GI0eOqDsfJyYmYsuWLerOvs7OzjA1NcWiRYtw7do1bN26FdOnT9d5HS8jaZBp3749hBBaQ3R0NID8prAvv/wSd+/eRWZmJnbv3o3atWtLWTIREVGhPD09ER8fj1q1aqFv375wd3fHiBEj0KFDBxw9erTIa8gUaNGiBVasWIEFCxagYcOGiI2N1QooAQEB2LZtG2JjY9G0aVO0aNEC3377LVxcXHS+PQ0aNMCBAwdw+fJltG3bVt3Q4OTkBABwcHBAdHQ0NmzYgHr16mH27NmYO3euzut4GYUoTc8kGUpNTUXFihWRkpICGxsbqcuhN0k57yPDs5aovMrMzMT169fh5ubGFngDV9x7XdK/3+X2rCUiIiKil2GQISIiItkqt6dfE8mesXH+RboKxomISOcYZIj0xcxM80qjRESkczy0RERE5ZKBn4tC0M17zCBDRETlivH/PxT7Oq8OS9IouK3Ci1ciLg0eWiLSl/T0/BvXAfmXSy9np18TlVcVKlSAhYUF7t+/DxMTExi9eEsBkj0hBDIyMpCcnAxbW1t1eC0LBhkifSrq7rtEVCSFQoFq1arh+vXruHnzptTlkB7Z2toWeyPokmCQISKicsfU1BSenp48vGTATExMXqklpgCDDBERlUtGRka8si+9FA88EhERkWwxyBAREZFsMcgQERGRbLGPDJG+GBkBfn7/jhMRkc4xyBDpi7k5sH+/1FUQERk0/kwkIiIi2WKQISIiItlikCHSl/R0wMEhf0hPl7oaIiKDxD4yRPr04IHUFRARGTS2yBAREZFsMcgQERGRbDHIEBERkWwxyBAREZFsMcgQERGRbPGsJSJ9MTICfH3/HSciIp1jkCHSF3NzIC5O6iqIiAwafyYSERGRbDHIEBERkWwxyBDpS0YG4OqaP2RkSF0NEZFBYh8ZIn0RArh5899xIiLSObbIEBERkWwxyBAREZFsMcgQERGRbDHIEBERkWwxyBAREZFs8awlIn1RKIB69f4dJyIinWOQIdIXCwvgwgWpqyAiMmg8tERERESyxSBDREREssUgQ6QvGRmAt3f+wFsUEBHpBfvIEOmLEMDFi/+OExGRzrFFhoiIiGSLQYaIiIhki0GGiIiIZItBhoiIiGSLQYaIiIhki2ctEemLQgG4uPw7TkREOscgQ6QvFhbAjRtSV0FEZNB4aImIiIhki0GGiIiIZItBhkhfnj0DmjbNH549k7oaIiKDxD4yRPqiUgHx8f+OExGRzrFFhoiIiGSLQYaIiIhki0GGiIiIZItBhoiIiGSLQYaIiIhki2ctEelT5cpSV0BEZNAYZIj0xdISuH9f6iqIiAxauT60lJeXh8mTJ8PNzQ3m5uZwd3fH9OnTIYSQujQiIiIqB8p1i8ycOXOwZMkSrFq1Ct7e3oiPj0dwcDAqVqyI0NBQqcsjIiIiiZXrIHPkyBF0794dQUFBAABXV1esXbsWJ06ckLgyohJ49gwIDMwf37EDMDeXth4iIgNUrg8ttWrVCnv27MHly5cBAGfPnsWhQ4cQWPDHoRBZWVlITU3VGIgkoVIBBw7kD7xFARGRXpTrFpnw8HCkpqaibt26MDY2Rl5eHmbOnIkBAwYU+ZqoqChERka+xiqJiIhIKuW6RWb9+vVYvXo11qxZg1OnTmHVqlWYO3cuVq1aVeRrIiIikJKSoh6SkpJeY8VERET0OpXrFpkJEyYgPDwc77//PgDAx8cHN2/eRFRUFAYPHlzoa5RKJZRK5essk4iIiCRSrltkMjIyYGSkWaKxsTFU7G9AREREKOctMt26dcPMmTPh7OwMb29vnD59GvPmzcPQoUOlLo2IiIjKgXIdZBYtWoTJkydj9OjRSE5OhpOTE0aOHIkpU6ZIXRpRyVhYSF0BEZFBK9dBxtraGvPnz8f8+fOlLoWo9CwtgfR0qasgIjJo5bqPDBEREVFxGGSIiIhIthhkiPQlMxMICsofMjOlroaIyCCV6z4yRLKWlwf8/vu/40REpHNskSEiIiLZYpAhIiIi2WKQISIiItlikCEiIiLZYpAhIiIi2WKQISIiItni6ddE+mJpCQghdRVERAaNLTJEREQkWwwyREREJFsMMkT6kpkJ9OmTP/AWBUREesEgQ6QveXnAxo35A29RQESkFwwyREREJFsMMkRERCRbDDJEREQkWwwyREREJFsMMkRERCRbDDJEREQkW7xFAZG+WFgAaWn/jhMRkc4xyBDpi0KRf78lIiLSGx5aIiIiItlikCHSl6wsYMiQ/CErS+pqiIgMEoMMkb7k5gKrVuUPublSV0NEZJAYZIiIiEi2GGSIiIhIthhkiIiISLYYZIiIiEi2GGSIiIhIthhkiIiISLZ4ZV8ifbGwAJKT/x0nIiKdY5Ah0heFAnBwkLoKIiKDxkNLREREJFsMMkT6kpUFhITkD7xFARGRXjDIEOlLbi7w3Xf5A29RQESkFwwyREREJFsMMkRERCRbDDJEREQkWwwyREREJFsMMkRERCRbDDJEREQkW7yyL5G+mJsD16//O05ERDrHIEOkL0ZGgKur1FUQERk0HloiIiIi2WKQIdKX7GxgwoT8ITtb6mqIiAwSgwyRvuTkAHPn5g85OVJXQ0RkkBhkiIiISLYYZIiIiEi2GGSIiIhIthhkiIiISLYYZIiIiEi2GGSIiIhItnhlXyJ9MTcHzp//d5yIiHSOQYZIX4yMAG9vqasgIjJoPLREREREslXmFpn09HQcOHAAt27dQvYLl18PDQ195cKIZC87G5g1K3980iTA1FTaeoiIDFCZgszp06fRpUsXZGRkID09Hfb29njw4AEsLCzg6OjIIEME5N+WIDIyf3zCBAYZIiI9KNOhpXHjxqFbt254/PgxzM3NcezYMdy8eRNNmjTB3LlzdVrg7du3MXDgQFSqVAnm5ubw8fFBfHy8TtdBRERE8lSmIHPmzBl88sknMDIygrGxMbKyslCzZk189dVXmDRpks6Ke/z4MVq3bg0TExPs2LEDFy9exDfffAM7OzudrYOIiIjkq0yHlkxMTGBklJ+BHB0dcevWLXh5eaFixYpISkrSWXFz5sxBzZo1sXLlSvU0Nzc3nS2fiIiI5K1MLTKNGjVCXFwcAMDPzw9TpkzB6tWrERYWhvr16+usuK1bt8LX1xd9+vSBo6MjGjVqhBUrVuhs+URERCRvZQoys2bNQrVq1QAAM2fOhJ2dHUaNGoX79+9j+fLlOivu2rVrWLJkCTw9PbFz506MGjUKoaGhWLVqVZGvycrKQmpqqsZAREREhkkhhBBSF1EUU1NT+Pr64siRI+ppoaGhiIuLw9GjRwt9zbRp0xBZcKbIc1JSUmBjY6O3Wom0pKcDVlb542lpgKWltPW8wDV8u9QlvNSN2UE6Wc6btK1EhiI1NRUVK1Z86d/vcn1BvGrVqqFevXoa07y8vHDr1q0iXxMREYGUlBT1oMs+O0SlYmYGnDiRP5iZSV0NEZFBKnFn38aNG2PPnj2ws7NDo0aNoFAoipz31KlTOimudevWuHTpksa0y5cvw8XFpcjXKJVKKJVKnayf6JUYGwNNm0pdBRGRQStxkOnevbs6IPTo0UNf9WgYN24cWrVqhVmzZqFv3744ceIEli9frtN+OERERCRfJQ4yU6dOLXRcn5o2bYpNmzYhIiICX375Jdzc3DB//nwMGDDgtayf6JVkZwMLFuSPjx3LK/sSEelBma4jExcXB5VKhebNm2tMP378OIyNjeHr66uT4gCga9eu6Nq1q86WR/Ta5OQAEyfmj48ezSBDRKQHZersGxISUmgn2tu3byMkJOSViyIiIiIqiTIFmYsXL6Jx48Za0xs1aoSLFy++clFEREREJVGmIKNUKnHv3j2t6Xfu3EGFCmU6WkVERERUamUKMp07d1Zfr6XAkydPMGnSJHTq1ElnxREREREVp0zNJ3PnzkW7du3g4uKCRo0aAci/I3aVKlXw008/6bRAIiIioqKUKchUr14d586dw+rVq3H27FmYm5sjODgY/fv3h4mJia5rJCIiIipUmTu0WFpaYsSIEbqshciwmJkB+/b9O05ERDpX5iCTmJiIffv2ITk5GSqVSuO5KVOmvHJhRLJnbAy0by91FUREBq1MQWbFihUYNWoUKleujKpVq2rcd0mhUDDIEBER0WtRpiAzY8YMzJw5E5999pmu6yEyHDk5QMF9wUaMANh/jIhI58oUZB4/fow+ffrouhYiw5KdDYwZkz8+ZAiDDBGRHpTpOjJ9+vRBbGysrmshIiIiKpUytch4eHhg8uTJOHbsGHx8fLROuQ4NDdVJcURERETFKVOQWb58OaysrHDgwAEcOHBA4zmFQsEgQ0RERK9FmYLM9evXdV0HERERUamVqY9MgezsbFy6dAm5ubm6qoeIiIioxMoUZDIyMjBs2DBYWFjA29sbt27dAgB8/PHHmD17tk4LJCIiIipKmYJMREQEzp49i/3798PsuUuv+/v74+eff9ZZcUSyplQC27blD0ql1NUQERmkMvWR2bx5M37++We0aNFC46q+3t7euHr1qs6KI5K1ChWAoCCpqyAiMmhlapG5f/8+HB0dtaanp6drBBsiIiIifSpTkPH19cX27dvVjwvCy/fff4+WLVvqpjIiucvJAaKj84ecHKmrISIySGU6tDRr1iwEBgbi4sWLyM3NxYIFC3Dx4kUcOXJE67oyRG+s7GwgODh/vE8f3qKAiEgPytQi06ZNG5w5cwa5ubnw8fFBbGwsHB0dcfToUTRp0kTXNRIREREVqkwtMgDg7u6OFStW6LIWIiIiolIpU5ApuG5MUZydnctUDBEREVFplCnIuLq6Fnt2Ul5eXpkLIiIiIiqpMgWZ06dPazzOycnB6dOnMW/ePMycOVMnhRERERG9TJmCTMOGDbWm+fr6wsnJCV9//TV69er1yoURERERvUyZO/sWpk6dOoiLi9PlIonkS6kE1q//d5yIiHSuTEEmNTVV47EQAnfu3MG0adPg6empk8KIZK9ChfzrxxARkd6UKcjY2tpqdfYVQqBmzZpYt26dTgojIiIiepkyBZm9e/dqBBkjIyM4ODjAw8MDFSro9GgVkXzl5gKbNuWP9+yZ30JDREQ6VaZv1vbt2+u4DCIDlJUF9O2bP56WxiBDRKQHZbpFQVRUFH744Qet6T/88APmzJnzykURERERlUSZgsyyZctQt25drene3t5YunTpKxdFREREVBJlCjJ3795FtWrVtKY7ODjgzp07r1wUERERUUmUKcjUrFkThw8f1pp++PBhODk5vXJRRERERCVRpt6HH374IcLCwpCTk4O3334bALBnzx5MnDgRn3zyiU4LJCIiIipKmYLMhAkT8PDhQ4wePRrZ2dkAADMzM3z22WeIiIjQaYFERERERSlTkFEoFJgzZw4mT56MhIQEmJubw9PTE0pehp3oX6amwMqV/44TEZHOvdKFLe7evYtHjx6hXbt2UCqVEEJoXfGX6I1lYgIMGSJ1FUREBq1MnX0fPnyIjh07onbt2ujSpYv6TKVhw4axjwwRERG9NmUKMuPGjYOJiQlu3boFCwsL9fR+/fohJiZGZ8URyVpuLrB9e/6Qmyt1NUREBqlMh5ZiY2Oxc+dO1KhRQ2O6p6cnbt68qZPCiGQvKwvo2jV/nLcoICLSizK1yKSnp2u0xBR49OgRO/wSERHRa1OmINO2bVv8+OOP6scKhQIqlQpfffUVOnTooLPiiIiIiIpTprbur776Ch07dkR8fDyys7MxceJEXLhwAY8ePSr0ir9ERERE+lCmFpn69evj8uXLaNOmDbp374709HT06tULp0+fhru7u65rJCIiIipUqVtkcnJy8M4772Dp0qX4/PPP9VETERERUYmUukXGxMQE586d00ctRERERKVSpkNLAwcOxP/93//puhYiw2JqCixenD/wFgVERHpRps6+ubm5+OGHH7B79240adIElpaWGs/PmzdPJ8URyZqJCRASInUVREQGrVRB5tq1a3B1dcX58+fRuHFjAMDly5c15uG9loiIiOh1KVWQ8fT0xJ07d7Bv3z4A+bckWLhwIapUqaKX4ohkLS8P+OOP/PG2bQFjY2nrISIyQKUKMkIIjcc7duxAenq6TgsiMhiZmUDBBSLT0oAXDsESEdGrK1Nn3wIvBhsiIiKi16lUQUahUGj1gWGfGCIiIpJKqQ8tDRkyRH1jyMzMTHz00UdaZy39+uuvuquQiIiIqAilCjKDBw/WeDxw4ECdFkNERERUGqUKMitXrtRXHURERESl9kqdfV+32bNnQ6FQICwsTOpSiIiIqBwo05V9pRAXF4dly5ahQYMGUpdCVDImJsBXX/07TkREOieLFpm0tDQMGDAAK1asgJ2dndTlEJWMqSkwYUL+wHstERHphSyCTEhICIKCguDv7//SebOyspCamqoxEBERkWEq94eW1q1bh1OnTiEuLq5E80dFRSEyMlLPVeVzDd/+WtbzKm7MDpK6BNnR1ftqpMpD/XtXAQDnq7hDZaS7WxTwfSUiyleuW2SSkpIwduxYrF69GmZmZiV6TUREBFJSUtRDUlKSnqskKpwyNwdbfxyPrT+OhzI3R+pyiIgMUrlukTl58iSSk5PVd9oGgLy8PBw8eBCLFy9GVlYWjF+4EZ9SqVRfsI+IiIgMW7kOMh07dsSff/6pMS04OBh169bFZ599phViiIiI6M1SroOMtbU16tevrzHN0tISlSpV0ppOREREb55y3UeGiIiIqDjlukWmMPv375e6BCIiIion2CJDREREsiW7Fhkiucg1Nsb81v3V40REpHsMMkR6kmNsgvltBkhdBhGRQeOhJSIiIpIttsgQ6YlCqODxIP/K0lcq14RQ8HcDEZGuMcgQ6YlZTjZ2/RACAPAatxHPTEt2mw0iIio5/kQkIiIi2WKQISIiItlikCEiIiLZYpAhIiIi2WKQISIiItlikCEiIiLZ4unXRHqSa2yMZc16qceJiEj3GGSI9CTH2ARRHYZKXQYRkUHjoSUiIiKSLbbIEOmJQqhQPfU+AOC2jQNvUUBEpAcMMkR6YpaTjUNLhwHgLQqIiPSFPxGJiIhIthhkiIiISLYYZIiIiEi2GGSIiIhIthhkiIiISLYYZIiIiEi2ePo1kZ7kGRnjx0ZB6nEiItI9BhkiPcmuYIIpnUdJXQYRkUHjoSUiIiKSLbbIEOmLELB/lgoAeGRuAygUEhdERGR4GGSI9MQ8JwunFg0AwFsUEBHpC4MMEZGBcA3fLnUJL3VjdpDUJZCBYR8ZIiIiki0GGSIiIpItBhkiIiKSLQYZIiIiki0GGSIiIpItnrVEpCd5RsbYWL+jepyIiHSPQYZIT7IrmODToHFSl0FEZNB4aImIiIhkiy0yRPoiBMxzsgAAz0yUvEUBEZEesEWGSE/Mc7KQ8G1vJHzbWx1oiIhItxhkiIiISLYYZIiIiEi2GGSIiIhIthhkiIiISLYYZIiIiEi2GGSIiIhItngdGSI9URkZYXud1upxIiLSPQYZIj3JqmCKkB4RUpdBRGTQ+DORiIiIZItBhoiIiGSLQYZIT8yzM3FjTlfcmNMV5tmZUpdDRGSQGGSIiIhIthhkiIiISLYYZIiIiEi2GGSIiIhIthhkiIiISLYYZIiIiEi2eGVfIj1RGRlhby1f9TgREekegwyRnmRVMMXQPtOkLoOIyKDxZyIRERHJVrkOMlFRUWjatCmsra3h6OiIHj164NKlS1KXRUREROVEuQ4yBw4cQEhICI4dO4Zdu3YhJycHnTt3Rnp6utSlEb2UeXYmLs57DxfnvcdbFBAR6Um57iMTExOj8Tg6OhqOjo44efIk2rVrJ1FVRCVnkZMldQlERAatXAeZF6WkpAAA7O3ti5wnKysLWVn//vFITU3Ve11EREQkDdkEGZVKhbCwMLRu3Rr169cvcr6oqChERka+xsreDK7h26Uu4aVuzA6SugQiInrNynUfmeeFhITg/PnzWLduXbHzRUREICUlRT0kJSW9pgqJiIjodZNFi8yYMWOwbds2HDx4EDVq1Ch2XqVSCaVS+ZoqIyIiIimV6yAjhMDHH3+MTZs2Yf/+/XBzc5O6JCIiIipHynWQCQkJwZo1a7BlyxZYW1vj7t27AICKFSvC3Nxc4uqIiqdSKHCsZn31OBER6V65DjJLliwBALRv315j+sqVKzFkyJDXXxBRKWSZKPH+B7OlLoOIyKCV6yAjhJC6BCIiIirHZHPWEhEREdGLGGSI9MQ8OxMnF36Akws/4C0KiIj0pFwfWiKSu0rPeGVpIiJ9YosMERERyRaDDBEREckWgwwRERHJFoMMERERyRaDDBEREckWz1oi0hOVQoGzVT3V40REpHsMMkR6kmWiRPfB30pdBhGRQeOhJSIiIpItBhkiIiKSLQYZIj0xy8nEoSVDcWjJUJjl8BYFRET6wD4yRHqiEECN1GT1OBER6R5bZIiIiEi2GGSIiIhIthhkiIiISLYYZIiIiEi2GGSIiIhItnjWEpGeCAVwuZKzepyIiHSPQYZITzJNzNB5+HdSl0FEZNB4aImIiIhki0GGiIiIZItBhkhPzHIyEfv9aMR+P5q3KCAi0hP2kSHSE4UAaj+8pR4nIiLdY4sMERERyRaDDBEREckWgwwRERHJFoMMERERyRY7+xIRkey4hm+XuoSXujE7SOoS3ggMMkR6IhTA3zaO6nEiItI9BhkiPck0MUObUT9IXQYRkUFjHxkiIiKSLQYZIiIiki0GGSI9UeZkYcuqcdiyahyUOVlSl0NEZJDYR4ZIT4yEQMO7iepxIiLSPbbIEBERkWwxyBAREZFsMcgQERGRbDHIEBERkWwxyBAREZFs8awlIj16aG4jdQlERAaNQYZIT56ZmqFJ6BqpyyAiMmg8tERERESyxSBDREREssUgQ6QnypwsrFsTjnVrwnmLAiIiPWEfGSI9MRICLZLOq8eJiEj32CJDREREssUgQ0RERLLFIENERESyxSBDREREssUgQ0RERLLFs5aI9CjDRCl1CUREBo1BhkhPnpmaod74X6Qug4jIoPHQEhEREckWgwwRERHJFg8tEemJMjcbSzbNAgCM6jkJWRVMJa6IiMjwMMgQ6YmRSoW3r8Wrx4mISPd4aImIiIhkSxZB5r///S9cXV1hZmaG5s2b48SJE1KXREREROVAuQ8yP//8M8aPH4+pU6fi1KlTaNiwIQICApCcnCx1aURERCSxch9k5s2bhw8//BDBwcGoV68eli5dCgsLC/zwww9Sl0ZEREQSK9dBJjs7GydPnoS/v796mpGREfz9/XH06FEJKyMiIqLyoFyftfTgwQPk5eWhSpUqGtOrVKmCv/76q9DXZGVlISsrS/04JSUFAJCamqrz+lRZGTpfpq7paru5raWXl52JgorysjKgEro7c0kX28r3tHzhe1o6b9K2vqkK9p8QovgZRTl2+/ZtAUAcOXJEY/qECRNEs2bNCn3N1KlTBQAOHDhw4MCBgwEMSUlJxWaFct0iU7lyZRgbG+PevXsa0+/du4eqVasW+pqIiAiMHz9e/VilUuHRo0eoVKkSFAqFXut9VampqahZsyaSkpJgY2MjdTl69aZs65uyncCbs61vynYCb862vinbCchrW4UQePr0KZycnIqdr1wHGVNTUzRp0gR79uxBjx49AOQHkz179mDMmDGFvkapVEKp1LzjsK2trZ4r1S0bG5ty/wHTlTdlW9+U7QTenG19U7YTeHO29U3ZTkA+21qxYsWXzlOugwwAjB8/HoMHD4avry+aNWuG+fPnIz09HcHBwVKXRkRERBIr90GmX79+uH//PqZMmYK7d+/irbfeQkxMjFYHYCIiInrzlPsgAwBjxowp8lCSIVEqlZg6darWoTFD9KZs65uyncCbs61vynYCb862vinbCRjmtiqEeNl5TURERETlU7m+IB4RERFRcRhkiIiISLYYZIiIiEi2GGSIiIhIthhkyoGDBw+iW7ducHJygkKhwObNm6UuSS+ioqLQtGlTWFtbw9HRET169MClS5ekLksvlixZggYNGqgvOtWyZUvs2LFD6rL0bvbs2VAoFAgLC5O6FJ2bNm0aFAqFxlC3bl2py9KL27dvY+DAgahUqRLMzc3h4+OD+Ph4qcvSOVdXV633VKFQICQkROrSdCovLw+TJ0+Gm5sbzM3N4e7ujunTp7/8HkYyIYvTrw1deno6GjZsiKFDh6JXr15Sl6M3Bw4cQEhICJo2bYrc3FxMmjQJnTt3xsWLF2FpaSl1eTpVo0YNzJ49G56enhBCYNWqVejevTtOnz4Nb29vqcvTi7i4OCxbtgwNGjSQuhS98fb2xu7du9WPK1QwvK/Qx48fo3Xr1ujQoQN27NgBBwcHJCYmws7OTurSdC4uLg55eXnqx+fPn0enTp3Qp08fCavSvTlz5mDJkiVYtWoVvL29ER8fj+DgYFSsWBGhoaFSl/fKDO9/oQwFBgYiMDBQ6jL0LiYmRuNxdHQ0HB0dcfLkSbRr106iqvSjW7duGo9nzpyJJUuW4NixYwYZZNLS0jBgwACsWLECM2bMkLocvalQoUKR93kzFHPmzEHNmjWxcuVK9TQ3NzcJK9IfBwcHjcezZ8+Gu7s7/Pz8JKpIP44cOYLu3bsjKCgIQH5L1Nq1a3HixAmJK9MNHloiyaSkpAAA7O3tJa5Ev/Ly8rBu3Tqkp6ejZcuWUpejFyEhIQgKCoK/v7/UpehVYmIinJycUKtWLQwYMAC3bt2SuiSd27p1K3x9fdGnTx84OjqiUaNGWLFihdRl6V12djb+97//YejQoeX+BsOl1apVK+zZsweXL18GAJw9exaHDh0ymB/QbJEhSahUKoSFhaF169aoX7++1OXoxZ9//omWLVsiMzMTVlZW2LRpE+rVqyd1WTq3bt06nDp1CnFxcVKXolfNmzdHdHQ06tSpgzt37iAyMhJt27bF+fPnYW1tLXV5OnPt2jUsWbIE48ePx6RJkxAXF4fQ0FCYmppi8ODBUpenN5s3b8aTJ08wZMgQqUvRufDwcKSmpqJu3bowNjZGXl4eZs6ciQEDBkhdmk4wyJAkQkJCcP78eRw6dEjqUvSmTp06OHPmDFJSUrBx40YMHjwYBw4cMKgwk5SUhLFjx2LXrl0wMzOTuhy9ev7Xa4MGDdC8eXO4uLhg/fr1GDZsmISV6ZZKpYKvry9mzZoFAGjUqBHOnz+PpUuXGnSQ+b//+z8EBgbCyclJ6lJ0bv369Vi9ejXWrFkDb29vnDlzBmFhYXBycjKI95RBhl67MWPGYNu2bTh48CBq1KghdTl6Y2pqCg8PDwBAkyZNEBcXhwULFmDZsmUSV6Y7J0+eRHJyMho3bqyelpeXh4MHD2Lx4sXIysqCsbGxhBXqj62tLWrXro0rV65IXYpOVatWTStse3l54ZdffpGoIv27efMmdu/ejV9//VXqUvRiwoQJCA8Px/vvvw8A8PHxwc2bNxEVFcUgQ1QaQgh8/PHH2LRpE/bv32+wHQiLolKpkJWVJXUZOtWxY0f8+eefGtOCg4NRt25dfPbZZwYbYoD8Ds5Xr17Ff/7zH6lL0anWrVtrXRbh8uXLcHFxkagi/Vu5ciUcHR3VnWENTUZGBoyMNLvEGhsbQ6VSSVSRbjHIlANpaWkav+quX7+OM2fOwN7eHs7OzhJWplshISFYs2YNtmzZAmtra9y9excAULFiRZibm0tcnW5FREQgMDAQzs7OePr0KdasWYP9+/dj586dUpemU9bW1lp9nCwtLVGpUiWD6/v06aefolu3bnBxccE///yDqVOnwtjYGP3795e6NJ0aN24cWrVqhVmzZqFv3744ceIEli9fjuXLl0tdml6oVCqsXLkSgwcPNsjT6YH8syhnzpwJZ2dneHt74/Tp05g3bx6GDh0qdWm6IUhy+/btEwC0hsGDB0tdmk4Vto0AxMqVK6UuTeeGDh0qXFxchKmpqXBwcBAdO3YUsbGxUpf1Wvj5+YmxY8dKXYbO9evXT1SrVk2YmpqK6tWri379+okrV65IXZZe/Pbbb6J+/fpCqVSKunXriuXLl0tdkt7s3LlTABCXLl2SuhS9SU1NFWPHjhXOzs7CzMxM1KpVS3z++eciKytL6tJ0QiGEgVzaj4iIiN44vI4MERERyRaDDBEREckWgwwRERHJFoMMERERyRaDDBEREckWgwwRERHJFoMMERERyRaDDBGVSXR0NGxtbaUuo0zat2+PsLCwN27dRIaIQYboDZGUlIShQ4fCyckJpqamcHFxwdixY/Hw4UOpS9Ob/fv3Q6FQ4MmTJxrTf/31V0yfPl2aoohIpxhkiN4A165dg6+vLxITE7F27VpcuXIFS5cuxZ49e9CyZUs8evSoyNdmZ2e/xkpfD3t7e1hbW0tdBhHpAIMM0RsgJCQEpqamiI2NhZ+fH5ydnREYGIjdu3fj9u3b+Pzzz9Xzurq6Yvr06Rg0aBBsbGwwYsQIAPmHkpydnWFhYYGePXsW2pKzZcsWNG7cGGZmZqhVqxYiIyORm5urfl6hUOD7779Hz549YWFhAU9PT2zdurXY2rOysvDpp5+ievXqsLS0RPPmzbF//3718zdv3kS3bt1gZ2cHS0tLeHt74/fff8eNGzfQoUMHAICdnR0UCgWGDBkCQPvwjqurK2bMmIFBgwbBysoKLi4u2Lp1K+7fv4/u3bvDysoKDRo0QHx8vPo1Dx8+RP/+/VG9enVYWFjAx8cHa9eu1ag9PT1dvcxq1arhm2++0dq+x48fY9CgQbCzs4OFhQUCAwORmJhY7D4houdIfbMnItKvhw8fCoVCIWbNmlXo8x9++KGws7MTKpVKCCGEi4uLsLGxEXPnzhVXrlwRV65cEceOHRNGRkZizpw54tKlS2LBggXC1tZWVKxYUb2cgwcPChsbGxEdHS2uXr0qYmNjhaurq5g2bZp6HgCiRo0aYs2aNSIxMVGEhoYKKysr8fDhwyLrHz58uGjVqpU4ePCguHLlivj666+FUqkUly9fFkIIERQUJDp16iTOnTsnrl69Kn777Tdx4MABkZubK3755Rf1DQHv3Lkjnjx5IoTQvrGli4uLsLe3F0uXLhWXL18Wo0aNEjY2NuKdd94R69evF5cuXRI9evQQXl5e6v30999/i6+//lqcPn1aXL16VSxcuFAYGxuL48ePq5c7atQo4ezsLHbv3i3OnTsnunbtKqytrTXW/e677wovLy9x8OBBcebMGREQECA8PDxEdnZ2yd5gojccgwyRgTt27JgAIDZt2lTo8/PmzRMAxL1794QQ+X/Ue/TooTFP//79RZcuXTSm9evXTyPIdOzYUSss/fTTT6JatWrqxwDEF198oX6clpYmAIgdO3YUWtvNmzeFsbGxuH37tsb0jh07ioiICCGEED4+Phph6XkFd5Z//PixxvTCgszAgQPVj+/cuSMAiMmTJ6unHT16VAAQd+7cKXRdQuSHqk8++UQIIcTTp0+FqampWL9+vfr5hw8fCnNzc/W6L1++LACIw4cPq+d58OCBMDc313gdERWtglQtQUT0eolS3Oje19dX43FCQgJ69uypMa1ly5aIiYlRPz579iwOHz6MmTNnqqfl5eUhMzMTGRkZsLCwAAA0aNBA/bylpSVsbGyQnJxcaB1//vkn8vLyULt2bY3pWVlZqFSpEgAgNDQUo0aNQmxsLPz9/fHee+9prKOknn9NlSpVAAA+Pj5a05KTk1G1alXk5eVh1qxZWL9+PW7fvo3s7GxkZWWpt/Pq1avIzs5G8+bN1cuwt7dHnTp11I8TEhJQoUIFjXkqVaqEOnXqICEhodTbQPQmYpAhMnAeHh5QKBSFhhEg/4+pnZ0dHBwc1NMsLS1LvZ60tDRERkaiV69eWs+ZmZmpx01MTDSeUygUUKlURS7T2NgYJ0+ehLGxscZzVlZWAIDhw4cjICAA27dvR2xsLKKiovDNN9/g448/LlX9z9elUCiKnFZQ69dff40FCxZg/vz58PHxgaWlJcLCwgyyczRRecbOvkQGrlKlSujUqRO+++47PHv2TOO5u3fvYvXq1ejXr5/6D3VhvLy8cPz4cY1px44d03jcuHFjXLp0CR4eHlqDkVHZvmoaNWqEvLw8JCcnay2zatWq6vlq1qyJjz76CL/++is++eQTrFixAgBgamoKIL9lSNcOHz6M7t27Y+DAgWjYsCFq1aqFy5cvq593d3eHiYmJxn57/PixxjxeXl7Izc3VmOfhw4e4dOkS6tWrp/OaiQwRgwzRG2Dx4sXIyspCQEAADh48iKSkJMTExKBTp06oXr26xuGgwoSGhiImJgZz585FYmIiFi9erHFYCQCmTJmCH3/8EZGRkbhw4QISEhKwbt06fPHFF2Wuu3bt2hgwYAAGDRqEX3/9FdevX8eJEycQFRWF7du3AwDCwsKwc+dOXL9+HadOncK+ffvg5eUFAHBxcYFCocC2bdtw//59pKWllbmWF3l6emLXrl04cuQIEhISMHLkSNy7d0/9vJWVFYYNG4YJEyZg7969OH/+PIYMGaIR6jw9PdG9e3d8+OGHOHToEM6ePYuBAweievXq6N69u85qJTJkDDJEbwBPT0/Ex8ejVq1a6Nu3L9zd3TFixAh06NABR48ehb29fbGvb9GiBVasWIEFCxagYcOGiI2N1QooAQEB2LZtG2JjY9G0aVO0aNEC3377LVxcXF6p9pUrV2LQoEH45JNPUKdOHfTo0QNxcXFwdnYGkN/aEhISAi8vL7zzzjuoXbs2vvvuOwBA9erVERkZifDwcFSpUgVjxox5pVqe98UXX6Bx48YICAhA+/btUbVqVfTo0UNjnq+//hpt27ZFt27d4O/vjzZt2qBJkyZa29ekSRN07doVLVu2hBACv//+u9YhOCIqnEKUpgcgERERUTnCFhkiIiKSLQYZIiIiki0GGSIiIpItBhkiIiKSLQYZIiIiki0GGSIiIpItBhkiIiKSLQYZIiIiki0GGSIiIpItBhkiIiKSLQYZIiIiki0GGSIiIpKt/wfacZjhvfEogAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Parámetros de la simulación\n",
    "iteraciones = 30\n",
    "N = 500\n",
    "orden_real = 3\n",
    "varianza_ruido = 0.001\n",
    "orden_max = 8\n",
    "C = 5.0\n",
    "\n",
    "ordenes_estimados = []\n",
    "aciertos = 0\n",
    "\n",
    "for _ in range(iteraciones):\n",
    "    # Generar señal sintética con 'orden_real' cosenos\n",
    "    x, _, _ = generar_muestra_coseno(N=N, num_cosenos=orden_real, varianza_ruido=varianza_ruido)\n",
    "    \n",
    "    # Estimar el orden\n",
    "    orden_est, _, _ = estimador_de_orden(x, orden_max=orden_max, C=C)\n",
    "    ordenes_estimados.append(orden_est)\n",
    "    \n",
    "    if orden_est == orden_real:\n",
    "        aciertos += 1\n",
    "\n",
    "# Resultados\n",
    "print(\"Orden real:\", orden_real)\n",
    "print(\"Tasa de acierto:\", aciertos / iteraciones)\n",
    "print(\"Órdenes estimados:\", ordenes_estimados)\n",
    "\n",
    "# Histograma de órdenes estimados\n",
    "plt.hist(ordenes_estimados, bins=range(1, orden_max+2), align=\"left\", rwidth=0.7)\n",
    "plt.axvline(orden_real, color=\"red\", linestyle=\"--\", label=\"Orden real\")\n",
    "plt.xlabel(\"Orden estimado\")\n",
    "plt.ylabel(\"Frecuencia\")\n",
    "plt.title(f\"Distribución de órdenes estimados en {iteraciones} corridas\")\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}