{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calcule numérique avec la méthode Monte-Carlo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Approximation numérique d'une surface avec la méthode Monte-Carlo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Conceptuellement, pour calculer la surface $\\mathcal A_1$ d'un objet $O_1$ avec la methode Monte-Carlo, il suffit:\n",
    "\n",
    "1. de placer cet objet $O_1$ entièrement dans une figure géométrique $O_2$ dont on connait la surface $\\mathcal A_2$ (par exemple un carré ou un rectangle)\n",
    "2. tirer aléatoirement un grand nombre de points dans cette figure  $O_2$ (tirage uniforme)\n",
    "3. compter le nombre de points $N_1$ tombés dans l'objet  $O_1$ dont on veut calculer la surface\n",
    "4. calculer le rapport $\\frac{N_1}{N}$ où $N$ est le nombre total de points tiré aléatoirement (en multipliant par 100 ce rapport on obtient le pourcentage de points tombés dans l'objet $O_1$)\n",
    "5. appliquer ce rapport à la surface $\\mathcal A_2$ de figure englobante $O_2$ (le carré, rectangle, ... dont on connait la surface) pour obtenir la surface $\\mathcal A_1$ recherchée: $\\mathcal A_1 \\simeq \\frac{N_1}{N} \\mathcal A_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "t = np.linspace(0., 2. * np.pi, 100)\n",
    "x = np.cos(t) + np.cos(2. * t)\n",
    "y = np.sin(t)\n",
    "\n",
    "N = 100\n",
    "rand = np.array([np.random.uniform(low=-3, high=3, size=N), np.random.uniform(low=-3, high=3, size=N)]).T\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(6, 6))\n",
    "ax.plot(rand[:,0], rand[:,1], '.k')\n",
    "ax.plot(x, y, \"-r\", linewidth=2)\n",
    "ax.plot([-3, -3, 3, 3, -3], [-3, 3, 3, -3, -3], \"-b\", linewidth=2)\n",
    "ax.set_axis_off()\n",
    "ax.set_xlim([-4, 4])\n",
    "ax.set_ylim([-4, 4])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le même principe peut être appliqué pour calculer un volume.\n",
    "\n",
    "Cette methode très simple est parfois très utile pour calculer la surface (ou le volume) de figures géometriques complexes. En revanche, elle suppose l'existance d'une procedure ou d'une fonction permettant de dire si un point est tombé dans l'objet $O_2$ ou non."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application au calcul d'intégrales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculer l'intégrale d'une fonction, revient à calculer la surface entre la courbe décrivant cette fonction et l'axe des abscisses (les surfaces au dessus de l'axe des abscisses sont ajoutées et les surfaces au dessous sont retranchées)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exemple written in Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing()\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The function to integrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return -x**2 + 3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Random points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N = 100000           # The number of random points\n",
    "\n",
    "x_lower_bound = -4.0\n",
    "x_upper_bound = 4.0\n",
    "y_lower_bound = -16.0\n",
    "y_upper_bound = 16.0\n",
    "\n",
    "random_points = np.array([np.random.uniform(low=x_lower_bound, high=x_upper_bound, size=N),\n",
    "                          np.random.uniform(low=y_lower_bound, high=y_upper_bound, size=N)]).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Numerical computation of the integral with Monte-Carlo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Points between f and the abscissa\n",
    "random_points_in_pos = np.array([p for p in random_points if 0 <= p[1] <= f(p[0])])\n",
    "random_points_in_neg = np.array([p for p in random_points if 0 >  p[1] >= f(p[0])])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ratio_pos = float(len(random_points_in_pos)) / float(N)\n",
    "ratio_neg = float(len(random_points_in_neg)) / float(N)\n",
    "print('Percentage of \"positive\" points between f and the abscissa: {:.2f}%'.format(ratio_pos * 100))\n",
    "print('Percentage of \"negative\" points between f and the abscissa: {:.2f}%'.format(ratio_neg * 100))\n",
    "\n",
    "s2 = (x_upper_bound - x_lower_bound) * (y_upper_bound - y_lower_bound)\n",
    "print(\"Box surface:\", s2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "s1 = ratio_pos * s2 - ratio_neg * s2\n",
    "print(\"Function integral (numerical computation using Monte-Carlo):\", s1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The actual integral value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = sp.symbols(\"x\")\n",
    "\n",
    "integ = sp.Integral(f(x), (x, x_lower_bound, x_upper_bound))\n",
    "sp.Eq(integ, integ.doit())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The error ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "actual_s1 = float(integ.doit())\n",
    "\n",
    "error = actual_s1 - s1\n",
    "print(\"Error ratio = {:.6f}%\".format(abs(error / actual_s1) * 100.))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Graphical illustration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_array = np.arange(x_lower_bound, x_upper_bound, 0.01)\n",
    "y_array = f(x_array)\n",
    "\n",
    "plt.axis([x_lower_bound, x_upper_bound, y_lower_bound, y_upper_bound])\n",
    "plt.plot(random_points[:,0], random_points[:,1], ',k')\n",
    "plt.plot(random_points_in_pos[:,0], random_points_in_pos[:,1], ',r')\n",
    "plt.plot(random_points_in_neg[:,0], random_points_in_neg[:,1], ',r')\n",
    "plt.hlines(y=0, xmin=x_lower_bound, xmax=x_upper_bound)\n",
    "plt.plot(x_array, y_array, '-r', linewidth=2)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application au calcul de la constante $\\pi$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exemple written in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Random points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N = 100000           # The number of random points\n",
    "\n",
    "lower_bound = 0.0\n",
    "upper_bound = 1.0\n",
    "\n",
    "random_points = np.random.uniform(low=lower_bound, high=upper_bound, size=(N, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Numerical computation of $\\pi$ with Monte-Carlo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Points in the circle, i.e. $x^2 + y^2 \\leq 1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...\n",
    "\n",
    "Compute $\\pi$ (with $r=1$):\n",
    "\n",
    "square_area = $r^2$\n",
    "circle_area = $\\pi * r^2$\n",
    "\n",
    "quarter_circle_area = $\\frac{\\pi * r^2}{4}$\n",
    "<=> quarter_circle_area = $\\pi$ * square_area / 4\n",
    "<=> $\\pi$ = quarter_circle_area / square_area * 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "random_points_in = np.array([p for p in random_points if p[0]**2 + p[1]**2 <= 1.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ratio = float(len(random_points_in)) / float(N)\n",
    "print('Percentage of points in the circle: {:.2f}%'.format(ratio * 100))\n",
    "\n",
    "s2 = (upper_bound - lower_bound)**2\n",
    "print(\"Box surface:\", s2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "s1 = ratio * s2\n",
    "pi = 4. * s1\n",
    "print(\"Pi (numerical computation using Monte-Carlo):\", pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The error ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "error = math.pi - pi\n",
    "print(\"Error ratio = {:.6f}%\".format(abs(error / math.pi) * 100.))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Graphical illustration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_array = np.arange(lower_bound, upper_bound, 0.01)\n",
    "y_array = np.sqrt(1. - x_array**2)\n",
    "\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "plt.axis('equal')\n",
    "plt.plot(random_points[:,0], random_points[:,1], ',k')\n",
    "plt.plot(random_points_in[:,0], random_points_in[:,1], ',r')\n",
    "plt.hlines(y=0, xmin=lower_bound, xmax=upper_bound)\n",
    "plt.plot(x_array, y_array, '-r', linewidth=2)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}