Commit c254c552 authored by Matthieu Boileau's avatar Matthieu Boileau
Browse files

Split mc2020.ipynb notebook in 3

parent 9eaaa1aa
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Earth movers\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"ein.tags": "worksheet-0",
"slideshow": {
"slide_type": "-"
}
},
"source": [
"First, some python initializations."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"autoscroll": false,
"ein.hycell": false,
"ein.tags": "worksheet-0",
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"from matplotlib import rcParams\n",
"from matplotlib import pyplot as plt\n",
"from ipywidgets import interact\n",
"import ipywidgets as widgets\n",
"rcParams['figure.figsize'] = (8., 6.) # Enlarge figure\n",
"rcParams['animation.html'] = 'html5' # to render animation in notebook"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solve the earth movers problem for 50-position samples."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from earth_movers import EarthMovers\n",
"\n",
"em = EarthMovers(50)\n",
"em.plot_ot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot histogram of distance for 2000-position samples."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"em2000 = EarthMovers(2000)\n",
"print(f\"Wasserstein distance: {em2000.get_wasserstein_distance():f}\")\n",
"em2000.plot_distance_histogram(bins=20)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Todo\n",
"\n",
"- Uniform random blue and red points on a square #\n",
"- Its optimal mathching, with p=1, n=500 #\n",
"- Histogram of matching length in d=1,2,3 #\n",
"- one dimensional matching for p=1.1 and p=0.9, comparison\n",
"- The scaling algorithm for local optimal matching\n",
"\n",
"PoT: <https://pot.readthedocs.io/en/stable/auto_examples/plot_OT_2D_samples.html>"
]
}
],
"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.7.4"
},
"name": "mc2020.ipynb",
"toc": {
"base_numbering": 1,
"nav_menu": {
"height": "120px",
"width": "252px"
},
"number_sections": false,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": null,
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Course of Sébastien Martineau: percolation"
]
},
{
"cell_type": "markdown",
"metadata": {
"ein.tags": "worksheet-0",
"slideshow": {
"slide_type": "-"
}
},
"source": [
"First, some python initializations."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"autoscroll": false,
"ein.hycell": false,
"ein.tags": "worksheet-0",
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"from matplotlib import rcParams\n",
"from matplotlib import pyplot as plt\n",
"from ipywidgets import interact\n",
"import ipywidgets as widgets\n",
"rcParams['figure.figsize'] = (8., 6.) # Enlarge figure\n",
"rcParams['animation.html'] = 'html5' # to render animation in notebook\n",
"\n",
"from percolation import PercolationRect, PercolationHex, percolation_vs_p"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Rectangular lattice"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"percorect = PercolationRect(20, 10)\n",
"interact(percorect.plot, p=widgets.FloatSlider(min=0., max=1., step=0.1, value=0.5));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Hexagonal lattice"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"percohex = PercolationHex(5, 5)\n",
"percohex.compute_clusters(0.5)\n",
"percohex.plot_clusters(add_cluster_id=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"percohex15 = PercolationHex(15, 15)\n",
"interact(percohex15.plot, p=widgets.FloatSlider(min=0., max=1., step=0.1, value=0.5));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Probability of crossing as a function of $p$\n",
"\n",
"Based on 50 simulations on a $30 \\times 30$ lattice."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"percolation_vs_p(30, 30, nsim=50)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Todo\n",
"\n",
"- Standard coupling on honeycomb lattice, \n",
"- Duality of honycomb lattice\n",
"- Standard coupling on square lattice and its dual #\n",
"- Monte Carlo for crossing probabilities on square lattice, threshold phenomena #\n",
"- Monte Carlo for crossing probabilities on honeycomb lattice\n",
"- Percolation on 4 regular trees\n",
"- Percolation on free group with 2 generators #\n",
"- Percolation on the dual of seven triangular tilling\n",
"\n",
"hex cells:\n",
"<https://stackoverflow.com/questions/46525981/how-to-plot-x-y-z-coordinates-in-the-shape-of-a-hexagonal-grid>\n",
"\n",
"graphs: networkx"
]
}
],
"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.7.4"
},
"name": "mc2020.ipynb",
"toc": {
"base_numbering": 1,
"nav_menu": {
"height": "120px",
"width": "252px"
},
"number_sections": false,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": null,
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}
......@@ -2,9 +2,14 @@
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"metadata": {
"ein.tags": "worksheet-0",
"slideshow": {
"slide_type": "-"
}
},
"source": [
"# A set of problems illustrations for Probability Master Class"
"# Course of Marielle Simon: simple random walk"
]
},
{
......@@ -38,14 +43,16 @@
"from ipywidgets import interact\n",
"import ipywidgets as widgets\n",
"rcParams['figure.figsize'] = (8., 6.) # Enlarge figure\n",
"rcParams['animation.html'] = 'html5' # to render animation in notebook"
"rcParams['animation.html'] = 'html5' # to render animation in notebook\n",
"\n",
"import srw # Import simple random walk module from srw.py"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## A simple 2D random walk\n"
"## Random walk on $\\mathbb{Z}^2$\n"
]
},
{
......@@ -73,9 +80,7 @@
},
"outputs": [],
"source": [
"import srw\n",
"\n",
"walk = srw.Walk2D(nstep=100)\n",
"walk = srw.Walk2D(nstep=100) # Create a 100-step random walk\n",
"anim = walk.generate_animation()\n",
"plt.close(anim._fig) # Close the initial figure to display only the animation figure\n",
"anim # Now play"
......@@ -168,49 +173,7 @@
}
},
"source": [
"## Course of Marielle Simon, simple random walk\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"autoscroll": false,
"ein.hycell": false,
"ein.tags": "worksheet-0",
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"from IPython.display import display, Math\n",
"\n",
"import numpy as np\n",
"import numpy.linalg as LA\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
"from sklearn.utils import check_random_state\n",
"from scipy.special import binom\n",
"import scipy.linalg as la\n",
"import multiprocessing as mp\n",
"mp.set_start_method('spawn', True) # see https://github.com/microsoft/ptvsd/issues/1443\n",
"from numba import jit"
]
},
{
"cell_type": "markdown",
"metadata": {
"ein.tags": "worksheet-0",
"slideshow": {
"slide_type": "-"
}
},
"source": [
"## Random walk on $\\mathbb{Z}$\n",
"\n",
" Consider the random walk on $\\mathbb{Z}$ with $0 < p < 1$, denoted by $(S_n)$. The chain is supposed to start from state 0."
"## Random walk on $\\mathbb{Z}$"
]
},
{
......@@ -222,10 +185,10 @@
}
},
"source": [
"1\\. Implement a function `random_walk_z` simulating the behaviour of the random walk for $n_{\\max}$ steps. #\n",
"Consider the random walk on $\\mathbb{Z}$ with $0 < p < 1$, denoted by $(S_n)$. The chain is supposed to start from state 0.\n",
"\n",
"2\\. Modify the function `random_walk_z` such that it further returns: #\n",
" - both the number of times the chain is returned to the initial state;\n",
"1\\. Implement a function `random_walk_z` simulating the behaviour of the random walk for $n_{\\max}$ steps, and represent it on a graph. Ensure that the function `random_walk_z`also returns: #\n",
" - the number of times the chain is returned to the initial state;\n",
" - the largest state reached by the chain."
]
},
......@@ -242,95 +205,56 @@
},
"outputs": [],
"source": [
"import numpy as np\n",
"from sklearn.utils import check_random_state\n",
"from scipy.special import binom\n",
"import multiprocessing as mp\n",
"mp.set_start_method('spawn', True) # see https://github.com/microsoft/ptvsd/issues/1443\n",
"from numba import jit\n",
"\n",
"@jit(nopython=True)\n",
"def find_first(item, vec):\n",
"def count_first(item: int, vec: np.ndarray):\n",
" \"\"\"\n",
" find_first: find the index of the first element in the array `vec` equal to the element `item`. \n",
"\n",
" :param item: elemtns against which the entries of `vec` are compared\n",
" :type item: int or double\n",
" :param vec: [description]\n",
" :type vec: array-like\n",
" :return: index of the first element in `vec` equal to `item`\n",
" :rtype: int\n",
" \"\"\" \n",
" Find the index of the first element in the array `vec` equal to the element `item`. \n",
" \"\"\" \n",
" c = 0\n",
" for i in range(len(vec)):\n",
" if item == vec[i]:\n",
" return i\n",
" return -1\n",
"\n",
"def random_walk_z(p, X0, n_max, random_state):\n",
" \"\"\" Simulate a simple 1D random walk in Z.\n",
"\n",
" Input:\n",
" ------\n",
" :param p:\n",
" Transition probability (:math:`0 < p <1`)\n",
" :type p:\n",
" float\n",
"\n",
" :param X0:\n",
" Initial state opf the chain.\n",
" :type X0:\n",
" int\n",
"\n",
" :param n_max:\n",
" Maximal number of time steps.\n",
" :type n_max:\n",
" int\n",
" c += 1\n",
" return c\n",
"\n",
" Output:\n",
" -------\n",
" :param random_state:\n",
" Random generator or seed to initialize it.\n",
" :type random_state:\n",
" None | int | instance of RandomState \n",
"\n",
"def random_walk_z(p, n_max, random_state):\n",
" \"\"\" Simulate a simple 1D random walk in Z.\n",
" \n",
" :returns:\n",
" - X (array-like) - trajectory of the chain\n",
" - Ti (:py:class:`int`) - return time to the initial state\n",
" - Ti (:py:class:`int`) - number of returns to the initial state\n",
" - state_max (:py:class:`int`) - farthest state reached by the chain (w.r.t the initial state)\n",
" \"\"\"\n",
"\n",
" rng = check_random_state(random_state)\n",
" Z = 2*rng.binomial(1, p, size=(n_max)) - 1\n",
" X = np.empty(shape=(n_max+1), dtype=float)\n",
" X[0] = X0\n",
" X[1:] = X0 + np.cumsum(Z)\n",
" X[0] = 0\n",
" X[1:] = np.cumsum(Z)\n",
"\n",
" Ti = find_first(0, X[1:]) + 1\n",
" Ti = count_first(0, X[1:])\n",
" id = np.argmax(np.abs(X))\n",
" state_max = X[id]\n",
"\n",
" return X, Ti, state_max"
]
},
{
"cell_type": "markdown",
"metadata": {
"ein.tags": "worksheet-0",
"slideshow": {
"slide_type": "-"
}
},
"source": [
"3\\. Simulate the random walk with $p = 3/4$, and display the histogram of the states reached by the chain. Do the same with $p=1/2$ and illustrate the central limit theorem stated in the lecture, ie: $\\lim_{n\\to\\infty} n^{-1/2}S_n$ is distributed as a standard normal random variable."
" \n",
" t = np.arange(0, n_max+1, 1)\n",
" plt.plot(t, X)\n",
" plt.show()\n",
" return Ti, state_max"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"autoscroll": false,
"ein.hycell": false,
"ein.tags": "worksheet-0",
"slideshow": {
"slide_type": "-"
}
},
"metadata": {},
"outputs": [],
"source": [
"# to do"
"random_walk_z(0.5, 1000, 500)"
]
},
{
......@@ -342,7 +266,7 @@
}
},
"source": [
"4\\. Assume now that two players $A$ and $B$ play heads or tails, where heads occur with probability $p$. Player $A$ bets $1$ euro on heads at each toss, and $B$ bets $1$ euro on tails. Assume that: \n",
"2\\. Assume now that two players $A$ and $B$ play heads or tails, where heads occur with probability $p$. Player $A$ bets $1$ euro on heads at each toss, and $B$ bets $1$ euro on tails. Assume that: \n",
"- the initial fortune of $A$ is $a \\in \\mathbb{N}$;\n",
"- the initial fortune of $B$ is $b\\in\\mathbb{N}$;\n",
"- the gain ends when a player is ruined.\n",
......@@ -363,176 +287,7 @@
},
"outputs": [],
"source": [
"# to do"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Standard coupling"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from percolation import PercolationRect, PercolationHex, percolation_vs_p"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Rectangular lattice"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"percorect = PercolationRect(20, 10)\n",
"interact(percorect.plot, p=widgets.FloatSlider(min=0., max=1., step=0.1, value=0.5));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Hexagonal lattice"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"percohex = PercolationHex(5, 5)\n",
"percohex.compute_clusters(0.5)\n",
"percohex.plot_clusters(add_cluster_id=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"percohex15 = PercolationHex(15, 15)\n",
"interact(percohex15.plot, p=widgets.FloatSlider(min=0., max=1., step=0.1, value=0.5));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Probability of crossing as a function of $p$\n",
"\n",
"Based on 50 simulations on a $30 \\times 30$ lattice."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"percolation_vs_p(30, 30, nsim=50)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Earth movers\n",
"\n",
"Solve the earth movers problem for 50-position samples."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from earth_movers import EarthMovers\n",
"\n",
"em = EarthMovers(50)\n",
"em.plot_ot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot histogram of distance for 2000-position samples."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},