#!/usr/bin/env python
# coding: utf-8

# > This is one of the 100 recipes of the [IPython Cookbook](http://ipython-books.github.io/), the definitive guide to high-performance scientific computing and data science in Python.
# 

# # 12.2. Simulating an elementary cellular automaton

# 1. We import NumPy and matplotlib.

# In[ ]:


import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')


# 2. We will use the following vector to obtain numbers written in binary representation.

# In[ ]:


u = np.array([[4], [2], [1]])


# 3. We write a function that performs one iteration on the grid, updating all cells at once according to the given rule in binary representation. The first step consists in stacking circularly shifted versions of the grid to get the LCR triplets of each cell (`y`). Then, we convert these triplets in 3-bit numbers (`z`). Finally, we compute the next state of every cell using the specified rule.

# In[ ]:


def step(x, rule_binary):
    """Compute a single stet of an elementary cellular
    automaton."""
    # The columns contains the L, C, R values
    # of all cells.
    y = np.vstack((np.roll(x, 1), x,
                   np.roll(x, -1))).astype(np.int8)
    # We get the LCR pattern numbers between 0 and 7.
    z = np.sum(y * u, axis=0).astype(np.int8)
    # We get the patterns given by the rule.
    return rule_binary[7-z]


# 4. We now write a function that simulates any elementary cellular automaton. First, we compute the binary representation of the rule (**Wolfram's code**). Then, we initialize the first row of the grid to random values. Finally, we apply the function `step` iteratively on the grid.

# In[ ]:


def generate(rule, size=80, steps=80):
    """Simulate an elementary cellular automaton given its rule
    (number between 0 and 255)."""
    # Compute the binary representation of the rule.
    rule_binary = np.array([int(_) 
                            for _ in np.binary_repr(rule, 8)],
                            dtype=np.int8)
    x = np.zeros((steps, size), dtype=np.int8)
    # Random initial state.
    x[0,:] = np.random.rand(size) < .5
    # Apply the step function iteratively.
    for i in range(steps-1):
        x[i+1,:] = step(x[i,:], rule_binary)
    return x


# 5. Now, we simulate and display 9 different automata.

# In[ ]:


plt.figure(figsize=(6, 6));
rules = [3, 18, 30, 
         90, 106, 110, 
         158, 154, 184]
for i, rule in enumerate(rules):
    x = generate(rule)
    plt.subplot(331+i)
    plt.imshow(x, interpolation='none', cmap=plt.cm.binary);
    plt.xticks([]); plt.yticks([]);
    plt.title(str(rule))


# It has been shown that Rule 110 is **Turing complete** (or **universal**): in principle, this automaton can simulate any computer program.

# > You'll find all the explanations, figures, references, and much more in the book (to be released later this summer).
# 
# > [IPython Cookbook](http://ipython-books.github.io/), by [Cyrille Rossant](http://cyrille.rossant.net), Packt Publishing, 2014 (500 pages).