Coverage for gpkit/solvers/cvxopt.py: 100%
58 statements
« prev ^ index » next coverage.py v7.4.0, created at 2024-01-05 22:09 -0500
« prev ^ index » next coverage.py v7.4.0, created at 2024-01-05 22:09 -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
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
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
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
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}