{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "887123a4",
   "metadata": {},
   "source": [
    "# Latent Ewald Summation (LES), block by block: reproducing the original implementation with `xnn`\n",
    "\n",
    "[LES (Cheng, *npj Comput Mater* **11**, 80, 2025)](https://doi.org/10.1038/s41524-025-01577-7)\n",
    "adds **long-range interactions to any short-range MLIP**: a small MLP maps each\n",
    "atom's invariant features to a hidden variable $q$ (eq 2; analogous to\n",
    "environment-dependent partial charges, but unconstrained), and an Ewald\n",
    "summation over the structure factor of $q$ supplies the long-range energy\n",
    "\n",
    "$$E_{lr} = \\frac{1}{V}\\sum_{0<k<k_c} \\frac{e^{-\\sigma^2k^2/2}}{k^2}\\,\\lvert S(k)\\rvert^2,\n",
    "\\qquad S(\\mathbf k) = \\sum_i q_i e^{i\\mathbf k\\cdot \\mathbf r_i}. \\tag{eqs 3–4}$$\n",
    "\n",
    "In xnn this is `xnn.common.models.LatentEwald`, a **model-agnostic wrapper**\n",
    "(like `ForceStressOutput`) that consumes the `\"node_features\"` every xnn model\n",
    "now exposes, so *all* GNNs and short-range models (MACE, NequIP, Allegro, CACE,\n",
    "SchNet, PhysNet, HDNNP, ANI) gain long-range physics the same way. This\n",
    "notebook checks every block against the **original**\n",
    "[`cace`](https://github.com/BingqingCheng/cace) `EwaldPotential` (the paper's\n",
    "reference implementation) and ends with a whole-model CACE-LR weight-transplant\n",
    "parity check, the same protocol as the other model series."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "460b1113",
   "metadata": {},
   "source": [
    "## 0. Setup: `float64` for exact comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f0e84aa9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-22T03:15:15.544559Z",
     "iopub.status.busy": "2026-09-22T03:15:15.544442Z",
     "iopub.status.idle": "2026-09-22T03:15:17.842388Z",
     "shell.execute_reply": "2026-09-22T03:15:17.841767Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xnn: 0.2.0 | cace (original): 0.1.0\n"
     ]
    }
   ],
   "source": [
    "# silence the expected warnings\n",
    "import logging, warnings\n",
    "logging.disable(logging.WARNING)\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import numpy as np\n",
    "import torch\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "torch.set_default_dtype(torch.float64)\n",
    "torch.manual_seed(0)\n",
    "\n",
    "import xnn, cace\n",
    "from xnn.common.models import EwaldSummation, LatentEwald\n",
    "print(\"xnn:\", xnn.__version__, \"| cace (original): 0.1.0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "261d4d0a",
   "metadata": {},
   "source": [
    "## Block 1: The Ewald summation of a hidden variable  ·  eqs 3–5\n",
    "\n",
    "`EwaldSummation` ports the original `cace.modules.EwaldPotential` math:\n",
    "reciprocal-space sum on a k-grid limited by $|k| \\le 2\\pi/\\mathrm{dl}$ with\n",
    "hemisphere symmetry factors and triclinic support, tinfoil boundary conditions\n",
    "(the $k=0$ term is omitted, so a net $q$ is fine), an optional Gaussian\n",
    "self-interaction removal, the $1/r^6$ dispersion kernel (paper eq 5), and the\n",
    "`erf`-converged real-space direct sum used for non-periodic structures.\n",
    "\n",
    "We compare against upstream in float64 through its dtype-safe orthorhombic\n",
    "reference loop (`compute_potential`); the production triclinic path is\n",
    "float32-only upstream (its k-grid is hard-cast to float32), so it is compared\n",
    "in float32 below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "78ef317a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-22T03:15:17.844021Z",
     "iopub.status.busy": "2026-09-22T03:15:17.843917Z",
     "iopub.status.idle": "2026-09-22T03:15:17.915208Z",
     "shell.execute_reply": "2026-09-22T03:15:17.914613Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "electrostatic kernel (p=1):  |xnn - original| = 2.220446049250313e-16\n",
      "dispersion kernel   (p=6):  |xnn - original| = 8.881784197001252e-16\n",
      "molecular direct sum     :  |xnn - original| = 0.0\n",
      "triclinic forward (float32): rel diff = 9.21e-08\n"
     ]
    }
   ],
   "source": [
    "from cace.modules.ewald import EwaldPotential\n",
    "\n",
    "rng = np.random.default_rng(1)\n",
    "pos = torch.tensor(rng.uniform(0, 6, (12, 3)))\n",
    "q = torch.tensor(rng.normal(size=(12, 4)))          # 4-dimensional hidden variable\n",
    "box = torch.tensor([6.0, 6.0, 6.0], dtype=torch.float64)\n",
    "\n",
    "up = EwaldPotential(dl=1.9, sigma=1.0, remove_self_interaction=False)\n",
    "mine = EwaldSummation(dl=1.9, sigma=1.0, remove_self_interaction=False)\n",
    "\n",
    "pot_up, _ = up.compute_potential(pos, q, box)\n",
    "print(\"electrostatic kernel (p=1):  |xnn - original| =\",\n",
    "      abs(float(mine.reciprocal(pos, q, torch.diag(box))) - float(pot_up.sum())))\n",
    "\n",
    "pot6, _ = EwaldPotential(dl=1.9, exponent=6, remove_self_interaction=False).compute_potential(pos, q, box)\n",
    "mine6 = EwaldSummation(dl=1.9, exponent=6, remove_self_interaction=False)\n",
    "print(\"dispersion kernel   (p=6):  |xnn - original| =\",\n",
    "      abs(float(mine6.reciprocal(pos, q, torch.diag(box))) - float(pot6.sum())))\n",
    "\n",
    "pot_rs, _ = up.compute_potential_realspace(pos, q)\n",
    "print(\"molecular direct sum     :  |xnn - original| =\",\n",
    "      abs(float(mine.realspace(pos, q)) - float(pot_rs.sum())))\n",
    "\n",
    "# float32: the upstream production (triclinic) forward path\n",
    "torch.set_default_dtype(torch.float32)\n",
    "cell32 = torch.tensor([[6.0, 0, 0], [1.2, 5.8, 0], [-0.7, 0.9, 6.3]], dtype=torch.float32)\n",
    "data = {\"positions\": pos.float(), \"cell\": cell32.unsqueeze(0), \"q\": q.float(),\n",
    "        \"batch\": torch.zeros(12, dtype=torch.long)}\n",
    "e_up = float(EwaldPotential(dl=1.9, sigma=1.0, remove_self_interaction=False,\n",
    "                            feature_key=\"q\")(data)[\"ewald_potential\"][0])\n",
    "e_x = float(EwaldSummation(dl=1.9, sigma=1.0, remove_self_interaction=False)\n",
    "            .reciprocal(pos.float(), q.float(), cell32))\n",
    "print(f\"triclinic forward (float32): rel diff = {abs(e_x - e_up)/abs(e_up):.2e}\")\n",
    "torch.set_default_dtype(torch.float64)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b13d65aa",
   "metadata": {},
   "source": [
    "### The reciprocal-space kernels\n",
    "\n",
    "$e^{-\\sigma^2 k^2/2}/k^2$ captures $1/r$ electrostatics; the $p=6$ kernel\n",
    "(eq 5) captures London dispersion. The smearing $\\sigma$ (default 1 Å) trades\n",
    "real-space error against the required k-cutoff (paper fig 7)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e7af8f06",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-22T03:15:17.916547Z",
     "iopub.status.busy": "2026-09-22T03:15:17.916471Z",
     "iopub.status.idle": "2026-09-22T03:15:18.336732Z",
     "shell.execute_reply": "2026-09-22T03:15:18.336068Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAE2CAYAAAB80vHuAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAp7tJREFUeJzs3Xd4k1X7wPFvVtPd0k1L2XvvvbcyBEUBRUHhh+BE5VVx4kRfVPBVUUG2A1AEBQEBFYrsPcqmBTroprtN2ub5/ZE2NLSFFtqm4/5c13MlOc+6nxRycuec5xyVoigKQgghhBBCCCGEKHVqWwcghBBCCCGEEEJUVZJ0CyGEEEIIIYQQZUSSbiGEEEIIIYQQooxI0i2EEEIIIYQQQpQRSbqFEEIIIYQQQogyIkm3EEIIIYQQQghRRiTpFkIIIYQQQgghyogk3UIIIYQQQgghRBmRpFsIIYQQQgghhCgjknSLSsnPz48pU6bYOoxKa/v27ahUKnbs2FEu5/vll19QqVQcOnSoXM5Xms6ePYtKpeL7778v1vaXLl1Cr9eza9euMo7szp05cwadTsfBgwdtHYoQohqyt7dnxowZltdxcXGoVCrmz59vs5hKiy2v5eLFi6hUKr777rtyP3dpsLe3Z+bMmcXaNiMjg8DAQL744osyjurOpaam4uPjU2n/HqJ0SdItRCHc3d2ZNm2azc6/ceNGVCoV//77r02PIUruP//5D3369KFXr17lfu7Q0FA+/vhj2rdvj0qlom/fvoVu16xZM8aNG2f1pVcIIYSoLObOnQvA1KlTbXL+jIwM3nvvPZo3b46DgwNNmjThww8/JCMjw7KNs7Mzr7zyCm+88QYpKSk2iVNUHJJ0C1ENDRw4EEVRikzKxJ05e/Ys69at49lnn7XJ+SdOnMj169f57rvv8PT0vOW2Tz/9NHv27GHnzp3lFJ0QQhTOy8sLRVGqxA+BVelaKqqMjAzmz5/P1KlT0ev15X7+tLQ0+vXrx2+//cbChQtJTEzk77//JjMzk82bN1ttO3nyZBITE1myZEm5xykqFq2tAxBCiKrim2++wcvLi3vuuccm5w8KCir2tl27dqVhw4Z888039OnTpwyjEkIIIUrP6tWruX79Oo8++qhNzv/6669z6dIlzp07h4eHBwABAQG8++67BbZ1d3dn+PDhfPPNNzz//PPlHaqoQKSlW1Q4y5Yto1OnTjg6OuLi4sLgwYOLfS9wcfddsWIFXbt2xdnZGS8vL0aNGsWpU6dITExEpVKRlJTEt99+i0qlKtBNN+9etKCgILp27Yq9vT0fffQRYO7e+/DDD+Pt7Y1er6dJkya8//77ZGdnW/bPyMjglVdeoUGDBjg4OFCvXj0mT55MWFgYAF9++SUjRowAoFevXpYY8u4JOnXqlKVMpVLh4OBAq1at+Oyzz1AUpVjHKOqe7uvXr/PSSy/RsGFD7O3tadiwIa+//jqpqanFir0kEhIS6Nu3L76+vuzbt89SvmrVKrp164aTkxNOTk7079+f3bt3W+3r5eXFlClT2L9/Pz179sTBwYE33niDffv2oVKpWL9+PUuXLrVcR8eOHQscI+963nzzTRo3boxer8fHx4dJkyYRFRVV4usB2LBhA7169UKrtf49M39cX331FXXr1sXe3p6uXbvatKW5X79+bNq0iZycHJvFIISouqKjo5kwYQLu7u64u7vz2GOPkZycXGC7ou6D/vzzz2nVqhXOzs74+/szatQojh49almf/x7mlStX0qhRI+zt7WnXrh1//PFHgfNkZWUxZ84cmjdvjr29PR4eHowdO5bLly8XesyffvqJFi1aoNPp2LhxY7FiKupaEhMTef7556lduzZ2dnYEBgby7LPPkpiYWOi5165dS7NmzdDr9bRu3Zo///yzBO+89TVPnjwZe3t7fvrpJ0v5rl27GDp0KO7u7tjb29O+fXtWr15tte/AgQPp2LEjV69eZeTIkbi6ujJ8+HAyMzNRqVTMnj2bv//+m/bt22Nvb0/jxo358ccfC8SgKApff/017dq1w8HBATc3N4YPH86pU6fu6Jo2bNhAvXr1qFu3rlV5/rj++OMPWrVqhb29Pc2aNWPFihV3dK6bpaWlsWjRIh555BFLwn07/fv35+zZs1y4cKFUYhCVlCJEBTJz5kxFr9crn3/+uRIdHa1EREQoTz31lGJvb68cOXLEsp2vr68yefLkO9p3xowZik6nUz788EPl6tWrSkJCgvLbb78pU6ZMsWzj5uamPPnkk4XGqNfrlXvuuUcZNWqUcvbsWeXKlSvKhg0blLCwMMXb21tp166dcvjwYSUxMVFZuXKl4ujoqDz44IOW/adNm6b4+voqQUFBSnp6uhIWFqYsXbpUefnlly3bbNiwQQGUXbt23fY9i4uLU7777jvF3t5emTdvXrGOsW3bNgVQ/vnnH0tZQkKC0rhxY6V+/frKpk2blMTEROXSpUvKnDlzlOXLlxc79sL8/PPPCqAcPHhQURRFOX/+vNKoUSOlRYsWSmhoqGW72bNnK1qtVvn444+VyMhIJSoqSnnppZcUOzs7Zffu3ZbtPD09lYEDByrDhw9XgoODlfDwcGXdunXK3r17FUAZO3as8sYbbyjXrl1Trly5ovTp00fx9PRUUlJSLMcwGAxKjx49lICAAOX3339XkpKSlODgYKVnz55Ko0aNlOTkZEVRFOXMmTMKoKxcufKW1xgZGakAygcffFBgXV5co0aNUmbNmqVERUUpoaGhyoMPPqjY2dkp+/bts9per9crwG2XBx54oMh4PD09lT59+twy5kWLFimAcujQoVtuJ4QQJZWRkaG0aNFCqVevnhIUFKQkJSUp69evV8aPH6/o9Xrl+eeft2wbGxurAFZ12DfffKPodDpl1apVSnJyshIXF6ds2LBBefjhhy3bXLhwQQGU++67T3nqqaeU8PBwJSwsTPm///s/Ra1WK5s3b7Zsm5OTowwbNkzx8vJSVq1apVy/fl25cOGCcs899yg1a9ZUoqOjrY45atQoZfr06cqVK1eUo0ePKnv27ClWTIVdS2ZmptKuXTulZs2ayp9//qkkJSUp27dvVwICApTWrVsr6enpVud+8MEHleeee04JCwtTIiMjlZEjRyoODg5KVFTULd/zvP0XLVqkKIq5Xu/Xr5/i7e1tVYeuXbtW0Wg0yrRp05SQkBDl+vXryoIFCxSdTmfZV1EUZcCAAUqLFi2Ue++9V9m1a5cSFxenfP/990pGRoblPXr00UeVS5cuKbGxscr//d//KSqVSgkODraKa8qUKYqTk5OycOFCJS4uTrl69ary6KOPKi4uLsq5c+cs2+n1euWll1665TUqiqL4+PhYfa/KkxfX8OHDlXHjximXLl1SoqOjlddff10BlCVLllht36VLl2LVtQ0aNLDs888//1jq+rFjxyo1atRQ7O3tlfbt2ys//fRTofEePHhQAZTvvvvuttcmqi5JukWFcerUKUWlUilvv/22VbnJZFJat26tDBs2zFJ2c9Jd3H2PHTumAMprr712y1hul3S7urpakrI8Tz31lKLVapVLly5Zlc+ZM8cq+W3atKkyfvz4W56/JEl3nilTpihNmzYt1jEKS7pnzpypqNVq5eTJk0WeozixFyZ/0r1jxw7Fw8NDGTx4sJKUlGTZ5tKlS4pGoym0wu3WrZtVAunp6ak4ODgosbGxVtvlJbf33nuvVfnhw4cVwPLjgaIoyldffVXgPVAURYmIiFD0er0yd+5cRVGKn3Tv27ev0Eo9f1wDBgywKjcYDEpAQIDSv39/q/LySro3btyoAMratWtvuZ0QQpTUt99+W+hn7DfffKMAt026x40bpzRr1uyW58hLMtu2bWtVbjKZlObNmyutW7e2lK1evVoBlJ9//tlq28TERMXd3V155ZVXrI7ZqlWrAucrTkyFXcvChQsVQPn999+ttt28ebMCKF988YXVuTt37my13ZUrVxRA+fTTT2957vxJ94ULF5QmTZoozZo1U0JCQizbGAwGxdfXt0B9pCjmH9a9vLyUrKwsRVHMSTdglbAryo3ktk6dOorRaLSUp6SkKI6OjsoLL7xgKdu9e3eB90NRFMVoNCr169dXHnnkEUtZcZLuzMxMBVCee+65Auvy4vL391cMBoPVuoEDByq+vr6Wa1OUO0u6V65cqQCKRqNRHn/8cSU8PFy5du2a8vTTTyuA8u233xaIKzw8XAGUN99885bXJqo26V4uKoyNGzeiKAoPPvigVblKpaJ///637IZb3H03bdoEwMMPP3xXsfbr1w8XFxersr/++ou2bdtSv359q/IxY8ZY1gO0adOG9evX89FHH3H+/Pk7Ov+SJUvo2rUrrq6uVl3HL168eEfHA/N707p1a1q2bFnkNncb+4oVKxg0aBAPPfQQf/zxB66urpZ1mzdvJicnp8DfEGDAgAHs3r3bqpt+9+7d8fLyKvQ8w4YNs3qdd00hISGWsg0bNuDl5VVgMDl/f3+aNWtW4m7feV0Eb/53kd/IkSOtXtvZ2XHPPfewa9cusrKyLOWZmZko5h9Fb7n88ssvJYrxZnnvf/7ujUIIURr++usvXFxcCnzGjho1qlj7t2nThjNnzjBjxgwOHz58y9tg8m6nyqNSqRg5ciQnTpwgLi4OMH/m29nZFfgcdnNzo2PHjgU+82/erqQx5ffXX39hZ2dXoG4aOnQozs7Olu8HeW7ernbt2ri4uFjVYbeye/duunbtSmBgIHv37qVevXqWdQcPHiQ6OrrQunbgwIHExcURHBxsKfP09KR79+6Fnmfw4MHodDrLa2dnZ+rUqVOgrgUKnE+n09G7d+8yqWuHDBmCnZ2dVdmoUaOIjo62urZ9+/YVq67N/93KZDIB5r/JwoULCQgIwM/Pjy+++IKWLVvy5ptvWrbJI3WtALmnW1QgeffRtmnTBq1Wi0ajQa1Wo1armT9/PqmpqRiNxrvaNyYmBjAPeHE3Cts/Pj4ePz+/AuV5ZXkV/4IFC5g0aRKffvopTZo0wd/fn6lTpxb7vugvv/ySyZMnM3r0aIKDg8nKyrKMlJo/KS2pmJiY274vdxv7Tz/9hE6nY9q0aQXue877G/bo0aPA3zDvvvj89wLeKtaaNWtavbazs0Ov11tVeFFRUcTFxaHVaq3Op1KpOHbsGPHx8cW6pjzu7u4Ahd6vmMfX17fQsqysrFvuV1byzlmjRo1yP7cQomqLj48v9DPPx8cHtfr2Xz9feuklZs+ezW+//UbHjh3x9PRkzJgxHDlypMC2RX225sUB5s98o9GIo6Njgc/87du3F/jML6yOKUlM+cXHxxd53b6+vpbvB3lursPAnLgVN2nbuHEj8fHxTJ06FTc3N6t1eXXtU089ZXkf8t6LvEaC/O9FSerawuLMO1+dOnUK1O3Lli0r17oWKPBel1Tej/3du3e3+h6jUqno3bs3MTExXLlyxWofqWsFSNItKpC8D7JLly6RnZ1NTk4OJpMJk8lk+bXx5l8uS7qvt7c3ABEREXcVa/5fdvN4eHgQHR1doDyvLC9GDw8PFixYYPnF9T//+Q+//vorffv2LfDraGFWrFhBt27deOWVVwgMDLR86IeGht7NJeHt7X3b9+VuY//xxx9p06YNffv2LXRwNIDjx48X+TfMP2hJYX+DPCqV6raxeHl5UbduXbKzs63Ol3euPXv23PYY+dWpUwfgloOwFfXvQ6fTWbX629vbWw2WV9SS9wXpTl27ds0qdiGEKC2enp6WH7rzi4mJKVZ9odPpePvttwkNDSU0NJR58+Zx6tQpevfuTXh4uNW2t6p786ZP9PLywtXVFaPRWOhn/s29twqrY0oSU34eHh7ExMRYBju9Oc6be20Vpw67lffee48JEyYwfvz4AgOI5Z1rxYoVlvfh5veif//+lu1Lo65VqVTExcUVWrenp6eX6Nr0ej2+vr53VNcCVtNpdu3atVh1bcOGDS37tGnTpsjz5v19b35fpK4VIEm3qECGDx8OUGD0zNLcN6/LVmGja+bn5OSEwWAoUQwDBgzg6NGjBX7hXLt2rWV9fmq1mubNm/PCCy8wY8YMQkJCLJW2k5MTQJEx3DwvZXR0NFu3bi1wDbc6xs2GDx/OiRMnijWa6K1iv5UaNWqwbds2OnfuzODBg61GY7333ntRq9V39Pe/EyNGjODy5cscOHCgVI7n5+dHgwYNbjnSfl43uzxGo5EtW7bQq1evW36xKSsHDhzA1dX1ll8ihBDiTvTv35/k5OQC3Yd/++23Eh+rbt26PP7443zyySekpaUVaFnOG1k8j6IobNiwgdatW1uSzBEjRpCcnMyWLVtKfP47iSm/AQMGYDQaLbe45dm6dSupqakFvh/cLa1Wy4oVK3jqqaeYNGkS//vf/yzrunTpgre3N2vWrCnVcxZlxIgRKIpSqufr2bPnLevarVu3FugZ+dtvv+Hj43PLW+iKIyAggO7du7Nnzx6r3oWKovDvv/9Ss2bNAsl13veMXr163dW5ReUmSbeoMNq0acN//vMf3nrrLebOncvVq1fJyMggODiYTz75hGeeeeau923Tpg3PP/88c+fO5aOPPiI8PJzr16+zYcMGpk6dajley5Yt2b9/v+XXyeKYNWsW7u7uPPDAAxw7dozk5GR+/PFH3nvvPe6//3569uwJmL+I/Pjjj1y+fBmDwUBwcDDr1q2jUaNGlm5cTZs2RavVsmXLFjIzM63OM3LkSIKCgvj+++9JS0vj6NGjjBkzhoEDB1ptd6tjFObVV1+lQYMG3HfffWzZsoXk5GRCQ0P56KOPLL+UFyf223FycmLDhg3ce++9jBgxwlIRN2rUiLfffpuPPvqI999/n8uXL5ORkcGZM2f4/PPPmTJlSrGOX1xPPvkkPXv2ZMyYMfz888/ExsaSlJTEgQMHeP755y3Tq5XEiBEj2LVrV5Hd/J2cnHjzzTct3c8mTpxIdHQ077//vtV25XVP9z///MO9996LRqO5q+MIIcTNHnvsMZo1a8bkyZPZvXs3KSkpbNiwgZ07dxb44bgwEydO5IsvvuDcuXMYDAauXLnCsmXLcHZ2pmPHjlbb5k2/FRERQUREBNOnT+fMmTPMmTPHss24ceMYPnw4kyZNYsWKFURFRZGSksKRI0d49dVX+eSTT0o1ppvfizZt2vDkk0+yfft2UlJS+Oeff5gyZQotW7bkiSeeuO25S0qlUvG///2PN998k+eff57Zs2cD5p5U33zzDZs2bWLq1KmcPXuWzMxMQkJC+P777xk8eHCpxtG7d2/+7//+j5deeokvvviC8PBw0tLSOHnyJB988AGvvPJKiY85YsQIQkNDraZ6y69t27Y8/vjjXL58mdjYWN5++222bdvGBx98YNUl/E7u6QaYN28e165dY9q0aURGRhIdHc2MGTM4efIkc+fOLdDS/ffff9O0aVMaNWpU4msVVUipDckmRClZtWqV0rt3b8XFxUVxdHRUWrZsqbz66qtKeHi4ZZvCpgwr7r6KoiiLFy9WOnTooDg4OCje3t7K6NGjlVOnTlnWnzx5UunWrZvi4OCgAFajQN881Ul+Fy9eVMaOHat4enoqOp1OadiwofLOO+9Yje556NAh5bHHHlPq1q2r2NvbK3Xq1FGmTZtWIMYFCxYo9erVUzQajdUUIEajUZk1a5ZSq1YtxcHBQenSpYuyY8cO5ZVXXlFu/i9d1DEKG71cURQlPj5eee6555Q6deoodnZ2SsOGDZXXX3/dMtVWcWO/2c1ThimKefqWKVOmKGq1Wlm4cKGlfN26dUq/fv0UNzc3xcHBQWnRooXy0ksvKVeuXLFs4+npWejfP2+U8HXr1hVYV9jfLTMzU/nwww+VVq1aKfb29kqNGjWUrl27Kl9++aWSlpamKErxRy9XFEU5e/ZsoSPU5o9r/vz5Su3atRW9Xq906tRJ+fvvv2973OIaO3ZskaOvvvfee4XGtGPHjlI7vxBC5Hft2jVl3Lhxiqurq+Lq6qo8/PDDSmJiYrGmDAsJCVFeeOEFpWnTpoq9vb3i7++vPPjgg8rRo0ct2+QfrXvJkiVK/fr1FTs7O6VNmzbKb7/9ViCe7OxsZf78+Ur79u0VR0dHxc3NTenYsaMyd+5cJTExscAxb1acmAq7FkUxT9/19NNPKwEBAYpWq1X8/f2V6dOnK/Hx8YVez80CAgKsRvouTFH7z58/X1GpVMpzzz2nmEwmRVEUZf/+/cro0aMVb29vS30/adIkq2sZMGCA0qFDhwLnyRsl/OYZYxTFPCJ4YTNnLF26VOnWrZvi7OysODs7K23atFHefvtty1RtilL8KcPS09MVDw8P5d133y0yrt9++01p3ry5YmdnpzRp0kRZunTpbY9bEvv27VMGDhyoODs7Kw4ODkq3bt0K1P2KoijXr19X9Hq9Mn/+/FI9v6h8VIpSyA0mQggh7siYMWNITk626u6/b98+unXrxrp164o9cm9Ze/TRRwkJCSlwb70QQlQWFy9epFGjRixatKjUe0OJiu3999/n22+/5eLFi5aeE5mZmTg4OPD2229bWvZt7dNPP2Xu3LlcuHDhliOui6pPupcLIUQp+vjjj9m5cydBQUG2DqVIZ86cYdWqVcyfP9/WoQghhBAl9tJLL6FSqfj2229tHUqR0tLS+Pjjj3n//fcl4RZob7+JEEKI4mrQoEGJB+Erb82aNbOaF1wIIYSoTBwcHLh69aqtw7glJyenQkfwF9WTtHQLIYQQQgghhBBlRO7pFkIIIYQQQgghyoi0dAshhBBCCCGEEGVEkm4hhBBCCCGEEKKMVJuB1EwmE5GRkbi4uBSYtF4IIYSoyBRFISUlBX9/f9TqqvV7udTPQgghKqvi1s/VJumOjIwkMDDQ1mEIIYQQdywsLIxatWrZOoxSJfWzEEKIyu529XO1Sbrz5scLCwvD1dXVxtEIIYQQxZecnExgYGCVnOtV6mchhBCVVXHr52qTdOd1WXN1dZVKXQghRKVUFbtfS/0shBCisrtd/Vy1bgwTQgghhBBCCCEqEEm6hRBCCCGEEEKIMiJJtxBCCCHuSkREBHv37iUmJsbWoQghhBAVjiTdQgghhLgjJpOJyZMn07BhQ6ZNm0bt2rV55ZVXbB2WEEIIUaFUm4HUhBBCCFG6vvrqK3799VeOHz9O48aN2b9/P7169aJTp06MGTPG1uEJIYQQFYK0dAshhBDijixdupQHH3yQxo0bA9ClSxcGDRrE0qVLbRyZEEIIUXFIS/cd2Hl4Pd8ceof3711Jg8CWtg5HCCGEKHc5OTmcPHmSJ5980qq8Q4cOfPfdd0XuZzAYMBgMltfJycllFmNFl23KxphjJMuUhTHHiNFkND8WUpZlysKkmMhRcjCZch8VEybFRLaSbVVW2KOiKCgoAChK7uNNr28uszzPfUQpuE/+bfIfx2q/ykJRQDHlLjlgUgATmHJy1+XkW5+7mEzmbRQwv0GK9eu8MiXf8wLlFCwn37p87/9NAVvHXvhFFf5aKWp9YcdSinhZyf6+otK58XlDvkdzaf5/poVvl1uu3PiXmv+f7r2tJtGt7diyC/4mknSXkCknh8+OvEmIPczaPIGlE3bh5HjrydCFEEKIqiYlJYXs7Gw8PDysyr28vLh+/XqR+82ZM4d33nmnrMMrM4qikJaVRpIxiWRDsuUx2ZhMkiGJtKw00rPTSc9KtzxmZGdYlWVkmV/nKDm2vhwhhKiW/C4dlKS7IlNrNPyn00e8dPRlzujhhZVDmP/YNhwdnGwdmhBCCFFu7OzsAMjIyLAqT09Pt6wrzKxZs3jxxRctr5OTkwkMDCybIG8hx5RDalYqiYZEEg2JJBmSLIvldb6EOi+pTjGmlEmyrEKFncYOO7UdOo3O8txOY4dOrUOn0aFVaVGr1GhUGtQqNWr1jedFPqpzt0WNSqWyOh9gKVOhsnpu2e6mskK3szyoCj1uAYoJstLBmAbGdMjOgKwMyMqE7Mwbr7Mzb5TlGAoepzyo1LmLxvq5WpVbpgGVKndRA6rc9yO3LP8judvcXF5g+7xtcsvzPWD1ft7puryim8sK2ea2xylsH1FqCvvvo0C2omAyKeTkLtm5jyblRplJUTApmLfL3d6kkFuet435dWVl+V+S+39IpbIuu/F5VLDMr0678gxVku470bPtMJ6KPM5nMT+y1z6FaSt68e7wH6kb0NTWoQkhhBDlwtHRES8vLyIiIqzKIyIiqFOnTpH76fV69Hp9mcS0J3IPF69fNLcq529lzn1MMaaYE+rcZPpuuj/bqe1w07vhaueKq94VNzs3XPWuOOmccNI54ah1xFHniKPWEQedg/l1/jKtA3qN3pxU5ybUqgJJUCWjKJCZCIlhkHgVkiMgNRpSoiE16sbz9Dhz4n0nVBrQu9xY7JzzvXY2v9bag84h36MetA6gs7d+1OpvrNfYgVoHGm3uow7U2kISUyFKR4Yxh4R0I9fTjFxPN3I9PYvEdCMJaUYS07MsZckZWaRkZpGSmU1KZjYZWWXXQ0avVWOv06DXqrHTqtFr1ei1Gstzu9zXep0avUZtfsy3Pm8brVqNTqNCq1GjVavQadRo1CpzmVqNVmMu06rN2+Qvz9tem1um06jQqM3P1WrQqMyvK9vnpSTdd2jiva+R84eJBTGrOGqfxeQtD/Cg6z1MGf0RWq28rUIIIaq+QYMG8fvvvzNr1iwAsrOz2bhxI/fff79N4tl4aSMbQjaUaB9HrSPuenfc9G646d2snucl0nmPrnaulkTbXmtfRldRwWUbICEU4i9CQggk5SbYeYm2MaV4x1GpwckHnL3B0RMcPMChBjh63PS8hnmxdzMn1DoHSYRFhaQoCmnGHGJTDMQkZxKbaiA25cYSk/uYkJtkG7Lv8IenXA46DS72WpzttbjY63C11+Jir8VFr8NJr8XRToODnQbH3MVep8HR7qZyndby3EGnQa2W/1tlRbLDu/DEsDeodaARn5x4j2s6NV+l/8nOxVt5qM7/cd+gp1FrZHB4IYQQVddbb71Fp06dmDx5MsOHD+f7778nPT2dmTNn2iSetj5tyVFyrFqT85476hxx0jnhrne/kVjbuaHT6GwSa4WmKOZW6Zgz5uQ6/5J49fat1I5e4B4IbrXA2Q9cfM2Pzr43njt5gVpTPtcjRClIycziWlImEYkZRFoW8+uopExiUwwlboW206hxd9RRw9GOGk7mR3dHOzzyPXfNTapd7LW45j4622vRSZ5RqaiUm4earKKSk5Nxc3MjKSkJV1fXUj12QnIMH6+bxnbTeYy5vxC1zlDo6dyXMUNex9u7ZqmeTwghRPVSlnXY3Tp9+jTz5s3j6tWrNGrUiJkzZ1K3bt1i71+Rr61aMKZD7BmIPg3RwRB9CmJOQ3p80fvYuYBXQ/BoAO61cxPs2jcSbTsZ50ZUPplZOYQlpHM5Pp0r8WlciU+3JNgRiRmkZGYX6zjOei3eLnq8nfV4u+Y+utxYvJz05kTbyQ4nO02l6yYtrBW3DpOkuxSduXKEz7a/xEFNLDm5/4F8snPomOlCB6/B9Oz1JP7+tcrk3EIIIaquqpyYVuVrq3CyMuDaCYg4DJFHIPIoxF+i0KmfVGpzUu3VGDwbgGdD8+LVCJy8pYu3qJSyc0xcjk/nYkwqV+LTuByfzuW4NK7Ep3EtObPomddyuTvqqOnmQIC7Pf7uDpalpps9PrlJtaOddCSuTqps0p2VlQWATley7mDlWamHRJ/h6+1vs8t4hrR8PT/qGbNok+FAffs2NG80igbt+uPlLl8whBBC3FpVTkyr8rXZlMlkbsGOOJy7HDG3ZBc28rqjF/i2AN+WuY/Nwbup+f5pISqhG8l1CuejUzkfncKF6FRC4lLJyik69XHWa6nr5UgdTyfqeDgS6OGIv7s5ya7p5oCTXhJqYa24dVil+pfz1ltvMW/ePFQqFS+//DJvvPGGrUMqVH3fZsx9ZA3GHCO/HVzBH+dWc0KJItROR6hdNnAY99ADdA02Emhwx8+uOZ7+PfFu2oMGDRrjZC/3lwkhhBCiBLKNcO0YXNkNV/ZC2D7ITCq4nZMPBHTIXdqBX2tw9in3cIUoLUkZWZyOTCY4MongyGTOXEsmJDYNY07hYw846DQ09HGmnpcTdT3NCXZdL0fqejrh4WQn3b1Fmag0Ld1Hjx7l/vvv5+DBgyiKQseOHdm0aRMtWrQo1v62/iU9xZjC5pOr2X5+A8czL5Outv4g8M/KpmNmJk0yNNTICkDj0Jgc7xY41WqNX4PW1Pergb1OBhwRQojqyNZ1WFmqytdWprIyzYn1lT3mJfyQeX7r/Oycwb9dviS7PbgGSNdwUWldTzNyPDyR4MhkTkWYk+yrCemFbutop6GRjzMNfVxo7OtMI19nGvm4EODuIKN0i1JTIVu6g4ODOXr0KD169KBevXoF1mdnZ7N7925iY2Np27YtDRs2tKwLCgpi9OjReHl5ATB69Gh27txZ7KTb1lzsXHiowxQe6jCFbFM2x2KO8s/Z39gfuZ+LxmgidVp+1zmDC0ACvtn/0ibpL1rHGHHdk80FgyfXNHVJdqxDlns9tN6NcPZvSq2aNant6Yibg7SOCyGEEFWWyWQe5CzkH7j0D1zdC9mZ1ts4ekLtblCnB9TpBr6tzPNOC1EJ5ZgUzkWlcOTqdY5cvc6xq4mExKUVum2tGg609Hejhb8rzf1daewrybWoWMrlk/jAgQP85z//ISYmhrNnz7Jy5coCSXdMTAwDBw4kJSWFJk2asGvXLmbOnMk777wDQEJCAjVq1LBs7+HhQUJCQnmEX+q0ai0d/TrR0a8TAOlZ6RyLOcbByD3suxrEmdTLRGu1bHXWstXZPAKoVlFoarxA68xgWl830DrKQK2jOcQrrlxQ/AhX+3NdH0i6Sx1M7nXRe9Wlhqcv/jUcqelmvg/FwU5ayoUQQohKIyniRpIdsgPS46zXu/hDvV65iXZ386Bn0ootKql0YzaHr1xnf0gCR65e53hYImnGgmMQ1PdyolUtc4Ld0t+N5v6uuDva2SBiIYqvXJLutLQ03nnnHfr27VvkfRIzZ85ErVZz+vRpHBwc2LZtG4MHD2bo0KF069YNPz8/jh8/btk+PDyczp07l0f4Zc5R50j3gO50D+jO851mkp6Vzqm4U5yIO8HxmOOciD1GgiGRU3o9p/R6fjQ3h+ORk0NLg5GWhihaGa7SzxCEW5wJ4oCLkKrYE654E6J4sUvxIl7nR6q9PwanWpjcAnF098Hb1R6vvKkMnPV4udjh6aRHI78MCiGEEOUrJ8vcgn3+T7iwFeLOW6+3c4a6PaF+P2jQT5JsUallZuVw5Mp19obEs/dSPMfDEwsMcuZkp6FtbXfa165B+9o1aBvoTg0nSbBF5VPu93SrVCpWrlzJhAkTLGVGoxF3d3c+/fRTpk+fbilv2bIlffv25csvv+Ty5ct06tSJX375BUVReOihhzh69CgBAQGFnsdgMGAwGCyvk5OTCQwMrJT3jCmKQkRqBCfjTnIi9gQnYk9wOuE02aaC8wUG5KhpYTDSLiOFVgYDTY1G9EX8hdMVPeGKF9cUT64pHkThwTXFk2jFg3QHP7Kda6J3dMfDWU8NJx0ejnbUcLLDw8mOGo65i5MODyc7HHQyz6AQQpSVqnzfc1W+tmJJi4eL2+D8Frj4FxiSb6xTqcG/vTnBrt8PanUCrSQconLKMSkcD08k6Hwsey7Fc+xqYoHBzvzd7Ola35OOdT1oV9udxr4u0hAkKrQKeU93US5evEhGRgYtW7a0Km/VqhUnT54EoG7dunz55Ze88MILqFQqvv766yITboA5c+ZYuqZXdiqViloutajlUot76t0DgCHHwNmEs5yKO8XJuJOcijvFleQrRGhMRDhq2epo7oqvValprPemmcqRppk5NE9NokFyNE6ZMTiqDDRWRdCYiIInzQGSIDXRnijFw5yUKx5E4sl5y2tzsp6EE3qtxpKMeziZk3M3By1uDjqrxTXv0V6Hm6MOF71WknUhhBDVh6JAzGlzkn3+Twg7gNU82Y5e0HgINBoM9fuAQ40iDyVERRebYiDofCw7zsey60IsielZVut9XfV0q+9JtwaedK3vSW0PR/leKKqkCpF0JyWZp7TIf882mO/bPnPmjOX12LFjGTt2bLGOOWvWLF588UXL67yW7qpCr9HTxrsNbbzbWMqSDElWSfjJuJMkZCZwOjOa03kbOYGzuxctPHrSyqUOLXXutEKPT0YKJEdgSoogJzECVUokWkMizqpMGqoiaUhkkbFkKHZcUzyITvcgKr0G0bF5SXkNTuU+xuJODgXvKVersCTibvkS8pvLzAm7Fld7HS72WlxyH/VatXw4CyGEqNiyMuHyvzcS7aSr1ut9W5kT7Sb3mFu21WrbxCnEXVIUheDIZLaejuafszGcjLCets7FXkvvRt70aOhF1/oe1PNyku9xolqoEEm3g4MDAKmpqVblKSkplnUlpdfr0ev1dx1bZeKmd6NHQA96BPQAzB98kWmR5iQ81pyEn0k4Q2pWKvujD7E/+pBlXx9HH1p7taZl3Xto5dWKFl4t0CpA8jVIjoDkyHyPkZAcbn5Mj8dBZaS+Kor6RBUZmwk119U1iMWDa0oNwnNqEJlTgyilBlGZHkRleHBBqUEG9iW6Zp1GhYu9Dme9NjcZ1+Ks1+Ga9zxfgu6sv5G0W5XbaWV0SyGEEKUrJRou/AnntpgHQ8vKN62R1h7q9TEn2o2HgFst28UpxF3KMSkcupzAn8HRbD0dRfh166nrWvi70reJN32b+NAu0B2tRn5UEtVPhUi669Wrh0ql4sqVK3Tt2tVSfuXKFerXr2/DyCo3lUpFgHMAAc4BDK07FIBsUzaXEi9ZtYZfTLxITHoM269uZ/vV7eZ9UdHAvQEtvVrSyqsVLQPb0KjGGHTqm6Ymy8qElMjc5DzyxvOUvOT8GqRGoTZl42mKx5N4mgKoc5ebGLUupOm9Sdb5kKDxJE7lSZTiQURODa5mu3HF6EqYwZEUowlFgawchYQ0Iwlpxrt6r/KS9hvJuw5ney2uuWVOevOjs96csFu9zrfeTisViRBCVEuKAtHBcG4znN8MEYet17v45ybZQ6Feb7BztE2cQpSC7BwTuy/Fs+nENbafiSY+3/cwe52a3o28GdTclz5NvPFxKVmDihBVUYVIut3c3OjRowdr1qyxdB8PCwtj9+7drFixwsbRVS1atZYmHk1o4tGEMY3HAOYpy07Hn7bqmh6ZFsnFxItcTLzI+ovrAXOX9mYezSyJeCvvVtRyroXKoz543OLHEVMOpMWZW8pTrt1oLc97nvdoTMUuOwW77BRqEEKdoo6n0aH41iTHuSZZDr5kOviSqvchWedNosaLOLUnMdQgyagm1ZBNcmYWKZnZpGZmk2LI9zwz2zKAR6ohm1RDwYHpSspOo85NyjU42ZkTeKfcpNzlpuTdKTeBd87d1vmmBF+6zgshRAWXbYDLu8yt2ee3QFKY9Xr/dtD4HmgyFPxay0jjolIzmRSOXL3O78cj2XTyGnGpNxJtNwcdA5r5MKSFH70becs0tULcpFyS7qioKLZv3255vWfPHgAaNmxoadn+9NNP6du3L5MmTaJdu3Z8++23dO/evdj3cIs756hzpKNfRzr6dbSUxWXEFbg/PMWYwrHYYxyLPWbZzl3vbknC8+4xd7Zztj6BWgMuvublVjKTb0rK81rN87q4X4O0WDBloUq8ijbxKlrAASh0mBlHL3CtaW5dqOEPrv7gUtNc5loLXGqSqXEmxZBDqiGblNzk3LzceJ5mzH00mJcUw43necl6ZpY5eTfmmHJb3u/kL2FNq1aZE3i7m1vYNZbkPC+Rd8ptoXeyy/c8X4Jvr5MEXgghSkVKNFz6C85tMs+fbcx3a5zWAer3NSfZjYaY6xshKrmzUcmsPxrJhuORRCTe6Dru4WTHva38uKdlTTrX80An3caFKFK5TBkWHBzMnDlzCpT369ePyZMnW16fPXuWpUuXEh8fT7t27ZgyZUqp3Zdd7ackuUsmxcTV5KtWSfjZhLNkmaxHoVSr1DSp0YR2Pu1o79ueDr4d8HLwKr1AcrIgJarorux5ZTmG2x8LQOeYm4jnT8rznvubvzA5+5p/OLiF7BwTaYYcUo3mVvTUm5LyNENuee5687ocUg1ZpBlyrJL5dGNOKbxR1jRqFU52mkK7xxfWdf5GMq/BRa/DSa+xtMTL9HBClL+qXIdV+GvLNpjnzr70N1z8G6JPWq939ssdBO1e6TYuqoyk9Cx+Px7BmkPhVoOhOdlpGNLSj5Ft/OnR0EsSbVHtFbcOK/d5um2lwlfqlZAxx8j56+c5GXeSk7EnORpzlPDU8ALbBboE0t7HnIC382lHHdc6ZZu0KQqkJxSSlOd1Zc9tOc9MLN7xVGrzlyrXmoUn5XmPdk6lEn6OSSHNmC9RN2SbE3pDljlRz8wizWhunb+RwN9I7m+0xOeUSpf5m6lVFNnKnpeY36r7fP4EXxJ4IYqnKtdhFe7aTDkQdRKu7Da3ZF/ZbT0IGkDNtjfuz67ZVkYbF1WCyaSwNySe1QfD2BIchTHb3ItPp1HRr4kPo9oF0L+pD/Y66TouRB5Jum9S4Sr1KiomPYYjMUc4Em1ezl8/j4L1PzEPew/a+7Sno19HOvt1pqF7Q9skXsZ0cxJe2H3mec9TokApZsuzvduNBLxAUp67OHqW6z19JpNCelZuS/pNLfD5k/W85D3FqjzHkvinGcyt9KX9aZGXwBfW6n6jm7wGZ73O3K2+kO72ecm/o50k8KLqqsp1mM2vzZgG107A1T1wZS+E7QdDsvU2zr7QoD80GAAN+oFTKfbgEsLG4lINrD4Yxk8HrlqNPN7E14WHOgUyul0AHk52NoxQiIpLku6b2LxSr6ZSjCkcizlmScRPxZ3CaLIeadzD3oPOfp3p5NeJLjW7UNuldsVJnkw55vvI8+4pL+ye8+RIyCrmTdwaO3Dxu5GIuweCW+6S99y+Yv77NJkUMrJyikzUrZL2TOvu81b3wWeWYQJvdyMZv5GQ35S039xKb9Uqb14kgRcVTVWuw8rt2hQF0uMh9hxcOw7Xjpkf486DYrLeVu8KgV2gXi9zou3bQgZBE1XOsbBEVuy5zMYT1ywDy7rotYxs68/YToG0CnCTulCI25Ck+yZV+QtLZWLIMRAcF8yRmCMcuHaAozFHyczJtNrG19GXLjW7mJNwvy7UdK7gA9EoirlVJP895YXdc54WW7zj2buBW23zvK3u+RPy3DJnn0r/5U9RFNKNOYUk6uaB6wok7YV0o0/Nt5T2p5hKBc52N+5rd7bXFTqAXf6p5fLmiM8/J7yzvRaNzAEvSkFVrsNK9driLkBCqPnzNi3GPHNGyjWIv2QuNyQVvp+zrznJrtPdvPi2vO1YHkJURplZOfxx4hor9l7mePiN/w9tAt15rGsd7m1VU0YeF6IEJOm+SVX+wlKZGXOMnIw7yYFrBzgQdYDjsccLDM4W6BJIZ7/OdPXvSrea3XDTu9ko2ruUbYTUqBv3lCdHmqeXSQwzPyaFQcb12x9Ho78pIc9NxvOSc9cA0Ohuf5wqQlFyW+AzrRPx1NyR5833wRfSnb6IUelNpfyJ6GinyTcHvC53HngtLnpzUp43L7xL3lzx9vm2y91HRp8XVbkOK9Vr+/lxCP711tu41YaaraFmmxuLi9/dnVeICu56mpEVe6+wYu9ly5zadho1w9vU5LFudWkb6G7bAIWopCTpvklV/sJSlWRkZ3As5hgHow6yP2o/wXHB5OS7p1qtUtPSsyXdA7rTw78HLb1aolVXiOnmS4chBZLCrRPxxHyPKdeA2/yXVanNA765BUKNOlCjLrjXMT93r2Pu1i4tOIWyJPCW7vM5pBQywnxegp83vdzNz5Mzsy0D0JSGvOnj8rem5yXpea3t+ZP0wpJ8Z70WrYwyW2lV5TqsVK/tnzlw7g9w8gEnb3D2Nrdi16gHng3Mn4c6h1KJW4jKICwhncX/hrL6YBgZWebvU/5u9jzStQ7jOgXi6Vw6swQJUV1J0n2TqvyFpSpLy0rjcPRh9l/bz57IPVxMvGi13sXOha41u9Ld35yEV/iu6HcrJ8vcSm5JysMh8Wq+5Dz89tOlqXXmFvG8RNwqKa8Ljh6Vvvt6RWDINrew5w1gl5yZZXltlajnT9hv3r6Uu8476DRWremuluQ877UOVwdt7qMut9t8bpmDDmc7LWrpLm8TVbkOq8rXJoStBEcmsTAohI0nrpGT24Wrhb8rT/ZpwL0t/eRHWCFKiSTdN5FKvWqISotib+RedkfuZm/kXpKN1iPM1nOrRw//HvQI6EEnv07oNdXsF1yTCdLjzAl44hW4fjn3Mfd5UhiYbjONmJ3LjVbxwlrKZQ7acpM3+nxeQp58U8Kemvs62dLynr/V/UYybyilVneVCktrel5Sbn40J+Z5iXz+Mtd8yby0tt+5qlyHVeVrE6K8HQtLZP728+w4d2McmZ4NvZjWpwE9GnrKrUpClDJJum8ilXrVk2PKITg+mN2Ru9kTsYcTcScw5RuB1kHrQJeaXehTqw99avXB29HbhtFWEKYc873k+ZPxvOT8+hXzPee34+QDHvVzl3rWjw41yvoKxB0wZpssSXlKvtZ0S3f4DHN5cqa5hT05w5zIp+Q+JmdkWUa2vVt5Leu3bFW3lJm7x+dfb6etnkl7Va7DqvK1CVFeToQnMm/bef7JTbbVKhjW2p8ne9enZUAlHQtHiEpAku6bSKVe9SUbk9l/bT+7I3bzb8S/RKdHW61v7tmcvrX60juwN809msuvvYXJyjB3V7dKxi/fSNBvnrv2Zg418iXk+ZYa9czz2sp7XmllZuXkJuQ3WtfNyXn+MvPz5NzkPv/6vHsJ71ZeF/mbW9rdHKwXVwdzAp+/zFmvrbT/76tyHVaVr02IsnYqIon528+z/UwMYE62R7erxbP9G1LXy8nG0QlR9UnSfROp1KsXRVE4f/08O8J2EBQexMm4kyj5BiDzcfChV61e9KnVh67+XXHQysA6t6Uo5tHVr1+G66GQEGKegich9/ntWsntXPK1jN+0uPhJQl7FZeWYCiTi5uT8xnNLC3shZamG29wWUQwatQpXe61VYn5zsn5z4u7moMPN0Ty6vC0T9opchyUnJ3PkyBGysrJo3bo1vr6+Jd6/ol6bEBXV+egUPvnzHFtPmxsY1CoY1TaAZwc0op4k20KUG0m6byKVevUWlxHHrvBdBIUHsSdyD+nZ6ZZ19hp7egT0YEDtAfSu1bvyTklma4ZUc0KeEGJe8ifmSeHcctR1rcONhNyrEXg2BM/cRyfP8roCUYHlmJTc+9qzSMpN3POS+KSMLMtjwaV0userVRRI0m+XtLs56PB20WOvu/vZAipqHfbWW2+xePFiGjVqhEajYe/evbz55pvMmjWr2MeoqNcmREUUnZzJZ1vP8/PhMEyK+ffq+9r48+yARjTwdrZ1eEJUO5J030QqdZHHmGPkUNQhdoTvYGfYTiLTIi3rtCotnWt2ZkDtAfSv3R8vBy8bRlqFZGWau6jntYrnXxKvgnKLrscONfIl4Q1uJOUe9WXqH1EsiqKQmWUqNClPLuJ5/uVuBqKbO6Y1D3YMvOtrqKh12KJFi3j44YdxcjK3rG3YsIGRI0eyb98+unTpUqxjVNRrE6IiSTVks3DnJRbtCrXcrjO0hR8zhzSmoY+LjaMTovqqkkl3ixYtMBjM0yF5e3uzd+/eYu8rlboojKIonLt+jr+u/sX2K9utpiRToaKtT1sG1B7AgNoDqOVSy4aRVmE5WeZR1eNDIOESxF2A+IvmJSnsFjuqzHOR50/E8xa3QFBXzwG3ROnLzMq5RUv6rZP2z8e1Y0gLv7uOoTLVYXq9ni+++IKpU6cWa/vKdG1ClLesHBOrDobx+fbzxKUaAehQpwav3duMDnVk8FIhbK1KJt0+Pj7s2bMHAK1WS926dYu9r1TqojguJ13mr6t/8dfVvzgZd9JqXVOPpgyuM5ihdYcS6Hr3LVeiGIzp5tbw+IsQfwHi85LyC5CZVPR+Gr05GfdsaE7IvZqAd2Pwagx2cq+bKD+KopTKveCVpQ4LCgqiT58+7N69m+7duxe6jcFgsPyADuZrCwwMrPDXJkR5CzofyzsbgrkUmwZAPS8nXhnahCEt/CrtoJBCVDU2SbqzsrKIiIjA29vb0tXsZhkZGSQlJeHr61viDwwXFxcaNmyIu7s7L774IiNGjCj2vpXlC4uoOKLSovjr6l/8ffVvDkUfspqOrIVnC4bWHcqQukOo6VzThlFWU4oC6fHmZDx/y3j8RXOSnmMsel+32rkJeJN8j03A0aP84heihMqrDjt58iTBwcG33KZ///74+PgUKI+Pj6dLly60aNGC3377rcj9Z8+ezTvvvFOgXOpnIczCr6fz/sYzbAk2D1Dq4WTHjIGNGN+5NjqN9OISoiIp16Q7MjKSr776iuXLlxMREcHKlSuZMGGC1TbZ2dk89dRTLF++HL1ej4uLC99++y3Dhw8HICIigj59+hR6/J07dxIQEEBISAjZ2dkcP36c6dOns2vXLpo1a1asGCXpFncjITOBf67+w5bLWzgQdcAqAW/r3Zah9YYyuM5gmQu8IjDlmO8Tz0vI486bl9hzkB5X9H5O3jcl4o3Buym41JSR1YXNlVcdtmbNGn799ddbbvPWW2/RvHlzq7KkpCQGDhyIRqNh27ZtuLgUfY+ptHQLUbjMrBwWBYXw1Y6LZGaZ0KhVTOxWlxmDGuFqr7N1eEKIQpRr0r1y5UquXLnClClTqFmzZqFJ95tvvsl3331HUFAQDRo04NNPP+XNN9/k9OnT1K9fn+zsbC5fvlzo8evWrYtWq7Uqmzp1Kl26dGHy5MnFilGSblFa4jPi2X5lO5svb+ZI9BHLVGQqVHTw7cDQukMZVHcQHvbSclrhpMVD3DlzAp6XiMedv/W943YuBVvGfZqCe125b1yUm4pchyUnJzNo0CAAtm7diptbyWaAqMjXJkR5+ftsNO9sOM2VePPsKl3qefDOfS1o6if/J4SoyGx2T7dKpSqQdJtMJnx9fXn++ed54403APN9brVq1eLxxx/n/fffL9E54uLi6NWrF19++SUDBgwo1j5SqYuyEJMew9bLW9lyeQvHY49byjUqDT0CejCi/gj6BvbFXmtvwyjFbRlSrVvE8x4TQooeWV3naO6W7t0MfJqBT3Pzo6u/tIyLUldR67Dk5GQGDx6Moih3lHDnHaMiXpsQ5SEqKZO3fjtlmW/b11XP68OaM6J1TblvW4hKoLh1mLbINaXo8uXLxMXFWQ2qolKp6NatGwcPHizWMXbt2sXjjz+OyWQiOjqaqVOn3jLhLqz7mhClzcfRhwnNJzCh+QQiUyPZenkrmy9v5nT8aYLCgwgKD8JJ58SgOoMYUX8EHf06olZJ62iFo3eGgPbmJb9so3lEdauW8XMQex6y0iHyqHmxOpZbbhLe9EYi7tMcnGT6OVH1DB8+nODgYD755BM2b95sKW/ZsiUtW7a0YWRCVGwmk8Kqg2HM2XSGFEM2WrWKyb3q8Wz/Rjjry+XruRCiHJXL/+rY2FgAPD09rcq9vLwICQkp1jE6dOjAli1b0Gg0+Pv7o9frb7n9nDlzCh2oRYiy4u/sz6SWk5jUchIhSSFsvLSRP0L+IDItkvUX17P+4np8HX0ZVn8YI+qPoGGNhrYOWdyO1i43ab5p7AhTjnnO8dgzEHMGYk6bH+MugCEJwvaZl/ycvK1bxL1zE3P7krcMClFR1K1bF39/f/755x+rcq1WK0m3EEUIiU1l1q8n2R+aAECbQHf++0BrmvjJfNtCVFXl0r384MGDdO7cmUOHDtGhQwdL+eTJkzlx4kSxW7tLQgZqERWBSTFxNOYoG0M28uflP0kxpljWNfVoyvD6wxlWfxheDtIKWiVkG8wDuFkS8bPmx+uXgSI+al1r3Ujs/VqBb0vzNGcaGTRH3FCVu2BX5WsTIr+sHBOLdoUwf/sFjNkmHHQaZg5pwqTuddGopSu5EJVRhepeXqtWLQCioqKsyqOioggICCiTc+r1+tu2hgtR1tQqNR18O9DBtwOvdn6VoPAgNl7aSFBEEGcTznI24SzzD8+nV61e3N/ofnoG9ESrlm5llZZWD74tzEt+xjRz1/T8reKxZyE5ApLDzcvFbTe219iZ7xf3bWU+ll9LczIuXdSFEKJSOnMtmZfWHOf0NfPtjr0aefHh6FYEejjaODIhRHkol2/3NWvWpHHjxvz5558MGzYMMM/XHRQUxLvvvlseIQhhc3qNnkF1BjGoziASMxP58/Kf/H7pd07EneCfsH/4J+wfvBy8GNFgBKMajqK+W31bhyxKi51T4feMZySak+/oYHMyHh1sXgzJEHXSvOTn7HcjAfdtaX7u2Qg08kONEEJURDkmhe92hfDp1vMYc0y4Oeh4c3hzHmgfIAOlCVGNlEr38szMTEsrdr169Zg3bx6jRo3CxcXFch/3jz/+yKRJk1iwYAHt2rXjww8/5MCBAwQHB5dLdzLpviYqqkuJl1h3YR0bQjaQkJlgKW/n047RDUczpO4QHHXyS3i1oSiQeAWiTuUm4SfNz6+HFr69Rm9uFffLbRX3bWl+7ihT1lUlVbkOq8rXJqq3sIR0Xvr5OAdy790e0NSHOQ+0wsdFZjQRoqoo1ynDdu/ezSOPPFKgfNy4cXz00UeW1z/88AMLFiwgPj6edu3a8cEHH1C/fvm05kmlLiq6LFMWQeFBrLuwjl0RuzApJgActA4MrTuU0Y1G09a7rfwyXl0ZUnNbw0/lJuS5SbkxtfDtXWqak2+/1lCzNdRsA+51ZDqzSqoq12FV+dpE9aQoCr8cDuedDadJNWTjaKfhreHNGdspUOpwIaoYm83TXVFJpS4qk5j0GH6/9DvrL67nSvIVS3lD94Y81OQhhtcfjoudjHJa7ZlM5lbxvAQ86qT5+fXLhW9v75abhLcxL36tzYO2qTXlGrYouapch1XlaxPVT3yqgVm/nrTMu92hTg0+e6gNdTydbByZEKIsSNJ9E6nURWWkKApHY47y64Vf+fPyn2TmZALm1u97693Lg00epIVni9scRVQ7hhSIPg1RJ+DacfNj9GkwZRXcVudo7pael4TXbGMeSV0rA1FWJFW5DqvK1yaqlz0X43h+9TFiUwzoNCpeGNSYJ3s3kJHJhajCJOm+iVTqorJLNiaz4dIGfj73M5eSLlnKW3q25KEmD8m93+LWso3mQdvykvBrJ8wt41lpBbdV68xziPvltojXbG2+V1zvXP5xC6Bq12FV+dpE9ZBjUvjfXxf4398XUBRo6OPM/LFtaRngZuvQhBBlTJLum0ilLqoKRVE4EnOENefWsO3KNrJyWy9ddC6MaDCCBxs/SMMaDW0cpagUTDmQEGJOxPOWqBOQcb2QjVXg2dCchPu3My8120giXk6qch1Wla9NVH3RyZk8v+oo+0LMg6U91LEW74xsiYOd3LYjRHUgSfdNpFIXVVFCZgK/XfyNn8//TFhKmKW8o29HHmn2CH0D+8q836JkFAWSwswt4ZZW8eOQcq2QjVXmkdP924F/e/OjX0vQOZR72FVdVa7DqvK1iapt5/lYXlx9jPg0I452Gj4Y3ZLR7WrZOiwhRDmSpPsmUqmLqsykmNh3bR8/n/uZf8L+IUfJAaCmU03GNhnLA40ewN3e3bZBisotNSY3ET8Kkccg8igkRxTcTq013xOePxH3aQ5au3IPuSqpynVYVb42UTVl55j4dNt5vt5hvtWrqZ8LXz3Sngbe0vNHiOpGku6bSKUuqouotCjWnFvDL+d/4brB3E1Yr9EzvP5wxjcdTxOPJjaOUFQZKdHm5NuyHIG02ILbafTmFvC8bun+7cGrMWikF0ZxVeU6rCpfm6h64lINPPPjEUt38ke61ObN4c2x10l3ciGqI0m6byKVuqhuDDkGNodu5oczP3A24aylvJNfJx5u+rB0PRelT1HMrd95SXjEEfNjZmLBbXWO5tHSA3JbwwM6gEd9mUe8CFW5DqvK1yaqlhPhiUxbeZjIpEyc7DR89EBrRrTxt3VYQggbkqT7JlKpi+oqb9qxH878wF9X/7Lqej6u6TjGNB6Dq538nxBlRFHM84ZH5ibgkcfMizGl4LYOHubku1YnqNXB/NyhRjkHXDFV5TqsKl+bqDp+PhTG6+tPYcw2Uc/LiYWPdqCRr4utwxJC2Jgk3TeRSl0Ic9fz1edW88v5X0g0JALgqHXk/kb3M6H5BAKcA2wboKgeTCaIv5ivRfywebC2HEPBbT0bQkBHqJW7+LYEja78Y7axqlyHVeVrE5WfMdvE+3+cZsXeKwAMbObDZ2Pb4mpf/T6HhBAFSdJ9E6nUhbghMzuTzaGbWXF6BRcTLwKgVqkZXGcwE1tMpKVXSxtHKKqdbCNEn4TwwxBxCMIPQcKlgttp7c1TlQV0NLeG1+oEboFVvlt6Va7DqvK1icotJiWTp384wsHL5vFRZgxsxHP9G6FWV+3PGyFE8UnSfROp1IUoSFEU9kTuYVnwMvZd22cp7+DbgUktJtG7Vm/UKrUNIxTVWnqCuRU8/NCNRLyw+8OdfMyt4AEdzI/+7cG+an3OV+U6rCpfm6i8TkUkMWX5IaKSM3HRa5k/ri0DmvnaOiwhRAUjSfdNpFIX4tbOJZxjefByNoduJlvJBqCua10ea/EYI+qPwF5rb+MIRbWnKJAQAuEHbyTiUSfBlH3ThirwbgqBnSGwi3nxbFCpW8Orch1Wla9NVE5/BkcxY9UxMrJyaOjjzMJHO1BfpgMTQhSiUifdP/zwA1lZWXh7ezNs2DCrdcePHyc0NJTu3bvj4+NT7GNKpS5E8USlRfHj2R/55dwvpGSZB7vysPdgfNPxjG86Hje9m40jFCKfrEyIOmGdiCdeLbidoyfU6gy1c5Nw/3agcyj/eO9QVa7DqvK1icpFURQWBoXw0ZazKAr0buzNlw+3k/u3hRBFqtRJ97Rp04iMjCQmJoZ9+250ef3ss8+YP38+7du35+DBg/zzzz80bty4WMeUSl2IkknLSuPXC7+y8vRKrqVdA8yDro1tMpZHmz+Kt6O3jSMUogipMeYkPOyAeYk8AtmZ1tuoteZ7wwO73GgRd624U/9U5TqsKl+bqDyM2SbeWH+SNYfCAXisWx3eGt4crUZusRJCFK1SJ90Ax44dY9q0aZak22g04u/vz9GjRwkMDGTu3LlcvHiRb7/9tljHk0pdiDuTbcpm6+WtLD61mPPXzwNgp7ZjVMNRTGo5iUCXQBtHKMRtZBvNreFh+83L1f2QGlVwO7dA6y7pvi1BUzHmsq/KdVhVvjZROSSmG5n2/WH2hSSgVsFbw5szqUc9W4clhKgEiluHlfjbREJCAsuWLePIkSM89dRTdO/evcA2586dY9myZcTGxtKuXTumTJmCXq8HIDU1lV9++aXQY48ZMwZn58LvmQkLC8PT05PAQPMX/AEDBrBu3bqShi+EKCGtWsu99e/lnnr3sCtiF4tOLOJY7DHWnF/D2gtrGVpvKJNbTqZRjUa2DlWIwmntbkw51u1p873hSWHm5DsvEY8+ZS5LCoNTa8376RzNg7PlJeG1OoKjh22vRQhRqkLj0nhi2UFC49JwstPw5cPt6de0+LcvCiFEcZQo6V61ahUvvfQSDz74ID/88ANDhw4tkHQfOHCAvn378tBDD9G2bVu++uor1qxZw99//41GoyEzM5MdO3YUevzhw4cXmXQbjUZ0uhv31Oh0OoxGY0nCF0LcBZVKRe9avekV0IvD0Yf57tR37I7YzR8hf/BHyB/0DezLlFZTaOPdxtahCnFrKhW41zYvrR80lxlSzSOlhx0wJ+HhByAzCS7vMi95vJtB7a5Qp7v50b22ba5BCHHXjoUl8sSygySkGQlwd2DxpI409ZPeFkKI0leipLtbt25cvHgRBwcHPv/880K3mTlzJkOHDmXZsmUAPPDAA9SrV481a9Ywfvx4vLy8LOtKolatWoSFhZGeno6joyOnTp2iXj3p+iNEeVOpVHT060hHv46cjj/N4pOL2XZlGzvCdrAjbAed/Tozrc00Ovl1snWoQhSf3hnq9zEvACYTxJ3LbQnPTcTjL0LsGfNyeKl5O9da5uQ7LxH3bgZquQdUiIrun3MxPPX9ETKycmgV4MbiSR3xcZFZOoQQZaNESXedOnVuuT4pKYl///2Xn376yVIWGBhIjx492LhxI+PHjy/WeX7//XeOHTtGbGwsy5Yto0ePHjRq1IiRI0fy4IMPMmjQID777DMWL15c5DEMBgMGg8HyOjk5uVjnFkIUX3PP5nza91NCk0JZemopGy5t4EDUAQ5EHaCDbwemt5lOZ7/OqCrxVE2imlKrwaeZeekwyVyWFgdX98HVvebl2nFIDodTv5gXAHs3CMxNwmt3g4D2oNXb7DKEEAX9cjicV9aeIMek0KuRF99M6ICTvmKM3yCEqJpK9ef40NBQFEUpkJzXqVOHkJCQYh/n0KFDhISE0KtXL3bs2MG1a+aRkxcvXsygQYMIDQ3lu+++Y9CgQUUeY86cObi5uVmWvHvBhRClr55bPd7t8S6bH9jM+Kbj0al1HI4+zJStU5i0ZRL7ru2jgo7ZKETxOXlBs+Ew5AP4v7/h1avw2O/Q9zWo3xd0TuYu6Rf+hL/egaVD4cAiW0ddbiIjI/nkk0/YsGGDrUMRolCKorBgx0Vm/nycHJPC6HYBLJ7YSRJuIUSZu+PRy1UqFStXrmTChAmWsr1799K9e3dOnjxJy5YtLeVPP/00u3fv5tixY3cdcHEV1tIdGBgoo6MKUQ6i0qJYcmoJa8+vxWgyj73Qzqcd09tMp2vNrtLyLaqmnGzzKOn5W8MfXm0ejO0uVfQRvk0mE/379+fw4cMMGTKkyAFTC1PRr01UDTkmhfc2nmbZnssAPNmnPq8MaYpaLfWREOLOldno5bfi5uYGwPXr163KExISLOvKi16vt4yYLoQoX35OfrzW5TUmt5zMklNL+OX8LxyNOcrUbVNp692W6W2m082/myTfomrRaM3dyQPaQ7enzKOkVxPvvfcenp6e9OnTx9ahCFGAITuHF1cf54+T5p6Tbw5vzuSeMi6QEKL8lGr38oYNG+Lg4MCpU6esyk+ePEmrVq1K81RCiErA18mXWV1msfmBzTzS7BHs1HYciz3Gk9ufZMLmCeyJ2CPdzkXVpVKZlypu165dLFmyhEWLqk9XelF5ZBhzmLL8EH+cvIadRs0X49tJwi2EKHelmnTb2dnxwAMPsHDhQjIyMgDYtm0bwcHBPPzww6V5KiFEJeLj6MOrnV9l8wObmdBsAnqNnhOxJ3hy+5M8/ufjHIk+YusQhRB3ICEhgQkTJrB48WI8PIo3h7nBYCA5OdlqEaIsJGdm8diS/ey6EIeDTsOSSZ0Y0cbf1mEJIaqhEt3THRwczJw5cwD44Ycf6N69O/Xq1aNfv35MnjwZgJiYGAYOHEhqaipNmjQhKCiIl156iXfffbdsrqCY5J4xISqO2PRYlpxawppzayz3fPcM6Mmz7Z6luWdzG0cnRMVTXnXY7t272bt37y23GTt2rGVw0vvuu4+GDRvy6aefAjB8+HDs7e1veU/37NmzeeeddwqUS/0sStP1NCOPLTnAyYgkXOy1LHu8Ex3qFO+HISGEKK7i1s8lSrqjoqLYvn17gfKGDRvStWtXy+vs7Gx27dpFfHw8bdq0oVGjRiUMv/RJ0i1ExROVFsW3J75l3YV15Cg5AAyqM4in2z5NA/cGNo5OiIqjvOqwLVu2FFrP5zd9+nQaNGjAtm3bGD58OG+99ZZlDJUlS5ag1Wp57LHHmDp1aqGxykCnoqzFJGcyYfF+zken4uFkx4onOtMyoHzHFhJCVA9lknRXZpJ0C1FxXU2+yoLjC9gUsgkFBbVKzfD6w5nWZhqBLjLdnxAVsQ4LDg5m6dKlVmXr1q1Dq9UyYsQIXn/9dWrUqHHb41TEaxOVV/j1dB75bj9X4tPxddXzw5QuNPRxsXVYQogqSpLum0ilLkTFd+H6Bb48+iV/h/0NgFal5YHGDzC19VR8HH1sHJ0QtlNZ6rDidC+/WWW5NlHxhcSm8sh3+7mWlEmghwM/TO5KbU9HW4clhKjCiluHlepAakIIcTca1WjE5/0/56dhP9HdvzvZSjarz63m3l/vZd7heSQbZcAlIYQQBV2MSeWhb/dxLSmTBt5O/Pxkd0m4hRAVhiTdQogKp6VXS74d9C1LhyylvU97DDkGlpxawr2/3svy4OUYc4y2DlEIUYj777+fkSNH2joMUc1cjElh3MJ9xKUaaOrnwponu+HnZm/rsIQQwkK6lwshKjRFUQgKD2L+kflcTLwIgL+TP8+0e4Zh9YehVslvh6Lqq8p1WFW+NlH2zAn3fuJSDTSr6coPU7rg4WRn67CEENWEdC8XQlQJKpWKPoF9+GXEL7zb/V18HHyITIvktX9fY9zGceyJ3GPrEIUQQtjAhWjrhPtHSbiFEBWUJN1CiEpBo9YwutFoNt6/kefbP4+zzpkzCWd4ctuTTN06lTPxZ2wdohBCiHJyITqF8YvMXcqb5ybcNSThFkJUUJJ0CyEqFQetA1NaTWHT/ZuY0GwCWrWWvdf28tDGh3h116tEpEbYOkQhhBBl6EbCbaR5bpdySbiFEBWZJN1CiEqphn0NXun8Cr+P+p176t0DwB8hfzBy3UjmH55PqjHVxhEKIYQobRdjbiTcLfwl4RZCVA6SdAshKrVAl0D+2/u/rBq+ii5+XTCajCw+tZhh64bxy/lfyDHl2DpEIYQQpeByXBoPL9ovCbcQotKRpFsIUSW08GzBosGL+F+//1HHtQ4JmQm8s/cdHtr4EPuv7bd1eEIIIe5CRGIGj3y3n5gUA018Xfh+chfcHSXhFkJUDpJ0CyGqDJVKRb/a/Vg3ch0vd3oZFzsXzl8/z5StU3j272e5nHTZ1iEKIYQooZjkTB5ZtI+IxAzqeznxvbRwCyEqGa2tAyjM+fPnMZlMODg4UKdOHUv5lStXyMjIsLxu1KgRGo3GFiEKISownUbHo80fZUT9ESw4voA159awI2wH/4b/y7im45jWZhpuejdbhymEEOI24lMNPPLdfi7Hp1OrhgM//F8XvF30tg5LCCFKRKUoimLrIG7WoUMHEhIS8PX1Zd++fZbyrl27Ehsbi06nA+Dff//Fy8urWMcs7sTlQoiqJyQxhE8OfcKuiF0AuOndmN5mOg81eQidWmfj6IS4vapch1XlaxN3Jykji4cX7SM4Mhk/V3vWPNmN2p6Otg5LCCEsiluHVcju5YcPH2bdunWFrlu9ejW7d+/m7NmzxU64hRDVW333+iwYuIBvBn5DA7cGJBmS+OjARzy04SEOXDtg6/CEEELcJNWQzaSlBwiOTMbL2Y4f/q+LJNxCiErrjrqX79mzhyNHjjB48GAaN25cYH1GRgZbt24lNjaWdu3a0aFDB8u6rKwsLl26VOhxGzRoYGnFLkzdunWZMGEC165do2vXrqxbtw5HR/kAFkIUT4+AHnSp2YW159fy5bEvuZh4kclbJzO4zmBmdpxJTeeatg5RCCGqvcysHKYsP8jRq4m4OehYObkLDbydbR2WEELcsRIl3UFBQTz77LM4ODiwf/9+Vq5cWSDpDg8Pp2/fvjg6OtK8eXNmzpzJxIkT+fzzzwGIiYlh1KhRhR7/r7/+IiAgoMjzr1q1CgCDwcDYsWNZsGABM2fOLMklCCGqOa1ay9imYxlabyhfHP2Cn8//zNYrWwkKD2JKqylMajkJvUbuFxRCCFvIzjHx3E9H2ReSgLNey4onOtOsptx2IISo3EqUdGu1WlasWEGbNm1QqVSFbvPiiy/i6enJrl27sLOzY/fu3fTs2ZNRo0bRr18/AgICOHv27F0FrdfrGTx4MGfOnLmr4wghqi83vRtvdH2DBxs/yIf7P+RIzBG+PPYl6y+u5+VOL9M3sG+Rn3NCCCFKn6IovL7uFFtPR2OnVbPosY60CXS3dVhCCHHXSnRPd/fu3WnTpk2R6zMzM/ntt9+YMmUKdnbmqRx69OhB27ZtWb16dbHPc/nyZS5fvkxmZiZnz54lKSmJnJwczp49y9mzZ9myZQvz5s2jb9++RR7DYDCQnJxstQghxM2aeDRh2dBlfNTrI3wcfAhPDee5f55j+l/TZYoxIYQoR//98xyrD4WhVsEX49vRrYGnrUMSQohSUapThl26dAmj0UizZs2syps1a8bp06eLfZwXXnjB0oo9atQo3n//ffr27cuoUaNQqVT4+vryn//8hwceeKDIY8yZM4d33nnnzi5ECFGtqFQqhtUfRt/Aviw8sZAVp1ewO2I3o6+N5rHmj/Fk6ydx1Mn4EUIIUVa+2xXC1zvMY/7Mub8VQ1r42TgiIYQoPaWadKekpADg7u5uVV6jRo0SJd1FjVxekm7ps2bN4sUXX7S8Tk5OJjAwsNj7CyGqHyedEy90eIHRDUfz8cGP+TfiX5acWsLGSxt5pfMrDKozSLqcCyFEKVt7OJz3/zA3trw8tAljO9W2cURCCFG6SnXKMAcHB+BG8p0nOTm53EcZ1+v1uLq6Wi1CCFEcdd3qsmDAAr7o/wW1nGsRkxHDSztfYvpf0wlLDrN1eEIIUWX8dSaal9eeAGBKz3pM79PAxhEJIUTpK9Wku379+qjVakJDQ63KQ0JCaNSoUWmeSgghypRKpaJvYF/W3beOaW2moVPrzF3Ofx/NN8e/wZhjtHWIQghRqR28nMBTPxwhx6Rwf7sAXru3mfQmEkJUSaWadLu4uNCvXz9++OEHS9nFixfZt28fI0eOLM1TCSFEubDX2vN026f5deSvdKnZBUOOga+OfcUDvz/A/mv7bR2eEEJUSmejkpm87CCGbBP9m/rw8ZjWqNWScAshqqYS3dMdHh7O+vXrLa+3bdtGYmIiLVq0oF+/fgB89tln9OrVi1GjRtGuXTuWLVvGkCFDuP/++0s1cCGEKE913eqyaNAiNoduZu6huVxOvsyUrVO4t969/KfTf/By8LJ1iEIIUSlEJGYwcckBkjOz6VinBl893B6dplTbgYQQokIp0SdcWlqaZdqup59+GhcXF86ePcu1a9cs27Ru3ZpTp07RpUsXUlJSeP/999mwYYN0FxJCVHoqlYp769/L76N+Z3zT8ahQsSl0EyPXjeSnsz+RY8qxdYhCCFGhJWVk8fjSA0QnG2js68ziiZ1wsNPYOiwhhChTKkVRFFsHUR6Sk5Nxc3MjKSlJBlUTQpSK4Lhg3t33LqfjzbMztPBswVvd3qK5Z3MbRyaqmqpch1XlaxPWDNk5TFxygH0hCfi66ln3VA/83R1sHZYQQtyx4tZh0pdHCCHuUAuvFvx474+81uU1nHXOBMcH8/AfD/PpoU/JyM6wdXhCCFFhmEwKL/9ygn0hCTjrtSyd1FkSbiFEtSFJtxBC3AWNWsP4puPZMHoDQ+sOJUfJYVnwMkb/Npo9kXtsHZ4QQlQIc7ee47djkWjVKr6e0J7m/tKrQQhRfUjSLYQQpcDLwYu5feby1YCv8HPyIyI1gie3Pcnr/75OYmaircMTokwdPnyY999/nw8++ICzZ8/aOhxRwXy/7wpf77gEwJz7W9GrkbeNIxJCiPIlSbcQQpSi3rV6s/6+9TzS7BFUqPj90u+MXD+SjSEbqSZDaIhq5j//+Q/9+vUjLi4OgAkTJhAUFGTjqERFsf10NG/9dgqAFwY25sGOgTaOSAghyp8MpCaEEGXkeOxxZu+ZzcXEiwD0COjBm13fJMA5wMaRicqmotZha9eu5aGHHmLPnj106dIFAKPRSExMDLVq1SrWMSrqtYm7dywskfEL95GRlcPYjoF89EArmc1GCFGlyEBqQghhY22827Bm+Bqebfcsdmo7dkfsZvRvo1kevJxsU7atwxPiri1YsIChQ4daEm4AOzu7Yifcouq6Ep/G5GUHycjKoU9jb94f3VISbiFEtSVJtxBClCGdRsfU1lP5ZeQvdPDtQEZ2Bp8c+oQJmyZw4foFW4cnxF05dOgQPXv25M8//2TWrFnMmzePixcv3nIfg8FAcnKy1SKqlqT0LB5fdpD4NCMt/F356pH26DTylVMIUX1J93IhhCgnJsXErxd+5bNDn5GSlYJWreXJ1k8yudVkdGqdrcMTFVh51WGbNm1i69att9zmmWeeoWHDhiiKglqtpmXLltSoUYMhQ4Zw+vRp1q5dy+rVq7nvvvsK3X/27Nm88847Bcqlfq4ajNkmJi45wN6QePzd7Fn/dA98XO1tHZYQQpSJ4tbPknQLIUQ5i0mP4b1977EjbAcATT2a8l6P92jq0dSmcYmKq7zqsL1797J///5bbjNmzBhL93EXFxcCAwM5efIkGo0GgCeffJIdO3Zw7ty5Qvc3GAwYDAbL6+TkZAIDA6V+rgIUReHVtSdZfSgMJzsNv0zvTrOa8jcVQlRdxa2fteUYkxBCCMDH0Yf/9fsfm0I3MefAHM4mnGX8xvFMaT2Fqa2motNIq7ewjW7dutGtW7dib9+iRQvq169vSbgBOnTowNKlS1EUpdB7ePV6PXq9vlTiFRXLwqAQVh8KQ62CLx9uLwm3EELkkhtshBDCBlQqFcPqD2P9fesZWHsg2Uo23xz/hrF/jCU4PtjW4QlRLOPHj2fPnj2kp6dbyrZv307r1q1l0KxqZsupKD7aYp6j/a3hzenX1MfGEQkhRMUhSbcQQtiQl4MXn/X9jLl95lJDX4ML1y/wyB+P8PmRzzHmGG0dnhC39NRTT9GiRQtat27N1KlT6dGjB0FBQXzzzTe2Dk2Uo5PhScxYfRRFgce61WFSj3q2DkkIISqUCnlP9+HDh1myZAlqtZrHH3+c9u3bA5CZmcmnn35KaGgoI0aMKHKQlsLIPd1CiIouITOBOfvnsOXyFgAauDXgvR7v0cq7lY0jE7ZWkeswRVEs93DXrFmT/v374+LiUuz9K/K1idu7lpTBfV/uJibFQJ/G3iye2BGtjFQuhKgmKu083YcPH+bZZ5+lRYsWeHl50a9fPyIjIwGYOHEiR48epWvXrrz44ots2rTJxtEKIUTp8bD3YG6fuczvOx9Pe08uJV1iwuYJzDs8T1q9RYWlUqno168f06ZN47777itRwi0qtzRDNpOXHSImxUBjX2e+eLidJNxCCFGICjeQWoMGDQgKCkKrNYe2efNmIiIicHZ25s8//yQqKgp7e3tcXV1ZvHgx9957r40jFkKI0jWgzgA6+Hbg44MfszFkI0tOLWFXxC4+7PmhjHAuhKgQckwKz686yulryXg527F4Yidc7WUQSCGEKEyJfo7Mysri559/pn///ri7u7NmzZpCt1uwYAGtWrXCz8+Pe+65h+DgG4MCxcbGMmbMmEKX2NhY3N3dLQn3zp07cXJyokOHDoSHhxMYGIi9vXmux+bNm3P58uU7vGwhhKjY3O3dmdNrDp/3+xwPew8uXL/A+D/Gs/DEQrJN2bYOTwhRzc3ZdIbtZ2Kw06pZ+FhHAj0cbR2SEEJUWCVq6f7yyy/ZvXs3r732GoMGDcJoLNjdcdGiRcycOZPly5fTrl073nvvPfr378+ZM2fw8PDAycmJcePGFXp8Jycny/PNmzczb948fv31V9RqNfb29lbzehoMBksCLoQQVVX/2v1p492G9/a9x19X/+KLo1+wM2wnH/T8gLpudW0dnhCiGvpx/1W++zcUgE8fbEP72jVsHJEQQlRsJRpILf+cmyqVipUrVzJhwgSrberXr8/999/PJ598AoDRaMTHx4fXXnuNl19+uVjnWbx4MevXr2f16tU4Opp/Oc3KysLf3599+/bRoEEDPvjgAyIjI/nqq6+KdUwZqEUIUZkpisLGkI3M2T+HlKwU7DX2zOgwg/FNx6NWyT2UVV1VrsOq8rVVRXsuxfHo4gPkmBReGtSYZwc0snVIQghhM2UykNrt5tyMiIggNDSUAQMGWMrs7Ozo3bs3//77b7HOsWPHDqZOnYpKpeKxxx5jzJgxHDx4EJ1Ox+zZs+nRowcDBw7km2++4cUXXyzyOAaDgeTkZKtFCCEqK5VKxYgGI/j1vl/pWrMrmTmZfHTgI6Zuncq11Gu2Dk8IUQ1ciU/jqR+OkGNSuK+tP8/0b2jrkIQQolIo1YHUIiIiAPD19bUq9/Hx4dixY8U6RuPGjVm9erVVWUBAAABPP/00gwcP5vLly3Tq1Al3d/cijzNnzhzeeeed4gcvhBCVgJ+TH98O+pbV51bz2aHP2B+1n/t/v59XO7/KyAYjb/vjqBBC3ImUzCymLD9EYnoWbWq58fEDreXzRgghiqlMRi9Xq60b0LVaLcXtxe7v78+YMWOKXN+oUSMaNbp9V6ZZs2ZZtYQnJycTGBhYrBiEEKIiU6vUjG86nu7+3Xn939c5HnucN3a/wV9X/2J299l42HvYOkQhRBWSY1KYseoYF2JS8XXVs/CxjtjrNLYOSwghKo1SvRHQx8cHgLi4OKvymJgYy7ryotfrcXV1tVqEEKIqqeNah+VDl/N8++fRqrX8E/YPD/z+AP9GFO92HiGEKI5Ptp7jr7Mx6LVqFj7aEV9XGchWCCFKolST7jp16uDn58euXbssZTk5OezevZsuXbqU5qmEEEIAGrWGKa2msGrYKhq6NyQuI47p26czZ/8cMrMzbR2eEKKSW380gq93XALgv2Na0ybQ3bYBCSFEJVSqSbdKpeLZZ5/lq6++4tChQxgMBt5++23S0tKYMmVKaZ5KCCFEPk08mvDTsJ94uOnDAPx49kfG/zGecwnnbByZEKKyOhaWyMtrTwAwvW8D7msbYOOIhBCicipR0r1z507c3d0tA5g9+eSTuLu789xzz1m2efXVV3niiSfo168fTk5OrF27lg0bNlCrVq1SDVwIIYQ1e609s7rM4uuBX+Np78nFxIuM/2M8y4OXY1JMtg5PCFGJRCVlMnXFIYzZJgY28+E/g5vYOiQhhKi0SjRPd3Z2NqmpqQXK9Xo9Dg4OVmWKomAwGLC3rxj3/cg8oEKI6iQhM4G397zNjrAdAHSt2ZX3e7yPr5PvLfcTFVNVrsOq8rVVVplZOYz9di/Hw5No7OvMr0/1wFlfJmPvCiFEpVYm83RrtVpLS3f+5eaEG8xdzStKwi2EENWNh70H/+v3P97s+ib2Gnv2XdvHAxseYPuV7bYOTQhRgSmKwitrT3A8PAl3Rx3fPdZJEm4hhLhLpXpPtxBCiIpDpVLxUJOHWDNiDc09m5NkSOKFHS/w1u63SM9Kt3V4QogK6Oudl/jtWCRatYoFj7SntqejrUMSQohKT5JuIYSo4uq51eP7e75nSqspqFCx7uI6Htr4EKfjT9s6NCFEBbL9dDRz/zQPvvj2yBZ0b+Bl44iEEKJqkKRbCCGqAZ1Gx/Ptn2fJkCX4OflxJfkKj2x6hJWnV1KCoT2EEFXUuagUnl91FEWBCV1r82jXOrYOSQghqgxJuoUQohrp6NeRX0b8woDaA8g2ZfPfg//lmb+fISEzwdahCSFs5HqakSkrDpJmzKFrfQ/eHtHC1iEJIUSVIkm3EEJUM256N+b1nccbXd7ATm1HUHgQY34fw4FrB2wdmhCinGXlmJj+w2HCEjII9HBgwSMd0Gnk66EQQpQm+VQVQohqSKVSMbbpWH4c9iP13eoTmxHLlK1T+OLoF2Sbsm0dnhCinLyzIZh9IQk42WlYPLETHk52tg5JCCGqHEm6hRCiGmvi0YSfhv3E/Y3uR0Fh4YmFPPHnE1xLvWbr0IQQZWzlvit8v+8qKhV8Pq4djX1dbB2SEEJUSZJ0CyFENeeoc+Sd7u8wt/dcnHXOHI05KnN6C1HF7bkUxzu/BwPwnyFNGNjc18YRCSFE1SVJtxBCCACG1hvKmhFraOXVihRjCi/seIH3972PIcdg69CEEKXoanw6T/1whGyTwn1t/Znep4GtQxJCiCpNkm4hhBAWgS6BLL9nOY+3fByA1edWM2HTBK4mX7VxZEKI0pCSmcWUFQdJTM+iTS03Pn6gNSqVytZhCSFElaa1dQAVSU5ODllZWbYOo8rTaDRotVqp5IWooHRqHS92eJGufl2Z9e8sziac5aGND/Fu93cZXHewrcMTQtwhk0nhhdXHOB+dio+Lnm8f7Yi9TmPrsIQQd0BRFLKzs8nJybF1KFVaaeUtFTLpzs7O5uDBg6jVajp16oRabW6QDwoKIiHhxlyy99xzD3q9vlTOmZqaSnh4OIqilMrxxK05OjpSs2ZN7OxklFQhKqruAd1ZM3wNLwe9zJGYI7y08yXGR49nZseZ2Gnk/64w279/P2vXriUmJgZ/f3/GjRtH69atbR2WKMTcrefYfiYGO62ahY91xM/N3tYhCSHugNFo5Nq1a6Snp9s6lGqhNPKWCpd0R0REMHr0aLRaLQkJCbi5ubFz507s7e15+eWXcXV1xdHREYC+ffuWStKdk5NDeHg4jo6OeHt7SwtsGVIUBaPRSGxsLKGhoTRq1Mjyo4oQouLxdfJl8ZDFfHH0C5acWsJPZ3/iROwJPunzCbVcatk6PGFjK1eu5PHHH+eFF16gb9++7Nu3j/bt27Nx40aGDh1q6/BEPuuPRvD1jksA/PeB1rQNdLdtQEKIO2IymQgNDUWj0eDv74+dnZ3kLmWkNPOWCpd0p6SksGjRItq0aYPJZKJVq1YcPXqUbt26ATB58mRq1apFx44dS62VOysrC0VR8Pb2xsHBoVSOKYrm4OCATqfjypUrGI1G7O3ll3YhKjKtWssLHV6gg28HXvv3NYLjg3lo40O83+N9+tfub+vwhA0tWbKEsWPHMnfuXAAmTZrE+fPnWb58uSTdFcixsEReXnsCgGl9GjCqXYCNIxJC3Cmj0YjJZCIwMNDSECnKTmnlLSVOutPT01m1ahVHjhxh4sSJdOrUqcA2YWFh/Pjjj8TGxtKuXTvGjRuHRqOx7L9169ZCjz148GCaNm1qeZ2cnEx2djbNmjUDoE+fPvz000+EhoaSkZHB7t278fb2LuklFEl+JSo/0rotROXTu1Zvfh7+MzODZnIi9gTP//M8jzV/jBkdZqBT62wdnrCB+vXrc+HCBRRFQaVSWbo89ujRw9ahiVxRSZlMXXEIY7aJgc18+M+QJrYOSQhRCuS7dPkpjfe6REdYv349DRs2JCgoiK+++opz584V2ObEiRO0bNmSvXv34uzszBtvvMHIkSMt90qnpqaybNmyQpfU1FTLcRISEhgzZgwLFy7E3d0dgI8//pj169dz/Phx+vXrx4IFC+7i0oUQQpRUTeeaLBuyjEebPwrAitMreHzL40SlRdk4MmEL8+bNo2HDhjRv3pzRo0fTrFkz7r33Xt54440i9zEYDCQnJ1stomxkZuUwdeUhYlIMNPZ1Zv64dmjU0sAghBDlrUQt3c2aNSM4OJgaNWqwfPnyQrd56aWX6NGjB+vXrwdgwoQJNGnShHXr1nH//ffj4+NjWVeUK1eu8Mgjj/DJJ5/QtWvXImO5fPlyScIXxbRr1y42btwIwPDhw+nVq9ctt1+4cCFHjhyxKmvcuDEvvvhimcUohLAdnUbHy51epoNvB978902Oxx5nzIYxzOk5h161bv15ISq2NWvW8Ouvv95ym7feeovmzZsDsHfvXjZs2MCECRNo1aoV3t7e/PTTT4wfP56OHTsWuv+cOXN45513Sj12YU1RFF7+5QQnwpOo4ajju8c64ayvcHcVCiHEXYmOjmb58uVERkbSrFkzJk6ceMsu4KdPn+Z///tfgfL33nuvVHtQ36xEn75Nmty6S1JKSgp///03K1eutJQ1bNiQrl278ttvv3H//fff9hyXLl2id+/eTJ8+naioKNavX0+XLl3w8PBg8+bNAFy+fJkPPviAn3/+ucjjGAwGDAaD5bX8kl48n3/+Oa+++irPPPMMKpWKwYMH89FHH/H8888Xuc/WrVuJiYnh4YcftpT5+/uXR7hCCBsaUHsAjUc0ZubOmZyOP83Tfz3Nk22eZFrraWjUMg1RZdSsWTNGjRp1y228vLwsz6dOncrEiRP55JNPAHjiiSeYMGECzz77LHv37i10/1mzZln9KJucnExgYODdBy+sLNhxid+PR6JVq1jwSAdqe8q9n0KIqiU0NJSuXbvSrl07evfuzZdffsnixYsJCgoqMvG+evUqixYt4quvvrIqL62xwopSqj95hoSEYDKZqFevnlV53j1fxREXF0enTp04dOgQhw4dAqBmzZo4ODiwbNkyVCoVvr6+/PLLL/Tp06fI41SnX9KNRiMhISFER0dbTXnWsmVLqy9Ht5OQkMCsWbOYN28e06ZNA6B27dq8/PLLPProo3h4eBS5b/PmzS37CCGqj0CXQFbes5L/Hvwvq8+t5pvj33Ay7iQf9/oYN72brcMTJdSqVStatWpVrG1zcnKIiIiwjLuSp2nTpvz1119F7qfX68v8y011t+10NJ9sNd8COHtkC7o18LRxREIIcUNkZCSXL1/GaDRayjw9PYtd/+R54403qFu3Ln/88QcajYapU6dSr149vvvuO5555pki91OpVOWet5Rq0p2RkQGAi4uLVbmrq2ux55Hr0qVLkd3Pb9ctPb/q8kv68uXLeemll/D29sZgMBAaGmpZd+TIEbZu3UpQUNAtj/HWW2/h7+/Ptm3byMzMZNy4cZZ1jzzyCM899xzbt2/noYceKvIYx44dY8aMGXh7e9O/f3/LaPNCiKrPTmPHG13foI13G97d+y67I3YzduNYPuv7Gc09m9s6PFFGNBoNbdq04ddff2XSpEloNBoMBgO///47HTp0sHV41dbZqGRmrDqKosCjXeswoWsdW4ckhBCAuUfzxIkTOXv2LLVr1+bUqVNkZWUB8NRTTzFz5kw+/vjjWx5j0KBBPPDAAyiKwoYNG5g9e7ZlwG4vLy+GDBnC77//fsukW1EU3njjDYxGIy1btmTcuHF3NQd3cZRq0u3s7AxAYmKiVfn169dxdXUtzVPd1t38kq4oChlZOaUcUfE46DTFHkV9x44d/N///R+rV69m9OjRKIrCmDFjuHDhAidOmKcGSU1NpW3btrc+Z+40aRcuXKBGjRqWgesAatSogZubGxcvXixyf7Vajb+/PwEBAYSEhNCvXz9mzJjBRx99VKzrEEJUDSMajKBxjca8sOMFwlLCeHTTo7zR9Q1GNxpt69BEGVm0aBH3338/TZs2pXnz5hw9ehQnJyd++uknW4dWLSWkGZmy/BBpxhy61ffkrRHyo5cQ1YGtcpeS5C0Gg4H777+fJk2asH37duzt7dm0aRPDhg1j//79dO7cmZiYmNvmLQEB5ikPo6OjSUlJKdDDul69eqxdu/aWx2jSpAkajQYnJyfef/99PvzwQ3bv3o2nZ9n1CirVpLthw4bY2dlx5swZevbsaSk/c+YMXbp0Kc1TlamMrByav/WnTc59+t0hONoV78/y3nvv8fjjjzN6tPkLrUqlYtSoUUyaNMmyTa9evW47EFqejIyMAr0U4PY9FebNm2f5DwAwdOhQRo0axZgxY4ocSEcIUTU18WjCquGreH3X6+wI38Fbe97ieOxxZnWZhV4jXYqrmvbt23P+/HmOHDlCdHQ0AQEBtGvXztLqIMqPMdvE9O8PE349gzqejix4pD06jUwpJER1YKvcpSR5y6+//sqVK1f4999/Lfdb33PPPej1eoKDg+ncuTM+Pj7F7vZ9pz2sO3TowPHjx9HpzFOdzpw5k5YtW/LGG2/w9ddfF+vcd6JUk257e3vuu+8+Fi9ezMSJE7Gzs2PPnj0cO3aMzz77rDRPVe0lJiayc+dOXn75ZavypKQkq19pfvzxx2J3L3d1dS3QSwFu31Mhf8INMHLkSJycnNi7d68k3UJUQ652rnze/3O+O/kdXx79krUX1nIm4Qzz+s7D31kGWaxq7OzsipxpRJQPRVGYvSGY/aEJOOu1fPdYR2o4lW1XSSGEKIn169czaNAgqyQ5MzMTo9FoyV1CQ0OL3b08LzcpaQ/rm0cod3JyYuTIkezYsaMEV1NyJUq6z507xxdffGF5vWLFCvbt20f37t0tI1d/9tln9OnTh44dO9KiRQs2b97Mc889R79+/Uo38jLkoNNw+t0hNjt3cYSEhJCTk1PgPvU///yTe++91/I6MDCw2N3LW7RoQVJSEpGRkZbRx8PDw0lJSaFFixbFvoacnByysrLIybFNF30hhO2pVWqmtp5KS8+WvLLrFU7Hn+ahjQ/xca+P6RHQw9bhCVGlrNx3hR/3X0Wlgv+Nb0sj34K91oQQVZetcpfi5i1gvo21b9++VmVbt25Fr9db8kQnJ6didy/39PTE19eX06dPW60/ffp0ifIWMCf/ZZ23lCjpdnJyomnTpgBWyXfNmjUtz2vVqsWpU6fYsmUL8fHxzJw5s9INqKJSqYrdVcJWatSoAcCJEycs86X++eefbN26lT179li2K0n38gEDBuDt7c1XX33FBx98AJj/znmDo4H516RXX32VJ598knbt2pGUlMTJkyetbieYN28e2dnZDB48uFSuVQhReXUP6M7q4at5cceLBMcHM337dJ5q+xRTW09FrZKur0Lcrd0X43hng/lL5ytDm9K/qa+NIxJClLfKkrucPHnS8jolJYXXXnuNqVOnWlq/S9K9HGDcuHGsXLmSF154ARcXF06fPl1g+uq8nPXZZ58FYNu2bfTt29fSvTwkJIQ1a9YwefLku77GWynRX6dWrVq3HAkuj4ODg+U+Y1E26tWrxz333MP06dM5dOgQCQkJ/PTTT3z55Zd3/COHg4MDS5YsYdy4cRw8eBCAPXv2sGrVKktreGpqKt9++y0DBw603Lc3e/ZsUlJSaNKkCRcuXOD06dN8++23lh8DhBDVm7+zP8vvWc5HBz7il/O/8NWxrwiOD2ZOzzk42znbOjwhKq2Q2FSe+uEIOSaF+9sF8GTv+rYOSQghCvX000/zwAMPMHr0aBo3bswvv/xC7dq1+e9//3vHx5w9eza7du2ibdu2dOzYkb/++ouxY8dazbj0xx9/ADeS7tOnT/P000/TunVrTCYTW7duZfDgwcyePfuuru92VEr+iZ2rsOTkZNzc3EhKSirQzz8zM5PQ0FDq1atX5ETqFVFmZiaLFy8mODgYHx8fHnrooVJJdKOiovj7778B6N+/P35+fpZ1aWlprFy5ksGDB1O//o3K/dixY5w6dQovLy86d+58yzm982KvjO+5EOLurLuwjvf3vY/RZKSeWz3+1+9/1HWra+uwKrxb1WGVXVW+trKUmG5k9II9hMal0TbQnVVTu2Jfgq6eQojKqTJ/h96xYwe///47JpOJXr16MXr0aNTqu+v1lp2dzV9//UVkZCTNmjUrMMbIpk2bAKxuv42JiWH37t1kZ2fTokWL2+ZPt3rPi1uHSdJN5f7HW1nJey5E9XUy9iQzdswgJj0GZ50zH/f+mN61ets6rAqtKiemVfnayoox28Sji/ezPzSBAHcH1j/dA28XmR1AiOpAvkOXv9JIuuWGOiGEEOWqlXcrVg9fTTufdqRmpfLMX8+w6MQiqslvwELcFUVReG3dSctI5UsmdZKEWwghKjhJuoUQQpQ7LwcvFg9ezEONH0JB4X9H/8dLO18iPavouTWFEPD1zkv8cjgctQq+fLgdTfxkpHIhhKjoJOkWQghhEzqNjje7vcnb3d5Gq9ay7co2Htn0CGHJYbYOTYgKadPJa/x3yzkAZo9sQd8mPjaOSAghRHFI0i2EEMKmxjQew9IhS/Fy8OJi4kXG/TGOPRF7br+jENXI8bBEXlh9DIBJ3evyWLe6No1HCCFE8UnSLYQQwuba+rRl9fDVtPZqTbIxmel/TWfpqaVyn7cQQERiBlNWHMKQbaJfE2/eHC5TcgohRGUiSbcQQogKwcfRh6VDlzK64WhMionPDn/GK0GvkJGdYevQhLCZVEM2k5cdJDbFQFM/F754uD0atcrWYQkhhCgBSbqFEEJUGHYaO97p/g6vd3kdrUrL5subeWzzY0SmRto6NCHKXXaOied+OsrZqBS8nPUsntQJZ73W1mEJIYQoIUm6hRBCVCgqlYpxTcexaPAiPOw9OJtwlvF/jOdYzDFbhyZEuVEUhTd/O8XfZ2PQa9V8N7EjAe4Otg5LCCHEHZCkWwghRIXU0a8jq4atokmNJiRkJvDEn0/w+6XfbR2WEOXiy78v8tOBMNQq+N/4drQNdLd1SEIIIe6QJN2igJycHP755x+2bNlS7H1MJhMHDx7kzz//JCoqqgyjE0JUJzWda7LinhUMqD2ALFMWr//7Op8d/owcU46tQxOizPx8KIxPt50HzFODDWnhZ+OIhBCi4jp//jzr168nOjq62PuEhYWxZcsWjh49WoaR3VAhk+4ffviBZs2aUa9ePd59911L+ZkzZ+jVqxe1atXiqaeewmg02jDKqum///0vDRs25LHHHmPSpEnF2ic2NpZOnToxatQo3n33XerXr8/nn39etoEKIaoNR50jn/X9jP9r9X8ALD21lBn/zCAtK83GkQlR+oLOxzLr15MATOvTQKYGE0KIIuzdu5dBgwYxZMgQRo8ezcGDB4u13+zZs2nSpAkffvghQ4YMoU+fPqSkpJRprBUu6U5KSmLHjh38/vvvbN68me+//56goCAAJkyYwLhx49i3bx+hoaF88803No62atqxYwfPP/98sbefMWMGJpOJCxcusHv3bpYvX86LL77IiRMnyjBKIUR1olapea79c3zU6yPs1HbsCN/BhE0TiEiNsHVoQpSaUxFJTP/+MNkmhVFt/Xl5SBNbhySEEBVWbGwsL7/8MkeOHCn2Pjt27ODdd99l8+bNBAUFcfbsWcLCwnjrrbfKMFK4oyEwExMTOXv2LA0bNsTLy6vQbcLDw4mLi6Nx48Y4OjoW+9hubm4sWrQIMCfgzs7OODs7Ex0dTUREBE8//TQAM2fO5OOPP+a55567k0uoMs6ePcubb77JwYMHiY6OtprTdtu2bbi6uhIaGnrLYwwYMAAXFxcAXn755RKdPz09nbVr1/LFF19Y/s4PPvggL7/8Mj/88AOtW7cu4RUJIUTRhtUfRqBLIM//8zwXEy8yfuN45vebT3vf9rYOTYi7EpaQzuPLDpJmzKF7A0/+O6YNapkaTAhRhWRmZjJ37lx+/vlnLl++bNVr+ZFHHuGTTz5h586dtzxGo0aNaNGiBQAjR44EzLlpca1YsYJOnTrRp08fADw8PHjiiSeYP38+n332GSpV2XzulijpvnDhAh999BGbNm0iKiqKlStXMmHCBKtt0tPTGTduHNu3b8ff35/o6Gi++OILS1flq1ev0rx580KPf/r0aWrXro3BYMDT05P09HSmTp1K+/btOXnyJL6+vpZta9asSUxMTAkvt2oJDw+nZ8+ejBkzhp07d5KZmcmDDz6I0Whk/vz5dOrUiWXLlt323uz27dtbku6SOnfuHAaDgTZt2liVt2nTRlq6hRBlorV3a34a9hPP/f0cZxLOMHnrZN7q+hajG422dWhC3JHEdCOTlh6wzMX9zaMdsNNWuM6IQghxV6ZMmcLhw4dZuHAhTZo04fPPP+fDDz/k66+/ZuDAgcTGxrJs2bJbHmPUqFGWpPtOHD9+nA4dOliVtW3blvj4eCIjIwkICLjjY99KiZLukydP0q1bN7744gucnJwK3ea1117j1KlTXL58GR8fH5YvX87kyZPp1KkTLVq0IDAwsMiBtvKOqdfriYqKIjo6mv/7v/9j1apV9O7dm/j4eMu2cXFxeHp6liT84lMUyEovm2Pfjs4RivkLy6xZs2jdurVVN/unn36a5557jqFDhwIwbdo0pk2bViahgrk3AkCNGjWsyj09PTl9+nSZnVcIUb35OfmxbOgy3tj9BtuubOOtPW9xMfEiL3Z4EY1aY+vwhCi2NEM2jy87yKXYNGq62bP08U642utsHZYQojKxVe5Sgrxl586drF69mlOnTtGkifnWmTfeeIO5c+fi6OhIw4YNAVi/fn1ZRQuYc5fC8hYwt5hXiKT7/vvvv+X6nJwcli1bxuuvv46Pjw8AEydOZPbs2Sxbtoy5c+eiUqlwdnYu8hj//vsvERERjB49Gm9vb5ycnEhISKBmzZo4Ozvzyy+/MGLECBYsWMCAAQOKPI7BYMBgMFheJycnF/9Cs9LhQ//ib1+aXosEu8J/0MjPYDDw22+/8e2331qV63Q6dLoblfXx48dL1L28pOzs7ABzD4f8UlNT0ev1d3RMIYQoDkedI5/0+YSvj3/NN8e/YcXpFYQmhfLf3v/F2a7oekaIisKQncO07w9z9Goibg46lj/RmZpuMhe3EKKEbJW7FDNvAVizZg29e/e2JNwA9vb2mEwmSz5x/fr1EnUvvxN2dnaF5i1AmeYud3RPd1EuXbpEUlISnTp1sirv3LlzsYdjb9euHUuXLmXKlCmAOdF/4oknUKlULF26lEcffZQJEyYwdOjQW97PPWfOHN555507v5gK7vz586SkpBT4R7dnzx569uxpeb13794y7V5ev359wDzsfv4u5mFhYVb/qYQQoiyoVWqebvs0Ddwa8MbuN9gVsYsJmybw1cCvCHAum1+rhSgNOSaFGauOsetCHI52GpY93onGvndWFwshREV36NAhunTpYlW2f/9+TCYTPXr0ACiX7uX169cnLCzMqiwsLAyNRkNgYOAdH/d2SjXpTkhIAMw3pOfn6enJ+fPni3UMJycnFi9ezMKFC9ForLsIduvWjYsXL2IymVCrb32v06xZs3jxxRctr5OTk4v/Ruoczb/c2IKueIPO5d3kn3/ggKtXr7JmzRrLQHRQ+t3LMzMz2bJlC507d8bf3x8/Pz/atm3L2rVrGT58uCWO/fv3M2PGjFI7rxBC3MrQekMJdAnkub+f41LSJR7+42G+7P8lrbxb2To0IQpQFIXXfj3J5lNR2GnULHy0I+1q17j9jkIIURhb5S7FzFvAnLvcPODZ3LlzGTBggCVHa9y4cal3L9+7dy9gziMB7rnnHl599VWSkpJwc3MD4Oeff6Z///6Vp6U7r1tz/m7dABkZGVZdnovj5oQ7v9sl3GDuHnDHb5xKVeyuErbStGlT6tevz0svvcScOXNISEjg9ddfZ8SIEYwdO/aOj7t3716io6MJDg7GYDBY/uEPHDgQZ2dn4uLiGD16ND///DNjxowB4JNPPmHo0KG4ubnRsmVLPv/8c7p168YDDzxQGpcqhBDF0sKrBT8M+4Fn/nqGc9fP8fifjzOn1xwG1Rlk69CEsFAUhTmbz7L6UBhqFfxvfFt6Nip8JhghhCiWSpC7DBs2jA8//JAePXrQuHFjli5dSlBQEIcPH77jY0ZGRnLgwAHS0tIA2LdvH9nZ2VZd0PN6Puf1/J08eTLffvstQ4YMYerUqezZs4cdO3awa9euu7zCWyvVoTHr1KkDmN+A/CIjIy3rROnQarX88ccf+Pn5MXXqVD777DOeeeYZVqxYcVfH3bZtG8uWLeP69ev06dOHZcuWsWzZMsuAaQ4ODtx3331WgwwMGDCAvXv3YjQa2b59OxMnTmTbtm23/OFECCHKgp+TH8vvWU7vWr0x5Bh4cceLLD211Go6RVF8p0+f5rnnnqNnz55s3bq10G1OnjzJE088wYABA5g2bRohISHlHGXlsmDHJRYGmd+jj+5vzdCWNW0ckRBClL1Zs2YxY8YMPv/8c6ZPn45er+fw4cPUrl37jo955coVli1bxs8//8x9993HqVOnWLZsGQcPHrRs0717d7p37255bW9vz7///svw4cP5888/cXR05NChQ3Ts2PGuru92VModfhNRqVSFThnWpk0bOnXqxHfffQeYuz/XrFmTefPmleko2reTnJyMm5sbSUlJuLq6Wq3LzMwkNDSUevXqYW9vb6MIqxd5z4UQZSnblM1/D/6Xn87+BMCYxmN4rctr6NSVc1ToW9VhZeXbb7/l888/Z+rUqbzwwguF1vmnT5+mS5cujB8/nhEjRrBy5Up27NjB0aNHiz0CrC2uzVa+3XmJOZvPAvD6vc34v971bRyREKKyke/Q5e9W73lx67AStXQnJyezb98+9u3bB5gHTtu3b5/Vr9offvghy5Yt4/333+ePP/5g9OjR1K1bl4kTJ5bkVEIIIcQd06q1vNblNV7t/CoqVPxy/hee3v40KcYUW4dWaYwdO5bTp0/fcnyO9957j9atW7Nw4UJGjBjBTz/9hLOzM59++mn5BVpJfLcrxJJwvziosSTcQghRjZQo6T537hwzZsxgxowZdOnShc2bNzNjxgxWr15t2WbYsGFs3ryZo0ePMnfuXNq2bUtQUBAODjIFhhBCiPL1SLNH+Lzf5zhoHdh7bS+PbX6MyFQbDZRZybi7u992m+3btzNixAjLa41Gw/Dhw9m2bVsZRlb5LPk3lPf/OAPA8wMa8dyARjaOSAghRHkq0UBqnTp1srRy38qgQYMYNEgGrhFCCGF7/Wr3Y9nQZTzz1zNcTLxoHtl8wJe09Gpp69AqtbS0NOLi/r+9ew+Lqs7/AP4eLjNyv4ogoiEKgqCSqah4ATQzWbMyXRNT0RRco8JLmvoTV03dFNvYUrNc8QaKeSELLwtquWK6YoqXQAIBQS4iMih3OL8/XGYjFBGZOczM+/U8Pk+cOefM+4s+feYz55zv926j28jt7e2RlZX1xOMqKysbTLgql8uVlrEtiDx7C389ch0A8J5vN3wwgg03EZG2adXZy4mIiNoiNys37BmzB3+J/wtSi1Mx/eijmc1HdBkhdjSV2bp1KyIjI5vcJyIiAp6ens06X3V1NQA0er7NwMAAVVVVTzxuzZo1itlkNd3Oc5lYHnsNABA83AmhI50VS34SEZH2YNNNRERawdbIFjtG78CC0wvwU85PCD0VitC+oZjac6pWNEKjRo2Cq6trk/t07dr854xNTEygr6+PoqKiBtuLiopgaWn5xOMWL16M0NBQxc9yuVyxRqsm+fqndMUt5bOHdsXCUS5a8e+MiIgaY9NNRERaw0jfCJ/7fo5159chOiUaGy5uQHZpNhYPWAw9Hc0uiZ07d36upVn+SFdXFx4eHo3WWD1//nyTV8tlMhlkMlmr5WhrBEFAREIawk+kAnjUcC8a3YMNNxGRFmvVdbqJiIjauvqZzT/q9xEkkGBf6j58cPIDlFWXiR1N7QQGBiImJga//vpoVu7ExETEx8cjMDBQ5GTiEAQBa+N+VTTcoSOd2XATERGbbiIi0j4SiQQBbgHY6LMRMl0ZTt8+jZnHZ6KovOjpB2uJ//znP/D29oa3tzeAR8uDeXt7Y926dYp9goODMXHiRPTp0wc9e/aEj48P5s+fjzfeeEOs2KKprROw9NBVbPnx0TKqS8e4IsSvOxtuIiLi7eVERKS9/Dr74euXv8bchLlIvpuMKXFTsHnEZnQ2bb3bsNVV9+7dsXbt2kbbbW1tFf+to6ODLVu2YMWKFcjOzoajoyOsra1VGbNNqKiuRUjUJRy/ng+JBPjkdQ9M6s9/Q0RE9AibbnosQRBa9O18XV0ddHR4AwURqY8+Nn2wc/ROBP8rGNml2Qj4IQD/8PsHerXvJXY0UZmZmSmucj+Nra1tg2Zcm9x7WIWZkReQlHUfUj0dbJzQB2N62Ykdi4hIa6hD38LuiBTKysoQHh6Onj17QiaTwd7eHgsXLkR5eXmTx9XV1WHx4sWwsrKCVCpFv379cP78eRWlJiJ6fo5mjtj16i64WbmhuLIYM47NwKnsU2LHojYuq6gM4zedRVLWfZi208OuGQPYcBMRqcjJkycxevRomJubw9TUFKNGjcLly5efelxCQgI8PDygr6+PDh06YM2aNUrPyqabFI4dO4bc3FzExMSgrKwMsbGxiI6ORkhISJPHffLJJ9i6dSu+//573Lt3DwMHDsQrr7yCwsJCFSUnInp+1gbW+Oeof8Lb3hsVtRV4/+T72JeyT+xY1Eadz7iH17/8N9LvPoS9uQG+DR6E/o5PXiqNiIha1/r16/HBBx/g9u3byMzMhK2tLUaOHIm7d+8+8Zj09HT4+/tj/PjxkMvl2LVrF1avXo3NmzcrNatEEARBqe/QRsjlcpiZmaGkpASmpqYNXquoqEBGRgYcHR3Rrl07kRK2TE1NDSIjI3HhwgXk5+fj93+dYWFhcHd3R01NTZPnkMlkT7wl49NPP8XatWsbrcNar7a2Fra2tnj//fexdOlSAEBVVRVsbGywZMkSLFiw4LHHqfPvnIg0W3VdNVYmrsTBtIMAgHc93sV7nu+JOiFWUzVM3anj2Hb/nInlh6+hpk5Az46m2DatHzqYspYRkfKp+2fo8+fPY//+/bh16xaqqqoU2729vTFv3jxUVlY2ebyenh709B7/hHRxcTEsLS1x4MABvP7664/dZ+HChYiJiUFGRoZi29y5c3HixAmkpKQ89pimfufNrWG80q3GqqqqMGTIEGzcuBGenp4YPnw4Tp06hYsXL+KFF16Avb09PvzwQ5ibmzf558aNG098jzt37sDCwuKJr6elpeHu3bsYNmyYYptUKsXAgQORmJjYquMlIlIFfR19rBi0AnN6zwEAbE3eiqX/Xorq2mqRk5HYqmrqsORgMpYcvIqaOgH+veywP2gQG24iomb44osvMHz4cOjr68Pf3x8lJSU4fPgwLCws4OLiggsXLjy1b1m2bNkTz3/nzh0AaLJ3SUxMxNChQxts8/HxQWpqKu7du9c6A30MTqT2GIIgoLym6eeYlcVAz6DZV1NWr16N/Px8XL58GSYmJgCA8vJyrFu3Dp999hkAICIiAhERES3KkpycjM2bN2PlypVP3Cc/Px8A0L59+wbbbWxskJaW1qL3JSISm0QiQXCfYNga2WJF4grE/haLwrJChA8Ph7HUWOx4JILc++UIibqE/2QWQyIB5r/sgjnDnbgkGBGJTqze5Vn6lrS0NISGhmLv3r0YN24cAOCVV16Bra0t/P398ac//QnAo6vKLVFbW4uQkBC4u7s3OQlofn4+BgwY0GCbjY2N4jVLS+U8JtQmm+68vDxs2bIFFRUV+POf/4zevXsDANasWYPffvtNsV94eLhSbkUrrynHgD0Dnr6jEvz89s8w1Dds1r7bt2/HvHnzFA03ANjb26OsrEzxc01NTYtuL8/MzMSYMWMwZswYhIaGPjVLXV1do5/5QYSI1N3r3V+HtYE15p2eh8Q7iZh2dBq+HPElbAxtxI5GKnTsWh4W7r+CkvJqGMv08Pc/94GfawexYxERARCvd3mWvmXPnj1wcXFRNNzAo2ZXX19f0bsIgtCi28sFQcDs2bNx5coVnDlz5om3n9d7XN8CQKm9S5u7vbysrAz+/v7Q0dGBTCaDr68v0tPTAQCHDx9G9+7d4eXlBS8vL+jr64ucVjxpaWnIysrCoEGDGmz/9ddf4ebmpvi5JbeXZ2VlwcfHB3379sWePXua/AdoZ/doltaCgoIG2wsKCrR2+Rgi0ixDOg3BP1/5J6zaWSGlOAVTfpiCjJKMpx9Iaq+iuhbLD1/F7J0XUVJejV6dzPB9iDcbbiKiZ5SQkNCob8nIyEBVVZWid2nJ7eWCIGDOnDk4fPgwEhIS4Ozs3GQOOzu7x/YtAJTau7S5K936+vo4ceKE4l788+fPIyUlBV27dgUA/Pbbb7h//z4mTJgAAwMDpWQw0DPAz2//rJRzN+e9m0MulwMAdHV1FdsqKioQFRWFWbNmKbY96+3l2dnZ8PHxQa9evbBv377HfrFRUVEBqVQKHR0dODk5wdbWFidPnsTw4cMVr589exZ//etfm/2+RERtWU+rntj16i4E/SsImfJMTI2bii9HfAl3a3exo5GSJGUVY9G3V5Ca/wAAMGtoV8x/2QVSvTZ3vYKItJxYvUtz+xbgUe/y+74FALZt2wYnJyf07dsXANC/f/9nvr187ty52L9/P+Lj4+Hu3rgmV1c/mo+lvqfx9vbGjh07GqztHR8fDzc3N5ibmz/Tez+LZ266T548iU2bNiEpKQkbNmzAa6+91mifffv2YdOmTSgsLISnpydWrVqFLl26AACKiorw0UcfPfbc69atg5WVFUxMTDBz5kzcuXMHZmZmGDFiBADg448/RkFBATIyMuDr64tjx46hf//+zzqEp5JIJM2+VUIszs7OMDIywhdffIFNmzbh4cOHmD17NoyMjPD++++36Jy5ubnw9fVFjx49sHPnTtTW1qK2thYAFDP13b59Gw4ODoiJicH48eOho6ODBQsWYMWKFRg4cCDc3d2xdOlSGBoaYvr06a02XiIisXUy6YQdo3dgzr/m4FrRNQQeC8RnPp9hUMdBTz+Y1MbDyhqsP56C7WdvQRAAKyMp1k/oDR8XPlJARG2TOvQunp6eOHz4MBYsWIDOnTsjOjoa69evx5EjR1p8zpCQEERFReHYsWNwdnZWNOz6+vqKBr/+WfGjR48CAIKDgxEREYF58+Zh/vz5OHv2LLZv346vvvrqOUfYtGdqujdu3IjDhw8jKCgIMTExKC0tbbTP3r17MWXKFPz973+Hp6cn1q5di6FDh+Lq1aswMTGBTCaDl5fXY88vk8kAPPqH4+Xlhfz8fOzatQvXr19H7969MXbsWMW+ZmZm2Lt3r1KabnVgbGyMr776CrNnz8bu3btRXV2NkSNHIi4ursV3AMTGxiI7OxvZ2dmNJkYrKiqCkZERJBIJZDJZg2+qQkNDUVVVhTlz5uDevXt46aWXEB8fr9Rvi4iIxGDZzhLfjPoGH5z8AOfunMNf4v+C1YNX49Wur4odjZ6TIAg4mVKAZYeuIef+owmJ3njRHsvGuMHCSCpyOiIi9bZy5UokJSWha9euaNeuHaytrbFz5074+fm16HyVlZWKRnnIkCENXlu/fj3mzp0L4NGqSr/XqVMnnDhxAvPmzYOrqys6dOiA8PBwvPPOOy3K0VzPtE53VVWVIrhEIsHOnTsREBDQYJ8ePXpg5MiRiluay8rKYGNjg9WrVzfrCuxvv/0GMzMzWFtbA3jU0NnZ2TVa73nevHkwMDDAqlWrmpVdU9fprqioQFZWFqytrZU2254yqPPvnIiourYaH5/5GEdvPfrmfFH/RZjsOllp76eOa1k3V1sY2407cnzyww38dPMuAMDe3ACfvOGBYc7tn3IkEZFqqftn6KysLNTV1aFLly5qM+lya6zT/UxXuv/4TcEf5eXlISUlBevXr1dsMzQ0xLBhw3Dq1Klm3/Y8cuRIODs748GDB0hKSsLp06chl8sVs2jfunULN27cwLlz5554jsrKygaz39U/A61p2rVr99QJA4iIqHXp6+pj3dB1sGxniT2/7sHa82tRVF6E9zzfU5sPEQSk5pfiHwlpOHIlF3UCINXVwfTBLyDErzuMZG1u2hsiIrXXuXNnsSOIolUrSnZ2NoDGM7/Z2triypUrzTqHk5MTTp48ifj4eEilUvj4+MDY2Bjl5eXw8vKCRCLBuHHjMHz4cBgbP3mt1DVr1mDFihUtHwwREVETdCQ6WNR/EawMrBBxKQJbk7fiXsU9LPVaCj0dNmxtlSAISEwvwo6zmTh2PQ/19/uN8bDDotE94GDZtp+LJCIi9dOqnwrq1zj744zXUqlUMSFXc5ibm+PNN99ssM3AwAAzZ85s9jkWL17cYH1puVwOBweHZh9PRET0NBKJBLN6zYJlO0usPLcS3978FsUVxVg3dB3a6anfbX+a7HZxGb6/cgcxF28jreCBYvsrPW0x17cb3O3NRExHRESarFWb7vrnsIuKihpsv3v3ruI1VZHJZIqJ2YiIiJRpvPN4WMgssPDHhUjITsDsE7MR4RcBU6lmPX+tTh5U1uBSVjHOpRfhzM27uHy7RPGaoVQXr3vaY+qgF+DcwUTElEREpA1atel2dHSElZUVEhMT4evrC+C/t3ElJip9RjgiIiIx+XXxw+aRmxGSEIKkgiRMPzodm0dsRntDTsalKtv/nYGElEKk5Zcit6ThWq8SCTDA0RL+vTpibJ+OMG2n/4SzEBERta5Wbbp1dHQwa9YsfPHFF5g0aRIcHR2xceNGFBYWIjAwsDXfioiIqM3pZ9sP21/ZjqB/BSG1OBUxqTGY02eO2LG0xrVcOX5MLVT83NGsHby6WsHLyQrDXdrDxoS3/BMRkeo9U9N99uzZBles58+fj7CwMEyYMAGffPIJAGDFihXIy8tDjx49YGxsDKlUir1796Jbt26tm1wJnmH1NHpO/F0TkaZysXTBztE7Ef1rNGb3mi12HK3y+ov26NvFAt1sjNHNxhjmhlxfm4g0Ez9Lq05r/K6fqen29PTE0aNHG203M/vf5CP6+vrYtm0bPv/8c9y/fx92dnbQ1dV97qDKVJ+vqqoKBgYGIqfRDmVlZQAaT7pHRKQJOpl0wvx+88WOoXUGOVljkJPYKYiIlKf+s3NZWRn7FhVpjb7lmZpuAwODZl+xNjY2bnJJr7ZET08PhoaGKCwshL6+PnR0dMSOpLEEQUBZWRkKCgpgbm7e5r+QISIiIiJqK3R1dWFubo6CggIAgKGhISQSicipNFNr9i1cSBSPlnyxs7NDRkYGMjMzxY6jFczNzRut505ERERERE2r/wxd33iTcrVG38Km+7+kUim6d++OqqoqsaNoPH19fV7hJiIiIiJqgfoLhjY2NqiurhY7jkZrrb6FTffv6OjooF07zmxKRERERERtm66uLi9kqQk+vExERERERESkJGy6iYiIiIiIiJREa24vr19fTS6Xi5yEiIjo2dTXLk1cl5X1mYiI1FVz67PWNN2lpaUAAAcHB5GTEBERtUxpaSnMzMzEjtGqWJ+JiEjdPa0+SwRN/Nr8Merq6pCbmwsTE5MWrWUnl8vh4OCA7OxsmJqaKiFh26XNYwe0e/wcO8fOsbcNgiCgtLQUHTt2hI6OZj0Zxvrccto8dkC7x8+xc+wce9vQ3PqsNVe6dXR00KlTp+c+j6mpaZv6i1YlbR47oN3j59g5dm3TFseuaVe467E+Pz9tHjug3ePn2Dl2bdMWx96c+qxZX5cTERERERERtSFsuomIiIiIiIiUhE13M8lkMixfvhwymUzsKCqnzWMHtHv8HDvHrm20eezqSpv/zrR57IB2j59j59i1jbqPXWsmUiMiIiIiIiJSNV7pJiIiIiIiIlISNt1ERERERERESsKmm4iIiIiIiEhJtGad7ueVlpaGkpISuLm5wcDAQOw4KlVYWIiUlBS4ubnB0tJS7DgqlZWVhfv376Nbt24wNDQUO45KlZaW4ubNm7CxsWmVNXTV0c2bN5Gfn48+ffrA2NhY7DhKl5ycjJKSkgbb2rdvDxcXF5ESqV5tbS1u3LgBAwMDODk5iR2HmoH1mfWZ9Vn7sD6zPqsdgZpUWFgoDBo0SDAzMxO6desmmJmZCQcOHBA7lkokJycLkyZNEjp06CAAEA4ePCh2JJX59ttvBVdXV8HBwUFwd3cXjIyMhFWrVokdSyXy8vKEgIAAwcrKSnjxxRcFMzMzYcCAAUJGRobY0VQqPT1dsLCwEAAIFy5cEDuOSgwePFhwcHAQBg8erPizfPlysWOpzKFDhwQ7OzvB0dFR8PDwEHx9fYWCggKxY9ETsD6zPrM+sz6zPmsHTajPvL38KYKDg1FRUYGcnBzcvHkTS5cuxeTJk5GTkyN2NKW7du0axowZg2vXrokdReUyMzNx4MABZGVlITk5GYcOHcLy5csRGxsrdjSlu337Nt58800UFhbi4sWLyM3NBQDMnTtX5GSqU11djUmTJmHWrFliR1G5wMBAnDlzRvEnLCxM7EgqcfbsWbz55psICwtDeno6rly5grCwMOTl5YkdjZ6A9Zn1mfWZ9VmbsD6rd31m092E4uJiHDx4EKGhoTAyMgIAhISEQCqVIioqSuR0yjdx4kRMnjxZbdfDex4ffvghevToofh5xIgR6N69O86cOSNiKtXo27cvxo0bB4lEAgAwNDTEgAEDtOKDbL0lS5bA0dERAQEBYkdROblcjgsXLuD27dtiR1GpsLAwDBs2rMEHuSFDhsDDw0PEVPQkrM+sz/VYn1mftQXrs3rXZz7T3YTk5GTU1taib9++im1SqRS9e/fGpUuXRExGqlZYWIjMzEx069ZN7Cgqk5SUBLlcjuTkZERFRWHLli1iR1KJ48ePY9++ffjll1+0rrABwObNm5GQkICMjAy88MILiIyMRJ8+fcSOpVQ1NTX48ccfsXbtWjx48AApKSmws7NDx44dxY5GT8D6TPVYn1mftQXrs3rXZzbdTbh37x4AwMrKqsF2KysrxWuk+erq6jBjxgzY2dlh8uTJYsdRmQ0bNiA1NRUpKSl4+eWX4e3tLXYkpcvPz8e0adMQFRUFc3NzrSvqQUFBeOONN2BoaIiHDx9iypQpGDduHK5evarRE9UUFRWhsrIS165dg7OzM2xtbZGWlob+/fsjKioK7du3Fzsi/QHrMwGsz6zP2oP1Wf3rM28vb4K+vj4AoKKiosH28vJySKVSMSKRigmCgKCgIJw7dw7fffed4jZGbbB7925cuHABOTk5KCkpwdixY8WOpHTvvfceXnzxRejq6uLMmTOKK2aXL19GSkqKyOmULyAgQDELsJGREcLDw5GZmYnExESRkylX/f/rjx8/jkuXLiEpKQm3bt1CTk4OQkNDRU5Hj8P6TKzPrM8A6zPrs/rgle4mdOnSBQCQk5MDBwcHxfacnBwMGzZMrFikIoIgYM6cOTh48CASEhLg5uYmdiRRmJiYICgoCOPHj0dxcTEsLCzEjqQ0FhYWyM3NxaJFiwAAZWVlAICIiAhkZ2drzaQl9Tp06AAAGv+8oIWFBUxNTTFu3DjFmC0tLTFhwgRERkaKnI4eh/VZu7E+P8L6zPrM+qw+2HQ3wc3NDR07dkRsbCy8vLwAAOnp6UhOTsaqVatETkfKNnfuXOzfvx/x8fFqN1nD83j48GGjKwZpaWkwMDDQ+CsJf3wu7urVq/Dw8MDXX3+Nl156SaRUqlFeXg6ZTAYdnf/dAHX8+HEAgLu7u1ixVEIikeDll19u9OHl9u3banXrmjZhfdZurM//w/rM+qzJNKk+s+lugo6ODtasWYOZM2fCysoKjo6OWLlyJby9veHv7y92PKW7d+8erl+/jvLycgDAjRs3YG1tDXt7ezg6OoqcTrkWLlyILVu2ICIiAnK5XDErqo2NDZydnUVOp1xhYWEoLS2Fn58fTExMcO7cOfztb3/Dxx9/zNs2NVhGRgamTZuGGTNmwNHREcnJyVi9ejUmTZqk8R9oAGDFihUYOHAgli1bhiFDhuDChQvYsWOHVsyErY5Yn1mfWZ9Zn7UF67Nm1GeJIAiC2CHautjYWERGRkIul2PgwIFYsGABTExMxI6ldD/99BMWL17caPuECRMQEhIiQiLVeeedd5Cent5o+6hRo7Bs2TIREqlOXV0doqKicOTIERQXF6NLly4ICAjAkCFDxI6mchkZGZgyZQq++eYbuLi4iB1H6a5fv45NmzYhNTUVdnZ28Pf3x/jx48WOpTLXr19HeHg4bt26BQcHB8yYMUMrJihSZ6zPDbE+sz5rC9Zn1md1q89suomIiIiIiIiUhLOXExERERERESkJm24iIiIiIiIiJWHTTURERERERKQkbLqJiIiIiIiIlIRNNxEREREREZGSsOkmIiIiIiIiUhI23URERERERERKwqabiERx5coVxMfHo66uTuwoRERE9F+sz0StT0/sAESkfXbv3o3/+7//Q4cOHdCtWzfs2LFD7EhERERaj/WZSDnYdBNpmYyMDCQnJ2Ps2LHN2j81NRVJSUkAgN69e8PV1fW5M+zfvx8HDx6Eh4cHnJ2dIQgCJBKJ4vW8vDycOnUKAODo6IgBAwY893sSERG1ZazPRJqLt5cTaZnTp09jzpw5zd7/hx9+QFBQEA4dOoSbN282eK20tBQHDhzA6dOnn3j8w4cPsXfv3ga3qfn4+GD9+vX49NNP4erq2qCgA0B+fj4OHTqERYsWYdOmTc3OSkREpK5Yn4k0F690E9FTderUCdHR0Yqfy8rKMG/ePBw6dAjAo2/Yhw0b9thjf/jhB6xatQoTJ05UbJs6dSrmz5+P2traBuet17t3b0RHRyMgIKB1B0JERKRBWJ+J1AOvdBMRfvnlF+zbtw+VlZXN2r+yshK9evVCamoq/Pz8mtz3yJEj8Pf3b7AtKioK9vb2CAwM5PNiRERET8D6TKQZ2HQTabnY2FgMHToU1dXVkMlkzTrGwsICwcHBMDExaXK/uro6xMXFNSrq27Ztw/Tp0/Huu+8iLi4Od+7caXF+IiIiTcT6TKQ52HQTabHt27cjICAA0dHRmDx5cquf/+eff0ZdXR28vLwU265evYqLFy9i2rRp6N+/P9zc3PhtOhER0e+wPhNpFjbdRFpqw4YNCA0NRVxcHF599VWlvMeRI0cwevRo6OrqKrZ98803GDFiBDp37gwAmDFjBrZt26aU9yciIlI3rM9EmodNN5EWunv3LhYuXIglS5Zg8ODBSnufPz4vVl1djd27d8PBwQHR0dGIjo6Gnp4eUlNTcebMGaXlICIiUgesz0SaibOXE2kha2trrF69GkFBQXBxcWn0TFdryM7Oxo0bNzBq1CjFttjYWNTU1ODBgweKmVUBwN3dHdu2bYO3t3er5yAiIlIXrM9EmolNN5GWmjp1KgRBwFtvvYWYmJhWL+xHjhzB4MGDYW5urti2bds2TJgwAZs3b26wb1xcHCZMmIDPP/8cxsbGrZqDiIhInbA+E2keNt1EWmzatGmKwr5//36MGTOm2ccePnwY5eXlyMzMRGlpKaKjo6Grq4u33noLQONb13Jzc3Hs2DEcOXKk0bn8/PwgkUiwb98+BAYGPv/AiIiI1BjrM5FmYdNNpGUcHR3x2muvKX6ePn06pFIpYmJi0K9fP9jY2DQ6pqSkBNHR0ejduzdcXV0BPPr2+/79+7C3twcAHDp0CFKpFG+99RbKy8tx8uRJhIeHK85x/fp1TJw4Eb6+vo3OL5VKsWTJEhQVFQEA8vLycOrUKWRmZsLJyalVx09ERNQWsT4TaS6JIAiC2CGIqO2Ki4tDZGQkAODtt9/G2LFjn3rMd999h9DQUNy8ebNF73n58mWsWbMGADBs2DAEBwe36DxERESaivWZSH2w6SaiVrd79248fPgQs2bNEjsKERER/RfrM5E42HQTERERERERKQnX6SYiIiIiIiJSEjbdRERERERERErCppuIiIiIiIhISdh0ExERERERESkJm24iIiIiIiIiJWHTTURERERERKQkbLqJiIiIiIiIlIRNNxEREREREZGSsOkmIiIiIiIiUhI23URERERERERK8v82/M4MNpl/zwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x320 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "k = torch.linspace(0.05, 2 * np.pi, 300)\n",
    "fig, ax = plt.subplots(1, 2, figsize=(10, 3.2))\n",
    "for sigma in (0.5, 1.0, 2.0):\n",
    "    ew = EwaldSummation(sigma=sigma)\n",
    "    ax[0].plot(k, ew._kfac(k ** 2), label=f\"$\\sigma$={sigma}\")\n",
    "    ew6 = EwaldSummation(sigma=sigma, exponent=6)\n",
    "    ax[1].plot(k, ew6._kfac(k ** 2), label=f\"$\\sigma$={sigma}\")\n",
    "ax[0].set_title(\"electrostatics kernel (p=1)\"); ax[0].set_yscale(\"log\")\n",
    "ax[1].set_title(\"dispersion kernel (p=6)\")\n",
    "for a in ax:\n",
    "    a.set_xlabel(\"k [1/Å]\"); a.legend()\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd32b23f",
   "metadata": {},
   "source": [
    "## Block 2: Physics checks\n",
    "\n",
    "* **Exact invariances**: rotation (including triclinic cells), translation,\n",
    "  and shifts by full lattice vectors. One deliberate improvement over\n",
    "  upstream: floating-point ties at the $|k| = k_c$ shell and in the grid size\n",
    "  are resolved consistently, so rotating a cell can never drop a k shell\n",
    "  (upstream truncates exactly and loses ~$10^{-4}$-level invariance for\n",
    "  rotated cells).\n",
    "* **Tinfoil boundary conditions**: for a *neutral* $q$, an isolated cluster in\n",
    "  a large periodic box converges to the molecular direct sum; a net $q$ adds\n",
    "  the constant background of the omitted $k=0$ term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2e594b27",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-22T03:15:18.338314Z",
     "iopub.status.busy": "2026-09-22T03:15:18.338231Z",
     "iopub.status.idle": "2026-09-22T03:15:18.357513Z",
     "shell.execute_reply": "2026-09-22T03:15:18.356933Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rotation diff (triclinic): 1.1102230246251565e-16\n",
      "lattice-shift diff       : 1.1102230246251565e-16\n",
      "isolated cluster, 30 Å box vs direct sum (neutral q): 7.723798462899722e-05\n"
     ]
    }
   ],
   "source": [
    "Q, _ = np.linalg.qr(np.random.default_rng(3).normal(size=(3, 3)))\n",
    "if np.linalg.det(Q) < 0: Q[:, 0] *= -1\n",
    "Qt = torch.tensor(Q)\n",
    "cell_t = torch.tensor([[6.0, 0, 0], [1.2, 5.8, 0], [-0.7, 0.9, 6.3]])\n",
    "ew = EwaldSummation()\n",
    "e0 = ew.reciprocal(pos, q, cell_t)\n",
    "print(\"rotation diff (triclinic):\", abs(float(e0) - float(ew.reciprocal(pos @ Qt.T, q, cell_t @ Qt.T))))\n",
    "print(\"lattice-shift diff       :\", abs(float(e0) - float(ew.reciprocal(pos + cell_t[1], q, cell_t))))\n",
    "\n",
    "qn = q - q.mean(dim=0, keepdim=True)      # neutral hidden variable\n",
    "cluster = torch.tensor(np.random.default_rng(2).uniform(14, 16, (6, 3)))\n",
    "qc = torch.tensor(np.random.default_rng(2).normal(size=(6, 1)));  qc -= qc.mean()\n",
    "ew_fine = EwaldSummation(dl=1.0)\n",
    "print(\"isolated cluster, 30 Å box vs direct sum (neutral q):\",\n",
    "      abs(float(ew_fine.realspace(cluster, qc))\n",
    "          - float(ew_fine.reciprocal(cluster, qc, torch.eye(3, dtype=torch.float64) * 30))))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90a73ea0",
   "metadata": {},
   "source": [
    "## Block 3: The latent charge head and the `node_features` contract  ·  eq 2\n",
    "\n",
    "Every xnn model now returns its **invariant per-atom features** under the\n",
    "``\"node_features\"`` output key (CACE: the concatenated symmetrized B features,\n",
    "exactly what the paper uses; MACE/NequIP: the scalar channels; Allegro: the\n",
    "environment-aggregated edge latents; SchNet/PhysNet: the feature vectors;\n",
    "HDNNP/ANI: the descriptors) and advertises the width as ``node_feature_dim``.\n",
    "`LatentEwald` puts the reference `cace-lr-fit` head on top: a bias-free\n",
    "[24, 12] MLP plus a parallel bias-free linear layer, mapping features to the\n",
    "``n_channels``-dimensional $q$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4186ddd6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-22T03:15:18.358976Z",
     "iopub.status.busy": "2026-09-22T03:15:18.358903Z",
     "iopub.status.idle": "2026-09-22T03:15:18.364529Z",
     "shell.execute_reply": "2026-09-22T03:15:18.364069Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LatentEwald around CACE | node_feature_dim = 384\n",
      "q head: Sequential(\n",
      "  (0): Linear(in_features=384, out_features=24, bias=False)\n",
      "  (1): SiLU()\n",
      "  (2): Linear(in_features=24, out_features=12, bias=False)\n",
      "  (3): SiLU()\n",
      "  (4): Linear(in_features=12, out_features=4, bias=False)\n",
      ") \n",
      "+ parallel Linear(in_features=384, out_features=4, bias=False)\n"
     ]
    }
   ],
   "source": [
    "from xnn.common.config import from_dict\n",
    "from xnn.common.models import build_model, ForceStressOutput\n",
    "\n",
    "cfg = from_dict({\"model\": {\"name\": \"cace\", \"cutoff\": 4.5, \"n_interactions\": 1, \"n_rbf\": 6,\n",
    "    \"extra\": {\"species\": [1, 8], \"n_atom_basis\": 2, \"n_radial_basis\": 8,\n",
    "              \"max_l\": 3, \"max_nu\": 3, \"avg_num_neighbors\": 9.0, \"embed_receiver_nodes\": True,\n",
    "              \"long_range\": {\"n_channels\": 4, \"dl\": 1.9, \"sigma\": 1.0}}}})\n",
    "model = build_model(cfg.model)          # the config hook wraps automatically\n",
    "print(type(model).__name__, \"around\", type(model.model).__name__,\n",
    "      \"| node_feature_dim =\", model.model.node_feature_dim)\n",
    "print(\"q head:\", model.q_net, \"\\n+ parallel\", model.q_linear)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2efcdc9",
   "metadata": {},
   "source": [
    "## Capstone: whole-model CACE-LR transplant vs the original\n",
    "\n",
    "The original composes `Cace` + `Atomwise` (SR energy) + `Atomwise` ($q$ head)\n",
    "+ `EwaldPotential` + `FeatureAdd`; xnn composes `CACE` + `LatentEwald`. With\n",
    "every weight copied across, the two are the same function (float32: the\n",
    "upstream Ewald path is float32-only; each ingredient is verified at ~1e-16 in\n",
    "float64 above and in `tests/test_cace.py`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4b679d45",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-22T03:15:18.366003Z",
     "iopub.status.busy": "2026-09-22T03:15:18.365936Z",
     "iopub.status.idle": "2026-09-22T03:15:18.545679Z",
     "shell.execute_reply": "2026-09-22T03:15:18.545141Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|dE|/|E|      = 6.36e-07\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max|dF|/max|F| = 1.09e-06\n",
      "max|dq|/max|q| = 2.09e-07\n",
      "E_lr: xnn 4.02507e+08 | original 4.02507e+08\n"
     ]
    }
   ],
   "source": [
    "torch.set_default_dtype(torch.float32)\n",
    "from cace.modules import BesselRBF as UpBessel, PolynomialCutoff as UpPoly, FeatureAdd\n",
    "from cace.modules.atomwise import Atomwise\n",
    "from cace.modules.forces import Forces\n",
    "from cace.representations import Cace as UpCace\n",
    "from cace.models.atomistic import NeuralNetworkPotential\n",
    "from xnn.common.data import structure_to_graph\n",
    "\n",
    "CUT, NRBF, NRB, NAB, LMAX, NU, T, NQ = 4.5, 6, 8, 2, 3, 3, 1, 4\n",
    "torch.manual_seed(7)\n",
    "rep = UpCace(zs=[1, 8], n_atom_basis=NAB, cutoff=CUT,\n",
    "             radial_basis=UpBessel(cutoff=CUT, n_rbf=NRBF, trainable=True),\n",
    "             cutoff_fn=UpPoly(cutoff=CUT, p=6), max_l=LMAX, max_nu=NU,\n",
    "             num_message_passing=T, type_message_passing=[\"M\", \"Ar\", \"Bchi\"],\n",
    "             n_radial_basis=NRB, avg_num_neighbors=9.0, embed_receiver_nodes=True)\n",
    "sr = Atomwise(n_layers=3, n_hidden=[32, 16], output_key=\"SR_energy\", add_linear_nn=True)\n",
    "qhead = Atomwise(n_layers=3, n_hidden=[24, 12], n_out=NQ, per_atom_output_key=\"q\",\n",
    "                 output_key=\"tot_q\", residual=False, add_linear_nn=True, bias=False)\n",
    "ep = cace.modules.EwaldPotential(dl=1.9, sigma=1.0, feature_key=\"q\",\n",
    "                                 output_key=\"ewald_potential\",\n",
    "                                 remove_self_interaction=False, aggregation_mode=\"sum\")\n",
    "eadd = FeatureAdd(feature_keys=[\"SR_energy\", \"ewald_potential\"], output_key=\"CACE_energy\")\n",
    "nnp = NeuralNetworkPotential(representation=rep,\n",
    "                             output_modules=[sr, qhead, ep, eadd,\n",
    "                                             Forces(energy_key=\"CACE_energy\",\n",
    "                                                    forces_key=\"CACE_forces\")])\n",
    "\n",
    "cfg = from_dict({\"model\": {\"name\": \"cace\", \"cutoff\": CUT, \"n_interactions\": T, \"n_rbf\": NRBF,\n",
    "    \"extra\": {\"species\": [1, 8], \"n_atom_basis\": NAB, \"n_radial_basis\": NRB,\n",
    "              \"max_l\": LMAX, \"max_nu\": NU, \"avg_num_neighbors\": 9.0,\n",
    "              \"embed_receiver_nodes\": True,\n",
    "              \"long_range\": {\"n_channels\": NQ, \"dl\": 1.9, \"sigma\": 1.0}}}})\n",
    "x = build_model(cfg.model)\n",
    "\n",
    "rng = np.random.default_rng(11)\n",
    "pos_np = rng.uniform(0, 6, (12, 3))\n",
    "g = structure_to_graph({\"pos\": pos_np, \"atomic_numbers\": [1, 8] * 6,\n",
    "                        \"cell\": np.eye(3) * 6.0, \"pbc\": [True] * 3}, CUT)\n",
    "cell = g.cell[0]\n",
    "data = {\"positions\": g.pos.clone().requires_grad_(True), \"atomic_numbers\": g.atomic_numbers,\n",
    "        \"edge_index\": g.edge_index, \"shifts\": g.cell_shifts.to(cell.dtype) @ cell,\n",
    "        \"unit_shifts\": g.cell_shifts.to(cell.dtype), \"batch\": g.batch,\n",
    "        \"cell\": cell.unsqueeze(0)}\n",
    "out_up = nnp(data, training=True)          # lazy-inits both Atomwise heads\n",
    "\n",
    "with torch.no_grad():\n",
    "    b = x.model\n",
    "    b.embed_sender.copy_(rep.node_embedding_sender.embedding_weights)\n",
    "    b.embed_receiver.copy_(rep.node_embedding_receiver.embedding_weights)\n",
    "    b.rbf.freqs.copy_(rep.radial_basis.bessel_weights * CUT)\n",
    "    b.radial_transform.weight.copy_(torch.stack(list(rep.radial_transform.weights)))\n",
    "    for t, (memory, ar, bchi) in enumerate(rep.message_passing_list):\n",
    "        xi = b.interactions[t]\n",
    "        xi.memory.memory_coef.copy_(torch.stack(list(memory.memory_coef)))\n",
    "        xi.message_ar.prefactor.copy_(torch.stack(list(ar.prefactor)))\n",
    "        xi.message_ar.inv_r0.copy_(torch.stack(list(ar.invr0)))\n",
    "        xi.message_bchi.h.weight.copy_(bchi.hnet[0].linear.weight)\n",
    "        xi.message_bchi.h.bias.copy_(bchi.hnet[0].linear.bias)\n",
    "    for j, dense in enumerate(sr.outnet):\n",
    "        b.readout_mlp[2 * j].weight.copy_(dense.linear.weight)\n",
    "        b.readout_mlp[2 * j].bias.copy_(dense.linear.bias)\n",
    "    b.readout_linear.weight.copy_(sr.linear_nn.linear.weight)\n",
    "    b.readout_linear.bias.copy_(sr.linear_nn.linear.bias)\n",
    "    for j, dense in enumerate(qhead.outnet):\n",
    "        x.q_net[2 * j].weight.copy_(dense.linear.weight)\n",
    "    x.q_linear.weight.copy_(qhead.linear_nn.linear.weight)\n",
    "\n",
    "E_up = out_up[\"CACE_energy\"].sum()\n",
    "F_up = -torch.autograd.grad(E_up, data[\"positions\"])[0]\n",
    "ox = ForceStressOutput(x)(g)\n",
    "print(f\"|dE|/|E|      = {abs(float(ox['energy']) - float(E_up))/abs(float(E_up)):.2e}\")\n",
    "print(f\"max|dF|/max|F| = {((ox['forces'].detach() - F_up).abs().max()/F_up.abs().max()).item():.2e}\")\n",
    "print(f\"max|dq|/max|q| = {((ox['latent_charges'].detach() - out_up['q']).abs().max()/out_up['q'].abs().max()).item():.2e}\")\n",
    "print(f\"E_lr: xnn {float(ox['energy_lr']):.6g} | original {float(out_up['ewald_potential']):.6g}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd6db2d0",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "| block | original `cace` | `xnn` | diff |\n",
    "|---|---|---|---|\n",
    "| latent charge head (eq 2) | `Atomwise` (bias-free MLP + linear) | `LatentEwald.q_net` + `q_linear` | 0 (same layers) |\n",
    "| reciprocal Ewald, p=1 (eqs 3–4) | `EwaldPotential.compute_potential` | `EwaldSummation.reciprocal` | ~1e-16 (float64) |\n",
    "| dispersion kernel, p=6 (eq 5) | same, `exponent=6` | same | ~1e-16 (float64) |\n",
    "| molecular direct sum | `compute_potential_realspace` | `EwaldSummation.realspace` | ~1e-16 (float64) |\n",
    "| production triclinic path | `compute_potential_triclinic` (float32-only) | same math, dtype-safe | ~1e-7 rel (float32) |\n",
    "| **whole CACE-LR model** | rep + 2×Atomwise + Ewald + FeatureAdd | `LatentEwald(CACE)` | **float32 round-off** |\n",
    "\n",
    "Documented deviations (both floating-point robustness fixes, not physics\n",
    "changes): the k-grid follows the input dtype (upstream hard-casts to float32,\n",
    "crashing float64), and exact ties at the $|k| = k_c$ shell resolve\n",
    "consistently so the energy is exactly rotation-invariant. Also note upstream's\n",
    "`remove_self_interaction=True` over-subtracts the self term for\n",
    "multi-dimensional $q$ (once per channel); xnn subtracts it once (identical\n",
    "for 1-D $q$ and for the default `remove_self_interaction=False`).\n",
    "\n",
    "Because `LatentEwald` only needs the `\"node_features\"` contract, the same\n",
    "wrapper (and the `model.extra[\"long_range\"]` config hook) works for **every**\n",
    "xnn model; see `tests/test_les.py::test_wraps_every_model`. Next:\n",
    "`les_water_train_test.ipynb` reproduces the paper's bulk-water experiment."
   ]
  }
 ],
 "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.13.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}