X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fcone.py;h=e9d0f1e643e3a0d40d8a754faf551ac317e36f37;hb=d27a3f72ac6b898534615de9111ace4082a4c55e;hp=ba5f51ea880ccdc2cc3344cb8b91022ff3e5b8cf;hpb=e041595c10751828f196db2cda86bd0f15a81191;p=sage.d.git diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py index ba5f51e..e9d0f1e 100644 --- a/mjo/cone/cone.py +++ b/mjo/cone/cone.py @@ -8,23 +8,12 @@ addsitedir(abspath('../../')) from sage.all import * -def drop_dependent(vs): - r""" - Return the largest linearly-independent subset of ``vs``. - """ - result = [] - m = matrix(vs).echelon_form() - for idx in range(0, m.nrows()): - if not m[idx].is_zero(): - result.append(m[idx]) - - return result - - def basically_the_same(K1,K2): r""" ``True`` if ``K1`` and ``K2`` are basically the same, and ``False`` - otherwise. + otherwise. This is intended as a lazy way to check whether or not + ``K1`` and ``K2`` are linearly isomorphic (i.e. ``A(K1) == K2`` for + some invertible linear transformation ``A``). """ if K1.lattice_dim() != K2.lattice_dim(): return False @@ -35,7 +24,7 @@ def basically_the_same(K1,K2): if K1.dim() != K2.dim(): return False - if lineality(K1) != lineality(K2): + if K1.lineality() != K2.lineality(): return False if K1.is_solid() != K2.is_solid(): @@ -65,7 +54,8 @@ def rho(K, K2=None): INPUT: - - ``K2`` -- another cone whose lattice has the same rank as this cone. + - ``K2`` -- another cone whose lattice has the same rank as this + cone. OUTPUT: @@ -118,7 +108,7 @@ def rho(K, K2=None): sage: set_random_seed() sage: K = random_cone(max_dim = 8) - sage: lineality(K) == lineality(rho(K)) + sage: K.lineality() == rho(K).lineality() True No matter which space we restrict to, the lineality should not @@ -126,9 +116,9 @@ def rho(K, K2=None): sage: set_random_seed() sage: K = random_cone(max_dim = 8) - sage: lineality(K) >= lineality(rho(K)) + sage: K.lineality() >= rho(K).lineality() True - sage: lineality(K) >= lineality(rho(K, K.dual())) + sage: K.lineality() >= rho(K, K.dual()).lineality() True If we do this according to our paper, then the result is proper:: @@ -136,11 +126,11 @@ def rho(K, K2=None): sage: set_random_seed() sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=False) sage: K_S = rho(K) - sage: P = rho(K_S.dual()).dual() - sage: P.is_proper() + sage: K_SP = rho(K_S.dual()).dual() + sage: K_SP.is_proper() True - sage: P = rho(K_S, K_S.dual()) - sage: P.is_proper() + sage: K_SP = rho(K_S, K_S.dual()) + sage: K_SP.is_proper() True :: @@ -148,11 +138,11 @@ def rho(K, K2=None): sage: set_random_seed() sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=False) sage: K_S = rho(K) - sage: P = rho(K_S.dual()).dual() - sage: P.is_proper() + sage: K_SP = rho(K_S.dual()).dual() + sage: K_SP.is_proper() True - sage: P = rho(K_S, K_S.dual()) - sage: P.is_proper() + sage: K_SP = rho(K_S, K_S.dual()) + sage: K_SP.is_proper() True :: @@ -160,11 +150,11 @@ def rho(K, K2=None): sage: set_random_seed() sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=True) sage: K_S = rho(K) - sage: P = rho(K_S.dual()).dual() - sage: P.is_proper() + sage: K_SP = rho(K_S.dual()).dual() + sage: K_SP.is_proper() True - sage: P = rho(K_S, K_S.dual()) - sage: P.is_proper() + sage: K_SP = rho(K_S, K_S.dual()) + sage: K_SP.is_proper() True :: @@ -172,47 +162,52 @@ def rho(K, K2=None): sage: set_random_seed() sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=True) sage: K_S = rho(K) - sage: P = rho(K_S.dual()).dual() - sage: P.is_proper() + sage: K_SP = rho(K_S.dual()).dual() + sage: K_SP.is_proper() True - sage: P = rho(K_S, K_S.dual()) - sage: P.is_proper() + sage: K_SP = rho(K_S, K_S.dual()) + sage: K_SP.is_proper() True - Test the proposition in our paper concerning the duals, where the - subspace `W` is the span of `K^{*}`:: + 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:: sage: set_random_seed() - sage: K = random_cone(max_dim = 8, solid=False, strictly_convex=False) - sage: K_W = rho(K, K.dual()) - sage: K_star_W_star = rho(K.dual()).dual() + sage: J = random_cone(max_dim = 8, solid=False, strictly_convex=False) + sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice()) + sage: K_W = rho(K, J) + sage: K_star_W_star = rho(K.dual(), J).dual() sage: basically_the_same(K_W, K_star_W_star) True :: sage: set_random_seed() - sage: K = random_cone(max_dim = 8, solid=True, strictly_convex=False) - sage: K_W = rho(K, K.dual()) - sage: K_star_W_star = rho(K.dual()).dual() + sage: J = random_cone(max_dim = 8, solid=True, strictly_convex=False) + sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice()) + sage: K_W = rho(K, J) + sage: K_star_W_star = rho(K.dual(), J).dual() sage: basically_the_same(K_W, K_star_W_star) True :: sage: set_random_seed() - sage: K = random_cone(max_dim = 8, solid=False, strictly_convex=True) - sage: K_W = rho(K, K.dual()) - sage: K_star_W_star = rho(K.dual()).dual() + sage: J = random_cone(max_dim = 8, solid=False, strictly_convex=True) + sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice()) + sage: K_W = rho(K, J) + sage: K_star_W_star = rho(K.dual(), J).dual() sage: basically_the_same(K_W, K_star_W_star) True :: sage: set_random_seed() - sage: K = random_cone(max_dim = 8, solid=True, strictly_convex=True) - sage: K_W = rho(K, K.dual()) - sage: K_star_W_star = rho(K.dual()).dual() + sage: J = random_cone(max_dim = 8, solid=True, strictly_convex=True) + sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice()) + sage: K_W = rho(K, J) + sage: K_star_W_star = rho(K.dual(), J).dual() sage: basically_the_same(K_W, K_star_W_star) True @@ -220,18 +215,16 @@ def rho(K, K2=None): if K2 is None: K2 = K - # First we project K onto the span of K2. This can be done with - # cones (i.e. without converting to vector spaces), but it's - # annoying to deal with lattice mismatches. + # 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) - V = K.lattice().vector_space() - - # Create the space W \times W^{\perp} isomorphic to V. - # First we get an orthogonal (but not normal) basis... - W_basis = drop_dependent(K2.rays()) - W = V.subspace_with_basis(W_basis) + # 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() # We've already intersected K with the span of K2, so every # generator of K should belong to W now. @@ -242,171 +235,6 @@ def rho(K, K2=None): -def lineality(K): - r""" - Compute the lineality of this cone. - - The lineality of a cone is the dimension of the largest linear - subspace contained in that cone. - - OUTPUT: - - A nonnegative integer; the dimension of the largest subspace - contained within this cone. - - REFERENCES: - - .. [Rockafellar] R.T. Rockafellar. Convex Analysis. Princeton - University Press, Princeton, 1970. - - EXAMPLES: - - The lineality of the nonnegative orthant is zero, since it clearly - contains no lines:: - - sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)]) - sage: lineality(K) - 0 - - However, if we add another ray so that the entire `x`-axis belongs - to the cone, then the resulting cone will have lineality one:: - - sage: K = Cone([(1,0,0), (-1,0,0), (0,1,0), (0,0,1)]) - sage: lineality(K) - 1 - - If our cone is all of `\mathbb{R}^{2}`, then its lineality is equal - to the dimension of the ambient space (i.e. two):: - - sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)]) - sage: lineality(K) - 2 - - Per the definition, the lineality of the trivial cone in a trivial - space is zero:: - - sage: K = Cone([], lattice=ToricLattice(0)) - sage: lineality(K) - 0 - - TESTS: - - The lineality of a cone should be an integer between zero and the - dimension of the ambient space, inclusive:: - - sage: set_random_seed() - sage: K = random_cone(max_dim = 8) - sage: l = lineality(K) - sage: l in ZZ - True - sage: (0 <= l) and (l <= K.lattice_dim()) - True - - A strictly convex cone should have lineality zero:: - - sage: set_random_seed() - sage: K = random_cone(max_dim = 8, strictly_convex = True) - sage: lineality(K) - 0 - - """ - return K.linear_subspace().dimension() - - -def codim(K): - r""" - Compute the codimension of this cone. - - The codimension of a cone is the dimension of the space of all - elements perpendicular to every element of the cone. In other words, - the codimension is the difference between the dimension of the - ambient space and the dimension of the cone itself. - - OUTPUT: - - A nonnegative integer representing the dimension of the space of all - elements perpendicular to this cone. - - .. seealso:: - - :meth:`dim`, :meth:`lattice_dim` - - EXAMPLES: - - The codimension of the nonnegative orthant is zero, since the span of - its generators equals the entire ambient space:: - - sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)]) - sage: codim(K) - 0 - - However, if we remove a ray so that the entire cone is contained - within the `x-y`-plane, then the resulting cone will have - codimension one, because the `z`-axis is perpendicular to every - element of the cone:: - - sage: K = Cone([(1,0,0), (0,1,0)]) - sage: codim(K) - 1 - - If our cone is all of `\mathbb{R}^{2}`, then its codimension is zero:: - - sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)]) - sage: codim(K) - 0 - - And if the cone is trivial in any space, then its codimension is - equal to the dimension of the ambient space:: - - sage: K = Cone([], lattice=ToricLattice(0)) - sage: K.lattice_dim() - 0 - sage: codim(K) - 0 - - sage: K = Cone([(0,)]) - sage: K.lattice_dim() - 1 - sage: codim(K) - 1 - - sage: K = Cone([(0,0)]) - sage: K.lattice_dim() - 2 - sage: codim(K) - 2 - - TESTS: - - The codimension of a cone should be an integer between zero and - the dimension of the ambient space, inclusive:: - - sage: set_random_seed() - sage: K = random_cone(max_dim = 8) - sage: c = codim(K) - sage: c in ZZ - True - sage: (0 <= c) and (c <= K.lattice_dim()) - True - - A solid cone should have codimension zero:: - - sage: set_random_seed() - sage: K = random_cone(max_dim = 8, solid = True) - sage: codim(K) - 0 - - The codimension of a cone is equal to the lineality of its dual:: - - sage: set_random_seed() - sage: K = random_cone(max_dim = 8, solid = True) - sage: codim(K) == lineality(K.dual()) - True - - """ - return (K.lattice_dim() - K.dim()) - - def discrete_complementarity_set(K): r""" Compute the discrete complementarity set of this cone. @@ -420,9 +248,7 @@ def discrete_complementarity_set(K): A list of pairs `(x,s)` such that, - * `x` is in this cone. * `x` is a generator of this cone. - * `s` is in this cone's dual. * `s` is a generator of this cone's dual. * `x` and `s` are orthogonal. @@ -720,7 +546,7 @@ def lyapunov_rank(K): sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z]) sage: lyapunov_rank(K) 19 - sage: K.lattice_dim()**2 - K.dim()*codim(K) + sage: K.lattice_dim()**2 - K.dim()*K.codim() 19 The Lyapunov rank should be additive on a product of proper cones @@ -827,10 +653,10 @@ def lyapunov_rank(K): sage: K = random_cone(max_dim=8) sage: actual = lyapunov_rank(K) sage: K_S = rho(K) - sage: P = rho(K_S.dual()).dual() - sage: l = lineality(K) - sage: c = codim(K) - sage: expected = lyapunov_rank(P) + K.dim()*(l + c) + c**2 + sage: K_SP = rho(K_S.dual()).dual() + sage: l = K.lineality() + sage: c = K.codim() + sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2 sage: actual == expected True @@ -863,25 +689,33 @@ def lyapunov_rank(K): sage: lyapunov_rank(K) == len(LL(K)) True + Test Theorem 3 in [Orlitzky/Gowda]_:: + + sage: set_random_seed() + sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True) + sage: L = ToricLattice(K.lattice_dim() + 1) + sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L) + sage: lyapunov_rank(K) >= K.lattice_dim() + True + """ - K_orig = K beta = 0 m = K.dim() n = K.lattice_dim() - l = lineality(K) + l = K.lineality() if m < n: - # K is not solid, project onto its span. + # K is not solid, restrict to its span. K = rho(K) # Lemma 2 beta += m*(n - m) + (n - m)**2 if l > 0: - # K is not pointed, project its dual onto its span. - # Uses a proposition from our paper, i.e. this is - # equivalent to K = rho(K.dual()).dual() + # 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()) # Lemma 3