heteromotility/hmsims module

Classes to generate simulated data sets for analysis, including 1. Random walk simulations 2. Biased random walk simulations 3. Levy and power flight simulations

class heteromotility.hmsims.RWSims(obj=100, length=100, bound_x=(-inf, inf), bound_y=(-inf, inf), speed_range=range(5, 10), bias_prop=0.5)

Bases: object

Simulates random walkers.

obj

int – number of walker objects to simulate.

length

int – number of time steps to simulate. i.e. length = total_time_points - 1

bound_x, bound_y

tuple, (min, max), optional – boundaries of X and Y dimensions for motion.

speed_range

iterable, optional – range of mean displacement sizes to simulate

bias_prop

float, optional – proportion of the time biased RWs move in the biased direction

rw_paths

dict – unbiased random walk simulation paths. keyed by mean displacement size, valued by dicts in cell_ids format

brw_paths

dict – biased random walk simulation paths. keyed by mean displacement size, valued by dicts in cell_ids format

random_walk(N, speed_mu=7.0, speed_sigma=None, origin=(0, 0), bound_x=(-inf, inf), bound_y=(-inf, inf))

Models a single random walk given an origin, length, and stepsize.

Parameters:
  • N (length of walk to model) –
  • speed_mu (float, optional.) – mean of Gaussian distribution speeds are sampled from
  • speed_sigma (float, optional.) – standard deviation of Gaussian distribution speeds are sampled from
  • origin (tuple.) – origin point, default = (0,0)
  • bound_y (bound_x,) – (min, max) allowed values for X and Y coordinates. Useful to confine model to a specified box.
Returns:

model_path
Return type:

list containing sequential tuples of XY coordinates

Notes

step_size = rate = c = sqrt(x_step^2 + y_step^2) = sqrt(2*step^2) x_step = y_step = rate/sqrt(2) (Random Walks in Biology, 1992, Howard C. Berg)

rwsim()

Generates simulated object paths using an unbiased random walk model

Parameters:
  • obj (integer specifying the number of object paths to be simulated) –
  • length (integer specifying the length of each path to be simulated) –
  • speed_range (tuple containing integers of step sizes to use when) –
  • models (building) –
  • bound_x (tuple, optional.) – (min, max) allowed values for X coordinates. Useful to confine model to a specified box.
  • bound_y (tuple, optional.) – (min, max) allowed values for Y coordinates. Useful to confine model to a specified box.
Returns:

rw_paths – keyed by speed, values with dicts in cell_ids format.
Return type:

dict
biased_random_walk(N, speed_mu=7.0, speed_sigma=None, bias_prop=0.5, origin=(0, 0), bound_x=(-inf, inf), bound_y=(-inf, inf))

Models a set of random walks biased to move in one direction.

Parameters:
  • N (int) – length of walk to be modeled
  • speed_mu (float) – mean displacement size to be modeled
  • bias_prop (float) – proportion of the time an object moves in the direction of bias
  • origin (array-like, len 2) – origin point for the start of the walk to be modeled
  • bound_y (bound_x,) – (min, max) allowed values for X and Y coordinates. Useful to confine model to a specified box.
Returns:

model_path – tuples, each containing an XY point.
Return type:

list

Notes

Randomly selects in which direction to bias movement. Object will move in the direction of bias a proportion of the time specified in bias_prop.

biased_rwsim()

Generates simulated object paths using a biased random walk model.

Parameters:
  • obj (int) – number of object paths to be simulated.
  • length (int) – length of each path to be simulated.
Returns:

brw_paths – keyed by mean displacement, valued with cell_ids format dicts.
Return type:

dict
class heteromotility.hmsims.LevySims(obj=100, length=100, bound_x=(-inf, inf), bound_y=(-inf, inf), p_scale=range(25, 30), alpha=2.0, levy=False, levy_scale=range(1, 5))

Bases: object

Simulated power and true ‘Levy’ flier walker models.

obj

int – number of walker objects to simulate.

length

int – number of time steps to simulate. i.e. length = total_time_points - 1

bound_x, bound_y

tuple, (min, max), optional – boundaries of X and Y dimensions for motion.

p_scale

