Skip to content

Instantly share code, notes, and snippets.

@bmorris3

bmorris3/piecewise_polynomial.ipynb

Created June 25, 2021 10:34
Show Gist options
  • Save bmorris3/198264ed8156299157cab0e9bce0c977 to your computer and use it in GitHub Desktop.
Save bmorris3/198264ed8156299157cab0e9bce0c977 to your computer and use it in GitHub Desktop.
Download ZIP
Kitzmann et al. (2020) T-P profile parameterization
Display the source blob
Display the rendered blob
Raw
piecewise_polynomial.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "unknown-fantasy",
"metadata": {},
"source": [
"# Kitzmann et al. (2020)'s finite element approach to the T-P profile\n",
"\n",
"* Paper ref: https://arxiv.org/abs/1910.01070\n",
"* Code ref: https://github.com/exoclime/Helios-r2/blob/master/helios_src/additional/piecewise_poly.cpp"
]
},
{
"cell_type": "code",
"execution_count": 561,
"id": "great-stable",
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import pymc3 as pm\n",
"\n",
"gl_0 = [0.0]\n",
"gl_1 = [-1.0, 1.0]\n",
"gl_2 = [-1.0, 0.0, 1.0]\n",
"gl_3 = [-1.0, -0.447214, 0.447214, 1.0]\n",
"gl_4 = [-1.0, -0.654654, 0.654654, 1.0]\n",
"gl_5 = [-1.0, -0.654654, 0.0, 0.654654, 1.0]\n",
"gl_6 = [-1.0, -7.65055 -0.285232, 0.285232, 0.765055, 1.0]\n",
"\n",
"quadrature_nodes = [gl_0, gl_1, gl_2, gl_3, gl_4, gl_5, gl_6]\n",
"\n",
"\n",
"class Element(object): \n",
" def __init__(self, edges, order):\n",
" self.reference_vertices = quadrature_nodes[order]\n",
"\n",
" self.nb_dof = len(self.reference_vertices)\n",
"\n",
" self.dof_values = np.zeros(self.nb_dof)\n",
" self.dof_vertices = np.zeros(self.nb_dof)\n",
"\n",
" for i in range(0, self.nb_dof):\n",
" self.dof_vertices[i] = self.referenceElementMap(self.reference_vertices[i], edges[0], edges[1])\n",
"\n",
" def lagrangeBase(self, r, i):\n",
" l = 1\n",
"\n",
" for j in range(0, self.nb_dof):\n",
" if (i != j):\n",
" l *= ((r - self.reference_vertices[j]) / \n",
" (self.reference_vertices[i] - self.reference_vertices[j]))\n",
" return l\n",
"\n",
" def getValue(self, x):\n",
" # coordinate on the reference element\n",
" r = self.realElementMap(x, self.dof_vertices[0], self.dof_vertices[-1])\n",
" \n",
" y = 0\n",
"\n",
" for i in range(0, self.nb_dof):\n",
" y += self.dof_values[i] * self.lagrangeBase(r, i) \n",
"\n",
" return y\n",
" \n",
" # maps the coordinate value r on the reference element [-1, +1] to the real element [x_l, x_r]\n",
" def referenceElementMap(self, r, x_l, x_r): \n",
" return x_l + (1.0 + r)/2.0 * (x_r - x_l)\n",
" \n",
" # maps the coordinate value x on the real element [x_l, x_r] to the reference element [-1, +1]\n",
" def realElementMap(self, x, x_l, x_r): \n",
" return 2.0 * (x - x_l) / (x_r - x_l) - 1.0\n",
" \n",
" \n",
"class PiecewisePolynomial(object):\n",
" def __init__(self, element_number, polynomial_order, domain_boundaries, dof_values):\n",
" self.nb_elements = 0\n",
" self.nb_edges = 0\n",
" self.elements = []\n",
" log_boundaries = [np.log10(domain_boundaries[0]), np.log10(domain_boundaries[1])]\n",
"\n",
" self.nb_elements = element_number\n",
"\n",
" self.nb_edges = self.nb_elements + 1\n",
" self.order = polynomial_order\n",
" if (polynomial_order < 1): order = 1\n",
" if (polynomial_order > 6): order = 6\n",
" self.createElementGrid(log_boundaries)\n",
" self.setDOFvalues(dof_values)\n",
"\n",
" \n",
" def createElementGrid(self, domain_boundaries): \n",
" domain_size = domain_boundaries[0] - domain_boundaries[1]\n",
" element_size = domain_size / self.nb_elements\n",
"\n",
" element_edges = np.zeros(self.nb_elements+1)\n",
"\n",
" element_edges[0] = domain_boundaries[0]\n",
" element_edges[-1] = domain_boundaries[1]\n",
"\n",
"\n",
" for i in range(1, self.nb_edges-1):\n",
" element_edges[i] = element_edges[i-1] - element_size\n",
"\n",
"\n",
" self.elements = []\n",
" for i in range(0, self.nb_elements):\n",
" edges = [element_edges[i], element_edges[i+1]]\n",
" self.elements.append(Element(edges, self.order))\n",
"\n",
" self.dof_vertices = []\n",
" for i in range(0, self.nb_elements):\n",
" for j in range(0, self.elements[i].nb_dof-1):\n",
" self.dof_vertices.append(self.elements[i].dof_vertices[j])\n",
"\n",
" self.dof_vertices.append(self.elements[-1].dof_vertices[-1])\n",
" self.nb_dof = len(self.dof_vertices)\n",
" \n",
" def setDOFvalues(self, values):\n",
" if len(values) != self.nb_dof:\n",
" raise ValueError(\"Passed vector length does not correspond to the number of dof!\\n\")\n",
" \n",
" self.dof_values = values\n",
"\n",
" # set the dof values in each element\n",
" self.global_dof_index = 0\n",
"\n",
" for i in range(0, self.nb_elements):\n",
" for j in range(0, self.elements[i].nb_dof):\n",
" self.elements[i].dof_values[j] = self.dof_values[self.global_dof_index]\n",
" self.global_dof_index += 1\n",
"\n",
" self.global_dof_index -= 1 # ; //elements share a common boundary\n",
" \n",
" def getValue(self, x):\n",
" # The validity check below hasn't yet been implemented, but is preserved in comments:\n",
" # {\n",
" # //check validity range\n",
" # if (x > dof_vertices.front() || x < dof_vertices.back())\n",
" # {\n",
" # std::cout << \"Requested x value outside of domain of the polynomial!\\n\";\n",
"\n",
" # return 0.0;\n",
" # }\n",
"\n",
" # first, we check if x is a global DOF\n",
" for i in range(0, self.nb_dof):\n",
" if (self.dof_vertices[i] == x):\n",
" return self.dof_values[i]\n",
"\n",
" element_index = 0\n",
" \n",
" # if not, find the element it is in\n",
" for i in range(0, len(self.elements)):\n",
" # Assumes pressure is decreasing from beginning of array towards the end:\n",
" if (self.elements[i].dof_vertices[0] > x) and (self.elements[i].dof_vertices[-1] < x): \n",
" element_index = i\n",
" break\n",
"\n",
" # get the value from the corresponding element\n",
" return self.elements[element_index].getValue(x)"
]
},
{
"cell_type": "markdown",
"id": "average-refrigerator",
"metadata": {},
"source": [
"Usage example: "
]
},
{
"cell_type": "code",
"execution_count": 563,
"id": "gorgeous-denver",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAB+aklEQVR4nOz9eXgd2X3fCX9OVd19v1jvxQ4Q3MlustmbWi2pxVZ3K1JLiu04UpJxFscavzN2JD2vPLHGieORH0dOpHks+bWTjMb2eIlj2WMtLcqtbkm0lpa61c0dXEES+46Lu+9L1Xn/uCAIEACJlQDR9XkePCQOqk796oKsb53fdoSUEhMTExMTk+VQttoAExMTE5PtjSkUJiYmJiZ3xRQKExMTE5O7YgqFiYmJicldMYXCxMTExOSuaFttwGZQW1sr29vbt9oMExMTkweKM2fOzEgp6+4c31FCIYR4EXhx165dnD59eqvNMTExMXmgEEIMLTW+o1xPUsoTUsqP+3y+rTbFxMTEZMewo4TCxMTExGTj2VFCIYR4UQjx5WQyudWmmJiYmOwYdpRQmK4nExMTk41nRwmFuaIwMTEx2Xi2fdaTEGIf8AmgFjgppfyvyx0rpTwBnDh27NgvrfV6E5Mv0d/3BQrFCSg70b5aovZHgoRPJfK+nyOZOEim4setJXjyPRq7f+4jc+dmz02TenUQPVHk5VaVL7ZBzOHAl0vzK26FFl8Hn3+1l/FEnrDfwa89v4ePHGlatY09PT2cPHmSZDKJz+fj+PHjHD58eK23bGJiYnJXNnVFIYT4EyHEtBDi0h3jLwgheoUQN4UQv363OaSUV6WUvwz8PHDsHtdb14piYvIlrl37DQrFcUCCJYv+D8sUjumUbEeYmHqcTCUIKGQqQb5/0sb1v/0GUBWJxNduoCeKfLtR4z/udhBzOkEIki4v/6lk5Xe/+0PGEnkkMJbI85mvXeQb58ZWZWNPTw8nTpzg1j0mk0lOnDhBT0/Pmu7ZxMTE5F5stuvpT4EX5g8IIVTgD4H3A/uBjwkh9gshDgkhvnXHV/3sOR8CfgycvNvF1huj6O/7AoaRXzinDdIf1unr/BCGalvws4q08cYPKgCkXh1Elg0A/nC3jYIqFh5rsTKzd+HqIV/W+fyrvauy8eTJk5TL5QVj5XKZkyfv+tGYmJiYrJlNdT1JKX8khGi/Y/gx4KaUsh9ACPEV4MNSys8BH1xmnm8C3xRC/B3wP5Y6RgjxceDjAK2trWuyt1CcWHJcD0LRFlzyZ5mKv3pMojg3NmUXSx/r9mIns2BsPJFf8tjlWG61ZMZlTExMNoutCGY3ASPzvh+dHVsSIcR7hBC/L4T4v4CXlztOSvllKeUxKeWxurpFFegrwm4LLTmuxsBWjC35M7eWqB7jv73aaCgsvRmUO5NaNBb2O1Zl43KrJTPTy8TEZLPYCqFY6nV72W32pJQ/kFL+Gynl/yyl/MNNtIvOrk+jKAsf3KIInpdUuvq/iaIXF/xME0WefE91UeZ9vh1hqX6c/+v1InZ94S1p5RK11xbGIxwWlV97fs+qbDx+/DgWi2XBmMVi4fjx46uax8TExGSlbIVQjAIt875vBsa3wI5FhBo/zN69v4PdFgYElF2oX7dgP61iLZ4j1PAmbi0GGLi1GM8cL85lPbmO1OP/mW5Uv433T1b436/nCeZyICW+bIp/ay3x6+97N01+BwJo8jv43M8cWnXW0+HDh3nxxRfnVhA+n48XX3zRzHoyMTHZNMRm75k9G6P4lpTy4Oz3GnAdOA6MAaeAfyKlvLxR1zx27Jg0mwKamJiYrA4hxBkp5aLs0s1Oj/0r4A1gjxBiVAjxi1LKCvArwKvAVeBvNlIkTExMTEw2ls3OevrYMuMvc5fAtImJiYnJ9mFHtfBYL/n8CInEbZeVNAxyp06R7h9l/EacUqGy+CQp6bn8dcbeeA09XQJgsljmzUSGkmEsODQdnWH0SrX2sFQxODMUYyK5uvTY1aLrOn19fRSLxXsfbGJiYrIEplDMo1iaJpcfvD1gGBT7B8hOJohP5pbMzapInWh0mL7e80xGbgIQLVe4nisg7kjwKuVzxMZHSc1MoymCwZkcfdPZTbwjKBQKpNNprl69Sjqd3tRrmZiY7ExMoZiHoRdRhHXReLlkIITAYlMX/UxTNJ7s/ihO1crp/pcYz4yT1w00IbAoC4UiGG7G6nAwefM6QkBbjZOxRI5SxVg070bhcrnYu3cvmqZx/fp1ZmZmNu1aJiYmOxNTKOYhpY5QLIvGy2WJZlMQytIV16rmobvxEbxGnjeuf4OJ3Ax2ZfFHKxSFxl17KGQyxCfGaK91oRswGs9t+L3Mx263s3fvXtxuN4ODg4yPb4tsZBMTkwcEUyjmIdERS3wklZKBZlm8mpjDkKjeJt4RehRfapLzUxfIlxdXYQP4Gxpx+HxM9t0g6NBw2zUGo5vrfgLQNI3u7m5qamoYHx9nYGAAw9i8lYyJicnOYUcJxXq7x94ZU0BUv9d1iWa5y0c1G7uwdjzN0+4OtGyMG7GrRPPRJQ8Pd++lUiwyPdhPe42TqVSRfElfk82rQVEUOjo6CIfDRKNRbt68ia5v/nVNTEwebHaUUKy3e6wQGpJ5D85Z95GuGyjq0m6n6omzf9q82JofZa9wYS9n+fHYj0mVFq8sXP4AvoZGIkMDhF0KUsJQbPNXFbcIh8O0tbWRTqe5fv36om60JiYmJvPZUUKxXoSwYBiled/PKoDOsvEJAGZ/JnUD6vejOQIcroBi6Lw2+hq58uIYRKh7DwjIjvQTdFkYnLl/QgFQV1dHV1cX+Xye3t5eSqXSvU8yMTF5W2IKxTwU1YZhFFnQ1kQRSGnc8kItidCqH6OsGKAoyNq9OIF34qCkl/jx2I8p6wvf2q12B/VtnSSnJ2m0lIhlyyTz9/fN3u/3093dTblcpre3l0KhcF+vb2Ji8mBgCsU8VNUJ0liweZFQVASSu8V9hbUa6Jal6kHC5kYG2gmkJnnS00GqlOKNiTcw5MJJ6to6sDocqNODgMHQfQhq34nH42H37t3ous7169fJ5ze3ANDExOTBwxSKeWiqG4BKZV5hmqpUhUJfXimU2foKY7ZyWwAy2AV2L42RmxytfYip3BQXIhcWnqeqhPfsQxbzeHMRBu6z++kWLpeL3bt3I6U0xcLExGQRplDMw2LxA1AuJ+bGhKqhYKBXlu+yKzQFxaZi5KpCoQrQhQJtT0ExTWcuye7Abm4mbtKX6Ftwrre2Hm99A870JOl0lkh6a1ptOJ1Odu/eDWCKhYmJyQJ2lFCsNz1WVR0oio1yOX57Tk1FEzqVe6SvKk4LRrYaY9CEoGJI8DRC7W6YusRhZzMhV4jz0+eZyS+sjm7avY+gy0Zl/OZ9qalYDofDsUAszJiFiYkJ7DChWG96LIDVWkexFLk9oKpoqkG5eA+hcFswChVkWceqKJRuBcSbj4FmRwy/weMNj+G0OHlj/A3yldtv7Ba7nebu3Xj1HH03BjCMzd0j5G7cEotbbiizmaCJicmOEoqNwGqrQ69kqFSqb/ZCs6ChIw15V7FQvdUeUZVkCasiKN6Kfms2aH4UshEssT7eEX4HFaPCmxNvLghu17S00RSuIzPax8jM2lZEG8UtsdB1nRs3bph1FiYmb3NMobgDuy0EQLE4AYCwaFiUauyhmFv+gam4LAiLgh4vYFcUioa8nWZb0wXeJhg7g09YONpwlEg+wpXolbnzhRAcOnoUTRpcOnNhmavcP5xOJ7t27aJUKnHjxg2zgtvE5G2MKRR3YLH4UFUX+cIoAMJiwSqqApHPLC8UQgi0gB09UcSBQCLJz3chtT4B0oCRN2nzttHubedq7CrTuem5Q5xeD217djM5NkZkGzTu83g8dHZ2ks/n6evrM3tDmZi8TdlRQrHeYPYtHI5WisUpdL2IsFhQjDJWu0Y+fffqZa3GgTQktllByc5/C7d7ofEwxAchOcqR+iN4LB7emnyLkn573gMH94HdxeWz5ymXtj4+4Pf7aW1tJZVKMTw8vNXmmJiYbAE7Sig2IpgN4HS2gTTI54cQFguyVMLptZJLlhZWbd+B4rGg2FXssWq2UPbO2ovGQ1XBGP4pmoTHQ49TrBQ5N31u7pB6n4Pgrv1EkjnGrm6PrcTr6uoIhULMzMwwMTGx1eaYmJjcZ3aUUGwUFksAiyVANteHsFpBgtOtolcMCtl7uJ/qXTgyFYyiTqpyh19fUaH1SSimYbKHgD3Avpp9DKeHGU2Pzh3W1VxHzhdmZmKC2PjYZt3mqmhqaiIYDDI2NkY8Hr/3CSYmJjsGUyiWweXqplJOUFYyADid1WZPmdjd3UGWOgcWTcWeKi8WCgBvGIKdMHUJCin2Bvfit/k5N31uzgXVXutCq20iozoZ771CKb+5GxutlPb2dlwuF4ODg2ZBnonJ24gdJRQbFaMAcDjaUBQrOTkCgCrLODxWUtG7PyCFpmBpcOLOVohllhGV5kdBqDDyJopQONZwjKJepCfSA4DXbqHWYyMTbAchGL7Ug9wGgWRFUejq6kJRFPr6+qhUKlttkomJyX1gRwnFRsUoABRFw+XaTZEZKiKHLBbx1jooZMoU83d/QFoaXQRUhdhMDn2pmIbVCeEjkByF+CABe4DuQDcDqYG5qu2OWhepioq3rZtcMsH0YP+672kjsFqtdHZ2UiwWGRwc3GpzTExM7gPibsHZB5Vjx47J06dPr3sewyjS8/I/J3v1DOUmA11zkhv+FSLT+8mnyriDNp78cBe7H29cdO6fnRvmpYEI1xxglYJfHTP46LFWXEfqb00OV78Jg69B7yuUU2P851Azr7hcJPQcDY5mGgr/kuvDQWyRfh4u9bOrNIqeTeKpqeXpj/4C+55+Zt33uFYmJyd58803uXHjBplMBp/Px/Hjxzl8+PCW2WRiYrI+hBBnpJTH7hzfUSuKjXQ9AQy+8h+JilMYXgkGUJujUv8dpFYtlMvEinz/L69x/c3JBed9dTLGZ2MxEkKyK2MwZRf8dofKX/2kn+y52boJRYFsBM7+GaRG+Y7LwdeskoRejUdM5Uc5k/pbonovjnKaltgV9EwCpCQ9E+E7X/4Drr72/Q25z7UwPT1NT08PmUw1hpNMJjlx4gQ9PT1bZpOJicnmsKOEYiNdTwDDhb8Cx+yW2NVO4tj9o3iabj8MKyWDN15a2BH2c/0TZBWYcChYDOjKGBRUwR90Wkm9Onj7wJ98CWY3NPq/fV5KysJfh2qbQnNf5eH0JTQWxigqpSKvfeXPN+Q+18LJkycXVWuXy2VOnjy5RRaZmJhsFjtKKDYa3Vd9EIoCs2oBFlcUzbEwPfTOTKixYvXhH7cpTDoU9icNXGXJlF2gJ+Ydm7ydEhtaIjAsLAmEmqGsWJa0Lx2dWXL8frDcqm2jVnMmJibbhx0lFBvtelKTs8sIKatiAQghUW2ZBce5g7YF3zfZbj/Yh5wKRQWOxSo05g1U/7xjfc1zf91TKtN4h1gIYSCUEtPWOgwW78Xqqaldy21tCMut2jZqNWdiYrJ92FFCsdGup1b7xxAlwCpQ49UHtaGrJAaenDtGsyo8+eGuBed9pjOEY/bBXlEF5wIqDQXJZ3qLeJ9vv33g8d8EiwOAd+XzPJovoMxLLpCGSjl1mJuuLmKzmyrdQlE13vGzH9uQ+1wLx48fx2JZuNJRVZX3vve9W2SRiYnJZrGjhGKjVxSd/+D/oN34ZyglBTUnUNIK9vSjqMY+AGxOjad+dteirKefbQzyhX0thIUKUlJW4D0pybOttdg654nY4Z+HF38ffC0cLZb4aEHypLQhEATtQR6tfY569QjT9gbOhp9B8wZBCNzBGg68+73Yvb4tq684fPgwL7744twKwuPxcPjwYWprt26VY2Jisjls+/RYIcRHgA8A9cAfSim/c69zNio99haF69fJnzuP98MvEk39EF3PURN8nv5zKewuC+2Hln84nkpmuZEt8HN1ASpXo8iygeNg7dw+2wu4+T1IT8KBnyFuFPne8Pc4XHuYfLaBqxMp/uGRJuyW6nnxyXFGLvVU97HYs2/D7nU99Pf3E4/HOXDgAHa7favNMTExWSVbkh4rhPgTIcS0EOLSHeMvCCF6hRA3hRC/frc5pJTfkFL+EvAvgH98j+tt6IriForTVbUlXyAQeBIpdVLpt6htdpNNFEnHlt8ytMVuRUcyUSlj3x0ACYXrcaS+hEA3PwqGDuNnCdgD1DvquZm4SVvQgZQwFL3dyiPQGKa2tY3oyNC26QfV3NyMoiiMjIxstSkmJiYbyGa7nv4UeGH+gBBCBf4QeD+wH/iYEGK/EOKQEOJbd3zVzzv1382etywbHaO4heJyAmBks1gsPny+oxSLk1g9g1gdGpP9yWW3L623atgUheF8CcWuYdvlx8iVKfYlFneitfug4QDM3IDsDN2BbnKVHFkjQtBlYWBm4X7aoV17cAWCjF27TC6Z2NB7XgtWq5VQKEQymTSzn0xMdhCbKhRSyh8BsTuGHwNuSin7pZQl4CvAh6WUF6WUH7zja1pU+U/At6WUZ5e7lhDi40KI00KI05FIZLnD1oTiqq4ojGz1Qe1ydeFwtpPOXCLYkqOUrxAdzSx9rhC02q2MFEtUDInmt2Fr81KJFygNpxef0Hi4GuAeeYuQK4RTc9Kf7Ke91kUsWyKZv929VigKbYceRrNaGew5R7m4/MrmflFfX4/NZmNsbOyuLdlNTEweHLYimN0EzPdNjM6OLcevAs8CPyeE+OXlDpJSfllKeUxKeayurm5jLJ1FsVoRFgtG5rYY+H2Poln8FI0zuIIVIiPpZXtAdThs6FIyUqh2h7U0urA0uihPZimN3SEwmhXCRyEzhYgP0uHrYCo3Rb1XIAQM3rGq0KxW2h9+BL1SZrDnHMYWb1mqKArhcJhcLme2Izcx2SFshVAsLgiYK2db4gdS/r6U8hEp5S9LKf/bJtp1VxS3a25FAdWmgTXBdyFQsPgvACXGr8eXfIuut2q4VJW+/O1iO2urB63WQWk0TXlq4cOf2m5wBGDsDG3uaq3FZH6ERp+dwWh20TUcbg+tBx4in0wyevUSW00wGMThcDA+Pm6uKkxMdgBbIRSjQMu875uBrd8g+h6objd6eqGrSNNcBINPI0QBZ/1lssk80bHsonOFEHQ6bEwUS2Rn96gQQmDr8KEF7BQHU5Snc/NPgJbHoJjGFR+i1l7LSHqEjhoX2aLOdHpx+3JffQMNXd0kJieY6r+5sTe/SoQQhEIhCoUCiURiS20xMTFZP1shFKeAbiFEhxDCCnwU+OYW2LEqFI8HI5tdVLdgs9URCDyBxZlC9VxiaiBJIbN4F7xdzmpF9o3c7Ye8UAS2XX40v43iQHKhWHjD4GuBiQu0OOpIlVJ4nWU0VSwKat+ioaOLQCjMVP9N4pNbq72BQACbzcbk5OS9DzYxMdnWbHZ67F8BbwB7hBCjQohflFJWgF8BXgWuAn8jpdwem0PfBcXtBskC99MtHI5WfL6juGri6OpFRq7FMO7YL9utqYRtVm7mChjz3DFCEdi6A7fFYnLe/M3HwKgQTld7Ok0XJmkNOhmO5ajcuR/3LE37DuAKBBi5fJFsYutiBEIIGhsbyWazpNNLBO1NTEweGDY76+ljUsqQlNIipWyWUv7x7PjLUsrdUsouKeXvbKYNG4Xq8QBgpFJL/tzt3oPffwh3fYRM9jzjNxenh+5x2ckbBsOzQe1bzIlF0E5xKHU7wO3wQ91enPFBAsLKRGaCjloXFV0yGl96pz1FUWk7fBSrw8HghbMUsktnY90PgsEgmqYxPT29ZTaYmJisnx3VwmMzUbxegEVxivl4vYcJ1u3HERwnMnma2PjC1UeTzYJHVbmSWfyQF4rA1uWfC3AXB5PVQHDoYVAsNGSiRAtRAk4Fl01lILq0+wlAs1joeLhaXDlw/gzl0t33+d4sVFWlpqaGRCJBubzYHWdiYvJgYArFClGsVhS7DX2ZFcUt/L5HaGg+iOocYvjm62ST82ISQrDP7SBarjBVXPzgFIrA3uXHEnJRnspRvJFAKjZoPER9IYMsJIkWorTVuJhMFsiXlk+FtTmddDz8CJVSkcHzZ9C3aH/r2tpapJREo9Etub6Jicn6MYViFShe77Kup/n4/cdo6jyM1Pq5eemHFOcVyXU5bNgUhUtLrCpuYWv1zhXlFa5GMfx7qHHWI+JDRAtROmpdSAmDd1lVADh9ftoOPUw+nWLo4rktaSDocDhwuVzEYnfWXZqYmDwomEKxClSfDz15b6EQQlBT8xjNux6hrPdzo+cklXL1jV5TBPtcdsaLJaKl5d/yLY0u7LsDGPkKhWtJROAIfr1CNHIFn8NC0GVdVHy3FN7aepr2HiATjTJy5dKW1DUEg0FyuRz5/PLiaGJisn0xhWIVqD4fslJBz9z7AQ1QV3+M1l2PUygMcv3CK1Qq1ZXFHpcdixBczOTuer4WsGPfXwMCClNeamQziakeMAw661zEc2USudJd5wCoaWqZrbEYZ+LGtRXZvpH4/X4As6bCxOQBZUcJxWZ1j72FOtts0FhFA7768BGa2t9BNjvCjUvfolIpYFUU9rsdjBRKzNxlVQGguiw4DtSguKz4S4exRqxkJy7QGnSiCOhfwaoCqjUWNS1tzAwPMTXQd+8TNhCbzYbL5TIbBZqYPKDsKKHYrO6xt7glFJVVvhmH2g4TbnkPmeQ0Ny9/k0olzV6XHZuicD5991UFgLCo2PcGcXQ04yw2kjpzFWulQtjvYCiaXbZz7Z2Ed+/F3xhmqu8G0dHhVd3DevF6vWSzWSpbFFQ3MTFZOztKKDYbYbFUez6t4c24qXMfjeFnySRS3LzyTWQ5ykG3g4liifHCvd1HQhF4u+pJtTsoZsvk3+ihw2YhXzKYTK2sa6wQgpb9B/HW1jF27Qrxifu3j4XX60VKaRbfmZg8gJhCsUpUv5/KGruitnR30RB6nkxcZ+D639EiJnCrKmdTuRUFmR2ag1K9k0yHRKRH8Y+lCaTKDERWXlQnFIXWww/jCgQZuXKR5PT9abHhcrlQFIVMZusKAE1MTNaGKRSrRAsEMDJZZOneq4ClaO5upqHxebIxF2P9P2SXGCBWLnMzd++iOCEEds1OPtSKI5TAoozTUpbkrkQpZFduj6KotD90FKfXz9DFC6RmNr9yWlEUnE4n2SVaoJiYmGxvTKFYJWogAEAlnljT+UIImnfX09D4XvKxENrkBRzFa5xNJimtoM7BptooqBZEXRd2ay/1e90oRYOJtyYoz6w8/VTVNDqOPILD7WGo5zzp2Mya7mc1OJ1OcrmVrZ5MTEy2D6ZQrBI1GARAX0fDPSEETbuDNIQfp5TYR3NiiEjyEqdjE/c816JYKBvlamsPKalXByh2eZksVSj2JSjciCPLKyusUzULHUeOYXM6Gbxwlkxsc6unnU4nhmFQLG5NSxETE5O1saOEYrPTYwEUmw3F6URfZ0sKIQTh7gD1TftwZ45Qn01zZqqHkcTd6xwsioWKUQG7F2p3w8x12oMw4NOo1DvQ40VyFyNUYisLcGtWKx1HH8VqdzBw4QyZ+OZVUNvtdgAKha3fstXExGTl7Cih2Oz02FtoNUEqsY1p4R3q8tHQ1sKu7GGMtI0fTtxgZuY1DGPpmIMQAoPZFUPoIRCCjvJNEIJRq8BxsAbFqlK4Eadwc2WrC4vVRuctsTh/etPak9ts1T05zBWFicmDxY4SivuFWlODkc1ibNCbcX2bl7buevaVdjERr+diIsL09CsUi5FFxwohbm8ca3VC3V4c6UHC9hL9M1kUpwX7/hqszR70WJFcT4TKCmIXFpt9Tiz6z53alJWFxWJBURRKa0wEMDEx2RpMoVgD2q04xQZ2RK0Ju3l8Xx11ej3npnaTKMLMzElS6UtIeXtVYEijKha3aDwEQmW3cZNMocJ0uoBQBNYmd3V1YVcp9CXI98YwincvdpsvFgPnT29KzMJisZgtx01MHjDETsxAOXbsmDx9+vSmzZ/4xjeY+tzvYiSTaOEw9Z/6JL4XX9yQuadjef7y8gTT+RwDxjmEdYJK0UtXROW94f+L3mKSeEUQFTU8d/Df8YHOD8DoGSqn/x++emGGNn2QZiL8kfWf8fAHPs6HHw5TmcpRGqkWulma3JQms6S/M4SeKKL6bXifb8d1pH7OhnKxQP/Z05QKOSqFAqe+9XXS0Rk8NbU8/dFfYN/Tz6zp3np6enjllVfI5XL4fD6OHz/O4cOHN+RzMzExWT9CiDNSymN3jmtbYcxmIYR4EXhx165dm3aN5IkTTP7W/4GcdTtVxseZ+Pe/CbAhYlEfdGBpsPPdniiu3G70YJCws5eR5jLn6MLOOXwaPGqb5mvnPgPAB5IJxIW/4mAlzBXZzjG1l/+t/F/4za9XgP+FjxxpQg3YKQ2lyLw+Tu7cNOjVFwQ9USTxtRsAc2JhsdnpfORRXvvLP+XCd7+NoVdXIumZCN/58h8ArFosenp6OHHixNxqIplMcuLECQBTLExMtjk7yvV0P4LZ07/3xTmRmLtuocD0731xw67xpzNxxv0aCZeKbaaWntLjJPBzg704hUFGr4YpnvPk+NLZL8EPPocqK7Qq05RRuSZbcIoSn+QrfP7VXgAUm4p9d4DC9ficSMzZXzZIvTq4YMxitXHj1BtzInGLSqnIa1/581Xf08mTJxe5nMrlMidPnlz1XCYmJveXHSUU94PKxNK1DsuNr4WxYhlDFUR8KgmXQmhS46J8mCg1NFslIOkrKngUyWR2EpKjAHhEgVYxzRWjnay0ERZRxhMLA9lGaulAsp5YnIm0XIwiHV19cd5yKctmR1kTk+3PjhKK+1FHoYVCqxpfC002CwCGKoh6VHJ2gSOnkcKHVYEOm0FSF7ye1ahz1IGvee7ch5Q+BJIfGg8xbNQR9jsWzK36bUteU/FaF415amqXPHa58bux3Cpvs1OZTUxM1s+OEor74Xqq/9QnwbrwoSrs9ur4BvGZzhC3HucFq0BXBI5KhffxMgANFthj07mYt/JU01NceOJfUVSqNnlEgXeql4hIH78hP87/97ndC+b2Pt+OsNzxa1cF9t0BypGFLc+f/ugvoFkXCotqsfL0R39h1fd0/PhxLBbLgjFN0zh+/Piq5zIxMbm/7Khg9v3A9+KLFK5fJ/H//m016ykU2tCsJ4Cfbaym33722ihTUsdeMmjKJpFuFxVFIaUb/LhYz7989DM0u5u5nuxj6ujP8nDPt6ktJdGlwjnrIzTtf45Grx0p5VxK7a2AderVwbmsJ8+zraguC8X+JLKoY232ALcD1q995c9Jz0Swuz0cePfxNWU93QpYnzx5kmQyicPh4LHHHjMD2SYmDwBmeuwqkVKS/Po3sLY043z00U25xp1kk0VevzDFtcA0Ids4P9v94oJaionMBKcu/gWVQoKDD/8Luv3dCCG4PJ7kwkiSrjoXj3UEF9Zf3IE0JKXBJOVIHq3Gga3Th1AWHh8dG2Hs6mWa9x8kGG5eZqaVMTAwQDwe59ChQ4tWGiYmJlvDcumxO8r1dD9iFEY6jSyXUWtX76dfK06vlTaHla5CjomKjTOphS6ikDvEc7UP02D1cSFygdfGXiNXznEg7ONgk5e+SJY3B2J37doqFIGt04+12UMlmqfQG0NWFrb/CIabcfoDTNzopVxaXxuOxsZGDMMgEllcfW5iYrK92FFCcT9iFPrsNqiq379p17gTIQSBBhehXI4OtcjVbJ7e7MIUXbtm5ylvF480PEI0H+W7Q99lJD3C4WY/h5p89EeyvN4Xvee2qdYmN7YuP0a6TP5qFKOkL7Cjed8B9EqF8WtX13VPDocDn8/H9PQ0xgraq5uYmGwdO0oo7suKIl9NN1Vdrk27xlL4G52ApKuQpslm5XQyy2RxXl2CaoNKkU5fJ8+2PYvb6uanEz/lrYm32Bty8nCLn6Fojh9ej1DW7/5gttQ6sO0JIIs6hStRjMLtWgq7y01D5y6S05Pr3h2voaGBSqVCdANboZiYmGw8O0oo7seKQhYKCFVBWBenk24mFquK028llyrxpNeJR1P5UTxNpjL7xm9xgF4CvYzH6uGZlmfYX7Of4fQw3xn6DvWBEo93BplMFfjulSmy9+j7pPls2PfVIHVJ/nIUPXtblOrbOnB4vIxdu0JlHQ3+vF4vTqeT6enN32HPxMRk7ewoobgfKwqpG6Comzb/3fDV2pEG5KYLvCfowZCSH8TSVAxZ3Z8CoJACQBEKB2oO8EzLMyhC4fsj3yfHAO/aXUO2WOE7VyaZydw9zqC6LDgO1CBUQeFqDD1dFQWhKDQfOESlXGb02uV13VNdXR35fJ50Or2ueUxMTDaPbS8UQohOIcQfCyH+9l7H3pf9KAQgt8anbnNqWO0a0bEMHkXhnQEP8UqFU6ksOKpbtJJb6MapcdTwbNuzdHg7uBa/xpXk6zyxy4EiBN+7MkVfJHPXayp2Dfu+GoRFoXAthp6siovD7aGxq5vU9BTxibE131MwGETTNHNVYWKyjdlUoRBC/IkQYloIcemO8ReEEL1CiJtCiF+/2xxSyn4p5S9upp2rQVityIqO3JIArMRT66Bc1ElM52i2WznodnIzV6DPsIFmh8zUorMsioVjjcd4KvwU+UqeU5Ef0dWUoc5j5c3+GG8NxNDvEuRWbCqO/UGEXaVwPU4lUQ2k17V14PQHGOu9SimfW/b8u6GqKjU1NSQSCXOfChOTbcpmryj+FHhh/oAQQgX+EHg/sB/4mBBivxDikBDiW3d81S+ecnnuy1aoDicARm5tD8b14nBZsbstREbSSEPysMdBg9XCm8ksSVcIUuOwTBps2B3mfW3vo95Zz9V4D5r7Op31Fm5OZ/jO5UlSheX3iRAWFce+GhSHRvF6gkq8gBCC1gOHABi+fHHN4llXV4eU0gxqm5hsUzZVKKSUPwLu3CrtMeDm7EqhBHwF+LCU8qKU8oN3fK3YHyGE+DjwH4Cjm/lmqnqrVcvGljSzE4Ckvs1LuaATn8ohhOCpgBtNCF5TGzFKOcgu37TPoTl4Z9M7OVp/lFghypT+FnuaS2RLOq9cmqT/Lq4ooSnY9wZRnBrFG1WxsDqcNO3dTy4RZ3qwf013Zbfb8Xg8zMzM3LXWw8TEZGtYsVAIIRxCiD0bcM0mYGTe96OzY8tdt0YI8d+AI0KIzyx3nJTyy1LKY1LKY3V1dRtg5tKofj8ogsqWvP0KJBJP0I7TayUynEbXDVyqyhM+F3FLgPPSDdGb95ypy9/F8bbjuC1u+rPnqK8bxWsX/LQ/xo9vzFCs6EueNycWrlmxSBQINIbxN4aZGri55i1Ua2trKRaLZDJ3j5mYmJjcf1YkFLMbAp0HXpn9/mEhxDfXeM2l+kgs+xoppYxKKX9ZStklpfzcGq+5YQhNQwsGKU8tjgVs+rWFMhdIb+jwUSnpREerD9ZWh40ut4srtiYiM0Og33u7Ua/VyzOtz7A/uJ/pwii68wJt9RVG4zlevjjBWGLpvbaFpmDfc3tloSeLNO3dj9XhZPjShTVVbfv9flRVZWZm9S3MTUxMNpeVrih+i6rLKAEgpTwPtK/xmqNAy7zvm4HxNc61JVjCYfRYHCObvb8XFgJJVSicXiveOgczoxlKswVxx3xOnP4wr5ds6JHrK5pSEQoHag/wnpb3oAiFifJZWkIzWFTBD3sjvN639Ori1sriVoCbvEHb4SPo5TIjl3pW7UJSVZVAIEAikUDXl17NmJiYbA0rFYqKlHKjnPKngG4hRIcQwgp8FFjr6mRLsLRUda40PHxfrytQFgSqG9qrtRNT/dXaCaui8ERDmJTNT894Hxgrf+DWOmp5tu1Z2rxtjOf7sHqv0FEvGIrm+LueCYaii0VRaAqOvUGEVaXQG8em2Gnau59MLMpU341V319NTQ26rpOYbZNiYmKyPVipUFwSQvwTQBVCdAsh/n/A6/c6SQjxV8AbwB4hxKgQ4hellBXgV4BXgavA30gp11e1dZ9R3W60+jqKN/vub5qsUOZWFABWu0Zdi4dUNE86Vk1ZDdutdDR0cKUAiclrq5reolh4tPFR3hF+BwU9z0TlFHuaszitGj+5GeX7vdNk7qjoFha1urJQBfneOP5gI4FwE9OD/aQiq6uNcLvdWK1WYrG1xTlMTEw2h5UKxa8CB4Ai8D+AJPDJe50kpfyYlDIkpbRIKZullH88O/6ylHL3bNzhd9Zo+5Zi696NkctRvo+rCiFUkMYCt05NsxurQ2OiL4E+28PpWGMrmsPHW6N9UC4sN92yNLmbeF/b+6hz1NGXuYTTf52DzQ4i6SJ/1zPOxdEklXn9ohSbin1PAAxJoTdOuHMvDo+X4cs9FHMrd88JIQgEAqRSKSqVu7cYMTExuX/cUyhm6x6+KaX8DSnlo7Nf/05Kufon0A7C0hRG9fvJX75831YVQlT3maouyqooiiDc7adc0IkMVdtg2FWFh1v3MqULBgfPrulaDs3B081Pc6TuCDP5GQYLb3C0U9IccHJxLMnLlyYZjd+uJVGcFuy7A8iSTqkvSevBhxFCMHjhLHrl3oH1WwQCAaSUpvvJxGQbcU+hkFLqQE4IYW5uPA8hBI5DBzEyWYo3Vu+PX9s1qz2mqr+S27h8NgIhF9GxDNnZFhu7g/UE/GHOzkyjJ9feYmNXYBfPtj2LU3NybuZNrO4B3rU7iCLgR9dn+PtrUyRy1boV1WPF1uVHz5SR40VaDz5EMZddVXD7lvvJFAoTk+3DSl1PBeDibM+l37/1tZmGPQhYwmEs4RCFS5fvS6W2MreiWPyG3tDhxWJXGb9edUEJITjWtp+s5uJa32morKPLq9XLe1vfy77gPoZSQ/TEf8RjXRYeaQsQy5b59qVJ3hqIkS/paEE7tjYvlVgBa9ZKePc+UjMRJvtWloUF1VTZVCplZj+ZmGwTVioUfwf8e+BHwJl5X297HEeOAJLcqVObfq3brqfFD1BVVWjaHaBUqDDVX01Qa3Q6aArv41JRUBp4bV3XVoTCwdqDvKflPQD8aOyHlLUhPnC4gd0NbvojGU5cqMYvqHNgaXBSnsjis9cRbGomMjhAbHxlKxu/349hGKRSqXXZbGJisjFoKzlISvlnm23Ig4rqduN46CFyZ85SuH4d++7dm3YtoVT3ljaMpX3+Lp+N2mYPM6Np3AE73loHD9c38nepdq7ErvDw5CVoPLguG2odtbyv7X2cj5znauwqU9kpHg09yu6GEBdGklwcS3JjOs2BkJdmj5XSQIrGvbso5nKMXr2E1eHAHQje9RputxtN00gkEgQCgXXZa2Jisn5WWpk9IITov/Nrs41bLfejKeBS2HbtwhIOk79wgcomVhbfzfV0i/o2Dw63lbHrCUqFCkGLRltDB9fszRRHz0BydN12WNRqGu2ToSfJlDOcHDrJdGGYd3bX8tyBBrx2C2eGE5zM55kulincSNK69zA2p5OhnnP3zIRSFAWv10sqlTJ7P5mYbANW6no6Bjw6+/U08PvAf98so9bKfdmPYhlcjz+G4nSSff31TavYrtYngmEsH28QiqB5X/UtfPRqHMOQHPI4KNfsptdSD/0/hNzG1Ck0e5p5rv05ahw1nJ0+y2ujr+GyGzy7v4Fn9tZhtWmckhUuDMUZuxil5eARAAbOnb5nmw+v10u5XCa3RV16TUxMbrMioZjtt3Tra0xK+UXgvZtr2oOFsFpxP/UUsqKTee3HyE3oYKso9xYKqBbiNe32k8+UmOxPErBoNDnsXAs+REWxwI3vQnFjmu85NAfvan4XR+qOEMlH+O7QdxnLjBHyOXjhYCNPHmyg0OigfzDBG29No7Tuo1QoMnjhLMZdgtVeb7Xq3IxTmJhsPSt1PR2d93VMCPHLgGeTbXvgUP1+XE+9Az2dIvPjnyA3uGhsTijkvUXIW+ugptlNfCJLfDLLfreDomploPndYFTgxqtQ2ri39flptK+Pv87pydOUjTItQSfvfbKV7v11uJIlrgwWuKE1MjQ6zcCFc8vWoFitVhwOh7lFqonJNmClrqf/c97X54BHgJ/fLKMeZCwNDbgee4xKJEL2Jz9BbmCKpxAKQrFgGCvrztrQ5sXltzFxM4m3KPFrGr26DbrfB+U8XH+l+ucGcSuNdm9gLwOpAb439D1m8jMIIWg+WMfhXTU8plpw+2oZszXykws3eOP1UwuqvOfj8XjIZDIYW7KboImJyS3ETgwWHjt2TJ4+fXpLbSj295M7dZrS0CDJb7+CPjmJFgpR/6lP4nvxxTXNOTH5Ejeu/zblShy7LUxn16cJNX74rudUyjr95yNM3Ejwp8kEP+iwkrcKDpTj/CfbdY7V1kMpCz/8T9VAt68Zjv8mHF7fe8BMfoa3Jt4iV8mxL7iP/mQ/f3Tqy+yaCrPX6ORw6jF6ExMk5QyiMUB07Cz2mQHsbi+vBx7nLaWdx11RDihj6OUiPp+P48ePc/jw4XXZZWJisjxCiDNSymN3jq/U9fQJIYRXVPkjIcRZIcRzG2/mzsHW2YmRSRP7879An5gAKamMjzPx73+T5IkTq55vYvIlrl37DcqVOACF4jjXrv0GE5Mv3fU8zaJSzuu8enWaIYegYBEgBJetQX652M3p3jfh7z4FyRFAVv888W+g52/Wcttz3EqjbfO28de9f81/eP0/0F8cxFNx8OzUozTlNd5ra+EwNTiHIuSKfq47dzFctHJg5Ac8UbrErsoAerm6ekomk5w4cYKenp512WViYrJ6Vup6+ldSyhTwHFAP/EvgdzfNqjWyVemxyxH7y/8Bd8QpZKHA9O99cdVz9fd9AcNY6CYyjDz9fV+457mnXxnk1cMObBVJKF6Za1U+aq3h122PLt7kqJyHk59dtY13ciuN9sdjP6Y8W/txNLsfC5a5YzptLTxkreOdwkFNaYa4JcCAqx2XXVBCXWhWuczJkyfXbZeJicnqWKlQ3NqV7h8A/4+U8sK8sW3DVqbHLkVlYmJV43ejUFz6nOXG55OJFRmr0RiptRBI6wQyt+Mml52dS5+0AfUWt5jJ364tmbDOcNZ5hRJVARVCUG9vo8FaT0Npht2Z69QVZygoDkb0AJO6m7yhzW3DsV1eAkxM3k6sVCjOCCG+Q1UoXhVCeAAzwngPtFBoyXF1DXt6221Lz7Xc+HzcQRu+nMGMTyXm0fDkDAKpqlg0lZfZ+9vXvGobl6PR1Xj776UaIpY4p9wXSanVFF0hBB5LtVrbgk5DKUKHMUmDkqaCQlLaGTQCRA0nLq9/w+wyMTFZGSsVil8Efh14VEqZAyxU3U8md6H+U59E2O0LBy0W3E+/k3xPz6rak3d2fRpFcSwYUxQHnV2fvue5T364i+NXClgqkokajaJVwVMwaEzpfCZQAcvCeVEtcPgfg74x6b2fOPoJ7Gr1c/ie/6c8ktmPFHDNPsioZYqCUeBS4nYvqgoKwxkrblGiSUkRUPLYKRMXXkT7Y7x+c4ZIevX7cpuYmKyNFfV6Ap4Ezksps0KIfwYcBb60eWbtDG5lN03/3hepTEyghULUfeLfYA2HKVy9RmUmiuvJJ1AcjnvMxFx2U9/N36VYmsZmbaBr17+9Z9YTwO7HG/kkYHlzkFc6LUz5FFrTkn+tuXi2493w4u9XYxLJUfA2wSP/HFx10PsydL0XbO71fAx8oPMDAHzp7Jf4IWfw2/z888iHGK6MEXdmmGqKM1acgJxAdfu5YW9nRndSkgEOOOOIQpZdXsk7n3kCa10r/ZEsg9EcAaeF7gY37TUuNHWl7zwmJiarZUXpsUKIHuAh4DDwF8AfAz8jpXz35pq3OoQQLwIv7tq165du3Kc9ItZKcWCA/NmzCE3D+dhjWJZxU91JqRQjEnmVYPCdOBwta7p22ZB8ZXyGwFSJzny15Ye35g6xSgzDwGsgBLQ/Df61XeteTA+M0nflCpFQhq7WPeyr2YdR0Rk4e4pCNkPb4SOMRWZQFIU9e/ZU7dcNhqJZrk9lSOTKWFRBe62L7no3fqd1U+w0MXk7sFx67EqF4qyU8qgQ4jeBMSnlH98a2wxj18t2qKNYCXoiQfanb6Ink9i6u3E8dBihqnc/Ry8yOfk1fL6juN171nztV2eSlCsG+ycrFDJlWvYH8QTvcJMVktD/g2pvqIaD0PQIKBv75i4NSebCJIOpIa41jhC0B3k09ChO4WDg7CnymTRabQMlCQ8//PCi86fTBW5OZxiJ5dANqHFb6a530xp0mqsME5NVsq46CiAthPgM8D8Bfze7ParlHueY3APV78fz7HFs3d0Ub9wg/Z3vUIndvWGfqtoQQkPX19d+o8FqIWEYhPcFsbk0Rq7GyMTv2N3W7oO9H4T6fTB1Ca59C/KJdV33ToQicHQE6LC38YR2dK4b7VB2mI6jj+Jwe4gODZBJxCkt0T+r3mPnHV21fPjhJo62+SnrBj/tj/H1c2OcGowRzZixDBOT9bJSofjHQJFqPcUk0AR8ftOsehshNA3n0SO43/0uZKVC+nvfI3/x4l1bf6iqk4q+vg61tVYNiSSBQdvBGmwOjeErS4iFokLrE9VYRSkLV0/A9LV1XftONL8d1WcjmHDzvuZn57rR/nT6TUKHD+L1+YiOjhAZG1l2DrtFZW+jlw8eDvPsvnqaAg76IxlevTzFyxcn6J1MUyibO+aZmKyFFbfwEEK0Ad1Syu8JIZyAKqXcVh3bHqQYxVLIUonc+fOUBgZRvR4cjzyCpb5+0XEzM9/HMErU1z+/5msVdIP/dyrGI14X+90OKmWdoYtRirnK0m4oqDYRHPwxpMbAG4a2p9Yd6L6Fni2TvzSDNezG0uzmZuImF2cuoikaB7z7uf6jc3gsKgcee4JAY3hFc5Yq1VhGXyRDLFtGEdAccNJR5yLktaMo264UyMRkS1lvjOKXgI8DQSlllxCiG/hvUsrjG2/q+nlQYhTLUZ6cJHf6DEY2i62zA/vhwyg229zPE4nT5PNDhEI/u67rfHUyRoPNwjsD1UbAetlg6FKUQrZM054AvrplsrEivTA6u/Vr86NQu7sa9F4nhRtx9EQR58N1CItKspjkrcm3iBfiqGMKjboTj1Bo2neAmqbVBdfj2RL9M1kGZ7IUKwYOq0J7jYuOWpcZADcxmWW9MYr/FXgKSAFIKW9QbeVhsglYGhvxvvA89r17KA4MkPr2tyn298/t9qZqbgyjtOIussvhs2gkK7fdMapFoe1QDXa3hdFrMeKTy7i36vbA/g+DsxaGXofeb29I7MLa7EEakvJkNf7is/k43nqcfcF9xI0EQ7YE0mNl7OplIkMDq5o74LLySFuAf3ikiae7awm6bPROpnn54iSvXJrg2mTKdE2ZmCzDSusoilLKkph9axRCaMC2azs7z/W01aasG6FpOB56CGtbG7mz58idOk2xrw/nkSNorqq7p1LJYLXa7jHT8vg0lZu5hX2eVE2h/VANw1dijN9IYOiSmqYl3Es2D+x5AWZuVFcXV16C0GFoPFyNa6wBxaGh1dgpT2WxhFwITUERCofqDlEKlbgyfYXrvij1ZQvG9avo5TKNu1a3R7miCFqCTlqCTgplnaFojoGZLGeHEpwbTtDos9NR46I54DCzpkxMZlmp6+k/AwngF4BfBf4X4IqU8jc21bo18qC7npaiNDRE/sIFjHwBpTVIunGMYMO7cTrb1zzn1Uye06ks/6ghiP2Oh6JhSMauxUlF89Q2e2jo8C4/UTkPI29CbKAqIK1Pgq9pTTYZuTK5izNYmz1Y5wnU8PAw05FpZJNkIDmAbTxPqOgj3N5N0579iHW6vpL5MgMzWYaiWbJFHU0RNAcctNWa8QyTtw/LuZ5WuqL4t8C/Bi4C/zPwMvBHG2eeyb2wtrVhCYcpXOul0HuZ/MwlsnEvjoNhhHVtPnbnrDjkDGORUCize29P3FSYGU1TKemEu/2IpR6YFgd0vgdqumHkp3DjOxBoq8YvbKvbCFFxWlB9tturitnrWa1WkHCk9giNrkbOqmfpG5kifT2NXirTcvAQyhpXMgA+h4WHW/w81OwjkikyOJNjOJZjMJrDpim01jhpCzqp89jWLUomJg8a9xQKIYQC9EgpDwL/9+abtHZ2kutpKYTFguPQQWydHaTOjZIfvEZyWGDftx/brq57FuvdiWO2eK6gG0tWxQghCHf70awKkeE05ZJOy/4g6nIuGV8TeD5SrbmY6Km2BGk4WHVHqSt9JwFLyEXhWoxKrIClthpQt1iqBpbLZVo8LdQ6ajllP8VI/1WS139MrpBh99HHUbX1lfcIIaj32Kn32DnWFmA8ma+6pyJZbkxlcFpVWoJO2mqc1LrX7vYzMXmQWKnr6S+Bz0gphzffpPWzE11PdzIz833K6RjuoRoqU9MoTif2AwewtrchVlg9nSxX+GYkwTv9Hjqcd3/oxSezjN9IYHdZaD1Qg8V2D1EqZmDsDMT6weqC8FGo6VpRdpSUkvzFGYQqcByoBSCVSnH9+nV2796N1+udO+5m4iZnr/2Y0lCE7sZ9PPLk+7DYlkjtXSdl3WA8URWN8UQeQ4LLptJW46I16CToMjOnTB581ut6CgGXhRBvAXOpMFLKD22QfRvCTl9RzMdi8VGyzeB+97upTE2Rv3iR3KlTFK9dxX7wIJaWlnu6SCyzglJewctCoNGFZlUZvRaj/3yEtgPV7Khlsbmh891QtxdG34LB12D6StUd5b17XyshBJZ6J8WhFEaujOK0zK0oKvM2ghJC0B3opv5oPT9x/j1Xr14k+t0I73rXP8TjDdzznlaDRVVoq3HRVuOiVDEYiVddU9cmUlwZT+G2a7QGnaZomOxIVrqiWLL5n5Tyhxtu0eJru4D/ApSAH0gp//Je57wdVhTZbB+JxFs0NLyIplWDvqXRMQqXLqEnk6g+L/b9++8qGCXD4K8nbxfdrYRCtszQpSh6xaB5TwBv7QrOk7K6shg7C6UM+Fqg6Sg4g8ufUjbInZtGa3Bia/NSLpe5cOECLS0tNDQ0LDrekAbnBt/i4qkfYBNWnnrnh2gNbf4LQ7GiMxrPMxzNMZkqICW47RotAQctQdM9ZfJgsaaCOyGEHfhlYBfVQPYfSynXvUmBEOJPgA8C07Oxj1vjL1BtX64CfySl/F0hxP8EJKSUJ4QQfy2l/Md3mfeBrsxeDcXSDDOR71JT827s9tuVylJKyiMjFK5cQU+mUL0e7Pv2YWltXeSSqhiSv5qMctTr4sAKhQKgXNIZuRIjny5R3+6lrmWFAWu9ApGr1fiFXqq6okIPg33pjKrCjTh6qoTzSD0SydmzZwmHw4TDy1dmT8RG+MGPvkYhl+XAw09xbN/TKOL+pLkWylXRGInnmEoWMCSzMQ0HLQEzEG6y/VmrUPw1UAZeA94PDEkpP7EBxrwLyAB/fksoZhsNXgfeB4wCp4CPAR8Gvi2lPC+E+B9Syn9yr/nfDisKwygxMfFVvN6H8Hj2L/q5lJLy6CiFy1fQk0kUlwv73j1YOzrmgt63hOKIx8lBj3N119cNxm4kSEXy+OochLv9KCutO6gUYfISTF+urjZquyH0UDWWMf+wWIHCjTj2vUE0n41z585RW1tLS8vdq7LzhRw//PHXmJgaora9g3c/+kG81ruk924CxYrOWDzPSDzPZDKPboBNU2ieXWk0eO2oZsqtyTZjrTGK/VLKQ7MT/DHw1kYYI6X8kRCi/Y7hx4CbUsr+2et9hapIjALNwHnuUkkuhPg41TYjtLa2boSZ2xpFsaKqLsrlxJI/F0JgbWnB0txMeWyc4rWr5M6cpXDlCrbdu7F1dmJo1V+/uoa3XEVVaNkbJOJMMz2Uopiv0LIviNW+grCXZoPmR6pdaScuwMx1iN6E2j3QeAisVdFS/TaEKtCjBTSfDUVR0O/SLPEWDruT5977Mc6d+QGXb5ziW5m/4LFHn2NXsHvV97lWbJpKZ52bzjo3pYrBRDLPaDzPUCxHXySLpgqa/NWVRqPPjlUzi/tMti/3+l89V7Yrpaxs8rK5CZjfHnQUeBz4feAPhBAfAE4sd7KU8svAl6G6othEO7cNFoufcjl+12OEEFibm7A2N1GemqZw9Qr5Cz0ULl/B6OxE+utRxNob+9W1erC5NMZ64/Sfj9CyN4jLv0K/vNUJbU9C48GqYESuwkzvnGAIqxPVb6cSL2CVXlRVxVjh9rGKovLIo8ep8zfy1vmT/OS1bzL+0CM82vw4Dm3lbraNwKrdDoTrhmQyVWA0lqsKRzSHIqDBa6c54KAp4MBpXXkqsYnJ/eBe/yIfEkKkZv8uAMfs9wKQUsqNXM8vpUJSSpnF3J97SSwWP4XCOIZRQVHu/XCxNNRjaainEo9T7O1l5uZN8rYpKrV+Knu70YLLB5fvhrfGge3hapvywYszNHT4qG1ehfjYPND+zmq9xWTPPMHYjebqphI1MDJlFEVZsVDcorX7AB5PkDOnv8fg6VNE09M81vEUTe61VY6vF1WpriSa/A4ek5JIpshovLraODUY59RgnKDLSnPAQXPAYTYsNNkW3PXpIqVce6nr6hkF5jufm4Hx+3j9Bw6LJQBIKpUkVmvNis/TAgG0J54gvSeJ1tuPOj5MemQIra4W2+7dWMLhFddi3MLmtNB5pI7x3gRTA0ny6RLhbj/qalwqdu88wbgIkV5UvReinejTEiHEqoUCINAY4h1Pv8iVMz+hv/cmP8p/h+7WQzxU/xAWZev235pf3He0NUAyV2Y0UV1p9Iwm6RlN4rKpVWEJOKj3mHENk61hO61xTwHdQogOYAz4KHDPwPXbmapQQKkcW5VQ3KJot2NtbaHu0B4cY6MUb9wg+5PXUZxObLu6sHZ2Lmhvfi9UVaFlf5CZ0QxTA0kK2TIt+4LYXat8GNu90P4UhB5CTF1GnRpBvzSB8AlksH11c83i9Pl56B3H8Vzw0T/Wy/X8GSL5CI82Pkqto3ZNc240PqcFn9PHgbCPfElnLJFnPJGnP1LdH1xTBSGfnbDfQdjnwGG9n+9xJm9ntkQohBB/BbwHqBVCjAL/YXYf7l8BXqWaHvsnUsrLW2Hfg4KmuVEUK+VSHFz3Pv5Osnr17dztsFcD3N3d1cD3jevkey5SuHwZS0srtu5dq3JL1Ta7cXgsjF6rxi1CXT4CjWsw0OaG1sdR9U7K1/shex6y02CLVoPeq2w8aLHb6T72JPYrbsZGbzJ2Y4TvF9LsrzvAvpp99y2NdiU4rCq76t3sqndT0Q2m0kXG4lXhGInlAQi6rDT5HYT9doIuq5l6a7JprHiHuweJt0N67C2qu90Vqa9/YdXnvpXI0J8v8tHQ4tWInkxSvHmT0uAgsqKjBgPYunZhbW1BaCt7v6iUdEZ742QTRXx1DkK7VumKujVPokihN8aINo0oT7HXEYdyDhyBai+pYMeqW5tPD/Yzdv0qkzJKKqxR623g0dCj9z2Ndi3Es6W51cZMprqPuN2iEPJVYx9mFpXJWlnXDncPGm8noUilekinrxAK/dyKAtrz+e5MkrKU/IM6/7LHyFKJ4uAgpf5+9GQKYbFgbW/D1tmJ6l/+vLnzpWRmJMP0UAqrXaN5XwCHe3UBWlkxyJ6ZYsiYRgva2dPdXa30nroE+ThYnFC/t5otZVl5n6fUzDTDl3qIFxNM15VQPA4O1x2my9+1Kvu2kkJZZzyRZyJZYCJZoFQxEAJq3TbCfjthn4OA2VLEZIWYQrFDyedHicVeo7b2WWy2ulWd+7eTMRrnbYV6LyqRCMW+PkojI2BItNoarB2dWFuaEZa7xyGyySKj1+LoZYP6di81Ta5VuUpyF6bpi41ga/Kwe/e8zYqSYzB1ubqPt6JBsLNan3GX9iDzKeayDF44SyadIOIvkvZLGp2NHGs8dt/TaNeLYUii2dKscOSJZavZ7Q6rQqO36qJq9NmxaWZsw2RpTKHYoeh6nsnJb+D1HcHj3rvi8/K6wd9Ora7P0y2MYpHSrVVGKo3QNKxtrVg7OtBqlg+qV8o64zcSpKMFXH4bTXsCWFYYkC3ciHPt5nVcu4J0dy9ROJeLwfRViPWBoYMnVBUMf+s9O9bqlQojVy6Smp4i49SZCOaxWKwcrT9Ks6d5RfZtR/IlnfFknsk7VhtBl5Wwz0HIbyfotJqbMpnMYQrFDmVi8iVuXP8s5UoCuy1MZ9enCTV++K7nfHUyxm/3jTNZqhCyavy7rjA/27i2GopKJEKxf4Dy6Eg1luHzYm1vx9rWhuJwkDxxgunf+yKViQm0UIj6T30S/dh7mexPoqiC8C4/k31J3nipj0ysiDto48kPd7H78cYF1zn3P17j/NULDIkZ3Iqddx99J49+6J2LDSoXqpXekd5qA0Kru7rHd+3uqluq52/g5Gere2X4muH4b8LhnweqcYvJvuvoVoVzsQucP/dDbNNFwkaQc91JzteO0+hq5BNHP8EHOj+wps9rq5CyutqYSBQYT+aJzsY2LKqg0Wcn5HMQ8tlx2W67L79xbozPv9rLeCJP2O/g157fw0eObE39icn9wRSKHcjE5Etcu/YbGEZ+bkxRHOzd+zvLisVXJ2N8uneEvHH79+5QBF/Y07JmsYBqLKM0OkppYIDKTBQElIaGif33/w7F4txxwm4n9NufxX78BUZ74wxdjHL9rUn0ym17NKvCM/9075xYnPrmj3n99BvYpZWIkqQiDDSp8Pwj711aLAAMA5LDMH0N0hPVYHd8CN74A6gUbh9nccCLvz8nFplYlJ9+9Sv0fP87RJ05ZvwltIpCMGXh/K4Eg0157Kqd33rHbz1wYjGfQllnOlWcW3HkStXWKB67RthvZyCS5d+/dJl8+XbLFIdF5XM/c8gUix3MckKxo1IjhBAvCiG+nEwmt9qU+0J/3xcWiASAYeS5eeNz6HpxyXM+1z+xQCQA8obkc/0T67JFWK3YOjvxHD+O9/0vYN+7l8RXv7pAJABkocD0732xWqD3UB2Dl2YWiARApWTwxkt9c9//8OyP0VBRENw6siIMfnj2x8sbpCgQaIc9L8CBj1S3aT39xwtFAqr7fZ/87Ny37mANfedOIXWdYNpG81S179RUsEhdshooL+gFvnT2S6v6fLYbdotKa42TJzpr+MiRJj5wOMQjbQHcdo0bUxn+47evLRAJgHxZ5/Ov9m6RxSZbyY4SCinlCSnlx30+31abcl8oFJd+uJfKESanvsFM9Aek01coFqcwjGpgc6xYXvKc5cbXgur14jh8GCORWPLnlYmq3UIR5FNLXzcTuy0wumEQNNyUhI4ubldmZ2RhqVMX4whUe0rll+mLlRy949rRub87SyodE06aIg58mdtumYnsBDfjN8mVcyuzYZvjc1jY0+jhmT31PN5ZQyxbWvK48UR+yXGTnc12qsw2WSV2W4hCcXGXE5u1EbdrD4XCKKnCLTERWCx+nrKk6Cu5iBGkiH0u0Ntk2/hWFlo4TGV8sX1a6PYOd+6gbYEozN2Ds9po0O+xsMtoJCuKROfajlVxCCu6rqOudK9wVx1kI4vHfQsD1g6vj3wyMfe9QOAuLPyv4rf5ORc5x7nIOXxWH2F3mJArRNAefKAL3yq6wUQij99hIZFfLOJh/4OVCWayMeyoFcXbjc6uT6MoC//jKoqDrl3/Gz7fwzQ0fJDGxp+hpubdeDz7URQb/7w2z2FxhYc5Qze9HJTn2Sv6+d9bFAxj6bfItVL7//lluKM4T9jt1H/qk3PfP/nhLjTrwn+GqkVw+D3NpIZTjH1/hHDNLsaNNIa47aJSpWBXayeXL18msczKZQHT12D3C9UW5/OxOKoB7VkmbvTSdeRRVG2hcFYUgzN7qisSu2rn1x/7dZ5vf57DtYexqlauxa7x9yN/z7f6v8WZqTNMZifRjXu3RN9O5Es637s6xWA0x79+ugP7HUV7DovKrz2/Z4usM9lKzGD2A87E5Ev0932BQnECuy10z6wnKSVfHx/mv/RfZ6Zcplub4RfqShz1VltsaBY/NmsdVmsNVmstmrbC3euWIH/xIvG//msyP/4J+tTUXNaT78UXFxx3/c3JRVlPXfuC5K5EyRYrJK0aN89cp3f0Glk1gUOTvPuRd7LvvQ8zNDREPp8nEAjQ0tKC1bpEcVkxA1e+Aa76qvvp7397yayn6NgIY1cvE2xuIR2Z5rWv/Dnp6Ayaz8WZ3Ym7Zj0V9SKT2UnGM+NMZiepyAqa0Gh0NdLsbqbR1YhF3boGhPciXSjz/d4IhZLOU921NPkdZtbT2xAz68lkEX8zGSNss/AOn41SKUqpNFP9KkeRszENRbFjtdXNikctFksAsYKeSFJKUidOoAYCuJ9+elV2GYUK+UtRhFXBsa8GVEEykmdmNEMxV0azqgTDLoIhF0KBqakpJiYmEELQ3NxMbW3tQvdP399XC/MOfKTa0nwJcskEN0+/iTsYpOOhR1bdPXc+uqEznZtmLDPGRHaCgl5AQaHeWU+zp5mwO4xN3T57aacLZb53dQrdgPfsqTP3+X4bs9Yd7kx2MI1WC1PFMoriwW4PYbdXYwdSVluXF0uRqnAUIxTy1T2lhNCwWmurX7Y6rJbaJVuH6NEoRr6A46HV7TYoDUnhRgIE2PcEEZbqA9vf4MTf4CQTLzAzmmF6MMXMSJpAo4vapnoCgQBDQ0MMDQ0Ri8VobW3F4XBANlpNiw0/vKxIGLrO8OUeLDYbbYceXpdIAKiKSsgdIuQOVesXClHGMmOMpcc4PXUaMSWod9bT5G4i7A5vaQV4vqTz/d4IugHP7qs3978wWZIdJRRCiBeBF3ft2rXVpjwQNNosDBWKJMsVfJbb/xSEqAa+LRY/uKpV0Lqeo1iMUJoVj3T6EqQBoWC1BLDaGuZWHYpipTwxAQIs8wLXK6E8lsHIlbHvCaDYFgep3QE77oCdQqbMzGia6FiG2HgWX72D9pZOUtkEo6OjXLlyhXA4TEP2KopqhfoDy15zqv8mpVyOzqOLYxPrRQhBraOWWkctD9U9RLwQZywzxmhmlLPTZzk3fY46Rx3NnmaaPc33daVhGJIf3ai6m95rioTJXTBdT29jMhWdr0/HOeZ1sW+1bTyMEqXSTHXVUZymVI6BNACBxRpAvziCprupOf4zKMrKHkBGoUK+Zwa1xo69y7+ic0qFCtGxLPHJLNKQeGsc+BptROKTxKMRHFNnad99ENeedy19fj7HtddfI9AYouXA4RXe/caQLCYZTY8ykhkhXUojEHPuqWZ3M1Z1cx/c54bjXJ1I885dtbTWODf1WiYPBqbryWQRbk3Fq6mMFkqrFgpFsWK3h7HbwwAYRoVSeYZScZpicZps7iZqTYDSxNewWoOzK456rNa6ZbvclsezIMDasvIAutWuEeryUdfiJjaRJTqWJRXN4/L7CDsSRHSda1FJw+gooVBoUSrt9NAAQggau3Yvc4XNw2fz4bP5OFB7gGQxyUh6hJH0CGemznBu6hwhV4gWbwshVwhtlZ2B70U0U+TqRJpd9W5TJEzuiSkUb3NabFauZguUDAPrOnzziqJhtzVitzViFIuo2RnUrjaEx0WxOE0mc42MvDLrqqrBZqvHZgthtdYghIKsGFSiebQaB8oadm7TrCr1bV5qmtzEJ7NERzNkJzO4jEbK7hCTk5MkEgna2trweKpCpFcqJCbG8Tc0YrGvvD35ZnBLNA7WHiReiDOcHmY4NcxYdgxNaDS5m2j1ttLgbNiQOo2zwwnsFoWHW/zrN95kx2MKxducFoeVy9k8o4Uync6N8Y8buRwCBYenFau3mk5pGOVZV9U0xeIk6fQV0unLCGHBZqtDyweQug17/dr7TQGomkJts4dgyEX89dNE4wGMpAu7aiVXjtLb20t9fT1NTU2kozMYuk6wqeXeE99HAvYAAXuAw7WHieQjDKeGGc2MMpQewq7aafO20eZtw2dbWweC6VSBSLrIYx0Bc4MjkxWxo4TCDGavnlqLhlNRGMoXN0woZLmaWivmBcgVxTIvs+ohDKNEsVgVjUJxgkysD90ok09fxV5uxGYLYbc3oihrs0lRFWo8GYItQRJ2P5GRNDJXTzYXZzgzRiKRwFYqomgaTu/2bPkiRDVmUe+s54hxhInsBEOpIW7Eb9Ab7yVgC9DmbaPV27qqIPjN6QwWVdBes4btaU3eluwooZBSngBOHDt27Je22pYHBSEEbQ4bvdkCRcPAts7UUIC5zn0s7yJRFCsORzMOR7V9RmZmmKIrgrAVKBTGyOUGAIHVGpwVjRAWS83q3C5GGaHZCDS68Nc7SUbyRIatpJIpIkMR0slBGhqC6IaBthH3vYmoijqXGVWoFBhJjzCUGuJ85Dw9kR7C7jAdvo57uqaklIwnC7QEnWjq9r5nk+3DjhIKk7XR4bBxNZtnKF9it2v9vnpxawe1VbSwUAwnLnsH9mAQKSXlcpRCYYJicYJ0+hLp9CUUxYbN1jgbRA/de7UhlDkbhCLwNzjx1TlIRjxMDbq5NNzHWDGB+tYF9hzYxYPSTNKu2ekOdNMd6CZRSDCYGmQoNcRoZhSn5qTd106HtwOnZXGQOlWoUKoY1HnMojqTlbOjhMJ0Pa2NGquGT9PozxU3Rihs1YeQUVhhd1eqz/RbmdpCiLmiPjiErhcpFicpFicoFMbJ54cAgdVWh90Wxm5vwmLxLp5Us8Md3V0XCMZkkELZSWw8yxvT52jbFaZ7byea9uD8t/Db/Txsf5hDtYcYz44zkBzgSvQKV6NXCblCdPo6aXQ1zq0ycqUKUN13wsRkpeyofy2m62ntdDltnE1lFxXfrQXFWX2TNTLZFZ8jLCpGobLkz1TVhtPZhtPZNm+1MU6hME4qdZ5U6jyq5sZhb8Zub8ZqnW3hYfdDIbH09RSBy28nXFeHw99G76V+bl4aYXwowoGHdtPYvLr9x7caVVFp8bTQ4mkhU8owkBxgMDXIeHYcp+Zkl38X7b52KnpVja2m28lkFewooTBZO50OG+dTOW7kihzzre+fhVAUVK8HfSVdXWdRXBYq8QKybMy17Vhy7nmrDa/3MJVKlkJxnEJhjEz2OpnMNRTFht3ejMMmsCViiFIOrIvdMBa7nXKxQHPYwxONhxkbmOZazw1O/fgCDaEGDjzcjcu3tWmza8FtdXOo7hAHag8wlhmjP9FPz0wPl6OXcSkN5CpeSpX6rTbT5AHCFAoTAByqQovdSl+uyMMeJ5qyvlx9NVhDeWIcKeWKAtCqzwqjUEkUsNStvABM01y4tW7cru7ZTKpJ8vlR8oVhcnoCUb6IfdSNI/QUNlvjgmI/u8tNJh4DQFEELV0NhFpruNrTx1DfKNHvxOna1UH73hBWx4P3X0URytwqI1lM0pfo41q0j97MJU4OjfKujsOEXeEHev8Mk/vDg/ev/y6YMYr1scdlZ6hQZCBfpHudsQpLYwOlwUH0eBwteO/aCNVtRXFoVKZyqxKK+VQzqVpxOFqR0qBYnCKfTFGIXSRvMxBCw24P43C0YreHcQdrSUxOkEsl51JkNYvGoUf20NoV4uLZa/Rev87U+Aydu9tpaPOiraEYcDvgs/k42nCUQ7WHSCZ+ymRmjNfHX8dtcdPt76bd177h1d8mO4cd5ah8u22FutE02CwENI1r2ZUHoZdDa2wEAeWRkRWfY2l0oWfLVGLrv74QCnZ7iEDTB2mUndRY9uBwtlMsThGL/ZiJya9jaAPoIk58YnTR+T6/lyff/QgHjnVSEhkuXrjIpdeHiAynMXRjiSs+GFhUC0ca9xPWnuSRusexqlbORc7x8sDLXIleobjMXusmb292lFAIIV4UQnw5mUxutSkPLPvcdhKVCuOF9e12p9hsWEJhSkNDSGNlD1at1oHi0CgNp5D6BjWrrNmFcPiwTw8S8D5CY+NHqKl5Dw57CyV9Es07yvjY14hF36RUii44VVVVOrs6OPbUYeraPETSo/Re7Of6qSkSUzke1Iaauxs86IagkPdzvPU4z7Q8Q429hsvRy7zc/zLnp8+Tr5h7Y5vcxuwea7IAQ0q+NhXHq6k8V7u+lVlpdIzsT36C6x1PYm1ZWZsMPVUifzWKVutYcQfZexIfqm5e1HQUQg/NDUupE5/ppf/CqzhrVXx1dWgWH05HB05nO6p6u1FipVJheHiYybFpSmmB31GPx++iscOHy//g1SR8/9o00WyJDx4OYbdU3WnJYpLeWC/D6WEUodDp62R3YPeS9RgmO5PlusfuqBWFyfpRhGC/28FUqUykVF7XXJamMKrXQ+Hq1RWfo3qtWJvcVGbylKdWnl57VwJtEGiH8fPVjYxmEUIlWLefhvCzlKNdOGwHUISFVOo8k5MvEY3+kHx+FCkNNE2js7OTvQd24wtZySpTxGIzDF6cYfhKlFJ+6dTe7crRtgAV3eDMUHxuzGfz8VjoMV5of4FWTyt9iT5eGXyFc9PnzBXG25xtLxRCiE4hxB8LIf52BcearqcNoNtpx6YoXEyv7+EghMC2dy96PEFpNbGKJjea30ZxKEUlukEPqNYnweKorizKC+ds3LUbTbMz05empua91Nd/ALdnH6VynFjsNSanvkkq1UOlkiEYDHLgwAHqQkHwZsmrURIzWW6emWZqIIX+gMQvfA4LB5t8DEVzXJ9KL/iZ2+rmWOMxXuioCkZ/op9XBl7hYuQiJX19LkmTB5MtEQohxJ8IIaaFEJfuGH9BCNErhLgphPh1ACllv5TyF1cyrxnM3hgsimC/y85YsbTuVYW1vR3V5yN/oQdZWdlbtxAC264AqttKsS+5IcFtLHboOg6VPPR9H/TbtlisNloOHKKQTjPeexWLxYvP+xCNDR8iGHwXVkuAdPoyU1MnmIn+AMOYYdeuXbS2tqLYK5RdMyjO6o57N09Pk5jK3cWQ7cOBsJew387ZoTiTycWfscvi4ljjMZ5vf56wO8y1+DVeHniZ3lgv+iras5g8+GzViuJPgRfmDwghVOAPgfcD+4GPCSH2r2ZSc0WxcexxObApChc2YFXhOHIEI5ulcOXKys9TBfbdARSXRuFmfGPcUK4aaH8aMlPQ/32YF2T31tZT195BbGyU6cH+WdsVHI4mamreTUPDh/B4DlAuJ4hGf0gk8m1criR79uzC4bCTMSJYaguoFsHY9TgDFyIUMusT2c1GCME7umrxOiz86HqE6dTSguy2unk89DjPtT1HraOWnpkeXh18lZH0yleJJg82WyIUUsofAbE7hh8Dbs6uIErAV4APr3ROIcTHgf8AHC2VzOXxerEogoNuBxPFEpPFdcYqGuqxdrRTuHaNSuzOX/vyCE3BvrcGzWejOJiiOJxCGutMvgh2QNtTkByF/r9f0LiwsWs3voZGJm9eJzq28CGoaS683sM0NnyIQOBJhKKRTJ4hmfoO4XCRujofqWycgjVCoNlKMV+h79w0EzcT6JXt646yagrv3VuP06byg+sRJpLLvxj4bD7e2fROnm56Gk3R+OnET/n+8PdJFs0Xs53OdopRNAHz/3eOAk1CiBohxH8DjgghPrPcyVLKL0spj0kpj9XVPVh9erYru512XKrK2VR23amgjocfRnE4yb7xxtx+FStBqAJbdwBLg5PyRJbCtRiyvE63R91uaH0CEiNw47tQqb5YCCFoOXAIT00tY1cvEx1d/MYshILT2U593fPU1b0Puy1ENtuLxXKOUChDpZJmIjaMtxkCjU5iE9X4RWJ6+7qj7BaV43sbcNs0ftgboS+Suevxja5G3tf2Ph6pf4R0Oc13h77L+enzlPXtvYIyWTvbSSiW6iMgpZRRKeUvSym7pJSfu+9WvY3RFMFDHgfRcoXB/DrrKqxWXE88jpHNkjt1alXnCkVga/dh6/JjZMvkLs5QSawzblG/DzreVXVD9b4MxWpAV1FU2h46gqe2jrFrl5nqv7nsFFZrLcHgUzQ0fACHswMhpgkG+7HZBxmbuEFGRmg56MdiVRnrjTN4cWbbZkc5rCrP7mugwWvnzf4YZ4ZiGHdZvQkh6PR38nz783T4OriRuMErg68wnhm/j1ab3C+2k1CMAvOT7ZsB81/dFtPpsBG0aJxL56is0+2j1dXhOHyY0sjoquIVt7DUOnAcqEFYFAq9cYr9SeR63Do1XbDrWShl4eq3ID0JVMWi/fARAqEwU/03Gbncg3GX4K2meQj4H6Wh4UV8vv0EAwYeby/xxJv0DfYQ7LAQ6vKTT5e5eXaayHB6/S60TcCqKbx7dx17Gj30Tmb47tUpMsW7C5tNtfFIwyMcbz2OTbXxk/Gf8Mb4GxQqG5CAYLJt2LKCOyFEO/AtKeXB2e814DpwHBgDTgH/REp5ebVzmwV3G8tUsczv9k9wMpoiUq7QZLPwmc4QP9u4tv2tsz/9KaWhYYxcjthf/AWViQm0UIj6T30S34sv3vN8aUjKYxlKExkUi4q1zUtxKEXq1UH0RBHVb8P7fDuuIyvskJpPQN/J6qoifBSmr8Lf/zYkR5nSOplq/hkc+99LIZ3mja99hXR0Bk9NLU9/9BfY9/Qzi6bT9TyZzFXi8SvMzESoVBqw2srEZv6IbMRHPtFEOtdCuOldzFxUycSKuIM2nvxwF7sfb1zlp7k5jMRy/LQ/igSOtga4NJbk86/2Mp7IE/Y7+LXn9/CRI00LzjGkwbXYNa5Gr2JVraSKKf7syp8xmZ2k0dXIJ45+gg90fmBrbshkRSxXcLclQiGE+CvgPUAtMAX8BynlHwsh/gHwRUAF/kRK+Ttrmd8Uio3lq5MxPnVtmNK8fyoORfCFPS1rEgup60x/4QvE/uK/w7yUWWG3E/rtz65ILAD0TIniQIpCb4zcuWmY1/ZDWBT8P9O9crGolGDoJ3Dpq9Dz1zCvXiApfbzl+Fku9QxgzEur1aw2nvv4rywpFgC6niORvMD163/D+PgppCHxeGeQRS+Z6T1Mje2jOLMHWyk4O5/CM/9077YRi0yxwpv9UV69PMk3zo1Tmlcj4rCofO5nDi0SC6hWeP/X8/+Vv+79ayry9udlV+381jt+yxSLbcy2qsyWUn5MShmSUlqklM1Syj+eHX9ZSrl7Nh6xJpEw2Xg+1z+xQCQA8obkc/0Ta5pPqCrJV15dIBIAslBg+ve+uOJ5VLcVx4EaCr2xBSIBIMsGqVcHV26UZoWuZ6rB7TuKynwiSf+lqwtEAqBSKvLaV/58eftUJzXBJ9Erb+F0xlC1Cul0DRWliL/1LRo73qJsTcybz+CNl/pWbvMm47ZpvHdvPT/ojSwQCYB8Wefzr/YueZ7P5uPk8MkFIgFQ0At86eyXNs1ek81jO8UoTLYpY8ukxy43vhL0ycklxysTqxMfoQiM9NJ26Ik1dELNTC09XF66vXg6OnPPKYulKazWIi5XHE0rUyy6yOS8eEOXKNqnF14ntr26twohiKSXtmk8sXwq7VRu6c9xMrv0791ke2MKhck9abJZVjW+ErRQaFXjd0NdpimfcGqUJ7OrCxz7mpccbnGlgcXzODxeKveo27HbqvekqjpudwyLJY+uW0gmGtGsC9NmbS6NXGp71QGF/Y4lxwMuK9HM0iLS6FrafVZjr9kwu0zuH6ZQmNyTz3SGcNyx451VVMfXSv2nPomw37E5kqYR/Gf/dNVzeZ9vX7x9qiZwHm2gOJQifyFSFYyVtC4//pvVnlDzUa08fSzEbn8Kh3r7Ia5qFnY98hi9b7xGdHRk2VqTzq5PA1UxEwJcrjQ2a4Zksp56/wxazTVAoloE3Y80MHAhwvTQBhQXbhC/9vweHJaFKyqbpvDCgUZevTzFj2/MkMwtXNV94ugnsKsLf78WxcLTzU9zavIUhty+RYgmizG3tDK5J7cC1p/rn2CsWKbWovJsjY+nAp41z3krYD39e1+kMjGB2tiI74XnUVxu8pcv4zhwYMVz3QpYL5X1pCeLlMYyFIdSlMYyWBpdWBqcCG2Zd6TDP1/98+Rnq9XbvmY4/puEmh7haO2f4XvzR0zE82Qc7bzjo/+KjiOPMnrtMmPXLhMdGybcvRd3cOFbc6ix2mDg2rXPoesRikUXkcg7CPl+junYBdTaGzhq0jz6xD9i92PtTPYliQynycSLNO0JYNvibVhvBazvzHr6B4dCXJtMcW0yzXAsR2vQyaEmHz6nZS5g/aWzX1qQ9dTua+dK9AolvcQToSdQlQdzx8C3G+Z+FCarpmJIvhVJIIEP1vmxrHN/7VtIwyD31ilKQ0NY29pwHnsEoW3MQ1JPlyiPZ6gkighFoNU7sTQ4UeyrmL9SgsmLMHUJhAINB6DhIGhWElOTTNy4RrlQwFtbR8Ou3Tjc9xZSKSWDgz9lYuJHWKwe9u75OTyeepKRPBM3ExiGJNTlI9DoWsfdby6Fsk7vZJreqTQVXdIccLA/7KXWvbRL8Gb8Juci52hwNvBU+ClTLLYR2yo9drOYt2f2L924cWOrzdnRTBfLvBpNsstp50m/e0Pnzl++TOHSZdSAH/dTT6G4Nu4haeTKlCay6NECUko0vw1LowvVt4rNhwopGDsD8UHQ7BA6DHV7MSREhgeJDA1gVCoEQmEaOndhddx7459odJCbfd/AMCSdHS/S0LCLclFnrDdONlnEV+cgtMuPutxKaBtwSzCuT6Up65I6j419IQ9NfgdCLHyZGEwOcmrqFCFXiHeE34Eitu99vZ14WwjFLcwVxf3hbCrL5Uye9wS9tNitGzp3eWyM7JtvgSJwPf44ljUEue+GUdKpTOcoT+eQZQPFoWGpd6LVOpZ3S91JdgbGzkJqDKwuaDwEtbupVHSmB/uJjg4jpSQYbqK+vfOegpHLRent/Vvy+RSNjc/R3n4IIQQzIxmmh1JY7RrN+wI43Bv7WW80Zd2gP5Ll2mSKbFHHY9fY3eCho9aFdd5n25fo4+z0Wdq97Tza+OgWWmxyC1MoTDYcXUq+HUmSMww+WOfHqW7sW6GeSpF9/Q30ZBL7/n3YDxxAKBt7DWlIKtE8lek8eqaEUARq0I6l3onqWeEDOTUB4+eqqbUW55xglMsVpgb7iI2NAhAIhalv78LmXF4wyuUMN278LYlEFK/3abq7H8Jms5FNFhm9FkcvGzR2+giGt68r6haGIRmO5eidShPNlNBUQWeti+4GDz5HNWPuSvQKl6OXOVR7iL3BvVtssYkpFCabQrJc4e8iSeptFo4HPYtcDOtFVirkz52j2D+AVhPE+cQTqO6NdXXdQs+WqURyVGbySF2iODS0WgdarQPFugI/+gLBcFQbD9bto1wxmB4aIDY2gpQGvrpG6to7cHqX3mCrUkkzNPRNIjMxbNYn6ejYhd/vp1LWGetNkIkX8NU5CHf7UTZYnDeLaKZI71Sa4WgOQ0K9x0ZXvZuWgIMz06cYTg/zrqZ30eBq2GpT39aYQmGyaVzPFngzmeGIx8lBz7398WuhNDxM7vQZkBLHkYexdXZuynUApG5QiRWoRPLo6RII0Hw2tFoHasCOuFfwPjVRDXgnR0G1QO0eqN9HGY2ZkSGioyMYlQquQJC6tnY8NXWLBLZYnGZi8lXiMQ1d7yYUChEOhwHmXFE2p4WW/cEtz4paDYWyTn8ky81IhkyhglVTaA1aGSm+hdUiea79OWzqKuJFJhuKKRQmm8qPYmmGCyWeq/FSv45CvLthZLNk33qLynQES6gR57FjKHdx42zINQsVKpE85Zk8sqQjVAUtWBUNxWO9+woqF6tmScUHqt8H2qH+ALrdT2xslJnhQcrFIlank9qWNgKhJtR5WV6p9CWSiQsUCl0kkxo+n4+Ojg40TSMTLzDaG0ca0Lw3gCdoX9qGbYqUkul0kZvTGUZiOTLlJKOlUxxq6OID3e/EbjEzobYCUyhMNpWSYfByJElFSj5Q58exSS4RKSXFGzcoXLwIQsHx8ENYOzo23OW11HWNVInyTB49XkDqEmFV0WrsaEE76t0CzMU0TF+DmevVPlLueqjfh/S2kpiZZmZkiHwyiaJpBEJN1DS3YHe5kdJgevrbVCvCjzE6OobVaqWrqwun00mpUGHkSoxCtkx9u5e6lrXXtWwlhbLOcCzHDwZPcyNxnT2eJ+iuDdNR4yLst6M9IO61ncDbQijM9NitJVau8EokSa1V49kaL8omPrz1TIbcW6eoRCJodbU4jx1D9Xo37XrzkbpETxSoRAvoiSJSShSbihq0o9U4UF3LrKj0MszcgOkrVfGwOKB2N9TtIZcvMzMyRGJqEqTEFQhS09yK1V0invgJweDT6Lqf/v5+KpUKbW1t1NTUYOgGY9cTpGby+OudhLr9KBtU13K/KRtlvt77d2TyFgLiIfIlA00VNAcctNW4CHntD+y9PSi8LYTiFuaKYuvoyxV4PZHhgMvBUd/mZuZIKSkNDJC/0IOslLHv3Yd9394NK9JbkQ0Vg0q8QCVWwEiWqqJhV1ED1ZWG4rIsXu1IWU2pnb4GyZFqXw9fSzVTylFHfGKc6NgI5Xwe1aKh+m7iq2sn1PQC5XKZvr4+MpkMDQ0NNDU1oSgKkeE000MpnF4rLfuDaA+o66Y31kvPTA/PND+DXnEzOJNlOJajrEtsmkJL0ElL0EGDxxSNzcAUCpP7xpuJDNdzBZ4OeGh3bH5g0igUyJ8/T2loGMXpxHHkCNbmxfskbDZzohEtYKSqoiGsKlrAhhawV2Madz7cimmIXK+6pSoFsLqhdhcyuIt0pkBsbITozE+pyEmCvvdT09SCt66B8YkJpqen8Xq9dHR0YLFYSEbyjF2PY7GqtB6seaCC3LcoG2W+1fctWjwtHGusPq90QzKRzDMUzTEWz1MxqqLRHHDQEnTSaK40NgxTKEzuG4aUfDeaIlqu8HyNjxrr/XlglaenyZ89h55MYmlswPHQQ6h+/3259p3IioGeKFKJz7qnDInQFFSfDS1gQ/XZFhb2GQYkhqquqdRYdczbBLW7SFFmfPQ7GJlOKjkVoSr46xsx7A5mEqkFcYtcqsTwlSgArftrcHq3d3HeUpyaPMVoepQPdX1oUXuPim4wkSwwEssxmshT0SVWTSHst9MScBLymTGN9WAKhcl9Ja8bfHsmiZSS929CMd5ySMOgePMmhctXkOUSts5O7AcPotzZqfY+InWJnpwnGhUDIQSK14rmt6H6bQt7ThXTVcGI9kEpQ0mUiVinqQl/CF1tIT4+RmJ6EqNSQReClC6xe7zs3ruPQCBAKV9h6FKUckmnZV/wgcuIGsuM8fr46zzT8gy1jtplj7u10hiJ5RlL5ClVDFQFGn0OmgMOmvwOM3tqlbwthMIMZm8vYuUKr8wkCWgq76vxod1H94BRLFK4fIXizRsIVcO2dw/2PXvua/xiKaSUGJny3GrDyFd3gVMcWnW14bfddlFJCelJytPnmI68SlBtw+FohZouDF8bqXSR+OQ48elJpuNJDEWhtaOT7v0HUDU7Q5ejFDJlmvYE8NdvbhrxRlLSS7zU9xIHaw6yr2bfis4xDEkkU2Q0nmM0nidb1AGodVsJ+6uiEXA9eKur+83bQihuYa4otg/D+SI/jKdps9t4OuDe9DTWO9FTKfI9FymPjaE47NgPHKim025wK5C1YhQqVdFIFG/HNRSB6quuNFSfjRIRotPfo1bbhS2VgPTsLoCuWgh2UXaFic/EuH7tKtGZGZw2G62tLfjrw6RjVoo5CHX5H4i2H7c40XeCRlfjmntAxbMlRuPVlUYsW91DxGlVaQo4CPsdNHhspotqCUyhMNkyrmTynEll70sm1HJUIhHyPT1UZqIobheOgwextLbed+G6G1KX6OkieqL6Zcy+FRe0IbLWqzSGP4TVH0DoeYgNQKyvWtQnBLgbIdjBaErhZl8flVyWoMOOIgT5jAWh+Gg/3Elj5/KunO3EyeGTWBQL72p+17rnypd0xpN5xuJ5JlMFKrpEVaDeY6fRZyfsc+Bzbk6R6IOGKRQmW8qtTKhHfS72upbeWvN+UB4bI3/xEnoyierzYj94EEtT07YSjFsY+Qp6skhk6u8pZ5PUqMersQ23pbri8FlRlAwiPlit/i6kQAiSai39aQu65qbG66WUSjJxc4pCpkR9R4iW/e346upX1P58q/jhyA8xpMEzrc9s6Ly6IZlOF5hIFphIFEjmqzvzuWwqjV47Yb+Deq8Nm/b2jG0sJxQPXv7cXZgXo9hqU0zu4DGfi5xucCqZxaEotN2HtNmlsDQ1oYXDlEdHKVy6RPYnr6P6fNgP7MfS3LytBENxaBiWItLI4Ws7gEMJoidL1V37RtMwCkJVUL0dqLV7UC05lPwwvvgAe7UpbkYKxBNe2nbto3XPu7hxboTI8BiFzHncfjt2jwdfXQPe+oYVbbJ0PzGksSm/C1URhHwOQj4HtEK2WKmKRjLPcCxHXySLEBB0WWn0VlcctW4b6ts8/dZcUZjcNyqG5LvRJPGyzvEaLw2b1BNqpUgpKQ8PU7hyBT2Vrq4w9u3D0tKybWIYicRpsrk+GhteRFVvrwBkWUdPlarCkbrtphJWFdVjRbVkkaVRBvovk85kCPlthJvaGI01EEtbcdUIkClyyQQAVocDb10D3rp6XD7/lt//twe+TcAe4InQE/ftmoYhmckWmUwWmEwWiGZLSAmaIqjz2mj02mnw2gk4lyii3CGYrieTbUFBN3g1miSvGzxX6yNo2fpFrZSS8shIVTCSKRS3C/vevVjb2xHq1rkgSqUYkch3cLl24fcv+r+7AKNYmRMNPVVClg0AhFVhLDlEPDVCjT1Lh7/C2LiDdMlH0/4GXE1hUtkKqcgU6VgUpES1WPDU1OKtrcdTW4uq3V9BrxgVXrr5EnuCezhYe/C+Xns+pYrBVKrAVKrqqkoXqhlqVk2hwWujYVY4bu2tsRN4WwiFmR77YJCt6LwaTaFLyfO1PrzbxB8spaQ8Nk7x2lUq0RiKw47joYewtrVtgS0Gkch30I0CDfXvR1FW56ozcmX0dKm66kiVmIpNM5GYxm230NXoYCYyQy6WoaU5j6fGCY2H0QOdZGIzpCLTpGYi6OUyCIHLH6Bl/8H7FtOI5CL8YPQHPBV+irA7fF+uuRJypQpTqeKceNxKwXVYFVoCTo61B7fYwvXztohRSClPACeOHTv2S1tti8nyuDSV40EPP4pn0LfRi4oQAmtzE9bmJspTUxSuXgV1a/6LGEYZRbHi8RxctUgAKE4LitOCpaGaZdaeC+IZCzIxMoFav4uWI48wdH4c6U2CGAchUDUNX30jvvpGpJTkkglSM9NkolE06/2LKRnSIGALUOOouW/XXAlOq0ZHrUZHbfUzTRfKTKWKTKcKGNvnn/GmsKNWFLcwXU8PBlLKHevr3a7M/8zNz9/kTpZbUWyPiN0GIYR4UQjx5WQyudWmmKwA8yF1/5n/mZufv8lK2VFCIaU8IaX8uM+39F7EJiYmJiarZ0cJhbmiMDExMdl4tn0wWwjxEeADQD3wh1LK7yx3rBnM3hl8dTLG5/onGCuWabJZ+ExniJ9tfDAzSpInTjD9e1+kMjGBFgpR/6lPknvUoL/vCxSKE9htITq7Pk2o8cNbbSo9PT288vJ3yOUzKLqNIN089+I72f144wZe5G/g5GfJRrtIGf+SihEkak3yJ7Vf51p4lE8c/QSdY05e+8qfk47O4Kmp5emP/gL7nt7YCu2dxjfOjfH5V3sZT+QJ+x382vN7+MiRjduTZUuC2UKIPwE+CExLKQ/OG38B+BKgAn8kpfzdeT8LAF+QUv7iXeY102MfcL46GePTvSPk56WROBTBF/a0PHBikTxxgol//5vIQmFuLPcOhdQ/1TFEeW5MURzs3fs7WyoWPT09vPSNb6IblduDhoIvt4cP/KN3b4xY9PwNnPg3ZPOPkaj8KpLb7c8LosiXQn/JeO4mT12qQZZv26FZbTz38V8xxWIZvnFujM987SL5sj435rCofO5nDq1aLLZbMPtPgRfmDwghVOAPgfcD+4GPCSH2zzvk383+fFnMGMWDz+f6JxaIBEDekHyuf2KLLFo707/3xQUiAZB+f2GBSAAYRp7+vi/cT9MWcfLkyYUiAaAYpO39vPFS3wZd5LNQzpOq/PMFIgFglzb+xfSHOXzVtUAkACqlIq995c83xoYdyOdf7V0gEgD5ss7nX+3dsGtsiVBIKX8ExO4Yfgy4KaXsl1KWgK8AHxZV/hPwbSnl2eXmFEJ8XAhxWghxOhKJbJ7xJpvKWLG8qvHtTGVisbjpyyyKCsWtFcLl4nqGWiQTK27QRUYB0Fm6g21dJYirsHTxZTo6szE27EDGE/lVja+F7RTMbgJG5n0/Ojv2q8CzwM8JIX55uZOllF+WUh6TUh6rq6vbXEtNNo2mZfo/LTe+ndFCoUVj6p2vR7PYbYuPvZ8stwpXdBvu4AYV2/maAVBZ+qEf0WJk7fqSP/PUPBjt0beCsH/pbszLja+F7SQUSyV1Synl70spH5FS/rKU8r/dd6tM7iuf6QzhuKNTp0MRfKZzax+ka6H+U59E3LEFq+fbdhS5UPQUxUFn16fvp2mLOH78OKpyR26LoeApdPLkh7s26CK/CRYHXu3PECx0yRVEkT+tf4mefVnEHf2/NKuNpz/6Cxtjww7k157fg+OOLV8dFpVfe37Phl1jO2U9jQIt875vBsa3yBaTLeJWwHonZD35XnwRYEHW065/+ElyB7Zf1tPhw4cBFmc9/aMNzHo6/PMAuE5+FqJ/sCjrqTc8xieO/jqdT5lZT6vhVsB6x2U9AQgh2oFv3cp6EkJowHXgODAGnAL+iZTy8mrnNlt4mJiYmKyebZX1JIT4K+ANYI8QYlQI8YtSygrwK8CrwFXgb9YiEiYmJiYmG8uWuJ6klB9bZvxl4OX7bI6JiYmJyV3YTsFsExMTE5NtiCkUJiYmJiZ3xRQKExMTE5O7siM3LhJCRIChNZxaC8tUA20/TFs3ngfFTjBt3SweFFs3y842KeWiiuUdKRRrRQhxeqnUsO2IaevG86DYCaatm8WDYuv9ttN0PZmYmJiY3BVTKExMTExM7oopFAv58lYbsApMWzeeB8VOMG3dLB4UW++rnWaMwsTExMTkrpgrChMTExOTu2IKhYmJiYnJXXlbCYUQ4lNCiMtCiEtCiL8SQtiFEEEhxHeFEDdm/wzMO/4zQoibQoheIcTzm2zbnwghpoUQl+aNrdo2IcQjQoiLsz/7fSHEUvt8bIatnxdCXBNC9Aghvi6E8G9XW+f97NNCCCmEqJ03tu1sFUL86qw9l4UQ/3mrbV3m9/+wEOKnQojzsztNPrbVds5eo0UI8X0hxNXZz+8Ts+Pb6v/W/7+9Mw25ogrj+O+PZqiVZWW2UkQGKmKapdFiC1FSWeCXKFCKSAsrPxSFEBYVaYtBQXtkq0sJWVCaQbZraLlVUubbItIG5fLBpZ4+nOf2Xm/vHb3me+fQfX4wzJkzM2f+d7jPPPecM/d5CnTmYVdm1hILKVveOqC7b88GxgHTgFu97lZgqpf7A8uBfYHjgLVAl07UdyYwBFhVVdewNmAJMIKUCOpN4MImaT0f6OrlqTlr9fqjSZGKvwMOyVUrcDawENjXt/uUrbWOzgWV6wCjgHfL1unXOBwY4uX9SakM+udmWwU6s7CrlupRkKLldlfKfdGDlBhpNDDD988ALvXyaGCmmW01s3XAN6S83p2CdZxHvCFtkg4HDjCzjy19Y56rOqdTtZrZAkuh4gE+ISWeylKrMx24Bah+myNHrROAe81sqx/zc9la6+g04AAv96I96VjZ93SDmS3z8iZSCoMjycy26unMxa5axlGY2XrgfuB7YAPwh5ktAA4zsw1+zAagj59SL4d3M2lU25Ferq1vNleRfslAhlolXQKsN7PlNbuy0wr0A86QtFjSIknDvD43rTcB90n6gWRnt3l9NjqVkqWdBCwmY9uq0VlNaXbVMo7CxyBHk7ppRwA9JV1ZdEoHdbm8S1xPW+maJU0GdgAvVqo6OKw0rZJ6AJOB2zva3UFd2fe1K3AQMBy4GZjtY865aZ0ATDKzo4FJwNNen4VOSfsBrwI3mdnGokM7qGua3no6y7arlnEUwHnAOjP7xcy2A3OB04CfvLuGrytd+xxyeDeq7Ufau6bV9U1B0ljgIuAK7/ZCflqPJ/1YWC6pza+7TFLfDLXi155riSXAX6SAcLlpHUuyKYA5tA/Tlq5T0j6kh++LZlbRmJ1t1dGZh13trcmY3BfgVGA1aW5CpHHJicB97DypNc3LA9h5suhbOnEy2695LDtPEDasjZRrfDjtE1mjmqT1AuAL4NCa47LTWrOvjfbJ7Oy0AuOBO73cjzTcoLK1dqDzS2Ckl88FluZwT73t54CHauqzsq0CnVnY1V7/oue8AHcAXwGrgOf9Jh8MvAN87eveVcdPJr1NsIZOeCOjRtvLpLmT7aRfBVfviTbgZP98a4FH8H/fN0HrN/4Q+9yXx3LVWrO/DXcUOWoFugEv+LWXAeeUrbWOztOBpf7wWgwMLVunX+N00tDLiqrv5qjcbKtAZxZ2FSE8giAIgkJaaY4iCIIg2APCUQRBEASFhKMIgiAICglHEQRBEBQSjiIIgiAopGvZAoKg2UiqvBoJ0Bf4E/jFt08xs22lCOsASSOBbWb2UclSghYmHEXQcpjZb8BgAElTgM1mdn9ZeiR1tfbAb7WMBDYDu+0oJHUxsz/3hrYggBh6CgLgnxj+iyQtlTS/KrzDu5KmS3rPcwUMkzTX8xjc5ccc6zkDZnjegFc8ptSu2r1H0iLgRkkXe+C/zyQtlHSYB4cbD0xSyvNwhqRnJY2p0r3Z1yM9n8FLwEpJXTyXwaeu6dqm3tDgf0U4iiBIoQ4eBsaY2VDgGeDuqv3bzOxM4DHgNeB6YCAwzoexAE4EnjCzQcBG4DqP3VPU7oFmdpaZPQB8AAw3s5OAmcAtZtbm15xuZoPN7P1dfI5TgMlm1p/0b+k/zGwYMAy4RtJxjd+aIIihpyCAFMplIPC2JwPrQgpRUWGer1cCq83DU0v6lhSY7XfgBzP70I97AbgBeGsX7c6qKh8FzPIeRzdSkq1GWWIpNwGkhDeDqnofvYAT9rDdoMUJRxEEqUex2sxG1Nm/1dd/VZUr2xUbqo2FUwn5XNTulqryw8CDZjbPJ7Cn1DlnBz4S4OHGu9VpT8BEM5tfp50g2G1i6CkI0sP/UEkjIIV7ljSgwTaOqZwPXE4aSlrTQLu9gPVeHltVv4mUGrNCGzDUy6OBfeq0Nx+Y4MNfSOonqefuf5wgaCccRRCknsEYYKqk5aQonac12MaXwFhJK4DewKP+mu3utjsFmCPpfeDXqvrXgcsqk9nAk8BZkpaQQudv+VdLiadI4amXSVoFPE6MIAR7SESPDYL/iL+d9IaZDSxbSxB0BtGjCIIgCAqJHkUQBEFQSPQogiAIgkLCUQRBEASFhKMIgiAICglHEQRBEBQSjiIIgiAo5G/oBBFGCmuK0QAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots()\n",
"\n",
"for dof_vals in pm.TruncatedNormal.dist(mu=1000, sigma=500, upper=8000, lower=800).random(size=(10, 7)): \n",
" pp = PiecewisePolynomial(\n",
" element_number=3, polynomial_order=2, \n",
" domain_boundaries=[1e2, 1e-3], dof_values=np.sort(dof_vals)[::-1]\n",
" )\n",
" log_p = np.linspace(1.8, -2.8, 100)\n",
" \n",
" ax.scatter(pp.dof_values, 10**np.asarray(pp.dof_vertices))\n",
" ax.semilogy([pp.getValue(x) for x in log_p], 10**log_p, alpha=0.4)\n",
" ax.set(xlabel='Temperature', ylabel='Pressure')\n",
"ax.invert_yaxis()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "copyrighted-italic",
"metadata": {},
"outputs": [],
"source": []
}
],
"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.7.10"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
@bmorris3
Copy link
Author

