docs/source/examples/Metropolis_algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Metropolis-Hastings Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In statistics and statistical physics, the Metropolis–Hastings algorithm is a Markov chain Monte Carlo (MCMC) method for obtaining a sequence of random samples from a probability distribution from which direct sampling is difficult. \n",
    "This sequence can be used to approximate the distribution (e.g. to generate a histogram) or to compute an integral (e.g. an expected value). \n",
    "Metropolis–Hastings and other MCMC algorithms are generally used for sampling from multi-dimensional distributions, especially when the number of dimensions is high."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Metropolis–Hastings algorithm can draw samples from any probability distribution $P(x)$, provided that we know a function $f(x)$ proportional to the density of $P$ and the values of $f(x)$ can be calculated. \n",
    "The requirement that $f(x)$ must only be proportional to the density, rather than exactly equal to it, makes the Metropolis–Hastings algorithm particularly useful, because calculating the necessary normalization factor is often extremely difficult in practice.\n",
    "\n",
    "The Metropolis–Hastings algorithm works by generating a sequence of sample values in such a way that, as more and more sample values are produced, the distribution of values more closely approximates the desired distribution $P(x)$. \n",
    "These sample values are produced iteratively, with the distribution of the next sample being dependent only on the current sample value (thus making the sequence of samples into a Markov chain). \n",
    "Specifically, at each iteration, the algorithm picks a candidate for the next sample value based on the current sample value. \n",
    "Then, with some probability, the candidate is either accepted (in which case the candidate value is used in the next iteration) or rejected (in which case the candidate value is discarded, and current value is reused in the next iteration)—the probability of acceptance is determined by comparing the values of the function $f(x)$ of the current and candidate sample values with respect to the desired distribution $P(x)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Monte Carlo methods vary, but tend to follow a particular pattern:\n",
    "\n",
    "* Define a domain of possible inputs\n",
    "* Generate inputs randomly from a probability distribution over the domain\n",
    "* Perform a deterministic computation on the inputs\n",
    "* Aggregate the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Intuition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First of all we need to include some useful (fundamental!) libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np # utility structs and math\n",
    "import pylab as plt # result visualization\n",
    "from time import time as now # time monitoring"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, consider a unitary circle inscribed in a unit square. \n",
    "\n",
    "The value of $\\pi$ can be approximated using a Monte Carlo method:\n",
    "\n",
    "* Draw a square, then inscribe a circle within it\n",
    "* Uniformly scatter a given number of points over the square\n",
    "* Count the number of points inside the circle, i.e. having a distance from the origin of less than $1$\n",
    "* The ratio of the inside-count and the total-sample-count is an estimate of the ratio of the two areas, $\\frac{\\pi}{4}$. \n",
    "* Multiply the result by 4 to estimate $\\pi$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets start from a simple square domain $[-1, 1]$ in both the axes.\n",
    "\n",
    "**Note:** In this way we can consider a unitary circle included in it, i.e\n",
    "\n",
    "$$\n",
    "r=1 \\rightarrow A = \\pi r^2 = \\pi\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUQAAAEvCAYAAAA92bhfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3daWyc13kv8P/DGa7iJoqUuGi35MSKV5KWZbtNFVvyEpPyTZHCDoo0LRpRKZxb4AK5idEAplzg4mZBPjRt0kgKgtv0Q9J8cULKsh1bge0kjWWTsmXLq2RLtsShJIqUKJLiMhye++GZ4SIOyVnOzDtz5v8DBiRnRuc9oyGfOfsjxhgQERGQ53UFiIgyBQMiEVEYAyIRURgDIhFRGAMiEVEYAyIRUZjf6wosprq62qxfv97rahCRY7q7uy8aY2quvT+jA+L69evR1dXldTWIyDEi8nG0+9llJiIKY0AkIgpjQCQiCmNAJCIKY0AkIgpjQCQiCmNAJCIKsxIQReRnInJBRI4v8LiIyA9F5KSIvCkijTauS0Rkk60W4v8D8MAijz8IYHP41gbg3y1dl4jIGis7VYwxL4vI+kWe8jCAnxs9nvsVEakUkTpjTK+N61OOmJrSmzHzbwAgMv+Wl6c3ohik6zelAcCZWT+fDd83j4i0iUiXiHT19fXFdZH169dDRFBfXw8A2Lt3L0Rk+tbd3Y3u7u459+3duxcAUF9fP31fU1MTAKCtrW3OcwOBADo7O+fct3///ki9p2+tra0AgNbW1jn3A8D+/fvn3NfZ2YlAIDDnvra2NgBAU1PT9H0uv6Z/fuIJlIigUgQrRfDmoUN466mncIsI7hbBThH84pFHgF/9Cv9QXY1HCwrwSGEhHv/0p4Ff/xr7H3oIj5SU4JHiYjxSVISBn/4U3f/0T3gkPx+P+nz4ggie/sd/BH7/e9wmgutFsEYEX77vPuDSJXzxwQeRz/cpa1/TeovnHYitnCrhFuJBY8yNUR57GsD/Ncb8IfzzYQDfNMZ0L1Zmc3OziWcvc1NTE7q7Fy2SvDQ+DgwPA1euAH19QH8/MDICTE7Of64IkJ8P+P36NT8f8Pn0/nhNTgLB4NyvU1Pzn+fzASUlQFUVUFMDVFQApaVAcXH816S0ERHEG8dEpNsY03zt/ek63OEsgDWzfl4NIGD7IgyGGWRsDBga0uB34YLerl6dCWhFRXpbvlwDUSr5/XpbSiikwfLCBeDj8N5/Y4DCQg2Qq1ZpkCwr0yCZSHCmjJaugNgB4Osi8ksAdwAYTMX4YVtb23STm9IoFAIuX9YW3/nzwMWLGhAjAaO4WG8VFd7Wcyk+n94KC+fWNRgELl0CentnWpZ+vwbJlSuBFSu0VRlL0CXrOjo6rJVlpcssIr8AsB1ANYDzANoB5AOAMeYnoh3+f4PORF8F8HfGmCX7wvF2mRNpOlOCxsaAgQHgzBngk0+0ZZWXNxP88vO9rmFqhULA6Ki2eqemNPjX1wPr1gHV1dr1prQIBALT45yxSmmX2RjzpSUeNwAes3Et8ogxM2N/p05pKxDQbm9VVeq7vZnG59PxxdJS/TkU0g+Is2f158pKYMOGmW42Z7pTpqGhwVpDiG18WtjkpP6R9/YCp09riwjQIFBbyzG02Xw+DXyRrvboKHDsmLYeCwq05bh6tX54FBR4W1dakFMBsaenx+sqZD9jdLzso4+ADz/Un/1+oLw888cAM0lk6ADQD5ZPPgFOnNAPkbVrgU2btGvNlmNGcSogdnd3xz2WQGETE9oSfPddbRUWFOgfbK51hVPB79eWIaAtxvPntcVdWgrccIO2HLm0J2G7d++2Vpa1dYipwEmVNLh0Sf84T5zQcbDycmDZMq9rlRvGxnR2HtDxxuuu0xlrDkWknNfrECmTBIMzrcH+fm0NrljB1mC6FRXpWOzUFNDTo0MU5eXaamxo0MdpSTY3ZDAg5pIrV3SG+IMPdFyrrEyXipC38vJmutSjo8Crr+p9GzboWGPkMYrq6NGj1spyKiDu27fP6ypkpuFh4J13gJMndX3g8uVcRJypIpMxoZAu4Tl5UscYb7pJ3zdKKaf+KiIb0yns6lXg/feB997TAFhby1nNbOHzacvQGF3z+cwzwMaNwJYt2q2maXV1ddbKciogclIlbGxMu8XvvKN/WCtXcnwwW4nMBMZAQJdDbd6s44yRReE5LhCwdywCmwsuGR/XINjRoRMmNTV6YzDMfpHAWFurqwI6O4HXX9deQI6LHDlmA5fduCAY1MmSyM6IFSs4Rui6yFbBUEi70ddfn7Oz0tl4/FdatLS0eF2F9DJGj6k6elRbhytWcFtYrvD5tPUfCmlv4L33gJtv1u40ewQJcyogdnZ2el2F9BkaArq7df1adTWXZuQqn08PkAgGtQt96hRwxx38fUiQU2OIkaPGnTY1pbtKDh3SLlNDg57fR7ktPx+oq9PA+OyzwJtv6vc5IJ5htaU41UI8ePCg11VIrcFBXbTb16fdJdfPHKT4RbZevv22Dqds26a/KxQTp1qIzgqFdPb40CGdVayvZzCkhfl82loUAX77Wx1jnpjwulYp09w8b24kYU61EJ00MAAcOaKHAHA9IcWjtFRP7j5xQo8fu+MODZS0IKdaiE4tuQkGdRzomWd033FdHYMhxS8vTyddCguB3/0OeOUVXbhPUTkVEJ1JMDU4qF2dd97RhbhlZV7XiLJdcbEOtZw5Azz9tGYWdER7e7u1srgwO9OcPQv88Y/6C8w9q5QKV6/qOZhbt+ppOjl4/uJCC7OdaiFmtakp7SK/9JImKGIwpFQpKdFu9JEjumohy5fn2Dwln5MqmWB8XH85z5zRLjLHCinV/H5dw3r6tA7R3HVX1h4W0dtrL8W7Uy1Emwmr0+byZR0vPH9ef0EZDCldRLSlODysi7kdGldMlFMBsampyesqxOfsWeC553RPMhfPkleqqnQx9/PP67FxWTYO39jYaK0sp7rMNhNWp9TUFHD8OPDWW7oPmVvvyGvFxdpafPVVXfva1JQ1i/9t5VMBHGshZoWxMeD3v9eAWFvLYEiZIzKu+PHHwOHD2pXOAjZPymdATKerV3Vx7IULuiaM44WUaSLjiqOjOrYdSZOawQ4cOGCtLKcCos2E1dYNDQEvvKAtRI4XUqZbvly7zM8/r6lqc4RTATFjd6oMDmownJriOXWUPcrKdM3iCy/kzAy0UwExI2eZBwb0U9bn0wXXRNlk2TLdJHD4sCa5ykA9PT3WynIqINpMWG1Ff79+uhYVcecJZa/iYu3ZvPiiLhXLMJxlzgaRYLhsWdbuACCaVlioS8ReeinjguKuXbusleVUQLSZsDopAwPaxSgt1YBI5IKCAg2KL7+suXwc5FRAtJmwOmEDAzMtQwZDck1hoWZ3fOklJ4OiUwHRZsLqhFy6pMGwpITBkNw1OyhaPFghUfv27bNWFs9DtGVkRBey5udzzJByw9iYLty+/35dt5hFeB5iKgWDeqirMQyGlDuKivT3/eWXdReWR8TiAbcMiMmamgJee027y1x0TbmmtFSzQv7hD1l/0CzgWEC0mbA6ZsePA6dOaUY8olxUVaUNgq4ubSBkMacCYtqdPq3H/tfW5mReCqJpK1cCH32kidHSrKWlxVpZTgVEmwmrl3TxIvCnPzFXMhGgDYLaWuDYMT0+LI06OzutleVUQEyb4WFdclBRoYtViUgbBjU1wH//d1pPyGltbbVWFgNivMbHdVbN59P1hkQ0o6BA9+2/+GLaDpg9ePCgtbKcCog2E1ZHFQppdrzhYZ5cQ7SQZcuAvDw9GX5iwuvaxMWpgJjynSrvvaepQjmjTLS45cv1UOTu7qxKWuVUQLSZsHqegQEdMK6tTd01iFxSU6Mzzyk+Hcfm7jSnAqLNhNVzBIM6o1xezhlloliJ6Ok4R46kdCeLzZPynQqIKfP229r857Y8ovhEskqmsOu8Z88ea2U5FRBtJqyeduGCBkSOGxIlprpa1yaePu11TZbkVEC0eZQ4AJ0h+9OfdIA4z6n/KqL0WrkSePXVjM/17NRfuc2E1QB0EmVsjOsNiZJVUKBH4736qvX9zh0dHdbKciog2kxYjUAA+OAD5lAmsqWqSg+U/fBDq8XazLZpJSCKyAMi8r6InBSRx6M8vl1EBkXkjfDtCRvXTZmxMe0qV1Wxq0xk06pVeirO4KC1IhsaGqyVlfRfu4j4APwIwIMAtgD4kohsifLU3xtjbg3f/jnZ66aMMcDrr+vX4mKva0PkFr9fh6COHNGdXxnGRvNnK4CTxpiPjDETAH4J4GEL5cbNSsLqQEAXk65YkXxZRDRfRYUe/nDypNc1mcdGQGwAcGbWz2fD913rThE5JiLPiMhnLFx3nqRnmUMh4OhRnVXm+YZEqVNdrWeJjo0lXdTu3bstVEjZCIjRIse1KzCPAlhnjLkFwL8C+PWChYm0iUiXiHT19fXFVZGkE1Z//LEuC+CsMlFq5efrbPP77yddVKbtVDkLYM2sn1cDmJMg2RhzxRgzHP7+EIB8EamOVpgxZr8xptkY01yTzhneiQkdO2RXmSg9qqv1hO2hoaSKybRZ5tcAbBaRDSJSAOBRAHMWBolIrYRTY4nI1vB103eCZCxOnAAmJ3ngK1G6+Hz693b8eFLFHD161FKFAH+yBRhjJkXk6wCeA+AD8DNjzNsi8rXw4z8B8EUA/yAikwBGATxqUpBAOeGE1Vev6pvC1iFRelVV6STmpz6VEVkrnUpUn7CuLs2cx0XYROk3OKgHp9xzT0KTmfX19QgEAks/cZacSFSfUMLqwUHdkcLWIZE3KiqAc+f0loB4g+FinAqICTl2DCgq4o4UIi9VVuqStwQWa9s8KT+3o0Bfn6YEWL7c65oQ5bZly7S39skncf/TJ5980lo1nAqIcSWsNkY/kcrLuQibKBOsWKFL34JBz6rgVECMK2H1hQu6faisLHUVIqLYFRZqmt8zZ5Z+boo4FRDjSlj93nvaTCeizFFZqSfUx3Fmos2VKE4FxJgTVg8N6SEO5eWprRARxae4WP8+L1705PJOBcSYnTqlq+Q5dkiUeUpK4trj3Nw8bzlhwnIvIAaD+p+dAaviiSiKigrN5ZzkHudEOBUQY9p1c/as7ln2J71rkYhSQUR7cKdOpf3STgXEJY8BMkZP16isTE+FiCgxVVXak4thCU57e7u1yzoVEJdMWH3xoi7+ZGoAoszm92tPLoZtedypkqgPPmAwJMoWFRV6CtUSQ2H19fXWLpk7AXF4WLcFsbtMlB1KSrRHt8QSnN7eXmuXdCogLpqw+vRpLrUhyjbFxdqzSxOnAuKCR4lHcjfwEAei7FJZqT27RZJRNTY2WrucUwFxwYTVly9rzpT8/PRWiIiSE+nRLdJtTjrb5ixOBcQF9fZy3SFRtiop0SGvBbS1tVm7lPsB0Rjgww+5b5koW5WVAT09C65JPHDggLVLORUQoyasHhoCRkb0aCEiyj55eToP0J/6RJ1OBcSoO1XOn2d6AKJsV1iYlnMSnYoUUWeZP/yQh8ASZbvycp1tjpJzpaenx9plnAqI8xJWX70KXLqkg7JElL38fl0pcunSvIc4yxyrvj6va0BEtvj9Ufc279q1y9olnAqIdXV1c+84dUoTYBNR9quo0L/pWI75S5BTAXFOwurxcV1/yIBI5IaCAh0GGxxM2SWcCohzjgHq79dV7ty7TOSOvDzNmDnLvn37rBUvMZ0y7ZHm5mYTT0YtEZk5NfvYMeDkSc31SkRuGB7WXt/27UkVIyLdxph5yVicaiHOce4cZ5eJXFNcrPuaZzXkxGIv0M2AGArpgQ5FRV7XhIhs8vn0JO2RkZQU71RAnO5eDw3pJwjHD4nclKKMfE4FxGmRgEhE7vH55uxrbmlpsVa0UwFxOmF1Xx8PcyByVUmJnlEQ1tnZaa1opwLitPPnmUyKyFWRiZWpKQBAa2urtaLdC4jBICdUiFyWl6dDYsPDAICDBw/aK9paSRmgvb1dxw+5IJvIbcakZGLFqYC4d+9e4MoVr6tBRKmWn5+SA2OdCoj19fW6rYcTKkRumzWxYnO3nVMBsbe3V2eYuUOFyG3FxcDAABAKRT8pP0FOBUQAOq5QUOB1LYgolUR0HHFiAnv27LFWrFMBceutt+o3nFAhyg3j41aLcyogHnn5Za+rQETpNDFhtTinAuL/euwxtg6JcoUxwPg4Ojo6rBXpVED8r//8z+nV60TkOJ8PGB2Nnm0zQU4FxEKALUSiXOH3A0NDaGhosFakUwGxFND/JCJyXyTHikVORY/f/eY3PPaLKFf4/dP7mW1xqoV48tgx3dJDRO7LzwdGRrB7925rRToVEL/zxBMMiES5wucDgkHs//GPrRXpTkCcmkIxwDFEolwigrtvv91ace4ExIkJGICzzEQ55p033rBWljsBcWoKu7/6Va9rQUTpJGI1iLkTEI3Bjnvv9boWRJROxqBu1SprxVkJiCLygIi8LyInReTxKI+LiPww/PibItJo47pzGINHvvQl68USUWY7/uab1spKOiCKiA/AjwA8CGALgC+JyJZrnvYggM3hWxuAf0/2uvNw/SFRTvr+975nrSwbLcStAE4aYz4yxkwA+CWAh695zsMAfm7UKwAqRaTOwrVnMCAS5aQf/OAH1sqyERAbAJyZ9fPZ8H3xPic5xqCx0X5PnIgym811JTYCYrT6XNtci+U5+kSRNhHpEpGuvr6+uCryrW9+M67nE1GWs9wztBEQzwJYM+vn1QACCTwHAGCM2W+MaTbGNNfU1MReCxF897vfjf35RJT9RPDss89aK85GQHwNwGYR2SAiBQAeBXDtiY0dAP4mPNu8DcCgMabXwrVniODo669bLZKIskCevdWDSe9zM8ZMisjXATwHwAfgZ8aYt0Xka+HHfwLgEIDPAzgJ4CqAv0v2uvNwhwpRTtpx333os9R1trLx1xhzCBr0Zt/3k1nfGwCP2bjWghgQiXKSzVFEd3aqiOC/fvELr2tBROkkwoAYld+PFw4f9roWRJRm//tb37JWljsBMT8fB376Uy7QJsolxuDx9nZrxbkTEEUwAgDBoNc1IaI0Wr1xo7Wy3AmI0OlrBkSiHDE5CRQUoOfcOWtFOhUQ937/+/qfRETuCwaB0lKrRToVED/V1ARMTHhdDSJKh2AQKCmxeoaBUwFx6z33AFNTXleDiNJhchIoK0N3d7e1Ip0KiBMAZ5mJcsXEBLBsGdra2qwV6VRAHAe4Y4UolxQX48CBA9aKcyog/vXf/i1biES5pKDAanFOBcR/3b/f6yoQUbqIAIWFVot0KiA2bdsG+HycWCHKBVNTug6xp8dakU4FxKNHjwLl5cD4uNdVIaJUCoUAvx8oLOQs86JWrQKuXvW6FkSUSqOjQHU1IIJdu3ZZK9apgFhXV6f/SVycTeS20VFt/FjmVEAMBAJAWRmX3hC5LhQCqqqsF+tUQNy7d6/ubRThxAqR68rKAAD79u2zVqRTAfHJJ5/UWeaqKmBszOvqEFEqRCZUSkoAgDtVlsSJFSJ3Xb0KrFw5PTQmFofI3AyI1dU8F5HIVSmaUAEcC4hdXV36DSdWiNw1NQUsX56Sop0KiNOWLdPk1aGQ1zUholQIT6gAQEtLi7VinQqIzc3N+k1ennabR0e9rRAR2RUM6v7l4uLpuzo7O60V71RAnGPVKgZEItdEJlRmaW1ttVa8uwGxpob5VYhcMzoK1NXNuevgwYPWincqILbPzs9aVcVxRCLXGDOvhWiTUwFx7969Mz/4/cDatcDQkGf1ISKLRkeBigrrmfZmcyog1tfXz71j3TqOIxK54soVIEpSemPxlHynAmJvb+/cO1as0K9MK0CU/aamgNraeXfvt3hSvlMBcZ7CQh2AHR72uiZElIyJCd27XFEx76E9e/ZYu4xTATFqwuoNGxgQibLd4KD+Lad4B5pTATHqUeI1NemvCBHZNTkJXDtHkAJOBcSoxwCVlOgSHJ5+Q5SdJieB/PwF9y93dHRYu5RTAXHBhNUbN3L5DVG2unJFV4z4fFEfbmpqsnYppwLiglat4gnaRNlqfBxYs2bBhxsaGqxdKjcCYlmZLubkKdpE2SUU0pZhZAldijkVEBdMWC0C3HADcOlSeitERMkZGAA2b9YxxDRwKiAumrB6zRrubSbKJsbocV9RdqfMtnv3bmuXdCogLpqwurBQP2nYSiTKDkNDujMlymLs2bhTJVEbN+qKd27lI8p8w8M61LUEzjInqrJSP3G4BIcos42P6xriGJJJHT161NplnQqIMSWsvuEGbuUjynSXLgGf+YyO+6eRUwExpoTVK1fqJ8/4eOorRETxC4V0Zcgiaw9nq7vmBO1kOBUQY0pY7fMBW7YAly+nvkJEFL9Ll4DrrgOKimJ6eiAQsHZppwJizNau1YkVLsEhyjwTExoQYzTnpPwk5WZALCrS/3AuwSHKLMPDekJVHInon3zySWuXdyogxpWwetMmLsEhyjSDgzqZ4hGnAmJcCauXL9euM1uJRJlheHhmaZxHnAqIcSesvvlmPfCBJ+EQeW9wEGhqinupTVdXl7UqOBUQ405YXV4OfPrTQH9/aipERLG5fFlPxI5hIXYqORUQE3LDDdpCnJz0uiZEuckYYGQEuOWWhHKmNDc3W6sKA2JxMXDTTcDFi17XhCg39ffrJGdVldc1SS4gikiViDwvIifCX6POlYvIaRF5S0TeEBF7Hf5rJJywetMmoKCAu1eI0i0U0iO+PJxZni3ZFuLjAA4bYzYDOBz+eSGfM8bcaoyx1769RsLHABUUALfdpodRElH6XLyowbC0NOEi2tvbrVVHEm5VARCR9wFsN8b0ikgdgBeNMZ+K8rzTAJqNMXH1S5ubm008M0gikngrMRQCnn1Wv1+2LLEyiCh2ExM6s9zaGvM2PVtEpDta4yzZFuIqY0wvAIS/rlzgeQbAb0WkW0RiOIHBAz4f0NjIPc5E6dLfD9x6a9LBsN5ivmb/Uk8QkRcARFsp+e04rnO3MSYgIisBPC8i7xljXl7gem0A2gBg7dq1cVzCgtpavQ0OLnlKLxElYXRUT53asCHponp7ey1USC3ZQjTG7DDG3Bjl9hsA58NdZYS/XligjED46wUATwHYusj19htjmo0xzTU1NXG9mKQTVovoWOLICA9+IEoVY3S8vrER8C/ZJkurZLvMHQC+Ev7+KwB+c+0TRGSZiJRFvgdwH4DjSV43KitHiVdVATfeCFyIGtuJKFn9/cD69cDq1VaKa2xstFIOkHxA/A6AnSJyAsDO8M8QkXoRORR+zioAfxCRYwBeBfC0MebZJK8blbWE1Vu26J7KK1fslEdEanRUe2KNjQktwo5m0WybcUoqIBpj+o0x9xpjNoe/DoTvDxhjPh/+/iNjzC3h22eMMf/HRsVTyu8H7ryTXWcimyJd5W3bdEOEJTGdlB8j7lRZSGWljiey60xkR1+fpgK21ZMLO3DggLWynAqINhNWA9A3r7qaS3GIknX1quZGv+UWr2uyKKcCos2E1QB0beIdd+gRYcGg3bKJcsXUlHaV77xTg2IGcyog2kxYPa28HGhuZteZKFEXLujKjZUL7dtITk9Pj7WynAqINhNWz7Fxo457cK8zUXyGh3WfcgoPb8iYWeackZcH3H67npk4MeF1bYiyQyiku77uvBPIz0/ZZXbt2mWtLKcCos2E1fMsW6bjiX19TDlAtBRjgHPnNE3HihVe1yZmTgVEmwmro1q3Tpv+586l9jpE2e7ChZm/lyziVEC0mbB6QTffDKxZw0kWooUMDOjhKFu3xp0wKhH79u2zVlZS5yGmWlrPQ4zHxARw+LAux4kjoTaR80ZG9O/igQcy+lzRVJ2HmJsKCoA//3MdSxwZ8bo2RJlhYkL3/2/fntZgKJb2RAMMiIkrLdU3/soVzjwThUI64Xj33Vk1iXItpwKizYTVMVmxArjrLv1F4CEQlKsiM8q33gqk+1Bny5wKiJ5Yt073Z547p78YRLnm/Hnguuv02DwPtLS0WCvLqYBoM2F1XLZs0d0snHmmXDMwoD2l5mZr5xvGq7Oz01pZTgVEz+Tl6S/E8uV6GjBRLrhyRQ9A+bM/8zQVQGtrq7WyGBBtyc/XmeeiIu55JvcNDekJUNu3Wz3sNREHDx60VpZTAdFmwuqEFBcDn/ucHnHEoEiuGhrSlRU7djiXndKpgJiWnSpLKSlhUCR3DQ8D4+PAvfc6FwwBxwKizYTVSZkdFC9d8ro2RHYMD+sulB07NMVGhrC5O82pgGgzYXXSIkGxoIATLZT9rlzRYHjvvRkVDAG7J+U7FRAzTkkJcM89uo3p4kWva0OUmMuXdePBzp0ZuXd/z5491spyKiDaTFhtTWSipbJSF7ASZZOBAV1as3Onk2OG13IqINo8StyqwkLgs58FVq3ijhbKHhcv6jKye+/Vvfs5wKmAaDNhtXUFBbqAdc0aoLeXe58pcxmjvZnych3yKSnxukaL6ujosFaWUwHRZsLqlPD7Nb/EjTdqS3F83OsaEc01OQn09Oge/c99TluIGc5mtk3v9tvkqrw8PXW7qgr44x91jLG83OtaEQGjozpmuHUrsHmzZ3uT49XQ0GBt6Y1TLcSssno1cP/9+kvX1+d1bSjXDQzoYcc7dwLXX581wdA2pwKizYTVaVFZCdx3n0629PRwXJHSLzJeWFqqx/6nKJl8tnAqIGbsLPNiCgt1suWmmziuSOkVGS9cv17HC7N0Jnn37t3WymKSqUxy9izHFSk9rl7VbaVbtwKbNuVcF5lJprLB6tXabRHRw2azObhT5rp0SSdQ7rsvqyZPFmJzlpkBMdNUVOgvakMDEAjo/lEiG4JB/Z0qK9MP3poar2tkxdGjR62V5dSyG5sJqz1VWKjJq9atA44c0fPnqquz/pOcPDQwoGcY3n675j9JQwL5bORUQMzonSrxEtEudHU18MYbwIcf6tpFj08npiwzMaHLuhoaNM1Flk6cLKaurs5aWU4FxKyfVImmqAjYtk1bi6+8oscw1dTwE54WZ4weOzc1pbmS161ztocRCASslcW/qmxRVwc89JDOCAYCelgnUTRjY/o7UlcHtAHI42AAAAjySURBVLToshpHgyFg96R8LrvJRn192locHtaFtD6f1zWiTGCM/m74fMAdd2g32eFAGJHI3/1Cy26c6jLbTFid0WpqdJbw3XeB48f1ANocOKuOFjE8DAwOag/illuy4lCGTORUQLSZsDrj5efrIRGrVwPHjumOg7IyLujONZEF1lVVmutk1Sqva5TVnBpDtJmwOmtUVWlu3J07dQa6p4fji7lgdFTHCaemgL/4Cz0oJEeDYTzDaktxqoVoM2F1VhHRscQdO3Sj/uuva2BcvjzjD/ekOI2P65rCkhLdA9/QwDFki5wKiDlPBKit1dZCIKCBMRDQwMj1i9ltYkKX0RQW6oTJ2rV64DChubnZ2mQq/0ddlJenY4t1ddpSfP11zZxWVaV/UJQ9gkENhH6/Lqxev17HjyklnAqIObHkJh4+n7Yk6uuBTz7RHS8DAzr54uCOBaeMjuqHmM+ns8bXXad5eSilnJpUsZmw2il+P7BxI9DaquNOhYXale7r46G0mSQU0g+syGTJHXcADz8M3HADg+Ei2tvbrZXFhdm56vJl4NQp4MQJ/UMsL9f1jJR+o6O6hhAANmzQ1uCKFTmxqNorObEwm+JQWQncdptmAOztBd55R1sm+fk61siZy9SamtIPpbExHcJobtZxXy6ojlt9fb21/cwMiLkuP1/HGdes0T/Qjz7Sk3UircaSErZUbIq0Bo3RCZJNm3i0W5J6e3utleVUQLSZsDrniOjynKYmze8SCAAffKB5XgANjKWlbDnGa2pKF8qPjOjPFRX6f7x6NZdCZSCnAqLNo8RzWkGBtl7Wr9cWTX8/cOaM3kIhnaQpL+cSnoUEg3pM28SELoGqq9OZ4upqjtOmQGNjo7WyOKlCsQuFdN9sb69OyIyMaMuytFT/0HO122eMfnAMDWmLsLBQP0waGnQ8lusGMw4nVSh5Pp+2cqqrdTJmeFiX7pw+rV1rEb0VF+vN1UAQCumhCqOjGgCN0cB32226hbKiInc/HDzQ1tZmbcldUi1EEfkrAHsB3ABgqzEmanNORB4A8C8AfAB+aoz5Tizls4WYRcbHdbLg8mUNkn19GjAigSESJLNtPV0wqK9jdHQmC6Lfr4Fv5Uqdra+o4HighzLpPMTjAP4SwILZnUTEB+BHAHYCOAvgNRHpMMa8k+S157GZsJriVFg4EySuv17vGxubOaevr09Tq/b3zwTJggK9+f3amvRqwiYU0qTtwaDexsdnHou8rpoaDX6lpRr82AJ0UlIB0RjzLqARehFbAZw0xnwUfu4vATwMwHpA5E6VDFNUpLfqal1sDOhEw9CQ3vr7tes5PKwty4mJmW737E/8/PyZW+Txa2+A/pvILfLz1JQGuUjAizwWuYYxWu6yZboesKREu7/l5foz1wXmFmNM0jcALwJoXuCxL0K7yZGfvwzg3xYpqw1AF4CutWvXmng0NjYaANO3rq4u09XVNee+9vZ2Y4wxdXV10/c1NjYaY4zZvXv3nOf29PSYjo6OOfft27fPGK3o9K2lpcUYY0xLS8uc+40xZt++fXPu6+joMD09PXPu271797z619XVGWOMaW9vz5nXlAeYIsCUA+a73/iGMWfOmG3V1eYmwGwDzNc2bTLm0CHzL/ffbx4CTCtgdgHm4oED5tXHHze7wvc9BJhDjz1mzKFDZgdg7gLMLYBp277dmE8+MX+9c6epAEwxYHx8n5x4TfEC0GWixJ8lxxBF5AUAtVEe+rYx5jfh57wI4BsmyhhieJzxfmPMV8M/fxk63vg/F70w4h9DJCKKRcJjiMaYHUle+yyANbN+Xg3AXt5AIiJL0nHazWsANovIBhEpAPAoAG4pIaKMk1RAFJEviMhZAHcCeFpEngvfXy8ihwDAGDMJ4OsAngPwLoBfGWPeTq7aRET2JTvL/BSAp6LcHwDw+Vk/HwJwKJlrERGlmlMHxBIRJYMBkYgojAGRiCiMAZGIKIwBkYgojAGRiCiMAZGIKCyjT8wWkT4AH8fxT6oBXExRdbJBLr9+vvbclcjrX2eMqbn2zowOiPESka5oG7ZzRS6/fr723HztgN3Xzy4zEVEYAyIRUZhrATHXj8zO5dfP1567rL1+p8YQiYiS4VoLkYgoYVkdEEXkr0TkbRGZEpEFZ5lE5AEReV9ETorI4+msYyqJSJWIPC8iJ8Jfly/wvNMi8paIvCEiWZ2TYan3UtQPw4+/KSKNXtQzFWJ47dtFZDD8Pr8hIk94Uc9UEJGficgFETm+wON23vdoiVay5QbNB/0pLJ7kygfgQwAbARQAOAZgi9d1t/T6vwfg8fD3jwP47gLPOw2g2uv6Wni9S76X0HM4nwEgALYBOOJ1vdP42rcDOOh1XVP0+j8LoBHA8QUet/K+Z3UL0RjzrjHm/SWeNp0G1RgzASCSBtUFDwP4j/D3/wHgf3hYl3SI5b18GMDPjXoFQKWI1KW7oing8u/xkowxLwMYWOQpVt73rA6IMWoAcGbWz2fD97lglTGmFwDCX1cu8DwD4Lci0i0ibWmrnX2xvJeuvt+xvq47ReSYiDwjIp9JT9UygpX3PakUAukQSxrUpYqIcl/WTK0v9vrjKOZuY0xARFYCeF5E3gt/4mabWN7LrH6/FxHL6zoK3ZI2LCKfB/BrAJtTXrPMYOV9z/iAaHI8Depir19EzotInTGmN9w9uLBAGYHw1wsi8hS0+5WNATGW9zKr3+9FLPm6jDFXZn1/SER+LCLVxphc2Ods5X3PhS6zy2lQOwB8Jfz9VwDMazGLyDIRKYt8D+A+AFFn6rJALO9lB4C/Cc86bgMwGBlWyHJLvnYRqRURCX+/Ffr33Z/2mnrDzvvu9exRkjNPX4B+MowDOA/gufD99QAOXTMD9QF0lu7bXtfb4utfAeAwgBPhr1XXvn7orOSx8O3tbH/90d5LAF8D8LXw9wLgR+HH38ICqw+y8RbDa/96+D0+BuAVAHd5XWeLr/0XAHoBBMN/83+fivedO1WIiMJyoctMRBQTBkQiojAGRCKiMAZEIqIwBkQiojAGRCKiMAZEIqIwBkQiorD/D5H0nrUA/+SxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.patches as patches # necessary for geometric drawing\n",
    "\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 5))\n",
    "ax.set_xlim(-1.1, 1.1)\n",
    "ax.set_ylim(-1.1, 1.1)\n",
    "\n",
    "rect = patches.Rectangle((-1, -1), 2, 2, \n",
    "                         linewidth=1, \n",
    "                         edgecolor='k', \n",
    "                         facecolor='none', \n",
    "                         linestyle='dashed')\n",
    "circle = plt.Circle((0, 0), 1, color='r', alpha=.25)\n",
    "\n",
    "ax.add_patch(rect)\n",
    "ax.add_patch(circle)\n",
    "_ = ax.set_xticks([-1, -.5, 0, .5, 1.])\n",
    "_ = ax.set_yticks([-1, -.5, 0, .5, 1.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can throw some points into our domain $[-1, 1]$, marking with a different color the points **inside** and **outside** the circle sub-domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUQAAAEvCAYAAAA92bhfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de2zc1Zk38O/j8T2+ZWInsZ2EJBBaUq6OGwLslpQQoCUJL6gLrVbd7pbGdNW+i5D6lmgr1WalV2/aqkLpbrtNUtBu9492+aOhTkihEAS03QK1A7lwTSCBxDaJYyeO7fg6c94/nhlfx/Zczsxv5sz3I43sGU/O70xm/PjcHzHGgIiIgByvK0BElC4YEImIQhgQiYhCGBCJiEIYEImIQhgQiYhCcr2uwGwqKyvN8uXLva4GETmmtbX1nDGmaurjaR0Qly9fjpaWFq+rQUSOEZGPIj3OLjMRUQgDIhFRCAMiEVEIAyIRUQgDIhFRCAMiEVEIAyIRUYiVgCgiT4rIWRE5OsPPRUR+IiLHReSwiNTZuC4RkU22Woj/AeCuWX7+BQCrQrcGAP9u6bpERNZY2alijHlFRJbP8pR7APzS6PHcr4pIhYhUG2M6bFyfskQwqDdjpt8AQGT6LSdHb0RRSNUnpRbAqQn3T4cem0ZEGkSkRURaOjs7Y7rI8uXLISKoqakBADQ1NUFExm6tra1obW2d9FhTUxMAoKamZuyxNWvWAAAaGhomPbe9vR179+6d9NiuXbvC9R67bd68GQCwefPmSY8DwK5duyY9tnfvXrS3t096rKGhAQCwZs2ascdcfk3/8v3vo1gEFSJYKILD+/fjyJ49uE4Et4hgowh+9cADwFNP4R8rK/Hl/Hw8UFCAbZ/+NPD009h19914oLgYDxQV4YHCQnT/4hdo/ed/xgN5efiyz4d7RfDMP/0T8Ic/4AYRXCmCpSL46h13AOfP40tf+ALy+D5l7GtabvG8A7GVUyXUQtxnjLk6ws+eAfD/jDF/DN0/AOC7xpjW2cqsr683sexlXrNmDVpbZy2SvDQ0BPT1ARcvAp2dQFcX0N8PjI5Of64IkJcH5Obq17w8wOfTx2M1OgqMjEz+GgxOf57PBxQXA34/UFUFlJcDJSVAUVHs16SUERHEGsdEpNUYUz/18VQd7nAawNIJ95cAaLd9EQbDNDI4CPT2avA7e1Zvly6NB7TCQr3Nn6+BKJlyc/U2l0BAg+XZs8BHob3/xgAFBRogFy3SIFlaqkEynuBMaS1VAbEZwLdF5NcAbgTQk4zxw4aGhrEmN6VQIABcuKAtvjNngHPnNCCGA0ZRkd7Ky72t51x8Pr0VFEyu68gIcP480NEx3rLMzdUguXAhsGCBtiqjCbpkXXNzs7WyrHSZReRXANYDqARwBkAjgDwAMMb8XLTD/2/QmehLAP7BGDNnXzjWLnM8TWeK0+Ag0N0NnDoFfPyxtqxycsaDX16e1zVMrkAAGBjQVm8wqMG/pga47DKgslK73pQS7e3tY+Oc0Upql9kY85U5fm4AfMvGtcgjxoyP/Z04oa1AQLu9fn/yu73pxufT8cWSEr0fCOgfiNOn9X5FBbBixXg3mzPdSVNbW2utIcQ2Ps1sdFR/yTs6gJMntUUEaBBYvJhjaBP5fBr4wl3tgQHg0CFtPebna8txyRL945Gf721daUZOBcS2tjavq5D5jNHxsg8/BD74QO/n5gJlZek/BphOwkMHgP5h+fhj4Ngx/SOybBlwxRXatWbLMa04FRBbW1tjHkugkOFhbQm+8462CvPz9Rc227rCyZCbqy1DQFuMZ85oi7ukBLjqKm05cmlP3LZu3WqtLGvrEJMh0ydVjJncq5x6Py2cP6+/nMeO6ThYWRkwb57XtcoOg4M6Ow/oeOPll+uMddp9SNzj9TrErNPUpJ/1xx/Xz7cxwCOP6Fh7aDG/d0ZGxluDXV3aGlywgK3BVCss1LHYYBBoa9MhirIybTXW1urPaU42N2QwICaBMRoMd+zQ+48/rsFwxw7g4Yc9bClevKgzxO+/r+NapaW6VIS8lZMz3qUeGABef10fW7FCxxrDP6OIDh48aK0spwLizp07va4CAA12jz+u3+/YMR4YH354vMWYUn19wNtvA8eP6/rA+fO5iDhdhSdjAgFdwnP8uI4xXnONvm+UVE6NIaYbYyZPIobX76bMpUvAe+8B776rAXDBAs5qZprwrP/gILByJbB6tXaraUxNTQ3a22PbCZwVY4jpNKkSHjOc6JFHUtRCHBzUbvHbb+u44MKFHB/MVCLaZTYGaG/X5VCrVuk4Y3hReJaLNRjOhs2FJAgHw/CYYTCoX3fs0MeTFrOHhjQINjfrhElVld4YDDNfODAuXqyrAvbuBd54Q3sBWa7J4iylU13mdGohpnSWeWREJ0vCOyMWLOAYoevCWwUDAe1GX3ll1s5KZ+LxXymxadMmr6swpqlp8mxyeKLFanfZGD2m6uBBbR0uWMBtYdnC59PWfyCgvYF33wWuvVa70+wRxM2pgLh3716vqzDJ1OBnNRj29gKtrbp+rbKSSzOylc+nB0iMjGgX+sQJ4MYb+XmIk1NjiOGjxp0WDOqukv37tctUW6vn91F2y8sDqqs1MD77LHD4sH6fBWyuRHGqhbhv3z6vq5BcPT26aLezU7tLrp85SLELb7186y0dTlm3Tj8rFBWnWojOCgR09nj/fp1VrKlhMKSZ+XzaWhQBfv97HWMeHva6VklTXz9tbiRuTrUQndTdDbz2mk5Zcz0hxaKkRE/uPnZMjx+78UYNlDQjp1qI6bLkxoqRER0H+t3vdN9xdTWDIcUuJ0cnXQoKgBdfBF59VRfuU0ROBURnEkz19GhX5+23dSFuaanXNaJMV1SkQy2nTgHPPKOZBR3R2NhorSwuzE43p08Df/qTfoC5Z5WS4dIl3R+9dq2eppOF5y/OtDDbqRZiRgsGtYv88su6nYXBkJKluFi70a+9pqsWMnx5js1T8jmpkg6GhvTDeeqUdpE5VkjJlpura1hPntQhmptvztjDIjo67KV4d6qFaDNhdcpcuKDjhWfO6AeUwZBSRURbin19upjboXHFeDkVENesWeN1FWJz+jTw3HO6J5mLZ8krfr8u5n7+eT02LsPG4evq6qyV5VSX2WbC6qQKBoGjR4EjR3QfMrfekdeKirS1+PrruvZ1zZqMWfxvK58K4FgLMSMMDgJ/+IMGxMWLGQwpfYTHFT/6CDhwQLvSGaChocFaWQyIqXTpki6OPXtW14RxvJDSTXhccWBAx7bDaVLT2O7du62V5VRAtJmw2rreXuCFF7SFyPFCSnfz52uX+fnnNVVtlnAqIKbtTpWeHg2GwSDPqaPMUVqqaxZfeCFrZqCdCohpOcvc3a1/ZX0+XXBNlEnmzdNNAgcOaJKrNNTW1matLKcCos2E1VZ0delf18JC7jyhzFVUpD2bl17SpWJphrPMmSAcDOfNy9gdAERjCgp0idjLL6ddUNyyZYu1spwKiNXpctZbd7d2MUpKNCASuSA/X4PiK69oLh8HORUQbSasjlt393jLkMGQXFNQoNkdX37ZyaDoVEC0mbA6LufPazAsLmYwJHdNDIoWD1aI186dO62VxfMQbenv14WseXkcM6TsMDioC7fvvFPXLWYQnoeYTCMjeqirMQyGlD0KC/Xz/sorugvLI2LxgFsGxEQFg8Bf/qLdZS66pmxTUqJZIf/4x4w/aBZwLCDaTFgdtaNHgRMnNCMeUTby+7VB0NKiDYQM5lRATLmTJ/XY/8WLszIvBdGYhQuBDz/UxGgptmnTJmtlORUQbSasntO5c8Cf/8xcyUSANggWLwYOHdLjw1Jo79691spyKiCmTF+fLjkoL9fFqkSkDYOqKuB//ielJ+Rs3rzZWlkMiLEaGtJZNZ9P1xsS0bj8fN23/9JLKTtgdt++fdbKciog2kxYHVEgoNnx+vp4cg3RTObNA3Jy9GT44WGvaxMTpwJi0neqvPuupgrljDLR7ObP10ORW1szKmmVUwHRZsLqabq7dcB48eLkXYPIJVVVOvOc5NNxbO5Ocyog2kxYPcnIiM4ol5VxRpkoWiJ6Os5rryV1J4vNk/KdCohJ89Zb2vzntjyi2ISzSiax6/zQQw9ZK8upgGgzYfWYs2c1IHLckCg+lZW6NvHkSa9rMienAqLNo8QB6AzZn/+sA8Q5Tv1XEaXWwoXA66+nfa5np37LbSasBqCTKIODXG9IlKj8fD0a7/XXre93bm5utlaWUwHRZsJqtLcD77/PHMpEtvj9eqDsBx9YLdZmtk0rAVFE7hKR90TkuIhsi/Dz9SLSIyJvhm7ft3HdpBkc1K6y38+uMpFNixbpqTg9PdaKrK2ttVZWwr/tIuID8FMAXwCwGsBXRGR1hKf+wRhzfej2L4leN2mMAd54Q78WFXldGyK35ObqENRrr+nOrzRjo/mzFsBxY8yHxphhAL8GcI+FcmNmJWF1e7suJl2wIPGyiGi68nI9/OH4ca9rMo2NgFgL4NSE+6dDj011k4gcEpHfichnLFx3moRnmQMB4OBBnVXm+YZEyVNZqWeJDg4mXNTWrVstVEjZCIiRIsfUFZgHAVxmjLkOwL8CeHrGwkQaRKRFRFo6OztjqkjCCas/+kiXBXBWmSi58vJ0tvm99xIuKt12qpwGsHTC/SUAJiVINsZcNMb0hb7fDyBPRCojFWaM2WWMqTfG1FelcoZ3eFjHDtlVJkqNyko9Ybu3N6Fi0m2W+S8AVonIChHJB/BlAJMWBonIYgmlxhKRtaHrpu4EyWgcOwaMjvLAV6JU8fn09+3o0YSKOXjwoKUKAbmJFmCMGRWRbwN4DoAPwJPGmLdE5Juhn/8cwJcA/KOIjAIYAPBlk4QEynEnrL50Sd8Utg6JUsvv10nMT30qLbJWOpWoPm4tLZo5j4uwiVKvp0cPTrnttrgmM2tqatDe3j73EyfIikT1cSWs7unRHSlsHRJ5o7wc+OQTvcUh1mA4G6cCYlwOHQIKC7kjhchLFRW65C2Oxdo2T8rP7ijQ2akpAebP97ommSPGpVBEUZk3T3trH38c8z997LHHrFXDqYAYU8JqY/QvUlkZF2FH6/Bh4P77gSNHvK4JuWjBAl36NjLiWRWcCogxJaw+e1a3D5WWJq9CLgkEgO3b9fvt29NyHypluIICTfN76tTcz00SpwJiTAmr331Xm+kUnd/8Bjh/Xr/v7gb27PG2PuSmigo9oT6GMxNtrkRxKiBGnbC6t1cPcSgrS26FXNHdDTz55Pi+08FB4Ikn9HEim4qK9Pfz3DlPLu9UQIzaiRO6Sp5jh9E5cGB6FzkQAF580Zv6kNuKi2Pa41xfP205YdyyLyCOjOh/dhqsis8YGzZMT7/q8+lCWiLbyss1l3OCe5zj4VRAjGrXzenTumc5N+Fdi9nD7we+/nVdrwno1wcf5B8VSg4R/YN74kTKL+1UQJzzGCBj9HSNiorUVMgl9903HgD9fuDee72tD7nN79eeXBRLcBobG61d1qmAOGfC6nPndPEnUwPEzucDHn1Uv9+2bXoXmsim3FztyUWxLY87VeL1/vsMhom49lrgqaeAa67xuiaUDcrL9RSqOYbCampqrF0yewJiX59uC2J3OTE8EYhSpbhYe3RzLMHp6OiwdkmnAuKsCatPnuRSG6JMU1SkPbsUcSogzniUeDh3Aw9xIMosFRXas5slGVVdXZ21yzkVEGdMWH3hguZMyctLbYWIKDHhHt0s3eaEs21O4FRAnFFHB9cdEmWq4mId8ppBQ0ODtUu5HxCNAT74gPuWiTJVaSnQ1jbjmsTdu3dbu5RTATFiwureXqC/X48WIqLMk5Oj8wBdyU/U6VRAjLhT5cwZpgcgynQFBSk5J9GpSBFxlvmDD3gILFGmKyvT2eYIBxO3tbVZu4xTAXFawupLl/RQ0+JibypERHbk5upKkfAhxRNwljlaTIhE5I7c3Ih7m7ds2WLtEk4FxOrq6skPnDihCbCJKPOVl+vvdDTH/MXJqYA4KWH10JCuP2RAJHJDfr4Og/X0JO0STgXESccAdXXpKnfuXSZyR06OZsycYOfOndaKl6hOmfZIfX29iSWjloiMn5p96BBw/LjmeiUiN/T1aa9v/fqEihGRVmPMtGQsTrUQJ/nkE84uE7mmqEj3NU9oyInFXqCbATEQ0AMdwjlAiMgNPp+epN3fn5TinQqIY93r3l79C8LxQyI3JSkjn1MBcUw4IBKRe3y+SfuaN23aZK1opwLiWMLqzk4e5kDkquJiPaMgZO/evdaKdiogjjlzhsmkiFwVnlgJBgEAmzdvtla0ewFxZIQTKkQuy8nRIbG+PgDAvn377BVtraQ00NjYqOOHXJBN5DZjkjKx4lRAbGpqAi5e9LoaRJRseXlJOTDWqYBYU1Oj23o4oULktgkTKzZ32zkVEDs6OnSGmTtUiNxWVAR0dwOBQOST8uPkVEAEoOMK+fle14KIkklExxGHh/HQQw9ZK9apgLj2+uv1G06oEGWHoSGrxTkVEF975RWvq0BEqTQ8bLU4pwLiI9/6FluHRNnCGGBoCM3NzdaKdCog/vd//dfY6nUicpzPBwwMRM62GSenAmIBwBYiUbbIzQV6e1FbW2utSKcCYgmg/0lE5L5wjhWLnIoeL/72tzz2iyhb5OaO7We2xakW4vFDh3RLDxG5Ly8P6O/H1q1brRXpVEDc/v3vMyASZQufDxgZwa6f/cxake4ExGAQRQDHEImyiQhu+exnrRXnTkAcHoYBOMtMlGXefvNNa2W5ExCDQWz9xje8rgURpZKI1SDmTkA0Brdv2OB1LYgolYxB9aJF1oqzEhBF5C4ReU9EjovItgg/FxH5Sejnh0WkzsZ1JzEGD3zlK9aLJaL0dvTwYWtlJRwQRcQH4KcAvgBgNYCviMjqKU/7AoBVoVsDgH9P9LrTcP0hUVb60Q9/aK0sGy3EtQCOG2M+NMYMA/g1gHumPOceAL806lUAFSJSbeHa4xgQibLSj3/8Y2tl2QiItQBOTbh/OvRYrM9JjDGoq7PfEyei9GZzXYmNgBipPlOba9E8R58o0iAiLSLS0tnZGVNFHv3ud2N6PhFlOMs9QxsB8TSApRPuLwHQHsdzAADGmF3GmHpjTH1VVVX0tRDBD37wg+ifT0SZTwTPPvusteJsBMS/AFglIitEJB/AlwFMPbGxGcDfhWab1wHoMcZ0WLj2OBEcfOMNq0USUQbIsbd6MOF9bsaYURH5NoDnAPgAPGmMeUtEvhn6+c8B7AfwRQDHAVwC8A+JXnca7lAhykq333EHOi11na1s/DXG7IcGvYmP/XzC9wbAt2xca0YMiERZyeYoojs7VUTw37/6lde1IKJUEmFAjCg3Fy8cOOB1LYgoxf7Po49aK8udgJiXh92/+AUXaBNlE2OwrbHRWnHuBEQR9APAyIjXNSGiFFqycqW1stwJiNDpawZEoiwxOgrk56Ptk0+sFelUQGz60Y/0P4mI3DcyApSUWC3SqYD4qTVrgOFhr6tBRKkwMgIUF1s9w8CpgLj2ttuAYNDrahBRKoyOAqWlaG1ttVakUwFxGOAsM1G2GB4G5s1DQ0ODtSKdCohDAHesEGWToiLs3r3bWnFOBcS//fu/ZwuRKJvk51stzqmA+K+7dnldBSJKFRGgoMBqkU4FxDXr1gE+HydWiLJBMKjrENvarBXpVEA8ePAgUFYGDA15XRUiSqZAAMjNBQoKOMs8q0WLgEuXvK4FESXTwABQWQmIYMuWLdaKdSogVldX638SF2cTuW1gQBs/ljkVENvb24HSUi69IXJdIAD4/daLdSogNjU16d5GEU6sELmutBQAsHPnTmtFOhUQH3vsMZ1l9vuBwUGvq0NEyRCeUCkuBgDuVJkTJ1aI3HXpErBw4djQmFgcInMzIFZW8lxEIlclaUIFcCwgtrS06DecWCFyVzAIzJ+flKKdCohj5s3T5NWBgNc1IaJkCE2oAMCmTZusFetUQKyvr9dvcnK02zww4G2FiMiukRHdv1xUNPbQ3r17rRXvVECcZNEiBkQi14QnVCbYvHmzteLdDYhVVcyvQuSagQGgunrSQ/v27bNWvFMBsXFifla/n+OIRK4xZloL0SanAmJTU9P4ndxcYNkyoLfXs/oQkUUDA0B5ufVMexM5FRBramomP3DZZRxHJHLFxYtAhKT0xuIp+U4FxI6OjskPLFigX5lWgCjzBYPA4sXTHt5l8aR8pwLiNAUFOgDb1+d1TYgoEcPDune5vHzajx566CFrl3EqIEZMWL1iBQMiUabr6dHf5STvQHMqIEY8SryqKvUVISK7RkeBqXMESeBUQIx4DFBxsS7B4ek3RJlpdBTIy5tx/3Jzc7O1SzkVEGdMWL1yJZffEGWqixd1xYjPF/HHa9assXYppwLijBYt4gnaRJlqaAhYunTGH9fW1lq7VHYExNJSXczJU7SJMksgoC3D8BK6JHMqIM6YsFoEuOoq4Pz51FaIiBLT3Q2sWqVjiCngVECcNWH10qXc20yUSYzR474i7E6ZaOvWrdYu6VRAnDVhdUGB/qVhK5Fc09npdQ2So7dXd6ZEWIw9EXeqxGvlSl3xzq185IrDh4H77weOHPG6Jvb19elQ1xw4yxyvigr9i8MlOOSCQADYvl2/377dreGgoSFdQxxFMqmDBw9au6xTATGqhNVXXcWtfOSG3/xmfAiouxvYs8fb+th0/jzwmc/ouH8KORUQo0pYvXCh/uUZGkp+hYiSpbsbePLJ8aVkg4PAE0/o45kuENCVIbOsPZyoesoJ2olwKiBGlbDa5wNWrwYuXEh+hYiS5cCB6V3kQAB48UVv6mPT+fPA5ZcDhYVRPb29vd3apZ0KiFFbtkwnVlwac6HssmHD9K1sPh9w223e1Mem4WENiFGadFJ+grIzIBYW6n84l+BQpvL7ga9/fbwVVVgIPPigPp7J+vr0hKoYEtE/9thj1i7vVECMKWH1FVdwCQ5ltvvuGw+Afj9w773e1seGnh6dTPGIUwExpoTV8+dr15mtRMpUPh/w6KP6/bZtM54GkzH6+saXxnnEqYAYc8Lqa6/V2TmehEOZ6tprgaeeAq65xuuaJK6nB1izJualNi0tLdaq4FRAjDlhdVkZ8OlPA11dyakQUSq4cCr8hQt6InYUC7GTyamAGJerrtIW4uio1zUhyk7GAP39wHXXxZUzpb6+3lpVGBCLirS7ce6c1zUhyk5dXTrJmQYz5AkFRBHxi8jzInIs9DXiXLmInBSRIyLypojY6/BPEXfC6iuuAPLzuXuFKNUCAT3iy8OZ5YkSbSFuA3DAGLMKwIHQ/Zl83hhzvTHGXvt2iriPAcrPB264wY1tT0SZ5Nw5DYYlJXEX0djYaK06EnerCoCIvAdgvTGmQ0SqAbxkjPlUhOedBFBvjImpX1pfX29imUESkfhbiYEA8Oyz+v28efGVQUTRGx7WmeXNm6PepmeLiLRGapwl2kJcZIzpAIDQ14UzPM8A+L2ItIpIFCcweMDnA+rquMeZKFW6uoDrr084GNZYzNecO9cTROQFAJFWSn4vhuvcYoxpF5GFAJ4XkXeNMa/McL0GAA0AsGzZshguYcHixXrr6ZnzlF4iSsDAgJ46tWJFwkV1dHRYqJCas4VojLndGHN1hNtvAZwJdZUR+np2hjLaQ1/PAtgDYO0s19tljKk3xtRXxbi+KuGE1SI6ltjfz4MfiJLFGB2vr6sDcudsk6VUol3mZgBfC33/NQC/nfoEEZknIqXh7wHcAeBogteNyMpR4n4/cPXVwNmIsZ2IEtXVBSxfDixZYqW4uro6K+UAiQfE7QA2isgxABtD9yEiNSKyP/ScRQD+KCKHALwO4BljzLMJXjciawmrV6/WPZUXL9opj4jUwID2xOrq4lqEHcms2TZjlFBANMZ0GWM2GGNWhb52hx5vN8Z8MfT9h8aY60K3zxhj/q+NiidVbi5w003sOhPZFO4qr1unGyIsieqk/Chxp8pMKip0PJFdZyI7Ojs1FbCtnlzI7t27rZXlVEC0mbAagL55lZVcikOUqEuXNDf6ddd5XZNZORUQbSasBqBrE2+8UY8IGxmxWzZRtggGtat8000aFNOYUwHRZsLqMWVlQH09u85E8Tp7VlduLJxp30Zi2trarJXlVEC0mbB6kpUrddyDe52JYtPXp/uUk3h4Q9rMMmeNnBzgs5/VMxOHh72uDVFmCAR019dNNwF5eUm7zJYtW6yV5VRAtJmwepp583Q8sbOTKQeI5mIM8MknmuJgwQKvaxM1pwKizYTVEV12mTb9P/kkudchynRnz47/vmQQpwKizYTVM7r2WmDpUk6yEM2ku1sPR1m7NuaEUfHYuXOntbISOg8x2VJ6HmIshoeBAwd0OU4MCbWJnNffr78Xd92V1ueKJus8xOyUnw/89V/rWGJ/v9e1IUoPw8O6/3/9+pQGQ7G0JxpgQIxfSYm+8RcvcuaZKBDQCcdbbsmoSZSpnAqINhNWR2XBAuDmm/WDwEMgKFuFZ5Svvx5I9aHOljkVED1x2WW6P/OTT/SDQZRtzpwBLr9cj83zwKZNm6yV5VRAtJmwOiarV+tuFs48U7bp7taeUn29tfMNY7V3715rZTkVED2Tk6MfiPnz9TRgomxw8aIegPJXf+VpKoDNmzdbK4sB0Za8PJ15LizknmdyX2+vngC1fr3Vw17jsW/fPmtlORUQbSasjktREfD5z+sRRwyK5KreXl1ZcfvtzmWndCogpmSnylyKixkUyV19fcDQELBhg3PBEHAsINpMWJ2QiUHx/Hmva0NkR1+f7kK5/XZNsZEmbO5Ocyog2kxYnbBwUMzP50QLZb6LFzUYbtiQVsEQsHtSvlMBMe0UFwO33abbmM6d87o2RPG5cEE3HmzcmJZ79x966CFrZTkVEG0mrLYmPNFSUaELWIkySXe3Lq3ZuNHJMcOpnAqINo8St6qgAPjc54BFi7ijhTLHuXO6jGzDBt27nwWcCog2E1Zbl5+vC1iXLgU6Orj3mdKXMdqbKSvTIZ/iYq9rNKvm5mZrZTkVEG0mrE6K3FzNL3H11dpSHBryukZEk42OAgakuRQAAAvoSURBVG1tukf/85/XFmKas5lt07v9NtkqJ0dP3fb7gT/9SccYy8q8rhURMDCgY4Zr1wKrVnm2NzlWtbW11pbeONVCzChLlgB33qkfus5OO2XaKoeyT3e3Hna8cSNw5ZUZEwxtcyog2kxYnRIVFcAdd+hkS1tbYuOKhw8D998PHDlir37kvvB4YUmJHvsfIZn81MaXy3OCTgXEtJ1lnk1BgU62XHNN/OOKgQCwfbt+v307J2woOuHxwuXLdbwwwkxyUxPwyCPjQdAYvZ8Ou2TDtm7daq0spwKizYTVKRUeV7z1Vl0Ee/FibP/+N78Z3yLY3Q3s2WO/juSWS5e0ZXjjjTpmGCGRvDH6cdyxYzwoPvKI3r9wIX1aityp4qolS7TbIqKHzUbzievuBp58UrdVAfr1iSd4sATN7Px5nUC5445ZJ09EgMcfBx5+WINgTo5+ffhhfTxdhhltzjIzIKab8nL9oNbWAu3t44FuJgcOTO8iBwLAiy8mr46UmUZG9DNVWqp/eKuq5vwn4aA4UToFQwA4ePCgtbKcCog2E1Z7qqBAk1fdeqvO/HV2ztxa3LBBt1ZN5PPpglqisO5uPWTks58d318fhXA3eaKJY4qucSogpvVOlViJaBf67rt1d0t7u3ZzpvL7ga9/fXwBbWEh8OCD+jjR8LBOnPj9wKZN2kXOie7XfuKY4cMPaxrycPc5nYJidXW1tbKcWpgtIlbPRksLhYXAunW6c+DVV3XCpapq8of6vvuAp5/WoOn3A/fe6119KT0Yoy3CYFBzJV92Wcz9XBFdGTZxzDDcfa6oSJ9uc3t7u7WyJJ0DSH19vYkl17KTAXGi4WFdZ/jOO/qJnLhM4vBh/eT+5Ce6hIey1+CgBsPly4Ebbkh4L7Ixk4Pf1Ptea2pqivm0fBFpNcZMS9PJgJiJOju1tdjXpwtpw2OInZ1RDZSTo4zRz4DPp8tpamvTK3IlSTy/9zMFRKe6zDYTVqe1qiqdJXznHeDoUR0gLy9nMMxmfX1ATw9wxRXAdddlxKEM6cipgGgzYXXay8vTxdxLlgCHDunAeWkpD4rINpcu6bpCv19znSxa5HWNMppTs8yxJKx2Zn+m36+5cTdu1JNz2tq0tUBuGxjQSbRgUJdn3Xln1gbDWIbV5uJUCzHahNVNTbr1KDxzFl5eUFGRXns0oyaiY4m3367bsd54QwPj/Plpf7gnxWhoSNcUFhfrHvja2unrUCluTrUQo5Ep+zPjIgIsXqythVtv1R0rM61fpMwyPKwnrff364TJ3XcDy5YxGAKor582NxI3p1qI0Zi4lmrHDr0B6bc/MyE5OTq2WF2tLcU33tBo7/frLhjKHCMjuoQmNxeor9elNBEOYiA7nFp2EwtjJq9tDgYdCYaRjI4CH38MvPmmdrlKS7MmaVDGGhjQP2I+n64rvfxyzctD09hcduNUlznaY4CybX8mcnOBlSuBzZt13KmgQLvSnZ08OzGdBAI6PhieLLnxRuCee4CrrmIwnEVjY6O1spxqIUbzl2Lq/szHH59+39mW4kQXLgAnTgDHjukvYllZ1Bv+ybKBAV1DCAArVmhrcMGCLPkgeiMrFmZHI1P2ZyZdRYVu67r6ah2sf/ttbZnk5elYIwfrkysY1D9Kg4M6hFFfr+O+XFAds5qaGmv7mbOuhRiW7vszUy48/f7hh8AHH4y3GouLs/w/xrJwa9AYnSC54gqgspL/xwng1r0ZxJKweurnL+s/jyK6bnHNGh3Eb28H3n9f87wAGhhLSthyjFUwqAvl+/v1fnm5/h8vWaIL6SmtOBUQbR4lntXy87X1sny5tmi6uoBTp/QWCOgkTVkZl/DMZGREj2kbHtalDNXVur+4spLjtElQV1dnrays7TJTHAIB3Tfb0aETMv392rIsKdFf9GxtZhujfzh6e7VFWFCgf0xqa3U8lusG005WdJkpyXw+beVUVupkTF+fLt05eVK71iJ6KyrSm6uBIBDQQxUGBjQAGqOB74YbdAtleXn2/nHwQENDg7XMewm1EEXkbwA0AbgKwFpjTMTmnIjcBWAHAB+AXxhjtkdTPluIGWRoSCcLLlzQINnZqQEjHBjCQTLT1tONjOjrGBgYX6iam6uBb+FCna0vL+d4oIfSaVLlKID7AMyY3UlEfAB+CmAjgNMA/iIizcaYtxO89jQ2E1ZTjAoKxoPElVfqY4OD4+f0dXZqatWurvEgmZ+vt9xcbU16NWETCOhunpERvQ0Njf8s/LqqqsZPKS8qYgvQUQkFRGPMO4BG6FmsBXDcGPNh6Lm/BnAPAOsB0WbCarKgsFBvlZW62BjQiYbeXr11dWnXs69PW5bDw+Pd7ol/8fPyxm/hn0+9Afpvwrfw/WBQg1w44IV/Fr6GMVruvHm6HrC4WLu/ZWV6n+sCs4sxJuEbgJcA1M/wsy9Bu8nh+18F8G+zlNUAoAVAy7Jly0ws6urqDICxW0tLi2lpaZn0WGNjozHGmOrq6rHH6urqjDHGbN26ddJz29raTHNz86THdu7caYxWdOy2adMmY4wxmzZtmvS4Mcbs3Llz0mPNzc2mra1t0mNbt26dVv/q6mpjjDGNjY1Z85pyAFMImDLA/OA73zHm1CmzrrLSXAOYdYD55hVXGLN/v9lx553mbsBsBswWwJzbvdu8vm2b2RJ67G7A7P/Wt4zZv9/cDpibAXMdYBrWrzfm44/N327caMoBUwQYH98nJ15TrAC0mAjxZ84xRBF5AcDiCD/6njHmt6HnvATgOybCGGJonPFOY8w3Qve/Ch1v/N+zXhjJPdyBiLJX3GOIxpjbE7z2aQBLJ9xfAsBe3kAiIktScdrNXwCsEpEVIpIP4MsAot9SQkSUIgkFRBG5V0ROA7gJwDMi8lzo8RoR2Q8AxphRAN8G8ByAdwA8ZYx5K7FqExHZl+gs8x4AeyI83g7gixPu7wewP5FrERElm1MHxBIRJYIBkYgohAGRiCiEAZGIKIQBkYgohAGRiCiEAZGIKCStT8wWkU4AH8XwTyoBnEtSdTJBNr9+vvbsFc/rv8wYUzX1wbQOiLESkZZIG7azRTa/fr727HztgN3Xzy4zEVEIAyIRUYhrATHbj8zO5tfP1569rL1+p8YQiYgS4VoLkYgobhkdEEXkb0TkLREJisiMs0wicpeIvCcix0VkWyrrmEwi4heR50XkWOjr/Bmed1JEjojImyKS0TkZ5novRf0k9PPDIlLnRT2TIYrXvl5EekLv85si8n0v6pkMIvKkiJwVkaMz/NzO+x4p0Uqm3KD5oD+F2ZNc+QB8AGAlgHwAhwCs9rrull7/DwFsC32/DcAPZnjeSQCVXtfXwuud872EnsP5OwACYB2A17yudwpf+3oA+7yua5Je/+cA1AE4OsPPrbzvGd1CNMa8Y4x5b46njaVBNcYMAwinQXXBPQD+M/T9fwL4Xx7WJRWieS/vAfBLo14FUCEi1amuaBK4/DmekzHmFQDdszzFyvue0QExSrUATk24fzr0mAsWGWM6ACD0deEMzzMAfi8irSLSkLLa2RfNe+nq+x3t67pJRA6JyO9E5DOpqVpasPK+J5RCIBWiSYM6VxERHsuYqfXZXn8MxdxijGkXkYUAnheRd0N/cTNNNO9lRr/fs4jmdR2EbknrE5EvAngawKqk1yw9WHnf0z4gmixPgzrb6xeRMyJSbYzpCHUPzs5QRnvo61kR2QPtfmViQIzmvczo93sWc74uY8zFCd/vF5GfiUilMSYb9jlbed+zocvschrUZgBfC33/NQDTWswiMk9ESsPfA7gDQMSZugwQzXvZDODvQrOO6wD0hIcVMtycr11EFouIhL5fC/397kp5Tb1h5333evYowZmne6F/GYYAnAHwXOjxGgD7p8xAvQ+dpfue1/W2+PoXADgA4Fjoq3/q64fOSh4K3d7K9Ncf6b0E8E0A3wx9LwB+Gvr5Ecyw+iATb1G89m+H3uNDAF4FcLPXdbb42n8FoAPASOh3/sFkvO/cqUJEFJINXWYioqgwIBIRhTAgEhGFMCASEYUwIBIRhTAgEhGFMCASEYUwIBIRhfx/L3kNGSCF828AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ax.scatter([-0.9, .7, -0.7], [-0.9, -0.8, .9], marker='x', color='b')\n",
    "ax.scatter([-0.5, .5, -0.3], [-0.8, -0.4, .4], marker='d', color='r')\n",
    "\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can repeat this procedure for many iterations, counting how many times the points fall into the circle, in relation to the total number of throws.\n",
    "We can mathematically prove that the ratio of the inside-count (*hints*) and the total-sample-count is an estimate of the ratio of the two areas, i.e\n",
    "\n",
    "$$\n",
    "\\frac{A_{\\mbox{circle}}}{A_{\\mbox{square}}} = \\frac{\\pi r^2}{l^2} = \\frac{\\pi}{4} \\rightarrow \\pi = 4 \\times \\frac{N_\\mbox{hints}}{N_{\\mbox{points}}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step1 - C++ like version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets start from a very naive implementation of this algorithm, following step-by-step the pseudo-code instructions.\n",
    "\n",
    "This scheme of programming is very close to the C++ language, in which we have a \"low-control\" of the execution and a capillar management of the objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pi with 100000 steps for 500 runs is 3.141358 in 37.608s\n"
     ]
    }
   ],
   "source": [
    "N_run = 500 # number of runs\n",
    "N = 100000 # number of MC events\n",
    "Nhits = 0.0 # number of points accepted\n",
    "pi = np.zeros(shape=(N_run,)) # values of pi\n",
    "\n",
    "tic = now() # start clock \n",
    "for I in range(N_run):\n",
    "    Nhits = 0.0\n",
    "    for i in range(N):\n",
    "        x = np.random.rand()*2 - 1\n",
    "        y = np.random.rand()*2 - 1\n",
    "        res = x*x + y*y\n",
    "        if res < 1:\n",
    "            Nhits += 1.0\n",
    "    pi[I] += 4. * Nhits/N\n",
    "\n",
    "toc = now()\n",
    "\n",
    "print ('pi with {:d} steps for {:d} runs is {:.6f} in {:.3f}s'.format(\n",
    "    N, N_run, np.mean(pi), toc-tic))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1b1a422f3d0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT8AAAEvCAYAAAAzcMYwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAZcElEQVR4nO3dfZAV9Z3v8feXJ2VmeBoZkDuDSwgEYQGDNatBI2uJeyPgFbKiIYkpzM0tyiTuze5NahncSmD/EOHeuNnU4iaFISuVjVKsirIqMcguC3ej6IgjaAbW0aCOwjxAeH4amO/+MQ074DAw5+nX0/15VU2dc7pPd39sjp/p3zndZ8zdERFJmx6hA4iIhKDyE5FUUvmJSCqp/EQklVR+IpJKKj8RSaVeoQMADB482EeMGBE6hsTQzp07ARgzZkzgJNIdvf76683uXtbRvFiU34gRI6iurg4dQ2JowYIFADz00EOBk0h3ZGbvX2heLMpP5EJUepIves9PRFJJ5Sexduedd3LnnXeGjiEJpGGvxNrevXtDR5CE0pGfiKSSjvxEIgcPHqSxsZGWlpbQUeQSFRcXU1FRQY8eXT+OU/mJ0FZ8DQ0NlJeX07dvX8wsdCS5iNbWVj766COam5sZMmRIl5dX+UmsTZ06tSDbaWxspLy8nKKiooJsT7LXo0cPhg4dyvvvv6/yk+T5/ve/X5DttLS00Ldv34JsS3Knd+/enDp1KqNl9YGHSERD3e4nm38zlZ/E2rRp05g2bVroGJJAGvZKzoyoej7n69zzdv0n1r1ryYycb0fSR0d+InKORx55hIkTJ9K/f3/69+/P5MmTef75zn+xbdq0iTvuuIPy8nLMjMcee6zT5y9evBgz4/77789qPdlQ+YnIOSoqKli6dClbt26lurqaW265hVmzZrFt27YLLnP48GHGjx/Pj3/844t+cPTKK6/w6KOPMnHixKzWky2Vn0g3N27cOMysw59FixZ1eX0zZ85k2rRpjBo1is985jM8+OCD9OvXj5dffvmCy0yfPp3Fixcze/bsTk84PnDgAF/96ldZsWIFgwYNyng9uaDyk1jr++nr6Pvp60LHiLVnnnkGgBdeeIHdu3fz8ccfU1RUxIoVK5g/fz6LFy+mpKSk05/Nmzd3uO7Tp0+zatUqDh8+zA033JB11nnz5jF79mxuueWWrNeVLX3gIbE24Po/Dbr9m2+++RPT7r77br71rW9x9OhRpk+f/on59957L/feey/Nzc3Mnj37E/O/+c1v8qUvfYkPP/yQr33ta+fM27hxY5czNjQ0YGbcdNNNlJSUUFdXx9GjR/n85z9P3759ue+++7j77rs7XUd5efk5j7dv387kyZM5fvw4JSUlrFmzhgkTJnQ5W3uPPvoodXV1/OIXv8hqPbmi8hPp5t58801GjhxJSUkJADU1NRQVFTFq1CgASktLKS0t7dI6x4wZQ01NDfv37+epp55i7ty5bNy4kfHjx2eUcefOnTzwwANs3ryZPn36ZLSOXFP5SaztebwKgCu/siTI9js7EisqKup0/uDBgzudP3z48IyO9M63bds2rrnmmrOPa2pqGD9+/Nn3zBYvXszixYs7Xce6deu46aabzj7u06fP2fKsrKzktdde40c/+hErVqzIKOPLL79Mc3PzOeV5+vRpNm3axE9/+lOOHDnCZZddltG6M6XyE+nmtm3bds6J4DU1NeeUYSbD3vO1trZy4sSJjDPOmjWLysrKc6Z9/etfZ/To0TzwwANBjgZVfiLdWGtrK2+99Rbz588/O+3dd99l8uTJZx93ddhbVVXFjBkzGD58OIcOHeLxxx9n48aN55zrt2zZMpYtW8aOHTuAtlNU6urqzmb64IMPqKmpobS0lKuuuoqBAwcycODAc7ZTXFxMaWnpOUeDF1tPLunTXpFu7N133+XIkSPnnDM3YcIEli5dyrp16zJa5549e7jnnnsYM2YMU6dO5bXXXmPdunXnHF02Nzef/bOiANXV1UyaNIlJkyZx7NgxFi5cyKRJk/jBD37QpW3naj2Xwtw95yvtqsrKStefruz+8nJ5Wwfv+eXj8rba2lrGjh2b8/VK/nX2b2dmr7t7ZUfzNOyVWCu++qaLP0kkAyo/ibV+1+pLDCQ/9J6fxFpry3FaW46HjiEJpCM/ibXGf1oEhDvPT5JL5SfdTj4+WFl+x5Wc/PD3Z78ZeGLFwIssIXGQzQe2GvaKAPuPt8Kpk6FjSBe1tLTQq1dmx3AqPxHg8W37aWzYjbecyOpoQgqntbWVhoYGBgwYkNHyGvaKADV7TvLIK018ZWILAy/vQe1h/QnL7qC4uJjBgwdntKzKT2KtZMKtBdtWzZ6T1OxpBPR3QtLgosNeM/u5mTWa2Vvtpv0/M9thZtvMbI2ZDWw3b4GZ1ZnZTjP7Qr6CSzqUTLi1oAUo6XEp7/k9Btx23rT1wHh3nwj8B7AAwMzGAXOAP4yW+Xsz65mztJI6p48e4PTRA6FjSAJdtPzcfROw77xpv3b3M38m/RWgIro/E1jl7ifc/XdAHaDvIJeMNT3zEE3PPBQ6hiRQLj7t/Z/Ama+PKAc+bDevPpomIhIrWZWfmf0VcAr45ZlJHTytw/MGzGyemVWbWXVTU1M2MUREuizj8jOzucDtwFf9v06MqgeGt3taBfBxR8u7+3J3r3T3yrKyskxjiIhkJKPyM7PbgPnAHe5+tN2stcAcM7vMzD4FjAZezT6miEhuXfQ8PzN7ArgZGGxm9cBC2j7dvQxYH10L+Yq73+fub5vZauC3tA2Hv+3up/MVXpKv36RP/mlIkVy4aPm5+5c7mHzBP+Hk7g8CD2YTSuSM4rFTQkeQhNK1vRJrpw42ceqgPhCT3FP5Saw1P/cwzc89HDqGJJDKT0RSSeUnIqmk8hORVFL5iUgq6fv8JNb6X/fF0BEkoVR+EmtFo64PHUESSsNeibWWvfW07K0PHUMSSOUnsbb3xWXsfXFZ6BiSQCo/EUkllZ+IpJLKT0RSSeUnIqmkU10k1gbcMCd0BEkolZ/EWt8Rnw0dQRJKw16JtZMN73Gy4b3QMSSBVH4Sa/s2LGffhuWhY0gCqfxEJJVUfiKSSio/EUkllZ+IpJJOdZFYGzhlbugIklAqP4m1yyvGho4gCaVhr8Ta8fpajtfXho4hCaTyk1jbv2kl+zetDB1DEkjlJyKppPITkVRS+YlIKqn8RCSVLlp+ZvZzM2s0s7faTSs1s/Vm9k50O6jdvAVmVmdmO83sC/kKLulQOnUepVPnhY4hCXQpR36PAbedN60K2ODuo4EN0WPMbBwwB/jDaJm/N7OeOUsrqdNn6Ej6DB0ZOoYk0EXLz903AfvOmzwTOHP+wUpgVrvpq9z9hLv/DqgDrstRVkmhY7tqOLarJnQMSaBMr/AY6u67Adx9t5kNiaaXA6+0e159NE0kIwd+swrQNzpL7uX6Aw/rYJp3+ESzeWZWbWbVTU1NOY4hItK5TMuvwcyGAUS3jdH0emB4u+dVAB93tAJ3X+7ule5eWVZWlmEMEZHMZFp+a4EzX7cxF3i23fQ5ZnaZmX0KGA28ml1EEZHcu+h7fmb2BHAzMNjM6oGFwBJgtZl9A/gAuAvA3d82s9XAb4FTwLfd/XSesouIZOyi5efuX77ArKkXeP6DwIPZhBI544ov3B86giSUvs9PYq33FRWhI0hC6fI2ibWjdVs4WrcldAxJIB35SawdfHUNAEWjrg+cRJJGR34ikkoqPxFJJZWfiKSSyk9EUkkfeEisDb79u6EjSEKp/CTWevXXdd+SHxr2Sqwdqd3EkdpNoWNIAunIT2Lt0BsvAFA8dkrgJJI0OvITkVRS+YlIKmnYmxIjqp4PHUEkVnTkJyKppCM/ibWyWQtCR5CEUvlJrPUsGhA6giSUhr0Sa4e3v8Th7S+FjiEJpPKTWFP5Sb6o/EQklVR+IpJKKj8RSSWVn4ikkk51kVgbctei0BEkoVR+Ems9el8eOoIklIa9EmuHtj7Poa26LllyT+UnsXZkx2aO7NgcOoYkkMpPRFJJ5SciqaTyE5FUyqr8zOwvzOxtM3vLzJ4ws8vNrNTM1pvZO9HtoFyFFRHJlYxPdTGzcuB/A+Pc/ZiZrQbmAOOADe6+xMyqgCpgfk7SSupc+ZUloSNIQmU77O0F9DWzXkAR8DEwE1gZzV8JzMpyGyIiOZdx+bn7R8APgQ+A3cABd/81MNTdd0fP2Q0MyUVQSacDW57mwJanQ8eQBMq4/KL38mYCnwL+G1BsZvd0Yfl5ZlZtZtVNTU2ZxpCEO/buqxx799XQMSSBshn23gr8zt2b3L0FeBq4AWgws2EA0W1jRwu7+3J3r3T3yrKysixiiIh0XTbl9wHwOTMrMjMDpgK1wFpgbvScucCz2UUUEcm9jD/tdfctZvYksBU4BbwBLAdKgNVm9g3aCvKuXAQVEcmlrL7Vxd0XAgvPm3yCtqNAkaxZr8tCR5CE0ldaSawNvfuvQ0eQhNLlbSKSSio/ibX9//4E+//9idAxJIFUfhJrx99/k+Pvvxk6hiSQyk9EUknlJyKppPITkVTSqS4Saz379g8dQRJK5SexVvbFB0JHkITSsFdEUknlJ7H2+397jN//22OhY0gCadgrsXbiox2hI0hC6chPRFJJ5SciqaTyE5FU0nt+Emu9+g0OHUESSuUnsTb4f3wvdARJKA17RSSVVH4Sa/teWs6+l5aHjiEJpGGvxNrJxvdCR5CEUvmJdGBE1fMF2c6uJTMKsh35JA17RSSVVH4ikkoa9kqs9S4tDx1BEkrlJ7F2xW1/FjqCJJSGvSKSSio/ibW9v/o79v7q70LHkATSsFdirWXfR6EjSELpyE9EUknlJyKplFX5mdlAM3vSzHaYWa2ZTTazUjNbb2bvRLeDchVWRCRXsj3y+zHwK3e/GrgGqAWqgA3uPhrYED0WyUifISPpM2Rk6BiSQBl/4GFm/YEpwL0A7n4SOGlmM4Gbo6etBDYC87MJKelVeuu80BEkobI58hsJNAH/YGZvmNnPzKwYGOruuwGi2yE5yCkiklPZlF8v4FrgJ+4+CThCF4a4ZjbPzKrNrLqpqSmLGJJkzf/8Q5r/+YehY0gCZVN+9UC9u2+JHj9JWxk2mNkwgOi2saOF3X25u1e6e2VZWVkWMSTJTh1q5tSh5tAxJIEyLj933wN8aGZjoklTgd8Ca4G50bS5wLNZJRQRyYNsr/D4M+CXZtYHeA/4Om2FutrMvgF8ANyV5TZERHIuq/Jz9xqgsoNZU7NZr4hIvunaXom1y8qvDh1BEkrlJ7E26I/vDR1BEkrX9opIKqn8JNaa1iymac3i0DEkgTTslVg7fexg6AiSUDryE5FUUvmJSCqp/EQklfSen8Ta5X9wTegIklAqP4m1gTd+OXQESSgNe0UklVR+EmsNqxfSsHph6BiSQBr2Sqz5qROhI0hC6chPRFJJ5SciqaTyE5FU0nt+Emt9P31d6AiSUCo/ibUB1/9p6AiSUBr2ikgqqfwk1vY8XsWexy/5z0GLXDKVn4ikkspPRFJJ5SciqaTyE5FU0qkuEmvFV98UOoIklMpPYq3ftTNCR5CE0rBXYq215TitLcdDx5AE0pGfxFrjPy0C4MqvLAkbRBJHR34ikkoqPxFJpazLz8x6mtkbZvZc9LjUzNab2TvR7aDsY4qI5FYujvy+A9S2e1wFbHD30cCG6LGISKxkVX5mVgHMAH7WbvJMYGV0fyUwK5ttSLqVTLiVkgm3ho4hCZTtp71/C/wl0K/dtKHuvhvA3Xeb2ZAstyEppuKTfMn4yM/Mbgca3f31DJefZ2bVZlbd1NSUaQxJuNNHD3D66IHQMSSBshn23gjcYWa7gFXALWb2j0CDmQ0DiG4bO1rY3Ze7e6W7V5aVlWURQ5Ks6ZmHaHrmodAxJIEyLj93X+DuFe4+ApgD/Iu73wOsBeZGT5sLPJt1ShGRHMvHeX5LgD8xs3eAP4kei4jESk4ub3P3jcDG6P5eYGou1isiki+6wkNEUklfbBADI6qeDx0htvpNmh46giSUyk9irXjslNARJKE07JVYO3WwiVMHdR6o5J7KT2Kt+bmHaX7u4dAxJIFUfiKSSio/EUkllZ+IpJLKT0RSSae6SKz1v+6LoSNIQqn8JNaKRl0fOoIklIa9Emste+tp2VsfOoYkkMpPYm3vi8vY++Ky0DEkgVR+IpJKKj8RSSWVn4ikkspPRFJJp7pIrA24YU7oCJJQKj+Jtb4jPhs6Ql4V4otsdy2ZkfdtdEca9kqsnWx4j5MN74WOIQmk8pNY27dhOfs2LA8dQxJI5SciqaTyE5FUUvmJSCqp/EQklXSqi8TawClzQ0eQhFL5SaxdXjE2dARJKA17JdaO19dyvL42dAxJIJWfxNr+TSvZv2ll6BiSQCo/EUmljMvPzIab2b+aWa2ZvW1m34mml5rZejN7J7odlLu4IiK5kc2R3yngu+4+Fvgc8G0zGwdUARvcfTSwIXosIhIrGZefu+92963R/UNALVAOzATOvEmzEpiVbUgRkVzLyakuZjYCmARsAYa6+25oK0gzG5KLbUg6lU6dFzqCJFTW5WdmJcBTwJ+7+0Ezu9Tl5gHzAK666qpsY0hC9Rk6MnQESaisPu01s960Fd8v3f3paHKDmQ2L5g8DGjta1t2Xu3ulu1eWlZVlE0MS7NiuGo7tqgkdQxIom097DVgB1Lr737SbtRY4c03SXODZzONJ2h34zSoO/GZV6BiSQNkMe28EvgZsN7Mzv5ofAJYAq83sG8AHwF3ZRRQRyb2My8/d/z9woTf4pma6XhGRQtAVHiKSSio/EUklfaWVxNoVX7g/dARJKJWfxFrvKypCR5CE0rBXYu1o3RaO1m0JHUMSSEd+EmsHX10DQNGo6wMnkaTRkZ+IpJLKT0RSSeUnIqmk8hORVNIHHhJrg2//bugIklAqP4m1Xv31dWeSHyq/Toyoej50hNQ7UrsJgOKxUwInkaRR+UmsHXrjBUDlJ7mnDzxEJJVUfiKSSio/EUkllZ+IpJI+8JBYK5u1IHQESSiVn8Raz6IBoSNIQmnYK7F2ePtLHN7+UugYkkAqP4k1lZ/ki8pPRFJJ5SciqaTyE5FUUvmJSCrpVBeJtSF3LQododsrxLcT7VoyI+/byDWVn8Raj96Xh44gCaXyk1g7tLXtqKXftd3vyCJNCvXdl7k8wtR7fhJrR3Zs5siOzaFjSALlrfzM7DYz22lmdWZWla/tiIhkIi/lZ2Y9gUeAacA44MtmNi4f2xIRyUS+3vO7Dqhz9/cAzGwVMBP4ba42oL+vISLZyNewtxz4sN3j+miaiEgs5OvIzzqY5uc8wWweMC96eNjMduYpy2CgOU/r7oq45ID4ZLnkHO8vvT3PUbrfPimAuGQ5m8OWdnnZP7jQjHyVXz0wvN3jCuDj9k9w9+XA8jxt/ywzq3b3ynxvp7vkgPhkiUsOiE+WuOSA+GTJV458DXtfA0ab2afMrA8wB1ibp22JiHRZXo783P2Umd0PvAj0BH7u7m/nY1siIpnI2xUe7v4C8EK+1t8FeR9aX6K45ID4ZIlLDohPlrjkgPhkyUsOc/eLP0tEJGF0eZuIpFK3KT8zu9zMXjWzN83sbTP76w6ec7WZvWxmJ8zse+2mDzezfzWz2mjZ77Sbt8jMPjKzmuhnep6zXHBZMys1s/Vm9k50OyhQjoLuk3bze5rZG2b2XIh9cpEcBd8nZrbLzLZH26sOtU86yRFinww0syfNbEf0//PkTPYJAO7eLX5oO3ewJLrfG9gCfO685wwB/gh4EPheu+nDgGuj+/2A/wDGRY8XtX9uAbJccFng/wJV0f0qYGmgHAXdJ+3m/x/gceC5dtMKtk8ukqPg+wTYBQzuYL0F3Sed5AixT1YC/yu63wcYmMk+cffuc+TnbQ5HD3tHP37ecxrd/TWg5bzpu919a3T/EFBLFlecZJmls2Vn0vaPS3Q7K1COLssmC4CZVQAzgJ+dN6tg++QiObos2yydKOg+yaVssphZf2AKsCJ63kl33x/N7tI+gW407IWzw5EaoBFY7+5bMljHCGASbb9xzrjfzLaZ2c8v6XA5yyydLDvU3XdDW2HT9hswRA4o8D4B/hb4S6D1vOkF3Sed5IDC7xMHfm1mr1vbFVFnFHqfXCgHFHafjASagH+I3pb4mZkVR/O6vE+6Vfm5+2l3/yxtV4xcZ2bju7K8mZUATwF/7u4Ho8k/AT4NfBbYDTyc7yzZ/ncUIEdB94mZ3Q40uvvrl5o9QI6Cv06AG939Wtq+HenbZjalC8sWIkeh90kv4FrgJ+4+CThC2xA3I92q/M6IDnU3Ardd6jJm1pu24vuluz/dbl0N0T9GK/Aobd9Ik9csnSzbYGbDorzDaPvNWPAcAfbJjcAdZrYLWAXcYmb/GM0r5D65YI4QrxN3/zi6bQTWtNtmQV8nF8oRYJ/UA/XtjhSfpK0MIYN90m3Kz8zKzGxgdL8vcCuw4xKXNdreJ6h19785b96wdg+/CLyV5yydLbsWmBvdnws8GyJHofeJuy9w9wp3H0HbpZD/4u73RLMLtk86yxHgdVJsZv3O3Af+e7ttFvJ1csEcAV4ne4APzWxMNGkq//U1eV3aJ2dW2C1+gInAG8C2aCf/IJp+H3BfdP9K2n47HAT2R/f7A5+n7X2LbUBN9DM9WuYXwPZo3lpgWJ6zdLhstMwVwAbgnei2NFCOgu6T89ZzM+d+ylqwfXKRHIV+nYwE3ox+3gb+KtDrpLMcBX+d0DbEro6WfwYYlMk+cXdd4SEi6dRthr0iIrmk8hORVFL5iUgqqfxEJJVUfiKSSio/EUkllZ+IpJLKT0RS6T8BHcLWIInF1t0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 5))\n",
    "ax.hist(pi, bins=10, align='mid')\n",
    "ax.axvline(np.mean(pi), linestyle='dashed', color='k', label=r'$\\pi$={:.3f}'.format(np.mean(pi)))\n",
    "ax.legend(loc='best', fontsize=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step2 - A more Pythonic version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case we can remember that the NumPy library provides a lot of functions for the random number generators able to perform a faster and more accurate computation.\n",
    "\n",
    "Moreover, is very (**very**) disadvantageous to use *for* loops in Python, and the use of pre-built functions is **always recommended**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pi with 100000 steps for 500 runs is 3.141406 in 1.066s\n"
     ]
    }
   ],
   "source": [
    "pi = np.empty(shape=(N_run,)) # values of pi\n",
    "\n",
    "tic = now() # start clock \n",
    "for I in range(N_run):\n",
    "    points = np.random.uniform(low=-1, high=1, size=(2, N))\n",
    "    areas = np.sum(points**2, axis=0)\n",
    "    Nhints = np.mean(areas < 1)\n",
    "    pi[I] = Nhints * 4\n",
    "\n",
    "toc = now()\n",
    "\n",
    "print ('pi with {:d} steps for {:d} runs is {:.6f} in {:.3f}s'.format(\n",
    "    N, N_run, np.mean(pi), toc-tic))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1b1a43645e0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT8AAAEvCAYAAAAzcMYwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAbS0lEQVR4nO3dfXRV9Z3v8fdXDJrwIEQSpAlOVBCxqMWVYtHqeMW5CjqCFS22dqHXu1i22tuZaVcFu1qYP0CcW9vpGtrpwuKV6VQZpj4xPtQqHQp3qmjECGJgjEoxgCRgeUxAHr7zRzbMIebx7HPO3vH3ea2Vdc7Zjx83+OG3z9n7xNwdEZHQnJR0ABGRJKj8RCRIKj8RCZLKT0SCpPITkSCp/EQkSCcnHQBgyJAhXlVVlXQMCcDGjRsBGDVqVMJJpBBef/31He5e1t68VJRfVVUVNTU1SceQAMyaNQuA+++/P+EkUghm9seO5qWi/EQKRaUnx+g9PxEJkspPgnLTTTdx0003JR1DUkCnvRKUnTt3Jh1BUqLLkZ+ZPWxmjWb2Vpvp3zSzjWa23sz+LmP6LDOrj+Zdk4/QIiJxdWfk9wiwAPinYxPM7H8Ak4EL3f2gmZVH088HpgGfBT4DvGRm57r7kTgh9+zZQ2NjI4cOHYqzGSmgoqIiysvLGThwYNJRRNrVZfm5+0ozq2oz+evAfHc/GC3TGE2fDCyJpr9vZvXAOODlbAPu2bOH7du3U1FRQXFxMWaW7aakQNydlpYWtmzZAqAClFTK9gOPc4HLzWy1mf3ezD4fTa8APshYriGalrXGxkYqKiooKSlR8fUSZkZJSQkVFRU0NjZ2vUIBTZgwgQkTJiQdQ1Ig2w88TgYGA18APg8sNbOzgfbaqd1vSzWzGcAMgDPPPLPDHR06dIji4uIsY0qSiouLU/dWxfe///2kI0hKZDvyawCe8FavAkeBIdH04RnLVQJb29uAuy9092p3ry4ra/fuk+M04uud9OcmaZZt+T0FXAVgZucCfYEdwDJgmpmdYmZnASOBV3MRVCQXJk6cyMSJE5OOISnQ5WmvmT0GXAkMMbMGYDbwMPBwdPnLx8B0b/1lIOvNbCnwNnAYuDvuJ70ibVXNfDbrdT9c39CtbWyaf13W+5DeocuRn7vf6u7D3L3I3SvdfZG7f+zut7n7GHe/2N1/l7H8XHc/x91Hufvz+Y0flp/+9KdceOGFDBw4kIEDBzJ+/Hiefbbz/4lXrlzJDTfcQEVFBWbGI4880uny8+bNw8y45557Ym1HJO10e1svUllZyQMPPMCaNWuoqanhqquuYsqUKaxdu7bDdfbt28eYMWP4yU9+0uUHR6+88goPPfQQF154YaztiPQGKr88Ov/88zGzdn/mzJnT4+1NnjyZiRMnMmLECM4991zmzp3LgAEDePnlji+jnDRpEvPmzWPq1KmcdFLHf9y7d+/mq1/9KosWLWLw4MFZb0ekt9Df4jx66qmnAHjuuefYtm0bW7dupaSkhEWLFnHvvfcyb948+vfv3+nPqlWr2t32kSNHWLJkCfv27ePSSy+NnXXGjBlMnTqVq666Kva20qz4nHEUnzMu6RiSAr32iw2uvPLKT0y75ZZb+MY3vkFzczOTJk36xPzbb7+d22+/nR07djB16tRPzP/617/Ol7/8ZT744AO+9rWvnTBvxYoVPc64fft2zIzLL7+c/v37U19fT3NzM1/84hcpLi7mrrvu4pZbbul0GxUVJ14jvm7dOsaPH8+BAwfo378/Tz75JBdccEGPs2V66KGHqK+v55e//GWs7fQGp13ypaQjSEr02vLrDd58803OPvts+vfvD0BtbS0lJSWMGDECgNLSUkpLS3u0zVGjRlFbW8uuXbt4/PHHmT59OitWrGDMmDFZZdy4cSP33Xcfq1atom/fvlltQ6Q36rXl19lIrKSkpNP5Q4YM6XT+8OHDsxrptbV27Vouuuii469ra2sZM2bM8ffM5s2bx7x58zrdxvPPP8/ll19+/HXfvn2Pl2d1dTWvvfYaP/7xj1m0aFFWGV9++WV27NhxQnkeOXKElStX8vOf/5z9+/dzyimnZLXtNPrw0ZkAnPGV+QknkaT12vLrDdauXXvCBbW1tbUnlGE2p71tHT16lIMHD2adccqUKVRXV58w7Y477mDkyJHcd999Gg3Kp5bKL0+OHj3KW2+9xb333nt82rvvvsv48eOPv+7pae/MmTO57rrrGD58OHv37uXRRx9lxYoVJ1zrt2DBAhYsWMCGDRuA1ktU6uvrj2favHkztbW1lJaWcuaZZzJo0CAGDRp0wn769etHaWnpCaPBrrYj0tvo0948effdd9m/f/8J18xdcMEFPPDAAzz/fHbXfn/44YfcdtttjBo1igkTJvDaa6/x/PPPnzC63LFjx/FfzwhQU1PD2LFjGTt2LC0tLcyePZuxY8fygx/8oEf7ztV2RNLCWu9KS1Z1dbV39Ksr6+rqGD16dIETSa7k488v1u1t3XzPT7e3fTqY2evuXt3ePJ32SlD6nXd51wtJEFR+EpQBF2tEJ630np8E5eihAxw9dCDpGJICGvlJUBr/dQ6g6/ykl4z80vChjPSc/twkzVJffkVFRbS0tCQdQ7LQ0tJCUVFR0jFE2pX68isvL2fLli00NzdrJNFLuDvNzc1s2bKF8vLypOOItCv17/kd+52vW7duTd1vApOOFRUVMXToUP3OXkmt1JcfcPxr20Xi6n/B1UlHkJToFeUnkisqPzkm9e/5ieTSkebdHGnenXQMSQGVnwSl6an7aXrq/qRjSAqo/EQkSCo/EQlSl+VnZg+bWaOZvdXOvO+YmZvZkIxps8ys3sw2mtk1uQ4sIpIL3Rn5PQJc23aimQ0H/gLYnDHtfGAa8NlonZ+ZWZ+cJBURyaEuy8/dVwIftTPrx8B3gczbLiYDS9z9oLu/D9QD+iWpkhoDxk5iwNhP/lpTCU9W1/mZ2Q3AFnd/08wyZ1UAr2S8boimtbeNGcAMQL8DQgqm3+grko4gKdHjDzzMrAT4HtDeL2+wdqa1e0Ouuy9092p3ry4rK+tpDJGsHN7TxOE9TUnHkBTI5tPec4CzgDfNbBNQCawxszNoHekNz1i2EtgaN6RIrux45kF2PPNg0jEkBXpcfu6+zt3L3b3K3atoLbyL3f1DYBkwzcxOMbOzgJHAqzlNLCKSA9251OUx4GVglJk1mNmdHS3r7uuBpcDbwG+Au939SK7CiojkSpcfeLj7rV3Mr2rzei4wN14sEZH80h0eIhIkfaWVBGXguBuTjiApofKToJSMuCTpCJISOu2VoBza2cChnQ1Jx5AUUPlJUHa+sICdLyxIOoakgMpPRIKk8hORIKn8RCRIKj8RCZIudZGgnHbptKQjSEqo/CQoxVWfSzqCpIROeyUoH29/j4+3v5d0DEkBlZ8E5aPlC/lo+cKkY0gKqPxEJEgqPxEJkspPRIKk8hORIOlSFwnKoCumJx1BUkLlJ0E5tXJ00hEkJXTaK0E50FDHgYa6pGNICqj8JCi7Vi5m18rFSceQFFD5iUiQVH4iEqTu/NLyh82s0czeypj2f81sg5mtNbMnzWxQxrxZZlZvZhvN7Jp8BRcRiaM7I79HgGvbTHsRGOPuFwL/CcwCMLPzgWnAZ6N1fmZmfXKWVkQkR7q81MXdV5pZVZtpv814+QowNXo+GVji7geB982sHhgHvJyTtCIxlU6YkXQESYlcXOf3v4B/iZ5X0FqGxzRE00RSoe/Qs5OOICkR6wMPM/secBj41bFJ7SzmHaw7w8xqzKymqakpTgyRbmvZVEvLptqkY0gKZD3yM7PpwPXABHc/VnANwPCMxSqBre2t7+4LgYUA1dXV7RakSK7t/sMSQN/oLFmWn5ldC9wL/Lm7N2fMWgY8amY/Aj4DjARejZ1SpMCqZj5bkP1smn9dQfYjn9Rl+ZnZY8CVwBAzawBm0/rp7inAi2YG8Iq73+Xu681sKfA2rafDd7v7kXyFFxHJVnc+7b21ncmLOll+LjA3TigRkXzTHR4iEiR9pZUE5fRr7kk6gqSEyk+CUnR6ZdIRJCV02itBaa5fTXP96qRjSApo5CdB2fPqkwCUjLgk4SSSNI38RCRIKj8RCZLKT0SCpPITkSDpAw8JypDrv510BEkJlZ8E5eSBZUlHkJTQaa8EZX/dSvbXrUw6hqSARn4SlL1vPAdAv9FXJJxEkqaRn4gESeUnIkFS+YlIkFR+IhIkfeAhQSmbMivpCJISKj8JSp+S05KOICmh014Jyr51L7Fv3UtJx5AUUPlJUFR+cozKT0SCpPITkSCp/EQkSF2Wn5k9bGaNZvZWxrRSM3vRzN6JHgdnzJtlZvVmttHMrslXcBGROLoz8nsEuLbNtJnAcncfCSyPXmNm5wPTgM9G6/zMzPrkLK1ITOU3z6H85jlJx5AU6LL83H0l8FGbyZOBxdHzxcCUjOlL3P2gu78P1APjcpRVJLaTik7lpKJTk44hKZDte35D3X0bQPRYHk2vAD7IWK4hmiaSCnvXPMveNc8mHUNSINcfeFg707zdBc1mmFmNmdU0NTXlOIZI+/ZvWMX+DauSjiEpkG35bTezYQDRY2M0vQEYnrFcJbC1vQ24+0J3r3b36rIyfbW4iBRWtuW3DJgePZ8OPJ0xfZqZnWJmZwEjgVfjRRQRyb0uv9jAzB4DrgSGmFkDMBuYDyw1szuBzcDNAO6+3syWAm8Dh4G73f1InrKLiGSty/Jz91s7mDWhg+XnAnPjhBIRyTd9pZUE5YyvzE86gqSEbm8TkSCp/CQou1c/we7VTyQdQ1JAp72SM1Uz03/xcMu7rRcfnHbJlxJOIknTyE9EgqTyE5EgqfxEJEh6z0+CYiefknQESQmVnwRl6C1/m3QESQmd9opIkFR+EpRd//EYu/7jsaRjSAqo/CQoB/74Jgf++GbSMSQFVH4iEiSVn4gESeUnIkHSpS4SlD7FA5OOICmh8pOglN14X9IRJCV02isiQVL5SVD+9PtH+NPvH0k6hqSATnslKAe3bEg6gqSERn4iEiSVn4gESeUnIkHSe34SlJMHDEk6gqRErPIzs78G/jfgwDrgDqAE+BegCtgE3OLuf4qVUiRHhvzld5KOICmR9WmvmVUA/weodvcxQB9gGjATWO7uI4Hl0WsRkVSJ+57fyUCxmZ1M64hvKzAZWBzNXwxMibkPkZz56KWFfPTSwqRjSApkfdrr7lvM7IfAZqAF+K27/9bMhrr7tmiZbWZWnqOsIrF93Phe0hEkJeKc9g6mdZR3FvAZoJ+Z3daD9WeYWY2Z1TQ1NWUbQ0QkK3FOe68G3nf3Jnc/BDwBXApsN7NhANFjY3sru/tCd6929+qysrIYMUREei5O+W0GvmBmJWZmwASgDlgGTI+WmQ48HS+iiEjuxXnPb7WZ/RpYAxwG3gAWAv2BpWZ2J60FeXMugorkQlFpRdIRJCViXefn7rOB2W0mH6R1FCiSOqdf+82kI0hK6PY2EQmSyk+CsvM3/8DO3/xD0jEkBXRvrwTl0Edbko4gKaGRn4gESeUnIkFS+YlIkPSenwSlb/nZSUeQlFD5SVBKr56RdARJCZ32ikiQVH4SlB3/9kN2/NsPk44hKaDTXgnK4b07ko4gKaGRn4gESeUnIkFS+YlIkPSenwTllIrzko4gKaHyk6AM/vPbk44gKaHTXhEJkspPgtL05DyanpyXdAxJAZ32SlCOtOxJOoKkhEZ+IhIklZ+IBEnlJyJB0nt+EpRT/+yipCNISqj8JCiDLrs16QiSErFOe81skJn92sw2mFmdmY03s1Ize9HM3okeB+cqrIhIrsR9z+8nwG/c/TzgIqAOmAksd/eRwPLotUgqbF86m+1LZycdQ1Ig6/Izs4HAFcAiAHf/2N13AZOBxdFii4EpcUOK5IofPogfPph0DEmBOCO/s4Em4P+Z2Rtm9gsz6wcMdfdtANFjeXsrm9kMM6sxs5qmpqYYMUREei5O+Z0MXAz8o7uPBfbTg1Ncd1/o7tXuXl1WVhYjhohIz8X5tLcBaHD31dHrX9NaftvNbJi7bzOzYUBj3JAin1ZVM5/N+z42zb8u7/vojbIe+bn7h8AHZjYqmjQBeBtYBkyPpk0Hno6VUCSHis8ZR/E545KOISkQ9zq/bwK/MrO+wHvAHbQW6lIzuxPYDNwccx8iOXPaJV9KOoKkRKzyc/daoLqdWRPibFdEJN90b68E5cNHZ/Lho7r0VFR+IhIolZ+IBEnlJyJBUvmJSJD0lVYSlH7nXZ50BEkJlZ8EZcDFuttBWum0V4Jy9NABjh46kHQMSQGN/CQojf86B4AzvjI/2SCSOI38RCRIKj8RCZLKT0SCpPITkSDpAw8JSv8Lrk46gqSEyk+CovKTY3TaK0E50rybI827k44hKaDyk6A0PXU/TU/dn3QMSQGVn4gESeUnIkHSBx6BKMSvSBTpTTTyE5EgaeQnQRkwdlLSESQlVH4SlH6jr0g6gqRE7NNeM+tjZm+Y2TPR61Ize9HM3okeB8ePKZIbh/c0cXhPU9IxJAVy8Z7ft4C6jNczgeXuPhJYHr0WSYUdzzzIjmceTDqGpECs8jOzSuA64BcZkycDi6Pni4EpcfYhIpIPcUd+fw98FziaMW2ou28DiB7LY+5DRCTnsi4/M7seaHT317Ncf4aZ1ZhZTVOT3oMRkcKKM/K7DLjBzDYBS4CrzOyfge1mNgwgemxsb2V3X+ju1e5eXVZWFiOGiEjPZV1+7j7L3SvdvQqYBvzO3W8DlgHTo8WmA0/HTimSIwPH3cjAcTcmHUNSIB/X+c0HlprZncBm4OY87EMkKyUjLkk6gqRETsrP3VcAK6LnO4EJudiuSK4d2tkAQNHplQknkaTp3l4Jys4XFrDzhQVJx5AUUPmJSJBUfiISJJWfiARJ5SciQdJXWklQTrt0WtIRJCVUfhKU4qrPJR1BUkKnvRKUj7e/x8fb30s6hqSAyk+C8tHyhXy0fGHSMSQFVH4iEiSVn4gESeUnIkFS+YlIkHSpiwRl0BXTu15IgqDyk6CcWjk66QiSEjrtlaAcaKjjQENd1wvKp57KT4Kya+Vidq1c3PWC8qmn8hORIKn8RCRIKj8RCZLKT0SCpEtdJCilE2YkHUFSQuUnQek79OykI0hK6LRXgtKyqZaWTbVJx5AUyLr8zGy4mf27mdWZ2Xoz+1Y0vdTMXjSzd6LHwbmLKxLP7j8sYfcfliQdQ1IgzsjvMPBtdx8NfAG428zOB2YCy919JLA8ei0ikipZl5+7b3P3NdHzvUAdUAFMBo5dQr8YmBI3pIhIruXkPT8zqwLGAquBoe6+DVoLEijPxT5ERHIpdvmZWX/gceCv3H1PD9abYWY1ZlbT1NQUN4aISI/EutTFzIpoLb5fufsT0eTtZjbM3beZ2TCgsb113X0hsBCgurra4+QQ6a7Tr7kn6QiSEnE+7TVgEVDn7j/KmLUMOPaNkdOBp7OPJ5JbRadXUnR6ZdIxJAXijPwuA74GrDOzYxdO3QfMB5aa2Z3AZuDmeBFFcqe5fjUAJSMuSTiJJC3r8nP3/w9YB7MnZLtdkXza8+qTgMpPdIeHiARK5SciQVL5iUiQVH4iEiR9pZUEZcj13046gqSEyk+CcvLAsqQjSErotFeCsr9uJfvrViYdQ1JAIz8Jyt43ngOg3+grEk4iSdPIT0SCpPITkSCp/EQkSCo/EQmSPvCQoJRNmZV0BEkJlZ8EpU/JaUlHKLiqmc/mfR+b5l+X933kmk57JSj71r3EvnUvJR1DUkAjvxQoxL/M0upY8fW/4OqEk0jSNPITkSCp/EQkSCo/EQmSyk9EgqQPPCQo5TfPSTrCp1KhPrTL5SU1Kj8JyklFpyYdQVJCp70SlL1rnmXvGl1aJCo/Ccz+DavYv2FV0jEkBfJWfmZ2rZltNLN6M5uZr/2IiGQjL+/5mVkf4KfAXwANwGtmtszd387VPnS/oojEka+R3zig3t3fc/ePgSXA5DztS0Skx/JVfhXABxmvG6JpIiKpkK9LXaydaX7CAmYzgBnRy31mtrGH+xgC7MgiW7fZA8ln6KY05EhDBuhmjj8+cH3iGQogDTlymqEb/0+29WcdzchX+TUAwzNeVwJbMxdw94XAwmx3YGY17l6d7fq5kIYMacmRhgxpyZGGDGnJkYYMHcnXae9rwEgzO8vM+gLTgGV52peISI/lZeTn7ofN7B7gBaAP8LC7r8/HvkREspG329vc/TnguXxtnxinzDmUhgyQjhxpyADpyJGGDJCOHGnI0C5z966XEhH5lNHtbSISpMTLz8xONbNXzexNM1tvZn/bzjLnmdnLZnbQzL6TMX24mf27mdVF634rY16pmb1oZu9Ej4MTyjHHzLaYWW30MylPGTpct8DHorMcBTkWGfP7mNkbZvZMNscijxm6fRxykcPMNpnZumhfNUkci04y9OhY5JS7J/pD6zWB/aPnRcBq4AttlikHPg/MBb6TMX0YcHH0fADwn8D50eu/A2ZGz2cCDySUY07msnnM0OG6BT4WneUoyLHImP83wKPAMxnTun0s8pih28chFzmATcCQdrZbsGPRSYYeHYtc/iQ+8vNW+6KXRdGPt1mm0d1fAw61mb7N3ddEz/cCdfz3nSSTgcXR88XAlIRydFvMDJ2tW8hj0eW63REnA4CZVQLXAb9oM6vbxyKPGXokbo5OFOxYpFHi5QfHTw1qgUbgRXdfncU2qoCxtP6LBDDU3bdBaznR+q9SEjkA7jGztWb2cDdOObPO0Mm6BT0WXaxbkGMB/D3wXeBom+k9OhZ5ygA9OA45yOHAb83sdWu9s+qYQh6LjjJAD49FrqSi/Nz9iLt/jtY7QcaZ2ZierG9m/YHHgb9y9z0py/GPwDnA54BtwIP5yhA3fwFyFORYmNn1QKO7v97dzAXO0KPjECdH5DJ3vxiYCNxtZlf0YN18Z+jxsciVVJTfMe6+C1gBXNvddcysiNbC+ZW7P5Exa7uZDYuWGUbrv1YFz+Hu26O/NEeBh2j9xpu8ZOhk3YIei47WLeCxuAy4wcw20fqNQleZ2T9H87I6FrnMkO1xyDIH7r41emwEnszYX6GORYcZ4hyLuBIvPzMrM7NB0fNi4GpgQzfXNWARUOfuP2ozexkwPXo+HXg6iRzH/nJFbgTeylOGztYt5LHocN1CHQt3n+Xule5eReutlb9z99ui2d0+FvnK0JPjEDeHmfUzswHHngP/M2N/BTkWnWXo6bHIKU/gU5bMH+BC4A1gbfQf/oNo+l3AXdHzM2j9soQ9wK7o+UDgi7S+l7AWqI1+JkXrnA4sB96JHksTyvFLYF00bxkwLE8Z2l03gWPRWY6CHIs227mSEz9p7faxyGOGbh+HHPx5nA28Gf2sB75X6GPRRYYeHYtc/ugODxEJUuKnvSIiSVD5iUiQVH4iEiSVn4gESeUnIkFS+YlIkFR+IhIklZ+IBOm/AETvuXRs2nTiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 5))\n",
    "ax.hist(pi, bins=10, align='mid')\n",
    "ax.axvline(np.mean(pi), linestyle='dashed', color='k', label=r'$\\pi$={:.3f}'.format(np.mean(pi)))\n",
    "ax.legend(loc='best', fontsize=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step3 - Further improvements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above example we have replaced the inner *for* loop with its corresponding \"NumPy version\".\n",
    "The same approach could be applied also for the outer *for* loop, taking care of the array dimensions of our object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pi with 100000 steps for 500 runs is 3.141567 in 1.066s\n"
     ]
    }
   ],
   "source": [
    "points = np.random.uniform(low=-1, high=1, size=(2, N_run, N))\n",
    "areas = points**2\n",
    "Nhints = np.sum(areas, axis=0) < 1\n",
    "pi = np.mean(Nhints, axis=1) * 4\n",
    "\n",
    "print ('pi with {:d} steps for {:d} runs is {:.6f} in {:.3f}s'.format(\n",
    "    N, N_run, np.mean(pi), toc-tic))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1b1a43c4820>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUcAAAEvCAYAAADW7gNcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAYkUlEQVR4nO3dfXBV9Z3H8fdXCJoQIkQSygZslkJRF1ScVItW1xF3W8AKXdHa1g463WGsusv2YdboTiv7hzzsaLudpd0Olq4MW5dlq6IrUlfRLFgtGjGCbmBFlypPeYACgQQM8N0/cmQSOCTkPnB+8XxeM5l77zn33PPhiB9+v3vOvTF3R0REujor6QAiIiFSOYqIxFA5iojEUDmKiMRQOYqIxFA5iojE6J90AIChQ4d6ZWVl0jEkRTZv3gzA2LFjE04i+fTGG280u3tZJtsGUY6VlZXU1tYmHUNS5L777gNg3rx5CSeRfDKz32e6bRDlKHKmqRSlJ3rPUUQkhspRUummm27ipptuSjqGBEzTakml3bt3Jx1BAqeRo4hIjD4xcty/fz+NjY20t7cnHUVOU0FBAeXl5ZSUlCQdRSQjwZfj/v37aWhooKKigsLCQsws6UjSA3enra2N7du3A6ggpU8KvhwbGxupqKigqKgo6ShymsyMoqIiKioq2LFjR5DlOGnSpKQjSOCCL8f29nYKCwuTjiEZKCwsDPatkB/84AdJR5DA9YkTMppK90367yZ9WZ8oR5Fcmzx5MpMnT046hgQs+Gm1yIkqq1dm/Rq73tnW7WttnT81631I39bjyNHMfmlmjWb2dqdlpWb2vJm9G90O6bTuPjPbYmabzeyL+QqeRj/96U+5+OKLKSkpoaSkhIkTJ7JyZfdFsWbNGm688UYqKiowMx599NFunz937lzMjHvuuafL8nnz5vG5z32OkpISysrK+PKXv8zbb799ilcR6ftOZ1r9KPClE5ZVA6vdfQywOnqMmV0E3Ar8SbTNz8ysX87SptyIESNYsGAB69evp7a2luuuu47p06ezYcOGU25z4MABxo0bx09+8pMeT2z97ne/45FHHuHiiy8+aV1NTQ133XUXr7zyCi+++CL9+/fn+uuvZ8+ePVn/uURC1GM5uvsa4MT/A6YBS6L7S4DpnZYvc/fD7v5/wBbg8hxl7XMuuugizCz2Z86cOb1+vWnTpjF58mRGjx7NZz/7WR588EEGDRrEq6++esptpkyZwty5c5kxYwZnnXXq/9z79u3jG9/4BosXL2bIkCEnrX/uuee44447GDduHOPHj2fp0qU0NTXx29/+ttd/DpG+INMTMsPcfSdAdFseLa8APuz0vG3RslRasWIFAM8++yw7d+5kx44dFBUVsXjxYu69917mzp1LcXFxtz9r166Nfe2jR4+ybNkyDhw4wJVXXpl11lmzZjFjxgyuu+6603p+S0sLx44diy3SvqDwM5dT+JnU/rstpyHXJ2Tirt3w2CeazQJmAZx//vm93tG111570rJbbrmFu+66i9bWVqZMmXLS+ttvv53bb7+d5uZmZsyYcdL6b3/723z1q1/lww8/5Jvf/GaXdTU1Nb3O2NDQgJlx9dVXU1xczJYtW2htbeULX/gChYWF3Hnnndxyyy3dvkZFRdd/WzZu3MjEiRM5dOgQxcXFPPnkk4wfP77X2Tp75JFH2LJlC0uXLj3tbWbPns2ll17KxIkTs9p3Us694i+6XZ+Lkz6nQyd+wpVpOTaY2XB332lmw4HGaPk2YGSn540AdsS9gLsvAhYBVFVVxRZoX/fWW28xatQoiouLAairq6OoqIjRo0cDUFpaSmlpaa9ec+zYsdTV1bF3714ef/xxZs6cSU1NDePGjcso4+bNm7n//vtZu3YtAwYMOK1tvvvd7/Lyyy/z8ssv06+f3lKWT6ZMy/FpYCYwP7p9qtPyx8zsR8AfAWOA17INGae7kVxRUVG364cOHdrt+pEjR2Y0UjzRhg0buOSSS44/rqurY9y4ccff+5s7dy5z587t9jVWrVrF1VdfffzxgAEDjpdrVVUVr7/+Oj/+8Y9ZvHhxRhlfffVVmpubu5Tr0aNHWbNmDT//+c85ePAgZ5999vF13/nOd1i2bBkvvfQSo0aNymifIdj1WDUAn/r6/ISTSKh6LEcz+zfgWmComW0DHqCjFJeb2beAD4CbAdz9HTNbDvwPcAS4292P5il78DZs2NDlQuO6urouZZnJtPpEx44d4/DhwxlnnD59OlVVVV2W3XHHHYwZM4b777+/y2hy9uzZLFu2jJqaGi644IKM9ynSF/RYju7+tVOsiv3kvrs/CDyYTahPgmPHjvH2229z7733Hl/23nvvdXmPrrfT6urqaqZOncrIkSNpaWnhscceo6ampsu1jgsXLmThwoVs2rQJ6LiUZ8uWLcczffDBB9TV1VFaWsr555/P4MGDGTx4cJf9DBw4kNLS0i6jybvvvpulS5eyYsUKhgwZwq5duwCOnzgS+aTRxwfz5L333uPgwYNdrhkcP348CxYsYNWqVRm95q5du7jtttsYO3YskyZN4vXXX2fVqlVdRqfNzc3Hf+0oQG1tLRMmTGDChAm0tbXxwAMPMGHCBH74wx/2at8/+9nPaGlpYdKkSQwfPvz4z0MPPZTRn0UkdPr4YJ6MGTMG967nmZYvX57Va/b06RaAOXPmdLmG8tprrz0pR0/i3m/t7WuI9HUqR0mlgRdc3fOTJNVUjpJKgy7T9YXSPb3nKKl0rP0Qx9oPJR1DAqaRo6RS43/MAXSdo5xanxg56mRA36T/btKXBV+OBQUFtLW1JR1DMtDW1kZBQUHSMUQyEnw5lpeXs337dlpbWzUS6SPcndbWVrZv3055eXnPG4gEKPj3HD/+tZ47duwI9jfZyckKCgoYNmxYkL+WVeR0BF+OwPFfCyBhO/lrvjbHPi8ExeOvTzqCBK5PlKNIrqkcpSfBv+cokg9HW/dxtHVf0jEkYCpHSaWmFfNoWjEv6RgSMJWjiEgMlaOISAyVo4hIDJWjiEgMXcojqTRowsm/ulekM5WjpNLAC69JOoIETtNqSaUj+5s4sr8p6RgSMJWjpFLzMw/T/MzDSceQgKkcRURiqBxFRGKoHEVEYqgcRURi6FIeSaWSy7+SdAQJnMpRUqlo9BVJR5DAaVotqdS+exvtu7clHUMCpnKUVNr93EJ2P7cw6RgSMJWjiEgMlaOISAyVo4hIDJWjiEgMXcojqXTulbcmHUECp3KUVCqsvDTpCBI4TasllT5qeJ+PGt5POoYETOUoqbRn9SL2rF6UdAwJmMpRRCSGylFEJIbKUUQkhspRRCRGVpfymNl3gL8EHNgI3AEUAf8OVAJbgVvc/Q9ZpRTJscHXzEw6ggQu45GjmVUAfw1Uufs4oB9wK1ANrHb3McDq6LFIUM4ZcSHnjLgw6RgSsGyn1f2BQjPrT8eIcQcwDVgSrV8CTM9yHyI5d2hbPYe21ScdQwKWcTm6+3bgIeADYCewz93/Cxjm7juj5+wEynMRVCSX9q5Zwt41S3p+oqRWNtPqIXSMEv8Y+CNgoJnd1ovtZ5lZrZnVNjU1ZRpDRCQvsplWXw/8n7s3uXs78ARwJdBgZsMBotvGuI3dfZG7V7l7VVlZWRYxRERyL5ty/AD4vJkVmZkBk4B64Gng41OBM4GnsosoInLmZXwpj7uvM7NfA+uBI8CbwCKgGFhuZt+io0BvzkVQEZEzKavrHN39AeCBExYfpmMUKRKs0kmzko4ggdP3OUoqDRg2KukIEjh9fFBSqW1rHW1b65KOIQHTyFFSad8ry4DkvxG8snpl3vexdf7UvO/jk0gjRxGRGCpHEZEYKkcRkRgqRxGRGDohI6l03hfvSTqCBE7lKKlUcN6IpCNI4DStllRq3bKO1i3rko4hAdPIUVJp/2tPAlA0+oqEk0ioNHIUEYmhchQRiaFyFBGJoXIUEYmhEzKSSkNv+F7SESRwKkdJpf4l+r1F0j1NqyWVDtav4WD9mqRjSMA0cpRUannzWQAGXnhNwkkkVBo5iojEUDmKiMRQOYqIxFA5iojE0AkZSaWy6fclHUECp3KUVOpXdG7SESRwmlZLKh3Y+AIHNr6QdAwJmMpRUknlKD1ROYqIxFA5iojEUDmKiMRQOYqIxNClPJJK5TfPSTqCBE7lKKl0VsE5SUeQwGlaLanUsn4lLetXJh1DAqZylFQ6uGktBzetTTqGBEzlKCISQ+UoIhJD5SgiEkPlKCISQ5fySCp96uvzk44ggdPIUUQkhspRUmnfuifYt+6JpGNIwLIqRzMbbGa/NrNNZlZvZhPNrNTMnjezd6PbIbkKK5Irbe+9Rtt7ryUdQwKW7cjxJ8Bv3P0C4BKgHqgGVrv7GGB19FhEpE/JuBzNrAS4BlgM4O4fufteYBqwJHraEmB6tiFFRM60bEaOo4Am4F/M7E0z+4WZDQSGuftOgOi2PAc5RUTOqGzKsT9wGfDP7j4BOEgvptBmNsvMas2stqmpKYsYIr1n/c/G+p+ddAwJWDbXOW4Dtrn7uujxr+koxwYzG+7uO81sONAYt7G7LwIWAVRVVXkWOUR6bdgtf590BAlcxiNHd98FfGhmY6NFk4D/AZ4GZkbLZgJPZZVQRCQB2X5C5q+AX5nZAOB94A46Cne5mX0L+AC4Oct9iOTc3t/+GwCDr/pawkkkVFmVo7vXAVUxqyZl87oi+Xbo92913FE5yinoEzIiIjFUjiIiMVSOIiIx9JVlkkr9CkuSjiCBUzlKKpV95f6kI0jgNK0WEYmhcpRU+sN/P8of/vvRpGNIwDStllQ6vH1T0hEkcBo5iojEUDmKiMRQOYqIxNB7jpJK/QcNTTqCBE7lKKk09MvfTzqCBE7TahGRGCpHSaU9LyxizwuLko4hAdO0WlLpo8b3k44ggdPIUUQkhspRRCSGylFEJIbec5RUKiitSDqCBE7lKKl03pf+KukIEjhNq0VEYqgcJZV2/+af2P2bf0o6hgRM02pJpfY925OOIIHTyFFEJIbKUUQkhspRRCSG3nOUVBpQPirpCBI4lWNKVFavTDpCUEqvn5V0BAmcptUiIjFUjpJKzf/5EM3/+VDSMSRgmlZLKh1paU46ggROI0cRkRgqRxGRGCpHEZEYes9RUunsiguSjiCBUzlKKg3509uTjiCB07RaRCSGylFSqenJuTQ9OTfpGBIwTasllY627U86ggROI0cRkRhZl6OZ9TOzN83smehxqZk9b2bvRrdDso8pInJm5WLkOBuo7/S4Gljt7mOA1dFjEZE+JatyNLMRwFTgF50WTwOWRPeXANOz2YdIPpzz6Us459OXJB1DApbtCZl/BP4WGNRp2TB33wng7jvNrDzLfYjk3OCrvpZ0BAlcxiNHM7sBaHT3NzLcfpaZ1ZpZbVNTU6YxRETyIptp9VXAjWa2FVgGXGdm/wo0mNlwgOi2MW5jd1/k7lXuXlVWVpZFDJHea1j+AA3LH0g6hgQs43J09/vcfYS7VwK3Ai+6+23A08DM6GkzgaeyTimSY37kMH7kcNIxJGD5uM5xPvBnZvYu8GfRYxGRPiUnn5Bx9xqgJrq/G5iUi9cVEUmKPiEjIhJDn62WVCr8zOVJR5DAqRwllc694i+SjiCB07RaRCSGylFSaddj1ex6TB/7l1NTOYqIxFA5iojEUDmKiMRQOYqIxNClPJJKAy+4OukIEjiVo6TSoMumJh1BAqdptaTSsfZDHGs/lHQMCZhGjpJKjf8xB4BPfV1fGiXxNHIUEYmhchQRiaFyFBGJoXIUEYmhEzKSSsXjr086ggRO5SiplKZyrKxemfd9bJ3/ybtuVNNqSaWjrfs42rov6RgSMJWjpFLTink0rZiXdAwJmMpRRCSGylFEJIbKUUQkhspRRCSGLuWRVBo0YUrSESRwKkdJpYEXXpN0BAmcptWSSkf2N3Fkf1PSMSRgKkdJpeZnHqb5mYeTjiEBUzmKiMRQOYqIxFA5iojEUDmKiMTQpTySSiWXfyXpCBI4laOkUtHoK5KOIIHTtFpSqX33Ntp3b0s6hgRM5SiptPu5hex+bmHSMSRgKkcRkRh6zzEAZ+J3fIhI72jkKCISQ+UoIhJD02pJpXOvvDXpCBK4jEeOZjbSzF4ys3oze8fMZkfLS83seTN7N7odkru4IrlRWHkphZWXJh1DApbNtPoI8D13vxD4PHC3mV0EVAOr3X0MsDp6LBKUjxre56OG95OOIQHLuBzdfae7r4/utwD1QAUwDVgSPW0JMD3bkCK5tmf1IvasXpR0DAlYTk7ImFklMAFYBwxz953QUaBAeS72ISJyJmVdjmZWDDwO/I277+/FdrPMrNbMapua9HX1IhKWrMrRzAroKMZfufsT0eIGMxserR8ONMZt6+6L3L3K3avKysqyiSEiknPZnK02YDFQ7+4/6rTqaWBmdH8m8FTm8UREkpHNdY5XAd8ENppZXbTsfmA+sNzMvgV8ANycXUSR3Bt8zcyenySplnE5uvvLgJ1i9aRMX1fkTDhnxIVJR5DA6eODkkqHttVzaFt90jEkYCpHSaW9a5awd82Snp8oqaVyFBGJoXIUEYmhchQRiaFyFBGJoe9zlFQqnTQr6QgSOJWjpNKAYaOSjiCBUzlKKrVt7fhQl77wNjfO1C+J2zp/6hnZD6gcJaX2vbIMUDnKqemEjIhIDJWjiEgMlaOISAyVo4hIDJ2QkVQ674v3JB1BAqdylFQqOG9E0hEkcJpWSyq1bllH65Z1SceQgGnkKKm0/7UnASgafUXCSSRUGjmKiMRQOYqIxFA5iojEUDmKiMTQCRlJpaE3fC/pCBI4lWM3ztTXMMmZ17+kLOkIEjhNqyWVDtav4WD9mqRjSMA0cpRUannzWQAGXnhNwkkkVBo5iojEUDmKiMRQOYqIxFA5iojE0AkZSaWy6fclHUECp3KUVOpXdG7SESRwmlZLKh3Y+AIHNr6QdAwJmMpRUknlKD1ROYqIxFA5iojEUDmKiMRQOYqIxNClPJJK5TfPSTqCBE7lKKl0VsE5SUeQwGlaLanUsn4lLev1ZcZyan125Khv6ZZsHNy0FoBBl01NOImEKm8jRzP7kpltNrMtZladr/2IiORDXsrRzPoBPwUmAxcBXzOzi/KxLxGRfMjXyPFyYIu7v+/uHwHLgGl52peISM7lqxwrgA87Pd4WLRMR6RPydULGYpZ5lyeYzQJmRQ8PmNnmLPY3FGjOYvszSVlzL+Ocv19wQ46j9KivHFMIMKstiF3cXc5PZ7qvfJXjNmBkp8cjgB2dn+Dui4BFudiZmdW6e1UuXivflDX3+kpOUNZ8yFfOfE2rXwfGmNkfm9kA4Fbg6TztS0Qk5/IycnT3I2Z2D/Ac0A/4pbu/k499iYjkQ94uAnf3Z4Fn8/X6J8jJ9PwMUdbc6ys5QVnzIS85zd17fpaISMros9UiIjGCKkczO8fMXjOzt8zsHTP7+5jnXGBmr5rZYTP7fqflI83sJTOrj7ad3WndHDPbbmZ10c+UhLOeclszKzWz583s3eh2SMBZgzqundb3M7M3zeyZTstyelzzmDO4Y2pmW81sY5SnttPyoP6u9pC198fV3YP5oeP6yOLofgGwDvj8Cc8pBz4HPAh8v9Py4cBl0f1BwP8CF0WP53R+bgBZT7kt8A9AdXS/GlgQcNagjmun9d8FHgOe6bQsp8c1jzmDO6bAVmBozOsG9Xe1h6y9Pq5BjRy9w4HoYUH04yc8p9HdXwfaT1i+093XR/dbgHry+KmcLLN2t+00YEl0fwkwPeCsOZdNVgAzGwFMBX5xwqqcHtc85sy5bLN2I6i/q7kWVDnC8alGHdAIPO/u6zJ4jUpgAh3/6nzsHjPbYGa/zMXwP9us3Ww7zN13Qkfh0/GvZKhZIbDjCvwj8LfAsROW5/y45iknhHdMHfgvM3vDOj7Z9rHg/q52kxV6eVyDK0d3P+rul9LxqZrLzWxcb7Y3s2LgceBv3H1/tPifgc8AlwI7gYeTzprtn7O38pQ1qONqZjcAje7+Ri5y9CRPOYM6ppGr3P0yOr5l624zuyYXmU4lT1l7fVyDK8ePufteoAb40uluY2YFdBTjr9z9iU6v1RAd8GPAI3R8a1CiWbvZtsHMhgNEt425SXnK/WW8bYDH9SrgRjPbSsc3QV1nZv8arcvbcc1lzgCPKe6+I7ptBJ7slCm4v6unyprJcQ2qHM2szMwGR/cLgeuBTae5rQGLgXp3/9EJ64Z3evgV4O2Es3a37dPAzOj+TOCpULOGdlzd/T53H+HulXR8ZPVFd78tWp3T45qvnKEdUzMbaGaDPr4P/HmnTKH9XT1l1oyOq+fwrFi2P8DFwJvAhij8D6PldwJ3Rvc/RccXW+wH9kb3S4Av0PF+wwagLvqZEm2zFNgYrXsaGJ5w1thto23OA1YD70a3pQFnDeq4nvA619L1LHBOj2secwZ1TIFRwFvRzzvA3wX8d7W7rL0+rvqEjIhIjKCm1SIioVA5iojEUDmKiMRQOYqIxFA5iojEUDmKiMRQOYqIxFA5iojE+H+jtbj/CT7/cwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 5))\n",
    "ax.hist(pi, bins=10, align='mid')\n",
    "ax.axvline(np.mean(pi), linestyle='dashed', color='k', label=r'$\\pi$={:.3f}'.format(np.mean(pi)))\n",
    "ax.legend(loc='best', fontsize=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extra - C++ version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " For sake of completeness we implement the same functions also in C++.\n",
    " \n",
    " >⚠️**Note:** C++ is always **faster** than Python but more complex to write!\n",
    " > Furthermore, it is not straightforward the visualization of the results: if you are interested into a graphic library (please do not spend a lot of time to re-implement the same things!) you can use [GnuPlot](http://www.gnuplot.info/), [OpenCV](https://docs.opencv.org/3.4/d0/d1e/classcv_1_1plot_1_1Plot2d.html), or the easier [cv-plot](https://github.com/Profactor/cv-plot) libraries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```c++\n",
    "// g++ metropolis.cpp -std=c++14 -O3 -o metropolis -fopenmp\n",
    "\n",
    "#include <random>\n",
    "#include <iostream>\n",
    "#include <iomanip>\n",
    "#include <numeric>\n",
    "#include <omp.h>\n",
    "\n",
    "int main()\n",
    "{\n",
    "  std :: random_device rd;\n",
    "  std :: mt19937 eng{rd()};\n",
    "  \n",
    "  int64_t N = 1e5, N_run = 1e4;\n",
    "  \n",
    "  double pi = 0.;\n",
    "  double start_time, run_time;\n",
    "    \n",
    "  std :: uniform_real_distribution < double > uniform_dist {-1., 1.};\n",
    "\n",
    "  omp_set_num_threads(omp_get_max_threads());\n",
    "  \n",
    "  start_time = omp_get_wtime();\n",
    "  #pragma omp parallel\n",
    "  {\n",
    "    #pragma omp for reduction (+ : pi)\n",
    "    for (int64_t I = 0; I < N_run; ++I)\n",
    "    {\n",
    "      int64_t Nhints = 0;\n",
    "            \n",
    "      #pragma omp parallel\n",
    "      {\n",
    "        #pragma omp for reduction (+ : Nhints)\n",
    "        for (int64_t i = 0; i < N; ++i)\n",
    "        {\n",
    "          const double x = uniform_dist(eng);\n",
    "          const double y = uniform_dist(eng);\n",
    "          const double res = x*x + y*y;\n",
    "          Nhints += res < 1 ? 1 : 0;\n",
    "        }\n",
    "      }\n",
    "      pi += 4.0 * (static_cast < double >(Nhints) / N);\n",
    "    }\n",
    "  }\n",
    "  \n",
    "  pi /= N_run;\n",
    "  run_time = omp_get_wtime() - start_time;\n",
    "  std :: cout << \"pi with \" << N << \" steps for \" << N_run << \" runs is \" << std :: setprecision(16) << pi << \" in \" << run_time << \" sec\" << std :: endl;\n",
    "    \n",
    "  return 0;\n",
    "}   \n",
    "```"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}