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

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 # CVXPY variable for the radius 

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

32 # CVXPY variable for the midpoint 

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

34 

35 # Minimize the radius 

36 objective = cp.Minimize(r) 

37 

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

39 constraints = [ 

40 cp.SOC( 

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

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

43 axis=1, 

44 ) 

45 ] 

46 

47 # Create and solve the problem 

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

49 problem.solve(**kwargs) 

50 

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

52 

53 

54def min_circle_cvx_3() -> cp.Problem: 

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

56 

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

58 

59 Returns: 

60 CVXPY Problem object ready to be solved 

61 """ 

62 # CVXPY variable for the radius 

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

64 # CVXPY variable for the midpoint 

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

66 # Parameter for the 3 points 

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

68 

69 # Minimize the radius 

70 objective = cp.Minimize(r) 

71 

72 # Second-order cone constraints for 3 points 

73 constraints = [ 

74 cp.SOC( 

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

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

77 axis=1, 

78 ) 

79 ] 

80 

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

82 return problem 

83 

84 

85def min_circle_cvx_2() -> cp.Problem: 

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

87 

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

89 

90 Returns: 

91 CVXPY Problem object ready to be solved 

92 """ 

93 # CVXPY variable for the radius 

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

95 # CVXPY variable for the midpoint 

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

97 # Parameter for the 2 points 

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

99 

100 # Minimize the radius 

101 objective = cp.Minimize(r) 

102 

103 # Second-order cone constraints for 2 points 

104 constraints = [ 

105 cp.SOC( 

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

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

108 axis=1, 

109 ) 

110 ] 

111 

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

113 return problem 

114 

115 

116if __name__ == "__main__": 

117 # Example with 3 points 

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

119 problem = min_circle_cvx_3() 

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

121 problem.solve(solver="CLARABEL") 

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

123 

124 # Performance testing with 2 points 

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

126 problem = min_circle_cvx_2() 

127 

128 def f() -> None: 

129 """Benchmark function for parameterized problem.""" 

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

131 problem.solve(solver="CLARABEL") 

132 

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

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

135 

136 def g() -> Circle: 

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

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

139 

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

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