{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 02.1\n",
"- Compute the following 1D integral via Monte Carlo \n",
"$$I = \\int_0^1 \\frac{\\pi}{2}\\cos(\\pi x/2) dx = 1$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"1. sampling a uniform distribution in $[0,1]$\n",
"2. using importance sampling (i.e. sampling a non-uniform probability in $[0,1]$)\n",
"\n",
"Show a picture of both your estimations of $I$ and their uncertainties with a large number of *throws* $M$ (e.g. $M\\ge 10^4$) as a function of the number of blocks, $N$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The taylor Expansion of the cosine is:\n",
"\n",
"$$\\cos(x) = 1 - \\frac{x^2}{2} + O(x^3).$$\n",
"\n",
"So I choose to sample the integral with distribution:\n",
"$$d(x) = 1 - x^2.$$\n",
"This must be normalized so that $\\int_0^1 d(x) \\,dx = 1$:\n",
"$$\\int_0^1 \\left(1 - x^2\\right) \\,dx = \\left(x - \\frac{x^3}{3}\\right)\\bigg|_0^1 = \\frac{2}{3} \\quad \\implies \\quad d(x) = \\frac{3}{2}\\left(1 - x^2\\right).$$\n",
"Inverting the cumulative distribution requires solving a cubic equation, so I resort to accept-reject to sample $d(x)$.\n",
"$$y = F(x) = \\int_0^x dx' d(x') \\rightarrow x = F^{-1}(y)$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"[]"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAAsTAAALEwEAmpwYAAAka0lEQVR4nO3deXRV9b3+8ffnZCQkEIaEOYAQwMhMAEVUnAEHtAUFVERtFRSttbba3lvbXvvrcLXWolJEBGepCiq1Vto64MCYMI8SQCAyzyQhCUm+vz+S5qYYJEDO3uecPK+1shb77J3k2SvhPNnT92vOOURERAACfgcQEZHQoVIQEZFKKgUREamkUhARkUoqBRERqRTtd4Az0bRpU9euXTu/Y4iIhJXs7Oy9zrmU6taFdSm0a9eOrKwsv2OIiIQVM9tyonU6fSQiIpVUCiIiUkmlICIilVQKIiJSSaUgIiKVPCkFM5tmZrvNbNUJ1puZTTSzHDNbYWa9vcglIiL/yasjhReAwd+yfgiQXvFxJ/BnDzKJiMhxPHlOwTn3qZm1+5ZNhgEvufJxvBeYWbKZtXDO7QhGnvU7j/C3FdsxMwJmBAyiooyYQIDoKCMuOoq46ADxMVEkxP77I5rE+GgaxEeTFB9DbLTOvIlI5AmVh9daAduqLOdWvPaNUjCzOyk/miAtLe20vlnO7jye+jiHM5lKol5MFI0SYmiYEEuT+rE0SYylcf1YUpLiSE2KJzUpjuYN42neMJ4G8TGn/41ERDwUKqVg1bxW7Vu2c24KMAUgMzPztN7Wr+regqu6X4VzDueg1DlKyxzHSssoKXUUlZRRVFJK4bEyCopLKCguJb+ohLyiEo4UlnD46DEOHT3GgYJjHCwoZn9BMVu3FrAvr4j84tJvfL/EuGhaJsfTulECrZLr0aZxPdIa16dtkwTaNkkgITZUfgwiUteFyrtRLtCmynJrYHuwv6mZYQYBjJgoiI+JOuOvWVBcwu7DRew+UsTOw4XsPHSU7QcL2X7wKLkHjpL11X4OF5b8x+e0aBhP+6b16ZCSSHqzRDqmJpKemkRKUtwZ5xERORWhUgqzgQlmNgPoDxwK1vWEYEuIjaZd02jaNa1/wm0OFRxjy/58tuwr4Ku9+Wzem8/Gvfm8s/RrjhT9X2E0qR9LlxZJdG7WgIyWDTinZQM6piYSE6XrGSISHJ6Ugpm9DgwCmppZLvALIAbAOTcZeB8YCuQABcBtXuTyS8OEGLonJNO9dfJ/vO6cY/eRIjbsyuPLXUdYv/MI63Yd4bVFWyg8VgZAbFSAs1sk0a11Q7q3SqZHm2Q6piYSFajuDJyIyKkxdyZXW32WmZnp6sIoqaVljs1781m9/RBrth9mRe4hVn19qPKoIjEumu6tG9I7rRF92jWid1ojGtbTxW0RqZ6ZZTvnMqtdp1IIT2Vljs378lm29SDLth1kydYDrNt5hNIyhxl0bpZE//aN6de+Cf3PakzTRF2fEJFyKoU6Ir+ohOXbDpK15QCLNu8ne8sBjh4rvxuqc7MkzuvQhAEdmnBuhya6TVakDlMp1FHHSstY+fUhFmzax/yN+1j81X4Kj5URFTB6tG7IhZ1SuKhTCt1bJ+uahEgdolIQAIpLyliy9QCfb9jLZzl7WZF7EOegUUIMF6SncOnZqVzUKYXkhFi/o4pIEKkUpFoH8ov5LGcvn6zfzdz1e9iXX0zAILNtYy7LSOXyjOa0/5Zba0UkPKkU5KTKyhzLcw/y0brd/GvtbtbuOAxAemoig7s258pzmnNOywaY6TSTSLhTKcgp27a/gH+t3cU/Vu9i4eZ9lDlo07geQ7u14OpuLenaSgUhEq5UCnJG9uUV8a+1u3h/5U6+yNlLSZkjrXEC1/RowbU9WtG5eZLfEUXkFKgUpNYcLCjmH6t38dcV25m3cR+lZY7OzZK4rlcrhvVsScvken5HFJGTUClIUOzNK+L9lTt4d9l2srccwAzObd+E7/RuxdBuLagfFypDa4lIVSoFCbot+/J5Z+l2Zi3NZcu+AhJioxjStQUjMlvTv31jXX8QCSEqBfGMc47sLQd4KzuX91bsIK+ohPZN6zMiszXD+7QmNSne74gidZ5KQXxxtLiU91fu4C9Z21i0eT/RAePyjGaM6pfGwI5NCegpahFfqBTEd5v25DFj8Tbeys5lf34xaY0TuKl/GiMy29C4vp6gFvGSSkFCRlFJKR+s2smrC7ay6Kv9xEYHuKZ7S8YOaEe31g39jidSJ6gUJCSt33mElxd8xawlX1NQXEqfto0YO6AdQ7o2J1qzy4kEjUpBQtrhwmO8mZXLS/O/Ysu+Alo2jOfWAe0Y2TeNhgka4luktqkUJCyUlTk+Wreb5z/fzPxN+0iIjeLGvm24/fz2tGmc4Hc8kYihUpCws2b7YaZ+tonZy7dT5hxDu7Vg3EUd6NpK1x1EzpRKQcLWjkNHeeGLr3h14Vbyikq4IL0p4y/qwHkdmuiBOJHTpFKQsHe48BivLtjKtC82s+dIET3bJDPh4o5cenaqykHkFKkUJGIUHitl5pJcJs/dyLb9R+nSPIn7Lk1n8DnN9TCcSA2pFCTilJSWMXv5dp7+OIdNe/JJT03k3kvTuapbC803LXISKgWJWKVljr+t3MFTH25gw+480lMT+cFl6Qzt2kJHDiIn8G2loCeEJKxFBYxre7Rkzv0X8vToXgBMeG0pQyd+xpzVOwnnP3pE/KBSkIgQCBhXd2/JB/dfyJ9G9qS4pIy7Xs5m2DNf8Mn63SoHkRpSKUhEiQoYw3q24h8/vJDHhndnf34xY6cvZuSUBWRvOeB3PJGQp2sKEtGKS8p4fdFWnvpoA3vzirk8oxkPDe5Mx1TNKy11ly40S52XX1TC9C828+zcTeQXl3Bj3zbcf1knmjXQpD9S96gURCrszy/mqY828MqCLUQHAnz/wrO468KzNJ+01Cm6+0ikQuP6sfzimnP48IFBXHp2KhM/3MCgxz9hxqKtlJaF7x9IIrXFs1Iws8Fmtt7Mcszs4WrWNzSzv5rZcjNbbWa3eZVN6p60Jgk8Pbo3s+4eQJtG9Xh41kqumvgZ8zbu9TuaiK88KQUziwKeAYYAGcAoM8s4brN7gDXOuR7AIOAPZqZ5GiWoeqc1Yub4ATw9uhdHCksY/dxC7nwpiy378v2OJuILr44U+gE5zrlNzrliYAYw7LhtHJBk5aObJQL7gRKP8kkdZlb+jMOHP7qIH1/Zmc9z9nL5E5/yvx+sI79Iv4JSt3hVCq2AbVWWcyteq+pp4GxgO7AS+IFzruz4L2Rmd5pZlpll7dmzJ1h5pQ6Kj4ninos78vGDg7i6ewsmfbKRix//hHeWfq2H36TO8KoUqhuE5vj/ZVcCy4CWQE/gaTNr8I1Pcm6Kcy7TOZeZkpJS2zlFaNYgnidu7MnM8QNo3jCe+/+yjBunLGDdzsN+RxMJOq9KIRdoU2W5NeVHBFXdBsxy5XKAzUAXj/KJfEOfto14++7z+c313fhy1xGumvg5//PXNRwpPOZ3NJGg8aoUFgPpZta+4uLxSGD2cdtsBS4FMLNmQGdgk0f5RKoVFTBG90/j4x8N4sa+bZg+bzOXPTGX91Zs1ykliUielIJzrgSYAMwB1gJvOOdWm9k4MxtXsdmjwAAzWwl8CDzknNP9gRISGtWP5TfXd+Ptu88nJSmOCa8tZcy0RbpLSSKOnmgWOUWlZY5XFmzhsTnrOVZaxn2XpvP9C84iNlrPgkp40BPNIrUoKmDcOqAd/3rgIi49O5XH5qzn6qc+0yisEhFUCiKnqXnDeCbd1IdpYzPJKyxh+OR5PPLuKl2IlrCmUhA5Q5d0acY/HriIW89rx8sLtnDFHz/lo3W7/I4lclpUCiK1IDEuml9eew6zxg8gKT6a21/I4v4ZS9mfX+x3NJFTolIQqUW90hrx3r0XcP9l6fxt5Q4uf2Iuf1uxw+9YIjWmUhCpZbHRAe6/rBPv3XsBrRrV457XlnD3q9nszSvyO5rISakURIKkc/MkZo0fwE8Gd+Zfa3brqEHCgkpBJIiiowLcPagjf7tvIGmNE7jntSVMeG2JrjVIyFIpiHggvVkSM8cP4MdXdmbO6p1c8cdP+XCt7lCS0KNSEPFIdFSAey7uyOwJA2maGMsdL2bx8MwV5GnOBgkhKgURj53dogHvTjifuwd14I2sbQx+8lMWf7Xf71gigEpBxBdx0VH8ZHAX3hx3HgEzbnx2Po9XjKUk4ieVgoiP+rRtzPs/uIDhfVrz9Mc5fPfP89i4J8/vWFKHqRREfJYYF83/Du/B5Jt7s3V/AVdP/JwZi7ZqvgbxhUpBJEQM7tqCD35wIb3Sknl41krGv7KEgwW6dVW8pVIQCSHNG8bzyh39+emQLny4bhdD/vQZCzft8zuW1CEqBZEQEwgYd13UgZnjBxAXHWDUcwv44z+/pEQXocUDKgWRENW9dTLv3XcB1/VsxZ8+3MDoqQvZeajQ71gS4VQKIiEsMS6aJ27syR9G9GDV14cYOvEzPl6/2+9YEsFUCiJh4Lt9WjN7wkBSk+K4bfpifvv3tXqmQYJCpSASJjqmJvLOPeczun8az87dxKgpC9hx6KjfsSTCqBREwkh8TBS/ub4bfxrZkzU7DnPVxM+Z++Uev2NJBFEpiIShYT1b8dd7B5KSGMfY6Yv44z+/pLRMD7vJmVMpiISpDinlp5Ou71V+d9LY6Ys0T4OcMZWCSBirFxvFH0b04Lff6cbCzfu5euJnLNt20O9YEsZUCiJhzswY1S+NWeMHYGbcMHk+ry3U2ElyelQKIhGia6uGvHfvQM7t0ISfvb2Sn7y1gsJjpX7HkjCjUhCJII3qxzJ9bF/uu6Qjb2bncsOz89l+ULetSs2pFEQiTFTAeOCKzky5pQ+b9uRzzVOfM3+jBtWTmlEpiESoK85pzrsTzqdR/Vhufn4hL3yxWdcZ5KQ8KwUzG2xm680sx8wePsE2g8xsmZmtNrO5XmUTiVQdUhJ5++4BXNw5lV/+dQ0/1nUGOQlPSsHMooBngCFABjDKzDKO2yYZmARc65w7BxjhRTaRSJcUH8OUW/pw36XpvJWdy41TFrDrsEZblep5daTQD8hxzm1yzhUDM4Bhx20zGpjlnNsK4JzTUJAitSQQMB64vBOTb+7Nhl1HuPbpz1mu5xmkGl6VQitgW5Xl3IrXquoENDKzT8ws28zGVPeFzOxOM8sys6w9ezTmi8ipGNy1BTPHDyA6EGDEs/N5Z+nXfkeSEONVKVg1rx1/xSsa6ANcBVwJ/NzMOn3jk5yb4pzLdM5lpqSk1H5SkQh3dosGzJ5wPr3aJHP/X5bx2Jx1lGncJKngVSnkAm2qLLcGtlezzQfOuXzn3F7gU6CHR/lE6pQmiXG8fEd/RvVrwzMfb2T8q9nkF5X4HUtCgFelsBhIN7P2ZhYLjARmH7fNu8AFZhZtZglAf2CtR/lE6pzY6AC/ub4bj1ydwT/X7GLEZD3oJh6VgnOuBJgAzKH8jf4N59xqMxtnZuMqtlkLfACsABYBU51zq7zIJ1JXmRm3D2zPtLF92bq/gOue+YIVuQf9jiU+snB+mCUzM9NlZWX5HUMkIqzfeYQ7XlzM3rwi/nhDT4Z0a+F3JAkSM8t2zmVWt05PNIsIAJ2bJ/HOPeeT0aIB419dwrNzN+oJ6DpIpSAilZomxvHa98/l6u4t+O3f1/Gzt1dxrLTM71jioWi/A4hIaImPiWLiyF60bZLAMx9vJPdAAZNu6k1SfIzf0cQDOlIQkW8IBIwfX9mF/x3enfkb9zFi8nx2HNKdSXWBSkFETuiGzDa8cFs/cg8c5bpnvmDN9sN+R5IgUymIyLcamN6Ut8afR8CMEZPn8dkGDS8TyU65FMysfsWopyJSR3Rp3oC37z6fNo0TuG36YmZm5/odSYLkpKVgZgEzG21mfzOz3cA6YEfFnAePmVl68GOKiN+aN4znjXHn0f+sxvzozeU883GOblmNQDU5UvgY6AD8FGjunGvjnEsFLgAWAL8zs5uDmFFEQkSD+Bimj+3HdT1b8tic9fz83VWUajC9iFKTW1Ivc84dO/5F59x+YCYw08x0r5pIHREbHeCJG3rSrGE8z87dxN4jxTw5sifxMTqrHAlOeqTw70IwsyfNrLohsKmuNEQkcgUCxk+HnM3Pr87gg9U7GfP8Ig4V6G0gEpzKheY8YLaZ1QcwsyvM7IvgxBKRcHDHwPZMHNWLpdsOcMOz8zXNZwSocSk45/4beB34xMw+B34EPBysYCISHq7t0ZLpY/uRe6CA70yax6Y9eX5HkjNQ41Iws0uB7wP5QApwn3Pus2AFE5HwMTC9Ka/feS6Fx0oZPnk+K3MP+R1JTtOpnD76L+DnzrlBwHDgL2Z2SVBSiUjY6d46mTfHnUe9mChGTpnPvI17/Y4kp+FUTh9d4pz7vOLfK4EhwK+DFUxEws9ZKYnMHD+Alsn1GDt9MXNW7/Q7kpyimjy8dqI7jnYAl37bNiJS9zRvGM8bd51XPi/DK9m8paefw0qNHl4zs3vNLK3qixVzLZ9nZi8CtwYlnYiEpUb1Y3n1e/0Z0KEpD765nOlfbPY7ktRQTUphMFAKvG5mO8xsjZltBjYAo4A/OudeCGJGEQlD9eOimXprJldkNONXf13DxA83aFiMMFCTJ5ovBt5xzk2qeHK5KXDUOXcwqMlEJOzFx0Qx6abe/GTmCp7455ccKTzGz4aejc44h66alMJ3gUfNrBnlg+EtA5ab2VJgnXOuNIj5RCTMRUcFeHx4D5Lionnus83kF5fy62FdCQRUDKHopKXgnPsegJndD6QDmyk/epgC7AdaBzGfiESAQMD45bXnkBAXzZ8/2cjR4lIeG96d6ChN6RJqTmWO5tuccz3+vWBmk4Af134kEYlEZsZDg7uQGBfNY3PWU1BcwlOjehMbrWIIJafy0zhsZn3+veCcywY61X4kEYlk91zckZ9fncGc1bu46+UsCo/pDHQoOZUjhduBV8xsDZANdAM0LKKInLI7BrYnPibAf729ijteXMxzYzJJiD2VtyMJllN5onkDMAB4H2gGrAWGBimXiES4m/q35fERPZi/cR9jpy0mr6jE70jCqR0pUHGn0ZsVHyIiZ2R4n9bERQe4/y/LGPP8Ql64vR8N4jVnl590hUdEfHVNj5Y8PaoXK3IPccvUhZqsx2cqBRHx3ZBuLfjzzX1Ys+MwNz2/gIMFxX5HqrNUCiISEi7PaMaUWzL5clceo59byIF8FYMfVAoiEjIu7pLKlFv6kLMnj9FTF7JfxeA5z0rBzAab2XozyzGzE07jaWZ9zazUzIZ7lU1EQsegzqlMHZPJpj15jH5ugYrBY56UgplFAc9QPjFPBjDKzDJOsN3vgTle5BKR0HRhpxSev7Uvm/fmqxg85tWRQj8gxzm3yTlXDMwAhlWz3b3ATGC3R7lEJEQNTG/K1Fsz2bw3n5um6hqDV7wqhVbAtirLuRWvVTKzVsD1wORv+0JmdqeZZZlZ1p49e2o9qIiEjgvSU3huTCYbK64x6K6k4POqFKobI/f42TaeBB462VDczrkpzrlM51xmSkpKbeUTkRB1YacUplYUw83P6zmGYPOqFHKBNlWWWwPbj9smE5hhZl8Bw4FJZnadJ+lEJKRd2CmFZ2/uw5c78xgzbSGHC1UMweJVKSwG0s2sfcXcziOB2VU3cM61d861c861A94C7nbOveNRPhEJcRd3SWXSTb1Zs+Mwt05bpLGSgsSTUnDOlQATKL+raC3whnNutZmNM7NxXmQQkfB3WUYznh7dm5W5h7ht+iIKilUMtc3CeSLtzMxMl5WV5XcMEfHYeyu2c9/rSzn3rCZMG9uX+JgovyOFFTPLds5lVrdOTzSLSNi5untL/nBDD+Zv2sddL2dTVKKJemqLSkFEwtL1vVrzu+90Y+6Xe7j3taUcKy3zO1JEUCmISNi6sW8av7r2HP6xZhcPvLGc0rLwPR0eKjT/nYiEtVsHtKOguJTff7COejEBfved7gQC1T0aJTWhUhCRsDd+UAeOHitl4ocbSIiN5hfXZGCmYjgdKgURiQg/vCyd/KISnv98M4lx0Tx4ZWe/I4UllYKIRAQz47+vOpuC4hKe/jiH+nHRjB/Uwe9YYUelICIRw8z49XXdKq8xJMZHc8u5bf2OFVZUCiISUaICxuMjepBfVMIj766iQXw0w3q2OvknCqBbUkUkAsVEBXh6dG/6tWvMj95YzkfrdvkdKWyoFEQkIsXHRDH11kwyWjZg/CtLWLhpn9+RwoJKQUQiVlJ8DC/c1o82jRP43otZrPr6kN+RQp5KQUQiWuP6sbx8Rz8a1Ith7PRFbN6b73ekkKZSEJGI16JhPV6+ox/Owc1TF7LzUKHfkUKWSkFE6oSzUhJ58fZ+HDp6jDHTNN/ziagURKTO6NqqIVPG9OGrvQXc8WIWR4s15PbxVAoiUqcM6NCUJ0f2ZMnWA9zz2hINuX0clYKI1DlDu7Xg0WFd+Wjdbh6auYJwnoGytumJZhGpk24+ty1784p48l8bSE2K5+EhXfyOFBJUCiJSZ/3g0nT2HCli8tyNpCbFcfvA9n5H8p1KQUTqLDPjf4Z1ZV9eMf/z3hqaJsVxbY+Wfsfyla4piEidFhUwnhzZk37tGvPgG8uZt3Gv35F8pVIQkTovPiaK58Zk0rZJAne9lM26nYf9juQblYKICNAwIYYXb+9H/bhobp22iK8PHvU7ki9UCiIiFVom1+OF2/tSUFzK2GmLOFRwzO9InlMpiIhU0aV5A569pQ9f7cvnzpezKCqpW089qxRERI4zoENTHh/Rg4Wb9/PgmysoK6s7D7fpllQRkWoM69mK7QcL+f0H62iZHM9Ph5ztdyRPqBRERE5g3EVn8fXBAp6du4k2jRK4+dy2fkcKOpWCiMgJmBm/vOYcth8s5JF3V9EyOZ5LujTzO1ZQeXZNwcwGm9l6M8sxs4erWX+Tma2o+JhnZj28yiYiciLRUQGeGtWLjJYNmPDaUlbmRvaUnp6UgplFAc8AQ4AMYJSZZRy32WbgIudcd+BRYIoX2URETqZ+XDTTbu1Lo4RYbn9xcUQ/w+DVkUI/IMc5t8k5VwzMAIZV3cA5N885d6BicQHQ2qNsIiInldognum39aWwuJTbpy/mSGFkPsPgVSm0ArZVWc6teO1E7gD+Xt0KM7vTzLLMLGvPnj21GFFE5Nt1apbEn2/uw8Y9edz9amRO0ONVKVg1r1V746+ZXUx5KTxU3Xrn3BTnXKZzLjMlJaUWI4qInNzA9Kb8v+u78tmGvTzy7qqIm6DHq7uPcoE2VZZbA9uP38jMugNTgSHOuX0eZRMROSU39k1jy74CJn2ykbOaJvL9C8/yO1Kt8epIYTGQbmbtzSwWGAnMrrqBmaUBs4BbnHNfepRLROS0PHhFZ4Z2a85v/r6WOat3+h2n1nhSCs65EmACMAdYC7zhnFttZuPMbFzFZo8ATYBJZrbMzLK8yCYicjoCAeOJG3rSvXUy989YxqqvI+NWVQvn82GZmZkuK0vdISL+2X2kkOufmUdJWRnv3jOQ5g3j/Y50UmaW7ZzLrG6dBsQTETkDqUnxPD82k7zCEr730mIKikv8jnRGVAoiImeoS/MGPDW6F2u2H+aHf1kW1qOqqhRERGrBJV2a8V9XZTBn9S4e+8d6v+OcNg2IJyJSS24/vx0b9+Tx50820jElke/2Cb+BGXSkICJSS8yMX117DgM6NOGns1aSvWW/35FOmUpBRKQWxUQFmHRTb1omx3PXy9nkHijwO9IpUSmIiNSy5IRYpt7al6KSMr73Yhb5ReFzR5JKQUQkCDqmJvL06N58uesID7wRPnckqRRERILkok4plXckPfmv8Bi9R3cfiYgE0e3nt2P9zsNM/CiHTs2TuLp7S78jfSsdKYiIBJGZ8eh1Xcls24gH31we8mMkqRRERIIsLjqKybf0oXFCLHe+lMWeI0V+RzohlYKIiAeaJsYxZUwm+wuKGf9KNsUloTlrm0pBRMQjXVs15LHhPcjacoBfzA7NWdt0oVlExEPX9GjJ2h2HmfTJRjJaNuSWc9v6Hek/6EhBRMRjP7qiM5d0SeVXs1ezaHNoDYWhUhAR8VhUwHhyZE/SGidw96vZbD941O9IlVQKIiI+aBAfw5QxmRQeK+POl7MoPFbqdyRApSAi4puOqYk8eWNPVn19mJ+9vTIkLjyrFEREfHRZRjN+eFknZi35mhfnfeV3HJWCiIjf7r2kI5ed3YxH/7aWBZv2+ZpFpSAi4rNAwPjjjT1o2ySBe15d4uuFZ5WCiEgISIqPYcotmRSVlDH+lWzfLjyrFEREQkTH1EQeH9GD5bmH+OXs1b5kUCmIiISQwV2bc8/FHZixeBuvL9rq+fdXKYiIhJgHLu/MhZ1S+MW7q1m69YCn31ulICISYqICxsSRPUltEMfdry5hX553Q22rFEREQlByQiyTb+7D/vxi7n19KSWl3gy1rVIQEQlRXVs15NfXdWXexn08/g9v5nhWKYiIhLARmW24qX8ak+du5INVO4P+/VQKIiIh7pFrMujRJpkH31zOpj15Qf1enpWCmQ02s/VmlmNmD1ez3sxsYsX6FWbW26tsIiKhLC46ikk39SYmyhj/yhIKikuC9r08KQUziwKeAYYAGcAoM8s4brMhQHrFx53An73IJiISDlol1+NPI3vx5e4j/PfbwZvK06sjhX5AjnNuk3OuGJgBDDtum2HAS67cAiDZzFp4lE9EJORd2CmlfETVpV/zWpAebPNqjuZWwLYqy7lA/xps0wrYUXUjM7uT8iMJ0tLSaj2oiEgom3BxR3J259E0MS4oX9+rUrBqXjv+2Kcm2+CcmwJMAcjMzPR/RgoREQ8FAsbEUb2C9/WD9pX/Uy7Qpspya2D7aWwjIiJB5FUpLAbSzay9mcUCI4HZx20zGxhTcRfSucAh59yO47+QiIgEjyenj5xzJWY2AZgDRAHTnHOrzWxcxfrJwPvAUCAHKABu8yKbiIj8H6+uKeCce5/yN/6qr02u8m8H3ONVHhER+SY90SwiIpVUCiIiUkmlICIilVQKIiJSyYI1foYXzGwPsOU0P70psLcW44QD7XPdoH2uG85kn9s651KqWxHWpXAmzCzLOZfpdw4vaZ/rBu1z3RCsfdbpIxERqaRSEBGRSnW5FKb4HcAH2ue6QftcNwRln+vsNQUREfmmunykICIix1EpiIhIpYgvBTMbbGbrzSzHzB6uZr2Z2cSK9SvMrLcfOWtTDfb5pop9XWFm88yshx85a9PJ9rnKdn3NrNTMhnuZLxhqss9mNsjMlpnZajOb63XG2laD3+2GZvZXM1tesc9hPdqymU0zs91mtuoE62v//cs5F7EflA/TvRE4C4gFlgMZx20zFPg75TO/nQss9Du3B/s8AGhU8e8hdWGfq2z3EeWj9Q73O7cHP+dkYA2QVrGc6nduD/b5Z8DvK/6dAuwHYv3Ofgb7fCHQG1h1gvW1/v4V6UcK/YAc59wm51wxMAMYdtw2w4CXXLkFQLKZtfA6aC066T475+Y55w5ULC6gfJa7cFaTnzPAvcBMYLeX4YKkJvs8GpjlnNsK4JwL9/2uyT47IMnMDEikvBRKvI1Ze5xzn1K+DydS6+9fkV4KrYBtVZZzK1471W3Cyanuzx2U/6URzk66z2bWCrgemExkqMnPuRPQyMw+MbNsMxvjWbrgqMk+Pw2cTflUviuBHzjnyryJ54taf//ybJIdn1g1rx1/D25NtgknNd4fM7uY8lIYGNREwVeTfX4SeMg5V1r+R2TYq8k+RwN9gEuBesB8M1vgnPsy2OGCpCb7fCWwDLgE6AD808w+c84dDnI2v9T6+1ekl0Iu0KbKcmvK/4I41W3CSY32x8y6A1OBIc65fR5lC5aa7HMmMKOiEJoCQ82sxDn3jicJa19Nf7f3OufygXwz+xToAYRrKdRkn28DfufKT7jnmNlmoAuwyJuInqv1969IP320GEg3s/ZmFguMBGYft81sYEzFVfxzgUPOuR1eB61FJ91nM0sDZgG3hPFfjVWddJ+dc+2dc+2cc+2At4C7w7gQoGa/2+8CF5hZtJklAP2BtR7nrE012eetlB8ZYWbNgM7AJk9TeqvW378i+kjBOVdiZhOAOZTfuTDNObfazMZVrJ9M+Z0oQ4EcoIDyvzTCVg33+RGgCTCp4i/nEhfGI0zWcJ8jSk322Tm31sw+AFYAZcBU51y1tzaGgxr+nB8FXjCzlZSfWnnIORe2Q2qb2evAIKCpmeUCvwBiIHjvXxrmQkREKkX66SMRETkFKgUREamkUhARkUoqBRERqaRSEBGRSioFERGppFIQEZFKKgWRWlQxX8MKM4s3s/oVY/p39TuXSE3p4TWRWmZmvwbiKR+ELtc591ufI4nUmEpBpJZVjMuzGCgEBjjnSn2OJFJjOn0kUvsaUz7BSxLlRwwiYUNHCiK1zMxmUz4rWHughXNugs+RRGosokdJFfFaxexmJc6518wsCphnZpc45z7yO5tITehIQUREKumagoiIVFIpiIhIJZWCiIhUUimIiEgllYKIiFRSKYiISCWVgoiIVPr/AsrsC+Z951IAAAAASUVORK5CYII=\n",
"text/plain": [
"