API Reference¶

def ball_with_counts(
    points: np.ndarray,
    pivot_rule: PivotRule = "heuristic",
    dynamic_qr: bool = True,
    verbose: bool = False,
) -> Ball:
    """Solve the smallest enclosing ball by the pivoting method, reporting the work done.

    Args:
        points: A ``(n, d)`` array with ``n >= 1``.
        pivot_rule: ``"heuristic"`` (the paper's own code) or ``"bland"`` (the
            rule Theorem 1 proves terminating).
        dynamic_qr: Carry the factorisation across pivots, repairing it in
            ``O(d r)`` as section 4 does. ``False`` refactorises from scratch at
            every pivot, in ``O(d r^2)`` -- the same answers, and the baseline the
            data structure is worth measuring against.
        verbose: If ``True``, print the phase, support size and radius per turn.

    Returns:
        The :class:`Ball`, including the support set and the pivot counts.

    Raises:
        ValueError: If ``points`` is not a finite, non-empty ``(n, d)`` array (see
            :func:`cvxball.solver._validate`), if ``pivot_rule`` is not one of the
            two rules, or if the iteration limit is reached -- which, termination
            being a theorem, means numerical trouble rather than a hard instance.
    """
    if pivot_rule not in ("heuristic", "bland"):
        raise ValueError(f"pivot_rule must be 'heuristic' or 'bland', got {pivot_rule!r}")  # noqa: TRY003

    pts = _validate(points)
    n = pts.shape[0]

    # Recentre, so that every tolerance below is governed by the extent of the
    # cloud rather than by its distance from an arbitrary origin, and difference
    # before squaring for the same reason the shipped solver does.
    shift = pts.mean(axis=0)
    pts = pts - shift
    extent = float(np.sqrt(np.einsum("ij,ij->i", pts, pts).max()))
    if extent == 0.0:
        # Every point is the same point; the ball is that point, radius zero.
        return Ball(0.0, shift.copy(), np.zeros(1, dtype=np.intp), 0, 0, 0)
    margin = _AFF_RTOL * extent

    # Initialisation, from Fig. 2: c is any point of S, and T the single point of
    # S farthest from it -- which is what makes B(c, T) enclose S to begin with.
    centre = pts[0].copy()
    offsets = pts - centre
    frame = _Frame(pts, int(np.argmax(np.einsum("ij,ij->i", offsets, offsets))), dynamic_qr)
    support = frame.support
    # Whether c is known to lie in aff(T). It is not, initially: c is one point of
    # the cloud and aff(T) is a different one. It becomes true exactly when a walk
    # runs to completion, since the centre is then the projection onto that hull.
    on_hull = False

    drops = insertions = 0
    limit = _MAX_ITER_PER_POINT * (n + 1)
    for iteration in range(limit):
        if on_hull:
            leaving = _leaving_point(frame.coefficients(centre), support, pivot_rule)
            if leaving is None:
                # c is in conv(T): by Lemma 1 this ball is SEB(S).
                radius = float(np.linalg.norm(pts[support] - centre, axis=1).max())
                if verbose:
                    print(f"[{iteration:4d}] optimal   support={len(support):3d} radius={radius:.12g}")
                return Ball(radius, centre + shift, np.array(support, dtype=np.intp), iteration, drops, insertions)
            if verbose:
                print(f"[{iteration:4d}] drop      support={len(support):3d} point={support[leaving]}")
            frame.remove(leaving)
            drops += 1
            on_hull = False
            continue

        # --- Walking phase ---------------------------------------------------
        # cc(T) is the orthogonal projection of c onto aff(T), by Lemma 3(i).
        direction = frame.direction_to_circumcentre(centre)
        direction_sq = float(direction @ direction)
        if np.sqrt(direction_sq) <= margin:
            # Already on the hull -- reached whenever the support has grown to
            # d + 1 points, whose affine hull is the whole space.
            on_hull = True
            continue

        offsets = pts - centre
        squared = np.einsum("ij,ij->i", offsets, offsets)
        radius_sq = float(squared[support].max())

        # Lemma 4: p can stop the walk only if <p - c, e> < <e, e>, i.e. only if
        # this denominator is positive. Requiring it to clear `margin * ||e||`
        # rather than zero is the stability threshold, and it is exactly the
        # right test: e is orthogonal to aff(T), so denom / ||e|| is p's distance
        # from aff(T) measured along e, and thus a lower bound on its full
        # distance from aff(T). Clearing the margin here therefore certifies that
        # T stays affinely independent when p joins it -- and the points it
        # excludes have a near-zero denominator, hence an enormous stopping
        # fraction, so they were never going to be the minimiser anyway.
        denominator = direction_sq - offsets @ direction
        stoppers = denominator > margin * np.sqrt(direction_sq)
        stoppers[support] = False

        fractions = np.full(n, np.inf)
        np.divide(
            np.maximum(radius_sq - squared, 0.0),
            2.0 * denominator,
            out=fractions,
            where=stoppers,
        )
        # Walk to the nearest stopper the support can actually take. A candidate
        # that fails `admits` is numerically inside aff(T): taking it would break
        # the affine independence Fig. 2's invariant rests on, so it is passed
        # over and the next-nearest considered. Usually the first one is fine.
        while True:
            shortest = float(fractions.min())
            if shortest >= 1.0:
                # Nothing stops the walk: the centre reaches cc(T) and lands on the hull.
                centre = centre + direction
                on_hull = True
                if verbose:
                    radius = float(np.linalg.norm(pts[support] - centre, axis=1).max())
                    print(f"[{iteration:4d}] arrive    support={len(support):3d} radius={radius:.12g}")
                break

            entering = _entering_point(fractions, shortest, pts, frame, pivot_rule)
            if not frame.admits(entering):
                fractions[entering] = np.inf
                continue

            centre = centre + shortest * direction
            frame.insert(entering)
            insertions += 1
            if verbose:
                radius = float(np.linalg.norm(pts[support] - centre, axis=1).max())
                print(f"[{iteration:4d}] insert    support={len(support):3d} radius={radius:.12g} step={shortest:.6g}")
            break

    raise ValueError(f"pivoting method did not converge in {limit} iterations")  # noqa: TRY003

