# Chaos game simulation over a
prescribed driver in Maxima CAS

authors: Krzysztof Le\'{s}niak, Nina Snigireva
date: 2nd June 2020 /// 2vi2020  

file name: Sym5disjunctiveChaosGame_2020-06-02.wxm

/* Clear all variables if needed */
kill(all)$

Numerical vs symbolic computation.
In principle, it is possible to run orbits for the chaos game as symbolic computations.  
However, for longer orbits it becomes impossible to manage them 
even by employing an automatic simplification in Maxima. 
Therefore, all values below are turned into bigfloats.  

Recall that code lines ending with "$" do not show the output of a computation 
and the lines ending with semicolon ";" show the output of a computation. 
(It is not advisable to print lengthy outputs like reading-in the driver from a file or the list of points in a computed orbit.)

## IFS

Define maps f1,f2,f3... comprising 
an iterated function system (IFS) fn=[f1,f2,f3,...].
The maps operate on vectors (1-column matrices).
(The IFSs below consist of affine maps given by matrices.)

The examples below arise in the investigation of the following problem:
Given an IFS F= {fi : i in I} and rotation R with R^^p = id, what is the relation between (semi) attractors of various IFSs.
Namely 
F~= {fi, R : i in I}, 
G= {R^^(-j) o fi o R^^j : i in I, j in Zp}, 
G~ =  {fi o R^^j : i in I, j in Zp}, 
G^ =  {R^^k o fi : i in I, k in Zp}.

I = {1,2,3} standard Sierpiński triangle, I={1} Sierpiński triangle by rotation, I={1,2,3,4} Barnsley fern, 
I={1,2,3,4,5} standard Sierpiński pentagon, I={1} Sierpiński pentagon by rotation.

### Equilateral Sierpiński triangle and its rotational symmetry

/* Standard encoding of the Sierpiński triangle with vertices B1, B2, B3 (third roots of unity).
Note that {fi : i=1,2,3} = {R^^(- (j-1)) o f1 o R^^(j-1) : j=1,2,3}, 
where R is the rotation by 120deg about the origin  O */

B1: matrix([1],[0]);
B1: bfloat(B1)$
f1(v) := (1/2) * (v + B1)$

B2: matrix([-1/2],[sqrt(3)/2]);
B2: bfloat(B2)$
f2(v) := (1/2) * (v + B2)$

B3: matrix([-1/2],[-sqrt(3)/2]);
B3: bfloat(B3)$
f3(v) := (1/2) * (v + B3)$

fn: [f1,f2,f3]$

/* Semiattractor of  F~= {f1,f2},
where  f2=R --- clockwise rotation by 120deg about the origin O, F={f1}, and
the IFS G={R^^(- (j-1)) o f1 o R^^(j-1) : j=1,2,3} encodes Sierpiński triangle */

R:  matrix(
 [-1/2,-sqrt(3)/2], 
 [sqrt(3)/2,-1/2]
);
R: bfloat(R)$
f2(v) := R . v$

B1: matrix([1],[0]);
B1: bfloat(B1)$
f1(v) := (1/2) * (v + B1)$
fn: [f1,f2]$

/* Attractor of  G~ = {f1 o R^^j : j=0,1,2},
where  R --- clockwise rotation by 120deg about the origin O, F={f1} and
G = {R^^(- (j-1)) o f1 o R^^(j-1) : j=1,2,3} encodes Sierpiński triangle with vertices at the third roots of unity: 
Bj=R^^(j-1)B1,  j=1,2,3 */

R:  matrix(
 [-1/2,-sqrt(3)/2], 
 [sqrt(3)/2,-1/2]
);
RR: R . R;
R: bfloat(R)$
RR: bfloat(RR)$

B1: matrix([1],[0]);
B1: bfloat(B1)$
f1(v) := (1/2) * (v + B1)$

f2(v) := (1/2) * (R . v + B1)$
f3(v) := (1/2)  * (RR . v + B1)$

fn: [f1,f2,f3]$

/* Attractor of G^ = {f1, R o f1, R^^2 o f1}, where
R --- 120deg rotation about the origin O, F={f1}.
Compare with the standard encoding of Sierpiński triangle 
via IFS G = {R^^(- (j-1)) o f1 o R^^(j-1) : j=1,2,3} */

R:  matrix(
 [-1/2,-sqrt(3)/2], 
 [sqrt(3)/2,-1/2]
);
RR: R . R;
R: bfloat(R); 
RR: bfloat(RR);

