Coverage for suppy\feasibility\_halfspaces\_ams_extrapolations.py: 93%

1from typing import List 

2import numpy as np 

3import numpy.typing as npt 

4from scipy import sparse 

5 

6try: 

7 import cupy as cp 

8 

9 NO_GPU = False 

10 

11except ImportError: 

12 NO_GPU = True 

13 cp = np 

14 

15from suppy.feasibility._halfspaces._ams_algorithms import SimultaneousAMSHalfspace 

16 

17 

18class ExtrapolatedLandweberHalfspace(SimultaneousAMSHalfspace): 

19 """ 

20 Extrapolated Landweber method for solving linear inequalities of the 

21 form Ax <= b. 

22 

23 Parameters 

24 ---------- 

25 A : npt.NDArray or sparse.sparray 

26 Matrix for linear systems 

27 b : npt.NDArray 

28 Bound for linear systems 

29 weights : None or List[float], optional 

30 The weights for the constraints, by default None. 

31 relaxation : float, optional 

32 The outer relaxation parameter, by default 1.0. 

33 proximity_flag : bool, optional 

34 A flag indicating whether to use proximity measures, by default True. 

35 

36 References 

37 ---------- 

38 - [1] https://doi.org/10.1007/s11075-025-02163-0 

39 - [2] https://doi.org/10.1007/978-3-642-30901-4 

40 """ 

41 

42 def __init__( 

43 self, 

44 A: npt.NDArray | sparse.sparray, 

45 b: npt.NDArray, 

46 relaxation: float = 1.0, 

47 weights: None | List[float] | npt.NDArray = None, 

48 proximity_flag: bool = True, 

49 ): 

50 super().__init__(A, b, 1.0, relaxation, weights, proximity_flag) 

51 self.a_i = self.A.row_norm(2, 2) 

52 self.weight_norm = self.weights / self.a_i 

53 self.sigmas = [] 

54 

55 def _project(self, x: npt.NDArray) -> np.ndarray: 

56 xp = cp if self._use_gpu else np 

57 p = self.map(x) 

58 res = self.b - p 

59 res_idx = res < 0 

60 if not xp.any(res_idx): 

61 self.sigmas.append(0) 

62 return x 

63 t = self.weight_norm[res_idx] * res[res_idx] 

64 t_2 = t @ self.A[res_idx, :] 

65 sig = (res[res_idx] @ t) / (t_2 @ t_2) 

66 self.sigmas.append(sig) 

67 x += sig * t_2 

68 

69 return x 

70 

71 

72class AdaptiveStepLandweberHalfspace(SimultaneousAMSHalfspace): 

73 """ 

74 Adaptive step-size Landweber algorithm for solving linear inequalities 

75 of 

76 the form Ax <= b. 

77 

78 Parameters 

79 ---------- 

80 A : npt.NDArray or sparse.sparray 

81 Matrix for linear systems 

82 b : npt.NDArray 

83 Bound for linear systems 

84 weights : None or List[float], optional 

85 The weights for the constraints, by default None. 

86 relaxation : float, optional 

87 The outer relaxation parameter, by default 1.0. 

88 proximity_flag : bool, optional 

89 A flag indicating whether to use proximity measures, by default True. 

90 

91 References 

92 ---------- 

93 - [1] https://doi.org/10.1515/jiip-2015-0082 

94 """ 

95 

96 def __init__( 

97 self, 

98 A: npt.NDArray | sparse.sparray, 

99 b: npt.NDArray, 

100 relaxation: float = 1.0, 

101 weights: None | List[float] | npt.NDArray = None, 

102 proximity_flag: bool = True, 

103 ): 

104 super().__init__(A, b, 1.0, relaxation, weights, proximity_flag) 

105 

106 def _project(self, x: npt.NDArray) -> np.ndarray: 

107 xp = cp if self._use_gpu else np 

108 p = self.map(x) 

109 res = self.b - p 

110 res_idx = res < 0 

111 

112 if not xp.any(res_idx): 

113 return x 

114 

115 u_t = (self.weights[res_idx] * res[res_idx]) @ self.A[res_idx, :] 

116 Au_t = self.A @ u_t # A*AT*res 

117 sig = (u_t @ u_t) / (Au_t @ (self.weights * Au_t)) 

118 x += sig * u_t 

119 

120 return x