From f514fa14d9f352a3aad72f08c07a400d25e1a77e Mon Sep 17 00:00:00 2001 From: Andrew Cumming Date: Fri, 29 Sep 2023 06:17:13 -0400 Subject: [PATCH] Adds some notes to monte carlo integration solutions --- montecarlo_solutions.ipynb | 75 +++++++++++++++++++++++++++----------- 1 file changed, 54 insertions(+), 21 deletions(-) diff --git a/montecarlo_solutions.ipynb b/montecarlo_solutions.ipynb index fdbbfa9..bebcd72 100644 --- a/montecarlo_solutions.ipynb +++ b/montecarlo_solutions.ipynb @@ -8,6 +8,14 @@ "# Monte Carlo Integration" ] }, + { + "cell_type": "markdown", + "id": "f71e7eb4", + "metadata": {}, + "source": [ + "Integrate $\\exp(-|x^3|)$ from $0$ to $\\infty$." + ] + }, { "cell_type": "code", "execution_count": 1, @@ -25,7 +33,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 22, "id": "5582478b", "metadata": {}, "outputs": [], @@ -33,19 +41,23 @@ "def integrate_MC(f, N, sampling = 'uniform'):\n", "\n", " if sampling == 'gaussian':\n", - " p = lambda x: np.exp(-x**2/2)/(np.sqrt(2*np.pi))\n", - " x = rng.normal(size = N) # Generate x values between +-xmax\n", + " # We need to normalize the Gaussian from 0 to infinity\n", + " p = lambda x: np.exp(-x**2/2) * np.sqrt(2/np.pi)\n", + " x = rng.normal(size = N) \n", " else:\n", - " xmax = 3\n", - " p = lambda x: np.ones_like(x) / (2*xmax)\n", - " x = xmax*(-1 + 2*rng.uniform(size = N))\n", + " # Generate x values between 0 and xmax\n", + " xmax = 10\n", + " p = lambda x: np.ones_like(x) / xmax\n", + " x = xmax*rng.uniform(size = N)\n", " \n", - " return np.sum(f(x)/p(x)) / N" + " # use np.mean to calculate the integral as an alternative to np.sum()/N\n", + " # also return the error in the mean, which is an error estimate for the integral\n", + " return np.mean(f(x)/p(x)), np.std(f(x)/p(x))/np.sqrt(N)" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 28, "id": "5c39e032", "metadata": { "scrolled": false @@ -55,14 +67,14 @@ "name": "stdout", "output_type": "stream", "text": [ - "Integral from scipy.quad = 1.78596 with error 5.49113e-09\n", - "Uniform: I = 1.78656 +- 0.0220357; frac err = 0.0123341\n", - "Gaussian: I = 1.78585 +- 0.00960208; frac err = 0.00537675\n" + "Integral from scipy.quad = 0.89298 with error 2.74557e-09 (135 function evaluations)\n", + "Uniform: I = 0.892806 +- 0.02417; frac err = 0.0270719; error estimate = 0.0248929\n", + "Gaussian: I = 0.893114 +- 0.00462334; frac err = 0.00517665; error estimate = 0.00463627\n" ] }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAus0lEQVR4nO3de3hU1b3/8c/cMplwSSA0XIMHpIgVEQTBQJUqUapUpWCLVhCVB1sbrYitlFq11QpeWqXWgJZjUY9QECpWgeJBquCFmyj94aVQqpYUJCglBHKZzGX9/kid0wgoQ/Zkr2S/X88zjys7O2u+yzCTz6y9Zo3PGGMEAABgMb/bBQAAAHwRAgsAALAegQUAAFiPwAIAAKxHYAEAANYjsAAAAOsRWAAAgPUILAAAwHpBtwv4rGQyqd27d6tNmzby+XxulwMAAI6BMUYHDx5Uly5d5Pc7Px9iXWDZvXu3CgsL3S4DAAAch7KyMnXr1s3xfq0LLG3atJFUP+C2bdu6XA0AADgWlZWVKiwsTP0dd5p1geXTy0Bt27YlsACQJMUSMT313lOSpPEnj1coEHK5IgBHk6nlHNYFFgD4rFgypgc2PyBJGnfSOAIL4EEEFgDWC/qDuvjEi1NtAN7DIx+A9bICWbr7q3e7XQYAF7EPCwAAsB6BBQAAWI/AAsB61bFqDV0wVEMXDFV1rNrtcgC4gDUsAJqFg7GDbpcAwEUEFgDWyw5ma9k3l6XaALyHwALAen6fXye0PcHtMgC4iDUsAADAesywALBeLBnTku1LJEmX9r5UIT873QJeQ2ABYL1YIqYZG2ZIki458RICC+BBBBYA1gv4AzrvhPNSbQDeQ2ABYJeKMql6X4NDYUkP9L6y/ovy99LvMydfyitsfG0AXENgAWCPijKpdLDk9OZwoRypZCOhBWjGCCwA7FG9rz6sjJkrdejtTJ+fbJeemVzfN4EFaLYILADs06G31KV/6suaeI2+8cw3JEnLxixTJBhxqTAAbiGwALCeMUZ7a/am2gC8h8ACwHrhQFiLL1qcagPwHgILAOsF/AH1ad/H7TIAuIit+QEAgPWYYQFgvVgypuXvL5ckjeo5ip1uAQ8isACwXiwR022v3SZJOv+E8wksgAcRWABYL+AP6KyuZ6XaALyHwALAeuFAWLOLZ7tdBgAXsegWAABYj8ACAACsR2ABYL2aeI1GPTNKo54ZpZp4jdvlAHABa1gAWM8Yo50Hd6baALyHwALAeuFAWE9e8GSqDcB7CCwArBfwBzSgYIDbZQBwEWtYAACA9ZhhAWC9eDKu1TtXS5JGdB+hoJ+nLsBreNQDsF5dok4/XPNDSdKG72wgsAAe1KhLQvfcc498Pp+mTJmSOlZbW6uSkhLl5+erdevWGjt2rMrLyxtbJwAP8/v8GtRxkAZ1HCS/jyvZgBcd98uUTZs26dFHH1W/fv0aHL/pppu0fPlyLV68WLm5ubr++us1ZswYvfbaa40uFoA3ZQezNe/r89wuA4CLjuulyqFDh3TFFVdo7ty5ateuXer4gQMH9Nhjj+mBBx7Queeeq4EDB2revHl6/fXXtX79eseKBgAA3nJcgaWkpESjRo1ScXFxg+ObN29WLBZrcLxPnz7q3r271q1bd8S+otGoKisrG9wAAAD+U9qXhBYuXKg333xTmzZtOux7e/bsUVZWlvLy8hoc79ixo/bs2XPE/mbOnKmf//zn6ZYBwENq47Uav2K8JOmpC59SdjDb5YoANLW0ZljKysp04403av78+crOduYJY/r06Tpw4EDqVlZW5ki/AFqOpElq2/5t2rZ/m5Im6XY5AFyQ1gzL5s2btXfvXp1++umpY4lEQmvXrtXDDz+sF154QXV1daqoqGgwy1JeXq5OnTodsc9wOKxwmK22ARxdOBDWo+c9mmoD8J60AsuIESO0devWBseuvvpq9enTR9OmTVNhYaFCoZBWr16tsWPHSpK2bdumnTt3qqioyLmqAXhKwB/Q0C5D3S4DgIvSCixt2rRR3759Gxxr1aqV8vPzU8cnTZqkqVOnqn379mrbtq1uuOEGFRUV6cwzz3SuagAA4CmObxf54IMPyu/3a+zYsYpGoxo5cqRmz57t9N0A8JB4Mq7Xd78uSRraZSg73QIe1OhH/csvv9zg6+zsbJWWlqq0tLSxXQOApPqt+UtWl0hia37Aq3jUA7Ce3+fXKfmnpNoAvIfAAsB62cFsLfzGQrfLAOAiXqoAAADrEVgAAID1CCwArFcbr9WEFRM0YcUE1cZr3S4HgAtYwwLAekmT1JaPt6TaALyHwALAelmBLM06Z1aqDcB7CCwArBf0BzWi+wi3ywDgItawAAAA6zHDAsB6iWRCb+59U5J0esHpCvgDLlcEoKkRWABYL5qI6poXrpFUvzV/jj/H5YoANDUCCwDr+Xw+nZh7YqoNwHsILACsFwlG9OzoZ90uA4CLWHQLAACsR2ABAADWI7AAsF5tvFaT/3eyJv/vZLbmBzyKNSwArJc0Sa3/aH2qDcB7CCwArJcVyNLMs2am2gC8h8ACwHpBf1Df6PkNt8sA4CLWsAAAAOsxwwLAeolkQu/96z1J0sntT2ZrfsCDCCwArBdNRHX58sslsTU/4FUEFgDW8/l86tKqS6oNwHsILACsFwlG9MKlL7hdBgAXsegWAABYj8ACAACsR2ABYL1oIqof/PkH+sGff6BoIup2OQBcwBoWANZLJBN6qeylVFu8qxnwHAILAOuFAiHdUXRHqg3AewgsAKwX8od0ae9L3S4DgItYwwIAAKzHDAsA6yVNUu9XvC9J6pnXU34fr7UAryGwALBebbxW33zum5L+vTV/iK35Aa8hsABoFtqF27ldAgAXEVgAWC8nlKO1l611uwwALuJCMAAAsB6BBQAAWI/AAsB60URU09ZO07S109iaH/AoAgsA6yWSCa34YIVWfLCifmt+AJ7DolsA1gsFQrrljFtSbQDeQ2ABYL2QP6QJX5ngdhkAXMQlIQAAYD1mWABYL2mS+qjqI0lS51ad2Zof8CACCwDr1cZr9fU/fF0SW/MDXkVgAdAsRIIRt0sA4CICCwDr5YRytPGKjW6XAcBFXAgGAADWI7AAAADrcUkIwPGrKJOq9znX3yfbj3i4LlGnGRtmSJJ+MuQnygpkOXefAJoFAguA41NRJpUOlmLVzvYbypFy8hsciifj+sPf/iBJuuWMWwgsgAcRWAAcn+p99WFlzFypQ2/n+s3Jl/IKGxwK+UO6YcANqTYA7yGwAGicDr2lLv0zehehQEjX9rs2o/cBwG4sugUAANZjhgWA9Ywx2h/dL0lqF24nn8/nckUAmhqBBYD1auI1Gr5ouCS25ge8iktCAADAesywALBeTihHWydudbsMAC5ihgUAAFiPwAIAAKzHJSEA1qtL1OnBzQ9Kkm4aeBM73QIexAwLAOvFk3E99d5Teuq9pxRPxt0uB4ALmGEBYL2QP6TJp05OtQF4D4EFgPVCgZB+cPoP3C4DgIu4JAQAAKzHDAsA6xljVBOvkSRFghG25gc8iBkWANariddoyIIhGrJgSCq4APAWAgsAALAel4QAWC8SjGjDdzak2gC8h8ACwHo+n49PaAY8jktCAADAesywALBeLBHTnL/MkSRdd9p1CgXYPA7wmrRmWObMmaN+/fqpbdu2atu2rYqKivSnP/0p9f3a2lqVlJQoPz9frVu31tixY1VeXu540QC8JZaMae7WuZq7da5iyZjb5QBwQVqBpVu3brrnnnu0efNmvfHGGzr33HN1ySWX6J133pEk3XTTTXr++ee1ePFirVmzRrt379aYMWMyUjgA7wj6gxp/8niNP3m8gn4mhgEvSuuRf9FFFzX4+u6779acOXO0fv16devWTY899pgWLFigc889V5I0b948nXzyyVq/fr3OPPNM56oG4ClZgSxNGzzN7TIAuOi4F90mEgktXLhQVVVVKioq0ubNmxWLxVRcXJw6p0+fPurevbvWrVt31H6i0agqKysb3AAAAP5T2oFl69atat26tcLhsL73ve9p6dKl+spXvqI9e/YoKytLeXl5Dc7v2LGj9uzZc9T+Zs6cqdzc3NStsLAw7UEAAICWLe3ActJJJ2nLli3asGGDrrvuOk2cOFHvvvvucRcwffp0HThwIHUrKys77r4AtEzVsWqd+sSpOvWJU1Udq3a7HAAuSHv1WlZWlnr16iVJGjhwoDZt2qRf//rXGjdunOrq6lRRUdFglqW8vFydOnU6an/hcFjhcDj9ygEAgGc0euO4ZDKpaDSqgQMHKhQKafXq1anvbdu2TTt37lRRUVFj7waAh0WCEa0Zt0Zrxq1ha37Ao9KaYZk+fbouuOACde/eXQcPHtSCBQv08ssv64UXXlBubq4mTZqkqVOnqn379mrbtq1uuOEGFRUV8Q4hAI3i8/nUPru922UAcFFagWXv3r268sor9dFHHyk3N1f9+vXTCy+8oPPOO0+S9OCDD8rv92vs2LGKRqMaOXKkZs+enZHCAQCAd6QVWB577LHP/X52drZKS0tVWlraqKIA4D/FEjHNe2eeJOnqU65ma37Ag9gyEoD1YsmYfvPWbyRJ408eT2ABPIjAAsB6QX9QY788NtUG4D088gFYLyuQpZ8N/ZnbZQBwUaPf1gwAAJBpBBYAAGA9AgsA61XHqjV4/mANnj+YrfkBj2INC4BmoSZe43YJAFxEYAFgvexgtlaOXZlqA/AeAgsA6/l9fnVt3dXtMgC4iDUsAADAesywALBeLBnTwr8ulCRd1ucyhfzsdAt4DYEFgPViiZju23SfJGnsl8cSWAAPIrAAsF7AH9CFPS5MtQF4D4EFgPXCgbDuPftet8sA4CIW3QIAAOsRWAAAgPUILACsVx2r1tkLz9bZC89ma37Ao1jDAqBZ2B/d73YJAFxEYAFgvexgtpZevDTVPi6fbHewon/LyZfyCp3vF8BhCCwArOf3+dWrXa/j++GcfCmUIz0z2dmipPp+SzYSWoAmQGAB0LLlFdaHiup9zvb7yfb6EFS9j8ACNAECCwDrxZIx/XHHHyVJl/S6JP2dbvMKCRVAM0dgAWC9WCKmn6/7uSTpwh4XsjU/4EEEFgDWC/gDOqfwnFQbgPcQWABYLxwI66FzH3K7DAAuYuM4AABgPQILAACwHoEFgPVq4jUauWSkRi4ZqZp4jdvlAHABa1gAWM8Yo91Vu1NtAN5DYAFgvXAgrN+P+n2qDcB7CCwArBfwB9S3Q1+3ywDgItawAAAA6zHDAsB68WRcKz9cKUn6+n99XUE/T12A1/CoB2C9ukSdpr8yXZJ0buG5BBbAg3jUA7Ce3+fXmZ3PTLUBeA+BBYD1soPZmnv+XLfLAOAiXqoAAADrEVgAAID1CCwArFcTr9HoZ0dr9LOj2Zof8CjWsACwnjFGfz/w91QbgPcQWABYLxwI63cjf5dqA/AeAgsA6wX8AZ3R6Qy3ywDgItawAAAA6zHDAsB68WRca/65RpI0vNtwdroFPIhHPQDr1SXqNOWlKZKkDd/ZQGABPIhHPQDr+X1+9f9S/1QbgPcQWABYLzuYrf+58H/cLgOAi3ipAgAArEdgAQAA1iOwALBebbxWly27TJctu0y18Vq3ywHgAtawALBe0iT1zr53Um0A3kNgAWC9rECWSkeUptoAvIfAAsB6QX9QZ3c72+0yALiINSwAAMB6zLAAsF4imdCGPRskSUM6DVHAH3C5IgBNjcACwHrRRFTfXfVdSfVb8+f4c1yuCEBTI7AAsJ7f59dJ7U5KtQF4D4EFgPWyg9lacvESt8sA4CJeqgAAAOsRWAAAgPUILACsVxuv1dUrr9bVK69ma37Ao1jDAsB6SZPUG+VvpNoAvIfAAsB6WYEs/XL4L1NtAN5DYAFgvaA/qJH/NdLtMgC4iDUsAADAesywALBeIpnQ//vk/0mS+nXox9b8gAcRWABYL5qI6so/XSmJrfkBryKwALCez+dT9zbdU20A3kNgAWC9SDCi5WOWu10GABex6BYAAFiPwAIAAKyXVmCZOXOmzjjjDLVp00YFBQUaPXq0tm3b1uCc2tpalZSUKD8/X61bt9bYsWNVXl7uaNEAvCWaiOr7L35f33/x+4omom6XA8AFaQWWNWvWqKSkROvXr9eqVasUi8V0/vnnq6qqKnXOTTfdpOeff16LFy/WmjVrtHv3bo0ZM8bxwgF4RyKZ0Cu7XtEru15RIplwuxwALkhr0e3KlSsbfP3444+roKBAmzdv1tlnn60DBw7oscce04IFC3TuuedKkubNm6eTTz5Z69ev15lnnulc5QA8IxQI6a5hd6XaALynUe8SOnDggCSpffv2kqTNmzcrFoupuLg4dU6fPn3UvXt3rVu37oiBJRqNKhr9vyneysrKxpQEoAUK+UMa3Wu022UAcNFxL7pNJpOaMmWKhg0bpr59+0qS9uzZo6ysLOXl5TU4t2PHjtqzZ88R+5k5c6Zyc3NTt8LCwuMtCQAAtFDHHVhKSkr09ttva+HChY0qYPr06Tpw4EDqVlZW1qj+ALQ8iWRCf/3XX/XXf/2VNSyARx3XJaHrr79ey5Yt09q1a9WtW7fU8U6dOqmurk4VFRUNZlnKy8vVqVOnI/YVDocVDoePpwwAHhFNRPWt578lia35Aa9Ka4bFGKPrr79eS5cu1Z///Gf16NGjwfcHDhyoUCik1atXp45t27ZNO3fuVFFRkTMVA/Acn8+ngkiBCiIFbM0PeFRaMywlJSVasGCB/vjHP6pNmzapdSm5ubmKRCLKzc3VpEmTNHXqVLVv315t27bVDTfcoKKiIt4hBOC4RYIRrf726i8+EUCLlVZgmTNnjiTpa1/7WoPj8+bN01VXXSVJevDBB+X3+zV27FhFo1GNHDlSs2fPdqRYAADgTWkFFmPMF56TnZ2t0tJSlZaWHndRAAAA/4nPEgJgvWgiqqkvT9XUl6eyNT/gUQQWANZLJBNa9Y9VWvWPVbytGfCoRu10CwBNIRQI6SdDfpJqA/AeAgsA64X8IV3e53K3ywDgIi4JAQAA6zHDAqBJ7aqo0f6qurR+JmmSKq/ZJUnqGOkqv69pX2u1a5WlrnmRJr1PAA0RWACvqCiTqvc5198n29P+kV0VNSr+1RrVxNJcOOurU5s+t0uSDv71TslkpX3fjREJBfTizcMJLYCLCCyAF1SUSaWDpVi1s/2GcqSc/GM+fX9VnWpiCc0a11+9Clof88/VJmpU8kr9+fOuK1J2oOmCw469hzRl0Rbtr6ojsAAuIrAAXlC9rz6sjJkrdejtXL85+VJeYdo/1qugtfp2zU3jJ3K14Yp1ad8PgJaDwAJ4SYfeUpf+blcBAGnjXUIAAMB6BBYA1qtL1OnWV2/Vra/eqrpEeu8wAtAyEFgAWC+ejOu5vz+n5/7+nOLJuNvlAHABa1gAWC/kD2nqwKmpNgDvIbAAsF4oENLVfa92uwwALuKSEAAAsB4zLACslzRJfVz9sSTpSzlfavKt+QG4j8ACwHq18VoVLymWJG34zgblhHJcrghAUyOwAGgWgj6ergAv4xkAgPVyQjl668q33C4DgIu4EAwAAKxHYAEAANbjkhAA69Ul6nTfpvskSbeccYuyAlkuVwSgqTHDAsB68WRci7Yt0qJti9iaH/AoZlgAWC/kD+m6065LtQF4D4EFgPVCgZC+3//7bpcBwEVcEgIAANZjhgWA9YwxOhg7KElqE2ojn8/nckUAmhqBBYD1auI1Gvb7YZLYmh/wKi4JAQAA6zHDAsB6kWBEb054UxKfKQR4FY98ANbz+XwK+Xg7M+BlXBICAADWY4YFgPViiZgeeushSdIPBvxAoQCzLYDXMMMCwHqxZEyPv/O4Hn/nccWSMbfLAeACZlgAWC/kD+mqU65KtQF4D4EFgPVCgZBuHnSz22UAcBGBBWjGdlXUaH9V3Reel/3JIfWStOPjQ6o1BzJf2FHs2HvItfsG0LwRWIBmaldFjYp/tUY1scQXnnuK7wMtD0s3Ltyid1wMLJIUCQXUrlVWWj9jjFHcxCXV78PC1vyA9xBYgGZqf1WdamIJzRrXX70KWn/uudmf5EpLpV9f1l+1HU5togqPrF2rLHXNi6T1MzXxGg1ZMEQSW/MDXkVgAZq5XgWt1bdr7uef5KsPNL2+1Frq8gXnAoCFCCwArBcJRvTa5a+l2gC8h8ACwHo+n09ts9q6WsNnFwzbspD58xzP5TfAVgQWAPgc7VplKRIKaMqiLQ2O27SQ+WgioYBevHk4oQUtAoEFgPViiZjmbp0rSZp86uQm3Zq/a15EL948/LC3j9u0kPlIduw9pCmLtmh/VR2BBS0CgQWA9WLJmOb8ZY4k6apTrmryzxLqmhc5/I8+C5mBJkVgAWC9oD+ocSeNS7UBeA+PfADWywpk6adn/tTtMgC4iE9rBgAA1iOwAAAA6xFYAFivOlatAU8O0IAnB6g6Vu12OQBcwBoWAM3Cpx9+CMCbCCwArJcdzNaLl76YagPwHgILAOv5fX51bNXR7TIAuIjAAgCN8cl2Z/vLyZfyCp3tE2gBCCwArBdLxPTUe09JksafPL7Jd7o9opx8KZQjPTPZ2X5DOVLJRkIL8BkEFgDWiyVjemDzA5KkcSeNsyOw5BXWB4vqfc71+cn2+gBUvY/AAnwGgQWA9YL+oC4+8eJU2xp5hQQLoIlY9MgHgCPLCmTp7q/e7XYZAFzExnEAAMB6BBYAAGA9AgsA61XHqjV0wVANXTCUrfkBj2INC4Bm4WDsoNslAHARgQWA9bKD2Vr2zWWpNgDvIbAAsJ7f59cJbU9wuwwALmINCwAAsB4zLACsF0vGtGT7EknSpb0vVchvwU63AJoUgQWA9WKJmGZsmCFJuuTESwgsgAcRWABYL+AP6LwTzku1AXgPgQWA9cKBsB742gNulwHARSy6BQAA1ks7sKxdu1YXXXSRunTpIp/Pp2effbbB940xuv3229W5c2dFIhEVFxfrb3/7m1P1AgAAD0o7sFRVVem0005TaWnpEb9/33336aGHHtIjjzyiDRs2qFWrVho5cqRqa2sbXSwAb6qJ12jE0yM04ukRqonXuF0OABekvYblggsu0AUXXHDE7xljNGvWLP30pz/VJZdcIkl68skn1bFjRz377LO67LLLGlctAE8yxmhvzd5UG4D3OLro9oMPPtCePXtUXFycOpabm6shQ4Zo3bp1Rwws0WhU0Wg09XVlZaWTJQFoAcKBsBZftDjVBuA9ji663bNnjySpY8eODY537Ngx9b3PmjlzpnJzc1O3wsJCJ0sC0AIE/AH1ad9Hfdr34W3NgEe5/i6h6dOn68CBA6lbWVmZ2yUBAADLOHpJqFOnTpKk8vJyde7cOXW8vLxc/fv3P+LPhMNhhcNM8QI4ulgypuXvL5ckjeo5ip1uAQ9ydIalR48e6tSpk1avXp06VllZqQ0bNqioqMjJuwLgIbFETLe9dptue+02xRIxt8sB4IK0Z1gOHTqkHTt2pL7+4IMPtGXLFrVv317du3fXlClT9Itf/EJf/vKX1aNHD912223q0qWLRo8e7WTdADwk4A/orK5npdoAvCftwPLGG2/onHPOSX09depUSdLEiRP1+OOP65ZbblFVVZWuvfZaVVRU6Ktf/apWrlyp7Oxs56oG4CnhQFizi2e7XQYAF6UdWL72ta997j4IPp9Pd955p+68885GFQYAAPAp198lBAAA8EUILACsVxOv0ahnRmnUM6PYmh/wKEff1gwAmWCM0c6DO1NtAN5DYAFgvXAgrCcveDLVBuA9BBYA1gv4AxpQMMDtMgC4iMAC/Nuuihrtr6pzu4xjtmPvIbdLAIAmQ2ABVB9Win+1RjWxhNulpCUSCqhdqyy3y8i4eDKu1Tvrd9Ae0X2Egn6eugCv4VEPSNpfVaeaWEKzxvVXr4LWbpdzzNq1ylLXvIjbZWRcXaJOP1zzQ0nShu9sILAAHsSjHvgPvQpaq2/XXLfLwGf4fX4N6jgo1QbgPQQWANbLDmZr3tfnuV0GABfxUgUAAFiPwAIAAKzHJSHANhVlUvU+Z/v8ZLuz/TWx2nitxq8YL0l66sKnlB3k098BryGwADapKJNKB0uxauf7DuVIOfnO99sEkiapbfu3pdoAvIfAAtikel99WBkzV+rQ29m+c/KlvEJn+2wi4UBYj573aKoNwHsILICNOvSWuvR3uwprBPwBDe0y1O0yALiIRbcAAMB6zLAAsF48Gdfru1+XJA3tMpSdbtPQ3D5zyiu7NyN9POoBWK8uUaeS1SWS2Jr/WLVrlaVIKKApi7a4XUpaIqGAXrx5OKEFh+FRD8B6fp9fp+Sfkmrji3XNi+jFm4c3u08gn7Joi/ZX1RFYcBgCCwDrZQeztfAbC90uo9npmhfhDz9aDF6qAAAA6xFYAACA9QgsAKxXG6/VhBUTNGHFBNXGa90uB4ALWMMCwHpJk9SWj7ek2gC8h8ACwHpZgSzNOmdWqg3AewgsAKwX9Ac1ovsIt8sA4CLWsAAAAOsxwwLAeolkQm/ufVOSdHrB6Qr4Ay5XBKCpEVgAWC+aiOqaF66RVL81f44/x+WKADQ1AgsA6/l8Pp2Ye2KqDcB7CCwArBcJRvTs6GfdLgOAi1h0CwAArEdgAQAA1iOwALBebbxWk/93sib/72S25gc8ijUsAKyXNEmt/2h9qg3AewgsAKyXFcjSzLNmptoAvIfAAsB6QX9Q3+j5DbfLAOAi1rAAAADrMcMCwHqJZELv/es9SdLJ7U9ma37AgwgsAKwXTUR1+fLLJbE1P+BVBBYA1vP5fOrSqkuqDcB7CCwArBcJRvTCpS+4XQYAF7HoFgAAWI8ZFmTErooa7a+qc7uMY7Zj7yG3SwAAfA4CCxy3q6JGxb9ao5pYwu1S0hIJBdSuFZuS2SiaiOpHa34kSbp/+P0KB8IuVwSgqRFY4Lj9VXWqiSU0a1x/9Spo7XY5x6xdqyx1zYu4XQaOIJFM6KWyl1Jt8a5mwHMILMiYXgWt1bdrrttloAUIBUK6o+iOVLvF+2S7s/3l5Et5hc72CTQxAgsA64X8IV3a+1K3y8i8nHwplCM9M9nZfkM5UslGQguaNQILANgir7A+WFTvc67PT7bXB6DqfQQWNGsEFgDWS5qk3q94X5LUM6+n/L4WvCNDXiHBAjgCAgsA69XGa/XN574p6d9b84fYmh/wGgILgGahXbid2yUAcBGBBYD1ckI5WnvZWrfLAOCiFnwhGAAAtBQEFgAAYD0CCwDrRRNRTVs7TdPWTlM0EXW7HAAuILAAsF4imdCKD1ZoxQcr6rfmB+A5LLoFYL1QIKRbzrgl1QbgPQQWANYL+UOa8JUJbpcBwEUEFgCAVXbsPeR2CWnhk96bBoEFgPWSJqmPqj6SJHVu1bllb83vYe1aZSkSCmjKoi1ul5KWSCigF28eTmjJMAILAOvVxmv19T98XRJb87dkXfMievHm4dpfVed2Kcdsx95DmrJoi/ZX1RFYMozAAqBZiAT5Y+AFXfMi/OHHERFYAFgvJ5SjjVdsdLsMAC7iQjAAALAeMyxAY1SUSdX7nOvvk+3O9QUALQiBBTheFWVS6WApVu1sv6EcKSff2T6bubpEnWZsmCFJ+smQnygrkOVyRc1QJsJwTr6UV+h8v8AREFiA41W9rz6sjJkrdejtXL/8EThMPBnXH/72B0nSLWfcQmBJR05+fQh+ZrLzfYdypJKN/HtFkyCwAI3VobfUpb/bVbRoIX9INwy4IdVGGvIK60OFk5cupfoZm2cm1/dLYGGzuyZAYGkGdlXUNLt9CQAnhQIhXdvvWrfLaL7yCgkVGcJmd00nY4GltLRU999/v/bs2aPTTjtNv/nNbzR48OBM3V2LtauiRsW/WqOaWPP6hNpIKKB2rZi2B9Cysdld08lIYFm0aJGmTp2qRx55REOGDNGsWbM0cuRIbdu2TQUFBZm4yxZrf1WdamIJzRrXX70KWrtdzjFrjtONsJcxRvuj+yVJ7cLt5PP5XK4I+D9sdtc0MhJYHnjgAU2ePFlXX321JOmRRx7R8uXL9bvf/U4//vGPM3GXx6y5Xl7pVdBafbvmulxNM8dbkJutmniNhi8aLomt+QGvcjyw1NXVafPmzZo+fXrqmN/vV3FxsdatW3fY+dFoVNFoNPX1gQMHJEmVlZVOl6bdFTW6+OFXVRtLOt53JmWH/AomalVZyavK41bxT2nuOVK8xtl+gxEpniVl4N8r/k91rFqJmvrLopWVlYqH4i5XBB08JEVN/X/599+sHDpYqWS0WocOVjr6d+XTv9vGGMf6bMA4bNeuXUaSef311xsc/9GPfmQGDx582Pl33HGHkcSNGzdu3LhxawG3srIyp6OFMcYY198lNH36dE2dOjX1dTKZ1L/+9S/l5+e32OvUlZWVKiwsVFlZmdq2bet2OU3Gq+OWvDt2r45b8u7YvTpuybtj/3TcO3fulM/nU5cuXTJyP44Hlg4dOigQCKi8vLzB8fLycnXq1Omw88PhsMLhcINjeXl5TpdlpbZt23rqH/WnvDpuybtj9+q4Je+O3avjlrw79tzc3IyO2/EPP8zKytLAgQO1evXq1LFkMqnVq1erqKjI6bsDAAAekJFLQlOnTtXEiRM1aNAgDR48WLNmzVJVVVXqXUMAAADpyEhgGTdunD7++GPdfvvt2rNnj/r376+VK1eqY8eOmbi7ZiccDuuOO+447FJYS+fVcUveHbtXxy15d+xeHbfk3bE31bh9xmTq/UcAAADOcHwNCwAAgNMILAAAwHoEFgAAYD0CCwAAsB6BpZHWrl2riy66SF26dJHP59Ozzz77uedfddVV8vl8h91OOeWU1DkzZ87UGWecoTZt2qigoECjR4/Wtm3bMjyS9GRi3P/pnnvukc/n05QpU5wvvpEyNfZdu3Zp/Pjxys/PVyQS0amnnqo33ngjgyNJTybGnUgkdNttt6lHjx6KRCI68cQTddddd2Xus0iOU7pjl6T58+frtNNOU05Ojjp37qxrrrlG+/Y1/PDNxYsXq0+fPsrOztapp56qFStWZGgExycT4547d67OOusstWvXTu3atVNxcbE2btyYwVEcn0z9zj+1cOFC+Xw+jR492tnCGylT466oqFBJSYk6d+6scDis3r17p/3vncDSSFVVVTrttNNUWlp6TOf/+te/1kcffZS6lZWVqX379vrWt76VOmfNmjUqKSnR+vXrtWrVKsViMZ1//vmqqqrK1DDSlolxf2rTpk169NFH1a9fP6fLdkQmxr5//34NGzZMoVBIf/rTn/Tuu+/qV7/6ldq1a5epYaQtE+O+9957NWfOHD388MN67733dO+99+q+++7Tb37zm0wN47ikO/bXXntNV155pSZNmqR33nlHixcv1saNGzV58uTUOa+//rouv/xyTZo0SW+99ZZGjx6t0aNH6+23387UMNKWiXG//PLLuvzyy/XSSy9p3bp1Kiws1Pnnn69du3ZlahjHJRNj/9SHH36oH/7whzrrrLOcLrvRMjHuuro6nXfeefrwww+1ZMkSbdu2TXPnzlXXrl3TKy4jn1DkUZLM0qVL0/qZpUuXGp/PZz788MOjnrN3714jyaxZs6aRFWaGk+M+ePCg+fKXv2xWrVplhg8fbm688UbnCs0Ap8Y+bdo089WvftXh6jLHqXGPGjXKXHPNNQ3OGzNmjLniiiucKDMjjmXs999/v+nZs2eDYw899JDp2rVr6utvf/vbZtSoUQ3OGTJkiPnud7/rWK1OcmrcnxWPx02bNm3ME0884USZGeHk2OPxuBk6dKj57//+bzNx4kRzySWXOFytc5wa95w5c0zPnj1NXV1do+phhsVljz32mIqLi3XCCScc9ZwDBw5Iktq3b99UZWXc0cZdUlKiUaNGqbi42KXKMu9IY3/uuec0aNAgfetb31JBQYEGDBiguXPnulil84407qFDh2r16tXavn27JOkvf/mLXn31VV1wwQVulemIoqIilZWVacWKFTLGqLy8XEuWLNGFF16YOmfdunWH/TsfOXKk1q1b19TlOuZYxv1Z1dXVisVizf757VjHfuedd6qgoECTJk1yqVJnHcu4n3vuORUVFamkpEQdO3ZU3759NWPGDCUSifTurFFxBw0ozVedu3btMoFAwCxatOio5yQSCTNq1CgzbNgwByrMDKfG/fvf/9707dvX1NTUGGNMi5xhOdrYw+GwCYfDZvr06ebNN980jz76qMnOzjaPP/64wxU7w6lxJxIJM23aNOPz+UwwGDQ+n8/MmDHD4Wqddaxjf/rpp03r1q1NMBg0ksxFF13U4BVmKBQyCxYsaPAzpaWlpqCgwOmSHeHUuD/ruuuuMz179kw97m3k1NhfeeUV07VrV/Pxxx8bY0yLmGEx5ovHfdJJJ5lwOGyuueYa88Ybb5iFCxea9u3bm5/97Gdp1cMMi4ueeOIJ5eXlfe6iq5KSEr399ttauHBh0xWWYUcad1lZmW688UbNnz9f2dnZ7hWXYUf7nSeTSZ1++umaMWOGBgwYoGuvvVaTJ0/WI4884k6hDjvauJ9++mnNnz9fCxYs0JtvvqknnnhCv/zlL/XEE0+4U6hD3n33Xd144426/fbbtXnzZq1cuVIffvihvve977ldWkalO+577rlHCxcu1NKlS5v94/6Lxn7w4EFNmDBBc+fOVYcOHVyu1jnH8jtPJpMqKCjQb3/7Ww0cOFDjxo3Trbfemv7zW7qJC0enNF51JpNJ06tXLzNlypSjnlNSUmK6detm3n//fYcqzAwnxr106VIjyQQCgdRNkvH5fCYQCJh4PJ6ByhvPqd959+7dzaRJkxocmz17tunSpYsTZTrOqXF369bNPPzwww2O3XXXXeakk05yosyMOJaxjx8/3lx66aUNjr3yyitGktm9e7cxxpjCwkLz4IMPNjjn9ttvN/369XOyXMc4Ne5P3X///SY3N9ds2rTJ6VId58TY33rrrcOe43w+X+o5bseOHRkcwfFx6nd+9tlnmxEjRjQ4Z8WKFUaSiUajx1wPMywuWbNmjXbs2HHE65jGGF1//fVaunSp/vznP6tHjx4uVJgZRxv3iBEjtHXrVm3ZsiV1GzRokK644gpt2bJFgUDApYqd83m/82HDhh321vXt27d/7tqm5uLzxl1dXS2/v+HTUCAQUDKZbKryMuJo45KUest2UVGRVq9e3eCcVatWqaioqGmKzIBjGbck3Xfffbrrrru0cuVKDRo0qElrzJQvGnufPn0Oe467+OKLdc4552jLli0qLCx0o+xGO5bf+bBhw7Rjx44Gj+vt27erc+fOysrKOvY7O+ZogyM6ePCgeeutt1Lp+YEHHjBvvfWW+cc//mGMMebHP/6xmTBhwmE/N378eDNkyJAj9nndddeZ3Nxc8/LLL5uPPvoodauurs7oWNKRiXF/lq1rWDIx9o0bN5pgMGjuvvtu87e//c3Mnz/f5OTkmKeeeiqjY0lHJsY9ceJE07VrV7Ns2TLzwQcfmGeeecZ06NDB3HLLLRkdS7rSHfu8efNMMBg0s2fPNn//+9/Nq6++agYNGmQGDx6cOue1114zwWDQ/PKXvzTvvfeeueOOO0woFDJbt25t8vEdTSbGfc8995isrCyzZMmSBs9vBw8ebPLxfZ5MjP2zbFzDkolx79y507Rp08Zcf/31Ztu2bWbZsmWmoKDA/OIXv0irNgJLI7300ktG0mG3iRMnGmPq/0EOHz68wc9UVFSYSCRifvvb3x6xzyP1J8nMmzcvs4NJQybG/Vm2BpZMjf355583ffv2NeFw2PTp0+eY/z81lUyMu7Ky0tx4442me/fuJjs72/Ts2dPceuutaU0TN4XjGftDDz1kvvKVr5hIJGI6d+5srrjiCvPPf/6zwTlPP/206d27t8nKyjKnnHKKWb58eRON6NhkYtwnnHDCEfu84447mm5gxyBTv/P/ZGNgydS4X3/9dTNkyBATDodNz549zd133532pX6fMZZtKQkAAPAZrGEBAADWI7AAAADrEVgAAID1CCwAAMB6BBYAAGA9AgsAALAegQUAAFiPwAIAAKxHYAEAANYjsAAAAOsRWAAAgPUILAAAwHr/H0JRdIfNffICAAAAAElFTkSuQmCC", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0lElEQVR4nO3de3wU9b3/8ffesrsBkkCUhECCEVHwCqJgkJaKUURQKdGqRYpooUcjClHRtEWLFmJtj6I2eOFBA1YoFhWqYqGSFjhWLgqlR20PoPArKZh4IwmSZLPZnd8fHOcY5LbJXvJNXs/HYx/97uzs7OfbuMN7Z77zHYdlWZYAAADixJnoAgAAQMdC+AAAAHFF+AAAAHFF+AAAAHFF+AAAAHFF+AAAAHFF+AAAAHFF+AAAAHHlTnQBhwuHw9q3b5+6dOkih8OR6HIAAMAJsCxLBw4cUFZWlpzOYx/baHPhY9++fcrOzk50GQAAoAUqKirUq1evY67T5sJHly5dJB0qPiUlJcHVAACAE1FbW6vs7Gz73/FjaXPh46tTLSkpKYQPALZgKKgX/vmCJOmm/jfJ4/IkuCIAR3IiQybaXPgAgCMJhoN6bMtjkqTrz7ie8AEYjPABwAhup1tX97nabgMwF99gAEZIciVp9rDZiS4DQBQwzwcAAIgrwgcAAIgrwgcAI9QF6zR0yVANXTJUdcG6RJcDoBUiCh+hUEgzZ85Ubm6u/H6/+vTpo4cffliWZdnrWJalBx54QD169JDf71d+fr527twZ9cIBdDwHggd0IHgg0WUAaKWIBpz+4he/0NNPP61FixbprLPO0rvvvqtJkyYpNTVVd955pyTp0Ucf1ZNPPqlFixYpNzdXM2fO1MiRI/WPf/xDPp8vJp0A0P753D69/t3X7TYAczmsrx+2OI4xY8YoIyNDCxYssJcVFBTI7/frhRdekGVZysrK0t1336177rlHklRTU6OMjAwtXLhQN9xww3E/o7a2VqmpqaqpqWGSMQAADBHJv98RnXYZOnSoysvLtWPHDknS3//+d7311lsaNWqUJGn37t2qrKxUfn6+/Z7U1FQNGTJEGzZsOOI2A4GAamtrmz0AAED7FdFpl/vvv1+1tbXq16+fXC6XQqGQZs+erfHjx0uSKisrJUkZGRnN3peRkWG/driSkhLNmjWrJbUD6ECC4aBe2vGSJOna06+Vx8kMp4CpIjry8fvf/16LFy/WkiVLtHXrVi1atEi/+tWvtGjRohYXUFxcrJqaGvtRUVHR4m0BaL+CoaDmbJqjOZvmKBgKJrocAK0Q0ZGPe++9V/fff789duOcc87Rv/71L5WUlGjixInKzMyUJFVVValHjx72+6qqqjRgwIAjbtPr9crr9bawfAAdhcvp0mW9L7PbAMwVUfioq6uT09n8YInL5VI4HJYk5ebmKjMzU+Xl5XbYqK2t1aZNm3TbbbdFp2IAHZLX5dVjA6ZLdZ9LVf88+orJ6VJadvwKAxCxiMLHVVddpdmzZysnJ0dnnXWW/va3v+mxxx7TLbfcIunQbXSnTZumn//85+rbt699qW1WVpbGjh0bi/oBdBTVFVLpYOl4E4x5kqXCzQQQoA2LKHw89dRTmjlzpm6//XZ98sknysrK0o9+9CM98MAD9jozZszQwYMHNWXKFFVXV2vYsGFatWoVc3wAaJ26zw8Fj3HzpZNOP/I6n+2QXpl8aF3CB9BmRTTPRzwwzweAI6mv2KQxf7pZ6nSyXr92tfxu/zdX2rdNem64NGWdlDUgzhUCHVsk/35HdOQDABLFkqVP3G4psF9t7DcTgAgRPgAYwetM0rK9H0sFC+R1cYUcYDLCBwAjuBxO9WsMSl16S1xqCxgtoknGAAAAWovwAcAIwXCTVnTupBX71isYZoZTwGScdgHQNlRXHLpE9iiCn/6PZp6cLn0wX5ef90Pu7QIYjPABIPFOYAIxl0P6VkamdMrFTK8OGI7wASDxTmACMa+keUydDrQLhA8AbcdJpzM5GNABMOAUAADEFeEDgBHqm+o1+pXRGv3KaNU31Se6HACtwGkXAEawLEt7Duyx2wDMRfgAYASvy6vnRz1vtwGYi/ABwAgup0sDuw9MdBkAooAxHwAAIK448gHACE3hJpXvKZckXZpzqdxOdl+Aqfj2AjBCY6hR96y7R5K06fubCB+Awfj2AjCC0+HUBRkX2G0A5iJ8ADCCz+1T2RVliS4DQBTw8wEAAMQV4QMAAMQV4QOAERqaGnTtq9fq2levVUNTQ6LLAdAKjPkAYISwFdb2/dvtNgBzET4AGMHr8urZy5612wDMRfgAYASX06WhWUMTXQaAKGDMBwAAiCuOfAAwQlO4SW/ve1uSNDRrKDOcAgbj2wvACI2hRhWWF0pienXAdHx7ARjB6XDqrPSz7DYAcxE+ABjB5/Zp6ZiliS4DQBTw8wEAAMQV4QMAAMQV4QOAERqaGjThjQma8MYEplcHDBdR+DjllFPkcDi+8SgsPDQCvaGhQYWFhUpPT1fnzp1VUFCgqqqqmBQOoGMJW2Ft+3Sbtn26jenVAcNFNOD0nXfeUSgUsp+///77uuyyy3TddddJkqZPn66VK1dq2bJlSk1N1R133KFx48bpr3/9a3SrBtDhJLmSNPeSuXYbgLkiCh8nn3xys+ePPPKI+vTpo+HDh6umpkYLFizQkiVLNGLECElSWVmZ+vfvr40bN+qiiy6KXtUAOhy3061Lcy5NdBkAoqDFYz4aGxv1wgsv6JZbbpHD4dCWLVsUDAaVn59vr9OvXz/l5ORow4YNUSkWAACYr8XzfKxYsULV1dW6+eabJUmVlZVKSkpSWlpas/UyMjJUWVl51O0EAgEFAgH7eW1tbUtLAtCOhcIhbf1kqyTp/O7ny+V0JbgiAC3V4iMfCxYs0KhRo5SVldWqAkpKSpSammo/srOzW7U9AO1TIBTQLatv0S2rb1EgFDj+GwC0WS0KH//617+0Zs0a/fCHP7SXZWZmqrGxUdXV1c3WraqqUmZm5lG3VVxcrJqaGvtRUVHRkpIAtHMOh0N9UvuoT2ofORyORJcDoBVadNqlrKxM3bt31+jRo+1lgwYNksfjUXl5uQoKCiRJ27dv1549e5SXl3fUbXm9Xnm93paUAaAD8bv9WjF2RaLLABAFEYePcDissrIyTZw4UW73/709NTVVt956q4qKitStWzelpKRo6tSpysvL40oXAABgizh8rFmzRnv27NEtt9zyjdcef/xxOZ1OFRQUKBAIaOTIkZo3b15UCgUAAO1DxOHj8ssvl2VZR3zN5/OptLRUpaWlrS4MAL6uoalBU/88VZL01Iin5HP7ElwRgJZq8aW2ABBPYSusjR9vtNsAzEX4AGCEJFeSSr5VYrcBmIvwAcAIbqdbY04dk+gyAERBiycZAwAAaAmOfAAwQigc0j+/+KckqX+3/kyvDhiM8AHACIFQQDeuvFGStOn7m5TsTE5wRQBaivABwAgOh0NZnbLsNgBzET4AGMHv9mv1tasTXQaAKGDAKQAAiCvCBwAAiCvCBwAjBEIB3fnnO3Xnn+9UIBRIdDkAWoExHwCMEAqH9JeKv9htcaUtYCzCBwAjeFwePZj3oN0GYC7CBwAjeJweXXv6tYkuA0AUMOYDAADEFUc+ABghbIW1q3qXJOnUtFPldPDbCTAV4QOAERqaGvTdV78r6X+nV/cwvTpgKsIHAGN09XZNdAkAooDwAcAIyZ5krb9hfaLLABAFnDQFAABxRfgAAABxRfgAYIRAKKD71t+n+9bfx/TqgOEIHwCMEAqH9MbuN/TG7jcOTa8OwFgMOAVgBI/LoxkXzrDbAMxF+ABgBI/TowlnTkh0GQCigNMuAAAgrjjyAcAIYSusjw9+LEnq0akH06sDBiN8ADBCQ1ODrnj5CklMrw6YjvABwBh+t//EVvxsx7FfT06X0rJbXxCAFiF8ADBCsidZm8dvPs5K6ZInWXpl8rHX8yRLhZsJIECCED4AtB9p2YdCRd3nR1/nsx2Hwknd54QPIEEIHwDal7RsQgXQxhE+ABihMdSoOZvmSJJ+POTHSnIlJbgiAC3FtWoAjNAUbtLLO1/WyztfVlO4KdHlAGiFiMPH3r17ddNNNyk9PV1+v1/nnHOO3n33Xft1y7L0wAMPqEePHvL7/crPz9fOnTujWjSAjsfj9GjqwKmaOnCqPE6mVwdMFtFpl/379+viiy/WJZdcoj/+8Y86+eSTtXPnTnXt2tVe59FHH9WTTz6pRYsWKTc3VzNnztTIkSP1j3/8Qz6fL+odANAxeFweTTl3SqLLABAFEYWPX/ziF8rOzlZZWZm9LDc3125blqW5c+fqpz/9qa655hpJ0vPPP6+MjAytWLFCN9xwQ5TKBgAAporotMurr76qCy64QNddd526d++ugQMHav78+fbru3fvVmVlpfLz8+1lqampGjJkiDZs2HDEbQYCAdXW1jZ7AMDhLMvSFw1f6IuGL2RZVqLLAdAKEYWPXbt26emnn1bfvn21evVq3Xbbbbrzzju1aNEiSVJlZaUkKSMjo9n7MjIy7NcOV1JSotTUVPuRnc0lcgC+qb6pXsNfHK7hLw5XfVN9ossB0AoRhY9wOKzzzz9fc+bM0cCBAzVlyhRNnjxZzzzzTIsLKC4uVk1Njf2oqKho8bYAAEDbF9GYjx49eujMM89stqx///56+eWXJUmZmZmSpKqqKvXo0cNep6qqSgMGDDjiNr1er7xebyRlAOiAkj3Jem/ie4kuA0AURHTk4+KLL9b27dubLduxY4d69+4t6dDg08zMTJWXl9uv19bWatOmTcrLy4tCuQAAwHQRHfmYPn26hg4dqjlz5uh73/ueNm/erOeee07PPfecJMnhcGjatGn6+c9/rr59+9qX2mZlZWns2LGxqB8AABgmovBx4YUXavny5SouLtZDDz2k3NxczZ07V+PHj7fXmTFjhg4ePKgpU6aourpaw4YN06pVq5jjA0CrNIYa9fiWxyVJ0wdNZ3p1wGAR39tlzJgxGjNmzFFfdzgceuihh/TQQw+1qjAA+LqmcJNe+OcLkqSpA6cSPgCDcWM5AEbwOD2afM5kuw3AXIQPAEbwuDy68/w7E10GgCjgrrYAACCuOPIBwAiWZdkzm/rdfjkcjgRXBKClOPIBwAj1TfUasmSIhiwZwvTqgOEIHwAAIK447QLACH63X5u+v8luAzAX4QOAERwOh5I9yYkuA0AUcNoFAADEFUc+ABghGArq6b8/LUm67bzb5HEx0RhgKsIHACMEw0HNf2++JOmH5/yQ8AEYjPABwAhup1s39b/JbgMwF99gAEZIciXpvsH3JboMAFHAgFMAABBXhA8AABBXhA8ARqgL1umcRefonEXnqC5Yl+hyALQC4QMAAMQVA04BGMHv9mvd9evsNgBzET4AGMHhcKibr1uiywAQBZx2AQAAccWRDwBGCIaCKvugTJI06axJzHAKGIzwAcAIwXBQT/3tKUnSTf1vInwABiN8ADCC2+lWQd8Cuw3AXHyDARghyZWknw39WaLLABAFDDgFAABxRfgAAABxRfgAYIS6YJ0GLx6swYsHM706YDjGfAAwRn1TfaJLABAFhA8ARvC5fVpVsMpuAzAX4QOAEZwOp3p27pnoMgBEAWM+AABAXHHkA4ARguGglv7PUknSDf1ukMfJDKeAqQgfAIwQDAX16DuPSpIK+hYQPgCDRXTa5Wc/+5kcDkezR79+/ezXGxoaVFhYqPT0dHXu3FkFBQWqqqqKetEAOh6X06Urc6/UlblXyuV0JbocAK0Q8ZGPs846S2vWrPm/Dbj/bxPTp0/XypUrtWzZMqWmpuqOO+7QuHHj9Ne//jU61QLosLwur37x7V8kugwAURBx+HC73crMzPzG8pqaGi1YsEBLlizRiBEjJEllZWXq37+/Nm7cqIsuuqj11QIAAONFfLXLzp07lZWVpVNPPVXjx4/Xnj17JElbtmxRMBhUfn6+vW6/fv2Uk5OjDRs2RK9iAABgtIiOfAwZMkQLFy7UGWecoY8//lizZs3St771Lb3//vuqrKxUUlKS0tLSmr0nIyNDlZWVR91mIBBQIBCwn9fW1kbWAwAdQl2wTle8fIUkaVXBKiV7khNcEYCWiih8jBo1ym6fe+65GjJkiHr37q3f//738vv9LSqgpKREs2bNatF7AXQs+wP7E10CgCho1SRjaWlpOv300/Xhhx8qMzNTjY2Nqq6ubrZOVVXVEceIfKW4uFg1NTX2o6KiojUlAWinfG6fll+9XMuvXs706oDhWhU+vvzyS3300Ufq0aOHBg0aJI/Ho/Lycvv17du3a8+ePcrLyzvqNrxer1JSUpo9AOBwTodTp3U9Tad1PU1OB5MzAyaL6LTLPffco6uuukq9e/fWvn379OCDD8rlcunGG29Uamqqbr31VhUVFalbt25KSUnR1KlTlZeXx5UuAADAFlH4+Pe//60bb7xRn3/+uU4++WQNGzZMGzdu1MknnyxJevzxx+V0OlVQUKBAIKCRI0dq3rx5MSkcQMcSDAf1hw//IEm65rRrmOEUMFhE4WPp0qXHfN3n86m0tFSlpaWtKgoADhcMBTVrw6HB6VfmXkn4AAzGvV0AGMHldOmS7EvsNgBzET4AGMHr8urJEU8mugwAUcCQcQAAEFcc+QAQW9UVUt3nx17nsx3xqQVAm0D4ABA71RVS6WApWHf8dT3JUnL6UV+ub6rX2BVjJUkrxq6Q392yWZUBJB7hA0Ds1H1+KHiMmy+ddPqx101Ol9Kyj/qyZVnad3Cf3QZgLsIHgNg76XQpa0CrNuF1efW70b+z2wDMRfgAYASX06WzTzo70WUAiAKudgEAAHHFkQ8ARmgKN2nV/1slSbrilCvkdrL7AkzFtxeAERpDjSr+r2JJ0ojsEYQPwGB8ewEYwelw6qIeF9ltAOYifAAwgs/t0/zL5ye6DABRwM8HAAAQV4QPAAAQV4QPAEb4anr1sSvGqr6pPtHlAGgFxnwAMIJlWfqo5iO7DcBchA8ARvC6vPrNyN/YbQDmInwAMILL6dKFmRcmugwAUcCYDwAAEFcc+QBghKZwk9b9e50kaXiv4cxwChiMby8AIzSGGjXtL9MkSZu+v4nwARiMby8AIzgdTg04eYDdBmAuwgcAI/jcPv32yt8mugwAUcDPBwAAEFeEDwAAEFeEDwBGaGhq0A2v36AbXr9BDU0NiS4HQCsw5gOAEcJWWB98/oHdBmAuwgcAIyS5klR6aandBmAuwgcAI7idbn2717cTXQaAKGDMBwAAiCuOfAAwQigc0qbKTZKkIZlD5HK6ElwRgJYifAAwQiAU0I/e/JGkQ9OrJzuTE1wRgJYifAAwgtPh1Bldz7DbAMzVqm/wI488IofDoWnTptnLGhoaVFhYqPT0dHXu3FkFBQWqqqpqbZ0AOjif26eXrn5JL139knxuX6LLAdAKLQ4f77zzjp599lmde+65zZZPnz5dr732mpYtW6Z169Zp3759GjduXKsLBQAA7UOLwseXX36p8ePHa/78+eratau9vKamRgsWLNBjjz2mESNGaNCgQSorK9Pbb7+tjRs3Rq1oAABgrhaFj8LCQo0ePVr5+fnNlm/ZskXBYLDZ8n79+iknJ0cbNmw44rYCgYBqa2ubPQDgcA1NDZq0apImrZrE9OqA4SIecLp06VJt3bpV77zzzjdeq6ysVFJSktLS0potz8jIUGVl5RG3V1JSolmzZkVaBoAOJmyF9W7Vu3YbgLkiCh8VFRW666679Oabb8rni86Ar+LiYhUVFdnPa2trlZ2dHZVtA2g/klxJ+tXwX9ltAOaKKHxs2bJFn3zyic4//3x7WSgU0vr16/XrX/9aq1evVmNjo6qrq5sd/aiqqlJmZuYRt+n1euX1eltWPYAOw+10a+QpIxNdBoAoiCh8XHrppXrvvfeaLZs0aZL69eun++67T9nZ2fJ4PCovL1dBQYEkafv27dqzZ4/y8vKiVzUAADBWROGjS5cuOvvss5st69Spk9LT0+3lt956q4qKitStWzelpKRo6tSpysvL00UXXRS9qgF0OKFwSP/92X9Lks496VymVwcMFvUZTh9//HE5nU4VFBQoEAho5MiRmjdvXrQ/BkAHEwgF9IM//kAS06sDpmt1+Fi7dm2z5z6fT6WlpSotLW3tpgHA5nA4lNMlx24DMBf3dgFgBL/br5XjVia6DABRwN2ZAABAXBE+AABAXBE+ABghEAro9jW36/Y1tysQCiS6HACtwJgPAEYIhUP6r73/ZbfFlbaAsQgfAIzgcXn08MUP220A5iJ8ADCCx+nR2NPGJroMAFHAmA8AABBXHPkAYIRQOKSd1TslSX3T+jK9OmAwwgcAIwRCAV332nWSmF4dMB3hA4ARHA6Huvu7220A5iJ8ADCC3+1X+ffKE10GgChgwCkAAIgrwgcAAIgrwgcAIwRCARWtLVLR2iKmVwcMR/gAYIRQOKQ3//Wm3vzXm4emVwdgLAacAjCCx+XRj4f82G4DMBfhA4ARPE6Pbux3Y6LLABAFnHYBAABxxZEPAEYIW2FVHKiQJGV3yZbTwW8nwFSEDwBGaGhq0JjlYyT97/TqHqZXB0xF+ABgjC6eLokuAUAUED4AGCHZk6y3v/92ossAEAWcNAUAAHFF+AAAAHFF+ABghMZQo37y1k/0k7d+osZQY6LLAdAKhA8ARmgKN+nVj17Vqx+9qqZwU6LLAdAKDDgFYASP06OiQUV2G4C5CB8AjOBxeTTp7EmJLgNAFHDaBQAAxBVHPgAYIWyF9Wndp5Kkk5NPZnp1wGCEDwBGaGhqUP5L+ZKYXh0wHeEDgDHcDnZZQHvANxmAEZI9yfrbD/6W6DIAREFEJ02ffvppnXvuuUpJSVFKSory8vL0xz/+0X69oaFBhYWFSk9PV+fOnVVQUKCqqqqoFw0AAMwVUfjo1auXHnnkEW3ZskXvvvuuRowYoWuuuUYffPCBJGn69Ol67bXXtGzZMq1bt0779u3TuHHjYlI4AAAwU0SnXa666qpmz2fPnq2nn35aGzduVK9evbRgwQItWbJEI0aMkCSVlZWpf//+2rhxoy666KLoVQ2gw2kMNerRdx6VJM24cIaSXEkJrghAS7X4WrVQKKSlS5fq4MGDysvL05YtWxQMBpWfn2+v069fP+Xk5GjDhg1H3U4gEFBtbW2zBwAcrincpBe3v6gXt7/I9OqA4SIecPree+8pLy9PDQ0N6ty5s5YvX64zzzxT27ZtU1JSktLS0pqtn5GRocrKyqNur6SkRLNmzYq4cAAdi8fp0W3n3Wa3AZgr4vBxxhlnaNu2baqpqdFLL72kiRMnat26dS0uoLi4WEVFRfbz2tpaZWdnt3h7ANonj8uj2wfcnugyAERBxOEjKSlJp512miRp0KBBeuedd/TEE0/o+uuvV2Njo6qrq5sd/aiqqlJmZuZRt+f1euX1eiOvHAAAGKnV8xOHw2EFAgENGjRIHo9H5eXl9mvbt2/Xnj17lJeX19qPAdDBWZal2sZa1TbWyrKsRJcDoBUiOvJRXFysUaNGKScnRwcOHNCSJUu0du1arV69Wqmpqbr11ltVVFSkbt26KSUlRVOnTlVeXh5XugBotfqmel38u4slMb06YLqIwscnn3yiH/zgB/r444+Vmpqqc889V6tXr9Zll10mSXr88cfldDpVUFCgQCCgkSNHat68eTEpHAAAmCmi8LFgwYJjvu7z+VRaWqrS0tJWFQUAh/O7/do6Yask7vECmI5vMAAjOBwOeRxcYgu0B60ecAoAABAJjnwAMEIwFNSTf3tSknTnwDvlcXEUBDAVRz4AGCEYDmrhBwu18IOFCoaDiS4HQCtw5AOAETxOj24+62a7DcBchA8ARvC4PLr7grsTXQaAKCB8AOiYPttx7NeT06U07jMFxALhA4ARLMtSk9Uk6dA8Hw6Ho2UbSk6XPMnSK5OPvZ4nWSrcTAABYoDwAcAI9U31GrJkiKRWTq+eln0oVNR9fvR1PttxKJzUfU74AGKA8AGg40nLJlQACUT4AGAEv9uvv974V7sNwFyEDwBGcDgcSklKSXQZAKKAScYAAEBcceQDgBGCoaDmvzdfkjT5nMlMrw4YjPABwAjBcFBP//1pSdLNZ91M+AAMRvgAYAS3063rz7jebgMwF99gAEZIciXppxf9NNFlAIgCBpwCAIC4InwAAIC4InwAMEJdsE4Dnx+ogc8PVF2wLtHlAGgFxnwAMMZXN5YDYDbCBwAj+Nw+rbl2jd0GYC7CBwAjOB1OZXTKSHQZAKKAMR8AACCuOPIBwAjBUFAv/PMFSdJN/W9ihlPAYIQPAEYIhoN6bMtjkqTrz7ie8AEYjPABwAhup1tX97nabgMwF99gAEZIciVp9rDZiS4DQBQw4BQAAMQV4QMAAMQV4QOAEeqCdRq6ZKiGLhnK9OqA4RjzAcAYB4IHEl0CgCggfAAwgs/t0+vffd1uAzBXRKddSkpKdOGFF6pLly7q3r27xo4dq+3btzdbp6GhQYWFhUpPT1fnzp1VUFCgqqqqqBYNoONxOpzqndJbvVN6y+ngjDFgsoi+wevWrVNhYaE2btyoN998U8FgUJdffrkOHjxorzN9+nS99tprWrZsmdatW6d9+/Zp3LhxUS8cAACYKaLTLqtWrWr2fOHCherevbu2bNmib3/726qpqdGCBQu0ZMkSjRgxQpJUVlam/v37a+PGjbrooouiVzmADiUYDuqlHS9Jkq49/Vp5nMxwCpiqVWM+ampqJEndunWTJG3ZskXBYFD5+fn2Ov369VNOTo42bNhwxPARCAQUCATs57W1ta0pCUA7FQwFNWfTHEnSNX2uIXwABmtx+AiHw5o2bZouvvhinX322ZKkyspKJSUlKS0trdm6GRkZqqysPOJ2SkpKNGvWrJaWAaCDcDlduqz3ZXYbgLlaHD4KCwv1/vvv66233mpVAcXFxSoqKrKf19bWKjs7u1XbBND+eF1ePfadxxJdBoAoaFH4uOOOO/T6669r/fr16tWrl708MzNTjY2Nqq6ubnb0o6qqSpmZmUfcltfrldfrbUkZABKtukKq+/zor3+2I361ADBGROHDsixNnTpVy5cv19q1a5Wbm9vs9UGDBsnj8ai8vFwFBQWSpO3bt2vPnj3Ky8uLXtUAEq+6QiodLB1vtlFPspScHp+aABghovBRWFioJUuW6A9/+IO6dOlij+NITU2V3+9Xamqqbr31VhUVFalbt25KSUnR1KlTlZeXx5UuQHtT9/mh4DFuvnTS6UdfLzldSmv9qdT6pnqNeWWMJOn1ca/L7/a3epsAEiOi8PH0009Lkr7zne80W15WVqabb75ZkvT444/L6XSqoKBAgUBAI0eO1Lx586JSLIA26KTTpawBMf8Yy7L0Sf0ndhuAuSI+7XI8Pp9PpaWlKi0tbXFRAHA4r8urZVcts9sAzMW9XQAYweV0qV+3fokuA0AUcIMEAAAQVxz5AGCEYDiolbtWSpJGnzqaGU4BgxE+ABghGApq5l9nSpIu73054QMwGOEDgBFcTpe+1fNbdhuAuQgfAIzgdXk1L5/L9oH2gAGnAAAgrggfAAAgrggfAIxQ31Sv0a+M1uhXRqu+qT7R5QBoBcZ8ADCCZVnac2CP3QZgLsIHACN4XV49P+p5uw3AXIQPAHG1t7pe+w82tui9Hp0qSfrnx19Gs6Rv8H32pU6T9OGnX6rBqlHXTknqmcZddIFoIXwAiJu91fXK/891qg+GEl3KMZ3l2K2VXumupdv0gVUjv8elNXcPJ4AAUUL4ABA3+w82qj4Y0tzrB+i07p0jem8o3KR3Pl0vSbrw5G/L5Yzd7sv3Waq0XHrihgF6P5yraS9u0/6DjYQPIEoIHwDi7rTunXV2z9SI3lMXrNNNf35AkrTp+5uU7EmORWmHOA4Fo9NO7qwGK7KQBOD4CB8AjOB0OHVBxgV2G4C5CB8AjOBz+1R2RVmiywAQBfx8AAAAcUX4AAAAccVpF8BgrZkzo7UOnwvjRHz4Scvn52hoatBNb9wkSXrhyhfkc/tavC0AiUX4AAyV6DkzDp8L40T5PS517ZQU8eeFrbC2799utwGYi/ABGKo1c2ZEw9fnwmg46ZwTfl9LZwv1urx69rJn7TYAcxE+AMO1ZM6MqPjaXBjKiv3nu5wuDc0aGvPPARB7hA8AOAGtGa+SKNyTBm0V4QOAEZrCTXp739uSpKFZQ+WO4fTqX9e1U5L8HpemvbgtLp8XTdyTBm0V4QOAERpDjSosL5R0aHr1eIWPnml+rbl7eMKuKmqpDz/5knvSoM0ifAAwgtPh1FnpZ9nteOqZ5ucfcCCKCB8AjOBz+7R0zNJElwEgCpjhFAAAxBXhAwAAxBXhA4ARGpoaNOGNCZrwxgQ1NDUkuhwArcCYDwBGCFthbft0m90GYC7CBwAjJLmSNPeSuXYbgLkIHwCM4Ha6dWnOpYkuA0AURDzmY/369brqqquUlZUlh8OhFStWNHvdsiw98MAD6tGjh/x+v/Lz87Vz585o1QsAAAwXcfg4ePCgzjvvPJWWlh7x9UcffVRPPvmknnnmGW3atEmdOnXSyJEj1dDAADEALRcKh/RO5Tt6p/IdhcKhRJcDoBUiPu0yatQojRo16oivWZaluXPn6qc//amuueYaSdLzzz+vjIwMrVixQjfccEPrqgXQYQVCAd2y+hZJh6ZXT3YmJ7giAC0V1Uttd+/ercrKSuXn59vLUlNTNWTIEG3YsOGI7wkEAqqtrW32AIDDORwO9Untoz6pfeRwOBJdDoBWiOqA08rKSklSRkZGs+UZGRn2a4crKSnRrFmzolkGgHbI7/ZrxdgViS4DQBQkfJKx4uJi1dTU2I+KiopElwQAAGIoquEjMzNTklRVVdVseVVVlf3a4bxer1JSUpo9AABA+xXV8JGbm6vMzEyVl5fby2pra7Vp0ybl5eVF86MAdDANTQ2a/KfJmvynyUyvDhgu4jEfX375pT788EP7+e7du7Vt2zZ169ZNOTk5mjZtmn7+85+rb9++ys3N1cyZM5WVlaWxY8dGs24AHUzYCmvjxxvtNgBzRRw+3n33XV1yySX286KiIknSxIkTtXDhQs2YMUMHDx7UlClTVF1drWHDhmnVqlXy+XzRqxpAh5PkSlLJt0rsNgBzRRw+vvOd78iyrKO+7nA49NBDD+mhhx5qVWEA8HVup1tjTh2T6DIAREHCr3YBAAAdCzeWA/7X3up67T/YmOgyTtiHn3yZ6BLiKhQO6Z9f/FOS1L9bf7mcrgRXBKClCB+ADgWP/P9cp/qgWfcM8Xtc6tqpY4x/CIQCunHljZKYXh0wHeEDkLT/YKPqgyHNvX6ATuveOdHlnLCunZLUM82f6DLiwuFwKKtTlt0GYC7CB/A1p3XvrLN7pia6DByB3+3X6mtXx/dDP9tx7NeT06W07PjUArQjhA8AOFxyuuRJll6ZfOz1PMlS4eY2HUBMGxvUkY7mdWSEDwA4XFr2oVBR9/nR1/lsx6FwUvd5mwwfXTslye9xadqL2xJdSkT8HpfW3D2cANLOET4AGCEQCujedfdKkn45/Jfyuryx/cC07DYZKk5UzzS/1tw93LgruKa9uE37DzYSPto5wgcAI4TCIf2l4i92W1xpe1w90/z8I442ifABwAgel0cP5j1otwGYi/ABwAgep0fXnn5tossAEAVMrw4AAOKKIx8AjBC2wtpVvUuSdGraqXI6+O0EmIrwAcAIDU0N+u6r35X0v9Ore5heHTAV4QOAMbp6uya6BABRQPgAYIRkT7LW37A+0WUAiAJOmgIAgLgifAAAgLgifAAwQiAU0H3r79N96+9TIBRIdDkAWoExH4iJvdX1xt1TAoeprjj+jdXiKBQO6Y3db0jSoZlOmV4dMBbhA1G3t7pe+f+5TvXBUKJLiYjf41LXTkmJLqNtqK6QSgdLwbpjr+dJPnT7+TjwuDyaceEMuw3AXIQPRN3+g42qD4Y09/oBOq1750SXc8K6dkriJlxfqfv8UPAYN1866fSjr5ecHrc7v3qcHk04c0JcPgtAbBE+EDOnde+ss3umJroMtMZJp0tZAxJdBYB2hvABwAhhK6yPD34sSerRqQfTqwMGI3wAMEJDU4OuePkKSUyv3t6ZNgCcU7aRI3wAMIbfzQ6+PevaKUl+j0vTXtyW6FIi4ve4tObu4QSQCBA+ABgh2ZOszeM3J7qMbzreJcdxHJRrup5pfq25e7hxl+lPe3Gb9h9sJHxEgPABAC2RnH7oUuNXJh97PU+yVLiZAHKCeqb5+Ue8AyB8AEBLpGUfChXHm4jtlcmH1iF8ADbCBwAjNIYaNWfTHEnSj4f8WEmuNjAhXFo2oQJoAcIH0BG1sanTT0RTuEkv73xZkjTjwhltI3wAaBHCB9DRtMGp00+Ex+nR1IFT7TYAcxE+DMBN2hBVbXDq9BPhcXk05dwpiS4DQBTELHyUlpbql7/8pSorK3Xeeefpqaee0uDBg2P1ce0WN2lDzDB1OoAEiUn4ePHFF1VUVKRnnnlGQ4YM0dy5czVy5Eht375d3bt3j8VHtlvcpA04xLIs7Q/slyR19XaVw+FIcEUAWiom4eOxxx7T5MmTNWnSJEnSM888o5UrV+o3v/mN7r///lh85Akz9RQGN2lDR1ffVK/hLw6XxPTqaHtMO92c6B+IUQ8fjY2N2rJli4qLi+1lTqdT+fn52rBhwzfWDwQCCgQC9vOamhpJUm1tbbRL077qel3967fUEAxHfdux5PM45Q41qLaWX3qIggNfSgHr0P/G4HsWK3XBOoXqD51+rK2tVZOnKcEVnQBD/7/GiXOHGpQUbtCdz7+d6FIi4vM49eodw5QVxQDy1b/blmUdf2Uryvbu3WtJst5+++1my++9915r8ODB31j/wQcftCTx4MGDBw8ePNrBo6Ki4rhZIeFXuxQXF6uoqMh+Hg6H9cUXXyg9Pd2oc7q1tbXKzs5WRUWFUlJSEl1OzHWk/nakvkr0tz3rSH2V6G+8WZalAwcOKCsr67jrRj18nHTSSXK5XKqqqmq2vKqqSpmZmd9Y3+v1yuv1NluWlpYW7bLiJiUlpUP8R/6VjtTfjtRXif62Zx2prxL9jafU1NQTWs8Z7Q9OSkrSoEGDVF5ebi8Lh8MqLy9XXl5etD8OAAAYJianXYqKijRx4kRdcMEFGjx4sObOnauDBw/aV78AAICOKybh4/rrr9enn36qBx54QJWVlRowYIBWrVqljIyMWHxcm+D1evXggw9+4xRSe9WR+tuR+irR3/asI/VVor9tmcOyTuSaGAAAgOiI+pgPAACAYyF8AACAuCJ8AACAuCJ8AACAuCJ8HENpaalOOeUU+Xw+DRkyRJs3bz7m+nPnztUZZ5whv9+v7OxsTZ8+XQ0NDfbrJSUluvDCC9WlSxd1795dY8eO1fbt22PdjRMS7b5+3SOPPCKHw6Fp06bFoPKWiUV/9+7dq5tuuknp6eny+/0655xz9O6778ayGycs2v0NhUKaOXOmcnNz5ff71adPHz388MMndk+HGIukr8FgUA899JD69Okjn8+n8847T6tWrWrVNuMt2v1ty/spKTZ/36+0tX1VLPraZvZTUbmhSzu0dOlSKykpyfrNb35jffDBB9bkyZOttLQ0q6qq6ojrL1682PJ6vdbixYut3bt3W6tXr7Z69OhhTZ8+3V5n5MiRVllZmfX+++9b27Zts6688korJyfH+vLLL+PVrSOKRV+/snnzZuuUU06xzj33XOuuu+6KcU9OTCz6+8UXX1i9e/e2br75ZmvTpk3Wrl27rNWrV1sffvhhvLp1VLHo7+zZs6309HTr9ddft3bv3m0tW7bM6ty5s/XEE0/Eq1tHFGlfZ8yYYWVlZVkrV660PvroI2vevHmWz+eztm7d2uJtxlMs+ttW91OWFZv+fqWt7ati0de2tJ8ifBzF4MGDrcLCQvt5KBSysrKyrJKSkiOuX1hYaI0YMaLZsqKiIuviiy8+6md88sknliRr3bp10Sm6hWLV1wMHDlh9+/a13nzzTWv48OFt4gttWbHp73333WcNGzYsNgW3Uiz6O3r0aOuWW25pts64ceOs8ePHR7HyyEXa1x49eli//vWvmy07vB+RbjOeYtHfw7WV/ZRlxa6/bXFfFYu+tqX9FKddjqCxsVFbtmxRfn6+vczpdCo/P18bNmw44nuGDh2qLVu22IfFdu3apTfeeENXXnnlUT+npqZGktStW7coVh+ZWPa1sLBQo0ePbrbtRItVf1999VVdcMEFuu6669S9e3cNHDhQ8+fPj21nTkCs+jt06FCVl5drx44dkqS///3veuuttzRq1KgY9ubYWtLXQCAgn8/XbJnf79dbb73V4m3GSyz6eyRtYT8lxba/bW1fFau+tqn9VKLTT1u0d+9eS5L19ttvN1t+7733WoMHDz7q+5544gnL4/FYbrfbkmT9x3/8x1HXDYVC1ujRo495ZCQeYtXX3/3ud9bZZ59t1dfXW5ZltZlfE7Hqr9frtbxer1VcXGxt3brVevbZZy2fz2ctXLgwJv04UbHqbygUsu677z7L4XBYbrfbcjgc1pw5c2LShxPVkr7eeOON1plnnmnt2LHDCoVC1p/+9CfL7/dbSUlJLd5mvMSiv4drK/spy4pdf9vivipWfW1L+ymOfETJ2rVrNWfOHM2bN09bt27VK6+8opUrV+rhhx8+4vqFhYV6//33tXTp0jhX2nrH62tFRYXuuusuLV68+BtJ3EQn8rcNh8M6//zzNWfOHA0cOFBTpkzR5MmT9cwzzySw8pY5kf7+/ve/1+LFi7VkyRJt3bpVixYt0q9+9SstWrQogZVH7oknnlDfvn3Vr18/JSUl6Y477tCkSZPkdLbPXWOk/TV5PyUdv7/taV91In/bNrWfinvcMUAgELBcLpe1fPnyZst/8IMfWFdfffUR3zNs2DDrnnvuabbst7/9reX3+61QKNRseWFhodWrVy9r165dUa27JWLR1+XLl1uSLJfLZT8kWQ6Hw3K5XFZTU1OsunNcsfrb5uTkWLfeemuzdebNm2dlZWVFr/gWiFV/e/Xq9Y3zyw8//LB1xhlnRK/4CLWkr1+pr6+3/v3vf1vhcNiaMWOGdeaZZ7Z6m7EWi/5+XVvaT1lWbPrbVvdVsfrbtqX9VPuM962UlJSkQYMGqby83F4WDodVXl6uvLy8I76nrq7uG78eXC6XJNmXH1qWpTvuuEPLly/Xn//8Z+Xm5saoBycuFn299NJL9d5772nbtm3244ILLtD48eO1bds2e91EiNXf9uKLL/7G5Yg7duxQ7969o1l+xGLV36OtEw6Ho1l+RFrS16/4fD717NlTTU1Nevnll3XNNde0epuxFov+Sm1zPyXFpr9tdV8Vq79tm9pPxT3uGGLp0qWW1+u1Fi5caP3jH/+wpkyZYqWlpVmVlZWWZVnWhAkTrPvvv99e/8EHH7S6dOli/e53v7N27dpl/elPf7L69Oljfe9737PXue2226zU1FRr7dq11scff2w/6urq4t6/r4tFXw/XFs6jfiUW/d28ebPldrut2bNnWzt37rQWL15sJScnWy+88ELc+3e4WPR34sSJVs+ePe1LbV955RXrpJNOsmbMmBH3/n1dpH3duHGj9fLLL1sfffSRtX79emvEiBFWbm6utX///hPeZiLFor9tdT9lWbHp7+Hayr4qFn1tS/spwscxPPXUU1ZOTo6VlJRkDR482Nq4caP92vDhw62JEyfaz4PBoPWzn/3M6tOnj+Xz+azs7Gzr9ttvb/aHl3TER1lZWfw6dRTR7uvh2soX+iux6O9rr71mnX322ZbX67X69etnPffcc3HqzfFFu7+1tbXWXXfdZeXk5Fg+n8869dRTrZ/85CdWIBCIY6+OLJK+rl271urfv7/l9Xqt9PR0a8KECdbevXsj2maiRbu/bXk/ZVmx+ft+XVvaV8Wir21lP+WwrDYwJSEAAOgwGPMBAADiivABAADiivABAADiivABAADiivABAADiivABAADiivABAADiivABAADiivABAADiivABAADiivABAADiivABAADi6v8DeyrgvNLFPp0AAAAASUVORK5CYII=", "text/plain": [ "
" ] @@ -76,40 +88,61 @@ " return np.exp(-np.abs(x**3))\n", "\n", "# Compute the integral using quad form comparison\n", - "I0, err = scipy.integrate.quad(f, -np.inf, np.inf)\n", - "print(\"Integral from scipy.quad = %g with error %g\" % (I0, err))\n", + "I0, err, info = scipy.integrate.quad(f, 0.0, np.inf, full_output = True)\n", + "print(\"Integral from scipy.quad = %g with error %g (%d function evaluations)\" % \n", + " (I0, err, info['neval']))\n", "\n", + "# Number of samples for the MC integration\n", "N = 10**4\n", "\n", + "# Now do the integration 1000 times for both uniform and Gaussian\n", "n_trials = 1000\n", "I_arr = np.zeros(n_trials)\n", "I2_arr = np.zeros(n_trials)\n", "\n", "for i in range(n_trials):\n", - " I_arr[i] = integrate_MC(f, N)\n", - " I2_arr[i] = integrate_MC(f, N, sampling = 'gaussian')\n", - "\n", + " I_arr[i], err = integrate_MC(f, N)\n", + " I2_arr[i], err2 = integrate_MC(f, N, sampling = 'gaussian')\n", + " \n", "I_mean = np.mean(I_arr)\n", "I_std = np.std(I_arr)\n", - "print(\"Uniform: I = %g +- %g; frac err = %g\" % (I_mean, I_std, I_std/I_mean))\n", + "print(\"Uniform: I = %g +- %g; frac err = %g; error estimate = %g\" % (I_mean, I_std, I_std/I_mean, err))\n", "\n", "I2_mean = np.mean(I2_arr)\n", "I2_std = np.std(I2_arr)\n", - "print(\"Gaussian: I = %g +- %g; frac err = %g\" % (I2_mean, I2_std, I2_std/I2_mean))\n", + "print(\"Gaussian: I = %g +- %g; frac err = %g; error estimate = %g\" % (I2_mean, I2_std, I2_std/I2_mean, err2))\n", "\n", "counts, bins = np.histogram(I_arr, bins=10, density = True)\n", "plt.stairs(counts, bins) \n", "counts, bins = np.histogram(I2_arr, bins=10, density = True)\n", - "plt.stairs(counts, bins) \n", + "plt.stairs(counts, bins)\n", "plt.plot([I0, I0], (0,1.05*max(counts)), \":\")\n", "plt.ylim((0,1.05*max(counts)))\n", "plt.show()" ] }, + { + "cell_type": "markdown", + "id": "b2344769", + "metadata": {}, + "source": [ + "Notes:\n", + "\n", + "- if you try different $N$ values, you'll see that the fractional error scales $\\propto 1/\\sqrt{N}$.\n", + "\n", + "- since our estimate for the integral is the mean value of the $N$ samples of $f(x)/p(x)$, another way to estimate the error in the integration is to calculate the error in the mean (or ``standard error''),\n", + "\n", + "$$\\sigma_I^2 = {1\\over N} \\mathrm{Var}(f) = {1\\over N} \\left( \\langle f^2\\rangle - \\langle f\\rangle^2\\right).$$\n", + "\n", + "The function `integrate_MC` in the code above does this by returning $\\sigma_I$ as `np.std(f(x)/p(x))/np.sqrt(N)`\n", + "\n", + "The values agree well with the standard deviation calculated from the histograms." + ] + }, { "cell_type": "code", "execution_count": null, - "id": "e3558622", + "id": "661d095e", "metadata": {}, "outputs": [], "source": []