Coverage for picos/modeling/file_out.py: 2.83%

1# ------------------------------------------------------------------------------ 

2# Copyright (C) 2012-2017 Guillaume Sagnol 

3# Copyright (C) 2019 Maximilian Stahlberg 

4# 

5# This file is part of PICOS. 

6# 

7# PICOS is free software: you can redistribute it and/or modify it under the 

8# terms of the GNU General Public License as published by the Free Software 

9# Foundation, either version 3 of the License, or (at your option) any later 

10# version. 

11# 

12# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY 

13# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 

14# A PARTICULAR PURPOSE. See the GNU General Public License for more details. 

15# 

16# You should have received a copy of the GNU General Public License along with 

17# this program. If not, see <http://www.gnu.org/licenses/>. 

18# ------------------------------------------------------------------------------ 

19 

20"""Functions for writing optimization problems to a file.""" 

21 

22# ------------------------------------------------------------------------------ 

23# TODO: This is a temporary solution using outdated code. Ideally, writing to 

24# a file should be integrated with the reformulation pipeline e.g. by 

25# distinguishing a SolutionStrategy and an ExportStrategy. 

26# ------------------------------------------------------------------------------ 

27 

28from itertools import chain 

29 

30import cvxopt 

31import numpy 

32 

33from ..apidoc import api_end, api_start 

34from ..constraints import ( 

35 AffineConstraint, 

36 LMIConstraint, 

37 RSOCConstraint, 

38 SOCConstraint, 

39) 

40 

41from ..solvers import CVXOPTSolver, CPLEXSolver, MOSEKSolver, GurobiSolver 

42from ..expressions import IntegerVariable, BinaryVariable, CONTINUOUS_VARTYPES 

43from ..expressions.vectorizations import FullVectorization 

44 

45_API_START = api_start(globals()) 

46# ------------------------------- 

47 

48 

49INFINITY = 1e16 #: A number deemed too large to appear in practice. 

50 

51 

52def write(picos_problem, filename, writer="picos"): 

53 r"""Write an optimization problem to a file. 

54 

55 :param ~picos.Problem P: The problem to write. 

56 

57 :param str filename: Path and name of the output file. The export format 

58 is inferred from the file extension. Supported extensions and their 

59 associated format are: 

60 

61 * ``'.cbf'`` -- Conic Benchmark Format. 

62 

63 This format is suitable for optimization problems involving second 

64 order and/or semidefinite cone constraints. This is a standard 

65 choice for conic optimization problems. Visit the website of 

66 `The Conic Benchmark Library <http://cblib.zib.de/>`_ or read 

67 `A benchmark library for conic mixed-integer and continuous 

68 optimization 

69 <http://www.optimization-online.org/DB_HTML/2014/03/4301.html>`_ 

70 by Henrik A. Friberg for more information. 

71 

72 * ``'.lp'`` -- `LP format 

73 <http://docs.mosek.com/6.0/pyapi/node022.html>`_. 

74 

75 This format handles only linear constraints, unless the writer 

76 ``'cplex'`` is used. In the latter case the extended `CPLEX LP format 

77 <https://www.ibm.com/support/knowledgecenter/SSSA5P_12.10.0/ 

78 ilog.odms.cplex.help/CPLEX/FileFormats/topics/LP.html>`_ 

79 is used instead. 

80 

81 * ``'.mps'`` -- `MPS format 

82 <http://docs.mosek.com/6.0/pyapi/node021.html>`_. 

83 

84 As the writer, you need to choose one of ``'cplex'``, ``'gurobi'`` 

85 or ``'mosek'``. 

86 

87 * ``'.opf'`` -- `OPF format 

88 <http://docs.mosek.com/6.0/pyapi/node023.html>`_. 

89 

90 As the writer, you need to choose ``'mosek'``. 

91 

92 * ``'.dat-s'`` -- `Sparse SDPA format 

93 <http://sdpa.indsys.chuo-u.ac.jp/sdpa/download.html#sdpa>`_. 

94 

95 This format is suitable for semidefinite programs. Second order 

96 cone constraints are stored as semidefinite constraints on an 

97 *arrow shaped* matrix. 

98 

99 :param str writer: The default ``'picos'`` denotes PICOS' internal 

100 writer, which can export to *LP*, *CBF*, and *Sparse SDPA* formats. 

101 If CPLEX, Gurobi or MOSEK is installed, you can choose ``'cplex'``, 

102 ``'gurobi'``, or ``'mosek'``, respectively, to make use of that 

103 solver's export function and get access to more formats. 

104 

105 .. warning:: 

106 

107 For problems involving a symmetric matrix variable :math:`X` 

108 (typically, semidefinite programs), the expressions involving 

109 :math:`X` are stored in PICOS as a function of :math:`svec(X)`, the 

110 symmetric vectorized form of :math:`X` (see `Dattorro, ch.2.2.2.1 

111 <http://meboo.convexoptimization.com/Meboo.html>`_), 

112 and are also exported in that form. As a result, using an external 

113 solver on a problem description file exported by PICOS will also 

114 yield a solution in this symmetric vectorized form. 

115 

116 The CBF writer tries to write symmetric variables :math:`X` in the 

117 section ``PSDVAR`` of the .cbf file. However, this is possible only 

118 if the constraint :math:`X \succeq 0` appears in the problem, and no 

119 other LMI involves :math:`X`. If these two conditions are not 

120 satisfied, then the symmetric vectorization of :math:`X` is used as 

121 a (free) variable of the section ``VAR`` in the .cbf file, as 

122 explained in the previous paragraph. 

123 

124 .. warning:: 

125 

126 This function is severly outdated and may fail or not function as 

127 advertised. 

128 """ 

