{ "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": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
xy
05.010.0
15.010.0
25.010.0
35.010.0
45.010.0
\n", "
" ], "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": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
center_of_mass_xcenter_of_mass_ymoment_of_inertia
05.010.00.0
\n", "
" ], "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": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
xy
09.49150618.354704
13.82271513.159536
25.3882530.309439
34.80751712.536099
41.88673310.880113
\n", "
" ], "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": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
center_of_mass_xcenter_of_mass_ymoment_of_inertia
9995.0156659.93420239.863907
\n", "
" ], "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\u001b[0m in \u001b[0;36m\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 }