diff --git a/docs/make_work.jl b/docs/make_work.jl index c6aaf3669a80..7a43d5645e15 100644 --- a/docs/make_work.jl +++ b/docs/make_work.jl @@ -16,10 +16,10 @@ import Documenter.Documents: Document, Page #= copied from Documenter (0.26.1) - Adds the loop to collect the AA/Nemo/Hecke doc files into the + Adds the loop to collect the AA/Nemo/Hecke doc files into the build tree of Oscar. -=# +=# function Selectors.runner(::Type{SetupBuildDirectory}, doc::Documents.Document) @info "SetupBuildDirectory: setting up build directory." @@ -219,7 +219,7 @@ makedocs(bib, "$(nemo)/series.md", "$(nemo)/puiseux.md"], ], - "Fields" => [ + "Fields" => [ "$(aa)/fields.md", "Rings/rational.md", "Number Fields" => [ @@ -234,7 +234,7 @@ makedocs(bib, "Local Fields" => [ "$(nemo)/padic.md", "$(nemo)/qadic.md"], - "$(nemo)/finitefield.md", + "$(nemo)/finitefield.md", ], "Groups" => [ "Groups/groups.md", "Groups/basics.md", @@ -274,6 +274,9 @@ makedocs(bib, "CommutativeAlgebra/ca_affine_algebras.md", "CommutativeAlgebra/ca_binomial_ideals.md", "CommutativeAlgebra/ca_invariant_theory.md"], + "Polyhedral Geometry" => ["PolyhedralGeometry/pg.md", + "PolyhedralGeometry/pg_polyhedra.md", + "PolyhedralGeometry/pg_linear_programs.md"], "References" => "references.md", "Index" => "manualindex.md", ] diff --git a/docs/oscar_references.bib b/docs/oscar_references.bib index 726cacdbd690..135edba4fa64 100644 --- a/docs/oscar_references.bib +++ b/docs/oscar_references.bib @@ -284,6 +284,27 @@ @article{EM19 URL = {https://doi.org/10.1112/jlms.12232}, } +@Book{JT13, + author = {Joswig, Michael and Theobald, Thorsten}, + title = {Polyhedral and Algebraic Methods in Computational Geometry}, + publisher = {Springer}, + year = {2013} +} + +@BOOK{Zie95, + title = {Lectures on polytopes}, + publisher = {Springer-Verlag}, + year = {1995}, + author = {Ziegler, Günter M.}, + volume = {152}, + pages = {x+370}, + series = {Graduate Texts in Mathematics}, + address = {New York}, + isbn = {0-387-94365-X}, + mrclass = {52Bxx}, + mrnumber = {MR1311028 (96a:52011)}, + mrreviewer = {Margaret M. Bayer} +} @inproceedings{BDLP19, TITLE = {{Computing integral bases via localization and Hensel lifting}}, diff --git a/docs/src/PolyhedralGeometry/pg.md b/docs/src/PolyhedralGeometry/pg.md new file mode 100644 index 000000000000..8ac9911620b2 --- /dev/null +++ b/docs/src/PolyhedralGeometry/pg.md @@ -0,0 +1,21 @@ +```@meta +CurrentModule = Oscar +``` + +```@setup oscar +using Oscar +``` + +```@contents +Pages = ["pg.md"] +``` + +# Introduction + +The polyhedral geometry part of OSCAR provides functionality for handling +- convex polytopes, unbounded polyhedra and cones +- linear programs + +General textbooks offering details on theory and algorithms include: +- [JT13](@cite) +- [Zie95](@cite) diff --git a/docs/src/PolyhedralGeometry/pg_cones.md b/docs/src/PolyhedralGeometry/pg_cones.md new file mode 100644 index 000000000000..4ec139fb7101 --- /dev/null +++ b/docs/src/PolyhedralGeometry/pg_cones.md @@ -0,0 +1,17 @@ +```@meta +CurrentModule = Oscar +``` + +```@setup oscar +using Oscar +``` + +```@contents +Pages = ["pg_cones.md"] +``` + +# Cones + +## Introduction + +## Constructions diff --git a/docs/src/PolyhedralGeometry/pg_linear_programs.md b/docs/src/PolyhedralGeometry/pg_linear_programs.md new file mode 100644 index 000000000000..dcefec606e93 --- /dev/null +++ b/docs/src/PolyhedralGeometry/pg_linear_programs.md @@ -0,0 +1,20 @@ +```@meta +CurrentModule = Oscar +``` + +```@setup oscar +using Oscar +``` + +```@contents +Pages = ["pg_linear_programs.md"] +``` + +# Linear programs + +## Introduction + +## Constructions + + + diff --git a/docs/src/PolyhedralGeometry/pg_polyhedra.md b/docs/src/PolyhedralGeometry/pg_polyhedra.md new file mode 100644 index 000000000000..59a73e64c4d6 --- /dev/null +++ b/docs/src/PolyhedralGeometry/pg_polyhedra.md @@ -0,0 +1,71 @@ +```@meta +CurrentModule = Oscar +``` + +```@setup oscar +using Oscar +``` + +```@contents +Pages = ["pg_polyhedra.md"] +``` + +# Polyhedra + +## Introduction + +A set $P \subseteq \mathbb{R}^n$ is called a (convex) polyhedron if it can be written as the intersection of finitely many closed affine halfspaces in $\mathbb{R}^n$. +That is, there exists a matrix $A$ and a vector $b$ such that +$$P = P(A,b) = \{ x \in \mathbb{R}^n \mid Ax \leq b\}.$$ +Writing $P$ as above is called an $H$-representation of $P$. + +When a polyhedron $P \subset \mathbb{R}^n$ is bounded, it is called a polytope and the fundamental theorem of polytopes states that it may be written as the convex hull of finitely many points. +That is $$P = \textrm{conv}(p_1,\ldots,p_N), p_i \in \mathbb{R}^n.$$ +Writing $P$ in this way is called a $V$-representation. +Polytopes are necessarily compact, i.e., they form convex bodies. + +Each polytope has a unique $V$-representation which is minimal with respect to inclusion (or cardinality). +Conversely, a polyhedron which is full-dimensional, has a unique minimal $H$-representation. +If the polyhedron is not full-dimensional, then there is no canonical choice of an $H$-representation. + + +## Construction + +Based on the definition of an $H$-representation, the constructor of `Polyhedron` can be used: + +```@docs +Polyhedron(A::Union{Oscar.MatElem,AbstractMatrix}, b) +``` + +The following defines a polytope, as a `Polyhedron`, via a $V$-representation by calling the following function: + +```@docs +convex_hull(::AnyVecOrMat) +``` + +### Computing convex hulls + +This is a standard triangle, defined via a (redundant) $V$-representation and its unique minimal $H$-representation: + +```@repl oscar +T = convex_hull([ 0 0 ; 1 0 ; 0 1; 0 1/2 ]) +facets_as_halfspace_matrix_pair(T) +``` + +## `Polyhedron` and `polymake`'s `Polytope` + +Many polyhedral computations are done through `polymake`. +This is visible in the structure `Polyhedron` via a pointer `pm_polytope` to the corresponding `polymake` object. +This allows to apply all functionality of `polymake` by calling a suitable `Polymake.jl` function on this pointer. + +To allow both `Oscar`'s and `polymake`'s functionality to be applicable to a polyhedron object, it can be converted back and forth: + +```@docs +Polyhedron(::Polymake.BigObject) +pm_polytope(::Polyhedron) +``` + +There are several design differences between `polymake` and `Oscar`. + +!!! warning + Polyhedra in `polymake` and `Polymake.jl` use homogeneous coordinates. The polyhedra in `Oscar` use affine coordinates. diff --git a/docs/src/index.md b/docs/src/index.md index 97e21daf2cbd..d727af7f3155 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,3 +1,6 @@ # Welcome to Oscar -Oscar is a new computer algebra system under development. +Oscar is a new computer algebra system. +While it is still under development it can already be used. + +Oscar features functions for groups, rings, and fields as well as linear and commutative algebra, number theory, algebraic and polyhedral geometry, and more. diff --git a/src/Polytopes/Cone/constructors.jl b/src/Polytopes/Cone/constructors.jl index 7d282a5c5b04..dfd3be8ef89c 100644 --- a/src/Polytopes/Cone/constructors.jl +++ b/src/Polytopes/Cone/constructors.jl @@ -5,10 +5,33 @@ ############################################################################### @doc Markdown.doc""" - Cone(Rays) + Cone(R [, L]) - A polyhedral cone, not necessarily pointed, defined by the positive hull - of the rays, Rays. +# Arguments +- `R::Matrix`: Rays generating the cone; encoded row-wise as representative vectors. +- `L::Matrix`: Generators of the Lineality space; encoded as row vectors. + +A polyhedral cone, not necessarily pointed, defined by the positive hull +of the rays, Rays. + +# Examples +To construct the positive orthant as a `Cone`, you can write: +```julia-repl +julia> R = [1 0; 0 1]; + +julia> PO = Cone(R) +A polyhedral cone in ambient dimension 2 +``` + +To obtain the upper half-space of the plane: +```julia-repl +julia> R = [0 1]; + +julia> L = [1 0]; + +julia> HS = Cone(R, L) +A polyhedral cone in ambient dimension 2 +``` """ struct Cone #a real polymake polyhedron pm_cone::Polymake.BigObject @@ -36,10 +59,21 @@ end """ - positive_hull(generators::Union{Oscar.MatElem,AbstractMatrix}) + positive_hull(generators) A polyhedral cone, not necessarily pointed, defined by the positive hull of the `generators`. Redundant rays are allowed in the generators. + +# Arguments +- `generators::Matrix`: Rays generating the cone; encoded row-wise as representative vectors. + +# Examples +```julia-repl +julia> R = [1 0; 0 1]; + +julia> PO = positive_hull(R) +A polyhedral cone in ambient dimension 2 +``` """ function positive_hull(generators::Union{Oscar.MatElem,AbstractMatrix}) # TODO: Filter out zero rows @@ -68,4 +102,3 @@ function Base.show(io::IO, C::Cone) end Polymake.visual(C::Cone; opts...) = Polymake.visual(pm_cone(C); opts...) - diff --git a/src/Polytopes/Cone/properties.jl b/src/Polytopes/Cone/properties.jl index 09b9acdf2054..f17245cf0f5e 100644 --- a/src/Polytopes/Cone/properties.jl +++ b/src/Polytopes/Cone/properties.jl @@ -20,9 +20,27 @@ Base.eltype(::Type{ConeRayIterator}) = Polymake.Vector{Polymake.Rational} Base.length(iter::ConeRayIterator) = nrays(iter.cone) """ - rays(C::Cone) + rays(C) -Returns the rays of a cone. +Return the rays of a cone. + +# Arguments +- `C::Cone`: A cone. + +# Examples +Here a cone is constructed from three rays. Calling `rays` reveals that one of these was redundant: +```julia-repl +julia> R = [1 0; 0 1; 0 2]; + +julia> PO = positive_hull(R); + +julia> collect(rays(PO)) +2-element Vector{Polymake.Vector{Polymake.Rational}}: + pm::Vector +1 0 + pm::Vector +0 1 +``` """ rays(C::Cone) = ConeRayIterator(C) @@ -37,30 +55,80 @@ rays(C::Cone) = ConeRayIterator(C) ############################################################################### """ - nrays(C::Cone) + nrays(C) + +Return the number of rays of the cone `C`. -Returns the number of rays of the cone `C`. +# Arguments +- `C::Cone`: A cone. + +# Examples +Here a cone is constructed from three rays. Calling `nrays` reveals that one of these was redundant: +```julia-repl +julia> R = [1 0; 0 1; 0 2]; + +julia> PO = positive_hull(R); + +julia> nrays(PO) +2 +``` """ nrays(C::Cone) = pm_cone(C).N_RAYS """ - dim(C::Cone) + dim(C) + +Return the dimension of a cone. + +# Arguments +- `C::Cone`: A cone. -Returns the dimension of a cone. +# Examples +The cone `C` in this example is 2-dimensional within a 3-dimensional ambient space. +```julia-repl +julia> C = Cone([1 0 0; 1 1 0; 0 1 0]); + +julia> dim(C) +2 +``` """ dim(C::Cone) = pm_cone(C).CONE_DIM """ - ambient_dim(C::Cone) + ambient_dim(C) + +Return the ambient dimension of a cone. + +# Arguments +- `C::Cone`: A cone. -Returns the ambient dimension of a cone. +# Examples +The cone `C` in this example is 2-dimensional within a 3-dimensional ambient space. +```julia-repl +julia> C = Cone([1 0 0; 1 1 0; 0 1 0]); + +julia> ambient_dim(C) +3 +``` """ ambient_dim(C::Cone) = pm_cone(C).CONE_AMBIENT_DIM """ - codim(C::Cone) + codim(C) Returns the codimension of a cone. + +# Arguments +- `C::Cone`: A cone. + +# Examples +The cone `C` in this example is 2-dimensional within a 3-dimensional ambient space. +```julia-repl +julia> C = Cone([1 0 0; 1 1 0; 0 1 0]); + +julia> codim(C) +1 +``` """ codim(C::Cone) = ambient_dim(C)-dim(C) @@ -68,16 +136,45 @@ codim(C::Cone) = ambient_dim(C)-dim(C) ## Boolean properties ############################################################################### """ - ispointed(C::Cone) + ispointed(C) + +Determine whether the cone is pointed, i.e. whether 0 is a face of the cone. + +# Arguments +- `C::Cone`: A cone. + +# Examples +A cone with lineality is not pointed, but a cone only consisting of a single ray is. +```julia-repl +julia> C = Cone([1 0], [0 1]); - Determine whether the cone is pointed, i.e. whether 0 is a face of the cone +julia> ispointed(C) +false + +julia> C = Cone([1 0]); + +julia> ispointed(C) +true +``` """ ispointed(C::Cone) = pm_cone(C).POINTED """ - isfulldimensional(C::Cone) + isfulldimensional(C) + +Determine whether the cone is full dimensional + +# Arguments +- `C::Cone`: A cone. + +# Examples +The cone `C` in this example is 2-dimensional within a 3-dimensional ambient space. +```julia-repl +julia> C = Cone([1 0 0; 1 1 0; 0 1 0]); - Determine whether the cone is full dimensional +julia> isfulldimensional(C) +false +``` """ isfulldimensional(C::Cone) = pm_cone(C).FULL_DIM @@ -86,9 +183,25 @@ isfulldimensional(C::Cone) = pm_cone(C).FULL_DIM ############################################################################### """ - rays_as_point_matrix(C::Cone) + rays_as_point_matrix(C) - Return the rays of a cone as rows in a matrix. +Return the rays of a cone as rows in a matrix. + +# Arguments +- `C::Cone`: A cone. + +# Examples +Here a cone is constructed from three rays. Calling `rays_as_point_matrix` reveals that one of these was redundant: +```julia-repl +julia> R = [1 0; 0 1; 0 2]; + +julia> PO = positive_hull(R); + +julia> rays_as_point_matrix(PO) +pm::Matrix +1 0 +0 1 +``` """ function rays_as_point_matrix(C::Cone) pm_cone(C).RAYS @@ -96,24 +209,70 @@ end """ - facets_as_point_matrix(C::Cone) + facets_as_point_matrix(C) + +Return the facets of a cone as rows of a matrix. -Returns the facets of a cone as rows of a matrix. +# Arguments +- `C::Cone`: A cone. + +# Examples +From this little example it is easy to see that the facets are displayed as their inside-pointing (w.r.t. the cone) normals. +```julia-repl +julia> R = [1 0; 1 1]; + +julia> C = positive_hull(R); + +julia> facets_as_point_matrix(C) +pm::Matrix +0 1 +1 -1 +``` """ facets_as_point_matrix(C::Cone) = pm_cone(C).FACETS """ - lineality_space(C::Cone) + lineality_space(C) + +Return a basis of the lineality space of a cone. + +# Arguments +- `C::Cone`: A cone. - Returns a basis of the lineality space of a cone. +# Examples +Three rays are used here to construct the upper half-plane. Actually, two of these rays point in opposite directions. +This gives us a 1-dimensional lineality. +```julia-repl +julia> UH = Cone([1 0; 0 1; -1 0]); + +julia> lineality_space(UH) +pm::Matrix +1 0 +``` """ lineality_space(C::Cone) = pm_cone(C).LINEALITY_SPACE """ - hilbert_basis(C::Cone) + hilbert_basis(C) + +Return the Hilbert basis of a pointed cone as the rows of a matrix. + +# Arguments +- `C::Cone`: A cone. + +# Examples +This (non-smooth) cone in the plane has a hilbert basis with three elements. +```julia-repl +julia> C = Cone([1 0; 1 2]) +A polyhedral cone in ambient dimension 2 - Returns the Hilbert basis of a pointed cone as the rows of a matrix. +julia> hilbert_basis(C) +pm::Matrix +1 0 +1 1 +1 2 +``` """ function hilbert_basis(C::Cone) if ispointed(C) diff --git a/src/Polytopes/Groups.jl b/src/Polytopes/Groups.jl index e1eb747c4d4c..8cd30f9a8ca5 100644 --- a/src/Polytopes/Groups.jl +++ b/src/Polytopes/Groups.jl @@ -47,9 +47,9 @@ function combinatorial_symmetries(P::Polyhedron) end """ - linear_symmetries(P::Polyhedron) + linear_symmetries(P::Polyhedron) - Get the group of linear symmetries on the vertices of a polyhedron. +Get the group of linear symmetries on the vertices of a polyhedron. """ function linear_symmetries(P::Polyhedron) if P.pm_polytope.BOUNDED diff --git a/src/Polytopes/LinearProgram.jl b/src/Polytopes/LinearProgram.jl index 988604e027ec..8a4f8704abd3 100644 --- a/src/Polytopes/LinearProgram.jl +++ b/src/Polytopes/LinearProgram.jl @@ -105,7 +105,7 @@ minimal_value(lp::LinearProgram) = lp.polymake_lp.MINIMAL_VALUE maximal_value(lp::LinearProgram) = lp.polymake_lp.MAXIMAL_VALUE """ - solve_lp(lp::LinearProgram) + solve_lp(lp::LinearProgram) Gives a pair `(m,v)` where the optimal value `m` of the objective function of lp is attained at `v` (if `m` exists). If the optimum @@ -129,7 +129,7 @@ end """ - primal_program(c, A, b) + primal_program(c, A, b) Constructs the primal linear program max{dot(c,x) | Ax<= b}. diff --git a/src/Polytopes/PolyhedralFan/constructors.jl b/src/Polytopes/PolyhedralFan/constructors.jl index c5422b85fb46..480da5122c6b 100644 --- a/src/Polytopes/PolyhedralFan/constructors.jl +++ b/src/Polytopes/PolyhedralFan/constructors.jl @@ -5,7 +5,7 @@ ############################################################################### @doc Markdown.doc""" - PolyhedralFan(Rays, Cones) + PolyhedralFan(Rays, Cones) A polyhedral fan formed from rays and cones made of these rays. The cones are given as an IncidenceMatrix, where the columns represent the rays and the rows @@ -35,7 +35,7 @@ function PolyhedralFan(Rays::Union{Oscar.MatElem,AbstractMatrix}, LS::Union{Osca end """ - pm_fan(PF::PolyhedralFan) + pm_fan(PF::PolyhedralFan) Get the underlying polymake object, which can be used via Polymake.jl. """ diff --git a/src/Polytopes/PolyhedralFan/properties.jl b/src/Polytopes/PolyhedralFan/properties.jl index 7de510f8459d..b388b33eeb68 100644 --- a/src/Polytopes/PolyhedralFan/properties.jl +++ b/src/Polytopes/PolyhedralFan/properties.jl @@ -20,7 +20,7 @@ Base.eltype(::Type{PolyhedralFanRayIterator}) = Polymake.Vector{Polymake.Rationa Base.length(iter::PolyhedralFanRayIterator) = nrays(iter.fan) """ - rays(PF::PolyhedralFan) + rays(PF::PolyhedralFan) Returns the rays of a polyhedral fan. """ @@ -28,7 +28,7 @@ rays(PF::PolyhedralFan) = PolyhedralFanRayIterator(PF) """ - maximal_cones(PF::PolyhedralFan) + maximal_cones(PF::PolyhedralFan) Returns the maximal cones of a polyhedral fan. """ @@ -36,7 +36,7 @@ Returns the maximal cones of a polyhedral fan. #TODO: should the documentation mention maximal_cones_as_incidence_matrix? # similarly for cone ray iterators and facet iterators? @doc Markdown.doc""" - maximal_cones(PF::PolyhedralFan, as = :cones) + maximal_cones(PF::PolyhedralFan, as = :cones) Returns an iterator over the maximal cones of the polyhedral fan `PF`. """ @@ -69,21 +69,21 @@ Base.length(iter::MaximalConeIterator) = nmaximal_cones(iter.PF) ############################################################################### """ - dim(PF::PolyhedralFan) + dim(PF::PolyhedralFan) Returns the dimension of a polyhedral fan. """ dim(PF::PolyhedralFan) = pm_fan(PF).FAN_DIM """ - nmaximal_cones(PF::PolyhedralFan) + nmaximal_cones(PF::PolyhedralFan) Returns the number of maximal cones in a polyhedral fan `PF`. """ nmaximal_cones(PF::PolyhedralFan) = pm_fan(PF).N_MAXIMAL_CONES """ - ambient_dim(PF::PolyhedralFan) + ambient_dim(PF::PolyhedralFan) Returns the ambient dimension of a polyhedral fan, which is the dimension of the embedding space. This is equal to the dimension of the fan if the fan is @@ -92,7 +92,7 @@ full-dimensional. ambient_dim(PF::PolyhedralFan) = pm_fan(PF).FAN_AMBIENT_DIM """ - nrays(PF::PolyhedralFan) + nrays(PF::PolyhedralFan) Returns the number of rays of a polyhedral fan. """ @@ -104,7 +104,7 @@ nrays(PF::PolyhedralFan) = pm_fan(PF).N_RAYS ############################################################################### """ - lineality_space(PF::PolyhedralFan) + lineality_space(PF::PolyhedralFan) Returns the lineality_space of a polyhedral fan. """ @@ -112,7 +112,7 @@ lineality_space(PF::PolyhedralFan) = pm_fan(PF).LINEALITY_SPACE """ - rays_as_point_matrix(PF::PolyhedralFan) + rays_as_point_matrix(PF::PolyhedralFan) Returns the rays of a polyhedral fan as rows of a matrix. """ @@ -120,7 +120,7 @@ rays_as_point_matrix(PF::PolyhedralFan) = pm_fan(PF).RAYS """ - maximal_cones_as_incidence_matrix(PF::PolyhedralFan) + maximal_cones_as_incidence_matrix(PF::PolyhedralFan) Returns the maximal cones of a polyhedral fan as an incidence matrix where the rows correspond to the maximal cones and the columns to the rays. @@ -133,21 +133,21 @@ end ## Boolean properties ############################################################################### """ - issmooth(PF::PolyhedralFan) + issmooth(PF::PolyhedralFan) Determine whether the fan is smooth. """ issmooth(PF::PolyhedralFan) = pm_fan(PF).SMOOTH_FAN """ - isregular(PF::PolyhedralFan) + isregular(PF::PolyhedralFan) Determine whether the fan is regular, i.e. the normal fan of a polytope. """ isregular(PF::PolyhedralFan) = pm_fan(PF).REGULAR """ - iscomplete(PF::PolyhedralFan) + iscomplete(PF::PolyhedralFan) Determine whether the fan is complete. """ diff --git a/src/Polytopes/PolyhedralFan/standard_constructions.jl b/src/Polytopes/PolyhedralFan/standard_constructions.jl index edb141730e5c..823476e8cf03 100644 --- a/src/Polytopes/PolyhedralFan/standard_constructions.jl +++ b/src/Polytopes/PolyhedralFan/standard_constructions.jl @@ -1,7 +1,7 @@ #TODO: inward/outward options? via polymake changes? """ - normal_fan(P::Polyhedron) + normal_fan(P::Polyhedron) Returns the normal fan of a polyhedron. """ @@ -12,7 +12,7 @@ function normal_fan(P::Polyhedron) end """ - face_fan(P::Polyhedron) + face_fan(P::Polyhedron) Returns the face fan of a polyhedron. """ diff --git a/src/Polytopes/Polyhedron/constructors.jl b/src/Polytopes/Polyhedron/constructors.jl index 62b1b2121fce..bdddecbbcbc8 100644 --- a/src/Polytopes/Polyhedron/constructors.jl +++ b/src/Polytopes/Polyhedron/constructors.jl @@ -4,35 +4,44 @@ ############################################################################### ############################################################################### +struct Polyhedron #a real polymake polyhedron + pm_polytope::Polymake.BigObject + boundedness::Symbol # Values: :unknown, :bounded, :unbounded +end + +@doc Markdown.doc""" + + Polyhedron(P::Polymake.BigObject) + +Construct a `Polyhedron` corresponding to a polymake `Polytope`. +""" +function Polyhedron(pm_polytope::Polymake.BigObject) + Polyhedron(pm_polytope, :unknown) +end @doc Markdown.doc""" - Polyhedron(A, b) + Polyhedron(A::Union{Oscar.MatElem,AbstractMatrix}, b) -The (metric) polyhedron defined by +The (convex) polyhedron defined by $$P(A,b) = \{ x | Ax ≤ b \}.$$ -see Def. 3.35 and Section 4.1. of Joswig, M. and Theobald, T. "Polyhedral and Algebraic Methods in Computational Geometry", Springer 2013. +see Def. 3.35 and Section 4.1. of [JT13] # Arguments - `A::Matrix`: Matrix corresponding to the linear coefficients of the inequalilites that describe P. - `b::Vector`: Vector corresponding to the constant term of the inequalilites that describe P. # Examples +The following lines define the square $[0,1]^2 \subset \mathbb{R}^2$: ```julia-repl -julia> A=[1 0; 0 1; -1 0 ; 0 -1]; -julia> b=[1,1,0,0]; -julia> Polyhedron([1 0;0 1;-1 0;0 -1],[1 1 0 0]) -A polyhedron of dimension 2 +julia> A = [1 0; 0 1; -1 0 ; 0 -1]; +julia> b = [1, 1, 0, 0]; +julia> Polyhedron(A,b) +A polyhedron in ambient dimension 2 ``` -""" struct Polyhedron #a real polymake polyhedron - pm_polytope::Polymake.BigObject - boundedness::Symbol # Values: :unknown, :bounded, :unbounded -end -function Polyhedron(pm_polytope::Polymake.BigObject) - Polyhedron(pm_polytope, :unknown) -end +""" function Polyhedron(A::Union{Oscar.MatElem,AbstractMatrix}, b) Polyhedron(Polymake.polytope.Polytope{Polymake.Rational}( INEQUALITIES = matrix_for_polymake(remove_zero_rows([b -A])), @@ -53,10 +62,44 @@ pm_polytope(P::Polyhedron) = P.pm_polytope @doc Markdown.doc""" convex_hull(V [, R [, L]]) -The polytope given as the convex hull of the columns of V. Optionally, rays (R) -and generators of the lineality space (L) can be given as well. +The polytope given as the convex hull of the rows of a set of points. + +see Def. 2.11 and Def. 3.1 of [JT13] + +# Arguments +- `V::Matrix`: Points whose convex hull is to be computed; encoded as row vectors. +- `R::Matrix`: Rays completing the set of points; encoded row-wise as representative vectors. +- `L::Matrix`: Generators of the Lineality space; encoded as row vectors. -see Def. 2.11 and Def. 3.1 of Joswig, M. and Theobald, T. "Polyhedral and Algebraic Methods in Computational Geometry", Springer 2013. +# Examples +The following lines define the square $[0,1]^2 \subset \mathbb{R}^2$: +```julia-repl +julia> Square = convex_hull([0 0; 0 1; 1 0; 1 1]) +A polyhedron in ambient dimension 2 +``` +To construct the positive orthant, rays have to be passed: +```julia-repl +julia> V = [0 0]; +julia> R = [1 0; 0 1]; +julia> PO = convex_hull(V, R) +A polyhedron in ambient dimension 2 +``` +The closed-upper half plane can be constructed by passing rays and a lineality space: +```julia-repl +julia> V = [0 0]; +julia> R = [0 1]; +julia> L = [1 0]; +julia> UH = convex_hull(V, R, L) +A polyhedron in ambient dimension 2 +``` +To obtain the x-axis in $\mathbb{R}^2$: +```julia-repl +julia> V = [0 0]; +julia> R = []; +julia> L = [1 0]; +julia> XA = convex_hull(V, R, L) +A polyhedron in ambient dimension 2 +``` """ function convex_hull(V::AnyVecOrMat; non_redundant::Bool=false) if !non_redundant @@ -119,4 +162,3 @@ function Base.show(io::IO, P::Polyhedron) end Polymake.visual(P::Polyhedron; opts...) = Polymake.visual(pm_polytope(P); opts...) - diff --git a/src/Polytopes/Polyhedron/properties.jl b/src/Polytopes/Polyhedron/properties.jl index e0b220e41205..0552d21a06f0 100644 --- a/src/Polytopes/Polyhedron/properties.jl +++ b/src/Polytopes/Polyhedron/properties.jl @@ -7,6 +7,14 @@ struct Polyhedra end struct Points end +@doc Markdown.doc""" + + Halfspaces + +Dummy type used for specifying the desired output format. +One halfspace `H(a,b)` is given by a vector `a` and a value `b` such that +$$H(a,b) = \{ x | ax ≤ b \}.$$ +""" struct Halfspaces end @@ -49,19 +57,67 @@ end Base.IteratorSize(::Type{PolyhedronFacePolyhedronIterator}) = Base.SizeUnknown() """ - faces(P::Polyhedron, face_dim::Int [, as::Type{T} = Polyhedron]) + faces(as, P, face_dim) -Returns the faces of `P` of dimension `face_dim` as an iterator over the type of object given -by `as`. Optional arguments for `as` include -* `Polyhedron` or `Polyhedra`: Returns for each face its realization as a polyhedron -""" -function faces(P::Polyhedron, face_dim::Int, as::Type{T} = Polyhedron) where {T} +Return the faces of `P` of dimension `face_dim` as an iterator over the type of object given +by `as`. +Optional arguments for `as` include +* `Polyhedron`/`Polyhedra`: Return the representation of a vertex as a point. + +# Arguments +- `as::Type{T}`: Element type of the returned iterator. +- `P::Polyhedron`: A polyhedron. +- `face_dim::Int`: Dimension of the desired faces. + +# Example +An `Array` containing the six sides of the 3-dimensional cube can be obtained via the following input: +```julia-repl +julia> F = faces(Polyhedron, cube(3), 2) +Oscar.PolyhedronFacePolyhedronIterator(A polyhedron in ambient dimension 3, 2) + +julia> collect(F) +6-element Array{Any,1}: + A polyhedron in ambient dimension 3 + A polyhedron in ambient dimension 3 + A polyhedron in ambient dimension 3 + A polyhedron in ambient dimension 3 + A polyhedron in ambient dimension 3 + A polyhedron in ambient dimension 3 +``` +""" +function faces(as::Type{T}, P::Polyhedron, face_dim::Int) where {T} if as == Polyhedron || as == Polyhedra PolyhedronFacePolyhedronIterator(P,face_dim-size(lineality_space(P),1)) else throw(ArgumentError("Unsupported `as` argument :" * string(as))) end end +""" + faces(P, face_dim) + +Return the faces of `P` of dimension `face_dim` as an iterator over `Polyhedron` objects. + +# Arguments +- `P::Polyhedron`: A polyhedron. +- `face_dim::Int`: Dimension of the desired faces. + +# Example +An `Array` containing the six sides of the 3-dimensional cube can be obtained via the following input: +```julia-repl +julia> F = faces(cube(3),2) +Oscar.PolyhedronFacePolyhedronIterator(A polyhedron in ambient dimension 3, 2) + +julia> collect(F) +6-element Array{Any,1}: + A polyhedron in ambient dimension 3 + A polyhedron in ambient dimension 3 + A polyhedron in ambient dimension 3 + A polyhedron in ambient dimension 3 + A polyhedron in ambient dimension 3 + A polyhedron in ambient dimension 3 +``` +""" +faces(P::Polyhedron, face_dim::Int) = faces(Polyhedron, P, face_dim) struct VertexPointIterator p::Polyhedron @@ -85,15 +141,38 @@ Base.eltype(::Type{VertexPointIterator}) = Polymake.Vector{Polymake.Rational} Base.length(iter::VertexPointIterator) = nvertices(iter.p) """ - vertices(P::Polyhedron, [,as::Type{T} = Points]) + vertices(as, P) -Returns an iterator over the vertices of a polyhedron `P` in the format defined by `as`. +Return an iterator over the vertices of a polyhedron `P` in the format defined by `as`. Optional arguments for `as` include * `Points`: Returns the representation of a vertex as a point. See also `vertices_as_point_matrix`. -""" -function vertices(P::Polyhedron, as::Type{T} = Points) where {T} + +# Arguments +- `as::Type{T}`: Element type of the returned iterator. +- `P::Polyhedron`: A polyhedron. + +# Examples +The following code computes the vertices of the Minkowski sum of a triangle and a square: +```julia-repl +julia> P = simplex(2) + cube(2); + +julia> collect(vertices(Points, P)) +5-element Array{Polymake.Vector{Polymake.Rational},1}: +pm::Vector +-1 -1 +pm::Vector +2 -1 +pm::Vector +2 1 +pm::Vector +-1 2 +pm::Vector +1 2 +``` +""" +function vertices(as::Type{T}, P::Polyhedron) where {T} if as == Points VertexPointIterator(P) else @@ -101,11 +180,58 @@ function vertices(P::Polyhedron, as::Type{T} = Points) where {T} end end +""" + vertices(P) + +Return an iterator over the vertices of a polyhedron `P` as points. + +See also `vertices_as_point_matrix`. +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +The following code computes the vertices of the Minkowski sum of a triangle and a square: +```julia-repl +julia> P = simplex(2) + cube(2); + +julia> collect(vertices(P)) +5-element Array{Polymake.Vector{Polymake.Rational},1}: +pm::Vector +-1 -1 +pm::Vector +2 -1 +pm::Vector +2 1 +pm::Vector +-1 2 +pm::Vector +1 2 +``` """ - `vertices_as_point_matrix(P::Polyhedron)` +vertices(P::Polyhedron) = vertices(Points, P) -Returns a matrix whose rows are the vertices of `P`. +""" + `vertices_as_point_matrix(P)` + +Return a matrix whose rows are the vertices of `P`. + +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +The following code computes the vertices of the Minkowski sum of a triangle and a square: +```julia-repl +julia> P = simplex(2) + cube(2); + +julia> vertices_as_point_matrix(P) +pm::Matrix +-1 -1 +2 -1 +2 1 +-1 2 +1 2 +``` """ function vertices_as_point_matrix(P::Polyhedron) decompose_vdata(pm_polytope(P).VERTICES).vertices @@ -134,31 +260,71 @@ Base.eltype(::Type{PolyhedronRayIterator}) = Polymake.Vector{Polymake.Rational} Base.length(iter::PolyhedronRayIterator) = nrays(iter.p) """ - nrays(P::Polyhedron) + nrays(P) + +Return the number of rays of a `Polyhedron`. + +# Arguments +- `P::Polyhedron`: A polyhedron. -Returns the number of rays of `P`. +# Examples +Reflecting the input, the upper half-plane indeed has one ray. +```julia-repl +julia> UH = convex_hull([0 0],[0 1],[1 0]); + +julia> nrays(UH) +1 +``` """ nrays(P::Polyhedron) = length(pm_polytope(P).FAR_FACE) """ - nvertices(P::Polyhedron) + nvertices(P) + +Return the number of vertices of the polyhedron `P`. + +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +The 3-cube's number of vertices can be obtained with this input: +```julia-repl +julia> C = cube(3); -Returns the number of vertices of `P`. +julia> nvertices(C) +8 +``` """ nvertices(P::Polyhedron) = pm_polytope(P).N_VERTICES - nrays(P) """ - rays(P::Polyhedron, [,as::Type{T} = Points]) + rays(as, P) -Returns minimal set of generators of the cone of unbounded directions of a polyhedron `P` (i.e. its rays) +Return minimal set of generators of the cone of unbounded directions of a polyhedron `P` (i.e. its rays) in the format defined by `as`. Optional arguments for `as` include * `Points`: Returns a vector representation of a ray. See also `rays_as_point_matrix`. -""" -function rays(P::Polyhedron, as::Type{T} = Points) where {T} + +# Arguments +- `as::Type{T}`: Element type of the returned iterator. +- `P::Polyhedron`: A polyhedron. + +# Examples +We can verify that the positive orthant of the plane is generated by the two rays in positive unit direction: +```julia-repl +julia> PO = convex_hull([0 0], [1 0; 0 1]); +julia> collect(rays(Points, PO)) +2-element Array{Polymake.Vector{Polymake.Rational},1}: + pm::Vector +1 0 + pm::Vector +0 1 +``` +""" +function rays(as::Type{T}, P::Polyhedron) where {T} if as == Points PolyhedronRayIterator(P) else @@ -166,6 +332,32 @@ function rays(P::Polyhedron, as::Type{T} = Points) where {T} end end +""" + rays(P) + + +Return minimal set of generators of the cone of unbounded directions of a polyhedron `P` (i.e. its rays) + as points. + +See also `rays_as_point_matrix`. + +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +We can verify that the positive orthant of the plane is generated by the two rays in positive unit direction: +```julia-repl +julia> PO = convex_hull([0 0], [1 0; 0 1]); +julia> collect(rays(PO)) +2-element Array{Polymake.Vector{Polymake.Rational},1}: + pm::Vector +1 0 + pm::Vector +0 1 +``` +""" +rays(P::Polyhedron) = rays(Points,P) + function rays_as_point_matrix(P::Polyhedron) decompose_vdata(pm_polytope(P).VERTICES).rays end @@ -206,23 +398,67 @@ Base.length(iter::PolyhedronFacetPolyhedronIterator) = nfacets(iter.p) # Base.eltype(::Type{PolyhedronFacetPolyhedronIterator}) = Polyhedron """ - nfacets(P::Polyhedron) + nfacets(P) + +Return the number of facets of the polyhedron `P`. + +# Arguments +- `P::Polyhedron`: A polyhedron. -Returns the number of facets of the polyhedron `P`. +# Examples +The number of facets of the 5-dimensional cross polytope can be retrieved via the following line: +```julia-repl +julia> nfacets(cross(5)) +32 +``` """ nfacets(P::Polyhedron) = pm_polytope(P).N_FACETS @doc Markdown.doc""" - facets(P::Polyhedron, as = :halfspaces) + facets(as::Type{T}, P::Polyhedron) -Returns the facets of the polyhedron `P` in the format defined by `as`. +Return the facets of the polyhedron `P` in the format defined by `as`. The allowed values for `as` are * `Halfspaces`: Returns for each facet the tuple `(A, b)` describing the halfspace `dot(A,x) ≤ b`. * `Polyhedron` or `Polyhedra`: Returns for each facet its realization as a polyhedron See also `facets_as_halfspace_matrix_pair`. -""" -function facets(P::Polyhedron, as::Type{T} = Halfspaces) where {T} + +# Arguments +- `as::Type{T}`: Element type of the returned iterator. +- `P::Polyhedron`: A polyhedron. + +# Examples +We can retrieve the six facets of the 3-dimensional cube this way: +```julia-repl +julia> C = cube(3); + +julia> collect(facets(Polyhedron, C)) +6-element Array{Any,1}: +A polyhedron in ambient dimension 3 +A polyhedron in ambient dimension 3 +A polyhedron in ambient dimension 3 +A polyhedron in ambient dimension 3 +A polyhedron in ambient dimension 3 +A polyhedron in ambient dimension 3 + +julia> collect(facets(Halfspaces, C)) +6-element Array{Tuple{Polymake.Vector{Polymake.Rational},Polymake.Rational},1}: +(pm::Vector +-1 0 0, 1) +(pm::Vector +1 0 0, 1) +(pm::Vector +0 -1 0, 1) +(pm::Vector +0 1 0, 1) +(pm::Vector +0 0 -1, 1) +(pm::Vector +0 0 1, 1) +``` +""" +function facets(as::Type{T}, P::Polyhedron) where {T} if as == Halfspaces PolyhedronFacetHalfspaceIterator(P) elseif as == Polyhedra || as == Polyhedron @@ -231,15 +467,66 @@ function facets(P::Polyhedron, as::Type{T} = Halfspaces) where {T} throw(ArgumentError("Unsupported `as` argument :" * string(as))) end end +@doc Markdown.doc""" + facets(P::Polyhedron) + +Return the facets of the polyhedron `P` as halfspaces. + +See also `facets_as_halfspace_matrix_pair`. + +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +We can retrieve the six facets of the 3-dimensional cube this way: +```julia-repl +julia> C = cube(3); + +julia> collect(facets(C)) +6-element Array{Tuple{Polymake.Vector{Polymake.Rational},Polymake.Rational},1}: +(pm::Vector +-1 0 0, 1) +(pm::Vector +1 0 0, 1) +(pm::Vector +0 -1 0, 1) +(pm::Vector +0 1 0, 1) +(pm::Vector +0 0 -1, 1) +(pm::Vector +0 0 1, 1) +``` +""" +facets(P::Polyhedron) = facets(Halfspaces, P) #TODO: how do underscores work in markdown? @doc Markdown.doc""" - `facets_as_halfspace_matrix_pair(P::Polyhedron)` + facets_as_halfspace_matrix_pair(P::Polyhedron) -Returns `(A,b)` such that $P=P(A,b)$ where +Return `(A,b)` such that $P=P(A,b)$ where $$P(A,b) = \{ x | Ax ≤ b \}.$$ + +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +```julia-repl +julia> C = cube(3); + +julia> facets_as_halfspace_matrix_pair(C) +(A = pm::SparseMatrix +(3) (0 -1) +(3) (0 1) +(3) (1 -1) +(3) (1 1) +(3) (2 -1) +(3) (2 1) +, b = pm::SparseVector +1 1 1 1 1 1) +``` """ function facets_as_halfspace_matrix_pair(P::Polyhedron) return decompose_hdata(pm_polytope(P).FACETS) @@ -256,31 +543,90 @@ end ## Scalar properties ############################################################################### """ - volume(P::Polyhedron) + volume(P) + +Return the (Euclidean) volume of a polyhedron. -Returns the (Euclidean) volume of a polyhedron. +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +```julia-repl +julia> C = cube(2); + +julia> volume(C) +4 +``` """ volume(P::Polyhedron) = (pm_polytope(P)).VOLUME """ - normalized_volume(P::Polyhedron) + normalized_volume(P) + +Return the (normalized) volume of a polyhedron. + +# Arguments +- `P::Polyhedron`: A polyhedron. -Returns the (normalized) volume of a polyhedron. +# Examples +```julia-repl +julia> C = cube(2); + +julia> normalized_volume(C) +8 +``` """ normalized_volume(P::Polyhedron) = factorial(dim(P))*(pm_polytope(P)).VOLUME """ - dim(P::Polyhedron) + dim(P) + +Return the dimension of a polyhedron. + +# Arguments +- `P::Polyhedron`: A polyhedron. -Returns the dimension of a polyhedron. +# Examples +```julia-repl +julia> V = [1 2 3; 1 3 2; 2 1 3; 2 3 1; 3 1 2; 3 2 1]; + +julia> P = convex_hull(V) +A polyhedron in ambient dimension 3 + +julia> dim(P) +2 +``` """ dim(P::Polyhedron) = Polymake.polytope.dim(pm_polytope(P)) """ - lattice_points(P::Polyhedron) + lattice_points(P) + +Return the integer points contained in a bounded polyhedron. -Returns the integer points contained in a bounded polyhedron. +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +```julia-repl +julia> S = 2 * simplex(2); + +julia> lattice_points(S) +6-element Array{Polymake.VectorAllocated{Polymake.Integer},1}: + pm::Vector +0 0 + pm::Vector +0 1 + pm::Vector +0 2 + pm::Vector +1 0 + pm::Vector +1 1 + pm::Vector +2 0 +``` """ function lattice_points(P::Polyhedron) if pm_polytope(P).BOUNDED @@ -294,16 +640,44 @@ end # scalable way to construct these iterators for so many functions """ - ambient_dim(P::Polyhedron) + ambient_dim(P) + +Return the ambient dimension of a polyhedron. -Returns the ambient dimension of a polyhedron. +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +```julia-repl +julia> V = [1 2 3; 1 3 2; 2 1 3; 2 3 1; 3 1 2; 3 2 1]; + +julia> P = convex_hull(V) +A polyhedron in ambient dimension 3 + +julia> ambient_dim(P) +3 +``` """ ambient_dim(P::Polyhedron) = Polymake.polytope.ambient_dim(pm_polytope(P)) """ - codim(P::Polyhedron) + codim(P) + +Return the codimension of a polyhedron. + +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +```julia-repl +julia> V = [1 2 3; 1 3 2; 2 1 3; 2 3 1; 3 1 2; 3 2 1]; + +julia> P = convex_hull(V) +A polyhedron in ambient dimension 3 -Returns the codimension of a polyhedron. +julia> codim(P) +1 +``` """ codim(P::Polyhedron) = ambient_dim(P)-dim(P) @@ -316,17 +690,56 @@ codim(P::Polyhedron) = ambient_dim(P)-dim(P) # Taylor: lineality space generators always look like [0, v] so # v is a natural output. """ - lineality_space(`P`::Polyhedron) + lineality_space(P) -Returns a matrix whose row span is the lineality space of a `P`. +Return a matrix whose row span is the lineality space of a polyhedron. + +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +```julia-repl +julia> UH = convex_hull([0 0],[0 1; 1 0; -1 0]); + +julia> lineality_space(UH) +pm::Matrix +1 0 +``` """ lineality_space(P::Polyhedron) = dehomogenize(pm_polytope(P).LINEALITY_SPACE) """ - recession_cone(P::Polyhedron) + recession_cone(P) + +Return the recession cone of `P`. + +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +```julia-repl +julia> P = Polyhedron([1 -2; -1 1; -1 0; 0 -1],[2,1,1,1]); -Returns the recession cone of `P`. +julia> collect(vertices(P)) +3-element Array{Polymake.Vector{Polymake.Rational},1}: + pm::Vector +0 -1 + pm::Vector +-1 0 + pm::Vector +-1 -1 + +julia> recession_cone(P) +A polyhedral cone in ambient dimension 2 + +julia> collect(rays(recession_cone(P))) +2-element Array{Polymake.Vector{Polymake.Rational},1}: + pm::Vector +1 1/2 + pm::Vector +1 1 +``` """ recession_cone(P::Polyhedron) = Cone(Polymake.polytope.recession_cone(pm_polytope(P))) @@ -334,48 +747,114 @@ recession_cone(P::Polyhedron) = Cone(Polymake.polytope.recession_cone(pm_polytop ## Boolean properties ############################################################################### """ - isfeasible(P::Polyhedron) + isfeasible(P) + +Check whether a polyhedron is feasible, i.e. non-empty. + +# Arguments +- `P::Polyhedron`: A polyhedron. - Check whether a polyhedron is feasible, i.e. non-empty. +# Examples +```julia-repl +julia> P = Polyhedron([1 -1; -1 1; -1 0; 0 -1],[-1,-1,1,1]) +A polyhedron in ambient dimension 2 + +julia> isfeasible(P) +false +``` """ isfeasible(P::Polyhedron) = pm_polytope(P).FEASIBLE """ - issmooth(P::Polyhedron) + issmooth(P) + +Check whether a polyhedron is smooth. + +# Arguments +- `P::Polyhedron`: A polyhedron. - Check whether a polyhedron is smooth. +# Examples +```julia-repl +julia> C = cube(8); + +julia> issmooth(C) +true +``` """ issmooth(P::Polyhedron) = pm_polytope(P).SMOOTH """ - isnormal(P::Polyhedron) + isnormal(P) + +Check whether a polyhedron is normal. + +# Arguments +- `P::Polyhedron`: A polyhedron. - Check whether a polyhedron is normal. +# Examples +TODO """ isnormal(P::Polyhedron) = pm_polytope(P).NORMAL """ - isbounded(P::Polyhedron) + isbounded(P) + +Check whether a polyhedron is bounded. + +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +```julia-repl +julia> P = Polyhedron([1 -3; -1 1; -1 0; 0 -1],[1,1,1,1]); - Check whether a polyhedron is bounded. +julia> isbounded(P) +false +``` """ isbounded(P::Polyhedron) = pm_polytope(P).BOUNDED """ - isfulldimensional(P::Polyhedron) + isfulldimensional(P) - Check whether a polyhedron is full dimensional. +Check whether a polyhedron is full dimensional. + +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +```julia-repl +julia> V = [1 2 3; 1 3 2; 2 1 3; 2 3 1; 3 1 2; 3 2 1]; + +julia> isfulldimensional(convex_hull(V)) +false +``` """ isfulldimensional(P::Polyhedron) = pm_polytope(P).FULL_DIM """ - f_vector(P::Polyhedron) + f_vector(P) + +Compute the vector`(f_1,f_2,...,f_(dim(P)-1))` where `f_i` is the number of faces of `P` of dimension `i`. -Computes the vector`(f_1,f_2,...,f_(dim(P)-1))` where `f_i` is the number of faces of `P` of dimension `i`. +# Arguments +- `P::Polyhedron`: A polyhedron. + +# Examples +Here we compute the f-vector of the 5-cube: +```julia-repl +julia> f_vector(cube(5)) +5-element Array{Int64,1}: + 32 + 80 + 80 + 40 + 10 +``` """ function f_vector(P::Polyhedron) f_vec=[length(collect(faces(P,i))) for i in 0:dim(P)-1] @@ -384,10 +863,29 @@ end """ - support_function(P::Polyhedron; convention = :max) + support_function(P; convention = :max) + +Produce a function `h(ω) = max{dot(x,ω) | x ∈ P}`. max may be changed +to min by setting convention = :min. + +# Arguments +- `P::Polyhedron`: A polyhedron. +- `convention::Symbol`: Convention to be applied. -Produces a function `h(ω) = max{dot(x,ω) | x ∈ P}`. max may be changed - to min by setting convention = :min. +# Examples +```julia-repl +julia> P = cube(3) + simplex(3); + +julia> φ = support_function(P); + +julia> φ([1,2,3]) +9 + +julia> ψ = support_function(P, convention = :min); + +julia> ψ([1,2,3]) +-6 +``` """ function support_function(P::Polyhedron; convention = :max) function h(ω::AbstractVector) @@ -396,5 +894,3 @@ function support_function(P::Polyhedron; convention = :max) end return h end - - diff --git a/src/Polytopes/Polyhedron/standard_constructions.jl b/src/Polytopes/Polyhedron/standard_constructions.jl index 4b0bd43b9e4a..e29652703457 100644 --- a/src/Polytopes/Polyhedron/standard_constructions.jl +++ b/src/Polytopes/Polyhedron/standard_constructions.jl @@ -5,9 +5,39 @@ ############################################################################### @doc Markdown.doc""" - orbit_polytope(V::AbstractVecOrMat, G::PermGroup) + orbit_polytope(V, G) -Construct the convex hull of the orbit of the point(s) in $V$ under the action of $G$. +Construct the convex hull of the orbit of one or several points under the action of a permutation group. + +# Arguments +- `V::AbstractVecOrMat`: Initial point(s). +- `P::PermGroup`: A permutation group. + +# Examples +This will construct the $3$-dimensional permutahedron: +```julia-repl +julia> V = [1 2 3]; + +julia> G = symmetric_group(3); + +julia> P = orbit_polytope(V, G) +A polyhedron in ambient dimension 3 + +julia> collect(vertices(P)) +6-element Vector{Polymake.Vector{Polymake.Rational}}: + pm::Vector +1 2 3 + pm::Vector +1 3 2 + pm::Vector +2 1 3 + pm::Vector +2 3 1 + pm::Vector +3 1 2 + pm::Vector +3 2 1 +``` """ function orbit_polytope(V::AbstractMatrix, G::PermGroup) if size(V)[2] != degree(G) @@ -23,11 +53,26 @@ function orbit_polytope(V::AbstractVector, G::PermGroup) end @doc Markdown.doc""" - cube(d [, u, l]) + cube(d [, l, u]) -Construct the $[-1,1]$-cube in dimension $d$. If $u$ and $l$ are given, the $[l,u]$-cube in dimension $d$ is returned. -""" cube(d) = Polyhedron(Polymake.polytope.cube(d)) -cube(d, u, l) = Polyhedron(Polymake.polytope.cube(d, u, l)) +Construct the $[-1,1]$-cube in dimension $d$. If $l$ and $u$ are given, the $[l,u]$-cube in dimension $d$ is returned. + +# Arguments +- `d::Int`: Dimension of the cube. +- `l::Rational`: Lower bound for each coordinate. +- `u::Rational`: Upper bound for each coordinate. + +# Examples +In this example the 5-dimensional unit cube is constructed to ask for one of its properties: +```julia-repl +julia> C = cube(5,0,1); + +julia> normalized_volume(C) +120 +``` +""" +cube(d) = Polyhedron(Polymake.polytope.cube(d)) +cube(d, l, u) = Polyhedron(Polymake.polytope.cube(d, u, l)) @@ -35,6 +80,30 @@ cube(d, u, l) = Polyhedron(Polymake.polytope.cube(d, u, l)) newton_polytope(poly) Compute the Newton polytope of the given polynomial `poly`. + +# Arguments +- `poly::Polynomial`: A multivariate polynomial. + +# Examples +```julia-repl +julia> S, (x, y) = PolynomialRing(ZZ, ["x", "y"]) +(Multivariate Polynomial Ring in x, y over Integer Ring, fmpz_mpoly[x, y]) + +julia> f = x^3*y + 3x*y^2 + 1 +x^3*y + 3*x*y^2 + 1 + +julia> NP = newton_polytope(f) +A polyhedron in ambient dimension 2 + +julia> collect(vertices(NP)) +3-element Array{Polymake.Vector{Polymake.Rational},1}: + pm::Vector +3 1 + pm::Vector +1 2 + pm::Vector +0 0 +``` """ function newton_polytope(f) exponents = reduce(hcat, Oscar.exponent_vectors(f))' @@ -44,10 +113,32 @@ end -""" - intersect(P::Polyhedron, Q::Polyhedron) +@doc Markdown.doc""" + intersect(P, Q) - Intersect two polyhedra. +Intersect two polyhedra. + +# Arguments +- `P::Polyhedron`: First polyhedron. +- `Q::Polyhedron`: Second polyhedron. + +# Examples +The positive orthant of the plane is the intersection of the two halfspaces with $x>0$ and $y>0$ respectively. +```julia-repl +julia> UH1 = convex_hull([0 0],[1 0],[0 1]); + +julia> UH2 = convex_hull([0 0],[0 1],[1 0]); + +julia> PO = intersect(UH1, UH2) +A polyhedron in ambient dimension 2 + +julia> collect(rays(PO)) +2-element Vector{Polymake.Vector{Polymake.Rational}}: + pm::Vector +1 0 + pm::Vector +0 1 +``` """ function intersect(P::Polyhedron, Q::Polyhedron) return Polyhedron(Polymake.polytope.intersection(pm_polytope(P), pm_polytope(Q))) @@ -55,9 +146,27 @@ end """ - minkowski_sum(P::Polyhedron, Q::Polyhedron) + minkowski_sum(P::Polyhedron, Q::Polyhedron) + +Minkowski sum of two polyhedra. + +# Arguments +- `P::Polyhedron`: First polyhedron. +- `Q::Polyhedron`: Second polyhedron. - Minkowski sum of two polyhedra. +# Examples +The Minkowski sum of a square and the 2-dimensional cross-polytope is an octagon: +```julia-repl +julia> P = cube(2); + +julia> Q = cross(2); + +julia> M = minkowski_sum(P, Q) +A polyhedron in ambient dimension 2 + +julia> nvertices(M) +8 +``` """ function minkowski_sum(P::Polyhedron, Q::Polyhedron; algorithm::Symbol=:standard) if algorithm == :standard @@ -65,7 +174,7 @@ function minkowski_sum(P::Polyhedron, Q::Polyhedron; algorithm::Symbol=:standard elseif algorithm == :fukuda return Polyhedron(Polymake.polytope.minkowski_sum_fukuda(pm_polytope(P), pm_polytope(Q))) else - throw(ArgumentError("Unknown minkowski sum `algorithm` argument :" * string(algorithm))) + throw(ArgumentError("Unknown minkowski sum `algorithm` argument:" * string(algorithm))) end end @@ -75,35 +184,114 @@ end #TODO: documentation + extend to different fields. """ - +(P::Polyhedron, Q::Polyhedron) + +(P::Polyhedron, Q::Polyhedron) + +Minkowski sum of two polyhedra. + +# Arguments +- `P::Polyhedron`: First polyhedron. +- `Q::Polyhedron`: Second polyhedron. + +# Examples +The Minkowski sum of a square and the 2-dimensional cross-polytope is an octagon: +```julia-repl +julia> P = cube(2); - Minkowski sum of two polyhedra. +julia> Q = cross(2); + +julia> M = minkowski_sum(P, Q) +A polyhedron in ambient dimension 2 + +julia> nvertices(M) +8 +``` """ +(P::Polyhedron, Q::Polyhedron) = minkowski_sum(P,Q) #TODO: extend to different fields -""" - *(k::Int, Q::Polyhedron) +@doc Markdown.doc""" + *(k::Int, Q::Polyhedron) + +Return the scaled polyhedron `kQ`. + +# Arguments +- `k::Int`: Scaling factor. +- `Q::Polyhedron`: A polyhedron. - Returns the scaled polyhedron `kQ`. +# Examples +Scaling an $n$-dimensional bounded polyhedron by the factor $k$ results in the volume being scaled by $k^n$. +This example confirms the statement for the 6-dimensional cube and $k = 2$. +```julia-repl +julia> C = cube(6); + +julia> SC = 2*C +A polyhedron in ambient dimension 6 + +julia> volume(SC)//volume(C) +64 +``` """ *(k::Int, P::Polyhedron) = Polyhedron(Polymake.polytope.scale(pm_polytope(P),k)) -""" - *(P::Polyhedron, k::Int) +@doc Markdown.doc""" + *(P::Polyhedron, k::Int) + +Return the scaled polyhedron `kP`. + +# Arguments +- `k::Int`: Scaling factor. +- `Q::Polyhedron`: A polyhedron. - Returns the scaled polyhedron `kP`. +# Examples +Scaling an $n$-dimensional bounded polyhedron by the factor $k$ results in the volume being scaled by $k^n$. +This example confirms the statement for the 6-dimensional cube and $k = 2$. +```julia-repl +julia> C = cube(6); + +julia> SC = C*2 +A polyhedron in ambient dimension 6 + +julia> volume(SC)//volume(C) +64 +``` """ *(P::Polyhedron,k::Int) = k*P -""" - +(P::Polyhedron, v::AbstractVector) +@doc Markdown.doc""" + +(P::Polyhedron, v::AbstractVector) - Returns the translation `P+v` of `P` by the vector `v`. +Return the translation `P+v` of `P` by the vector `v`. + +# Arguments +- `P::Polyhedron`: A polyhedron. +- `v::AbstractVector`: A vector of the same dimension as the ambient space of `P`. + +# Examples +We construct a polyhedron from its $V$-description. Shifting it by the right vector reveals that its inner geometry +corresponds to that of the 3-simplex. +```julia-repl +julia> P = convex_hull([100 200 300; 101 200 300; 100 201 300; 100 200 301]); + +julia> v = [-100, -200, -300]; + +julia> S = P + v +A polyhedron in ambient dimension 3 + +julia> collect(vertices(S)) +4-element Vector{Polymake.Vector{Polymake.Rational}}: + pm::Vector +0 0 0 + pm::Vector +1 0 0 + pm::Vector +0 1 0 + pm::Vector +0 0 1 +``` """ function +(P::Polyhedron,v::AbstractVector) if ambient_dim(P) != length(v) @@ -114,19 +302,98 @@ function +(P::Polyhedron,v::AbstractVector) end -""" - +(v::AbstractVector,P::Polyhedron) +@doc Markdown.doc""" + +(v::AbstractVector,P::Polyhedron) + +Return the translation `v+P` of `P` by the vector `v`. - Returns the translation `v+P` of `P` by the vector `v`. +# Arguments +- `P::Polyhedron`: A polyhedron. +- `v::AbstractVector`: A vector of the same dimension as the ambient space of `P`. + +# Examples +We construct a polyhedron from its $V$-description. Shifting it by the right vector reveals that its inner geometry +corresponds to that of the 3-simplex. +```julia-repl +julia> P = convex_hull([100 200 300; 101 200 300; 100 201 300; 100 200 301]); + +julia> v = [-100, -200, -300]; + +julia> S = v + P +A polyhedron in ambient dimension 3 + +julia> collect(vertices(S)) +4-element Vector{Polymake.Vector{Polymake.Rational}}: + pm::Vector +0 0 0 + pm::Vector +1 0 0 + pm::Vector +0 1 0 + pm::Vector +0 0 1 +``` """ +(v::AbstractVector,P::Polyhedron) = P+v @doc Markdown.doc""" - simplex(d[,n]) + simplex(d[,n]) Construct the simplex which is the convex hull of the standard basis vectors -along with the origin in R^$d$, optionally scaled by $n$. +along with the origin in $\mathbb{R}^d$, optionally scaled by $n$. + +# Arguments +- `d::Int`: Dimension of the simplex (and its ambient space). +- `n::Scalar`: Scaling factor. + +# Examples +Here we take a look at the facets of the 7-simplex and a scaled 7-simplex: +```julia-repl +julia> s = simplex(7) +A polyhedron in ambient dimension 7 + +julia> collect(facets(s)) +8-element Vector{Tuple{Polymake.Vector{Polymake.Rational}, Polymake.Rational}}: + (pm::Vector +-1 0 0 0 0 0 0, 0) + (pm::Vector +0 -1 0 0 0 0 0, 0) + (pm::Vector +0 0 -1 0 0 0 0, 0) + (pm::Vector +0 0 0 -1 0 0 0, 0) + (pm::Vector +0 0 0 0 -1 0 0, 0) + (pm::Vector +0 0 0 0 0 -1 0, 0) + (pm::Vector +0 0 0 0 0 0 -1, 0) + (pm::Vector +1 1 1 1 1 1 1, 1) + +julia> t = simplex(7, 5) +A polyhedron in ambient dimension 7 + +julia> collect(facets(t)) +8-element Vector{Tuple{Polymake.Vector{Polymake.Rational}, Polymake.Rational}}: + (pm::Vector +-1 0 0 0 0 0 0, 0) + (pm::Vector +0 -1 0 0 0 0 0, 0) + (pm::Vector +0 0 -1 0 0 0 0, 0) + (pm::Vector +0 0 0 -1 0 0 0, 0) + (pm::Vector +0 0 0 0 -1 0 0, 0) + (pm::Vector +0 0 0 0 0 -1 0, 0) + (pm::Vector +0 0 0 0 0 0 -1, 0) + (pm::Vector +1 1 1 1 1 1 1, 5) +``` """ simplex(d::Int64,n) = Polyhedron(Polymake.polytope.simplex(d,n)) simplex(d::Int64) = Polyhedron(Polymake.polytope.simplex(d)) @@ -134,24 +401,75 @@ simplex(d::Int64) = Polyhedron(Polymake.polytope.simplex(d)) @doc Markdown.doc""" - cross(d[,n]) + cross(d[,n]) + +Construct a $d$-dimensional cross polytope around origin with vertices located at $\pm e_i$ for each unit vector $e_i$ of $R^d$, scaled by $n$. -Construct a $d$-dimensional cross polytope around origin with vertices located at $\pm e_i$ for each unit vector $e_i$ of $R^d$. -If $n$ is not given, construct the unit cross polytope around origin. +# Arguments +- `d::Int`: Dimension of the cross polytope (and its ambient space). +- `n::Scalar`: Scaling factor. + +# Examples +Here we print the facets of a non-scaled and a scaled 3-dimensional cross polytope: +```julia-repl +julia> C = cross(3) +A polyhedron in ambient dimension 3 + +julia> collect(facets(C)) +8-element Vector{Tuple{Polymake.Vector{Polymake.Rational}, Polymake.Rational}}: + (pm::Vector +1 1 1, 1) + (pm::Vector +-1 1 1, 1) + (pm::Vector +1 -1 1, 1) + (pm::Vector +-1 -1 1, 1) + (pm::Vector +1 1 -1, 1) + (pm::Vector +-1 1 -1, 1) + (pm::Vector +1 -1 -1, 1) + (pm::Vector +-1 -1 -1, 1) + +julia> D = cross(3, 2) +A polyhedron in ambient dimension 3 + +julia> collect(facets(D)) +8-element Vector{Tuple{Polymake.Vector{Polymake.Rational}, Polymake.Rational}}: + (pm::Vector +1 1 1, 2) + (pm::Vector +-1 1 1, 2) + (pm::Vector +1 -1 1, 2) + (pm::Vector +-1 -1 1, 2) + (pm::Vector +1 1 -1, 2) + (pm::Vector +-1 1 -1, 2) + (pm::Vector +1 -1 -1, 2) + (pm::Vector +-1 -1 -1, 2) +``` """ cross(d::Int64,n) = Polyhedron(Polymake.polytope.cross(d,n)) cross(d::Int64) = Polyhedron(Polymake.polytope.cross(d)) @doc Markdown.doc""" - archimedean_solid(s) + archimedean_solid(s) Construct an Archimedean solid with the name given by String `s` from the list below. The polytopes are realized with floating point numbers and thus not exact; Vertex-facet-incidences are correct in all cases. # Arguments -- `s::String`: the name of the desired Archimedean solid +- `s::String`: The name of the desired Archimedean solid. Possible values: @@ -185,13 +503,27 @@ exact; Vertex-facet-incidences are correct in all cases. Regular polytope with 30 square, 20 hexagonal and 12 decagonal facets. "snub_dodecahedron" : Snub Dodecahedron. - Regular polytope with 80 triangular and 12 pentagonal facets. The - vertices are realized as floating point numbers. This is a chiral - polytope. + Regular polytope with 80 triangular and 12 pentagonal facets. + The vertices are realized as floating point numbers. + This is a chiral polytope. + +# Examples +```julia-repl +julia> T = archimedean_solid("cuboctahedron") +A polyhedron in ambient dimension 3 + +julia> sum([nvertices(F) for F in faces(T, 2)] .== 3) +8 + +julia> sum([nvertices(F) for F in faces(T, 2)] .== 4) +6 + +julia> nfacets(T) +14 +``` """ archimedean_solid(s::String) = Polyhedron(Polymake.polytope.archimedean_solid(s)) - @doc Markdown.doc""" upper_bound_theorem(d::Int, n::Int) diff --git a/src/Polytopes/Serialization.jl b/src/Polytopes/Serialization.jl index 801f2cd78c91..02e5635e5493 100644 --- a/src/Polytopes/Serialization.jl +++ b/src/Polytopes/Serialization.jl @@ -1,6 +1,6 @@ ############################################################################## """ - save_cone(Cone, String) + save_cone(Cone, String) Save a cone to a file in JSON format. The first argument is the cone, the second argument is the filename. @@ -11,7 +11,7 @@ function save_cone(C::Cone, filename::String) end """ - load_cone(String) + load_cone(String) Load a cone stored in JSON format, given the filename as input. """ @@ -25,7 +25,7 @@ end ############################################################################## """ - save_polyhedron(Cone, String) + save_polyhedron(Cone, String) Save a polyhedron to a file in JSON format. The first argument is the polyhedron, the second argument is the filename. @@ -36,7 +36,7 @@ function save_polyhedron(P::Polyhedron, filename::String) end """ - load_polyhedron(String) + load_polyhedron(String) Load a polyhedron stored in JSON format, given the filename as input. """ @@ -51,7 +51,7 @@ end ############################################################################## """ - save_polyhedralfan(PolyhedralFan, String) + save_polyhedralfan(PolyhedralFan, String) Save a polyhedral fan to a file in JSON format. The first argument is the polyhedral fan, the second argument is the filename. @@ -62,7 +62,7 @@ function save_polyhedralfan(PF::PolyhedralFan, filename::String) end """ - load_polyhedralfan(String) + load_polyhedralfan(String) Load a polyhedral fan stored in JSON format, given the filename as input. """ diff --git a/src/Polytopes/helpers.jl b/src/Polytopes/helpers.jl index ba7a5321d47e..290b30a1ae2d 100644 --- a/src/Polytopes/helpers.jl +++ b/src/Polytopes/helpers.jl @@ -71,6 +71,8 @@ stack(A::AbstractMatrix, B::AbstractMatrix) = [A; B] stack(A::AbstractMatrix, B::AbstractVector) = isempty(B) ? A : [A; B'] stack(A::AbstractVector, B::AbstractMatrix) = isempty(A) ? B : [A'; B] stack(A::AbstractVector, B::AbstractVector) = isempty(A) ? B : [A'; B'] +stack(A::AbstractVector,nothing) = A' +stack(nothing,B::AbstractVector) = B' #= function stack(A::Array{Polymake.Vector{Polymake.Rational},1}) if length(A)==2 @@ -85,7 +87,7 @@ end =# """ - decompose_vdata(A::AbstractMatrix) + decompose_vdata(A::AbstractMatrix) Given a (homogeneous) polymake matrix split into vertices and rays and dehomogenize. """ diff --git a/test/Polytopes/Group.jl b/test/Polytopes/Group.jl index 851bde027ea5..f2519fda9d47 100644 --- a/test/Polytopes/Group.jl +++ b/test/Polytopes/Group.jl @@ -7,7 +7,7 @@ P = convex_hull(M) @test ambient_dim(P) == 4 - F = facets(P, Polyhedra) + F = facets(Polyhedra, P) #@test nvertices.(F) == [6, 6, 4, 6, 4, 4, 6, 4, 6, 6, 4, 6, 6, 4] #Since different convex hull algorithms will result in differnt vertex orders, this test may fail. #Best way to handle this is to use "prefer" option in polymake, which is not available in OSCAR yet diff --git a/test/Polytopes/Polyhedron.jl b/test/Polytopes/Polyhedron.jl index f196082ece09..36229b5aeb1f 100644 --- a/test/Polytopes/Polyhedron.jl +++ b/test/Polytopes/Polyhedron.jl @@ -7,7 +7,7 @@ Q1 = convex_hull(pts, [1 1]) Q2 = convex_hull(pts, [1 1], [1 1]) C0 = cube(2) - C1 = cube(2, 1, 0) + C1 = cube(2, 0, 1) Pos=Polyhedron([-1 0 0; 0 -1 0; 0 0 -1],[0,0,0]) point = convex_hull([0 1 0]) s = simplex(2) diff --git a/test/Polytopes/extended.jl b/test/Polytopes/extended.jl index d4591c82edb3..ff2de197100a 100644 --- a/test/Polytopes/extended.jl +++ b/test/Polytopes/extended.jl @@ -7,7 +7,7 @@ Q1 = convex_hull(pts, [1 1]) Q2 = convex_hull(pts, [1 1], [1 1]) C0 = cube(2) - C1 = cube(2, 1, 0) + C1 = cube(2, 0, 1) @testset "(de)homogenize" begin dehomogenize, homogenize = Oscar.dehomogenize, Oscar.homogenize @@ -63,8 +63,8 @@ @test isnormal(C0) @test isfeasible(C0) @test isfulldimensional(C0) - @test minkowski_sum(C0,C0) == cube(2,2,-2) - @test minkowski_sum(C0,C0; algorithm=:fukuda) == cube(2,2,-2) + @test minkowski_sum(C0,C0) == cube(2,-2,2) + @test minkowski_sum(C0,C0; algorithm=:fukuda) == cube(2,-2, 2) @test intersect(C0,C0) == C0 end