X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fcone.py;h=f3543a147ad8da3c5000015f3c53e837781180e5;hb=fbaecc56ec029d6f813d76e26bd8891a41416bf0;hp=821b95958e81f1e9a93dcf176f774a4c82e9317f;hpb=35b57bd6de2e1db25309b7cf20b39a5675bdb44f;p=sage.d.git diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py index 821b959..f3543a1 100644 --- a/mjo/cone/cone.py +++ b/mjo/cone/cone.py @@ -91,34 +91,41 @@ def _basically_the_same(K1, K2): -def _rho(K, K2=None): +def _restrict_to_space(K, W): r""" - Restrict ``K`` into its own span, or the span of another cone. + Restrict this cone a subspace of its ambient space. INPUT: - - ``K2`` -- another cone whose lattice has the same rank as this - cone. + - ``W`` -- The subspace into which this cone will be restricted. OUTPUT: - A new cone in a sublattice. + A new cone in a sublattice corresponding to ``W``. - EXAMPLES:: + EXAMPLES: + + When this cone is solid, restricting it into its own span should do + nothing:: sage: K = Cone([(1,)]) - sage: _rho(K) == K + sage: _restrict_to_space(K, K.span()) == K True + A single ray restricted into its own span gives the same output + regardless of the ambient space:: + sage: K2 = Cone([(1,0)]) - sage: _rho(K2).rays() + sage: K2_S = _restrict_to_space(K2, K2.span()).rays() + sage: K2_S N(1) in 1-d lattice N sage: K3 = Cone([(1,0,0)]) - sage: _rho(K3).rays() + sage: K3_S = _restrict_to_space(K3, K3.span()).rays() + sage: K3_S N(1) in 1-d lattice N - sage: _rho(K2) == _rho(K3) + sage: K2_S == K3_S True TESTS: @@ -127,8 +134,7 @@ def _rho(K, K2=None): sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8) - sage: K_S = _rho(K) - sage: K_S.is_solid() + sage: _restrict_to_space(K, K.span()).is_solid() True And the resulting cone should live in a space having the same @@ -136,22 +142,22 @@ def _rho(K, K2=None): sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8) - sage: K_S = _rho(K, K.dual() ) - sage: K_S.lattice_dim() == K.dual().dim() + sage: K_P = _restrict_to_space(K, K.dual().span()) + sage: K_P.lattice_dim() == K.dual().dim() True This function should not affect the dimension of a cone:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8) - sage: K.dim() == _rho(K).dim() + sage: K.dim() == _restrict_to_space(K,K.span()).dim() True Nor should it affect the lineality of a cone:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8) - sage: K.lineality() == _rho(K).lineality() + sage: K.lineality() == _restrict_to_space(K, K.span()).lineality() True No matter which space we restrict to, the lineality should not @@ -159,79 +165,71 @@ def _rho(K, K2=None): sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8) - sage: K.lineality() >= _rho(K).lineality() + sage: S = K.span(); P = K.dual().span() + sage: K.lineality() >= _restrict_to_space(K,S).lineality() True - sage: K.lineality() >= _rho(K, K.dual()).lineality() + sage: K.lineality() >= _restrict_to_space(K,P).lineality() True If we do this according to our paper, then the result is proper:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8) - sage: K_S = _rho(K) - sage: K_SP = _rho(K_S.dual()).dual() + sage: K_S = _restrict_to_space(K, K.span()) + sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual() sage: K_SP.is_proper() True - sage: K_SP = _rho(K_S, K_S.dual()) + sage: K_SP = _restrict_to_space(K_S, K_S.dual().span()) sage: K_SP.is_proper() True Test the proposition in our paper concerning the duals and restrictions. Generate a random cone, then create a subcone of - it. The operation of dual-taking should then commute with rho:: + it. The operation of dual-taking should then commute with + _restrict_to_space:: sage: set_random_seed() sage: J = random_cone(max_ambient_dim = 8) sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice()) - sage: K_W_star = _rho(K, J).dual() - sage: K_star_W = _rho(K.dual(), J) + sage: K_W_star = _restrict_to_space(K, J.span()).dual() + sage: K_star_W = _restrict_to_space(K.dual(), J.span()) sage: _basically_the_same(K_W_star, K_star_W) True """ - if K2 is None: - K2 = K - - # First we project K onto the span of K2. This will explode if the - # rank of ``K2.lattice()`` doesn't match ours. - span_K2 = Cone(K2.rays() + (-K2).rays(), lattice=K.lattice()) - K = K.intersection(span_K2) - - # Cheat a little to get the subspace span(K2). The paper uses the - # rays of K2 as a basis, but everything is invariant under linear - # isomorphism (i.e. a change of basis), and this is a little - # faster. - W = span_K2.linear_subspace() + # First we want to intersect ``K`` with ``W``. The easiest way to + # do this is via cone intersection, so we turn the subspace ``W`` + # into a cone. + W_cone = Cone(W.basis() + [-b for b in W.basis()], lattice=K.lattice()) + K = K.intersection(W_cone) # We've already intersected K with the span of K2, so every # generator of K should belong to W now. - W_rays = [ W.coordinate_vector(r) for r in K.rays() ] + K_W_rays = [ W.coordinate_vector(r) for r in K.rays() ] - L = ToricLattice(K2.dim()) - return Cone(W_rays, lattice=L) + L = ToricLattice(W.dimension()) + return Cone(K_W_rays, lattice=L) def discrete_complementarity_set(K): r""" - Compute the discrete complementarity set of this cone. - - The complementarity set of a cone is the set of all orthogonal pairs - `(x,s)` such that `x` is in the cone, and `s` is in its dual. The - discrete complementarity set is a subset of the complementarity set - where `x` and `s` are required to be generators of their respective - cones. + Compute a discrete complementarity set of this cone. - For polyhedral cones, the discrete complementarity set is always - finite. + A discrete complementarity set of `K` is the set of all orthogonal + pairs `(x,s)` such that `x \in G_{1}` and `s \in G_{2}` for some + generating sets `G_{1}` of `K` and `G_{2}` of its dual. Polyhedral + convex cones are input in terms of their generators, so "the" (this + particular) discrete complementarity set corresponds to ``G1 + == K.rays()`` and ``G2 == K.dual().rays()``. OUTPUT: A list of pairs `(x,s)` such that, * Both `x` and `s` are vectors (not rays). - * `x` is a generator of this cone. - * `s` is a generator of this cone's dual. + * `x` is one of ``K.rays()``. + * `s` is one of ``K.dual().rays()``. * `x` and `s` are orthogonal. REFERENCES: @@ -635,8 +633,8 @@ def lyapunov_rank(K): sage: set_random_seed() sage: K = random_cone(max_ambient_dim=8) sage: actual = lyapunov_rank(K) - sage: K_S = _rho(K) - sage: K_SP = _rho(K_S.dual()).dual() + sage: K_S = _restrict_to_space(K, K.span()) + sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual() sage: l = K.lineality() sage: c = K.codim() sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2 @@ -671,7 +669,7 @@ def lyapunov_rank(K): if m < n: # K is not solid, restrict to its span. - K = _rho(K) + K = _restrict_to_space(K, K.span()) # Non-solid reduction lemma. beta += (n - m)*n @@ -680,10 +678,59 @@ def lyapunov_rank(K): # K is not pointed, restrict to the span of its dual. Uses a # proposition from our paper, i.e. this is equivalent to K = # _rho(K.dual()).dual(). - K = _rho(K, K.dual()) + K = _restrict_to_space(K, K.dual().span()) # Non-pointed reduction lemma. beta += l * m beta += len(LL(K)) return beta + + + +def is_lyapunov_like(L,K): + r""" + Determine whether or not ``L`` is Lyapunov-like on ``K``. + + We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle + L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs + `\left\langle x,s \right\rangle` in the complementarity set of + ``K``. It is known [Orlitzky]_ that this property need only be + checked for generators of ``K`` and its dual. + + INPUT: + + - ``L`` -- A linear transformation or matrix. + + - ``K`` -- A polyhedral closed convex cone. + + OUTPUT: + + ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``, + and ``False`` otherwise. + + .. WARNING:: + + If this function returns ``True``, then ``L`` is Lyapunov-like + on ``K``. However, if ``False`` is returned, that could mean one + of two things. The first is that ``L`` is definitely not + Lyapunov-like on ``K``. The second is more of an "I don't know" + answer, returned (for example) if we cannot prove that an inner + product is zero. + + REFERENCES: + + .. [Orlitzky] M. Orlitzky. The Lyapunov rank of an + improper cone (preprint). + + EXAMPLES: + + todo. + + TESTS: + + todo. + + """ + return all([(L*x).inner_product(s) == 0 + for (x,s) in discrete_complementarity_set(K)])