Finished implementing all random walk functions.

Author: rbanerjee20

fdm-devito-notebooks/03_diffu/diffu_rw.ipynb

@@ -592,6 +592,7 @@
     "where $\\alpha = \\frac{1}{2}$.\n",
     "\n",
     "## Implementation of multiple walks\n",
+    "<div id=\"diffu:randomwalk:1D:multiplewalks\"></div>\n",
     "\n",
     "Our next task is to implement an ensemble of walks (for statistics,\n",
     "see the section [Statistical considerations](#diffu:randomwalk:1D:EVar)) and also provide data from\n",
@@ -647,7 +648,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 216,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1060,7 +1061,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 74,
+   "execution_count": 211,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1120,7 +1121,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 72,
+   "execution_count": 212,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1146,12 +1147,12 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "A vectorized version is desired. We follow the ideas from the section [Playing around with some code](#diffu:randomwalk:1D:code1). In fact, the only change to our Devito code is an extra `Dimension` added to our `TimeFunction`, representing the number of space dimensions. This can be represented as follows:"
+    "A vectorized version is desired. We follow the ideas from the section [Playing around with some code](#diffu:randomwalk:1D:code1). We need to add an extra `Dimension` to our `TimeFunction`, representing the number of space dimensions. We also need to generate `N*d` random values, which we can do using the `uniform` function from before and specifying the `size` argument as `N*d`. This can be represented as follows:"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 73,
+   "execution_count": 219,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1165,8 +1166,8 @@
     "    r = TimeFunction(name='r', dimensions=(x_d,d_d), shape=(N+1,d))\n",
     "    r.data[0] = x0\n",
     "    \n",
-    "    rs = random.uniform(0, 1, size=N)\n",
-    "    steps = Function(name='steps', dimensions=(x_d,), shape=(N+1,))\n",
+    "    rs = random.uniform(0, 1, size=N*d).reshape(N,d)\n",
+    "    steps = Function(name='steps', dimensions=(x_d,d_d), shape=(N+1,d))\n",
     "    steps.data[:N] = np.where(rs <= p, -1, 1)\n",
     "    \n",
     "    # The next value is computed using the previous value and the steps function at that point\n",
@@ -1192,100 +1193,54 @@
     "\n",
     "## Multiple random walks in any number of space dimensions\n",
     "\n",
