{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Escape Time Computation\n\nMean Exit Time For a SDE \n\n\\begin{align}\begin{equation}\begin{pmatrix}{dx \\ dy}\\end{pmatrix} = \begin{pmatrix}1 & 0 \\\n 0 & 1\\end{pmatrix}dW_t \\end{equation}\\end{align}\n\n\\ \n\nBy using two approaches:\n\n- kinetic Monte Carlo Simulations\n- solving the Boundary Value Problem with Finite Elements\n\nWe will illustrate that numerically computed Mean Exit Times for this 2D\nexample apparoximate the theoretical value of the exit time that is\n\\ $0.5$\\ \\* Oksendal, Bernt. Stochastic differential equations:\nan introduction with applications. Springer Science & Business Media,\n2013.\n\n## Dependencies\n\n- numpy pip install numpy\n- matplotlib pip install matplotlib\n- tqdm pip install tqdm\n- shapely pip install Shapely\n- joblib pip install joblib\n- seaborn pip install seaborn\n- fenics (for solving the Boundary Value Problem)\n [Installation_Instruction](https://fenicsproject.org/download/archive/)_\n\n.. code:: ipython3\n\n import numpy as np\n import matplotlib.pyplot as plt\n from tqdm import tqdm\n from shapely.geometry import Polygon, Point, MultiPoint\n from joblib import Parallel, delayed\n import seaborn as sns\n \n plt.rc('xtick', labelsize=20)\n plt.rc('ytick', labelsize=20)\n plt.rc('axes', labelsize=21)\n plt.rc('font', family='serif')\n plt.rc('font', family='serif')\n plt.rcParams['image.cmap'] = 'Spectral'\n np.random.seed(2)\n rng = np.random.default_rng(5)\n\n### Define Equations For the Circle\n\n.. code:: ipython3\n\n def equation(theta,r):\n x = np.cos(theta)*r\n y = np.sin(theta)*r\n return x,y\n\n## The user chooses the number of points to sample on the circle\n\n.. code:: ipython3\n\n N = int(input(\"Enter Number of Points to Create the Circle: \"))\n T = np.linspace(0,360,N)*np.pi/180\n Circle = np.array(equation(T,1)).T\n\n\n.. parsed-literal::\n\n Enter Number of Points to Create the Circle: 1000\n\n\n.. code:: ipython3\n\n fig = plt.figure(figsize=(4,4))\n ax = fig.add_subplot(111)\n ax.scatter(Circle[:,0],Circle[:,1],c='k',s=5)\n ax.set_xlabel('x')\n ax.set_ylabel('y')\n \n ##The Circle is saved and will be used for the Finite Difference Code as well. \n np.savetxt('Circle_Boundary.csv',Circle,delimiter=',')\n\n\n\n![](\"file://output_8_0.png\")
\n\n\n### The Shapely Library is used to find based on the data sampled on the circle which are inside the circle and which are outside.\n\nThis will help as to identify when a trajectory escapes from the circle.\nA grid of points is created and we test if the points are inside or\noutside of the circle\n\n.. code:: ipython3\n\n poly = Polygon(Circle)\n X,Y = np.meshgrid(np.linspace(-1.1,1.1,200),np.linspace(-1.1,1.1,200))\n meshpoint = np.c_[X.reshape(-1,1),Y.reshape(-1,1)]\n meshpoint_ = MultiPoint(meshpoint)\n\nTest if a point in the grid is inside or outside of the limit cycle.\n\n.. code:: ipython3\n\n judge_list = []\n \n for i in tqdm(range(len(meshpoint))):\n result = poly.contains(meshpoint_.geoms[i])\n judge_list.append(result)\n\n\n.. parsed-literal::\n\n 100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 40000/40000 [00:01<00:00, 20398.44it/s]\n\n\nSimulate a trajectory until it escapes from the circle starting from the\norigin.\n\n.. code:: ipython3\n\n t = 0.0\n dt=1e-5\n traj =[]\n x = np.zeros(2)\n traj.append(x)\n while True:\n traj.append(x)\n if poly.contains(Point(x)) == False:\n break\n traj.append(x)\n x = x+ np.sqrt(dt) * np.random.randn(2)\n traj = np.array(traj)\n\n### We visualize the points inside and outside of the circle based on Shapely and also a trajectory (starting from the origin) integrated with Euler-Maruyama until it escapes.\n\n.. code:: ipython3\n\n judge_list = np.array(judge_list)\n \n inside_pts = meshpoint[judge_list == True]\n outside_pts = meshpoint[judge_list == False]\n fig = plt.figure(figsize=(4,4))\n ax = fig.add_subplot(111)\n ax.scatter(inside_pts[:, 0], inside_pts[:, 1], color='r', label='Interior',s=2,zorder=0)\n ax.scatter(outside_pts[:, 0], outside_pts[:, 1], color='b', label='Exterior',s=2,zorder=0)\n ax.scatter(Circle[:,0],Circle[:,1],c='k',s=5)\n ax.plot(traj[:,0],traj[:,1],'g-',label='Trajectory',linewidth=0.2)\n ax.plot(0,0,'ko',markersize=5)\n ax.set_xlim(-1.1,1.1)\n ax.set_ylim(-1.1,1.1)\n ax.set_xticks([-1,0,1])\n ax.set_yticks([-1,0,1])\n ax.set_xlabel('x')\n ax.set_ylabel('y')\n ax.legend(fontsize=15,frameon=True,loc = 'upper right')\n plt.tight_layout()\n\n\n\n![](\"file://output_16_0.png\")
\n\n\n### kinetic Monte Carlo run Experiments - Parallel\n\nInstructions: \\* To compute the escape time with kinetic Monte Carlo the\nuser specifies the number of experiments (trajectories) to simulate and\nalso the time step (dt).\n\nSuggestions: \\* Use a time step dt~1e-5 \\* Run about 10,000 experiments\nto get a good sample of the statistics (this might involve running for\n~20 mins)\n\n.. code:: ipython3\n\n n = int(input(\"Enter Number of Experiments: \")) # For 10,000 experiments 10,000 (needs ~22mins in 24 cores)\n dt: float = float(input(\"Select Time Step: \")) # Time step length (suggested time step ~1e-5)\n \n def integrate_until_exit(random_seed):\n np.random.seed(random_seed)\n \n t = 0.0\n x = np.zeros(2)\n \n while True:\n if poly.contains(Point(x)) == False:\n # if np.linalg.norm(x) >= 1.0: \n ### For this toy example we can use also the above as a stopping condition.\n break\n x += np.sqrt(dt) * np.random.randn(2)\n t += dt\n \n return t\n \n \n fpts = Parallel(n_jobs=-1, verbose=1)(\n delayed(integrate_until_exit)(i) for i in (range(0,n)))\n \n print('MFPT', np.mean(fpts), '+/-', np.std(fpts))\n \n fpts = np.array(fpts)\n\n\n.. parsed-literal::\n\n Enter Number of Experiments: 100\n Select Time Step: 1e-3\n\n\n.. parsed-literal::\n\n [Parallel(n_jobs=-1)]: Using backend LokyBackend with 12 concurrent workers.\n [Parallel(n_jobs=-1)]: Done 26 tasks | elapsed: 1.0s\n\n\n.. parsed-literal::\n\n MFPT 0.5049199999999969 +/- 0.34225363927939706\n\n\n.. parsed-literal::\n\n [Parallel(n_jobs=-1)]: Done 77 out of 100 | elapsed: 1.4s remaining: 0.4s\n [Parallel(n_jobs=-1)]: Done 100 out of 100 | elapsed: 1.6s finished\n\n\n.. code:: ipython3\n\n fig = plt.figure(figsize=(4*2,4))\n ax = fig.add_subplot(121)\n sns.kdeplot(fpts,fill=True);\n\n\n\n![](\"file://output_19_0.png\")
\n\n\n### Imports the Code that Solves the Boundary Value Problem with Finite Elements\n\n.. code:: ipython3\n\n import BVP_Stationary as BV\n\n\n.. parsed-literal::\n\n Solving linear variational problem.\n MFPT: 0.49996583014439444\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.7.9"
}
},
"nbformat": 4,
"nbformat_minor": 0
}