{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "korean-liabilities", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 96, "id": "metallic-olympus", "metadata": {}, "outputs": [], "source": [ "from scipy.special import gamma\n", "\n", "#volume of n-dimensional ball\n", "def vol(n,R):\n", " return np.pi**(n/2.)/gamma(n/2.+1)*R**n\n" ] }, { "cell_type": "code", "execution_count": 171, "id": "auburn-duncan", "metadata": {}, "outputs": [], "source": [ "#radius\n", "L=0.5\n", "\n", "#number of cells for Riemann sums\n", "N=2**20\n", "\n", "\n", "#volume estimated with riemann sums\n", "m_riemann=np.zeros(ndim)\n", "#volume estimated by MC sampling\n", "m_sample=np.zeros((ndim,2))\n", "\n", "ndim=10\n", "\n", "for D in range(1,ndim+1):\n", " \n", " #number of cells per dimension\n", " x=N**(1./D)\n", " x1=np.ceil(x).astype(int)\n", " #number of cells\n", " N1=x1**D\n", "\n", " #cell centers\n", " cell_centers=-L+np.arange(x1)*2*L/x1+L/x1\n", " \n", " m=0 #Riemann sum\n", " shape=np.repeat(x1,D)\n", "\n", " for n in range(N1):\n", " \n", " idx=np.unravel_index(n, shape)\n", " y=cell_centers[np.asarray(idx)]\n", " \n", " #add cell to Riemann sum if |y|**2 < r**2\n", " m=m+(np.sum(np.square(y))" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#volume of n-dimensional ball of radius 0.5\n", "x=np.arange(1,ndim+1)\n", "plt.plot(x,np.abs(vol(x,L)-m_riemann)/(vol(x,L)),label='Riemann')\n", "plt.plot(x,np.abs(vol(x,L)-m_sample[:,0])/(vol(x,L)),label='Sample')\n", "plt.yscale('log')\n", "plt.ylabel('error in integral',size=15)\n", "plt.xlabel('D',size=15)\n", "plt.xticks(size=15)\n", "plt.yticks(size=15)\n", "plt.legend(fontsize=15)\n", "plt.grid()" ] }, { "cell_type": "code", "execution_count": 176, "id": "electronic-appendix", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 256 4096 65536 1048576]\n" ] } ], "source": [ "#dimension\n", "\n", "D=4\n", "\n", "N_vec=2**(4*np.arange(2,6))\n", "\n", "print(N_vec)\n", "\n", "#volume estimated with riemann sums\n", "m_riemann=np.zeros(len(N_vec))\n", "#volume estimated by MC sampling\n", "m_sample=np.zeros((len(N_vec),2))\n", "\n", "\n", "for i in range(len(N_vec)):\n", " \n", " N=N_vec[i]\n", " \n", " #number of cells per dimension\n", " x=N**(1./D)\n", " x1=np.ceil(x).astype(int)\n", " #number of cells\n", " N1=x1**D\n", "\n", " #cell centers\n", " cell_centers=-L+np.arange(x1)*2*L/x1+L/x1\n", " \n", " m=0 #Riemann sum\n", " shape=np.repeat(x1,D)\n", "\n", " for n in range(N1):\n", " \n", " idx=np.unravel_index(n, shape)\n", " y=cell_centers[np.asarray(idx)]\n", " \n", " #add cell to Riemann sum if |y|**2 < r**2\n", " m=m+(np.sum(np.square(y))" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(N_vec,np.abs(vol(D,L)-m_riemann),label='Riemann')\n", "plt.plot(N_vec,np.abs(vol(D,L)-m_sample[:,0]),label='Sampling')\n", "plt.xscale('log')\n", "plt.yscale('log')\n", "#plt.plot(N_vec,1./N_vec**(0.5),label)\n", "#plt.plot(N_vec,1./N_vec**(0.25))\n", "#plt.ylabel('error in integral',size=15)\n", "#plt.xlabel('D',size=15)\n", "#plt.xticks(size=15)\n", "#plt.yticks(size=15)\n", "plt.legend(fontsize=15)\n", "plt.grid()\n", "\n", "print(m_riemann)\n", "print(m_sample)\n", "print(vol(D,L))" ] }, { "cell_type": "code", "execution_count": 156, "id": "colored-dietary", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.9999845595017082\n", "1.0003537533306635\n", "1.9999679752502466\n", "2.00233385393464\n", "2.9999506830251526\n", "2.9998787967965876\n", "3.9999373318961835\n", "3.999560467228359\n", "4.999935496966371\n", "5.003683958855469\n", "5.99994728737391\n", "5.998498952469344\n", "7.000798791233446\n", "6.997226212216778\n", "8.079964307047268\n", "7.999139901382353\n", "8.802107735641766\n", "9.004354593649584\n", "6.527005160137841\n", "9.999578851362283\n", "6.569491041688399\n", "10.99711667811596\n", "6.557606750716823\n", "12.001457055152214\n", "49.817052701742696\n", "13.00469270878783\n", "71.89463865065876\n", "13.999227549358515\n", "103.22705098116226\n", "14.99954757724035\n" ] } ], "source": [ "# mean and std of normal\n", "mu=0\n", "sigma=1\n", "\n", "# define cells in [-5*sigma,5*sigma]^D\n", "L=5*sigma\n", "\n", "N=2**20\n", "\n", "#volume estimated with riemann sums\n", "m_riemann=np.zeros(ndim)\n", "#volume estimated by MC sampling\n", "m_sample=np.zeros(ndim)\n", "\n", "ndim=15\n", "\n", "for D in range(1,ndim+1):\n", "\n", " #number of cells per dimension\n", " x=N**(1./D)\n", " x1=np.ceil(x).astype(int)\n", " #number of cells\n", " N1=x1**D\n", "\n", " #cell centers\n", " cell_centers=-L+np.arange(x1)*2*L/x1+L/x1\n", " \n", " m=0 #Riemann sum\n", " \n", " shape=np.repeat(x1,D)\n", "\n", " for n in range(N1):\n", " \n", " idx=np.unravel_index(n, shape)\n", " y=cell_centers[np.asarray(idx)]\n", " \n", " #add norm(mu,sigma,y)*|y|**2 to Riemann sums\n", " m=m+np.exp(-0.5*np.sum(np.square(y)))*1./(2*np.pi)**(D/2.)*np.sum(np.square(y))\n", " \n", " #multiply by cell volume\n", " m_riemann[D-1]=m*(2*L/x1)**D\n", " print(m_riemann[D-1])\n", " \n", " samples=np.random.normal(0,1,(N,D))\n", " \n", " m_sample[D-1]=np.mean(np.sum(np.square(samples),axis=1))\n", "\n", " print(m_sample[D-1])\n", "\n", " " ] }, { "cell_type": "code", "execution_count": 157, "id": "downtown-falls", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEPCAYAAAC6Kkg/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABKw0lEQVR4nO3dd3hUZfbA8e9JD4EACUXpHUQQgYggKoiAiIJlLdgVfyKIZe0N266uvXdEAVeEBURFilgRFZRmoTdpoYcO6cn5/fEOGEISJskkU3I+z5MnzJ07956EJGfedl5RVYwxxpjiCPN3AMYYY4KPJQ9jjDHFZsnDGGNMsVnyMMYYU2yWPIwxxhSbJQ9jjDHFZsnDGGNMsVnyMMYYU2xBmTxEpImIvC8iE/0dizHGVEQSKCvMReQD4Hxgu6q2yXO8D/AqEA6MUNVn8jw3UVUv8eb6NWrU0EaNGvk2aGOMCXELFixIUdWa+Y9H+COYQowC3gA+PHRARMKBN4FeQDIwT0Qmq+rS4l68UaNGzJ8/30ehGmNMxSAi6ws6HjDdVqo6C9iV73AnYLWq/qWqmcA44AJvrykig0RkvojM37Fjhw+jNcaYii1gkkch6gIb8zxOBuqKSKKIvAO0F5EHC3uxqg5X1SRVTapZ86hWlzHGmBIKpG4rr6nqTmCwN+eKSD+gX7Nmzco2KGOMqUACveWxCaif53E9zzGvqeoXqjqoatWqPg3MGGMqskBPHvOA5iLSWESigAHA5OJcQET6icjwvXv3lkmAxhhTEQVM8hCRscAcoKWIJIvIjaqaDdwKzACWAeNVdUlxrmstD2OM8b2AGfNQ1SsKOT4NmFbS63o75rFv3z62b99OVlZWSW9lKqjIyEhq1apFfHy8v0MxptwETPIoK6r6BfBFUlLSTYWds2/fPrZt20bdunWJjY1FRMoxQhPMVJW0tDQ2bXJDcZZATKDYm5bFgvW7+HXtLm46owk1Kkf79Pohnzy8aXls376dunXrUqlSpfILzIQEEaFSpUrUrVuXzZs3W/IwfrNtXzpz1+5i3rpdzF27ixXb9qMKEWFCt+Y1qdHMkkexeNPyyMrKIjY2thyjMqEmNjbWujxNuVFV1u1MZd7aXcz1JIsNu1IBqBQVTocG1Tm3zfGc0rg67etXJzYq3OcxhHzy8JZ1VZnSsJ8fU5ZycpVlW/Yxb92hlsVuUg5kAFC9UiSnNErg2i4NOaVRAq3rxBMZXvZzoUI+edgiQWNMsMnIzuHP5L3MXetaFQvX72Z/RjYAdavFcnqzRE5pnECnRgk0rVmZsLDyf/MS8snDm24rY4zxp9TMbOat283ctTuZt3Y3vyfvITM7F4DmtSrT7+Q6dGqUwCmNE6hbLTC62EM+eVQkjz/+OE888cThx7Vr1yYpKYn//Oc/nHTSSQDMnDmTs846i0WLFtGmTZvCLmWMKUPpWTks3LCbX9bsZPaanfy+cQ/ZuUp4mNCmTjzXdm7IKY0TOKVRAglxUf4Ot0AhnzwqWrdV1apV+fLLLwFYt24djz76KL169WLZsmUkJCTQoUMH5syZQ9OmTf0cqTEVR1ZOLn8m72H26p3M+Wsn89fvJjM7lzCBtnWr8n9nNKFL00SSGlYnLjo4/iwHR5SlUNG6rSIiIujcuTMAnTt3plGjRnTp0oUvv/ySK6+8kvj4+MPPG2PKRk6usnTzPmavSWH2mp3MW7eL1MwcAE44Pp6rT23IaU0T6dQkgfiYSD9HWzIBU57ElI127doBsHGjq2w/c+ZMRITFixcfPic3N5dnnnmGZs2aER0dTYsWLRg9evQR1+nevTuXXHIJI0eOpHHjxlSuXJlrrrmGjIwM5s6dS6dOnahcuTLdu3dnw4YNR7z2gQceoG3btlSuXJl69epx1VVXsXXr1iPOadSoEffccw8vv/wy9erVo3r16gwYMIA9e/YcPudQ7DNnzuTSSy+lcuXKNGnShLfeesuX3zJjii03V1m+dR8jf17LTR/Op/2/vqLfGz/x9PTlJO9O5R8d6vHWVR1Y+Egvpt9xBo/2a03P1rWDNnFABWh5VHSH/pA3bty40HNuu+02Ro8ezaOPPkqHDh34+uuvGThwIImJiZx//vmHz/vll19ISUnh9ddfZ8OGDdx5553Exsby66+/ct999xEXF8ftt9/OoEGDDnedgVuE+dBDD1GnTh127NjBiy++SI8ePVi8eDFhYX+/fxk/fjwnnXQSw4cPJzk5mbvuuouHHnroqORw0003cd111zFo0CDGjh3L0KFDSUpKolOnTr76thlTJFVlbcpBZq9x3VC/rNnJzoOZANRPiOXcNsdzWrNEOjdJpHZ8jJ+jLRshnzxKOubxxBdLWLp5X9kEdQyt68TzWL8TS/z67Gw3pW/9+vXceuutnHzyyVxwQcEbMK5evZq3336bkSNHct111wHQs2dPtmzZwhNPPHFE8jhw4ACff/45h4pMzpw5k/fee48ffviBM888E4DNmzczdOhQUlNTD6/Y/+CDDw5fIycnhy5dulCvXj1++umnw68DVyPqs88+IyLC/VguXbqUcePGHZU8rrjiCoYNGwa4FtEXX3zBpEmTLHmYcvHerL94/6e1bN2XDkDt+GjObFGTLk0T6dIkkfoJFaNSRcgnj4o25rFz504iI/9uCicmJjJv3jyiowsuTfDtt98SFhbGRRdddDjpAJx99tmMHTuWnJwcwsPd6tSkpCTyVidu1qwZUVFRnH766UccA5dEDv17+vTp/Pvf/2bJkiXs2/d3Ql65cuURyeOss846nDgAWrdufbhYZd6vqXfv3of/HRkZSfPmzUlOTvbyO2RMyX3++yaemraMLk0SubVHM05rmkjjGnEVcpFoyCePkirNO39/qlq1Kt988w05OTn88ccf3HPPPVx55ZX8/PPPR3QRHZKSkkJOTg6FlazfsmUL9erVA6BatWpHPBcVFUWVKlWOuG5UlJtWmJ7u3pXNmzeP/v37c9FFF/HAAw9Qq1YtRITOnTsfPueQgq6vqmRkZByRPAo6L/+1jPG1xZv2ct/EP+nUKIHRAzsRFVGxh4wteYSYiIgIkpKSADj11FOJjY3l2muvZcKECVx++eVHnZ+QkEBEREShyaVWrVqliufTTz+lZs2a/O9//zv87mz9+vWluqYx5W3ngQxu/u8CEuKiePOqDhU+cYAlj5B39dVX8+yzz/Lss88WmDx69OhBTk4Oe/fupVevXj6/f1paGpGRkUc068eMGePz+xhTVrJychn68UJ2HMhg4uAu1Kzi2+q0wcqSR4gTER566CGuuuoqvv3228PjF4e0bNmSwYMHM2DAAO677z6SkpJIT09nyZIlrFy5khEjRpTq/r169eKVV17hn//8J/369WP27Nl89NFHpbqmMeXpqanL+OWvXbx4aTtOqlfN3+EEjJBve9ke5nD55ZfTvHlznnvuuQKff/PNN3nkkUf48MMP6du3L9dffz1Tp049YjC7pPr27cuzzz7LJ598Qv/+/fnhhx+YMmVKqa9rTHkYP38jo2avY2DXxvyjYz1/hxNQRFX9HUO5SEpK0vnz5xf43LJlyzjhhBPKOSITauznKLT8tmE3l7/7C0mNqvPhwE5ElEOZ80AkIgtUNSn/8Yr53TDGmCJs35/O4I8WUCs+mjeu7FBhE0dRbMzDGGPyyMzOZchHC9mblsWkIV0Dtqqtv1nyMMaYPB7/YgkL1u/m9Sva07qO7UlfGGuLGWOMx5hf1/PxrxsY0r0p/drV8Xc4Ac2ShzHGAPPX7eLxyUvo1qIm9/Ru6e9wAl5QdluJSBzwFpAJzFRVW3VmjCmxLXvTGPzRQupWi+W1Ae0J98Oe4MEmYFoeIvKBiGwXkcX5jvcRkRUislpEHvAcvhiYqKo3Af3LPVhjTMhIz8ph8H8XkJaZzfBrk6haKXj32ChPAZM8gFFAn7wHRCQceBM4F2gNXCEirYF6wEbPaTnlGKMxJoSoKg9/upg/kvfy0uUn06J2FX+HFDQCJnmo6ixgV77DnYDVqvqXqmYC44ALgGRcAoEivgYRGSQi80Vk/o4dO8oibGNMEBs9ex2fLEzmjrObc86Jx/k7nKASMMmjEHX5u4UBLmnUBSYB/xCRt4EvCnuxqg5X1SRVTapZs2bZRupnjz/+OCJC8+bNC3y+efPmiAiPP/54mdx/+PDhfPbZZz69ZkZGBi+88ALt27cnLi6OSpUqccopp/Diiy+Slpbmk3sc2v7WVDxz1uzk31OX0fOE2txxdsG/N6ZwQTlgrqoHgRu8ObekOwkGo5iYGNauXcv8+fMPl2UHt6fGunXriIkpu+0whw8fTps2bbjwwgt9cr20tDR69+7NokWL+Oc//3l4w6k5c+bw7LPPEhERwR133OGTe5mKJ3l3KkM/XkijxEq8fHk7wmyAvNgCPXlsAurneVzPc8xrFWknwbi4ODp06MC4ceOOSB7jxo2jR48eLFiwwI/RFc+wYcNYuHAhv/76K23atDl8vGfPngwdOpTly5eX6vppaWnExsaWNkwThNIyc7j5vwvIys7lvWuTqBJjA+QlEejdVvOA5iLSWESigAHA5OJcoKJV1R0wYADjx4/nUMFLVWX8+PEMGDCgwPPHjx9P27ZtiY6Opn79+jz88MNHbEc7atQoRIRFixbRq1cv4uLiaNWqFZMmTTp8Tvfu3VmwYAGjR49GRBARRo0adfj5ESNGcOKJJxIdHU3Dhg0Lre57SGpqKu+++y6DBw8+InEckpCQwGmnnQa4nQ4HDhxIkyZNiI2NpUWLFgwbNozMzMzD569btw4RYcyYMVx77bVUq1aNfv36FXr/Y31PTPBSVe7/5E+WbtnHa1e0p0nNyv4OKWgFTPIQkbHAHKCliCSLyI2qmg3cCswAlgHjVXVJca6rql+o6qDCtlkNNRdffDHbtm3jp59+AuDHH39kx44dXHzxxUed+9VXX3H55ZfToUMHPv/8c2677TZeeOEFbr311qPOvfLKK+nfvz+ffvopzZs3Z8CAAYf3DX/rrbdo1aoVffv2Zc6cOcyZM4fzzjsPgOeff54hQ4Zw4YUXMmXKFIYMGcIjjzzCG2+8UejXsGDBAg4ePEifPn0KPeeQlJQUEhISeOmll/jyyy+59957GTlyJLfddttR595zzz1UqVKFCRMm8NBDDxV4veJ8T0zwee/Hv5j8x2bu6d2Ss1qVbpfMii5guq1U9YpCjk8DppX0uiUe85j+AGxdVNLbls5xbeHcZ0r00mrVqtGnTx/GjRvHGWecwbhx4+jTp0+Be5Q/+uijdO/endGjRwMc/mP94IMPMmzYsMN7lwPceeedDBw4EICOHTtSu3ZtpkyZwuDBg2ndujVxcXHUrFmTzp07H37Nvn37eOKJJxg2bBiPPfYY4DaHSk1N5cknn2TIkCFHbU4FsGmT65ls0KDBMb/etm3b8sILLxx+3LVrV+Li4hg4cCCvv/764T3VATp37sybb75Z5PWK8z0xwWXWyh08M305fdsexy3dm/o7nKAXMC2PslLRWh7guq4mTpxIRkYGEydOLLDLKicnh4ULF3LppZcecfzyyy8nNzeXOXPmHHG8d+/eh/+dmJhIrVq1Drc8CjNnzhwOHjzIpZdeSnZ29uGPHj16sG3btmO+Pu/WtYVRVV555RVat25NbGwskZGRXHXVVWRkZLBhw4Yjzj3UGipMcb8nJnis33mQ28b+RovaVXj+knZe/WyZogVMy6OslLjlUcJ3/oGgf//+/N///R8PP/wwBw8eLLB/PyUlhaysLGrXrn3E8UOPd+06cslNtWrVjngcFRVFenp6kXGkpKQAcOKJJxb4/MaNG2nYsOFRx+vWrQvAhg0baNGiRZH3eOWVV7j33nu5//776datG9WrV2fevHkMHTr0qPjyf60FxVuc74kJDgczshn04QJEYPg1ScRFh/yfvXIR8t/FijTb6pC4uDjOP/98Xn75ZS699FLi4uKOOqdGjRpERkayffv2I45v27YNcIPSpXXoGlOmTCnwD3fLlgUXn0tKSiIuLo4ZM2bQs2fPIu8xYcIELrnkEp566qnDx5YuXVrgucd6t1ke3xNTvlSVeyb8wart+/lw4Kk0SKzk75BCRsh3W1VUQ4YMoV+/fgwePLjA58PDw+nYsSMTJkw44vj48eMJCwujS5cuxbpfQS2RLl26EBsby+bNm0lKSjrqo0qVgktBxMbGcvPNN/P2228XmAj27NlzuAspLS2N6OjoI54fM6ZkdTJ9/T0x/vfm96uZvngrD/U9gdOb1/B3OCEl5FseFWmRYF7du3ene/fuRZ7zxBNPcM4553DDDTcwYMAAFi1axCOPPMJNN91U7IHhVq1aMWPGDGbMmEFiYiKNGzcmMTGRxx9/nDvuuIP169dz5plnkpuby8qVK/n+++/59NNPC73ek08+ydy5c+natSt33nknXbt2BeDXX3/l9ddf54EHHqBLly706tWL1157jVNPPZWmTZsyZswYVq9eXazYy+p7YsqOqpKWlcPOA5nsTs1k58FMdh/MZNfBv/+dciCTb5dv48KT63Dj6Y39HXLICfnkURG7rbzVu3dvxo0bx5NPPsmYMWOoVasWd999N0888USxrzVs2DA2bNjAZZddxr59+xg5ciTXX3899913H3Xq1OHll1/mxRdfJCYmhhYtWnD55ZcXeb3Y2Fi++eYbXn/9dT766COeecaNQZ144oncd9993HzzzYCbHbVjxw6GDRsGuKnKr732WpHrOIriy++J8V5OrrI3LYtdBzPYdfDIz4eSwc6DLlHsOuD+nZGdW+C1IsKEhLgoEuKiuOjkuvzn4rY2QF4G5NBislCXlJSk8+fPL/C5ZcuWccIJJ5RzRCbU2M/R0XJylZ0HM9i+L4Nt+9LZvt993rYvgx373edt+9JJOZBBbiF/iipHR5AQF0X1uCgS46KoXimKxMouOSRUijryubgo4mMiLFn4kIgsUNWk/McLbXmIyA7A68yiqgG54qaidlsZU5Zyc5WdBzM9CSHdkxwy2LY/ne15kkTKgUxyCsgKNSpHUatKDLXio2l9fDy14qMP//FPjIumelzk4c/REUevBTL+V1S31ZsUI3kEKuu2MsY3cnKVeyf8wew1O9lxIKPApJAYF0XNKtHUjo+h1XFVqB0fQ60q0dSKj6F2fAy146OpUTmayHCbqxPsCk0eqvp4OcZhjAlwn/++iUm/beKcE2vTrFZlT2JwrYfa8THUrBxNVIQlhYoi5AfMjTGll5Gdw4tfreTEOvG8fVVHK2FuvE8eItIFuBFoARy1MYSqdvJhXD7j7ZiHqtogmymxUJ94MuaXDWzak8bTF7e1xGEALxcJikgvYBZuP43TgR3AAaAdkAgsLqsAS8ub2laRkZE+25nOVExpaWlERobmvhD707N44/vVnNY0kTNsoZ3x8LaD8l/Aq8ChynKPqGoPXCskC5jp+9DKT61atdi0aROpqakh/w7S+JaqkpqayqZNm6hVKyAnHJbaiB/XsutgJvf3aWWtc3OYt91WrYFhQC5uBlYcgKquF5HHgSeAD8siwPIQHx8PwObNm8nKyvJzNCbYREZGUrt27cM/R6Ek5UAGI378i75tj6Nd/Wr+DscEEG+TRzoQpqoqIluApsCPnuf24bqzglp8fHxI/vIbUxpvfLea9Oxc7u5dcBFLU3F5mzz+AFoCXwPfAg+KyCYgE9el5addk4wxZWXDzlTG/Lqey5Lq0dS2azX5eDvm8Qp/Lxh8CDiI2xr2e6AWMNTnkflIRdvD3BhfeenrFYSJcMfZRe+pYiomr5KHqk5T1Tc9/94EdMS1RE4GmqnqgjKLsJQq4k6CxpTW0s37+PyPzdzQtTHHVT1qZr4xx04eIhIjIitFpM+hY+qsUtU/VTWzbEM0xpS352Ysp0p0BEO62V7fpmDHTB6qmg5Uw820MsaEuF/+2snMFTu45axmVK0UmmtXTOl5O+YxBrihLAMxxvifqvLsl8s5Lj6G609r5O9wTADzdrbVBuAyEZkHTAe2cWTFXVXVt30dnDGmfH21dBu/bdjDMxe3JSbSSqGbwnmbPF70fD4eN1ienwKWPIwJYtk5uTw/YwVNasZxScegX7plypi3s63CjvFRrm9RRKSJiLwvIhPL877GhLJJCzexevsB7u3dkgjbb8McQ7n/hIjIByKyXUQW5zveR0RWiMhqEXmgqGuo6l+qemPZRmpMxZGelcPL36ykXf1q9GlznL/DMUHAq24rETmziKdzcSVKVqhqhheXGwW8QZ5aWCISjtu5sBeQDMwTkclAOPB0vtcPVNXt3sRtjPHOh3PWsWVvOi9e1s6KHxqveDvmMZMjB8iFo7eoTReREcBdqppT2IVUdZaINMp3uBOwWlX/AhCRccAFqvo0cL6XMR5FRAYBgwAaNGhQ0ssYE9L2pmXx5vdrOLNFTU5raiXXjXe87bbqiZtx9Q7QF0jyfH4X2Ahchmsh3ISrdVVcdT3XOSTZc6xAIpIoIu8A7UXkwcLOU9Xhqpqkqkk1a9YsQVjGhL7hs9awNy2L+86x4ofGe962PG4FRhewr/kMT0n261W1n4hEANcDD/sswgKo6k5gsDfneruToDEV0fZ96bz/01r6t6tDm7pWwsd4z9uWR2/gp0Ke+xk4y/PvWbjpvMW1Caif53E9z7FSs9pWxhTu1W9XkZ2j3N3bih+a4vE2eewC+hfyXH/P8wCVgJKUr50HNBeRxiISBQwAJpfgOkexqrrGFGxtykHGzdvIlac2oGFinL/DMUHG2+TxHHCriEwWkZtE5ELP5ym4cuzPes47C5cICiUiY4E5QEsRSRaRG1U1G9c1NgNYBoxX1SUl+YLys5aHMQV74asVREeEcVuP5v4OxQQhr8Y8VPUNz+ZPDwJv4abQ5gC/Af9Q1U89pz6N2yCqqGtdUcjxacA0L+P2mo15GHO0Rcl7mfrnFm7r0YyaVaL9HY4JQl4vElTVT1W1ExCDG9eIUdVOeRIHqpqiqvvKIM4Ss5aHMUd79svlVK8UyaAzm/g7FBOkirXCXNzqoTq4PcyDYocYG/Mw5kg/rUrhp9UpDD2rGVVirOS6KRmvk4eI3IKbAbUe+BG3kyAiMklE/lkm0fmAtTyM+Zuq8tyM5dStFsvVnRv6OxwTxLxKHiJyL/AS8B7QA7fC/JCZwOU+j8xHrOVhzN+mLdrKn8l7ubNXCyu5bkrF25bHUOBRVX0M1+rIawUQsJPEreVhjJOVk8sLX62gZe0qXNS+0AIOxnjF2+RxHLCgkOdyCZLxD2MqsvHzN7I25SD3ntOS8DArfmhKx9vksRroVshzZwJLfROOMaYspGXm8Oo3q0hqWJ2zT6jl73BMCPC2ttUrwFsikgkc2oCplojcCNyFK4gYkGydhzHwwc9r2b4/gzev6mAl141PeLuT4AhcscP7gUMrv6cBrwKPq+rHZRNe6dmYh6no9qRm8s4Pazi7VS1OaZTg73BMiPC25YGqPu8pg94FqIGrZzVHVW0akzEB7O2ZaziQkc29fazkuvEdb3cSvBaY6imF/lW+5xKA81X1wwJfbIzxmy170xg1ex0Xta9Lq+Pi/R2OCSHeDpiPxK0qL0hjz/MBydZ5mIrsla9XoQp39QrY2fQmSHmbPIoaYUvE7WEekGzMw1RUq7fvZ8KCjVzduSH1qlfydzgmxBTabSUiFwAX5Dn0iIjsyHdaDHAGxyjDbowpf8/PWEGlqAhu7WEzDY3vFTXmUQtom+dxU9xiwbwycWMgT/o4LmNMKSzcsJsZS7ZxV68WJMRF+TscE4IKTR6q+h6ulhUi8j0wRFWXl1dgxpjiy8zOZezcDbz+3SpqVI7ixtMb+zskE6K83QzqrGOfZYzxl9xc5Ys/N/PiVyvZsCuVUxsn8Gi/1sRFez0b35hi8fonS0TqAOcD9Ti6lpWq6v2+DMxXbIW5CWWqysyVO3juyxUs27KPE46PZ9QNp9CtRU1bSW7KlKjqsU8SuQgYi9t+djtHbzWrqhrQW5IlJSXp/Pnz/R2GMT6zcMNunp2+nF/X7qJBQiXu7t2CfifVIcyKHhofEpEFqpqU/7i3LY//4AbGr1fVXT6NzBhTLKu37+f5GSuYsWQbNSpH8a8LTmTAKQ2IiijWxqDGlIq3yaM+cJslDmP8Z/OeNF79ZhUTFmykUlQEd/VqwY2nN7ZxDeMX3v7UzcZtO/tNGcZijCnAntRM3pq5hlGz14HC9ac1ZuhZTUmsHO3v0EwF5m3yuAsYIyIHgK+BPflPUNVUH8ZlTIWXlpnDBz+v5Z0fXGHDi9vX485ezW21uAkI3iaPPz2fRwKFjbDbhsjG+EBWTi7/m7eR175dxfb9GfQ8oRb3ntOKlsdV8XdoxhzmbfIYSOFJwy9E5ELgPCAeeF9Vvyr6FcYEttxcZdriLbz41UrWphwkqWF13ryqg+3BYQKSt4sER/nypiLyAW7NyHZVbZPneB/cBlPhwAhVfaaImD4DPhOR6sAL5CsVb0ww+WlVCs9+uZxFm/bSsnYV3r8uiR6tatlaDROw/DVNYxTwBnB4DxARCQfeBHoBycA8EZmMSyRP53v9QFXd7vn3MM/rjAk62/alc/f4P/hpdQp1q8Xy4qXtuLB9XcJtrYYJcEVV1Z2LW9exVETmcYxuK1Xt5O1NVXWWiDTKd7gTsFpV//Lcfxxwgao+jWul5I9PgGeA6aq6sJCvYRAwCKBBgwbehmdMuXnr+9XMXbuLR85vzdWdGxAdYUOHJjgU1fJYAqTl+XdZj3nUBTbmeZwMnFrE+bcBPYGqItJMVd/Jf4KqDgeGg1th7sNYjSm13Fxl+uKtnNWqphUwNEGnqKq6N+T59/XlEk0xqOprwGvHOs9qW5lAtWDDbrbvz6Bv2+P9HYoxxRZI9Qw24VayH1LPc6xUbCdBE6imLdpCVEQYZ59Q29+hGFNsgZQ85gHNRaSxiEQBA4DJpb2o7WFuAlFurjJ90Va6tahJZSsvYoKQX5KHiIwF5gAtRSRZRG5U1WzgVmAGsAwYr6pL/BGfMWXtt4172Lovnb5t82/OaUxw8MtbHlW9opDj04BpPr7XF8AXSUlJN/nyusaUxrRFW4gKty4rE7wCqdvKmApBVZm+aAtnNK9BfEykv8MxpkRCPnnYmIcJNH8k72Xz3nTOtVlWJoh5lTxEJFJE7hGR2SKyQUS25/8o60BLymZbmUAzbdEWIsOFXtZlZYKYt2MeLwM3A1OA7zl6G9qAZes8TCBRVaYt2kLXZjWoWsm6rEzw8jZ5XAo8oKovlmUwZcEGzE0gWbRpL8m707j97Ob+DsWYUvF2zEP4e08PY0wJTVu0lYgwoXdr67Iywc3b5PEeUOD0WmOMdw51WZ3WrAbVKkX5OxxjSsXbbqttwFUi8j0Fb0Orqvq2LwPzFRvzMIFiyeZ9bNiVyi3dm/o7FGNKzdvk8YrncwOgWwHPKxCQycPGPEygmL54C+FhQu8TbVW5CX7e7iQY8utBjClLrstqK12aJJIQZ11WJvhZUjCmHCzfup+1KQc512pZmRDhdW0rEamGW+txOpAA7AJ+BIar6p6yCM4XbMzDBIJpi7YQJnCOdVmZEOHtCvOmwCLgX0AcsMHz+V/An57nA5KtMDf+pqpMXbSFUxsnUqNytL/DMcYnirPCfA/QWVUPb9AkInVxVXBfAi7weXTGhICV2w7w146D3NDVtpo1ocPbMY/uwKN5EweA5/G/gLN8HJcxIWPaoi2IwDkn2sJAEzq8TR4KhBdxDfVNOMaEnmmLttCpUQK1qsT4OxRjfMbb5PE98G8RaZj3oOfxv4BvfR2YMaFg1bb9rNp+gL5Wft2EGG+Txz+BaGCViPwiIp+LyBxgFRAF3FVG8ZWa7edh/Gn64q2IQJ82NsvKhBavkoeqrgNaAbcDS4BIYCluz/ETPM8HJJttZfxp2qItJDWsTu1467IyocXrdR6qmgm84/kwxhzDmh0HWL51P4/1a+3vUIzxOVthbkwZmb5oC2BdViY0WfIwpoxMW7SVDg2qcXzVWH+HYozPWfIwpgysSznI0i37bJaVCVmWPIwpA9MWuy6rcy15mBB1zOQhIjEi8p6IdC6PgLwhIieIyDsiMlFEhvg7nmLbvxXGXgHzRvg7ElNGpi/aysn1q1G3mnVZmdB0zOShqunAAMAncw1F5AMR2S4ii/Md7yMiK0RktYg8cIyYlqnqYOAyoKsv4io3yQtgeHdYMQ2+exIyU/0dkfGxDTtTWbRpL32t/LoJYd52W32H7+pXjQL65D0gIuHAm8C5QGvgChFpLSJtRWRKvo9antf0B6biCjMGh98/hpHnQngk9H0B0nbDn//zd1TGx6Yf6rJqY11WJnR5u87jTWCEiMTh/lhvI189K1Vd6s2FVHWWiDTKd7gTsFpV/wIQkXHABar6NHB+IdeZDEwWkanAxwWdIyKDgEEADRo08Ca8spGTDV8Ng1/fhsZnwiWjoFIC/PZf+OVt6Hg9iPgvPuNT0xZt4aR6VamfUMnfoRhTZrxteXwJ1MOVIfkG+BO3v8ciYLHnc2nUBTbmeZzsOVYgEekuIq+JyLsU0fJQ1eGqmqSqSTVr1ixliCWUugs+usgljlOHwNWfQlyiSxadh0LKClhjpcFCRfLuVP5I3muzrEzI87blEVAl11V1JjDTm3P9upPg1sUw7grYvw0ueAvaX3Xk8ydeBF8/4lofzXqWf3zG56Yv2grAubYw0IQ4r5KHqv5QxnFsAurneVzPcyx4LfkUPrsFYqrCDdOhXsejz4mIglNugu+fhB0roGbL8o/Tn3JzQHPdGFCImLZ4CyfWiadhYpy/QzGmTBVrnYeInCoid4vIU57Pp/oojnlAcxFpLCJRuNldk31x4XIvjJibC9/+CyZcD7XbwKCZBSeOQ5JugIgY1/qoSHKyYcyl8HIb+Kus35uUj8170vhtwx7rsjIVgrd7mMeJyDRgDvA0MNDzebaITBURr0cGRWSs5zotRSRZRG5U1Wxchd4ZwDJgvKouKebXUtj9yq8ke/pe103144vQ/hq4fgpUOUb3RVwNOOky+GOcGx+pKL5/yo31iMCHF7iEm5Pt76hKZfpi12VlycNUBN62PJ4DugCXAzGqejxu3ccAz/Fnvb2hql6hqseraqSq1lPV9z3Hp6lqC1VtqqpPFe/LKPJ+5dPySFkF750Nq79x03D7vw4R0d699tQhkJ0GC0aVaYgBY/lU+Okl6HAd3LYA2l/tEu6ovrBng7+jK7Hpi7ZwwvHxNK5hXVYm9HmbPP4B3K+qE1Q1F0BVc1V1AvAAcGlZBRgUVn4F7/WAtF1w7efQ6abiTb2t3RqadIe570FOVpmFGRB2roFPh8DxJ8O5z0FUHFzwBvzjfdi2FN45HZZ+7u8oi23r3nTmr99NXxsoNxWEt8mjKkdOpc1rIxDvm3B8r0y7rVTdO+aPL4Pqjdz4RqPTS3atzkNh/+ag/MPptcxUGH8thIXBZR9CZJ6iBW0vgcE/QkJTd86UOyErzX+xFtOXnoWBfU+yLitTMXibPP4Ahogc+Xba83iI5/mAVGbdVpkHYeINrq++zcUwcAZUK8VCxGY9IbEZzHnTJaVQowpT74JtS+DiEVC94dHnJDR238fTbof5H7jW3PZl5R9rCUxbtJWWtavQtGZlf4diTLnwNnk8BJwDLBeRZ0TkThF5Gje43dvzfEAqk5bH7vXw/jmw5DPo+YTrcokq5WrisDA4dTBsXggb5/okzICyYCT8MRa63Q/Ni1jTEhEFvf8NV38CB3fA8LPcWFAAJ9Tt+9KZt35X4AyU52S7ltvntwb0980EN2/3MP8OaA/8hhvfeApXlHAh0EFVvy+zCEvJ5y2PtbNcYcM9G+CqiXD6P31XWuTkK926kF/e8s31AsWmBTD9fte66na/d69p1hMG/wwNToUv7nCtvPRymDFXAl8u2YoqgVEIMTcHPr/Ftdx++y8s+8LfEZkQ5U1J9mgReRiIVNUBntlQlTyfr/S2plXQU4Vf3oEPL4S4mjDo+6LfQZdEVJyrc7VsclDPOjrCwZ0w/jqofBxc/J5rYXmrSm1XzuXsx2DpZDeYvnFe2cVaQtMWbaF5rco0r13Fv4Hk5sLk21yxzbMehtpt4csHXRerMT7mTUn2DOBhoFqZRxOostLh86Hw5f3Q4hz4v28gsWnZ3KvTIEBg7vCyuX55ys2BSf8HB7bBZaNdMcjiCguDM+5yYyEAI/vATy+7P5TlIW03/DneDeK/chIsGH3E0zv2ZzB37S7/b/qUmwtT/gm/j4HuD0K3++C8F2FfMvzwnH9jMyHJ27eBvwIdyjKQslLqMY99W2DUee6Xstv9cPkYiCnDyWVV60HrC2DBh5BxoOzuUx5+eBbWfOem5NYt5Y9P/VPg5h+h1fnwzePw0cWuZlhZ2JsMvw6H0f3h+WYw6SbY8IvrUvzidvj+P4fHEmYs2Uquv7usVGHaPbBwNJxxz99dgw1OhZOvhjlvuPI3pvRU4c8J1h0IiHoxoCYip+DKnr9K4SXZA3pXo6SkJJ0/f37xXzjuKljzPVz0DrTu7/vACrJxHrzf0y027HRT+dzT11Z+BR9fCidfBRe86btxIVX3R3L6/RBdxf2/lLaopCpsX+oWLy6fClt+d8drtIBW57mEVacDaI57d//bR+7r6vcqV41cwJa96Xx7VzfEH2X1Vd33Yu670PUON4EjbxwHdsAbHeH4dnDtZCv9Xxrpe93425JPQcLdmq7GZ/g7qjInIgtUNemo414mj7x9BAW+QFXDSx5e2Stx8ti3xS3+q32i74Mqyntnuy6TW+cXb5wgEOxeB+92g6r14cavSj8TrSDbl8GEG2DHMje1t8cjbqaWt3JzYOOvnoQxxcUMUO8UlzBangc1Wxz9OlXXopr5NJmNutNxxTVc170t95zjh6KWqm6fmDlvuHVC5zxVcHKYNwKm3u1mBba9pPzjDAWbFsDEgbBno2vZLZoA6Xtg0A9QtdDdI0JCaZPHdcc6R1VHH+scfypx8vCXRRPhkxvhiv9Byz7HPj9QZKXDB71h1zq4eSYkNCnDe6XBjIfczKK6Hd0fx4TGRZ+/5nuXMFZOh9SdEB4Fjbt5Esa5x65FdsjC/5L7xR0szalP1LUTadGsuW++Jm+puu67n1+BTjfDuc8W3qrIzYERZ7s3QrfNdy024x1VN/vx68egcm245H1o0Nl1A77XA2qdANdP9b4UURAqcfIQkWjgHmCKqgbsYsDC5NnP46ZVq1b5Oxzv5WTBq+3cwsHrfFJguHxMvg0WfggDxkKrvuVzz6Wfw+e3AQrnv3zku+vUXbByhmtdrPkOslIhuiq06A0t+7ourxKOYT33xhvclvJvYqrWRK7+pHxL6n/3JMx6HpIGwnkvHbs7atMC15rt4mmhmGM7uNNNe175pWuJXvDGkZM+lnwGE66DpBvh/Jf8FmZZK23L4yDQtxz29SgzQdfyADer6JvHYcjs8u82K4mF/4XJt8Lpd0HPx8r33rvXwyf/B8lzXaHF2m1cC2P9bDdWUaWOS2atzoOGpxevi6ug2x3MJOmpb3i0QwbXrb0PcjLhinHQsIuPvqAizHwWZv4HOlwL57/qfbfmF3e4/6PBPwbWz9O2Je5NR6vz4dSb3ZR1f1v3s/t5Sk2B3k+6WZAFJeivHoHZr8GFb7t1WiGosOThbWf6XIJ0tlVQ63AdRFYKjkWDW/5wM34ad4Mew8r//tUbwg3T4Iy74bcx8OUDcDAFTr8Tbvoe7lrqpq427VHqxAHw1dKt5OQqHbv0gP/72q39+fCCsq9N9uOLLnG0u7J4iQPcepmYqjD1nsBZeZ66C8Ze4YpifvsEvHqym+mWneGfeHJz3NTm0ee72ms3fu0SWmEtu7Mfg0ZnuBX9W4KuY6ZUvP3Juw+4RURuFZEmnv09KuX9KMsgK6xKCdDuCjc18MAOf0dTuLTd8L9rIDbBjTuE+WnuRHgknP0oDJ0Lty2Eob/A2Y+4acI+nmU0bdFWGiRU4sQ68a4o5o1fQZ2T3YLIstrY6+fXXC21tpe5LpTiTqSolAC9noANs91CQn/LyXaVA/ZvcXvfDJwBNZrD9Hvh9Y5uVlt57vGyfyv890K310ybS+DmWe7/tCjhEXDJSKiUCP+7ukLtyVOcdR5NgdeAVcA+YH++D1MWTh0MORmuNlQgys2FSTfDvs2uUm7lmv6OyM2SKqtFnMDe1Cx+Xp3CuW2P+3t6bqUEN3Wz1Xmu1TPjYd8uZJzzltvv/sSLXRdJSRP0yVe7GWVfDYO0Pb6LryS+eQz+munGqeoluYHo66fC1ZPcJmmfD4W3OsPiSWW/KHTVN/B2V0ie76aWXzzc+4kFlWvCZf91yeeT/3OtlwrA2+QxELjB8zGwkA9TFmq2gGa93F4f/mrKF+WnF2HVDDjnP24hXwXw1dKtZOcqfdvkW1UeGesSaKeb3fTZTwa62WelNfc9mPEgnNDf/VELjyj5tcLCXPdd6k73Dttf/hjnvkedBrkxqkNEoNnZrqvx8o8gLMK1Tt4900188HV3W04WfP0ojPkHVK7ltlVof3XxW6r1OrrFsGu+hZlP+zbGQKWqIf0B9AOGN2vWTIPWqm9UH4tX/e1jf0dypNXfqj5WVXXijaq5uf6OptzcMHKunvb0t5pb2Necm6v606vu/+yDc1VTd5X8ZvPed9f5+ArV7MySXye/qfeoPl5NdfPvvrumt5IXqP6rpuoHfY/9NeVkq/4+TvWVk9z3YUQv1b9m+SaOXetU3zvbXXfyHaqZqaW7Xm6u6me3uOstm+KTEAMBMF8L+NtarE5TEWktIteIyEMicpznWDMRCdiJ41pe29CWpaY9oGYrN3AeKAOdezbCxBtdXP1erTArl/emZfHjqh30zdtllZ8IdL3djf8kz4MP+rjvV3Et/NANxDY/By4d6cZ0fOWsh10//dS7y69OGMCB7W5soHItV+/sWF9TWDi0u9wtlj3/Zfd9HH2+K1C6aUHJ41g6Gd49w63XuGQk9HvFtRxLQwT6vuh2yfx0MKSsLt31ApxXyUNEKovIeGAxMAL4N1DH8/R/gHKel1nBiEDnIbD1T1j/s7+jcd1nE65zTf7L/xsYUyvLybfLtpGVo97t3dH2Etd/v28LjOgJW/70/ka/fwyTb4emZ7uuMF8vQoutBr3+7ZLb7x/59tqFyc50BSZTd8GAMW5cw1vhkW5Ny+0LofdT7nfhvR6ufNC2YhT2zkp3CXP8NW7Xyptnuc3cfCUyxv1OhEW4JBns9emK4G3L4yXgNOBsoAqQ9y3XNCCIlkAHqZMud7OZymomT3HMeMi967vwTTc7pgKZtmgLdarGcHL9at69oPEZMPBL9w56ZF+3UPFY/pwAn90CTbq5P7J5t+v1pXYDoEEXt3q6PGYJffkAbJjjZood365k14iMhdNuhTv+cK2ntbPg7dPgk5tg55qiX5uyyiXxeSOgy61udldRFQlKqloDuOQDSFnh1q8ESm+Bj3mbPC4G7le36VP+qQTrgQL2FDU+FRnr3nktnwq71vovjj/+5375TrvNVf+tQPanZzFrZQrntj2+eEUQa7d2ZfyrNYAxl7rB4sIsngSfDoJGp7tV+qXtSimKiBs8T9/rpgCXpQWjYP77rg6ZL+prRVdxZefv+MMVhFz2BbzZyS2E3Lvp6PP/GOfqre3bBFeOd6vsfbDep1BNz3L11pZMCo51WiXgbfKIBXYW8lwVjk4opiyc8n/uHay/9vrYtsT9cjbsCmc/7p8Y/Oi75dvJzMktWfn1+DowcDo0PA0+vRlmvXD0O9Klk91Uz/qnutXqZVFQMr/aJ7rp4AtGlW4MoSgbfnULE5v2gJ6P+/bah9au3PG7e3P12xh4rT18+ZBbJJpxAD4d4r7ndU6GwT+5PXnKw+l3ulXzXz0C634qn3uWI2+Txzzg2kKeuwSY7ZtwTJHij3fz/Bf+F9L3le+90/e6PtyYeDfAWJrpokFq6p9bOC4+hvb1q5fsAjFV4apP3CK/7/4NU+/6exHc8mluSmrdjnDVBIiu7LvAj6X7A67o35S7fL9GYd9mN75QtZ7ryimrBaRVjoO+z8NtC6DtpfDr227zrrdPgz/Gukq4104u3wq4Im5NTkITmHC9+16EEG+TxyPAxSLyDfB/uLLsfUXkv7g9zct9wNyzyn2+iJxf3vf2q85DIHO/W31bXlRdH/zu9XDpKLc9bAVzICObH1buoE+b4wgLK8XMsogouOhd9650/gcuIS/5zA0kH98Orp5Y/lVvY+JdN86W310LxFey0v8eNL5iLMSWMOkWR/WGbizull9d8cuwcFdY9KyH/POGJyberVfJTHX/x9mZ5R9DGfEqeajqj7jB8mjgDdyA+RNAE6Cnqnq9sbSIfCAi20Vkcb7jfURkhYisFpEHvLjU/cB4b+8bMup2cIOcv75TPitZU3e50vDLp0Dvf7tulwpo+A9ryMjO5eIOPnjnGhbmum/6vuAWWE64zo2LXD3JtU78oc0/XI2mb//luntKS9W1rDYtgIvfdaXLy1PNFu6Nzu2/QeMzy/fe+dVq5RJa8jy32DNEeL3OQ1V/VtUzgHigHlBFVbuqanHnjo4i3+wsEQkH3gTOBVoDV3jWlLQVkSn5PmqJSC9gKbC9mPcODZ2HwJ71sGJa2d5n5QxXHmLp53DWMOh8S9neL0Al707l3Vl/0b9dHU6qV813F+50Ewz42M2ku+YzN33WXw4NnmcedLOvSmvu8L+3bj6hX+mvF+xOvMhNMpk3An4f6+9ofKLY7ThVTQPSSnpDVZ0lIo3yHe4ErFbVvwBEZBxwgao+DRzVLSUi3YE4XKJJE5FpqnrUSicRGQQMAmjQoEFJQw48Lc+Dqg3ctN2y+MVM3+um4/72EdRq7frgSzq1MgQ8PX05IvDAua18f/GW57qPQFCzpdvv4+dXoMM1rtZUSaydBV8+6PZL6eZNJ0IFcfbjsPl3t5Vx7dZB/zsVKPub1gXyLsFN9hwrkKo+rKr/xO2r/l5BicNz3nBVTVLVpJo1A6Bgn6+ER7gy0et/dj+MvvTXTHjrNLdI7fS7XK2fIP8hL425a3cx9c8tDO7WlDrVynDabKA4816Ir+tmR5Wkou3u9a6ycGIzN7YTbFsol6UjKvBeE/QVeIP6f1ZVR6nqlKLOEZF+IjJ879695RVW+ehwDURV9t2iwcyD7g/Ghxe4RWkDv3IbOoXw9prHkpOrPPHFEupUjeHmM8uuSm9Aia4MfZ6GbYtcF0txZKa6Fd+5Oa47roQ7NIa0wxV4t8Ckm4K6Am+gJI9NQP08j+t5jpnCxFSFk6+CxZ+4UtClsX6OK0c97z03rnHzjxWmQm5RJi7YyJLN+7j/3FbERvlpjxJ/OKG/K4vy/VPe/2ypuhLq2xbDP0ZAjWZlG2MwO1SBd/U3MPMZf0dTYoGSPOYBzUWksYhEAQMAn2zcHRKFEQtz6s2Qmw3z3i/Z67PS3b4OI891W7VeN8W96yyPxWkBbn96Fs/PWEFSw+r0b1fn2C8IJSJuzUR2uitX7o2fX3Wrqc9+1E2RNUXreL0r/T7rOVgx3d/RlEi5Jw8RGQvMAVqKSLKI3Kiq2cCtwAxgGTBeVZf46H6h2W0FbsOjlue6sg/F3Tdi0wK3R8Ls190P8pDZrg6TAeCN71aTciCTR/u1Ll4pklCR2NSV/fjzf8deHb3qG/jmcTej6PQ7yyW8oJe3Au+kQceuyxWAyj15qOoVqnq8qkaqaj1Vfd9zfJqqtlDVpqrqs11qQrrlAW7abupOWOTlkpfsTPjuKRjRCzL2w9WfuHLU5b0wLYCtSznIBz+v5ZKO9Xw7NTfYnH6Xq8c19R5XQbkgO9fAxIFQu43bga8iJtqSyl+BN/OgvyMqlkDptjIl1egM94v7y9vHrt65dTGM6OGayiddBrfMgWY9yyfOIPLUtGVEhYdx3zkt/R2Kf0VVgj7Pwo5lBU/MyNgPY69wq7gHjKlQpfl9ploDuOR92LE86CrwhnzyCOluK/Ds9XELbF8Ka38o+JycbPjxRRje3Q2ADvgYLnrHv4vSAtRPq1L4euk2hvZoRq34MiqFHkxa9YUWfdzAbt5qtYf2rt+52q3krm6FtUusaQ/oMcxNfvnf1fDXD0GRREI+eYR8txW40hJxNWFOAaWfd6yED3q7shOt+rqaP63OK/8Yg0B2Ti7/mrKE+gmxDOxaBvs8BKtzn3UTKr56+O9js56DFVNdTawm3fwXW6g4/S63Gn/9z/Bhf3jjFPf7nLbb35EVKuSTR8i3PMD1nSbd6OokHdr6MjcX5rzpttrc9ZfbEvXS0RCX6N9YA9jYuRtYue0AD/dtTUxkBZqaeyzVG7k/bks+hTXfw7IpMPNpaHelK+duSk/EFW+8axlc6OkVmPEgvHgCfDYUkhcEXGtENMACKitJSUk6f/58f4dRdg5sh5dPhA7XuRITnw9172Ja9HF7jFcpwR4UFcie1EzOemEmrY6L5+ObTq2YM6yKkpXu6pxprpugUaMF3DC97HY5NG7b4vkfwJ/jIeugq/SQdKPbTKscx5dEZIGqJh113JJHCPnsFrcTnYS5j3OfcQsJ7Q/hMT0+eQkfzlnH1NvP4ITjbWV0gVZ9A2P+AXG1XNma8twboyJL3+emTM//wI1tRld1WwgnDXQVe8uYJY+KkDy2Lobh3dwWpv3fgGr1j/0aw+rt+znnlR+5/JT6/Oeitv4OJ7D99hHUae92IDTlSxU2/OKSyNLPICfT7eqZNNAVSC2jUkIVNnmISD+gX7NmzW5atWqVv8Mpewd3uq05rbXhFVXlupHz+G3Dbmbe053EyhW3lpcJIgdTXCJfMBJ2r4NKNVy9u443+HzmW4VNHodUiJaHKbbvlm9j4Kj5DDvvBP7vjCb+DseY4snNhTXfudbIyumuddK8lxsbad7LJ9v+FpY8Kt5G1MZ4ZGbn8uSUZTSpEce1XRr5Oxxjii8sDJr3dB97k2HBaFj4IYy9HKrWh47XQftry2Tr6JCfqmtMYT6cs46/Ug7yyPmtiYqwXwUT5KrWgx4Pw52L4bIPIaEJfPckvNzatU58LORbHnnGPPwdigkgOw9k8Oq3q+jWoiZntarl73CM8Z3wSGh9gftIWQULR0O9Tj6/Tci/3aoQK8xNsb349UpSM3N45PwT/B2KMWWnRnPo/aTb5MvHQj55GJPfsi37GDd3A9d0bkizWlZN2JiSsORhKhRV5V9fLKVqbCR39mzh73CMCVqWPEyFMmPJVub8tZO7erWgaqVIf4djTNCy5GEqjPSsHJ6atoyWtatwRacG/g7HmKAW8smjQlTVNV754Oe1bNyVxqP9WhMRHvI/+saUqZD/DQr12VaqSnZOLulZORzIyGZvWhZpmTnk5laMygHe2r4vnTe+W02v1rXp2qyGv8MxJuiF/DqP0tqxP4O9aVlk5eSSlZNLZnYumTm5ZOUomdm5h49nHPp3vuczPceycnLJzPearJxcsnOV7BwlOzfX8znfv484xz2Xk6Nk5TmnMJHhQkxEONGRYUTn/RwRRswR/3afDz1f0HNx0RHUqx5Lg4RKJMRFBV3J8udmrCArJ5eH+9rUXGN8wZLHMTw9fRmTFm469omFiAwXIsPDiIoIc589/44Ic8cjw4XwMCHCc7xSuHsuIkyICBciwsI8n9057rkCjnkeh4mQnatkZOWSnp1zxOeM7BzSPZ8zsnPZk5pJRrZLfOlZ7lhGVg7p2bnkFJGUKkWF0yChEvUTKlG/eiUaJMTSILESDRIqUa96pYDbSOnP5D1MXJDMzd2a0KiG7bNtjC9Y8jiGq05tQLcWNYn2/PHPnwgiI8R99hx3x1xSiAoPC7p36Idk5xyZVPanZ5O8O5UNu9zHxl2pbNiZyk+rUkjLyjnitbWqRNMgwZNMPJ9dsomldpUYwsLK73uiqjzxxVJqVI7m1rOsyoAxvmLJ4xg6Nkygo28rHAeFiPAwIsJdd9UhLY87ekGdqpJyIJMNu1JdctnpSS67U/l17S4+/X3TEbtnRkWEUa96rKfFUokWx1XhtKaJNKkRVyaJdvIfm1mwfjfP/qMtVWJsaq4xvmLJw5SKiFCzSjQ1q0TTsWH1o57PzM5l054011LxtFg2elowv23Yzb70bACOi4/htKaJnNasBqc1TaROtdhSx5aWmcMz05dzYp14LuloG2MZ40tBmTxEpDvwb2AJME5VZ/ozHlO4qIgwGteIo3EBYw2qyoZdqfy8eiez16Tww8odTPrNjS81rhHnkknTGnRpmkhCXFSx7/3urDVs2ZvOqwPaE16OXWXGVATlnjxE5APgfGC7qrbJc7wP8CoQDoxQ1WeKuIwCB4AYILkMwzVlSERomBhHw8Q4rjy1Abm5yopt+5m9ZiezV6fw+e+bGfPrBgBaHx/PaU0T6dqsBqc0TqBydNE/upv3pPHOD2s476Tj6dQ4oTy+HGMqlHLfSVBEzsT94f/wUPIQkXBgJdALlwzmAVfgEsnT+S4xEEhR1VwRqQ28pKpXHeu+tpNg8MnOyeXPTXuZvTqF2Wt2Mn/9bjKzc4kIE9rVr0bXpol0aVqDDg2rER1x5Ayv28f+xowlW/n27m7Uq17JT1+BMcEvYHYSVNVZItIo3+FOwGpV/QtARMYBF6jq07hWSmF2A4VuOi0ig4BBAA0aWDmKYBMRHkaHBtXp0KA6t/ZoTnpWDgvW72b2mhR+Xr2TN75fzWvfrSYmMoxTGiXQpWkiXZvWICM7l8l/bOb2Hs0scRhTRgJlzKMusDHP42Tg1MJOFpGLgXOAasAbhZ2nqsOB4eBaHr4I1PhPTGQ4XZvVoGuzGtx7DuxLz2LuX7v4eU0Ks1fv5LkvVwArADcAP7h7U/8GbEwIC5TkUSyqOgmY5M25tpNg6IqPiaRn69r0bO32Z96xP4M5f+1k7tqdnNvmeCpFBeWPtzFBIVB+uzYBeedS1vMcM8ZrNatE079dHfq3q+PvUIwJeYFSGHEe0FxEGotIFDAAmOyLC4d6YURjjPGHck8eIjIWmAO0FJFkEblRVbOBW4EZwDJgvKou8dH9rCS7Mcb4WLlP1fUXm6prjDHFV9hU3UDptjLGGBNEQj55WLeVMcb4XsgnDxswN8YY3wv55GEtD2OM8b2QTx7W8jDGGN+rMLOtRGQHsN7fceRTA0jxdxBeCqZYIbjiDaZYIbjiDaZYITDjbaiqNfMfrDDJIxCJyPyCpsAFomCKFYIr3mCKFYIr3mCKFYIr3pDvtjLGGON7ljyMMcYUmyUP/xru7wCKIZhiheCKN5hiheCKN5hihSCK18Y8jDHGFJu1PIwxxhSbJY9yJiL1ReR7EVkqIktE5A5/x+QNEQkXkd9EZIq/YymKiFQTkYkislxElolIF3/HVBQRudPzc7BYRMaKSIy/Y8pLRD4Qke0isjjPsQQR+VpEVnk+V/dnjIcUEuvznp+FP0XkUxGp5scQj1BQvHmeu1tEVERq+CM2b1jyKH/ZwN2q2hroDAwVkdZ+jskbd+DK5Qe6V4EvVbUV0I4AjllE6gK3A0mq2gYIx+1lE0hGAX3yHXsA+FZVmwPfeh4HglEcHevXQBtVPQlYCTxY3kEVYRRHx4uI1Ad6AxvKO6DisORRzlR1i6ou9Px7P+6PW13/RlU0EakHnAeM8HcsRRGRqsCZwPsAqpqpqnv8GtSxRQCxIhIBVAI2+zmeI6jqLGBXvsMXAKM9/x4NXFieMRWmoFhV9SvPfkEAv+B2KQ0IhXxvAV4G7gMCekDakocfiUgjoD3wq59DOZZXcD/MuX6O41gaAzuAkZ4uthEiEufvoAqjqpuAF3DvMLcAe1X1K/9G5ZXaqrrF8++tQG1/BlMMA4Hp/g6iKCJyAbBJVf/wdyzHYsnDT0SkMvAJ8E9V3efveAojIucD21V1gb9j8UIE0AF4W1XbAwcJnC6Vo3jGCi7AJb06QJyIXO3fqIpH3XTNgH6HDCAiD+O6jMf4O5bCiEgl4CHgUX/H4g1LHn4gIpG4xDFGVSf5O55j6Ar0F5F1wDigh4h85N+QCpUMJKvqoZbcRFwyCVQ9gbWqukNVs4BJwGl+jskb20TkeADP5+1+jqdIInI9cD5wlQb22oSmuDcSf3h+3+oBC0XkOL9GVQhLHuVMRATXJ79MVV/ydzzHoqoPqmo9VW2EG8z9TlUD8t2xqm4FNopIS8+hs4GlfgzpWDYAnUWkkufn4mwCeIA/j8nAdZ5/Xwd87sdYiiQifXBdrv1VNdXf8RRFVRepai1VbeT5fUsGOnh+rgOOJY/y1xW4BvcO/nfPR19/BxVCbgPGiMifwMnAf/wbTuE8LaSJwEJgEe73MaBWGIvIWGAO0FJEkkXkRuAZoJeIrMK1np7xZ4yHFBLrG0AV4GvP79o7fg0yj0LiDRq2wtwYY0yxWcvDGGNMsVnyMMYYU2yWPIwxxhSbJQ9jjDHFZsnDGGNMsVnyMBWeiLTxVDDt7nmsInKrf6MqHhG53hN3ZX/HYiqGCH8HYEwA6gKs9XcQxTQVF3dAL4QzocOShzH5qOov/o6huFR1B64opDHlwrqtTIUjIreIyEYROSgiXwDH53v+iG4rEZnp2WDqBhFZKyIHROS/IhItIp1EZK7n2EwRaZDvWjEi8pznfhki8kf+igIisk5EXvBsDJUsIrtFZFzejYtEJNJzzgbPdTZ7NjeK8jx/VLeViNQQkdEislNEUj3xJRX33sYUxFoepkLxlLx+E3gH+AzoBnzgxUs7AzVw5U8a4PZcSANOBZ7DVfB9DVdeJO8GPxOBTsBjwBrgMmCyiCSp6u95zrsM+BMYhCuI9xKutMotnucfBK7CVQleCxwH9MVtIFWYz4BmwD1ACnAv8L2ItFfV1cW4tzFHU1X7sI8K8wHMBabnO/Yerqx4d89jBW7N8/xMYA9QNc+x8Z7zzsxz7BbPsUqex2d7HnfLd79ZwIQ8j9fhEktEnmOvAFvzPJ4CvFjE13W9516VPY/75L83EIfr2nq3OPe2D/so6MO6rUyF4dmtrwNHV4H1piz+fFXdm+fxaiAT+CnfMXB7c4ArGrgV+FlEIg594LZuPaL7CPhe/97xDlw14Fqe8v0AvwPXi8h9InKSpwpvUTrh9mH54dABVT2IS0KnF/PexhzFuq1MRVID182Tf/8Jb/aj2JPvcSawX1Vz8x0DiMlzv+OArAKul+PF9QWI9rz+SdxOjrcAzwKbROR5VX21kHgL22djG5BQzHsbcxRLHqYiScH90a6V73j+x76yC9iED/b4VtV03A5zj4pIc2Aw8IqIrFDVLwt4yRYK/rpqU/C+2cYUi3VbmQrD0zXzG27r17wuLqNbfotreRxQ1fn5P0p6UVVdhRsEzwBaF3Lar7iupzMPHfBsc3oeR3a1GVMi1vIwFc1/gEki8jbwKW62VZ+iX1JiXwMzcBsRPQssAeJxm1TFqOqD3l5IRD4FFuCSXxpwCe73d1ZB56vqDBGZDfxPRB4AduISTizwfEm/IGMOsZaHqVBU9VPcdNt+uKms7YEy2cFNVRXXqvkA+CcukbyLWwle3Hf/s3HdXx/jBvw7Av84RgvmQlwCewWYgBvH6KFHTtM1pkRsJ0FjjDHFZi0PY4wxxWbJwxhjTLFZ8jDGGFNsljyMMcYUmyUPY4wxxWbJwxhjTLFZ8jDGGFNsljyMMcYUmyUPY4wxxfb/CQGciEgBNJIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(np.arange(1,D+1),np.divide(np.abs(m_riemann-np.arange(1,D+1)),np.arange(1,D+1)),label='Riemann')\n", "plt.plot(np.arange(1,D+1),np.divide(np.abs(m_sample-np.arange(1,D+1)),np.arange(1,D+1)),label='Monte Carlo')\n", "#plt.plot(np.arange(1,D+1),N**(-1./np.arange(1,D+1)))\n", "plt.xlabel('dimension',size=15)\n", "plt.ylabel('error on integral',size=15)\n", "plt.legend(fontsize=15)\n", "plt.yscale('log')\n" ] }, { "cell_type": "code", "execution_count": null, "id": "equipped-incident", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "dangerous-yeast", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.7" } }, "nbformat": 4, "nbformat_minor": 5 }