iterable – scaling values for power fliers. 25 = 5 px step mean.

alpha

float, optional – alpha exponent for power distribution.

levy

boolean, optional – simulate Levy fliers.

levy_scale

iterable, optional – scaling values for Levy fliers.

power_paths

dict – keyed by scale factor, valued with dicts in cell_ids format.

levy_paths

dict – keyed by scale factor, valued with dicts in cell_ids format.

power_flier(N, x_min=0.1, alpha=2.0, scale_coeff=25, origin=(0, 0), bound_x=(-inf, inf), bound_y=(-inf, inf))

Models a power flier path.

Parameters:
  • N (int) – number of steps in the flight
  • x_min (int) – integer value, minimum step size
  • scale_coeff (float) – scaling coefficient for step sizes
  • alpha (float) – alpha exponent for power flier model, see Notes
  • origin (array-like) – initial site of walk, len 2
  • bound_y (bound_x,) – (min, max) allowed values for X and Y coordinates.
Returns:

model_path – tuples of XY coordinates
Return type:

list

Notes

Directions are chosen randomly, as in a random walk Step sizes however are chosen according to the distribution:

P(l_j) ~ l_j**(-a) where a is a parameter in the range 1 < u <= 3

Defaults to a = 2 as the optimal parameter for foraging behavior see : Viswanathan 1999, Nature

Obtaining Levy distributed variables :

To transform uniformly distributed random variable into another distribution, use inverse cum. dist. func. (CDF)

if F is CDF corresponding to probability density of variable f and u is a uniformly distributed random variable on [0,1]

x = F^-1(u) for pure power law distribution P(l) = l**-a F(x) = 1 - ( x / x_min ) where x_min is a lower bound on the random variable

F^-1(u) = x_min(1 - u)^(-1/a)

See: (Devroye 1986, Non-uniform random variable generation)

or

X = F^-1(U) = c / ( phi^-1(1-U/2) )**2 + mu

c : scaling parameter Phi(x) : CDF of Gaussian distribution mu : location

for a pure Levy distribution

see http://www.math.uah.edu/stat/special/Levy.html

power_sim()

Generates simulated object paths using a levy flight model with a power law distribution of displacements

Parameters:
  • obj (int) – number of object paths to be simulated
  • length (int) – length of each path to be simulated
  • scale_range (iterable) – scaling factors to simulate
  • bound_y (bound_x,) – (min, max) allowed values for X and Y coordinates.
  • alpha (float) – alpha exponent for power flier model, see Notes
Returns:

  • levy_paths (dict keyed by obj id, containing lists of tuples, each tuple)
  • specifying the XY position of a simulated object at a given timepoint – levy_paths = { 0: [ (x1,y1), (x2,y2)…(x_n, y_n) ] 1: [ (x1,y1), (x2,y2)…(x_n, y_n) ] … }
  • Distribution with scale factor = 25 has mean = 5
levy_flier(N, scale_coeff=1, origin=(0, 0), bound_x=(-inf, inf), bound_y=(-inf, inf))

Defines a levy flight path beginning at the origin, taking N steps, with step sizes pulled from a true Levy distribution

Parameters:
  • N (number of steps in the flight) –
  • step_min (integer value, minimum step size) –
  • origin (initial site of walk) –
  • = parameter for probability distribution of step sizes (u) –
Returns:

model_path
Return type:

a list containing tuples of XY coordinates

Notes

scale = 2, median disp = 4.4 scale = ~2.3 has a median displacement of 5 units

Directions are chosen randomly, as in a random walk Step sizes however are chosen according to the distribution:

P(l_j) ~ l_j**(-a) where a is a parameter in the range 1 < u <= 3

Defaults to a = 2 as the optimal parameter for foraging behavior see : Viswanathan 1999, Nature

Obtaining Levy distributed random variables:

To transform uniformly distributed random variable into another distribution, use inverse cum. dist. func. (CDF)

if F is CDF corresponding to probability density of variable f and u is a uniformly distributed random variable on [0,1]

x = F^-1(u) for pure power law distribution P(l) = l**-a F(x) = 1 - ( x / x_min ) where x_min is a lower bound on the random variable

