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")
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
# 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)
Check if a connected component contains cycles
# Check if cycles
list(nx.simple_cycles(nx.subgraph(g,ccs[1])))
Yep! One between 10 and 5
# 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)