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:`DetRootN`.""" 

21 

22import operator 

23from collections import namedtuple 

24 

25import cvxopt 

26import numpy 

27 

28from .. import glyphs 

29from ..apidoc import api_end, api_start 

30from ..constraints import DetRootNConstraint 

31from .data import convert_and_refine_arguments, convert_operands, cvx2np 

32from .exp_affine import AffineExpression, ComplexAffineExpression 

33from .expression import Expression, refine_operands, validate_prediction 

34 

35_API_START = api_start(globals()) 

36# ------------------------------- 

37 

38 

39class DetRootN(Expression): 

40 r"""The :math:`n`-th root of the determinant of an :math:`n\times n` matrix. 

41 

42 :Definition: 

43 

44 For an :math:`n \times n` positive semidefinite hermitian matrix :math:`X`, 

45 this is 

46 

47 .. math:: 

48 

49 \sqrt[n]{\det X}. 

50 

51 .. warning:: 

52 

53 When you pose a lower bound on the :math:`n`-th root of a determinant of 

54 the matrix :math:`X`, then PICOS enforces positive semidefiniteness 

55 :math:`X \succeq 0` through an auxiliary constraint during solution 

56 search. 

57 """ 

58 

59 # -------------------------------------------------------------------------- 

60 # Initialization and factory methods. 

61 # -------------------------------------------------------------------------- 

62 

63 @convert_and_refine_arguments("x") 

64 def __init__(self, x): 

65 """Construct a :class:`DetRootN`. 

66 

67 :param x: The matrix concerned. Must be hermitian by definition. 

68 :type x: ~picos.expressions.ComplexAffineExpression 

69 """ 

70 if not isinstance(x, ComplexAffineExpression): 

71 raise TypeError("Can only form the determinant of an affine " 

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

73 elif not x.square: 

74 raise TypeError("Can't take the determinant of non-square {0}." 

75 .format(x.string)) 

76 elif not x.hermitian: 

77 raise NotImplementedError("Taking the n-th root of the determinant " 

78 "of {0} is not supported as {0} is not necessarily hermitian." 

79 .format(x.string)) 

80 

81 self._x = x 

82 

83 Expression.__init__(self, "n-th Root of a Determinant", 

84 glyphs.power(glyphs.det(x.string), glyphs.div(1, x.shape[0]))) 

85 

86 # -------------------------------------------------------------------------- 

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

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

89 

90 def _get_refined(self): 

91 if self._x.constant: 

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

93 elif len(self._x) == 1: 

94 return self._x.renamed(self._symbStr) 

95 else: 

96 return self 

97 

98 Subtype = namedtuple("Subtype", ("diag", "complex")) 

99 

100 def _get_subtype(self): 

101 return self.Subtype(self.n, self._x.complex) 

102 

103 def _get_value(self): 

104 value = self._x._get_value() 

105 

106 det = numpy.linalg.det(cvx2np(value)) 

107 

108 if det < 0: 

109 raise ArithmeticError("Cannot evaluate {}: {} is negative." 

110 .format(self.string, glyphs.eq(glyphs.det(self.x.string), det))) 

111 

112 return cvxopt.matrix(det**(1.0 / self._x.shape[0])) 

113 

114 def _get_mutables(self): 

115 return self._x._get_mutables() 

116 

117 def _is_convex(self): 

118 return False 

119 

120 def _is_concave(self): 

121 return True 

122 

123 def _replace_mutables(self, mapping): 

124 return self.__class__(self._x._replace_mutables(mapping)) 

125 

126 def _freeze_mutables(self, freeze): 

127 return self.__class__(self._x._freeze_mutables(freeze)) 

128 

129 # -------------------------------------------------------------------------- 

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

131 # -------------------------------------------------------------------------- 

132 

133 @convert_operands(scalarRHS=True) 

134 @refine_operands() 

135 def __mul__(self, other): 

136 if isinstance(other, AffineExpression): 

137 if not other.constant: 

138 raise NotImplementedError("You may only multiply a nonconstant " 

139 "PICOS n-th root of a determinant with a constant term.") 

140 

141 root = DetRootN(other.value * self._x) 

142 root._typeStr = "Scaled " + root._typeStr 

143 root._symbStr = glyphs.clever_mul(self.string, other.string) 

144 return root 

145 else: 

146 return NotImplemented 

147 

148 @convert_operands(scalarRHS=True) 

149 @refine_operands() 

150 def __rmul__(self, other): 

151 if isinstance(other, AffineExpression): 

152 mean = self.__mul__(other) 

153 # NOTE: __mul__ always creates a fresh expression. 

154 mean._symbStr = glyphs.clever_mul(other.string, self.string) 

155 return mean 

156 else: 

157 return NotImplemented 

158 

159 # -------------------------------------------------------------------------- 

160 # Methods and properties that return modified copies. 

161 # -------------------------------------------------------------------------- 

162 

163 @property 

164 def x(self): 

165 """The matrix concerned.""" 

166 return self._x 

167 

168 # -------------------------------------------------------------------------- 

169 # Methods and properties that describe the expression. 

170 # -------------------------------------------------------------------------- 

171 

172 @property 

173 def n(self): 

174 """Diagonal length of :attr:`x`.""" 

175 return self._x.shape[0] 

176 

177 # -------------------------------------------------------------------------- 

178 # Constraint-creating operators, and _predict. 

179 # -------------------------------------------------------------------------- 

180 

181 @classmethod 

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

183 assert isinstance(subtype, cls.Subtype) 

184 

185 if relation == operator.__ge__: 

186 if issubclass(other.clstype, AffineExpression) \ 

187 and other.subtype.dim == 1: 

188 return DetRootNConstraint.make_type( 

189 diag=subtype.diag, complex=subtype.complex) 

190 

191 return NotImplemented 

192 

193 @convert_operands(scalarRHS=True) 

194 @validate_prediction 

195 @refine_operands() 

196 def __ge__(self, other): 

197 if isinstance(other, AffineExpression): 

198 return DetRootNConstraint(self, other) 

199 else: 

200 return NotImplemented 

201 

202 

203# -------------------------------------- 

204__all__ = api_end(_API_START, globals())