{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6ef1e92d",
   "metadata": {},
   "source": [
    "# Example Testcase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "aa8366d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# GPU support\n",
    "\n",
    "# This library can use CuPy for GPU acceleration.\n",
    "# We recommend installing CuPy via conda-forge:\n",
    "# conda install -c conda-forge cupy cuda-version=12.4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a8ca13c",
   "metadata": {},
   "source": [
    "The following example shows an optimization problem, where we place 16 cubic magnets around a ring, and we would like to determine their optimal orientations, such that the magnetic field homogeneity in Bx is optimized. It is important to note, that the field can only be optimized for the set of measurement points, where we estimate the field, here called the sensor points. Using the function sensorring, the psoitions of a ring of sensors can be generated. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a40ce0da",
   "metadata": {},
   "source": [
    "Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "aa79db1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cupy as cp\n",
    "import matplotlib.pyplot as plt\n",
    "from magoptlib.genetic_algorithm import genetic_algorithm_gpu\n",
    "from magoptlib.sh_fitting import precompute_sh_coeffs, get_m_l_vals\n",
    "from magoptlib.sh_gpu import get_orientations_gpu\n",
    "from magoptlib.synthetic_measurement import simulate_magnetic_measurement\n",
    "from magoptlib.pos import magnetring_gpu, sensorring"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "306635e6",
   "metadata": {},
   "source": [
    "Generate synthetic measurement data\n",
    "\n",
    "The library was created to optimize the orientations of real magnets. However, it also offers the possibility to generate synthetic measurement datasets. This is possible via the Magpylib Python library.\n",
    "\n",
    "The sampling is done using the Fibonacci lattice method, which offers an even distribution on the sphere. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c6a5f4a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_theta = 18 # Number of theta for the Fibonacci lattice sampling\n",
    "n_phi   = 36 # Number of phi for the Fibonacci lattice sampling\n",
    "radius = 20.0 # Radius of the measurement sphere\n",
    "n_magnets = 16 # Number of magnets in the Halbach ring\n",
    "\n",
    "measurements = [\n",
    "    simulate_magnetic_measurement(n_theta, n_phi, radius)\n",
    "    for _ in range(n_magnets)\n",
    "]\n",
    "\n",
    "Bx = np.stack([m[0] for m in measurements], axis=0)  # (16, 648)\n",
    "By = np.stack([m[1] for m in measurements], axis=0)  # (16, 648)\n",
    "Bz = np.stack([m[2] for m in measurements], axis=0)  # (16, 648)\n",
    "\n",
    "# DO NOT stack these:\n",
    "theta_values = measurements[0][3]   # (648,)\n",
    "phi_values   = measurements[0][4]   # (648,)\n",
    "\n",
    "assert Bx.shape == By.shape == Bz.shape == (n_magnets, n_theta * n_phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cba4478d",
   "metadata": {},
   "source": [
    "Initialization\n",
    "\n",
    "Since the optimization is done by using the SH coefficients, we first have to map the magnetic field data to spherical harmonics. We first have to initialize the maximum order of the approximation, and therefore get the degree and order of the spherical harmoincs. \n",
    "\n",
    "- $l (l=0, 1, 2, ...)$: is the degree, which determines the complexity of the function in terms of the number of nodal lines\n",
    "- $m (m=-l,...,+l)$ is the order, specifies the azimuthal variation within a given degree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "47196e83",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ============================================================\n",
    "# Initialization\n",
    "# ============================================================\n",
    "\n",
    "max_sh_order = 24 # Maximum spherical harmonic order for the fitting\n",
    "\n",
    "# Get the (l,m) values for the spherical harmonics\n",
    "m_values, l_values = get_m_l_vals(max_sh_order)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a728c7ad",
   "metadata": {},
   "source": [
    "In this example, the initial orientations are given by the theoretical Halbach configuration. \n",
    "\n",
    "```{image} images/halbach_orient.png\n",
    ":width: 50%\n",
    ":align: center\n",
    ":alt: Halbach orientation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c5734087",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial angles for the magnets\n",
    "init_angles = cp.linspace(0, 4*np.pi, n_magnets, endpoint=False, dtype=np.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fa6d9d16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# based on the number of magnets and their initial angles, get their positions and orientations\n",
    "positions_cpu, init_orientations_gpu = magnetring_gpu(n_magnets, init_angles, radius = 40)\n",
    "\n",
    "# Generate the sensor positions on a ring where we want to optimize the field\n",
    "points_of_interest = sensorring(24, radius=20)\n",
    "\n",
    "# Map the magnetic field measurement data to spherical harmonics and compute the SH coefficients\n",
    "sh_coeff_X, sh_coeff_Y, sh_coeff_Z = precompute_sh_coeffs(Bx, By, Bz, theta_values, phi_values, n_magnets, radius, max_sh_order)\n",
    "\n",
    "# Convert the positions, and the SH coefficients to cp arrays\n",
    "positions_gpu = cp.asarray(positions_cpu)           # (8,3)\n",
    "shX_gpu = cp.asarray(sh_coeff_X, dtype=cp.float64)\n",
    "shY_gpu = cp.asarray(sh_coeff_Y, dtype=cp.float64)\n",
    "shZ_gpu = cp.asarray(sh_coeff_Z, dtype=cp.float64)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa2f2a9c",
   "metadata": {},
   "source": [
    "## Genetic Algorithms\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03b0f82f",
   "metadata": {},
   "source": [
    "First, the set of possible angles has to be given that the algorithm can choose from. Naturally, the more possible angles to choose from, the more possible combinations, which we have to take into account when setting the max number of iterations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "8112f0cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Give the possible angles the algorithm can choose from\n",
    "possible_angles_cpu = np.linspace(-15*np.pi/180, 15*np.pi/180, 5, dtype=np.float64)\n",
    "possible_angles_gpu = cp.asarray(possible_angles_cpu, dtype=cp.float64)  # Convert to cp array"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56601f75",
   "metadata": {},
   "source": [
    "Before starting the algorithm, a few parameters have to be set. The population size sets the number of individuals in a generation, i.e. sets the number of potential solutions examined in one iteration. Choosing a too big number could result in performance issues in parallel execution, depending on the specific GPU memory.\n",
    "\n",
    "Setting the num_gens sets the maximum number of iterations. More iterations will result in a longer runtime, but the probability of finding the best optimal configuration is also higher.\n",
    "\n",
    "The mut_rate sets the probability of mutations of the parents for genetic algorithms, it introduces randomness thus helps to avoid premature convergence. The mut_rate is not used for ACROMUSE, since it uses adaptive mutation rates. \n",
    "\n",
    "Since GA is a heuristic algorithm, finding the global optimum is not guaranteed. Therefore, it is recommended to restart the algorithm multiple times and choose the best solution as the final optimal solution. A great rule of thumb is to set the number of restarts to be 10. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "609f866e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best fitness: 0.08083379807215738\n",
      "sum of Bx: 2601.082654525003\n",
      "sum of By: 538.6248074558314\n",
      "sum of Bz: 3.7416517047157973e-06\n",
      "Run 1/1 done, best fitness = 8.083e-02\n",
      "Best angles (rad): [ 0.         -0.26179939 -0.26179939 -0.26179939 -0.13089969  0.26179939\n",
      "  0.26179939  0.26179939  0.         -0.26179939 -0.26179939 -0.26179939\n",
      " -0.13089969  0.26179939  0.26179939  0.26179939]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAHHCAYAAABXx+fLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABYLklEQVR4nO3deVxU9f4/8NcwAzMsgjiMuCGmIIIbimspLpmK+3ZbTcoSI5dWNa3MNSyt7Hcl83tLsa5Laouk3gpLM9RcEELB8KqIYC44LIMLMMx8fn8ok3MBYWCGA8zr+XjM4945c86Z9zlxm9f9nM8iE0IIEBEREdkRB6kLICIiIqptDEBERERkdxiAiIiIyO4wABEREZHdYQAiIiIiu8MARERERHaHAYiIiIjsDgMQERER2R0GICIiIrI7DEBEVKuKi4vRqVMntGzZElevXq217500aRJcXFxw7ty5WvtOIqq7GICIqiAmJgYymcz0cnBwgKurK0JCQvDPf/4TRqPRJt+blJSERYsW4a+//qryMQMHDjSrtVGjRujduze++OILm9RoqaioKKSlpeHrr7+Gt7d3hftt2rQJEydOBACcO3cOzZo1q/TcISEhmDlzZrmfbdy4Eb6+voiIiLC45qtXr2LevHkICgqCi4sLXF1dERwcjHnz5uHs2bMWn4+IpKeQugCi+mT+/Plo3bo1hBC4fv06du7cidmzZ0MIgdmzZ1v9+5KSkrB48WKMGjUKLVq0qPJxLVq0wNtvvw0hBLKzs7F9+3aEh4dDq9XilVdesXqdVZWTk4P3338fM2bMQJ8+fe6775EjR9CrVy8AwNGjR9G7d+9y98vLy0N6ejo2bNiAEydOoG/fvuXu5+rqin/+85945JFH8OOPP2LYsGFVqvno0aMYOXIkbt68iSeffBLTp0/HrVu3cPjwYXz44YdYtWoVDAZDlc5FRHUHAxCRBSZMmIAePXqY3r/++uvw8vLCf/7zH5sEoOry9PTECy+8YHo/Z84cdOvWDcuXL8dLL70EBwdpGn8//fRTGAwGzJ8/v9J9jxw5ghUrVgC4fwBavXo1Fi9eXKXvHzJkCPr06YMPP/ywSgEoPz8fEyZMgEqlwqFDh+Dv72/2+ZkzZxAeHl6l7yaiuoWPwIhqwNHREUIItGzZ0my7EAIff/wxgoKCoFQq0axZM7z44ou4ceOG2X779u3Dww8/DI1GA09PT/Tp0wdLlixBcXExYmJi8OyzzwIAevbsaXqkFRMTY3Gdzs7OCAsLg1arRXZ2tmn7M888Azc3tzL7z5w5EzKZzGybTCbDjBkz8OWXX6JLly5wdnZGx44dsXPnzirXsWnTJgwZMuS+j74AoKioCMnJyaaweb8A9OSTT+L777/H999/X6UannjiCcTFxZndh4r83//9Hy5duoR169aVCT8A0L59e8THx5tt27BhA7p06QKVSoXmzZtjxowZyM/PN9tn4MCB6NSpEw4fPowhQ4bA1dUVrVu3xtKlS02PU69duwZHR0fMmDGjzPeePHkSMpkMa9euNW3766+/MHXqVHh7e0OpVKJz587YtGmT2XGLFi2CTCZDcnIyIiIi4OXlBZVKhZSUFABAdnY2nn/+eTRp0gRKpRIhISHo0aNHmb8F4M4/y5CQEDg7O0OtVmPy5Mm4cuWK2T6W/M3k5eXh5Zdfho+PD5ydnREYGIhXXnnF7PFvWloaJk2ahCZNmsDZ2Rm9evXCjz/+WOZcRFUiiKhSGzZsEADE7t27RWZmprh48aI4duyYeOyxx0SLFi3EuXPnzPafPn26cHBwEE899ZT44IMPxMyZM4VKpRITJkww7ZOUlCSUSqXo1q2beP/998WKFSvEiBEjhIODg7h8+bL4888/xZQpUwQAMX/+fLF27Vqxdu1a8eeff9631gEDBoiOHTuW2f7qq68KAOLWrVumbeHh4cLV1bXMvjNmzBD/+68HAMLJyUloNBoxb948sXz5cuHv7y8UCkWZ6y/P1atXBQDx2WefVbiPr6+vAHDfV3h4eIXHAxAzZsy4bx0XL14UAMS2bdsqrblfv36iZcuWle5X6v333xcAxKBBg8SKFStEZGSkcHJyEj179hTFxcWm/QYMGCAcHR2FUqkUzz//vFi1apUYPHiwACA2btxo2m/kyJFCo9EIvV5v9j3z5s0TTk5OIicnRwhx5962bt1aeHt7izfeeEO8//77YsiQIQKA2Llzp+m4d955RwAQrq6uIiwsTHz00Udi+fLlIisrS9y4cUMEBgYKlUolZs6cKVatWiUmT54sXFxcyvwtREVFCQBi9OjRYuXKlWLOnDnCw8NDhISECKPRaNqvqn8zBQUFolOnTsLJyUlMnz5dfPjhh+K5554TLi4uYsOGDUIIIVJTU0Xjxo1Fu3btxMKFC0VUVJTo2bOnkMvl4sSJE1X+Z0RUigGIqApKA1B5r0mTJon09HTTvgcPHhQAxOeff252js2bNwsAIikpSQghxJw5c4SDg4MoKCgw2+/kyZPixo0bZt977NixKtdaXgDS6/WiQ4cOolu3bmbbLQ1Affv2FdnZ2aZtp06dEgDEu+++W2ldsbGxAoC4cOFChfukpKSIxMREMX36dBEeHi4SExPFv/71L9GtWzeRmJgoEhMTRUZGRoXHVyUACSGEn5+feO211yrdr0mTJmL06NGV7ieEEFqtVqhUKjF69GizELBt2zYBQKxdu9a0bcCAAUKtVos//vjDtE2v14tmzZqJoUOHmrZt3bpVABB79uwxbTMajcLX19csTE+bNk1oNBpx+fJls5qGDh0qgoODTe9LA9C6devK1L9s2TIBQPz4449m2//3byEjI0MoFArx9ttvm+136NAhAUB89913pm1V/ZtZuHChACD+85//mJ0zKyvLFG4eeeQR0aFDB9P/NkrvWVBQkBg3blyZ6yGqDPsAEVngo48+gp+fH4A7j2lSU1Oxdu1a9O3bFydOnEDz5s2xbds2uLu745FHHkFWVpbp2K5duwIAEhMT0bVrV8jlchiNRvz8888YO3asab9OnTrVuM6SkhJkZWWZ/nPNmjX473//iz179tTovN27d4eXl5fpfWBgIGQyGTIyMio99urVq1CpVGjdunWF+wQFBQEA/vzzT7zxxhsIDg7Gt99+ixEjRiA4OLhGtd/L39+/SkPwdTod3N3dq3TOuLg4FBYWYtasWWaPjCZNmoSWLVvi+++/N+uX1axZM3Tp0sX0XqFQoH379mb3csyYMXB3d8fmzZsRFhYGADh48CAyMjLw8ccfA7jzuHXHjh0YP3686Z93qb59+2L58uUoKSmBQvH3v+4nTJhQpv4dO3agW7duGDp06H2v85tvvoHBYMCjjz5q9l0+Pj5o1KgREhMTzf6eq/I38/XXXyMkJATDhw83+66WLVuiZcuWyMnJwd69e/H2228jNzcXubm5Zte4d+/e+9ZMVB4GICIL9OvXz6wT9MSJEzFx4kR07NgRa9aswfLly3HmzBnodLoKf+gzMzMBANOnT8fGjRsxbtw4+Pr6om/fvhg8eDCefPJJuLq61qjOtLQ0+Pj4mN43adIEu3btqvTHzVIODg5wdnZGYWFhpftqtVp4enqW258EAK5cuYKSkhIUFRXh+PHjaNeuHbKysrBnzx68/vrryMrKgqurKzw9PWtcd+PGjaHVaivdr1GjRtDpdFU6Z3p6OgCYAnIpmUyGtm3b4vz585Wew83NzSwYODs7Y9KkSdi2bRtu374NZ2dnbN68GV5eXhgxYgSAO/12cnNzsX79eqxfv77c816+fNns76E8Z86cwbhx4yqt8cyZMxBCoHPnzuV+Xvr3XZHy/mbOnTuHSZMmVXjM2bNnIYTAkiVLsGTJknLPSWQpBiCiGgoKCoK3tzeSk5MBAEajES1atMC6devK3b99+/YAgDZt2iA1NRVbt27FL7/8gvj4eGzduhUffPABjh8/Xm7n5Kpq3bo1oqOjAQCffPIJfvzxR9y+fbva57ufigJNefvdb98+ffqY/fiX3icAePzxxwEA4eHh1eoE/r+EEFWqOyAgAAkJCVXev/Tc1VXed0yZMgXr169HbGwsJk6ciO3bt+Pxxx+Ho6MjAJg6TT/99NN49NFHyz3vvS0wFdHr9VCpVJXuZzQaIZfL8d1335X7eWVBCyh7nZXd39JrfP311zFgwIBKz09UFQxARDUkhEBRUZHpEUObNm3w66+/YuDAgZWGmMaNG+OFF14wPRr55JNPMGPGDHz77bd4+umnIZfLAcDiiRYbNWqEUaNGAbgz4qhfv3546qmnsH//ftPcOgCgUqlw+/Zti37gq0utVpcZDXWv9evX49atW/jqq69w8+ZNTJ06FUlJSdi0aRNWrlwJoGo/rlWRm5uLpk2bVrrfiBEjsHDhQsTGxpo91ilPmzZtANxprWjbtq1puxACZ8+erfYjvNDQUPj6+mLz5s1o1KgRrl+/jilTppg+12g0cHV1hdFoNP0zr45mzZrh0qVLle7Xpk0bGAwGdOjQoUxrV3X5+vri9OnT9/0cAJRKZY2ukehebDckqqGvvvoKeXl5GDJkCABg3LhxKCwsxNy5c8vse+3aNVy8eBEAsH37dly/ft3s89LHa0VFRQBgGi5+5syZatfn5uaGXbt2oXHjxhgzZgwuXLhg+qx169YwGo2m1ivgziOVEydOVPv7KuLt7Y2bN29WOKv14MGDMWrUKNy4cQOTJk3CqFGj4OjoaNo+atQoUz+qmkpLS6tSAIqMjIRarUZkZCRSU1PLfJ6dnY1//OMfAIBHHnkESqUSH3/8sVlg3bJlCy5fvozRo0dXq1aZTIannnoKP/zwAz755BMEBgaiZ8+eps/lcjlGjRqFrVu34tChQ2WOP378eJW+p3///ti/f7/ZUPbCwsIyM12PHTsWMpkMr776apkJIAsKCpCWlmbJ5QEARo0ahePHjyMuLs5s++XLl3HmzBk0b94cvXv3xv/7f/+v3Jm3q3qNRPdiCxCRBb755hscP34cQggUFBQgMTERO3bsQK9evfD8888DAIYPH44nnngCa9euRWJiIsLCwuDq6oqEhATs3LkTX3/9NVq3bo1NmzZh2rRpeOyxxxAUFIT8/HzExMSgadOmpr4YPXv2hLu7O+bMmYOMjAzcunULo0aNqnC244q0atUK33//PUJDQzFixAgcOnQIjRs3xqOPPop33nkHjz/+OCIiInD16lV8/vnn922pqa7SeXz27duHp556qsL9SmdYBoBDhw5V+FinVGpqqlk4OXfuHHbs2AEA5fYrOXv2LC5cuFDpTNTAnUdHO3bswNixY9G9e3c89thj6NatG4qLi/HHH38gNjbW9GhRrVZj0aJFmD9/PgYPHowRI0YgPT0dn332Gbp3747nnnuu0u+ryNNPP413330Xu3fvxrvvvlvm8/feew/79+/HwIED8eSTT6JLly7Q6XT46aefkJ6ejsuXL1f6HXPmzMGOHTsQGhqKKVOmwGAw4MsvvyzTdykwMBBvvPEGoqKi0LVrV4wfPx5qtRonT57EN998g6ioKAQEBFh0ffPnz8fXX3+NUaNGYerUqfD398fp06exdetWfPnll2jfvj2io6MxcOBABAcHY8qUKWjfvj2uXbuGXbt2wcXFBb///rtF30nEYfBEVVDeMHilUinat28vFixYIG7evGm2v8FgEP/85z9F9+7dhbOzs2jcuLHo1auXWL58uWnulkOHDoknn3xStGjRQjg6OgofHx8RHh5uNqReCCH27NkjOnbsKJRKpWjVqpXYv3//fWutaB4gIYT47rvvhIODgxg0aJBpXpqdO3eazh8UFCSio6NFZGRkucPgyxti7urqet+5ee7VsWNHs+Hb/+vs2bPC29vb9F6j0Yjz58/f95ylQ7vLe5XnvffeEzKZTFy7dq1KNQshRHp6uoiMjBR+fn5CqVQKNzc30b17d/HWW2+Vqe9f//qX6Nixo3BychLe3t4iMjLS9M+8VEX/jEaOHCl8fX3LraFHjx7CwcFBXLx4sdzPs7KyxPTp04WPj4/p72ns2LHi66+/Nu1Teq/uHZZ+rx9++EF06tRJODo6ilatWok333xTPPfcc0KpVJbZd9OmTeLBBx8Ubm5uws3NTQQHB4sFCxaIrKws0z6W/M1cunRJPPPMM8Lb21solUoREBAg3njjDbNaT58+LZ588knRrFkz4eTkJB544AHx+OOPi59//rnc6yG6H5kQNeixR0RkgRUrVmDp0qW4ePEi1Gp1rX+/0WhE165d0apVK/znP/+p9e+vj5544gn8/vvvplFuRA0F+wARUa0pXYl9+fLlknz/pk2bcOrUKUkXhK3L/rezvU6nw969e9G/f3+JKiKyHbYAEVGtWrp0KZYvX47k5GSzoe62lp+fj86dO6N9+/acOK8cN27cQGBgIEaPHo2AgADcunULGzZsQGZmJo4dO2aVCTqJ6hIGICKqVXq9Ht26dYNMJkN8fDw8PDxs/p1GoxFjxozBL7/8gpMnT6Jdu3Y2/876pqioCGPHjkViYiJyc3Ph6uqKvn37YvHixWajzogaCgYgIiIisjvsA0RERER2hwGIiIiI7A4nQqyA0WjEX3/9hUaNGtl8iQAiIiKyDnF3otoWLVrcd6FcBqAK/PXXX1Zbd4iIiIhqV2ZmJlq1alXh5wxAFWjUqBGAOzfQ3d1d4mqIiIioKnQ6HXx8fEy/4xVhAKpA6WMvd3d3BiAiIqJ6prLuK+wETURERHaHAYiIiIjsDgMQERER2R0GICIiIrI7DEBERERkdxiAiIiIyO4wABEREZHdkTQACSGwdOlS+Pj4QKlUIjg4GHFxcfc9JiYmBv3794dGo4GHhwdCQ0MRHx9f4f5xcXFwcnLC6tWrrVw9ERER1VeSBqDo6GisXr0aa9asQXJyMsLCwjBmzBikp6dXeMyRI0cwbtw47Nq1CwcOHICfnx9GjBiBy5cvl9k3OTkZkydPhpOTky0vg4iIiOoZSQPQunXrMHfuXIwdOxYBAQGIioqCv78/YmJiKjxm7dq1eO2119C7d2907doV69atQ2FhIQ4fPmy2X1ZWFkaNGoVPPvkEXl5eNr4SIiIiqk8kC0CFhYVISUlBSEiI2fZ+/fohISGhyufJy8uDXq+Hp6enaZtOp8OIESPwyiuvYOLEiVU6T1FREXQ6ndmLiIiIGibJApBWq4UQosw6W2q1GteuXavyeRYuXAh/f3+EhoYCAPR6PSZMmICBAwfilVdeqfJ5oqKi4OHhYXpxJXgiIqKGS/JRYApF9ddjXbVqFbZs2YJt27ZBLpcDAObPn49GjRpZ3Ol5/vz5yM/PN70yMzOrXdf96A1GZGhvQnujyCbnJyIiospJFoDUajVkMhlyc3PNtmu1Wmg0mkqPX7p0KVasWIG4uDgEBwebtp89exZ79uyBi4sLVCoVVCoVMjIyMGfOHDRr1qzC8ymVStPK77ZcAf6Vr5IwYOV+fJt4ySbnJyIiospVv/mlhlQqFYKCghAfH4+HH37YtP3w4cMYOXJkhccVFxfjhRdewL59+3Dw4EEEBASYfR4dHY2CggKzbQ8//DDCw8Px7LPPWvciqqGVpwsAIDPnlsSVEBER2S/JAhAARERE4M0330S3bt0QEBCA7du3IzU1Fdu2bQMALFmyBEuWLMG5c+fg6+sL4E6YuXLlCr766isolUpcuHABwJ1Haa1atULLli3LfI+joyOaNm0Kf3//Wru2ivg0cQYAZObelrgSIiIi+yVpAJo1axZycnIQGRmJ7OxsBAYGIjY21hRUjEYjDAYDhBCmY0onPezdu7fZuXx9fU1hqC5r3YQtQERERFKTiXvTBZnodDp4eHggPz/fqv2BLly/iYGr9kPl6IDTS4ZDJpNZ7dxERET2rqq/35KPArM3LRo7QyYDCvVGZHMkGBERkSQYgGqZk8IBzd1VAIDMHPYDIiIikgIDkAR87vYDysplPyAiIiIpMABJoDQAXdQyABEREUmBAUgCPqVzAbEFiIiISBIMQBIwzQXEPkBERESSYACSQOkjMLYAERERSYMBSAKlkyFezi+E3mCUuBoiIiL7wwAkAY2bEk4KBxiMApfzCqUuh4iIyO4wAEnAwUGGVp6la4LxMRgREVFtYwCSiA9XhSciIpIMA5BEWrMjNBERkWQYgCRSOhT+IofCExER1ToGIInwERgREZF0GIAkwvXAiIiIpMMAJJHSFqDrN4pxq7hE4mqIiIjsCwOQRDxcHOGuUgDgkhhERES1jQFIQqYlMdgPiIiIqFYxAEmIq8ITERFJgwFIQlwVnoiISBoMQBLiZIhERETSYACSUCv2ASIiIpIEA5CE7p0MUQghcTVERET2gwFIQqUrwt8sNiD3ll7iaoiIiOwHA5CEVI5yNG2kBMDHYERERLWJAUhipR2hLzIAERER1RoGIIn5cCQYERFRrWMAkpiPJ+cCIiIiqm0MQBJrxVXhiYiIah0DkMRacy4gIiKiWscAJLHSPkCX8m7DYORcQERERLWBAUhizdxVcJTLoDcIXNEVSl0OERGRXZA0AAkhsHTpUvj4+ECpVCI4OBhxcXH3PSYmJgb9+/eHRqOBh4cHQkNDER8fb7bPypUr0atXL3h6ekKtViMsLAynTp2y5aVUm9xBhhaNSztC8zEYERFRbZA0AEVHR2P16tVYs2YNkpOTERYWhjFjxiA9Pb3CY44cOYJx48Zh165dOHDgAPz8/DBixAhcvnzZtM/Ro0fx3HPP4eeff8YPP/wAuVyOsLAwFBbWzRaWe5fEICIiItuTCQkXoercuTMmT56MefPmmbZ16dIF48ePx+LFi6t0Dr1eD1dXV2zduhUTJkwod5+MjAy0adMGCQkJ6N69e5XOq9Pp4OHhgfz8fLi7u1fpmOqa/81JbDl6EbMH++HVoQE2/S4iIqKGrKq/34parMlMYWEhUlJSEBISYra9X79+SEhIqPJ58vLyoNfr4enpWeE+165dA4D77lNUVISioiLTe51OV+Uaasqnyd1HYLmcC4iIiKg2SPYITKvVQghRJp2p1WpTYKmKhQsXwt/fH6GhoeV+LoTAokWLMHjwYDzwwAMVnicqKgoeHh6ml4+PT5VrqCk+AiMiIqpdko8CUyiq3wi1atUqbNmyBdu2bYNcLi93n9mzZyM5ORkxMTH3Pdf8+fORn59vemVmZla7LktxOQwiIqLaJdkjMLVaDZlMhtzcXLPtWq0WGo2m0uOXLl2Kjz/+GHFxcQgODi7zucFgQGRkJPbu3Yv9+/dX2qKjVCqhVCotugZrKZ0M8aquCIV6A1SO5Yc5IiIisg7JWoBUKhWCgoLKDGE/fPgwunXrVuFxxcXFmDp1KtavX4+DBw+iZ8+eZfbR6XQYPXo0jh8/jkOHDqFdu3ZWr9+aPF0c4ep0J/RksR8QERGRzUn6CCwiIgKrVq1CbGws0tLSsGzZMqSmpiI8PBwAsGTJEigUCmRkZJiOefjhh/Hbb7/hq6++glKpxIULF3DhwgVkZWUBAAoKCtC3b19cv34dX375JQoLC037XLlyRZLrrIxMJuNjMCIiolok2SMwAJg1axZycnIQGRmJ7OxsBAYGIjY2Fv7+/gAAo9EIg8GAe0fql7YY9e7d2+xcvr6+uHDhArRaLVJTUwEAnTp1MttnwIAB2L9/vw2vqPpaebrgzysFyGJHaCIiIpuTdB6guqw25wECgCXfp2L9wXREhLbFghGBNv8+IiKihqiqv9+SjwKjO0rnArqoZQsQERGRrTEA1RGmuYDYB4iIiMjmGIDqCFMnaPYBIiIisjkGoDqileedR2C6whLk39JLXA0REVHDxgBUR7gqFfBycwLAx2BERES2xgBUh7TimmBERES1ggGoDuFkiERERLWDAagO8bnbDygzh8thEBER2RIDUB1SuijqRT4CIyIisikGoDqEj8CIiIhqBwNQHVI6GWJW7m0YjVyhhIiIyFYYgOqQ5o1VcJABxSVGZN8okrocIiKiBosBqA5xlDuguUdpR2g+BiMiIrIVBqA6hh2hiYiIbI8BqI4pXRWeQ+GJiIhshwGojuGq8ERERLbHAFTHcFV4IiIi22MAqmMYgIiIiGyPAaiOKe0DdFlXiOISo8TVEBERNUwMQHWMxk0JlaMDhAD+ymNHaCIiIltgAKpjZDIZWrEjNBERkU0xANVBpavCcy4gIiIi22AAqoNamzpC8xEYERGRLTAA1UGlI8EytDclroSIiKhhYgCqgzq19AAAHE3P4arwRERENsAAVAd1b+0JVyc5tDeLkXpZJ3U5REREDQ4DUB3kpHBA33ZqAMCvZ7IlroaIiKjhYQCqo/r7awAAv/2XAYiIiMjaGIDqqND2dwJQQkYubhaVSFwNERFRw8IAVEe1UbvAp4kz9AaBw+e0UpdDRETUoDAA1VEymYyPwYiIiGxE8gAkhMDSpUvh4+MDpVKJ4OBgxMXF3feYmJgY9O/fHxqNBh4eHggNDUV8fLzZPvn5+QgPD4enpydcXV0xfPhwnD9/3paXYnWhdwPQgf9el7gSIiKihkXyABQdHY3Vq1djzZo1SE5ORlhYGMaMGYP09PQKjzly5AjGjRuHXbt24cCBA/Dz88OIESNw+fJl0z7PPfccUlJSsGfPHhw+fBgODg4YNWoUDAZDbVyWVTzop4bcQYb06zeRyWUxiIiIrEYmhJB0pr3OnTtj8uTJmDdvnmlbly5dMH78eCxevLhK59Dr9XB1dcXWrVsxYcIEZGdno1mzZjh06BB69+4NALh+/Tq8vb2xd+9eDBo0qNJz6nQ6eHh4ID8/H+7u7tW7OCuYtPYQjmfkYvn4Tniqt69kdRAREdUHVf39lrQFqLCwECkpKQgJCTHb3q9fPyQkJFT5PHl5edDr9fD09AQAJCYmwmg0onv37qZ9vLy8EBAQUOF5i4qKoNPpzF51QelosAOcD4iIiMhqJA1AWq0WQogyCU2tVuPatWtVPs/ChQvh7++P0NBQAEB2djacnZ3h6OhY5fNGRUXBw8PD9PLx8bHwamyjNAAdOqtFicEocTVEREQNg+R9gABAoVBU+9hVq1Zhy5Yt2LZtG+RyebXPOX/+fOTn55temZmZ1a7Jmjq39EBjF0cUFJUgKTNP6nKIiIgaBEkDkFqthkwmQ25urtl2rVYLjUZT6fFLly7FihUrEBcXh+DgYNN2jUaDmzdvQq/XV/m8SqUS7u7uZq+6QO4gw0N+XgD4GIyIiMhaJA1AKpUKQUFBZYawHz58GN26davwuOLiYkydOhXr16/HwYMH0bNnT7PPg4ODIYTAoUOHTNtycnKQlpZ23/PWVQM4HJ6IiMiqJH8EFhERgVWrViE2NhZpaWlYtmwZUlNTER4eDgBYsmQJFAoFMjIyTMc8/PDD+O233/DVV19BqVTiwoULuHDhArKysgAATZs2xbhx4zBz5kwcOXIEycnJiIiIgJ+fHwYOHCjFZdZI//Z3WoCSs/KQd6tY4mqIiIjqv+p3vrGSWbNmIScnB5GRkcjOzkZgYCBiY2Ph7+8PADAajTAYDLh3tH5pi1HpEPdSvr6+uHDhAgDg888/x6xZszBs2DAUFxejf//+2L17d436G0mluYcz/Ju64b/XbiD+7HWM6tJC6pKIiIjqNcnnAaqr6so8QKWW7krF5/HpeKyHD96b1EXqcoiIiOqkejEPEFWdaT6g/2aDmZWIiKhmGIDqiV5tmsBJ4YDL+YU4e+2G1OUQERHVawxA9YSzkxy9H2gCgKPBiIiIaooBqB4xrQ7P+YCIiIhqhAGoHikdDn8kXYtCff1Z1Z6IiKiuYQCqRwK8G8HbXYlCvRHHL+RWfgARERGViwGoHpHJZOjv//doMCIiIqoeBqB6pr8/1wUjIiKqKQageqa/vwYyGfDnlQJc0xVKXQ4REVG9xABUzzRxdULnlh4AOByeiIiouhiA6iE+BiMiIqoZBqB6qHQ+oPiz12E0clkMIiIiSzEA1UPdfT3h6iRHzs1ipPylk7ocIiKieocBqB5ylDugb7u7j8E4HJ6IiMhiDED11IC7s0L/yn5AREREFmMAqqce8rsTgJIy86A3GCWuhoiIqH5hAKqnHvByRSOVAsUlRpy5WiB1OURERPUKA1A9JZPJ0KXVnfmAkrPyJa6GiIiofmEAqsc6t2wMgAGIiIjIUgxA9VhXUwtQnrSFEBER1TMMQPVY57sBKO1KAQr1BomrISIiqj8YgOqxlo2d0cTVCSVGgT+vsCM0ERFRVTEA1WPmHaHzpC2GiIioHmEAque6tORIMCIiIksxANVzXVo1BgCcZAAiIiKqMgageq60I/R/rxXgVnGJxNUQERHVDwxA9Zy3uwre7koYBbgyPBERURUxADUApY/B2A+IiIioahiAGoC/O0LnSVsIERFRPcEA1ACU9gNiR2giIqKqYQBqAEofgZ2/fhO6Qr20xRAREdUDDEANQBNXJ7TydAYAnGIrEBERUaUkDUBCCCxduhQ+Pj5QKpUIDg5GXFxcpccZDAbs27cPTk5OSEpKKvN5Tk4Opk+fjhYtWsDNzQ09e/bEzp07bXAFdYdpRuhLDEBERESVkTQARUdHY/Xq1VizZg2Sk5MRFhaGMWPGID09vcJjDhw4AIVCgcGDB0OvL/9xz2OPPYZTp07hq6++wtGjRzFy5EhMnDgRR44csdWlSI4TIhIREVWdpAFo3bp1mDt3LsaOHYuAgABERUXB398fMTExFR7To0cPnD59Gnv37i33c71ej59//hlvvPEG+vfvj6CgICxatAht2rRBfHy8ja5EeqUjwf7gSDAiIqJKKaT64sLCQqSkpCAkJMRse79+/ZCQkFDhcS4uLujQoQNUKlW5nzs6OqJ79+54//33ERgYCD8/P1y8eBFXrlzB4MGDKzxvUVERioqKTO91uvo1qWCnu4/AsnJvI+dmMZq4OklcERERUd0lWQuQVquFEALu7u5m29VqNa5du1ajc+/evRvZ2dkICAjAgAEDMGjQIPz73/9Gt27dKjwmKioKHh4eppePj0+Naqht7ipHtPVyBcD5gIiIiCoj+SgwhcL6jVDz5s3Dgw8+iL/++gtTp06Fh4cH5s6diwsXLlR4zPz585Gfn296ZWZmWr0uW+N8QERERFVT4wCk1+uRmJiIGzduWHScWq2GTCZDbm6u2XatVguNRlPten777Tf8+9//xj//+U94e3sjPDwcR48ehbu7OxYtWlThcUqlEu7u7mav+sa0JAZHghEREd2XxQHo9ddfx7BhwwAAt2/fRo8ePRASEoLWrVvj2LFjVT6PSqVCUFBQmY7Jhw8fvu+jqsrk5+dDCIHCwkLTNoVCAW9vb+Tk5FT7vPWBaSg8H4ERERHdl8UBaOfOnYiIiDD99ytXruDMmTOIjIzEvHnzLDpXREQEVq1ahdjYWKSlpWHZsmVITU1FeHg4AGDJkiVQKBTIyMgwHVNSUoK8vDxTJ+WCggLk5eWhpKQEANC/f3+0bNkS48ePx4EDB5CSkoKVK1fihx9+wLPPPmvp5dYrHVu4w0EGXNUV4aqusPIDiIiI7JTFAejSpUvo2rUrACA+Ph6PPPII/Pz8MGXKFCQmJlp0rlmzZuG1115DZGQkOnfujO3btyM2Nhb+/v4AAKPRCIPBACGE6Zj4+Hh4enqaaggNDYWnp6epJcnDwwO//PILmjVrhkcffRR9+vTBN998g507d2L8+PGWXm694uKkgH/TRgDYD4iIiOh+ZOLedFEFfn5+iI6OxrBhw9ChQwfMmjULM2bMQEJCAoYNG4br16/bqtZapdPp4OHhgfz8/HrVH+j17X9gR0IWZg/2w6tDA6Quh4iIqFZV9ffb4hagZ599Fk8//TT69OmDrKwsTJo0CQDw3XffISgoqPoVk1V05ZIYRERElbJ4DPqbb74JjUaDkydP4sMPP4S3tzcMBgMyMjLw6quv2qJGskDne5bEEEJAJpNJWxAREVEdVK1JeEo7QZeSy+X44osvrFIQ1UyHZo2gcJBBe7MYl/Juo5Wni9QlERER1TkWPwL78ssvMWvWLNP7F154AR4eHhgyZAiuXLli1eLIcipHOTo0Z0doIiKi+7E4AK1evRodOnQAAOzduxdffPEFPvjgAzg7O+O1116zeoFkuc4tGwNgPyAiIqKKWPwI7MyZM6ZFRX/++Wc8/PDDeP755xESEoKwsDCrF0iW69LKA1uOckJEIiKiiljcAtSoUSPcvHkTwJ05efr06QPgzszOpdtJWn/PCJ0PC2c5ICIisgsWtwANHz4cs2fPRr9+/XD48GGsXbsWAHDgwAG0bdvW6gWS5dp7N4KTwgEFhSXI0N5Cm7urxBMREdEdFrcAffDBB2jTpg3+85//YOXKlejUqRMMBgPWrl2Lxx9/3BY1koUc5Q4Ian5n8qc/+BiMiIioDItbgDw9PbF582azbXK5HElJSdaqiaygaysPJGXm4WRWPsYGt5S6HCIiojrF4hYgAPj1118xbdo0DBs2DFlZWQCAFStWYP/+/dasjWqgdEJEjgQjIiIqy+IAtHHjRgwbNgw3btzAvn37cOvWrTsncnBAVFSU1Quk6ildEuPUpXwYjOwITUREdC+LA9D777+PTz/9FFu2bIFcLjdtHzp0qMWrwZPttNW4wcVJjlvFBpzPviF1OURERHWKxQHo/Pnz6Nu3b5ntcrkcBQUFVimKak7uIEOnFndagf7gjNBERERmLA5Abdq0walTp8ps37Nnj2mGaKobSucDOsmRYERERGYsHgU2Z84czJo1C3q9HgBw9OhRbNmyBVFRUfj000+tXiBVX+fSCRHZEZqIiMiMxQFo6tSpMBgMmDdvHoqKijBlyhR4e3tj5cqVeOaZZ2xQIlVXl7sjwVL/0kFvMMJRXq1Bf0RERA2OxQEIAKZNm4Zp06bh+vXrMBqNaNq0qbXrIitoo3aBytEBhXoj/sq7DV81Z4QmIiICqjkPUCkvLy+GnzpMJpPBy00JALh+o0jiaoiIiOoOi1uAtFotFi5ciEOHDiE/v2zfkvPnz1ulMLIOTSMlsnJvI7ugWOpSiIiI6gyLA9AzzzyDlJQUPPHEE9BoNJDJZLaoi6yktAUomy1AREREJhYHoP379+Onn34qdy4gqns0je4+AitgACIiIiplcR8gjUYDR0dHW9RCNsA+QERERGVZHIBmz56Njz76yBa1kA2UtgBlswWIiIjIxOJHYJ999hlOnz6N1NTUcvv/nDhxwiqFkXVo3JwAsAWIiIjoXhYHoEmTJrHjcz1iagFiACIiIjKxOAAtWrTIBmWQrZj6ABUUQwjB8EpERIRq9AGSy+W4ePFime0//vgj2rZta5WiyHpKA9BtvQE3iw0SV0NERFQ3WByAhBDlbvf29saVK1dqXBBZl6tSARcnOQAOhSciIipV5UdgX3zxBYA7yyt8/fXXUKvVps+KioqwceNGBAYGWr9CqjEvNyUu5txC9o0itPHiemBERERVDkAvvfQSgDstQIsWLYKDw9+NRyqVCh07dkR0dLT1K6Qa0zS6E4DYAkRERHRHlR+B5ebmIjc3FwMGDMCZM2dM73Nzc3H58mXs3bsXXbt2tbgAIQSWLl0KHx8fKJVKBAcHIy4urtLjDAYD9u3bBycnJyQlJZW7T2ZmJqZPn47WrVvDyckJgwcPtri+hsDr7lB4jgQjIiK6w+JRYPv27bNqAdHR0Vi9ejXWr1+PDh06ICYmBmPGjEFqaioeeOCBco85cOAABgwYcN/znjlzBgMHDsSECROwadMmeHt7o6CgwKq11xdcDoOIiMhclQPQ4MGDsXXrVjz++OP33e+XX36xqIB169Zh7ty5GDt2LAAgKioKu3fvRkxMDBYvXlzuMT169MDp06dx6dIlDBkypNx9pk+fjpdffhlz5861qJ6GiAuiEhERmatyAPL19YVCoajWY66KFBYWIiUlBSEhIWbb+/Xrh4SEhAqPc3FxQYcOHaBSqcr9PCMjA/v370dwcDCCgoJw+fJltG7dGvPnz68wwBUVFaGo6O+AoNPpqnFFddPfy2EUS1wJERFR3VDlABQdHQ0XFxerrgOm1WohhIC7u7vZdrVajePHj1f7vMnJyZDL5XBzc8Onn34KNzc37NixA08++SSaNGmCoUOHljkmKiqqwhan+o4LohIREZmrcifoRo0amU2AuHPnTty8edMqRSgUFndFuq/8/Hy4uLhg6dKlCA0NRffu3fHuu+9iwIAB2LhxY7nHzJ8/H/n5+aZXZmamVWuSEhdEJSIiMlflAPS/EyA+8cQTuHTpUo2+XK1WQyaTITc312y7VquFRqOp9nnd3d1x8+ZNFBebP/Lx9/fH9evXyz1GqVTC3d3d7NVQaO5pAapoIksiIiJ7YvFM0KWs8UOqUqkQFBSE+Ph4s+2HDx9Gt27dqn3ekJAQCCGwd+9es+2nTp1C+/btq33e+qr0EVhRiREFRSUSV0NERCS9agcga4mIiMCqVasQGxuLtLQ0LFu2DKmpqQgPDwcALFmyBAqFAhkZGaZjSkpKkJeXZ+qoXFBQgLy8PJSU3Plxb9myJZ566ilMmzYN3333HU6ePIk33ngDJ0+exOzZs2v/IiXm7CSHm/LOY0YOhSciIrKgE7RMJiuzkrg1VhafNWsWcnJyEBkZiezsbAQGBiI2Nhb+/v4AAKPRCIPBYNbiFB8fj0GDBpneh4aGArgzR9HAgQMBAP/617+wePFivPTSS7h69Sp69uyJX375xXRee+Pl5oQbRSXILihCW42b1OUQERFJSiaq+CzLwcEB7u7upiUw8vLyzN6XysnJsX6VEtDpdPDw8EB+fn6D6A/0j08P4diFXEQ/2R0juzSXuhwiIiKbqOrvd5VbgDZs2GCVwkgapskQCwolroSIiEh6VQ5ApX1yqH4yLYdxg5MhEhERSd4JmmrH3y1A7ARNRETEAGQn/m4BYgAiIiJiALITXA6DiIjobwxAdsLLzQkAH4EREREB1QhAFy9ehNFoLLP91q1buHr1qlWKIuu7txM0l8MgIiJ7Z3EAeuCBB5CVlVVm+9GjRxESEmKVosj6Sh+BFRuM0N3mchhERGTfLA5AFbUeuLq6llnUlOoOlaMcjVR3Zj3IZj8gIiKyc1WeB2jJkiUA7ix/sXr1ajRu3Nj0WVFREXbs2IGePXtavUCyHo2bEgWFd5bD8GvK5TCIiMh+VTkAffvttwDutAD9+OOPcHJyMn2mUqnQv39/LFy40PoVktV4NVLi/PWbHAlGRER2r8oBKDExEQDw7LPPYtWqVVCr1TYrimxDw8kQiYiIAFSjD9Cnn35qtgBqTk4Ovv32W5w7d86qhZH1cTJEIiKiO6rcAlRq5syZuHTpEvbs2YO8vDx06dIF2dnZkMlk2LlzJ4YNG2aLOskKOBcQERHRHRa3AMXFxWHWrFkA7vQLksvlyM3NxapVq/DOO+9YvUCyHrYAERER3WFxALp27Rr8/f0BAAcPHsSQIUPg4uKCYcOGITU11eoFkvX8vRwGV4QnIiL7Vq2JEBMTE1FSUoK4uDg89NBDAIDr16/DxcXF6gWS9XBFeCIiojss7gP00ksvYcqUKfDy8kJRUREmTpwIANiyZQu6du1q9QLJekofgWlvFsFoFHBwkElcERERkTQsDkARERFo06YNTp48ibFjx8LDwwMGgwEuLi548803bVEjWYn6bidovUEg/7Yenq5OlRxBRETUMFkcgABg6NCheOSRR3D58mUYDAbI5XKsWLHC2rWRlSkVcng4OyL/th7XbxQxABERkd2yuA+QTqfDM888AxcXF7Ru3do0/8/o0aOxfPlyqxdI1sWh8ERERNUIQK+99hpOnjyJb7/9FgrF3w1IU6ZMwbZt26xaHFlfaT8gLohKRET2zOJHYLt378b27dvx0EMPQSb7uxNt165dORt0PcCRYERERNV8BHbvSvClrl69arZAKtVNf0+GyLmAiIjIflkcgEJDQ7Fp0ybTe5lMhpKSErz33nsYMGCAVYsj62MLEBERUTUega1atQoDBgzAiRMnUFJSgnnz5iElJQXXr19HfHy8LWokK+JyGERERNVoAQoKCsLJkyfRu3dvDBs2DIWFhRg/fjwSExMRGBhoixrJijRsASIiIqrePEDNmjXD4sWLrV0L1YK/1wNjACIiIvtlUQtQcXExtFotAEAIgTlz5sDLywsajQbLli2zSYFkXX8vh1EMo1FIXA0REZE0qtwC9PPPP2PixIkoKCjAE088gf79++Pf//433nrrLeh0Orz//vtwdXXFK6+8Yst6qYZKl8MwGAVybxVDfbdFiIiIyJ5UOQC98847GDx4MJ577jl89tlnePXVV7Fu3TpMnjwZANCkSRN8+umnDEB1nKPcAZ4ujsi9pcf1GwxARERkn6r8CCwxMRELFizAyJEj8dlnn+H27dvo3Lmz6fNBgwZxIsR6gkPhiYjI3lU5ABUWFsLLywsAoFaroVQq4ezsbPrczc0NxcWWTa4nhMDSpUvh4+MDpVKJ4OBgxMXFVXqcwWDAvn374OTkhKSkpPvuGxcXBycnJ6xevdqi2hoyDoUnIiJ7V+UAJISAXC4323bvUhjVER0djdWrV2PNmjVITk5GWFgYxowZg/T09AqPOXDgABQKBQYPHgy9Xn/f8ycnJ2Py5Mmcofp/sAWIiIjsnUXD4KdNmwYXFxcAgF6vx4wZM+Dm5gYAuHXrlsVfvm7dOsydOxdjx44FAERFRWH37t2IiYmpcJh9jx49cPr0aVy6dAlDhgyp8NxZWVkYNWoUPvnkE7z22msW19aQsQWIiIjsXZUDUHh4uNn7p59+2uy9h4cHpkyZUuUvLiwsREpKCkJCQsy29+vXDwkJCRUe5+Ligg4dOkClUlW4j06nw4gRI/DKK69g4sSJVQpARUVFKCr6OxDodLoqXEX9xBYgIiKyd1UOQBs2bLDqF2u1Wggh4O7ubrZdrVbj+PHj1T6vXq/HhAkTMHDgQItGpEVFRdnN5I6lLUDZbAEiIiI7ZfFSGNamUFRrMuoKzZ8/H40aNbK40/P8+fORn59vemVmZlq1rrrE6+5cQGwBIiIieyVZAFKr1ZDJZMjNzTXbrtVqodFoqn3es2fPYs+ePXBxcYFKpYJKpUJGRgbmzJmDZs2aVXicUqmEu7u72auh+ns5DMtG7RERETUU1m1+sYBKpUJQUBDi4+Px8MMPm7YfPnwYI0eOrPZ5o6OjUVBQYLbt4YcfRnh4OJ599tlqn7chaXr3EVjOzSIYjAJyh5qN5iMiIqpvJAtAABAREYE333wT3bp1Q0BAALZv347U1FRs27YNALBkyRIsWbIE586dg6+vLwCgpKQEN27cMHVSLigoQF5eHtzc3KBQKNCyZcsy3+Po6IimTZvC39+/9i6uDmvi6gSZDDAKIOdmsalPEBERkb2QNADNmjULOTk5iIyMRHZ2NgIDAxEbG2sKKkajEQaDAUL8vWhnfHw8Bg0aZHofGhoKANi3bx8GDhxYq/XXVwq5A5q4OEF7sxjXbxQxABERkd2RiXvTBZnodDp4eHggPz+/QfYHGvbRAaRdLcAXU3shtH31+1wRERHVJVX9/ZZ8FBhJg5MhEhGRPWMAslMcCk9ERPaMAchOsQWIiIjsGQOQneJyGEREZM8YgOzU3y1AnAyRiIjsDwOQnWILEBER2TMGIDv193IYDEBERGR/GIDsVOkjsJxbxSgxGCWuhoiIqHYxANmpJq5OcJAB4u5yGERERPaEAchOyR1kaOJ6tx8QH4MREZGdYQCyY5wMkYiI7BUDkB3jUHgiIrJXDEB2TMOh8EREZKcYgOwYl8MgIiJ7xQBkxzgZIhER2SsGIDvm1ehOJ2i2ABERkb1hALJjGjcVALYAERGR/WEAsmNsASIiInvFAGTHSkeB5d7SQ8/lMIiIyI4wANkxTxcnyB1kAAAt5wIiIiI7wgBkxxwcZFC78jEYERHZHwYgO8eh8EREZI8YgOxc6WSIXBCViIjsCQOQnWMLEBER2SMGIDvH5TCIiMgeMQDZOS+3O52g2QJERET2hAHIzrEFiIiI7BEDkJ3TsA8QERHZIQYgO+dlagHiRIhERGQ/GIDsXGkLUP5tPYpKDBJXQ0REVDsYgOych7MjFFwOg4iI7IykAUgIgaVLl8LHxwdKpRLBwcGIi4ur9DiDwYB9+/bByckJSUlJZT5fuXIlevXqBU9PT6jVaoSFheHUqVM2uIL6z8FBZpoLiB2hiYjIXkgagKKjo7F69WqsWbMGycnJCAsLw5gxY5Cenl7hMQcOHIBCocDgwYOh1+vL3efo0aN47rnn8PPPP+OHH36AXC5HWFgYCgsLbXUp9VrpSLAr+bw/RERkHyQNQOvWrcPcuXMxduxYBAQEICoqCv7+/oiJianwmB49euD06dPYu3dvhfts374d06dPR/fu3dGzZ09ER0cjKysLqampNriK+s9X7QIASL9+U+JKiIiIaodCqi8uLCxESkoKQkJCzLb369cPCQkJFR7n4uKCDh06QKVSVfm7rl27BgDw9PSsXrENnF9TNwDA2Ws3JK6EiIiodkgWgLRaLYQQcHd3N9uuVqtx/Phxq32PEAKLFi3C4MGD8cADD1S4X1FREYqK/u4Do9PprFZDXWcKQNkMQEREZB8kHwWmUNg2g82ePRvJycn3fawGAFFRUfDw8DC9fHx8bFpXXdJOcycAnbt2A0IIiashIiKyPckCkFqthkwmQ25urtl2rVYLjUZT4/MbDAZERERg9+7d2L9/f6WBZv78+cjPzze9MjMza1xDffGAlyscZICusATZHAlGRER2QLIApFKpEBQUhPj4eLPthw8fRrdu3Wp0bp1Oh9GjR+P48eM4dOgQ2rVrV+kxSqUS7u7uZi97oXKUw6fJnY7Q7AdERET2QNJHYBEREVi1ahViY2ORlpaGZcuWITU1FeHh4QCAJUuWQKFQICMjw3RMSUkJ8vLyTH10CgoKkJeXh5KSEtP7vn374vr16/jyyy9RWFiICxcu4MKFC7hy5UrtX2Q9ce9jMCIiooZOsk7QADBr1izk5OQgMjIS2dnZCAwMRGxsLPz9/QEARqMRBoPBrF9KfHw8Bg0aZHofGhoKANi3bx8GDhwIrVZrGu7eqVMns+8bMGAA9u/fb+Orqp/8mrrhlz+v4Vw2h8ITEVHDJxPs9VounU4HDw8P5Ofn28XjsG3HMjH362T08/PCv5/vLXU5RERE1VLV32/JR4FR3dCuqSsA9gEiIiL7wABEAAA/TSMAwBVdIW4UlUhcDRERkW0xABEAwMPF0bQoKjtCExFRQ8cARCbtNHwMRkRE9oEBiExKl8Q4xyUxiIiogWMAIhMuikpERPaCAYhMuCgqERHZCwYgMimdDfqi9hb0BqPE1RAREdkOAxCZNPdQwdVJjhKjQIaWM0ITEVHDxQBEJjKZDO3YD4iIiOwAAxCZMS2KyjXBiIioAWMAIjMcCUZERPaAAYjMlLYAMQAREVFDxgBEZvzuLop6LvsGjEYhcTVERES2wQBEZnzVrlA4yHCr2IArukKpyyEiIrIJBiAy4yh3gK/aBQAfgxERUcPFAERlsCM0ERE1dAxAVMbfQ+EZgIiIqGFiAKIy2AJEREQNHQMQlVEagNgCREREDRUDEJXR9u4jsOs3ipF3q1jiaoiIiKyPAYjKcFMq0NxDBYCtQERE1DAxAFG52A+IiIgaMgYgKhcXRSUiooaMAYjK1Y4tQERE1IAxAFG5/LgoKhERNWAMQFSu0j5Ambm3UKg3SFwNERGRdTEAUbm83JzgrlJACCD9OvsBERFRw8IAROWSyWQcCUZERA0WAxBViAGIiIgaKgYgqhAXRSUiooaKAYgqxBYgIiJqqCQPQEIILF26FD4+PlAqlQgODkZcXFylxxkMBuzbtw9OTk5ISkoq83l+fj7Cw8Ph6ekJV1dXDB8+HOfPn7fBFTRcpQHo/PWbMBiFxNUQERFZj+QBKDo6GqtXr8aaNWuQnJyMsLAwjBkzBunp6RUec+DAASgUCgwePBh6vb7cfZ577jmkpKRgz549OHz4MBwcHDBq1CgYDBzSXVWtPF3gpHBAcYkRl3JvS10OERGR1UgegNatW4e5c+di7NixCAgIQFRUFPz9/RETE1PhMT169MDp06exd+/ecj/Pzs7Gt99+i+joaPTt2xddunTBF198gbS0NBw4cMBGV9LwyB1kaOvlCgA4m10gcTVERETWI2kAKiwsREpKCkJCQsy29+vXDwkJCRUe5+Ligg4dOqBdu3blfp6YmAij0Yju3bubtnl5eSEgIKDC8xYVFUGn05m9iEtiEBFRwyRpANJqtRBCwN3d3Wy7Wq3GtWvXqn3e7OxsODs7w9HRscrnjYqKgoeHh+nl4+NT7e9vSEqXxDh3jZMhEhFRwyH5IzAAUCgUkp9z/vz5yM/PN70yMzOtXlN9ZGoB4lB4IiJqQKyfPCygVqshk8mQm5trtl2r1UKj0VT7vBqNBjdv3oRerzdrBbrfeZVKJZRKZbW/s6G6d1FUIQRkMpnEFREREdWcpC1AKpUKQUFBiI+PN9t++PBhdOvWrdrnDQ4OhhAChw4dMm3LyclBWlpajc5rj9pqXCGTAfm39dDeLJa6HCIiIquQ/BFYREQEVq1ahdjYWKSlpWHZsmVITU1FeHg4AGDJkiVQKBTIyMgwHVNSUoK8vDxTR+WCggLk5eWhpKQEANC0aVOMGzcOM2fOxJEjR5CcnIyIiAj4+flh4MCBtX6N9ZnKUY5Wns4A2BGaiIgaDskD0KxZs/Daa68hMjISnTt3xvbt2xEbGwt/f38AgNFohMFggBB/T8QXHx8PT09PdO3aFQAQGhoKT09Ps5akzz//HF27dsWwYcPQp08fFBQUYPfu3Tbpb9TQ3fsYjIiIqCGQiXuTBZnodDp4eHggPz+/zCg1e7N8dyr+9Vs6nn2oDd4Z3VHqcoiIiCpU1d9vyVuAqO5rxxYgIiJqYBiAqFKla4KdYwAiIqIGggGIKlUagP7KL8TNohKJqyEiIqo5BiCqVGMXJ3i5OQEAzmdzRmgiIqr/OCSKqqStxg3Xb+Tgg7g0+DZxkbocIrsgk8mgcpTDxUkOZ0c5VHf/0/nuNpWjHCpHBzhYaYJSf283uDjxZ4HsA//SqUqCmrvjaHoO9qdlS10KEdmIl5sSG6f2RMcWHlKXQmRzDEBUJTMH+6GpuxKFxQapSyGyGwYhUKg34rbegNvFd1/6e/7z7n+3hoJCPa7fKMLj637Hv8J7oE9btVXOS1RXcR6gCnAeICKyJ7pCPZ7feBxH03PgpHDAmie6YWjHZlKXRWQxzgNERERV5q5yxBdTe+GRIG8Ulxjxwr8TsO1YptRlEdkMAxAREQG4s/bf2qe649EerWAUwNyvk/Hpr+ekLovIJhiAiIjIRCF3wHsTu2D6gLYAgBX/+RPv7jkNo5G9JahhYQAiIiIzMpkM88MCsWBEBwDA/x04jzk7kqE3GCWujMh6GICIiKhcEaHtsOofXSF3kOHrE1l44csEFOo5EpQaBo4CqwBHgRER3bE39SpmbD6BohIjWjdxgaaRssbndHGSY8YgPw63J6ur6u83A1AFGICIiP52ND0Hz208hoJC660HqFQ4YP0zPfGQn5fVzknEAFRDDEBEROayC4pw4mKuVc711bFM/PLnNagc74SgB9sxBJF1MADVEAMQEZHtFJUY8MKXCdiXlg2VowM2PNMLfdvxcRjVHCdCJCKiOkupkGPt5BAMDNCgUG/E1JhjOHJeK3VZZEcYgIiISBIqRzk+nRyC0PYa3NYb8GzMMRxNz5G6LLITDEBERCQZlaMc//d0CPr7e+FWsQHPbDiKYxcYgsj2GICIiEhSKkc5/jWlB/r53Q1B648iIYMhiGyLAYiIiCRXGoIebKfGzWIDwtcfQ0KGdUacEZWHo8AqwFFgRES173axAc/GHMXv53PgplRg6kNtIHeo+f9X92rkhGEdm8HLreaTOFLdxmHwNcQAREQkjVvFJXhmg/U7RMsdZHjIzwvjgltgaMdmcFMqrHp+qhsYgGqIAYiISDo3i0qwPj4dl3WFNT6XEEDqZR3+yMwzbVM5OmBIoDfGBrfEgPYaOCnYI6ShYACqIQYgIqKGJf36TcQm/YWdSZdw/vpN03YPZ0eM6NwcQ4O84VrHWoVaejqjZWNnqcuoVxiAaogBiIioYRJC4NQlHb5LuoTv//gL1wqKpC6pQg4yYM6wDnhhQFvIZDKpy6kXGIBqiAGIiKjhMxgFjpzX4rukSzhxMQ/GOvSTWGIQuJhzCwDwSJA3Vv2jKzycHSWuqu5jAKohBiAiIpKSEAJbj2XinZ0pKDYY4at2wdqnQhDUgr9J98O1wIiIiOoxmUyGJ3q1xo7IvmjZ2BkZ2lsY/8lB7EjIkrq0BoEBiIiIqA7r0qoxds/uh4EBGhSVGPH69j8w/5uTKNQbpC6tXpM8AAkhsHTpUvj4+ECpVCI4OBhxcXH3Paa4uBgvv/wymjZtCmdnZzz00ENISEgw2+fixYt4/PHH0bRpU3h4eGDAgAH47bffbHkpRERENtHYxQnrw3vilSHtIZMBW45exD8+PYzMu32EyHKSj/eLjo7G6tWrsX79enTo0AExMTEYM2YMUlNT8cADD5R7zIIFCxAbG4utW7eiWbNmWLlyJYYPH45z587B3d0dRqMRjzzyCDp06IDdu3dDqVRi7dq1GD58OFJTU+Hr61vLV0lERFQzDg4yvDTEH8GtG+OlrYk4eSkfo/4Zj/cmdkGnlvWzX1BjFyfJJqSUvBN0586dMXnyZMybN8+0rUuXLhg/fjwWL15cZv+SkhJoNBp8+umneOyxxwAABoMBGo0GH3zwAZ599lmcPXsW/v7+SE5ORufOnQEARqMRKpUKW7ZswcSJEyuti52giYiorsrKvYUZm07gj6x8qUupkXfHd8aTvVtb9Zz1ohN0YWEhUlJSEBISYra9X79+ZR5plTp//jzy8vLMjpHL5ejTp4/pmFatWqFly5Z45513cPnyZQDAsWPHoFQq8eCDD9roaoiIiGpHK08XbHuhL555sA3clAooFQ718iWXMIVI+ghMq9VCCFEmoanVahw/frzcY7KzswGg3GOuXbsGAFCpVNi3bx/69esHX19fDB48GGlpafjxxx/RvHnzcs9bVFSEoqK/J8PS6XTVvi4iIiJbUyrkWDSmIxaN6Sh1KfWS5J2gAUChsDyH3e+Y4uJizJgxAxEREcjIyMCYMWOgUCjwyiuvIDc3t9xjoqKi4OHhYXr5+PhYXBMRERHVD5IGILVaDZlMViaUaLVaaDSaco8p3X6/Y7766iukpKRgyZIlaN68OV588UUkJSUhKysLH3/8cbnnnT9/PvLz802vzMzMml4eERER1VGSBiCVSoWgoCDEx8ebbT98+DC6detW7jFt27aFh4eH2TFGoxFHjhwxHZOfn4+ioiKUlJSY9nF1dUXjxo2Rk5NT7nmVSiXc3d3NXkRERNQwSf4ILCIiAqtWrUJsbCzS0tKwbNkypKamIjw8HACwZMkSKBQKZGRkALjz6Gvq1KlYsGAB9u/fj9OnT2PWrFkAYBrdNWrUKOj1ejz66KM4cuQIkpOTMWfOHPz3v//F5MmTpblQIiIiqjMknwdo1qxZyMnJQWRkJLKzsxEYGIjY2Fj4+/sDuNO6YzAYcO9o/aioKOj1evzjH/9AQUEBunfvjh9++AGenp4AgDZt2uCXX37BW2+9hZEjR6KkpAQ9evTAvn370KtXL0muk4iIiOoOyecBqqs4DxAREVH9Uy/mASIiIiKSAgMQERER2R0GICIiIrI7DEBERERkdxiAiIiIyO4wABEREZHdYQAiIiIiu8MARERERHZH8pmg66rS+SF1Op3ElRAREVFVlf5uVzbPMwNQBQoKCgAAPj4+EldCREREliooKICHh0eFn3MpjAoYjUb89ddfaNSoEWQymdXOq9Pp4OPjg8zMTC6xUQt4v2sX73ft4z2vXbzftas691sIgYKCArRo0QIODhX39GELUAUcHBzQqlUrm53f3d2d/+OpRbzftYv3u/bxntcu3u/aZen9vl/LTyl2giYiIiK7wwBEREREdocBqJYplUq88847UCqVUpdiF3i/axfvd+3jPa9dvN+1y5b3m52giYiIyO6wBYiIiIjsDgMQERER2R0GICIiIrI7DEBERERkdxiAapEQAkuXLoWPjw+USiWCg4MRFxcndVkNitFoRGJiIry8vPDdd9+ZfVZcXIyXX34ZTZs2hbOzMx566CEkJCRIU2g9t2vXLjz88MNo3rw53Nzc0KNHD8TGxpo+5722rvj4eDzyyCNo0aIFnJ2dERgYiPfee8+01hHvt+1kZWWhVatWGDdunGkb77d17d+/HzKZrMwrODgYgO3uNwNQLYqOjsbq1auxZs0aJCcnIywsDGPGjEF6errUpTUIGRkZUCgU6N69O7RabZnPFyxYgNjYWGzduhUJCQlo3749hg8fzgVvq+HYsWMIDQ3F9u3bceTIEQwdOhQTJ05EUlISAN5razt58iR69eqFLVu2ICkpCa+//jrefvttbNy4EQDvt63odDqMGDECRUVFZtt5v23jjz/+QHp6uum1Z88eADa834JqTadOncSKFSvMtnXu3FksXLhQoooaluLiYnH69Glx+vRpAUB8++23ps/0er1o3Lix2Lp1q2lbSUmJ8PT0FOvXr5eg2oanTZs24oMPPuC9riXBwcHi7bff5v22keLiYjFkyBAxe/ZsER4eLsaOHSuE4L9LbGHfvn0CgMjNzS3zmS3vN1uAaklhYSFSUlIQEhJitr1fv35sOrUSR0dHdOjQAR06dCjz2fnz55GXl2d2/+VyOfr06cP7bwV6vR55eXnw9PTkvbax4uJibN26FVlZWXjqqad4v23k+eefh5ubGz766COz7bzfttOqVSs0adIEPXr0wPr16wHY9n5zMdRaotVqIYQos5ibWq3G8ePHJarKfmRnZwNAuff/2rVrUpTUoHz44YdwcHDA+PHjkZKSAoD32hZWrFiBBQsWoEmTJti0aRMCAgJw8OBBALzf1rRw4UKkpaVh3759ZVYT579LrK99+/bYtWsXWrZsiVu3bmH37t2YNm0a5HI5/Pz8ANjmfjMA1TKFgrdcSrz/1vfVV19h8eLF+Pbbb9G4cWPTdt5r65s2bRqGDx+Oo0eP4tFHH8Xnn3+O5s2bA+D9tpbvv/8emzdvxuHDh+Hs7Fzhfrzf1tOiRQu0aNHC9P7BBx/EpUuXsHbtWnzwwQcAbHO/+QislqjVashkMuTm5ppt12q10Gg0ElVlP0rvMe+/da1fvx7Tpk3D9u3bMWzYMAC817akVqsRHByMiIgIPP300/i///s/3m8rO3fuHDIyMuDj4wOVSgWVSoUvv/wS33//PVQqFRo1agSA99vW/P39kZuba9O/bwagWqJSqRAUFIT4+Hiz7YcPH0a3bt0kqsp+tG3bFh4eHmb332g04siRI7z/1SCEwNtvv425c+fihx9+wMiRI02f8V7XDp1OB6PRyPttZeHh4Th58iSSkpJMrzFjxmDQoEFISkpCQEAA77eVGQyGMtuOHTuGTp062fTvm214tSgiIgJvvvkmunXrhoCAAGzfvh2pqanYtm2b1KU1CEaj0WxY5M2bN5GXlwcXFxc4OTlh6tSpWLBgAXx9feHt7Y01a9YAACZOnChVyfXWlClT8NNPP2Hr1q1o0aIFLly4YPqsTZs2vNdW9uijj6Jv377o3bs33N3dsWvXLmzevBlff/01FAoF77cVeXp6wtPT02ybh4cHhBCmARa839b1zDPPwMfHB+PHj4dKpcLmzZuxa9cuHDt2zLZ/3zUaQ0YWMRqN4p133hEtWrQQjo6OokuXLuKHH36QuqwGIz09XQAo89qwYYMQQojCwkIxc+ZM4eXlJZRKpejbt684evSotEXXU76+vuXe69J/pfBeW9dHH30kevbsKZo0aSLc3NxEr169xM6dO02f837b1r3D4IXg/ba2DRs2iB49eojGjRsLpVIp+vTpIw4ePGj63Fb3WybE3alEiYiIiOwE+wARERGR3WEAIiIiIrvDAERERER2hwGIiIiI7A4DEBEREdkdBiAiIiKyOwxAREREZHcYgIiIiMjuMAARUa3Q6XRYtGgRunTpAldXV7i7u6NTp06IjIzEiRMnpC6vUi+//DIGDhxotm3//v2QyWRmS4EQUf3AtcCIyOauXLmC0NBQuLi44K233kLXrl1x48YNnDp1CuvXr8err76K/fv3S12mxfr06YP09HS0atVK6lKIyEJsASIim3vxxRdhMBhw6NAhPProowgICEBISAjCw8Px66+/4sMPPwQA5OfnY9q0adBoNPDw8MCwYcOQlpZmOs/AgQMxefJkzJw5E15eXlCr1XjppZdgNBpN+xQVFWHOnDlo0aIF3Nzc8NBDD+HIkSOmz5955hmEhYVhyZIl8Pf3h5OTE9LT07F06VI88MADcHZ2hlqtxsiRI3HmzBkAQExMDD7++GP8+uuvkMlkkMlkiImJwZ9//okHHngAWVlZpvNv2LAB7du3h5OTEwICAvDFF1+Y3QuZTIZly5ZhzJgxcHV1RevWrbFx40ab3HciqhgDEBHZVG5uLnbu3Im33noLLi4u5e7TvXt3CCEwcuRIXLp0Cd9//z0OHDgAtVqNUaNGobi42LTvjh070KRJE+zfvx9r167F2rVrsX37dtPnU6dOxYEDB7Bp0yb8/vvv6N27N8LCwqDVak37/Pjjj7hy5Qp27NiBEydOoFmzZmjatClWrlyJhIQE7N69G/n5+fjHP/4BAJg0aRKeffZZ9O7dG+np6UhPT8ekSZPKXMeuXbsQGRmJl156CQkJCZg9ezYiIiLwww8/mO33ySefYNy4cTh27BimTJmCadOmITMzs0b3mYgsVOPlVImI7uPIkSMCgDhx4sR99/v555+Fu7u7uHXrlmnbjRs3hIODgzhw4IAQQogBAwaI2bNnmx3Xt29f07azZ88KBwcHkZmZafrcaDSKZs2aiS+++EIIcWdl71GjRlVa93fffScAiBs3bgghhHjppZfEgAEDzPZJTEwUAER6eroQQoh+/fqJF1980WyfF198UfTv39/0HoD45ptvTO9v374tHBwczLYRke2xBYiIaoVSqTT997i4OKhUKtNr6NChSEhIQEFBATw9PU3b1Wo1jEYjMjIyTMfKZDKz82o0GuTn5wMATpw4AaPRCD8/P9M5nJ2dcfXqVbNzyOXyMvWdPn0aERER6NKlC1q0aIEnn3wSAKDX66t8jSkpKejdu7fZtl69eiElJcVs273XoFKp4ObmZroGIqod7ARNRDbl5+cHBwcHHD16FEFBQQCABx98EElJSQCAqKgoZGRkwGg0olmzZvjll1/KnKN58+YVnv/eMGE0GiGXy3H8+HEoFOb/evPy8qrwHJcuXULv3r3Rs2dPvPXWW2jdujVOnz6NqVOnWnKp1fa/oY6IbI8BiIhsqkmTJhgxYgRWrFiBCRMmwN3dHa6urujQoQMAwNPTExkZGejSpQuuXLkCR0dHtGvXrlrf1blzZxgMBmRnZ2PQoEFVPu7gwYMoLCzEjz/+aApOhYWFZvuoVCoUFRXd9zxBQUH4/fffMWXKFNO2I0eOoGPHjhZcBRHVBj4CIyKb++STT3D79m307t0bX3zxBU6fPo3U1FTs2rULR48eBQAMHToUISEhGDlyJL777jukpaVh7969mDJlCpKTk6v0PUFBQRg/fjwmT56MTZs24c8//8Rvv/2G2bNnY8+ePRUeFxAQgJKSEmzYsAGnT5/G119/jblz55rt07FjR5w4cQKxsbE4duwYzp07V+Y88+bNw/r16xEdHY1Tp04hOjoa69evxxtvvGHB3SKi2sAWICKyOR8fHyQmJmLVqlV4//33TeHBz88PI0eORHh4OORyOX766ScsWLAAL774IrRaLVq0aIEhQ4ZYNM/O5s2bsXjxYrz55pv466+/0LRpU4SGhiIwMLDCY7p27YqPP/4YCxcuREFBAXr16mUapVXq8ccfx08//YTJkydDLpdj48aNaN26tdl5Ro8ejTVr1uC9997DK6+8gjZt2mDt2rUYMWKEhXeMiGxNJoQQUhdBREREVJv4CIyIiIjsDgMQERER2R0GICIiIrI7DEBERERkdxiAiIiIyO4wABEREZHdYQAiIiIiu8MARERERHaHAYiIiIjsDgMQERER2R0GICIiIrI7DEBERERkd/4/MXOE1CcNS/YAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best angles (rad): [ 0.         -0.26179939 -0.26179939 -0.26179939 -0.13089969  0.26179939\n",
      "  0.26179939  0.26179939  0.         -0.26179939 -0.26179939 -0.26179939\n",
      " -0.13089969  0.26179939  0.26179939  0.26179939]\n",
      "Best fitness: 8.083380e-02\n"
     ]
    }
   ],
   "source": [
    "# GA parameters\n",
    "pop_size   = 100\n",
    "num_gens   = 50\n",
    "mut_rate   = 0.3\n",
    "\n",
    "# Penalty multiplier for By and Bz\n",
    "alpha      = 1e8\n",
    "\n",
    "# Number of restarts\n",
    "n_restarts = 1\n",
    "\n",
    "\n",
    "# Preallocate storage\n",
    "d = positions_cpu.shape[0]   # number of magnets\n",
    "best_angles     = np.zeros((n_restarts, d))\n",
    "best_fitness    = np.zeros(n_restarts)\n",
    "Bx_best         = np.zeros(n_restarts)\n",
    "By_best         = np.zeros(n_restarts)\n",
    "Bz_best         = np.zeros(n_restarts)\n",
    "fitness_history = np.zeros((n_restarts, num_gens))\n",
    "\n",
    "\n",
    "for run in range(n_restarts):\n",
    "\n",
    "    (best_angles[run],\n",
    "        best_fitness[run],\n",
    "        _,\n",
    "        _,\n",
    "        _,\n",
    "        fitness_history[run]\n",
    "    ) = genetic_algorithm_gpu(\n",
    "        positions_cpu, \n",
    "        init_orientations_gpu,\n",
    "        l_values, m_values,\n",
    "        sh_coeff_X, sh_coeff_Y, sh_coeff_Z, points_of_interest,\n",
    "        possible_angles_cpu,\n",
    "        population_size=pop_size,\n",
    "        generations=num_gens,\n",
    "        mutation_rate=mut_rate,\n",
    "        objective='homogeneity', # homogeneity or maxBx\n",
    "        parent_selection='ACROMUSE_adaptive',\n",
    "        alpha=alpha\n",
    "    )\n",
    "    print(f\"Run {run+1}/{n_restarts} done, best fitness = {best_fitness[run]:.3e}\")\n",
    "    print(f\"Best angles (rad): {best_angles[run]}\")\n",
    "\n",
    "best_run = np.argmin(best_fitness) # min fitness wins\n",
    "# best_run = np.argmax(best_fitness) # max fitness wins\n",
    "\n",
    "\n",
    "# Plot the results\n",
    "plt.plot(fitness_history[best_run])\n",
    "plt.xlabel(\"Generation\")\n",
    "plt.ylabel(\"Best Fitness\")\n",
    "plt.title(f\"Best Run (#{best_run+1}) Convergence\")\n",
    "plt.show()\n",
    "\n",
    "# Print the best angles and best fitness values\n",
    "print(\"Best angles (rad):\", best_angles[best_run])\n",
    "print(f\"Best fitness: {best_fitness[best_run]:.6e}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "magoptlib",
   "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.10.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}