Coverage for gpkit/solvers/cvxopt.py: 100%

58 statements  

« prev     ^ index     » next       coverage.py v7.4.0, created at 2024-01-07 22:56 -0500

1"Implements the GPkit interface to CVXOPT" 

2import numpy as np 

3from cvxopt import spmatrix, matrix 

4from cvxopt.solvers import gp 

5from gpkit.exceptions import UnknownInfeasible, DualInfeasible 

6 

7 

8# pylint: disable=too-many-locals,too-many-statements,too-many-branches 

9# pylint: disable=invalid-name 

10def optimize(*, c, A, k, meq_idxs, use_leqs=True, **kwargs): 

11 """Interface to the CVXOPT solver 

12 

13 Definitions 

14 ----------- 

15 "[a,b] array of floats" indicates array-like data with shape [a,b] 

16 n is the number of monomials in the gp 

17 m is the number of variables in the gp 

18 p is the number of posynomial constraints in the gp 

19 

20 Arguments 

21 --------- 

22 c : floats array of shape n 

23 Coefficients of each monomial 

24 A : floats array of shape (n, m) 

25 Exponents of the various free variables for each monomial. 

26 k : ints array of shape p+1 

27 k[0] is the number of monomials (rows of A) present in the objective 

28 k[1:] is the number of monomials present in each constraint 

29 

30 Returns 

31 ------- 

32 dict 

33 Contains the following keys 

34 "success": bool 

35 "objective_sol" float 

36 Optimal value of the objective 

37 "primal_sol": floats array of size m 

38 Optimal value of free variables. Note: not in logspace. 

39 "dual_sol": floats array of size p 

40 Optimal value of the dual variables, in logspace. 

41 """ 

42 log_c = np.log(np.array(c)) 

43 A = A.tocsr() 

44 maxcol = A.shape[1]-1 

45 lse_mons, lin_mons, leq_mons = [], [], [] 

46 lse_posys, lin_posys, leq_posys = [], [], [] 

47 constraint_hashes = set() 

48 for i, n_monomials in enumerate(k): 

49 start = sum(k[:i]) 

50 mons = range(start, start+k[i]) 

51 A_m = A[mons, :].tocoo() 

52 chash = hash((c[i], tuple(A_m.data), tuple(A_m.row), tuple(A_m.col))) 

53 if chash in constraint_hashes: 

54 continue # already got it 

55 if i: # skip cost posy 

56 constraint_hashes.add(chash) 

57 if use_leqs and start in meq_idxs.all: 

58 if start in meq_idxs.first_half: 

59 leq_posys.append(i) 

60 leq_mons.extend(mons) 

61 elif i != 0 and n_monomials == 1: 

62 lin_mons.extend(mons) 

63 lin_posys.append(i) 

64 else: 

65 lse_mons.extend(mons) 

66 lse_posys.append(i) 

67 if leq_mons: 

68 A_leq = A[leq_mons, :].tocoo() 

69 log_c_leq = log_c[leq_mons] 

70 kwargs["A"] = spmatrix([float(r) for r in A_leq.data]+[0], 

71 [int(r) for r in A_leq.row]+[0], 

72 [int(r) for r in A_leq.col]+[maxcol], tc="d") 

73 kwargs["b"] = matrix(-log_c_leq) 

74 if lin_mons: 

75 A_lin = A[lin_mons, :].tocoo() 

76 log_c_lin = log_c[lin_mons] 

77 kwargs["G"] = spmatrix([float(r) for r in A_lin.data]+[0], 

78 [int(r) for r in A_lin.row]+[0], 

79 [int(r) for r in A_lin.col]+[maxcol], tc="d") 

80 kwargs["h"] = matrix(-log_c_lin) 

81 k_lse = [k[i] for i in lse_posys] 

82 A_lse = A[lse_mons, :].tocoo() 

83 log_c_lse = log_c[lse_mons] 

84 F = spmatrix([float(r) for r in A_lse.data]+[0], 

85 [int(r) for r in A_lse.row]+[0], 

86 [int(r) for r in A_lse.col]+[maxcol], tc="d") 

87 g = matrix(log_c_lse) 

88 try: 

89 solution = gp(k_lse, F, g, **kwargs) 

90 except ValueError as e: 

91 raise DualInfeasible() from e 

92 if solution["status"] != "optimal": 

93 raise UnknownInfeasible("solution status " + repr(solution["status"])) 

94 la = np.zeros(len(k)) 

95 la[lin_posys] = list(solution["zl"]) 

96 la[lse_posys] = [1.] + list(solution["znl"]) 

97 for leq_posy, yi in zip(leq_posys, solution["y"]): 

98 if yi >= 0: 

99 la[leq_posy] = yi 

100 else: # flip it around to the other "inequality" 

101 la[leq_posy+1] = -yi 

102 return {"status": solution["status"], 

103 "objective": np.exp(solution["primal objective"]), 

104 "primal": np.ravel(solution["x"]), 

105 "la": la}