Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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 

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

10 """Interface to the CVXOPT solver 

11 

12 Definitions 

13 ----------- 

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

15 n is the number of monomials in the gp 

16 m is the number of variables in the gp 

17 p is the number of posynomial constraints in the gp 

18 

19 Arguments 

20 --------- 

21 c : floats array of shape n 

22 Coefficients of each monomial 

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

24 Exponents of the various free variables for each monomial. 

25 k : ints array of shape p+1 

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

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

28 

29 Returns 

30 ------- 

31 dict 

32 Contains the following keys 

33 "success": bool 

34 "objective_sol" float 

35 Optimal value of the objective 

36 "primal_sol": floats array of size m 

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

38 "dual_sol": floats array of size p 

39 Optimal value of the dual variables, in logspace. 

40 """ 

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

42 lse_mons, lin_mons, leq_mons = [], [], [] 

43 lse_posys, lin_posys, leq_posys = [], [], [] 

44 for i, n_monomials in enumerate(k): 

45 start = sum(k[:i]) 

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

47 if use_leqs and start in meq_idxs.all: 

48 if start in meq_idxs.first_half: 

49 leq_posys.append(i) 

50 leq_mons.extend(mons) 

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

52 lin_mons.extend(mons) 

53 lin_posys.append(i) 

54 else: 

55 lse_mons.extend(mons) 

56 lse_posys.append(i) 

57 A = A.tocsr() 

58 maxcol = A.shape[1]-1 

59 if leq_mons: 

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

61 log_c_leq = log_c[leq_mons] 

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

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

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

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

66 if lin_mons: 

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

68 log_c_lin = log_c[lin_mons] 

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

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

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

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

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

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

75 log_c_lse = log_c[lse_mons] 

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

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

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

79 g = matrix(log_c_lse) 

80 try: 

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

82 except ValueError as e: 

83 raise DualInfeasible() from e 

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

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

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

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

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

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

90 if yi >= 0: 

91 la[leq_posy] = yi 

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

93 la[leq_posy+1] = -yi 

94 return dict(status=solution["status"], 

95 objective=np.exp(solution["primal objective"]), 

96 primal=np.ravel(solution["x"]), 

97 la=la)