Coverage for gpkit/solvers/mosek

1"""Implements the GPkit interface to MOSEK (version >= 9)

2 through python-based Optimizer API"""

3import numpy as np

4import mosek # pylint: disable=import-error

5from ..exceptions import (UnknownInfeasible, InvalidLicense,

6 PrimalInfeasible, DualInfeasible)

8def optimize(*, c, A, k, p_idxs, **kwargs):

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

10 # pylint: disable=invalid-name

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

19 Arguments

20 ---------

21 c : floats array of shape n

22 Coefficients of each monomial

23 A : gpkit.small_classes.CootMatrix, 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 (rows of A) present in each constraint

28 p_idxs : ints array of shape n.

29 sel = p_idxs == i selects rows of A

30 and entries of c of the i-th posynomial

31 fi(x) = c[sel] @ exp(A[sel,:] @ x).

32 The 0-th posynomial gives the objective function, and the

33 remaining posynomials should be constrained to be <= 1.

35 Returns

36 -------

37 dict

38 Contains the following keys

39 "status": string

40 "optimal", "infeasible", "unbounded", or "unknown".

41 "primal" np.ndarray or None

42 The values of the ``m`` primal variables.

43 "la": np.ndarray or None

44 The dual variables to the ``p`` posynomial constraints, when

45 those constraints are represented in log-sum-exp ("LSE") form.

46 """

47 #

48 # Initial transformations of problem data.

49 #

50 # separate monomial constraints (call them "lin"),

51 # from those which require

52 # an LSE representation (call those "lse").

53 #

54 # NOTE: the epigraph of the objective function always gets an "lse"

55 # representation, even if the objective is a monomial.

56 #

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

58 lse_posys = [0]

59 lin_posys = []

60 for i, val in enumerate(k[1:]):

61 if val > 1:

62 lse_posys.append(i+1)

63 else:

64 lin_posys.append(i+1)

65 if lin_posys:

66 A = A.tocsr()

67 lin_posys_set = frozenset(lin_posys)

68 lin_idxs = [i for i, p in enumerate(p_idxs) if p in lin_posys_set]

69 lse_idxs = np.ones(A.shape[0], dtype=bool)

70 lse_idxs[lin_idxs] = False

71 A_lse = A[lse_idxs, :].tocoo()

72 log_c_lse = log_c[lse_idxs]

73 A_lin = A[lin_idxs, :].tocoo()

74 log_c_lin = log_c[lin_idxs]

75 else:

76 log_c_lin = None # A_lin won't be referenced later,

77 A_lse = A # so no need to define it.

78 log_c_lse = log_c

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

80 n_lse = sum(k_lse)

81 p_lse = len(k_lse)

82 lse_p_idx = []

83 for i, ki in enumerate(k_lse):

84 lse_p_idx.extend([i] * ki)

85 lse_p_idx = np.array(lse_p_idx)

86 #

87 # Create MOSEK task. Add variables, bound constraints, and conic

88 # constraints.

89 #

90 # Say that MOSEK's optimization variable is a block vector,

91 # [x;t;z], where ...

92 # x is the user-defined primal variable (length m),

93 # t is an aux variable for exponential cones (length 3 * n_lse),

94 # z is an epigraph variable for LSE terms (length p_lse).

95 #

96 # The variable z[0] is special,

97 # because it's the epigraph of the objective function

98 # in LSE form. The sign of this variable is not constrained.

99 #

100 # The variables z[1:] are epigraph terms for "log",

101 # in constraints that naturally write as

102 # LSE(Ai @ x + log_ci) <= 0. These variables need to be <= 0.

103 #

104 # The main constraints on (x, t, z) are described

105 # in next comment block.

106 #

107 env = mosek.Env()

108 task = env.Task(0, 0)

109 m = A.shape[1]

110 msk_nvars = m + 3 * n_lse + p_lse

111 task.appendvars(msk_nvars)

112 # "x" is free

113 task.putvarboundlist(np.arange(m), [mosek.boundkey.fr] * m,

114 np.zeros(m), np.zeros(m))

115 # t[3 * i + i] == 1, other components are free.

116 bound_types = [mosek.boundkey.fr, mosek.boundkey.fx, mosek.boundkey.fr]

117 task.putvarboundlist(np.arange(m, m + 3*n_lse), bound_types * n_lse,

118 np.ones(3*n_lse), np.ones(3*n_lse))

119 # z[0] is free; z[1:] <= 0.

120 bound_types = [mosek.boundkey.fr] + [mosek.boundkey.up] * (p_lse - 1)

121 task.putvarboundlist(np.arange(m + 3*n_lse, msk_nvars), bound_types,

122 np.zeros(p_lse), np.zeros(p_lse))

123 # t[3*i], t[3*i + 1], t[3*i + 2] belongs to the exponential cone

124 task.appendconesseq([mosek.conetype.pexp] * n_lse, [0.0] * n_lse,

125 [3] * n_lse, m)

126 #

127 # Exponential cone affine constraints (other than t[3*i + 1] == 1).

128 #

129 # For each i in {0, ..., n_lse - 1}, we need

130 # t[3*i + 2] == A_lse[i, :] @ x + log_c_lse[i] - z[lse_p_idx[i]].

131 #

132 # For each j from {0, ..., p_lse - 1}, the "t" should also satisfy

133 # sum(t[3*i] for i where i == lse_p_idx[j]) <= 1.

134 #

135 # When combined with bound constraints ``t[3*i + 1] == 1``, the

136 # above constraints imply

137 # LSE(A_lse[sel, :] @ x + log_c_lse[sel]) <= z[i]

138 # for ``sel = lse_p_idx == i``.

139 #

140 task.appendcons(n_lse + p_lse)

141 # Linear equations between (x,t,z).

142 # start with coefficients on "x"

143 rows = list(A_lse.row)

144 cols = list(A_lse.col)

145 vals = list(A_lse.data)

146 # add coefficients on "t"

147 rows += list(range(n_lse))

148 cols += (m + 3*np.arange(n_lse) + 2).tolist()

149 vals += [-1.0] * n_lse

150 # add coefficients on "z"

151 rows += list(range(n_lse))

152 cols += [m + 3*n_lse + lse_p_idx[i] for i in range(n_lse)]

153 vals += [-1.0] * n_lse

154 task.putaijlist(rows, cols, vals)

155 cur_con_idx = n_lse

156 # Linear inequalities on certain sums of "t".

157 rows, cols, vals = [], [], []

158 for i in range(p_lse):

159 sels = np.nonzero(lse_p_idx == i)[0]

160 rows.extend([cur_con_idx] * sels.size)

161 cols.extend(m + 3 * sels)

162 vals.extend([1] * sels.size)

163 cur_con_idx += 1

164 task.putaijlist(rows, cols, vals)

165 # Build the right-hand-sides of the [in]equality constraints

166 type_constraint = [mosek.boundkey.fx] * n_lse + [mosek.boundkey.up] * p_lse

167 h = np.concatenate([-log_c_lse, np.ones(p_lse)])

168 task.putconboundlist(np.arange(h.size, dtype=int), type_constraint, h, h)

169 #

170 # Affine constraints, not needing the exponential cone

171 #

172 # Require A_lin @ x <= -log_c_lin.

173 #

174 if log_c_lin is not None:

175 task.appendcons(log_c_lin.size)

176 rows = cur_con_idx + np.array(A_lin.row)

177 task.putaijlist(rows, A_lin.col, A_lin.data)

178 type_constraint = [mosek.boundkey.up] * log_c_lin.size

179 con_indices = np.arange(cur_con_idx, cur_con_idx + log_c_lin.size)

180 h = -log_c_lin #pylint: disable=invalid-unary-operand-type

181 task.putconboundlist(con_indices, type_constraint, h, h)

182 cur_con_idx += log_c_lin.size

183 #

184 # Constrain choice variables to discrete choices

185 #

186 choicevaridxs = kwargs.get("choicevaridxs", {})

187 if choicevaridxs:

188 n_choicevars = 0

189 for var, idx in choicevaridxs.items():

190 choices = sorted(var.choices)

191 m_choices = len(choices) - 1 # the first option is the default

192 choiceidxs = np.arange(msk_nvars + n_choicevars,

193 msk_nvars + n_choicevars + m_choices)

194 n_choicevars += m_choices

195 task.appendvars(m_choices)

196 task.putvartypelist(choiceidxs,

197 [mosek.variabletype.type_int]*m_choices)

198 task.putvarboundlist(choiceidxs, [mosek.boundkey.ra]*m_choices,

199 np.zeros(m_choices), np.ones(m_choices))

200 task.appendcons(m_choices)

201 for i in range(m_choices - 1):

202 # each larger choice requires those before it

203 task.putarow(cur_con_idx + i, choiceidxs[i:i+2], [1.0, -1.0])

204 task.putconbound(cur_con_idx + i, mosek.boundkey.lo, 0.0, 0.0)

205 cur_con_idx += 1

206 base = np.log(choices[0])

207 logdiffs = np.diff(np.log(choices)).tolist()

208 task.putarow(cur_con_idx, choiceidxs.tolist() + [idx],

209 logdiffs + [-1]) # choices are summed

210 task.putconbound(cur_con_idx, mosek.boundkey.fx, -base, -base)

211 cur_con_idx += 1

212 #

213 # Set the objective function

214 #

215 task.putclist([int(m + 3*n_lse)], [1])

216 task.putobjsense(mosek.objsense.minimize)

217 #

218 # Set solver parameters, and call .solve().

219 #

220 verbose = kwargs.get("verbose", True)

221 if verbose:

222 def streamprinter(text):

223 "Stream printer for output from mosek."

224 print(text)

225

226 env.set_Stream(mosek.streamtype.log, streamprinter)

227 task.set_Stream(mosek.streamtype.log, streamprinter)

228 task.putintparam(mosek.iparam.infeas_report_auto, mosek.onoffkey.on)

229 task.putintparam(mosek.iparam.log_presolve, 0)

230

231 try:

232 task.optimize()

233 except mosek.Error as e: # pragma: no cover

234 if e.errno in [mosek.rescode.err_missing_license_file,

235 mosek.rescode.err_license_version,

236 mosek.rescode.err_license_expired]:

237 raise InvalidLicense() from e

238 raise e

239

240 if verbose:

241 task.solutionsummary(mosek.streamtype.msg)

242 #

243 # Recover the solution

244 #

245 if choicevaridxs:

246 sol = mosek.soltype.itg

247 optimal = mosek.solsta.integer_optimal

248 else:

249 sol = mosek.soltype.itr

250 optimal = mosek.solsta.optimal

251 msk_solsta = task.getsolsta(sol)

252 if msk_solsta == mosek.solsta.prim_infeas_cer:

253 raise PrimalInfeasible()

254 if msk_solsta == mosek.solsta.dual_infeas_cer:

255 raise DualInfeasible()

256 if msk_solsta != optimal: # pragma: no cover

257 raise UnknownInfeasible("solution status: ", msk_solsta)

258

259 # recover primal variables

260 x = [0.] * m

261 task.getxxslice(sol, 0, m, x)

262 # recover binary variables

263 # xbin = [0.] * (n_choicevars)

264 # task.getxxslice(sol, msk_nvars, msk_nvars + n_choicevars, xbin)

265 # wrap things up in a dictionary

266 solution = {"status": "optimal", "primal": np.array(x),

267 "objective": np.exp(task.getprimalobj(sol))}

268 # recover dual variables for log-sum-exp epigraph constraints

269 # (skip epigraph of the objective function).

270 if choicevaridxs: # no dual solution

271 solution["la"] = []

272 solution["nu"] = []

273 else:

274 z_duals = [0.] * (p_lse - 1)

275 task.getsuxslice(mosek.soltype.itr, m + 3*n_lse + 1, msk_nvars, z_duals)

276 z_duals = np.array(z_duals)

277 z_duals[z_duals < 0] = 0

278 # recover dual variables for the remaining user-provided constraints

279 if log_c_lin is None:

280 solution["la"] = z_duals

281 else:

282 aff_duals = [0.] * log_c_lin.size

283 task.getsucslice(mosek.soltype.itr, n_lse + p_lse, cur_con_idx,

284 aff_duals)

285 aff_duals = np.array(aff_duals)

286 aff_duals[aff_duals < 0] = 0

287 # merge z_duals with aff_duals

288 merged_duals = np.zeros(len(k))

289 merged_duals[lse_posys[1:]] = z_duals # skipping the cost

290 merged_duals[lin_posys] = aff_duals

291 solution["la"] = merged_duals[1:]

292

293 task.__exit__(None, None, None)

294 env.__exit__(None, None, None)

295 return solution

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

142 statements