F^-1(u) = x_min(1 - u)^(-1/a)

See: (Devroye 1986, Non-uniform random variable generation)

or

X = F^-1(U) = c / ( phi^-1(1-U/2) )**2 + mu

c : scaling parameter Phi(x) : CDF of Gaussian distribution mu : location

for a pure Levy distribution

References

[1] http://www.math.uah.edu/stat/special/Levy.html

levy_sim()

Generates simulated object paths using a levy flight model with a true Levy distribution of displacements

Parameters:
  • obj (integer specifying the number of object paths to be simulated) –
  • length (integer specifying the length of each path to be simulated) –
Returns:

levy_paths – keyed by scale factor, values are cell_ids format dicts
Return type:

dict
class heteromotility.hmsims.fBmsims(obj=100, length=100, H=0.8, bound_x=(-inf, inf), bound_y=(-inf, inf), scale_range=range(1, 5))

Bases: object

Simulates fractal Brownian motion walkers.

obj

int – number of objects to simulate.

length

int – length of simulated tracks.

H

float – Hurst coefficient for fractal Brownian process

bound_x, bound_y

tuple, optional. – (min, max) allowed values for X and Y coordinates.

scale_range

iterable – scaling values for fBm walkers. 0.277 = 5 px step mean.

fbm_paths

dict – keyed by scale factor, values are cell_ids format dicts.

fbm_mat(gridN: int = 100, T: int = 100, H: float = 0.99) → numpy.ndarray

Generates a Cholesky decomposition of a covariance matrix

Gamma = (R(i,j), i, j) for i,j = 0…n

where n is a grid size parameter

Parameters:
  • gridN (int) – size of grid to generate for the covariance matrix
  • T (int) – length of displacement series to generate
  • H (float, (0,1)) – hurst parameter (0,1)
Returns:

X – series of displacements following fBm for the given Hurst index
Return type:

np.ndarray
fbm_series(gridN: int = 100, T: int = 100, H: float = 0.99) → numpy.ndarray

Generates a series of displacements simulating a fractal Brownian process with a given Hurst parameter of length T using the Cholesky decomposition method [1].

Parameters:
  • gridN (int) – size of grid to generate for the covariance matrix
  • T (int) – length of displacement series to generate
  • H (float, (0,1)) – hurst parameter (0,1)
Returns:

u – series of a fBm process.
Return type:

np.ndarray

Notes

  1. Generate a matrix Gamma = [ R(t_i,t_j), for i,j = 0…n ]
  2. take the cholesky decomposition of Gamma
  3. Find Sigma = cholesky_Gamme**0.5
  4. Generate a vector v or length n with Gaussian distributed values
  5. vector u = Sigma * v represents a fBm series

References

  1. Dieker, Simulation of fractional Brownian motion, 2004
fbm_walker(N, scale=1, H=0.99, origin=(0, 0), bound_x=(-inf, inf), bound_y=(-inf, inf))
fbm_sim()

Generates simulated object paths using a fractal Brownian motion model

Parameters:
  • obj (integer specifying the number of object paths to be simulated) –
  • length (integer specifying the length of each path to be simulated) –
  • H (Hurst coefficient for generation of fBm displacement series) –
  • scale_range (range of scales to use) –
Returns:

fbm_paths
Return type:

dict keyed by scale, values are cell_ids format dicts.

Notes

scale = 0.277 is equivalent to mean displacement = 5 units

class heteromotility.hmsims.MultiSim(obj=1000, tau=20)

Bases: object

Simulates objects that switch motion model state.

obj

int – number of object paths to be simulated.

tau

int – length of each state to be simulated.

pwr2fbm

dict. – power walkers switching to fractal Brownian motion after tau. keyed by object id, values are lists of coordinate tuples.

pwr2urw

dict. – power walkers switching to unbiased random walks after tau. keyed by object id, values are lists of coordinate tuples.

fbm2urw

dict. – fractal Brownian motion switching to unbiased random walks after tau. keyed by object id, values are lists of coordinate tuples.

join_paths(ids1, ids2)
power_fbm(obj=1000, tau=20)
power_urw(obj=1000, tau=20)
fbm_rw(obj=1000, tau=20)