ConwayLife.com

Code: Select all

import math
import numpy as np
from scipy.ndimage import convolve
from sys import argv
from PIL import Image

# Simulation method modified from https://scipython.com/blog/simulating-the-belousov-zhabotinsky-reaction/

# Side length of base rhombus (stored as an nXn grid)
n = 200 if len(argv) < 2 else int(argv[1])
# Value for all Belousov-Zhabotinsky reaction coefficients. 
param = 1 if len(argv) < 3 else float(argv[2])
# Number of interations to run.
niter = 500 if len(argv) < 4 else int(argv[3])

# Reaction parameters.
alpha, beta, gamma = param, param, param

# Shape of hexagonal neighborhood out to radius 3.
FILTER_STR = '''
1111000
1111100
1111110
1111111
0111111
0011111
0001111
'''
FILTER_INT = [[int(c) for c in row] for row in FILTER_STR.split()]
FILTER = np.array(FILTER_INT)/sum(sum(row) for row in FILTER_INT)

def update(p,arr):
    """Update arr[p] to arr[q] by evolving in time."""

    # Count the average amount of each species in the neighborhood around each cell
    # by convolution with FILTER
    q = (p+1) % 2
    s = np.zeros((3, arr.shape[2], arr.shape[3]))
    for k in range(3):
        s[k] = convolve(arr[p,k], FILTER, mode='wrap')

    # Apply the reaction equations
    arr[q,0] = s[0] + s[0]*(alpha*s[1] - gamma*s[2])
    arr[q,1] = s[1] + s[1]*(beta*s[2] - alpha*s[0])
    arr[q,2] = s[2] + s[2]*(gamma*s[0] - beta*s[1])

    # Ensure the species concentrations are kept within [0,1].
    np.clip(arr[q], 0, 1, arr[q])
    return arr

def equivalence(i, j, n):
  """Find the equivalent coordinates in the base rhombus by rotating through symmetries."""
  rot = 0
  while True:
    if ((i // n) + (j //n)) % 3 == 0:
      return i % n, j % n 
    i, j = -j - 1, i - j
    rot += 1

# Initialize the base rhombus with random amounts of A, B and C.
data = np.random.random(size=(2, 3, n, n))

# Build a 3nX3n patch based on 120-degree rhombic packing.
# This can be tiled with simple toroidal wrapping at the boundaries.
arr = np.zeros((2, 3, n * 3, n * 3))
for i in range(arr.shape[2]):
  for j in range(arr.shape[2]):
    ti, tj = equivalence(i, j, n)
    arr[:, :, i, j] = data[:, :, ti, tj]

# run the iteration
for i in range(0, niter):
  arr = update(i % 2, arr)

# Create an image from just the top nXn square, equivalent to the base rhombus
disp = np.concatenate((arr[(i + 1) % 2].transpose(2, 1, 0)[:n, :n, :], np.full((n, n, 1), 1, dtype=np.uint8)), axis=2)
img = Image.fromarray((disp * 255).astype(np.uint8), mode='RGBA')
w, h = img.size
img = img.transform((int(w * 3), int(h * math.sqrt(3))),
                    Image.Transform.AFFINE, (0.5, 0.5 / math.sqrt(3), -0.5 * w, 0, 1 / math.sqrt(3), 0),
                    resample=Image.Resampling.BICUBIC) 
img.save('/Users/calla/Downloads/zb3a.png')
img.show()