1# ------------------------------------------------------------------------------ 

2# Copyright (C) 2019 Maximilian Stahlberg 

3# Based on the original picos.expressions module by Guillaume Sagnol. 

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"""Implements :class:`Entropy` and :class:`NegativeEntropy`.""" 

21 

22import math 

23import operator 

24from collections import namedtuple 

25 

26import cvxopt 

27import numpy 

28 

29from .. import glyphs 

30from ..apidoc import api_end, api_start 

31from ..caching import cached_selfinverse_unary_operator, cached_unary_operator 

32from ..constraints import KullbackLeiblerConstraint 

33from .data import convert_and_refine_arguments, convert_operands, cvx2np 

34from .exp_affine import AffineExpression 

35from .expression import Expression, refine_operands, validate_prediction 

36 

37_API_START = api_start(globals()) 

38# ------------------------------- 

39 

40 

41class Entropy(Expression): 

42 r"""Entropy or negative relative entropy of an affine expression. 

43 

44 Negative relative entropy is also known as the perspective of the logarithm. 

45 

46 :Definition: 

47 

48 Let :math:`x` be an :math:`n`-dimensional real affine expression. 

49 

50 1. If no additional expression :math:`y` is given, this is the entropy 

51 

52 .. math:: 

53 

54 -\sum_{i = 1}^n \operatorname{vec}(x)_i 

55 \log(\operatorname{vec}(x)_i). 

56 

57 2. If an additional affine expression :math:`y` of same shape as :math:`x` 

58 is given, this is the negative relative entropy (or logarithmic 

59 perspective) 

60 

61 .. math:: 

62 

63 -\sum_{i = 1}^n \operatorname{vec}(x)_i 

64 \log\left( 

65 \frac{\operatorname{vec}(x)_i}{\operatorname{vec}(y)_i} 

66 \right) 

67 &= -\sum_{i = 1}^n \operatorname{vec}(x)_i 

68 \left[ 

69 \log(\operatorname{vec}(x)_i) - \log(\operatorname{vec}(y)_i) 

70 \right] \\ 

71 &= \sum_{i = 1}^n \operatorname{vec}(x)_i 

72 \left[ 

73 \log(\operatorname{vec}(y)_i) - \log(\operatorname{vec}(x)_i 

74 \right] \\ 

75 &= \sum_{i = 1}^n \operatorname{vec}(x)_i 

76 \log\left( 

77 \frac{\operatorname{vec}(y)_i}{\operatorname{vec}(x)_i} 

78 \right). 

79 

80 .. warning:: 

81 

82 When you pose a lower bound on this expression, then PICOS enforces 

83 :math:`x \geq 0` through an auxiliary constraint during solution search. 

84 When an additional expression :math:`y` is given, PICOS enforces 

85 :math:`y \geq 0` as well. 

86 """ 

87 

88 # -------------------------------------------------------------------------- 

89 # Initialization and factory methods. 

90 # -------------------------------------------------------------------------- 

91 

92 @convert_and_refine_arguments("x", "y", allowNone=True) 

93 def __init__(self, x, y=None): 

94 """Construct an :class:`Entropy`. 

95 

96 :param x: The affine expression :math:`x`. 

97 :type x: ~picos.expressions.AffineExpression 

98 :param y: An additional affine expression :math:`y`. If necessary, PICOS 

99 will attempt to reshape or broadcast it to the shape of :math:`x`. 

100 :type y: ~picos.expressions.AffineExpression 

101 """ 

102 if not isinstance(x, AffineExpression): 

103 raise TypeError("Can only take the elementwise logarithm of a real " 

104 "affine expression, not of {}.".format(type(x).__name__)) 

105 

106 if y is not None: 

107 if not isinstance(y, AffineExpression): 

108 raise TypeError("The additional parameter y must be a real " 

109 "affine expression, not {}.".format(type(x).__name__)) 