129 # HACK: This method abuses the internal problem representation of CVXOPT. 

130 # TODO: Add a proper method that transforms the problem into the canonical 

131 # form required here, and, if applicable, make also CVXOPT use it. 

132 

133 from .strategy import NoStrategyFound 

134 

135 # We first prepare the problem so it becomes exportable 

136 COMMERCIAL_WRITERS = ['cplex', 'mosek', 'gurobi'] 

137 if writer in COMMERCIAL_WRITERS: 

138 P = picos_problem.prepared(solver=writer) 

139 else: 

140 if picos_problem.continuous: 

141 # This should ensure that all quadratics have been cast as SOC 

142 # constraints. 

143 P = picos_problem.prepared(solver='cvxopt', assume_conic=True) 

144 else: 

145 try: 

146 P = picos_problem.prepared(solver=None) 

147 except NoStrategyFound as error: 

148 # FIXME: This will typically happen for integer linear programs 

149 # if no IP solver is available. 

150 raise RuntimeError("You try to export a problem for which no " 

151 "solver is available") from error 

152 

153 assert 'Linear' in P.type, "only LINEAR integer problems can be " \ 

154 "exported with the default writer" 

155 

156 # TODO: Need proper way to detect exponential cone constraints. This catches 

157 # only the Geometric Programs. 

158 if P.numberLSEConstraints: 

159 raise NotImplementedError("It is not possible (yet) to export GPs or " 

160 "problems with Exponential Cone Constraints") 

161 

162 # automatic extension recognition 

163 if not (any(filename.endswith(ext) for ext in 

164 (".mps", ".opf", ".cbf", ".ptf", ".lp", ".dat-s"))): 

165 if writer == "gurobi": 

166 if (P.numberConeConstraints + P.numberQuadConstraints) == 0: 

167 filename += ".lp" 

168 else: 

169 filename += ".mps" 

170 elif writer == "cplex": 

171 filename += ".lp" 

172 elif writer == "mosek": 

