Coverage for picos/modeling/file_in.py: 4.14%

362 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-03-26 07:46 +0000

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 reading optimization problems from a file.""" 

21 

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

23# TODO: This file is using ad-hoc parser implementations. Use a proper parser 

24# generation library and a unified 'read' function for different formats. 

25# ------------------------------------------------------------------------------ 

26 

27import warnings 

28 

29import cvxopt 

30import numpy 

31 

32from ..apidoc import api_end, api_start 

33from ..expressions import AffineExpression, Constant 

34from .problem import Problem 

35 

36_API_START = api_start(globals()) 

37# ------------------------------- 

38 

39 

40def _spmatrix(*args, **kwargs): 

41 """Create a CVXOPT sparse matrix. 

42 

43 A wrapper around :func:`cvxopt.spmatrix` that converts indices to 

44 :class:`int`, if necessary. 

45 

46 Works around PICOS sometimes passing indices as :class:`numpy:numpy.int64`. 

47 """ 

48 try: 

49 return cvxopt.spmatrix(*args, **kwargs) 

50 except TypeError as error: 

51 # CVXOPT does not like NumPy's int64 scalar type for indices, so attempt 

52 # to convert all indices to Python's int. 

53 if str(error) == "non-numeric type in list": 

54 newargs = list(args) 

55 

56 for argNum, arg in enumerate(args): 

57 if argNum in (1, 2): # Positional I, J. 

58 newargs[argNum] = [int(x) for x in args[argNum]] 

59 

60 for kw in "IJ": 

61 if kw in kwargs: 

62 kwargs[kw] = [int(x) for x in kwargs[kw]] 

63 

64 return cvxopt.spmatrix(*newargs, **kwargs) 

65 else: 

66 raise 

67 

68 

69def _break_cols(mat, sizes): 

70 """Help with the reading of CBF files.""" 

71 n = len(sizes) 

72 I, J, V = [], [], [] 

73 for i in range(n): 

74 I.append([]) 

75 J.append([]) 

76 V.append([]) 

77 cumsz = numpy.cumsum(sizes) 

78 import bisect 

79 

80 for i, j, v in zip(mat.I, mat.J, mat.V): 

81 block = bisect.bisect(cumsz, j) 

82 I[block].append(i) 

83 V[block].append(v) 

84 if block == 0: 

85 J[block].append(j) 

86 else: 

87 J[block].append(j - cumsz[block - 1]) 

88 return [ 

89 _spmatrix(V[k], I[k], J[k], (mat.size[0], sz)) 

90 for k, sz in enumerate(sizes) 

91 ] 

92 

93 

94def _break_rows(mat, sizes): 

95 """Help with the reading of CBF files.""" 

96 n = len(sizes) 

97 I, J, V = [], [], [] 

98 for i in range(n): 

99 I.append([]) 

100 J.append([]) 

101 V.append([]) 

102 cumsz = numpy.cumsum(sizes) 

103 import bisect 

104 

105 for i, j, v in zip(mat.I, mat.J, mat.V): 

106 block = bisect.bisect(cumsz, i) 

107 J[block].append(j) 

108 V[block].append(v) 

109 if block == 0: 

110 I[block].append(i) 

111 else: 

112 I[block].append(i - cumsz[block - 1]) 

113 return [ 

114 _spmatrix(V[k], I[k], J[k], (sz, mat.size[1])) 

115 for k, sz in enumerate(sizes) 

116 ] 

117 

118 

119def _block_idx(i, sizes): 

120 """Help with the reading of CBF files.""" 

121 # if there are blocks of sizes n1,...,nk and i is 

122 # the index of an element of the big vectorized variable, 

123 # returns the block of i and its index inside the sub-block. 

124 cumsz = numpy.cumsum(sizes) 

125 import bisect 

126 

127 block = bisect.bisect(cumsz, i) 

128 return block, (i if block == 0 else i - cumsz[block - 1]) 

129 

130 

131def _read_cbf_block(P, blocname, f, parsed_blocks): 

132 """Help with the reading of CBF files.""" 

133 if blocname == "OBJSENSE": 

134 objsense = f.readline().split()[0].lower() 

135 P.set_objective(objsense, P.objective.normalized.function) 

136 return None 

137 elif blocname == "PSDVAR": 

138 n = int(f.readline()) 

139 vardims = [] 

140 XX = [] 

141 for i in range(n): 

142 ni = int(f.readline()) 

143 vardims.append(ni) 

144 Xi = P.add_variable("X[" + str(i) + "]", (ni, ni), "symmetric") 

145 XX.append(Xi) 

146 P.add_constraint(Xi >> 0) 

147 return vardims, XX 

148 elif blocname == "VAR": 

149 Nscalar, ncones = [int(fi) for fi in f.readline().split()] 

150 tot_dim = 0 

151 var_structure = [] 

152 xx = [] 

153 for i in range(ncones): 

154 lsplit = f.readline().split() 

155 tp, dim = lsplit[0], int(lsplit[1]) 

156 tot_dim += dim 

157 var_structure.append(dim) 

158 if tp == "F": 

159 xi = P.add_variable("x[" + str(i) + "]", dim) 

160 elif tp == "L+": 

161 xi = P.add_variable("x[" + str(i) + "]", dim, lower=0) 

162 elif tp == "L-": 

163 xi = P.add_variable("x[" + str(i) + "]", dim, upper=0) 

164 elif tp == "L=": 

165 xi = P.add_variable("x[" + str(i) + "]", dim, lower=0, upper=0) 

166 elif tp == "Q": 

167 xi = P.add_variable("x[" + str(i) + "]", dim) 

168 P.add_constraint(abs(xi[1:]) < xi[0]) 

169 elif tp == "QR": 

170 xi = P.add_variable("x[" + str(i) + "]", dim) 

171 P.add_constraint(abs(xi[2:]) ** 2 < 2 * xi[0] * xi[1]) 

172 xx.append(xi) 

173 if tot_dim != Nscalar: 

174 raise Exception("VAR dimensions do not match the header") 

175 return Nscalar, var_structure, xx 

176 elif blocname == "INT": 

177 n = int(f.readline()) 

178 ints = {} 

179 for k in range(n): 

180 j = int(f.readline()) 

181 i, col = _block_idx(j, parsed_blocks["VAR"][1]) 

182 ints.setdefault(i, []) 

183 ints[i].append(col) 

184 x = parsed_blocks["VAR"][2] 

185 for i in ints: 

186 if len(ints[i]) == x[i].size[0]: 

187 x[i].vtype = "integer" 

188 else: 

189 x.append( 

190 P.add_variable( 

191 "x_int[" + str(i) + "]", len(ints[i]), "integer" 

192 ) 

193 ) 

194 for k, j in enumerate(ints[i]): 

195 P.add_constraint(x[i][j] == x[-1][k]) 

196 return x 

197 elif blocname == "CON": 

198 Ncons, ncones = [int(fi) for fi in f.readline().split()] 

199 cons_structure = [] 

200 tot_dim = 0 

201 for i in range(ncones): 

202 lsplit = f.readline().split() 

203 tp, dim = lsplit[0], int(lsplit[1]) 

204 tot_dim += dim 

205 cons_structure.append((tp, dim)) 

206 if tot_dim != Ncons: 

207 raise Exception("CON dimensions do not match the header") 

208 return Ncons, cons_structure 

209 elif blocname == "PSDCON": 

210 n = int(f.readline()) 

211 psdcons_structure = [] 

212 for i in range(n): 

213 ni = int(f.readline()) 

214 psdcons_structure.append(ni) 

215 return psdcons_structure 

216 elif blocname == "OBJACOORD": 

217 n = int(f.readline()) 

218 J = [] 

219 V = [] 

220 for i in range(n): 

221 lsplit = f.readline().split() 

222 j, v = int(lsplit[0]), float(lsplit[1]) 

223 J.append(j) 

224 V.append(v) 

225 return _spmatrix(V, [0] * len(J), J, (1, parsed_blocks["VAR"][0])) 

226 elif blocname == "OBJBCOORD": 

227 return float(f.readline()) 

228 elif blocname == "OBJFCOORD": 

229 n = int(f.readline()) 

230 Fobj = [ 

231 _spmatrix([], [], [], (ni, ni)) for ni in parsed_blocks["PSDVAR"][0] 

232 ] 

233 for k in range(n): 

234 lsplit = f.readline().split() 

235 j, row, col, v = ( 

236 int(lsplit[0]), 

237 int(lsplit[1]), 

238 int(lsplit[2]), 

239 float(lsplit[3]), 

240 ) 

241 Fobj[j][row, col] = v 

242 if row != col: 

243 Fobj[j][col, row] = v 

244 return Fobj 

245 elif blocname == "FCOORD": 

246 n = int(f.readline()) 

247 Fblocks = {} 

248 for k in range(n): 

249 lsplit = f.readline().split() 

250 i, j, row, col, v = ( 

251 int(lsplit[0]), 

252 int(lsplit[1]), 

253 int(lsplit[2]), 

254 int(lsplit[3]), 

255 float(lsplit[4]), 

256 ) 

257 if i not in Fblocks: 

258 Fblocks[i] = [ 

259 _spmatrix([], [], [], (ni, ni)) 

260 for ni in parsed_blocks["PSDVAR"][0] 

261 ] 

262 Fblocks[i][j][row, col] = v 

263 if row != col: 

264 Fblocks[i][j][col, row] = v 

265 return Fblocks 

266 elif blocname == "ACOORD": 

267 n = int(f.readline()) 

268 J = [] 

269 V = [] 

270 I = [] 

271 for k in range(n): 

272 lsplit = f.readline().split() 

273 i, j, v = int(lsplit[0]), int(lsplit[1]), float(lsplit[2]) 

274 I.append(i) 

275 J.append(j) 

276 V.append(v) 

277 return _spmatrix( 

278 V, I, J, (parsed_blocks["CON"][0], parsed_blocks["VAR"][0]) 

279 ) 

280 elif blocname == "BCOORD": 

281 n = int(f.readline()) 

282 V = [] 

283 I = [] 

284 for k in range(n): 

285 lsplit = f.readline().split() 

286 i, v = int(lsplit[0]), float(lsplit[1]) 

287 I.append(i) 

288 V.append(v) 

289 return _spmatrix(V, I, [0] * len(I), (parsed_blocks["CON"][0], 1)) 

290 elif blocname == "HCOORD": 

291 n = int(f.readline()) 

292 Hblocks = {} 

293 for k in range(n): 

294 lsplit = f.readline().split() 

295 i, j, row, col, v = ( 

296 int(lsplit[0]), 

297 int(lsplit[1]), 

298 int(lsplit[2]), 

299 int(lsplit[3]), 

300 float(lsplit[4]), 

301 ) 

302 if j not in Hblocks: 

303 Hblocks[j] = [ 

304 _spmatrix([], [], [], (ni, ni)) 

305 for ni in parsed_blocks["PSDCON"] 

306 ] 

307 Hblocks[j][i][row, col] = v 

308 if row != col: 

309 Hblocks[j][i][col, row] = v 

310 return Hblocks 

311 elif blocname == "DCOORD": 

312 n = int(f.readline()) 

313 Dblocks = [ 

314 _spmatrix([], [], [], (ni, ni)) for ni in parsed_blocks["PSDCON"] 

315 ] 

316 for k in range(n): 

317 lsplit = f.readline().split() 

318 i, row, col, v = ( 

319 int(lsplit[0]), 

320 int(lsplit[1]), 

321 int(lsplit[2]), 

322 float(lsplit[3]), 

323 ) 

324 Dblocks[i][row, col] = v 

325 if row != col: 

326 Dblocks[i][col, row] = v 

327 return Dblocks 

328 else: 

329 raise Exception("unexpected block name") 

330 

331 

332def import_cbf(filename): 

333 """Create a :class:`~picos.Problem` from a CBF file. 

