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

58 statements

, created at 2022-07-28 12:35 -0400

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

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

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

10 """Interface to the CVXOPT solver

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

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

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 A = A.tocsr()

43 maxcol = A.shape[1]-1

44 lse_mons, lin_mons, leq_mons = [], [], []

45 lse_posys, lin_posys, leq_posys = [], [], []

46 constraint_hashes = set()

47 for i, n_monomials in enumerate(k):

48 start = sum(k[:i])

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

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

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

52 if chash in constraint_hashes:

53 continue # already got it

54 if i: # skip cost posy

56 if use_leqs and start in meq_idxs.all:

57 if start in meq_idxs.first_half:

58 leq_posys.append(i)

59 leq_mons.extend(mons)

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

61 lin_mons.extend(mons)

62 lin_posys.append(i)

63 else:

64 lse_mons.extend(mons)

65 lse_posys.append(i)

66 if leq_mons:

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

68 log_c_leq = log_c[leq_mons]

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

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

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

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

73 if lin_mons:

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

75 log_c_lin = log_c[lin_mons]

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

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

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

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

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

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

82 log_c_lse = log_c[lse_mons]

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

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

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

86 g = matrix(log_c_lse)

87 try:

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

89 except ValueError as e:

90 raise DualInfeasible() from e

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

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

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

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

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

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

97 if yi >= 0:

98 la[leq_posy] = yi

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

100 la[leq_posy+1] = -yi

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

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

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

104 la=la)