173 if (P.numberConeConstraints + P.numberQuadConstraints 

174 + P.numberSDPConstraints) == 0: 

175 filename += ".lp" 

176 elif P.numberQuadConstraints == 0: 

177 filename += ".cbf" 

178 else: 

179 filename += ".mps" 

180 elif writer == "picos": 

181 assert not P.numberQuadConstraints, \ 

182 "Expected no quadratic constraints after call to 'prepared'." 

183 if (P.numberConeConstraints + P.numberSDPConstraints) == 0: 

184 filename += ".lp" 

185 elif P.numberConeConstraints == 0: 

186 filename += ".dat-s" 

187 else: 

188 filename += ".cbf" 

189 else: 

190 raise Exception("unexpected writer") 

191 

192 if writer == "cplex": 

193 cpl = CPLEXSolver(P) 

194 cpl._load_problem() 

195 cpl.int.write(filename) 

196 

197 elif writer == "mosek": 

198 if filename.endswith(".cbf"): 

199 # This ensures that all quadratics are converted to SOC 

200 P = picos_problem.prepared(solver='cvxopt') 

201 msk = MOSEKSolver(P) 

202 msk._load_problem() 

203 msk.int.writedata(filename) 

204 

205 elif writer == "gurobi": 

206 grb = GurobiSolver(P) 

207 grb._load_problem() 

208 grb.int.write(filename) 

209 

210 elif writer == "picos": 

211 if filename[-3:] == ".lp": 

212 _write_lp(P, filename) 

213 elif filename[-6:] == ".dat-s": 

214 _write_sdpa(P, filename) 

215 elif filename[-4:] == ".cbf": 

216 _write_cbf(P, filename) 

217 else: 

218 raise Exception("unexpected file extension") 

219 else: 

220 raise Exception("unknown writer") 

221 

222 

223def _write_lp(P, filename): 

224 """Write the problem to a file in LP format.""" 

225 # add extension 

226 if filename[-3:] != ".lp": 

227 filename += ".lp" 

228 # check lp compatibility 

229 if ( 

230 P.numberConeConstraints 

231 + P.numberQuadConstraints 

232 + P.numberLSEConstraints 

233 + P.numberSDPConstraints 

234 ) > 0: 

235 raise Exception("the picos LP writer only accepts (MI)LP") 

236 # open file 

237 f = open(filename, "w") 

238 f.write("\\* file " + filename + " generated by picos*\\\n") 

239 

240 # HACK: This abuses the internal problem representation of CVXOPT. 

241 if P.continuous: 

242 localCvxoptInstance = CVXOPTSolver(P) 

243 else: 

244 localCvxoptInstance = CVXOPTSolver(P.continuous_relaxation()) 

245 localCvxoptInstance.import_variable_bounds = False 

246 localCvxoptInstance._load_problem() 

247 cvxoptVars = localCvxoptInstance.internal_problem() 

248 

249 # variable names 

250 varnames = {} 

251 i = 0 

252 for name, v in P.variables.items(): 

253 # full vectorization is used, name variables as x, x[j] or x[j,k] 

254 if isinstance(v._vec, FullVectorization): 

255 for ind in range(v.dim): 

256 if v.size == (1, 1): 

257 varnames[i] = name 

258 elif v.size[1] == 1: 

259 varnames[i] = name + "(" + str(ind) + ")" 

260 else: 

261 k, j = divmod(ind, v.size[0]) 

262 varnames[i] = name + "(" + str(j) + "," + str(k) + ")" 

263 varnames[i] = varnames[i].replace("[", "(") 

264 varnames[i] = varnames[i].replace("]", ")") 

265 i += 1 

266 else: 

267 for ind in range(v.dim): 

268 varnames[i] = 'vec_' + name + "(" + str(ind) + ")" 

269 i += 1 

270 

271 # affexpr writer 

272 def affexp_writer(constraint_name, indices, coefs): 

273 s = "" 

274 s += constraint_name 

