diff --git a/fractals/README.md b/fractals/README.md new file mode 100644 index 000000000000..ddb196192f91 --- /dev/null +++ b/fractals/README.md @@ -0,0 +1,44 @@ +# Fractals + +A fractal is a geometric figure that is *self-similar*: zooming into a piece of +it reveals a copy of the whole. Fractals show up in mathematics, physics, +computer graphics and even biology (coastlines, ferns, snowflakes). + +This directory collects small, self-contained fractal generators. They fall +into two groups: + +- **Visual demos** that open a window (via `turtle`) or produce an image + (via `matplotlib`/`PIL`). Run these directly to see the picture. +- **Pure-computation** generators whose output can be checked with `doctest`, + so they run in CI without a display. + +## Contents + +| File | Fractal | Output | Notes | +| ---- | ------- | ------ | ----- | +| [`barnsley_fern.py`](barnsley_fern.py) | Barnsley fern | matplotlib (optional) | Iterated function system; deterministic with a seed | +| [`julia_sets.py`](julia_sets.py) | Julia sets | matplotlib | Complex-plane escape-time fractal | +| [`koch_snowflake.py`](koch_snowflake.py) | Koch snowflake | matplotlib | Line-segment subdivision | +| [`mandelbrot.py`](mandelbrot.py) | Mandelbrot set | PIL image | Complex-plane escape-time fractal | +| [`sierpinski_carpet.py`](sierpinski_carpet.py) | Sierpinski carpet | text | Integer arithmetic, fully doctested | +| [`sierpinski_triangle.py`](sierpinski_triangle.py) | Sierpinski triangle | turtle | Recursive midpoint subdivision | +| [`vicsek.py`](vicsek.py) | Vicsek fractal | turtle | Recursive cross pattern | + +## Running + +```bash +# text fractal – prints to the terminal +python fractals/sierpinski_carpet.py + +# image fractal – opens a matplotlib window (needs matplotlib) +python fractals/barnsley_fern.py + +# turtle fractal – opens a drawing window (needs a display) +python fractals/vicsek.py +``` + +## Further reading + +- Benoit B. Mandelbrot, *The Fractal Geometry of Nature* (1982) +- Michael Barnsley, *Fractals Everywhere* (1988) +- diff --git a/fractals/barnsley_fern.py b/fractals/barnsley_fern.py new file mode 100644 index 000000000000..4497873cf946 --- /dev/null +++ b/fractals/barnsley_fern.py @@ -0,0 +1,147 @@ +""" +The Barnsley fern is a fractal that resembles the black spleenwort fern. It was +described by the British mathematician Michael Barnsley in his 1988 book +*Fractals Everywhere* and is a classic example of an iterated function system +(IFS). + +An IFS builds a fractal by repeatedly applying a small set of affine +transformations, each chosen at random with a fixed probability. Starting from +the point ``(0, 0)`` the fern uses four transformations: + +=============== =========================================== ============ +Transformation Effect Probability +=============== =========================================== ============ +Stem collapse onto the y-axis 1% +Successive leaf the main self-similar copy of the fern 85% +Left leaflet a smaller rotated/reflected copy 7% +Right leaflet another smaller rotated/reflected copy 7% +=============== =========================================== ============ + +Because the whole picture is produced by chance the doctests below seed Python's +random generator so that the results are reproducible. Plotting the points with +matplotlib is optional and only happens when the module is run directly. + +Reference: https://en.wikipedia.org/wiki/Barnsley_fern +""" + +import random + +# Each row is (a, b, c, d, e, f) for the affine map +# x' = a*x + b*y + e +# y' = c*x + d*y + f +# and the running cumulative probabilities used to pick a transformation. +TRANSFORMATIONS: tuple[tuple[float, float, float, float, float, float], ...] = ( + (0.00, 0.00, 0.00, 0.16, 0.00, 0.00), # stem + (0.85, 0.04, -0.04, 0.85, 0.00, 1.60), # successive smaller leaflets + (0.20, -0.26, 0.23, 0.22, 0.00, 1.60), # left-hand leaflet + (-0.15, 0.28, 0.26, 0.24, 0.00, 0.44), # right-hand leaflet +) +CUMULATIVE_PROBABILITIES: tuple[float, ...] = (0.01, 0.86, 0.93, 1.00) + + +def transform(point: tuple[float, float], index: int) -> tuple[float, float]: + """ + Apply the affine transformation ``index`` to ``point`` and return the image. + + >>> transform((0.0, 0.0), 0) + (0.0, 0.0) + >>> transform((1.0, 1.0), 1) + (0.89, 2.41) + >>> transform((2.0, 3.0), 3) + (0.54, 1.68) + >>> transform((0.0, 0.0), 4) + Traceback (most recent call last): + ... + IndexError: index must be in range 0..3, got 4 + """ + if not 0 <= index < len(TRANSFORMATIONS): + msg = f"index must be in range 0..3, got {index}" + raise IndexError(msg) + a, b, c, d, e, f = TRANSFORMATIONS[index] + x, y = point + new_x = round(a * x + b * y + e, 12) + new_y = round(c * x + d * y + f, 12) + return (new_x, new_y) + + +def choose_transformation(sample: float) -> int: + """ + Map a value ``sample`` from ``[0, 1)`` to a transformation index using the + cumulative probabilities of the fern. + + >>> choose_transformation(0.0) + 0 + >>> choose_transformation(0.5) + 1 + >>> choose_transformation(0.9) + 2 + >>> choose_transformation(0.97) + 3 + """ + for index, threshold in enumerate(CUMULATIVE_PROBABILITIES): + if sample < threshold: + return index + return len(CUMULATIVE_PROBABILITIES) - 1 + + +def generate_fern( + iterations: int, seed: int | None = None +) -> list[tuple[float, float]]: + """ + Generate ``iterations`` points of the Barnsley fern, starting at ``(0, 0)``. + + Passing a ``seed`` makes the (otherwise random) output reproducible, which is + what keeps the doctests deterministic. + + >>> points = generate_fern(5, seed=0) + >>> len(points) + 5 + >>> points[0] + (0.0, 0.0) + >>> points # doctest: +NORMALIZE_WHITESPACE + [(0.0, 0.0), (0.0, 1.6), (0.064, 2.96), + (0.1728, 4.11344), (0.3114176, 5.089512)] + + Every fern point lives inside the well known bounding box. + + >>> cloud = generate_fern(2000, seed=42) + >>> all(-2.182 <= x <= 2.6558 for x, _ in cloud) + True + >>> all(0.0 <= y <= 9.9984 for _, y in cloud) + True + >>> generate_fern(0) + Traceback (most recent call last): + ... + ValueError: iterations must be positive, got 0 + """ + if iterations <= 0: + msg = f"iterations must be positive, got {iterations}" + raise ValueError(msg) + rng = random.Random(seed) + point = (0.0, 0.0) + points = [point] + for _ in range(iterations - 1): + index = choose_transformation(rng.random()) + point = transform(point, index) + points.append(point) + return points + + +if __name__ == "__main__": + import doctest + + doctest.testmod() + + try: + import matplotlib.pyplot as plt + except ImportError: + print("matplotlib is required to plot the fern (pip install matplotlib).") + else: + fern_points = generate_fern(100_000, seed=0) + xs = [x for x, _ in fern_points] + ys = [y for _, y in fern_points] + plt.figure(figsize=(4, 8)) + plt.scatter(xs, ys, s=0.2, color="forestgreen") + plt.axis("off") + plt.title("Barnsley fern") + plt.show() diff --git a/fractals/sierpinski_carpet.py b/fractals/sierpinski_carpet.py new file mode 100644 index 000000000000..f3b2a0a2e6b6 --- /dev/null +++ b/fractals/sierpinski_carpet.py @@ -0,0 +1,130 @@ +""" +The Sierpinski carpet is a plane fractal first described by Wacław Sierpiński +in 1916. It is a two-dimensional generalisation of the Cantor set and a close +relative of the Sierpinski triangle. + +Construction + Start from a filled square. Divide it into a 3x3 grid of nine equal + sub-squares and remove the central one. Then apply the same procedure + recursively to each of the eight remaining sub-squares, forever. + +A convenient way to decide whether a single cell of the ``3**n x 3**n`` grid is +filled (part of the carpet) or empty (a hole) is to look at the base-3 digits +of its row and column indices: the cell is a hole if and only if, at some +level, both the row digit and the column digit are equal to ``1`` (the centre +of that 3x3 block). + +This module builds the carpet purely with integer arithmetic, so every +function is deterministic and can be verified with doctests -- no plotting or +turtle graphics required. + +Reference: https://en.wikipedia.org/wiki/Sierpi%C5%84ski_carpet +""" + + +def is_filled(row: int, col: int) -> bool: + """ + Return ``True`` when the cell at (``row``, ``col``) belongs to the carpet + and ``False`` when it falls inside one of the removed central squares. + + The result is independent of the fractal depth: a cell is a hole as soon as + any pair of matching base-3 digits equals ``(1, 1)``. + + >>> is_filled(0, 0) + True + >>> is_filled(1, 1) # the very first central square is removed + False + >>> is_filled(4, 4) # centre of the centre block -> still a hole + False + >>> is_filled(0, 4) + True + >>> [is_filled(1, col) for col in range(3)] + [True, False, True] + + Negative coordinates make no sense for a grid index. + + >>> is_filled(-1, 0) + Traceback (most recent call last): + ... + ValueError: row and col must be non-negative, got (-1, 0) + """ + if row < 0 or col < 0: + msg = f"row and col must be non-negative, got ({row}, {col})" + raise ValueError(msg) + while row > 0 or col > 0: + if row % 3 == 1 and col % 3 == 1: + return False + row //= 3 + col //= 3 + return True + + +def generate_carpet(depth: int, filled: str = "#", hole: str = " ") -> list[str]: + """ + Build the Sierpinski carpet of the given ``depth`` as a list of strings. + + A depth of ``0`` is a single filled cell; each extra level multiplies the + side length by three. + + >>> generate_carpet(0) + ['#'] + >>> for line in generate_carpet(1): + ... print(line) + ### + # # + ### + >>> for line in generate_carpet(2, filled="X", hole="."): + ... print(line) + XXXXXXXXX + X.XX.XX.X + XXXXXXXXX + XXX...XXX + X.X...X.X + XXX...XXX + XXXXXXXXX + X.XX.XX.X + XXXXXXXXX + >>> generate_carpet(-1) + Traceback (most recent call last): + ... + ValueError: depth must be non-negative, got -1 + """ + if depth < 0: + msg = f"depth must be non-negative, got {depth}" + raise ValueError(msg) + size = 3**depth + return [ + "".join(filled if is_filled(row, col) else hole for col in range(size)) + for row in range(size) + ] + + +def count_filled_cells(depth: int) -> int: + """ + Return how many cells are filled in a carpet of the given ``depth``. + + Each level keeps eight of the nine sub-squares, so the count is ``8**depth``. + Verifying this closed form against a brute-force scan is a nice sanity check. + + >>> [count_filled_cells(depth) for depth in range(4)] + [1, 8, 64, 512] + >>> all( + ... count_filled_cells(depth) + ... == sum(line.count("#") for line in generate_carpet(depth)) + ... for depth in range(4) + ... ) + True + """ + if depth < 0: + msg = f"depth must be non-negative, got {depth}" + raise ValueError(msg) + return 8**depth + + +if __name__ == "__main__": + import doctest + + doctest.testmod() + + for carpet_line in generate_carpet(3): + print(carpet_line)