{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Solving Partial Differential Equation with Q-Exponential Process\n", "\n", "## Introduction\n", "In this example, we introduce a Bayesian solver to PDE by modeling its solution as a QEP $u\\sim \\textrm{q-EP}(0, \\mathcal C)$. When choosing a differentiable kernel $\\mathcal C$, the derivatives of the PDE solution, $\\frac{\\partial^{\\boldsymbol \\alpha}}{\\partial {\\bf x}^{\\boldsymbol \\alpha}}$, are also QEPs based on the preservation property of linear combination. Assuming a QEP prior, $\\textrm{q-EP}(0, \\tilde{\\mathcal C})$, for the extended function $\\tilde u = (u, \\frac{\\partial}{\\partial {\\bf x}}, \\cdots, \\frac{\\partial^k}{\\partial {\\bf x}^k})$, we define the likelihood by penalizing the discrepancy between the left-hand side $P(\\tilde u)$ and the right-hand size $h({\\bf x})$ of PDE and approximate the posterior using sparse variational Bayes. See more details in our [NIPS2025 paper](https://nips.cc/virtual/2025/loc/san-diego/poster/117572)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using the cpu device...\n" ] } ], "source": [ "import torch\n", "import qpytorch\n", "import random\n", "import numpy as np\n", "import tqdm\n", "from matplotlib import pyplot as plt\n", "\n", "%matplotlib inline\n", "\n", "# Setting manual seed for reproducibility\n", "seed=2025\n", "random.seed(seed)\n", "np.random.seed(seed)\n", "torch.manual_seed(seed)\n", "torch.cuda.manual_seed(seed)\n", "torch.cuda.manual_seed_all(seed)\n", "# set device\n", "device = 'cuda' if torch.cuda.is_available() else 'cpu'\n", "print('Using the '+device+' device...')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Burgers' Equation\n", "We consider the following Burgers' equation with $\\nu=0.1$:\n", "\n", "\\begin{align}\n", "\\frac{\\partial}{\\partial t} u + u \\frac{\\partial}{\\partial x} u - \\nu \\frac{\\partial^2}{\\partial x^2}u &=0, \\quad (x, t)\\in (-1, 1)\\times (0,1], \\\\\n", "u(x, 0) &= -\\sin(\\pi x), \\quad x \\in (-1, 1), \\\\\n", "u(-1, t) &= u(1,t) = 0, \\quad t \\in (0,1] .\n", "\\end{align}\n", "\n", "We first create $25\\times 25$ mesh using the following code." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# generate mesh points\n", "def sampled_pts_grid(N_domain, N_boundary, domain, time_dependent = False):\n", " x1l, x1r = domain[0]\n", " x2l, x2r = domain[1]\n", " \n", " N_pts = int(torch.sqrt(torch.tensor(N_domain + N_boundary)).item()) - 2\n", " xx = torch.linspace(x1l, x1r, N_pts + 2)\n", " yy = torch.linspace(x2l, x2r, N_pts + 2)\n", " XX, YY = torch.meshgrid(xx, yy, indexing='ij')\n", " \n", " if not time_dependent:\n", "\n", " XX_int = XX[1:-1, 1:-1]\n", " YY_int = YY[1:-1, 1:-1]\n", " \n", " XXv_bd = torch.cat((XX[:-1, 0], XX[0, :-1], XX[:-1, -1], XX[-1, :-1]))\n", " YYv_bd = torch.cat((YY[:-1, 0], YY[0, :-1], YY[:-1, -1], YY[-1, :-1]))\n", " \n", " else:\n", "\n", " XX_int = XX[1:-1, 1:]\n", " YY_int = YY[1:-1, 1:]\n", " \n", " XXv_bd = torch.cat((XX[-1,1:], XX[:, 0], XX[0, 1:]))\n", " YYv_bd = torch.cat((YY[-1,1:], YY[:, 0], YY[0, 1:]))\n", " \n", " # vectorized (x,y) coordinates\n", " XXv_int = XX_int.flatten().unsqueeze(1)\n", " YYv_int = YY_int.flatten().unsqueeze(1)\n", " \n", " XXv_bd = XXv_bd.unsqueeze(1)\n", " YYv_bd = YYv_bd.unsqueeze(1)\n", " \n", " X_domain = torch.cat((XXv_int, YYv_int), dim=1)\n", " X_boundary = torch.cat((XXv_bd, YYv_bd), dim=1)\n", " return X_domain, X_boundary" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now on the mesh with $N=N_d+N_b$ points of which $N_d$ are in the domain and $N_b$ are on the boundary, we have the extended function `mu` as an $N\\times D$ array where $D=1+kd$ with max $k=2$ order derivatives and $d=2$ physical dimensions. Define the PDE class to specify the left-hand function `lhs_f`, the right-hand function `rhs_f`, the boundary condition `bdy_g`, and extract the solution/derivatives in the domain and on the boundary respectively." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# PDE class\n", "class Burgers(object):\n", " def __init__(self, alpha=1.0, nu=0.02, domain=np.array([[-1, 1], [0, 1]])):\n", " self.alpha = alpha\n", " self.nu = nu\n", " self.rhs_f = lambda x: torch.zeros(x.shape[0])\n", " self.bdy_g = lambda x: (-torch.sin(torch.pi*x[:,:-1].squeeze()) * (x[:,-1]==0) + 0*((x[:,0]==1)+(x[:,0]==-1)))\n", " self.domain = domain\n", " self.dim = self.domain.shape[0]\n", "\n", " def sampled_pts(self, N_domain, N_boundary, sampled_type='random'):\n", " Xd, Xb = sampled_pts_grid(N_domain, N_boundary, self.domain, time_dependent=True)\n", " if sampled_type == 'random':\n", " Xd += torch.randn(*Xd.shape) * 1e-2\n", " elif sampled_type == 'grid':\n", " pass\n", " self.X_domain = Xd\n", " self.X_boundary = Xb\n", " self.Nd, self.Nb = self.X_domain.shape[0], self.X_boundary.shape[0]\n", " self.rhs = self.rhs_f(self.X_domain)\n", " self.bdy = self.bdy_g(self.X_boundary)\n", " \n", " def lhs_f(self, u_):\n", " u_d, u_b, du, d2u_x = self.extract_solution(u_)[0]\n", " lhs = torch.cat([du[:,-1] + self.alpha*u_d*du[:,:-1].squeeze() - self.nu*d2u_x[:,:self.dim-1].squeeze(),u_b],-1)\n", " return lhs\n", " \n", " def extract_solution(self, mu, var=None):\n", " u_d, u_b = mu[...,:self.Nd,0], mu[...,self.Nd:,0]\n", " du = mu[...,:self.Nd,1:1+self.dim]\n", " d2u_x = mu[...,:self.Nd,1+self.dim:]\n", " if var is None: var = torch.ones_like(mu)\n", " v_d, v_b = var[...,:self.Nd,0], var[...,self.Nd:,0]\n", " dv = var[...,:self.Nd,1:1+self.dim]\n", " d2v_x = var[...,:self.Nd,1+self.dim:]\n", " return [u_d, u_b, du, d2u_x], [v_d, v_b, dv, d2v_x]\n", "\n", " def plot_solution(self, u, ax=None, **kwargs):\n", " import matplotlib.pyplot as plt\n", " x, t = self.X_domain[:,0], self.X_domain[:,1]\n", " N_pts = int(np.sqrt(self.Nd+self.Nb))-2\n", " x = x.reshape((N_pts,N_pts+1))\n", " t = t.reshape((N_pts,N_pts+1))\n", " if ax is None:\n", " ax = plt.gca()\n", " ctf = ax.contourf(x, t , u, **kwargs)\n", " clb = plt.colorbar(ctf, ax=ax)\n", " clb.ax.tick_params(labelsize = ax.yaxis.get_tick_params()['labelsize'])\n", "\n", "# generate PDE\n", "alpha = 1.0\n", "nu = 0.1\n", "domain=np.array([[-1, 1], [0, 1]])\n", "N_dom, N_bdy = 625, 100\n", "burgers = Burgers(alpha, nu, domain)\n", "burgers.sampled_pts(N_dom, N_bdy, sampled_type='grid')\n", "burgers_X = torch.cat([burgers.X_domain,burgers.X_boundary]).type(torch.float).to(device)\n", "lims = torch.from_numpy(domain.T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reference solution\n", "As a reference, we solve the Burgers' equation using highly-resolved finite difference method with Cole-Hopf transformation." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAHRCAYAAACxcxlEAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVIRJREFUeJzt3QmczfX+x/GPbYxdMrcIY08lQpckkRR/kiyRqxupW7fQwk2RrdIebaiuLq66lYwt1CWyRRSTrSxlG2MJN8sMM5Yx/8fne/3mnpk5M3P28zvn93o+HufxO37nnN/5zTn4vee7fL6FMjMzMwUAAMAmCof7BAAAAFwRTgAAgK0QTgAAgK0QTgAAgK0QTgAAgK0QTgAAgK0QTgAAgK0QTgAAgK0QTgAAgK0QToAo0LdvXylUqJBUr15d7GDq1KnmfPS2Z8+ecJ8OgAhDOEHYLVu2LOtC5slNL3yevL5o0aJSoUIFqVGjhtx8883y5JNPysyZM+Xs2bMBPX9dAeKLL76QXr16SZ06daR06dLmvcuXLy/169eXu+++W15//XXZuHFjQN8XAKIV4QRRKyMjQ44dO2Z+c1+5cqW89dZb0r17d6lSpYqMGTNGzp8/7/d7/Pbbbyb4dO7cWT777DP59ddf5dSpU+a9T5w4IT/99JMkJCTIkCFD5LrrrpNt27ZJpHINgXofAIKlaNCODPjgkUcekUcffTTf52i48PT1qampJqBs2rRJlixZIosXL5YjR47IiBEjZN68eTJ//nyJi4vz6bvSFpjbbrtNNm/ebP7cqFEjuf/++00IKVOmjJw8eVK2bt0qK1askAULFpiw4qRuJr0BgC8IJ7CVP/zhD6YrJNCv/7//+z95+umn5eeff5Z7771XfvzxR/n++++lS5cu8s0330hMTIzX7zVp0qSsYKKh5MMPP5TChbM3RmqrysMPPyxnzpyRTz/91HT1AADyR7cOHOXqq6+WVatWmVYOpfcnTJjg07Hmzp1rtjq+ZNy4cbmCiavixYubloTLL7/cxzMHAOcgnMBxSpQoIR999JEZO6HeeOMNOXfunNfHSUpKMtuKFSsGpEVEu6BeeeUVad68uRnIq4FGu7B0nIx2P/lKx9zkNZg4J53to89z7ZKxXn/LLbdk7dP7+Q1U9nS2jnaxDR8+3IRF/QxjY2PNOfz5z3+Wb7/91qtz3b59u/zlL38x+/Wzu+yyy0zL2Jo1azz6nADYB+EEjnTNNdeY8SLqwIED8sMPP3h9DKsrSAfF/v77736dj3YzXXnllTJ06FBzMdVxMjqmZf/+/WaGUadOnaRbt26Snp4u0WLRokVSu3ZtefHFF2XDhg1mTI52f+3du1c+/vhjadmypQwYMEAuXLhQ4LFmz54tjRs3Nl1r+nr97A4fPixz5syRm266SaZPnx6SnwlAYBBO4Fht27bNuq+zebylF0NrKrH+xq4tH77QAHLrrbeakKQtATp+ZeHChbJu3TqZNm2aNGzY0Dxv1qxZYRlkesUVV5ixNZMnT87ap/d1n+vtrrvu8viYGkY0cOmg4WLFiplp3kuXLjXjgD744AMz/Vtpl5sGtvzoe//pT38yLSXjx4834e67776T0aNHm5YYnTn10EMPmVYaAJGBAbGwFf1td8uWLfkOeNVbIFjhQu3YscPr1+usIO0e0t/sNTjowFq94Opv/M2aNTOtM0WKFCnwOE888YRpKbEG2T7wwANZjzVp0kR69OhhBvTqxVtbAPr06WP+HCoaHnSQ8dGjR7P2aXjwZ+CyhgVt3dDPR7usbr/99qzH/vjHP5raMNrioQOYtdvtvvvuM5+nO4mJieZz0s+/bNmyWftvuOEG0zKjA6A1BGlrjIYgAPZHywls5b333pNrr702z9vEiRMD9l6XXnpp1n0rHHijadOm5rd8vXir48ePm7CiF15t7ShXrpy56Grg0Non7mhriXZJqPbt22cLJhYdP6EtFTrwVmnrQCTT1hGrG01bnFyDieWSSy6Rv//97+a+hr+Cvnf9fFyDiUVbVCpXruxz6xiA8CCcwLG0kqslJSXFp2M8+OCDpltBu2K0tokrDSRff/21CStaOfbf//53rtdrMTPtdlDugolFB3laY2RcXxOJtNaMJb+fuUWLFnLVVVflek1OGlobNGjg9jHtJrNmZu3atcuPswYQSoQT2MqoUaPMGI68bjqOIFBcA4m737o9pQNZ9Tf3//znP7J69Wozrbh3797ZisUdPHhQ7rjjjlwXWdcuLO0Kyo/1+OnTpyP6Qmv9zDqgWAvWefIz//LLL3kuO1CvXr18j6Ezn/wJoABCjzEncCzXMRTWBcwf2r2j04D1ZtFxEI899pgpY6+tHTpORae8WtOYXWf5FDSWxrVGir+zg8LJOnf9zK2uqoJ+Zg2m2vWmg15zKlmyZL7HsOrPRHJrE+A0tJzAsXT6rmvrRzC0adPGdO1Y4UdbAHSmijtWYHEKp/28ADxHOIFjaWiw6MyQYKlUqZJ07Ngx68+6OKC7Fhutl5KfQ4cOuX1dQVwr1xZUMySvgbuBZJ27doMVtPii9TNrkNFBsgCcgXACR9JxD7oQoKpatapcf/31QX0/a8ZIzhYD1+m4a9euLXCWi9WNUbNmTY/f23Wgbn6zkrS7RQNDsFs6rJ9Zx5Dk1YqU82fWAcW+rH8EIDIRTuA4aWlppm6GjmNQf/vb3woc++CO9XpPaEE1i2uwaN26dVYtFNciZ+5K5VstPa6v8YS2OFjl9V3PI6fPPvss359JC5pZtJJrIIrf5fczayE1rXOS8zUAoh/hBI6iFzvtwrHGm7Rq1UoeeeQRn47VtWtXU3+joK4QXWfGaqWpVq1a1tRWq0VF139RX331lfzzn//M9XptYejXr1/W+j9a0t1bujqytVjhzp07cz2ug3RHjBhRYPeUxd0xvKkPY7VUaQ0Y67NxpaXsdTVnq1vK1+8IQGRitg6iusKsBgftyti0aZO5CGrrg9U6oBVEExISsoqoeWvfvn3Sv39/efrpp01lWA0AOrBWWyp0DZxt27bJjBkz5Msvv8zqFnnzzTdzdY/oPj03PU8NIbrgXc+ePc1x9BhaIdXq/rCqxXpLZwl98cUXptVIW150SraGJC25r+/99ttvS1xcnGmRyavMuwYrnR6dnJxszknv689rteLoTJqctV7yoqFEpwlr8OrQoYMMHDjQfIalSpUywVEXQLSmS2vLlj/VaAFEoEwgzJYuXappwdxGjRrl1+s9ucXFxWW++OKLmefOnfPrvDt37uzxe5YrVy5z2rRpeR4rMTExs3Llyvkeo2vXrplpaWluX9+nTx/znPj4+Dzf47HHHsvz2NWqVcv8+eefzev1z3o8dyZOnJjnMaZMmZL1PL1v7d+9e7fbYy1cuDCzbNmy+f7M/fv3z8zIyHD7+oLO1ZvPBoC90HKCqKXdAfqbvJaRj4+PN+uv6Lo3WgwtEIMrdcVb7Q7RRfpWrVplaploq4K2Ruj4DC2Pr7/xa3l2LcpWsWLFPI+lrRh6LC1Nbx1Xi63pa7SFRxf805YFf2jriB7r/fffNy0x2k2krSHaraStE67l/POi3SvaQqJl+/UYOoi2oBk3edHPRWcuvfXWW6Z1SVtKdCyLHl+/p7/+9a9BnUUFwL4KaULx5gXr1683TeM6il5vuqKq8vIwWbQpW5uY9T9knTaoRZf0P0vdZw3iAwAg2tjhepqRkSHvvPOOGZyuvyzosh633HKLPPfcc1nLR0REONFl0XVQXU6+fJhaoVOraeoHojMYdJCc/vapt7p165rR+oGo3AkAgN2E+3p64cIF6d69u1l8VMPLrbfeao6zYsUKKVGihFkJXQewh4W3/UCvvPJK5ogRIzK/+OKLzIMHD2YWL17c9Of6onfv3ll96a79/wMHDvSoLxkAgEgV7uvppEmTzGN16tTJPHToUNb+hIQEs7927dp+j83zldctJzlp37r2E3t7GF0ITUf7a30JreHgumaGHk8LY2l/ti4pX9CaIwAARLpQX0+vvvpq2bp1q2k50VYcV507dzYz/HRGY7du3cQxdU50+XhtUtKBbzkX8ypevLgZ/Kd9YdY0TAAAEJjr6e7du00w0e4b1+U1LNrdo+bNmyfhELZwsnHjRrNt3Lix28et/VqfAgAABO56ar1GZxS6q/UU7mtw2MKJNj0pbYpyx9q/d+/ekJ4XAACRJMmH66ndr8Fhq3OitSCsRczc0UqRKiUlJd/jaH+a6zof2rSlfWtas4El2QEgcunYC70G6DIPrqtrB4pWctYqxYE615zXHO1S0Zsdr6epAboGB0vEF2F7+eWXzXxsAEB00qUi8voN359gUq1aKTly5EJAjqf1QawLvmXUqFGmxggiKJzoF6m0CqY71mJqBa3VMXToUBk0aFC2BcO06uUV456RwiX+t4qq01SrfDTcpxAVkg7kXdU13GL2+V/lNpjK7PNrIqDjpVTN/lu4E104ky67xz7v8ZpN3tAWEw0my9b+QUqX9u+zTk3NlNbNDpsQVbZs2az9oWg18fV6GqhrcNSFEw0QSst9u2Pt17Lj+cmr2UyDiVPDSfUqunBbaP5RRDs7/x0qEmvvcFIkhnDi1+cXSzixBLOLXoNJ6TL+dhn9t/VFg4lrOLHz9bRagK7BUTcgtmHDhmabmJjo9nFrf4MGDUJ6XgBgB2X2Eu4QvOup9RpdxV3X2fLkNY4IJ+3btzcDnFauXGmWuXelA1x1brUuxa7LqQPhbYUCAPtq78P1tEaNGmbtnLS0NFmwYEGuY2rxNeXvgqO2DSe6ymq9evXM2BBXlSpVkl69epl+v0cffTTbyqZDhgyRI0eOyL333kt1WCAC8Vs/YP/r6aCL4zX1Oa6hZtasWaY6bO3atU2l2IgYc6IJ64UXXsj6szUNS5dit4wYMSKr4pwuIqTLv2t53Zx0qfQ1a9bIzJkzzQduLVSkzUx16tSRcePG+fpzORa/6QNAZAj39bRfv36maqyWr9fXWAv/LV++3FSO/fjjj01J/IhoOdEEtnbt2qybtQaA6z59jicqVqxolokeOHCg+VL0A9LZNo899pjZz4rEAJyMFqjoFu7raeHChWXGjBkyduxYU0tm/vz5snnzZrOWzrp166RZs2YSLn4v/Gc3J0+elHLlyknV90bbeqZFsNByEnh7kuPEjoon2Xe2DhfVwEmJd+6snYz0dNn50jBzkQ30LBjrWrHup8v8nq2TmnJBrr/mt6Ccp1OFbUAsAACAO4QToAC0RgFAaBFOoggXUdgBXToA/EU4AQAbI+zBiQgnAADAVggnAADAVggnUYLxJgCAaEE4ATxA+AOA0CGcAAgYBm8GB58rnIZwAgAAbIVwAgAAbIVwEgUYDwEAiCaEEwAAYCuEEwCIAAyKhZMQTgAP0X2WPy6eAAKFcAIAAGyFcBLh+G0eABBtCCcAECHoOoNTEE4AAICtEE4AAICtEE4ALzDGxz26GwAEEuEkgnGhBABEI8IJAEQQWqngBIQTAABgK4QTAABgK4QTAABgK4STCMVgWD57AIhWhBMAfmGAZujxmSPaEU4AAICtEE4AAICtEE4AAICtEE4iEINhASA6pKWlyciRI6Vu3boSGxsrlStXln79+sn+/fs9PsbUqVOlUKFCBd6mTZuW7XV9+/bN9/nvv/++hEvRsL0zAMCvQbEp8YX4BCNYenq6tGnTRtasWSOVKlWSzp07y549e2TKlCkyf/58s79mzZoFHqd27drSp08ft4+dOHFC5syZY+7fdNNNbp/Trl07ufzyy3Ptv/LKKyVcCCeAj61Xe5LjHP/ZMWsE8N2YMWNMAGnevLksWrRISpcubfaPGzdOBg8ebFpQli1bVuBxNHTkFTzee+89E05atGiRZ9B55plnpHXr1rb6KunWAQAgxM6ePSvjx4839ydMmJAVTNSgQYOkQYMGsnz5clm/fr1f7/Pxxx+b7Z///GeJJIQTAABCbNWqVabLpVatWtKoUaNcj3fv3t1s582b5/N77N69W1avXi0xMTHSo0cPiSR060QYBsMCQOTbuHGj2TZu3Njt440v7t+0aZPfrSYdO3aUSy65JM/nzZo1S2bOnCkZGRlSo0YN6dSpk9SrV0/CiXACABGKQbGRKykpyWyrVKni9vEqF/fv3bs36F067777brY/P/300/LII4/I22+/LUWLhicm0K0DAEAAnDx5MtvtzJkzeT43NTXVbEuWLOn28VKlSpltSkqKT+fy/fffy44dO6RChQqm5cQd7U7S6cL6vNOnT8uuXbvM+Jfy5cvLxIkT5amnnpJwoeUE8BEzdoDINyelocRmFvPrGOmp50RkkVStWjXb/lGjRsno0aMlHD6+2GqiY010zIk7jz/+eLY/a5fOo48+Kq1atTLdSjpgVwfn5vy5QoGWEwA+YRoxkN2+ffvMIFfrNnTo0Dw/Imt2jrZYuHPq1Kn//jsrU8brj/n8+fMyffp0n2fpXHPNNXLnnXea4yxZskTCgZaTCMJgWACwr7Jly5qbJ6pVq2a2ycnJbh9Pvrg/Pj7e6/PQmimHDx82dU1uvPFG8UWdOnXM9uDBgxIOtJwAQASjBSsyNWzY0GwTExPdPp54cb/WO/G1S+fee+/1+fyOHTuWbexLqBFOAAAIMa3YWq5cOdm5c6ds2LAh1+MJCQlmq9N6vaEDbefOnetXONGBvAsWLMh3qnOwEU4AAAgxHaQ6YMAAc79///5ZY0ys8vWbNm0yA1ObNGmStV8HqGr9kfzGsmjNEh3HcsMNN2R1zbizbds2+eijj3LNKDpy5Ijcc889ZvyMtu5oiAoHxpwAABAGw4cPl8WLF5sqrhokWrZsaeqarF27VuLi4mTy5MnZnn/06FHZvn17vuNAPK1tcujQIbnvvvvMjJ3rr7/evN+BAwdMuXydvqx1Vj7//HOzOnE40HISIRgMa09O/V4Y5wD4LzY2VpYuXSojRoww9U50gT4NJ3379jVjTjxZkdiVhpZvvvlGihUrJj179sz3uXXr1pUnnnjCrDy8efNmmTFjhqxbt86EJJ0CrS03+pxwKZSZmZkpUUQL32g/XtX3RkvhErESLZx6EYwE4VqduHiS+9oFoUA4sZeU+PD8dhtsGenpsvOlYWZarqezYLy9Vgxfc7vElva/zsmYGxYF5TydipYTAIhwhEVEG8IJAACwFcIJAACwFcJJBGC8CQDASQgngJ8IjwAQWIQTAF5h8KU98b0gmhBOAACArRBOAACArRBOAACArRBObI7BlgAApyGcAECUYFAsogXhBAgAWrgAIHAIJwA8xm/mAEKBcAIAAGyFcGJjdBUAAJyIcAIAUYSuNzg6nKSlpcnIkSOlbt26EhsbK5UrV5Z+/frJ/v37vT7W119/LR07dpS4uDgpVqyYXHrppXL77bfL7NmzfT09AADgpHCSnp4ubdq0kRdeeEFSU1Olc+fOUrVqVZkyZYo0atRIdu3a5fGx3nrrLRNEvvrqKxN0unXrJvXq1ZPFixdL165d5dlnn/XlFAEAgJPCyZgxY2TNmjXSvHlz2bFjh0yfPl3Wrl0rY8eOlSNHjpgWFE/oc5955hnTWrJ06VJZtWqVfPbZZ2a7bNkyKV68uLz88stehR0gXBgjBABhCidnz56V8ePHm/sTJkyQ0qVLZz02aNAgadCggSxfvlzWr19f4LE00Jw5c8a0wrRq1SrbYzfffLO0a9dOMjMzZd26dd6eJoAAYywDANuGE23VOHHihNSqVct04eTUvXt3s503b16Bx9KWEU/oGBSn4bdwAIBTeR1ONm7caLaNGzd2+7i1f9OmTQUeq2nTplK+fHn55ptvTGuLqxUrVsjChQulTp060rJlS29PEwAci1YuRLqi3r4gKSnJbKtUqeL2cWv/3r17CzxWuXLl5B//+If86U9/kltuuUVuvPFG8/rk5GRZvXq1tGjRQqZNmyYxMTF5HkO7hfRmOXnypLc/EgAAiORworNzVMmSJd0+XqpUKbNNSUnx6Hg6I0dn6vTo0cN0GVnKli1rZvFcccUV+b5eB8w+99xzXvwEAADAzsJehE1n+LRt29YMgNWuIA0/utVBslpHRcNLfoYOHWrGwFi3ffv2hezcAQCADcKJNTvn9OnTbh8/deqU2ZYpU6bAY+l04b/97W9y3XXXyYwZM+Taa681LS+6TUhIMPsXLFhgWlbyG1SrrSyut0jHYNjIFa3fHWMYANg6nFSrVs1sdVyIO9b++Pj4Ao/10UcfmW2XLl2kcOHsp1KkSJGsVhMdHAsA8ByBEo4KJw0bNjTbxMREt49b+7XeSUGsIKMDY92x9h87dszb0wQAAE4JJzqDRkPDzp07ZcOGDbke1+4Y1alTpwKPdfnll5ttXkXWfvjhB7OtXr26t6cJAACcEk50Wu+AAQPM/f79+2eNMVHjxo0zg1m12muTJk2y9mtFWV0vRwevurrrrrvM9l//+pfMnz8/22Nz586VTz75xHT3aLcPAABwBq+nEqvhw4ebhfm0FolVJE3rmmg5el1ZePLkydmef/ToUdm+fbscPHgwVzi5++67zWBYbWm5/vrrpUaNGrJ79+6s1pQXX3xRrrzySn9+RgAAEO1TiWNjY81CfSNGjDD1TubMmWPCSd++fc2Yk5o1a3p0nEKFCplFA7UQm04l/vXXX2X27NmyZ88e6dChg5mlM2zYMHGSaJ3tAQBA0OuclChRQp5//nkTKLRCq7aKTJkyxW3l2NGjR5sF/KZOneo2oOgqxlq+Xge+njt3zqxWrFOI27dv7+vpAWETbQGTWR+Ri+/O/tLS0kxNr7p165pf/CtXrmyuifv37/fqODo2U6+ned22bdvm9nUZGRny5ptvmhIeel3X3g8tirp161aJuG4dAADgn/T0dFNwdM2aNVKpUiXp3Lmz6TmYMmWKGYep+z3tibD06dPH7X53s2IvXLhghlZoj4Wuc9exY0czDEMntmgDgfaQ6Bp44UA4AQAgDMaMGWMCSPPmzWXRokVZRU7HjRsngwcPNi0oWqzUG+56KPKi40M1mOjY0ZUrV8pll11m9s+cOVO6d+8uvXv3Ni0oRYsWdV75egAAnObs2bNmJquaMGFCVjBRgwYNMrXCdLjD+vXrg3YOGoLUa6+9lhVMVLdu3eTOO+80wzZ05mw4EE4AAAgxXehW14OrVauWNGrUKNfj3bt3N9t58+YF5f11Vqy2iug4E+3OCfX7F4RuHRuJtoGUAAD3Nm7caLaNGzd2+3jji/u1dpg3Xn/9dVMkVdedu+aaa0ydMB3kmtf7169fX4oVKxaw9w8UwgkARPmMnZT4QuE+DeSQlJRktu5muLru1zId3hgyZEi2Pz/55JPy7rvvmvEroXj/QKFbBwgCWsEA5zl58mS2m5bZyEtqaqrZaq0wd0qVKmW2KSkpHr23jhGZNWuWCROnT5+WLVu2mLEreg4PPvhgrrEjgX7/QKPlBECeqJOBaLf0t7pSNLW4X8c4f0pDyCKpWrVqtv2jRo0ydb5C4Z133sn2Z+3SGTt2rFk65qGHHpKnn37aTFWOFIQTAAACYN++fVK2bNmsP+u4j7xYs3O0lcOdUxfXrStTpoxf5/TAAw+YJWd0CRmtoWItpBuq9/cV3ToAAASABhPXW37hpFq1amabnJzs9vHki/vj4+P9OiddPFdnBCnX9e1C9f6+IpzYBGMUAMA5GjZsaLa6Hp07iRf3a70Tf+nSMK7jSFzfX8em6LIxwXx/XxBOACDKMXbIflq0aGFKyuu03w0bNuR6PCEhwWw7derk1/v89NNPpktHB77q+BNLjRo15KqrrjJr+2ip+mC9v68IJwAAhFhMTIwMGDDA3O/fv3/WGA+rcuumTZukVatW0qRJk6z9WlFWA8bQoUOzHevLL7+Ub775Jtd76DF07RxdeFdn7Oh7utLZPNb048OHD2ft11k/X3zxhdSuXTtsg2gZEAsAQBjoQNXFixfL6tWrzfo2LVu2NFOB165dawqn6do3rnRRPm0FcR07or7//nt57rnnzPgQ7a7RVpJdu3aZrpnz589L69at5ZVXXsn1/lr7RIONrq+joefWW28176Fl87Vy7McffxyWdXUULSdAkET6OCK6AoDgio2NNSv/jhgxwgSKOXPmmHDSt29fEyw8XZG4Xbt2JmjoIFwti69dMrouzk033SSTJk0yAUjDhrvBsjNmzDBTjitXrmxWQt68ebNZW2fdunXSrFkzCZdCmdreE0W08I3241V9b7QULhErkSLSL2Rwb09y7rLRgVI8KXsTbaARTqJLJFaJzUhPl50vDTNr0LhO0Q3ktaLF3AFStJT/dU5WdR4flPN0KlpObIBgAgDA/xBOAMABaAlDJCGcAAAAWyGcAAAAWyGcAAAAWyGcAAAAWyGchBkzdQAAyI5wAgAAbIVwAgAOwXRiRArCCRBEkdptx0UMQDgRTgAAgK0QTgAAgK0QTsIoUpv8AQAIJsIJAACwFcIJADgIg50RCQgnAADAVggnAADAVggnAADAVggnAADAVggnYcI0YgAA3COcAEFGEIXdMGMHdkc4AZANFy4A4UY4AQAAtkI4AQAAtkI4AQAAtkI4CQMGSAIAkDfCCQA4EAOfYWeEEwAAYCuEEwAAYCuEEwAAYCuEEwAAYCuEkxBjpg4AwJKWliYjR46UunXrSmxsrFSuXFn69esn+/fv9/hDOn78uHzyySfSq1cvqVGjhsTExEiZMmWkWbNm8vbbb8u5c+fcvq5v375SqFChPG/vv/9+2L6oomF7Z8BhoXRPcpzYHTM4gNBJT0+XNm3ayJo1a6RSpUrSuXNn2bNnj0yZMkXmz59v9tesWbPA47zxxhvy4osvmkBx3XXXmVBy5MgRWbVqlXz//feSkJAgCxculJIlS7p9fbt27eTyyy/Ptf/KK6+UcCGcAIBDaRhNiS8U7tNwrDFjxpgA0rx5c1m0aJGULl3a7B83bpwMHjzYtKAsW7aswOOUKlVKhgwZIv3795dq1apl7f/ll1+kbdu28u2335r3eumll9y+/plnnpHWrVuLndCtAwBAiJ09e1bGjx9v7k+YMCErmKhBgwZJgwYNZPny5bJ+/foCjzV06FB59dVXswUTVadOHXnllVfM/U8//VQiCeEEAIAQ0y6XEydOSK1ataRRo0a5Hu/evbvZzps3z6/3adiwodkeOHBAIgndOgAAhNjGjRvNtnHjxm4fb3xx/6ZNm/x6n127dpmtuzElllmzZsnMmTMlIyPDDKjt1KmT1KtXT8KJcAIAQIglJSWZbZUqVdw+XuXi/r179/r1PjpbR+lg27y8++672f789NNPyyOPPGJeW7RoeGIC3TohxDRiAIheJ0+ezHY7c+ZMns9NTU0127xm0JQqVcpsU1JSfD4fnQq8ePFiKV++vBn0mpN2J+lzduzYIadPnzatLDr+RZ8/ceJEeeqppyRcaDkBAAdz+oydpAMVpXCJWL+OcSEt3WyrVq2abf+oUaNk9OjREg4rV66Uxx9/3Ewvnjx5sqmfkpM+7kq7dB599FFp1aqV6VbSAbs6ODfnzxUKtJwAABAA+/btM4NcrZvOosmLNTtHWyzcOXXqlNlqMTVvbdmyxXTj6Iwg7Zrp0qWLV6+/5ppr5M4775Tz58/LkiVLJBxoOQEAIADKli1rbp6wpv0mJye7fTz54v74+HivzmH37t1y++23y7Fjx0yrzcCBA8UXOg1ZHTx4UMKBlhMAAELMmuKbmJjo9vHEi/u13omnNEjcdtttZqtdNtqt5CsNN65jX0KNcAIAQIi1aNFCypUrJzt37pQNGzbkejwhIcFsdVqvp2FCy9Dr8e6//3558803fT43Hci7YMGCfKc6BxvhJESYqQMAsOjifAMGDDD3tey8NcbEKl+/adMmMzC1SZMmWft1gKrWH8k5lkXHrXTs2FE2b94sPXr0kEmTJpmBsPnZtm2bfPTRR7lmFOmaPPfcc48ZP6OtOxqiwoExJ0CI2H3xPxb9cy6nz9gJl+HDh5upvqtXrzZjPFq2bGnqmqxdu1bi4uLMLBtXR48ele3bt+caB/Lss8/Kd999J0WKFDF1SR544AG37zd16tSs+4cOHZL77rvPdP9cf/315v20iqyWy9fpy1pn5fPPPy8w5AQL4QQAgDCIjY2VpUuXyssvvyyffPKJzJkzRypUqCB9+/aVF154Ic8CbXmND9EKr3qcvLiGk7p168oTTzxhFh7UFpf//Oc/Urx4cbNfu5I0tFxyySUSLoUyMzMzvX1RWlqa+TA/++wzU+VOP8z27dubD/OKK67w+iR0iWhdnEiXdNbkplOnNEV27drV6yIwWvhG+/Gqvjfa77nrgUS3DlQgW06KJ8UE9EOl5cTZ7NhykpGeLjtfGmam5Xo6CyYc1wqtc7LvkdFBOU+n8rrlJD09Xdq0aWPSVqVKlcxcag0XU6ZMkfnz55v9NWvW9Ph4X331lVngSAOPDry54YYbTILTJPfBBx+EtUIdwuO2y7dl+/PXh8K7xgMAwObhZMyYMSaANG/eXBYtWpRVSEYH8AwePFj69esny5Yt8+hYOiBHW0e0peTrr7+WG2+8MeuxCxcu5DnFCs4JJjn3EVTgiXI7/zvI70St4nxgQLSHE602p6OFldbft4KJ0hK3//znP2X58uVmQI3rCOO86Gu0JUZXQ3QNJqpw4cJmkA6cHUwKeg5hBTkDSc59BBQg8ng1lXjVqlWmT61WrVpmwaCctHtGzZs3r8Bj6TQlHWOiXUAdOnSQaMZ4k8AEk7xeZ93gLBo8XG/5PQ9AFLecbNy4Md+iLNZ+nZ9dEO360a4bbTHR+v2zZs0y4UdHG9evX1969uwZ1pHCiDy0qkQ/gkbwMJ0YERtOdGaOymt6k7Vf52kX5OeffzZb7RrSud06jiXnvG2tkHfLLbd4c4qIQLR6INhhhO4dIIq7dVJTU822ZMmSbh+3avBrARdP52V/+OGHZmCszs3+/fffTYGZe++919zXlRT379+f73G0up1OCXO9IXIEM5gQeuCKVhcgcoStfL126Sjt0tEpw7169TLdOFoARkvq/vGPfzTjWyZOnJjvcbTeis5Vt25Vq1YN0U8AfxEeEOowQUABojCcWLNztI6/O9baADo12NNj6fbuu+/O9bguXKR09k9+dI0BDTHWTQfawv5CFUwIQMiJgAJE2ZiTatWqmW1ycrLbx6398fHxBR7Leo4e013t/urVq5vt4cOH8z2OltvVm10xUyc3AgPCHSAYgwJEUcuJrlCo8iqOZu1v0KBBgceypiJbY09y0jEnyrWWCiKf04OJXcOqE0vX04KSmxP/HiAKwokunazjOnbu3CkbNmzI9bjOrlG6aFBBdArxpZdealZG1EGwOVndOe7qqQDecHogQt4IKEAUhJOYmBgZMGCAud+/f/+sMSZW+Xqtb9KqVats1WG1omy9evXM2BBXuqyzVojVdQf1WK6zbHQJaV09Ubt7Hn74YX9+PtgIIQGeIDAA8HptneHDh5vwsHr1arNysNYo0boma9eulbi4OJk8eXK25x89etS0jBw8eDDXsXRRP10uWo+ns3R00T99vtY80WJsL774ojRt2pRvKQqEO5jo+1PqHu4w/gSIgqnEsbGxJlCMGDHC1DuZM2eOCSd9+/Y1Y068WZG4WLFi8uWXX8qrr74qFStWNOXsdTVibX3REvjDhg3z9vRgQ+EOJogc4Wo1obUGsJdCmdqvEkW0e8jUO3lvtBQuERvu07HtAEgnBhO7tJzsSY4LyHGKJ8VItA2EDHdIYJFAkZT43LMnwyUjPV12vjTMlIkoW7asba8VF9LSZd8jo4Nynk4VtiJsTkAwsU8wsVtQgv2CiV3OAQDhBACycXpAsUsrGpyNlhM4qpXCrucFAPgfwgkCjgCASG+tsNv5AE5DOIHjEJ7gCQIKED6EEwQUF354ixAAICfCSZA4faYOEA0ITkB4EE7gyFaTcJ4rwfV/uPjbEzN2EG6EE8DhuBDljwAFhB7hBI5rNYnkcwYAJyCcAAiLSGqRiKRzBaIB4QSOboGI5HMHgGhFOAkCBjwC+aMlAkB+CCcA4AECFRA6hBOI07tFouFnAAKNWVwIJ8IJgJCK5BaISD532FNaWpqMHDlS6tatK7GxsVK5cmXp16+f7N+/3+tjHTt2TB5//HGJj4+X4sWLm+0TTzwhx48fz/M1GRkZ8uabb8q1114rJUqUkLi4OOnRo4ds3bpVwolwAp9FU4tDNP0sACJDenq6tGnTRl544QVJTU2Vzp07S9WqVWXKlCnSqFEj2bVrl8fHOnr0qDRt2lTeeecdKVq0qNx1111SpkwZefvtt6VZs2by+++/53rNhQsX5O6775ZBgwZJcnKydOzYUa655hpJSEiQ66+/Xr7//nsJF8IJgJCJhpaHaPgZYA9jxoyRNWvWSPPmzWXHjh0yffp0Wbt2rYwdO1aOHDliWlA8pS0kv/76q3Tt2lW2b99ujrVlyxYZOHCgObYGkJwmT54ss2fPljp16si2bdtMKFm2bJnMmDFDTp8+Lb1795bz589LOBBO4JNobGmIxp8JgD2dPXtWxo8fb+5PmDBBSpcunfXYoEGDpEGDBrJ8+XJZv359gcc6ePCgfPrppxITEyMTJ040LSeW119/3XTVfPzxx3L48OFsrxs3bpzZvvbaa3LZZZdl7e/WrZvceeedJuzMnTtXwoFwEmBMIwaiv8Uhmn4WhMeqVavkxIkTUqtWLdOFk1P37t3Ndt68eQUe69///rfpomnZsmW2kKF07EmnTp3M2JIvv/wya//u3bvNuBIdZ6LdOf68fzAQTuA1WhgAZ2DGTvBs3LjRbBs3buz28cYX92/atCkox7JeU79+fSlWrJhf7x8MhBPABcELnqL1BP5ISkoy2ypVqrh9vMrF/Xv37g3KsQL5/sHwv44pwANcvAPX/bcnOc4xvxlzIYcTnDx5MleXit7c0dk5qmTJkm4fL1WqlNmmpKQU+L6+HCuQ7x8MhBPATQD7+lA9PhfAAWL2xUiR2Bi/jpGRfsFsdRqwq1GjRsno0aP9OrZTEU4ABFU0t5roz3ailvvfjOE8+/btk7Jly2b9Oa9WE2XNztEpu+6cOnXKbLVWSUF8OVYg3z8YCCfwGF06AJA3DSau4SQ/1apVM1stfuZO8sX9WuU1GMcK5PsHAwNiA4hpxNGDIBYY0dxq4qSfEYHXsGFDs01MTHT7eOLF/VrvJBjHsl6jhdrOnTvn1/sHA+EEHuFiDTgT04mDo0WLFlKuXDnZuXOnbNiwIdfjCQkJZqs1SgrSvn17KVy4sKxcuTJXobUzZ86YWiVFihSRDh06ZO2vUaOGXHXVVWZtnwULFvj1/sFAOAEAP9F6Am9pNdcBAwaY+/37988a42FVbt20aZO0atVKmjRpkrVfK8rWq1dPhg4dmu1YlSpVkl69epmqs48++mi2kvNDhgwxpfDvvfde+cMf/pDtdVZJe32Oa6iZNWuWfPHFF1K7dm2z3k84MOYEBaLVBAACb/jw4bJ48WJZvXq1Wd9GK7xqXZG1a9eakvO69k3Oxf103RwtV5/TW2+9ZdbpmTlzpgkwunDfTz/9ZLpt9NhWqXpXunaPVo3V9XX0Nbfeeqt5Dy2br5VjteS9ayn8UKLlBMgDocw/TmtNcNrPC//FxsbK0qVLZcSIEabeyJw5c0w46du3rxnzUbNmTY+PVbFiRbOKsC70py0oGji0PP5jjz1m9leoUCHXa7QrSBf504UGK1euLPPnz5fNmzebtXXWrVtnVjMOl0KZmZmhqcIUwiI42o9X9b3RUrhEbEjfO1oHxDr5Ih3Meif+FmErnuRfbYZgjydw4sU6WqcVp8QXCsv7ZqSny86XhpmLrKezYLy9VtQa9pIUiY217Xk6FS0nyJeTgwkAIDwIJwESra0mADznxNYiIBgIJ8gTrSYAFNOJEWqEEyAfBDR4i9YTwH+EE7jFRRn+4AINwB+EEwAIMMIZ4B/CCXKh1QQAEE6EE8ChM7wY5BhctJ4AviOcRMFFBgCAaEI4QTZ06eTGZwLQ0obQIpwACCi6M/gsAH8RTgAAgK0QTpCF7gsAgB0QTgAgiOjmArxHOPETM3UAAAgswgkMunTyx+cDAKFDOAEQMHRhRPfnQuE+hArhBAAA2ArhBHRZAABshXACACEQLV07QCgQThyOgZ4AALshnPiBacTOQpADgNAgnAAICLot+IyAQCGcOBgtAQC8xXRihALhBABCiBYmoGCEE8CB+O0XgJ0RThyKLh17YFA1AORGOAGAEKNrB8gf4cRH/MbrTLQ4ucfFFkAgEU4ciAssAH8wZgnBRjgBgDCgtQkIQjhJS0uTkSNHSt26dSU2NlYqV64s/fr1k/3794s/fvnlFylRooQUKlRI2rZt69exAACIZqtWrZIOHTpIhQoVpHTp0tK0aVOZNm2a18dZv369jB49Wm688UYpX768xMTESNWqVeXee++VTZs2uX3Nnj17zLU6r9vll1/u889V1JcXpaenS5s2bWTNmjVSqVIl6dy5sznJKVOmyPz5883+mjVr+nRCDz30kJw5wwJZwUKXDgBEh5kzZ0rPnj3lwoULcvPNN0vFihVlyZIl0qdPHxMo3njjDY+Oc/78ebn++uvNfQ05GlBKlSolP/74o/zrX/+SGTNmmG337t3dvv6yyy6T9u3b59pfrly50IaTMWPGmADSvHlzWbRokUlraty4cTJ48GDTgrJs2TKvj/uPf/zDvE4Dyt///ndfTg0AIqpr50St4uE+DUSg33//3VxrMzIyTEjp2rWr2f/bb7/JTTfdJGPHjpU77rhDWrdu7dHx/vjHP8qzzz5rXlOkSBGzT0OP9pC8+OKL5r30WBqAcqpXr55MnTo1vN06Z8+elfHjx5v7EyZMyAomatCgQdKgQQNZvny5aSLyhn6gTz31lNx2223Sq1cvb08LAADH+PDDD+XkyZOm58IKJlYrxmuvvWbua0DxRNGiReX77783x7KCiSpcuLC88MILcuWVV0pKSoosWLBAQqWwL/1bJ06ckFq1akmjRo1yPW41+8ybN8+r4z7++ONmHMvEiRPF7iJ1GjFdOnyOwcDATiD0FlwMCu66Wjp27GjGgi5evNgMw/CHjh3RRgd14MABsW042bhxo9k2btzY7ePW/rwG0Ljz5ZdfyvTp02XYsGFSu3Ztb08JACJWpIY7phOH18Z8rsU6mLV+/fommOzYscPv99q1a5fZ5jXAVXs+Ro0aZYZkaA9IQkKC6WXxh9djTpKSksy2SpUqbh+39u/du9ej4506dUoeffRR02z09NNPe3s6AAA4ysmTJ00PRkHX4nXr1plrsdXy4Ytvv/3WDNPQwONu0Kvatm2bPP/889n2VatWzQyk1dlDIWk5SU1NNduSJUu6fVxH+Crtn/LE8OHDzYf3/vvvmx/eWzqzR78o1xtyo0sHAIIr57UoWDNPUy9ehwN5LXZHfwYdCKuefPJJMzvXVfHixeWRRx4xE1m09USf/91335mpzdqQ0a5dO48bKgIyWydQNNW98847ct9993k8ojinl19+WZ577rmAnxsAhAqzdsKnzL5MKRKT6dcxMs7+9/VaF8SVdnVo7RB3unTpIlu3bvXqfaZNm+ZzS4S3dBZQ7969Te0xfc+cLSNKw0rOcaI33HCDGQ+jr/3kk0/kpZdekg8++CD44cSanXP69Ok8u2lUmTJlCpxX/Ze//MUUe/F0LrY7Q4cONbOELJrccv4FcTpaTQAg+Pbt2ydly5bN1rKQl927d8v27du9Ov7pi9dd11myus/1Pb29FudFW0S0bpkOudCw4W3Pho4h1XCycOFCn97f63Ci/UgqOTnZ7ePW/vj4+HyPo8/bsGGDGWBz9913Z3vs+PHjZqv9XFaLSl51U/TLz+8vAAAAoaAhwV1QcEevf/68T7ly5cy4E72WXn311T5fi9155plnZNKkSeYX/a+//tptbZOC1KlTx2wPHjwoIQknDRs2NNvExES3j1v7PR2Ac+jQIXNzR0OK1kyxk0idRozAt0Z9faie4z/WSJ1pYkd07cDba/GKFSvMNTdnODl37pxs2bLFTCfWJWa8oTVSXn31VfnDH/5ggomvPRHHjh3LNvYl6ANiW7RoYRLbzp073SY/nUKkOnXqlO9xqlevLpmZmW5vS5cuNc+59dZbs/bBN3TpAAgWphOHT8eOHbNdc11pd4xOI9b16TSgeEpbS3TWrA630O4Y7dLxlVatza/sSMDDifY7DRgwwNzv379/Vr+WVb5e65u0atVKmjRpkrVfK8pqeVsdHwIAcI+WKHjqwQcfNN07c+fOlVmzZmXtP3z4sAwZMsTc1+VkctJrsd5yLtKrIeevf/2rGc+itceuu+46j8KMTiPOSc9Hu4asnBCy2To6/Vcrz61evdr0K7Vs2dJMF1q7dq3ExcXJ5MmTsz3/6NGjZuCPr31PAADgf3SBPr3W9ujRw1SJ1fGZl156qbk265AInSjibhasNQhXu35cA43OrtG1dGrUqGFm17ibYXPXXXeZm0UXA9TCazqMQ7uP9PU///xzVmDRgmw6Kylk4USbibTrRafx6mjcOXPmmA+qb9++pg5/XkVhEFp06QBA9OrWrZsZd2ItxqtVWXX8ifZu6MrEntIZP1ZF182bN5tbXsMxXMOJzrjVBgkd4qGLAOsSNPpnXetHZ/tot5KvCmVG2YAOnUqsY2KqvjdaCpfwvK8tGgfEEk6CK1ADYvckx/n0uuJJ3hctDPRYAbohgiOSVipOiS8UtGNnpKfLzpeGmVkpns6C8fZacd2fX5QiMf5dKzLOpsuGj54Nynk6lddjTpwskoIJIkeo/14xiBGA3RFOohStJnzGABCpCCcAfEKXTvBE0mdLSxyCgXACAABshXAShejSAQBEMsIJANhQJHXtAIFGOAEAALZCOImyacR06QAAIh3hBABsiq4dOBXhBADgF6YTI9AIJ1GELh0+81DhN3oAwUQ4AQAbIwjCiQgnAADAVggnUYIuHQBAtCCcRNE0YgDRia4dOA3hBADgN2bsIJAIJ1GALh0AQDQhnABABKBrB05COAEAALZCOIlwdOmEH98BAAQW4QSAV+heABBshJMCMI0YgF0QDOEUhBMAQEAwnRiBQjiJYIx1AABEI8IJAEQQunbgBIQTAABgK4STCEWXDgAgWhFO8sFMHQB2RNcOoh3hBAAA2ArhJALRpQPArphOjEAgnABABKJrB9GMcAIEgFNas7ggAvayatUq6dChg1SoUEFKly4tTZs2lWnTpnl9nKlTp0qhQoXyvN1zzz15vvann36Su+++W+Li4qREiRJy7bXXyltvvSUXLlzw+ecq6vMrERZOuQgiOGhyB6LHzJkzpWfPniYE3HzzzVKxYkVZsmSJ9OnTRzZt2iRvvPGG18ds2LChXHfddbn2N2vWzO3zv/vuO7n11lslLS3NBKPq1avLihUr5Mknn5TVq1fL9OnTTbjxFuEEACK4JetEreLhPg2Ewe+//y79+vWTjIwME1K6du1q9v/2229y0003ydixY+WOO+6Q1q1be3Xcu+66S0aPHu3Rc8+dOye9e/c2wWTcuHEmkKjU1FS5/fbbZcaMGaZVp2/fvl7/fHTr5IFpxAgl/r4B8MaHH34oJ0+elM6dO2cFE3XZZZfJa6+9Zu5rQAmm2bNny+7du01rixVMlHYvjR8/3q9zIJxEELp0AABqwYIFZtu9e3fJqWPHjhIbGyuLFy+W9PR0Ccc5NG7cWGrWrClbtmyRPXv2eH1swgkARDA7DlJmbFPwbdy4MSsE5BQTEyP169c3wWTHjh1eHXf9+vXy1FNPycMPPyyjRo2S5cuX+3QOrvt1/Iu3GHMCAEAEOXnypJw4ccLcr1Klitvn6P5169bJ3r17pUGDBh4fe/78+eZmef7556VVq1ZmYKt2GblKSkoq8ByUnoO3aDmJEHTpAID9Q4Pr7cyZ4LRqpaamZt0vWbKk2+eUKlXKbFNSUjw6ZqVKlcxA2B9//NEEn0OHDskXX3wh9erVM60nOrhWB9+6O49AnYMrWk5gO93LJpptwkn3TYUAsmPWju/K7j4jRYt6P9XV1fnz/w0hVatWzbZfu0XymvnSpUsX2bp1q1fvM23aNDNdNxjatWtnbpayZctKp06d5JZbbpEmTZqYVpjPP/9cevXqJaFAOIEtg4nrfUIKgEiwb98+c1G3FC+e9zRvneWyfft2r45/+vTprNkwrvtc39Ny6tQpsy1Tpoz4Q9/rsccekwEDBsjChQuzhRN97NixY1nnFchzoFsnAqZ1OqVLxzWYeLIfyEvMtuSsGxAqGhJcb/mFkw0bNkhmZqZXt9YXa5boscuVK2fuJye7/ztu7Y+Pj/f756pTp47ZHjx4MNv+atWqBe0cCCewhYICiD5OSIEncgYSpwQUu83aYcZOcDVs2NBsExMT3RZH0ym8Op24bt26fr+Xto64jiHx5Bxc93szINdCOEFYeRs6CCjIS34tJU4JKHCOjh07mm1CQkKux3S2jU4jbtu2rQko/tIKtO6mDOd3DjqwdteuXWZKs5a09xbhxOaiuUvH16BBKwpcedp9QzcPosmDDz5ounfmzp0rs2bNytp/+PBhGTJkiLk/ePDgXK/T2Td6279/f7b9L7/8shw9ejRXC8xzzz1nytDrgn73339/rkG9NWrUMPVO3nzzzWxjTfr375/nOXiCAbEIi0C0gDBgFr60iOhrztZzX5ch0jFrxzkqVKggkydPlh49epgKrToe5dJLLzVVYY8fPy6DBg1yu66ONQhXg4erYcOGmSBy/fXXm1lHOhVax8UcOHDAtL58/PHHcsUVV2R7TbFixcx+baHR99NaKDq+ZOXKlWZ8ip6XLkLoC1pOEHKB7pqhq8d5/G0FoZsH0aBbt25mBWCdAqzdKF9++aXUrl1bpk6d6vWaNiNHjjQrG+uMI22N+eabb0z9Eq0UqyHFdf0eVzfeeKP88MMP5lx+/fVXUxtFg5MuBOjrisSKlhMbz9SJti6dYIYIO7Si6Pf19aF6Eq3sMuAyUMEimltQ4BwtWrSQr776yuPn66wfd7TVxFfXXHON23En/qDlBCERqtYNWlGiVzDGjERjC4pdQiTgD8IJoi4wMGA2+gQzRDBQNniYTgxfEU5sKlq6dGjJQKQEh2hsRQEiFeEEQUMwQaQhoAD2QDhB1Har2OEcEHmiIaAw7gSRjnBiQ9HSpQNEakhgHAoQXoQTm04jjlS0VgAA/EU4AYAo7OKhaweRjHBiM5HcpWPHVhM7nhPgJEwnhi8IJwgIQgCitcXCTucCOAXhBH4jmCDaRWpAoWsHkYpwAgAAbIVwYiORON6EVhM4RaS2ngCRiHByEdOIozeYRMp5wv4iMaDQtYNIRDgBYCt2DwB2Pz87YsYOQhZO0tLSZOTIkVK3bl2JjY2VypUrS79+/WT//v0eH+P48ePyySefSK9evaRGjRoSExMjZcqUkWbNmsnbb78t586dE6eItC4dWiMAALYKJ+np6dKmTRt54YUXJDU1VTp37ixVq1aVKVOmSKNGjWTXrl0eHeeNN96Q3r17y/Tp0+WSSy6Rrl27StOmTWXjxo3yxBNPmPc4ffq0L6eIICKYwOkirfWErh04IpyMGTNG1qxZI82bN5cdO3aYcLF27VoZO3asHDlyxLSgeKJUqVIyZMgQ2bNnjyQmJspnn30mS5Yskc2bN0u1atXk22+/Ne8FAHYTaQEFiOpwcvbsWRk/fry5P2HCBCldunTWY4MGDZIGDRrI8uXLZf369QUea+jQofLqq6+aIOKqTp068sorr5j7n376qUS7SOrSidRWk0g9bwBwIq/DyapVq+TEiRNSq1Yt04WTU/fu3c123rx5fp1Yw4YNzfbAgQN+HQeBwwUewRZprRGRdL507SCqw4mOB1GNGzd2+7i1f9OmTX6dmDVu5fLLL5dgYxoxACcEFCBqw0lSUpLZVqlSxe3j1v69e/f6dWI6W0fpYFuEH60mAPzBdGJ4o6hXzxYxs3NUyZIl8xzkqlJSUsRX77//vixevFjKly8vzzzzTL7PPXPmjLlZTp48KZEkEsabEEyAgltPztZz/wsbgCgowrZy5Up5/PHHpVChQjJ58mRTPyU/L7/8spQrVy7rplOaETgEEyB6uncYd4KoDSfW7Jy86o+cOnXKbLWYmre2bNliunF0RpB263Tp0sWjGT86QNe67du3z+v3hTMQtOwtEi7uAGzarWNN+01Odv8fibU/Pj7eq+Pu3r1bbr/9djl27JiMHj1aBg4c6NHrihcvbm6RyO5dOlzMAe/QvQOEqeXEmuKrRdPcsfZrvRNPHTx4UG677Taz1S6dUaNGeXtaCDCCCRCdLUB07SAqw0mLFi3M2I6dO3fKhg0bcj2ekJBgtp06dfLoeNpS0q5dO3O8+++/X958800JJaYR50YwAaI7oABRF050cb4BAwaY+/37988aY6LGjRtn6pu0atVKmjRpkrVfK8rWq1fPjA9xpeNWOnbsaMrV9+jRQyZNmmQGwjqB3bt0APiHgJIb04kRtDEnavjw4Waq7+rVq02p+ZYtW5q6Jrq+TlxcnJll4+ro0aOyfft2023j6tlnn5XvvvtOihQpIkWLFpUHHnjA7ftNnTrVl9OED2g1AaKfdu2cqBWZY/XgDD6Fk9jYWFm6dKmZxvvJJ5/InDlzpEKFCtK3b1+zUnFeBdrcdemojIwMc5y8EE4ARCIGyAIhrnNSokQJef755+XXX381RdC0VWTKlClug4nOvsnMzMwVMvTPur+gW7Sxa5cOrSYAEFlWrVolHTp0MA0EWuqjadOmMm3aNK+PU716dTOsIr9bzZo1s71mz549+T7fn+VnfGo5AQBEdusJXTuRb+bMmdKzZ0+5cOGC3HzzzVKxYkVZsmSJ9OnTx4z/fOONNzw+li7aq0Mw3Fm+fLkJIjqEw53LLrtM2rdvn2u/Tp7xlaPDCTN1nNdqoj9nwkn3i1YCQKT4/fffpV+/fmZYhIaUrl27mv2//fab3HTTTTJ27Fi54447pHXr1h4dL68go8HH6hH585//7PY5OuEl0MMvbFe+PtrZtUsHCKdon9kS7T8fQu/DDz80a8lpVXUrmFitGK+99pq5rwHFX9oSo8M2rrjiCmnTpo2ECuEEjmk1ARB+TCcOjAULFmR1x+SkJTp04orOqk1PT/frfT7++GOz/dOf/iSFC4cuMhBOHI5gAji39YRqsZFr48aNZtu4cWO39cjq169vgsmOHTt8fo+0tDSZPXu2uX/vvffm+TztStLK7g899JA89dRTphirrpHnD0ePOQEAINKcPHnSLHSr8irdofvXrVtnapB5s5yMKy0TkpKSYl6f3zG2bdtmZu/mXIdvxowZZvaQL2g5cfB4E1pNgNCyY+sJAhsaXG9aZiMYUlNTs+6XLFnS7XNKlSplthoufPXRRx/lOxBWF9195JFHZNmyZab1RH9mLayqU5uTkpLM0jQajnxBywkAOHhqsdOnFMfsOCBFC8f4dYzCF/7bhVG1atVs+7WrQ+t8udOlSxfZunWrV+8zbdo0n1sivHX48GH5+uuvzTgTHW/iTqVKlWTixInZ9t1www1mPEzv3r1NcdWXXnpJPvjgA6/f37HhxOnTiGk1sZc9yXHhPgUAftq3b5+ULVs2W8tCXnbv3m2WdfHG6dOnzVaLrbnuc31Pi7XuXZkyZcQXn332mZw/f15uu+02qVy5stevHzZsmAknCxcu9On96dZxaJcOnCkl3hkLa9qd07t3onXGjoYE11t+4WTDhg0eVUjPdLlZNUv02FaBs+Rk93+XrP3x8fF+zdLJbyBsfnTdPZVzTT1PEU4ciFYTAK6YtRN5GjZsaLaJiblLQZw7d062bNliphPXrVvX62PrDJ8ffvjBjGdxraHiDWvtPGvsi7cIJw5DMAHswemtJ/BPx44dzVan7eY0f/58M424bdu2JqD42mqi42Jcu5C8oVVr85rq7AnCSQjQpQPAHQIKfPXggw+a7p25c+fKrFmzsg1kHTJkiLk/ePBgt6Xm9bZ///48j/2vf/0r31k6lkmTJplpxDnp+TzzzDPmfv/+/cUXjh0Q60S0mgTX14fqSTTTGR00/0cvp8/aiTQVKlSQyZMnS48ePUyVWB2Pcumll5qqsMePH5dBgwa5XVfHGoSrXT/urF69Wnbt2mVWFNaWl/xoiNHCa1oDRbuPdB2en3/+OSuwaEE2bX3xBeEEAMLIblOLETm6desmK1askDFjxsiaNWtMVdarr75aBgwYYFYm9oXVpdOrVy8pUqRIvs/9y1/+InFxcWZw76JFi0xFWf2zjlPR+icFhZv8ODKchHIasV26dGg1AYDo06JFC/nqq688fr7O+smP1i3JWbskL1rLRG/BwJgTByCYAPZml7Enoey2i9bpxAgMwgkA2IBdAgpgB4STKO/SodUEABBpCCdRjGACRBY7tJ4wIwt2QDgBAAC24rhwEqqZOuHu0qHVBIhMdmg9AcLNceEEAOwu3AGFrh2EG+EkCtFqAiASMJ0YeSGcRGGXDoDIF+7WEyCcCCdRhlYTAIFA1w7CiXASRQgmQHSh9QRORTgJMLp0AADwj6PCSSgX/As1Wk2A6BTO1hO6dhAujgonQLB8faieIz7cE7WKh/sUHCmau3eYsQN3CCdR0KVDqwkAIJoQTgAgAoSr9YSuHYQD4STC0WoCAIg2hJMAYZYOgGCL5rEngCPDSTTO1KHVBEAo0LWDUHNMOAHsak9yXEjfLyW+UEjfD4FF6wmcgHASoWg1ARAtmE6MnAgnAcB4EwChROsJoh3hJALRagIg1Bh3glAinABABKL1BNGMcBJhXTq0mtiPU0rXWyhhDyDYHBFOomUaMcEEQDhbT+jaQag4IpwAloSTjfkwbOhsvSrhPgWEGTN24IpwEiFdOrSaAHCHsSeIRoQTAIhwoQwodO0gFAgnEYBWEwCAkxBObN6lQzAB4Am6dxBNoj6cRMtMHQCwC7p2wu/UqVPy0UcfycCBA6VZs2ZSvHhxKVSokIwePdqv486bN09atWolZcuWNbfWrVvLggUL8n3NTz/9JHfffbfExcVJiRIl5Nprr5W33npLLly44PN5FPX5lQg6Wk0AeNt6Eskzn3TGDgtTeuaXX36R++67L6CfvwaKJ598UooWLSpt27Y1gWfRokVyxx13yLvvvisDBgzI9ZrvvvtObr31VklLS5OmTZtK9erVZcWKFeY4q1evlunTp5vQ5K2obzkJBtbSgVMLsFkoxGZfdO84Q5kyZeSBBx6Q999/X9avXy/PP/+8X8fbvn27/O1vfzOBRMPFV199JXPmzJENGzbIpZdeasLGr7/+mu01586dk969e5tgMm7cOFm7dq0JIxqcmjdvLjNmzJB//vOfPp0P4cSmaDUBYGd07YRXrVq15MMPP5SHH35YGjduLMWKFfPreG+//bZkZGTIX//6VxMsLHXr1pVnn31Wzp8/b57javbs2bJ7925p2LChCS+W0qVLy/jx4839sWPH+nQ+hBMbIpgA8AetJ/CWNa6ke/fuuR6z9ul4FE9fo4GpZs2asmXLFtmzZ4/X50M48RJdOpHLjtVh9yTHheV97divH8ljJYBIdvz4cUlKSjL3GzVqlOvxqlWrSsWKFWXv3r1y8uTJrP0bN27MCiLuWPs3bdrk9TlFdTiJxJk6tJoAiJTWE7p2okPSxWByySWXSKlSpdw+p0qV//7yoAEl5+usxzx5jWNn62RmZprthbR0OX/qTMCPn556ToIptZDvU68Q+u/O379j+vfUVxnp/v1dyTj7338rvjp/PvD/vgpfOBvwYzrZ+fO+//0K1d+jXMdL/1+r3oUz6dn+Xw+G85lnRS4E4Bgi2VoVlA4u1Zvdpaammm3JkiXzfI4VWlJSUjx+nbvXODacWB/C/kGvyP4gHH+VBNeYIB/f2RZFyDGBADkcgk/y29D8v16uXLmAHjMmJkYuv/xyWXZoWkCOp4NAtfvD1ahRo/KsO9KlSxfZunWrV+8xbdo0M13XCaIunFSuXFn27dtnpln5MrcauelvA/qPTj9XLcqD8OM7sR++k8DTFhMNJvr/eqDFxsaamSZnz54N2LnmvObk12qi763Td71x+vRpCQYNVgUdX4u+Kb22ur7u2LFjeb7O3WscG04KFy6cZ/8X/GNVDIR98J3YD99JYAW6xSRnQNFbOGj9ELuoVq2a2WrQ0EDhbtxJcvJ/xzDFx8dne52+Rh9r0KCBR6/xVFQPiAUAAPkrX758VkD58ccfcz2ureZHjx41IcP1F1Stb6ISExPdHtfa7y64FIRwAgCAw3Xs2NFsExIScj1m7evUqZPHr9GQs2vXLqlfv74pae8twgkKpP2mOrArEkadOwXfif3wnSAS1KtXz9z2788+ZeTxxx+XIkWKmHL4a9asydqvpehffPFFs96OPifnoN4aNWqYeidvvvlm1n7tGurfv7+5P3jwYJ/Os1BmMOdoAQCAoOjSpYscPHjQ3D9w4IDpfrniiiuyxl1WqlTJlJh3ZQ3a1QG5OVs0NGAMGjTIBJHbbrvNzGjShf907Zx33nnHrICcky7up4sE6nN0dWTt+lm5cqU5L60c+/nnn/s0OYVwAgBABKpevXq+Bc40KOQsHZ9fOLFK1L/++utZY0+0YuyQIUPMysR5+emnn0zr+rJly0yria77o4sSakuLTlLxBeEEAADYCmNOAACArRBOkMuqVaukQ4cOUqFCBVNkRysSamVCb02dOtU0IeZ1u+eee/j0L9L+2pEjR5rlybXughad6tevX65Ba57QugPanKpNujpIU7dPPPGEWdwLof9OtOk8v38H27Zt42sBor0IG/wzc+ZM6dmzp1y4cEFuvvlmsxLlkiVLpE+fPmZlyTfeeMPrY+pc+Ouuuy7Xfh08BZH09HRp06aNGSGvA9g6d+5s+omnTJki8+fPN/t16XFPaC2C5s2by6+//mpec9ddd5n+4Lffflu++uor+e6770zoROi+E4v+Gwp1kTEgYulsHUD95z//ySxbtqzO3sqcOXNm1ody6NChzNq1a5v9S5cu9fjDmjJlinnNqFGj+IDz8eyzz5rPqXnz5pkpKSlZ+8eOHWv2t2rVyuPPr3fv3uY1Xbt2zTx37lzW/oEDB5r9ffr04bsI8XcSHx9vXgPAc/yLQZZXX33V/CfauXPnXJ/KrFmzzGN33HGHx58Y4aRgZ86cySxXrpz5bBMTE3M93qBBA/PYunXrCjzWgQMHMgsXLpwZExNjAqWr9PT0zLi4uMwiRYpk/vbbbx5/h04UyO9EEU4A7zHmBFkWLFhgtjo3PSetBKj97osXLzZN3gjc+J4TJ06YqXc6ZS8n67vQ6X0F+fe//22641q2bCmXXXZZtsd07IlWd8zIyJAvv/ySry9E3wkA3zDmBFm0yp9q3Lhxrk9Fi/FoGeJ169bJjh07vForYf369fLUU0+ZVVt1iXLty2/VqhWffAGfuet+He8TiGNNnjzZo2M5WSC/E1daO2Lnzp0mKF5zzTWmgFZcXFwAzhiIPoQTGBoc9LdFldeqzrpfw4kW/fEmnOgAQr1Znn/+eRNOpk+fnus3fKdJSkoq8DNX+RVaCsaxnCxYn6MWsnL15JNPyrvvvmtmAAHIjm4dGKmpqVmfRMmSJd1+KtYy2ikpKR59ajrLYfTo0abSoAafQ4cOyRdffGHWdVi+fLmpOKjdDE5mfe6B+MwDeSwnC/TneOedd8qsWbNMmDl9+rRs2bLFlAg/c+aMPPjggzJ37twAnj0QHWg5iSLaTLx161avXqP1S7SOSTC0a9fO3Cy61LaOe7jlllukSZMmphVG113o1atXUN4fsANdk8SVdumMHTvWhPSHHnpInn76aTNVGcD/EE6iiK6VsH37dq9eo7/JKS225rpPg0ROumaCKlOmjF/nqe/12GOPyYABA2ThwoWODifW5259D/585oE8lpOF6nPUtUeGDx9u/s1qDRVflpUHohXdOlFkw4YNOjXcq1vr1q3NazWMWMWgkpOT3R7f2q8VR/1Vp04ds7VW1HSqatWqBewzD+SxnCxUn6MuiKYzgpTT/x0AORFOkK2Sq0pMTMz1qZw7d870let0Yi3n7S8tse7af+9U+X3mrvs9GYAcyGM5WSg/R/4dAO4RTpCtlolKSEjI9anobButb9K2bVsTUAJRJj+/6ZpO0aJFC9NipVNMteUrJ+u70LE6BWnfvr35bXzlypVy+PDhbI/p4Euty1GkSBGzbhJC853kR5cV0C4dHXir408AuPChcBscVr5eK4rmV77+yiuvNLfk5ORs+1966aXMI0eOZNt39uzZzNGjR5tjlShRItdrnFwq/cYbb8xMTU0tsFT6u+++az7vZ555Js/y9d26dctWvv6xxx6jfH0YvpMFCxZkLlmyJNfxN27cmHnVVVeZY+l3AyA7wgmySUhIMCXQCxUqlHnLLbdkdu/ePbN8+fLmP9FBgwa5/bT0Mb3t3r071/7ixYtntmjRIvOee+7J7NChQ2blypXN/tjY2GwByMnS0tIymzVrZj6XSpUqZfbo0SPrz1pyfufOndmer2sV5bVOjobBWrVqmcd127Nnz8z69eubP9epU8cEUITuO7H2awn7O++80/w7aNq0aWbRokXN/tatW2eePn2arwTIgW4dZNOtWzdZsWKFmQKs9Um01Hnt2rVl6tSpZvqjN3S5eV3ZeN++faaWwzfffGOasB9++GHTXN61a1c+fRHTTbZ06VIZMWKE+XzmzJljamL07dvXjG/wZvVbXUX6+++/l4EDB8rZs2dl9uzZpsaMzo7S/axIHNrvRP8daZE1HXCuZfG1S0hXjL7ppptk0qRJZjmIEiVK8O8AyKGQJpScOwEAAMKFlhMAAGArhBMAAGArhBMAAGArhBMAAGArhBMAAGArhBMAAGArhBMAAGArhBMAAGArhBMAAGArhBMAAGArhBMAAGArhBMAAGArhBMAACB28v++8Tc3UKbCqAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Solve Burgers' equation using finite difference method.\n", "def u_truth(t, x, nu=0.02):\n", " [Gauss_pts, weights] = np.polynomial.hermite.hermgauss(80)\n", " temp = x - np.sqrt(4 * nu * t) * Gauss_pts\n", " val1 = weights * np.sin(np.pi * temp) * np.exp(-np.cos(np.pi * temp) / (2 * np.pi * nu))\n", " val2 = weights * np.exp(-np.cos(np.pi * temp) / (2 * np.pi * nu))\n", " return - val1.sum(axis=1) / val2.sum(axis=1)\n", "\n", "def solve_Burgers(n, nu=0.02):\n", " t = np.linspace(0, 1, n+2)\n", " x = np.linspace(-1, 1, n+2)\n", " tt, xx = np.meshgrid(t, x)\n", " tx = np.vstack([tt.flatten(), xx.flatten()])\n", " u = u_truth(tx[0].reshape(-1, 1), tx[1].reshape(-1, 1), nu)\n", " return np.stack([tt, xx, u.reshape(n+2, n+2)])\n", "\n", "N_pts = int(np.sqrt(burgers.Nd+burgers.Nb))-2\n", "txu = solve_Burgers(N_pts, nu)\n", "txu = txu[:, 1:-1, 1:]\n", "u_FD = txu[2]\n", "\n", "# plot\n", "plt.figure(figsize=(6,5))\n", "plt.tick_params(axis='both', which='major', labelsize=15)\n", "plt.contourf(txu[1], txu[0], u_FD)#, 50)#, cmap='inferno')\n", "clb = plt.colorbar()\n", "clb.ax.tick_params(labelsize=15)\n", "plt.title(\"FD Solution\", fontsize=20)\n", "\n", "None" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setting up the model\n", "We model the extended function $\\tilde u$ using some variational distribution $\\textrm{q-ED}(\\mu, \\Sigma)$. With nonlinear left-hand side function $P$, $P(\\tilde u)$ no longer follows q-exponential distribution. QPyTorch has `PDESolver` model to propagate the variational distribution $\\textrm{q-ED}(\\mu, \\Sigma)$ through the PDE by linearizing $P$:\n", "\n", "\\begin{align}\n", "P(\\tilde{u}) &\\approx P(\\tilde{u}_0) + \\nabla P(\\tilde{u}_0)(\\tilde{u}-\\tilde{u}_0) \\sim \\textrm{q-ED}(\\mu_*, \\Sigma_*), \\\\\n", "\\mu_* &= P(\\tilde{u}_0) + \\nabla P(\\tilde{u}_0)(\\mu-\\tilde{u}_0), \\quad \\Sigma_* = P(\\tilde{u}_0) \\Sigma P(\\tilde{u}_0)^{\\mathsf T} + \\delta I.\n", "\\end{align}\n", "where the reference point can be chosen as $\\tilde{u}_0=\\mu$, and $\\delta>0$ is a small nugget added to make the propagated covariance postive definite.\n", "\n", "The `PDESolver` class inherits from `ApproximateQEP` to approximate the posterior with sparse variatinoal Bayes. To faciliate modeling derivatives and differential equations, QPyTorch defines new methods including `Matern32KernelGrad`, `Matern52KernelGradGrad`, `RQKernelGrad`, and `RQKernelGradGrad` which are not available in GPyTorch." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# QEP solver\n", "POWER = 1.0\n", "class QEPsolver(qpytorch.models.PDESolver):\n", " def __init__(self, pde_lhs, power=torch.tensor(1.0), num_inducing=256):\n", " self.power = power\n", " input_dims = burgers.dim\n", " output_dims = input_dims*2\n", " \n", " # inducing_points = lims[0] + torch.rand(num_inducing, input_dims) * lims.diff(dim=0)\n", " inducing_points = burgers_X[torch.randperm(burgers_X.size(0))[:num_inducing]]\n", " batch_shape = torch.Size([output_dims])\n", " variational_distribution = qpytorch.variational.NaturalVariationalDistribution(\n", " num_inducing_points=num_inducing,\n", " batch_shape=batch_shape,\n", " power=self.power\n", " )\n", " variational_strategy = qpytorch.variational.MultitaskVariationalStrategy(\n", " self,\n", " inducing_points,\n", " variational_distribution,\n", " learn_inducing_locations=True,\n", " # jitter_val = 1.0e-4\n", " )\n", " \n", " super().__init__(pde_lhs, variational_strategy)\n", " # self.mean_module = qpytorch.means.ConstantMeanGradGrad()\n", " self.mean_module = qpytorch.means.LinearMeanGradGrad(input_dims)\n", " # self.base_kernel = qpytorch.kernels.RBFKernelGradGrad(ard_num_dims=input_dims, interleaved=False)\n", " # self.base_kernel = qpytorch.kernels.RQKernelGradGrad(ard_num_dims=input_dims, interleaved=False)\n", " self.base_kernel = qpytorch.kernels.Matern52KernelGradGrad(ard_num_dims=input_dims, eps=1e-4, interleaved=False)\n", " self.covar_module = qpytorch.kernels.ScaleKernel(self.base_kernel)\n", " self.covar_module.base_kernel.register_constraint(\"raw_lengthscale\", qpytorch.constraints.Interval(1e-2, 1))\n", " \n", " def forward(self, x):\n", " N = x.shape[-2]\n", " mean_x = self.mean_module(x) # ... x N x (2D+1)\n", " covar_x = self.covar_module(x)# ... x N(2D+1) x N(2D+1)\n", " mean_x = mean_x[...,:-1]\n", " covar_x = covar_x[...,:-N,:-N]\n", " return qpytorch.distributions.MultitaskMultivariateQExponential(mean_x, covar_x, power=self.power, interleaved=False)\n", " \n", " def predict(self, testX):\n", " with torch.no_grad():\n", " output = self(testX)\n", " if type(likelihood) is qpytorch.likelihoods.QExponentialLikelihood:\n", " output = output.to_data_uncorrelated_dist()\n", " predictions = likelihood(output)\n", " pred_m = predictions.mean\n", " pred_v = predictions.variance\n", " if pred_m.ndim == 3: pred_m = pred_m.mean(0)\n", " if pred_v.ndim == 3: pred_v = pred_v.mean(0)\n", " return pred_m, pred_v\n", "\n", "# likelihood\n", "likelihood = qpytorch.likelihoods.QExponentialLikelihood(power=torch.tensor(POWER), noise_constraint=qpytorch.constraints.Interval(1e-2,1.0))\n", "likelihood.noise = torch.tensor(0.02)\n", "# define the model\n", "model = QEPsolver(pde_lhs=burgers.lhs_f, power=torch.tensor(POWER))\n", "model.covar_module.base_kernel.lengthscale = torch.tensor([.25,.5] if POWER>1.0 else [0.15, 0.2])\n", "# set device\n", "model = model.to(device)\n", "likelihood = likelihood.to(device)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Training\n", "\n", "We use the variational inference with natural gradient descent as in this [NGD example](../04_Variational_and_Approximate_QEPs/Natural_Gradient_Descent.ipynb)." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "b8c12647c15d455e8d6ef14f3c14d369", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/200 [00:00" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Set into eval mode\n", "model.eval()\n", "likelihood.eval()\n", "\n", "# prepare quantites to plot\n", "pred1_m, pred1_v = output.mean, output.variance\n", "N_pts = int(np.sqrt(burgers.Nd+burgers.Nb))-2\n", "\n", "# Initialize plots\n", "fig, axes = plt.subplots(1,2, figsize=(14,5))\n", "\n", "# solution\n", "u = burgers.extract_solution(pred1_m)[0][0].detach().cpu()\n", "axes.flat[0].tick_params(axis='both', which='major', labelsize=15)\n", "burgers.plot_solution(u.reshape((N_pts,N_pts+1)), ax=axes.flat[0])#, levels=50)\n", "induc_pts = model.variational_strategy.inducing_points.detach().cpu()\n", "axes.flat[0].autoscale(False)\n", "axes.flat[0].scatter(induc_pts[:,0],induc_pts[:,1], marker='x')\n", "axes.flat[0].set_title('Solution (q='+str(round(POWER,1))+')', fontsize=20)\n", "# uncertainty\n", "v = burgers.extract_solution(pred1_v)[0][0].detach().cpu()\n", "axes.flat[1].tick_params(axis='both', which='major', labelsize=15)\n", "burgers.plot_solution(v.reshape((N_pts,N_pts+1)), ax=axes.flat[1])#, levels=50)\n", "axes.flat[1].set_title('Uncertainty (q='+str(round(POWER,1))+')', fontsize=20)\n", "\n", "None" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One can try for longer time to obtain solution with improved precision. Refer to the [NIPS2025 paper](https://nips.cc/virtual/2025/loc/san-diego/poster/117572)." ] } ], "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.10.18" } }, "nbformat": 4, "nbformat_minor": 4 }