Coverage for picos/modeling/file_in.py: 3.88%
361 statements
« prev ^ index » next coverage.py v7.6.12, created at 2025-04-12 07:53 +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, IntegerVariable,
34 RealVariable, SymmetricVariable)
35from .problem import Problem
37_API_START = api_start(globals())
38# -------------------------------
41def _spmatrix(*args, **kwargs):
42 """Create a CVXOPT sparse matrix.
44 A wrapper around :func:`cvxopt.spmatrix` that converts indices to
45 :class:`int`, if necessary.
47 Works around PICOS sometimes passing indices as :class:`numpy:numpy.int64`.
48 """
49 try:
50 return cvxopt.spmatrix(*args, **kwargs)
51 except TypeError as error:
52 # CVXOPT does not like NumPy's int64 scalar type for indices, so attempt
53 # to convert all indices to Python's int.
54 if str(error) == "non-numeric type in list":
55 newargs = list(args)
57 for argNum, arg in enumerate(args):
58 if argNum in (1, 2): # Positional I, J.
59 newargs[argNum] = [int(x) for x in args[argNum]]
61 for kw in "IJ":
62 if kw in kwargs:
63 kwargs[kw] = [int(x) for x in kwargs[kw]]
65 return cvxopt.spmatrix(*newargs, **kwargs)
66 else:
67 raise
70def _break_cols(mat, sizes):
71 """Help with the reading of CBF files."""
72 n = len(sizes)
73 I, J, V = [], [], []
74 for i in range(n):
75 I.append([])
76 J.append([])
77 V.append([])
78 cumsz = numpy.cumsum(sizes)
79 import bisect
81 for i, j, v in zip(mat.I, mat.J, mat.V):
82 block = bisect.bisect(cumsz, j)
83 I[block].append(i)
84 V[block].append(v)
85 if block == 0:
86 J[block].append(j)
87 else:
88 J[block].append(j - cumsz[block - 1])
89 return [
90 _spmatrix(V[k], I[k], J[k], (mat.size[0], sz))
91 for k, sz in enumerate(sizes)
92 ]
95def _break_rows(mat, sizes):
96 """Help with the reading of CBF files."""
97 n = len(sizes)
98 I, J, V = [], [], []
99 for i in range(n):
100 I.append([])
101 J.append([])
102 V.append([])
103 cumsz = numpy.cumsum(sizes)
104 import bisect
106 for i, j, v in zip(mat.I, mat.J, mat.V):
107 block = bisect.bisect(cumsz, i)
108 J[block].append(j)
109 V[block].append(v)
110 if block == 0:
111 I[block].append(i)
112 else:
113 I[block].append(i - cumsz[block - 1])
114 return [
115 _spmatrix(V[k], I[k], J[k], (sz, mat.size[1]))
116 for k, sz in enumerate(sizes)
117 ]
120def _block_idx(i, sizes):
121 """Help with the reading of CBF files."""
122 # if there are blocks of sizes n1,...,nk and i is
123 # the index of an element of the big vectorized variable,
124 # returns the block of i and its index inside the sub-block.
125 cumsz = numpy.cumsum(sizes)
126 import bisect
128 block = bisect.bisect(cumsz, i)
129 return block, (i if block == 0 else i - cumsz[block - 1])
132def _read_cbf_block(P, blocname, f, parsed_blocks):
133 """Help with the reading of CBF files."""
134 if blocname == "OBJSENSE":
135 objsense = f.readline().split()[0].lower()
136 P.set_objective(objsense, P.objective.normalized.function)
137 return None
138 elif blocname == "PSDVAR":
139 n = int(f.readline())
140 vardims = []
141 XX = []
142 for i in range(n):
143 ni = int(f.readline())
144 vardims.append(ni)
145 Xi = SymmetricVariable("X[" + str(i) + "]", (ni, ni))
146 XX.append(Xi)
147 P.add_constraint(Xi >> 0)
148 return vardims, XX
149 elif blocname == "VAR":
150 Nscalar, ncones = [int(fi) for fi in f.readline().split()]
151 tot_dim = 0
152 var_structure = []
153 xx = []
154 for i in range(ncones):
155 lsplit = f.readline().split()
156 tp, dim = lsplit[0], int(lsplit[1])
157 tot_dim += dim
158 var_structure.append(dim)
159 if tp == "F":
160 xi = RealVariable("x[" + str(i) + "]", dim)
161 elif tp == "L+":
162 xi = RealVariable("x[" + str(i) + "]", dim, lower=0)
163 elif tp == "L-":
164 xi = RealVariable("x[" + str(i) + "]", dim, upper=0)
165 elif tp == "L=":
166 xi = RealVariable("x[" + str(i) + "]", dim, lower=0, upper=0)
167 elif tp == "Q":
168 xi = RealVariable("x[" + str(i) + "]", dim)
169 P.add_constraint(abs(xi[1:]) < xi[0])
170 elif tp == "QR":
171 xi = RealVariable("x[" + str(i) + "]", dim)
172 P.add_constraint(abs(xi[2:]) ** 2 < 2 * xi[0] * xi[1])
173 xx.append(xi)
174 if tot_dim != Nscalar:
175 raise Exception("VAR dimensions do not match the header")
176 return Nscalar, var_structure, xx
177 elif blocname == "INT":
178 n = int(f.readline())
179 ints = {}
180 for k in range(n):
181 j = int(f.readline())
182 i, col = _block_idx(j, parsed_blocks["VAR"][1])
183 ints.setdefault(i, [])
184 ints[i].append(col)
185 x = parsed_blocks["VAR"][2]
186 for i in ints:
187 if len(ints[i]) == x[i].size[0]:
188 x[i] = IntegerVariable(
189 x[i].name, x[i].shape, lower=x[i]._lower, upper=x[i]._upper)
190 else:
191 x.append(
192 IntegerVariable("x_int[" + str(i) + "]", len(ints[i]))
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 (
419 "PSDVAR" not in seen_blocks
420 ):
421 raise Exception("missing PSDVAR block")
422 if lsplit[0] in ("OBJACOORD", "ACOORD", "HCOORD") and (
423 "VAR" not in seen_blocks
424 ):
425 raise Exception("missing VAR block")
426 if lsplit[0] in ("DCOORD", "HCOORD") and (
427 "PSDCON" not 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())