275 s += " : " 

276 start = True 

277 for (i, v) in zip(indices, coefs): 

278 if v > 0 and not (start): 

279 s += "+ " 

280 s += "%.12g" % v 

281 s += " " 

282 s += varnames[i] 

283 # not the first term anymore 

284 start = False 

285 if not (coefs): 

286 s += "0.0 " 

287 s += varnames[0] 

288 return s 

289 

290 print("writing problem in " + filename + "...") 

291 

292 # objective 

293 if P.objective.direction == "max": 

294 f.write("Maximize\n") 

295 # max handled directly 

296 cvxoptVars["c"] = -cvxoptVars["c"] 

297 else: 

298 f.write("Minimize\n") 

299 I = cvxopt.sparse(cvxoptVars["c"]).I 

300 V = cvxopt.sparse(cvxoptVars["c"]).V 

301 

302 f.write(affexp_writer("obj", I, V)) 

303 f.write("\n") 

304 

305 f.write("Subject To\n") 

306 bounds = {} 

307 # equality constraints: 

308 Ai, Aj, Av = (cvxoptVars["A"].I, cvxoptVars["A"].J, cvxoptVars["A"].V) 

309 ijvs = sorted(zip(Ai, Aj, Av)) 

310 del Ai, Aj, Av 

311 itojv = {} 

312 lasti = -1 

313 for (i, j, v) in ijvs: 

314 if i == lasti: 

315 itojv[i].append((j, v)) 

316 else: 

317 lasti = i 

318 itojv[i] = [(j, v)] 

319 ieq = 0 

320 for i, jv in itojv.items(): 

321 J = [jvk[0] for jvk in jv] 

322 V = [jvk[1] for jvk in jv] 

323 if len(J) == 1: 

324 # fixed variable 

325 b = cvxoptVars["b"][i] / V[0] 

326 bounds[J[0]] = (b, b) 

327 else: 

328 # affine equality 

329 b = cvxoptVars["b"][i] 

330 f.write(affexp_writer("eq" + str(ieq), J, V)) 

331 f.write(" = ") 

332 f.write("%.12g" % b) 

333 f.write("\n") 

334 ieq += 1 

335 

336 # inequality constraints: 

337 Gli, Glj, Glv = (cvxoptVars["Gl"].I, cvxoptVars["Gl"].J, cvxoptVars["Gl"].V) 

338 ijvs = sorted(zip(Gli, Glj, Glv)) 

339 del Gli, Glj, Glv 

340 itojv = {} 

341 lasti = -1 

342 for (i, j, v) in ijvs: 

343 if i == lasti: 

344 itojv[i].append((j, v)) 

345 else: 

346 lasti = i 

347 itojv[i] = [(j, v)] 

348 iaff = 0 

349 for i, jv in itojv.items(): 

350 J = [jvk[0] for jvk in jv] 

351 V = [jvk[1] for jvk in jv] 

352 b = cvxoptVars["hl"][i] 

353 f.write(affexp_writer("in" + str(iaff), J, V)) 

354 f.write(" <= ") 

355 f.write("%.12g" % b) 

356 f.write("\n") 

357 iaff += 1 

358 

359 # variable bounds 

360 # retrieve as a dictionary {index -> (lo,up)} 

361 i_var = 0 

362 for var in P.variables.values(): 

363 for ind, lo in var.bound_dicts[0].items(): 

364 (current_lo, current_up) = bounds.get( 

365 i_var + ind, (-INFINITY, INFINITY)) 

366 new_lo = max(lo, current_lo) 

367 bounds[i_var + ind] = (new_lo, current_up) 

368 for ind, up in var.bound_dicts[1].items(): 

369 (current_lo, current_up) = bounds.get( 

370 i_var + ind, (-INFINITY, INFINITY)) 

371 new_up = min(up, current_up) 

372 bounds[i_var + ind] = (current_lo, new_up) 

