circle_packing

pleat.circle_packing

Discrete circle packing for half-edge graphs.

Implements Collins-Stephenson radii iteration on a triangulated half-edge graph, producing a circle packing whose combinatorics match the input. Two public entry points:

• :func:pack_euclidean: euclidean packing with prescribed boundary radii.
• :func:pack_hyperbolic: maximal circle packing in the Poincaré disk.

Both functions require triangulated, simply-connected (disk topology) input and raise :class:ValueError otherwise. Numerical non-convergence raises :class:ConvergenceError.

The output is a fresh :class:EuclideanPositionHEG (copy of the input by default; pass copy_graph=False to mutate in place) with:

• v['pos']: euclidean center of the circle (complex for hyperbolic output, 2D real array for euclidean output).
• v['radius']: euclidean radius of the circle.

The "natural euclidean" mode (maximal packing transferred to euclidean geometry) is obtained by composition::

G_natural = pack_hyperbolic(G).convert_to_euclidean()

ConvergenceError

Bases: RuntimeError

Raised when the radii iteration fails to reach the requested tolerance.

boundary_angles_from_positions

boundary_angles_from_positions(
    G: EuclideanPositionHEG,
) -> dict[Vertex, float]

Sum the incident-triangle angles at each boundary vertex from v['pos'].

For a flat triangulated tiling these angles automatically satisfy Σ_bdy θ = π(F − 2 V_int) (Gauss-Bonnet), making the result a valid input for :func:pack_euclidean in angle-prescribed mode. Suitable for matching a packing's combinatorics-shape to an existing flat tiling.

Parameters:

Name	Type	Description	Default
G	EuclideanPositionHEG	Half-edge graph with v['pos'] populated on boundary vertices (and their interior triangle neighbors).	required

Returns:

Type	Description
dict[Vertex, float]	Mapping from each boundary vertex to its total incident-triangle angle.

Source code in pleat/circle_packing.py

def boundary_angles_from_positions(G: EuclideanPositionHEG) -> dict[Vertex, float]:
    """Sum the incident-triangle angles at each boundary vertex from ``v['pos']``.

    For a flat triangulated tiling these angles automatically satisfy
    ``Σ_bdy θ = π(F − 2 V_int)`` (Gauss-Bonnet), making the result a valid
    input for :func:`pack_euclidean` in angle-prescribed mode. Suitable for
    matching a packing's combinatorics-shape to an existing flat tiling.

    Args:
        G: Half-edge graph with ``v['pos']`` populated on boundary vertices
            (and their interior triangle neighbors).

    Returns:
        Mapping from each boundary vertex to its total incident-triangle angle.
    """
    result: dict[Vertex, float] = {}
    for v in G.border_vertices():
        p_v = _position_as_2d(v["pos"])
        total = 0.0
        for h in v.outgoing_iter():
            if h.face is None:
                continue
            p_u = _position_as_2d(h.dest["pos"])
            p_w = _position_as_2d(h.nex.dest["pos"])
            vu = p_u - p_v
            vw = p_w - p_v
            cos_a = float(np.dot(vu, vw) / (np.linalg.norm(vu) * np.linalg.norm(vw)))
            cos_a = max(-1.0, min(1.0, cos_a))
            total += float(np.arccos(cos_a))
        result[v] = total
    return result

pack_euclidean

pack_euclidean(
    G: EuclideanPositionHEG,
    boundary_radii: (
        float
        | Mapping[Vertex, float]
        | Callable[[Vertex], float]
        | None
    ) = None,
    *,
    boundary_angles: (
        float
        | Mapping[Vertex, float]
        | Callable[[Vertex], float]
        | str
        | None
    ) = None,
    alpha: Vertex | None = None,
    beta: Vertex | None = None,
    tol: float = 1e-10,
    max_iter: int = 10000,
    copy_graph: bool = True
) -> EuclideanPositionHEG

Compute a euclidean circle packing with prescribed boundary data.

The boundary can be pinned in one of two ways (mutually exclusive):

• boundary_radii: fix each boundary vertex's radius; iterate only interior vertices toward angle sum 2π. This is the classical euclidean boundary-value problem.
• boundary_angles: fix each boundary vertex's total incident-triangle angle; iterate every vertex (boundary and interior). For a flat input tiling pass "from_positions" to read the angles off the input positions. The angles must satisfy Gauss-Bonnet: Σ_bdy θ = π(F − 2 V_int).

If neither is given, defaults to boundary_radii=1.0.

Parameters:

