Coverage for picos/constraints/con_vecnorm.py: 65.99%

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# ------------------------------------------------------------------------------ 

19 

20"""Implementation of :class:`VectorNormConstraint`.""" 

21 

22from collections import namedtuple 

23 

24from .. import glyphs 

25from ..apidoc import api_end, api_start 

26from .constraint import Constraint, ConstraintConversion 

27 

28_API_START = api_start(globals()) 

29# ------------------------------- 

30 

31 

32class VectorNormConstraint(Constraint): 

33 """Bound on a (generalized) vector :math:`p`-norm.""" 

34 

35 class Conversion(ConstraintConversion): 

36 """Bound on a (generalized) :math:`p`-norm constraint conversion.""" 

37 

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 

43 

44 n, num, den, relation = subtype 

45 p = float(num) / float(den) 

46 

47 assert p >= 0.0 

48 

49 if relation == Constraint.LE: 

50 assert p >= 1.0 

51 

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 

66 

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) 

80 

81 @classmethod 

82 def convert(cls, con, options): 

83 """Implement :meth:`~.constraint.ConstraintConversion.convert`.""" 

84 from ..expressions import GeometricMean 

85 from ..modeling import Problem 

86 

87 norm = con.norm 

88 relation = con.relation 

89 rhs = con.rhs 

90 

91 p = norm.p 

92 m = len(norm.x) 

93 

94 P = Problem() 

95 

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) 

164 

165 return P 

166 

167 def __init__(self, norm, relation, rhs): 

168 """Construct a :class:`VectorNormConstraint`. 

169 

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 

178 

179 assert isinstance(norm, Norm) 

180 assert isinstance(rhs, AffineExpression) 

181 assert relation in self.LE + self.GE 

182 assert len(rhs) == 1 

183 

184 assert norm.q == norm.p, \ 

185 "Can't create a p-norm constraint for a (p,q)-norm." 

186 

187 if relation == self.LE: 

188 assert norm.convex, \ 

189 "Upper bounding a p-norm requires p s.t. the norm is convex." 

190 

191 if relation == self.GE: 

192 assert norm.concave, \ 

193 "Lower bounding a p-norm requires p s.t. the norm is concave." 

194 

195 self.norm = norm 

196 self.relation = relation 

197 self.rhs = rhs 

198 

199 super(VectorNormConstraint, self).__init__(norm._typeStr) 

200 

201 Subtype = namedtuple("Subtype", ("argdim", "num", "den", "relation")) 

202 

203 def _subtype(self): 

204 return self.Subtype( 

205 len(self.norm.x), self.norm.pnum, self.norm.pden, self.relation) 

206 

207 @classmethod 

208 def _cost(cls, subtype): 

209 return subtype.argdim + 1 

210 

211 def _expression_names(self): 

212 yield "norm" 

213 yield "rhs" 

214 

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) 

220 

221 if self.norm.p < 1: 

222 return glyphs.and_(str, glyphs.ge(self.norm.x.string, 0)) 

223 else: 

224 return str 

225 

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 

231 

232 

233# -------------------------------------- 

234__all__ = api_end(_API_START, globals())