B1: matrix([1],[0]);
B1: bfloat(B1);
f1(v) := (1/2) * (v + B1)$

f2(v) := R . f1(v)$
f3(v) := RR . f1(v)$

fn: [f1,f2,f3]$

### Sierpiński triangle enriched with 40 deg clockwise rotation

/* Semiattractor of  F~= {fi : i =1,2,3,4},
where  F={f1,f2,f3} encodes the Sierpiński triangle with vertices 2*B1, 2*B2, 2*B3; 
f4=R --- rot by 40deg about the origin O */

alpha: 2*%pi/9$
c: cos(alpha)$
s: sin(alpha)$
R:  matrix(
 [c, s], 
 [-s,c]
);
R: bfloat(R);
f4(v) := R . v$

B1: matrix([-1/4],[sqrt(3)/4]);
B1: bfloat(B1);
f1(v) := (1/2) * v + B1$

B2: matrix([1/2],[0]);
B2: bfloat(B2);
f2(v) := (1/2) * v + B2$

B3: matrix([-1/4],[-sqrt(3)/4]);
B3: bfloat(B3);
f3(v) := (1/2) * v + B3$

fn: [f1,f2,f3,f4]$

/* Attractor of  G~= {fi o R^^j : i =1,2,3; j=0,1,2,3,4,5,6,7,8},
where  F={f1,f2,f3} encodes Sierpiński triangle with vertices 2*B1, 2*B2, 2*B3; 
R --- rot by 40deg about the origin O */

alpha: 2*%pi/9$
c: cos(alpha)$
s: sin(alpha)$
R:  matrix(
 [c, s], 
 [-s,c]
)$
R2: R^^2$
R3: R^^3$
R4: R^^4$
R5: R^^5$
R6: R^^6$
R7: R^^7$
R8: R^^8$
R: bfloat(R)$
R2: bfloat(R2)$
R3: bfloat(R3)$
R4: bfloat(R4)$
R5: bfloat(R5)$
R6: bfloat(R6)$
R7: bfloat(R7)$
R8: bfloat(R8)$

B1: matrix([-1/4],[sqrt(3)/4]);
B1: bfloat(B1);
f1(v) := (1/2) * v + B1$

B2: matrix([1/2],[0]);
B2: bfloat(B2);
f2(v) := (1/2) * v + B2$

B3: matrix([-1/4],[-sqrt(3)/4]);
B3: bfloat(B3);
f3(v) := (1/2) * v + B3$

f4(v) := f1(R . v)$
f5(v) := f2(R . v)$
f6(v) := f3(R . v)$

f7(v) := f1(R2 . v)$
f8(v) := f2(R2 . v)$
f9(v) := f3(R2 . v)$

f10(v) := f1(R3 . v)$
f11(v) := f2(R3 . v)$
f12(v) := f3(R3 . v)$

f13(v) := f1(R4 . v)$
f14(v) := f2(R4 . v)$
f15(v) := f3(R4 . v)$

f16(v) := f1(R5 . v)$
f17(v) := f2(R5 . v)$
f18(v) := f3(R5 . v)$

f19(v) := f1(R6 . v)$
f20(v) := f2(R6 . v)$
f21(v) := f3(R6 . v)$

f22(v) := f1(R7 . v)$
f23(v) := f2(R7 . v)$
f24(v) := f3(R7 . v)$

f25(v) := f1(R8 . v)$
f26(v) := f2(R8 . v)$
f27(v) := f3(R8 . v)$

fn: [f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18,f19,f20,f21,f22,f23,f24,f25,f26,f27]$

### Sierpiński pentagon and its rotational symmetry 

/* Semiattractor of  F~= {f1,f2},
where  f2=R --- rot by alpha=72deg about the origin O, F={f1}, and
the IFS G={R^^(- (j-1)) o f1 o R^^(j-1) : j=1,2,3,4,5} encodes the Sierpiński pentagon */

a: (3-sqrt(5))/2$
A1:  matrix(
 [a, 0], 
 [0, a]
)$
A1: bfloat(A1);
B1: matrix([1-a],[0])$
B1: bfloat(B1);
f1(v) := A1 . v + B1$

alpha: 2*%pi/5$
c: cos(alpha)$
s: sin(alpha)$
R:  matrix(
 [c, s], 
 [-s,c]
)$
R: bfloat(R);
f2(v):= R . v$

fn: [f1,f2]$