-    "As we did in 1D, we extend one single walk to a number of walks (`num_walks`\n",
-    "in the code).\n",
-    "\n",
-    "### Scalar code\n",
-    "\n",
-    "As always, we start with implementing the scalar case:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 49,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def random_walksdD(x0, N, p, num_walks=1, num_times=1,\n",
-    "                   random=random):\n",
-    "    \"\"\"Simulate num_walks random walks from x0 with N steps.\"\"\"\n",
-    "    d = len(x0)\n",
-    "    position  = np.zeros((N+1, d))   # Accumulated positions\n",
-    "    position2 = np.zeros((N+1, d))   # Accumulated positions**2\n",
-    "    # Histogram at num_times selected time points\n",
-    "    pos_hist = np.zeros((num_walks, num_times, d))\n",
-    "    pos_hist_times = [(N//num_times)*i for i in range(num_times)]\n",
-    "\n",
-    "    for n in range(num_walks):\n",
-    "        num_times_counter = 0\n",
-    "        current_pos = np.array(x0, dtype=float)\n",
-    "        for k in range(N):\n",
-    "            if k in pos_hist_times:\n",
-    "                pos_hist[n,num_times_counter,:] = current_pos\n",
-    "                num_times_counter += 1\n",
-    "            # current_pos corresponds to step k+1\n",
-    "            for i in range(d):\n",
-    "                r = random.uniform(0, 1)\n",
-    "                if r <= p:\n",
-    "                    current_pos[i] -= 1\n",
-    "                else:\n",
-    "                    current_pos[i] += 1\n",
-    "            position [k+1,:] += current_pos\n",
-    "            position2[k+1,:] += current_pos**2\n",
-    "    return position, position2, pos_hist, np.array(pos_hist_times)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Vectorized code\n",
-    "\n",
-    "Significant speed-ups can be obtained by vectorization. We get rid of\n",
-    "the loops in the previous function and arrive at the following vectorized\n",
-    "code."
+    "As we did in 1D, we extend one single walk to a number of walks (`num_walks` in the code). We can take a look at our function for multiple random walks (`random_walks1D`) from the [Implementation of multiple walks](#diffu:randomwalk:1D:multiplewalks) section. The changes we make here are similar to those we made for `random_walkdD_vec`, in that we add an extra `Dimension` to represent the number of space dimensions."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 76,
+   "execution_count": 218,
    "metadata": {},
    "outputs": [],
    "source": [
     "def random_walksdD_vec(x0, N, p, num_walks=1, num_times=1, random=random):\n",
-    "    \"\"\"Vectorized version of random_walks1D; no loops.\"\"\"\n",
+    "    \"\"\"Vectorized version of random_walksdD.\"\"\"\n",
     "    d = len(x0)\n",
     "    position  = np.zeros((N+1, d))  # Accumulated positions\n",
     "    position2 = np.zeros((N+1, d))  # Accumulated positions**2\n",
-    "    walks = np.zeros((num_walks, N+1, d))  # Positions of each walk\n",
-    "    walks[:,0,:] = x0\n",
     "    # Histogram at num_times selected time points\n",
     "    pos_hist = np.zeros((num_walks, num_times, d))\n",
     "    pos_hist_times = [(N//num_times)*i for i in range(num_times)]\n",
     "\n",
-    "    r = random.uniform(0, 1, size=N*num_walks*d)\n",
-    "    steps = np.where(r <= p, -1, 1).reshape(num_walks, N, d)\n",
-    "    walks[:,1:,:] = x0 + np.cumsum(steps, axis=1)\n",
-    "    position  = np.sum(walks,    axis=0)\n",
-    "    position2 = np.sum(walks**2, axis=0)\n",
-    "    pos_hist[:,:,:] = walks[:,pos_hist_times,:]\n",
+    "    # Create and initialise our TimeFunction r\n",
+    "    x_d = Dimension(name='x_d')\n",
+    "    t_d = Dimension(name='t_d')\n",
+    "    d_d = Dimension(name='d_d')\n",
+    "    r = TimeFunction(name='r', dimensions=(x_d, t_d, d_d), shape=(N+1, num_walks, d))\n",
+    "    # Setting each walk's first element to x0\n",
+    "    r.data[0] = x0\n",
+    "    \n",
+    "    steps = Function(name='steps', dimensions=(x_d,t_d,d_d), shape=(N+1,num_walks,d))\n",
+    "        \n",
+    "    # Populating steps with -1 if value in rs <= p at that point and 1 otherwise\n",
+    "    rs = random.uniform(0, 1, size=N*num_walks*d).reshape(num_walks, N, d)\n",
+    "    for n in range(num_walks):\n",
+    "        steps.data[:N, n] = np.where(rs[n] <= p, -1, 1)\n",
+    "    \n",
+    "    # Creating and applying operator that contains equation for adding steps\n",
+    "    eq = Eq(r.forward, r + steps)\n",
+    "    op = Operator(eq)\n",
+    "    op.apply()\n",
+    "    \n",
+    "    # Summing over data to give positions\n",
+    "    position = np.sum(r.data, axis=1)               # Accumulated positions\n",
+    "    position2 = np.sum(r.data**2, axis=1)           # Accumulated positions**2\n",
+    "    \n",
+    "    for k in range(num_walks):\n",
+    "        pos_hist[k,:] = r.data[pos_hist_times,k]\n",
+    "    \n",
     "    return position, position2, pos_hist, np.array(pos_hist_times)"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 77,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
   {
    "cell_type": "markdown",
    "metadata": {},