Coverage for src / cvxball / solver.py: 100%
40 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-12 10:49 +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: dict[str, 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 r * np.ones(points.shape[0]),
51 points - x, # Broadcasting handles this automatically
52 axis=1,
53 )
54 ]
56 problem = cp.Problem(objective=objective, constraints=constraints)
57 problem.solve(**kwargs)
59 # Ensure the problem was solved successfully
60 if r.value is None or x.value is None:
61 raise ValueError("Optimization failed to find a solution") # noqa: TRY003
63 return float(r.value[0]), x.value
66def min_circle_clarabel(points: np.ndarray, verbose: bool = False) -> tuple[float, np.ndarray]:
67 """Compute the smallest enclosing circle for a set of points using Clarabel directly.
69 This function solves the same convex optimisation problem as
70 :func:`min_circle_cvx` but bypasses CVXPY and calls the Clarabel solver
71 directly. The problem is assembled in Clarabel's standard form::
73 minimise (1/2) z' P z + q' z
74 subject to A z + s = b, s ∈ K
76 where the decision vector is ``z = [r, x₁, …, x_d]`` (radius followed by
77 the d centre coordinates), the objective is to minimise *r* (so ``P = 0``,
78 ``q = e₀``), and the feasible set is a product of *n* second-order cones.
80 For each point ``p_i`` we require ``[r, p_i - x] in Q^{d+1}``, which gives
81 one SOC block of dimension ``d + 1`` per point.
83 Args:
84 points: A numpy array of shape ``(n, d)`` where *n* is the number of
85 points and *d* is the ambient dimension.
86 verbose: If ``True``, print Clarabel's iteration log. Defaults to
87 ``False``.
89 Returns:
90 A tuple ``(radius, center)`` where *radius* is the optimal enclosing
91 radius (float) and *center* is a numpy array of shape ``(d,)``.
93 Raises:
94 ValueError: If Clarabel does not return a ``Solved`` status.
96 Example:
97 >>> import numpy as np
98 >>> from cvxball.solver import min_circle_clarabel
99 >>> points = np.array([[0, 0], [1, 0], [0, 1]])
100 >>> radius, center = min_circle_clarabel(points)
101 """
102 n, d = points.shape
103 n_vars = 1 + d # decision vector: [r, x_1, ..., x_d]
105 # --- Objective: minimise r -----------------------------------------------
106 p_mat = sp.csc_matrix((n_vars, n_vars))
107 q = np.zeros(n_vars)
108 q[0] = 1.0
110 # --- Constraints: one SOC block of size (d+1) per point ------------------
111 # We need b - a_mat @ z = s where s in K.
112 # For point i the desired slack is s = [r, p_i - x], so:
113 # row i*(d+1) : b = 0, a_mat col 0 = -1 (gives s_0 = r)
114 # row i*(d+1)+j : b = p_i[j], a_mat col j = +1 (gives s_j = p_ij - x_j)
115 total_rows = n * (d + 1)
117 # Entries for the r column (column 0): -1 at each block's first row
118 r_rows = np.arange(n) * (d + 1)
120 # Entries for the x columns (columns 1..d): +1 at each block's inner rows
121 x_row_offsets = np.arange(n)[:, None] * (d + 1) + np.arange(1, d + 1)[None, :] # (n, d)
122 x_rows = x_row_offsets.ravel()
123 x_cols = np.tile(np.arange(1, d + 1), n)
125 all_rows = np.concatenate([r_rows, x_rows])
126 all_cols = np.concatenate([np.zeros(n, dtype=np.intp), x_cols])
127 all_vals = np.concatenate([-np.ones(n), np.ones(n * d)])
129 a_mat = sp.csc_matrix((all_vals, (all_rows, all_cols)), shape=(total_rows, n_vars))
131 b = np.zeros(total_rows)
132 b[x_rows] = points.ravel()
134 # --- Cones: n SOC cones each of dimension (d+1) --------------------------
135 cones = [clarabel.SecondOrderConeT(d + 1) for _ in range(n)] # ty: ignore[unresolved-attribute]
137 # --- Solve ---------------------------------------------------------------
138 settings = clarabel.DefaultSettings.default() # ty: ignore[unresolved-attribute]
139 settings.verbose = verbose
141 solver = clarabel.DefaultSolver(p_mat, q, a_mat, b, cones, settings) # ty: ignore[unresolved-attribute]
142 solution = solver.solve()
144 if solution.status != clarabel.SolverStatus.Solved: # ty: ignore[unresolved-attribute]
145 raise ValueError(f"Clarabel did not converge: status = {solution.status}") # noqa: TRY003
147 return float(solution.x[0]), np.asarray(solution.x[1:])