Coverage for src/cvxcla/_events.py: 100%

1"""Critical-lambda event scans for the Critical Line Algorithm. 

2 

3Along a critical-line segment ``w(lam) = r_alpha + lam * r_beta`` two families of 

4events end the segment's validity: a *box* event (a free weight reaching a bound, 

5or a blocked weight's multiplier changing sign) and an *inequality-row* event (an 

6inactive ``G w <= h`` row's slack reaching zero, or an active row's multiplier 

7changing sign). Both reduce to the same ``-intercept / slope`` critical-lambda 

8ratio, computed here as pure functions and stacked by :func:`segment_events` into 

9the ``(n + p, 4)`` matrix the generic path tracer scans. 

10""" 

11 

12from __future__ import annotations 

13 

14import numpy as np 

15from numpy.typing import NDArray 

16 

17from ._kkt import Segment 

18 

19 

20def event_ratios( 

21 r_alpha: NDArray[np.float64], 

22 r_beta: NDArray[np.float64], 

23 gamma: NDArray[np.float64], 

24 delta: NDArray[np.float64], 

25 free_in: NDArray[np.bool_], 

26 at_upper: NDArray[np.bool_], 

27 at_lower: NDArray[np.bool_], 

28 lower: NDArray[np.float64], 

29 upper: NDArray[np.float64], 

30) -> NDArray[np.float64]: 

31 """Critical lambda for every candidate box event, as an ``(n, 4)`` matrix. 

32 

33 Along the segment ``w(lam) = r_alpha + lam * r_beta`` a free weight can 

34 reach a box bound (columns 0/1, "moves to a bound") and a blocked weight's 

35 multiplier can change sign so it re-enters the free set (columns 2/3, 

36 "leaves a bound"). Entries with no event are ``-inf``. 

37 

38 A free weight moves with even a tiny slope, so given a long enough lam 

39 range it still crosses a bound; filtering slopes at the classification 

40 tolerance would miss such crossings and let weights drift out of bounds. 

41 Only slopes at floating-point noise level are excluded: below 

42 ``sqrt(machine epsilon)`` a slope is indistinguishable from solve noise, and 

43 the huge ratios it would produce only amplify rounding errors. 

44 

45 Args: 

46 r_alpha: Segment intercept ``w(0)``. 

47 r_beta: Segment slope ``dw/dlam``. 

48 gamma: Multiplier gradient for the alpha system. 

49 delta: Multiplier gradient for the beta system. 

50 free_in: Mask of assets in the reduced solve. 

51 at_upper: Mask of assets blocked at their upper bound. 

52 at_lower: Mask of assets blocked at their lower bound. 

53 lower: Per-asset lower bounds. 

54 upper: Per-asset upper bounds. 

55 

56 Returns: 

57 The ``(n, 4)`` matrix of critical lambdas. 

58 """ 

59 ns = len(r_alpha) 

60 eps = np.sqrt(np.finfo(np.float64).eps) 

61 # 4 columns = the 4 event types; extra unused columns are harmless. 

62 l_mat = np.full((ns, 4), -np.inf) # pragma: no mutate 

63 

64 # Precompute each event mask exactly once. The <,> vs <=,>= choice at 

65 # the eps boundary is numerically irrelevant — a slope/derivative 

66 # landing exactly on +/-sqrt(machine-eps) never occurs with real 

67 # data — so those boundary comparisons are marked no-mutate. 

68 beta_down = free_in & (r_beta < -eps) # pragma: no mutate 

69 beta_up = free_in & (r_beta > +eps) # pragma: no mutate 

70 delta_down = at_upper & (delta < -eps) # pragma: no mutate 

71 delta_up = at_lower & (delta > +eps) # pragma: no mutate 

72 

73 # Columns 0,1 are "moves to a bound" (free->blocked) and 2,3 are 

74 # "leaves a bound" (blocked->free); the next-free update only tests 

75 # dirchg >= 2, so swapping a column *within* a group (0<->1 or 2<->3) is 

76 # behaviourally identical and marked no-mutate. Crossing the 1<->2 group 

77 # boundary IS exercised by the frontier tests. 

78 l_mat[beta_down, 0] = (upper[beta_down] - r_alpha[beta_down]) / r_beta[beta_down] # pragma: no mutate 

79 l_mat[beta_up, 1] = (lower[beta_up] - r_alpha[beta_up]) / r_beta[beta_up] 

80 l_mat[delta_down, 2] = -gamma[delta_down] / delta[delta_down] # pragma: no mutate 

81 l_mat[delta_up, 3] = -gamma[delta_up] / delta[delta_up] 

82 return l_mat 

83 

84 