And a jax friendly version:

import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
from jax import numpy as jnp
from jax.ops import index_update

gl_0 = [0.0]
gl_1 = [-1.0, 1.0]
gl_2 = [-1.0, 0.0, 1.0]
gl_3 = [-1.0, -0.447214, 0.447214, 1.0]
gl_4 = [-1.0, -0.654654, 0.654654, 1.0]
gl_5 = [-1.0, -0.654654, 0.0, 0.654654, 1.0]
gl_6 = [-1.0, -7.65055 -0.285232, 0.285232, 0.765055, 1.0]

quadrature_nodes = [gl_0, gl_1, gl_2, gl_3, gl_4, gl_5, gl_6]


class Element(object): 
    def __init__(self, edges, order):
        self.reference_vertices = quadrature_nodes[order]

        self.nb_dof = len(self.reference_vertices)

        self.dof_values = []
        self.dof_vertices = []

        for i in range(0, self.nb_dof):
            self.dof_vertices.append(self.referenceElementMap(self.reference_vertices[i], edges[0], edges[1]))

    def lagrangeBase(self, r, i):
        l = 1

        for j in range(0, self.nb_dof):
            if (i != j):
                l *= ((r - self.reference_vertices[j]) / 
                      (self.reference_vertices[i] - self.reference_vertices[j]))
        return l

    def getValue(self, x):
        # coordinate on the reference element
        r = self.realElementMap(x, self.dof_vertices[0], self.dof_vertices[-1])
        
        y = 0

        for i in range(0, self.nb_dof):
            y += self.dof_values[i] * self.lagrangeBase(r, i)  

        return y
            
    # maps the coordinate value r on the reference element [-1, +1] to the real element [x_l, x_r]
    def referenceElementMap(self, r, x_l, x_r): 
        return x_l + (1.0 + r)/2.0 * (x_r - x_l)
    
    # maps the coordinate value x on the real element [x_l, x_r] to the reference element [-1, +1]
    def realElementMap(self, x, x_l, x_r): 
        return 2.0 * (x - x_l) / (x_r - x_l) - 1.0
       
        