373 i_var += var.dim 

374 

375 f.write("Bounds\n") 

376 for i in range(P.numberOfVars): 

377 if i in bounds: 

378 bl, bu = bounds[i] 

379 else: 

380 bl, bu = -INFINITY, INFINITY 

381 if bl == -INFINITY and bu == INFINITY: 

382 f.write(varnames[i] + " free") 

383 elif bl == bu: 

384 f.write(varnames[i] + (" = %.12g" % bl)) 

385 elif bl < bu: 

386 if bl == -INFINITY: 

387 f.write("-inf <= ") 

388 else: 

389 f.write("%.12g" % bl) 

390 f.write(" <= ") 

391 f.write(varnames[i]) 

392 if bu == INFINITY: 

393 f.write("<= +inf") 

394 else: 

395 f.write(" <= ") 

396 f.write("%.12g" % bu) 

397 f.write("\n") 

398 

399 # general integers 

400 f.write("Generals\n") 

401 i_var = 0 

402 for name, v in P.variables.items(): 

403 if isinstance(v, IntegerVariable): 

404 for ind in range(v.dim): 

405 f.write(varnames[i_var + ind] + "\n") 

406 i_var += v.dim 

407 

408 # binary variables 

409 f.write("Binaries\n") 

410 i_var = 0 

411 for name, v in P.variables.items(): 

412 if isinstance(v, BinaryVariable): 

413 for ind in range(v.dim): 

414 f.write(varnames[i_var + ind] + "\n") 

415 i_var += v.dim 

416 

417 f.write("End\n") 

418 print("done.") 

419 f.close() 

420 

421 

422def _write_sdpa(P, filename): 

423 """Write the problem to a file in Sparse SDPA format.""" 

424 # HACK: This abuses the internal problem representation of CVXOPT. 

425 localCvxoptInstance = CVXOPTSolver(P) 

426 localCvxoptInstance._load_problem() 

427 cvxoptVars = localCvxoptInstance.internal_problem() 

428 

429 dims = {} 

430 dims["s"] = [int(numpy.sqrt(Gsi.size[0])) for Gsi in cvxoptVars["Gs"]] 

431 dims["l"] = cvxoptVars["Gl"].size[0] 

432 dims["q"] = [Gqi.size[0] for Gqi in cvxoptVars["Gq"]] 

433 G = cvxoptVars["Gl"] 

434 h = cvxoptVars["hl"] 

435 

436 # handle the equalities as 2 ineq 

437 if cvxoptVars["A"].size[0] > 0: 

438 G = cvxopt.sparse([G, cvxoptVars["A"]]) 

439 G = cvxopt.sparse([G, -cvxoptVars["A"]]) 

440 h = cvxopt.matrix([h, cvxoptVars["b"]]) 

441 h = cvxopt.matrix([h, -cvxoptVars["b"]]) 

442 dims["l"] += 2 * cvxoptVars["A"].size[0] 

443 

444 for i in range(len(dims["q"])): 

445 G = cvxopt.sparse([G, cvxoptVars["Gq"][i]]) 

446 h = cvxopt.matrix([h, cvxoptVars["hq"][i]]) 

447 

448 for i in range(len(dims["s"])): 

449 G = cvxopt.sparse([G, cvxoptVars["Gs"][i]]) 

450 h = cvxopt.matrix([h, cvxoptVars["hs"][i]]) 

451 

452 # Remove the lines in A and b corresponding to 0==0 

453 JP = list(set(cvxoptVars["A"].I)) 

454 IP = range(len(JP)) 

455 VP = [1] * len(JP) 

456 

457 # is there a constraint of the form 0==a(a not 0) ? 

458 if any([b for (i, b) in enumerate(cvxoptVars["b"]) if i not in JP]): 

459 raise Exception("infeasible constraint of the form 0=a") 

460 

461 from cvxopt import sparse, spmatrix 

