Coverage for src / min_circle / cvx.py: 100%

1"""CVXPY implementation for finding the minimum enclosing circle. 

2 

3This module provides functions to find the minimum enclosing circle 

4using convex optimization with the CVXPY library. 

5""" 

6 

7import statistics 

8import timeit 

9from typing import Any 

10 

11import cvxpy as cp 

12import numpy as np 

13 

14from .utils.circle import Circle 

15 

16 

17def min_circle_cvx(points: np.ndarray, **kwargs: Any) -> Circle: 

18 """Find the minimum enclosing circle using convex optimization. 

19 

20 Uses CVXPY to formulate and solve the minimum enclosing circle problem 

21 as a second-order cone program (SOCP). 

22 

23 Args: 

24 points: Array of 2D points with shape (n, 2) 

25 **kwargs: Additional keyword arguments to pass to the solver 

26 

27 Returns: 

28 Circle object containing the center and radius of the minimum enclosing circle 

29 

30 Examples: 

31 >>> import numpy as np 

32 >>> points = np.array([[0.0, 0.0], [2.0, 0.0]]) 

33 >>> circle = min_circle_cvx(points, solver="CLARABEL") 

34 >>> round(circle.radius, 6) 

35 1.0 

36 >>> np.round(circle.center, 6) 

37 array([1., 0.]) 

38 """ 

39 # CVXPY variable for the radius 

40 r = cp.Variable(name="Radius") 

41 # CVXPY variable for the midpoint 

42 x = cp.Variable(points.shape[1], name="Midpoint") 

43 

44 # Minimize the radius 

45 objective = cp.Minimize(r) 

46 

47 # Second-order cone constraints ensuring all points are within the circle 

48 constraints = [ 

49 cp.SOC( 

50 r * np.ones(points.shape[0]), 

51 points - cp.outer(np.ones(points.shape[0]), x), 

52 axis=1, 

53 ) 

54 ] 

55 

56 # Create and solve the problem 

57 problem = cp.Problem(objective=objective, constraints=constraints) 

58 problem.solve(**kwargs) 

59 

60 return Circle(radius=float(r.value), center=x.value) 

61 

62 

63def min_circle_cvx_3() -> cp.Problem: 

64 """Create a CVXPY problem for finding the minimum circle for 3 points. 

65 

66 Creates a parameterized problem that can be reused with different point sets. 

67 

68 Returns: 

69 CVXPY Problem object ready to be solved 

70 """ 

71 # CVXPY variable for the radius 

72 r = cp.Variable(name="Radius") 

73 # CVXPY variable for the midpoint 

74 x = cp.Variable((1, 2), name="Midpoint") 

75 # Parameter for the 3 points 

76 p = cp.Parameter((3, 2), "points") 

77 

78 # Minimize the radius 

79 objective = cp.Minimize(r) 

80 

81 # Second-order cone constraints for 3 points 

82 constraints = [ 

83 cp.SOC( 

84 cp.hstack([r, r, r]), 

85 p - cp.vstack([x, x, x]), 

86 axis=1, 

87 ) 

88 ] 

89 

90 problem = cp.Problem(objective=objective, constraints=constraints) 

91 return problem 

92 

93 

94def min_circle_cvx_2() -> cp.Problem: 

95 """Create a CVXPY problem for finding the minimum circle for 2 points. 

96 

97 Creates a parameterized problem that can be reused with different point sets. 

98 

99 Returns: 

100 CVXPY Problem object ready to be solved 

101 """ 

102 # CVXPY variable for the radius 

103 r = cp.Variable(name="Radius") 

104 # CVXPY variable for the midpoint 

105 x = cp.Variable((1, 2), name="Midpoint") 

106 # Parameter for the 2 points 

107 p = cp.Parameter((2, 2), "points") 

108 

109 # Minimize the radius 

110 objective = cp.Minimize(r) 

111 

112 # Second-order cone constraints for 2 points 

113 constraints = [ 

114 cp.SOC( 

115 cp.hstack([r, r]), 

116 p - cp.vstack([x, x]), 

117 axis=1, 

118 ) 

119 ] 

120 

121 problem = cp.Problem(objective=objective, constraints=constraints) 

122 return problem 

123 

124 

125if __name__ == "__main__": # pragma: no cover 

126 # Example with 3 points 

127 p = np.array([[0, 0], [0.0, 2.0], [4.0, 4.0]]) 

128 problem = min_circle_cvx_3() 

129 problem.param_dict["points"].value = p 

130 problem.solve(solver="CLARABEL") 

131 print(problem.var_dict["Radius"].value) 

132 

133 # Performance testing with 2 points 

134 p = np.random.randn(2, 2) 

135 problem = min_circle_cvx_2() 

136 

137 def f() -> None: 

138 """Benchmark function for parameterized problem.""" 

139 problem.param_dict["points"].value = p 

140 problem.solve(solver="CLARABEL") 

141 

142 results = timeit.repeat(f, number=1, repeat=10) 

143 print(f"Parameterized problem mean time: {statistics.mean(results):.6f} seconds") 

144 

145 def g() -> Circle: 

146 """Benchmark function for direct problem formulation.""" 

147 return min_circle_cvx(p, solver="CLARABEL") 

148 

149 results = timeit.repeat(g, number=1, repeat=10) 

150 print(f"Direct problem mean time: {statistics.mean(results):.6f} seconds")