110 elif x.shape != y.shape: 

111 y = y.reshaped_or_broadcasted(x.shape) 

112 

113 if y.is1: 

114 y = None 

115 

116 self._x = x 

117 self._y = y 

118 

119 if y is None: 

120 typeStr = "Entropy" 

121 if len(x) == 1: 

122 symbStr = glyphs.neg(glyphs.mul(x.string, glyphs.log(x.string))) 

123 else: 

124 symbStr = glyphs.neg( 

125 glyphs.dotp(x.string, glyphs.log(x.string))) 

126 else: 

127 typeStr = "Logarithmic Perspective" 

128 if len(x) == 1: 

129 symbStr = glyphs.mul( 

130 x.string, glyphs.log(glyphs.div(y.string, x.string))) 

131 else: 

132 symbStr = glyphs.dotp(x.string, 

133 glyphs.sub(glyphs.log(y.string), glyphs.log(x.string))) 

134 

135 Expression.__init__(self, typeStr, symbStr) 

136 

137 # -------------------------------------------------------------------------- 

138 # Abstract method implementations and method overridings, except _predict. 

139 # -------------------------------------------------------------------------- 

140 

141 def _get_refined(self): 

142 if self._x.constant and (self._y is None or self._y.constant): 

143 return AffineExpression.from_constant(self.value, 1, self._symbStr) 

144 else: 

145 return self 

146 

147 Subtype = namedtuple("Subtype", ("argdim", "y")) 

148 

149 def _get_subtype(self): 

150 return self.Subtype(len(self._x), self._y is not None) 

151 

152 def _get_value(self): 

153 x = numpy.ravel(cvx2np(self._x._get_value())) 

154 

155 if self._y is None: 

156 s = -numpy.dot(x, numpy.log(x)) 

157 else: 

158 y = numpy.ravel(cvx2np(self._y._get_value())) 

159 s = numpy.dot(x, numpy.log(y / x)) 

160 

161 return cvxopt.matrix(s) 

162 

163 @cached_unary_operator 

164 def _get_mutables(self): 

165 if self._y is None: 

166 return self._x._get_mutables() 

167 else: 

168 return self._x._get_mutables().union(self._y.mutables) 

169 

170 def _is_convex(self): 

171 return False 

172 

173 def _is_concave(self): 

174 return True 

175 

176 def _replace_mutables(self, mapping): 

177 return self.__class__(self._x._replace_mutables(mapping), 

178 None if self._y is None else self._y._replace_mutables(mapping)) 

179 

180 def _freeze_mutables(self, freeze): 

181 return self.__class__(self._x._freeze_mutables(freeze), 

182 None if self._y is None else self._y._freeze_mutables(freeze)) 

183 

184 # -------------------------------------------------------------------------- 

185 # Python special method implementations, except constraint-creating ones. 

186 # -------------------------------------------------------------------------- 

187 

188 @cached_selfinverse_unary_operator 

189 def __neg__(self): 

190 return NegativeEntropy(self._x, self._y) 

191 

192 @convert_operands(scalarRHS=True) 

193 @refine_operands() 

194 def __add__(self, other): 

195 if isinstance(other, AffineExpression): 

196 other = Entropy(other, other * math.e) 

197 

198 if isinstance(other, Entropy): 

199 if self._y is None and other._y is None: 

