...
 
Commits (20)
......@@ -26,7 +26,7 @@ doctest:
- >
pacman -Syu --noconfirm
m2r python-cvxopt python-matplotlib python-networkx python-numpy
python-pip python-sphinx
python-pip python-sphinx python-sphinx_rtd_theme
- pip install sphinx-automodapi
script:
- sphinx-build -b doctest doc doctest
......@@ -91,8 +91,9 @@ pdfdoc:
before_script:
- >
pacman -Syu --noconfirm
git glpk graphviz make m2r python-cvxopt python-matplotlib python-networkx
python-numpy python-pip python-scipy python-sphinx texlive-most
git glpk graphviz make m2r python-cvxopt python-matplotlib
python-networkx python-numpy python-pip python-scipy python-sphinx
python-sphinx_rtd_theme texlive-most
- pip install sphinx-automodapi swiglpk
script:
- make -C doc latexpdf
......
......@@ -21,7 +21,19 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
## [Unreleased]
### Added
- Support for MOSEK 9.
- Support for ECOS 2.0.7.
- Support for multiple subsystems with
:func:`partial_trace <picos.tools.partial_trace>`.
- Possibility to use :func:`picos.sum` to sum over the elements of a single
multidimensional expression.
### Fixed
- Upgrading the PyPI package via pip.
- A regression that rendered the Kronecker product unusable.
- Noisy exception handling in :func:`spmatrix <picos.tools.spmatrix>`.
- Shape detection for matrices given by string.
- The ``hotstart`` option when solving with CPLEX.
## [1.2.0] - 2019-01-11
......
......@@ -181,5 +181,6 @@ def setup(app):
(u'∑', em('\\sum{}')),
(u'∈', em('\\in{}')),
(u'…', '\\ldots{}'),
(u'∀', em('\\forall{}'))
(u'∀', em('\\forall{}')),
(u'λ', em('\\lambda{}'))
]
......@@ -1109,20 +1109,30 @@ class AffinExp(Expression):
"""
Partial trace, see also :func:`the partial_trace tool
<picos.tools.partial_trace>`
Ref. J. Maziero, Computing partial traces and reduced density matrices, Int. J. Mod. Phys. C 28, 1750005 (2017)
"""
T = [list,tuple]
sz = self.size
<<<<<<< HEAD
dim_ = dim
=======
#If dim is None, we assume a bipartite system
>>>>>>> upstream/master
if dim is None:
if type(k) in T:
if len(k) > 1:
raise UserWarning("For partial_trace over mutliple system, dim is required")
if sz[0] == sz[1] and (sz[0]**0.5) == int(sz[0]**0.5) \
and (sz[1]**0.5) == int(sz[1]**0.5):
dim = (int(sz[0]**0.5), int(sz[1]**0.5))
else:
raise Exception("The default parameter dim=None assumes X is a "
"{} matrix".format(glyphs.size(*[gylphs.cubed("n")]*2)))
"{} matrix".format(glyphs.size(*[glyphs.cubed("n")]*2)))
else:
dim_ = list(dim)
# checks if dim is a list (or tuple) of lists (or tuples) of two
# integers each
T = [list,tuple]
if type(dim) in T and all([type(d) in T and len(d) == 2 for d in dim]) \
and all([type(n) is int for d in dim for n in d]):
dim = [d for d in zip(*dim)]
......@@ -1161,7 +1171,11 @@ class AffinExp(Expression):
for i,partition in enumerate(sorted(k)):
self = self.partial_trace(k=partition-i,dim=dim_)
del dim_[partition-i]
<<<<<<< HEAD
self.string = glyphs.ptrace("".join(['{},'.format(p) for p in k[:-1]]+[str(k[-1])]), str_)
=======
self.string = glyphs.ptrace(",".join([str(p) for p in k]), str_)
>>>>>>> upstream/master
return self
if type(k) is int:
......
......@@ -695,7 +695,7 @@ class Problem(object):
* an ``int`` :math:`n`, in which case the variable is an
:math:`n`-dimensional vector,
* or a ``tuple`` :math:`(n, m)`, in which case the variable is a
:math:`n × m` matrix.
:math:`n \times{} m` matrix.
:type size: int or tuple
:param str vtype: Domain of the variable.
......@@ -1977,7 +1977,9 @@ class Problem(object):
duals = []
for cstNum, cst in enumerate(realP.constraints):
if self.constraints[cstNum].is_complex():
if isinstance(cst, LMIConstraint):
if cst.dual is None:
duals.append(None)
elif isinstance(cst, LMIConstraint):
n = int(cst.dual.size[0] / 2.)
if cst.dual.size[1] != cst.dual.size[0]:
raise Exception('Dual must be square matrix')
......@@ -3443,8 +3445,9 @@ class Problem(object):
def as_real(self):
"""
Returns a modified copy of the problem, where hermitian :math:`n × n`
matrices are replaced by symmetric :math:`2n × 2n` matrices.
Returns a modified copy of the problem, where hermitian
:math:`n \times{} n` matrices are replaced by symmetric
:math:`2n \times{} 2n` matrices.
"""
import copy
real = Problem()
......
......@@ -181,7 +181,7 @@ class CPLEXSolver(Solver):
self.int.variables.add(
lb=lowerBounds, ub=upperBounds, types=types, names=names)
# Map PICOS indices to CPLEX indices.
# Map PICOS indices to CPLEX names.
for localIndex in range(varDim):
picosIndex = picosVar.startIndex + localIndex
self._cplexVarName[picosIndex] = names[localIndex]
......@@ -215,7 +215,7 @@ class CPLEXSolver(Solver):
quadratic expression.
:returns: :class:`SparseTriple <cplex.SparseTriple>` mapping a pair of
CPLEX variable indices to scalar constants.
CPLEX variable names to scalar constants.
"""
import cplex
......@@ -474,7 +474,7 @@ class CPLEXSolver(Solver):
def _import_constraint(self, picosConstraint):
# Import constraint based on type and keep track of the corresponding
# CPLEX constraint and auxiliary variable indices.
# CPLEX constraint and auxiliary variable names.
if isinstance(picosConstraint, AffineConstraint):
self._cplexLinConNames[picosConstraint] = \
self._import_linear_constraint(picosConstraint)
......@@ -732,8 +732,8 @@ class CPLEXSolver(Solver):
if picosVar.is_valued():
for localIndex in range(picosVar.dim):
picosIndex = picosVar.startIndex + localIndex
cplexName = _cplexVarName[picosIndex]
indices.append(cplexName)
cplexName = self._cplexVarName[picosIndex]
names.append(cplexName)
values.append(picosVar.valueAsMatrix[localIndex])
if names:
self.int.MIP_starts.add(
......
......@@ -20,7 +20,7 @@
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# This file implements the MOSEK solver through its low level Optimizer API.
# This file implements the MOSEK 8/9 solver through its low level Optimizer API.
# The low level API is tedious to interface, but is currently much faster than
# the high level Fusion API, which would be the prefered interface otherwise.
#-------------------------------------------------------------------------------
......@@ -126,9 +126,8 @@ class MOSEKSolver(Solver):
@classmethod
def _get_environment(cls):
import mosek
if not hasattr(cls, "mosekEnvironment"):
import mosek
cls.mosekEnvironment = mosek.Env()
return cls.mosekEnvironment
......@@ -138,38 +137,60 @@ class MOSEKSolver(Solver):
MOSEKSolver instances. (The MOSEK documentation states that "[a]ll tasks in
the program should share the same environment.")"""
@staticmethod
def _streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
@classmethod
def _get_major_version(cls):
if not hasattr(cls, "mosekVersion"):
import mosek
cls.mosekVersion = mosek.Env.getversion()
@staticmethod
def _status_string(statusCode):
from mosek import solsta
return cls.mosekVersion[0]
ver = property(lambda self: self.__class__._get_major_version())
"""The major version of the available MOSEK library."""
@classmethod
def _get_solution_status_map(cls):
if not hasattr(cls, "solstaMap"):
from mosek import solsta
cls.solstaMap = {
solsta.unknown: "unknown",
solsta.optimal: "optimal",
solsta.prim_feas: "primal feasible",
solsta.dual_feas: "dual feasible",
solsta.prim_and_dual_feas: "feasible",
solsta.prim_infeas_cer: "primal infeasible",
solsta.dual_infeas_cer: "dual infeasible",
solsta.prim_illposed_cer: "primal illposed",
solsta.dual_illposed_cer: "dual illposed",
solsta.integer_optimal: "integer optimal"
}
if cls._get_major_version() < 9:
cls.solstaMap.update({
solsta.near_optimal: "near optimal",
solsta.near_prim_feas: "near primal feasible",
solsta.near_dual_feas: "near dual feasible",
solsta.near_prim_and_dual_feas: "near feasible",
solsta.near_prim_infeas_cer: "near primal infeasible",
solsta.near_dual_infeas_cer: "near dual infeasible",
solsta.near_integer_optimal: "near integer optimal"
})
return cls.solstaMap
@classmethod
def _status_string(cls, statusCode):
try:
return {
solsta.unknown: "unknown",
solsta.optimal: "optimal",
solsta.prim_feas: "primal feasible",
solsta.dual_feas: "dual feasible",
solsta.prim_and_dual_feas: "feasible",
solsta.near_optimal: "near optimal",
solsta.near_prim_feas: "near primal feasible",
solsta.near_dual_feas: "near dual feasible",
solsta.near_prim_and_dual_feas: "near feasible",
solsta.prim_infeas_cer: "primal infeasible",
solsta.dual_infeas_cer: "dual infeasible",
solsta.near_prim_infeas_cer: "near primal infeasible",
solsta.near_dual_infeas_cer: "near dual infeasible",
solsta.prim_illposed_cer: "primal illposed",
solsta.dual_illposed_cer: "dual illposed",
solsta.integer_optimal: "integer optimal",
solsta.near_integer_optimal: "near integer optimal"
}[statusCode]
return cls._get_solution_status_map()[statusCode]
except KeyError:
return "unknown"
@staticmethod
def _streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
@staticmethod
def _low_tri_indices(rowCount):
"""
......@@ -588,7 +609,7 @@ class MOSEKSolver(Solver):
# Domains:
# - "": linear,
# - "co": conic,
# - "nl": nonlinear
# - "nl": nonlinear (MOSEK < 9)
# - "qo": quadratic
# Tolerances:
# - "pfeas": primal feasibility,
......@@ -598,7 +619,7 @@ class MOSEKSolver(Solver):
# Omitted tolerances:
# - "infeas": infeasibility,
# - "near_rel": (not actually a tolerance)
for domain in ("", "co", "nl", "qo"):
for domain in ("", "co", "qo") + (("nl",) if self.ver < 9 else ()):
for tolerance in ("pfeas", "dfeas", "mu_red", "rel_gap"):
parameter = "msk_dpar_intpnt_{}{}tol_{}".format(
domain, "_" if domain else "", tolerance).upper()
......
This diff is collapsed.
......@@ -52,7 +52,9 @@ def sum(lst, it=None, indices=None):
This is a replacement for Python's :func:`sum` that produces sensible string
representations when summing PICOS expressions.
:param lst: A list of :class:`expressions <picos.expressions.Expression>`.
:param lst: A list of :class:`expressions <picos.expressions.Expression>`,
or a single affine expression whose elements shall be summed.
:type lst: list or tuple or AffinExp
:param it: DEPRECATED
:param indices: DEPRECATED
......@@ -66,6 +68,8 @@ def sum(lst, it=None, indices=None):
<Quadratic Expression: x[0]·x[1] + x[1]·x[2] + x[2]·x[3] + x[3]·x[4]>
>>> picos.sum(e)
<Quadratic Expression: ∑(x[i]·x[i+1] : i ∈ [0…3])>
>>> picos.sum(x) # The same as (x|1).
<1×1 Affine Expression: ∑(x)>
"""
from .expressions import Expression, AffinExp
......@@ -73,6 +77,15 @@ def sum(lst, it=None, indices=None):
warn("Arguments 'it' and 'indices' to picos.sum are ignored.",
DeprecationWarning)
if isinstance(lst, Expression):
if isinstance(lst, AffinExp):
theSum = (lst|1)
theSum.string = glyphs.sum(lst.string)
return theSum
else:
raise TypeError("PICOS doesn't know how to sum over a single {}."
.format(type(lst).__name__))
if not all(isinstance(expression, Expression) for expression in lst):
return builtin_sum(lst)
......@@ -483,7 +496,7 @@ def partial_transpose(exp, dims_1=None, subsystems = None, dims_2=None):
def partial_trace(X, k=1, dim=None):
r"""Partial trace of an Affine Expression, with respect to the ``k`` th subsystem for a tensor product of dimensions ``dim``.
r"""Partial trace of an Affine Expression, with respect to the ``k``th subsystem or to every subsystems in ``k`` for a tensor product of dimensions ``dim``.
If ``X`` is a matrix
:class:`AffinExp <picos.expressions.AffinExp>`
that can be written as :math:`X = A_0 \otimes \cdots \otimes A_{n-1}`
......@@ -491,11 +504,11 @@ def partial_trace(X, k=1, dim=None):
of respective sizes ``dim[0] x dim[0]``, ... , ``dim[n-1] x dim[n-1]`` (``dim`` is a list of ints if all matrices are square),
or ``dim[0][0] x dim[0][1]``, ...,``dim[n-1][0] x dim[n-1][1]`` (``dim`` is a list of 2-tuples if any of them except the ``k`` th one is rectangular),
this function returns the matrix
:math:`Y = \operatorname{trace}(A_k)\quad A_0 \otimes \cdots A_{k-1} \otimes A_{k+1} \otimes \cdots \otimes A_{n-1}`.
:math:`Y = \operatorname{trace}(A_k)\quad A_0 \otimes \cdots A_{k-1} \otimes A_{k+1} \otimes \cdots \otimes A_{n-1}` if ``k`` is an ``int`` or a ``list`` of one element, :math:`Y = \operatorname{trace}(A_{k_0,k_1,\ldots,k_m}) \quad A_0 \otimes \cdots A_{k_0-1} \otimes A_{k_0+1} \otimes \cdots \otimes A_{k_m-1} \otimes A_{k_m+1} \otimes \cdots \otimes A_{n-1}` if ``k`` is a list of subsystems.
The default value ``dim=None`` automatically computes the size of the subblocks,
assuming that ``X`` is a :math:`n^2 \times n^2`-square matrix
with blocks of size :math:`n \times n`.
with blocks of size :math:`n \times n`, in this case only bipartite subsystem is assume, hence ``k`` as to be either ``0`` or ``1``.
**Example:**
......@@ -503,7 +516,9 @@ def partial_trace(X, k=1, dim=None):
>>> import cvxopt as cvx
>>> P = pic.Problem()
>>> X = P.add_variable('X',(4,4))
>>> Y = P.add_variable('Y',(8,8))
>>> X.value = cvx.matrix(range(16),(4,4))
>>> Y.value = cvx.matrix(range(64),(8,8))
>>> print(X) #doctest: +NORMALIZE_WHITESPACE
[ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
[ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
......@@ -517,6 +532,10 @@ def partial_trace(X, k=1, dim=None):
>>> print(pic.partial_trace(X,0)) #doctest: +NORMALIZE_WHITESPACE
[ 1.00e+01 1.80e+01]
[ 1.20e+01 2.00e+01]
>>> #Partial trace in a multiple subsystem scenario:
>>> print(pic.partial_trace(Y,k=[0,2],dim=[2,2,2])) #doctest: +NORMALIZE_WHITESPACE
[ 9.00e+01 1.54e+02]
[ 9.80e+01 1.62e+02]
"""
return X.partial_trace(k, dim)
......@@ -1343,7 +1362,10 @@ def retrieve_matrix(mat, exSize=None):
forcedSize = int(forcedSize[0]), int(forcedSize[1])
size = forcedSize
else:
size = exSize
if is_integer(exSize):
size = (exSize, exSize)
else:
size = exSize
if size is None:
raise ValueError("The size of the target matrix was not specified.")
......@@ -2155,48 +2177,33 @@ def is_realvalued(x):
isinstance(x, np.int64) or
isinstance(x, np.int32))
def spmatrix(*args,**kwargs):
def spmatrix(*args, **kwargs):
"""
A wrapper around :func:`cvxopt.spmatrix` that converts indices to
:class:`int`, if necessary.
It works around PICOS sometimes passing indices as `numpy.int64`.
"""
try:
return cvx.spmatrix(*args,**kwargs)
except TypeError as ex:
print(type(ex))
print(str(ex))
print('non-numeric' in str(ex))
if 'non-numeric' in str(ex):#catch exception with int64 bug of cvxopt
size_tc = {}
if len(args)>0:
V = args[0]
elif 'V' in kwargs:
V = kwargs['V']
else:
V = []
if len(args)>1:
I = args[1]
elif 'I' in kwargs:
I = kwargs['I']
else:
I = []
if len(args)>2:
J = args[2]
elif 'J' in kwargs:
J = kwargs['J']
else:
J = []
if len(args) > 3:
size_tc['size'] = args[3]
elif 'size' in kwargs:
size_tc['size'] = kwargs['size']
if len(args) > 4:
size_tc['tc'] = args[4]
elif 'tc' in kwargs:
size_tc['tc'] = kwargs['tc']
return cvx.spmatrix(V, [int(i) for i in I], [int(j) for j in J],**size_tc)
return cvx.spmatrix(*args, **kwargs)
except TypeError as error:
# CVXOPT does not like NumPy's int64 scalar type for indices, so attempt
# to convert all indices to Python's int.
if str(error) == "non-numeric type in list":
newargs = list(args)
for argNum, arg in enumerate(args):
if argNum in (1, 2): # Positional I, J.
newargs[argNum] = [int(x) for x in args[argNum]]
for kw in "IJ":
if kw in kwargs:
kwargs[kw] = [int(x) for x in kwargs[kw]]
return cvx.spmatrix(*newargs, **kwargs)
else:
raise
except Exception as ex:
print(type(ex))
import pdb;pdb.set_trace()
def kron(A,B):
......