update notebook

usnistgov · Nov 9, 2020 · 6b33621 · 6b33621
1 parent 665ed7c
commit 6b33621
Showing 1 changed file with 71 additions and 0 deletions.
diff --git a/benchmarks/benchmark8.ipynb b/benchmarks/benchmark8.ipynb
@@ -969,6 +969,77 @@
     "\n"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Appendix\n",
+    "\n",
+    "The following code is given as a resource for those implementing the benchmark to help count the particles in Part (b) and Part (c). One possibility as outlined here is to use [`scipy.ndimage.label`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.label.html#scipy.ndimage.label). This returns a tuple with a labeled array of unique features as the first item and a count of features as the second item. In the case of the phase field the values will require thresholding to 0 or 1 before using `scipy.ndimage.label`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 138,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Actual number of particles: 14\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAADKCAYAAAC11LviAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAEKdJREFUeJzt3X2sJXV5wPHvI+zuZaGoC8vChkQKhdpKZEtXkLZa01VWNm3TxheoosWoa4OxiS1ohKZC04QatDRotVmb1r5YUSBgSQoXWUqlGmkWukWUZQG7lrhuWWWl8iIqPP3jzMXj8dx75p6XmTlzvp/kZmd+d2bPc+fMPOeZ329mTmQmkqTp95y6A5AkjYcJXZJawoQuSS1hQpekljChS1JLmNAlqSVM6JLUEiMl9IjYEhH3RMR9EXHRuIKSJC3f0Ak9Ig4FPga8EngRcFZEnDquwCRJy3PwCOueBtyVmfsAIuIaYAtw12IrrIxVOcehI7wknPTiJ0Zaf1i7715dy+tKs2Dcx3XbjtfvcuBbmbl20HKjJPT1wMNd8/uBE3sXioitwFaAOVZzemwa4SVhfn7nSOsPa/P6DbW8rjQLxn1ct+14vSWv+XqZ5UZJ6ADP9Myv7F0gM7cB2wAOjzWteHDM/N7l7Xxt27kkNdMoCX0fcGTX/NqirZWWm8QXW9fkLk3ewjE3a8fbKFe53AG8JCKOioiDgdcC28cTliRpuYau0DPzsYh4F/CvwArgHzPz38YWWY9RKuQmmdXKQdLkjdSHnpk3ADeMKRZJ0ghGHRSVpMaZ1TNgb/2vyfzena3pRpLUDCZ0SWqJqUnom9dvaOVplFW61N7ju2pTk9AlSUszoUtqjFGr9Fmv9E3okhpl2KQ8y4l8gQldklpi6q5D37x+gwOJ0gzorrh7j3mr8f6s0CWpJaauQpc0e6zIy5nKCn3WR7IlqZ+pTOiSpJ801Ql9UpV6lWcAnmlIGpepTuiSpB9pxaDoUpc3DfN/SNI0skKXpJZoRYXebZyV9iRvYvKMQNK4tS6hj9tC4h1XYjeRS5oUu1wkqSVM6CWN41JGq3NJk2RCl6SWsA99mZZziaQVuaQqmdBHYMKW1CQDu1wiYi4ibomIByNid0RcVLQfHxFfLNo+FRFzkw9XkrSYsn3oH8jME4BTgLMjYgPw18ClmXkSsAc4fzIhSpLKGJjQM/N7mfm5YvpJ4AFgHXAycHOx2FXAlkkFKUkabFlXuUTEOuClwD3AgczM4lf7gaPHHJskaRlKD4pGxCrgauDioumZnkVWLrLeVmArwByrhwhRklRGqQo9IlYC1wI3ZuYn6FTkz+taZC2wr9+6mbktMzdm5sYVrBoxXEnSYspc5bIauAG4PTMvA8jM7wP3RcSmYrFzgO0Ti1KSNFCZLpfTgFcAL4iItxRt1wFvA/4hIv4KuAt4S//VJUlVGJjQM/M2WLSv5IyxRiNJGprPcpGkljChS1JLmNAlqSVM6JLUEiZ0SWoJH58rTbne5/L7WOfZZYUuSS1hQpemWL9vzRr0TVpqLxO6JLWEfejSFLIKVz9W6NIU2rx+g4Of+gkmdElqiUq7XE568RPMz/c/VbTakKTRWKFLUks0ZlC0e5DHal0qx2NF3RpZoc/v3ekoviQtUyMTuiRp+Rqd0K3SJam8Rid0SVJ5jU/o9qdLUjmNucpFmoR+xYBXhqitGl+hS5LKMaGrteyq06wxoUtSSywroUfEhRFxTzF9RETcFBG7i3/XTCbEDgdHJWlppRN6RPwy8IaupsuB6zLzJOA64JLxhiZJWo5SCT0ijgSuAH6vq3kT8Oli+ipgy3hDk0bT72oWr3BRmw28bDEiAvg74D3A/3b96ojM/A5AZj466S6XpQzqivEgnl2+95olZSr0dwNfzMzbetqzZ35lv5UjYmtE7IiIHfu//fQQIUqSyihzY9FPA2dGxJuAFcCxEXE7cCAiDsvMxyLiucAj/VbOzG3ANoCNp8z1fgiU1l1pLXdw1EfzSpoFAyv0zHxXZv5sZr6QTr/5/Zn5MuBW4OxisXOA7ZMLU5I0yCi3/l8IfDIi3gvsAd44logGGPXSxd71rdgltcWyEnpm7gFOLqb3A2dOIKa+JnUN+sL/a2KXNO28U1SSWsKELkktYUKXpJbweeiF+b077UfX0HzuupqgkQl94UDwYVxqorL7pfc/qGp2uUhSSzSmQreCUdONcsZota4qWKFLUktUWqHvvnt1qeqkrr7zum4yWurvtZqrn2M5mhaN6XKZNcMMrC0wyU83707WpNjlIkktYYVesXGcvlvhSerHCl2SWiIyh/7OiWU7PNbk6bGp9PJVD0ZNuuKdxN9jlT5ZVeyDvoca5Ja85s7M3DhoObtcKjDJpGD3i6QFdrlIUkuY0CWpJUzoktQSjU7oVfYLT+K15vfu9C7DKTfJfXDz+g2OfWisGj8oOulH6bblgHJwVFKjK3RJUnmNr9AXjLtSt5JVWe57mhZW6JLUElOX0EetbhyIktRWU9Pl0q1fQp71x8zO0t9al+5tvNzuF98fVaFUhR4RKyLiQxHxQEQ8FBHPj4jjI+KLEbE7Ij4VEXOTDlaStLiyFfpHgW8CJ3a1XQtcmpnzEXEZcD7w52OOr7QmVkCTvuRS9end32b9DFHNMLBCj4ijgV8CLskCnQ+Ck4Gbi8WuArZMLEpJ0kBlKvSTgQRujYhjgB3Ae4AD+aNn7+4Hjp5MiFLzWY2rCcok9KOA3cDvAD8ELgcuAZ7pWW5lv5UjYiuwFWCO1cPGOdUm2fViIpG0oMyg6AHg8cx8KjOfBq4HjgWe17XMWmBfv5Uzc1tmbszMjStYNXLAkqT+yiT0LwAvj4jjivmzirb7Ip79+qFzgO1jj05LsjqX1G1gQs/M/wPeCnw2Ir5KpwvmcuBtwJ9GxP3A8UWbJKkmpS5bzMxbgFN6mh8Azhh7RBrIylxSP1N3678kqb+pvPV/FlmVSxrEhF4B7xSVVAW7XCSpJUzoU8IqX9IgJnRJagkTuiS1hAldklrChC5JLWFCnxJehy5pEBO6JLWECb0Cm9dvsMKWNHEm9JLm9+4c+VrwYZK6HwaSyjKhS1JL+CyXJfSryP12d0lNZYUuSS1hhd7HcvvK5/futEqXJmi5x+SsHo8m9C6jDHourDurO5I0KcMcl93rzNIxaZeLJLWEFbqkRhrXI6Nn6ezZCl2SWsIKveAXSEjN4LE4PBN6xRZO+wbttOM+PfT6eU2DSSbzWeh6sctFklqiVIUeEb8LXAisBO4GzgNWAZ8Ejge+BrwhMx+ZTJiTt3n9hkpP9aqoEgb9Pb2/b3PlIs2CgRV6RKwD3g+ckZknAQ8D7wIuB64r2q4DLplgnJKkAcp0uawEDgUOK+b3Ad8HNgGfLtquAraMPToNZdgnQ47jiZJS07V5Hx/Y5ZKZD0XEFcC9EXENsA54PXBpZn6nWObRiFgz2VAnr+yA5VLrSlJdynS5PBf4TeAMYJ5On/mvAdmz6MpF1t8aETsiYscPeGrEcCVJiykzKPoq4N7MvJdOlf4Y8E7gQEQclpmPFUm/74BoZm4DtgEcHmt6PwQaabmVepOq83GcTvqwMam8Jj03pkwf+teAl3V1qWwEdgG3AmcXbecA28cfniSprDJ96HdFxEeAL0XE08BOYCswB3wyIt4L7AHeOMlA61D3p62k8Rv3cd2kPFHqOvTMvBK4sqf5u8CZY49IkjQU7xSVpJbwWS4NMo7LJcdxx2uTTiE1W6q+Y7ttrNAlqSWs0Gs2ajXi81ikcmbh2DCh12RSp5Xd15D7xbqaRqPcsb3U/zcL7HKRpJawQq9YFQM+/R7k7xdcaNqMo1KftX3cCl2SWsIKfUbMWqWi9hh0prnYsrPICl2SWsIKvUJV3zAxC1+Kq9nivrw0K3RJagkTuiS1hAldklrChC5JLWFCl6SW8CoXSeqy1NVoTb/KxgpdklrCCl2N0aRvT9fsKXOfSNPv7bBCl6SWsEJXbZaqiKzWVaXl3sXd/b0DTWJCVy2WcwA1/TRX06tt319ql4sktYQJvUJVV5ib129oVVU7v3dn6yoqaZxM6JLUEiZ0STOrbWexDopWbNzfaL7Ua0iaLVboktQSkZnVvVjEfuBx4FuVvejwjqT5cU5DjGCc42ac4zUNcb4gM9cOWqjShA4QETsyc2OlLzqEaYhzGmIE4xw34xyvaYmzDLtcJKklTOiS1BJ1JPRtNbzmMKYhzmmIEYxz3IxzvKYlzoEq70OXJE2GXS6S1BKVJfSI2BIR90TEfRFxUVWvO0hEzEXELRHxYETsXogtIi6JiIcjYlfx888NiPW2iNjTFdMfRcQREXFTEftNEbGm5hhP6YpvV0Q8UMR9XkQc6Gq/s6b4To2Iu7vm+26/iHhORHy4aP/PiDi15jgvKLblroi4MSLWFu3HRcRTPdv8RTXGuej7HBEXF8f/PRFxVo0x/kvP9no8Io4rfvdMz+/OrCrOscjMif8AhwJfB46mc3fq7cCpVbx2idjmgFcV04cA/wVsAC4BLqg7vp5YbwM29rT9DfCOYvodwJV1x9kT31bgCuA84CM1x/Ih4NvAPYO2H/Bm4J+K6RcBd9Uc5yuB1cX0RcAVxfRx3cs1YHv2fZ+BlwP/DhwEHAPsBlbUEWPP79cUsawq5h+rY1uO66eqCv00OgfEvsz8IXANsKWi115SZn4vMz9XTD8JPACsqzeqZdkEfLqYvoqGbFeAiDgY+APgg3XHApCZfwj8Yk/zYtvv2fbM/AoQEXFsXXFm5i2Z+UQx+2U6xVGtFtmei9kEXJ2ZT2fmN4GvAKdPLLhCiRjfDXw8M5+adCxVqCqhrwce7prfTwN2yF4RsQ54KXBH0XRhRNwfEddHxPoaQ1uQwDXFaeuVRcI8IjO/A5CZj9KpOJriTcDnM/Mbxfwbiu35uYj4+ToD67LY9mvyPnsusL1r/oRiu94ZEb9VV1Bd+r3PjdueEfE84I3Ax7qa54rYvxwRb68ptKFVOSj6TM/8ygpfe6CIWAVcDVxcHOB/lpnrgJPonCr+RZ3xFc7KzOOAX6Bz2rqVTpLv1ojtGhEHAe8BPlA0fYpO8jwR+DjwiZpC67XU9mvcPhsR5wNHAH9bND0EHF5s13OBj0bE8+uKj6Xf56Ztz98H/j4zH+tqO6yI/dXABVWOR4xDVQl9H53nJSxYW7Q1QkSsBK4FbszMT0CnK6b4N4HPAD9XW4CFrpieAG6gE9OBiDgMICKeCzxSX4Q/5hw63WwPAmTmU8W2hE6X24m1RfbjFtt+jdtnI+LNdM56XpOZTwMUXRg/KKbvBXYBx9cV4xLvc6O2Z0T8FPBW4Mru9q5j7BvAF4AXVh/d8KpK6HcAL4mIo4pugtfy46eMtYmI1XSS4+2ZeVlX+6YiVoDXA1+qI76ueOYi4hXF9Argt4uYbgXOLhY7hwZs14h4DvA+oHt7vjwiDilmXwPsqCO2PhbbftuLeYoq7dDM/Fr14XVExFbg7XTO0h7tan9xRBxTTJ9AJ5nfV0+US77P24HXRcRBRbynAv9RR4yFdwKfycxnC6CIOCEifqaYXgv8ClDL1VhDq2r0FfgNOgMhu4E/rns0uCuuVwBP0alsFn4uAz5M58qc+4DrgbU1x3kI8Hngv4uYPkjnA3ktcHOxXW+uO84i1tcBn+1pex+wp9i+24Hja4jrT4C7gSfpJJpfXWz70bka4y+L9p3AaTXHuadr++0CdhXLvrqYv7+I88ya41z0fQbeX+y7XwV+vcYYVwP/AxzTs+ypxbIPFLnq3Kr30VF/vFNUklrCO0UlqSVM6JLUEiZ0SWoJE7oktYQJXZJawoQuSS1hQpekljChS1JL/D/klfMYaYnqyQAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x7f23c00f4fd0>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "import scipy.ndimage\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "def make_particles(shape, n_particles, max_radius, seed=99):\n",
+    "    \"\"\"Make an array with random particles\n",
+    "    \n",
+    "    Args:\n",
+    "      shape: shape of array\n",
+    "      n_particles: number of particles in arrau\n",
+    "      max_radius: maximum possible radius of particles\n",
+    "      seed: random seed\n",
+    "      \n",
+    "    Returns: numpy array mask with particles\n",
+    "    \"\"\"\n",
+    "    np.random.seed(seed)\n",
+    "        \n",
+    "    def xy_grid(shape):\n",
+    "        \"\"\"Generate a grid of x, y values in a single array\n",
+    "        \n",
+    "        The returned array has shape (len(shape),) + shape + (1,). The last axis is\n",
+    "        for the operations with particles.\n",
+    "        \"\"\"\n",
+    "        linspace = lambda x: np.linspace(0, x - 1, x) + 0.5\n",
+    "        return np.array(np.meshgrid(*map(linspace, shape), indexing=\"ij\"))[..., None]\n",
+    "    \n",
+    "    # first axis is dimensions\n",
+    "    # second and third axes are shape of domain\n",
+    "    # last axis is the number of particles. \n",
+    "    centers = (np.random.random((len(shape), n_particles)) * np.array(shape)[:, None])[:, None, None]\n",
+    "    radii = (np.random.random(n_particles) * max_radius)[None, None]\n",
+    "    return numpy.any(((xy_grid(shape) - centers)**2).sum(0) < radii**2, axis=-1)\n",
+    "\n",
+    "arr = make_particles((100, 200), 20, 10)\n",
+    "plt.imshow(arr)\n",
+    "print('Actual number of particles:', scipy.ndimage.label(arr)[1])"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,