Ggthjy

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 Spacenvia {} 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]

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.