462 

463 PP = spmatrix(VP, IP, JP, (len(IP), cvxoptVars["A"].size[0])) 

464 cvxoptVars["A"] = PP * cvxoptVars["A"] 

465 cvxoptVars["b"] = PP * cvxoptVars["b"] 

466 c = cvxoptVars["c"] 

467 # ------------------------------------------------------------# 

468 # make A,B,and blockstruct. # 

469 # This code is a modification of the conelp function in SMCP # 

470 # ------------------------------------------------------------# 

471 Nl = dims["l"] 

472 Nq = dims["q"] 

473 Ns = dims["s"] 

474 if not Nl: 

475 Nl = 0 

476 

477 P_m = G.size[1] 

478 

479 P_b = -c 

480 P_blockstruct = [] 

481 if Nl: 

482 P_blockstruct.append(-Nl) 

483 for i in Nq: 

484 P_blockstruct.append(i) 

485 for i in Ns: 

486 P_blockstruct.append(i) 

487 

488 # write data 

489 # add extension 

490 if filename[-6:] != ".dat-s": 

491 filename += ".dat-s" 

492 

493 # open file 

494 f = open(filename, "w") 

495 f.write('"file ' + filename + ' generated by picos"\n') 

496 if P.options.verbosity >= 1: 

497 print("writing problem in " + filename + "...") 

498 f.write(str(P.numberOfVars) + " = number of vars\n") 

499 f.write(str(len(P_blockstruct)) + " = number of blocs\n") 

500 # bloc structure 

501 f.write(str(P_blockstruct).replace("[", "(").replace("]", ")")) 

502 f.write(" = BlocStructure\n") 

503 # c vector (objective) 

504 f.write(str(list(-P_b)).replace("[", "{").replace("]", "}")) 

505 f.write("\n") 

506 # coefs 

507 for k in range(P_m + 1): 

508 if k != 0: 

509 v = sparse(G[:, k - 1]) 

510 else: 

511 v = +sparse(h) 

512 

513 ptr = 0 

514 block = 0 

515 # lin. constraints 

516 if Nl: 

517 u = v[:Nl] 

518 for i, j, value in zip(u.I, u.I, u.V): 

519 f.write( 

520 "{0}\t{1}\t{2}\t{3}\t{4}\n".format( 

521 k, block + 1, j + 1, i + 1, -value 

522 ) 

523 ) 

524 ptr += Nl 

525 block += 1 

526 

527 # SOC constraints 

528 for nq in Nq: 

529 u0 = v[ptr] 

530 u1 = v[ptr + 1:ptr + nq] 

531 tmp = spmatrix( 

532 u1.V, [nq - 1 for j in range(len(u1))], u1.I, (nq, nq) 

533 ) 

534 if not u0 == 0.0: 

535 tmp += spmatrix(u0, range(nq), range(nq), (nq, nq)) 

536 for i, j, value in zip(tmp.I, tmp.J, tmp.V): 

537 f.write( 

538 "{0}\t{1}\t{2}\t{3}\t{4}\n".format( 

539 k, block + 1, j + 1, i + 1, -value 

540 ) 

541 ) 

542 ptr += nq 

543 block += 1 

544 

545 # SDP constraints 

546 for ns in Ns: 

547 u = v[ptr:ptr + ns ** 2] 

548 for index_k, index in enumerate(u.I): 

549 j, i = divmod(index, ns) 

550 if j <= i: 

551 f.write( 

552 "{0}\t{1}\t{2}\t{3}\t{4}\n".format( 

553 k, block + 1, j + 1, i + 1, -u.V[index_k] 

554 ) 

555 ) 

556 ptr += ns ** 2 

557 block += 1 

558 

559 f.close() 

560 

561 

562# This is an ugly function from a former version of picos 

563# It separates a linear constraint between 'plain' vars and matrix 'bar' 

