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();