30:T621,viewof s = { let sigma_choices = [ function s(i) { if (i == 0) { return [0, 1]; } else if (i == 1) { return [0, 2]; } else if (i == 2) { return [0]; } else { return i.map(s).flat(); } }, function s(i) { if (i == 0) { return [0, 0, 1]; } else if (i == 1) { return [0, 0, 2]; } else if (i == 2) { return [0]; } else { return i.map(s).flat(); } }, function s(i) { if (i == 0) { return [0, 1]; } else if (i == 1) { return [2]; } else if (i == 2) { return [0]; } else { return i.map(s).flat(); } }, function s(i) { if (i == 0) { return [0, 0, 1]; } else if (i == 1) { return [0, 2]; } else if (i == 2) { return [0]; } else { return i.map(s).flat(); } }, function s(i) { if (i == 0) { return [0, 1, 0, 2]; } else if (i == 1) { return [2]; } else if (i == 2) { return [0]; } else { return i.map(s).flat(); } } ]; let input = Inputs.radio(sigma_choices, { value: sigma_choices[0], format: (s) => tex`\Large \sigma(0) = ${s(0) .toString() .replace(/,/g, "")}, \: \sigma(1) = ${s(1) .toString() .replace(/,/g, "")}, \: \sigma(2) = ${s(2) .toString() .replace(/,/g, "")}` }); d3.select(input).selectAll("label").style("display", "block"); return input; }31:T72e,## Strings and substitutions When we say that ${tex`\sigma`} is a *substitution* on the three symbols ${tex`0,1,`} and ${tex`2`}, we mean that ${tex`\sigma`} is a function that accepts one of those symbols and returns a string of one or more of them back. In our working example, selected by the radio buttons above, we have - ${tex`\sigma(0) = ${s(0).toString().replace(/,/g,'')}`}, - ${tex`\sigma(1) = ${s(1).toString().replace(/,/g,'')}`}, and - ${tex`\sigma(2) = ${s(2).toString().replace(/,/g,'')}`}. This extends naturally via concatenation to a finite or even an infinite string of these symbols. For example: ${tex.block`\sigma(0102) = ${[0, 1, 0, 2] .map((a) => s(a).toString().replace(/,/g, "")) .toString() .replace(/,/g, "\\:")} = ${s([0,1,0,2]).toString().replace(/,/g, '')}.`} We can iterate this operation to produce a sequence of these strings. If the starting point of our iteration is just the single character 0, the first few terms in our sequence are:
${(function() { let i = 0; let ii = ` ${i}`; for (let cnt = 0; cnt < 6; cnt++) { i = s(i); ii = ii + "
\n" + " " + i.toString().replace(/,/g, ""); } return ii; })()}
This sequence converges to an infinitely long string that is fixed under application of ${tex`\sigma`}; let's call that string ${tex`W`}. The strings that we generate in this process are often called *words*, which is why we use the symbol ${tex`W`} above. We represent an arbitrary word (finite or infinite) with a ${tex`w`} and denote the symbol in the ${tex`i^{\text{th}}`} position by ${tex`w_i`}. If we want to call these things "words", then it might make more sense to use letters as our symbols than to use numbers. Numbers are easier to translate into positions in space, though, which is our task in the next section.32:T736,## Linear algebra and projection The mathematical level needs to step up a bit to describe the details of the projection. In particular, we need to use a fair amount of linear algebra, including real and complex eigenvalues and eigenvectors, along with change of bases using those tools. Of course, we'll probably need a matrix, as well. In particular, we'll use the *incidence matrix* ${tex`M`} for the substitution ${tex`\sigma.`} If the entries of ${tex`M`} are denoted ${tex`m_{ij}`}, then ${tex.block`m_{ij} = \# \text{ occurrences of } i \text{ in } \sigma(j).`} For our working example, we have ${tex.block`M = ${math.parse(M.toString()).toTex()}.`} It turns out that multiplication by the matrix ${tex`M`} is related to application of the substitution ${tex`\sigma`} via a function ${tex`L`} that maps words into ${tex`R^3.`} ${tex`L`} is defined by ${tex.block`L(w_1 w_2 \cdots w_k) = \vec{e}_{w_1} + \vec{e}_{w_2} + \cdots + \vec{e}_{w_k}.`} ${tex`L`} relates ${tex`M`} to ${tex`\sigma`} via the formula ${tex.block`L(\sigma(w)) = ML(w).`} It's actually pretty easy to understand this formula, if you think of matrix multiplication in the following way: The product ${tex`M\vec{e}_i`} yields exactly the ${tex`i^{\text{th}}`} column of ${tex`M`}. Of course, those columns are chosen exactly to tell us how many of each symbol appears in ${tex`\sigma(i)`}. Thus, ${tex`L(\sigma(w))`} and ${tex`ML(w)`} yield exactly the sum of the same of basis vectors. The relationship between ${tex`\sigma`} and ${tex`M`} provided by the map ${tex`L`} is called a *semi-conjugacy* and it allows us to relate application of the substitution operator to multiplication by a matrix. This is great, since linear algebra is a well understood topic with lots of tools to help us understand it. In particular, the eigenvalues of the matrix ${tex`M`} are33:T609,## An alternative projection Since the Rauzy fractal lives in a plane, it's natural to try to generate a two-dimensional picture of it. One simple way to do so would be to rotate the 3D image we already have so that it is parallel to the ${tex`xy`}-plane and then drop the third coordinate. There is a much more direct approach, though, that takes advantage of the complex structure of the eigenvectors and leads to a much simpler computation and code. As it turns out, if we take the transpose of the matrix ${tex`M`} and compute its complex eigenvectors, then they are orthogonal to the one-dimensional real eigenspace of ${tex`M.`} This suggests that we might project the staircase onto one of those, grab the scalar component, and plot the real and imaginary components in the plane. Now, there's a very simple formula for the projection of a vector ${tex`\vec{v}`} onto the vector ${tex`u`} that takes on an even simpler form when ${tex`\vec{u}`} is normalized so that its magnitude is one, namely: ${tex.block`\text{proj}_{\vec{u}}\vec{v} = \frac{\vec{v}\cdot\vec{u}}{\vec{u}\cdot{\vec{u}}}\vec{u} = (\vec{v}\cdot\vec{u})\vec{u}.`} We just want the scalar component, namely the ${tex`\vec{v}\cdot\vec{u}`}. Thus, once we have the stair case and matrix ${tex`M`}, all we need to do is compute a complex eigenvector of ${tex`M^T`}, normalize it to obtain ${tex`\vec{u}`} and then compute the dot product ${tex`\vec{u}\cdot\vec{v}`} where ${tex`\vec{v}`} ranges all the points in the staircase. Here's the result for our working example:34:T65c,// We delay evaluation so that other cells render quicker per this forum post: // https://talk.observablehq.com/t/rauzy-fractals-and-delayed-evaluation/7453 Promises.delay(delay).then(() => { // It's just a few lines to compute the points in the fractal. let M1 = math.transpose(M); let eigenInfo = math.eigs(M1); let u = math.column(eigenInfo.vectors, 1); // Math.js computes eigenvectors using an iterative algorithm // that has a random component. This next ensures consistency. if (u.get([0, 0]).re > 0) { u = math.multiply(-1, u); } let pts = staircase.map((pt) => math.dot(u, pt)); // Set up the scales let [xmin, xmax] = d3.extent(pts.map((z) => z.re)); let [ymin, ymax] = d3.extent(pts.map((z) => z.im)); let w = d3.min([800, width]); let h = (w * (ymax - ymin)) / (xmax - xmin); let pad = 20; let x_scale = d3 .scaleLinear() .domain([xmin, xmax]) .range([pad, w - pad]); let y_scale = d3 .scaleLinear() .domain([ymin, ymax]) .range([h - pad, pad]); // Setup the canvas and draw the points let canvas = d3.create("canvas").attr("width", w).attr("height", h); let context = canvas.node().getContext("2d"); pts.forEach(function (pt, i) { context.fillStyle = d3.schemeCategory10[W[i]]; context.fillRect(x_scale(pt.re), y_scale(pt.im), 3, 3); }); // If axes are desired, we could use draw_axes defined in // import { draw_axes } from "@mcmcclur/iteratedfunctionsystem-class" // draw_axes( // context, // [ // [xmin, xmax], // [ymin, ymax] // ], // x_scale, // y_scale // ); return canvas.node(); })35:T6f6,// 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]), 0 ) .flat(), //.toArray().flat(), orientation: v }; }36:[{"id":"ccd9227d4c4353220b3cac2982f5892b6e065e8c2ef8cc2f9ac5bfbc31ae665a52bdfc3235c8b52c3b1bfbf56787a4294174976aaf7ce27d079426ad4e2ee78f","href":"https://static.observableusercontent.com/files/ccd9227d4c4353220b3cac2982f5892b6e065e8c2ef8cc2f9ac5bfbc31ae665a52bdfc3235c8b52c3b1bfbf56787a4294174976aaf7ce27d079426ad4e2ee78f","name":"rauzy0.png","type":"image/png","lastModified":1671805879842.637,"size":323716},{"id":"00c4226fb6c5e4deea1655ec656a7f05b9749116e7ff4d519d468305b7320774ab410c255cb3ccf7af5dfbf3f7668d618c06ed6993d7b8e4883a44f1e52480a9","href":"https://static.observableusercontent.com/files/00c4226fb6c5e4deea1655ec656a7f05b9749116e7ff4d519d468305b7320774ab410c255cb3ccf7af5dfbf3f7668d618c06ed6993d7b8e4883a44f1e52480a9","name":"rauzy1.png","type":"image/png","lastModified":1671805879740.051,"size":260223},{"id":"5702ded5ffa0219d535cde9c653a1951f7b37a8cfed9a6b9e24ee70c1e49c924b1ab687e25c03e7e5780a28de371d139811ae92a4d90044f02da2e1e6677b9eb","href":"https://static.observableusercontent.com/files/5702ded5ffa0219d535cde9c653a1951f7b37a8cfed9a6b9e24ee70c1e49c924b1ab687e25c03e7e5780a28de371d139811ae92a4d90044f02da2e1e6677b9eb","name":"rauzy2.png","type":"image/png","lastModified":1671805879606.287,"size":251753},{"id":"be6b2ef21580021544861acd2a98ec97373cf279a2e76d666e1b9312ed63943c4fafaca05f680c8ba2e6316a2fe7e98dd5b64befa35fffd57e4e0639d909cd30","href":"https://static.observableusercontent.com/files/be6b2ef21580021544861acd2a98ec97373cf279a2e76d666e1b9312ed63943c4fafaca05f680c8ba2e6316a2fe7e98dd5b64befa35fffd57e4e0639d909cd30","name":"rauzy3.png","type":"image/png","lastModified":1671805879834.237,"size":513495},{"id":"d895dd4e3ed8bfa9b5adcb6ac6f7f96b2e173eae3655c7c1c38d2320778ce322b97b0a6554393baf9ee90b1ff013425b2bbfe7f345eb15909c2d7f3c69a4155e","href":"https://static.observableusercontent.com/files/d895dd4e3ed8bfa9b5adcb6ac6f7f96b2e173eae3655c7c1c38d2320778ce322b97b0a6554393baf9ee90b1ff013425b2bbfe7f345eb15909c2d7f3c69a4155e","name":"rauzy4.png","type":"image/png","lastModified":1671805879811.355,"size":336247}] 2e:{"title":"Rauzy Fractals","cells":[{"id":0,"value":"# Rauzy Fractals\n\nIn 1982, [Gerard Rauzy](https://fr.wikipedia.org/wiki/G%C3%A9rard_Rauzy) described an intriguing construction in discrete mathematics that leads to beautiful fractal objects, like the following:","pinned":false,"hidden":false,"mode":"md"},{"id":1718,"value":"// These images are all computed in the code below\n// but included here as FileAttachments for speed.\npics[\n [\"0,1,0,2\", \"0,0,1,0,0,2\", \"0,1,2\", \"0,0,1,0,2\", \"0,1,0,2,2\"].indexOf(\n [0, 1].map(s).flat().toString()\n )\n]","pinned":false,"hidden":false,"mode":"ojs"},{"id":1278,"value":"The picture is determined by a substitution ${tex`\\sigma`}, which you can choose from the following:","pinned":false,"hidden":false,"mode":"md"},{"id":605,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":1283,"value":"If you'd like to know what a substitution is and how it might determine that picture, then read on!","pinned":false,"hidden":false,"mode":"md"},{"id":542,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":647,"value":"## An infinite staircase through 3D\n\nWe're now going to use the infinite word ${tex`W`} that's defined as a fixed point of ${tex`\\sigma`} in the previous section to define a discrete subset of ${tex`\\mathbb R^3`} that we call a *staircase*. To do so, let\n${tex.block`\n\\vec{e}_0 = \\begin{bmatrix}1\\\\0\\\\0\\end{bmatrix}, \\:\\: \n\\vec{e}_1 = \\begin{bmatrix}0\\\\1\\\\0\\end{bmatrix}, \\text{ and } \\: \n\\vec{e}_2 = \\begin{bmatrix}0\\\\0\\\\1\\end{bmatrix}\n`}\ndenote the standard unit basis vectors in ${tex`\\mathbb R^3`}. The ${tex`k^{\\text{th}}`} element ${tex`\\vec{s}_k`} of our staircase is then\n${tex.block`\\vec{s}_k = \\sum_{i=0}^k \\vec{e}_{W_i}.`}\nWe might think of this as produced by a sequence of steps, each of which in the North, East, or Up directions.\nThe first 55 points in the staircase look something like so:","pinned":false,"hidden":false,"mode":"md"},{"id":559,"value":"// We delay evaluation so that other cells render quicker per this forum post:\n// https://talk.observablehq.com/t/rauzy-fractals-and-delayed-evaluation/7453\nPromises.delay(delay).then(function () {\n let to_show = gatheredPoints55\n .map(function (A, i) {\n return A.map(function (pt) {\n return create_sphere(pt, 0.2, { color: d3.schemeCategory10[i] });\n });\n })\n .flat();\n\n let n = 55;\n let v = math.column(eigen_basis, 0).toArray().flat();\n to_show.push(\n create_indexedLineSet([\n staircase.slice(0, n),\n [[0, 0, 0], v.map((x) => 34 * x)]\n ])\n );\n return show_x3d(to_show, {\n viewpoint: {\n position: \"40.95545843876255 9.059054609016371 6.442387308246059\",\n orientation:\n \"0.6406938865933923 0.5547318003272355 0.5308332821023546,2.0095131485497375\",\n centerOfRotation: \"0 0 0\"\n }\n });\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":997,"value":"

\n Figure 1: The staircase associated with ${tex`\\sigma`}\n

","pinned":false,"hidden":false,"mode":"html"},{"id":730,"value":"The color of each point indicates the direction of the last step. Note that the same symbol ${tex`0`}, ${tex`1`}, or ${tex`2`} never appears too many times in a row in the word ${tex`W`}; how many times it can occur depends on ${tex`\\sigma`} but there's always a bound. As a result, the staircase never meanders too far away from an average line, as drawn in the figure.","pinned":false,"hidden":false,"mode":"md"},{"id":1760,"value":"## The Rauzy fractal in 3D\n\nWe're now ready to describe how the Rauzy fractal is generated at an intuitive level at least. The idea is to project the points in the staircase along the average line and onto a plane. When we do so, we get the following image of the Rauzy fractal in 3D space:","pinned":false,"hidden":false,"mode":"md"},{"id":342,"value":"// We delay evaluation so that other cells render quicker per this forum post:\n// https://talk.observablehq.com/t/rauzy-fractals-and-delayed-evaluation/7453\nPromises.delay(delay).then(() =>\n show_x3d(\n gatheredProjectedPoints.map((pts, i) =>\n create_pointSet(pts, {\n color: d3.color(d3.schemeCategory10[i]).darker()\n })\n ),\n // Include this to check that we're really looing orthogonally at the object:\n // .concat([create_arrow([0, 0, 0], viewpoint.position)]),\n {\n viewpoint: viewpoint\n }\n )\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1001,"value":"

\n Figure 2: The Rauzy fractal in 3D\n

","pinned":false,"hidden":false,"mode":"html"},{"id":1309,"value":"We obviously need to be more precise about exactly what that \"average line\" is, what plane we project onto, and how we accomplish all that. That's the topic of the next section.","pinned":false,"hidden":false,"mode":"md"},{"id":744,"value":"$32","pinned":false,"hidden":false,"mode":"md"},{"id":767,"value":"{\n let evals = eigenInfo.values.toArray();\n let real_eval = evals[real_pos];\n\n let complex_pos = real_pos == 0 ? 1 : 0;\n let re = evals[complex_pos].re;\n let im = math.abs(evals[complex_pos].im);\n\n return tex`\\begin{aligned}\n \\lambda_0 &\\approx ${d3.format(\"0.3f\")(real_eval)} \\text{ and} \\\\\n \\lambda_{\\pm} &\\approx ${d3.format(\"0.3f\")(re)} \\pm ${d3.format(\"0.3f\")(\n im\n )}\\,i.\n \\end{aligned}`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":795,"value":"The one real eigenvalue tells us there is a one-dimensional subspace that is invariant under the action of ${tex`M`} and that, furthermore, ${tex`M`} is expansive on that subspace, since the eigenvalue is larger than ${tex`1`}. The pair of complex conjugate eigenvalues tell us there is a two dimensional subspace that is also invariant that involves some kind of rotational behavior. The fact that ${tex`|\\lambda_{\\pm}|<1`} tells us that ${tex`M`} is eventually contractive on that subspace, i.e. ${tex`M^k`} is contractive on the subspace for sufficiently large ${tex`k`}.\n\nWe can get more information on the geometric action of ${tex`M`} from the eigenvectors associated with the eigenvalues:","pinned":false,"hidden":false,"mode":"md"},{"id":827,"value":"{\n let v0_tex = math\n .parse(\n math\n .column(eigen_basis, 0)\n .toString()\n .replace(real_regexp, d3.format(\"0.3f\"))\n )\n .toTex();\n let vp_tex = math\n .parse(\n math\n .column(eigen_basis, 1)\n .toString()\n .replace(real_regexp, d3.format(\"0.3f\"))\n )\n .toTex();\n let vm_tex = math\n .parse(\n math\n .column(eigen_basis, 2)\n .toString()\n .replace(real_regexp, d3.format(\"0.3f\"))\n )\n .toTex();\n return tex.block`\\vec{v}_0 \\approx ${v0_tex} \\: \\text{ and } \\: \\vec{v}_{\\pm} \\approx ${vp_tex} \\pm ${vm_tex}i.`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":897,"value":"Specifically, the vector ${tex`\\vec{v}_0`} spans the invariant one-dimensional subspace of ${tex`\\mathbb R^3`} on which ${tex`M`} is expansive. This one-dimensional subspace also happens to be exactly the line drawn in the illustration of the staircase above; that is the direction along which we should project the staircase to obtain the 3D image of the Rauzy fractal. The real and imaginary parts of the complex eigenvector ${tex`\\vec{v}_{\\pm}`} span the invariant two dimensional subspace on which ${tex`M`} is eventually contractive. That is the subspace onto which we should project the staircase to obtain the 3D image of the Rauzy fractal.\n\nTo actually perform the projection, we'll use multiplication by a projection matrix. To form that matrix, we can use the column vectors just described as a change of basis matrix applied to a simple projection onto the ${tex`yz`}-plane. Thus, our projection matrix ${tex`P`} is","pinned":false,"hidden":false,"mode":"md"},{"id":938,"value":"{\n let S = math\n .parse(eigen_basis.toString().replace(real_regexp, d3.format(\"0.3f\")))\n .toTex();\n return tex.block`\\begin{aligned}\n P &\\approx ${S} \\begin{bmatrix} 0 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & 0 & 1\\end{bmatrix} ${S}^{-1} \\\\\n &\\approx ${math\n .parse(P.toString().replace(real_regexp, d3.format(\"0.3f\")))\n .toTex()}.\n\\end{aligned}\n`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1318,"value":"The image shown in figure 2 is obtained by multiplying each of the points generated in the construction of the staircase by this matrix ${tex`P`}.","pinned":false,"hidden":false,"mode":"md"},{"id":1234,"value":"$33","pinned":false,"hidden":false,"mode":"md"},{"id":1570,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":1625,"value":"This process is equivalent to *valuation*, as defined by equation 1.3 in [reference 2](https://observablehq.com/d/28d955eb226a309d#References).","pinned":false,"hidden":false,"mode":"md"},{"id":1766,"value":"## Fractal structure\n\nIt's hard to avoid the observation that the Rauzy fractal appears to display some sort of self-similarity. We'll explore that possibility in a forthcoming notebook.","pinned":false,"hidden":false,"mode":"md"},{"id":1584,"value":"## References\n\nI learned about Rauzy fractals primarily from the following sources:\n\n 1. [Substitutions, Rauzy fractals and tilings](https://www.cambridge.org/core/books/abs/combinatorics-automata-and-number-theory/substitutions-rauzy-fractals-and-tilings/5777998A6799720E4468DD829550B02C)
\n Berthe, Siegel, and Thuswaldner
\n - [PDF](https://www.irif.fr/~berthe/Articles/Chapter5.pdf) An accessible PDF version of that paper\n - [Introduction to Rauzy fractals](http://www.irisa.fr/symbiose/people/asiegel/Canton/Siegel.pdf) a gentler presentation by one of the same authors.\n 2. [Self-affine tiling via substitution dynamical systems and Rauzy fractals](https://msp.org/pjm/2002/206-2/p11.xhtml)
\n by Sirvent and Wang","pinned":false,"hidden":false,"mode":"md"},{"id":1011,"value":"## Code","pinned":false,"hidden":false,"mode":"md"},{"id":1333,"value":"### Substitution and the staircase","pinned":false,"hidden":false,"mode":"md"},{"id":1336,"value":"// Ultimately, everything depends on the substitution function sigma,\n// which is denoted with an s in the code.\n// s, though, is defined in the radio bitton at the top of the notebook.\n// You can see how it's defined here:\n\n[0, 1, 2].map((i) => s(i).toString().replace(/,/g, \"\"))","pinned":true,"hidden":false,"mode":"ojs"},{"id":88,"value":"// W is supposed to be an infinite word.\n// This is just a real long approximation.\n// It's obtained by iteration of s\nW = {\n let i = [0];\n while (i.length < 40000) {\n i = s(i);\n }\n return i;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":92,"value":"// Use the word W to define the stair case.\nstaircase = {\n let basis = [\n [1, 0, 0],\n [0, 1, 0],\n [0, 0, 1]\n ];\n let accumulated = [basis[W[0]]];\n for (let i = 1; i < W.length; i++) {\n let last = accumulated.slice(-1)[0];\n let next = basis[W[i]];\n accumulated.push([last[0] + next[0], last[1] + next[1], last[2] + next[2]]);\n }\n return accumulated;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1343,"value":"### M and P used to define the 3D image","pinned":false,"hidden":false,"mode":"md"},{"id":1348,"value":"// Use the substitution s to define M.\nM = {\n let M = math.zeros(3, 3);\n for (let j = 0; j < 3; j++) {\n let counts = _.countBy(s(j));\n for (let i = 0; i < 3; i++) {\n M.set([i, j], counts[i] || 0);\n }\n }\n return M;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1350,"value":"// Here's what M looks like\nM.toArray()","pinned":true,"hidden":false,"mode":"ojs"},{"id":769,"value":"// Compute the eigenvalues and eigenvectors of M\neigenInfo = math.eigs(M)","pinned":true,"hidden":false,"mode":"ojs"},{"id":250,"value":"// Convert the eigenvectors to a change of basis matrix.\neigen_basis = {\n let real_evec = math.column(eigenInfo.vectors, real_pos);\n let leading_term = real_evec.get([0, 0]);\n real_evec = math.re(math.divide(real_evec, leading_term));\n let complex_pos = eigenInfo.values\n .toArray()\n .map((x) => typeof x)\n .indexOf(\"object\");\n let complex_evec = eigenInfo.vectors.columns()[complex_pos].toArray();\n leading_term = complex_evec[0][0];\n complex_evec = math.divide(complex_evec, leading_term);\n\n let eigen_basis = math.identity(3);\n eigen_basis.subset(math.index([0, 1, 2], 0), real_evec);\n eigen_basis.subset(math.index([0, 1, 2], 1), math.re(complex_evec));\n eigen_basis.subset(math.index([0, 1, 2], 2), math.im(complex_evec));\n\n return eigen_basis;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":774,"value":"// Find the position of the real eigenvalue.\n// Used to help construct the eigen_basis, as well as for\n// the display of the eigenvalues in the exposition\nreal_pos = eigenInfo.values\n .toArray()\n .map((x) => typeof x)\n .indexOf(\"number\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":291,"value":"// Compute the projection matrix.\nP = math.multiply(\n eigen_basis,\n [\n [0, 0, 0],\n [0, 1, 0],\n [0, 0, 1]\n ],\n math.inv(eigen_basis)\n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1366,"value":"## Set up the 3D images","pinned":false,"hidden":false,"mode":"md"},{"id":335,"value":"// Project the points in the long staircase and gather them according to\n// the last step taken. These points are used to generate the 3D fractal.\ngatheredProjectedPoints = {\n let pts = staircase.map((pt) => math.multiply(P, pt).toArray());\n let gatheredPoints = [[], [], []];\n W.forEach((address, idx) => gatheredPoints[address].push(pts[idx]));\n return gatheredPoints;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1204,"value":"// Gather just a few points to display the staircase.\ngatheredPoints55 = {\n let gatheredPoints = [[], [], []];\n W.slice(0, 55).forEach((address, idx) =>\n gatheredPoints[address].push(staircase[idx])\n );\n return gatheredPoints;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":671,"value":"$35","pinned":true,"hidden":false,"mode":"ojs"},{"id":1382,"value":"## A formatter","pinned":false,"hidden":false,"mode":"md"},{"id":925,"value":"// This regular expression is used to identify strings that look like\n// real numbers for formatting in the exposition.\nreal_regexp = /([0-9]*\\.?[0-9]+)/g","pinned":true,"hidden":false,"mode":"ojs"},{"id":1737,"value":"delay = 0","pinned":true,"hidden":false,"mode":"ojs"},{"id":950,"value":"## Imports","pinned":false,"hidden":false,"mode":"md"},{"id":411,"value":"plotx3d_canvas_style","pinned":true,"hidden":false,"mode":"ojs"},{"id":293,"value":"import {\n show_x3d,\n create_pointSet,\n create_indexedLineSet,\n create_sphere,\n create_arrow,\n style as plotx3d_canvas_style\n} from \"@mcmcclur/plotx3d\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":1008,"value":"import { AffineFunction, DigraphIFS } from \"@mcmcclur/digraphifs-class\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":413,"value":"import { cross, normalize, rotation_matrix } from \"@mcmcclur/x3dom-utilities\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":4,"value":"math = require(\"mathjs@11\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":415,"value":"mlm = require(\"ml-matrix@6\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":1720,"value":"pics = Promise.all([\n FileAttachment(\"rauzy0.png\").image(),\n FileAttachment(\"rauzy1.png\").image(),\n FileAttachment(\"rauzy2.png\").image(),\n FileAttachment(\"rauzy3.png\").image(),\n FileAttachment(\"rauzy4.png\").image()\n // \"blank\"\n])","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W36","stdlib":"1","titling":"auto"}