564# variables J and V denote the sparse indices/values of the constraints for the 

565# whole (s-)vectorized vector (without offset) 

566def _separate_linear_cons_plain_bar_vars(J, V, idx_sdp_vars): 

567 # sparse values of the constraint for 'plain' variables 

568 jj = [] 

569 vv = [] 

570 # sparse values of the constraint for the next svec bar variable 

571 js = [] 

572 vs = [] 

573 mats = [] 

574 offset = 0 

575 if idx_sdp_vars: 

576 idxsdpvars = [ti for ti in idx_sdp_vars] 

577 nextsdp = idxsdpvars.pop() 

578 else: 

579 return J, V, [] 

580 for (j, v) in zip(J, V): 

581 if j < nextsdp[0]: 

582 jj.append(j - offset) 

583 vv.append(v) 

584 elif j < nextsdp[1]: 

585 js.append(j - nextsdp[0]) 

586 vs.append(v) 

587 else: 

588 while j >= nextsdp[1]: 

589 mats.append( 

590 devectorize( 

591 cvxopt.spmatrix( 

592 vs, 

593 js, 

594 [0] * len(js), 

595 (nextsdp[1] - nextsdp[0], 1) 

596 )).T) 

597 js = [] 

598 vs = [] 

599 offset += (nextsdp[1] - nextsdp[0]) 

600 try: 

601 nextsdp = idxsdpvars.pop() 

602 except IndexError: 

603 nextsdp = (float('inf'), float('inf')) 

604 if j < nextsdp[0]: 

605 jj.append(j - offset) 

606 vv.append(v) 

607 elif j < nextsdp[1]: 

608 js.append(j - nextsdp[0]) 

609 vs.append(v) 

610 while len(mats) < len(idx_sdp_vars): 

611 mats.append( 

612 devectorize( 

613 cvxopt.spmatrix( 

614 vs, 

615 js, 

616 [0] * len(js), 

617 (nextsdp[1] - nextsdp[0], 1) 

618 )).T) 

619 js = [] 

620 vs = [] 

621 nextsdp = (0, 1) # doesnt matter, it will be an empt matrix anyway 

622 return jj, vv, mats 

623 

624 

625def devectorize(vec): 

626 """Create a matrix from a symmetric vectorization.""" 

627 from ..expressions.vectorizations import SymmetricVectorization 

628 v = vec.size[0] 

629 n = int((1 + 8 * v) ** 0.5 - 1) // 2 

630 return SymmetricVectorization((n, n)).devectorize(vec) 

631 

632 

633def _write_cbf(P, filename, uptri=False): 

634 """Write the problem to a file in Sparse SDPA format. 

635 

636 :param bool uptri: Whether upper triangular elements of symmetric 

637 matrices are specified. 

638 """ 

639 # write data 

640 # add extension 

641 if filename[-4:] != ".cbf": 

642 filename += ".cbf" 

643 

644 # parse variables 

645 semidef_vars = set() 

646 for cons in P.constraints: 

647 cs = P.constraints[cons] 

648 if isinstance(cs, LMIConstraint): 

649 if cs.semidefVar: 

650 semidef_vars.add(cs.semidefVar) 

651 

652 NUMVAR_SCALAR = int( 

653 sum( 

654 [var.dim 

655 for var in P.variables.values() 

656 if var not in semidef_vars] 

657 ) 

658 ) 

659 

660 ind = 0 

661 indices = [] 

662 start_indices = {} 

663 for v in P.variables.values(): 

664 indices.append((ind, ind + v.dim, v)) 

665 start_indices[v] = ind 

666 ind += v.dim 

667 

668 indices = sorted(indices) 

669 idxsdpvars = [ 

670 (si, ei) for (si, ei, v) in indices[::-1] if v in semidef_vars] 

671 # search if some semidef vars are implied in other semidef constraints 

672 PSD_not_handled = [] 

673 for c in P.constraints.values(): 

