Import packages

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
from math import *

And this is just to turn off the warnings for the notebook (since some of the plot functions will spit out various warnings):

import warnings
warnings.filterwarnings("ignore")

Set up our CA

Some basic parameters for our model

g = nx.DiGraph() # Make an empty graph that will be the phase space
L = 4 #9 #10     # Grid size (start with something small)

Choose a rule set for our CA

We'll start with Rule 50 and then try a few others

# Rule 50
neighrule = {
        (0,0,0):0,
        (0,0,1):1,
        (0,1,0):0,
        (0,1,1):0,
        (1,0,0):1,
        (1,0,1):1,
        (1,1,0):0,
        (1,1,1):0
        }

# Rule 13
#neighrule = {
#        (0,0,0):1,
#        (0,0,1):0,
#        (0,1,0):1,
#        (0,1,1):1,
#        (1,0,0):0,
#        (1,0,1):0,
#        (1,1,0):0,
#        (1,1,1):0
#        }

# Rule 30
#neighrule = {
#        (0,0,0):0,
#        (0,0,1):1,
#        (0,1,0):1,
#        (0,1,1):1,
#        (1,0,0):1,
#        (1,0,1):0,
#        (1,1,0):0,
#        (1,1,1):0
#        }

Useful Functions: Functions for running the model

We'll start with two functions to convert between decimal integers and binary model configurations

# Takes a configuration and returns the corresponding integer
def config2int(config):
    return int(''.join(map(str, config)),2) #maps the config->strings, joins them, and then converts to int from binary

# Takes an integer and converts it to a configuration (list of cell states)
def int2config(x):
    return [1 if x & 2**i > 0 else 0 for i in range(L - 1, -1, -1)]

Next, a function to run the model--take the current configuration and map to next timestep

# Function to update the CA one timestep
def update(config):
    nextconfig = [0]*L
    for x in range(L):
        nextconfig[x] = neighrule[(config[(x - 1) % L],config[x],config[(x + 1) % L])]
    return nextconfig

Map out the phase space!

# Go through every possible configuration and add an edge linking it to where it goes next
for x in range(2**L):
    g.add_edge(x, config2int(update(int2config(x))))

Plot the phase space network

# Plot each connected component of the phase space
ccs = [cc for cc in nx.connected_components(g.to_undirected())]
n = len(ccs)
w = ceil(sqrt(n))
h = ceil(n / w)
for i in range(n):
    plt.subplot(h, w, i + 1)
    nx.draw_networkx(nx.subgraph(g, ccs[i]), with_labels = True)
plt.show()

Highlight an attracting sub-component

# Plot each component with the attracting subcomponent highlighted
subg = nx.subgraph(g, ccs[1])
attr = set().union(*nx.attracting_components(subg))

pos=nx.spring_layout(subg) # positions for all nodes
nx.draw_networkx_nodes(subg, pos, nodelist = attr, node_color='r', node_size=500, alpha=0.8)
nx.draw_networkx_nodes(subg,pos, nodelist = (set(subg.nodes()) - attr), node_color='#00b4e9', node_size=500, alpha=0.8)
nx.draw_networkx_edges(subg,pos,width=2.0)
nx.draw_networkx_labels(subg,pos)

{1: Text(0.5169937444849606, -0.12428797082992608, '1'),
 10: Text(0.7038321131882188, 0.4015101274979579, '10'),
 2: Text(-0.12066314943124436, -0.13618536138051654, '2'),
 5: Text(0.2559792897492209, 0.28220023311008946, '5'),
 4: Text(0.5243166229707061, 0.7789484885741887, '4'),
 7: Text(-1.0, -0.36407846672195343, '7'),
 8: Text(-0.4324368257049425, 0.008427034707937249, '8'),
 11: Text(-0.3119194206851483, 0.5080032501560819, '11'),
 13: Text(-0.3457099471580606, -0.5649093511375037, '13'),
 14: Text(0.20960757258628734, -0.7896279839763553, '14')}

Check if a connected component contains cycles

# Check if cycles
list(nx.simple_cycles(nx.subgraph(g,ccs[1])))

[[10, 5]]

Yep! One between 10 and 5

Simulate the model to observe the dynamics we've mapped

# Plot the CA behavior starting at an interesting IC from the above
# Run the model for a few steps and plot
steps = 20
output = np.zeros([steps,L])
output[0,:] = int2config(11)
for i in range(1,steps):
    output[i,:] = update(output[i-1,:])
plt.cla()
plt.imshow(output)

<matplotlib.image.AxesImage at 0xa14fba510>