You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
520 lines
48 KiB
Plaintext
520 lines
48 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Generacion de muestras"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"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": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"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": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"\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": 5,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Orden real: 3\n",
|
|
"Orden estimado: 2\n",
|
|
"Varianzas: [2.323137726399167, 1.8383255140487047, 1.837929702038553, 1.8368502482272016, 1.8355685769211567, 1.8287296226777754, inf, inf]\n",
|
|
"Criterios: [np.float64(426.95937009636503), np.float64(315.4275552824338), np.float64(320.8198881079817), np.float64(326.02614158247326), np.float64(331.177142348153), np.float64(334.8107652653657), inf, inf]\n",
|
|
"Parámetros óptimos: [{'amplitud': np.float64(0.07807501614245879), 'frecuencia': np.float64(0.053583028276843), 'fase': np.float64(1.4717432359599454)}, {'amplitud': np.float64(0.9879446102229781), 'frecuencia': np.float64(9.435587905228454), 'fase': np.float64(0.37640582199551625)}]\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": 6,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Orden real: 3\n",
|
|
"Tasa de acierto: 0.16666666666666666\n",
|
|
"Órdenes estimados: [4, 6, 7, 5, 3, 6, 3, 1, 5, 3, 5, 2, 4, 1, 4, 1, 6, 1, 4, 3, 6, 5, 4, 3, 8, 4, 6, 7, 4, 8]\n"
|
|
]
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAHHCAYAAACRAnNyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABQl0lEQVR4nO3deXxM5/4H8M9km+wbCUEWkiBiKYl9CRXSWGqppcq1VSnRWLqJtkiV0LqW1q2tt+iCopRaInbX2sRWW0liS20hSGSRZeb5/eGXU2MSYsw4x/i8X6/zcpZnzvmemTH5zrONSgghQERERKRAFnIHQERERFQaJipERESkWExUiIiISLGYqBAREZFiMVEhIiIixWKiQkRERIrFRIWIiIgUi4kKERERKRYTFXrh5OfnY+rUqdiyZYvcoRjF7du3ERsbi0OHDskdChGR4jBRUbBJkyZBpVI9l2u1bt0arVu3lrZ37doFlUqF1atXP5frP0ylUmHSpEmlHh87dix+/vlnNG7c+LnEM3DgQPj5+Znk3EII9O/fH7t27UL9+vWf+Xx+fn4YOHDgswf2Enj0Pa8EfP3Mw5M+w4o9z8/4FxkTledkyZIlUKlU0mJra4tKlSohIiICX3/9Ne7du2eU61y9ehWTJk3CsWPHjHI+pVm5ciV+++03bN68Ga6urnKH88y+/PJLXLx4EWvXroWNjY3c4Zid06dPY9KkSbh48aLcobzU8vLy8Pbbb6N27dpwcXGBo6Mj6tWrhzlz5qCwsFCv/N27dzF06FB4eHjAwcEBbdq0wZEjR2SInJTASu4AXjaff/45qlatisLCQly/fh27du3C6NGjMXPmTKxfvx5169aVyn766acYN27cU53/6tWriI2NhZ+fH1555ZUyPy4hIeGprmNKeXl5sLLSf2sKIfD3339j8+bN8PHxkSEy47p//z6KioqwadMms0i6lOj06dOIjY1F69at9WrFlPSeN3d5eXk4deoUOnToAD8/P1hYWGD//v0YM2YMDh06hGXLlklltVotOnbsiOPHj+PDDz9E+fLl8e2336J169Y4fPgwAgMDZbyTsintM4wMw2fyOYuMjERoaKi0HRMTgx07dqBTp054/fXXcebMGdjZ2QEArKysTP5mz83Nhb29vaK+zdva2pa4X6VSYezYsc85GtOxtbXFJ5988lSPycnJgYODg4kierko6T1v7tzd3XHw4EGdfe+++y5cXFwwd+5czJw5ExUrVgQArF69Gvv378eqVavQo0cPAECvXr1QvXp1TJw4USepkUNp/we1Wi0KCgpga2tb6mcYGYZNPwrw6quv4rPPPsOlS5fw008/SftLar/cunUrWrRoAVdXVzg6OqJGjRoYP348gAf9Sho2bAgAGDRokNTMtGTJEgAP2uRr166Nw4cPo1WrVrC3t5ceW1p7vUajwfjx41GxYkU4ODjg9ddfR1pamk6Z0trVSzrn/fv3MWnSJFSvXh22trbw8vJC9+7dkZqaKpUpqX336NGjiIyMhLOzMxwdHdG2bVu9D77i5rV9+/Zh7NixUrVxt27dcPPmTb34SvLbb7+hdu3asLW1Re3atbF27doSy2m1WsyePRvBwcGwtbVFhQoVMGzYMNy5c6dM19mxYwdatmwJBwcHuLq6okuXLjhz5oxOmeLX//Tp03jrrbfg5uaGFi1aAHhQu/TFF1+gSpUqsLe3R5s2bXDq1KkSr3X37l2MHj0a3t7eUKvVCAgIwPTp06HVaqUyFy9ehEqlwowZM7Bw4UL4+/tDrVajYcOGSExM1DvnX3/9hR49esDd3R22trYIDQ3F+vXrdcoUFhYiNjYWgYGBsLW1Rbly5dCiRQts3br1ic9PWWIGgBUrViAkJAROTk5wdnZGnTp1MGfOHAAP3g89e/YEALRp00b6/7Br1y4ApffLWrlyJWJjY1G5cmU4OTmhR48eyMzMRH5+PkaPHg1PT084Ojpi0KBByM/P14ln8eLFePXVV+Hp6Qm1Wo1atWph3rx5evf3NK/f+fPn0bNnT7i7u8Pe3h5NmjTBxo0b9cp98803CA4Ohr29Pdzc3BAaGlqmP+r5+fmYOHEiAgICoFar4e3tjY8++kjv3lQqFUaOHCn9H1Gr1QgODkZ8fPwTr1Ga4lquu3fvSvtWr16NChUqoHv37tI+Dw8P9OrVC+vWrdOLqySbN29GWFiY9L5o2LCh3nOxatUqhISEwM7ODuXLl0e/fv1w5coVnTIDBw6Eo6MjUlNT0aFDBzg5OaFv374A/nk+fv75ZwQHB0OtVkvPRUmfYXv37kXDhg1ha2sLf39/LFiwoMTYy/oeSkpKQkREBMqXLw87OztUrVoVgwcPfuJz8yJijYpC/Otf/8L48eORkJCAd955p8Qyp06dQqdOnVC3bl18/vnnUKvVSElJwb59+wAAQUFB+PzzzzFhwgQMHToULVu2BAA0a9ZMOkdGRgYiIyPx5ptvol+/fqhQocJj45oyZQpUKhU+/vhjpKenY/bs2QgPD8exY8ekmp+y0mg06NSpE7Zv344333wTo0aNwr1797B161acPHkS/v7+pd53y5Yt4ezsjI8++gjW1tZYsGABWrdujd27d+t1qn3vvffg5uaGiRMn4uLFi5g9ezZGjhyJX3755bHxJSQk4I033kCtWrUQFxeHjIwMDBo0CFWqVNErO2zYMCxZsgSDBg1CdHQ0Lly4gLlz5+Lo0aPYt28frK2tS73Otm3bEBkZiWrVqmHSpEnIy8vDN998g+bNm+PIkSN6TRQ9e/ZEYGAgpk6dCiEEAGDChAn44osv0KFDB3To0AFHjhxB+/btUVBQoPPY3NxchIWF4cqVKxg2bBh8fHywf/9+xMTE4Nq1a5g9e7ZO+WXLluHevXsYNmwYVCoVvvzyS3Tv3h3nz5+X7unUqVNo3rw5KleujHHjxsHBwQErV65E165d8euvv6Jbt24AHiRacXFxGDJkCBo1aoSsrCwkJSXhyJEjaNeuXanPT1lj3rp1K/r06YO2bdti+vTpAIAzZ85g3759GDVqFFq1aoXo6Gh8/fXXGD9+PIKCggBA+rc0cXFxsLOzw7hx45CSkoJvvvkG1tbWsLCwwJ07dzBp0iQcPHgQS5YsQdWqVTFhwgTpsfPmzUNwcDBef/11WFlZ4ffff8eIESOg1WoRFRUllSvr63fjxg00a9YMubm5iI6ORrly5bB06VK8/vrrWL16tfRcL1q0CNHR0ejRowdGjRqF+/fv488//8ShQ4fw1ltvlXqvWq0Wr7/+Ovbu3YuhQ4ciKCgIJ06cwKxZs3Du3Dn89ttvOuX37t2LNWvWYMSIEXBycsLXX3+NN954A5cvX0a5cuUe+7wCQEFBAbKyspCXl4ekpCTMmDEDvr6+CAgIkMocPXoUDRo0gIWF7vfoRo0aYeHChTh37hzq1KlT6jWWLFmCwYMHIzg4GDExMXB1dcXRo0cRHx8vPRfF/3cbNmyIuLg43LhxA3PmzMG+fftw9OhRnabYoqIiREREoEWLFpgxYwbs7e2lYzt27MDKlSsxcuRIlC9fvtRO9ydOnED79u3h4eGBSZMmoaioCBMnTizx87cs76H09HTpfOPGjYOrqysuXryINWvWPPE1eCEJei4WL14sAIjExMRSy7i4uIj69etL2xMnThQPv0SzZs0SAMTNmzdLPUdiYqIAIBYvXqx3LCwsTAAQ8+fPL/FYWFiYtL1z504BQFSuXFlkZWVJ+1euXCkAiDlz5kj7fH19xYABA554zu+//14AEDNnztQrq9VqpXUAYuLEidJ2165dhY2NjUhNTZX2Xb16VTg5OYlWrVpJ+4qf4/DwcJ3zjRkzRlhaWoq7d+/qXfdhr7zyivDy8tIpl5CQIAAIX19fad///vc/AUD8/PPPOo+Pj48vcX9J1/H09BQZGRnSvuPHjwsLCwvRv39/aV/x69+nTx+dx6enpwsbGxvRsWNHnfscP368AKDzWkyePFk4ODiIc+fO6Zxj3LhxwtLSUly+fFkIIcSFCxcEAFGuXDlx+/Ztqdy6desEAPH7779L+9q2bSvq1Kkj7t+/L+3TarWiWbNmIjAwUNpXr1490bFjx8c+FyUpa8yjRo0Szs7OoqioqNRzrVq1SgAQO3fu1DtW2nu+du3aoqCgQNrfp08foVKpRGRkpM7jmzZtqvO+EEKI3NxcvetERESIatWqSdtP8/qNHj1aABD/+9//pH337t0TVatWFX5+fkKj0QghhOjSpYsIDg4u9XkozY8//igsLCx0zi+EEPPnzxcAxL59+6R9AISNjY1ISUmR9h0/flwAEN98802Zrrd8+XIBQFpCQ0PFn3/+qVPGwcFBDB48WO+xGzduFABEfHx8qee/e/eucHJyEo0bNxZ5eXk6x4qf64KCAuHp6Slq166tU2bDhg0CgJgwYYK0b8CAAQKAGDdunN61AAgLCwtx6tSpEo89+hlma2srLl26JO07ffq0sLS01PmMF6Js76G1a9c+8e+JOWHTj4I4Ojo+dvRPcZa/bt06vSrwslKr1Rg0aFCZy/fv3x9OTk7Sdo8ePeDl5YVNmzY99bV//fVXlC9fHu+9957esdKG6Gk0GiQkJKBr166oVq2atN/LywtvvfUW9u7di6ysLJ3HDB06VOd8LVu2hEajwaVLl0qN7dq1azh27BgGDBgAFxcXaX+7du1Qq1YtnbKrVq2Ci4sL2rVrh1u3bklLSEgIHB0dsXPnzideZ+DAgXB3d5f2161bF+3atSvxeX333Xd1trdt24aCggK89957Ovc5evRovceuWrUKLVu2hJubm06s4eHh0Gg02LNnj0753r17w83NTdourpU7f/48gAdzvuzYsQO9evXCvXv3pPNlZGQgIiICycnJUvW5q6srTp06heTk5FKfj5KUNWZXV1fk5OSUqSnpafTv31+nRqxx48YQQuhVqzdu3BhpaWkoKiqS9j1cy5iZmYlbt24hLCwM58+fR2ZmJoCne/02bdqERo0aSU1+wIPPiaFDh+LixYs4ffo0gAfPxd9//11iM93jrFq1CkFBQahZs6bOc/3qq68CgN57OTw8XKfms27dunB2dpbeH0/Spk0bbN26FatWrcK7774La2tr5OTk6JTJy8uDWq3We2xxv4+8vLxSz79161bcu3cP48aN0+snUvxcJyUlIT09HSNGjNAp07FjR9SsWbPEZrXhw4eXeL2wsDC9z4dHaTQabNmyBV27dtUZBBAUFISIiAi98mV5DxX/LdiwYUOJo6bMDRMVBcnOztZJCh7Vu3dvNG/eHEOGDEGFChXw5ptvYuXKlU+VtFSuXPmpOhE+2sNepVIhICDAoOGeqampqFGjxlN1EL558yZyc3NRo0YNvWNBQUHQarV6fWYeHRFU/If3cf1HipOYkkYUPHrt5ORkZGZmwtPTEx4eHjpLdnY20tPTn3id0u7n1q1beh/cVatWLVOsHh4eOklGcazx8fF6cYaHhwOAXqxPeu5SUlIghMBnn32md86JEyfqnPPzzz/H3bt3Ub16ddSpUwcffvgh/vzzz1Kfm6eNecSIEahevToiIyNRpUoVDB48+Jn6S5T2HBQnrt7e3nr7tVqt9McDAPbt24fw8HCp75GHh4fUD6y43NO8fpcuXSr1vfLwuT7++GM4OjqiUaNGCAwMRFRUlNQk/DjJyck4deqU3nNdvXp1AE9+fwAP3iNl7ZtVoUIFhIeHo0ePHpg3bx46deqEdu3a4fr161IZOzu7Evuh3L9/XzpemuK+brVr1y61zOP+D9asWVPvC42VlVWJzb+A/v/Nkty8eRN5eXll+mwByvYeCgsLwxtvvIHY2FiUL18eXbp0weLFi8vUf+dFxD4qCvH3338jMzNTp632UXZ2dtizZw927tyJjRs3Ij4+Hr/88gteffVVJCQkwNLS8onXedp+JWXxuNqQssRkbKVdU/x//45npdVq4enpiZ9//rnE4x4eHka5TrFnec20Wi3atWuHjz76qMTjxX+Qij3puStOij/44IMSvw0CkN7DrVq1QmpqKtatW4eEhAR89913mDVrFubPn48hQ4Y8c8yenp44duwYtmzZgs2bN2Pz5s1YvHgx+vfvj6VLl5Z6/icp7Tl40nOTmpqKtm3bombNmpg5cya8vb1hY2ODTZs2YdasWQbXgpZFUFAQzp49iw0bNiA+Ph6//vorvv32W0yYMAGxsbGlPk6r1aJOnTqYOXNmiccfTc6M/X+rR48e+OSTT7Bu3ToMGzYMwIPa0mvXrumVLd5XqVIlg65lKLVarddfppixP0/L+h4qnozz4MGD+P3337FlyxYMHjwY//73v3Hw4EE4OjoaNS65MVFRiB9//BEASv3wL2ZhYYG2bduibdu2mDlzJqZOnYpPPvkEO3fuRHh4uNFnOXy02l4IgZSUFJ35Xtzc3HR67Re7dOmSTnONv78/Dh06hMLCwsd2Nn2Yh4cH7O3tcfbsWb1jf/31FywsLPQ+TA3h6+sLQP9+Aehd29/fH9u2bUPz5s2f+oOq+Dql3U/58uWfOPz44Vgffn5v3ryp983W398f2dnZUm3Esyq+nrW1dZnO6e7ujkGDBmHQoEHIzs5Gq1atMGnSpMcmKk8Ts42NDTp37ozOnTtDq9VixIgRWLBgAT777DMEBAQ811k/f//9d+Tn52P9+vU6NQ+PNp88zevn6+tb6nvl4XMBgIODA3r37o3evXujoKAA3bt3x5QpUxATE1PqcFl/f38cP34cbdu2lWWG1OJmnIdrpV555RX873//g1ar1UkQDh06BHt7e73k+mHFzVInT54s9Uvfw/8Hi5u4ip09e1bnOTUGDw8P2NnZlemzpazvoWJNmjRBkyZNMGXKFCxbtgx9+/bFihUrHvv/60XEph8F2LFjByZPnoyqVatKQ99Kcvv2bb19xZO6FVf5Ff+RKylxMMQPP/yg029m9erVuHbtGiIjI6V9/v7+OHjwoM6IhQ0bNug1ybzxxhu4desW5s6dq3ed0r6RWVpaon379li3bp1Oc9ONGzewbNkytGjRAs7OzobensTLywuvvPIKli5dqvOhuXXrVqkfQLFevXpBo9Fg8uTJeucpKip67HP/8HUeLnfy5EkkJCSgQ4cOT4w1PDwc1tbW+Oabb3Set0dH8BTHeuDAgRJ/F+nu3bs6/SvKwtPTE61bt8aCBQtK/Nb78DDwjIwMnWOOjo4ICAh4YvV0WWN+9PwWFhZSAm2q/w+PU1zb8PBrkpmZicWLF+uUe5rXr0OHDvjjjz9w4MABaV9OTg4WLlwIPz8/qX/Eo8+FjY0NatWqBSHEY/sw9OrVC1euXMGiRYv0juXl5ek1Qxrq1q1bJf4f/+677wBAZ26pHj164MaNGzojWG7duoVVq1ahc+fOJfZfKda+fXs4OTkhLi5OaioqVnz90NBQeHp6Yv78+Trvxc2bN+PMmTPo2LGjYTdZCktLS0REROC3337D5cuXpf1nzpzRe4+X9T10584dvefz0b8F5oQ1Ks/Z5s2b8ddff6GoqAg3btzAjh07sHXrVvj6+mL9+vWPnSjo888/x549e9CxY0f4+voiPT0d3377LapUqSJ1tvP394erqyvmz58PJycnODg4oHHjxmVqSy2Ju7s7WrRogUGDBuHGjRuYPXs2AgICdIZQDxkyBKtXr8Zrr72GXr16ITU1FT/99JPecOP+/fvjhx9+wNixY/HHH3+gZcuWyMnJwbZt2zBixAh06dKlxBi++OILaf6YESNGwMrKCgsWLEB+fj6+/PJLg+6rJHFxcejYsSNatGiBwYMH4/bt29LcFNnZ2VK5sLAwDBs2DHFxcTh27Bjat28Pa2trJCcnY9WqVZgzZ440UVVJvvrqK0RGRqJp06Z4++23peHJLi4uZfp9EA8PD3zwwQeIi4tDp06d0KFDBxw9ehSbN29G+fLldcp++OGHWL9+PTp16oSBAwciJCQEOTk5OHHiBFavXo2LFy/qPeZJ/vOf/6BFixaoU6cO3nnnHVSrVg03btzAgQMH8Pfff+P48eMAgFq1aqF169YICQmBu7s7kpKSsHr1aowcOfKx5y9rzEOGDMHt27fx6quvokqVKrh06RK++eYbvPLKK1IfjldeeQWWlpaYPn06MjMzoVarpTkqjK19+/ZSDc+wYcOQnZ2NRYsWwdPTUyepe5rXb9y4cVi+fDkiIyMRHR0Nd3d3LF26FBcuXMCvv/4q1Ti0b98eFStWRPPmzVGhQgWcOXMGc+fORceOHR/b7+1f//oXVq5ciXfffRc7d+5E8+bNodFo8Ndff2HlypXYsmWLThJhqJ9++gnz58+XOsXfu3cPW7ZswdatW9G5c2edmo0ePXqgSZMmGDRoEE6fPi3NTKvRaB7bjAUAzs7OmDVrFoYMGYKGDRtK8w8dP34cubm5WLp0KaytrTF9+nQMGjQIYWFh6NOnjzQ82c/PD2PGjHnm+31UbGws4uPj0bJlS4wYMQJFRUXSZ8vD/bbK+h5aunQpvv32W3Tr1g3+/v64d+8eFi1aBGdn5zJ92XnhPP+BRi+n4qGzxYuNjY2oWLGiaNeunZgzZ47OEOBijw5P3r59u+jSpYuoVKmSsLGxEZUqVRJ9+vTRG8a5bt06UatWLWFlZaUzVDksLKzUIYylDdVcvny5iImJEZ6ensLOzk507NhRZ4hdsX//+9+icuXKQq1Wi+bNm4ukpCS9cwrxYOjdJ598IqpWrSqsra1FxYoVRY8ePXSGHuORoX1CCHHkyBEREREhHB0dhb29vWjTpo3Yv39/ic/xo0P2iu+lpCGqj/r1119FUFCQUKvVolatWmLNmjViwIABesNQhRBi4cKFIiQkRNjZ2QknJydRp04d8dFHH4mrV68+8Trbtm0TzZs3F3Z2dsLZ2Vl07txZnD59WqdM8etf0nB0jUYjYmNjhZeXl7CzsxOtW7cWJ0+eLHGo+L1790RMTIwICAgQNjY2onz58qJZs2ZixowZ0jDc4uHJX331ld61Sno9UlNTRf/+/UXFihWFtbW1qFy5sujUqZNYvXq1VOaLL74QjRo1Eq6ursLOzk7UrFlTTJkyRWfob2nKEvPq1atF+/bthaenp7CxsRE+Pj5i2LBh4tq1azrnWrRokahWrZo0FLT4fVDae37VqlU6jy/tfVXS67N+/XpRt25dYWtrK/z8/MT06dOlYfkXLlyQyj3N65eamip69OghXF1dha2trWjUqJHYsGGDTpkFCxaIVq1aiXLlygm1Wi38/f3Fhx9+KDIzM5/4XBcUFIjp06eL4OBgoVarhZubmwgJCRGxsbE6jwcgoqKi9B5f2vQED0tMTBQ9e/YUPj4+Qq1WCwcHB9GgQQMxc+ZMUVhYqFf+9u3b4u233xblypUT9vb2Iiws7KmG4q5fv140a9ZM+v/VqFEjsXz5cp0yv/zyi6hfv75Qq9XC3d1d9O3bV/z99986ZQYMGCAcHBxKvEZpz0fxsUf/z+zevVuEhIQIGxsbUa1aNTF//ny9z/ji2J/0Hjpy5Ijo06eP9Hx6enqKTp06iaSkpDI/Ry8SlRBG6mFIREREZGTso0JERESKxUSFiIiIFIuJChERESkWExUiIiJSLCYqREREpFhMVIiIiEixXugJ37RaLa5evQonJydZpn8mIiKipyeEwL1791CpUqVSf0up2AudqFy9etUov/NCREREz19aWlqpv05d7IVOVIqnhk5LSzPK770QlVlODlD8K65XrwJP+CFBIiL6R1ZWFry9vR/7Ew/FXuhEpbi5x9nZmYkKPV8P/9y9szMTFSIiA5Sl2wY70xIREZFiMVEhIiIixXqhm36IZGNlBQwY8M86ERGZBD9hiQyhVgNLlsgdBdELTaPRoLCwUO4wyASsra1h+XBfvmfARIWIiJ4rIQSuX7+Ou3fvyh0KmZCrqysqVqz4zPOcMVEhMoQQQG7ug3V7e4ATDhKVWXGS4unpCXt7e07YaWaEEMjNzUV6ejoAwMvL65nOx0SFyBC5uYCj44P17GwOTyYqI41GIyUp5cqVkzscMhE7OzsAQHp6Ojw9PZ+pGYijfoiI6Lkp7pNib28vcyRkasWv8bP2Q2KiQkREzx2be8yfsV5jJipERESkWLImKn5+flCpVHpLVFSUnGERERGZ3JIlS+Dq6ip3GAbbtWsXVCqVyUdvyZqoJCYm4tq1a9KydetWAEDPnj3lDIuIiEhPWloaBg8ejEqVKsHGxga+vr4YNWoUMjIy5A7NrMmaqHh4eKBixYrSsmHDBvj7+yMsLEzOsIiIiHScP38eoaGhSE5OxvLly5GSkoL58+dj+/btaNq0KW7fvl3qYwsKCp5jpE+mtHieRDF9VAoKCvDTTz9h8ODB7GRFymdpCfTo8WAx0uyLRKRcUVFRsLGxQUJCAsLCwuDj44PIyEhs27YNV65cwSeffCKV9fPzw+TJk9G/f384Oztj6NChAB409fj4+MDe3h7dunUrsSZm3bp1aNCgAWxtbVGtWjXExsaiqKhIOq5SqfDdd9+hW7dusLe3R2BgINavX//Y2EuLZ+/evWjZsiXs7Ozg7e2N6Oho5OTkSI/78ccfERoaCicnJ1SsWBFvvfWWNDfKcyUU4pdffhGWlpbiypUrpZa5f/++yMzMlJa0tDQBQGRmZj7HSImIyFB5eXni9OnTIi8vT/dAdnbpy9OUzc0tW9mnkJGRIVQqlZg6dWqJx9955x3h5uYmtFqtEEIIX19f4ezsLGbMmCFSUlJESkqKOHjwoLCwsBDTp08XZ8+eFXPmzBGurq7CxcVFOs+ePXuEs7OzWLJkiUhNTRUJCQnCz89PTJo0SSoDQFSpUkUsW7ZMJCcni+joaOHo6CgyMjJKjb+keFJSUoSDg4OYNWuWOHfunNi3b5+oX7++GDhwoPS4//73v2LTpk0iNTVVHDhwQDRt2lRERkZKx3fu3CkAiDt37pR43VJfayFEZmZmmf9+q/7/xmUXEREBGxsb/P7776WWmTRpEmJjY/X2Z2ZmwtnZ2ZThEb1Q/MZtlDuEJ7o4raPcIZAM7t+/jwsXLqBq1aqwtbX958DjatI7dAA2PvSednD4Z2boR4WFAbt2/bPt4QHcuqVf7in+9B06dAhNmjTB2rVr0bVrV73js2bNwtixY3Hjxg14enrCz88P9evXx9q1a6Uyb731FjIzM7Hxoft48803ER8fL3VGDQ8PR9u2bRETEyOV+emnn/DRRx/h6tWrAB7UqHz66aeYPHkyACAnJweOjo7YvHkzXnvttRLjLymeIUOGwNLSEgsWLJD27d27F2FhYcjJydF9bf5fUlISGjZsiHv37sHR0RG7du1CmzZtcOfOnRI7BZf6WgPIysqCi4tLmf5+K6Lp59KlS9i2bRuGDBny2HIxMTHIzMyUlrS0tOcUIRERveye5nt9aGiozvaZM2fQuHFjnX1NmzbV2T5+/Dg+//xzODo6Sss777yDa9euIfehxKxu3brSuoODA5ydnZ/YJPNoPMePH8eSJUt0rhUREQGtVosLFy4AAA4fPozOnTvDx8cHTk5OUv/Ry5cvl/FZMA5FTKG/ePFieHp6omPHx3/DUqvVUKvVzykqosfIyeEU+kTGlJ1d+rFH+4E97o+yxSPfvy9eNDikYgEBAVCpVDhz5gy6deumd/zMmTNwc3ODh4eHtM/BgM+E7OxsxMbGonv37nrHHq6RsLa21jmmUqmg1Wofe+5H48nOzsawYcMQHR2tV9bHxwc5OTmIiIhAREQEfv75Z3h4eODy5cuIiIh47p1xZU9UtFotFi9ejAEDBsDKSvZwiIhIDk/zh91UZUtRrlw5tGvXDt9++y3GjBkj/Y4N8OAHFn/++Wf079//sQNBgoKCcOjQIZ19Bw8e1Nlu0KABzp49i4CAgGeO+UkaNGiA06dPl3qtEydOICMjA9OmTYO3tzeAB00/cpC96Wfbtm24fPkyBg8eLHcoREREJZo7dy7y8/MRERGBPXv2IC0tDfHx8WjXrh0qV66MKVOmPPbx0dHRiI+Px4wZM5CcnIy5c+ciPj5ep8yECRPwww8/IDY2FqdOncKZM2ewYsUKfPrpp0a/n48//hj79+/HyJEjcezYMSQnJ2PdunUYOXIkgAe1KjY2Nvjmm29w/vx5rF+/XuoX87zJnqi0b98eQghUr15d7lCIiIhKFBgYiKSkJFSrVg29evWCv78/hg4dijZt2uDAgQNwd3d/7OObNGmCRYsWYc6cOahXrx4SEhL0EpCIiAhs2LABCQkJaNiwIZo0aYJZs2bB19fX6PdTt25d7N69G+fOnUPLli1Rv359TJgwAZUqVQLwYJ6zJUuWYNWqVahVqxamTZuGGTNmGD2OslDMqB9DPE2vYSKjUngfFY76IaV63EgQMi9mNeqHiIiIqCRMVIiIiEixOMyGyBCWlg8moSpeJyIik2CiQmQIW1vdmTKJiMgk2PRDRETP3Qs8joPKyFivMRMVIiJ6bopnVc0t7bd6yGwUv8aPzqT7tNj0Q2SInBzA0/PBenq64oYnEymVpaUlXF1dpd+msbe3f+yMrvTiEUIgNzcX6enpcHV1heUz9uNjokJkKH4jJDJIxYoVAeCJP6RHLzZXV1fptX4WTFSIiOi5UqlU8PLygqenJwoLC+UOh0zA2tr6mWtSijFRISIiWVhaWhrtjxmZL3amJSIiIsViokJERESKxUSFiIiIFIt9VIgMYWEBhIX9s05ERCbBRIXIEHZ2wK5dckdBRGT2+FWQiIiIFIuJChERESkWExUiQ+TkAB4eD5acHLmjISIyW+yjQmSoW7fkjoCIyOyxRoWIiIgUi4kKERERKRYTFSIiIlIsJipERESkWExUiIiISLE46ofIEBYWQGjoP+tERGQSTFSIDGFnByQmyh0FEZHZ41dBIiIiUiwmKkRERKRYTFSIDJGbC/j5PVhyc+WOhojIbLGPCpEhhAAuXfpnnYiITII1KkRERKRYTFSIiIhIsZioEBERkWIxUSEiIiLFYqJCREREisVRP0SGUKmAWrX+WSciIpNgokJkCHt74NQpuaMgIjJ7bPohIiIixWKiQkRERIole6Jy5coV9OvXD+XKlYOdnR3q1KmDpKQkucMierzcXCA4+MHCKfSJiExG1j4qd+7cQfPmzdGmTRts3rwZHh4eSE5Ohpubm5xhET2ZEMDp0/+sExGRSciaqEyfPh3e3t5YvHixtK9q1aoyRkRERERKImvTz/r16xEaGoqePXvC09MT9evXx6JFi0otn5+fj6ysLJ2FiIiIzJesicr58+cxb948BAYGYsuWLRg+fDiio6OxdOnSEsvHxcXBxcVFWry9vZ9zxERERPQ8qYSQr4HdxsYGoaGh2L9/v7QvOjoaiYmJOHDggF75/Px85OfnS9tZWVnw9vZGZmYmnJ2dn0vMRACAnBzA0fHBenY24OAgbzyP8Bu3Ue4QnujitI5yh0BEMsnKyoKLi0uZ/n7LWqPi5eWFWsWze/6/oKAgXL58ucTyarUazs7OOgsRERGZL1k70zZv3hxnz57V2Xfu3Dn4+vrKFBFRGalUQPH7lFPoExGZjKyJypgxY9CsWTNMnToVvXr1wh9//IGFCxdi4cKFcoZF9GT29sDFi3JHQURk9mRt+mnYsCHWrl2L5cuXo3bt2pg8eTJmz56Nvn37yhkWERERKYTsP0rYqVMndOrUSe4wiIiISIFkn0Kf6IWUlwc0bPhgycuTOxoiIrMle40K0QtJqwWKf5NKq5U3FiIiM8YaFSIiIlIsJipERESkWExUiIiISLGYqBAREZFiMVEhIiIixeKoHyJDlS8vdwRERGaPiQqRIRwcgJs35Y6CiMjssemHiIiIFIuJChERESkWExUiQ+TlAa1bP1g4hT4RkcmwjwqRIbRaYPfuf9aJiMgkWKNCREREisVEhYiIiBSLiQoREREpFhMVIiIiUiwmKkRERKRYHPVDZCh7e7kjICIye0xUiAzh4ADk5MgdBRGR2WPTDxERESkWExUiIiJSLCYqRIa4fx/o2PHBcv++3NEQEZkt9lEhMoRGA2za9M86ERGZBGtUiIiISLGYqBAREZFiMVEhIiIixWKiQkRERIrFRIWIiIgUi4kKERERKRaHJxMZwsEBEELuKIiIzB5rVIiIiEixmKgQERGRYjFRITLE/ftAz54PFk6hT0RkMkxUiAyh0QCrVz9YOIU+EZHJMFEhIiIixWKiQkRERIrFRIWIiIgUi4kKERERKZasicqkSZOgUql0lpo1a8oZEhERESmI7DPTBgcHY9u2bdK2lZXsIREREZFCyJ4VWFlZoWLFinKHQfR07O2B7Ox/1omIyCRk76OSnJyMSpUqoVq1aujbty8uX75catn8/HxkZWXpLESyUKke/N6Pg8ODdSIiMglZa1QaN26MJUuWoEaNGrh27RpiY2PRsmVLnDx5Ek5OTnrl4+LiEBsbK0OkZC78xm2UO4Qnujito9whkALxvUsvK1lrVCIjI9GzZ0/UrVsXERER2LRpE+7evYuVK1eWWD4mJgaZmZnSkpaW9pwjJnrApqgQMzbOwoyNs2BTVCh3OEREZkv2PioPc3V1RfXq1ZGSklLicbVaDbVa/ZyjItJnqdWgx8ntAIDP2g0HYC1vQEREZkr2PioPy87ORmpqKry8vOQOhYiIiBRA1kTlgw8+wO7du3Hx4kXs378f3bp1g6WlJfr06SNnWERERKQQsjb9/P333+jTpw8yMjLg4eGBFi1a4ODBg/Dw8JAzLCIiIlIIWROVFStWyHl5IiIiUjhF9VEhIiIiehgTFSIiIlIsRQ1PJnpR5Fmr0eC9n6V1IiIyDSYqRIZQqXDb3kXuKIiIzB6bfoiIiEixWKNCZACbokJ8uuM7AMAXrw5BgRVnpiUiMgXWqBAZwFKrQf+jG9H/6EZYajVyh0NEZLaYqBAREZFiMVEhIiIixWKiQkRERIrFRIWIiIgUi4kKERERKRYTFSIiIlIszqNCZID71jZo8e5/pXUiIjINJipEBhAqC/ztUkHuMIiIzB6bfoiIiEixWKNCZABrTSE+2PMjAGBGq3+h0JJT6BMRmQJrVIgMYKXRYNgfazDsjzWw0nAKfSIiU2GiQkRERIrFRIWIiIgUi4kKERERKRYTFSIiIlIsJipERESkWExUiIiISLE4jwqRAe5b26Dd4P9I60REZBpMVIgMIFQWSPbwlTsMIiKzx6YfIiIiUiyDa1RycnKwe/duXL58GQUFBTrHoqOjnzkwIiWz1hQi6sBKAMB/mvbiFPpERCZiUKJy9OhRdOjQAbm5ucjJyYG7uztu3boFe3t7eHp6MlEhs2el0WD0vuUAgAWN3mCiQkRkIgY1/YwZMwadO3fGnTt3YGdnh4MHD+LSpUsICQnBjBkzjB0jERERvaQMSlSOHTuG999/HxYWFrC0tER+fj68vb3x5ZdfYvz48caOkYiIiF5SBiUq1tbWsLB48FBPT09cvnwZAODi4oK0tDTjRUdEREQvNYP6qNSvXx+JiYkIDAxEWFgYJkyYgFu3buHHH39E7dq1jR0jERERvaQMqlGZOnUqvLy8AABTpkyBm5sbhg8fjps3b2LhwoVGDZCIiIheXgbVqISGhkrrnp6eiI+PN1pARERERMU4My2RAfKtrPF6/5nSOhERmUaZE5UGDRpg+/btcHNzQ/369aFSqUote+TIEaMER6RUWgtL/OlVXe4wiIjMXpkTlS5dukCtVgMAunbtaqp4iIiIiCRlTlQmTpxY4jrRy8haU4hBSesBAItDX+fMtEREJmLQqJ/ExEQcOnRIb/+hQ4eQlJRkUCDTpk2DSqXC6NGjDXo80fNkpdFg/K7FGL9rMaw0GrnDISIyWwYlKlFRUSVO7HblyhVERUU99fkSExOxYMEC1K1b15BwiIiIyEwZlKicPn0aDRo00Ntfv359nD59+qnOlZ2djb59+2LRokVwc3MzJBwiIiIyUwYlKmq1Gjdu3NDbf+3aNVhZPd2I56ioKHTs2BHh4eGGhEJERERmzKB5VNq3b4+YmBisW7cOLi4uAIC7d+9i/PjxaNeuXZnPs2LFChw5cgSJiYllKp+fn4/8/HxpOysr6+kCJyIioheKQYnKjBkz0KpVK/j6+qJ+/foAHvyicoUKFfDjjz+W6RxpaWkYNWoUtm7dCltb2zI9Ji4uDrGxsYaEbBC/cRuf27UMdXFaR7lDIJIV/58SmTeDmn4qV66MP//8E19++SVq1aqFkJAQzJkzBydOnIC3t3eZznH48GGkp6ejQYMGsLKygpWVFXbv3o2vv/4aVlZW0JQwkiImJgaZmZnSwl9qJiIiMm8GT6Hv4OCAoUOHGnzhtm3b4sSJEzr7Bg0ahJo1a+Ljjz+GpaWl3mPUarU06RyRnPKtrPFmn6nSOhERmYbBiUpycjJ27tyJ9PR0aLVanWMTJkx44uOdnJxQu3ZtnX0ODg4oV66c3n4ipdFaWOKgD4fTExGZmkGJyqJFizB8+HCUL18eFStW1PndH5VKVaZEhYiIiOhJDEpUvvjiC0yZMgUff/yxUYPZtWuXUc9HZCpWmiL0OR4PAFhe7zUUWfKHyImITMGgT9c7d+6gZ8+exo6F6IVhrSnC5K3zAQCra4czUSEiMhGDRv307NkTCQkJxo6FiIiISIdBXwMDAgLw2Wef4eDBg6hTpw6srXVHPURHRxslOCIiInq5GZSoLFy4EI6Ojti9ezd2796tc0ylUjFRISIiIqMwKFG5cOGCseMgIiIi0mNQH5ViBQUFOHv2LIqKiowVDxEREZHEoEQlNzcXb7/9Nuzt7REcHIzLly8DAN577z1MmzbNqAESERHRy8ugRCUmJgbHjx/Hrl27dH5QMDw8HL/88ovRgiNSqgIrawzqMRGDekxEAafQJyIyGYP6qPz222/45Zdf0KRJE51ZaYODg5Gammq04IiUSmNhiZ3+DeUOg4jI7BlUo3Lz5k14enrq7c/JydFJXIiIiIiehUGJSmhoKDZu3ChtFycn3333HZo2bWqcyIgUzEpThB4ntqHHiW2w0rAzORGRqRjU9DN16lRERkbi9OnTKCoqwpw5c3D69Gns379fb14VInNkrSnCjE2zAQAba7TgFPpERCZiUI1KixYtcOzYMRQVFaFOnTpISEiAp6cnDhw4gJCQEGPHSERERC8pg78G+vv7Y9GiRcaMhYiIiEiHQYlK8bwppfHx8TEoGCIiIqKHGZSo+Pn5PXZ0j0ajMTggIiIiomIGJSpHjx7V2S4sLMTRo0cxc+ZMTJkyxSiBERERERmUqNSrV09vX2hoKCpVqoSvvvoK3bt3f+bAiIiIiIw6prJGjRpITEw05imJFKnAyhojuoyT1omIyDQMSlSysrJ0toUQuHbtGiZNmoTAwECjBEakZBoLS2yq2ULuMIiIzJ5BiYqrq6teZ1ohBLy9vbFixQqjBEZERERkUKKyY8cOnUTFwsICHh4eCAgIgJUVZ+gk82ep1SDi3AEAwJbqTaGxsJQ5IiIi82RQVtG6dWsjh0H0YrEpKsS366YBAILGrEaeDRMVIiJTMGgK/bi4OHz//fd6+7///ntMnz79mYMiIiIiAgxMVBYsWICaNWvq7Q8ODsb8+fOfOSgiIiIiwMBE5fr16/Dy8tLb7+HhgWvXrj1zUERERESAgYmKt7c39u3bp7d/3759qFSp0jMHRURERAQY2Jn2nXfewejRo1FYWIhXX30VALB9+3Z89NFHeP/9940aIBEREb28DEpUPvzwQ2RkZGDEiBEoKCgAANja2uLjjz9GTEyMUQMkIiKil5dBiYpKpcL06dPx2Wef4cyZM7Czs0NgYCDUarWx4yNSpEJLK3zQYbS0TkREpvFMn7DXr1/H7du30apVK6jVaggh9GasJTJHRZZWWF0nXO4wiIjMnkGdaTMyMtC2bVtUr14dHTp0kEb6vP322+yjQkREREZjUKIyZswYWFtb4/Lly7C3t5f29+7dG/Hx8UYLjkipLLUatElNRJvURFhqNXKHQ0Rktgxq+klISMCWLVtQpUoVnf2BgYG4dOmSUQIjUjKbokIsXh0LgFPoExGZkkE1Kjk5OTo1KcVu377NDrVERERkNAYlKi1btsQPP/wgbatUKmi1Wnz55Zdo06aN0YIjIiKil5tBTT9ffvkl2rZti6SkJBQUFOCjjz7CqVOncPv27RJnrCUiIiIyhEE1KrVr18a5c+fQokULdOnSBTk5OejevTuOHj0Kf39/Y8dIREREL6mnrlEpLCzEa6+9hvnz5+OTTz4xRUxEREREAAyoUbG2tsaff/5piliIiIiIdBjU9NOvXz/897//feaLz5s3D3Xr1oWzszOcnZ3RtGlTbN68+ZnPS2RqhZZW+Kzdu/is3bucQp+IyIQM+oQtKirC999/j23btiEkJAQODg46x2fOnFmm81SpUgXTpk1DYGAghBBYunQpunTpgqNHjyI4ONiQ0IieiyJLK/zYoJPcYRARmb2nSlTOnz8PPz8/nDx5Eg0aNAAAnDt3TqfM0/zWT+fOnXW2p0yZgnnz5uHgwYNMVIiIiOjpEpXAwEBcu3YNO3fuBPBgyvyvv/4aFSpUeOZANBoNVq1ahZycHDRt2rTEMvn5+cjPz5e2s7Kynvm6RIaw0GrQ6O9TAIA/qgRDa8GZaYmITOGpEhUhhM725s2bkZOT80wBnDhxAk2bNsX9+/fh6OiItWvXolatWiWWjYuLQ2xs7DNdj8gY1EWFWLF8PABOoU9EZEoGdaYt9mjiYogaNWrg2LFjOHToEIYPH44BAwbg9OnTJZaNiYlBZmamtKSlpT3z9YmIiEi5nqpGRaVS6fVBeZo+KSWxsbFBQEAAACAkJASJiYmYM2cOFixYoFdWrVbzt4SIiIheIk/d9DNw4EApWbh//z7effddvVE/a9asMTggrVar0w+FiIiIXl5PlagMGDBAZ7tfv37PdPGYmBhERkbCx8cH9+7dw7Jly7Br1y5s2bLlmc5LRERE5uGpEpXFixcb9eLp6eno378/rl27BhcXF9StWxdbtmxBu3btjHodIiIiejHJOqWmMWa3JSIiIvPFub+JDFBkaYmprQdJ60REZBpMVIgMUGhpjYWN35A7DCIis/dM86gQERERmRJrVIgMYKHVoPaNVADAyQr+nEKfiMhEmKgQGUBdVIj1P4wFwCn0iYhMiU0/REREpFhMVIiIiEixmKgQERGRYjFRISIiIsViokJERESKxUSFiIiIFIvDk4kMUGRpidnN+0jrRERkGkxUiAxQaGmN2S36yh0GEZHZY9MPERERKRZrVIgMoBJaBNxKAwCklPeGUDHnJyIyBSYqRAawLSzA1u+jABRPoW8rc0REROaJXwOJiIhIsZioEBERkWIxUSEiIiLFYqJCREREisVEhYiIiBSLiQoREREpFocnExmgyNISCxp1l9aJiMg0mKgQGaDQ0hpxbQbLHQYRkdlj0w8REREpFmtUiAygElpUzroJALji7MEp9ImITISJCpEBbAsLsHf+2wA4hT4RkSnxayAREREpFhMVIiIiUiwmKkRERKRYTFSIiIhIsZioEBERkWIxUSEiIiLF4vBkIgNoLCzxQ/2O0joREZkGExUiAxRYWWNC++Fyh0FEZPbY9ENERESKxRoVIkMIAfe8LADAbTtnQKWSOSAiIvPERIXIAHaF+TjyTV8AnEKfiMiU2PRDREREiiVrohIXF4eGDRvCyckJnp6e6Nq1K86ePStnSERERKQgsiYqu3fvRlRUFA4ePIitW7eisLAQ7du3R05OjpxhERERkULI2kclPj5eZ3vJkiXw9PTE4cOH0apVK5miIiIiIqVQVGfazMxMAIC7u3uJx/Pz85Gfny9tZ2VlPZe4iIiISB6KSVS0Wi1Gjx6N5s2bo3bt2iWWiYuLQ2xs7HOOjIiInie/cRvlDuGJLk7raJTzvEz3aijFJCpRUVE4efIk9u7dW2qZmJgYjB07VtrOysqCt7f38wiPSIfGwhKra7eV1omIyDQUkaiMHDkSGzZswJ49e1ClSpVSy6nVaqjV6ucYGVHJCqys8UHHMXKHQURk9mRNVIQQeO+997B27Vrs2rULVatWlTMcIiIiUhhZE5WoqCgsW7YM69atg5OTE65fvw4AcHFxgZ2dnZyhET2eELArfNCxO89azSn0iYhMRNZ5VObNm4fMzEy0bt0aXl5e0vLLL7/IGRbRE9kV5uPMrB44M6uHlLAQEZHxyd70Q0RERFQa/tYPERERKRYTFSIiIlIsJipERESkWExUiIiISLGYqBAREZFiKWJmWqIXjdbCAhtrNJfWiYjINJioEBkg38oGUV1j5A6DiMjs8asgERERKRYTFSIiIlIsJipEBrAruI+L0zvh4vROsCu4L3c4RERmi4kKERERKRYTFSIiIlIsJipERESkWExUiIiISLGYqBAREZFiMVEhIiIixeLMtEQG0FpYYEe1UGmdiIhMg4kKkQHyrWwwuOckucMgIjJ7/CpIREREisVEhYiIiBSLiQqRAewK7uP0zDdweuYbnEKfiMiE2EeFyED2hflyh0BEZPZYo0JERESKxUSFiIiIFIuJChERESkWExUiIiJSLCYqREREpFgc9UNkAK1KhYPetaV1IiIyDSYqRAbIt1bjzbemyR0GEZHZY9MPERERKRYTFSIiIlIsJipEBrAruI/DX7+Fw1+/xSn0iYhMiH1UiAxULi9L7hCIiMwea1SIiIhIsZioEBERkWIxUSEiIiLFYqJCREREisVEhYiIiBSLo36IDKBVqXC8YqC0TkREpiFrjcqePXvQuXNnVKpUCSqVCr/99puc4RCVWb61Gl0GzEKXAbOQb62WOxwiIrMla6KSk5ODevXq4T//+Y+cYRAREZFCydr0ExkZicjISDlDICIiIgV7ofqo5OfnIz8/X9rOyuLMoCQP28L72PbdCABA+JBvcd/aVuaIiIjM0wuVqMTFxSE2NlbuMMyO37iNcofwRBendZQ7BB0qAVTJSpfWiYjINF6o4ckxMTHIzMyUlrS0NLlDIiIiIhN6oWpU1Go11GqOsCAiInpZvFA1KkRERPRykbVGJTs7GykpKdL2hQsXcOzYMbi7u8PHx0fGyIiIiEgJZE1UkpKS0KZNG2l77NixAIABAwZgyZIlMkVFRERESiFrotK6dWsIwSET9OIRKuBcOR9pnYiITOOF6kxLpBT3rW3Rfsi3codBRGT22JmWiIiIFIuJChERESkWExUiA9gW3kfCdyOQ8N0I2BbelzscIiKzxT4qRAZQCaB6xmVpnYiITIM1KkRERKRYTFSIiIhIsZioEBERkWIxUSEiIiLFYqJCREREisVRP0QGECrgb2dPaZ2IiEyDiQqRAe5b26LF8O/lDoOIyOyx6YeIiIgUi4kKERERKRYTFSIDqAvzsW7pGKxbOgbqwny5wyEiMlvso0JkAAshUO96srRORESmwRoVIiIiUiwmKkRERKRYTFSIiIhIsZioEBERkWIxUSEiIiLF4qgfIgNl2DnLHQIRkdljokJkgDwbW4REL5M7DCIis8emHyIiIlIsJipERESkWExUiAygLszHimXjsGLZOE6hT0RkQuyjQmQACyHQJO2ktE5ERKbBGhUiIiJSLCYqREREpFhMVIiIiEixmKgQERGRYjFRISIiIsXiqB8iA+Vaq+UOgYjI7DFRITJAno0tao39Ve4wiIjMHpt+iIiISLGYqBAREZFisemHyADqogLMWzsVADC823jkW9nIHBERkXliokJkAAutFq+eT5LWiYjINNj0Q0RERIqliETlP//5D/z8/GBra4vGjRvjjz/+kDskIiIiUgDZE5VffvkFY8eOxcSJE3HkyBHUq1cPERERSE9Plzs0IiIikpnsicrMmTPxzjvvYNCgQahVqxbmz58Pe3t7fP/993KHRkRERDKTNVEpKCjA4cOHER4eLu2zsLBAeHg4Dhw4IGNkREREpASyjvq5desWNBoNKlSooLO/QoUK+Ouvv/TK5+fnIz8/X9rOzMwEAGRlZZkkPm1+rknOa0zGuPeX5T4B492rpuA+iiPS5OdCK4w38oev6dN5We71ZblPgPeqNKb4G1t8TiHEkwsLGV25ckUAEPv379fZ/+GHH4pGjRrplZ84caIAwIULFy5cuHAxgyUtLe2JuYKsNSrly5eHpaUlbty4obP/xo0bqFixol75mJgYjB07VtrWarW4ffs2ypUrB5VKZfJ4n0VWVha8vb2RlpYGZ2dnucMxqZflXl+W+wRennt9We4T4L2aoxfpPoUQuHfvHipVqvTEsrImKjY2NggJCcH27dvRtWtXAA+Sj+3bt2PkyJF65dVqNdRq3V+sdXV1fQ6RGo+zs7Pi30DG8rLc68tyn8DLc68vy30CvFdz9KLcp4uLS5nKyT4z7dixYzFgwACEhoaiUaNGmD17NnJycjBo0CC5QyMiIiKZyZ6o9O7dGzdv3sSECRNw/fp1vPLKK4iPj9frYEtEREQvH9kTFQAYOXJkiU095kStVmPixIl6TVfm6GW515flPoGX515flvsEeK/myFzvUyVEWcYGERERET1/ss9MS0RERFQaJipERESkWExUiIiISLGYqBAREZFiMVExsT179qBz586oVKkSVCoVfvvtN7lDMom4uDg0bNgQTk5O8PT0RNeuXXH27Fm5wzKJefPmoW7dutKkSk2bNsXmzZvlDsvkpk2bBpVKhdGjR8sditFNmjQJKpVKZ6lZs6bcYZnMlStX0K9fP5QrVw52dnaoU6cOkpKS5A7LqPz8/PReU5VKhaioKLlDMzqNRoPPPvsMVatWhZ2dHfz9/TF58uSy/Y7OC0ARw5PNWU5ODurVq4fBgweje/fucodjMrt370ZUVBQaNmyIoqIijB8/Hu3bt8fp06fh4OAgd3hGVaVKFUybNg2BgYEQQmDp0qXo0qULjh49iuDgYLnDM4nExEQsWLAAdevWlTsUkwkODsa2bdukbSsr8/x4vHPnDpo3b442bdpg8+bN8PDwQHJyMtzc3OQOzagSExOh0Wik7ZMnT6Jdu3bo2bOnjFGZxvTp0zFv3jwsXboUwcHBSEpKwqBBg+Di4oLo6Gi5w3tm5vk/UUEiIyMRGRkpdxgmFx8fr7O9ZMkSeHp64vDhw2jVqpVMUZlG586ddbanTJmCefPm4eDBg2aZqGRnZ6Nv375YtGgRvvjiC7nDMRkrK6sSf2PM3EyfPh3e3t5YvHixtK9q1aoyRmQaHh4eOtvTpk2Dv78/wsLCZIrIdPbv348uXbqgY8eOAB7UJi1fvhx//PGHzJEZB5t+yCQyMzMBAO7u7jJHYloajQYrVqxATk4OmjZtKnc4JhEVFYWOHTsiPDxc7lBMKjk5GZUqVUK1atXQt29fXL58We6QTGL9+vUIDQ1Fz5494enpifr162PRokVyh2VSBQUF+OmnnzB48GDF/4CtIZo1a4bt27fj3LlzAIDjx49j7969ZvMlmTUqZHRarRajR49G8+bNUbt2bbnDMYkTJ06gadOmuH//PhwdHbF27VrUqlVL7rCMbsWKFThy5AgSExPlDsWkGjdujCVLlqBGjRq4du0aYmNj0bJlS5w8eRJOTk5yh2dU58+fx7x58zB27FiMHz8eiYmJiI6Oho2NDQYMGCB3eCbx22+/4e7duxg4cKDcoZjEuHHjkJWVhZo1a8LS0hIajQZTpkxB37595Q7NKJiokNFFRUXh5MmT2Lt3r9yhmEyNGjVw7NgxZGZmYvXq1RgwYAB2795tVslKWloaRo0aha1bt8LW1lbucEzq4W+edevWRePGjeHr64uVK1fi7bffljEy49NqtQgNDcXUqVMBAPXr18fJkycxf/58s01U/vvf/yIyMhKVKlWSOxSTWLlyJX7++WcsW7YMwcHBOHbsGEaPHo1KlSqZxWvKRIWMauTIkdiwYQP27NmDKlWqyB2OydjY2CAgIAAAEBISgsTERMyZMwcLFiyQOTLjOXz4MNLT09GgQQNpn0ajwZ49ezB37lzk5+fD0tJSxghNx9XVFdWrV0dKSorcoRidl5eXXkIdFBSEX3/9VaaITOvSpUvYtm0b1qxZI3coJvPhhx9i3LhxePPNNwEAderUwaVLlxAXF8dEhaiYEALvvfce1q5di127dpll57zH0Wq1yM/PlzsMo2rbti1OnDihs2/QoEGoWbMmPv74Y7NNUoAHHYhTU1Pxr3/9S+5QjK558+Z6UwecO3cOvr6+MkVkWosXL4anp6fU0dQc5ebmwsJCt8uppaUltFqtTBEZFxMVE8vOztb5VnbhwgUcO3YM7u7u8PHxkTEy44qKisKyZcuwbt06ODk54fr16wAAFxcX2NnZyRydccXExCAyMhI+Pj64d+8eli1bhl27dmHLli1yh2ZUTk5Oen2MHBwcUK5cObPre/TBBx+gc+fO8PX1xdWrVzFx4kRYWlqiT58+codmdGPGjEGzZs0wdepU9OrVC3/88QcWLlyIhQsXyh2a0Wm1WixevBgDBgww2+HmwIORiFOmTIGPjw+Cg4Nx9OhRzJw5E4MHD5Y7NOMQZFI7d+4UAPSWAQMGyB2aUZV0jwDE4sWL5Q7N6AYPHix8fX2FjY2N8PDwEG3bthUJCQlyh/VchIWFiVGjRskdhtH17t1beHl5CRsbG1G5cmXRu3dvkZKSIndYJvP777+L2rVrC7VaLWrWrCkWLlwod0gmsWXLFgFAnD17Vu5QTCorK0uMGjVK+Pj4CFtbW1GtWjXxySefiPz8fLlDMwqVEGYydR0RERGZHc6jQkRERIrFRIWIiIgUi4kKERERKRYTFSIiIlIsJipERESkWExUiIiISLGYqBAREZFiMVEhIj1LliyBq6ur3GEYpHXr1hg9evRLd20ic8VEhcgMpKWlYfDgwahUqRJsbGzg6+uLUaNGISMjQ+7QTGbXrl1QqVS4e/euzv41a9Zg8uTJ8gRFREbHRIXoBXf+/HmEhoYiOTkZy5cvR0pKCubPn4/t27ejadOmuH37dqmPLSgoeI6RPh/u7u5wcnKSOwwiMhImKkQvuKioKNjY2CAhIQFhYWHw8fFBZGQktm3bhitXruCTTz6Ryvr5+WHy5Mno378/nJ2dMXToUAAPmnp8fHxgb2+Pbt26lVgTs27dOjRo0AC2traoVq0aYmNjUVRUJB1XqVT47rvv0K1bN9jb2yMwMBDr169/bOz5+fn44IMPULlyZTg4OKBx48bYtWuXdPzSpUvo3Lkz3Nzc4ODggODgYGzatAkXL15EmzZtAABubm5QqVQYOHAgAP3mFz8/P3zxxRfo378/HB0d4evri/Xr1+PmzZvo0qULHB0dUbduXSQlJUmPycjIQJ8+fVC5cmXY29ujTp06WL58uU7sOTk50jm9vLzw73//W+/+7ty5g/79+8PNzQ329vaIjIxEcnLyY58TInqE3D82RESGy8jIECqVSkydOrXE4++8845wc3MTWq1WCCGEr6+vcHZ2FjNmzBApKSkiJSVFHDx4UFhYWIjp06eLs2fPijlz5ghXV1fh4uIinWfPnj3C2dlZLFmyRKSmpoqEhATh5+cnJk2aJJUBIKpUqSKWLVsmkpOTRXR0tHB0dBQZGRmlxj9kyBDRrFkzsWfPHpGSkiK++uoroVarxblz54QQQnTs2FG0a9dO/PnnnyI1NVX8/vvvYvfu3aKoqEj8+uuv0g/OXbt2Tdy9e1cIof/Dib6+vsLd3V3Mnz9fnDt3TgwfPlw4OzuL1157TaxcuVKcPXtWdO3aVQQFBUnP099//y2++uorcfToUZGamiq+/vprYWlpKQ4dOiSdd/jw4cLHx0ds27ZN/Pnnn6JTp07CyclJ59qvv/66CAoKEnv27BHHjh0TERERIiAgQBQUFJTtBSYiwUSF6AV28OBBAUCsXbu2xOMzZ84UAMSNGzeEEA/+aHft2lWnTJ8+fUSHDh109vXu3VsnUWnbtq1eMvTjjz8KLy8vaRuA+PTTT6Xt7OxsAUBs3ry5xNguXbokLC0txZUrV3T2t23bVsTExAghhKhTp45OMvSw4l8mv3Pnjs7+khKVfv36SdvXrl0TAMRnn30m7Ttw4IAAIK5du1bitYR4kDS9//77Qggh7t27J2xsbMTKlSul4xkZGcLOzk669rlz5wQAsW/fPqnMrVu3hJ2dnc7jiOjxrOSqySEi4xFP8SPooaGhOttnzpxBt27ddPY1bdoU8fHx0vbx48exb98+TJkyRdqn0Whw//595Obmwt7eHgBQt25d6biDgwOcnZ2Rnp5eYhwnTpyARqNB9erVdfbn5+ejXLlyAIDo6GgMHz4cCQkJCA8PxxtvvKFzjbJ6+DEVKlQAANSpU0dvX3p6OipWrAiNRoOpU6di5cqVuHLlCgoKCpCfny/dZ2pqKgoKCtC4cWPpHO7u7qhRo4a0febMGVhZWemUKVeuHGrUqIEzZ8489T0QvayYqBC9wAICAqBSqUpMNoAHfyzd3Nzg4eEh7XNwcHjq62RnZyM2Nhbdu3fXO2ZrayutW1tb6xxTqVTQarWlntPS0hKHDx+GpaWlzjFHR0cAwJAhQxAREYGNGzciISEBcXFx+Pe//4333nvvqeJ/OC6VSlXqvuJYv/rqK8yZMwezZ89GnTp14ODggNGjR5tl52MipWNnWqIXWLly5dCuXTt8++23yMvL0zl2/fp1/Pzzz+jdu7f0h7gkQUFBOHTokM6+gwcP6mw3aNAAZ8+eRUBAgN5iYWHYx0j9+vWh0WiQnp6ud86KFStK5by9vfHuu+9izZo1eP/997Fo0SIAgI2NDYAHNTvGtm/fPnTp0gX9+vVDvXr1UK1aNZw7d0467u/vD2tra53n7c6dOzplgoKCUFRUpFMmIyMDZ8+eRa1atYweM5G5YqJC9IKbO3cu8vPzERERgT179iAtLQ3x8fFo164dKleurNNcU5Lo6GjEx8djxowZSE5Oxty5c3WafQBgwoQJ+OGHHxAbG4tTp07hzJkzWLFiBT799FOD465evTr69u2L/v37Y82aNbhw4QL++OMPxMXFYePGjQCA0aNHY8uWLbhw4QKOHDmCnTt3IigoCADg6+sLlUqFDRs24ObNm8jOzjY4lkcFBgZi69at2L9/P86cOYNhw4bhxo0b0nFHR0e8/fbb+PDDD7Fjxw6cPHkSAwcO1EnaAgMD0aVLF7zzzjvYu3cvjh8/jn79+qFy5cro0qWL0WIlMndMVIhecIGBgUhKSkK1atXQq1cv+Pv7Y+jQoWjTpg0OHDgAd3f3xz6+SZMmWLRoEebMmYN69eohISFBLwGJiIjAhg0bkJCQgIYNG6JJkyaYNWsWfH19nyn2xYsXo3///nj//fdRo0YNdO3aFYmJifDx8QHwoLYkKioKQUFBeO2111C9enV8++23AIDKlSsjNjYW48aNQ4UKFTBy5MhniuVhn376KRo0aICIiAi0bt0aFStWRNeuXXXKfPXVV2jZsiU6d+6M8PBwtGjRAiEhIXr3FxISgk6dOqFp06YQQmDTpk16TWREVDqVeJpeeERERETPEWtUiIiISLGYqBAREZFiMVEhIiIixWKiQkRERIrFRIWIiIgUi4kKERERKRYTFSIiIlIsJipERESkWExUiIiISLGYqBAREZFiMVEhIiIixWKiQkRERIr1f6BsHll9ObzzAAAAAElFTkSuQmCC",
|
|
"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"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Segundo intento del detector de orden\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Generador de muestras\n",
|
|
"\n",
|
|
"Esta sección se describe los cambios y mejoras implementadas en el generador de señales cardiorespiratorias sintéticas para simulaciones con radar UWB. La función ha sido adaptada para generar señales más realistas, controlables y útiles para métodos no lineales de análisis espectral.\n",
|
|
"\n",
|
|
"- Se modificaron los intervalos de frecuencia fundamentales para que reflejen rangos fisiológicos realistas para mediciones con radar UWB: la frecuencia cardíaca se estableció entre 0.6 y 1.1 Hz, lo que corresponde aproximadamente a 36 a 66 latidos por minuto, y la frecuencia respiratoria entre 0.15 y 0.35 Hz, equivalente a 9 a 21 respiraciones por minuto. Los centros de estos intervalos se utilizan como semillas para la inicialización de métodos no lineales, siendo 0.85 Hz para la cardíaca y 0.25 Hz para la respiratoria.\n",
|
|
"- Se definió que el número de componentes de cada señal corresponde al fundamental más los armónicos especificados por el usuario. Esto permite simular señales más complejas y con mayor riqueza espectral, manteniendo control sobre la cantidad de información armónica incluida.\n",
|
|
"- Los armónicos no son múltiplos exactos de la frecuencia fundamental. Se introduce un factor de dispersión aleatoria de ±5% para cada armónico, de modo que las frecuencias no sean exactos múltiplos enteros del fundamental, simulando así la inarmonicidad típica de señales fisiológicas reales.\n",
|
|
"- Las amplitudes de los componentes fueron ajustadas para reflejar la relación típica en señales UWB. La señal respiratoria es significativamente mayor que la cardíaca, con amplitudes generadas en el intervalo de 1.0 a 2.0 para respiratoria y de 0.01 a 0.1 para cardíaca, manteniendo la proporción de orden de magnitud observada en mediciones reales.\n",
|
|
"- Se mantiene la generación aleatoria de fases y offsets para cada componente dentro de rangos definidos, lo que asegura variabilidad y realismo en la forma de la señal.\n",
|
|
"- Se añadió la posibilidad de introducir ruido aditivo gaussiano con varianza especificada por el usuario. Además, la función calcula automáticamente la relación señal a ruido (SNR) resultante en decibelios, lo que permite evaluar la calidad de la señal generada en condiciones controladas.\n",
|
|
"- Los valores por defecto para la frecuencia de muestreo y la duración de la señal son 200 Hz y 30 segundos, respectivamente, parámetros adecuados para simulaciones de señales cardiorespiratorias con radar UWB.\n",
|
|
"- Se incluyen diccionarios de salida detallados: uno con los parámetros utilizados (frecuencias, amplitudes, fases, offsets, varianza y SNR), otro con los límites de cada parámetro para posibles optimizaciones, y otro con las semillas de inicialización. Además, se entrega metadata con información de la simulación, como número de componentes, frecuencia de muestreo, duración y factor de dispersión aplicado, facilitando la reutilización y análisis posterior de la señal generada."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 8,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"\n",
|
|
"def generar_muestras_cardiorespiratorias(N_card=3, N_resp=2, fs=200, T=30, sigma2=0.01, epsilon=0.05, seed=None):\n",
|
|
" \"\"\"\n",
|
|
" Generador de señales cardiorespiratorias sintéticas para radar UWB.\n",
|
|
"\n",
|
|
" Args:\n",
|
|
" N_card (int): número de componentes de la señal cardíaca (fundamental + armónicos).\n",
|
|
" N_resp (int): número de componentes de la señal respiratoria (fundamental + armónicos).\n",
|
|
" fs (float): frecuencia de muestreo [Hz].\n",
|
|
" T (float): duración de la señal [s].\n",
|
|
" sigma2 (float): varianza del ruido aditivo gaussiano.\n",
|
|
" epsilon (float): factor de dispersión para armónicos (±5% por defecto).\n",
|
|
" seed (int or None): semilla para reproducibilidad.\n",
|
|
"\n",
|
|
" Returns:\n",
|
|
" signal (np.array): señal sintetizada con ruido.\n",
|
|
" params (dict): diccionario con parámetros de frecuencias, amplitudes, fases.\n",
|
|
" bounds (dict): límites de cada parámetro.\n",
|
|
" seeds (dict): centros de los intervalos para inicialización no lineal.\n",
|
|
" metadata (dict): información de la simulación.\n",
|
|
" \"\"\"\n",
|
|
"\n",
|
|
" if seed is not None:\n",
|
|
" np.random.seed(seed)\n",
|
|
"\n",
|
|
" t = np.arange(0, T, 1/fs)\n",
|
|
"\n",
|
|
" # Intervalos de frecuencias fundamentales\n",
|
|
" f_card_interval = np.array([0.6, 1.1])\n",
|
|
" f_resp_interval = np.array([0.15, 0.35])\n",
|
|
" f_card_center = f_card_interval.mean()\n",
|
|
" f_resp_center = f_resp_interval.mean()\n",
|
|
"\n",
|
|
" # Amplitudes\n",
|
|
" A_card_interval = np.array([0.01, 0.1])\n",
|
|
" A_resp_interval = np.array([1.0, 2.0])\n",
|
|
" A_card_center = A_card_interval.mean()\n",
|
|
" A_resp_center = A_resp_interval.mean()\n",
|
|
"\n",
|
|
" # Fases y offsets\n",
|
|
" phase_card = np.random.uniform(0, 2*np.pi, N_card)\n",
|
|
" phase_resp = np.random.uniform(0, 2*np.pi, N_resp)\n",
|
|
"\n",
|
|
" offset_card = np.random.uniform(-0.1,0.1)\n",
|
|
" offset_resp = np.random.uniform(-0.1,0.1)\n",
|
|
"\n",
|
|
" # Frecuencia fundamental\n",
|
|
" f_card0 = np.random.uniform(*f_card_interval)\n",
|
|
" f_resp0 = np.random.uniform(*f_resp_interval)\n",
|
|
"\n",
|
|
" # Generar armónicos con dispersión\n",
|
|
" f_card = [f_card0]\n",
|
|
" for n in range(2, N_card+1):\n",
|
|
" delta = np.random.uniform(-epsilon, epsilon)\n",
|
|
" f_card.append(n*f_card0*(1+delta))\n",
|
|
"\n",
|
|
" f_resp = [f_resp0]\n",
|
|
" for n in range(2, N_resp+1):\n",
|
|
" delta = np.random.uniform(-epsilon, epsilon)\n",
|
|
" f_resp.append(n*f_resp0*(1+delta))\n",
|
|
"\n",
|
|
" f_card = np.array(f_card)\n",
|
|
" f_resp = np.array(f_resp)\n",
|
|
"\n",
|
|
" # Amplitudes para cada componente\n",
|
|
" A_card = np.random.uniform(A_card_interval[0], A_card_interval[1], N_card)\n",
|
|
" A_resp = np.random.uniform(A_resp_interval[0], A_resp_interval[1], N_resp)\n",
|
|
"\n",
|
|
" # Generar señal\n",
|
|
" signal_card = np.sum([A_card[i]*np.cos(2*np.pi*f_card[i]*t + phase_card[i]) for i in range(N_card)], axis=0) + offset_card\n",
|
|
" signal_resp = np.sum([A_resp[i]*np.cos(2*np.pi*f_resp[i]*t + phase_resp[i]) for i in range(N_resp)], axis=0) + offset_resp\n",
|
|
"\n",
|
|
" signal_clean = signal_card + signal_resp\n",
|
|
"\n",
|
|
" # Ruido y SNR\n",
|
|
" noise = np.random.normal(0, np.sqrt(sigma2), len(t))\n",
|
|
" signal_noisy = signal_clean + noise\n",
|
|
"\n",
|
|
" snr = 10*np.log10(np.var(signal_clean)/sigma2)\n",
|
|
"\n",
|
|
" # Preparar diccionarios de salida\n",
|
|
" params = {\n",
|
|
" 'f_card': f_card,\n",
|
|
" 'f_resp': f_resp,\n",
|
|
" 'A_card': A_card,\n",
|
|
" 'A_resp': A_resp,\n",
|
|
" 'phase_card': phase_card,\n",
|
|
" 'phase_resp': phase_resp,\n",
|
|
" 'offset_card': offset_card,\n",
|
|
" 'offset_resp': offset_resp,\n",
|
|
" 'sigma2': sigma2,\n",
|
|
" 'SNR_dB': snr\n",
|
|
" }\n",
|
|
"\n",
|
|
" bounds = {\n",
|
|
" 'f_card': f_card_interval,\n",
|
|
" 'f_resp': f_resp_interval,\n",
|
|
" 'A_card': A_card_interval,\n",
|
|
" 'A_resp': A_resp_interval,\n",
|
|
" 'phase': [0, 2*np.pi],\n",
|
|
" 'offset': [-0.1, 0.1]\n",
|
|
" }\n",
|
|
"\n",
|
|
" seeds = {\n",
|
|
" 'f_card_center': f_card_center,\n",
|
|
" 'f_resp_center': f_resp_center,\n",
|
|
" 'A_card_center': A_card_center,\n",
|
|
" 'A_resp_center': A_resp_center\n",
|
|
" }\n",
|
|
"\n",
|
|
" metadata = {\n",
|
|
" 'N_card': N_card,\n",
|
|
" 'N_resp': N_resp,\n",
|
|
" 'fs': fs,\n",
|
|
" 'T': T,\n",
|
|
" 'epsilon': epsilon\n",
|
|
" }\n",
|
|
"\n",
|
|
" return signal_noisy, params, bounds, seeds, metadata"
|
|
]
|
|
}
|
|
],
|
|
"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
|
|
}
|