Visualizing Attractors#
An attractor is a set of values to which a numerical system tends to evolve. An attractor is called a strange attractor if the resulting pattern has a fractal structure. This notebook shows how to calculate and plot two-dimensional attractors of a variety of types, using code and parameters primarily from Lázaro Alonso, François Pacull, Jason Rampe, Paul Bourke, and James A. Bednar.
Clifford Attractors#
For example, a Clifford Attractor is a strange attractor defined by two iterative equations that determine the x,y locations of discrete steps in the path of a particle across a 2D space, given a starting point (x0,y0) and the values of four parameters (a,b,c,d):
At each time step, the equations define the location for the following time step, and the accumulated locations show the areas of the 2D plane most commonly visited by the imaginary particle.
It’s easy to calculate these values in Python using Numba. First, we define the iterative attractor equation:
import numpy as np, pandas as pd
from numba import jit
from math import sin, cos, sqrt, fabs
@jit
def Clifford(x, y, a, b, c, d, *o):
return sin(a * y) + c * cos(a * x), \
sin(b * x) + d * cos(b * y)
We then evaluate this equation 10 million times, creating a set of x,y coordinates visited. The @jit here and above is optional, but it makes the code 50x faster.
n=10_000_000
@jit
def trajectory_coords(fn, x0, y0, a, b=0, c=0, d=0, e=0, f=0, n=n):
x, y = np.zeros(n), np.zeros(n)
x[0], y[0] = x0, y0
for i in np.arange(n-1):
x[i+1], y[i+1] = fn(x[i], y[i], a, b, c, d, e, f)
return x,y
def trajectory(fn, x0, y0, a, b=0, c=0, d=0, e=0, f=0, n=n):
x, y = trajectory_coords(fn, x0, y0, a, b, c, d, e, f, n)
return pd.DataFrame(dict(x=x,y=y))
%%time
df = trajectory(Clifford, 0, 0, -1.3, -1.3, -1.8, -1.9)
CPU times: user 1.08 s, sys: 80.4 ms, total: 1.16 s
Wall time: 1.16 s
df.tail()
| x | y | |
|---|---|---|
| 9999995 | 1.857884 | 1.395837 |
| 9999996 | 0.375264 | -0.205512 |
| 9999997 | -1.326024 | -2.301316 |
| 9999998 | 0.423711 | 2.867012 |
| 9999999 | -0.981134 | 1.060115 |
We can now aggregate these 10,000,000 continuous coordinates into a discrete 2D rectangular grid with Datashader, counting each time a point fell into that grid cell:
import datashader as ds
/home/runner/work/examples/examples/attractors/envs/default/lib/python3.11/site-packages/dask/dataframe/__init__.py:31: FutureWarning:
Dask dataframe query planning is disabled because dask-expr is not installed.
You can install it with `pip install dask[dataframe]` or `conda install dask`.
This will raise in a future version.
warnings.warn(msg, FutureWarning)
%%time
cvs = ds.Canvas(plot_width = 700, plot_height = 700)
agg = cvs.points(df, 'x', 'y')
print(agg.values[190:195,190:195],"\n")
[[ 37 41 26 24 24]
[ 28 30 42 27 26]
[111 53 29 37 28]
[145 177 113 63 44]
[ 78 94 136 140 96]]
CPU times: user 391 ms, sys: 7.91 ms, total: 399 ms
Wall time: 398 ms
A small portion of that grid is shown above, but it’s difficult to see the grid’s structure from the numerical values. To see the entire array at once, we can turn each grid cell into a pixel, using a greyscale value from white to black:
ds.transfer_functions.Image.border=0
ds.tf.shade(agg, cmap = ["white", "black"])
As you can see, the most-visited areas of the plane have an interesting structure for this set of parameters. To explore further, let’s wrap up the above aggregation and shading commands into a function so we can apply them more easily:
from datashader.colors import inferno, viridis
def dsplot(fn, vals, n=n, cmap=viridis, label=True):
"""Return a Datashader image by collecting `n` trajectory points for the given attractor `fn`"""
lab = ("{}, "*(len(vals)-1)+" {}").format(*vals) if label else None
df = trajectory(fn, *vals, n=n)
cvs = ds.Canvas(plot_width = 300, plot_height = 300)
agg = cvs.points(df, 'x', 'y')
img = ds.tf.shade(agg, cmap=cmap, name=lab)
return img
And let’s load some colormaps that we can use for subsequent plots:
from colorcet import palette
palette["viridis"]=viridis
palette["inferno"]=inferno
We can now use these colormaps with a pre-selected set of Clifford attractor parameter values (stored in a separate YAML-format text file) to show a wide variety of trajectories that these equations can form:
import yaml
vals = yaml.load(open("data/strange_attractors.yml","r"), Loader=yaml.FullLoader)
def args(name):
"""Return a list of available argument lists for the given type of attractor"""
return [v[1:] for v in vals if v[0]==name]
def plot(fn, vals=None, **kw):
"""Plot the given attractor `fn` once per provided set of arguments."""
vargs=args(fn.__name__) if vals is None else vals
return ds.tf.Images(*[dsplot(fn, v[1:], cmap=palette[v[0]][::-1], **kw) for v in vargs]).cols(4)
plot(Clifford)
0, 0, -1.3, -1.3, -1.8, -1.9 | 0, 0, -1.4, 1.6, 1.0, 0.7 | 0, 0, 1.7, 1.7, 0.6, 1.2 | 0, 0, 1.7, 0.7, 1.4, 2.0 |
0, 0, -1.7, 1.8, -1.9, -0.4 | 0, 0, 1.1, -1.32, -1.03, 1.54 | 0, 0, 0.77, 1.99, -1.31, -1.45 | 0, 0, -1.9, -1.9, -1.9, -1.0 |
0, 0, 0.75, 1.34, -1.93, 1.0 | 0, 0, -1.32, -1.65, 0.74, 1.81 | 0, 0, -1.6, 1.6, 0.7, -1.0 | 0, 0, -1.7, 1.5, -0.5, 0.7 |
Here the values shown are the arguments for the first call to Clifford(x, y, a, b, c, d), with each subsequent call using the x,y location of the previous call.
Randomly sampling the parameter space typically yields much less dramatic patterns, such as all trajectory locations being on a small number of points:
import numpy.random
numpy.random.seed(21)
num = 4
rvals=np.c_[np.zeros((num,2)), numpy.random.random((num,4))*4-2]
plot(Clifford, vals=[["kbc"]+list(rvals[i]) for i in range(len(rvals))], label=True)
0.0, 0.0, -1.8051004767634908, -0.8435613608407326, 0.8838653873249194, -1.9135350003362008 | 0.0, 0.0, -1.176308938930244, -1.7969069732184928, -0.7909124241415326, 0.6556411784987999 | 0.0, 0.0, -0.7675424270524847, 0.33436510487439497, -1.7217161815495978, 1.4696179359723534 | 0.0, 0.0, -1.467037922993009, -1.2875013537620066, -0.01628180069537688, 1.4547985782328943 |
If you wish, Datashader could easily be used to filter out such uninteresting examples, by applying a criterion to the aggregate array before shading and showing only those that remain (e.g. rejecting those where 80% of the pixel bins are empty).
De Jong attractors#
A variety of other sets of attractor equations have been proposed, such as these from Peter de Jong:
@jit
def De_Jong(x, y, a, b, c, d, *o):
return sin(a * y) - cos(b * x), \
sin(c * x) - cos(d * y)
plot(De_Jong)
0, 0, -1.244, -1.251, -1.815, -1.908 | 0, 0, 1.7, 1.7, 0.6, 1.2 | 0, 0, 1.4, -2.3, 2.4, -2.1 | 0, 0, -2.7, -0.09, -0.86, -2.2 |
0, 0, -0.827, -1.637, 1.659, -0.943 | 0, 0, -2.24, 0.43, -0.65, -2.43 | 0, 0, 2.01, -2.53, 1.61, -0.33 | 0, 0, 1.4, 1.56, 1.4, -6.56 |
Svensson attractors#
From Johnny Svensson:
@jit
def Svensson(x, y, a, b, c, d, *o):
return d * sin(a * x) - sin(b * y), \
c * cos(a * x) + cos(b * y)
plot(Svensson)
0, 0, 1.5, -1.8, 1.6, 0.9 | 0, 0, -1.78, 1.29, -0.09, -1.18 | 0, 0, -0.91, -1.29, -1.97, -1.56 | 0, 0, 1.4, 1.56, 1.4, -6.56 |
Bedhead Attractor#
From Ivan Emrich and Jason Rampe:
@jit
def Bedhead(x, y, a, b, *o):
return sin(x*y/b)*y + cos(a*x-y), \
x + sin(y)/b
plot(Bedhead)
1, 1, -0.81, -0.92 | 1, 1, -0.64, 0.76 | 1, 1, 0.06, 0.98 |