674 if isinstance(c, LMIConstraint) and c not in semidef_vars: 

675 for v in (c.lhs - c.rhs)._coefs: 

676 if v in semidef_vars: 

677 start_index = start_indices[v] 

678 idx = (start_index, start_index + v.dim) 

679 if idx in idxsdpvars: 

680 PSD_not_handled.append(v) 

681 NUMVAR_SCALAR += idx[1] - idx[0] 

682 idxsdpvars.remove(idx) 

683 

684 barvars = bool(idxsdpvars) 

685 

686 # find integer variables, retrieve bounds and put 0-1 bounds on binaries 

687 ints = [] 

688 ind = 0 

689 varbounds_lo = {} 

690 varbounds_up = {} 

691 

692 for k, var in P.variables.items(): 

693 if isinstance(var, BinaryVariable): 

694 for relind, absind in enumerate(range(ind, ind + var.dim)): 

695 ints.append(absind) 

696 clb = var.bound_dicts[0].get(relind, -INFINITY) 

697 cub = var.bound_dicts[1].get(relind, INFINITY) 

698 varbounds_lo[absind] = max(0.0, clb) 

699 varbounds_up[absind] = min(1.0, cub) 

700 

701 elif (isinstance(var, IntegerVariable) or 

702 isinstance(var, CONTINUOUS_VARTYPES)): 

703 for relind, absind in enumerate(range(ind, ind + var.dim)): 

704 if isinstance(var, IntegerVariable): 

705 ints.append(absind) 

706 clb = var.bound_dicts[0].get(relind, -INFINITY) 

707 cub = var.bound_dicts[1].get(relind, INFINITY) 

708 varbounds_lo[absind] = clb 

709 varbounds_up[absind] = cub 

710 

711 else: 

712 raise Exception("variable type not handled by _write_cbf()") 

713 ind += var.dim 

714 

715 if barvars: 

716 ints, _, mats = _separate_linear_cons_plain_bar_vars( 

717 ints, [0.0] * len(ints), idxsdpvars 

718 ) 

719 if any([bool(mat) for mat in mats]): 

720 raise Exception( 

721 "semidef vars with integer elements are not supported" 

722 ) 

723 

724 # open file 

725 f = open(filename, "w") 

726 f.write("#file " + filename + " generated by picos\n") 

727 print("writing problem in " + filename + "...") 

728 

729 f.write("VER\n") 

730 f.write("1\n\n") 

731 

732 f.write("OBJSENSE\n") 

733 if P.objective.direction == "max": 

734 f.write("MAX\n\n") 

735 else: 

736 f.write("MIN\n\n") 

737 

738 # VARIABLEs 

739 

740 if barvars: 

741 f.write("PSDVAR\n") 

742 f.write(str(len(idxsdpvars)) + "\n") 

743 for si, ei in idxsdpvars: 

744 ni = int(((8 * (ei - si) + 1) ** 0.5 - 1) / 2.0) 

745 f.write(str(ni) + "\n") 

746 f.write("\n") 

747 

748 # bounds 

749 cones = [] 

750 conecons = [] 

751 Acoord = [] 

752 Bcoord = [] 

753 iaff = 0 

754 offset = 0 

755 for si, ei, v in indices: 

756 if v in semidef_vars and not (v in PSD_not_handled): 

757 offset += ei - si 

758 else: 

759 if all(varbounds_lo[ind] == 0 for ind in range(si, ei)): 

760 cones.append(("L+", ei - si)) 

761 elif all(varbounds_up[ind] == 0 for ind in range(si, ei)): 

762 cones.append(("L-", ei - si)) 

763 else: 

764 cones.append(("F", ei - si)) 

765 if any(varbounds_lo[ind] != 0 for ind in range(si, ei)): 

766 for ind in range(si, ei): 

767 l = varbounds_lo[ind] 

768 if l - 1 > -INFINITY: