Walking on lines, lattices and volumes

Random walk on line.

Random walk on lattice.

Random walk on volume (3+ environment).

Simple random walks simulation in one dimention

In [2]:
import numpy as np
import random

def random_walk(n=100):
    ys = []
    current = 0
    for _ in range(n):
        if random.random() > 0.5:
            current += 1
        else:
            current -= 1
        ys.append(current)
    return ys

Load plotting tools

In [3]:
%matplotlib inline

from matplotlib import pylab

Plot like random walks look like

Each walk has a different color

Histogram of visited places.

In [4]:
def print_random_walks():
    all_ys = []
    pylab.figure(figsize=(10, 10))
    for _ in range(50):
        ys = random_walk(1000)
        pylab.plot(ys)
        all_ys += ys

    pylab.figure(figsize=(10, 10))
    pylab.hist(all_ys, normed=1, bins=10)
    
print_random_walks()

Some visualization questions

How to represent it in 2D?

How to represnt it in 3D?

How to represent a random walk on line and using hypergraphs in the same way?

1 Dimensional random walk on a hypergraph

I would think about it in this way:

[o(o{o][o)o}..

It's very easy to model random walks on Markov Chains, however - with indefinite number of states (possible in random walk on line) can get tricky.

How to approximate random walk on line with Markov Chain?

In [5]:
import pykov


def prepare_random_walk_transition_matrix(N):
    matrix = np.zeros((N, N))

    for i in range(1, N - 1):
        matrix[i][i + 1] = 0.5
        matrix[i][i - 1] = 0.5

    matrix[0][1] = 1
    matrix[-1][-2] = 1
    return matrix


def transition_matrix_to_pykov_chain(matrix):
    chain = pykov.Chain()
    
    for i, row in enumerate(matrix):
        for j, column in enumerate(row):
            chain[(i, j)] = column
    return chain


def markov_chain_random_walk(n, matrix_size=1000):
    matrix = prepare_random_walk_transition_matrix(matrix_size)
    chain = transition_matrix_to_pykov_chain(matrix)
    return chain.steady().items()

results = markov_chain_random_walk(100)
In [6]:
xs, ys = zip(*results)
pylab.plot(xs, ys)
Out[6]:
[<matplotlib.lines.Line2D at 0x7f0b65cc5b38>]

We see here that steady states aren't the good way of approaching this problem

Random walks usually start somewhere and that should matter.

This doesn't look like the thing I wanted to see. We need a different approach.

In [7]:
from hypergraph.markov_diffusion import (create_markov_matrix_model_nodes,
                                         create_markov_matrix_model_hyper_edges)

from matplotlib import pyplot as plt

# todo - move to some library
def plot_transition_matrix(matrix):
    conf_arr = matrix

    norm_conf = []
    for i in conf_arr:
        a = 0
        tmp_arr = []
        a = sum(i, 0)
        for j in i:
            tmp_arr.append(float(j)/float(a))
        norm_conf.append(tmp_arr)


    fig = plt.figure(figsize=(12, 10))
    plt.clf()
    ax = fig.add_subplot(111)
    ax.set_aspect(1)
    res = ax.imshow(np.array(norm_conf), cmap=plt.cm.gray, 
                    interpolation='nearest')

    width = len(conf_arr)
    height = len(conf_arr[0])

    for x in range(width):
        for y in range(height):
            ax.annotate("%.2f" % (conf_arr[x][y]), xy=(y, x), 
                        horizontalalignment='center',
                        verticalalignment='center')

    cb = fig.colorbar(res)
    index = range(1, len(matrix) + 1)
    plt.xticks(range(width), index)
    plt.yticks(range(height), index)
    
In [8]:
matrix_size = 10
markov_matrix = prepare_random_walk_transition_matrix(matrix_size)
plot_transition_matrix(markov_matrix)
In [9]:
def plot_binned(all_ys, N):
    binned_xs = range(0, len(all_ys), N)
    norm = 1 / sum(all_ys)
    binned_ys = [sum(all_ys[x:x + N]) * norm for x in range(0, len(all_ys), N)]
    plt.bar(binned_xs, binned_ys, width=N)
    return binned_xs, binned_ys
In [10]:
matrix_size = 200
markov_matrix = prepare_random_walk_transition_matrix(matrix_size)
In [11]:
from hypergraph.diffusion_engine import DiffusionEngine
In [12]:
def simulate_walk(markov_matrix, t_max=1000, starting_point=None):
    
    engine = DiffusionEngine(markov_matrix, t_per_walker=t_max)
    
    if starting_point is None:
        starting_point = int(len(markov_matrix) / 2)
    
    frequencies, states = engine.simulate_plain(t_max, starting_point)

    frequencies = [(node, frequency) for node, frequency in frequencies]
    connected_elements = [freq[0] + 1 for freq in frequencies]


    missing_nodes = set(range(1, len(markov_matrix) + 1)) - set(connected_elements)
    for missing_node in missing_nodes:
        frequencies.append((missing_node, 0))
                
    frequencies.sort(key=lambda x: x[0])
    xs, ys = zip(*frequencies)

    ys = np.array(ys, dtype='float')
    ys /= sum(ys)

    return xs, ys
In [14]:
def illustrate_walk(markov_matrix, N=2, starting_point=None):
    if starting_point is None:
        starting_point = int(len(markov_matrix) / 2.0)
    xs = None
    all_ys = None
    for x in range(N):
        xs, ys = simulate_walk(markov_matrix, starting_point)
        if all_ys is None:
            all_ys = ys
        else:
            all_ys += ys
        print('%s/%s' % (x + 1, N))

    plt.figure(figsize=(10, 10))
    plt.bar(xs, all_ys, color='crimson', label='Simulated')
    
    return xs, all_ys
In [16]:
xs, random_walk_ys = illustrate_walk(markov_matrix, 2)
1/2
2/2
In [17]:
plot_binned(random_walk_ys, 5)
Out[17]:
(range(0, 200, 5),
 [0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.034999999999999983,
  0.23999999999999988,
  0.15499999999999994,
  0.074999999999999969,
  0.099999999999999964,
  0.069999999999999965,
  0.089999999999999955,
  0.14499999999999993,
  0.089999999999999955,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0])

This actually looks better

Let's think now a little bit about line hypergraphs

In [19]:
from hypergraph.hypergraph_models import HyperGraph

def create_line_hypergraph(N, cardinality):
    HG = HyperGraph()
    nodes = range(N)
    hyper_edges = [range(i, i + cardinality) for i in range(N + 1 - cardinality)]
    HG.add_nodes_from(nodes)
    HG.add_edges_from(hyper_edges)
    return HG

HG = create_line_hypergraph(200, 3)
tiny_HG = create_line_hypergraph(20, 3)
In [20]:
# similar to line
HG = create_line_hypergraph(200, 2)
tiny_HG = create_line_hypergraph(20, 2)
hypergraph_line_walk_matrix = create_markov_matrix_model_nodes(HG)
tiny_hypergraph_line_walk_matrix = create_markov_matrix_model_nodes(tiny_HG)
plot_transition_matrix(tiny_hypergraph_line_walk_matrix)

xs, line_edges_ys = illustrate_walk(hypergraph_line_walk_matrix)
1/2
2/2
In [29]:
plot_binned(line_edges_ys, N=2)
Out[29]:
(range(0, 200, 2),
 [0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0032000000000000002,
  0.0080000000000000002,
  0.0058000000000000005,
  0.017999999999999999,
  0.046800000000000001,
  0.064400000000000013,
  0.10680000000000002,
  0.16300000000000003,
  0.20879999999999999,
  0.15059999999999998,
  0.11219999999999999,
  0.060000000000000019,
  0.029999999999999999,
  0.0091999999999999998,
  0.0035999999999999999,
  0.0042000000000000006,
  0.0044000000000000003,
  0.001,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0])

Node model first

In [22]:
HG = create_line_hypergraph(200, 3)
tiny_HG = create_line_hypergraph(10, 3)
hypergraph_line_walk_matrix = create_markov_matrix_model_nodes(HG)
tiny_hypergraph_line_walk_matrix = create_markov_matrix_model_nodes(tiny_HG)
plot_transition_matrix(tiny_hypergraph_line_walk_matrix)
In [34]:
xs, line_edges_ys = illustrate_walk(hypergraph_line_walk_matrix)
1/50
2/50
3/50
4/50
5/50
6/50
7/50
8/50
9/50
10/50
11/50
12/50
13/50
14/50
15/50
16/50
17/50
18/50
19/50
20/50
21/50
22/50
23/50
24/50
25/50
26/50
27/50
28/50
29/50
30/50
31/50
32/50
33/50
34/50
35/50
36/50
37/50
38/50
39/50
40/50
41/50
42/50
43/50
44/50
45/50
46/50
47/50
48/50
49/50
50/50
In [36]:
plot_binned(line_edges_ys, N=10)
Out[36]:
(range(0, 200, 10),
 [0.00038000000000000029,
  0.001720000000000001,
  0.0042400000000000024,
  0.014280000000000006,
  0.020560000000000012,
  0.041100000000000025,
  0.066000000000000059,
  0.10048000000000007,
  0.12228000000000011,
  0.17572000000000018,
  0.17098000000000019,
  0.10466000000000009,
  0.065160000000000051,
  0.051240000000000036,
  0.028760000000000018,
  0.015420000000000008,
  0.010320000000000006,
  0.0053800000000000028,
  0.0013200000000000009,
  0.0])

Hyper edge model next

In [21]:
hypergraph_line_walk_edges_matrix = create_markov_matrix_model_hyper_edges(HG)
tiny_hypergraph_line_walk_edges_matrix = create_markov_matrix_model_hyper_edges(tiny_HG)
plot_transition_matrix(tiny_hypergraph_line_walk_edges_matrix)
In [37]:
xs, line_hyper_edges_ys = illustrate_walk(hypergraph_line_walk_edges_matrix)
1/50
2/50
3/50
4/50
5/50
6/50
7/50
8/50
9/50
10/50
11/50
12/50
13/50
14/50
15/50
16/50
17/50
18/50
19/50
20/50
21/50
22/50
23/50
24/50
25/50
26/50
27/50
28/50
29/50
30/50
31/50
32/50
33/50
34/50
35/50
36/50
37/50
38/50
39/50
40/50
41/50
42/50
43/50
44/50
45/50
46/50
47/50
48/50
49/50
50/50
In [39]:
plot_binned(line_hyper_edges_ys, N=10)
Out[39]:
(range(0, 198, 10),
 [0.00021999999999999998,
  0.0007000000000000001,
  0.0013600000000000001,
  0.0037000000000000002,
  0.010120000000000001,
  0.019180000000000003,
  0.044900000000000002,
  0.060000000000000012,
  0.090980000000000033,
  0.15752000000000008,
  0.18332000000000009,
  0.14888000000000007,
  0.11822000000000005,
  0.069660000000000027,
  0.029320000000000006,
  0.019279999999999999,
  0.0184,
  0.015560000000000003,
  0.0068999999999999999,
  0.0017800000000000003])

What about the hypergraphs of cardinality 4 (even) and 5 (odd)?

In [85]:
HG = create_line_hypergraph(200, 5)
tiny_HG = create_line_hypergraph(20, 5)
In [86]:
hypergraph_line_walk_matrix = create_markov_matrix_model_nodes(HG)
tiny_hypergraph_line_walk_matrix = create_markov_matrix_model_nodes(tiny_HG)
plot_transition_matrix(tiny_hypergraph_line_walk_matrix)
In [87]:
xs, line_edges_ys = illustrate_walk(hypergraph_line_walk_matrix)
1/50
2/50
3/50
4/50
5/50
6/50
7/50
8/50
9/50
10/50
11/50
12/50
13/50
14/50
15/50
16/50
17/50
18/50
19/50
20/50
21/50
22/50
23/50
24/50
25/50
26/50
27/50
28/50
29/50
30/50
31/50
32/50
33/50
34/50
35/50
36/50
37/50
38/50
39/50
40/50
41/50
42/50
43/50
44/50
45/50
46/50
47/50
48/50
49/50
50/50
In [88]:
xs, ys = plot_binned(line_edges_ys, N=10)
In [89]:
import numpy
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Define some test data which is close to Gaussian
data = line_hyper_edges_ys
bin_edges = ys


# Define model function to be used to fit to the data above:
def gauss(x, *p):
    A, mu, sigma = p
    return A*numpy.exp(-(x-mu)**2/(2.*sigma**2))

# p0 is the initial guess for the fitting coefficients (A, mu and sigma above)
p0 = [1., 100, 100.]

coeff, var_matrix = curve_fit(gauss, xs, ys, p0=p0)
# Get the fitted curve
hist_fit = gauss(xs, *coeff)

plt.plot(xs, ys, 'o', label='Test data')
plt.plot(xs, hist_fit, label='Fitted data')

# Finally, lets get the fitting parameters, i.e. the mean and standard deviation:
print('Fitted mean = ', coeff[1])
print('Fitted standard deviation = ', coeff[2])

plt.show()
Fitted mean =  102.79737442
Fitted standard deviation =  44.1545677506
In [100]:
import numpy as np


def square_grid_matrix(n):
    adjacency_matrix = np.zeros((n * n, n * n))
    for i in range(n ** 2):
        for j in range(n ** 2):
            if j == i + 1 and j // n == i // n:
                adjacency_matrix[i, j] = 1
            elif j == i - 1 and j // n == i // n:
                adjacency_matrix[i, j] = 1
            elif j == i + n:
                adjacency_matrix[i, j] = 1
            elif j == i - n:
                adjacency_matrix[i, j] = 1
        adjacency_matrix[i] /= sum(adjacency_matrix[i])
    return adjacency_matrix
In [101]:
small_square_matrix = square_grid_matrix(4)
plot_transition_matrix(small_square_matrix)
square_matrix = square_grid_matrix(20)
In [121]:
xs, line_edges_ys = illustrate_walk(square_matrix, starting_point=210)
1/50
2/50
3/50
4/50
5/50
6/50
7/50
8/50
9/50
10/50
11/50
12/50
13/50
14/50
15/50
16/50
17/50
18/50
19/50
20/50
21/50
22/50
23/50
24/50
25/50
26/50
27/50
28/50
29/50
30/50
31/50
32/50
33/50
34/50
35/50
36/50
37/50
38/50
39/50
40/50
41/50
42/50
43/50
44/50
45/50
46/50
47/50
48/50
49/50
50/50
In [122]:
def one_to_two_d(xs, ys):
    z_data = ys

    print(np.sqrt(len(xs)))
    x_data, y_data = np.meshgrid(np.arange(np.sqrt(len(xs))),
                                 np.arange(np.sqrt(len(xs))))
    x_data = x_data.flatten()
    y_data = y_data.flatten()
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')
    ax.bar3d( x_data,
          y_data,
          np.zeros(len(z_data)),
          1, 1, z_data )
    
    plt.show()
In [124]:
one_to_two_d(xs, line_edges_ys)
20.0