/* Attractor of  G~= {f1 o R^^j : j=0,1,2,3,4},
where  R --- rot by alpha=72deg about the origin O, F={f1}, and
the IFS G={R^^(- (j-1)) o f1 o R^^(j-1) : j=1,2,3,4,5} encodes the Sierpiński pentagon */

alpha: 2*%pi/5$
c: cos(alpha)$
s: sin(alpha)$
R:  matrix(
 [c, s], 
 [-s,c]
);
R: bfloat(R);

a: (3-sqrt(5))/2$
A1:  matrix(
 [a, 0], 
 [0, a]
);
A1: bfloat(A1);
B1: matrix([1-a],[0]);
B1: bfloat(B1);
f1(v) := A1 . v + B1$

f2(v) := f1(R . v)$
f3(v) := f1(R^^2 . v)$
f4(v) := f1(R^^3 . v)$
f5(v) := f1(R^^4 . v)$

fn: [f1,f2,f3,f4,f5]$

/*  Sierpiński pentagon by Field & Golubitsky "Symmetry in Chaos" : 
attractor of G^ = {f1, R o f1, R^^2 o f1, R^^3 o f1, R^^4 o f1}, 
R --- alpha=72deg rotation about Origin, F={f1}.
Compare with the standard encoding of Sierpiński pentagon 
via IFS G = {R^^(- (j-1)) o f1 o R^^(j-1) : j=1,2,3,4,5} 
as in Schlicker & Dennis (1995) "Sierpinski n-gons," Pi Mu Epsilon Journal Vol.10
and in L. Riddle http://larryriddle.agnesscott.org/ifs/pentagon/pentagon.htm */

alpha: 2*%pi/5$
c: cos(alpha)$
s: sin(alpha)$
R:  matrix(
 [c, s], 
 [-s,c]
);
R: bfloat(R);

a: (3-sqrt(5))/2$
A1:  matrix(
 [a, 0], 
 [0, a]
);
A1: bfloat(A1);
B1: matrix([1-a],[0]);
B1: bfloat(B1);
f1(v) := A1 . v + B1$

f2(v) := R . f1(v)$
f3(v) := R^^2 . f1(v)$
f4(v) := R^^3 . f1(v)$
f5(v) := R^^4 . f1(v)$

fn: [f1,f2,f3,f4,f5]$

If one really needs more powerful trig simplification, then package "ntrig" could be of help.

load(ntrig)$

/* Angle alpha=72deg is not well handled in pure Maxima; "ntrig" resolves the issue */
alpha: 2*%pi/5$
c: cos(alpha)$
s: sin(alpha)$
R:  matrix(
 [c, s], 
 [-s,c]
);
M: R^^5;
trigsimp(M);

### Barnsley fern enriched with 180 deg rotation

/* Barnsley fern IFS F={f1,f2,f3,f4} */
A1:  matrix(
 [0, 0], 
 [0, 0.16]
);
A1: bfloat(A1);
B1: matrix([0],[0]);
B1: bfloat(B1);
f1(v) := A1 . v + B1$
A2:  matrix(
 [0.85, 0.04], 
 [-0.04, 0.85]
);
A2: bfloat(A2);
B2: matrix([0],[1.6]);
B2: bfloat(B2);
f2(v) := A2 . v + B2$
A3:  matrix(
 [0.2, -0.26], 
 [0.23, 0.22]
);
A3: bfloat(A3);
B3: matrix([0],[1.6]);
B3: bfloat(B3);
f3(v) := A3 . v + B3$
A4:  matrix(
 [-0.15, 0.28], 
 [0.26, 0.24]
);
A4: bfloat(A4);
B4: matrix([0],[0.44]);
B4: bfloat(B4);
f4(v) := A4 . v + B4$
fn: [f1,f2,f3,f4]$

/* Attractor of  G~= {fi o R^^j : i =1,2,3,4; j=0,1} 
= {f1,f2,f3,f4,f5=f1 o R,f6=f2 o R, f7=f3 o R, f8=f4 o R},
where  F={f1,f2,f3,f4} encodes Barnsley fern, 
R --- rot by alpha=180deg about the origin O */

alpha: %pi$
c: cos(alpha)$
s: sin(alpha)$
R:  matrix(
 [c, s], 
 [-s,c]
);
R: bfloat(R);

A1:  matrix(
 [0, 0], 
 [0, 0.16]
);
A1: bfloat(A1);
B1: matrix([0],[0]);
B1: bfloat(B1);
f1(v) := A1 . v + B1$

