Simpson's Rule on the Logarithmic Scale

simpson-log-scale.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "39e0a9b5-6b36-44d8-8329-d200a2a2c5ad",
   "metadata": {},
   "source": [
    "Recently at [work](https://monetate.com/) I have been prototyping an idea that requires numerical integration of probabilities and began to suspect I was facing issues with floating point arithmetic in certain circumstances.  To work around these issues, the standard approach is to work with on the [logarithmic scale](https://en.wikipedia.org/wiki/Log_probability) where quite small probabilities are represented by reasonably-sized negative numbers.  This post shows how [Simpsons's rule](https://en.wikipedia.org/wiki/Simpson%27s_rule) for numerical integration can be adapted to work with log scale data.  While the idea is straightforward, I made a few mistakes in the details in my first attempted implementation, so I want to record the correct approach for my future reference. \n",
    "\n",
    "Recall that (a basic form of) Simpson's rule works by dividing the interval of integration $[a, b]$ into many subintervals of length $h = \\frac{b - a}{n}$ with endpoints $x = a + i \\cdot h$ for $i = 0, \\ldots, n$.  On sliding subintervials $[x_i, x_{i + 2}]$, the integrand, $f$, is approximated using a [quadratic interpolation](https://en.wikipedia.org/wiki/Polynomial_interpolation) of the points $(x_i, f(x_i))$, $(x_{i + 1}, f(x_{i + 1}))$, and $(x_{i + 2}, f(x_{i + 2}))$.  Integral over $[x_i, x_{i + 2}]$ of the resulting parabola can be evaluated exactly and used as an approximation of the integral of $f$ if $h$ is small enough.  Combining the integrals over these subintervals with appropriate weights gives the [composite Simpson's rule](https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_1/3_rule),\n",
    "\n",
    "$$\n",
    "\\int_a^b f(x)\\ dx \\approx \\frac{h}{3} \\left(f(x_0) + 4 \\sum_{i = 1}^{\\lfloor n / 2\\rfloor} f(x_{2 i - 1}) + 2 \\sum_{i = 1}^{\\lfloor n / 2\\rfloor - 1} f(x_{2 i}) + f(x_n) \\right).\n",
    "$$\n",
    "\n",
    "We can rewrite the right hand side as\n",
    "\n",
    "$$\n",
    "\\int_a^b f(x)\\ dx \\approx \\sum_{i = 0}^n w_i \\cdot f(x_i)\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "w_i = \\begin{cases}\n",
    "    \\frac{h}{3} & i = 0 \\\\\n",
    "    \\frac{4 \\cdot h}{3} & 0 < i < n \\text{ is odd} \\\\\n",
    "    \\frac{2 \\cdot h}{3} & 0 < i < n \\text{ is even} \\\\\n",
    "    \\frac{h}{3} & i = n\n",
    "\\end{cases}.\n",
    "$$\n",
    "\n",
    "If we want to approximate the logarithm of the integral using values of the function $\\log f$, we can rewrite the above as\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\log \\int_a^b f(x)\\ dx\n",
    "    & \\approx \\log \\sum_{i = 0}^n w_i \\cdot f(x_i) \\\\\n",
    "    & = \\log \\sum_{i = 0}^n \\exp(\\log w_i + \\log f(x_i)) \\\\\n",
    "    & = \\operatorname{LSE} (\\log w_0 + \\log f(x_0), \\ldots, \\log w_n + \\log f(x_n)).\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Here $\\text{LSE}$ is the [LogSumExp](https://en.wikipedia.org/wiki/LogSumExp) function, which is [frequently used](https://en.wikipedia.org/wiki/LogSumExp#log-sum-exp_trick_for_log-domain_calculations) for increasing the accuracy of log-probability calculations and is [available](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.logsumexp.html) in SciPy.  In this case, the log weights are\n",
    "\n",
    "$$\n",
    "\\log w_i = \\begin{cases}\n",
    "    \\log h - \\log 3 & i = 0 \\\\\n",
    "    \\log 4 + \\log h - \\log 3 & 0 < i < n \\text{ is odd} \\\\\n",
    "    \\log 2 + \\log h - \\log 3 & 0 < i < n \\text{ is even} \\\\\n",
    "    \\log h - \\log 3 & i = n\n",
    "\\end{cases}.\n",
    "$$\n",
    "\n",
    "We now implement this log-Simpson's rule in Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d1bb986b-3b87-4976-bd25-2ae2f7716d03",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3fa628ce-7e46-4b72-b232-41b963d62690",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.special import logsumexp\n",
    "from scipy.stats import norm\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4e464158-757b-43d4-9fb8-20f0b2c7faf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(color_codes=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "63073ccf-ec47-466f-ae97-8ed2053ceeea",
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_simpson(log_y, dx, axis=-1):\n",
    "    log_w = np.full(log_y.shape[axis], np.log(dx) - np.log(3.0))\n",
    "    log_w[1:-1:2] += np.log(4.0)\n",
    "    log_w[2:-1:2] += np.log(2.0)\n",
    "\n",
    "    shape = [1] * len(log_y.shape)\n",
    "    shape[axis] = log_y.shape[axis]\n",
    "    log_w = np.reshape(log_w, shape)\n",
    "\n",
    "    return logsumexp(log_w + log_y, axis=axis)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44b55634-53ec-4536-a788-ea6b8f6db60c",
   "metadata": {},
   "source": [
    "We validate this approximation on the log PDF ($\\varphi$) and CDF ($\\Phi$) of the normal distribution.  First we define an evenly spaced grid of 300 points on which to test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "191658ba-eac4-4bc4-8a2a-6ce1f769ade7",
   "metadata": {},
   "outputs": [],
   "source": [
    "Z_MAX = norm.isf(1e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "68377be6-4f36-4270-991d-437c6cfffe2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 300"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3a1ede7b-b274-4d3e-baad-a57200cbc7b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "z = np.linspace(-Z_MAX, Z_MAX, N)\n",
    "dz = z[1] - z[0]\n",
    "\n",
    "log_φz = norm.logpdf(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb99d49c-e60c-43da-a28e-c4698900034f",
   "metadata": {},
   "source": [
    "Since $\\Phi(0) = \\frac{1}{2}$,\n",
    "\n",
    "$$\n",
    "\\log \\int_{-\\infty}^0 \\varphi(z)\\ dz = \\log \\Phi(0) = \\log \\frac{1}{2},\n",
    "$$\n",
    "\n",
    "which is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a452d2f9-2460-4e77-a2e3-e547cd333c06",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.6931471805599453"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.log(0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3175d84a-5e68-4dc9-b291-bb29acb2bd0e",
   "metadata": {},
   "source": [
    "The approximation of this quantity using our log-Simpson's rule is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9e96160d-2bf7-4cd5-882a-90ab953147e3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.6889293664385137"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log_simpson(np.where(z < 0, log_φz, -np.inf), dz)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10c16686-4b8d-4e51-8b34-3fbfccec8452",
   "metadata": {},
   "source": [
    "which is not too far off the true value.  Note that here we use $-\\infty$ as the log integrand for $z \\geq 0$ because $\\log 0 = -\\infty$.\n",
    "\n",
    "To gain further confidence in our approximation, we can approximate\n",
    "\n",
    "$$\n",
    "\\log \\int_{-\\infty}^{z_i} \\varphi(z)\\ dz = \\log \\Phi(z_i)\n",
    "$$\n",
    "\n",
    "for each of our grid points $z_i$ for $i = 2, \\ldots, n$ as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "91f8d7ea-00f2-487b-aaa9-be345e883e14",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_grid = np.tril(np.ones((z.size, z.size)) * z)[2:]\n",
    "log_φz_grid = np.where(z_grid == 0, -np.inf, norm.logpdf(z_grid))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "8c6e3868-b12f-495b-afbd-5cff711d3ae4",
   "metadata": {},
   "outputs": [],
   "source": [
    "cum_logΦ_approx = log_simpson(log_φz_grid, dz)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92dbc9db-530a-4226-aa5a-9cf5b20b7820",
   "metadata": {},
   "source": [
    "We see that, once we have accumulated a reasonable number of intervals, our approximation is reasonable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "93c015ed-3f26-4175-b1d3-7731e19a6634",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAAG1CAYAAAAcMztGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABTKklEQVR4nO3dd3gVZeL28e8p6SEhoUMoGnqv0ouiSBFF1EUE1LXuqyvrulZ01Z9tC6zr2leBVVER1wYCSl9BRKQJIk1aIEAgJCEJSU6def9IiEZ62pxyf66LCzJzzpk7D4ecm5lnZmymaZqIiIiICHarA4iIiIgEChUjERERkRIqRiIiIiIlVIxERERESqgYiYiIiJRQMRIREREpoWIkIiIiUkLFSERERKSEipGIiIhICafVAYKRaZoYRuhfMNxut4XF9xmINPbW0dhbQ+NunXAZe7vdhs1mO+vjVIzKwTBMsrMLrI5RpZxOO0lJceTlFeLzGVbHCSsae+to7K2hcbdOOI19cnIcDsfZi5EOpYmIiIiUUDESERERKaFiJCIiIlJCxUhERESkhCZfVxHDMPD7fVbHKDfDsOFyOfB43Pj9oX+2ghUcDid2u/5vIiISSFSMKplpmuTlZVNUdNzqKBV29KgdwwjtsxSsFhMTT0JC8jmdQioiIlVPxaiSnShF8fFJREZGBfUHnsNh096iKmKaJh6Pm+PHcwBITKxlcSIREQEVo0plGP7SUhQfn2B1nApzOu0hf10LK0VGRgFw/HgONWok6bCaiEgA0E/iSuT3+4GfP/BEzubEeyWY56OJiISSsChGhmHw4osv0r9/fzp37sztt9/O/v37q2x7wXz4TKqX3isiIoElLIrRq6++yvvvv8/TTz/NBx98gGEY3HbbbXg8HqujiYiISAAJ+WLk8XiYPn06EydOZNCgQbRu3Zp//vOfZGRksHDhQqvjiYiISAAJ+WK0bds2CgoK6N27d+myhIQE2rZty5o1ayxMFjwKCo4zeHBfRo4cgs93fnNhNm36no0bv6+0LIcOHaRfv+6sX7+20l5TRETkhJA/Ky0jIwOABg0alFlet27d0nXl4XSe3CkNI3Tmi5yY+mKzweLFC0lKSiY7O4uvvlrK4MFDzvl17rrrNiZNeoJOnTpXTdAQ4XDYSt9TDkfZ36X6aOyrn2maGCYUFHnJL3DjdRXh93jw+dwYXg+m14Pf68H0eTB8PkzDj9/vxzT8GH4/hmFglnyNYQBG8e9m8S/TNKHkV/GfS5bx83IwyzyueJlR/ByKH3qa9L/6Xk6/zvaLlWe+CIp5ij+e62VTzFM89Fc5TvEsh8OO32+cx3aqlrPOBXQcOsq67Vu25WpSVFQEQGRkZJnlUVFR5Obmlus17XYbSUlxJy13uRwcPWov8yEX7BwOO/Pnz6F3775kZBxizpxPuPzyoef1GnZ75Y3HLz+4QmGMDcOG3W4nMTGW6OjoMusSEmIsSiUa+5MZXg+Febnk5+RQlJeLp/A4HpcLr9uF1+XG53HjL/ll+DwYXg/4PGB4sfu92E0fdqP4d4fpw2n6cOAnAh8RtuLfHTYTG8UfTCH/4SSn5c/7nqjIUcTGnfw5Wx1C/r134sPG4/GU+eBxu93ExJTvh59hmOTlFZ603ONxl9wKxCxz/R/TNPF4rbseUGSE/bzPfrLZisvHrl27+PHHzdxww43k5+fx178+w+7de2jSpCkAPp+Pt96ayhdfzOXYsRyaNbuQ3/3ubnr06EW/ft0BeOaZJ1m3bi233HIH1113JS+++DpduxavO3ToYJllHo+HN998jf/9bwmZmUeIiYmle/eLuO++h0hKSir5Xw34/UZIXGPJ7zcxDIPc3EKKioov9+Bw2ElIiCEvr6j0+5XqEU5jb3hcFB47Sn5WFkU5WbjysvEXHANXHg53Pg5fEQ7DRaThJsp0EWnzl3l+ZMmvcrFx6l0Xv+Az7Xhx4sNR8suJYXNgYsO02TGxY9rsYLP94s8nlhcvw2bDtNl+3qDt59/Nkt9tv/hz6WNsZR9/8o9PW5k/madeddLjfv0SZslS269XnOHL0w6c7Qzb4udtnfw0Gw6nHb/PCJD9RRDTMJXaHnB7Cir1dRMSYs5pb3DIF6MTh9COHDlCkyZNSpcfOXKEVq1alft1T/WhfKqrRJumyV/eXc/OA+XbO1UZmqck8si4rudVjk7s9Z07dw4xMbH06tUHt9vNlCl/Zfbsj7nnnvsAeOGFKfzvf0v4058eomXL1sydO5uHHrqPt956n9mzv+Sqq4YyceKfGD58JPn5eWfd7quvvsjKlSuYNOkJGjRoyM6dP/Hcc//HO+9M5w9/+FO5vv9g8OsyXbwsNMpfMAr2sTeK8inIyiA/6yhFRw/iy83ELMrF6ckj0ldArFFAlM0LQHTJrzMq+dFhmDaKzEiKiMJji8KwOzHsEZglv2zOCHBEYHNGYo+IxO6MLPNne0QkjogoHJFROCIicUZG4YyKIjo2hlp1alLgMjBszuLXsAX/HuFg4HTaSUqKIyenIKDe81ZmCfli1Lp1a+Lj41m9enVpMcrLy2PLli2MHz++ekIE6dQjn8/HggXz6ddvAFFR0URFRXPRRb354ot53HHH3fj9PubNm8299z7AxRdfCsCdd94NQEFBAU2aNAMgPj6e+Pj4cypGbdq05eKLB9OpUxcA6tdvQI8eF7F7986q+SZFgpThKiDvwC7yDu2n6NhRbMcOEuXOIsp/nDhcANQo+XWSkp9JbtPJcTOGQnscHmc8/sgaGFEJ2KJr4IyJxxlbg6j4GkTH1SA2IZG4hBokRDor/fpbTqedmKQ4XDkFmAH04SzhKeSLUWRkJOPHj2fKlCkkJyfTqFEjJk+eTP369Rky5NwnEZeXzWbjkXFdg+5QGsA336wkOzurzGTrSy+9nG++WcGyZYtp1uwCvF4v7dp1KPO8E+WoPC6/fDhr1qzmtddeYv/+fezbt5d9+9Lo2LFzuV9TJJgZruPk7t3KsYNpuLMPE52/n3hvFlF4cQBJJb9+LdeIoZAYCpw1cUfXwh6bhDO+JlEJycTUTKZGrTokJNagdmTIfwyInJew+BcxceJEfD4fjz32GC6Xix49ejBt2jQiIiKqZfs2m42oSEe1bKsyzZs3B4BHH33gpHWzZ3/Mffc9XOFtnLiNygmTJz/HsmVLGDZsBP36DaBFi9uYOfNdjhw5XOFtiQQy0zTxZqWTlbaTgoz9OLN2Eus6QgwunEDtUzwn2x9Hlr0W/qgEjBr1cSY3JD6pFgkNUqiTXJOoiOD7uSNitbAoRg6HgwceeIAHHjj5A15OLScnm5UrVzB8+Eiuv35cmXWzZr3PvHlzsNlsOJ1Otm37kebNW5Suv+OOmxk8+DLGjCn7vBNFtLDw5wl1+/fvK/1zbu4xZs/+hP/7v+fK7KXau3cPsbGxlfr9iVjJNE1MVz55+3eRu3M9RtZ+YoqOEEch8UD8rx5/xJ9AlrMuZmwyJDUmpn4zatVvQMO6yTQNgbMzRQJJWBQjOX8LFszH7/czfvxNpXOFTrjxxlv44ou5zJnzCddcM4Y333yNmjWTuOCCVObOnc3u3Tt57LEnAYiJiWXv3j3k5h6jVq3aNGjQkA8/nEnjxk3JzT3Gm2++VnqYLy6ueC7SihVf0apVG9xuNx99NIsdO7bRtm37ah4Bkcpjmib+7P0c3bWN/PSd1MjeSpyRjwNI/sXjPKaDA0ZtXFHJeBKbEd0wlbopTWjSoBapTu39EakOKkZySvPnf06PHj1PKkUAjRql0L//QBYu/IKPPpqLw+Fg8uS/cPx4Ps2bt2Ty5H+VPu/668fx/vvvkJa2h7/97Z889thT/OtfU7j55rE0atSYiRPv44EH/gCA0+nk6af/yssvv8CNN15PQkICXbt2584772bGjLdwuVzVOAIiFeMvOEbW5m8o2L+dyGNpJBjHiANOXJnFMCHHiOeQMwV30oUkNmxC/eat6VAnUTcXFrGQzTTNQLl0QdDw+w2ys0++voLX6yEr6xC1ajUgIqLcV/gIGE6nPaBO3wxFp3rPBOrps+GgImNv+n0U7N9KzsavsGftJt6bjf0X/cZr2tnrr4c7tj5mgzYkN+9I05TaxETp/6d6z1snnMY+OTlO1zESEalK/uPZHN26Du/2FcQXpuPAoOaJlTbY569NRkxzYus3o16bLnRMqYNTtxsRCWgqRiIi58g0DfzHMjiyeTX2nSuI82bzy9MCCoxIfrKn4m3QiUYtW9MytQntNDlaJKioGImInIFpmvhzDnF440oidi8nxp9fetFEw7SRYdTkQHxbopr3pGXr5vRP0hmUIsFMxUhE5BSMwmMc/X4FbFtEjC+PhJLlHtNBulGbIzU7kdyxP+1bNKRVEF6nTEROTcVIRKSE6fdxaM0Sjq74hLiiDE7cZrq4DNUhs1ZXancaQLvUOnTS6fMiIUnFSETCnufwLrJWfkLU0e1E4COO4tPpD/qTSUvsSt2uF9OhRX0iVIZEQp6KkYiEJdN1nNx1X1Cwcx0J7ozSq03nG9FsjupMdPvBdG3fjDYx1XPrIBEJDCpGIhI2TNPEd2ALR35YRVT6GqJMNwkUT6L+3t+cotRBXDykD8MjHCF/TRcROTUVIxEJeabPgyt9K9nfziEhb1fpROqDvppsjetBSqee9Gl/ITHRztKL3YlIeFIxkrMqKDjOlVdeTmxsHJ9+Oh+nM/DfNs8++ySHDh3k5ZffqLJt7N69i4yMQ/Tp06/KtiEVYxoGud8vxlj/CRGGiwSKr0C9wdccW0on2vYbyKhav75lq4iEs8D/hBPLLV68kKSkZLKzs/jqq6Vl7nwfqP7wh/sxDH+VbuOhh/7I0KEjVIwCkFF4jGPffIx/7zqijUIcQK4Rw06aYWs/jD492upWHCJySvrJIGc1b94cevXqQ0bGIWbP/iQoilF8fNXvBdBtBgOPP/sA2T+swL5jGZGmmwiKr0a9NqoXDXoPZ0DLetjtukGriJyerlUvZ7R37x62bNlMjx49GTToEtavX8u+fWml66+9diRvvTWV++77PZdc0pcxY0Yxd+5npevnz/+cq68ezpw5n3LVVUO57LIBPPLI/Rw9mlnmNV5++QXGj7+OESMGs2HDOvx+P7NmvcfYsaO55JI+jB07ms8++6j0Of/859+57LIBZGQcAqCoqIjrrx/NQw/9ESg+lPb7398BwPr1axk4sCdffbWs5PX6MnHi7zh8OIMXXpjM0KGDuOKKy3j77Wmlr+/xeHjllX9x3XVXMmhQL4YNu4Q///lhcnJySjNnZBziP/95s3Q7x48f529/e5YrrriUyy8fyMSJv2Pbti2V/Dcip2K6CziycBqFHz1K9PYviTTd7PPVYkGN0RSM+CtX/PYWurWur1IkImelYlQNTNPE9Lqt+1WBPRvz5s0hJiaWXr36MGDAxTidTmbP/rjMY95+exrt23fkrbfeY/To6/j7359jyZKFpetzcrL58MOZPPXUX3nxxdc4cuQw9933e3w+X+ljPvnkQ/7wh/v5xz9eol27Drz88gu89dY0fvvbO3j77Q8YPfo3/Otf/+DDD98H4K67JlKnTh3+/vfnAHjppedxuYqYNOmJU34ffr+fd96ZzhNPPMOLL77GTz/t4OabbyAiIpI33nibUaOu4c03X2PXrp0AvPrqi/zvf0uZNOkJPvjgUx599EnWrVvDO+9MB+DNN9+hbt16XH/9eJ57bjKmafLAAxM5ePAAf/vbC7zxxtu0a9eB//f/bmXHjm3lHn85MyM/k6y5L3Ds7YnE7F0BwGZPCisSRlDj6j9z7dgrSW2cbHFKEQkmOpRWxUzTpHDOsxiHd1qWwVGvBTFXTsJmO7//Lft8PhYsmE+/fgOIioomKiqaiy7qzRdfzOOOO+4mKioKgB49enHLLcV7TZo0acaWLZv58MOZpYfcfD4fjz32f7Ru3QaAxx9/mnHjrmXdujX07NkbgF69+tKjR0+geLL3p5/+l3vu+SNDhgwFoHHjJhw6dIAZM97iuuvGEhUVzeOPP82dd/6Wv/zlKb74Yi4vvPAqiYk1T/v93Hbb72jdui0A3br1YMuWzdx110RsNhsTJtzMW29NZffunaSmNqdNm7ZcfPFgOnXqAkD9+g3o0eMidu8u/ntMSkrCbrcTExNDQkIia9d+x+bNPzBv3mISEhIBuPPOu/nhh438978f8OijT57X2MuZGfmZ5PywHNuPC4g0PQAc8tdkW50hdB80kN614yxOKCLBSsWoGtgIzt3333yzkuzsrDJzii699HK++WYFy5YtZujQEQB07dqtzPPat+/EN998Xfp1bGxcaSkCaNq0GTVqJLBr187SYpSS0rh0fVraXnw+Hx07di7zup07d+PDD2eSk5NNcnItWrduy4033sJ//vMmv/nNWLp27X7G7+eX24iJiaFBg4alZTEqKhoAr9cLwOWXD2fNmtW89tpL7N+/j3379rJvX9pJmU7YsWMbpmlyzTVXlFnu8Xhwu91nzCXnzjQMcld/hu2Hz4mkeE/obm8dttUbysCLe3FlLRUiEakYFaMqZrPZiLlyEvg81oVwRp733iIoPowG8OijD5y0bvbsj0uL0a9P3zcMP3b7z0dpT3V6v2H4cTh+fsyJvU8ApzvyZ5rGSa+3fftWHA4HGzasw+v1EhFx+qsU/zrHmcZk8uTnWLZsCcOGjaBfvwG0aHEbM2e+y5Ejh0/5eMMwiIuLY9q0d09ad6ZMcm5Mn4fC9fM5/uNXxHqL53n95K3HwcROdBg8gjENEi1OKCKhQsWoGthsNoiIOvsDA0hOTjYrV65g+PCRXH/9uDLrZs16n3nz5pQeVtq6tewE4x9+2ETLlq1Lv87Ly+XAgXQaNUoBiq//U1BQUOYxv9SsWTOcTiebNn1PixatSpdv3LiBWrVqUaNG8eX5PvvsY9au/Y4XXniVxx57kGnT/s3vfvf7Cn/vubnHmD37E/7v/54rs7ds7949xMbGln79y2J14YXNKSgowOv1csEFF5Yu/9vfnqF58xZcc82YCucKV97MPeQs+DcxhRnEAoVGBF9HD6LDsBFc0bim1fFEJMSoGMkpLVgwH7/fz/jxN9GkSbMy62688Ra++GIun31WPAl78eIFtG3bnosu6sWKFf9j+fJl/O1v/yzznKeffpx7730An8/H88//lfbtO9K5c9dTbjsuLp6rrhrN1Kn/JiEhkTZt2rF69So+/fQj7rjjbmw2G+np+3nllRf47W/voEuXbvzhD/fzzDNP0KdPv9Me7jpXcXHxxMfHs2LFV7Rq1Qa3281HH81ix45ttG3bvvRxMTExpKfvJzs7i549e9OiRUueeOIR7r33AerWrcenn/6X+fM/5/nnX65QnnDlz9rPscVvEpm7jxiK72H2P1sv2g66jNGtGpZrL6iIyNmoGMkpzZ//OT169DypFAE0apRC//4DWbjwC5xOJ8OHj2T58v/xyisvkJLSmKee+gu9e/ct85whQ4bxwAN/wOv10LfvAP7wh/vP+MF2zz33kZhYk9dee4mcnGxSUhrzxz8+yJVXXo3f7+fppx+nSZNm3HDDhNLXX7x4Ac888wRvvTWzQt+70+nk6af/yssvv8CNN15PQkICXbt2584772bGjLdwuVxER0dz7bXX88orL7B79y7efnsm//znq7z66r94/PGHKSoqolmzC3n22cl069ajQnnCjek6zrHvl8KmOUTiw2fa2exvirfTtfymVxucDp1MKyJVx2bqKnXnze83yM4++V5KXq+HrKxD1KrVgIiISAuSVS6n037WG2lee+1Ihg27gltvvfOU6+fP/5znnvs/vv56bVVEDHqnes84nfbS+3WF241MXdu+xrXibRxm8ST4rd6G7Gl2DcMu7kBCbNX/mwrnsbeSxt064TT2yclxZea2no72GImI5fyZe8haPZeYg+twAAd8SeyI7UKXq0ZxUb2Esz5fRKSyqBiJiGVM0+D4qo9g83xiSpYt9Xamdv/rGNlB84hEpPqpGEmFfPTR52dcP3z4SIYPH1lNaSRYmKaJL+17jn47h/i8PQCsczcjt3F/Bl8+gPgYXeJARKyhYiQi1co0fOQtewv7rq+JB3ymnQX2QXS7+koGpdS0Op6IhDkVIxGpFqZp4tuzluxVnxFbcADDtLHc04boDkO4tn9HnW0mIgFBxagK6EQ/OVfh8l4xDT95S6dj372SWMBtOvkyYgiDRg8npU681fFEREqpGFUih8MBgMfjJjIyuK50LdbweIrvo+ZwhOY/RdM08e3/gayVHxGXvw/DtLHU3Z64zkMY268dDrv2EolIYAnNn8YWsdsdxMTEc/x48b2cIiOjgvqsGsOw4feHxx6N6maaJh6Pm+PHc4iJiS9zb7lQYZoGBV+/h7l1CXGAx3TwZcTlDBo9nEbaSyQiAUrFqJIlJCQDlJajYGa32zGM0L7gl9ViYuJL3zOhxJ+9n+zls4g+shmA5a42RHQcyvUDOmgvkYgENBWjSmaz2UhMrEWNGkn4/T6r45Sbw2EjMTGW3NxC7TWqIg6HMyT3FLk2LsCzeibRgN+0Mc8cRO9rr+bChrpQo4gEvpAvRocOHWLy5MmsXr0aj8dDx44defjhh2nRokWVbtdut2O3B+9tQZxOO9HR0RQV+UP+MvFSOUxPETnffkrEtoXYgI2eJhxudDGjh/cnJirkf9SISIgIvf+u/oLH4+GOO+4gMzOT119/nffff5+4uDhuuukmsrOzrY4nEjJ86ZvJffdPRGxbCMBiT2ciL7mLa66+WKVIRIJKSP/EWrt2LTt27GD58uXUq1cPgMmTJ9OzZ0+WLl3Ktddea3FCkeDn3rOBosWv4DB9HPYnsCGmN5eMuZpaidFWRxMROW8hXYxatGjBG2+8UVqKgNI5HXl5eVbFEgkJxvFs8pdNw37oRxzAD54UDradwOhBLTTBWkSCVkgXozp16jBw4MAyy2bMmIHL5aJv374Vem2nM7R/8DtKrkLs0NWIq10wjL1RcIzsz/5CRGEmftPGGl8r6g39Lde3bmB1tAoJhrEPRRp362jsTxbUxSg9PZ3Bgwefdv2qVatITv75VOhFixbxj3/8g5tvvplWrVqVe7t2u42kpLhyPz+YJCTEnP1BUiUCcexN0yBv4zIyFr5HlDeXLH8c8+Ku5nc3D6F+rdD5NxGIYx8ONO7W0dj/zGYG8T0JvF4v+/btO+36Zs2alV6NeubMmTz99NNceeWVPPfccxU6TdrvN8jLKyr384OBw2EnISGGvLwi/H6dlVadAnXsTdMkf+k0/NuXA5Dlj2d9k/GMGtGTSKfD4nSVI1DHPtRp3K0TTmOfkBBzTnvGgnqPUUREBKmpqWd93OTJk5k6dSq//e1veeihhyrlatThcgq732+EzfcaaAJp7E1PETnL3ydi9woM08Z8V1eaDryK0Z2bAKH37yGQxj6caNyto7H/WVAXo3NxohQ99NBD3HLLLVbHEQk6voPbOL7gZSK8xwGY4+9H3+uu0wUbRSQkhXQxWr16NVOnTmXChAmMHDmSzMzM0nWxsbHExYXOnAiRquA7mkb+/H/iNNwc8ddgdfQARt5wJTXjdZNkEQlNIV2M5s6dCxSfiTZjxowy637/+99zzz33WBFLJOCZfh+u7+dRtH4uEaaXn7z12Nh0HOOHtiMiROYTiYicSlBPvraK32+QnV1gdYwq5XTaSUqKIyenQMedq5nVY2+aBscXvgJp6wDY7atLVrfbufii5pUyPy+QWT324Urjbp1wGvvk5LjQn3wtIpXL9Lk59r93caatw2fa+djdl14jruSS1DpWRxMRqRYqRiICgD9zD/lfvoizKAeA2f6BXH79NTSuG29xMhGR6qNiJCIY+Znkzf0HTu9xcvyxLI/sz1VjR5NUQ5OsRSS8qBiJhDlf+mZyFk8j2nucdF8SS2vfwK2juhAdqR8PIhJ+9JNPJIy5N3yOZ83HRAPZ/jg2N7uB3w3tppvAikjYUjESCVPuXWvwrPkYgOWuVkR0G821fVqG/JlnIiJnomIkEmZMn5ui1R/h/XExdmC5qzUJg26kf8eGVkcTEbGcipFIGCm+RtGrkL4RO7DOcyENL7+ZLq3qWx1NRCQgqBiJhAnTNMlbOQt7+ka8pp133Zdw+agRtGqSZHU0EZGAoWIkEgaM41nkL3ode+ZPAMz29WfUmKtoUq+GxclERAKLipFIiDN9bvLnPY899wBu08lSsycjxl1P3ZoxVkcTEQk4KkYiIcz0echZ+AYRuQfIN6J51341t98wSBduFBE5DRUjkRDlz95P/hcvElGQid+0Mdd5OXeOvZjEuEiro4mIBCwVI5EQZLoLyJ//TxyF2RwzYljqvJixN4ygRqxKkYjImagYiYQY01NE1vyXiSrM5qg/no9ixnDX9b2Ii46wOpqISMBTMRIJIf7MveR/+S+iinLwmnYWx43g92N6ExOlf+oiIudCPy1FQoTpOk7+Fy/gcB3jqD+eFXGXM+H6y3QzWBGR86A7RYqEANPrJufLV3C4jnHEX4NP4m/gN9cPUykSETlP+qkpEuR+efaZ17SzOHY4d/3mIpUiEZFy0E9OkSBmeoo4Pv8FHIVZ5PhjWRI9hAljLtOcIhGRctKhNJEgZRp+chZPxV6YRZY/nplRYxkzdgSx0SpFIiLlpZ+gIkHIOJ5F/oJXiMjajd+08WXkEO4e24tYnZIvIlIhKkYiQcb0+zj+xb+w5+yjyIhgoWMQN4wdSnyMSpGISEXpUJpIkDn+7UfYcvZx3IhiGtdwxdhrSdAVrUVEKoX2GIkECdPnpvDr92HHVwDM8ffn1nEDqRmvG8KKiFQWFSORIGCaJoVL/o2Rth6AxZ4ujPjNKOrWjLE4mYhIaNGhNJEg4Nn+NUbaenymnTcLL6PD1TeTUjfe6lgiIiFHxUgkwHl3r6FgxQwAFrg6M+SqoTRvlGhxKhGR0KRDaSIBzLNxPu7VHxIBbPfWp/mQ39D+glpWxxIRCVnaYyQSoHyZe3F99xEAS4vacrzXXXRvU9/iVCIioU3FSCQA+QvzyPriVWymwUZPE/ydRzOoWxOrY4mIhDwdShMJML4je9g7/wUiC3PIN6LZ03gk4wekWh1LRCQsqBiJBBDT6yL/i39hK8zhsD+B5QkjufGKbthsNqujiYiEBR1KEwkg+d98iK0gmyx/HB9GXcfYay/G6dA/UxGR6hJWP3HXrl1LmzZtWL16tdVRRMowDT+Fq2Zh274UgC/tg5g4thcxUdqpKyJSncLmp25+fj4PPvgghmFYHUXkJK7lb+HfsQKApd5O3Pz/riM+0o7Pp/eriEh1Cps9Rk8++SSNGze2OobISXzpm/HtWIFhwjsFA2g36hYa16thdSwRkbAUFsVo9uzZbNiwgUmTJlkdRaQMoyCH3CXTAFjubkPHwUNp3TTJ4lQiIuEr5A+lpaen8+yzz/Lqq68SFxdXaa/rdIZ2p3SUTPh1aOJvlfEd3kX+588T6cknyx+Hv+OVXNw1RWNvIY29NTTu1tHYnyyoi1F6ejqDBw8+7fqVK1fywAMPMGbMGLp37056enqlbNdut5GUVHklK5AlJOju7VXB8BSRtuAl7J58DviSWN/oN9xzTQ/s9p9Py9fYW0djbw2Nu3U09j8L6mJUr1495s+ff9r1H3zwAUVFRdxzzz2Vul3DMMnLK6zU1ww0DoedhIQY8vKK8Ps1Abiy5X01A/N4Nkf98Xwacy0PjuxFbm7xe0pjbx2NvTU07tYJp7FPSIg5pz1jQV2MIiIiSE09/RWBP/nkE44cOULPnj0BME0TgNtvv51Ro0bx1FNPlXvb4XK2kN9vhM33Wh1M08T942K8Py7GDnxh9OV313bHYbedNM4ae+to7K2hcbeOxv5nQV2MzmbGjBn4fL7Srw8fPsyECRN45pln6Nu3r4XJJFx51n2Gd/1s7MA3nlaMuPYKkmpEWR1LRERKhHQxatSoUZmvHQ4HUHwIrlatWlZEkjDmz9qHe/3n2IDZhV1pffkYmtbXafkiIoFE09BFqoFp+MldOh0bBhs8TYntegXd29SzOpaIiPxKSO8x+rWUlBS2b99udQwJM6a7gPyFrxCRsxeX6WR7nSHc2v8Cq2OJiMgphFUxEqlupmlSuOhlbIe24jadfG67lHFX9cRus539ySIiUu10KE2kCvl2rcY4uBWP6eDfrmEMu2YksdH6/4iISKBSMRKpIkbhMfJWvAfA4qIOXDFiIA1qhceFQUVEgpX+6ypSBfzZ6eTPnUKkN59Mfw3iu4+gU/PaVscSEZGz0B4jkUpm+r0ULHgRh+sYh/0JfF37Oob3bW51LBEROQcqRiKVzL1pAbb8I+QaMbznGM11V/bBpsnWIiJBQcVIpBL5s/ZTtG4OAPPd3bn16u7EROmItYhIsNBPbJFK4juwhYIv/4XT8LDLW5cOg4eRUife6lgiInIetMdIpBKYniIKF7+G3e9mh7c+2y+8gd7tG1odS0REzpOKkUglKFo3B5s7nyP+GiyMv5rRl3W0OpKIiJSDipFIBfkyfsLzw0IAFvh6cefVnYhw6p+WiEgw0hwjkQrwpW2gYOErOPDzo6cRA0cMJTkh2upYIiJSTvpvrUg5me4CCpdNxW76+MGTQmbHm2h3YS2rY4mISAWoGImUU9Ha2dg8BRzy1WRl0ihG9GthdSQREakgFSORcvBn/IT3x8UALDB6cvtVHbDbdRFHEZFgpzlGIufJt38TBQtewo7BZk8KFw8fQs34KKtjiYhIJdAeI5HzYHoKKVzyb+yGlx89jchoO472mlckIhIyVIxEzoNr/VxsngIy/IksS7iSkQNbWR1JREQqkYqRyDnyZx/AXXK9ooW+Htx+VUecDv0TEhEJJZpjJHIO/Ed2kz93Cg7Txw5vffoMGULtxBirY4mISCXTf3dFzsL0eShY+BIOXyF7fbXZlXo9XVrVtTqWiIhUARUjkbPw/LgYW2EOOf5YPo0cxdWXdrA6koiIVBEVI5EzMIryKFz7OQAL3F245crORDgdFqcSEZGqojlGIqdhFOSQN+dvOP1FHPLVpFmfIaTUjbc6loiIVCHtMRI5BdM0KPzyBRz5GRwzYllR80oG92hidSwREaliKkYip+Db9R1mVhpFRgRTPVdwzZV9sdt0yw8RkVCnYiTyK6bhI//bjwBY6mrHqGE9dMsPEZEwoWIk8gum4ef4kqlEFB4l34jGaHUJXVrUsTqWiIhUExUjkV9wr3gL9nyL37Sx0D6Qay5tZ3UkERGpRipGIiX8R/fi3b4Cw4S3CgYxcOQIoiJ0ar6ISDhRMRIpcbxkXtE6z4U069GfCxokWJxIRESqm4qRCODdvQb7wc34TRs/xPfmij7NrI4kIiIWUDGSsOfdvYaixa8BsNLTmmtH9sHp0D8NEZFwpJ/+EtZMdwFFX03HhsF37gtx9PgNjWrHWR1LREQsEhbFaNq0aQwePJiOHTsyevRovv32W6sjSYBwb/wSm7eIg76arE0azpCLmlkdSURELBTyxejVV1/l5Zdf5k9/+hNz5syhc+fO/L//9//Yv3+/1dHEYkZRHq5NCwBY5O3KLVe0w27X1a1FRMJZSBejwsJC3nzzTe6//36GDx9Os2bNePTRR2ncuDHr1q2zOp5YyHQXkPf5FByGh3RfEh0GXEKdmjFWxxIREYs5rQ5QldatW0dRUREjRowoXeZwOJgzZ46FqcRqpmlSuPAlHMf2cdyIYnXSFdzUuZHVsUREJACEdDHas2cPiYmJbN++nRdeeIG9e/fSvHlz/vjHP9K1a9cKvbbTGdI723CUnJXlCMGzs7z7NmEc2obHdPCmaxj33tifiAC6kGMoj32g09hbQ+NuHY39yWymaZpWhyiv9PR0Bg8efNr1f/jDH/j3v/9N/fr1+dOf/kTDhg2ZNWsWn332GZ999hmpqanl2q5pmth0p/WgZJome6c+hHlkF8uK2tJ45G0M6dnU6lgiIhIgKm2PkdvtxmazERkZWVkveVb16tVj/vz5p12/ZMkSXC4XkyZNYuDAgQC0a9eODRs28O677/LEE0+Ua7uGYZKXV1iu5wYLh8NOQkIMeXlF+P2G1XEqjWvrcswju/CYDg7U7cvVLWqRk1NgdawyQnXsg4HG3hoad+uE09gnJMSc056xchej1atXs2TJEtavX8+uXbtwuVwAREdHk5qaSpcuXbj00kvp2bNneTdxVhEREWfc6/Pjjz8C0KpVq9JlNpuN1NRU0tPTK7Rtny+030An+P1GyHyv3l3fUbRsGjZgubc91w3vit9vAoG50zSUxj7YaOytoXG3jsb+Z+dVjLxeL7NmzeI///kPBw4cIDExkXbt2jFy5EgSExMxTZO8vDzS09OZM2cOM2bMoGHDhtxyyy2MGTOGiIiIqvo+Tql79+7YbDa+//57hg4dChQfStm5cye9e/eu1ixiLdN1vORCjiarXM1J7HsNtRN1FpqIiJR1XsVoyJAheL1eRo0axbBhw2jXrt0ZH79582a+/PJLXn/9daZPn87SpUsrFPZ8NWzYkGuuuYZnnnmGmJgYmjRpwowZM0hPT+eGG26o1ixiLffmhdh8Lg74ktiQPIwHuja2OpKIiASg8ypGd955J6NHjz7neUTt27enffv2TJw4kU8++aRcASvqySef5OWXX+axxx4jNzeXtm3bMn36dC688EJL8kj1Mz2FFG1ciBNY6unEzSPaYNfkeREROYWgPivNKn6/QXZ2YE3YrWxOp52kpDhycgqC+riz6feRN/8F7Ic2c8iXyN7u9zHkosA+Cy1Uxj4YaeytoXG3TjiNfXJy3DlNvq7QhQsuuugi7rvvPtxud0VeRqTKuL6aiv3QZjymg69jB3Np9yZWRxIRkQBWoWKUl5fH/Pnzue222ygoOHkPyrZt25g5c2ZFNiFSbv7Mvfh2fovftDG94BIuv2Kw7oUmIiJnVOFLXf7ud78jLS2NW265hby8vDLrduzYwVNPPVXRTYiUS+G64lu/rPdcQIsevWlUO87iRCIiEugqXIxSU1N57733yMzM5KabbiI7O7sycolUiP9oGuxbD8CGqO6M6N3M2kAiIhIUKuXmKI0bN+a9996jsLCQG2+8kczMzMp4WZFyMfIyyZ/3PADfe5owclgfIkL83nYiIlI5Ku3TokGDBrz77rsAjB8/nkOHDlXWS4ucM9MwKPzynzjcuRzyJZJ+wShaNq5pdSwREQkSlfrf6Dp16vDuu+8SFxfHDTfcQFpaWmW+vMhZ+fauxTx2kAIjkhn+EVw1uIPVkUREJIhUqBh9+eWXJ91ao2bNmrzzzjvUq1ePV155pULhRM6HaZoUrPkcgOXu1lx1WWdio6v3NjQiIhLcKlSMmjVrRu3atU9aHh8fz/Tp0xk8eDAxMboflVQP767vcOTux206yWnQl26t6lgdSUREgsx53RLkfMTGxmqPkVQb38GtFC17EzvwjbcN1w3thE23/RARkfN0XnuMhg8fzmeffYbH4znn53g8Hj7++GOGDx9+3uFEzoXpdVO46FXspo9NnsbE9LyG5IRoq2OJiEgQOq89RldffTV/+ctfePbZZ7nkkkvo3bs37dq1IyUlpfSQWWFhIenp6WzevJlvvvmGZcuWERERwa233lol34CId/tybO58Mv01WBo7jEnddNsPEREpn/MqRrfffjtjx47lo48+4tNPP2X27NmlhyscDgcAfr8fKJ4I26JFC+655x6uvfZa4uPjKzm6CJiGj4L183ECy1xtGXdVexx2XbNIRETK57znGMXHx3PzzTdz8803k56ezoYNG9i9ezfHjh0Dis9Ku/DCC+ncuTONGzeu7LwiZbg2LcLpyiHfiCaqdX8ubJhgdSQREQliFZp8nZKSQkpKSmVlETkv3l3f4f3uQ2zAcn8nrr64ldWRREQkyFXKMQev18vRo0fPa1K2SEWYruMUfTUNGyYrXK1o3P8q4nTNIhERqaAKFaMtW7Zw66230qVLF/r370+XLl248cYb2bhxY2XlEzklz9Zl2HxuDviS+CH5Mvp0aGB1JBERCQHlPpS2du1abrvtNrxeL126dKFevXocOXKEdevWMX78eF5//XX69u1bmVlFADD9Xgo3LsQJ/M/djvGXt9Y1i0REpFKUqxh5PB4effRRkpKS+Pe//03Lli1L1+3YsYM77riDP//5z8yfP5/oaF1PRipX0fp5OD355Phjqd25P43q6IxHERGpHOU6lLZq1SrS0tJ46qmnypQigJYtW/LMM89w8OBB/ve//1VGRpFS3h0r8W/4DICv6c7Ivs2tDSQiIiHlvPcY/fjjjyxatIiYmBiSkpL48ccfT3pMzZo1iYuLY9myZTRp0oS2bdtWSlgJb4Yrn6IVb2EDlha1pfWwkURFOqyOJSIiIeS8i9E111xTOp/juuuuO+NjZ8+ezZw5c9i6dWv50on8gnfLMmx+L/t9yeysdxlXtjj5BsYiIiIVcd7F6LnnnmPBggV89913PP7446d8jGmaPPPMM3Tt2pURI0ZUOKSI6fdSuGkRTuArdztuuKylJlyLiEilO+9iNHr0aGJjY1m+fDlNmzala9euJz1mw4YNFBYWcvXVV6sYSaVwbVpYOuG6Vse+NKgVZ3UkEREJQeWafN2/f3/q1KnDI488QlpaWpl1aWlpPPTQQ9SuXZuBAwdWSkgJb949a/Gu+QiAr43OjOynCdciIlI1ynW6flxcHE8//TR33303I0aMoFu3bqXXMVq7di2mafLSSy/pxrFSYaa7gKJlb2LD5GtXS5oOuoKYqArdyUZEROS0yv0JM3DgQN566y1eeOEF1q5di9/vx26307lzZ+69914uuuiiyswpYcq77StsPjcHfTXZUPNSHtYVrkVEpApV6L/e3bt3591338XtdpObm0tCQoIu6CiVxjT8FG5ciAP4ytWGG0a3xq4J1yIiUoUq5ZhEVFQUdevWrYyXEinl2bESh+sY+UY00a370qx+gtWRREQkxFXoJrIiVcWf8ROuFe8A8I2vLaMGtTzLM0RERCpOxUgCjunzULjoZeymjx88KST1GkVCbKTVsUREJAyoGEnA8e5cBUW5ZPnjWBI9lEFdU6yOJCIiYULFSAKKaZoUfr8AgBXu1vxmSDscdr1NRUSkelRo8nXr1q3PeluGqKgo6tWrR8+ePbntttto0qRJRTYpIc53YAuOvIO4TSeuxr1o1STJ6kgiIhJGKlSM7r77bpYsWcLOnTvp378/TZs2BWDv3r18/fXXtGzZkp49e7Jv3z4++eQT5s2bx3vvvUfr1q0rJfy5KCgoYMqUKSxatAiXy0WXLl14+OGHSU1NrbYMcm6MghyOL3kDJ7DWk8pVg9tbHUlERMJMhYpR3bp1ycnJ4YsvvqBx48Zl1qWlpTFhwgSaN2/OQw89xN69exkzZgzPP/88b7zxRoVCn4+nn36aTZs28eKLL5KYmMjkyZO57bbb+PLLL4mKiqq2HHJmpmlQuOBFnO5cMvyJFLYZSd2aMVbHEhGRMFOhyRvTpk1j3LhxJ5UigKZNmzJu3LjSEtSsWTOuv/56NmzYUJFNnrfFixczduxYunbtSmpqKvfeey8HDx5k586d1ZpDzsy//wfMo3soMiJ413s5Q/tV315FERGREypUjDIyMnA4HKdd73A4OHToUOnXKSkpeDyeimzyvNWqVYv58+eTlZWFx+Pho48+ombNmprrFGCKNi0C4Ft3Cwb160hstO6HJiIi1a9Cnz7Nmzdn5syZXHXVVdSuXbvMuszMTGbOnEnz5j/fCX3//v0nPa6qPfvsszz44IP06dMHh8NBbGws06dPp0aNGhV6XacztM+UcjjsZX6vSv5jGXBwM4YJP8V25v7uKWF9Jlp1jr2UpbG3hsbdOhr7k9lM0zTL++TVq1dz++2343A4uPTSS0snX6elpbF48WJ8Ph9Tp06lZ8+euN1uBg8ezIABA3juuecqJXx6ejqDBw8+7fpVq1bxxRdfsHDhQu68805iY2N58803+eGHH5g1axYNGpTvhqSmaZ71bDw5N6bhJ23G0xjpP/CjpxENxz5Kt9b1rI4lIiJhqkLFCGDLli28+OKLfPvtt7hcLqD4FP3evXtzzz330K5du0oJeiper5d9+/addv2xY8cYN24cS5cupWHDhqXPGTZsGBdffDGPPvpoubbr9xvk5RWV67nBwuGwk5AQQ15eEX6/UWXbKfz6XdybFuIxHcxLuJ5bJlxeZdsKFtU19nIyjb01NO7WCaexT0iIOac9YxWeyNG2bVtef/11DMMgKysLKJ7XY6+GQyERERFnPO1+2rRp1KpVq7QUnXhO27ZtSUtLq9C2fb7QfgOd4PcbVfa9+nMO4t60EID3CvpzzTV9wmZcz0VVjr2cmcbeGhp362jsf1Zp7aWoqIi8vDzy8vIoKgqMvSn169cnJyeHI0eOlC4zDIOdO3fSrFkz64IJAJ4flwAU3w+tXW8a1Y6zOJGIiIS7Cu8x2rRpE5MnT2b9+vUYRnHbtNvtdOvWjQceeIAOHTpUOGR5XXzxxTRu3JiJEyfyyCOPEB8fz/Tp0zl06BA33nijZbkETK8L9/avcQDfGW25pf8FVkcSERGpWDHauHEjEyZMICIigmuvvbb0sNauXbuYN28e48ePZ8aMGXTs2LFSwp6v2NhY3nnnHf7+979z991343a76dChAzNnziQlRTcmtVLRD0tw+N0c8deg1UW9SYiNtDqSiIhIxSZf33zzzRw4cID333+fOnXqlFl39OhRxo4dS0pKCv/5z38qHDSQ+P0G2dkFVseoUk6nnaSkOHJyCir9uLPv4FYK5k7GjsE8ox/X3HYLESF++YPzUZVjL2emsbeGxt064TT2yclx5zT5ukKfRhs3bmTMmDEnlSKA2rVr85vf/Ibvv/++IpuQEGN6XRQuehU7BuvczUgdOFKlSEREAkaFPpHsdjt+v/+06w3DqJaz0yR4eHd+i82dT6a/Bt/UGEr31nWtjiQiIlKqQq2lS5cuvPfeexw4cOCkdQcPHuT999+na9euFdmEhBDTNCn6YTEAX7tbcc0lrXShTBERCSgVmnx93333MW7cOIYNG8Zll11Wegr8nj17WLJkCQ6Hgz/96U+VkVNCgHFkF/Zj6XhNOwUNutOqSZLVkURERMqoUDFq27YtH374IS+88AJLly4tvX5RTEwM/fv359577y1zrzQJX6bPQ+5XM4gANniaMfKS9lZHEhEROUmFr2PUokULXnnlFQzDIDs7G4Dk5GTNLZIyXCveJuJYGgVGJJlNLmNwnXirI4mIiJzkvIrRwYMHz+lxGRkZZb7+5S05JPz4j6bh+2klftPGjKJB3HpxF6sjiYiInNJ5FaNLLrmkXJNlt27det7PkdDh2bIMgI2epjTrchHJCdEWJxIRETm18ypGzz33nM4ikvNieopw7/gGB7DObMPvejW1OpKIiMhpnVcxGj16dFXlkBBVtHUFDsPDYX8CbS/qRWx0hNWRRERETkszpKXK+I+m4fnuQwDW047B3XV/OhERCWwqRlIlTJ+HwoUv4TB9bPE0JKXfFUQ4HVbHEhEROSMVI6kSvr3r4fhRco0YlsQMpVd7nZkoIiKBT8VIqkTh5qUAfONqyZUXt8OuSfsiIhIEVIyk0hnHMrAf2YFh2sis3ZV2FyRbHUlEROScqBhJpTu2+lMAtnobMvTizrrEg4iIBA0VI6lUnm1fEZG2GsOEtFp9SG2UaHUkERGRc6ZiJJXGKDyG6+sZAMwv6kyfwQMtTiQiInJ+VIyk0nh3rMRm+Ejz1SK32WAa19WNYkVEJLioGEmlME2Tws3F90T7xt2KUQNSLU4kIiJy/lSMpFL4D23DWXgUl+kkukVP6ifHWh1JRETkvKkYSYWZho+cr4tv/bHBewEj+reyOJGIiEj5qBhJhblX/5foY3twmREUpg6mVmK01ZFERETKRcVIKsR/7CDeHxYAMMvVj4sHdLU4kYiISPmpGEmFeLeuAGCzJ4UGnfuRGBdpcSIREZHyUzGScjMNH4XbiovReqMVQ3s2sTiRiIhIxagYSbl59m4kwnucfCOaZt36EBcdYXUkERGRClExknIxfW5yVxafibbJaM7gHk0tTiQiIlJxKkZSLkUrZhBTdJg8I5rozsOJjnRaHUlERKTCVIzkvPmPpuH/6WsM08bH/kvoe5GuWyQiIqFBxUjOm3vbcgC+9zSlbc/eRDgdFicSERGpHCpGcl5Mvxf39m8A+NHRmoGdG1qcSEREpPKoGMl5ce/ZgNNfxDEjhlYX9dHeIhERCSkqRnLOTJ+b3G8+BmCz2ZIBnVMsTiQiIlK5QqYYPf744zz88MMnLV+1ahWjR4+mU6dODB06lHnz5lmQLjQUff0esa7iM9Hiug4jwhkybx8REREgBIqRYRg8//zzzJo166R1u3bt4s4776R///588sknXHfddTz44IOsWrXKgqTBzX80Df+O5RgmfGZeQu9uLa2OJCIiUumC+uIzu3bt4tFHHyUtLY2GDU+eBPz222/TqlUr/vjHPwKQmprKli1bmDp1Kr17967uuEHNvf1rADZ6mtK2X1/tLRIRkZAU1J9u3377LampqcydO5eUlJPnu6xdu/akAtSrVy/WrVuHaZrVFTPomYaPopIz0bY6W9O/YwOLE4mIiFSNoN5jNG7cuDOuz8jIoH79+mWW1a1bl6KiInJyckhOTi73tp0hvsfE4bCX/l64fzORvgLyjWja9OpHdFRQv20C3i/HXqqXxt4aGnfraOxPFrCfcOnp6QwePPi061etWnXWYuNyuYiMjCyz7MTXHo+n3NnsdhtJSXHlfn4wqREfSdqqT4kDtthacO0lrXUYrZokJMRYHSFsaeytoXG3jsb+ZwFbjOrVq8f8+fNPuz4xMfGsrxEVFXVSATrxdUxM+d8EhmGSl1dY7ucHA4fDTkJCDAcWvkdcXhouM4IaXYZxPL/I6mgh78TY5+UV4fcbVscJKxp7a2jcrRNOY5+QEHNOe8YCthhFRESQmppaoddo0KABR44cKbPsyJEjxMbGUqNGjQq9ts8X2m8gAG/2Qdzr5mAD5pv9GdetdVh834HC7zc03hbR2FtD424djf3PQvqYSPfu3fnuu+/KLPv222/p2rUrdntIf+uVIuf7/2HDZKu3ARf0vhSnjkGLiEiIC+lPugkTJrBp0yamTJnCrl27mD59Ol9++SW33Xab1dECnmmaZG9YBsBWRxv6tK9/lmeIiIgEv5AuRi1atODVV1/lq6++YtSoUfz3v/9l8uTJuobROXAf2EGkKxuX6aRZj/7aWyQiImEhYOcYna8ZM2accvmAAQMYMGBANacJfgdWziMZ2GZeQJ9OTa2OIyIiUi20G0BO4t61huSs7zFMiGx7iU7PFxGRsKFPPCnDdB2n8KvpAKw0OtGjfy+LE4mIiFQfFSMpw7PrO5y+IjL8iST0vY7ICIfVkURERKqNipGUceyH5QCsN1pxed/mFqcRERGpXipGUsqXl0ls3l4ME5I69CU2OsLqSCIiItVKxUhK7V+9FIA9Rn0G9mlvcRoREZHqp2IkAPjzs4jfU1yMXA27E6e9RSIiEoZUjATTNMla8DrRuNnnr03bwcOtjiQiImIJFSPBd2ALMdk/4TEd7G52LQnxsVZHEhERsYSKkZC1cQUA6zyp9OvXxeI0IiIi1lExCnOm34fjwAYA/E26k1QjyuJEIiIi1lExCnP7N35HNG7yjBguGtDP6jgiIiKWUjEKc3kbFwNwOKEttZI0t0hERMKbilEYy1j/FU29u/GbNlJ6D7U6joiIiOVUjMKU6TqOc91MADbH9aLuBS0sTiQiImI9FaMwlfPjN0SZLjL8iaRccp3VcURERAKCilGYOvbjNwCkxbanWcNki9OIiIgEBhWjMJSXdZRaRWkANOkxyNowIiIiAUTFKAxt/XoZdpvJYVsdWrRpbnUcERGRgKFiFGYKC10kHSo+jGZv1h2bzWZxIhERkcChYhRm9iyaRQN7DoVE07TfMKvjiIiIBBQVozDizjlMw4yvADjS/CocMQkWJxIREQksKkZhZO+3S3HaDPYaDWjVf4jVcURERAKOilGYMAwT9q8HwNekB5ERDosTiYiIBB4VozCxadMOGpKJYUKrPoOsjiMiIhKQVIzCgGma7F+3AoDc2CbE1tQFHUVERE5FxSgMbN2TSWvPDwDUbNfb4jQiIiKBS8UoDGSs+IT6jlzc9lhqtBtgdRwREZGApWIU4vbvTqOL6zsA7Bddjy0qzuJEIiIigUvFKMTtXl18in5mRCOSOvS3Oo6IiEhAUzEKYUdzi6iZvRWA2Na9dfsPERGRs1AxCmHLV2+nmTMTgLodNOlaRETkbFSMQlShy0vOtrXYbSbuGo2wx9eyOpKIiEjAUzEKUV9tPEgnxy4AarTobnEaERGR4BBSxejxxx/n4YcfPmn5xx9/zMiRI+ncuTNDhgzhjTfewO/3W5Cwevj8BvvWraRVxCFMm52IFn2sjiQiIhIUQqIYGYbB888/z6xZs05aN2fOHJ544gnGjx/PnDlzuPfee/n3v//Na6+9ZkHS6rHmx3SG2L4BwNn+cuyJ9SxOJCIiEhycVgeoqF27dvHoo4+SlpZGw4YNT1o/c+ZMRo0axZgxYwBo0qQJe/bs4b///S+///3vqztulTNNk92rV9DecRy3swa1ul9ldSQREZGgEfR7jL799ltSU1OZO3cuKSkpJ62///77ufXWW8sss9vt5ObmVlfEarVlbw71XcVzi6Ja9sEWEW1xIhERkeAR9HuMxo0bd8b13bp1K/N1fn4+M2fOpH//il3s0OkMzE656Ls0xkQcACC+Rbdy53Q47GV+l+qjsbeOxt4aGnfraOxPFtDFKD09ncGDB592/apVq0hOPvc7xRcUFHDXXXfhdrt58MEHy53LbreRlBR4t9bYczCX3H07iE90Q2Qsddp0wuao2F9xQkJMJaWT86Wxt47G3hoad+to7H8W0MWoXr16zJ8//7TrExMTz/m1MjMzufPOO0lPT2fatGmnPOx2rgzDJC+vsNzPryofLtpOu8h0ACIat+dYnhtwl+u1HA47CQkx5OUV4fcblZhSzkZjbx2NvTU07tYJp7FPSIg5pz1jAV2MIiIiSE1NrfDr7Nq1i9tuuw3DMHjvvfdo0aJFhV/T5wusN1BOvps1mw/wUI3dADgad6qUjH6/EXDfa7jQ2FtHY28Njbt1NPY/C/mDivv37+emm24iJiaGDz74oFJKUSBavG4//SN/pJajAFtcMs4LdVFHERGR8xXQe4wqw6RJk/B4PDz//PM4nU4yMzNL19WpU8fCZJWnyO1j7YZd3B+7GYCoi67F5oyyOJWIiEjwCelidPjwYb777jsArrrq5Ov5bN++vbojVYmvNx2iOXuItnmxJzfG2byX1ZFERESCUkgVoxkzZpT5ul69eiFTfk7HMEwWrd3PyIiDADgv7I7NFvJHSEVERKqEPkGD3IafMsnOLaRVRAYAzpT2FicSEREJXipGQW7Rmv00dmQRY/NAZCz22hdYHUlERCRoqRgFsb0ZeexIz6VN5CEAnI3aYrPrr1RERKS89CkaxBatSQdMetUovgWIQ4fRREREKkTFKEgdO+7mu62H6RixjyRfJjijcDbranUsERGRoKZiFKSWrT+AYfgZlfADAJEdL8cek2BxKhERkeCmYhSEPF4/yzYcoE3EQWqZ2RAVR2THoVbHEhERCXoqRkHo2y2HOV7kpWPcUQAiLuiOLTLW4lQiIiLBT8UoyJhm8QUdATqUFCNHw9ZWRhIREQkZKkZBZmtaDgcyC0iM8BPvKj5N39FAxUhERKQyqBgFmUVrivcWDU/1gmliS6iLPS7J4lQiIiKhQcUoiBzOLmTjriwAutXMBsCpvUUiIiKVRsUoiJyYW9TpwmSiMrcB4GjQyspIIiIiIUXFKEgUuLx8/UPxnKKRFxZi5KQXX9SxSSeLk4mIiIQOFaMgsWLjITxeg0Z14qh/aDkAEa0HYouOtziZiIhI6FAxCgJ+w2DJuuLDaFe2ceA/tA1sDiI7Xm5xMhERkdCiYhQENuw4Slaem/iYCNrHHAbA2aQj9vhaFicTEREJLSpGQWDJunQABnVpCEf3AGCv18LKSCIiIiFJxSjA7T9ynO37j2G32RjUuRH+I7sBcNS90OJkIiIioUfFKMCd2FvUtWVtajqKMAuywWbDUaeZtcFERERCkIpRACtwefn2xwwABndLwX9kFwD2pBRsEdFWRhMREQlJKkYBbMXGQ3h8Bil14mnZuCaGDqOJiIhUKRWjAGUYJss2FB9GG9ytETabDf/hnQDYVYxERESqhIpRgNq0O4vMYy7iop30alcfI/8o/oyfAHA2bGtxOhERkdCkYhSgTky67t+xIVERDrzbVwAmjkZtsSfUsTaciIhIiFIxCkCHsgr4cU82NmBQ10aYhlFSjCCi1QBrw4mIiIQwFaMAtHT9AQA6Na9N3Zox+A9tKz5NPyoO5wXdLE4nIiISulSMAkyR28fKHw4BcEm3RgD4M4vPRnOmdMDmiLAsm4iISKhTMQow32zOwOXxUz85lrbNkgEwsov3INmTU6yMJiIiEvJUjAKIaZosXX/iFP0U7DYbAEZ28TJHciPLsomIiIQDFaMAsiUth0NZhURFOujTvj4ApuHDOFZ8aM2epD1GIiIiVUnFKIAsWVu8Z6hf+wbERDkBMHKPgOGDiGhsNWpZGU9ERCTkqRgFiMxjRWzceRT4edI1/HwYzZ7UEJtNf10iIiJVSZ+0AWLZhgOYQLtmSTSoFVe63MgpmV+kw2giIiJVLmSK0eOPP87DDz982vWmaXLrrbcyYcKEakx17rbvywFgcLfGZZbrjDQREZHqE/TFyDAMnn/+eWbNmnXGx7399tt8/fXX1ZTq/N1wWUtuHdGGTs1/nkdkmib+zD2AipGIiEh1cFodoCJ27drFo48+SlpaGg0bNjzt47Zv384rr7xC586dqy/ceUptmEhqw8Qyy4zs9OIrXjsicdRrblEyERGR8BHUe4y+/fZbUlNTmTt3Likpp96j4na7uf/++5k4cSIXXHBBNSesGN++jQA4GrXB5oy0OI2IiEjoC+o9RuPGjTvrYyZPnkzdunUZP348jzzySKVt2+ms+k5ZuL+4GEVd0KVatvdLDoe9zO9SfTT21tHYW0Pjbh2N/ckCthilp6czePDg065ftWoVycnJZ3yN5cuX8/nnnzNnzhxsJVeRrgx2u42kpLizP7AC/IV55GTsBKB2xz44E6p2e6eTkBBjyXZFY28ljb01NO7W0dj/LGCLUb169Zg/f/5p1ycmJp52HUB2djaTJk3iySefpF69epWazTBM8vIKK/U1f82zcx1g4qjVmHx/NOQUVOn2fs3hsJOQEENeXhF+v1Gt2w53GnvraOytoXG3TjiNfUJCzDntGQvYYhQREUFqamq5n//VV1+RmZnJpEmTmDRpEgAejwfDMOjSpQvz5s0744Tts/H5qvYN5D12BABbUkqVb+tM/H7D0u2HM429dTT21tC4W0dj/7OALUYVddlll9G1a9cyy6ZMmUJGRgZTpkyhbt26FiU7N+bxLADs8boNiIiISHUJ2WIUHx9PfHx8mWVxcXFER0fTtGlTi1KdO6OkGNlUjERERKqNpqEHKLNAe4xERESqW8jsMZoxY8ZZH/PXv/61GpJUDuN4NqA9RiIiItVJe4wCkOkpAnfxWWj2+DNfkkBEREQqj4pRADIKivcWERmLLVLXlhAREakuKkYBSGekiYiIWEPFKAAZ+SfOSNNhNBERkeqkYhSAtMdIRETEGipGAUjXMBIREbGGilEAMksmX2uPkYiISPVSMQpApt8LgL1mA4uTiIiIhJeQucBjKIkecAtGzgEctQP/1iUiIiKhRMUoADmSU3Akp1gdQ0REJOzoUJqIiIhICRUjERERkRIqRiIiIiIlVIxERERESqgYiYiIiJRQMRIREREpoWIkIiIiUkLFSERERKSEipGIiIhICRUjERERkRIqRiIiIiIlVIxERERESqgYiYiIiJSwmaZpWh0i2JimiWGE/rA5HHb8fsPqGGFJY28djb01NO7WCZext9tt2Gy2sz5OxUhERESkhA6liYiIiJRQMRIREREpoWIkIiIiUkLFSERERKSEipGIiIhICRUjERERkRIqRiIiIiIlVIxERERESqgYiYiIiJRQMRIREREpoWIkIiIiUkLFSERERKSEipGIiIhICRUjOWdr166lTZs2rF692uooYeHQoUPcd9999O3blx49enDrrbfy008/WR0rJBmGwYsvvkj//v3p3Lkzt99+O/v377c6Vlg4duwYjz/+OAMGDKBr166MHTuWtWvXWh0rrOzZs4cuXbrwySefWB0lIKgYyTnJz8/nwQcfxDAMq6OEBY/Hwx133EFmZiavv/4677//PnFxcdx0001kZ2dbHS/kvPrqq7z//vs8/fTTfPDBBxiGwW233YbH47E6Wsi777772LBhA88//zwff/wxbdq04dZbb2X37t1WRwsLXq+X+++/n8LCQqujBAwVIzknTz75JI0bN7Y6RthYu3YtO3bsYMqUKXTo0IEWLVowefJkCgsLWbp0qdXxQorH42H69OlMnDiRQYMG0bp1a/75z3+SkZHBwoULrY4X0tLS0li5ciVPPvkk3bt354ILLuDPf/4zdevW5fPPP7c6Xlh46aWXiI+PtzpGQFExkrOaPXs2GzZsYNKkSVZHCRstWrTgjTfeoF69eqXL7Pbif655eXlWxQpJ27Zto6CggN69e5cuS0hIoG3btqxZs8bCZKEvKSmJN954gw4dOpQus9ls2Gw2vc+rwZo1a5g1axZ//etfrY4SUFSM5IzS09N59tln+fvf/05cXJzVccJGnTp1GDhwYJllM2bMwOVy0bdvX4tShaaMjAwAGjRoUGZ53bp1S9dJ1UhISGDgwIFERkaWLluwYAFpaWn079/fwmShLy8vjwcffJDHHnvspPd+uHNaHUCsk56ezuDBg0+7fuXKlTzwwAOMGTOG7t27k56eXo3pQtvZxn7VqlUkJyeXfr1o0SL+8Y9/cPPNN9OqVavqiBg2ioqKAMp8OANERUWRm5trRaSwtX79eh555BGGDBnCoEGDrI4T0p588km6dOnCyJEjrY4ScFSMwli9evWYP3/+add/8MEHFBUVcc8991RjqvBwtrFPTEws/fPMmTN5+umnufLKK3nwwQerI15YiY6OBornGp34M4Db7SYmJsaqWGFn8eLF3H///XTt2pUpU6ZYHSekffbZZ6xdu1bzuE7DZpqmaXUICUyXXHIJR44cISIiAgDTNCkqKiIqKopRo0bx1FNPWZww9E2ePJmpU6fy29/+loceegibzWZ1pJCzadMmrrvuOhYtWkSTJk1Kl48dO5ZWrVrx5JNPWhcuTLz77rs8++yzDB06lL/97W8n7b2TyjVhwgTWr19fZpwLCwuJjIykZ8+eTJ061cJ01tMeIzmtGTNm4PP5Sr8+fPgwEyZM4JlnntE8l2pwohQ99NBD3HLLLVbHCVmtW7cmPj6e1atXlxajvLw8tmzZwvjx4y1OF/pOXCZhwoQJPProoyr/1WDKlCm4XK4yy4YMGcLEiRO58sorLUoVOFSM5LQaNWpU5muHwwEUHwaqVauWFZHCxurVq5k6dSoTJkxg5MiRZGZmlq6LjY3VRPhKFBkZyfjx45kyZQrJyck0atSIyZMnU79+fYYMGWJ1vJC2Z88ennvuOS677DLuvPNOjh49WrouOjqaGjVqWJgudP3ybNdfqlWr1mnXhRMVI5EANHfuXKB4r92MGTPKrPv973+veV+VbOLEifh8Ph577DFcLhc9evRg2rRppYeRpWosWLAAr9fLokWLWLRoUZl1V199tU4jF0tojpGIiIhICV3HSERERKSEipGIiIhICRUjERERkRIqRiIiIiIlVIxERERESqgYiYiIiJRQMRIREREpoWIkIiIiUkLFSERERKSEbgkiImFtz549DB069JTratSowdq1a6s5kYhYScVIRMJafHw8f//738ss2717N6+//jp9+/a1KJWIWEX3ShMR+YX9+/dzww030LhxY6ZNm0ZMTIzVkUSkGqkYiYiUOHToEOPGjSM5OZm33nqL+Ph4qyOJSDVTMRIRATIzMxk/fjxRUVG888471KxZ0+pIImIBFSMRCXs5OTlMmDABv9/Pu+++S61atayOJCIW0eRrEQlr+fn53HrrrbhcLt577z2VIpEwp2IkImGroKCA22+/naysLN577z3q1atndSQRsZgOpYlI2Lr77rtZvHgxd911F82aNSuzrmXLlrRp08aaYCJiGRUjEQlLpmnStWtXCgsLT7n+8ccfZ9y4cdWcSkSspmIkIiIiUkL3ShMREREpoWIkIiIiUkLFSERERKSEipGIiIhICRUjERERkRIqRiIiIiIlVIxERERESqgYiYiIiJRQMRIREREpoWIkIiIiUkLFSERERKSEipGIiIhIif8Pr3D7j3vooa8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(z, norm.logcdf(z), label=\"Actual\");\n",
    "ax.plot(z[2:], cum_logΦ_approx, label=\"Approximate\");\n",
    "\n",
    "ax.set_xlabel(\"$z$\");\n",
    "ax.set_ylabel(r\"$\\log \\Phi(z)$\");\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd5934ef-e3a9-4b07-b769-0e4007da5c70",
   "metadata": {},
   "source": [
    "This post is available as a Jupyter notebook [here]()."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6c047641-625e-492f-a470-9c34d79fcd45",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Last updated: Fri Jun 06 2025\n",
      "\n",
      "Python implementation: CPython\n",
      "Python version       : 3.12.5\n",
      "IPython version      : 8.29.0\n",
      "\n",
      "matplotlib: 3.9.2\n",
      "numpy     : 1.26.4\n",
      "scipy     : 1.14.1\n",
      "seaborn   : 0.13.2\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -n -u -v -iv"
   ]
  }
 ],
 "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.12.5"
  },
  "nikola": {
   "date": "2025-06-06",
   "title": "Simpson's Rule on the Log Scale"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}

Is this kind of machine learning related alg?