85def ineq_event_ratios( 

86 r_alpha: NDArray[np.float64], 

87 r_beta: NDArray[np.float64], 

88 eta_alpha: NDArray[np.float64], 

89 eta_beta: NDArray[np.float64], 

90 active_ineq: NDArray[np.bool_], 

91 g: NDArray[np.float64], 

92 h: NDArray[np.float64], 

93) -> NDArray[np.float64]: 

94 """Critical lambda for every inequality-row event, as a ``(p, 4)`` matrix. 

95 

96 The row analogue of :func:`event_ratios`. Along the segment an *inactive* 

97 row ``i`` becomes active when its slack ``s_i(lam) = g_i w(lam) - h_i`` 

98 rises to zero from the feasible (negative) side (column 0); an *active* 

99 row releases when its multiplier ``eta_i(lam)`` falls to zero (column 1). 

100 Both are affine in ``lam``, so the critical lambda is the same 

101 ``-intercept / slope`` ratio used for the box events, with the same 

102 ``sqrt(machine eps)`` slope floor: a slope below noise level is 

103 indistinguishable from solve round-off and would only produce a huge, 

104 rounding-dominated ratio. Entries with no event are ``-inf``. Columns 2 

105 and 3 are unused (kept so the block stacks onto the ``(n, 4)`` box block). 

106 

107 Args: 

108 r_alpha: Segment intercept ``w(0)``. 

109 r_beta: Segment slope ``dw/dlam``. 

110 eta_alpha: Affine inequality-multiplier intercept (length ``p``). 

111 eta_beta: Affine inequality-multiplier slope (length ``p``). 

112 active_ineq: Boolean mask (length ``p``) of the active inequality rows. 

113 g: Inequality-constraint matrix ``G`` of ``G w <= h`` (``(p, n)``). 

114 h: Inequality-constraint right-hand side ``h`` (length ``p``). 

115 

116 Returns: 

117 The ``(p, 4)`` matrix of critical lambdas. 

118 """ 

119 p = g.shape[0] 

120 l_mat = np.full((p, 4), -np.inf) # pragma: no mutate 

121 if p == 0: 

122 return l_mat 

123 

124 eps = np.sqrt(np.finfo(np.float64).eps) 

125 inactive = ~active_ineq 

126 

127 # Enter: an inactive row's slack rises to zero. The slope/intercept split 

128 # comes straight from the affine weights; the slope sign mirrors the box 

129 # "moves to a bound" event (decreasing lam must raise the slack). 

130 s_alpha = g @ r_alpha - h 

131 s_beta = g @ r_beta 

132 enter = inactive & (s_beta < -eps) # pragma: no mutate 

133 l_mat[enter, 0] = -s_alpha[enter] / s_beta[enter] 

134 

135 # Release: an active row's non-negative multiplier falls back to zero, 

136 # the row analogue of a blocked multiplier changing sign. 

137 release = active_ineq & (eta_beta > +eps) # pragma: no mutate 

138 l_mat[release, 1] = -eta_alpha[release] / eta_beta[release] 

139 return l_mat 

140 

141 

142def segment_events( 

143 segment: Segment, 

144 lower: NDArray[np.float64], 

145 upper: NDArray[np.float64], 

146 g: NDArray[np.float64], 

147 h: NDArray[np.float64], 

148) -> NDArray[np.float64]: 

149 """Return the ``(n + p, 4)`` matrix of candidate critical lambdas for ``segment``. 

150 

151 The first ``n`` rows are the box events of :func:`event_ratios` (a free weight 

152 reaching a bound, a blocked multiplier changing sign); the trailing ``p`` rows 

153 are the inequality-row events of :func:`ineq_event_ratios` (an inactive row's 

154 slack reaching zero, an active row's multiplier changing sign). The generic 

155 tracer treats the two blocks uniformly; ``CLA.step`` decodes a row index 

156 ``>= n`` as a row event. 

157 

158 Args: 

159 segment: The critical-line segment to scan. 

160 lower: Per-asset lower bounds. 

161 upper: Per-asset upper bounds. 

162 g: Inequality-constraint matrix ``G`` of ``G w <= h`` (``(p, n)``). 

163 h: Inequality-constraint right-hand side ``h`` (length ``p``). 

164 

165 Returns: 

166 The stacked ``(n + p, 4)`` event matrix. 

167 """ 

168 box = event_ratios( 

169 segment.r_alpha, 

170 segment.r_beta, 

171 segment.gamma, 

172 segment.delta, 

173 segment.free_in, 

174 segment.at_upper, 

175 segment.at_lower, 

176 lower, 

177 upper, 

178 ) 

179 ineq = ineq_event_ratios( 

180 segment.r_alpha, segment.r_beta, segment.eta_alpha, segment.eta_beta, segment.active_ineq, g, h 

181 ) 

182 return np.vstack([box, ineq])