Coverage for gpkit\solvers\cvxopt.py : 0%
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
8# pylint: disable=too-many-locals
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 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)