A2:  matrix(
 [0.85, 0.04], 
 [-0.04, 0.85]
);
A2: bfloat(A2);
B2: matrix([0],[1.6]);
B2: bfloat(B2);
f2(v) := A2 . v + B2$

A3:  matrix(
 [0.2, -0.26], 
 [0.23, 0.22]
);
A3: bfloat(A3);
B3: matrix([0],[1.6]);
B3: bfloat(B3);
f3(v) := A3 . v + B3$

A4:  matrix(
 [-0.15, 0.28], 
 [0.26, 0.24]
);
A4: bfloat(A4);
B4: matrix([0],[0.44]);
B4: bfloat(B4);
f4(v) := A4 . v + B4$

f5(v) := f1(R . v)$
f6(v) := f2(R . v)$
f7(v) := f3(R . v)$
f8(v) := f4(R . v)$

fn: [f1,f2,f3,f4,f5,f6,f7,f8]$

/* Semiattractor of  F~= {fi : i =1,2,3,4,5},
where  F={f1,f2,f3,f4} encodes Barnsley fern, 
f5=R --- rot by alpha=180deg about the origin O */

alpha: %pi$
c: cos(alpha)$
s: sin(alpha)$
R:  matrix(
 [c, s], 
 [-s,c]
);
R: bfloat(R);
f5(v) := R . v;

A1:  matrix(
 [0, 0], 
 [0, 0.16]
);
A1: bfloat(A1);
B1: matrix([0],[0]);
B1: bfloat(B1);
f1(v) := A1 . v + B1$

A2:  matrix(
 [0.85, 0.04], 
 [-0.04, 0.85]
);
A2: bfloat(A2);
B2: matrix([0],[1.6]);
B2: bfloat(B2);
f2(v) := A2 . v + B2$

A3:  matrix(
 [0.2, -0.26], 
 [0.23, 0.22]
);
A3: bfloat(A3);
B3: matrix([0],[1.6]);
B3: bfloat(B3);
f3(v) := A3 . v + B3$

A4:  matrix(
 [-0.15, 0.28], 
 [0.26, 0.24]
);
A4: bfloat(A4);
B4: matrix([0],[0.44]);
B4: bfloat(B4);
f4(v) := A4 . v + B4$

fn: [f1,f2,f3,f4,f5]$

### Flat Sierpiński triangle in 3d enriched with 90 deg Oy rotation

/* Semiattractor of  F~= {fi : i =1,2,3,4},
where  F={f1,f2,f3} encodes flat Sierpiński triangle in 3d with vertices B1,B2,B3, 
f4=R --- rot by 90deg about the Oy axis */

B1: matrix([1], [0], [0]);
B1: bfloat(B1)$
f1(v) := (1/2) * (v + B1)$

B2: matrix([-1/2], [sqrt(3)/2], [0]);
B2: bfloat(B2)$
f2(v) := (1/2) * (v + B2)$

B3: matrix([-1/2], [-sqrt(3)/2], [0]);
B3: bfloat(B3)$
f3(v) := (1/2) * (v + B3)$

alpha: %pi/2$
c: cos(alpha)$
s: sin(alpha)$
R: matrix([c, 0, s], [0,1,0], [-s,0,c]);
R: bfloat(R)$
f4(v) := R . v$

fn: [f1,f2,f3,f4]$

/* Attractor of  G~= {fi o R^^j : i =1,2,3; j=0,1,2,3},
where  F={f1,f2,f3} encodes flat Sierpiński triangle in 3d with vertices B1,B2,B3, 
R --- rot by 90deg about the Oy axis */

B1: matrix([1], [0], [0])$
B1: bfloat(B1)$
f1(v) := (1/2) * (v + B1)$

B2: matrix([-1/2], [sqrt(3)/2], [0])$
B2: bfloat(B2)$
f2(v) := (1/2) * (v + B2)$

B3: matrix([-1/2], [-sqrt(3)/2], [0])$
B3: bfloat(B3)$
f3(v) := (1/2) * (v + B3)$

alpha: %pi/2$
c: cos(alpha)$
s: sin(alpha)$
R: matrix([c, 0, s], [0,1,0], [-s,0,c])$
R: bfloat(R);

f4(v) := f1(R . v)$
f5(v) := f2(R . v)$
f6(v) := f3(R . v)$

