GitHub_collection_pykan/docs/Examples/Example_6_PDE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5d904dee",
   "metadata": {},
   "source": [
    "# Example 6: Solving Partial Differential Equation (PDE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d568912",
   "metadata": {},
   "source": [
    "We aim to solve a 2D poisson equation $\\nabla^2 f(x,y) = -2\\pi^2{\\rm sin}(\\pi x){\\rm sin}(\\pi y)$, with boundary condition $f(-1,y)=f(1,y)=f(x,-1)=f(x,1)=0$. The ground truth solution is $f(x,y)={\\rm sin}(\\pi x){\\rm sin}(\\pi y)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0e2bc449",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "pde loss: 5.92e+00 | bc loss: 7.98e-02 | l2: 3.07e-02 : 100%|█| 20/20 [00:18<00:\n"
     ]
    }
   ],
   "source": [
    "from kan import KAN, LBFGS\n",
    "import torch\n",
    "import matplotlib.pyplot as plt\n",
    "from torch import autograd\n",
    "from tqdm import tqdm\n",
    "\n",
    "dim = 2\n",
    "np_i = 21 # number of interior points (along each dimension)\n",
    "np_b = 21 # number of boundary points (along each dimension)\n",
    "ranges = [-1, 1]\n",
    "\n",
    "model = KAN(width=[2,2,1], grid=5, k=3, grid_eps=1.0, noise_scale_base=0.25)\n",
    "\n",
    "def batch_jacobian(func, x, create_graph=False):\n",
    "    # x in shape (Batch, Length)\n",
    "    def _func_sum(x):\n",
    "        return func(x).sum(dim=0)\n",
    "    return autograd.functional.jacobian(_func_sum, x, create_graph=create_graph).permute(1,0,2)\n",
    "\n",
    "# define solution\n",
    "sol_fun = lambda x: torch.sin(torch.pi*x[:,[0]])*torch.sin(torch.pi*x[:,[1]])\n",
    "source_fun = lambda x: -2*torch.pi**2 * torch.sin(torch.pi*x[:,[0]])*torch.sin(torch.pi*x[:,[1]])\n",
    "\n",
    "# interior\n",
    "sampling_mode = 'random' # 'radnom' or 'mesh'\n",
    "\n",
    "x_mesh = torch.linspace(ranges[0],ranges[1],steps=np_i)\n",
    "y_mesh = torch.linspace(ranges[0],ranges[1],steps=np_i)\n",
    "X, Y = torch.meshgrid(x_mesh, y_mesh, indexing=\"ij\")\n",
    "if sampling_mode == 'mesh':\n",
    "    #mesh\n",
    "    x_i = torch.stack([X.reshape(-1,), Y.reshape(-1,)]).permute(1,0)\n",
    "else:\n",
    "    #random\n",
    "    x_i = torch.rand((np_i**2,2))*2-1\n",
    "\n",
    "# boundary, 4 sides\n",
    "helper = lambda X, Y: torch.stack([X.reshape(-1,), Y.reshape(-1,)]).permute(1,0)\n",
    "xb1 = helper(X[0], Y[0])\n",
    "xb2 = helper(X[-1], Y[0])\n",
    "xb3 = helper(X[:,0], Y[:,0])\n",
    "xb4 = helper(X[:,0], Y[:,-1])\n",
    "x_b = torch.cat([xb1, xb2, xb3, xb4], dim=0)\n",
    "\n",
    "steps = 20\n",
    "alpha = 0.1\n",
    "log = 1\n",
    "\n",
    "def train():\n",
    "    optimizer = LBFGS(model.parameters(), lr=1, history_size=10, line_search_fn=\"strong_wolfe\", tolerance_grad=1e-32, tolerance_change=1e-32, tolerance_ys=1e-32)\n",
    "\n",
    "    pbar = tqdm(range(steps), desc='description')\n",
    "\n",
    "    for _ in pbar:\n",
    "        def closure():\n",
    "            global pde_loss, bc_loss\n",
    "            optimizer.zero_grad()\n",
    "            # interior loss\n",
    "            sol = sol_fun(x_i)\n",
    "            sol_D1_fun = lambda x: batch_jacobian(model, x, create_graph=True)[:,0,:]\n",
    "            sol_D1 = sol_D1_fun(x_i)\n",
    "            sol_D2 = batch_jacobian(sol_D1_fun, x_i, create_graph=True)[:,:,:]\n",
    "            lap = torch.sum(torch.diagonal(sol_D2, dim1=1, dim2=2), dim=1, keepdim=True)\n",
    "            source = source_fun(x_i)\n",
    "            pde_loss = torch.mean((lap - source)**2)\n",
    "\n",
    "            # boundary loss\n",
    "            bc_true = sol_fun(x_b)\n",
    "            bc_pred = model(x_b)\n",
    "            bc_loss = torch.mean((bc_pred-bc_true)**2)\n",
    "\n",
    "            loss = alpha * pde_loss + bc_loss\n",
    "            loss.backward()\n",
    "            return loss\n",
    "\n",
    "        if _ % 5 == 0 and _ < 50:\n",
    "            model.update_grid_from_samples(x_i)\n",
    "\n",
    "        optimizer.step(closure)\n",
    "        sol = sol_fun(x_i)\n",
    "        loss = alpha * pde_loss + bc_loss\n",
    "        l2 = torch.mean((model(x_i) - sol)**2)\n",
    "\n",
    "        if _ % log == 0:\n",
    "            pbar.set_description(\"pde loss: %.2e | bc loss: %.2e | l2: %.2e \" % (pde_loss.cpu().detach().numpy(), bc_loss.cpu().detach().numpy(), l2.detach().numpy()))\n",
    "\n",
    "train()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2246bab",
   "metadata": {},
   "source": [
    "Plot the trained KAN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "02e2a0ba",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABcy0lEQVR4nO3dd1RU19oG8OfM0AZB6XZREUQBFURsscaS2EuikZLEJMauiMZuLLHmRgRL7MmNIFhijcbeEhUVUIwKiIgKgiBV+tT9/cGd80mihjIw7f2tdddd6wrM5s5snrPbuznGGAMhhBCiQgJ1N4AQQojuoXAhhBCichQuhBBCVI7ChRBCiMpRuBBCCFE5ChdCCCEqR+FCCCFE5ShcCCGEqByFCyGEEJWjcCGEEKJyFC6EEEJUjsKFEEKIylG4EEIIUTkKF0IIISpH4UIIIUTlDNTdAEK0AWMM2dnZKCwshJmZGaytrcFxnLqbRYjGopELIe+Ql5eH4OBgODo6wtbWFi1atICtrS0cHR0RHByMvLw8dTeREI3E0U2UhLzZmTNnMHr0aBQXFwMoG70oKUctpqamOHToEAYOHKiWNhKiqShcCHmDM2fOYPDgwWCMQaFQvPXrBAIBOI7DyZMnKWAIeQ2FCyF/k5eXhyZNmqCkpOSdwaIkEAggEonw/PlzWFhY1HwDCdECtOZCyN/88ssvKC4urlCwAIBCoUBxcTH27NlTwy0jRHvQyIWQ1zDG4OjoiKSkJFSma3Ach5YtW+LRo0e0i4wQULgQUk5WVhZsbW2r9f3W1tYqbBEh2ommxQh5TWFhYbW+v6CgQEUtIUS7UbgQ8pqcnJxqfb+5ubmKWkKIdqNwIXovNzcXu3btQt++fdGtWzcIBFXrFgKBAEuWLMGNGzcqtV5DiC6iNReil8RiMU6fPo2wsDCcOnUKMpkM77//Pry9vZGSkoL58+dXOiD69u2L+Ph4pKWloVWrVvD19YWPjw9atWpVQ78FIZqLwoXoDcYYIiIiEBYWhkOHDiE3NxcdOnSAt7c3xo4di/r16wOo3jkXc3NzXL58GSEhITh06BAKCwvRtWtX+Pr6YuzYsbTYT/QGhQvReQkJCQgPD0d4eDiePn2KJk2awNvbG5988gnatm37xu+p7An933//HQMGDCj3b8XFxTh27BhCQ0Nx5swZCAQCfPjhh/Dz88OQIUNgYmKi0t+TEE1C4UJ0UmZmJg4ePIiwsDBERUWhbt26GDVqFHx8fNC9e/cKratUtLbY4cOH/xEsf5eRkYH9+/cjJCQEUVFRqFevHj7++GP4+vqiR48eVV7nIURTUbgQnVFcXIwTJ04gPDwcZ8+eBcdx+OCDDzBu3DgMGjQIIpGo0j8zLy8Pe/bswcaNG/H48WP+f3dwcMCMGTPw2WefoV69epX6mfHx8QgNDUVoaCiePXuGZs2awcfHB35+fmjTpk2l20iIJqJwIVpNLpfjzz//xN69e3H06FEUFBTAy8sLPj4++Oijj1S2xsEYQ05ODgoKCmBubg4rK6tqn8RXKBS4fv06QkJCcODAAeTl5cHDwwN+fn745JNP0KBBA5W0nRB1oHAhWun+/fv8OkpaWhpatmwJb29vjBs3Dg4ODupuXqWJxWKcPHkSoaGhOHHiBORyOQYMGABfX1+MGDECderUUXcTCakUCheiNV68eIF9+/YhPDwcf/31F6ysrPDxxx/D29sbXl5eOlPTKycnBwcPHkRISAiuXbuGOnXqYNSoUfDz80Pfvn0hFArV3URC/hWFC9FoBQUFOH78OPbu3YtLly7ByMgIgwcPhre3NwYMGAAjIyN1N7FGJSUlYe/evQgJCcGjR4/QsGFDeHt7w9fXF+3bt9eZQCW6h8KFaByZTIYLFy4gPDwcx44dQ0lJCXr06AFvb2+MHDlSL+9MYYwhMjISoaGhCA8PR1ZWFlxdXfmDmk2aNFF3Ewkph8KFaATGGG7fvo3w8HDs378fmZmZcHZ25s+jNGvWTN1N1BhSqRRnz55FSEgIjh07BrFYjN69e8PPzw+jR49G3bp11d1EQihciHo9e/aMX5h/+PAh7OzsMHbsWHh7e6NDhw407fMv8vPzcejQIYSGhuLSpUswNjbGsGHD4Ofnh4EDB8LQ0FDdTSR6isKF1Lrc3FwcOXIEe/fuxbVr12Bqaorhw4dj3Lhx6Nu3LwwMDNTdRK2UkpKC8PBwhISE4P79+7CxscEnn3wCX19fndrwQLQDhQupFRKJBKdPn0Z4eDhOnjwJmUyGPn36wMfHB8OGDYOZmZm6m6gzGGP466+/EBISgrCwMLx48QKOjo7w9fWFr68vWrZsqe4mEj1A4UJqDGMMN2/exN69e3Ho0CHk5OSgffv2GDduHMaOHYuGDRuqu4k6Ty6X4+LFiwgNDcWhQ4dQVFSE7t27w9fXF2PGjIGVlZW6m0h0FIULUbnExESEh4cjLCwMT548QePGjTFu3DiMGzcOLi4u6m6e3ioqKsKxY8cQEhKCs2fPQigUYvDgwfDz88PgwYNhbGys7iYSHULhQlQiKyuLLxQZGRkJc3NzjBw5Ej4+PlSYUQOlp6dj3759CA0NRXR0NCwsLDBmzBj4+vpWuLAnIe9C4UKqrKSkBCdPnkRYWBjOnj0LABgwYADGjRuHIUOGVKlQJKl9cXFxfCHN5ORkNG/eHD4+PvD19YWzs7O6m0e0FIULqRSFQoE///wTYWFhOHLkCPLz89GpUyd4e3vj448/ho2NjbqbSKpIoVDg6tWrCAkJwcGDB/Hq1St4enrC19cX48aNg52dnbqbSLQIhQupkNjYWISFhSE8PBypqalo0aIFf8DR0dFR3c0jKlZaWoqTJ08iJCQEv//+OxQKBQYMGAA/Pz8MHz4cpqam6m4i0XAULuSt0tPTsX//foSFheHu3buwtLTERx99BB8fH3Tu3JnOTeiJ7OxsHDhwACEhIYiIiICZmRlGjx4NX19f9OnThwppkjeicCHlFBYW4vjx4wgLC8PFixdhYGCAwYMHY9y4cRg4cCDtKNJzjx8/5gtpJiYmolGjRvD29oafnx/atWun7uYRDULhQiCTyXDp0iWEhYXh2LFjKC4uRvfu3eHj44ORI0fC0tJS3U0kGoYxhlu3biEkJAT79u1DdnY23Nzc4OfnB29vbzRu3FjdTSRqRuGipxhjiImJ4QtFZmRkwMnJib9wy97eXt1NJFpCIpHgzJkzCA0NxbFjxyCRSNC3b1/4+vpi9OjRMDc3V3cTiRpQuOiZ5ORk/sKtuLg42NraYuzYsRg3bhw8PDxoHYVUy6tXr3Do0CGEhITg8uXLEIlEGD58OPz8/NC/f38qpKlHKFz0QF5eHo4ePYq9e/fizz//hEgkwrBhw+Dt7Y2+fftShyc1Ijk5GWFhYQgJCUFsbCxsbW0xbtw4+Pr6wtPTkx5kdByFi46SSCQ4e/YswsPDceLECUgkEvTp0wfe3t4YPnw4TVWQWqOcgg0NDUVYWBjS09PRunVr/qKzFi1aqLuJpAZQuOgQ5W2Fe/fuxcGDB5GTkwM3Nzd4e3tj7NixaNSokbqbSPScTCbDxYsXERISgsOHD6O4uBjvvfce/Pz88PHHH9PmER1C4aIDHj9+zF+49fjxYzRs2JAvFOnm5qbu5hHyRoWFhTh69ChCQ0Nx7tw5GBgYYMiQIfD19cWgQYNo27uWo3DRUtnZ2fj111+xd+9e3Lp1C2ZmZhg5ciS8vb3Rs2dPOthGtMqLFy+wb98+hISE4M6dO7C0tCxXSJPWZ7QPhYsWUZbkCA8Px+nTp8EYQ//+/eHt7Y0hQ4ZQSQ6iEx48eIDQ0FDs3bsXKSkpaNGiBXx8fODn5wcnJyd1N49UEIWLFmCMYebMmdi/fz9fTHDcuHEYM2YMbG1t1d08QmqEQqHAH3/8gdDQUBw8eBD5+fnw8vLC+vXr8d5776m7eeRfULhoCYlEAo7jIBQK6a4NoneUf6YYY+A4jqbJtACFCyGEEJUzUHcDtI2+ZzE9MRLqA9QHKoLCpQpOnTqFmJgYvfqQMcbQt29fdOnSRd1NIRrg5MmTuH37tt71gf79+6Nr167qbopWoHCpgrNnz6J+/fqVvgJWoVBo7XzxH3/8gZs3b1K4EADA6dOnUb9+fbRp0wZyuRxyuRxGRkZa+dmuqMuXLyMiIoLCpYIoXKqA4zj06NGjQjtWGGPIy8vDrl27EBERgdatW2Py5Mlo2rSpVnXEgoIC5ObmqrsZRENwHIdu3brh0qVLuHXrFurUqYM9e/bodFmh/Px86gOVQOFSw4qKijBjxgwcOnQIjDGcOnUKly9fRlhYGJo3b65VAUPI6ziOw82bN3HhwgXUr18fBQUFOh0upHJoT2sNYoxh7969OHLkCBhjfPXhmJgYLFmyBGKxuNbbI5VKUVpaqveLsqT6hEIh6tevD6DsIerVq1dqbhHRJBQuNSg3Nxc7duyAXC6Hubk5Nm3ahLZt2wIATpw4gUuXLtXaH3nGGM6fP4+RI0eif//+2L9/P+Ryea28NtFNHMehYcOGAACxWIzs7GyVvwZjDGlpafj9999x584dyGQylb8GqRkULjWEMYbLly/j4cOHAIDBgwfDx8cHCxcuhKGhIcRiMbZv3w6JRFIrbbl37x4mTpyICxcuICoqCjNmzKjVcCO6qUGDBgDKqh2/fPlSpT9bWeW7X79+GDZsGHr37o2VK1fWSp8h1UfhUkMUCgWOHj3K76Lx8/ODUCjEwIED4eHhAQC4evUq7t69W+N/4OVyOTZu3IgXL17w/1tBQQHWrFmDwsLCGn1totsaNmwIjuPAGEN6erpKP8uvXr3CrFmzEB8fD4VCgYKCAnz//ff8NDPRbBQuNSQzMxMREREAgJYtW6Jjx47gOA516tTBp59+CoFAgKKiIhw6dKjG25KUlITTp08DABwcHODl5QUAiIyMxNWrV6mjkiqrX78+X4E7PT1dZT+XMYbjx4/j1q1bAIAWLVrA2NgYpaWlWLduHfLy8lT2WqRmULjUAOU0lLKz9ezZE3Xr1gVQNk89cOBAfq769OnTNbq9kTGG06dP8/Phfn5+mDNnDgwMDCCVSvHrr79CoVDU2OsT3WZjY8NvVElLS1PZz5VIJNizZw/kcjlEIhG2bt2KoUOHAgDu3buHU6dO0UORhqNwqSFXr16FTCaDQCBAr169yv1bgwYN0Lt3bwDAkydPEB0dXWMdRSqV8qOWunXrYvDgwejevTtatmwJoOxwZGZmZo28NtF9FhYW/FUPGRkZKntQSUhIQGRkJACgU6dO6NGjB6ZMmQITExPI5XKEhoZCKpWq5LVIzaBwqQESiQTXr18HAFhaWsLd3b3ceRaBQIChQ4dCKBSW++NfE1JTU/HXX38BANq2bQsHBwdYWFigb9++AMouaarJcCO6zdzcnD/b8vLlS5X8wVfubCwoKAAADB8+HCYmJvDy8oK7uzsA4MaNG0hMTKz2a5GaQ+FSAzIyMvhdYk5OTvwUmBLHcfDy8uJ32ly5coXvSKrEGMPt27f5abeePXvCxMQEHMehf//+MDAwgFwux8WLF1X+2kQ/iEQi/t77nJwclJaWVvtnymQynD9/HkBZePXr1w8cx0EkEmHEiBEAgLy8PJw9e5YeijQYhUsNiI2NRU5ODgDAy8vrjXeB29nZ8QvrSUlJiI+Pr5G2XLt2DQqFAgYGBujRowdf26xDhw6ws7MDAFy/fh3FxcU18vpEtxkaGvIX1uXn56tk9+HLly8RExMDAGjdujUcHBwA/P96pXKkdPr06RqdGmOM4cmTJ3j69CmFWBVQuKgYYwxRUVGQy+UQCATo3LnzG79OIBDwT2QlJSU1smurtLQUUVFRAABbW1v+ACdQFm7t2rUDADx+/BjJyckqfW2iH14/pV9cXFztU/rKzTDKdcBu3bpBJBLx/+7o6Mh/ju/cuaPSTQR/V1paCn9/f/Tu3RsbN25UyahMn1C4qJhcLkd0dDSAsiG9i4vLG+uHcRyHrl278k9hly9fVvnp44yMDDx+/BhA2fScjY0N/29CoRDdu3cHABQWFuLOnTv0dEaq5PVT+llZWdX+edevX+c3w/Ts2bPcv5mYmKBPnz4AgOzsbNy6datGPreMMRw7dgxnzpxBcnIy9u3bR4c3K4nCRcVevXrFT3E1bdoUjRo1euvX2tvbw9HREUDZ9kpVn3COj4/nnyQ9PT35LaPA/6/7GBoagjHGb0AgpLKU4SKXy5GRkVGtnyWTyXDz5k0AQL169dC+fft/PJz16dMHhoaGUCgUuHz5crVe720yMzOxZs0aSCQSmJiYYMGCBfxxAlIxFC4qlpyczHcwV1dXfpvmm4hEIv5uiKysLNy/f19l7VAu5iun5zw9Pf/RSVu3bs2vu0RHR6OkpERlr0/0A8dxaNSoEX9KPy0trVojiZycHMTFxQEoO/D7ps0w7dq146fiIiIiUFRUVPVf4A2UBWeV/XH48OEYOHCgSl9DH1C4qBBjDA8ePODnZv++BflNunfvDoFAAJlMptJ1F4VCgTt37gAAzMzM0KZNm398jZWVFT9//fTp0xqdvya6q0GDBjAwKLu9o7qfocTERH69xd3dvdx6i5K1tTXat28PoGwzzNOnT6v1mn8nFotx7NgxMMZQt25dfPPNNzAyMlLpa+gDChcVU65dGBoavnFI/zrlri3lVs4bN26obF63sLCQ3w7dqFGjfzwBAoCBgQE6deoEoGynjypHTkR/2NjY8DsiqzNyYYwhJiaGv4pC+dn8O6FQyF/UV1hYqPJzWo8fP+Z3q3l4eLx13ZS8G4WLCkmlUty9exdA2eFJ5RbKd2nQoAGcnJwAAA8fPqz2nLVSWloaX37G2dkZderU+cfXcByHTp06QSgUQqFQIDIykhb1SaVZWFjAzMwMQFl9sape5cAY4zfDmJiYoF27du/cDGNkZATGGK5du6ayzy1jDOfOnUN+fj4AYNCgQTRqqSIKFxXKycnBkydPAADNmzcvtzvrbYyNjfnzLjk5OXjw4IFK2pKQkMDPRbdr1w4CwZvf6jZt2sDCwgJA2boLldQglVWnTh3Uq1cPQNlCeFVH32KxGPfu3QNQtnW+efPmb/1aZ2dnft0lKipKZeuFUqkUJ0+eBFC223PAgAE0aqkiChcVevr0KV8g0tXV9Y2HJ9+ka9euEAgEkMvluHHjRrWfwpRnBRQKBQQCAdzc3N76tXZ2dnydsUePHtXIhU9Et5mYmPAPUjk5OVU+kJuZmYlnz54BAFq1asVPF7+JlZUVXFxcAJTV53v+/HmVXvPvnj59itu3bwMoeyhT7uYklUfhoiLKP+jK+WJlDaR/w3Ec3Nzc+NHDrVu3qj16UCgU/PqJmZkZHB0d3/r0ZWJiwrc1MzOT6jWRSjMwMODX9AoKCqpcDj8xMZH/3nbt2vGbBN5EKBSiS5cu/Guq4l4kxhhu3LjBt2HAgAEwMTGp1s/UZxQuKqTcnWVsbAxXV9cKD6cbNmzIr888fPiw2gfRiouLkZCQAKDsvg1lDbO36dSpEziOg1Qqxe3bt2ndhVQKx3Fo3LgxgLJT7VX5/DLG8Ndff0Emk/EbXf7tNTt37gwDAwMoFAr+bEx1MMZw5coVMMZgZGSEXr160ZRYNVC4qEhpaSk/X2xtbY0WLVpU+HtNTEz42ymzsrL4YKiqzMxMfjHfwcHhjYv5SspzA8rzOLdu3aL7XUilNWnSBEDZmsXrN55WhrJ6t7GxcYV2aLVp04afOouKiqr2TsuioiK+zH+jRo3euH2fVByFi4pkZWVVeL74TTp37syPHqKioqo1enj69ClfZblt27b8TYFv06xZM76SwP379+nqY1IpHMehSZMm/EHK1NTUSn9+xWIxYmNjAZRtbW7WrNm/fo+dnR2/JpKQkFDtEX9ycjJ/ZqZdu3aV7sOkPAoXFXn06BFf2r59+/blSq38G47j0L59e347Z3VGD4wxxMbG8tMLykXPdzE3N4erqyuAsvtfqIglqayGDRvyayRVWVzPzs7mH86aN2/Or0G+i7GxMTw9Pfnvr05lccYY7ty5w++w7Nq1678+lJF3o3BRAcYY7t69y/9B79ixY6V/RtOmTdG0aVMAZaOH6lSXVW5nNjY2hpOT079OLwgEAn47dHFxMe7du0frLqRS7Ozs+MXvlJSUSn9+nj17xj+cubi4VPjhzMvLCxzHQSaTVfuc1vXr1/kD0MqZBFJ1FC4qoDyACACmpqZo27ZtpT+YderU4Uvgp6en8+dlKksikfAn8y0tLfm58HfhOA4eHh58EUtVLI4S/WJpaclX+E5NTa1UhW/laFu5ZqLsB//m7yP+mzdvVnnEX1JSwl9PYWdnB2dn5yr9HPL/KFxUoLCwkB8tNGrUiB+BVIaySjFQ9kGvagn8vLw8flqradOmFZpeAMruyVBe+kRFLEllmZmZwdraGkDZZV+VvftEuZhvaGhYqXIrqhrxp6enIykpCUDZRgHl70KqjsJFBZ49e4bU1FQAZYcnlU9wlaEcPSinFm7evFmlcElNTeVvwXRycqrwQU4bGxt+d0xiYqLKDqUR/WBsbMyfdcnJyeHLp1SEVCrlz2VZWFhUaqelmZkZX8QyLS2tykUsHzx4wAdTp06d3nnGhlQMhUs1KUvbK08ld+nS5a2lVv6Ng4MDX9Li7t27lT7pzBhDYmIi/9SoXKSvCAMDA/7ysPz8fJUXAyS6TSgU8ju8CgsL+crGFZGfn8+PGpo1a1ahsklKyvMuQNmIvyrntBhjiIyM5K+nUK7jkOqhcKkmxhgiIiLAGOPrhFX1g2lhYcGXwH99NFQZ9+/fB2MMQqEQbdq0qXBbOI5D9+7d+WKAf/zxB4ULqRR7e3sAZduKK/PZTUlJ4bcROzs7V3i0DYDfQKP8nqqsF8pkMty6dQtA2QVlFdlhSf4dhUs1FRYW8guBDRo0QKtWrar8s4RCIb/uUlhYWOldWwqFgj8rYGZmxtcMq6g2bdrwUxvXrl2rchkPon84jkPz5s3BcRwUCgWePXtW4c9ufHw8v8b3tkrI7+Lo6Mhfevf6LEJF5ebm8heUNW/e/J23x5KKo3CppoSEBH6e193dvVoHrziOg6enJ1/SQvk0VVFFRUV4/PgxgLKyL8oF+oqytrbmp8aePXuGv/76i0YvpMKaNWvGr1VUdLejsiafcrTt5uZW6XD5exHLyo74Hz9+zF8x7u7uTvXEVITCpRoYYzh58iSKi4vBcRz69u1b5fUWJWdnZ1hZWQGofEmLypR9eROO4/Dhhx9CIBBAIpHgzJkzlfp+ot8aNmzIlxF68uRJhbYFy+VyfqeYubl5lUb+QqGQX3fJz8+vVBFL5R0yyoKzdL5FdShcqiEvLw/Hjh0DUHb/RL9+/ar9wbS1teU72OtPVBWRlJRUqbIvf8dxHLp06cJvKjh9+jR/sI2Qf2NlZcU/GCUnJ1fowaioqAiPHj0CULaNX/nZqwzl51Y54r9+/XqFv/f1c10ikahCV5OTiqFwqYZr167xHaNPnz4Vqof0b4yMjPiSFjk5OYiPj6/QU5jyIJpcLgfHcZXaKfa6hg0bonfv3gDKStocP36cpsZIhZiZmfHrFenp6RXajvzixQt+tN26dWt+5FNZLi4u/NmUmzdvVvicTVFREX97bIMGDSq9TknejsKliuRyOQ4cOACZTAZDQ0OMHTu22lNiSp07d4ZAIIBMJqvU7pfXy7686w6XdxEIBBg/fjxMTU0hl8uxbNkyHD16lColk39lYGDAn1HJy8ur0Kj74cOHfD0vNze3KvchOzs7fqdlQkIC0tLSKvR9z58/5w8du7q68jdqkuqjcKmivLw8PHnyBEKhEE5OTujatatKhtPKEvh169YFANy4caNCpTTEYjFfuM/a2rpCZV/e9vqdO3fGuHHjIBAIkJGRgRUrVlS74izRfRzH8VO6JSUl/1oAlTGGmJgY/sbUDh06VLkPGRoa8ptR8vLyKnROS1kTUDmV7OXlpbIHRELhUmVWVlY4fPgw/vvf/2LevHkqfeJp1KgRPzyPi4ur0B/27OxsvjPb29tXa9eaoaEhVq5ciS+//BK2trb47rvvKr3zjOgnZaFUhUKBhISEd/6BVygU/JSUmZlZtet5vffee/y6yx9//FGh71FeK25gYECHJ1WMwqWKOI6Dra0tRo8ejY8++kilH0qRSMSvu2RmZlaolHhycnKVqsq+CcdxqFevHtatW4fjx49j0KBB1OnIv+I4Dg4ODjAyMgKAf/3cFhUV8V/TsGFD/oxVVV/bzc2N3xBw7do1frrtbcRiMb/d39rami4HUzEKl2riOK5G/vB2796dvzxMWQr8bRhjuH//fqWryr4Lx3EQiUTo0KEDTRWQCmvcuDFfLPXhw4eQSqVv/drU1FR+baRNmzZ8deOqsrW15W90TUxMRGJi4ju/PjU1lb/1tU2bNvxBTKIa9FdDA3EcB3d3d36q7dq1a+/spAAQExMDoGwx39XVlUYaRC1ev+bh6dOnb90xptzdqLz11N3dvdoPMUKhEO+//z6AslHR5cuX3/pQprwcTFmFonv37tUa7ZN/onDRUE2aNIGTkxMAIDY29p07b8RiMV9V1traGs2bN6+NJhLyD8bGxvzaSVZW1jura0dFRUGhUEAoFKJjx47VfiDiOA69evXiR0Bnz55952aYq1evQqFQlCvaSlSHwkVDmZiYoEuXLgDKOum77nfJysriS9C0aNGC7v4mavP6GauSkhI8fPjwjZ9bmUzG1+SrV6+eyi7ncnR05NdOoqKikJKS8savKykp4Q9b2tjYVKnsDHk3ChcN1qtXLwiFQsjlcly+fPmtX5eYmMjf4dKhQwd+QZWQ2qZcWBcKhWCM8aVd/u71O+9btmxZrcX814lEInzwwQf8a1y8ePGN4ZaYmMi/focOHWi9pQZQuGgojuPQoUMHfvfL1atX+fnp1ylrI8lkMr7wJSHq5OTkxJ/Tun37NuRy+T++5uHDh/ydLx07dlRpschBgwbB1NQUjDEcPnz4H+uVjDFcvnyZ70/9+/eny8FqAIWLBrOzs0PHjh0BlD1pKXe2vO716sl16tSh4T1Ru4YNG/J3u8TFxf2jPp3yDiSJRAKO49CtWzeVvbZy5KTcMXnjxg2+RJOSTCbD2bNnAZT1md69e1OfqQEULhpMKBSif//+4DgOxcXFuHDhwj+G+Lm5ufzUQ7NmzVRS34yQ6lBuYQeAjIyMN/5xv3LlCoCySsienp4qPyc2evRocByHvLw87N+/v1y/efbsGV9WqU2bNvzGGaJaFC4aTLn7RXlu4NSpU/8oyJeQkIAXL14AKJteqGyZfUJUTTka4TgOYrGYv6lVKSMjgz+Z7+TkxI9yVPn6I0eO5KeUw8LC+D7CGMOxY8f4NcqhQ4dWuVgmeTcKFw3XvHlz/mDYvXv3+JsmAfDXEYvFYnAch549e9Lwnqgdx3Hw8vLiH3QuXbrEbwlmjOHWrVv8ekuPHj0gEolU3gZ7e3uMGjUKQNndMrt27YJCoUB2djZ++eUXMMZgYWGBkSNHUp+pIRQuGs7Q0BAjRowAx3EoKirC4cOH+adAsViMCxcuAAAsLCzooiOiMRwcHPjppujo6HIjh99++w1yuRyGhoYYMGBAjby+QCDApEmTYGtrC8YYNm/ejPPnz2Pz5s38A9qgQYPQunXrGnl9QuGi8TiOw8CBA/mtmkeOHEFGRgaAsikx5fRCu3btaL2FaAxTU1P069cPQFl9POWW4NTUVJw/fx5A2ahc1estShzHoW3btpg6dSoEAgGysrIwZswYrF27FgqFAtbW1pg1axbtEqtBFC5aoHHjxhgyZAiAspIa+/btg1wux8GDB1FQUACO4zB06FA630I0yrBhwyASiaBQKLB3716Ulpbi0KFDfD2xoUOH1uiBX4FAAH9/f4waNQocxyE/Px8SiQQikQiLFy+mWydrGMW2FuA4DuPHj8fBgweRm5uLTZs2oXHjxggLCwNQdoPe4MGDqaMQjaE8p+Xl5YUrV67g6tWr2LZtG3788UcwxlCvXj34+PjU+GfW3Nwc27dvh5ubG06cOAEzMzN89dVX+Oijj6ggaw2j/3e1gHLvvp+fHziOQ1paGiZMmIC0tDRwHAdvb2+aEiMaRyQSYdq0aTAyMoJYLMa8efPw+PFjAMCYMWNq5UwWx3GwtLTE4sWLceXKFZw+fRqffPIJFamsBRQuWkIgEGD27Nno2bMngLLFfKDs9rxp06bRqIVoHI7jMHjwYPj5+fHXdgNA+/btsWDBAgiFwlpti4mJCQwNDamv1BKaFtMSysvJfv75Z2zZsgWRkZFwc3PDjBkzUL9+feowRCMZGxvjhx9+QMuWLXHhwgU4OTkhICAA9vb29JnVcRQuVcAYQ1xcnNoW0IcOHYrBgwdDIBAgPT0d6enpNf6aT5484Q9zEqK8j8XY2LhCX9+3b1/07NkTQqEQ2dnZyM7OruEWqh71gcqhcKkCT09PXL9+nb+gSx8oFAp07dpV3c0gGqJTp064du0a7ty5o+6m1BqFQqHSOmi6jmPvuj+X/IO+/99FUxmE+gD1gYqgcCGEEKJytFuMEEKIylG4EEIIUTla0NcSjDEwxsBxHM35Er31+iw+9QPNRiMXLRETEwNTU1O92qFGyN/duXMHAoFAr3apaSsKF0IIISpH4UIIIUTlKFwIIYSoHIULIYQQlaNwIYQQonIULoQQQlSOwoUQQojKUbgQQghROQoXQgghKkfhQgghROUoXAghhKgchQshhBCVo3AhhBCichQuhBBCVI7ChRBCiMpRuBBCCFE5ChdCCCEqR+FCCCFE5ShcCCGEqByFCyGEEJWjcCGEEKJyFC6EEEJUjsKFEEKIylG4EEIIUTkKF0IIISpH4UIIIUTlKFy0AGMMubm55f6bEH2j/PwDoH6gBShcNFheXh6Cg4Ph6OiI999/HxKJBO+//z4cHR0RHByMvLw8dTeRkBr3ej/o168fAKBfv37UDzQcxyj+NdKZM2cwevRoFBcXA0C5pzSO4wAApqamOHToEAYOHKiWNhJS06gfaC8KFw105swZDB48GIwxKBSKt36dQCAAx3E4efIkdSyic6gfaDcKFw2Tl5eHJk2aoKSk5J0dSkkgEEAkEuH58+ewsLCo+QYSUguoH2g/WnPRML/88guKi4sr1KEAQKFQoLi4GHv27KnhlhFSe6gfaD8auWgQxhgcHR2RlJRUqZ0wHMehZcuWePToET8PTYi2on6gGyhcNEhWVhZsbW2r9f3W1tYqbBEhtY/6gW6gaTENUlhYWK3vLygoUFFLCFGfnJycan0/9QPNYKDuBpD/Z2ZmVq3vX79+PXr06AFPT0+0aNGCpgaIVigpKUF8fDwePHiA2NhY3Lt3r1o/z9zcXEUtI9VB02IaRDnX/Pjx40p/r5WVFfr27YvY2FjI5XJYWVmhY8eO6NSpEzw9PdG6dWsIhcIaaDUhlZOfn4+4uDg8ePAADx484NdW6tWrBxcXF7Rt2xZTp05FcnJypX4urbloFhq5aBCO4zBkyBAEBwdX+vuWLl2KGTNmoKioCDExMYiKikJUVBR++OEHSCQSmJmZwcPDA56enujUqRNcXV1haGhYQ78JIf8vOzubH5U8ePAAKSkpAAA7Ozu0bdsWH3zwAVxcXNCwYUM+FB4/foxZs2ZVakGfMYYZM2ZQsGgIGrloCJlMhsDAQGzbtg0pKSmQy+UV6lj/tr9fIpHg3r17fNjcvn0bRUVFMDY2Rvv27eHp6QlPT0+4u7tDJBLVwG9G9AljDC9evOCDJDY2FhkZGQCAJk2a8COTtm3bvnPRvrLnXADAwMAAkZGR6NChgyp+FVJNFC4aIDk5GTNnzsSDBw8wZ84cNG3aFEOHDq3wyeTff/8dAwYMqNBryeVyxMfH82ETFRWF3NxcCIVCuLi48CMbDw8P1KtXT1W/ItFRjDE8e/as3MgkLy+Pn6Jq27YtHyh169at1M+uzAl9ABg1ahQ4joOfnx+GDBlCIxg1o3BRsyNHjuDbb7+FjY0NgoOD0a5dOwAVr6l0+PDhCgfLmzDGkJSUhKioKERGRiIqKgrp6ekAACcnJ3Tq1AkdO3aEp6cn7Ozsqvw6RDfI5XIkJibyQRIXF4eioiIYGBjA0dGRDxJnZ2eYmppW+/Uq0w/69OmDsLAwHD9+HO7u7pg2bRqd1lcjChc1KSwsxLfffotjx45h5MiRWL58OerUqVPua/Ly8rBnzx5s3Lix3CK/g4MDZsyYgc8++0zlowvGGNLS0hAZGYno6GhERkbi6dOnAIBmzZrxIxtPT080adKEng51nFgsRkJCAj8yefjwIcRiMUxMTODs7MyPTBwdHWFkZFQjbahsP4iJicGmTZvAGMP06dPh7u5eI+0i70bhogYxMTHw9/dHTk4OVq5ciWHDhr3z6xljyMnJQUFBAczNzWFlZVWrf9SzsrL4oImKisLDhw/BGIOdnR0fNJ6ennBwcOCnKIh2Kioq4ndyxcbGIjExEXK5HGZmZuWmuFq2bFnruw8r0w/y8vKwZcsW3L59G0OHDoWPjw9tYKllFC61SC6XY8eOHdiwYQPc3NwQFBSEpk2bqrtZlZafn4/bt2/zazb37t2DXC5HvXr1+KDx9PRE27ZtafuzhsvLy+OnuB48eIBnz56BMQYrKys+SFxcXNC0aVOtG6UyxnDy5EmEhISgadOmCAgIQKNGjdTdLL1B4VJLMjIyEBAQgJs3b2Ly5MmYOXMmDAx0Yyd4aWlpue3PMTExKC0thUgkgru7Oz+6cXNzg4mJibqbq7cYY8jMzCy3+J6WlgYAaNCgAVxcXPhAqV+/vtaFyds8efIEgYGByM7Oxpdffom+ffvqzO+myShcasH58+cxb948GBsbIzAwEF26dFF3k2qUVCpFbGwsP40WHR2NgoICGBoaol27duW2P1e3KgF5O8YYUlNT+VFJbGwssrKyAAD29vblprmsrKzU3NqaJRaL8dNPP+H8+fPo1q0bJk2a9I81TqJaFC41qLS0FGvWrEFoaCj69euHdevW6eXuFblcjkePHvEjm8jISGRnZ0MgEKBNmzZ82HTs2FHn/8jVJLlcjqdPn5Y7Y5Kfnw+BQAAHBwd+ZNKmTRu9DfWIiAhs3boVIpEIs2bNgrOzs7qbpLMoXGpIQkICZs6ciWfPnmHRokXw9vamofj/MMaQnJzMj2yioqLw/PlzAGU7gF5ft2nYsKGaW6u5pFIpEhMT+SCJi4tDSUkJjIyM4OTkxI9MWrduTdORr8nKykJQUBDi4+MxZswYjB49mtYGawCFi4oxxhAaGorVq1ejefPmCA4OhpOTk7qbpfFevHiB6OhoPmwSExMBAI0bN+aDplOnTrC3t9fbkC4tLUV8fDw/MklISIBUKoVIJEKbNm34Ka5WrVrRzqh/IZfLcejQIRw4cADOzs7w9/eHjY2NupulUyhcVCg3Nxfz58/H+fPn4efnh/nz59MTYxXl5ubyYRMZGYm4uDgoFApYW1vzQdOxY0c4OTnp7FNnQUEBYmNj+f88fvwYCoUCdevWLbeTq3nz5rQFvIri4uIQFBSEkpISTJ48GV27dlV3k3QGhYuKREREYPbs2RCLxVi3bh369eun7ibplMLCQsTExPCHO+/evQupVApzc3N07NiRrwDt4uKitU/t2dnZ5dZLlFWBbWxsyu3katy4sd6O3mpCYWEhtm/fjuvXr6Nfv3744osvYGxsrO5maT0Kl2qSyWTYsGEDtm/fji5dumD9+vWoX7++upul88RiMf766y/+cOft27dRUlICExMTtG/fnt/+3KFDB40cPTLGkJ6eXu6MibLAY+PGjcuFSXVuZSQVwxjDxYsXsXv3blhbWyMgIAAtWrRQd7O0GoVLNSQnJ8Pf3x/3799HQEAAJkyYoLNTNJpOLpcjLi6u3CaBV69eQSgUwtXVlQ8bDw+PShdQVAXlJobXz5jk5uaC4zi0aNGCn+Jq06aNXu4o1BSpqanYsGEDUlJS4Ofnh8GDB9MosYooXKro6NGj+Pbbb2FlZYXg4GC0b99e3U0ir1EoFEhKSioXNhkZGeA4Dq1bty53kVpNLOTK5XIkJSXxo5K4uDgUFhZCKBT+o8AjnbfQLFKpFHv37sVvv/0GDw8PTJ06lQK/CihcKqmwsBBLly7F0aNHMWLECCxfvlxvzwxoE8YYnj9/zk+jRUZG8msa9vb25WqkVWVNQyKRICEhgR+VxMfHQywWw9jYGK1bt+anuRwdHWk+X0vcuXMHmzZtAsdxmD59Ot0TU0kULpVw9+5dzJw5Ezk5OVixYgVGjBih7iaRasjMzCx3r01CQgIYY2jQoEG5szYODg7/CJvi4uJy974nJCRALpejTp065U6+Ozg40FSpFsvLy8OmTZsQExNDBTAricKlAhQKBXbs2IHAwEC4uroiKCgIzZo1U3eziIopC3Iqp9Lu378PuVwOS0tLuLm5wc7ODkZGRsjJyeELPFpaWvJh4uLigmbNmtEcvY5hjOHEiRMICQmBvb09Zs2aRQUwK4DC5V9kZGRg9uzZuHHjBiZNmgR/f3+dKThJ3i4zMxN37tzB2bNncevWLTx58gT5+fngOA7m5uZwdXVF79690bdvX7Rr146muvRAUlISNmzYgOzsbHz11Vfo06cPPUi8A4XLO1y4cAFz586FkZERAgMD6YCVjlIWeHz9jElmZiaAsgvSlCOTVq1aIT09nT9rEx0djcLCQhgZGfEFOTt16oQOHTrQIr2OKi0txe7du3Hx4kUqgPkvKFzeoLS0FGvXrkVISAj69euHtWvXwtLSUt3NIiqiUCjw9OnTctuCXy/w2LZtW/4/5ubmb/05crkcDx8+LHeRWk5ODoRCIdq0acNXEejYsSN9fnTM9evXsXXrVtSpUwezZs1C69at1d0kjUPh8jePHj3CjBkz8PTpUyxatAg+Pj409NVyMpmML/Co3MlVXFwMQ0NDODk58YvvrVu3hkgkqvLrMMbw9OnTcldEK+9LadWqVbkroumgrfbLzMzEhg0b8OjRI74AJpXh+X8ULv/DGENYWBhWrVqFZs2aITg4mJ5GtFRpaSkePnxYrsCjRCKBSCSCs7Mzv/ju4OBQY/e+K7148YKvjxYVFYWkpCQAQJMmTcqFDW0E0E5yuRwHDx7Er7/+SgUw/4bCBWXbDefPn49z587B19cXCxYs0MiSIeTNCgsLy937/vjxY8jlcpibm/+jwKO6twVnZ2eX25EWHx8PhUIBGxubcmdtHB0d6SlYi7xeAHPKlCk6fyFgReh9uNy8eROzZs2igpNaJDc3t1xNruTkZDDGYG1tXa4mV5MmTTR+NFBQUIA7d+7wZ23++usvyGQy1K1bFx4eHnzguLi40C5FDVdYWIht27YhIiIC/fv3x/jx4/V6F6HehotMJkNwcDC2bt2Kzp07IzAwkObBNRBjDC9fviy3+P7ixQsAQKNGjcqNTGxtbTU+TP5NaWkp7t27x49s7ty5wxfkdHd350c27du3p9G1BmKM4cKFC9i9ezdsbW0REBCA5s2bq7tZaqGX4ZKSkgJ/f3/cu3cPs2bNwtdff6326RJShjGGlJSUctuCs7OzwXEc7O3ty13Vqw87sGQyGWJjY8tVEsjPz4eBgQHc3Nz4sPHw8HjnzjZSu54/f44NGzbg+fPn+PTTTzFo0CCtf/CpLL0Ll+PHj2Px4sWwsrJCUFAQ1QtSM2WBR+WFWA8ePOALPLZq1apctWA6T1C2jToxMZEPmsjISGRmZoLjODg7O5crW2Ntba3u5uo1qVSK0NBQnDhxAh4eHpg2bRrq1aun7mbVGr0Jl6KiIixduhRHjhzB8OHDsWLFCio4qQYSiQSPHj3iRyXx8fEoLS2FkZERnJ2d+fMlTk5ONO1TAcqRnjJooqOj+YKcLVq0KHdFNJUsUY/bt29j06ZNEAgEmDFjht5UUNeLcImPj8eUKVOQmZmJFStWYOTIkepukt6JiYnB/v378ejRI8hkMpiamv6jwCMtWKtGRkYGoqKi+LM2jx49AgA0bNgQXbp0werVq/Vuikbd8vLysHHjRty9exfDhw+Hn5+fzr8HehEuMpkMOTk5sLS0pIqmapKfn4/09HTUrVsX5ubmMDU11fnOpSnkcjlKS0tRWloKhUJBN1uqCWMML168gFQqhb29vbqbU+P0IlwIIYTUrlqdh9D3HNOEJ3V6D+g9UDd6D9SvNt6DWp/kvnTpEh48eKARH7DawhhD9+7d4eHhoe6mAABfhkTf3oP27dvD2dlZ3U0BAFy+fBmxsbF69x5069YN7u7u6m4KACA6Olpv+0FtlLaq9XC5cuUKbG1t0apVq9p+6VollUqRkZGBxo0b4+bNm7hz547GhMudO3dgYWGBJk2aIDMzE1ZWVjq/mH7v3j08fPhQY8Llzz//hI2NDRwcHNTdlBrFGEN6ejrq16+PqKgoxMTEaEy4vN4PdJlEIkFBQQGsra1x//59PHz4UDfDheM4eHl5wcvLq7ZfulYwxlBaWort27fj6NGj+P777+Hq6opXr16pu2nltG3bFlKpFKdOnUL//v0xbNgwnd7sUFxcjMLCQnU3g8dxHDp16oROnTqpuyk1gjGG/Px8hISE4PDhw1i9ejVcXFw0rh8odyvqKqlUin379uHWrVuYNGkSWrZsiaKiolp5baqMp0KMMeTl5WHZsmXYunUrUlJSsGzZMmRnZ6u7af+QmZmJ7du3IyMjA/v378eJEycgk8nU3SyiAxQKBe7du4cpU6bgxx9/RFpaGgIDA5Gfn6/upukVmUyG3377DSdOnEBaWhp27NhRq+FO4aIijDEkJyfD398fhw4d4osP+vn5aeSpXGtrawwdOhRGRkaQSqUIDw/HyZMnKWBIlTHGIBaLsX//fkyaNAnR0dFQKBSwsrLCBx98oNdFHGubXC7HhQsXcPDgQchkMhgbG2PEiBG1+reIwkUFGGO4ffs2Jk2ahKtXr4IxhgYNGuCHH37A2LFjNXI9QygU4oMPPoCPjw8fMGFhYThz5gzkcrm6m0e0DGMMubm5WLlyJVatWoWcnBwIBAJ07twZO3bswOeff07hUksUCgVu3LiBkJAQSCQSGBoawtvbGz179qzVzQua91dPy8hkMpw4cQKrV6/mp79cXFywatUquLq6avROFKFQiEGDBkGhUCAsLAwSiQQhISEQCoXo378/FfMkFcIYQ2JiIpYvX47o6GgwxiASifDZZ5/hiy++gLm5uUb3A13CGMO9e/ewa9culJSUQCgUYuTIkRg4cGCt92cKlypijKGkpAQ7duzg30iBQIA+ffpg6dKlaNSokVZ0KAMDAwwZMgRyuRz79++HWCzGL7/8AgMDA/Tp04cChryTQqFAREQEli9fztc0a9SoEebPn4++fftCKBRqRT/QBYwxPH78GFu3bkV+fj44jsOAAQMwYsQItWzWoXCpAsYYMjMzsWrVKpw6dQpyuRxGRkbw9fXFjBkzYGZmplUdysDAAMOGDYNMJsPBgwdRWlqKn3/+GYaGhujRowfdiEjeSCaT4fjx4/j++++Rl5cHAHB3d8eyZcvg5OSkVX1A2zHGkJqais2bN/NVsrt37w5vb+8av8r7bShcKokxhvj4eCxZsgR37twBAFhYWGD27Nn46KOP1PZGVpehoSFGjRoFqVSKo0ePori4GDt37oSBgQG6du1KAUN4yoX7n376CTt27EBpaSm/hjd//nzY2NhQsNQixhiys7OxZcsWpKSkAAA6dOiAL7/8EiKRSG3vBYVLJSgUCly5cgVLly5FamoqAMDe3h4rVqxAt27dtP4PsKGhIcaMGcOvIxUVFWHbtm0QCoXw8vLS+t+PVB9jDEVFRQgMDMSBAwcgk8lgZGSEzz//HBMnTqSCpLWMMYaCggJs27YNCQkJAAAnJydMmjRJ7WtdFC4VwBiDVCrFwYMH8cMPP/DzmZ06dcJ3330HBwcHnelQhoaGGDduHGQyGU6dOoXCwkJs3boVhoaG8PDw0Jnfk1QeYww5OTlYuXIlzpw5A4VCATMzM/j7+2PMmDFaO2rXVsoD2z///DNiYmIAAE2aNMHUqVNhbW2t9r5Kj6L/gjGG4uJiBAUFYeXKlfwVsyNHjsTmzZt1KliAspPjyvWj/v37g+M45OfnY/Pmzbh7967eF/zTV8o5/Tlz5uD06dNQKBSwsbHBqlWrMG7cOAoWNZBIJAgLC8Off/4JxhhsbGwwbdo0NG7cWCP+JlG4vINy7/7SpUuxa9cuSCQSiEQiTJ06FStWrICVlZVGvImqxnEcjI2N8dlnn6FPnz7gOA55eXnYtGkTHjx4QAGjZ5S7kGbNmoWIiAgwxtC0aVOsX78eAwYMoB2FaiCTyXDs2DGcOXMGjDHUrVsXkydPRqtWrTTmbxKFy1soL/aZPXs2jh49CrlcDgsLCyxfvhxTpkxR60JZbeA4DiYmJvjiiy/Qo0cPcByHnJwcbNy4EfHx8RQweoIxhvv372PmzJm4d+8eAKB169YIDg6Gl5eXTvcBTSWXy3Hu3DkcPnwYcrkcIpEIX331Fdq1a6dR7weFyxson9SmTZuGP/74A4wxNG7cGIGBgRg5cqRGnrivCRzHwdTUFBMmTEDXrl3BcRwyMzMRHByMxMREChgdxxhDZGQk/P398fjxYwCAh4cHgoOD0aZNG436Q6YvFAoFrl+/jr1790IqlfJT2F26dNG4DTea1RoNoDzhOm3aNNy9exdA2e6LLVu2oGfPnhr3BtY0juNQp04dTJw4EZ06dQLHccjIyMCGDRuQlJREAaOjFAoFrl69ijlz5iA1NRUcx6FHjx4IDAyEvb09BYsaMMZw584d7N69mz99P3r0aPTr108jpyb16y/lv2CM4ebNm5g+fToePXoEAOjYsSN+/PFHjS/lUpM4joO5uTkmT57M30nz4sULBAUFITk5mQJGxygUCly8eBHz5s3Dy5cvwXEcBg4ciHXr1qF+/fp62w/UiTGGuLg4bNu2DQUFBRAIBBg0aBCGDRumsTMpFC7/o1Ao8McffyAgIADPnz8Hx3Ho2bMngoOD0bx5c73vUBzHoV69epgyZQrat28PAHj+/Dk2bNiA58+fU8DoCIVCgbNnz2Lx4sV88cnhw4djxYoVsLS01Pt+oA6MMSQlJWHz5s3IyckBx3Ho1asXxo4dq9G79ChcUNahzp8/j2+++QYZGRngOA4ffPABfvjhBzRo0IA61P9wHAdLS0tMmzYNLi4uAIBnz54hKCgIL168oIDRcgqFAqdPn8bSpUuRl5cHoVCIjz76CEuWLFH7gTx9pdwCvmnTJv5vk5eXFz7//HOYmJiou3nvpPfholAocOrUKSxYsADZ2dkQCAQYOXIkVq9erbNbjauD4zhYW1tjxowZ/JXBSUlJCAoKQkZGBgWMllIGy/Lly/Hq1SsIhUJ88sknmD9/Pp26VxPGGF6+fIlNmzbxZV3atWuHr7/+GnXq1NH490Svw0WhUODkyZNYsmQJ/6Q2duxYLF26lJ7U3oHjONja2mLmzJlwdHQEADx69AjBwcHIysqigNEyCoUC586d44PFwMAAPj4+mDNnjs5vuddUymoImzdvRmJiIoCyLeCTJ09GvXr1tOI90dtwUQbLt99+yz+p+fj4YOHChVrxVKBuHMehfv36mDlzJlq2bAkAiI+Px6ZNm5Cbm0sBoyUYY7h8+TKWLl3K94Nx48Zh1qxZFCxqwhhDfn4+tm7ditjYWABA8+bNMW3aNK0qCqqX4fJ6sCjLuXz66aeYO3cudahK4DgOjRo1gr+/P5o1awYAuHfvHrZs2YL8/HwKGA3HGMO1a9fw7bff8iP3Tz75BAEBARo/n6+rGGMoLCzEtm3b+HphjRs3xsyZM9GwYUOt+tukd+GiUChw5swZLF26tFywzJ49GyYmJlr15mkCjuPQpEkT+Pv7o1GjRgCAO3fuYPv27SgqKqKA0VDKMxOLFy9GVlYWBAIBRo0axQcL9YPap6xjuGvXLkRGRoIxhvr162PGjBlo2rSp1r0nehUuCoUCFy5cwJIlS/gpAF9fX+pQ1cRxHJo3bw5/f3/Ur18fjDHcuHEDu3fvRmlpKQWMhlHeSbRgwQKkp6eD4zgMHjwY8+bNo8V7NXm9wvG1a9fAGIO1tTWmT5+utcVx9SZcGGP4888/sWjRIuTm5kIoFMLb2xtz5syhYFEBjuPQqlUrTJ8+HVZWVmCM4cqVKwgJCYFEIlF388j/MMbw9OlTzJ8/H8+ePQMA9O3bF4sWLaK1RjVRXr62Z88eXL58GYwxWFhYYOrUqXB2dtba90QvwoUxhlu3bpXbbvzRRx/hm2++oWBRIY7j0LZtW0ydOhV169YFYwxnz57FwYMHIZVK1d08vafc2rpo0SI8fPgQANClSxcsW7ZMa3Yg6RrGGCQSCfbu3Yvz58/zFY6nTJmicYUoK0vnw0VZK2zu3Ln8IaThw4djwYIFNAVQAziOQ4cOHfD111/D1NQUcrkcx48fx8mTJyGXy9XdPL3FGMOrV6+wdOlS3L59GwDg5uaGlStXatUOJF2iDJbw8PByl69NmjQJ7u7uWv+e6HS4MMaQkJCAOXPm8CVdPvjgAyxZsgRmZmZa/+ZpKoFAgC5dumD8+PEwNjaGVCrFvn37cOnSJSgUCnU3Ty8VFxdj7dq1uHLlCgCgZcuWWLNmjcZcLKWPpFIpDhw4gN9//x1yubxcgVhdKJCr/b/BWzDGkJKSgm+++YYvF96zZ08sX74cdevWpQ5Vw4RCIfr06YNPPvkEhoaGEIvF+Pnnn3Hz5k1a4K9lEokEW7ZswW+//QbGGBo0aIDVq1dr7UKxLpBIJDhw4AB+++03yOVymJqa4quvvtLI0vlVpRu/xd8o55bnzZuHBw8eAAC8vLywZs0aKulSi4RCIQYPHoxhw4ZBKBSiuLgYO3bswL179yhgaolMJkNISAhCQkL4C+9WrFiB9u3bUz9QE4lEgoMHD+L48eP8ZV9ffvklunfvrjPBAuhguCjnlr/99lvcunULAODi4oK1a9fCzs6OOlQtMzQ0xMcff4z+/fvz1yVv2bKFLhurBQqFAr/99hs2b94MqVQKkUiE+fPn47333qN+oCbKYDl27Fi5YOnRo4dG3slSHToVLowxlJSUYPXq1bhw4QIAoEWLFvjPf/6DZs2aUYdSEyMjI/j5+aF79+7gOI4vxpeamkoBU0OUW+/Xrl2LkpISGBoaYvr06RgyZIhOPR1rE7FYjAMHDpQLli+++AI9e/bUuWABdCxcJBIJNm3ahKNHj/Jzy+vWrYOTkxMFixpxHMff892hQwcAQEpKCjZt2kSFLmsAYwx3794tVy/s008/ha+vr8ZeLKXLlOdY9u3b94+psF69eulksAA6FC5yuRyhoaH4+eef+bnllStXwsPDg4JFA3Acx+/fd3JyAgAkJCTgxx9/xKtXryhgVIQxhidPnmDRokX86fuhQ4diypQpGn2xlK5SBktYWBi/HV8ZLLo6YlHSiXBRFqIMDg6GVCqFqakpFi5ciF69elGwaBDlXTDTp09H06ZNAQB3797Frl27UFxcTAFTTYwxZGRkYOHChfwOyV69evFlXUjtUpZ02bNnD7/d2NTUFBMmTND5YAF0IFwYY4iIiMB3332HoqIifm55+PDhNLesgTiO46u82tnZgTGG69evY8+ePVQmphqUG1mWLVvGV9Pt0KEDli1bBgsLC3rIqmXKIpQ//fQTzp49yx+QnDhxok4u3r+JVv/1ZYwhLi4OixYt4u/79vHxwWeffUZzyxqM4zi0bNkSU6dOhYWFBRhjuHDhAvbv309lYqqopKTkH4ckV61ahfr161Ow1DJl2fwdO3bg0qVLYIzB3NwckydPRrdu3fTmoVdrf0vl3dLz589HSkoKOI7Dhx9+CH9/f5pb1gIcx8HNzQ0TJ05EnTp1IJfL8dtvv+HYsWOQyWTqbp5WedMhyVWrVqFly5YULLVMOYLcsmULX93YwsIC06dPh5eXl94EC6Cl4cIYQ25uLhYtWsQfkuzcuTO+/fZbKuuiRTiOg5eXF7744guYmJhAJpPhwIEDOH36NNUhqyDlIck9e/ZALpejXr16WL58OTp06ED9oJYxxpCVlYWgoCD+PhYbGxvMnDkTHh4eehUsAKCVc0fFxcVYvXo1rl27BgBwdnbG6tWrYW1tTR1KywgEAvTq1QslJSX45ZdfIJVKERoaCmNjY/Tt21cv5qar6m2HJHv06EH9oJYxxpCWlobNmzcjISEBANCgQQNMnz4drVu31sv3Q+vCRSKRYPPmzTh+/DgYY2jUqBHWrl1LhyS1mFAoxMCBA1FaWop9+/ZBLBbjv//9L4yMjNCjRw+9e+KrCOV9OX8/JDl06FD6/6uWKbd/b9q0CcnJyQCAZs2aYfr06WjRooXe/l3SqnBRnmX573//C7lcDktLS6xcuRKurq56+wbqCgMDAwwbNgxisRiHDx9GcXExdu7cCUNDQ50q5qcKyiuKly1bxh+S/Pzzz+mQpBooNxVt3rwZGRkZAABHR0dMmzZN7ytOa80nUXmWJSgoCBKJBCKRCAsWLKApAB1iaGiIjz76CBKJBCdOnEBRURG2bdsGoVCoM2XIq4sxhsTERCxatIi/n2jEiBGYPHkybWSpZQqFAtHR0di+fTtyc3PBcRxcXV0xdepUuiMHWrKgzxjD1atXy51lmTZtGp1l0UFGRkYYN24cBg4cCIFAgIKCAvz444+Ijo7W+7tglPP6CxYswJMnTwAAffr0wdy5cyESidTcOv0il8tx5coVbNq0iQ+Wzp07w9/fn4LlfzT+L7OyTtLrZ1l8fX0xfvx4mgLQUUZGRvD19UW/fv3AcRzy8/OxZcsW3L59W29P8St3SC5ZsgT3798HAHTs2BHLli2j+4lqmVQqxYkTJ7Bz504UFRVBIBDg/fffx5QpU+i66NdodLgwxvD48WPMmzcPaWlpfJ0kOsui2ziOg4mJCT7//HP07dsXHMfxZwf0MWCUp71XrVqFiIgIAGXz+itXroStrS39MaslynIu+/btQ1hYGMRiMb9WOH78eLo2/W80NlyUhyTnzp2LxMREAGU3SS5evBh16tShN1HHKQPmiy++QJ8+ffi7YDZt2oSoqCi9miITi8UIDAzEqVOnwBhD48aNsWbNGjRv3pz6QS1hjKGoqAi7d+/G8ePHIZPJYGxsDB8fH3zyyScwNjam9+JvNDJcXr9J8u7duwAADw8PrFq1CpaWlvQm6gllqf4vv/ySD5hXr15h8+bNiIyM1IuAkUgk2LlzJ/bv3w+FQgErKyt89913cHFxoX5QS5RTkps3b8alS5f4OmETJkzA4MGDYWhoSO/FG2hcuDDGkJ2djYULF+LGjRsAyg5Jfv/992jQoAG9iXrm9YB5//33y63BRERE6HTAyGQyhIeHY9euXZDJZDAzM8OSJUvQpUsX6ge1RLmJYv369fypeysrK0yfPl2n72JRBY1aEX+9rIuyAF+LFi2wfv16mgLQY8qA+eKLLyAUCnHu3DkUFBRg69atKC0tRe/evXWuk8vlchw+fBhBQUEQi8UwMTHBN998gwEDBtAOyVrCGOPvHHr+/DkAoHHjxpgyZYrenrqvDI0JF8YYcnJysHjxYly4cAGMMTRt2hTr16+nN5KUW+Q3MDDA6dOnUVRUhJ07d6KwsBCDBg2CoaGhupupEnK5HMePH8e6dev40/dTpkzB6NGjdS5ENZVcLkdUVBR27dqFnJwcAGWbKKZOnYomTZrQ36MK0IhweX0q7OLFi/yi5Q8//IB27drRG0kAlAWMsbEx/Pz8YGxsjOPHj0MsFiM0NBTZ2dkYO3as1u/YkclkOHr0KNauXYuioiIYGBhgwoQJfKiSmieVSnHu3DmEh4ejuLgYHMfBw8MDEydOhJWVlVZ/vmqT2j+tjDG8ePEC8+fP5wtRNmnSBOvXr0fHjh3pjSTlcBwHIyMjfPLJJzAzM8P+/fshFotx4sQJpKWlYfz48WjUqJFWfm4kEgnCwsIQHByMkpISGBgYYPz48Zg4cSJtva8Fyi3fysrcMpkMQqEQffr0gZ+fH+1SrSS1hotyTnPBggX8rjB7e3v88MMPcHd3pzeSvBHHcTA0NMSwYcNgZWWFn3/+Ga9evUJ0dDSeP3+OsWPHolu3bjAyMtKKzxBjjK9EEBYWBolEAkNDQ3z55ZeYPHkyjI2N1d1EnccYQ2ZmJnbv3o3o6GgwxmBsbIxRo0Zh6NChWvNZ0iRqCxeFQoHr169jyZIlfCVRJycn/Oc//6FtlqRChEIhevToATs7O+zcuRNPnjxBRkYGtmzZgitXruDDDz9E27ZtUadOHXU39a0UCgUSEhLw/fff48aNG1AoFBCJRJgyZQo+++wzGrHUAoVCgdjYWOzevZv/W1SvXj18/vnn6N69O61zVZFawkUikSAkJATBwcF49eoVAKBTp050MIxUmkAggLOzMxYvXoz9+/fjypUrEIvFuHv3Lu7fvw9bW1t07NgRTZo0UXdT/6GkpAS7d+/Gf//7X2RnZwMAbGxsMHfuXAwaNIjWWGqBVCrFqVOncODAARQWFgIoK5c/ceJEODk50c68alDLpzcnJwe//PILXy78ww8/xOLFi6ngG6kSjuNgaWmJCRMmoHPnzjhy5Aji4+Mhk8mQnp6Oly9fwt7eXt3N/IeSkhIcP34c2dnZ4DgO7dq1w4IFC9CuXTv6o1ZL8vPzcerUKRQWFkIgEMDT0xPjx4+nsjoqoJZwqV+/PubPn4+lS5fC29sbX375JUQiEb2ZpMo4joOBgQHc3d3h4uKCuLg4XL9+HbGxsRp7xaylpSWmTJmCNWvW4OOPP4avry8sLCyoH9QiKysr+Pn5YefOnRgwYACGDh0KExMTeg9UoNbDRVmMsnXr1pg7dy6aNWvGXwuqq5KTk1GvXj11N6OclJQUnTkX8iampqZ4//330a1bNwiFQty7d0+j1l4YY0hKSkKrVq0wb948NGrUCMnJyfycvy5KSUlB3bp11d2Mcp4/fw57e3v4+fmhfv36SElJUXeTalRGRkat9YNaD5d27dohOjqaLxseGRlZ202odYwxdOzYUd3N4LVq1QpxcXFISkpSd1NqDWMMzs7O6m4Gz83NDbdv38aDBw/U3ZRawxiDh4eHupvBa9WqFeLj4/l+EBsbq+YW1bza7Accq8X65fpWKv3vNGGoTe8BvQfqRu+B+tXGe1Cr4UIIIUQ/aN4qJyGEEK1H4UIIIUTl9CJcZDIZMjIyIJVK1d0UvZWXl4eHDx/q9P0rmq6oqAgvX76k90BNlHfDPHv2TN1NqRV6ES6JiYkYN24cPDw8cOTIEb1fzKtNUqkUP//8Mz7//HPs27cPxcXF6m6S3vrjjz/wwQcf4IMPPsC9e/fU3Ry9kpubi++++w7Tp0/HH3/8oRd/g/RmQb+oqAjLli3D4cOHMWzYMHz33XcwMzNTd7N0WmpqKtavX49nz57hs88+w9ChQzVip5A+S05Oxpw5c/DgwQP4+/vjyy+/1MgDprrk9u3b2LRpEwQCAWbMmIH27duru0m1Qm/CRem3337D4sWLYWlpiQ0bNsDd3V3dTdI5jDFcuHABO3fuhI2NDWbPno2WLVuqu1nkf2QyGTZu3IidO3eiS5cuWLduHezs7NTdLJ0jlUoREhKCkydPomPHjpg2bZrGHSKtSXoXLkDZSeFZs2bhr7/+gr+/PyZOnEiVT1WkqKgIP/74I65du4Z+/frhq6++gomJibqbRd4gIiICc+fOhVQqxZo1a9CnTx91N0lnPH/+HBs2bMDz58/x6aefYtCgQXo3atfLcAH+/+ntxx9/hJeXFzZs2ID69euru1laLS4uDoGBgSgqKsLUqVPRvXt3dTeJ/Ivc3FwsXLgQly5dgo+PD+bOnUv3x1QDYwznz5/HTz/9BDs7O8yaNQvNmzdXd7PUQm/DRenmzZuYNWsWxGIx1q1bh379+qm7SVpHoVDg4MGD2LdvH5ydnREQEABbW1t1N4tUEGMMYWFhWLduHZo3b47AwEC0atVK3c3SOoWFhdi6dStu3LiB/v37Y/z48Xod1HofLkDZNtn58+fj3Llz8PHxwcKFC2kqp4IyMzOxYcMGxMXFYezYsfj4449pilFLJSQkICAgACkpKViwYAHGjh2rd1M5VRUXF4egoCCUlpZi8uTJ6NKli7qbpHYULv/DGEN4eDhWrlyJZs2aITg4GK1bt1Z3szTa9evXsWXLFohEIgQEBKBt27bqbhKpptLSUqxbtw7h4eHo168fVq5cCQsLC3U3S2PJ5XIcPHgQv/76K9q0aYOZM2fCxsZG3c3SCBQuf/Po0SPMnDkTT548wcKFC+Hr60tPb39TWlqKn376CWfPnkW3bt0wZcoU2tatY86fP49FixbBxMQE//nPf+Dl5aXuJmmcly9fIigoCI8ePcKYMWMwevRo2tb9GgqXNxCLxVi7di327NmD999/H+vWrYOlpaW6m6URnjx5gvXr1yMzMxNfffUV+vXrR+Gro9LT0zF37lxERkZi4sSJmDZtGl29/D/Xrl3Dtm3bUKdOHcyaNYtmOd6AwuUdLly4gLlz58LIyAjr169Ht27d1N0ktWGM4cSJE/jll1/QtGlTzJ49WyPvpSeqJZfLsWvXLmzcuBGurq5Yv369Xr/vpaWl2L17Ny5evIju3btj4sSJGnUJnSahcPkXGRkZmDNnDiIiIjBx4kTMmjVL757eXr16hY0bNyI6OhpDhw7Fp59+qtO3WJJ/unv3LmbPno28vDwsX74cgwcPVneTal1SUhI2bNiA7OxsTJgwAb1796ZR+ztQuFSAQqHAzp07sX79eri6uiIoKAjNmjVTd7NqRUxMDIKCgqBQKODv769RNwmS2lVQUIDly5fjxIkTGDFiBJYsWaIXT+3KUXtISAjs7e0xa9YsNGrUSN3N0ngULpVw9+5d+Pv7Izs7GytWrMCIESPU3aQaI5PJEBoaiqNHj6JDhw6YOXMmrTsRMMZw7NgxrFixAjY2NggMDISrq6u6m1Vj8vLysGnTJsTExGDYsGHw8fHRu5mLqqJwqaSioiIsXboUR44cwfDhw7FixQqd2ymVlpaG9evX4+nTp/j0008xbNgwGv6TcpKTkxEQEIC4uDgEBARg/PjxOrdT6vbt29i8eTM4jsP06dPRoUMHdTdJq1C4VNGxY8ewZMkSWFlZISgoSCc+eIwxXLp0CTt27ICVlRVmz54NBwcHdTeLaCipVIrg4GDs2rUL3bp1w7p163SiMoNUKkVoaChOnDgBDw8PTJs2DfXq1VN3s7QOhUs1pKSkYObMmbh//z5mzZqFr7/+WmtPpxcVFWHbtm34888/8f7772PChAlUpYBUyPXr1zF37lzI5XKsWbMGvXv3VneTqiw1NRUbNmxASkqK3hacVBUKl2qSyWQICgrCtm3b0LlzZwQGBmpdAcz4+HisX78eRUVFmDJlCt577z11N4lomZycHCxcuBCXL1+Gr68vvvnmG62qq6W8JuKnn36CtbU1AgIC0KJFC3U3S6tRuKhIREQEZs+erVUFMBUKBX799Vfs27cPTk5OCAgIoHs9SJUxxrB37158//33aNGiBQIDA7ViWrWwsBDbtm1DREQEFZxUIQoXFcrLy8O8efNw/vx5+Pr6YsGCBRo7tZSVlYUNGzYgNjYWY8aMwZgxY7R2So9olocPHyIgIACpqalYuHAhPv74Y42dWlIWnCwpKcHkyZPRtWtXdTdJZ1C4qJjy6W316tWwt7dHcHAwnJyc1N2sciIiIrBlyxYYGxsjICAALi4u6m4S0TGlpaVYs2YN9u/fj/79+2PlypUatSgul8vx66+/4uDBg3B2doa/vz8VnFQxCpcakpCQgJkzZ+LZs2dYuHAhfHx81P70JhaLsXv3bpw9exZdu3bF1KlTdW4bNdEs586dw6JFi2Bqaor//Oc/6NSpk7qbhMzMTAQFBeHhw4cYO3YsRo0aRaP2GkDhUoNKS0uxdu1ahISEoF+/fli3bp3aypc/ffoUP/zwA16+fImvvvoK/fv3V3vYEf2Qnp6Ob775BtHR0Zg0aRKmTp2qtj/m169fx9atW2FqaopZs2bB2dlZLe3QBxQuteD8+fOYN28ejI2NERgYWKsXCTHGcPLkSfzyyy9o3Lgx5syZo9eFB4l6yOVy7NixA5s3b4abmxvWr1+Pxo0b19rrK0ftFy5cQLdu3TBp0iS9KF2jThQutSQjIwMBAQG4efMmJk2aBH9//xovI5Gfn4+NGzciKioKQ4YMwaeffgojI6MafU1C3iUmJgazZ8/Gq1evaq0A5pMnTxAYGIjs7Gx89dVX6NOnD43aawGFSy2Sy+XYuXMnAgMD4ebmhqCgIDRt2rRGXismJgbBwcGQy+WYMWMGPD09a+R1CKmsgoICLFu2DCdPnsTIkSOxZMkSmJqaqvx1lAUnQ0ND0bRpUwQEBFDByVpE4aIGd+/excyZM5GTk4PvvvsOw4cPV9nPlslk2Lt3L44cOUIFJ4nGUhbAXL58Oezs7BAYGKjSXYt5eXnYvHkz7ty5g6FDh8LHx4euiahlFC5qUlhYiKVLl+Lo0aMYOXIkli9f/tY5YMYYsrOzUVhYCDMzM1hbW79xWJ+WlobAwEA8efIEfn5+GD58OA3/iUZ7+vQpZs+ezZ+N+fzzz99aALOi/SAmJgYbN24EAEyfPh3u7u41+juQt2BErY4cOcJcXV1Z79692d27d8v9W25uLgsKCmIODg4MAP8fBwcHFhQUxHJzcxljjCkUCnbhwgU2duxYNnnyZPbo0SM1/CaEVI1EImHff/89a926Nfviiy9YZmZmuX+vaD+QSCTs559/ZqNGjWLfffcd/78T9aCRiwZITk6Gv78/7t+/j9mzZ2PChAk4d+4cRo8ejeLiYgBlT21Kyqc1U1NT7N27F4mJifjzzz/Rt29ffP311xpbFYCQd1EWwFQoFFizZg169eqFM2fOVKgf7NixA7dv30ZycjL8/PwwZMgQGrWrGYWLhpDJZNiwYQO2b98Oe3t7XL58GUBZ/a+3UXaeXr16YcWKFejRo0dtNJWQGpOdnY2FCxfiypUr6Nq1K/bs2QOgYv1g6NChCAoKooKTGoLCRcOcOXMGgwYNemdn+jtTU1Okpqaq7YAmIarEGMP27dsxZcoUVPTPE8dxEIlE1A80iG5dHacD4uPjKxUsAFBSUsI/4RGi7TiOg1gsrnCwAGWBRP1As9DIRYMwxuDo6IikpKRKdSyO49CyZUs8evSI5pmJ1qN+oBsoXDRIVlZWta6JzcrKgrW1tQpbREjto36gG2haTIMUFhZW6/sLCgpU1BJC1If6gW6gcNEg1S1/b25urqKWEKI+1A90A4WLBrG2toaDg0Ol54s5joODgwOsrKxqqGWE1B7qB7qBwkWDcByH6dOnV+l7Z8yYQYuYRCdQP9ANtKCvYfLy8tCkSROUlJRUaEuyQCCASCTC8+fPaX8/0RnUD7QfjVw0jIWFBQ4dOgSO495awE9JIBCA4zgcPnyYOhTRKdQPtB+FiwYaOHAgTp48CZFIBI7j/jHMV/5vIpEIv//+OwYMGKCmlhJSc6gfaDcKFw01cOBAPH/+HEFBQWjZsmW5f2vZsiWCgoKQmppKHYroNOoH2ovWXLQAYww5OTkoKCiAubk5rKysaNGS6B3qB9qFwoUQQojK0bQYIYQQlaNwIYQQonIULoQQQlSOwoUQQojKUbgQQghROQoXQgghKkfhQgghROUoXAghhKgchQshhBCVo3AhhBCichQuhBBCVI7ChRBCiMpRuBBCCFE5ChdCCCEq93/a7CHBsUQMWQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 500x400 with 7 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot(beta=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64d2573b",
   "metadata": {},
   "source": [
    "Fix the first layer activation to be linear function, and the second layer to be sine functions (caveat: this is quite sensitive to hypreparams)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e2e78752",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best value at boundary.\n",
      "r2 is 0.9969676978399866\n",
      "Best value at boundary.\n",
      "r2 is 0.9983639008937205\n",
      "Best value at boundary.\n",
      "r2 is 0.9974491732032462\n",
      "Best value at boundary.\n",
      "r2 is 0.9978791881996706\n",
      "r2 is 0.9723468700787765\n",
      "r2 is 0.9844055428126749\n"
     ]
    }
   ],
   "source": [
    "for i in range(2):\n",
    "    for j in range(2):\n",
    "        model.fix_symbolic(0,i,j,'x')\n",
    "        \n",
    "for i in range(2):\n",
    "    model.fix_symbolic(1,i,0,'sin')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fae3f32",
   "metadata": {},
   "source": [
    "After setting all to be symbolic, we further train the model (affine parameters are still trainable). The model can now reach machine precision!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "308b72af",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "pde loss: 1.37e-16 | bc loss: 3.89e-18 | l2: 7.38e-18 : 100%|█| 20/20 [00:07<00:\n"
     ]
    }
   ],
   "source": [
    "train()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35985ae9",
   "metadata": {},
   "source": [
    "Print out the symbolic formula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f0ec310e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0.5 \\sin{\\left(3.14159 x_{1} - 3.14159 x_{2} + 7.85398 \\right)} - 0.5 \\sin{\\left(3.14159 x_{1} + 3.14159 x_{2} + 1.5708 \\right)}$"
      ],
      "text/plain": [
       "0.5*sin(3.14159*x_1 - 3.14159*x_2 + 7.85398) - 0.5*sin(3.14159*x_1 + 3.14159*x_2 + 1.5708)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula, var = model.symbolic_formula(floating_digit=5)\n",
    "formula[0]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}