Clean up the fractals directory

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)
+- <https://en.wikipedia.org/wiki/Fractal>

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()

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)