Coverage for src/cvxball/solver.py: 100%
40 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-07-02 08:03 +0000
1"""Convex utilities for computing the minimum enclosing circle/ball.
3Provides two solvers for the smallest enclosing ball problem:
5- :func:`min_circle_cvx`: uses CVXPY to model and then dispatch to a backend
6 solver (default: CLARABEL).
7- :func:`min_circle_clarabel`: bypasses CVXPY and calls the Clarabel solver
8 directly, which removes the CVXPY canonicalisation overhead.
9"""
11from typing import Any
13import clarabel
14import cvxpy as cp
15import numpy as np
16import scipy.sparse as sp
19def min_circle_cvx(points: np.ndarray, **kwargs: Any) -> tuple[float, np.ndarray]:
20 """Compute the smallest enclosing circle for a set of points using convex optimization.
22 This function solves the convex optimization problem to find the minimum radius
23 circle that contains all the given points. It uses a second-order cone constraint
24 to enforce that all points lie within the circle.
26 Args:
27 points: A numpy array of shape (n, d) where n is the number of points
28 and d is the dimension of the space.
29 **kwargs: Additional keyword arguments to pass to the solver.
30 Common options include 'solver' to specify which CVXPY solver to use.
32 Returns:
33 A tuple containing:
34 - The radius of the minimum enclosing circle (float)
35 - The center coordinates of the circle (numpy.ndarray)
37 Example:
38 >>> import numpy as np
39 >>> from cvxball.solver import min_circle_cvx
40 >>> points = np.array([[0, 0], [1, 0], [0, 1]])
41 >>> radius, center = min_circle_cvx(points, solver="CLARABEL")
42 """
43 # cvxpy variable for the radius
44 r = cp.Variable(shape=1, name="Radius")
45 # cvxpy variable for the midpoint
46 x = cp.Variable(points.shape[1], name="Midpoint")
47 objective = cp.Minimize(r)
48 constraints: list[cp.Constraint] = [
49 cp.SOC(
50 # Elementwise broadcast of the scalar radius across all points.
51 # `cp.multiply` (not `*`) avoids CVXPY's deprecated `*`-as-matmul
52 # path, which is ambiguous when n == 1 ((1,) * (1,) -> dot product).
53 cp.multiply(r, np.ones(points.shape[0])), # type: ignore[attr-defined] # cvxpy re-exports atoms via star-import; stubs don't expose them
54 points - x, # Broadcasting handles this automatically
55 axis=1,
56 )
57 ]
59 problem = cp.Problem(objective=objective, constraints=constraints)
60 problem.solve(**kwargs) # type: ignore[no-untyped-call] # cvxpy's Problem.solve is unannotated
62 # Ensure the problem was solved successfully
63 if r.value is None or x.value is None:
64 raise ValueError("Optimization failed to find a solution") # noqa: TRY003
66 return float(r.value[0]), x.value
69def min_circle_clarabel(points: np.ndarray, verbose: bool = False) -> tuple[float, np.ndarray]:
70 """Compute the smallest enclosing circle for a set of points using Clarabel directly.
72 This function solves the same convex optimisation problem as
73 :func:`min_circle_cvx` but bypasses CVXPY and calls the Clarabel solver
74 directly. The problem is assembled in Clarabel's standard form::
76 minimise (1/2) z' P z + q' z
77 subject to A z + s = b, s ∈ K
79 where the decision vector is ``z = [r, x₁, …, x_d]`` (radius followed by
80 the d centre coordinates), the objective is to minimise *r* (so ``P = 0``,
81 ``q = e₀``), and the feasible set is a product of *n* second-order cones.
83 For each point ``p_i`` we require ``[r, p_i - x] in Q^{d+1}``, which gives
84 one SOC block of dimension ``d + 1`` per point.
86 Args:
87 points: A numpy array of shape ``(n, d)`` where *n* is the number of
88 points and *d* is the ambient dimension.
89 verbose: If ``True``, print Clarabel's iteration log. Defaults to
90 ``False``.
92 Returns:
93 A tuple ``(radius, center)`` where *radius* is the optimal enclosing
94 radius (float) and *center* is a numpy array of shape ``(d,)``.
96 Raises:
97 ValueError: If Clarabel does not return a ``Solved`` status.
99 Example:
100 >>> import numpy as np
101 >>> from cvxball.solver import min_circle_clarabel
102 >>> points = np.array([[0, 0], [1, 0], [0, 1]])
103 >>> radius, center = min_circle_clarabel(points)
104 """
105 n, d = points.shape
106 n_vars = 1 + d # decision vector: [r, x_1, ..., x_d]
108 # --- Objective: minimise r -----------------------------------------------
109 p_mat = sp.csc_matrix((n_vars, n_vars))
110 q = np.zeros(n_vars)
111 q[0] = 1.0
113 # --- Constraints: one SOC block of size (d+1) per point ------------------
114 # We need b - a_mat @ z = s where s in K.
115 # For point i the desired slack is s = [r, p_i - x], so:
116 # row i*(d+1) : b = 0, a_mat col 0 = -1 (gives s_0 = r)
117 # row i*(d+1)+j : b = p_i[j], a_mat col j = +1 (gives s_j = p_ij - x_j)
118 total_rows = n * (d + 1)
120 # Entries for the r column (column 0): -1 at each block's first row
121 r_rows = np.arange(n) * (d + 1)
123 # Entries for the x columns (columns 1..d): +1 at each block's inner rows
124 x_row_offsets = np.arange(n)[:, None] * (d + 1) + np.arange(1, d + 1)[None, :] # (n, d)
125 x_rows = x_row_offsets.ravel()
126 x_cols = np.tile(np.arange(1, d + 1), n)
128 all_rows = np.concatenate([r_rows, x_rows])
129 all_cols = np.concatenate([np.zeros(n, dtype=np.intp), x_cols])
130 all_vals = np.concatenate([-np.ones(n), np.ones(n * d)])
132 a_mat = sp.csc_matrix((all_vals, (all_rows, all_cols)), shape=(total_rows, n_vars))
134 b = np.zeros(total_rows)
135 b[x_rows] = points.ravel()
137 # --- Cones: n SOC cones each of dimension (d+1) --------------------------
138 cones = [clarabel.SecondOrderConeT(d + 1) for _ in range(n)] # ty: ignore[unresolved-attribute]
140 # --- Solve ---------------------------------------------------------------
141 settings = clarabel.DefaultSettings.default() # ty: ignore[unresolved-attribute]
142 settings.verbose = verbose
144 solver = clarabel.DefaultSolver(p_mat, q, a_mat, b, cones, settings) # ty: ignore[unresolved-attribute]
145 solution = solver.solve()
147 if solution.status != clarabel.SolverStatus.Solved: # ty: ignore[unresolved-attribute]
148 raise ValueError(f"Clarabel did not converge: status = {solution.status}") # noqa: TRY003
150 return float(solution.x[0]), np.asarray(solution.x[1:])