]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
1 # Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we
2 # have to explicitly mangle our sitedir here so that "mjo.cone"
4 from os
.path
import abspath
5 from site
import addsitedir
6 addsitedir(abspath('../../'))
11 def random_cone(min_dim
=None, max_dim
=None, min_rays
=None, max_rays
=None):
13 Generate a random rational convex polyhedral cone.
15 Lower and upper bounds may be provided for both the dimension of the
16 ambient space and the number of generating rays of the cone. Any
17 parameters left unspecified will be chosen randomly.
21 - ``min_dim`` (default: random) -- The minimum dimension of the ambient
24 - ``max_dim`` (default: random) -- The maximum dimension of the ambient
27 - ``min_rays`` (default: random) -- The minimum number of generating rays
30 - ``max_rays`` (default: random) -- The maximum number of generating rays
35 A new, randomly generated cone.
39 It's hard to test the output of a random process, but we can at
40 least make sure that we get a cone back::
42 sage: from sage.geometry.cone import is_Cone
43 sage: K = random_cone()
44 sage: is_Cone(K) # long time
49 def random_min_max(l
,u
):
51 We need to handle four cases to prevent us from doing
52 something stupid like having an upper bound that's lower than
53 our lower bound. And we would need to repeat all of that logic
54 for the dimension/rays, so we consolidate it here.
56 if l
is None and u
is None:
57 # They're both random, just return a random nonnegative
59 return ZZ
.random_element().abs()
61 if l
is not None and u
is not None:
62 # Both were specified. Again, just make up a number and
63 # return it. If the user wants to give us u < l then he
64 # can have an exception.
65 return ZZ
.random_element(l
,u
)
67 if l
is not None and u
is None:
68 # In this case, we're generating the upper bound randomly
69 # GIVEN A LOWER BOUND. So we add a random nonnegative
70 # integer to the given lower bound.
71 u
= l
+ ZZ
.random_element().abs()
72 return ZZ
.random_element(l
,u
)
74 # Here we must be in the only remaining case, where we are
75 # given an upper bound but no lower bound. We might as well
77 return ZZ
.random_element(0,u
)
79 d
= random_min_max(min_dim
, max_dim
)
80 r
= random_min_max(min_rays
, max_rays
)
83 rays
= [L
.random_element() for i
in range(0,r
)]
85 # We pass the lattice in case there are no rays.
86 return Cone(rays
, lattice
=L
)
91 Compute the Lyapunov (or bilinearity) rank of this cone.
93 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
95 1. The dimension of the Lie algebra of the automorphism group of the
98 2. The dimension of the linear space of all Lyapunov-like
99 transformations on the cone.
103 A closed, convex polyhedral cone.
107 An integer representing the Lyapunov rank of the cone. If the
108 dimension of the ambient vector space is `n`, then the Lyapunov rank
109 will be between `1` and `n` inclusive; however a rank of `n-1` is
110 not possible (see the first reference).
114 In the references, the cones are always assumed to be proper. We
115 do not impose this restriction.
123 The codimension formula from the second reference is used. We find
124 all pairs `(x,s)` in the complementarity set of `K` such that `x`
125 and `s` are rays of our cone. It is known that these vectors are
126 sufficient to apply the codimension formula. Once we have all such
127 pairs, we "brute force" the codimension formula by finding all
128 linearly-independent `xs^{T}`.
132 1. M.S. Gowda and J. Tao. On the bilinearity rank of a proper cone
133 and Lyapunov-like transformations, Mathematical Programming, 147
136 2. G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
137 optimality constraints for the cone of positive polynomials,
138 Mathematical Programming, Series B, 129 (2011) 5-31.
142 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`::
144 sage: positives = Cone([(1,)])
145 sage: lyapunov_rank(positives)
147 sage: quadrant = Cone([(1,0), (0,1)])
148 sage: lyapunov_rank(quadrant)
150 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
151 sage: lyapunov_rank(octant)
154 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one::
156 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
157 sage: lyapunov_rank(L31)
160 Likewise for the `L^{3}_{\infty}` cone::
162 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
163 sage: lyapunov_rank(L3infty)
166 The Lyapunov rank should be additive on a product of cones::
168 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
169 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
170 sage: K = L31.cartesian_product(octant)
171 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
174 Two isomorphic cones should have the same Lyapunov rank. The cone
175 ``K`` in the following example is isomorphic to the nonnegative
176 octant in `\mathbb{R}^{3}`::
178 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
179 sage: lyapunov_rank(K)
182 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
185 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
186 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
191 The Lyapunov rank should be additive on a product of cones::
193 sage: K1 = random_cone(0,10,0,10)
194 sage: K2 = random_cone(0,10,0,10)
195 sage: K = K1.cartesian_product(K2)
196 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
199 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
202 sage: K = random_cone(0,10,0,10)
203 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
207 V
= K
.lattice().vector_space()
209 xs
= [V(x
) for x
in K
.rays()]
210 ss
= [V(s
) for s
in K
.dual().rays()]
212 # WARNING: This isn't really C(K), it only contains the pairs
213 # (x,s) in C(K) where x,s are extreme in their respective cones.
214 C_of_K
= [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
216 matrices
= [x
.column() * s
.row() for (x
,s
) in C_of_K
]
218 # Sage doesn't think matrices are vectors, so we have to convert
219 # our matrices to vectors explicitly before we can figure out how
220 # many are linearly-indepenedent.
222 # The space W has the same base ring as V, but dimension
223 # dim(V)^2. So it has the same dimension as the space of linear
224 # transformations on V. In other words, it's just the right size
225 # to create an isomorphism between it and our matrices.
226 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
230 Convert a matrix to a vector isomorphically.
234 vectors
= [phi(m
) for m
in matrices
]
236 return (W
.dimension() - W
.span(vectors
).rank())