Rauzy Fractals / Mark McClure

I'm a professor of mathematics at the University of North Carolina Asheville. I've also done a fair amount of consulting work over the years focusing on scientific and data visualization. I'm happy to be part of the first Observable Ambassador's cohort.

Workspace

Public

Edited

Mar 26, 2024

1 fork

3 stars

// Ultimately, everything depends on the substitution function sigma,

// which is denoted with an s in the code.

// s, though, is defined in the radio bitton at the top of the notebook.

// You can see how it's defined here:

[0, 1, 2].map((i) => s(i).toString().replace(/,/g, ""))

// W is supposed to be an infinite word.

// This is just a real long approximation.

// It's obtained by iteration of s

W = {

let i = [0];

while (i.length < 40000) {

i = s(i);

}

return i;

}

// Use the word W to define the stair case.

staircase = {

let basis = [

[1, 0, 0],

[0, 1, 0],

[0, 0, 1]

];

let accumulated = [basis[W[0]]];

for (let i = 1; i < W.length; i++) {

let last = accumulated.slice(-1)[0];

let next = basis[W[i]];

accumulated.push([last[0] + next[0], last[1] + next[1], last[2] + next[2]]);

}

return accumulated;

}

// Use the substitution s to define M.

M = {

let M = math.zeros(3, 3);

for (let j = 0; j < 3; j++) {

let counts = _.countBy(s(j));

for (let i = 0; i < 3; i++) {

M.set([i, j], counts[i] || 0);

}

return M;

}

// Here's what M looks like

M.toArray()

// Compute the eigenvalues and eigenvectors of M

eigenInfo = math.eigs(M)

// Convert the eigenvectors to a change of basis matrix.

eigen_basis = {

let real_evec = math.column(eigenInfo.vectors, real_pos);

let leading_term = real_evec.get([0, 0]);

real_evec = math.re(math.divide(real_evec, leading_term));

let complex_pos = eigenInfo.values

.toArray()

.map((x) => typeof x)

.indexOf("object");

let complex_evec = eigenInfo.vectors.columns()[complex_pos].toArray();

leading_term = complex_evec[0][0];

complex_evec = math.divide(complex_evec, leading_term);

let eigen_basis = math.identity(3);

eigen_basis.subset(math.index([0, 1, 2], 0), real_evec);

eigen_basis.subset(math.index([0, 1, 2], 1), math.re(complex_evec));

eigen_basis.subset(math.index([0, 1, 2], 2), math.im(complex_evec));

return eigen_basis;

}

// Find the position of the real eigenvalue.

// Used to help construct the eigen_basis, as well as for

// the display of the eigenvalues in the exposition

real_pos = eigenInfo.values

.toArray()

.map((x) => typeof x)

.indexOf("number")

// Compute the projection matrix.

P = math.multiply(

eigen_basis,

[

[0, 0, 0],

[0, 1, 0],

[0, 0, 1]

math.inv(eigen_basis)

)

// Project the points in the long staircase and gather them according to

// the last step taken. These points are used to generate the 3D fractal.

gatheredProjectedPoints = {

let pts = staircase.map((pt) => math.multiply(P, pt).toArray());

let gatheredPoints = [[], [], []];

W.forEach((address, idx) => gatheredPoints[address].push(pts[idx]));

return gatheredPoints;

}

// Gather just a few points to display the staircase.

gatheredPoints55 = {

let gatheredPoints = [[], [], []];

W.slice(0, 55).forEach((address, idx) =>

gatheredPoints[address].push(staircase[idx])

);

return gatheredPoints;

}

// Compute an orthogonal viewpoint for the Rauzy fractal in 3D. Based on

// https://observablehq.com/@mcmcclur/understanding-x3d-viewpoints#orientation_string

// Code is a bit of mess since it combines mathjs and ml-matrix.

viewpoint = {

// Set up look, right, and up vectors.

let v1 = math.column(eigen_basis, 1);

let v2 = math.column(eigen_basis, 2);

let u = normalize(math.cross(v1, v2).toArray().flat());

let look = normalize(u.flat());

let right = normalize(cross(look, [0, 0, 1]));

let up = cross(right, look);

// Build a matrix with those vectors as columns

let M = new mlm.Matrix([right, up, look.map((x) => -x)]).transpose();

// Compute the eigenvalue decomposition of the matrix

let eig = new mlm.EigenvalueDecomposition(M);

// Find the position of the 1 in the list of eigenvalues

let tolerance = 10 ** -8;

let b1 = eig.realEigenvalues.map((e) => Math.abs(e - 1) < tolerance);

let b2 = eig.imaginaryEigenvalues.map((e) => Math.abs(e) < tolerance);

let one_position = (b1 && b2).indexOf(true);

// The corresponding eigenvector forms the first three entries of our

// orientation vector.

let v = eig.eigenvectorMatrix.transpose().to2DArray()[one_position];

// Find the position of eigenvalue with positive imaginary part.

let rot_position = eig.imaginaryEigenvalues

.map((y) => y > 10 ** -8)

.indexOf(true);

// The argument of that eigenvalue should be the rotation we need.

let rot = Math.atan2(

eig.imaginaryEigenvalues[rot_position],

eig.realEigenvalues[rot_position]

);

// Push that final value on to our vector and return.

v.push(rot);

return {

position: math

.column(

u.map((x) => [-4 * x]),

)

.flat(), //.toArray().flat(),

orientation: v

};

}

// This regular expression is used to identify strings that look like

// real numbers for formatting in the exposition.

real_regexp = /([0-9]*\.?[0-9]+)/g

delay = 0

plotx3d_canvas_style

import {

show_x3d,

create_pointSet,

create_indexedLineSet,

create_sphere,

create_arrow,

style as plotx3d_canvas_style

} from "@mcmcclur/plotx3d"

import { AffineFunction, DigraphIFS } from "@mcmcclur/digraphifs-class"

import { cross, normalize, rotation_matrix } from "@mcmcclur/x3dom-utilities"

math = require("mathjs@11")

mlm = require("ml-matrix@6")

pics = Promise.all([

FileAttachment("rauzy0.png").image(),

FileAttachment("rauzy1.png").image(),

FileAttachment("rauzy2.png").image(),

FileAttachment("rauzy3.png").image(),

FileAttachment("rauzy4.png").image()

// "blank"

])

Purpose-built for displays of data

Observable is your go-to platform for exploring data and creating expressive data visualizations. Use reactive JavaScript notebooks for prototyping and a collaborative canvas for visual data exploration and dashboard creation.

Learn more