200 entropy = Entropy(self._x.vec // other._x.vec) 

201 elif self._y is not None and other._y is None: 

202 one = AffineExpression.from_constant(1.0, (other.n, 1)) 

203 entropy = Entropy( 

204 self._x.vec // other._x.vec, self._y.vec // one) 

205 elif self._y is None and other._y is not None: 

206 one = AffineExpression.from_constant(1.0, (self.n, 1)) 

207 entropy = Entropy( 

208 self._x.vec // other._x.vec, one // other._y.vec) 

209 else: 

210 entropy = Entropy( 

211 self._x.vec // other._x.vec, self._y.vec // other._y.vec) 

212 

213 entropy._symbStr = glyphs.clever_add(self.string, other.string) 

214 

215 return entropy 

216 else: 

217 return NotImplemented 

218 

219 @convert_operands(scalarRHS=True) 

220 @refine_operands() 

221 def __sub__(self, other): 

222 if isinstance(other, AffineExpression): 

223 other = NegativeEntropy(other, other / math.e) 

224 

225 if isinstance(other, NegativeEntropy): 

226 return self + (-other) 

227 else: 

228 return NotImplemented 

229 

230 @convert_operands(scalarRHS=True) 

231 @refine_operands() 

232 def __radd__(self, other): 

233 if isinstance(other, (AffineExpression, Entropy)): 

234 entropy = self.__add__(other) 

235 # NOTE: __add__ always creates a fresh expression. 

236 entropy._symbStr = glyphs.clever_add(other.string, self.string) 

237 return entropy 

238 else: 

239 return NotImplemented 

240 

241 # -------------------------------------------------------------------------- 

242 # Methods and properties that return expressions. 

243 # -------------------------------------------------------------------------- 

244 

245 @property 

246 def x(self): 

247 """The expression :math:`x`.""" 

248 return self._x 

249 

250 @property 

251 def y(self): 

252 """The additional expression :math:`y`, or :obj:`None`.""" 

253 return self._y 

254 

255 # -------------------------------------------------------------------------- 

256 # Methods and properties that describe the expression. 

257 # -------------------------------------------------------------------------- 

258 

259 @property 

260 def n(self): 

261 """Length of :attr:`x`.""" 

262 return len(self._x) 

263 

264 # -------------------------------------------------------------------------- 

265 # Constraint-creating operators, and _predict. 

266 # -------------------------------------------------------------------------- 

267 

268 @classmethod 

269 def _predict(cls, subtype, relation, other): 

270 assert isinstance(subtype, cls.Subtype) 

271 

272 if relation == operator.__ge__: 

273 if issubclass(other.clstype, AffineExpression) \ 

274 and other.subtype.dim == 1: 

275 return KullbackLeiblerConstraint.make_type( 

276 argdim=subtype.argdim) 

277 

278 return NotImplemented 

279 

280 @convert_operands(scalarRHS=True) 

281 @validate_prediction 

282 @refine_operands() 

283 def __ge__(self, other): 

284 if isinstance(other, AffineExpression): 

285 return KullbackLeiblerConstraint(-self, -other) 

286 else: 

287 return NotImplemented 

288 

289 

290class NegativeEntropy(Expression): 

291 r"""Negative or relative entropy of an affine expression. 

292 

293 Relative entropy is also known as the Kullback-Leibler divergence. 

294 

295 :Definition: 

296 

297 Let :math:`x` be an :math:`n`-dimensional real affine expression. 

298 

299 1. If no additional expression :math:`y` is given, this is the negative 

300 entropy 

301 

302 .. math:: 

303 

304 \sum_{i = 1}^n \operatorname{vec}(x)_i 

305 \log(\operatorname{vec}(x)_i). 

306 

307 2. If an additional affine expression :math:`y` of same shape as :math:`x` 

308 is given, this is the relative entropy (or Kullback-Leibler divergence) 

309 

310 .. math:: 

311 

312 \sum_{i = 1}^n \operatorname{vec}(x)_i 

313 \log\left( 

314 \frac{\operatorname{vec}(x)_i}{\operatorname{vec}(y)_i} 

315 \right) 

316 = \sum_{i = 1}^n \operatorname{vec}(x)_i 

317 \left[ 

318 \log(\operatorname{vec}(x)_i) - \log(\operatorname{vec}(y)_i) 

319 \right]. 

320 

321 .. warning:: 

322 

323 When you pose an upper bound on this expression, then PICOS enforces 

324 :math:`x \geq 0` through an auxiliary constraint during solution search. 

325 When an additional expression :math:`y` is given, PICOS enforces 

326 :math:`y \geq 0` as well. 

327 """ 

328 

329 # -------------------------------------------------------------------------- 

330 # Initialization and factory methods. 

331 # -------------------------------------------------------------------------- 

332 

333 @convert_and_refine_arguments("x", "y", allowNone=True) 

334 def __init__(self, x, y=None): 

335 """Construct a :class:`NegativeEntropy`. 

336 

337 :param x: The affine expression :math:`x`. 

338 :type x: ~picos.expressions.AffineExpression 

339 :param y: An additional affine expression :math:`y`. If necessary, PICOS 

340 will attempt to reshape or broadcast it to the shape of :math:`x`. 

341 :type y: ~picos.expressions.AffineExpression 

342 """ 

343 if not isinstance(x, AffineExpression): 

344 raise TypeError("Can only take the elementwise logarithm of a real " 

345 "affine expression, not of {}.".format(type(x).__name__)) 

346 

347 if y is not None: 

348 if not isinstance(y, AffineExpression): 

349 raise TypeError("The additional parameter y must be a real " 

350 "affine expression, not {}.".format(type(x).__name__)) 

351 elif x.shape != y.shape: 

352 y = y.reshaped_or_broadcasted(x.shape) 

353 

354 if y.is1: 

355 y = None 

356 

357 self._x = x 

358 self._y = y 

359 

360 if y is None: 

361 typeStr = "Negative Entropy" 

362 if len(x) == 1: 

363 symbStr = glyphs.mul(x.string, glyphs.log(x.string)) 

364 else: 

365 symbStr = glyphs.dotp(x.string, glyphs.log(x.string)) 

366 else: 

367 typeStr = "Relative Entropy" 

368 if len(x) == 1: 

369 symbStr = glyphs.mul( 

370 x.string, glyphs.log(glyphs.div(x.string, y.string))) 

371 else: 

372 symbStr = glyphs.dotp(x.string, 

373 glyphs.sub(glyphs.log(x.string), glyphs.log(y.string))) 

374 

375 Expression.__init__(self, typeStr, symbStr) 

376 

377 # -------------------------------------------------------------------------- 

378 # Abstract method implementations and method overridings, except _predict. 

379 # -------------------------------------------------------------------------- 

380 

381 def _get_refined(self): 

382 if self._x.constant and (self._y is None or self._y.constant): 

383 return AffineExpression.from_constant(self.value, 1, self._symbStr) 

384 else: 

385 return self 

386 

387 Subtype = namedtuple("Subtype", ("argdim", "y")) 

388 

389 def _get_subtype(self): 

390 return self.Subtype(len(self._x), self._y is not None) 

391 

392 def _get_value(self): 

393 x = numpy.ravel(cvx2np(self._x._get_value())) 

394 

395 if self._y is None: 

396 s = numpy.dot(x, numpy.log(x)) 

397 else: 

398 y = numpy.ravel(cvx2np(self._y._get_value())) 

399 s = numpy.dot(x, numpy.log(x / y)) 

400 

401 return cvxopt.matrix(s) 

402 

403 @cached_unary_operator 

404 def _get_mutables(self): 

405 if self._y is None: 

406 return self._x._get_mutables() 

407 else: 

408 return self._x._get_mutables().union(self._y.mutables) 

409 

410 def _is_convex(self): 

411 return True 

412 

413 def _is_concave(self): 

414 return False 

415 

416 def _replace_mutables(self, mapping): 

417 return self.__class__(self._x._replace_mutables(mapping), 

418 None if self._y is None else self._y._replace_mutables(mapping)) 

419 

420 def _freeze_mutables(self, freeze): 

421 return self.__class__(self._x._freeze_mutables(freeze), 

422 None if self._y is None else self._y._freeze_mutables(freeze)) 

423 

424 # -------------------------------------------------------------------------- 

425 # Python special method implementations, except constraint-creating ones. 

426 # -------------------------------------------------------------------------- 

427 

428 @cached_selfinverse_unary_operator 

429 def __neg__(self): 

430 return Entropy(self._x, self._y) 

431 

432 @convert_operands(scalarRHS=True) 

433 @refine_operands() 

434 def __add__(self, other): 

435 if isinstance(other, AffineExpression): 

436 other = NegativeEntropy(other, other / math.e) 

437 

438 if isinstance(other, NegativeEntropy): 

439 if self._y is None and other._y is None: 

440 negent = NegativeEntropy(self._x.vec // other._x.vec) 

441 elif self._y is not None and other._y is None: 

442 one = AffineExpression.from_constant(1.0, (other.n, 1)) 

443 negent = NegativeEntropy( 

444 self._x.vec // other._x.vec, self._y.vec // one) 

445 elif self._y is None and other._y is not None: 

446 one = AffineExpression.from_constant(1.0, (self.n, 1)) 

447 negent = NegativeEntropy( 

448 self._x.vec // other._x.vec, one // other._y.vec) 

449 else: 

450 negent = NegativeEntropy( 

451 self._x.vec // other._x.vec, self._y.vec // other._y.vec) 

452 

453 negent._symbStr = glyphs.clever_add(self.string, other.string) 

454 

455 return negent 

456 else: 

457 return NotImplemented 

458 

459 @convert_operands(scalarRHS=True) 

460 @refine_operands() 

461 def __radd__(self, other): 

462 if isinstance(other, (AffineExpression, NegativeEntropy)): 

463 negent = self.__add__(other) 

464 # NOTE: __add__ always creates a fresh expression. 

465 negent._symbStr = glyphs.clever_add(other.string, self.string) 

466 return negent 

467 else: 

468 return NotImplemented 

469 

470 @convert_operands(scalarRHS=True) 

471 @refine_operands() 

472 def __sub__(self, other): 

473 if isinstance(other, AffineExpression): 

474 other = Entropy(other, other * math.e) 

475 

476 if isinstance(other, Entropy): 

477 return self + (-other) 

478 else: 

479 return NotImplemented 

480 

481 # -------------------------------------------------------------------------- 

482 # Methods and properties that return expressions. 

483 # -------------------------------------------------------------------------- 

484 

485 @property 

486 def x(self): 

487 """The expression :math:`x`.""" 

488 return self._x 

489 

490 @property 

491 def y(self): 

492 """The additional expression :math:`y`, or :obj:`None`.""" 

493 return self._y 

494 

495 # -------------------------------------------------------------------------- 

496 # Methods and properties that describe the expression. 

497 # -------------------------------------------------------------------------- 

498 

499 @property 

500 def n(self): 

501 """Length of :attr:`x`.""" 

502 return len(self._x) 

503 

504 # -------------------------------------------------------------------------- 

505 # Constraint-creating operators, and _predict. 

506 # -------------------------------------------------------------------------- 

507 

508 @classmethod 

509 def _predict(cls, subtype, relation, other): 

510 assert isinstance(subtype, cls.Subtype) 

511 

512 if relation == operator.__le__: 

513 if issubclass(other.clstype, AffineExpression) \ 

514 and other.subtype.dim == 1: 

515 return KullbackLeiblerConstraint.make_type( 

516 argdim=subtype.argdim) 

517 

518 return NotImplemented 

519 

520 @convert_operands(scalarRHS=True) 

521 @validate_prediction 

522 @refine_operands() 

523 def __le__(self, other): 

524 if isinstance(other, AffineExpression): 

525 return KullbackLeiblerConstraint(self, other) 

526 else: 

527 return NotImplemented 

528 

529 

530# -------------------------------------- 

531__all__ = api_end(_API_START, globals())