Name	Type	Description	Default
G	EuclideanPositionHEG	Triangulated, simply-connected disk EuclideanPositionHEG.	required
boundary_radii	float | Mapping[Vertex, float] | Callable[[Vertex], float] | None	Per-vertex boundary radii. Accepts a positive scalar (uniform), Mapping[Vertex, float], or callable Vertex -> float.	None
boundary_angles	float | Mapping[Vertex, float] | Callable[[Vertex], float] | str | None	Per-vertex boundary angle sums in (0, 2π). Accepts a scalar (uniform), Mapping, callable, or the string "from_positions" to infer angles from G's positions.	None
alpha	Vertex | None	Optional anchor vertex (placed at origin). Defaults to the interior vertex with maximum graph-distance to the boundary.	None
beta	Vertex | None	Optional second anchor (placed on +x axis through alpha).	None
tol	float	Max permitted angle defect over iterated vertices.	1e-10
max_iter	int	Maximum Collins-Stephenson iterations.	10000
copy_graph	bool	If True (default), the input graph is copied.	True

Returns:

Type	Description
EuclideanPositionHEG	EuclideanPositionHEG with EuclideanGeometry, v['pos'] as 2D real
EuclideanPositionHEG	array centers, and v['radius'] as positive floats.

Raises:

Type	Description
ValueError	input not triangulated/disk; both modes specified; boundary_angles violate Gauss-Bonnet; invalid spec values.
ConvergenceError	iteration did not reach tol within max_iter.

Source code in pleat/circle_packing.py

def pack_euclidean(
    G: EuclideanPositionHEG,
    boundary_radii: float | Mapping[Vertex, float] | Callable[[Vertex], float] | None = None,
    *,
    boundary_angles: float | Mapping[Vertex, float] | Callable[[Vertex], float] | str | None = None,
    alpha: Vertex | None = None,
    beta: Vertex | None = None,
    tol: float = 1e-10,
    max_iter: int = 10_000,
    copy_graph: bool = True,
) -> EuclideanPositionHEG:
    """Compute a euclidean circle packing with prescribed boundary data.

    The boundary can be pinned in one of two ways (mutually exclusive):

    - ``boundary_radii``: fix each boundary vertex's radius; iterate only
      interior vertices toward angle sum 2π. This is the classical
      euclidean boundary-value problem.
    - ``boundary_angles``: fix each boundary vertex's total incident-triangle
      angle; iterate every vertex (boundary and interior). For a flat input
      tiling pass ``"from_positions"`` to read the angles off the input
      positions. The angles must satisfy Gauss-Bonnet:
      ``Σ_bdy θ = π(F − 2 V_int)``.

    If neither is given, defaults to ``boundary_radii=1.0``.

    Args:
        G: Triangulated, simply-connected disk EuclideanPositionHEG.
        boundary_radii: Per-vertex boundary radii. Accepts a positive scalar
            (uniform), Mapping[Vertex, float], or callable Vertex -> float.
        boundary_angles: Per-vertex boundary angle sums in (0, 2π). Accepts
            a scalar (uniform), Mapping, callable, or the string
            ``"from_positions"`` to infer angles from ``G``'s positions.
        alpha: Optional anchor vertex (placed at origin). Defaults to the
            interior vertex with maximum graph-distance to the boundary.
        beta: Optional second anchor (placed on +x axis through alpha).
        tol: Max permitted angle defect over iterated vertices.
        max_iter: Maximum Collins-Stephenson iterations.
        copy_graph: If True (default), the input graph is copied.

    Returns:
        EuclideanPositionHEG with EuclideanGeometry, ``v['pos']`` as 2D real
        array centers, and ``v['radius']`` as positive floats.

    Raises:
        ValueError: input not triangulated/disk; both modes specified;
            boundary_angles violate Gauss-Bonnet; invalid spec values.
        ConvergenceError: iteration did not reach ``tol`` within ``max_iter``.
    """
    _validate_triangulated_disk(G)

    if boundary_radii is not None and boundary_angles is not None:
        raise ValueError("pack_euclidean: specify either `boundary_radii` or `boundary_angles`, not both.")
    if boundary_radii is None and boundary_angles is None:
        boundary_radii = 1.0

    P = G.copy() if copy_graph else G
    boundary_verts = P.border_vertices()
    boundary_set = set(boundary_verts)
    interior_verts = [v for v in P.vertices if v not in boundary_set]

    targets: dict[Vertex, float]
    update_verts: list[Vertex]
    radii: dict[Vertex, float]

    if boundary_angles is not None:
        if isinstance(boundary_angles, str):
            if boundary_angles == "from_positions":
                boundary_angles = boundary_angles_from_positions(P)
            else:
                raise ValueError(f"Unknown boundary_angles string: {boundary_angles!r}. " f"Expected 'from_positions'.")
        angle_spec = _resolve_boundary_angles(boundary_angles, boundary_verts)
        n_faces = len(P.faces)
        n_int = len(interior_verts)
        expected_sum = float(np.pi * (n_faces - 2 * n_int))
        actual_sum = float(sum(angle_spec.values()))
        if not np.isclose(actual_sum, expected_sum, rtol=1e-9, atol=1e-9):
            raise ValueError(
                f"pack_euclidean: boundary_angles sum to {actual_sum:.10f}, but "
                f"Gauss-Bonnet requires Σ θ = π(F − 2·V_int) = π·({n_faces} − 2·{n_int}) "
                f"= {expected_sum:.10f}. Defect = {actual_sum - expected_sum:.3e}."
            )
        targets = {v: 2.0 * np.pi for v in interior_verts}
        targets.update(angle_spec)
        update_verts = list(targets.keys())
        radii = {v: 1.0 for v in P.vertices}
    else:
        radii_spec = _resolve_boundary_radii(boundary_radii, boundary_verts)
        targets = {v: 2.0 * np.pi for v in interior_verts}
        update_verts = interior_verts
        boundary_mean = float(np.mean(list(radii_spec.values()))) if radii_spec else 1.0
        radii = {}
        for v in P.vertices:
            if v in boundary_set:
                radii[v] = radii_spec[v]
            else:
                radii[v] = boundary_mean

    all_vertices = list(P.vertices)
    idx_of = {v: i for i, v in enumerate(all_vertices)}
    r_arr = np.fromiter((radii[v] for v in all_vertices), dtype=float, count=len(all_vertices))
    update_idx, U_arr, W_arr, valid_arr, degrees_arr = _build_flower_arrays(all_vertices, update_verts)
    targets_arr = np.fromiter((targets[v] for v in update_verts), dtype=float, count=len(update_verts))

    _, final_defect, converged = _iterate_with_superstep(
        r_arr,
        update_idx,
        U_arr,
        W_arr,
        valid_arr,
        degrees_arr,
        targets_arr,
        _euclidean_angle_sums_vec,
        _bowers_stephenson_euclidean_update_vec,
        tol,
        max_iter,
    )
    if not converged:
        raise ConvergenceError(
            f"pack_euclidean: did not converge in {max_iter} iterations; "
            f"final max angle defect = {final_defect:.3e}"
        )
    radii = {v: float(r_arr[i]) for v, i in idx_of.items()}

    # Choose anchors
    if alpha is None:
        alpha = _choose_alpha(P)
    if beta is None:
        beta = _choose_beta(alpha)

    # Layout
    pos = _layout_euclidean(P, radii, alpha, beta)

    # Write back onto graph
    for v in P.vertices:
        v["pos"] = pos[v]
        v["radius"] = radii[v]
    P.geometry = EuclideanGeometry
    P.recompute_lengths_and_angles()
    return P

