02_separation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Separation\n",
    "\n",
    "Let's think about the essential pieces here:\n",
    "- a periodic 2d domain\n",
    "    - we need to map positions back to the central coordinate range\n",
    "- a group of particles\n",
    "    - particle positions need to be tracked\n",
    "    - particles need to be moved around\n",
    "    - we need to diagnose center of mass\n",
    "    - we need to diagnose moment of inertia\n",
    "- a random number generator\n",
    "    - we want to initialize it in a reproducible way (see rule 6 from [Sandve et al. (2013)](https://doi.org/10.1371/journal.pcbi.1003285))\n",
    "    - we want to draw arrays of random numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tech preamble"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The spatial domain\n",
    "\n",
    "Has info on its size (`length_x`, `length_y`), has a method to return these sizes, has a method to normalize given positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class PeriodicSpace:\n",
    "    def __init__(self, length_x=10, length_y=20):\n",
    "        self.length_y = length_y\n",
    "        self.length_x = length_x\n",
    "    \n",
    "    def get_sizes(self):\n",
    "        return self.length_x, self.length_y\n",
    "        \n",
    "    def normalize_positions(self, x, y):\n",
    "        return np.mod(x, self.length_x), np.mod(y, self.length_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_space = PeriodicSpace()\n",
    "\n",
    "# ensure that positions within the base interval\n",
    "# are not changed\n",
    "lx, ly = test_space.get_sizes()\n",
    "centered_x, centered_y = lx / 2, ly / 2\n",
    "assert (centered_x, centered_y) == test_space.normalize_positions(x=centered_x, y=centered_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A group of particles\n",
    "\n",
    "- track state of all particles (positions `x`, `y`)\n",
    "- method to move particles\n",
    "- methods to diagnose\n",
    "  - center of mass\n",
    "  - moment of inertia\n",
    "- method to wrap diagnostics in pandas dataframe\n",
    "- method to wrap positions in a pandas dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Particles:\n",
    "    def __init__(\n",
    "        self,\n",
    "        rng=np.random.RandomState(),\n",
    "        space=PeriodicSpace(),\n",
    "        x=None, y=None,\n",
    "        step_length=0.5\n",
    "    ):\n",
    "        self.rng = rng\n",
    "        self.space = space\n",
    "        self.x, self.y = x, y\n",
    "        self.step_length = step_length\n",
    "        self.steps_done = 0\n",
    "\n",
    "    def move(self):\n",
    "        self.x += self.step_length * self.rng.normal(size=self.x.shape)\n",
    "        self.y += self.step_length * self.rng.normal(size=self.y.shape)\n",
    "        \n",
    "        self.x, self.y = self.space.normalize_positions(self.x, self.y)\n",
    "        \n",
    "        self.steps_done += 1\n",
    "\n",
    "    def center_of_mass(self):\n",
    "        return self.x.mean(), self.y.mean()\n",
    "    \n",
    "    def moment_of_inertia(self):\n",
    "        return self.x.var() + self.y.var()\n",
    "\n",
    "    def diagnostics(self):\n",
    "        com = self.center_of_mass()\n",
    "        mi = self.moment_of_inertia()\n",
    "        return pd.DataFrame(\n",
    "            {\n",
    "                \"center_of_mass_x\": com[0],\n",
    "                \"center_of_mass_y\": com[1],\n",
    "                \"moment_of_inertia\": mi\n",
    "            },\n",
    "            index=[self.steps_done, ],        \n",
    "        )\n",
    "    \n",
    "    def positions(self):\n",
    "        return pd.DataFrame(\n",
    "            {\n",
    "                \"x\": self.x,\n",
    "                \"y\": self.y\n",
    "            }\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random number generator\n",
    "\n",
    "There is already an object with all we need: `numpy.random.RandomState`.  So we'll use this for our RNG."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Drive the experiment\n",
    "\n",
    "Parameters we want to be able to control:\n",
    "- number of particles\n",
    "- number of steps to run the model\n",
    "- size of the spatial domain\n",
    "- step size in the updates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_random_walk(\n",
    "    length_x=10,\n",
    "    length_y=20,\n",
    "    number_particles=100,\n",
    "    number_steps=100,\n",
    "    step_length=0.5\n",
    "):\n",
    "    # initialize spatial domain\n",
    "    space = PeriodicSpace(\n",
    "        length_x=length_x,\n",
    "        length_y=length_y\n",
    "    )\n",
    "    \n",
    "    # initialize random number generator\n",
    "    rng = np.random.RandomState()\n",
    "    \n",
    "    # initialize group of particles and register spatial domain and rng\n",
    "    particles = Particles(\n",
    "        space=space,\n",
    "        rng=rng,\n",
    "        x=np.ones((number_particles, )) * length_x / 2.0,\n",
    "        y=np.ones((number_particles, )) * length_y / 2.0,        \n",
    "        step_length=step_length\n",
    "    )\n",
    "    \n",
    "    # track initial positions\n",
    "    initial_positions = particles.positions()\n",
    "    \n",
    "    # initialize dataframe with diagnostics\n",
    "    diags = particles.diagnostics()\n",
    "    \n",
    "    # run `number_steps` updates and accumulate diagnostics\n",
    "    for step in range(1, number_steps):\n",
    "        particles.move()\n",
    "        \n",
    "        diags = diags.append(\n",
    "            particles.diagnostics(),\n",
    "            ignore_index=True\n",
    "        )\n",
    "    \n",
    "    # track final positions\n",
    "    final_positions = particles.positions()\n",
    "        \n",
    "    return diags, initial_positions, final_positions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "    diags,\n",
    "    initial_positions,\n",
    "    final_positions\n",
    ") = run_random_walk(number_particles=100, number_steps=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data and visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>5.0</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>5.0</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>5.0</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>5.0</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     x     y\n",
       "0  5.0  10.0\n",
       "1  5.0  10.0\n",
       "2  5.0  10.0\n",
       "3  5.0  10.0\n",
       "4  5.0  10.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>center_of_mass_x</th>\n",
       "      <th>center_of_mass_y</th>\n",
       "      <th>moment_of_inertia</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   center_of_mass_x  center_of_mass_y  moment_of_inertia\n",
       "0               5.0              10.0                0.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>9.491506</td>\n",
       "      <td>18.354704</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>3.822715</td>\n",
       "      <td>13.159536</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>5.388253</td>\n",
       "      <td>0.309439</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>4.807517</td>\n",
       "      <td>12.536099</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>1.886733</td>\n",
       "      <td>10.880113</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          x          y\n",
       "0  9.491506  18.354704\n",
       "1  3.822715  13.159536\n",
       "2  5.388253   0.309439\n",
       "3  4.807517  12.536099\n",
       "4  1.886733  10.880113"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>center_of_mass_x</th>\n",
       "      <th>center_of_mass_y</th>\n",
       "      <th>moment_of_inertia</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>999</td>\n",
       "      <td>5.015665</td>\n",
       "      <td>9.934202</td>\n",
       "      <td>39.863907</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     center_of_mass_x  center_of_mass_y  moment_of_inertia\n",
       "999          5.015665          9.934202          39.863907"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(initial_positions.head(5))\n",
    "display(diags.head(1))\n",
    "\n",
    "display(final_positions.head(5))\n",
    "display(diags.tail(1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'init' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-8-16352240deee>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mdisplay\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mdisplay\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdiags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'init' is not defined"
     ]
    }
   ],
   "source": [
    "display(init)\n",
    "display(diags)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_positions.plot.scatter(x=\"x\", y=\"y\")\n",
    "final_positions.plot.scatter(x=\"x\", y=\"y\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diags.plot();"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:climnum-exercise]",
   "language": "python",
   "name": "conda-env-climnum-exercise-py"
  },
  "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}