class PiecewisePolynomial(object):
    def __init__(self, element_number, polynomial_order, domain_boundaries, dof_values):
        self.nb_elements = 0
        self.nb_edges = 0
        self.elements = []
        log_boundaries = [jnp.log10(domain_boundaries[0]), jnp.log10(domain_boundaries[1])]

        self.nb_elements = element_number
        self.dof_vertices = []
        self.nb_edges = self.nb_elements + 1
        self.order = polynomial_order
#         if (polynomial_order < 1): order = 1
#         if (polynomial_order > 6): order = 6
        self.createElementGrid(log_boundaries)
        self.setDOFvalues(dof_values)

        
    def createElementGrid(self, domain_boundaries): 
        domain_size = domain_boundaries[0] - domain_boundaries[1]
        element_size = domain_size / self.nb_elements

        element_edges = []

        element_edges.append(domain_boundaries[0])

        for i in range(1, self.nb_edges-1):
            element_edges.append(element_edges[i-1] - element_size)
            
        element_edges.append(domain_boundaries[1])
        
        for i in range(0, self.nb_elements):
            edges = [element_edges[i], element_edges[i+1]]
            self.elements.append(Element(edges, self.order))

        for i in range(0, self.nb_elements):
            for j in range(0, self.elements[i].nb_dof-1):
                self.dof_vertices.append(self.elements[i].dof_vertices[j])

        self.dof_vertices.append(self.elements[-1].dof_vertices[-1])
        self.nb_dof = len(self.dof_vertices)
        
    def setDOFvalues(self, values):
        if len(values) != self.nb_dof:
            raise ValueError("Passed vector length does not correspond to the number of dof!\n")
        
        self.dof_values = values

        # set the dof values in each element
        self.global_dof_index = 0

        for i in range(0, self.nb_elements):
            for j in range(0, self.elements[i].nb_dof):
                self.elements[i].dof_values.append(self.dof_values[self.global_dof_index])
                self.global_dof_index += 1

            self.global_dof_index -= 1 # ; //elements share a common boundary
    
    def __call__(self, x_vector):
        x_lowers = jnp.array([self.elements[i].dof_vertices[-1] for i in range(len(self.elements))])
        x_uppers = jnp.array([self.elements[i].dof_vertices[0] for i in range(len(self.elements))])
        element_bools = jnp.where((x_vector < x_uppers[:, None]) & (x_vector > x_lowers[:, None]), True, False).T

        element_vals = jnp.array([[self.elements[i].getValue(x_vector[j]) for i in range(len(self.elements))]
                            for j in range(len(x_vector))])

        values = jnp.sum(
            jnp.where(element_bools, element_vals, 0), 
            axis=1
        )

        return values
    
from jax import jit


def piecewise_poly(log_p, domain_boundaries, dof_values, element_number, polynomial_order):
    pp = PiecewisePolynomial(
        element_number=element_number, polynomial_order=polynomial_order, 
        domain_boundaries=jnp.array(domain_boundaries), dof_values=jnp.sort(dof_values)[::-1]
    )
    return pp(jnp.asarray(log_p))

ppj = jit(piecewise_poly, static_argnums=(3, 4))

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment