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

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 """

11 Definitions

12 -----------

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

14 n is the number of monomials in the gp

15 m is the number of variables in the gp

16 p is the number of posynomial constraints in the gp

18 Arguments

19 ---------

20 c : floats array of shape n

21 Coefficients of each monomial

22 A : gpkit.small_classes.CootMatrix, of shape (n, m)

23 Exponents of the various free variables for each monomial.

24 k : ints array of shape p+1

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

26 k[1:] is the number of monomials (rows of A) present in each constraint

27 p_idxs : ints array of shape n.

28 sel = p_idxs == i selects rows of A

29 and entries of c of the i-th posynomial

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

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

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

34 Returns

35 -------

36 dict

37 Contains the following keys

38 "status": string

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

40 "primal" np.ndarray or None

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

42 "la": np.ndarray or None

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

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

45 """

46 #

47 # Initial transformations of problem data.

48 #

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

50 # from those which require

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

52 #

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

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

55 #

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

57 lse_posys = [0]

58 lin_posys = []

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

60 if val > 1:

61 lse_posys.append(i+1)

62 else:

63 lin_posys.append(i+1)

64 if lin_posys:

65 A = A.tocsr()

66 lin_posys_set = frozenset(lin_posys)

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

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

69 lse_idxs[lin_idxs] = False

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

71 log_c_lse = log_c[lse_idxs]

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

73 log_c_lin = log_c[lin_idxs]

74 else:

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

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

77 log_c_lse = log_c

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

79 n_lse = sum(k_lse)

80 p_lse = len(k_lse)

81 lse_p_idx = []

82 for i, ki in enumerate(k_lse):

83 lse_p_idx.extend([i] * ki)

84 lse_p_idx = np.array(lse_p_idx)

85 #

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

87 # constraints.

88 #

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

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

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

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

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

94 #

95 # The variable z[0] is special,

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

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

98 #

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

100 # in constraints that naturally write as

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

102 #

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

104 # in next comment block.

105 #

106 env = mosek.Env()

107 task = env.Task(0, 0)

108 m = A.shape[1]

109 msk_nvars = m + 3 * n_lse + p_lse

110 task.appendvars(msk_nvars)

111 # "x" is free

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

113 np.zeros(m), np.zeros(m))

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

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

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

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

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

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

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

121 np.zeros(p_lse), np.zeros(p_lse))

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

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

124 [3] * n_lse, m)

125 #

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

127 #

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

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

130 #

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

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

133 #

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

135 # above constraints imply

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

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

138 #

139 task.appendcons(n_lse + p_lse)

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

141 # start with coefficients on "x"

142 rows = list(A_lse.row)

143 cols = list(A_lse.col)

144 vals = list(A_lse.data)

145 # add coefficients on "t"

146 rows += list(range(n_lse))

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

148 vals += [-1.0] * n_lse

149 # add coefficients on "z"

150 rows += list(range(n_lse))

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

152 vals += [-1.0] * n_lse

153 task.putaijlist(rows, cols, vals)

154 cur_con_idx = n_lse

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

156 rows, cols, vals = [], [], []

157 for i in range(p_lse):

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

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

160 cols.extend(m + 3 * sels)

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

162 cur_con_idx += 1

163 task.putaijlist(rows, cols, vals)

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

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

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

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

168 #

169 # Affine constraints, not needing the exponential cone

170 #

171 # Require A_lin @ x <= -log_c_lin.

172 #

173 if log_c_lin is not None:

174 task.appendcons(log_c_lin.size)

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

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

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

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

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

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

181 cur_con_idx += log_c_lin.size

182 #

183 # Constrain choice variables to discrete choices

184 #

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

186 if choicevaridxs:

187 n_choicevars = 0

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

189 choices = sorted(var.choices)

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

191 choiceidxs = np.arange(msk_nvars + n_choicevars,

192 msk_nvars + n_choicevars + m_choices)

193 n_choicevars += m_choices

194 task.appendvars(m_choices)

195 task.putvartypelist(choiceidxs,

196 [mosek.variabletype.type_int]*m_choices)

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

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

199 task.appendcons(m_choices)

200 for i in range(m_choices - 1):

201 # each larger choice requires those before it

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

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

204 cur_con_idx += 1

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

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

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

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

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

210 cur_con_idx += 1

211 #

212 # Set the objective function

213 #

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

215 task.putobjsense(mosek.objsense.minimize)

216 #

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

218 #

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

220 if verbose:

221 def streamprinter(text):

222 "Stream printer for output from mosek."

223 print(text)

224

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

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

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

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

229

230 try:

231 task.optimize()

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

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

234 mosek.rescode.err_license_version,

235 mosek.rescode.err_license_expired]:

236 raise InvalidLicense() from e

237 raise e

238

239 if verbose:

240 task.solutionsummary(mosek.streamtype.msg)

241 #

242 # Recover the solution

243 #

244 if choicevaridxs:

245 sol = mosek.soltype.itg

246 optimal = mosek.solsta.integer_optimal

247 else:

248 sol = mosek.soltype.itr

249 optimal = mosek.solsta.optimal

250 msk_solsta = task.getsolsta(sol)

251 if msk_solsta == mosek.solsta.prim_infeas_cer:

252 raise PrimalInfeasible()

253 if msk_solsta == mosek.solsta.dual_infeas_cer:

254 raise DualInfeasible()

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

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

257

258 # recover primal variables

259 x = [0.] * m

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

261 # recover binary variables

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

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

264 # wrap things up in a dictionary

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

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

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

268 # (skip epigraph of the objective function).

269 if choicevaridxs: # no dual solution

270 solution["la"] = []

271 solution["nu"] = []

272 else:

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

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

275 z_duals = np.array(z_duals)

276 z_duals[z_duals < 0] = 0

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

278 if log_c_lin is None:

279 solution["la"] = z_duals

280 else:

281 aff_duals = [0.] * log_c_lin.size

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

283 aff_duals)

284 aff_duals = np.array(aff_duals)

285 aff_duals[aff_duals < 0] = 0

286 # merge z_duals with aff_duals

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

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

289 merged_duals[lin_posys] = aff_duals

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

291

292 task.__exit__(None, None, None)

293 env.__exit__(None, None, None)

294 return solution

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

142 statements