pack_hyperbolic

pack_hyperbolic(
    G: EuclideanPositionHEG,
    boundary_x_radii: (
        float
        | Mapping[Vertex, float]
        | Callable[[Vertex], float]
    ) = 1.0,
    *,
    alpha: Vertex | None = None,
    beta: Vertex | None = None,
    tol: float = 1e-10,
    max_iter: int = 10000,
    copy_graph: bool = True
) -> EuclideanPositionHEG

Compute a hyperbolic circle packing in the Poincaré disk.

Boundary vertices are assigned the prescribed x-radii (where x = 1 - exp(-2h), h being the hyperbolic radius); interior x-radii are iterated via Collins-Stephenson to satisfy hyperbolic angle sum = 2*pi. Setting boundary_x_radii = 1.0 produces a maximal packing (boundary horocycles).

The output stores the euclidean (Poincaré-disk) representation::

v['pos']    — euclidean center of the circle (complex, inside unit disk)
v['radius'] — euclidean radius of the circle (float)

Two hyperbolic disks are hyperbolically tangent iff their euclidean representations are euclidean tangent, so adjacent circles satisfy |c_u − c_v| = r_u + r_v. The intrinsic x-radius of each vertex is recoverable via :func:_x_radius_from_euclidean (and equals 1 for boundary horocycles).

Parameters:

Name	Type	Description	Default
G	EuclideanPositionHEG	Triangulated, simply-connected disk EuclideanPositionHEG.	required
boundary_x_radii	float | Mapping[Vertex, float] | Callable[[Vertex], float]	Per-vertex boundary x-radii in (0, 1]. Accepts a scalar (uniform), Mapping, or callable. Defaults to 0.5; pass 1.0 for the maximal packing.	1.0
alpha	Vertex | None	Optional anchor vertex (placed at Poincaré disk origin). Must be interior. Defaults to the interior vertex with max graph-distance to the boundary.	None
beta	Vertex | None	Optional second anchor (placed on +x axis through alpha).	None
tol	float	Max angle defect over interior vertices.	1e-10
max_iter	int	Max Collins-Stephenson iterations.	10000
copy_graph	bool	If True (default), return a new EHEG; otherwise mutate.	True

