Coverage for picos/modeling/file_out.py: 2.83%
636 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-02-15 14:21 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-02-15 14:21 +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 writing optimization problems to a file."""
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# ------------------------------------------------------------------------------
28from itertools import chain
30import cvxopt
31import numpy
33from ..apidoc import api_end, api_start
34from ..constraints import (
35 AffineConstraint,
36 LMIConstraint,
37 RSOCConstraint,
38 SOCConstraint,
39)
41from ..solvers import CVXOPTSolver, CPLEXSolver, MOSEKSolver, GurobiSolver
42from ..expressions import IntegerVariable, BinaryVariable, CONTINUOUS_VARTYPES
43from ..expressions.vectorizations import FullVectorization
45_API_START = api_start(globals())
46# -------------------------------
49INFINITY = 1e16 #: A number deemed too large to appear in practice.
52def write(picos_problem, filename, writer="picos"):
53 r"""Write an optimization problem to a file.
55 :param ~picos.Problem P: The problem to write.
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:
61 * ``'.cbf'`` -- Conic Benchmark Format.
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.
72 * ``'.lp'`` -- `LP format
73 <http://docs.mosek.com/6.0/pyapi/node022.html>`_.
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.
81 * ``'.mps'`` -- `MPS format
82 <http://docs.mosek.com/6.0/pyapi/node021.html>`_.
84 As the writer, you need to choose one of ``'cplex'``, ``'gurobi'``
85 or ``'mosek'``.
87 * ``'.opf'`` -- `OPF format
88 <http://docs.mosek.com/6.0/pyapi/node023.html>`_.
90 As the writer, you need to choose ``'mosek'``.
92 * ``'.dat-s'`` -- `Sparse SDPA format
93 <http://sdpa.indsys.chuo-u.ac.jp/sdpa/download.html#sdpa>`_.
95 This format is suitable for semidefinite programs. Second order
96 cone constraints are stored as semidefinite constraints on an
97 *arrow shaped* matrix.
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.
105 .. warning::
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.
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.
124 .. warning::
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.
133 from .strategy import NoStrategyFound
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
153 assert 'Linear' in P.type, "only LINEAR integer problems can be " \
154 "exported with the default writer"
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")
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")
192 if writer == "cplex":
193 cpl = CPLEXSolver(P)
194 cpl._load_problem()
195 cpl.int.write(filename)
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)
205 elif writer == "gurobi":
206 grb = GurobiSolver(P)
207 grb._load_problem()
208 grb.int.write(filename)
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")
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")
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()
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
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
290 print("writing problem in " + filename + "...")
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
302 f.write(affexp_writer("obj", I, V))
303 f.write("\n")
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
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
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
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")
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
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
417 f.write("End\n")
418 print("done.")
419 f.close()
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()
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"]
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]
444 for i in range(len(dims["q"])):
445 G = cvxopt.sparse([G, cvxoptVars["Gq"][i]])
446 h = cvxopt.matrix([h, cvxoptVars["hq"][i]])
448 for i in range(len(dims["s"])):
449 G = cvxopt.sparse([G, cvxoptVars["Gs"][i]])
450 h = cvxopt.matrix([h, cvxoptVars["hs"][i]])
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)
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")
461 from cvxopt import sparse, spmatrix
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
477 P_m = G.size[1]
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)
488 # write data
489 # add extension
490 if filename[-6:] != ".dat-s":
491 filename += ".dat-s"
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)
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
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
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
559 f.close()
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
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)
633def _write_cbf(P, filename, uptri=False):
634 """Write the problem to a file in Sparse SDPA format.
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"
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)
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 )
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
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)
684 barvars = bool(idxsdpvars)
686 # find integer variables, retrieve bounds and put 0-1 bounds on binaries
687 ints = []
688 ind = 0
689 varbounds_lo = {}
690 varbounds_up = {}
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)
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
711 else:
712 raise Exception("variable type not handled by _write_cbf()")
713 ind += var.dim
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 )
724 # open file
725 f = open(filename, "w")
726 f.write("#file " + filename + " generated by picos\n")
727 print("writing problem in " + filename + "...")
729 f.write("VER\n")
730 f.write("1\n\n")
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")
738 # VARIABLEs
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")
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:
769 Acoord.append((iaff, ind - offset, 1.0))
770 Bcoord.append((iaff, -l))
771 iaff += 1
772 if any(varbounds_up[ind] != 0 for ind in range(si, ei)):
773 for ind in range(si, ei):
774 u = varbounds_up[ind]
775 if u + 1 < INFINITY:
776 Acoord.append((iaff, ind - offset, -1.0))
777 Bcoord.append((iaff, u))
778 iaff += 1
779 if iaff:
780 conecons.append(("L+", iaff))
782 f.write("VAR\n")
783 f.write(str(NUMVAR_SCALAR) + " " + str(len(cones)) + "\n")
784 for tp, n in cones:
785 f.write(tp + " " + str(n) + "\n")
787 f.write("\n")
789 # integers
790 if ints:
791 f.write("INT\n")
792 f.write(str(len(ints)) + "\n")
793 for i in ints:
794 f.write(str(i) + "\n")
795 f.write("\n")
797 # constraints
798 psdcons = []
799 isdp = 0
800 Fcoord = []
801 Hcoord = []
802 Dcoord = []
803 ObjAcoord = []
804 ObjBcoord = []
805 ObjFcoord = []
807 # dummy constraint for the objective
808 dummy_cons = P.objective.function >= 0
809 setattr(dummy_cons, "dummycon", None)
811 for cons in chain((dummy_cons,), P.constraints.values()):
812 if isinstance(cons, LMIConstraint):
813 v = cons.semidefVar
814 if v is not None and v not in PSD_not_handled:
815 continue
817 # get sparse indices
818 if isinstance(cons, AffineConstraint):
819 expcone = cons.lhs - cons.rhs
820 if hasattr(cons, "dummycon"):
821 conetype = "0" # Dummy type for the objective function.
822 elif cons.is_equality():
823 conetype = "L="
824 elif cons.is_increasing():
825 conetype = "L-"
826 elif cons.is_decreasing():
827 conetype = "L+"
828 else:
829 assert False, "Unexpected constraint relation."
830 elif isinstance(cons, SOCConstraint):
831 expcone = (cons.ub) // (cons.ne[:])
832 conetype = "Q"
833 elif isinstance(cons, RSOCConstraint):
834 expcone = (cons.ub1) // (0.5 * cons.ub2) // (cons.ne[:])
835 conetype = "QR"
836 elif isinstance(cons, LMIConstraint):
837 if cons.is_increasing():
838 expcone = cons.rhs - cons.lhs
839 conetype = None
840 elif cons.is_decreasing():
841 expcone = cons.lhs - cons.rhs
842 conetype = None
843 else:
844 assert False, "Unexpected constraint relation."
845 else:
846 assert False, "Unexpected constraint type."
848 ijv = []
849 for var, fact in expcone._coefs.items():
850 if not isinstance(fact, cvxopt.base.spmatrix):
851 fact = cvxopt.sparse(fact)
852 sj = start_indices[var]
853 ijv.extend(zip(fact.I, fact.J + sj, fact.V))
854 ijvs = sorted(ijv)
856 itojv = {}
857 lasti = -1
858 for (i, j, v) in ijvs:
859 if i == lasti:
860 itojv[i].append((j, v))
861 else:
862 lasti = i
863 itojv[i] = [(j, v)]
865 if conetype:
866 if conetype != "0":
867 dim = expcone.size[0] * expcone.size[1]
868 conecons.append((conetype, dim))
869 else:
870 dim = expcone.size[0]
871 psdcons.append(dim)
873 if conetype:
874 for i, jv in itojv.items():
875 J = [jvk[0] for jvk in jv]
876 V = [jvk[1] for jvk in jv]
877 J, V, mats = _separate_linear_cons_plain_bar_vars(
878 J, V, idxsdpvars)
879 for j, v in zip(J, V):
880 if conetype != "0":
881 Acoord.append((iaff + i, j, v))
882 else:
883 ObjAcoord.append((j, v))
884 for k, mat in enumerate(mats):
885 for row, col, v in zip(mat.I, mat.J, mat.V):
886 if conetype != "0":
887 Fcoord.append((iaff + i, k, row, col, v))
888 else:
889 ObjFcoord.append((k, row, col, v))
890 if uptri and row != col:
891 if conetype != "0":
892 Fcoord.append((iaff + i, k, col, row, v))
893 else:
894 ObjFcoord.append((k, col, row, v))
895 constant = expcone._const
896 if not (constant is None):
897 constant = cvxopt.sparse(constant)
898 for i, v in zip(constant.I, constant.V):
899 if conetype != "0":
900 Bcoord.append((iaff + i, v))
901 else:
902 ObjBcoord.append(v)
903 else:
904 for i, jv in itojv.items():
905 col, row = divmod(i, dim)
906 if not (uptri) and row < col:
907 continue
908 J = [jvk[0] for jvk in jv]
909 V = [jvk[1] for jvk in jv]
910 J, V, mats = _separate_linear_cons_plain_bar_vars(
911 J, V, idxsdpvars)
912 if any([bool(m) for m in mats]):
913 raise Exception("SDP cons should not depend on PSD var")
914 for j, v in zip(J, V):
915 Hcoord.append((isdp, j, row, col, v))
917 constant = expcone._const
918 if not (constant is None):
919 constant = cvxopt.sparse(constant)
920 for i, v in zip(constant.I, constant.V):
921 col, row = divmod(i, dim)
922 if row < col:
923 continue
924 Dcoord.append((isdp, row, col, v))
926 if conetype:
927 if conetype != "0":
928 iaff += dim
929 else:
930 isdp += 1
932 if iaff > 0:
933 f.write("CON\n")
934 f.write(str(iaff) + " " + str(len(conecons)) + "\n")
935 for tp, n in conecons:
936 f.write(tp + " " + str(n))
937 f.write("\n")
939 f.write("\n")
941 if isdp > 0:
942 f.write("PSDCON\n")
943 f.write(str(isdp) + "\n")
944 for n in psdcons:
945 f.write(str(n) + "\n")
946 f.write("\n")
948 if ObjFcoord:
949 f.write("OBJFCOORD\n")
950 f.write(str(len(ObjFcoord)) + "\n")
951 for (k, row, col, v) in ObjFcoord:
952 f.write("{0} {1} {2} {3}\n".format(k, row, col, v))
953 f.write("\n")
955 if ObjAcoord:
956 f.write("OBJACOORD\n")
957 f.write(str(len(ObjAcoord)) + "\n")
958 for (j, v) in ObjAcoord:
959 f.write("{0} {1}\n".format(j, v))
960 f.write("\n")
962 if ObjBcoord:
963 f.write("OBJBCOORD\n")
964 v = ObjBcoord[0]
965 f.write("{0}\n".format(v))
966 f.write("\n")
968 if Fcoord:
969 f.write("FCOORD\n")
970 f.write(str(len(Fcoord)) + "\n")
971 for (i, k, row, col, v) in Fcoord:
972 f.write("{0} {1} {2} {3} {4}\n".format(i, k, row, col, v))
973 f.write("\n")
975 if Acoord:
976 f.write("ACOORD\n")
977 f.write(str(len(Acoord)) + "\n")
978 for (i, j, v) in Acoord:
979 f.write("{0} {1} {2}\n".format(i, j, v))
980 f.write("\n")
982 if Bcoord:
983 f.write("BCOORD\n")
984 f.write(str(len(Bcoord)) + "\n")
985 for (i, v) in Bcoord:
986 f.write("{0} {1}\n".format(i, v))
987 f.write("\n")
989 if Hcoord:
990 f.write("HCOORD\n")
991 f.write(str(len(Hcoord)) + "\n")
992 for (i, j, row, col, v) in Hcoord:
993 f.write("{0} {1} {2} {3} {4}\n".format(i, j, row, col, v))
994 f.write("\n")
996 if Dcoord:
997 f.write("DCOORD\n")
998 f.write(str(len(Dcoord)) + "\n")
999 for (i, row, col, v) in Dcoord:
1000 f.write("{0} {1} {2} {3}\n".format(i, row, col, v))
1001 f.write("\n")
1003 print("done.")
1004 f.close()
1007# --------------------------------------
1008__all__ = api_end(_API_START, globals())