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 

43 from . import GeometricMeanConstraint as GM 

44 

45 n, num, den, relation = subtype 

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

47 

48 assert p >= 0.0 

49 

50 if relation == Constraint.LE: 

51 assert p >= 1.0 

52 

53 if p == 1.0: 

54 yield ("var", RV.make_var_type(dim=n, bnd=0), 1) 

55 yield ("con", AC.make_type(dim=n, eq=False), 2) 

56 yield ("con", AC.make_type(dim=1, eq=False), 1) 

57 elif p == float("inf"): 

58 yield ("con", AC.make_type(dim=n, eq=False), 2) 

59 else: # 1 < p < inf 

60 yield ("var", RV.make_var_type(dim=n, bnd=0), 2) 

61 yield ("con", AC.make_type(dim=1, eq=False), 2*n) 

62 yield ("con", GM.make_type(argdim=num), n) 

63 yield ("con", AC.make_type(dim=1, eq=False), 1) 

64 else: 

65 assert relation == Constraint.GE 

66 assert p <= 1.0 

67 

68 if p == 1.0: 

69 yield ("con", AC.make_type(dim=n, eq=False), 1) 

70 yield ("con", AC.make_type(dim=1, eq=False), 1) 

71 elif p == float("-inf"): 

72 yield ("con", AC.make_type(dim=n, eq=False), 1) 

73 elif p > 0: 

74 yield ("var", RV.make_var_type(dim=n, bnd=0), 1) 

75 yield ("con", GM.make_type(argdim=den), n) 

76 yield ("con", AC.make_type(dim=1, eq=False), 1) 

77 else: # -inf < p < 0 

78 yield ("var", RV.make_var_type(dim=n, bnd=0), 1) 

79 yield ("con", GM.make_type(argdim=(num + den)), n) 

80 yield ("con", AC.make_type(dim=1, eq=False), 1) 

81 

82 @classmethod 

83 def convert(cls, con, options): 

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

85 from ..expressions import GeometricMean, RealVariable 

86 from ..modeling import Problem 

87 

88 norm = con.norm 

89 relation = con.relation 

90 rhs = con.rhs 

91 

92 p = norm.p 

93 m = len(norm.x) 

94 

95 P = Problem() 

96 

97 if relation == Constraint.LE: 

98 if p == 1: 

99 v = RealVariable('__v', m) 

100 P.add_constraint(norm.x[:] <= v) 

101 P.add_constraint(-norm.x[:] <= v) 

102 P.add_constraint(v.sum <= rhs) 

103 elif p == float('inf'): 

104 P.add_constraint(norm.x <= rhs) 

105 P.add_constraint(-norm.x <= rhs) 

106 else: 

107 x = RealVariable('__x', m) 

108 v = RealVariable('__v', m) 

109 amb = norm.pnum - norm.pden 

110 b = norm.pden 

111 oneamb = '|1|(' + str(amb) + ',1)' 

112 oneb = '|1|(' + str(b) + ',1)' 

113 for i in range(m): 

114 if amb > 0: 

115 if b == 1: 

116 vec = (v[i]) // (rhs * oneamb) 

117 else: 

118 vec = (v[i] * oneb) // (rhs * oneamb) 

119 else: 

120 # TODO: This shouldn't be possible? 

121 if b == 1: 

122 vec = v[i] 

123 else: 

124 vec = (v[i] * oneb) 

125 P.add_constraint(norm.x[i] <= x[i]) 

126 P.add_constraint(-norm.x[i] <= x[i]) 

127 P.add_constraint(x[i] <= GeometricMean(vec)) 

128 P.add_constraint(v.sum <= rhs) 

129 else: 

130 if p == 1: 

131 P.add_constraint(norm.x >= 0) 

132 P.add_constraint(norm.x.sum >= rhs) 

133 elif p == float('-inf'): 

134 P.add_constraint(norm.x >= rhs) 

135 elif p >= 0: 

136 v = RealVariable('__v', m) 

137 bma = -(norm.pnum - norm.pden) 

138 a = norm.pnum 

139 onebma = '|1|(' + str(bma) + ',1)' 

140 onea = '|1|(' + str(a) + ',1)' 

141 for i in range(m): 

142 if a == 1: 

143 vec = (norm.x[i]) // (rhs * onebma) 

144 else: 

145 vec = (norm.x[i] * onea) // (rhs * onebma) 

146 P.add_constraint(v[i] <= GeometricMean(vec)) 

147 P.add_constraint(rhs <= v.sum) 

148 else: 

149 v = RealVariable('__v', m) 

150 b = abs(norm.pden) 

151 a = abs(norm.pnum) 

152 oneb = '|1|(' + str(b) + ',1)' 

153 onea = '|1|(' + str(a) + ',1)' 

154 for i in range(m): 

155 if a == 1 and b == 1: 

156 vec = (norm.x[i]) // (v[i]) 

157 elif a > 1 and b == 1: 

158 vec = (norm.x[i] * onea) // (v[i]) 

159 elif a == 1 and b > 1: 

160 vec = (norm.x[i]) // (v[i] * oneb) 

161 else: 

162 vec = (norm.x[i] * onea) // (v[i] * oneb) 

163 P.add_constraint(rhs <= GeometricMean(vec)) 

164 P.add_constraint(v.sum <= rhs) 

165 

166 return P 

167 

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

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

170 

171 :param ~picos.expressions.Norm norm: 

172 Left hand side expression. 

173 :param str relation: 

174 Constraint relation symbol. 

175 :param ~picos.expressions.AffineExpression rhs: 

176 Right hand side expression. 

177 """ 

178 from ..expressions import AffineExpression, Norm 

179 

180 assert isinstance(norm, Norm) 

181 assert isinstance(rhs, AffineExpression) 

182 assert relation in self.LE + self.GE 

183 assert len(rhs) == 1 

184 

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

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

187 

188 if relation == self.LE: 

189 assert norm.convex, \ 

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

191 

192 if relation == self.GE: 

193 assert norm.concave, \ 

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

195 

196 self.norm = norm 

197 self.relation = relation 

198 self.rhs = rhs 

199 

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

201 

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

203 

204 def _subtype(self): 

205 return self.Subtype( 

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

207 

208 @classmethod 

209 def _cost(cls, subtype): 

210 return subtype.argdim + 1 

211 

212 def _expression_names(self): 

213 yield "norm" 

214 yield "rhs" 

215 

216 def _str(self): 

217 if self.relation == self.LE: 

218 str = glyphs.le(self.norm.string, self.rhs.string) 

219 else: 

220 str = glyphs.ge(self.norm.string, self.rhs.string) 

221 

222 if self.norm.p < 1: 

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

224 else: 

225 return str 

226 

227 def _get_slack(self): 

228 if self.relation == self.LE: 

229 return self.rhs.safe_value - self.norm.safe_value 

230 else: 

231 return self.norm.safe_value - self.rhs.safe_value 

232 

233 

234# -------------------------------------- 

235__all__ = api_end(_API_START, globals())