f7(v) := f1(R^^2 . v)$
f8(v) := f2(R^^2 . v)$
f9(v) := f3(R^^2 . v)$

f10(v) := f1(R^^3 . v)$
f11(v) := f2(R^^3 . v)$
f12(v) := f3(R^^3 . v)$

fn: [f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12]$

## Driver

Define a driver along which the chaos game orbit will be generated.
The driver is precomputed and stored in a separate file (as comma separated numbers 0,1,2,...) or explicitly in the script.  
To compute the driver one can, for example, use the algorithm implemented in Python by E. Bach: https://git.sagemath.org/sage.git/tree/src/sage/combinat/debruijn_sequence.pyx 

/* Load driver from a file and print its length */
driver: read_list("debruijndriver.txt", comma)$
length(driver);

/* Direct driver: de Bruijn word on the cycle 
over 3-letter alphabet {0,1,2} for strings of length up to 4  */
driver:  [0, 0, 0, 0,  1, 0, 0, 0,  2, 0, 0, 1,  1, 0, 0, 1,  2, 0, 0, 2,  
    1, 0, 0, 2,  2, 0, 1, 0,  1, 0, 2, 0,  1, 1, 1, 0,  1, 1, 2, 0,  
    1, 2, 1, 0,  1, 2, 2, 0,  2, 0, 2, 1,  1, 0, 2, 1,  2, 0, 2, 2,  
    1, 0, 2, 2,  2, 1, 1, 1,  1, 2, 1, 1,  2, 2, 1, 2,  1, 2, 2, 2,  
    2, 0, 0, 0]$
length(driver);

## Chaos game

Define the hopping procedure over IFS functions fn.
Choose a starting point (either in 2d or in 3d depending on the IFS).
Iterate IFS maps, fn, over a prescribed driver to generate the orbit.

/* Orbit hopping; 
v --- vector representing a point, 
(symbol+1) --- the index of a map fn chosen to act on v. 
Note that indices in Maxima's list  are counted from 1 
while symbols in a driver are counted from 0 */
g(v,symbol) := fn[symbol+1](v)$

/* Starting point for IFS on 2d (two dimensional space) */
startpnt: matrix([0],[0])$
startpnt: bfloat(startpnt)$

/* Starting point for IFS on 3d (three dimensional space) */
startpnt: matrix([0],[0],[0])$
startpnt: bfloat(startpnt)$

/* Main procedure: iteration and dumping chaos game orbit.
Points are stored as d-lists (d= dimension of vector/point v); 
orbit is a list of d-lists */ 
location: startpnt$
orbit: [list_matrix_entries(location)]$
for i in driver do
    block(
    location: g(location,i),
    orbit: append(orbit, [list_matrix_entries(location)])
)$

## Drawing orbit

/* Call package "draw" dedicated to plotting. 
(draw: Gnuplot by M.R. Riotorto)
http://riotorto.users.sourceforge.net/Maxima/gnuplot/index.html */
load(draw)$

/* Ancillary variables */
picture: orbit$
frameL: -1$
frameR: 1.5$
frame: 1.2$

/* Scene for 2d plot.
Set x,y tick marks;  range of x,y axis; 
color, ptype - representing symbol, and size of points;
ptype= filled_diamant, circle, dot */
scene: gr2d(
xtics  = [frameL, 0.5, frameR],
ytics  = [frameL, 0.5, frameR],
xrange = [frameL, frameR],
yrange = [frameL, frameR],
color = red,
point_type = filled_diamant,
point_size = 0.2,
points(picture)
)$

/* Scene for 3d plot. 
Set x,y,z tick marks; range of x,y,z axis; 
viewing angles of coordinate system;
color, ptype - representing symbol, and size of points;
ptype= filled_diamant, circle, dot */
scene: gr3d(
xtics  = [-frameL,0.25,frameR],
ytics  = [-frameL,0.25,frameR],
ztics  = [-frame,0.5,frame], 
xrange = [-frameL,frameR],
yrange = [-frameL,frameR],
zrange = [-frame,frame],
view = [50,20],
color = red,
point_type = filled_diamant,
point_size = 0.25,
points(picture)
)$

/* Plot into file; file types: 'png, 'eps, 'jpg; 
dimensions of image in pixels */
draw(file_name = "AttractorFromChaosGame",
    terminal   = 'jpg,
    dimensions  = [600,600],
    scene)$

/* Plot directly to Maxima output 
(Not recommended for large orbits) */
draw(
    dimensions  = [500,500],
    scene)$