Coverage for src/cvxball/solver.py: 100%
43 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-26 05: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 _build_soc_program(
70 points: np.ndarray,
71) -> tuple[sp.csc_matrix, np.ndarray, sp.csc_matrix, np.ndarray, list[Any]]:
72 """Assemble the Clarabel second-order-cone program for the enclosing ball.
74 The problem is written 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.
90 Returns:
91 A tuple ``(p_mat, q, a_mat, b, cones)`` of the objective quadratic
92 ``P``, the objective linear term ``q``, the constraint matrix ``A``,
93 the constraint right-hand side ``b``, and the list of *n* second-order
94 cones — the exact positional arguments Clarabel's ``DefaultSolver``
95 expects.
96 """
97 n, d = points.shape
98 n_vars = 1 + d # decision vector: [r, x_1, ..., x_d]
100 # --- Objective: minimise r -----------------------------------------------
101 p_mat = sp.csc_matrix((n_vars, n_vars))
102 q = np.zeros(n_vars)
103 q[0] = 1.0
105 # --- Constraints: one SOC block of size (d+1) per point ------------------
106 # We need b - a_mat @ z = s where s in K.
107 # For point i the desired slack is s = [r, p_i - x], so:
108 # row i*(d+1) : b = 0, a_mat col 0 = -1 (gives s_0 = r)
109 # row i*(d+1)+j : b = p_i[j], a_mat col j = +1 (gives s_j = p_ij - x_j)
110 total_rows = n * (d + 1)
112 # Entries for the r column (column 0): -1 at each block's first row
113 r_rows = np.arange(n) * (d + 1)
115 # Entries for the x columns (columns 1..d): +1 at each block's inner rows
116 x_row_offsets = np.arange(n)[:, None] * (d + 1) + np.arange(1, d + 1)[None, :] # (n, d)
117 x_rows = x_row_offsets.ravel()
118 x_cols = np.tile(np.arange(1, d + 1), n)
120 all_rows = np.concatenate([r_rows, x_rows])
121 all_cols = np.concatenate([np.zeros(n, dtype=np.intp), x_cols])
122 all_vals = np.concatenate([-np.ones(n), np.ones(n * d)])
124 a_mat = sp.csc_matrix((all_vals, (all_rows, all_cols)), shape=(total_rows, n_vars))
126 b = np.zeros(total_rows)
127 b[x_rows] = points.ravel()
129 # --- Cones: n SOC cones each of dimension (d+1) --------------------------
130 cones = [clarabel.SecondOrderConeT(d + 1) for _ in range(n)] # ty: ignore[unresolved-attribute]
132 return p_mat, q, a_mat, b, cones
135def min_circle_clarabel(points: np.ndarray, verbose: bool = False) -> tuple[float, np.ndarray]:
136 """Compute the smallest enclosing circle for a set of points using Clarabel directly.
138 This function solves the same convex optimisation problem as
139 :func:`min_circle_cvx` but bypasses CVXPY and calls the Clarabel solver
140 directly. The second-order-cone program is assembled by
141 :func:`_build_soc_program`; this function then solves it and extracts the
142 optimal radius and centre.
144 Args:
145 points: A numpy array of shape ``(n, d)`` where *n* is the number of
146 points and *d* is the ambient dimension.
147 verbose: If ``True``, print Clarabel's iteration log. Defaults to
148 ``False``.
150 Returns:
151 A tuple ``(radius, center)`` where *radius* is the optimal enclosing
152 radius (float) and *center* is a numpy array of shape ``(d,)``.
154 Raises:
155 ValueError: If Clarabel does not return a ``Solved`` status.
157 Example:
158 >>> import numpy as np
159 >>> from cvxball.solver import min_circle_clarabel
160 >>> points = np.array([[0, 0], [1, 0], [0, 1]])
161 >>> radius, center = min_circle_clarabel(points)
162 """
163 p_mat, q, a_mat, b, cones = _build_soc_program(points)
165 # --- Solve ---------------------------------------------------------------
166 settings = clarabel.DefaultSettings.default() # ty: ignore[unresolved-attribute]
167 settings.verbose = verbose
169 solver = clarabel.DefaultSolver(p_mat, q, a_mat, b, cones, settings) # ty: ignore[unresolved-attribute]
170 solution = solver.solve()
172 if solution.status != clarabel.SolverStatus.Solved: # ty: ignore[unresolved-attribute]
173 raise ValueError(f"Clarabel did not converge: status = {solution.status}") # noqa: TRY003
175 return float(solution.x[0]), np.asarray(solution.x[1:])