Separation¶

Let's think about the essential pieces here:

• a periodic 2d domain
• we need to map positions back to the central coordinate range
• a group of particles
• particle positions need to be tracked
• particles need to be moved around
• we need to diagnose center of mass
• we need to diagnose moment of inertia
• a random number generator
• we want to initialize it in a reproducible way (see rule 6 from Sandve et al. (2013))
• we want to draw arrays of random numbers

Tech preamble¶

In [1]:
%matplotlib inline
import numpy as np
import pandas as pd


The spatial domain¶

Has info on its size (length_x, length_y), has a method to return these sizes, has a method to normalize given positions.

In [2]:
class PeriodicSpace:
def __init__(self, length_x=10, length_y=20):
self.length_y = length_y
self.length_x = length_x

def get_sizes(self):
return self.length_x, self.length_y

def normalize_positions(self, x, y):
return np.mod(x, self.length_x), np.mod(y, self.length_y)

In [3]:
test_space = PeriodicSpace()

# ensure that positions within the base interval
# are not changed
lx, ly = test_space.get_sizes()
centered_x, centered_y = lx / 2, ly / 2
assert (centered_x, centered_y) == test_space.normalize_positions(x=centered_x, y=centered_y)


A group of particles¶

• track state of all particles (positions x, y)
• method to move particles
• methods to diagnose
• center of mass
• moment of inertia
• method to wrap diagnostics in pandas dataframe
• method to wrap positions in a pandas dataframe
In [4]:
class Particles:
def __init__(
self,
rng=np.random.RandomState(),
space=PeriodicSpace(),
x=None, y=None,
step_length=0.5
):
self.rng = rng
self.space = space
self.x, self.y = x, y
self.step_length = step_length
self.steps_done = 0

def move(self):
self.x += self.step_length * self.rng.normal(size=self.x.shape)
self.y += self.step_length * self.rng.normal(size=self.y.shape)

self.x, self.y = self.space.normalize_positions(self.x, self.y)

self.steps_done += 1

def center_of_mass(self):
return self.x.mean(), self.y.mean()

def moment_of_inertia(self):
return self.x.var() + self.y.var()

def diagnostics(self):
com = self.center_of_mass()
mi = self.moment_of_inertia()
return pd.DataFrame(
{
"center_of_mass_x": com[0],
"center_of_mass_y": com[1],
"moment_of_inertia": mi
},
index=[self.steps_done, ],
)

def positions(self):
return pd.DataFrame(
{
"x": self.x,
"y": self.y
}
)


Random number generator¶

There is already an object with all we need: numpy.random.RandomState. So we'll use this for our RNG.

Drive the experiment¶

Parameters we want to be able to control:

• number of particles
• number of steps to run the model
• size of the spatial domain
• step size in the updates
In [5]:
def run_random_walk(
length_x=10,
length_y=20,
number_particles=100,
number_steps=100,
step_length=0.5
):
# initialize spatial domain
space = PeriodicSpace(
length_x=length_x,
length_y=length_y
)

# initialize random number generator
rng = np.random.RandomState()

# initialize group of particles and register spatial domain and rng
particles = Particles(
space=space,
rng=rng,
x=np.ones((number_particles, )) * length_x / 2.0,
y=np.ones((number_particles, )) * length_y / 2.0,
step_length=step_length
)

# track initial positions
initial_positions = particles.positions()

# initialize dataframe with diagnostics
diags = particles.diagnostics()

# run number_steps updates and accumulate diagnostics
for step in range(1, number_steps):
particles.move()

diags = diags.append(
particles.diagnostics(),
ignore_index=True
)

# track final positions
final_positions = particles.positions()

return diags, initial_positions, final_positions

In [6]:
(
diags,
initial_positions,
final_positions
) = run_random_walk(number_particles=100, number_steps=1000)


Data and visualization¶

In [7]:
display(initial_positions.head(5))
display(diags.head(1))

display(final_positions.head(5))
display(diags.tail(1))

x y
0 5.0 10.0
1 5.0 10.0
2 5.0 10.0
3 5.0 10.0
4 5.0 10.0
center_of_mass_x center_of_mass_y moment_of_inertia
0 5.0 10.0 0.0
x y
0 9.491506 18.354704
1 3.822715 13.159536
2 5.388253 0.309439
3 4.807517 12.536099
4 1.886733 10.880113
center_of_mass_x center_of_mass_y moment_of_inertia
999 5.015665 9.934202 39.863907
In [8]:
display(init)
display(diags)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-16352240deee> in <module>
----> 1 display(init)
2 display(diags)

NameError: name 'init' is not defined
In [ ]:
initial_positions.plot.scatter(x="x", y="y")
final_positions.plot.scatter(x="x", y="y");

In [ ]:
diags.plot();