def min_circle_active_set(
    points: np.ndarray,
    verbose: bool = False,
    maintain: bool | None = None,
) -> tuple[float, np.ndarray]:
    """Compute the smallest enclosing circle with an active-set method.

    An active-set QP method on the dual of the enclosing-ball problem, in place of
    handing a cone program to a conic solver.  It
    maintains a *support set* of points held on the ball's boundary and repeatedly

    1. centres the ball on that support set by solving one small linear system
       (:func:`_face_weights`),
    2. shrinks the support when it cannot hold — either because a weight would turn
       negative, or because the set has become affinely dependent
       (:func:`_affine_null_space`) — moving as far as non-negativity allows, and
    3. adds the farthest point that is still outside the ball.

    Each subproblem is a ``k x k`` solve with ``k <= d``, so the cost per iteration
    is driven by the dimension rather than by the number of points, and the method
    stops at an exact vertex of the dual feasible set instead of at an
    interior-point tolerance.

    Args:
        points: A numpy array of shape ``(n, d)`` where *n* is the number of
                points and *d* is the ambient dimension.
        verbose: If ``True``, print the support size and radius per iteration.
                 Defaults to ``False``.
        maintain: Whether to carry the support's factorisation across iterations
                 and repair it, rather than rebuilding it each time.  ``None``,
                 the default, decides on the ambient dimension: repairing costs
                 ``O(d r)`` against ``O(d r^2)`` to rebuild, which is decisive by
                 ``d = 250`` and a net loss below ``d = 100`` where SciPy's
                 per-update overhead exceeds what it saves (see
                 :data:`_MAINTAIN_MIN_DIM`).  Both settings compute the same ball;
                 pass one explicitly only to measure the difference.

    Returns:
        A tuple ``(radius, center)`` where *radius* is the optimal enclosing
        radius (float) and *center* is a numpy array of shape ``(d,)``.

    Raises:
        ValueError: If ``points`` is not a finite, non-empty ``(n, d)`` array
                    (see :func:`_validate`), or if the iteration limit is reached
                    without the optimality certificate holding.

    Example:
        Three points forming a right triangle, whose smallest enclosing circle is
        the one on its hypotenuse.  Both values are pinned to full precision here,
        where the cone program in ``experiments/clarabel_ball.py`` can only pin its
        centre to three decimals — this method stops at an exact vertex of the
        dual feasible set, so on an input whose answer is exactly representable it
        returns that answer bit-for-bit.

        >>> import numpy as np
        >>> from cvxball import min_circle_active_set
        >>> points = np.array([[0, 0], [1, 0], [0, 1]])
        >>> radius, center = min_circle_active_set(points)
        >>> radius == 2**0.5 / 2
        True
        >>> center
        array([0.5, 0.5])
    """
    pts = _validate(points)
    n = pts.shape[0]

    # Work in coordinates centred on the cloud, undoing the shift on the way out.
    # Every quantity below -- the subproblem's Gram matrix, the squared distances, the
    # rounding floor -- is then governed by the *extent* of the cloud rather than by
    # its distance from an arbitrary origin.  Without this, a cloud of extent 1 sitting
    # a million units out gets a rounding floor of the same order as its own radius,
    # and the optimality test below then accepts a ball that is visibly too small.
    shift = pts.mean(axis=0)
    pts = pts - shift

    # Recentring fixes the origin but not the magnitude, and this method squares
    # everything it touches: the Gram matrix, the squared distances, the noise floor.
    # Squaring halves the usable exponent range, so a cloud of extent 1e-160 -- whose
    # coordinates are ordinary doubles -- has a Gram matrix of order 1e-320, deep in
    # the subnormals, and the solve for the circumcentre then overflows to infinity.
    # The far end fails too, and worse: past 1e+154 the squares saturate and the
    # method returns a radius of zero rather than raising.
    #
    # So normalise the extent to order one.  The factor is a *power of two*, which
    # makes the rescaling exact in binary floating point: not one bit of the answer
    # moves, so the method keeps returning representable answers bit-for-bit (see the
    # example above), and the whole iteration runs where doubles are well behaved.
    # Measuring the extent as a largest coordinate rather than a largest norm is what
    # keeps the measurement itself in range -- a norm would already have squared.
    largest = float(np.abs(pts).max())
    if largest == 0.0:
        # Every point is the same point, so the ball is that point with radius zero.
        return 0.0, shift
    exponent = int(np.frexp(largest)[1])
    pts = np.ldexp(pts, -exponent)

    sq_norms: np.ndarray = np.einsum("ij,ij->i", pts, pts)

    # Warm start with one point and whatever sits farthest from it: that puts both
    # ends of a near-diameter into the support straight away, which is usually
    # where the optimum keeps them.
    seed = int(np.argmax(_sq_dist(pts, pts[0])))
    # Rounding floor of a squared distance, set by the magnitude of the coordinates
    # that go into it rather than by any fixed constant.
    noise_floor = _FEAS_NOISE * float(np.finfo(np.float64).eps) * float(sq_norms.max())

    if maintain is None:
        maintain = pts.shape[1] >= _MAINTAIN_MIN_DIM
    support = _MaintainedFace(pts, seed) if maintain else _RebuiltFace(pts, seed)
    free = support.support
    weights = np.ones(1)

    iteration_limit = _MAX_ITER_PER_POINT * (n + 1)
    for iteration in range(iteration_limit):
        face = support.face
        null_space = support.null_space()

        if null_space.size:
            # Affinely dependent support: g is linear on the null space, so follow
            # steepest descent there until some weight is driven down to zero.
            centre = face.T @ weights
            # grad_i = ||x||^2 - ||p_i - x||^2, less the constant ||x||^2: null
            # directions sum to zero, so dropping it leaves every projection alone
            # while avoiding the cancellation of ||x||^2 against ||p_i||^2.
            gradient = -_sq_dist(face, centre)
            descent = -(null_space @ (null_space.T @ gradient))
            # Rescale to a unit-max direction.  The gradient carries units of length^2,
            # so comparing it against the dimensionless weight tolerance below would
            # make the whole method scale-dependent: on a cloud of extent 1e-20 no
            # component would ever look binding.  Only the direction matters here --
            # `alpha` cancels any positive factor -- so normalising costs nothing.
            largest = float(np.abs(descent).max())
            step = descent / largest if largest > 0.0 else descent
        else:
            target = support.circumcentre_weights()
            if target.min() >= -_DROP_TOL:
                # The subproblem solution is feasible: take it, then test the KKT
                # certificate and, if it fails, free the most violated point.
                weights = np.maximum(target, 0.0)
                weights /= weights.sum()
                centre = face.T @ weights
                dist_sq = _sq_dist(pts, centre)
                radius_sq = float(dist_sq[free].max())
                radius = float(np.sqrt(max(radius_sq, 0.0)))
                if verbose:
                    print(f"[{iteration:4d}] support={len(free):3d} radius={radius:.12g}")

                worst = int(np.argmax(dist_sq))
                if dist_sq[worst] <= radius_sq * (1.0 + _FEAS_RTOL) + noise_floor:
                    # Undo the exact power-of-two normalisation, then the shift.
                    return float(np.ldexp(radius, exponent)), np.ldexp(centre, exponent) + shift

                keep = weights > _DROP_TOL
                _shrink(support, keep)
                support.insert(worst)
                weights = np.append(weights[keep], 0.0)
                continue

            step = target - weights

        # Longest step along `step` that keeps every weight non-negative.  A step
        # of zero is impossible: the newly freed point always has step > 0, so the
        # support strictly shrinks here and `g` strictly decreases.
        binding = step < -_DROP_TOL
        ratios = np.where(binding, -weights / np.where(binding, step, -1.0), np.inf)
        alpha = float(ratios.min())
        if not np.isfinite(alpha):
            # Nothing blocks the step, so the support cannot shrink.  Unreachable in
            # exact arithmetic; bail out rather than propagate a non-finite weight.
            raise ValueError("active-set method stalled: no support point blocks the step")  # noqa: TRY003
        weights = np.maximum(weights + alpha * step, 0.0)

        keep = weights > _DROP_TOL
        _shrink(support, keep)
        weights = weights[keep]
        weights /= weights.sum()
        if verbose:
            print(f"[{iteration:4d}] support={len(free):3d} drop step={alpha:.6g}")

    raise ValueError(f"active-set method did not converge in {iteration_limit} iterations")  # noqa: TRY003

