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# ------------------------------------------------------------------------------
20"""Functions for reading optimization problems from a file."""
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# ------------------------------------------------------------------------------
27import warnings
29import cvxopt
30import numpy
32from ..apidoc import api_end, api_start
33from ..expressions import AffineExpression, Constant
34from .problem import Problem
36_API_START = api_start(globals())
37# -------------------------------
40def _spmatrix(*args, **kwargs):
41 """Create a CVXOPT sparse matrix.
43 A wrapper around :func:`cvxopt.spmatrix` that converts indices to
44 :class:`int`, if necessary.
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)
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]]
60 for kw in "IJ":
61 if kw in kwargs:
62 kwargs[kw] = [int(x) for x in kwargs[kw]]
64 return cvxopt.spmatrix(*newargs, **kwargs)
65 else:
66 raise
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
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 ]
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
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 ]
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
127 block = bisect.bisect(cumsz, i)
128 return block, (i if block == 0 else i - cumsz[block - 1])
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")
332def import_cbf(filename):
333 """Create a :class:`~picos.Problem` from a CBF file.
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``.
340 :returns: A tuple ``(P, x, X, params)`` where
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()
353 try:
354 f = open(filename, "r")
355 except IOError:
356 filename += ".cbf"
357 f = open(filename, "r")
359 line = f.readline()
360 while not line.startswith("VER"):
361 line = f.readline()
363 ver = int(f.readline())
364 if ver != 1:
365 warnings.warn("CBF file has a version other than 1.")
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 ]
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])
436 f.close()
437 # variables
438 if "VAR" in parsed_blocks:
439 Nvars, varsz, x = parsed_blocks["VAR"]
440 else:
441 x = None
443 if "INT" in parsed_blocks:
444 x = parsed_blocks["INT"]
446 if "PSDVAR" in parsed_blocks:
447 psdsz, X = parsed_blocks["PSDVAR"]
448 else:
449 X = None
451 # objective
452 obj_constant = parsed_blocks.get("OBJBCOORD", 0)
453 bobj = Constant("bobj", obj_constant)
454 obj = Constant("bobj", obj_constant)
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]
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]
473 P.set_objective(P.objective.normalized.direction, obj)
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]
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))
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]
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]
508 consexp[i] += (
509 "e_" + str(row) + "(" + str(szcons[i]) + ",1)"
510 ) * row_exp
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")
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", {})
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))
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]
551 for exp in consexp:
552 P.add_constraint(exp >> 0)
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 }
565 return P, x, X, params
568# --------------------------------------
569__all__ = api_end(_API_START, globals())