Returns:

Type	Description
EuclideanPositionHEG	EuclideanPositionHEG with PoincareDiskModel geometry.

Source code in pleat/circle_packing.py

def pack_hyperbolic(
    G: EuclideanPositionHEG,
    boundary_x_radii: float | Mapping[Vertex, float] | Callable[[Vertex], float] = 1.0,
    *,
    alpha: Vertex | None = None,
    beta: Vertex | None = None,
    tol: float = 1e-10,
    max_iter: int = 10_000,
    copy_graph: bool = True,
) -> EuclideanPositionHEG:
    """Compute a hyperbolic circle packing in the Poincaré disk.

    Boundary vertices are assigned the prescribed x-radii (where x = 1 - exp(-2h),
    h being the hyperbolic radius); interior x-radii are iterated via
    Collins-Stephenson to satisfy hyperbolic angle sum = 2*pi. Setting
    ``boundary_x_radii = 1.0`` produces a maximal packing (boundary horocycles).

    The output stores the **euclidean** (Poincaré-disk) representation::

        v['pos']    — euclidean center of the circle (complex, inside unit disk)
        v['radius'] — euclidean radius of the circle (float)

    Two hyperbolic disks are hyperbolically tangent iff their euclidean
    representations are euclidean tangent, so adjacent circles satisfy
    ``|c_u − c_v| = r_u + r_v``. The intrinsic x-radius of each vertex is
    recoverable via :func:`_x_radius_from_euclidean` (and equals 1 for
    boundary horocycles).

    Args:
        G: Triangulated, simply-connected disk EuclideanPositionHEG.
        boundary_x_radii: Per-vertex boundary x-radii in (0, 1]. Accepts a
            scalar (uniform), Mapping, or callable. Defaults to 0.5; pass
            ``1.0`` for the maximal packing.
        alpha: Optional anchor vertex (placed at Poincaré disk origin).
            Must be interior. Defaults to the interior vertex with max
            graph-distance to the boundary.
        beta: Optional second anchor (placed on +x axis through alpha).
        tol: Max angle defect over interior vertices.
        max_iter: Max Collins-Stephenson iterations.
        copy_graph: If True (default), return a new EHEG; otherwise mutate.

    Returns:
        EuclideanPositionHEG with PoincareDiskModel geometry.
    """
    _validate_triangulated_disk(G)

    P = G.copy() if copy_graph else G

    boundary_verts = P.border_vertices()
    radii_spec = _resolve_boundary_radii(boundary_x_radii, boundary_verts)
    for v, x in radii_spec.items():
        if not (0.0 < x <= 1.0):
            raise ValueError(f"pack_hyperbolic: boundary x-radius must be in (0, 1]; got {x} for vertex {v}.")
    boundary_set = set(boundary_verts)

    # Init: boundary at spec, interior uniformly at 0.5 (mid-range)
    x_radii: dict[Vertex, float] = {}
    for v in P.vertices:
        if v in boundary_set:
            x_radii[v] = radii_spec[v]
        else:
            x_radii[v] = 0.5

    interior_verts = [v for v in P.vertices if v not in boundary_set]
    all_vertices = list(P.vertices)
    idx_of = {v: i for i, v in enumerate(all_vertices)}
    x_arr = np.fromiter((x_radii[v] for v in all_vertices), dtype=float, count=len(all_vertices))
    update_idx, U_arr, W_arr, valid_arr, degrees_arr = _build_flower_arrays(all_vertices, interior_verts)
    targets_arr = np.full(len(interior_verts), 2.0 * np.pi, dtype=float)

    _, final_defect, converged = _iterate_with_superstep(
        x_arr,
        update_idx,
        U_arr,
        W_arr,
        valid_arr,
        degrees_arr,
        targets_arr,
        _hyperbolic_angle_sums_vec,
        _bowers_stephenson_hyperbolic_update_vec,
        tol,
        max_iter,
        max_value=1.0,
    )
    if not converged:
        raise ConvergenceError(
            f"pack_hyperbolic: did not converge in {max_iter} iterations; "
            f"final max angle defect = {final_defect:.3e}"
        )
    x_radii = {v: float(x_arr[i]) for v, i in idx_of.items()}

    if alpha is None:
        alpha = _choose_alpha(P)
    if beta is None:
        beta = _choose_beta(alpha)

    centers, eucl_radii = _layout_hyperbolic(P, x_radii, alpha, beta)

    for v in P.vertices:
        v["pos"] = centers[v]
        v["radius"] = eucl_radii[v]
    P.geometry = PoincareDiskModel
    P.recompute_lengths_and_angles()
    return P