Coverage for picos/constraints/con_vecnorm.py: 65.99%
147 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) 2018 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"""Implementation of :class:`VectorNormConstraint`."""
22from collections import namedtuple
24from .. import glyphs
25from ..apidoc import api_end, api_start
26from .constraint import Constraint, ConstraintConversion
28_API_START = api_start(globals())
29# -------------------------------
32class VectorNormConstraint(Constraint):
33 """Bound on a (generalized) vector :math:`p`-norm."""
35 class Conversion(ConstraintConversion):
36 """Bound on a (generalized) :math:`p`-norm constraint conversion."""
38 @classmethod
39 def predict(cls, subtype, options):
40 """Implement :meth:`~.constraint.ConstraintConversion.predict`."""
41 from ..expressions import RealVariable as RV
42 from . import AffineConstraint as AC, GeometricMeanConstraint as GM
44 n, num, den, relation = subtype
45 p = float(num) / float(den)
47 assert p >= 0.0
49 if relation == Constraint.LE:
50 assert p >= 1.0
52 if p == 1.0:
53 yield ("var", RV.make_var_type(dim=n, bnd=0), 1)
54 yield ("con", AC.make_type(dim=n, eq=False), 2)
55 yield ("con", AC.make_type(dim=1, eq=False), 1)
56 elif p == float("inf"):
57 yield ("con", AC.make_type(dim=n, eq=False), 2)
58 else: # 1 < p < inf
59 yield ("var", RV.make_var_type(dim=n, bnd=0), 2)
60 yield ("con", AC.make_type(dim=1, eq=False), 2*n)
61 yield ("con", GM.make_type(argdim=num), n)
62 yield ("con", AC.make_type(dim=1, eq=False), 1)
63 else:
64 assert relation == Constraint.GE
65 assert p <= 1.0
67 if p == 1.0:
68 yield ("con", AC.make_type(dim=n, eq=False), 1)
69 yield ("con", AC.make_type(dim=1, eq=False), 1)
70 elif p == float("-inf"):
71 yield ("con", AC.make_type(dim=n, eq=False), 1)
72 elif p > 0:
73 yield ("var", RV.make_var_type(dim=n, bnd=0), 1)
74 yield ("con", GM.make_type(argdim=den), n)
75 yield ("con", AC.make_type(dim=1, eq=False), 1)
76 else: # -inf < p < 0
77 yield ("var", RV.make_var_type(dim=n, bnd=0), 1)
78 yield ("con", GM.make_type(argdim=(num + den)), n)
79 yield ("con", AC.make_type(dim=1, eq=False), 1)
81 @classmethod
82 def convert(cls, con, options):
83 """Implement :meth:`~.constraint.ConstraintConversion.convert`."""
84 from ..expressions import GeometricMean
85 from ..modeling import Problem
87 norm = con.norm
88 relation = con.relation
89 rhs = con.rhs
91 p = norm.p
92 m = len(norm.x)
94 P = Problem()
96 if relation == Constraint.LE:
97 if p == 1:
98 v = P.add_variable('__v', m)
99 P.add_constraint(norm.x[:] <= v)
100 P.add_constraint(-norm.x[:] <= v)
101 P.add_constraint((1 | v) <= rhs)
102 elif p == float('inf'):
103 P.add_constraint(norm.x <= rhs)
104 P.add_constraint(-norm.x <= rhs)
105 else:
106 x = P.add_variable('__x', m)
107 v = P.add_variable('__v', m)
108 amb = norm.pnum - norm.pden
109 b = norm.pden
110 oneamb = '|1|(' + str(amb) + ',1)'
111 oneb = '|1|(' + str(b) + ',1)'
112 for i in range(m):
113 if amb > 0:
114 if b == 1:
115 vec = (v[i]) // (rhs * oneamb)
116 else:
117 vec = (v[i] * oneb) // (rhs * oneamb)
118 else:
119 # TODO: This shouldn't be possible?
120 if b == 1:
121 vec = v[i]
122 else:
123 vec = (v[i] * oneb)
124 P.add_constraint(norm.x[i] <= x[i])
125 P.add_constraint(-norm.x[i] <= x[i])
126 P.add_constraint(x[i] <= GeometricMean(vec))
127 P.add_constraint((1 | v) <= rhs)
128 else:
129 if p == 1:
130 P.add_constraint(norm.x >= 0)
131 P.add_constraint((1 | norm.x) >= rhs)
132 elif p == float('-inf'):
133 P.add_constraint(norm.x >= rhs)
134 elif p >= 0:
135 v = P.add_variable('__v', m)
136 bma = -(norm.pnum - norm.pden)
137 a = norm.pnum
138 onebma = '|1|(' + str(bma) + ',1)'
139 onea = '|1|(' + str(a) + ',1)'
140 for i in range(m):
141 if a == 1:
142 vec = (norm.x[i]) // (rhs * onebma)
143 else:
144 vec = (norm.x[i] * onea) // (rhs * onebma)
145 P.add_constraint(v[i] <= GeometricMean(vec))
146 P.add_constraint(rhs <= (1 | v))
147 else:
148 v = P.add_variable('__v', m)
149 b = abs(norm.pden)
150 a = abs(norm.pnum)
151 oneb = '|1|(' + str(b) + ',1)'
152 onea = '|1|(' + str(a) + ',1)'
153 for i in range(m):
154 if a == 1 and b == 1:
155 vec = (norm.x[i]) // (v[i])
156 elif a > 1 and b == 1:
157 vec = (norm.x[i] * onea) // (v[i])
158 elif a == 1 and b > 1:
159 vec = (norm.x[i]) // (v[i] * oneb)
160 else:
161 vec = (norm.x[i] * onea) // (v[i] * oneb)
162 P.add_constraint(rhs <= GeometricMean(vec))
163 P.add_constraint((1 | v) <= rhs)
165 return P
167 def __init__(self, norm, relation, rhs):
168 """Construct a :class:`VectorNormConstraint`.
170 :param ~picos.expressions.Norm norm:
171 Left hand side expression.
172 :param str relation:
173 Constraint relation symbol.
174 :param ~picos.expressions.AffineExpression rhs:
175 Right hand side expression.
176 """
177 from ..expressions import AffineExpression, Norm
179 assert isinstance(norm, Norm)
180 assert isinstance(rhs, AffineExpression)
181 assert relation in self.LE + self.GE
182 assert len(rhs) == 1
184 assert norm.q == norm.p, \
185 "Can't create a p-norm constraint for a (p,q)-norm."
187 if relation == self.LE:
188 assert norm.convex, \
189 "Upper bounding a p-norm requires p s.t. the norm is convex."
191 if relation == self.GE:
192 assert norm.concave, \
193 "Lower bounding a p-norm requires p s.t. the norm is concave."
195 self.norm = norm
196 self.relation = relation
197 self.rhs = rhs
199 super(VectorNormConstraint, self).__init__(norm._typeStr)
201 Subtype = namedtuple("Subtype", ("argdim", "num", "den", "relation"))
203 def _subtype(self):
204 return self.Subtype(
205 len(self.norm.x), self.norm.pnum, self.norm.pden, self.relation)
207 @classmethod
208 def _cost(cls, subtype):
209 return subtype.argdim + 1
211 def _expression_names(self):
212 yield "norm"
213 yield "rhs"
215 def _str(self):
216 if self.relation == self.LE:
217 str = glyphs.le(self.norm.string, self.rhs.string)
218 else:
219 str = glyphs.ge(self.norm.string, self.rhs.string)
221 if self.norm.p < 1:
222 return glyphs.and_(str, glyphs.ge(self.norm.x.string, 0))
223 else:
224 return str
226 def _get_slack(self):
227 if self.relation == self.LE:
228 return self.rhs.safe_value - self.norm.safe_value
229 else:
230 return self.norm.safe_value - self.rhs.safe_value
233# --------------------------------------
234__all__ = api_end(_API_START, globals())