/
mosekglue.py
275 lines (232 loc) · 10.2 KB
/
mosekglue.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
"""
Copyright 2015 Enzo Busseti
[ Modified by Balasubramanian Narasimhan ]
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
"""
import settings as s
import numpy as np
import scipy.sparse as sp
import mosek
def mosek_intf(A, b, G, h, c, dims, offset, solver_opts, verbose = False):
data = {s.A : A,
s.B : b,
s.G : G,
s.H : h,
s.C : c,
s.OFFSET : offset,
s.DIMS : dims}
with mosek.Env() as env:
with env.Task(0, 0) as task:
kwargs = sorted(solver_opts.keys())
if "mosek_params" in kwargs:
self._handle_mosek_params(task, solver_opts["mosek_params"])
kwargs.remove("mosek_params")
if kwargs:
raise ValueError("Invalid keyword-argument '%s'" % kwargs[0])
if verbose:
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
import sys
sys.stdout.write(text)
sys.stdout.flush()
env.set_Stream(mosek.streamtype.log, streamprinter)
task.set_Stream(mosek.streamtype.log, streamprinter)
# size of problem
numvar = len(c) + sum(dims[s.SOC_DIM])
numcon = len(b) + dims[s.LEQ_DIM] + sum(dims[s.SOC_DIM]) + \
sum([el**2 for el in dims[s.SDP_DIM]])
# otherwise it crashes on empty probl.
if numvar == 0:
result_dict = {s.STATUS: s.OPTIMAL}
result_dict[s.PRIMAL] = []
result_dict[s.VALUE] = 0. + data[s.OFFSET]
result_dict[s.EQ_DUAL] = []
result_dict[s.INEQ_DUAL] = []
return result_dict
# objective
task.appendvars(numvar)
task.putclist(np.arange(len(c)), c)
task.putvarboundlist(np.arange(numvar, dtype=int),
[mosek.boundkey.fr]*numvar,
np.zeros(numvar),
np.zeros(numvar))
# SDP variables
if sum(dims[s.SDP_DIM]) > 0:
task.appendbarvars(dims[s.SDP_DIM])
# linear equality and linear inequality constraints
task.appendcons(numcon)
if A.shape[0] and G.shape[0]:
constraints_matrix = sp.bmat([[A], [G]])
else:
constraints_matrix = A if A.shape[0] else G
coefficients = np.concatenate([b, h])
row, col, el = sp.find(constraints_matrix)
task.putaijlist(row, col, el)
type_constraint = [mosek.boundkey.fx] * len(b)
type_constraint += [mosek.boundkey.up] * dims[s.LEQ_DIM]
sdp_total_dims = sum([cdim**2 for cdim in dims[s.SDP_DIM]])
type_constraint += [mosek.boundkey.fx] * \
(sum(dims[s.SOC_DIM]) + sdp_total_dims)
task.putconboundlist(np.arange(numcon, dtype=int),
type_constraint,
coefficients,
coefficients)
# cones
current_var_index = len(c)
current_con_index = len(b) + dims[s.LEQ_DIM]
for size_cone in dims[s.SOC_DIM]:
row, col, el = sp.find(sp.eye(size_cone))
row += current_con_index
col += current_var_index
task.putaijlist(row, col, el) # add a identity for each cone
# add a cone constraint
task.appendcone(mosek.conetype.quad,
0.0, # unused
np.arange(current_var_index,
current_var_index + size_cone))
current_con_index += size_cone
current_var_index += size_cone
# SDP
for num_sdp_var, size_matrix in enumerate(dims[s.SDP_DIM]):
for i_sdp_matrix in range(size_matrix):
for j_sdp_matrix in range(size_matrix):
coeff = 1. if i_sdp_matrix == j_sdp_matrix else .5
task.putbaraij(current_con_index,
num_sdp_var,
[task.appendsparsesymmat(size_matrix,
[max(i_sdp_matrix,
j_sdp_matrix)],
[min(i_sdp_matrix,
j_sdp_matrix)],
[coeff])],
[1.0])
current_con_index += 1
# solve
task.putobjsense(mosek.objsense.minimize)
task.optimize()
if verbose:
task.solutionsummary(mosek.streamtype.msg)
return format_results(task, data)
def choose_solution(task):
"""Chooses between the basic and interior point solution.
Parameters
----------
task : mosek.Task
The solver status interface.
Returns
-------
soltype
The preferred solution (mosek.soltype.*)
solsta
The status of the preferred solution (mosek.solsta.*)
"""
import mosek
def rank(status):
# Rank solutions
# optimal > near_optimal > anything else > None
if status == mosek.solsta.optimal:
return 3
elif status == mosek.solsta.near_optimal:
return 2
elif status is not None:
return 1
else:
return 0
solsta_bas, solsta_itr = None, None
if task.solutiondef(mosek.soltype.bas):
solsta_bas = task.getsolsta(mosek.soltype.bas)
if task.solutiondef(mosek.soltype.itr):
solsta_itr = task.getsolsta(mosek.soltype.itr)
# As long as interior solution is not worse, take it
# (for backward compatibility)
if rank(solsta_itr) >= rank(solsta_bas):
return mosek.soltype.itr, solsta_itr
else:
return mosek.soltype.bas, solsta_bas
def format_results(task, data):
"""Converts the solver output into standard form.
Parameters
----------
task : mosek.Task
The solver status interface.
data : dict
Information about the problem.
Returns
-------
dict
The solver output in standard form.
"""
import mosek
# Map of MOSEK status to CVXPY status.
# taken from:
# http://docs.mosek.com/7.0/pythonapi/Solution_status_keys.html
STATUS_MAP = {mosek.solsta.optimal: s.OPTIMAL,
mosek.solsta.prim_infeas_cer: s.INFEASIBLE,
mosek.solsta.dual_infeas_cer: s.UNBOUNDED,
mosek.solsta.near_optimal: s.OPTIMAL_INACCURATE,
mosek.solsta.near_prim_infeas_cer: s.INFEASIBLE_INACCURATE,
mosek.solsta.near_dual_infeas_cer: s.UNBOUNDED_INACCURATE,
mosek.solsta.unknown: s.SOLVER_ERROR}
soltype, solsta = choose_solution(task)
if solsta in STATUS_MAP:
result_dict = {s.STATUS: STATUS_MAP[solsta]}
else:
result_dict = {s.STATUS: s.SOLVER_ERROR}
# Callback data example:
# http://docs.mosek.com/7.1/pythonapi/The_progress_call-back.html
# Retrieving double information items:
# http://docs.mosek.com/7.1/pythonapi/Task_getdouinf_.html#@generated-ID:5ef16e0
# http://docs.mosek.com/7.1/pythonapi/Double_information_items.html
result_dict[s.SOLVE_TIME] = task.getdouinf(mosek.dinfitem.optimizer_time)
result_dict[s.SETUP_TIME] = task.getdouinf(mosek.dinfitem.presolve_time)
result_dict[s.NUM_ITERS] = task.getintinf(mosek.iinfitem.intpnt_iter)
if result_dict[s.STATUS] in s.SOLUTION_PRESENT:
# get primal variables values
result_dict[s.PRIMAL] = np.zeros(task.getnumvar(), dtype=np.float)
task.getxx(soltype, result_dict[s.PRIMAL])
# get obj value
result_dict[s.VALUE] = task.getprimalobj(soltype) + \
data[s.OFFSET]
# get dual
y = np.zeros(task.getnumcon(), dtype=np.float)
task.gety(soltype, y)
# it appears signs are inverted
result_dict[s.EQ_DUAL] = -y[:len(data[s.B])]
result_dict[s.INEQ_DUAL] = \
-y[len(data[s.B]):len(data[s.B])+data[s.DIMS][s.LEQ_DIM]]
return result_dict
def _handle_mosek_params(task, params):
if params is None:
return
import mosek
def _handle_str_param(param, value):
if param.startswith("MSK_DPAR_"):
task.putnadouparam(param, value)
elif param.startswith("MSK_IPAR_"):
task.putnaintparam(param, value)
elif param.startswith("MSK_SPAR_"):
task.putnastrparam(param, value)
else:
raise ValueError("Invalid MOSEK parameter '%s'." % param)
def _handle_enum_param(param, value):
if isinstance(param, mosek.dparam):
task.putdouparam(param, value)
elif isinstance(param, mosek.iparam):
task.putintparam(param, value)
elif isinstance(param, mosek.sparam):
task.putstrparam(param, value)
else:
raise ValueError("Invalid MOSEK parameter '%s'." % param)
for param, value in params.items():
if isinstance(param, str):
_handle_str_param(param.strip(), value)
else:
_handle_enum_param(param, value)