Vectorized N-Dimensional Random Walk in Cartesian Coordinates

Question

I have written a random-walk routine that I hope to build upon in the future. Before doing so, I was hoping to get some critique.

I believe the implementation is correct. I noticed that many other implementations found online use for-loops or modules I am unfamiliar with; my goal was to vectorize the walk in n-dimensional space with the optional use of boundary conditions.

The main idea is to generate an n-dimensional array of random numbers according to the desired distribution. As of now, only the 'normal' distribution is implemented. If a threshold is not set, then the average of the distribution is used as a threshold. Numbers greater than this threshold are taken in the positive direction, whereas numbers less than this threshold are taken in the negative direction. Should the number exactly equal this threshold, then no step is taken. The initial steps array (called base initially consists of all zeros; the indices corresponding to the positive and negative steps are used to mask this array with the respective step vectors (magnitude and direction).

If edge_type is not None, then the boundary conditions corresponding to edge_type will be used. If edge_type='bounded', then the steps at the boundaries will be zero. If edge_type='pacman', then the steps at the boundaries will be of magnitude max_edge - min_edge and taken to be in the direction away from the respective edge.

import numpy as np
import matplotlib.pyplot as plt

class IndicialDistributions():

    """
    Steps are taken in a positive or negative direction according
    to a random number distribution. The methods of this class
    return the indices for positive and negative steps given
    a distribution type.

    As of now, only the 'normal' distribution type is implemented.
    """

    def __init__(self, nshape):
        """
        nshape              :   type <int / tuple / array>
        """
        self.nshape = nshape

    def get_normal_indices(self, mu, sigma, threshold=None):
        """
        mu                  :   type <int / float>
        sigma               :   type <int / float>
        threshold           :   type <int / float> or None
        """
        if threshold is None:
            threshold = mu
        random_values = np.random.normal(mu, sigma, size=self.nshape)
        pos = (random_values > threshold)
        neg = (random_values < threshold)
        return pos, neg

    def get_binomial_indices(self, p_success, threshold=None):
        """
        p_success           :   type <float>
        threshold           :   type <int / float> or None
        """
        raise ValueError("not yet implemented")

    @property
    def indicial_function_mapping(self):
        res = {}
        res['normal'] = self.get_normal_indices
        res['binomial'] = self.get_binomial_indices
        return res

    def dispatch_indices(self, distribution, **kwargs):
        """
        distribution        :   type <str>
        """
        available_keys = list(self.indicial_function_mapping.keys())
        if distribution not in available_keys:
            raise ValueError("unknown distribution: {}; available distributions: {}".format(distribution, available_keys))
        f = self.indicial_function_mapping[distribution]
        pos, neg = f(**kwargs)
        return pos, neg

class BoundaryConditions():

    """
    The methods of this class account for steps taken at the
    n-dimensional edges.

    As of now, 'bounded' and 'pacman' edges are implemented.
    """

    def __init__(self, steps, movements, max_edge, min_edge):
        """
        steps               :   type <array>
        movements           :   type <array>
        max_edge            :   type <int / float>
        min_edge            :   type <int / float>
        """
        self.steps = steps
        self.movements = movements
        self.max_edge = max_edge
        self.min_edge = min_edge

    def get_maximum_edge_indices(self):
        indices = (self.movements >= self.max_edge)
        return indices

    def get_minimum_edge_indices(self):
        indices = (self.movements <= self.min_edge)
        return indices

    def apply_maximum_bounded_edge(self):
        indices = self.get_maximum_edge_indices()
        self.steps[indices] = 0

    def apply_minimum_bounded_edge(self):
        indices = self.get_minimum_edge_indices()
        self.steps[indices] = 0

    def apply_pacman_edges(self):
        max_indices = self.get_maximum_edge_indices()
        min_indices = self.get_minimum_edge_indices()
        self.steps[max_indices] = self.min_edge - self.max_edge
        self.steps[min_indices] = self.max_edge - self.min_edge

    def apply_to_dimension(self, edge_type):
        """
        edge_type           :   type <str>
        """
        if edge_type is not None:
            if edge_type == 'bounded':
                self.apply_maximum_bounded_edge()
                self.apply_minimum_bounded_edge()
            elif edge_type == 'pacman':
                self.apply_pacman_edges()
            else:
                raise ValueError("unknown edge_type: {}; available edge_type = 'bounded', 'pacman', or None".format(edge_type))

class CartesianRandomWalker():

    """
    This class has methods to perform a random walk in n-dimensional space
    with the optional use of boundary conditions.
    """

    def __init__(self, initial_position, nsteps, edge_type=None, max_edges=(), min_edges=()):
        """
        initial_position    :   type <tuple / list / array>
        nsteps              :   type <int>
        edge_type           :   type <str>
        max_edges           :   type <tuple / list / array>
        min_edges           :   type <tuple / list / array>
        """
        self.initial_position = initial_position
        self.nsteps = nsteps
        self.edge_type = edge_type
        self.max_edges = max_edges
        self.min_edges = min_edges
        self.ndim = len(initial_position)
        self.nshape = (self.ndim, nsteps)
        self.base = np.zeros(self.nshape).astype(int)
        # self.boundary_crossings = 0
        self.current_position = np.array(initial_position)

    def __repr__(self):
        if np.all(self.base == 0):
            string = 'Initial Position:\t{}'.format(self.initial_position)
        else:
            string = 'Initial Position:\t{}\nNumber of Steps:\t{}\nCurrent Position:\t{}'.format(self.initial_position, self.nsteps, self.current_position)
        return string

    @property
    def position(self):
        return tuple(self.current_position)

    @property
    def movement(self):
        return np.cumsum(self.base, axis=1)

    def initialize_steps(self, distribution, delta_steps, **kwargs):
        """
        distribution        :   type <str>
        delta_steps         :   type <tuple / list / array>
        """
        pos, neg = IndicialDistributions(self.nshape).dispatch_indices(distribution, **kwargs)
        self.base[pos] = delta_steps[0]
        self.base[neg] = delta_steps[1]

    def apply_boundary_conditions(self):
        if self.edge_type is not None:
            for idx in range(self.ndim):
                max_edge, min_edge = self.max_edges[idx], self.min_edges[idx]
                steps = self.base[idx]
                movements = self.movement[idx] + self.initial_position[idx]
                BC = BoundaryConditions(steps, movements, max_edge, min_edge)
                BC.apply_to_dimension(self.edge_type)
                self.base[idx, :] = BC.steps

    def update_positions(self, distribution, delta_steps=(1, -1), **kwargs):
        """
        distribution        :   type <str>
        delta_steps         :   type <tuple / list / array>
        """
        self.initialize_steps(distribution, delta_steps, **kwargs)
        self.apply_boundary_conditions()
        delta_position = self.movement[:, -1]
        self.current_position += delta_position

    def view(self, ticksize=7, labelsize=8, titlesize=10):
        """
        ticksize            :   type <int>
        labelsize           :   type <int>
        titlesize           :   type <int>
        """
        if self.ndim == 1:
            raise ValueError("not yet implemented")
        elif self.ndim == 2:
            fig, ax = plt.subplots()
            x_movement = self.movement[0] + self.initial_position[0]
            y_movement = self.movement[1] + self.initial_position[1]
            ax.scatter(*self.initial_position, color='k', label='Initial Position', marker='x', s=100)
            ax.scatter(*self.current_position, color='k', label='Current Position', marker='.', s=100)
            ax.plot(x_movement, y_movement, color='r', alpha=1/3, label='Random Walk')
            ax.grid(color='k', linestyle=':', alpha=0.3)
            ax.set_xlabel('X', fontsize=labelsize)
            ax.set_ylabel('Y', fontsize=labelsize)
            ax.tick_params(axis='both', labelsize=ticksize)
            if self.edge_type is None:
                title = r'${}-$D Random Walk in Cartesian Space'.format(self.ndim)
            elif self.edge_type in ('bounded', 'pacman'):
                title = '${}-$D Random Walk in Cartesian Space\nvia {} Boundary Conditions'.format(self.ndim, self.edge_type.title())
            ax.set_title(title, fontsize=titlesize)
            plt.subplots_adjust(bottom=0.2)
            fig.legend(loc='lower center', mode='expand', fancybox=True, ncol=3, fontsize=labelsize)
            plt.show()
            plt.close(fig)
        elif self.ndim == 3:
            raise ValueError("not yet implemented")
        else:
            raise ValueError("invalid ndim: {}; can only view 1 <= ndim <= 3".format(self.ndim))

As of now, only the 2-dimensional case is viewable. I can implement something similar for the 1-dimensional and 3-dimensional cases, but I am more concerned with the methods rather than the graph (for now). That said, one can run this algorithm in 10-dimensional space without the graph.

Here is an example call:

np.random.seed(327) ## reproduce random results

## initial position
# pos = (50, 50, 50, 50, 50, 50, 50, 50, 50, 50) ## 10-D
pos = (50, 50) ## 2-D
## number of steps to travel
nsteps = 100

## random number distribution
## average = 50, spread=10
## positive step if random number > 50
## negative step if random number < 50
## no step if random number = 0
distribution = 'normal'
kwargs = {'mu' : 50, 'sigma' : 10} # 'threshold' : 50}

## boundary conditions
max_edges = np.array([60 for element in pos])
min_edges = np.array([40 for element in pos])
edge_type = None
# edge_type = 'pacman'
# edge_type = 'bounded'

RW = CartesianRandomWalker(pos, nsteps, edge_type, max_edges, min_edges)
RW.update_positions(distribution, **kwargs)
print(RW)
RW.view()

Here is an example output from the 2-D case:

Initial Position:   (50, 50)
Number of Steps:    100
Current Position:   [36 58]

And here is an output from the 10-D case:

Initial Position:   (50, 50, 50, 50, 50, 50, 50, 50, 50, 50)
Number of Steps:    100
Current Position:   [36 58 52 58 38 58 42 78 28 48]

Answer 1

As a first comment, instead of specifying your types in the comments of each method, you can use the typing module for Python 3.5 and above https://docs.python.org/3/library/typing.html Also, I think the comments on variables should also describe what they represent.

Then, if I were you, I wouldn't try creating a very generic code, it usually leads to unnecessary over-engineering and more complex code to maintain. If you need more things, then you'll adapt your code to it as soon as you need. You usually first make it work with your new needs anyhow and then try to refactor it to make it better.

However, I'll try to give my point of view of how I'd change your code so it can be adapted more easily to what you said about multiple distributions.

IndicalDistributions

The get_normal_indices and the get_binomial_indices functions don't really belong to this class. IndicalDistributions should not know about all the possible distributions there can be. And you should definitely not store a dictionary with all of them in here, it leads to very messy and hard to maintain code.

• I think the easiest way to fix this, is to have separate implementations for each kind of distribution and choose the class when needed in your code, so you would end up removing IndicalDistributions and having something like NormalIndicalDistribution and BinomialIndicalDistribution classes separately. When you'll have to create BinomialIndicalDistribution, you'll see what are the common parts of NormalIndicalDistribution and BinomialIndicalDistribution, and create some kind of abstraction to do the common stuff (maybe strategy or template method pattern).

If you want to specify strings to choose the correct distribution class like you do in dispatch_indices, you can just have a function create_indical_distribution(distribution, **kwargs) that's just a bunch of ifs that return the correct object constructed. This is usually called a factory method.

But again, for now just create one NormalIndicalDistribution class and then you'll see what happens when you need the binomial one.

BoundaryConditions

You said nothing about having multiple edge types apart from these 2, so I would not touch this and if I would, the approach would be similar to the class before.

CartesianRandomWalker

self.base is a reserved name for the Python language, name it something else, it can give you very weird bugs.

In view, I think it'd be cleaner if you'd do something like

if self.ndim not in (1, 2, 3):
  raise ValueError(f"invalid ndim: {self.ndim}; can only view 1 <= ndim <= 3")
if self.ndim in (1, 2):
  raise ValueError("not yet implemented")

# Your implementation goes here 

However, I think the implementation for each dimension should go to a separate class and maybe have some common utilities for each dimension, it'd be much easier to maintain. Similar approach to what I explained before.

Usually, instantiating classes in the middle of other classes is not a very good idea, the objects should be passed in __init__ (I mean IndicalDistributions and BoundaryConditions). They should be constructed before and passed to CartesianRandomWalker, this will let you use more kinds of indical distributions and edge types, since you just pass whichever you want in the construction and that's it.

• For IndicalDistributions, it should be no problem, just remove the distribution string and pass the correct object when constructing.

• However, the BoundaryConditions object depends on parameters that you can only know in apply_boundary_conditions. This is a bit tricky. Honestly, since you didn't say anything about having more than these two edge types, I'd leave as is, otherwise, look up Builder Pattern.