Coverage for src/cvxcla/_events.py: 100%
30 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-11 09:58 +0000
1"""Critical-lambda event scans for the Critical Line Algorithm.
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 ``CLA.event_matrix`` into
9the ``(n + p, 4)`` matrix the generic path tracer scans.
10"""
12from __future__ import annotations
14import numpy as np
15from numpy.typing import NDArray
18def event_ratios(
19 r_alpha: NDArray[np.float64],
20 r_beta: NDArray[np.float64],
21 gamma: NDArray[np.float64],
22 delta: NDArray[np.float64],
23 free_in: NDArray[np.bool_],
24 at_upper: NDArray[np.bool_],
25 at_lower: NDArray[np.bool_],
26 lower: NDArray[np.float64],
27 upper: NDArray[np.float64],
28) -> NDArray[np.float64]:
29 """Critical lambda for every candidate box event, as an ``(n, 4)`` matrix.
31 Along the segment ``w(lam) = r_alpha + lam * r_beta`` a free weight can
32 reach a box bound (columns 0/1, "moves to a bound") and a blocked weight's
33 multiplier can change sign so it re-enters the free set (columns 2/3,
34 "leaves a bound"). Entries with no event are ``-inf``.
36 A free weight moves with even a tiny slope, so given a long enough lam
37 range it still crosses a bound; filtering slopes at the classification
38 tolerance would miss such crossings and let weights drift out of bounds.
39 Only slopes at floating-point noise level are excluded: below
40 ``sqrt(machine epsilon)`` a slope is indistinguishable from solve noise, and
41 the huge ratios it would produce only amplify rounding errors.
43 Args:
44 r_alpha: Segment intercept ``w(0)``.
45 r_beta: Segment slope ``dw/dlam``.
46 gamma: Multiplier gradient for the alpha system.
47 delta: Multiplier gradient for the beta system.
48 free_in: Mask of assets in the reduced solve.
49 at_upper: Mask of assets blocked at their upper bound.
50 at_lower: Mask of assets blocked at their lower bound.
51 lower: Per-asset lower bounds.
52 upper: Per-asset upper bounds.
54 Returns:
55 The ``(n, 4)`` matrix of critical lambdas.
56 """
57 ns = len(r_alpha)
58 eps = np.sqrt(np.finfo(np.float64).eps)
59 # 4 columns = the 4 event types; extra unused columns are harmless.
60 l_mat = np.full((ns, 4), -np.inf) # pragma: no mutate
62 # Precompute each event mask exactly once. The <,> vs <=,>= choice at
63 # the eps boundary is numerically irrelevant — a slope/derivative
64 # landing exactly on +/-sqrt(machine-eps) never occurs with real
65 # data — so those boundary comparisons are marked no-mutate.
66 beta_down = free_in & (r_beta < -eps) # pragma: no mutate
67 beta_up = free_in & (r_beta > +eps) # pragma: no mutate
68 delta_down = at_upper & (delta < -eps) # pragma: no mutate
69 delta_up = at_lower & (delta > +eps) # pragma: no mutate
71 # Columns 0,1 are "moves to a bound" (free->blocked) and 2,3 are
72 # "leaves a bound" (blocked->free); the next-free update only tests
73 # dirchg >= 2, so swapping a column *within* a group (0<->1 or 2<->3) is
74 # behaviourally identical and marked no-mutate. Crossing the 1<->2 group
75 # boundary IS exercised by the frontier tests.
76 l_mat[beta_down, 0] = (upper[beta_down] - r_alpha[beta_down]) / r_beta[beta_down] # pragma: no mutate
77 l_mat[beta_up, 1] = (lower[beta_up] - r_alpha[beta_up]) / r_beta[beta_up]
78 l_mat[delta_down, 2] = -gamma[delta_down] / delta[delta_down] # pragma: no mutate
79 l_mat[delta_up, 3] = -gamma[delta_up] / delta[delta_up]
80 return l_mat
83def ineq_event_ratios(
84 r_alpha: NDArray[np.float64],
85 r_beta: NDArray[np.float64],
86 eta_alpha: NDArray[np.float64],
87 eta_beta: NDArray[np.float64],
88 active_ineq: NDArray[np.bool_],
89 g: NDArray[np.float64],
90 h: NDArray[np.float64],
91) -> NDArray[np.float64]:
92 """Critical lambda for every inequality-row event, as a ``(p, 4)`` matrix.
94 The row analogue of :func:`event_ratios`. Along the segment an *inactive*
95 row ``i`` becomes active when its slack ``s_i(lam) = g_i w(lam) - h_i``
96 rises to zero from the feasible (negative) side (column 0); an *active*
97 row releases when its multiplier ``eta_i(lam)`` falls to zero (column 1).
98 Both are affine in ``lam``, so the critical lambda is the same
99 ``-intercept / slope`` ratio used for the box events, with the same
100 ``sqrt(machine eps)`` slope floor: a slope below noise level is
101 indistinguishable from solve round-off and would only produce a huge,
102 rounding-dominated ratio. Entries with no event are ``-inf``. Columns 2
103 and 3 are unused (kept so the block stacks onto the ``(n, 4)`` box block).
105 Args:
106 r_alpha: Segment intercept ``w(0)``.
107 r_beta: Segment slope ``dw/dlam``.
108 eta_alpha: Affine inequality-multiplier intercept (length ``p``).
109 eta_beta: Affine inequality-multiplier slope (length ``p``).
110 active_ineq: Boolean mask (length ``p``) of the active inequality rows.
111 g: Inequality-constraint matrix ``G`` of ``G w <= h`` (``(p, n)``).
112 h: Inequality-constraint right-hand side ``h`` (length ``p``).
114 Returns:
115 The ``(p, 4)`` matrix of critical lambdas.
116 """
117 p = g.shape[0]
118 l_mat = np.full((p, 4), -np.inf) # pragma: no mutate
119 if p == 0:
120 return l_mat
122 eps = np.sqrt(np.finfo(np.float64).eps)
123 inactive = ~active_ineq
125 # Enter: an inactive row's slack rises to zero. The slope/intercept split
126 # comes straight from the affine weights; the slope sign mirrors the box
127 # "moves to a bound" event (decreasing lam must raise the slack).
128 s_alpha = g @ r_alpha - h
129 s_beta = g @ r_beta
130 enter = inactive & (s_beta < -eps) # pragma: no mutate
131 l_mat[enter, 0] = -s_alpha[enter] / s_beta[enter]
133 # Release: an active row's non-negative multiplier falls back to zero,
134 # the row analogue of a blocked multiplier changing sign.
135 release = active_ineq & (eta_beta > +eps) # pragma: no mutate
136 l_mat[release, 1] = -eta_alpha[release] / eta_beta[release]
137 return l_mat