Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_lsq/bvls.py: 5%

1"""Bounded-variable least-squares algorithm.""" 

2import numpy as np 

3from numpy.linalg import norm, lstsq 

4from scipy.optimize import OptimizeResult 

5 

6from .common import print_header_linear, print_iteration_linear 

7 

8 

9def compute_kkt_optimality(g, on_bound): 

10 """Compute the maximum violation of KKT conditions.""" 

11 g_kkt = g * on_bound 

12 free_set = on_bound == 0 

13 g_kkt[free_set] = np.abs(g[free_set]) 

14 return np.max(g_kkt) 

15 

16 

17def bvls(A, b, x_lsq, lb, ub, tol, max_iter, verbose, rcond=None): 

18 m, n = A.shape 

19 

20 x = x_lsq.copy() 

21 on_bound = np.zeros(n) 

22 

23 mask = x <= lb 

24 x[mask] = lb[mask] 

25 on_bound[mask] = -1 

26 

27 mask = x >= ub 

28 x[mask] = ub[mask] 

29 on_bound[mask] = 1 

30 

31 free_set = on_bound == 0 

32 active_set = ~free_set 

33 free_set, = np.nonzero(free_set) 

34 

35 r = A.dot(x) - b 

36 cost = 0.5 * np.dot(r, r) 

37 initial_cost = cost 

38 g = A.T.dot(r) 

39 

40 cost_change = None 

41 step_norm = None 

42 iteration = 0 

43 

44 if verbose == 2: 

45 print_header_linear() 

46 

47 # This is the initialization loop. The requirement is that the 

48 # least-squares solution on free variables is feasible before BVLS starts. 

49 # One possible initialization is to set all variables to lower or upper 

50 # bounds, but many iterations may be required from this state later on. 

51 # The implemented ad-hoc procedure which intuitively should give a better 

52 # initial state: find the least-squares solution on current free variables, 

53 # if its feasible then stop, otherwise, set violating variables to 

54 # corresponding bounds and continue on the reduced set of free variables. 

55 

56 while free_set.size > 0: 

57 if verbose == 2: 

58 optimality = compute_kkt_optimality(g, on_bound) 

59 print_iteration_linear(iteration, cost, cost_change, step_norm, 

60 optimality) 

61 

62 iteration += 1 

63 x_free_old = x[free_set].copy() 

64 

65 A_free = A[:, free_set] 

66 b_free = b - A.dot(x * active_set) 

67 z = lstsq(A_free, b_free, rcond=rcond)[0] 

68 

69 lbv = z < lb[free_set] 

70 ubv = z > ub[free_set] 

71 v = lbv | ubv 

72 

73 if np.any(lbv): 

74 ind = free_set[lbv] 

75 x[ind] = lb[ind] 

76 active_set[ind] = True 

77 on_bound[ind] = -1 

78 

79 if np.any(ubv): 

80 ind = free_set[ubv] 

81 x[ind] = ub[ind] 

82 active_set[ind] = True 

83 on_bound[ind] = 1 

84 

85 ind = free_set[~v] 

86 x[ind] = z[~v] 

87 

88 r = A.dot(x) - b 

89 cost_new = 0.5 * np.dot(r, r) 

90 cost_change = cost - cost_new 

91 cost = cost_new 

92 g = A.T.dot(r) 

93 step_norm = norm(x[free_set] - x_free_old) 

94 

95 if np.any(v): 

96 free_set = free_set[~v] 

97 else: 

98 break 

99 

100 if max_iter is None: 

101 max_iter = n 

102 max_iter += iteration 

103 

104 termination_status = None 

105 

106 # Main BVLS loop. 

107 

108 optimality = compute_kkt_optimality(g, on_bound) 

109 for iteration in range(iteration, max_iter): # BVLS Loop A 

110 if verbose == 2: 

111 print_iteration_linear(iteration, cost, cost_change, 

112 step_norm, optimality) 

113 

114 if optimality < tol: 

115 termination_status = 1 

116 

117 if termination_status is not None: 

118 break 

119 

120 move_to_free = np.argmax(g * on_bound) 

121 on_bound[move_to_free] = 0 

122 

123 while True: # BVLS Loop B 

124 

125 free_set = on_bound == 0 

126 active_set = ~free_set 

127 free_set, = np.nonzero(free_set) 

128 

129 x_free = x[free_set] 

130 x_free_old = x_free.copy() 

131 lb_free = lb[free_set] 

132 ub_free = ub[free_set] 

133 

134 A_free = A[:, free_set] 

135 b_free = b - A.dot(x * active_set) 

136 z = lstsq(A_free, b_free, rcond=rcond)[0] 

137 

138 lbv, = np.nonzero(z < lb_free) 

139 ubv, = np.nonzero(z > ub_free) 

140 v = np.hstack((lbv, ubv)) 

141 

142 if v.size > 0: 

143 alphas = np.hstack(( 

144 lb_free[lbv] - x_free[lbv], 

145 ub_free[ubv] - x_free[ubv])) / (z[v] - x_free[v]) 

146 

147 i = np.argmin(alphas) 

148 i_free = v[i] 

149 alpha = alphas[i] 

150 

151 x_free *= 1 - alpha 

152 x_free += alpha * z 

153 x[free_set] = x_free 

154 

155 if i < lbv.size: 

156 on_bound[free_set[i_free]] = -1 

157 else: 

158 on_bound[free_set[i_free]] = 1 

159 else: 

160 x_free = z 

161 x[free_set] = x_free 

162 break 

163 

164 step_norm = norm(x_free - x_free_old) 

165 

166 r = A.dot(x) - b 

167 cost_new = 0.5 * np.dot(r, r) 

168 cost_change = cost - cost_new 

169 

170 if cost_change < tol * cost: 

171 termination_status = 2 

172 cost = cost_new 

173 

174 g = A.T.dot(r) 

175 optimality = compute_kkt_optimality(g, on_bound) 

176 

177 if termination_status is None: 

178 termination_status = 0 

179 

180 return OptimizeResult( 

181 x=x, fun=r, cost=cost, optimality=optimality, active_mask=on_bound, 

182 nit=iteration + 1, status=termination_status, 

183 initial_cost=initial_cost)