def min_circle_fgk(points: np.ndarray, verbose: bool = False) -> tuple[float, np.ndarray]:
    """Compute the smallest enclosing ball, in this package's solver signature.

    Args:
        points: A numpy array of shape ``(n, d)`` where *n* is the number of
                points and *d* is the ambient dimension.
        verbose: If ``True``, print one line per pivot step. Defaults to ``False``.

    Returns:
        A tuple ``(radius, center)``, matching :func:`cvxball.min_circle_active_set`.

    Raises:
        ValueError: If ``points`` is not a finite, non-empty ``(n, d)`` array.

    Example:
        The right triangle whose smallest enclosing circle is the one on its
        hypotenuse.

        The values are rounded here, as they are for the cone program, and the
        reason is worth stating because it is easy to assume otherwise: this
        method terminates at an exact *combinatorial* configuration -- the
        support set ``{(1, 0), (0, 1)}`` -- but the centre it reports is not the
        exact circumcentre of that set. It is the running sum of the walks that
        got there, so it lands a few ulp out (four, on this input). The shipped
        active-set method solves afresh for the centre of its final support and
        so returns ``sqrt(2) / 2`` bit-for-bit; the difference is one of
        arithmetic, not of which ball the two methods identify.

        >>> import numpy as np
        >>> from cvxball import min_circle_fgk
        >>> radius, center = min_circle_fgk(np.array([[0, 0], [1, 0], [0, 1]]))
        >>> round(radius, 12)
        0.707106781187
        >>> np.round(center, 12)
        array([0.5, 0.5])
    """
    ball = ball_with_counts(points, verbose=verbose)
    if verbose:
        print(f"fgk: iterations={ball.iterations} drops={ball.drops} insertions={ball.insertions}")
    return ball.radius, ball.centre