The Chaos Equations
Dear all,
These are tough times. The world seems to be crashing all around us. And while there's not a lot we can do save from making sure the people that matter the most to us are feel safe and supported, perhaps the nerdier among us can take some small comfort in the beauty of maths.
Enter the chaos equations.
The chaos equations are a set of surprisingly simple rules for updating points on a screen over time. When you make small changes to the way these rules are formulated, you get radically different and beautiful animations out the other side. Before I go into any explanation, check out the video below (weed optional).
Disclaimer: the original idea is ripped straight from this guy , who incidentally does a better job than me at making the animations look just gorgeous. However, his source code is written in pure c++: nobody has time for that. I wanted to write a really simply Python version that anybody can run in a regular python environment.
Song credit Colored World by Kolby Wade .
To create these animations we need to keep track of three variables, \(x, y\) and \(t\) . The first two, \(x\) and \(y\) will represent the coordinates of a single point on the screen, and the third \(t\) will represent time as the animation progresses. We start at some time \(t_1\) , which usually pays to be between roughly \(-2\) and \(2\) . We then discretely increase time across \(N\) time steps to \(t_N\) , also usually between about \(-2\) and \(2\) . At the \(n\) -th time step, we iterate an equation \(M\) times to get \(M\) pairs of coordinates \((x_i, y_i)\) , which we can plot in a variety of colours. The update equation is as follows:
with the initial condition that
The choice of equations \(f(\cdot )\) and \(g(\cdot )\) are what gives rise to the different animations. Any update equations will do, although just simple polynomials seem to create great results. The only other things to choose are the times to run between and the number of iterations to perform.
For those more computationally, minded, an equivalent formulation in pseudo-code would be
for t from t0 to tn:
for i from 1 to M:
if i == 1:
x_i, y_i = t, t
else:
x_i, y_i = f(x_i, y_i, t), g(x_i, y_i, t)
That's literally it! There's nothing more to it. Just give a different colour to each iteration, plot and watch the results.
In order to implement this in Python, I used the popular library Pygame. I wanted to make it easy to explore these animations so I added an axis and the ability to pan and zoom. Check out the sporadically docstring-ed code below.
import pygame
from typing import Union
import numpy as np
from pygame import gfxdraw
import sys
import re
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
from matplotlib import cm
import matplotlib.colors as mpl_colors
colors = {**mpl_colors.CSS4_COLORS, **mpl_colors.TABLEAU_COLORS}
def get_cmap(cmap: str, vmin: float=0, vmax: float=1, alpha: float=1):
cmap = plt.get_cmap(cmap)
cNorm = Normalize(vmin=vmin, vmax=vmax)
scalarMap = cm.ScalarMappable(norm=cNorm, cmap=cmap)
return lambda value: [int(col * 255) for col in scalarMap.to_rgba(value)[:-1]] + [int(alpha * 255)]
def hex_to_RGB(hex_code: str) -> list:
"""
Convert a hex color code into an RGB triplet
Parameters
----------
hex_code A string specifying a hex color code
Returns
-------
RGB A 3-element RGB vector
"""
hex_code = hex_code.lstrip('#')
return [int(hex_code[i:i + 2], 16) for i in (0, 2, 4)]
def is_valid_hex(hex_code: str) -> bool:
"""
Determine whether a string is a valid hex color code
Parameters
----------
hex_code A string
Returns
-------
True is hex_code is a valid hex color code. Else False
"""
match = re.search(r'^#(?:[0-9a-fA-F]{3}){1,2}$', hex_code)
if match:
return True
else:
return False
def prepare_color(color: Union[str, tuple], alpha=1) -> pygame.Color:
"""
Convert a color, which could be in any form, into a usuable RGB tuple
Parameters
----------
color A color in any form
Returns
-------
RGB A (r, g, b, a) tuple
"""
alpha = int(255 * alpha)
if isinstance(color, str):
# if it is a simple hex color code, convert to RGB and return
if is_valid_hex(color):
return pygame.Color(*tuple(hex_to_RGB(color) + [alpha]))
# if it is a preset matplotlib color, return the RGB code associated
elif color in colors:
return pygame.Color(*tuple(hex_to_RGB(colors[color]) + [alpha]))
else:
raise ValueError('{} is not a valid color'.format(color))
elif isinstance(color, tuple):
# RGB tuple - add alpha
if len(color) == 3:
return pygame.Color(*tuple(list(color) + [alpha]))
# nothing to do
elif len(color) == 4:
return pygame.Color(*color)
else:
raise ValueError('Invalid color tuple. Must be RGB(A) but has length {}'.format(len(color)))
elif isinstance(color, pygame.Color):
return color
else:
raise ValueError('color must be a string, tuple or pygame color')
def color_mapper(cmap: str, vmin: float=0, vmax: float=1, alpha: float=1):
"""
Return a function that maps floats between vmin and vmax based on a chosen
color map
Parameters
----------
cmap A valid matplotlib color map string
vmin The minimum float value
vmax The maximum float value
alpha Alpha
Returns
-------
A function that takes a single value and outputs a color
"""
alpha = int(255 * alpha)
cmap = plt.get_cmap(cmap)
cNorm = Normalize(vmin=vmin, vmax=vmax)
scalarMap = cm.ScalarMappable(norm=cNorm, cmap=cmap)
def mapper(value):
"""
This is the function that gets returned
Parameters
----------
value A number between vmin and vmax
Returns
-------
An RGB color
"""
out = scalarMap.to_rgba(value)
if isinstance(out, tuple):
return tuple([255 * out[i] for i in range(3)] + [alpha])
elif isinstance(out, np.ndarray):
out[:, :-1] *= 255
out[:, 3] = alpha
return out
return mapper
class Screen:
"""
Screen class.
"""
def __init__(self, size: tuple=(1920, 1080), xlim: tuple=(0, 1), ylim: tuple=(0, 1), background='black', frame_rate=None, axis=False):
"""
Initialise a Screen object
Parameters
----------
size Screen size in pixels (length, height)
xlim The x-limits. (left, right)
ylim The y-limits. (bottom, top)
screen_color The screen color
"""
self.canvas = pygame.display.set_mode(size)
self.background = prepare_color(background)
self.canvas.fill(self.background)
self.size = self.canvas.get_size()
self.xlim = np.array(xlim).astype(float)
self.ylim = np.array(ylim).astype(float)
self.width = abs(self.xlim[1] - self.xlim[0])
self.height = abs(self.ylim[1] - self.ylim[0])
self.CLICKED = False
self.last_mouse_position = None
self.DRAW_MODE = False
self.raw_lines = []
self.interpolated_lines = []
self.current_line = []
self.n_lines = 0
self.frame_rate = frame_rate
self.clock = pygame.time.Clock()
pygame.font.init()
self.font = pygame.font.SysFont('Comic Sans MS', 30)
self.axis = axis
def draw_cricle(self, pos: tuple, color: Union[str, tuple]=(255, 255, 255), radius: int=5):
if isinstance(color, str):
color = pygame.Color(color)
# pygame.draw.circle(self.canvas, color, self.xy_to_pix(pos), radius)
pygame.gfxdraw.filled_circle(self.canvas, *self.xy_to_pix(pos), radius, color)
def xy_to_pix(self, pos: tuple):
px = int(np.round(self.size[0] * (pos[0] - self.xlim[0]) / (self.xlim[1] - self.xlim[0]), 0))
py = int(np.round(self.size[1] * (self.ylim[1] - pos[1]) / (self.ylim[1] - self.ylim[0]), 0))
return px, py
def pix_to_xy(self, pos: tuple):
x = self.xlim[1] + (self.xlim[0] - self.xlim[1]) * (1 - pos[0] / self.size[0])
y = self.ylim[0] + (self.ylim[1] - self.ylim[0]) * (1 - pos[1] / self.size[1])
return x, y
def xy_array_to_pix(self, pos: np.ndarray):
px = np.round(self.size[0] * (pos[0, :] - self.xlim[0]) / (self.xlim[1] - self.xlim[0]), 0).astype(int)
py = np.round(self.size[1] * (self.ylim[1] - pos[1, :]) / (self.ylim[1] - self.ylim[0]), 0).astype(int)
return np.array([px, py]).T
def plot(self, x: np.ndarray, y: np.ndarray, color: Union[str, tuple, pygame.Color]=(255, 255, 255), width: int=1):
"""
Plot a curve, similar to pyplot.plot
Parameters
----------
x numpy array (n, ) holding x values to plot at
y numpy array (n, ) holding y values to plot
color the line color
width the line width in pixels
"""
pygame.draw.lines(self.canvas, prepare_color(color), False, self.xy_array_to_pix(np.array([x, y])), width)
def add_text(self, text: str, pos: tuple):
text = self.font.render(text, False, (255, 255, 255))
self.canvas.blit(text, self.xy_to_pix(pos))
def get_mouse_position(self):
return self.pix_to_xy(pygame.mouse.get_pos())
def handle_zooming(self, event):
"""
Handle zooming, triggered by a mouse roll
Parameters
----------
event The zoom event
"""
# the zoom factor
zoom = 0.05
# zoom in
if event.button == 4:
self.ylim = (1 - zoom) * self.ylim + zoom * (self.ylim[1] + self.ylim[0]) / 2
self.xlim = (1 - zoom) * self.xlim + zoom * (self.xlim[1] + self.xlim[0]) / 2
# zoom out
if event.button == 5:
self.ylim = (1 + zoom) * self.ylim - zoom * (self.ylim[1] + self.ylim[0]) / 2
self.xlim = (1 + zoom) * self.xlim - zoom * (self.xlim[1] + self.xlim[0]) / 2
self.width = abs(self.xlim[1] - self.xlim[0])
self.height = abs(self.ylim[1] - self.ylim[0])
def handle_panning(self):
"""
Handle panning, triggered by dragging the mouse
"""
if self.CLICKED:
pos_change = self.last_mouse_position - np.array(self.get_mouse_position())
self.xlim += pos_change[0]
self.ylim += pos_change[1]
self.width = abs(self.xlim[1] - self.xlim[0])
self.height = abs(self.ylim[1] - self.ylim[0])
self.last_mouse_position = np.array(self.get_mouse_position())
def handle_drawing(self):
if self.CLICKED:
self.current_line.append(self.get_mouse_position())
if len(self.raw_lines) > self.n_lines:
del self.raw_lines[-1]
self.raw_lines.append(np.array(self.current_line))
self.last_mouse_position = np.array(self.get_mouse_position())
def plot_axis(self):
pad_left = 0.1
pad_right = 0.1
pad_top = 0.1
pad_bottom = 0.2
n_x_ticks = 10
n_y_ticks = 10
x0 = self.xlim[0] + pad_left * self.width
x1 = self.xlim[1] - pad_right * self.width
y0 = self.ylim[0] + pad_bottom * self.height
y1 = self.ylim[1] - pad_top * self.height
# x axis
self.plot(np.array([x0, x1]), np.array([y0, y0]))
# y axis
self.plot(np.array([x0, x0]), np.array([y0, y1]))
# x-ticks
for x in np.linspace(x0, x1, n_x_ticks):
self.plot(np.array([x, x]), np.array([y0 - 0.005 * self.height, y0 + 0.005 * self.height]))
self.add_text(f'{x:.2g}', (x - 0.01 * self.width, y0 - 0.01 * self.height))
# y-ticks
for y in np.linspace(y0, y1, n_y_ticks):
self.plot(np.array([x0 - 0.005 * self.width, x0 + 0.005 * self.width]), np.array([y, y]))
self.add_text(f'{y:.2g}', (x0 - 0.04 * self.width, y + 0.01 * self.height))
# axis labels
self.add_text('x', ((x1 + x0) / 2, y0 - (pad_bottom * self.height) / 4))
self.add_text('y', (x0 - pad_left * self.width / 1.5, (y0 + y1) / 2))
def update(self):
for event in pygame.event.get():
# if we exit
if event.type == pygame.QUIT:
sys.exit()
# handle zooming
if event.type == pygame.MOUSEBUTTONDOWN:
self.handle_zooming(event)
if event.type == pygame.KEYDOWN:
# enter draw mode
if event.key == 100:
if not self.DRAW_MODE:
self.DRAW_MODE = True
else:
self.DRAW_MODE = False
# handle panning
if pygame.mouse.get_pressed()[0]:
if not self.DRAW_MODE:
self.handle_panning()
else:
self.handle_drawing()
self.CLICKED = True
else:
if self.CLICKED:
if self.DRAW_MODE:
self.current_line = []
self.n_lines += 1
self.CLICKED = False
pressed_keys = np.argwhere(pygame.key.get_pressed())
# if we hit ctrl+z, remove last line
if 122 in pressed_keys and 306 in pressed_keys:
if len(self.raw_lines) > 0:
del self.raw_lines[-1]
# quit if we hit esc
if 27 in pressed_keys:
sys.exit()
# draw lines
for line in self.raw_lines:
if line.shape[0] > 1:
self.plot(line[:, 0], line[:, 1])
if self.axis:
self.plot_axis()
pygame.display.update()
self.canvas.fill(self.background)
if self.frame_rate is not None:
self.clock.tick(self.frame_rate)
FRAME_RATE = 30 # frames per second
duration = 30 # animation duration in seconds
I = 1j
presets = {1: {'t0': -0.42,
't1': -0.40,
'x0': -1.03,
'x1': 0.28,
'y0': -0.55,
'y1': 0.40,
'x_': lambda x, y, t: -y ** 2 - x * t + y,
'y_': lambda x, y, t: x ** 2 - x * y + t,
'n_iters': 150,
'trail_length': 10},
2: {'t0': 0.0249,
't1': 0.025,
'x0': -0.4,
'x1': 0.65,
'y0': -0.5,
'y1': 0.7,
'x_': lambda x, y, t: -x ** 2 + x * t + y,
'y_': lambda x, y, t: x ** 2 - y ** 2 - t ** 2 - x * y + y * t -x + y,
'n_iters': 350,
'trail_length': 20},
3: {'t0': 0.400,
't1': 0.400001,
'x0': -0.75,
'x1': 1.7,
'y0': -1.29,
'y1': 0.45,
'x_': lambda x, y, t: x ** 2 + -x * t + y * t - x,
'y_': lambda x, y, t: - y ** 2 - t ** 2 - x * y - y * t - x * t - y,
'n_iters': 130,
'trail_length': 10},
4: {'t0': 0.021,
't1': 0.0215,
'x0': -0.4,
'x1': 0.6,
'y0': -0.5,
'y1': 0.6,
'x_': lambda x, y, t: -x ** 2 + x * t + y,
'y_': lambda x, y, t: x ** 2 - y ** 2 - t ** 2 - x * y + y * t - x + y,
'n_iters': 400,
'trail_length': 10},
5: {'t0': -0.2,
't1': 0,
'x0': -1,
'x1': 1,
'y0': -1,
'y1': 1,
'x_': lambda x, y, t: -x ** 2 + x * t * y + y,
'y_': lambda x, y, t: x ** 2 - y ** 2 - t ** 2 - x * y + y * t - x + y,
'n_iters': 350,
'trail_length': 10},
6: {'t0': 0.10810076,
't1': 0.10810077,
'x0': -1.5,
'x1': 1.3,
'y0': 0,
'y1': 2,
'x_': lambda x, y, t: -t ** 2 - x * y + t,
'y_': lambda x, y, t: - x * y + x * t + y + t,
'n_iters': 800,
'trail_length': 10},
7: {'t0': -1.189,
't1': -1.188,
'x0': -0.8,
'x1': -0.2,
'y0': -2,
'y1': 1,
'x_': lambda x, y, t: t ** -2 + x ** 2 - t ** 2,
'y_': lambda x, y, t: - x * y + y ** 2 - t** 2,
'n_iters': 150,
'trail_length': 10},
8: {'t0': -0.67,
't1': -0.5,
'x0': -1,
'x1': 1,
'y0': -1,
'y1': 1,
'x_': lambda x, y, t: t ** 2 - x * y + x * t,
'y_': lambda x, y, t: y ** 2 - x ** 2 + t ** 2,
'n_iters': 150,
'trail_length': 10},
9: {'t0': 0.19501,
't1': 0.19502,
'x0': -0.9,
'x1': -0.5,
'y0': -0.4,
'y1': -0.1,
'x_': lambda x, y, t: t - x * y + x * t - x * 2,
'y_': lambda x, y, t: y ** 2 - x ** 2 + t ** 2,
'n_iters': 500,
'trail_length': 10},
10: {'t0': -0.788,
't1': -0.786,
'x0': -1.8,
'x1': 0.85,
'y0': -1.2,
'y1': 1.8,
'x_': lambda x, y, t: t - x * y + x * t - x ** 2,
'y_': lambda x, y, t: t ** 2 - y ** 2 + t * x,
'n_iters': 350,
'trail_length': 10},
11: {'t0': -0.00437,
't1': -0.00434,
'x0': -1,
'x1': 1,
'y0': -1,
'y1': 1,
'x_': lambda x, y, t: y + x + t,
'y_': lambda x, y, t: y ** 2 - x ** 2 + t - y,
'n_iters': 500,
'trail_length': 10},
12: {'t0': -0.0001,
't1': 0.035,
'x0': -2,
'x1': 0.2,
'y0': -2,
'y1': 1.2,
'x_': lambda x, y, t: t * x + t ** 2 - y ** 2,
'y_': lambda x, y, t: x ** 2 - t ** 2 + y + x - t,
'n_iters': 300,
'trail_length': 10},
13: {'t0': -0.6,
't1': -0.37,
'x0': -1.5,
'x1': 1.5,
'y0': -1.5,
'y1': 1.5,
'x_': lambda x, y, t: x * t - y * x + y * t - x ** 2 + y ** 2,
'y_': lambda x, y, t: y ** 2 - t ** 2 + x ** 2 - y * t + x - y,
'n_iters': 300,
'trail_length': 10},
14 : {'t0': 0,
't1': 0.9,
'x0': -0.5,
'x1': 1.5,
'y0': -0.08,
'y1': 1.3,
'x_': lambda x, y, t: x * y - t * x + y ** 3,
'y_': lambda x, y, t: np.cos(x ** 0.5),
'n_iters': 200,
'trail_length': 10},
15: {'t0': 1.7605,
't1': 1.760502,
'x0': -2,
'x1': 2,
'y0': -2,
'y1': 2,
'x_': lambda x, y, t: np.cos(t * x) - np.sin(y ** 2),
'y_': lambda x, y, t: x **2 + x * y - t * np.sin(x),
'n_iters': 100,
'trail_length': 10},
16: {'t0': 0.144,
't1': 0.146,
'x0': -0.6,
'x1': 1.5,
'y0': -0.7,
'y1': 2.1,
'x_': lambda x, y, t: np.cos(t * x) - np.sin(y ** 2),
'y_': lambda x, y, t: x ** 2 + x * y - t * np.sin(x),
'n_iters': 60,
'trail_length': 10},
17: {'t0': 0.15,
't1': 0.3,
'x0': -0.5,
'x1': 1,
'y0': -0.8,
'y1': 0.8,
'x_': lambda x, y, t: - y + t + x * y + t ** 2 - y * t,
'y_': lambda x, y, t: + t ** 2 - x * y - t + x ** 2 + x ** 2,
'n_iters': 599,
'trail_length': 10},
18 : {'t0': 1.1,
't1': 1.4,
'x0': -5,
'x1': 5,
'y0': -5,
'y1': 5,
'x_': lambda x, y, t:- y + t + x * y + t ** 2 - y * t ,
'y_': lambda x, y, t:+ t ** 2 - x * y - t + x ** 2 + x ** 2,
'n_iters': 599,
'trail_length': 10},
19: {'t0': -0.53,
't1': -0.4,
'x0': -1,
'x1': 1.5,
'y0': -2,
'y1': 1.6,
'x_': lambda x, y, t: - x ** 2 + x ** 2 + x * y + x * t + y ** 2 ,
'y_': lambda x, y, t: - y + t ** 2 + x * t + x * t + y * t ,
'n_iters': 599,
'trail_length': 10},
20: {'t0': -0.2,
't1': 0.15,
'x0': -0.5,
'x1': 0.5,
'y0': -0.5,
'y1': 0.5,
'x_': lambda x, y, t:+ x ** 2 - y - x ** 2 + x - t ,
'y_': lambda x, y, t: + x * y + y * t - t ** 2 + x + x ** 2 ,
'n_iters': 500,
'trail_length': 10},
20: {'t0': -0.8,
't1': -0.3,
'x0': -0.35,
'x1': 1.15,
'y0': -0.2,
'y1': 1.4,
'x_': lambda x, y, t: - t ** 2 * y - y ** 2 * t - x ** 2 * t - y * t - t ** 2 * y ,
'y_': lambda x, y, t: - t - t ** 2 * x + y - x ** 2 * y - t ** 2 * y ,
'n_iters': 500,
'trail_length': 10}
}
def run(n, color_map='hsv'):
preset = presets[n]
t0, t1 = preset['t0'], preset['t1']
screen = Screen((1920, 1080), xlim=(preset['x0'], preset['x1']), ylim=(preset['y0'], preset['y1']), frame_rate=FRAME_RATE, axis=True)
next_x, next_y = preset['x_'], preset['y_']
n_iters = preset['n_iters']
trail_length = preset['trail_length']
# animation time in seconds
n_frames = duration * FRAME_RATE
ts = np.linspace(t0, t1, n_frames)
output = np.zeros((n_frames, n_iters, 2))
output[:, 0, 0] = ts
output[:, 0, 1] = ts
for i in range(n_iters-1):
output[:, i+1, 0] = next_x(output[:, i, 0], output[:, i, 1], ts)
output[:, i+1, 1] = next_y(output[:, i, 0], output[:, i, 1], ts)
mappers = [get_cmap(color_map, vmin=0, vmax=n_iters, alpha=1 - i / trail_length) for i in range(trail_length)]
colors = [[mappers[i](j) for j in range(n_iters)] for i in range(trail_length)]
for j, t in enumerate(ts):
x1, x0, y1, y0 = screen.xlim[1], screen.xlim[0], screen.ylim[1], screen.ylim[0]
if j % 100 == 0:
print(x1, x0, y1, y0)
for i, points in enumerate(output[j-trail_length:j, :, :]):
for (X, Y), color in zip(points, colors[i]):
if all([X < x1, X > x0, Y < y1, Y > y0]):
screen.draw_cricle((X, Y), radius=2, color=color)
screen.add_text(f't = {t:.9g}', (x0, y1))
screen.update()
run(20, 'CMRmap_r')
' Accent_r, Blues, Blues_r, BrBG, BrBG_r, BuGn, BuGn_r, BuPu, BuPu_r, CMRmap, CMRmap_r, Dark2, Dark2_r, GnBu, GnBu_r, Greens, Greens_r, Greys, Greys_r, OrRd, OrRd_r, Oranges, Oranges_r, PRGn, PRGn_r, Paired, Paired_r, Pastel1, Pastel1_r, Pastel2, Pastel2_r, PiYG, PiYG_r, PuBu, PuBuGn, PuBuGn_r, PuBu_r, PuOr, PuOr_r, PuRd, PuRd_r, Purples, Purples_r, RdBu, RdBu_r, RdGy, RdGy_r, RdPu, RdPu_r, RdYlBu, RdYlBu_r, RdYlGn, RdYlGn_r, Reds, Reds_r, Set1, Set1_r, Set2, Set2_r, Set3, Set3_r, Spectral, Spectral_r, Wistia, Wistia_r, YlGn, YlGnBu, YlGnBu_r, YlGn_r, YlOrBr, YlOrBr_r, YlOrRd, YlOrRd_r, afmhot, afmhot_r, autumn, autumn_r, binary, binary_r, bone, bone_r, brg, brg_r, bwr, bwr_r, cividis, cividis_r, cool, cool_r, coolwarm, coolwarm_r, copper, copper_r, cubehelix, cubehelix_r, flag, flag_r, gist_earth, gist_earth_r, gist_gray, gist_gray_r, gist_heat, gist_heat_r, gist_ncar, gist_ncar_r, gist_rainbow, gist_rainbow_r, gist_stern, gist_stern_r, gist_yarg, gist_yarg_r, gnuplot, gnuplot2, gnuplot2_r, gnuplot_r, gray, gray_r, hot, hot_r, hsv, hsv_r, inferno, inferno_r, jet, jet_r, magma, magma_r, nipy_spectral, nipy_spectral_r, ocean, ocean_r, pink, pink_r, plasma, plasma_r, prism, prism_r, rainbow, rainbow_r, seismic, seismic_r, spring, spring_r, summer, summer_r, tab10, tab10_r, tab20, tab20_r, tab20b, tab20b_r, tab20c, tab20c_r, terrain, terrain_r, twilight, twilight_r, twilight_shifted, twilight_shifted_r, viridis, viridis_r, winter, winter_r'