334 

335 The created problem contains one (multidimensional) variable for each cone 

336 specified in the section ``VAR`` of the .cbf file, and one 

337 (multidimmensional) constraint for each cone specified in the sections 

338 ``CON`` and ``PSDCON``. 

339 

340 :returns: A tuple ``(P, x, X, params)`` where 

341 

342 - ``P`` is the imported picos :class:`~picos.Problem` object, 

343 - ``x`` is a list of multidimensional variables representing the scalar 

344 variables found in the file, 

345 - ``X`` is a list of symmetric variables representing the positive 

346 semidefinite variables found in the file, and 

347 - ``params`` is a dictionary containing PICOS cosntants used to define 

348 the problem. Indexing is with respect to the blocks of variables as 

349 defined in the sections ``VAR`` and ``CON`` of the file. 

350 """ 

351 P = Problem() 

352 

353 try: 

354 f = open(filename, "r") 

355 except IOError: 

356 filename += ".cbf" 

357 f = open(filename, "r") 

358 

359 line = f.readline() 

360 while not line.startswith("VER"): 

361 line = f.readline() 

362 

363 ver = int(f.readline()) 

364 if ver != 1: 

365 warnings.warn("CBF file has a version other than 1.") 

366 

367 structure_keywords = ["OBJSENSE", "PSDVAR", "VAR", "INT", "PSDCON", "CON"] 

368 data_keywords = [ 

369 "OBJFCOORD", 

370 "OBJACOORD", 

371 "OBJBCOORD", 

372 "FCOORD", 

373 "ACOORD", 

374 "BCOORD", 

375 "HCOORD", 

376 "DCOORD", 

377 ] 

378 

379 structure_mode = True # still parsing structure blocks 

380 seen_blocks = [] 

381 parsed_blocks = {} 

382 while True: 

383 line = f.readline() 

384 if not line: 

385 break 

386 lsplit = line.split() 

387 if lsplit and lsplit[0] in structure_keywords: 

388 if lsplit[0] == "INT" and ("VAR" not in seen_blocks): 

389 raise Exception("INT BLOCK before VAR BLOCK") 

390 if lsplit[0] == "CON" and not ( 

391 "VAR" in seen_blocks or "PSDVAR" in seen_blocks 

392 ): 

393 raise Exception("CON BLOCK before VAR/PSDVAR BLOCK") 

394 if lsplit[0] == "PSDCON" and not ( 

395 "VAR" in seen_blocks or "PSDVAR" in seen_blocks 

396 ): 

397 raise Exception("PSDCON BLOCK before VAR/PSDVAR BLOCK") 

398 if lsplit[0] == "VAR" and ( 

399 "CON" in seen_blocks or "PSDCON" in seen_blocks 

400 ): 

401 raise Exception("VAR BLOCK after CON/PSDCON BLOCK") 

402 if lsplit[0] == "PSDVAR" and ( 

403 "CON" in seen_blocks or "PSDCON" in seen_blocks 

404 ): 

405 raise Exception("PSDVAR BLOCK after CON/PSDCON BLOCK") 

406 if structure_mode: 

407 parsed_blocks[lsplit[0]] = _read_cbf_block( 

408 P, lsplit[0], f, parsed_blocks 

409 ) 

410 seen_blocks.append(lsplit[0]) 

411 else: 

412 raise Exception("Structure keyword after first data item") 

413 if lsplit and lsplit[0] in data_keywords: 

414 if "OBJSENSE" not in seen_blocks: 

415 raise Exception("missing OBJSENSE block") 

416 if not ("VAR" in seen_blocks or "PSDVAR" in seen_blocks): 

417 raise Exception("missing VAR/PSDVAR block") 

418 if lsplit[0] in ("OBJFCOORD", "FCOORD") and not ( 

419 "PSDVAR" in seen_blocks 

420 ): 

421 raise Exception("missing PSDVAR block") 

422 if lsplit[0] in ("OBJACOORD", "ACOORD", "HCOORD") and not ( 

423 "VAR" in seen_blocks 

424 ): 

425 raise Exception("missing VAR block") 

426 if lsplit[0] in ("DCOORD", "HCOORD") and not ( 

427 "PSDCON" in seen_blocks 

428 ): 

429 raise Exception("missing PSDCON block") 

430 structure_mode = False 

431 parsed_blocks[lsplit[0]] = _read_cbf_block( 

432 P, lsplit[0], f, parsed_blocks 

433 ) 

434 seen_blocks.append(lsplit[0]) 

435 

436 f.close() 

437 # variables 

438 if "VAR" in parsed_blocks: 

439 Nvars, varsz, x = parsed_blocks["VAR"] 

440 else: 

441 x = None 

442 

443 if "INT" in parsed_blocks: 

444 x = parsed_blocks["INT"] 

445 

446 if "PSDVAR" in parsed_blocks: 

447 psdsz, X = parsed_blocks["PSDVAR"] 

448 else: 

449 X = None 

450 

451 # objective 

452 obj_constant = parsed_blocks.get("OBJBCOORD", 0) 

453 bobj = Constant("bobj", obj_constant) 

454 obj = Constant("bobj", obj_constant) 

455 

456 aobj = {} 

457 if "OBJACOORD" in parsed_blocks: 

458 obj_vecs = _break_cols(parsed_blocks["OBJACOORD"], varsz) 

459 aobj = {} 

460 for k, v in enumerate(obj_vecs): 

461 if v: 

462 aobj[k] = Constant("c[" + str(k) + "]", v) 

463 obj += aobj[k] * x[k] 

464 

465 Fobj = {} 

466 if "OBJFCOORD" in parsed_blocks: 

467 Fbl = parsed_blocks["OBJFCOORD"] 

468 for i, Fi in enumerate(Fbl): 

469 if Fi: 

470 Fobj[i] = Constant("F[" + str(i) + "]", Fi) 

471 obj += Fobj[i] | X[i] 

472 

473 P.set_objective(P.objective.normalized.direction, obj) 

474 

475 # cone constraints 

476 bb = {} 

477 AA = {} 

478 FF = {} 

479 if "CON" in parsed_blocks: 

480 Ncons, structcons = parsed_blocks["CON"] 

481 szcons = [s for tp, s in structcons] 

482 

483 b = parsed_blocks.get("BCOORD", _spmatrix([], [], [], (Ncons, 1))) 

484 bvecs = _break_rows(b, szcons) 

485 consexp = [] 

486 for i, bi in enumerate(bvecs): 

487 bb[i] = Constant("b[" + str(i) + "]", bi) 

488 consexp.append(Constant("b[" + str(i) + "]", bi)) 

489 

490 A = parsed_blocks.get("ACOORD", _spmatrix([], [], [], (Ncons, Nvars))) 

491 Ablc = _break_rows(A, szcons) 

492 for i, Ai in enumerate(Ablc): 

493 Aiblocs = _break_cols(Ai, varsz) 

494 for j, Aij in enumerate(Aiblocs): 

495 if Aij: 

496 AA[i, j] = Constant("A[" + str((i, j)) + "]", Aij) 

497 consexp[i] += AA[i, j] * x[j] 

498 

499 Fcoords = parsed_blocks.get("FCOORD", {}) 

500 for k, mats in Fcoords.items(): 

501 i, row = _block_idx(k, szcons) 

502 row_exp = AffineExpression.zero() 

503 for j, mat in enumerate(mats): 

504 if mat: 

505 FF[i, j, row] = Constant("F[" + str((i, j, row)) + "]", mat) 

506 row_exp += FF[i, j, row] | X[j] 

507 

508 consexp[i] += ( 

509 "e_" + str(row) + "(" + str(szcons[i]) + ",1)" 

510 ) * row_exp 

511 

512 for i, (tp, sz) in enumerate(structcons): 

513 if tp == "F": 

514 continue 

515 elif tp == "L-": 

516 P.add_constraint(consexp[i] < 0) 

517 elif tp == "L+": 

518 P.add_constraint(consexp[i] > 0) 

519 elif tp == "L=": 

520 P.add_constraint(consexp[i] == 0) 

521 elif tp == "Q": 

522 P.add_constraint(abs(consexp[i][1:]) < consexp[i][0]) 

523 elif tp == "QR": 

524 P.add_constraint( 

525 abs(consexp[i][2:]) ** 2 < 2 * consexp[i][0] * consexp[i][1] 

526 ) 

527 else: 

528 raise Exception("unexpected cone type") 

529 

530 DD = {} 

531 HH = {} 

532 if "PSDCON" in parsed_blocks: 

533 Dblocks = parsed_blocks.get( 

534 "DCOORD", 

535 [_spmatrix([], [], [], (ni, ni)) for ni in parsed_blocks["PSDCON"]], 

536 ) 

537 Hblocks = parsed_blocks.get("HCOORD", {}) 

538 

539 consexp = [] 

540 for i, Di in enumerate(Dblocks): 

541 DD[i] = Constant("D[" + str(i) + "]", Di) 

542 consexp.append(Constant("D[" + str(i) + "]", Di)) 

543 

544 for j, Hj in Hblocks.items(): 

545 i, col = _block_idx(j, varsz) 

546 for k, Hij in enumerate(Hj): 

547 if Hij: 

548 HH[k, i, col] = Constant("H[" + str((k, i, col)) + "]", Hij) 

549 consexp[k] += HH[k, i, col] * x[i][col] 

550 

551 for exp in consexp: 

552 P.add_constraint(exp >> 0) 

553 

554 params = { 

555 "aobj": aobj, 

556 "bobj": bobj, 

557 "Fobj": Fobj, 

558 "A": AA, 

559 "b": bb, 

560 "F": FF, 

561 "D": DD, 

562 "H": HH, 

563 } 

564 

565 return P, x, X, params 

566 

567 

568# -------------------------------------- 

569__all__ = api_end(_API_START, globals())