11"""Core calibration routines for fairlex.
22
3- This module contains implementations of leximin‐ style calibration for survey
3+ This module contains implementations of leximin- style calibration for survey
44weights. Two variants are provided:
55
66* ``leximin_residual`` minimises the worst absolute deviation between the
7- calibrated and target margins (a ``min– max`` problem). It is akin to
7+ calibrated and target margins (a ``min- max`` problem). It is akin to
88 solving a Chebyshev approximation on the residuals. While this drives
99 margin errors down, it can lead to large deviations from the original
1010 weights if the margin targets are difficult to meet within bounds.
3838
3939import numpy as np
4040
41+ # Constants
42+ EXPECTED_MATRIX_DIMENSIONS = 2
43+
4144try :
4245 # SciPy is used for linear programming; HiGHS is fast and reliable.
4346 from scipy .optimize import linprog # type: ignore
4447except Exception : # pragma: no cover
45- linprog = None # type: ignore
48+ linprog = None
4649
4750
4851@dataclass
@@ -72,7 +75,9 @@ class CalibrationResult:
7275 message : str
7376
7477
75- def _validate_inputs (A : np .ndarray , b : np .ndarray , w0 : np .ndarray ) -> tuple [np .ndarray , np .ndarray , np .ndarray ]:
78+ def _validate_inputs (
79+ A : np .ndarray , b : np .ndarray , w0 : np .ndarray
80+ ) -> tuple [np .ndarray , np .ndarray , np .ndarray ]:
7681 """Validate and coerce input arrays to ensure they have compatible shapes.
7782
7883 Parameters
@@ -87,7 +92,7 @@ def _validate_inputs(A: np.ndarray, b: np.ndarray, w0: np.ndarray) -> tuple[np.n
8792 Returns
8893 -------
8994 (A, b, w0) : tuple of ndarrays
90- Validated and dtype‐ coerced versions of the inputs.
95+ Validated and dtype- coerced versions of the inputs.
9196
9297 Raises
9398 ------
@@ -97,8 +102,8 @@ def _validate_inputs(A: np.ndarray, b: np.ndarray, w0: np.ndarray) -> tuple[np.n
97102 A = np .asarray (A , dtype = float )
98103 b = np .asarray (b , dtype = float )
99104 w0 = np .asarray (w0 , dtype = float )
100- if A .ndim != 2 :
101- raise ValueError (f"A must be two‐ dimensional, got shape { A .shape } " )
105+ if A .ndim != EXPECTED_MATRIX_DIMENSIONS :
106+ raise ValueError (f"A must be two- dimensional, got shape { A .shape } " )
102107 m , n = A .shape
103108 if b .shape != (m ,):
104109 raise ValueError (f"b must be of shape { (m ,)} , got { b .shape } " )
@@ -107,7 +112,12 @@ def _validate_inputs(A: np.ndarray, b: np.ndarray, w0: np.ndarray) -> tuple[np.n
107112 return A , b , w0
108113
109114
110- def _solve_lp (c , A_ub , b_ub , bounds ):
115+ def _solve_lp (
116+ c : np .ndarray ,
117+ A_ub : np .ndarray ,
118+ b_ub : np .ndarray ,
119+ bounds : list [tuple [float | None , float | None ]],
120+ ) -> "scipy.optimize.OptimizeResult" : # type: ignore[name-defined] # noqa: F821
111121 """Solve a linear programming problem using SciPy HiGHS.
112122
113123 This helper centralises the call to ``scipy.optimize.linprog`` and
@@ -220,7 +230,74 @@ def leximin_residual(
220230 x = res .x
221231 w = x [:n ]
222232 epsilon = x [- 1 ]
223- return CalibrationResult (w = w , epsilon = epsilon , t = None , status = res .status , message = res .message )
233+ return CalibrationResult (
234+ w = w , epsilon = epsilon , t = None , status = res .status , message = res .message
235+ )
236+
237+
238+ def _setup_weight_fair_constraints (
239+ A : np .ndarray ,
240+ b : np .ndarray ,
241+ w0 : np .ndarray ,
242+ epsilon_opt : float ,
243+ * ,
244+ min_ratio : float ,
245+ max_ratio : float ,
246+ slack : float ,
247+ ) -> tuple [np .ndarray , np .ndarray , list [tuple [float | None , float | None ]]]:
248+ """Set up constraints for the weight-fair stage of calibration.
249+
250+ Returns
251+ -------
252+ A_ub : ndarray
253+ Inequality constraint matrix.
254+ b_ub : ndarray
255+ Inequality constraint right hand side.
256+ bounds : list
257+ Variable bounds.
258+ """
259+ m , n = A .shape
260+
261+ # Variables: w (n) and t (1)
262+ # Bounds: w within [w0*min_ratio, w0*max_ratio], t >= 0
263+ bounds = [(w0 [i ] * min_ratio , w0 [i ] * max_ratio ) for i in range (n )] + [(0 , None )]
264+
265+ # Build inequality constraints
266+ # Residual constraints: +/- (A_j w - b_j) <= epsilon_opt + slack
267+ # We'll build 2*m inequalities of the form A_j w + 0*t <= b_j + epsilon_opt + slack
268+ # and -A_j w + 0*t <= -b_j + epsilon_opt + slack
269+ total_constraints = 2 * m + 2 * n # residual constraints + weight change bounds
270+ A_ub = np .zeros ((total_constraints , n + 1 ))
271+ b_ub = np .zeros (total_constraints )
272+
273+ # Residual constraints
274+ for j in range (m ):
275+ # A_j w <= b_j + epsilon_opt + slack
276+ A_ub [2 * j , :n ] = A [j ]
277+ A_ub [2 * j , - 1 ] = 0.0
278+ b_ub [2 * j ] = b [j ] + epsilon_opt + slack
279+ # -A_j w <= -b_j + epsilon_opt + slack
280+ A_ub [2 * j + 1 , :n ] = - A [j ]
281+ A_ub [2 * j + 1 , - 1 ] = 0.0
282+ b_ub [2 * j + 1 ] = - b [j ] + epsilon_opt + slack
283+
284+ # Weight change bounds: for each i, w_i - w0_i <= t * w0_i and -(w_i - w0_i) <= t * w0_i
285+ offset = 2 * m
286+ for i in range (n ):
287+ # w_i - w0_i - t * w0_i <= 0 -> 1*w_i - w0_i* t <= w0_i
288+ row = np .zeros (n + 1 )
289+ row [i ] = 1.0
290+ row [- 1 ] = - w0 [i ]
291+ A_ub [offset + 2 * i ] = row
292+ b_ub [offset + 2 * i ] = w0 [i ]
293+ # -w_i + w0_i - t * w0_i <= 0 -> -1*w_i - w0_i* t <= -w0_i
294+ row = np .zeros (n + 1 )
295+ row [i ] = - 1.0
296+ row [- 1 ] = - w0 [i ]
297+ A_ub [offset + 2 * i + 1 ] = row
298+ b_ub [offset + 2 * i + 1 ] = - w0 [i ]
299+
300+ return A_ub , b_ub , bounds
224301
225302
226303def leximin_weight_fair (
@@ -280,64 +357,40 @@ def leximin_weight_fair(
280357 if return_stages :
281358 return stage1 , stage1
282359 return stage1
360+
283361 # Set up the second stage: minimise t subject to residual constraints and weight change bounds
284362 A , b , w0 = _validate_inputs (A , b , w0 )
285- m , n = A .shape
286- epsilon_opt = stage1 . epsilon
363+ n = A .shape [ 1 ]
364+
287365 # Variables: w (n) and t (1)
288366 # Objective: minimise t
289367 c = np .zeros (n + 1 )
290368 c [- 1 ] = 1.0
291- # Bounds: w within [w0*min_ratio, w0*max_ratio], t >= 0
292- bounds = [(w0 [i ] * min_ratio , w0 [i ] * max_ratio ) for i in range (n )] + [(0 , None )]
293- # Build inequality constraints
294- # Residual constraints: +/- (A_j w - b_j) <= epsilon_opt + slack
295- # We'll build 2*m inequalities of the form A_j w + 0*t <= b_j + epsilon_opt + slack
296- # and -A_j w + 0*t <= -b_j + epsilon_opt + slack
297- total_constraints = 2 * m + 2 * n # residual constraints + weight change bounds
298- A_ub = np .zeros ((total_constraints , n + 1 ))
299- b_ub = np .zeros (total_constraints )
300- # Residual constraints
301- for j in range (m ):
302- # A_j w <= b_j + epsilon_opt + slack
303- A_ub [2 * j , :n ] = A [j ]
304- A_ub [2 * j , - 1 ] = 0.0
305- b_ub [2 * j ] = b [j ] + epsilon_opt + slack
306- # -A_j w <= -b_j + epsilon_opt + slack
307- A_ub [2 * j + 1 , :n ] = - A [j ]
308- A_ub [2 * j + 1 , - 1 ] = 0.0
309- b_ub [2 * j + 1 ] = - b [j ] + epsilon_opt + slack
310- # Weight change bounds: for each i, w_i - w0_i <= t * w0_i and -(w_i - w0_i) <= t * w0_i
311- offset = 2 * m
312- for i in range (n ):
313- # w_i - w0_i - t * w0_i <= 0 -> 1*w_i - w0_i* t <= w0_i
314- row = np .zeros (n + 1 )
315- row [i ] = 1.0
316- row [- 1 ] = - w0 [i ]
317- A_ub [offset + 2 * i ] = row
318- b_ub [offset + 2 * i ] = w0 [i ]
319- # -w_i + w0_i - t * w0_i <= 0 -> -1*w_i - w0_i* t <= -w0_i
320- row = np .zeros (n + 1 )
321- row [i ] = - 1.0
322- row [- 1 ] = - w0 [i ]
323- A_ub [offset + 2 * i + 1 ] = row
324- b_ub [offset + 2 * i + 1 ] = - w0 [i ]
369+
370+ # Set up constraints using helper function
371+ A_ub , b_ub , bounds = _setup_weight_fair_constraints (
372+ A , b , w0 , stage1 .epsilon , min_ratio = min_ratio , max_ratio = max_ratio , slack = slack
373+ )
374+
325375 res = _solve_lp (c , A_ub , b_ub , bounds )
326376 if not res .success :
327377 stage2 = CalibrationResult (
328378 w = np .full_like (w0 , np .nan ),
329- epsilon = epsilon_opt ,
379+ epsilon = stage1 . epsilon ,
330380 t = np .nan ,
331381 status = res .status ,
332382 message = res .message ,
333383 )
334384 if return_stages :
335385 return stage1 , stage2
336386 return stage2
387+
337388 x = res .x
338389 w = x [:n ]
339390 t_opt = x [- 1 ]
340- stage2 = CalibrationResult (w = w , epsilon = epsilon_opt , t = t_opt , status = res .status , message = res .message )
391+ stage2 = CalibrationResult (
392+ w = w , epsilon = stage1 .epsilon , t = t_opt , status = res .status , message = res .message
393+ )
341394 if return_stages :
342395 return stage1 , stage2
343396 return stage2
0 commit comments