diff --git a/docs/src/PolyhedralGeometry/Polyhedra/auxiliary.md b/docs/src/PolyhedralGeometry/Polyhedra/auxiliary.md index adebd185c2f3..4567224be2f0 100644 --- a/docs/src/PolyhedralGeometry/Polyhedra/auxiliary.md +++ b/docs/src/PolyhedralGeometry/Polyhedra/auxiliary.md @@ -64,7 +64,7 @@ is_smooth(P::Polyhedron{QQFieldElem}) is_very_ample(P::Polyhedron{QQFieldElem}) lattice_points(P::Polyhedron{QQFieldElem}) lattice_volume(P::Polyhedron{QQFieldElem}) -normalized_volume(P::Polyhedron{T}) where T<:scalar_types +normalized_volume(P::Polyhedron) polarize(P::Polyhedron{T}) where T<:scalar_types project_full(P::Polyhedron{T}) where T<:scalar_types print_constraints(A::AnyVecOrMat, b::AbstractVector; trivial::Bool = false) diff --git a/docs/src/PolyhedralGeometry/Polyhedra/constructions.md b/docs/src/PolyhedralGeometry/Polyhedra/constructions.md index 777d26eef733..e6992c1bd70f 100644 --- a/docs/src/PolyhedralGeometry/Polyhedra/constructions.md +++ b/docs/src/PolyhedralGeometry/Polyhedra/constructions.md @@ -17,8 +17,8 @@ from other objects in OSCAR. ### Intersecting halfspaces: $H$-representation ```@docs -polyhedron(::Type{T}, A::AnyVecOrMat, b::AbstractVector) where T<:scalar_types -polyhedron(::Type{T}, I::Union{Nothing, AbstractCollection[AffineHalfspace]}, E::Union{Nothing, AbstractCollection[AffineHyperplane]} = nothing) where T<:scalar_types +polyhedron(::Oscar.scalar_type_or_field, A::AnyVecOrMat, b::AbstractVector) +polyhedron(::Oscar.scalar_type_or_field, I::Union{Nothing, AbstractCollection[AffineHalfspace]}, E::Union{Nothing, AbstractCollection[AffineHyperplane]} = nothing) ``` The complete $H$-representation can be retrieved using [`facets`](@ref facets) @@ -54,7 +54,7 @@ true ### Computing convex hulls: $V$-representation ```@docs -convex_hull(::Type{T}, ::AnyVecOrMat; non_redundant::Bool=false) where T<:scalar_types +convex_hull(::Oscar.scalar_type_or_field, ::AnyVecOrMat; non_redundant::Bool=false) ``` This is a standard triangle, defined via a (redundant) $V$-representation and diff --git a/docs/src/PolyhedralGeometry/cones.md b/docs/src/PolyhedralGeometry/cones.md index 92effbd2bc7e..38ea2b4f2690 100644 --- a/docs/src/PolyhedralGeometry/cones.md +++ b/docs/src/PolyhedralGeometry/cones.md @@ -27,7 +27,7 @@ reason for keeping cones as a distinct type. ## Construction ```@docs -positive_hull(::Type{T}, ::Union{Oscar.MatElem, AbstractMatrix}) where T<:scalar_types +positive_hull(::Oscar.scalar_type_or_field, ::Union{Oscar.MatElem, AbstractMatrix}) cone_from_inequalities cone_from_equations secondary_cone(SOP::SubdivisionOfPoints{T}) where T<:scalar_types diff --git a/docs/src/PolyhedralGeometry/intro.md b/docs/src/PolyhedralGeometry/intro.md index 672a0b804e77..f185306a969b 100644 --- a/docs/src/PolyhedralGeometry/intro.md +++ b/docs/src/PolyhedralGeometry/intro.md @@ -72,6 +72,13 @@ any of the corresponding notions below. ### Vectors +There are two specialized `Vector`-like types, `PointVector` and `RayVector`, which commonly are returned by functions from Polyhedral Geometry. These can also be manually constructed: + +```@docs +point_vector +ray_vector +``` + While `RayVector`s can not be used do describe `PointVector`s (and vice versa), matrices are generally allowed. @@ -98,6 +105,15 @@ Type | A `RayVector` corresponds to... ### Halfspaces and Hyperplanes +Similar to points and rays, there are types `AffineHalfspace`, `LinearHalfspace`, `AffineHyperplane` and `LinearHyperplane`: + +```@docs +affine_halfspace +linear_halfspace +affine_hyperplane +linear_hyperplane +``` + These collections allow to mix up affine halfspaces/hyperplanes and their linear counterparts, but note that an error will be produced when trying to convert an affine description with bias not equal to zero to a linear description. diff --git a/docs/src/PolyhedralGeometry/subdivisions_of_points.md b/docs/src/PolyhedralGeometry/subdivisions_of_points.md index 3b2f3b9a0825..4e4e9039af95 100644 --- a/docs/src/PolyhedralGeometry/subdivisions_of_points.md +++ b/docs/src/PolyhedralGeometry/subdivisions_of_points.md @@ -34,8 +34,8 @@ vector, but if it does, it is called `regular`. ```@docs -subdivision_of_points(::Type{T}, Points::Union{Oscar.MatElem,AbstractMatrix}, Incidence::IncidenceMatrix) where T<:scalar_types -subdivision_of_points(::Type{T}, Points::Union{Oscar.MatElem,AbstractMatrix}, Weights::AbstractVector) where T<:scalar_types +subdivision_of_points(::Oscar.scalar_type_or_field, Points::Union{Oscar.MatElem,AbstractMatrix}, Incidence::IncidenceMatrix) +subdivision_of_points(::Oscar.scalar_type_or_field, Points::Union{Oscar.MatElem,AbstractMatrix}, Weights::AbstractVector) ``` From a subdivision of points one can construct the @@ -50,6 +50,6 @@ is_regular(SOP::SubdivisionOfPoints) maximal_cells min_weights n_maximal_cells(SOP::SubdivisionOfPoints) -points(SOP::SubdivisionOfPoints) +points(SOP::SubdivisionOfPoints{T}) where T<:scalar_types secondary_cone ``` diff --git a/src/AlgebraicGeometry/Surfaces/K3Auto.jl b/src/AlgebraicGeometry/Surfaces/K3Auto.jl index e2806317765a..c1a00dde74fa 100644 --- a/src/AlgebraicGeometry/Surfaces/K3Auto.jl +++ b/src/AlgebraicGeometry/Surfaces/K3Auto.jl @@ -1594,7 +1594,7 @@ function span_in_S(L, S, weyl) if N==0 M = zero_matrix(QQ, 0, rank(S)) else - M = reduce(vcat, (matrix(QQ, 1, rank(S), v.a) for v in spanC)) + M = linear_equation_matrix(spanC) end k, K = kernel(M) gensN = transpose(K)[1:k,:] diff --git a/src/AlgebraicGeometry/ToricVarieties/Proj/constructors.jl b/src/AlgebraicGeometry/ToricVarieties/Proj/constructors.jl index 90f7fe73a8f1..f7ac0894c643 100644 --- a/src/AlgebraicGeometry/ToricVarieties/Proj/constructors.jl +++ b/src/AlgebraicGeometry/ToricVarieties/Proj/constructors.jl @@ -74,14 +74,14 @@ function _proj_and_total_space(is_proj::Bool, E::Vector{T}) where T <: Union{Tor end end - total_rays_gens = vcat([modified_ray_gens[ray] for ray in rays(v)], [vcat(RayVector(zeros(Int64, dim(v))), l[i]) for i in eachindex(E)]) + total_rays_gens = vcat([modified_ray_gens[ray] for ray in rays(v)], [vcat(ray_vector(zeros(Int64, dim(v))), l[i]) for i in eachindex(E)]) return normal_toric_variety(polyhedral_fan(total_rays_gens, IncidenceMatrix(new_maximal_cones))) end function _m_sigma(sigma::Cone{QQFieldElem}, pol_sigma::Cone{QQFieldElem}, D::Union{ToricDivisor, ToricLineBundle})::RayVector{QQFieldElem} - ans = RayVector(zeros(QQFieldElem, dim(sigma))) + ans = ray_vector(zeros(QQFieldElem, dim(sigma))) coeff = coefficients(isa(D, ToricDivisor) ? D : toric_divisor(D)) dual_ray = QQFieldElem[] diff --git a/src/PolyhedralGeometry/Cone/constructors.jl b/src/PolyhedralGeometry/Cone/constructors.jl index ade7eae79831..31b0d06d989a 100644 --- a/src/PolyhedralGeometry/Cone/constructors.jl +++ b/src/PolyhedralGeometry/Cone/constructors.jl @@ -62,25 +62,25 @@ julia> HS = positive_hull(R, L) Polyhedral cone in ambient dimension 2 ``` """ -function positive_hull(f::Union{Type{T}, Field}, R::AbstractCollection[RayVector], L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) where T<:scalar_types - parent_field, scalar_type = _determine_parent_and_scalar(f, R, L) - inputrays = remove_zero_rows(unhomogenized_matrix(R)) - if isnothing(L) || isempty(L) - L = Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, _ambient_dim(R)) - end - - if non_redundant - return Cone{scalar_type}(Polymake.polytope.Cone{_scalar_type_to_polymake(scalar_type)}(RAYS = inputrays, LINEALITY_SPACE = unhomogenized_matrix(L),), parent_field) - else - return Cone{scalar_type}(Polymake.polytope.Cone{_scalar_type_to_polymake(scalar_type)}(INPUT_RAYS = inputrays, INPUT_LINEALITY = unhomogenized_matrix(L),), parent_field) - end +function positive_hull(f::scalar_type_or_field, R::AbstractCollection[RayVector], L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) + parent_field, scalar_type = _determine_parent_and_scalar(f, R, L) + inputrays = remove_zero_rows(unhomogenized_matrix(R)) + if isnothing(L) || isempty(L) + L = Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, _ambient_dim(R)) + end + + if non_redundant + return Cone{scalar_type}(Polymake.polytope.Cone{_scalar_type_to_polymake(scalar_type)}(RAYS = inputrays, LINEALITY_SPACE = unhomogenized_matrix(L),), parent_field) + else + return Cone{scalar_type}(Polymake.polytope.Cone{_scalar_type_to_polymake(scalar_type)}(INPUT_RAYS = inputrays, INPUT_LINEALITY = unhomogenized_matrix(L),), parent_field) + end end # Redirect everything to the above constructor, use QQFieldElem as default for the # scalar type T. positive_hull(R::AbstractCollection[RayVector], L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) = positive_hull(QQFieldElem, R, L; non_redundant=non_redundant) cone(R::AbstractCollection[RayVector], L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) = positive_hull(QQFieldElem, R, L; non_redundant=non_redundant) -cone(f::Union{Type{T}, Field}, R::AbstractCollection[RayVector], L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) where T<:scalar_types = positive_hull(f, R, L; non_redundant=non_redundant) -cone(f::Union{Type{T}, Field}, x...) where T<:scalar_types = positive_hull(f, x...) +cone(f::scalar_type_or_field, R::AbstractCollection[RayVector], L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) = positive_hull(f, R, L; non_redundant=non_redundant) +cone(f::scalar_type_or_field, x...) = positive_hull(f, x...) function ==(C0::Cone{T}, C1::Cone{T}) where T<:scalar_types @@ -115,16 +115,16 @@ julia> rays(C) [1, 1] ``` """ -function cone_from_inequalities(f::Union{Type{T}, Field}, I::AbstractCollection[LinearHalfspace], E::Union{Nothing, AbstractCollection[LinearHyperplane]} = nothing; non_redundant::Bool = false) where T<:scalar_types - parent_field, scalar_type = _determine_parent_and_scalar(f, I, E) - IM = -linear_matrix_for_polymake(I) - EM = isnothing(E) || isempty(E) ? Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(IM, 2)) : linear_matrix_for_polymake(E) - - if non_redundant - return Cone{scalar_type}(Polymake.polytope.Cone{_scalar_type_to_polymake(scalar_type)}(FACETS = IM, LINEAR_SPAN = EM), parent_field) - else - return Cone{scalar_type}(Polymake.polytope.Cone{_scalar_type_to_polymake(scalar_type)}(INEQUALITIES = IM, EQUATIONS = EM), parent_field) - end +function cone_from_inequalities(f::scalar_type_or_field, I::AbstractCollection[LinearHalfspace], E::Union{Nothing, AbstractCollection[LinearHyperplane]} = nothing; non_redundant::Bool = false) + parent_field, scalar_type = _determine_parent_and_scalar(f, I, E) + IM = -linear_matrix_for_polymake(I) + EM = isnothing(E) || isempty(E) ? Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(IM, 2)) : linear_matrix_for_polymake(E) + + if non_redundant + return Cone{scalar_type}(Polymake.polytope.Cone{_scalar_type_to_polymake(scalar_type)}(FACETS = IM, LINEAR_SPAN = EM), parent_field) + else + return Cone{scalar_type}(Polymake.polytope.Cone{_scalar_type_to_polymake(scalar_type)}(INEQUALITIES = IM, EQUATIONS = EM), parent_field) + end end @doc raw""" @@ -152,11 +152,11 @@ julia> dim(C) 1 ``` """ -function cone_from_equations(f::Union{Type{T}, Field}, E::AbstractCollection[LinearHyperplane]; non_redundant::Bool = false) where T<:scalar_types - parent_field, scalar_type = _determine_parent_and_scalar(f, E) - EM = linear_matrix_for_polymake(E) - IM = Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(EM, 2)) - return cone_from_inequalities(f, IM, EM; non_redundant = non_redundant) +function cone_from_equations(f::scalar_type_or_field, E::AbstractCollection[LinearHyperplane]; non_redundant::Bool = false) + parent_field, scalar_type = _determine_parent_and_scalar(f, E) + EM = linear_matrix_for_polymake(E) + IM = Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(EM, 2)) + return cone_from_inequalities(f, IM, EM; non_redundant = non_redundant) end cone_from_inequalities(x...) = cone_from_inequalities(QQFieldElem, x...) diff --git a/src/PolyhedralGeometry/Cone/properties.jl b/src/PolyhedralGeometry/Cone/properties.jl index c432b40358d5..9fcb39394b45 100644 --- a/src/PolyhedralGeometry/Cone/properties.jl +++ b/src/PolyhedralGeometry/Cone/properties.jl @@ -4,17 +4,17 @@ ############################################################################### ############################################################################### -rays(as::Type{RayVector{T}}, C::Cone) where T<:scalar_types = lineality_dim(C) == 0 ? _rays(as, C) : _empty_subobjectiterator(as, C) -_rays(as::Type{RayVector{T}}, C::Cone) where T<:scalar_types = SubObjectIterator{as}(C, _ray_cone, _nrays(C)) +rays(as::Type{RayVector{T}}, C::Cone{T}) where T<:scalar_types = lineality_dim(C) == 0 ? _rays(as, C) : _empty_subobjectiterator(as, C) +_rays(as::Type{RayVector{T}}, C::Cone{T}) where T<:scalar_types = SubObjectIterator{as}(C, _ray_cone, _nrays(C)) -_ray_cone(::Type{RayVector{T}}, C::Cone, i::Base.Integer) where T<:scalar_types = RayVector{T}(coefficient_field(C).(view(pm_object(C).RAYS, i, :))) +_ray_cone(U::Type{RayVector{T}}, C::Cone{T}, i::Base.Integer) where T<:scalar_types = ray_vector(coefficient_field(C), view(pm_object(C).RAYS, i, :))::U _vector_matrix(::Val{_ray_cone}, C::Cone; homogenized=false) = homogenized ? homogenize(pm_object(C).RAYS, 0) : pm_object(C).RAYS _matrix_for_polymake(::Val{_ray_cone}) = _vector_matrix -rays(::Type{RayVector}, C::Cone{T}) where T<:scalar_types = rays(RayVector{T}, C) -_rays(::Type{RayVector}, C::Cone{T}) where T<:scalar_types = _rays(RayVector{T}, C) +rays(::Type{<:RayVector}, C::Cone{T}) where T<:scalar_types = rays(RayVector{T}, C) +_rays(::Type{<:RayVector}, C::Cone{T}) where T<:scalar_types = _rays(RayVector{T}, C) @doc raw""" rays(C::Cone) @@ -111,13 +111,13 @@ julia> RML.lineality_basis ``` """ rays_modulo_lineality(C::Cone{T}) where T<:scalar_types = rays_modulo_lineality(NamedTuple{(:rays_modulo_lineality, :lineality_basis), Tuple{SubObjectIterator{RayVector{T}}, SubObjectIterator{RayVector{T}}}}, C) -function rays_modulo_lineality(as::Type{NamedTuple{(:rays_modulo_lineality, :lineality_basis), Tuple{SubObjectIterator{RayVector{T}}, SubObjectIterator{RayVector{T}}}}}, C::Cone) where T<:scalar_types +function rays_modulo_lineality(::Type{NamedTuple{(:rays_modulo_lineality, :lineality_basis), Tuple{SubObjectIterator{RayVector{T}}, SubObjectIterator{RayVector{T}}}}}, C::Cone{T}) where T<:scalar_types return ( rays_modulo_lineality = _rays(C), lineality_basis = lineality_space(C) ) end -rays_modulo_lineality(as::Type{RayVector}, C::Cone) = _rays(C) +rays_modulo_lineality(::Type{<:RayVector}, C::Cone) = _rays(C) @doc raw""" @@ -147,8 +147,11 @@ function faces(C::Cone{T}, face_dim::Int) where T<:scalar_types return SubObjectIterator{Cone{T}}(C, _face_cone, size(Polymake.polytope.faces_of_dim(pm_object(C), n), 1), (f_dim = n,)) end -function _face_cone(::Type{Cone{T}}, C::Cone, i::Base.Integer; f_dim::Int = 0) where T<:scalar_types - return Cone{T}(Polymake.polytope.Cone{_scalar_type_to_polymake(T)}(RAYS = pm_object(C).RAYS[collect(Polymake.to_one_based_indexing(Polymake.polytope.faces_of_dim(pm_object(C), f_dim)[i])), :], LINEALITY_SPACE = pm_object(C).LINEALITY_SPACE)) +function _face_cone(::Type{Cone{T}}, C::Cone{T}, i::Base.Integer; f_dim::Int = 0) where T<:scalar_types + R = pm_object(C).RAYS[collect(Polymake.to_one_based_indexing(Polymake.polytope.faces_of_dim(pm_object(C), f_dim)[i])), :] + L = pm_object(C).LINEALITY_SPACE + PT = _scalar_type_to_polymake(T) + return Cone{T}(Polymake.polytope.Cone{PT}(RAYS = R, LINEALITY_SPACE = L), coefficient_field(C)) end function _ray_indices(::Val{_face_cone}, C::Cone; f_dim::Int = 0) @@ -156,8 +159,11 @@ function _ray_indices(::Val{_face_cone}, C::Cone; f_dim::Int = 0) return IncidenceMatrix([collect(f[i]) for i in 1:length(f)]) end -function _face_cone_facet(::Type{Cone{T}}, C::Cone, i::Base.Integer) where T<:scalar_types - return Cone{T}(Polymake.polytope.Cone{_scalar_type_to_polymake(T)}(RAYS = pm_object(C).RAYS[collect(pm_object(C).RAYS_IN_FACETS[i, :]), :], LINEALITY_SPACE = pm_object(C).LINEALITY_SPACE), coefficient_field(C)) +function _face_cone_facet(::Type{Cone{T}}, C::Cone{T}, i::Base.Integer) where T<:scalar_types + R = pm_object(C).RAYS[collect(pm_object(C).RAYS_IN_FACETS[i, :]), :] + L = pm_object(C).LINEALITY_SPACE + PT = _scalar_type_to_polymake(T) + return Cone{T}(Polymake.polytope.Cone{PT}(RAYS = R, LINEALITY_SPACE = pm_object(C).LINEALITY_SPACE), coefficient_field(C)) end _ray_indices(::Val{_face_cone_facet}, C::Cone) = pm_object(C).RAYS_IN_FACETS @@ -439,11 +445,11 @@ julia> f = facets(Halfspace, c) -x₂ + x₃ ≦ 0 ``` """ -facets(as::Type{<:Union{LinearHalfspace{T}, Cone{T}}}, C::Cone) where T<:scalar_types = SubObjectIterator{as}(C, _facet_cone, nfacets(C)) +facets(as::Type{<:Union{LinearHalfspace{T}, Cone{T}}}, C::Cone{T}) where T<:scalar_types = SubObjectIterator{as}(C, _facet_cone, nfacets(C)) -_facet_cone(::Type{LinearHalfspace{T}}, C::Cone, i::Base.Integer) where T<:scalar_types = linear_halfspace(coefficient_field(C), -pm_object(C).FACETS[[i], :]) +_facet_cone(U::Type{LinearHalfspace{T}}, C::Cone{T}, i::Base.Integer) where T<:scalar_types = linear_halfspace(coefficient_field(C), -pm_object(C).FACETS[[i], :])::U -_facet_cone(::Type{Cone{T}}, C::Cone, i::Base.Integer) where T<:scalar_types = Cone{T}(Polymake.polytope.facet(pm_object(C), i-1), coefficient_field(C)) +_facet_cone(::Type{Cone{T}}, C::Cone{T}, i::Base.Integer) where T<:scalar_types = Cone{T}(Polymake.polytope.facet(pm_object(C), i-1), coefficient_field(C)) _linear_inequality_matrix(::Val{_facet_cone}, C::Cone) = -pm_object(C).FACETS @@ -455,7 +461,7 @@ _incidencematrix(::Val{_facet_cone}) = _ray_indices facets(C::Cone{T}) where T<:scalar_types = facets(LinearHalfspace{T}, C) -facets(::Type{Halfspace}, C::Cone{T}) where T<:scalar_types = facets(LinearHalfspace{T}, C) +facets(::Type{<:Halfspace}, C::Cone{T}) where T<:scalar_types = facets(LinearHalfspace{T}, C) facets(::Type{Cone}, C::Cone{T}) where T<:scalar_types = facets(Cone{T}, C) @@ -477,7 +483,7 @@ julia> lineality_space(UH) """ lineality_space(C::Cone{T}) where T<:scalar_types = SubObjectIterator{RayVector{T}}(C, _lineality_cone, lineality_dim(C)) -_lineality_cone(::Type{RayVector{T}}, C::Cone, i::Base.Integer) where T<:scalar_types = RayVector{T}(coefficient_field(C).(view(pm_object(C).LINEALITY_SPACE, i, :))) +_lineality_cone(U::Type{RayVector{T}}, C::Cone{T}, i::Base.Integer) where T<:scalar_types = ray_vector(coefficient_field(C), view(pm_object(C).LINEALITY_SPACE, i, :))::U _generator_matrix(::Val{_lineality_cone}, C::Cone; homogenized=false) = homogenized ? homogenize(pm_object(C).LINEALITY_SPACE, 0) : pm_object(C).LINEALITY_SPACE @@ -501,7 +507,7 @@ x₃ = 0 """ linear_span(C::Cone{T}) where T<:scalar_types = SubObjectIterator{LinearHyperplane{T}}(C, _linear_span, size(pm_object(C).LINEAR_SPAN, 1)) -_linear_span(::Type{LinearHyperplane{T}}, C::Cone, i::Base.Integer) where T<:scalar_types = linear_hyperplane(coefficient_field(C), view(pm_object(C).LINEAR_SPAN, i, :)) +_linear_span(U::Type{LinearHyperplane{T}}, C::Cone{T}, i::Base.Integer) where T<:scalar_types = linear_hyperplane(coefficient_field(C), view(pm_object(C).LINEAR_SPAN, i, :))::U _linear_equation_matrix(::Val{_linear_span}, C::Cone) = pm_object(C).LINEAR_SPAN @@ -530,7 +536,7 @@ function hilbert_basis(C::Cone{QQFieldElem}) return SubObjectIterator{PointVector{ZZRingElem}}(C, _hilbert_generator, size(pm_object(C).HILBERT_BASIS_GENERATORS[1], 1)) end -_hilbert_generator(::Type{PointVector{ZZRingElem}}, C::Cone, i::Base.Integer) = PointVector{ZZRingElem}(view(pm_object(C).HILBERT_BASIS_GENERATORS[1], i, :)) +_hilbert_generator(T::Type{PointVector{ZZRingElem}}, C::Cone{QQFieldElem}, i::Base.Integer) = point_vector(ZZ, view(pm_object(C).HILBERT_BASIS_GENERATORS[1], i, :))::T _generator_matrix(::Val{_hilbert_generator}, C::Cone; homogenized=false) = homogenized ? homogenize(pm_object(C).HILBERT_BASIS_GENERATORS[1], 0) : pm_object(C).HILBERT_BASIS_GENERATORS[1] @@ -586,4 +592,4 @@ Base.in(v::AbstractVector, C::Cone) = Polymake.polytope.contains(pm_object(C), v Compute a point in the relative interior point of `C`, i.e. a point in `C` not contained in any facet. """ -relative_interior_point(C::Cone{T}) where T<:scalar_types = PointVector{T}(coefficient_field(C), view(Polymake.common.dense(pm_object(C).REL_INT_POINT), :)) # broadcast_view +relative_interior_point(C::Cone{T}) where T<:scalar_types = point_vector(coefficient_field(C), view(Polymake.common.dense(pm_object(C).REL_INT_POINT), :))::PointVector{T} # broadcast_view diff --git a/src/PolyhedralGeometry/PolyhedralComplex/constructors.jl b/src/PolyhedralGeometry/PolyhedralComplex/constructors.jl index a0bd21e680e8..e508b354c753 100644 --- a/src/PolyhedralGeometry/PolyhedralComplex/constructors.jl +++ b/src/PolyhedralGeometry/PolyhedralComplex/constructors.jl @@ -75,36 +75,36 @@ julia> lineality_dim(PC) 1 ``` """ -function polyhedral_complex(f::Union{Type{T}, Field}, +function polyhedral_complex(f::scalar_type_or_field, polyhedra::IncidenceMatrix, vr::AbstractCollection[PointVector], far_vertices::Union{Vector{Int}, Nothing} = nothing, L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false - ) where T<:scalar_types - parent_field, scalar_type = _determine_parent_and_scalar(f, vr, L) - points = homogenized_matrix(vr, 1) - LM = isnothing(L) || isempty(L) ? zero_matrix(QQ, 0, size(points, 2)) : homogenized_matrix(L, 0) - - # Rays and Points are homogenized and combined and - # If some vertices are far vertices, give them a leading 0 - if !isnothing(far_vertices) - points[far_vertices,1] .= 0 - end - - if non_redundant - return PolyhedralComplex{scalar_type}(Polymake.fan.PolyhedralComplex{_scalar_type_to_polymake(scalar_type)}( - VERTICES = points, - LINEALITY_SPACE = LM, - MAXIMAL_CONES = polyhedra, - ), parent_field) - else - return PolyhedralComplex{scalar_type}(Polymake.fan.PolyhedralComplex{_scalar_type_to_polymake(scalar_type)}( - POINTS = points, - INPUT_LINEALITY = LM, - INPUT_CONES = polyhedra, - ), parent_field) - end + ) + parent_field, scalar_type = _determine_parent_and_scalar(f, vr, L) + points = homogenized_matrix(vr, 1) + LM = isnothing(L) || isempty(L) ? zero_matrix(QQ, 0, size(points, 2)) : homogenized_matrix(L, 0) + + # Rays and Points are homogenized and combined and + # If some vertices are far vertices, give them a leading 0 + if !isnothing(far_vertices) + points[far_vertices,1] .= 0 + end + + if non_redundant + return PolyhedralComplex{scalar_type}(Polymake.fan.PolyhedralComplex{_scalar_type_to_polymake(scalar_type)}( + VERTICES = points, + LINEALITY_SPACE = LM, + MAXIMAL_CONES = polyhedra, + ), parent_field) + else + return PolyhedralComplex{scalar_type}(Polymake.fan.PolyhedralComplex{_scalar_type_to_polymake(scalar_type)}( + POINTS = points, + INPUT_LINEALITY = LM, + INPUT_CONES = polyhedra, + ), parent_field) + end end # default scalar type: `QQFieldElem` polyhedral_complex(polyhedra::IncidenceMatrix, @@ -114,11 +114,11 @@ polyhedral_complex(polyhedra::IncidenceMatrix, non_redundant::Bool = false) = polyhedral_complex(QQFieldElem, polyhedra, vr, far_vertices, L; non_redundant=non_redundant) -function polyhedral_complex(f::Union{Type{T}, Field}, v::AbstractCollection[PointVector], vi::IncidenceMatrix, r::AbstractCollection[RayVector], ri::IncidenceMatrix, L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) where T<:scalar_types - vr = [unhomogenized_matrix(v); unhomogenized_matrix(r)] - far_vertices = collect((size(v, 1) + 1):size(vr, 1)) - polyhedra = hcat(vi, ri) - return polyhedral_complex(f, polyhedra, vr, far_vertices, L; non_redundant = non_redundant) +function polyhedral_complex(f::scalar_type_or_field, v::AbstractCollection[PointVector], vi::IncidenceMatrix, r::AbstractCollection[RayVector], ri::IncidenceMatrix, L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) + vr = [unhomogenized_matrix(v); unhomogenized_matrix(r)] + far_vertices = collect((size(v, 1) + 1):size(vr, 1)) + polyhedra = hcat(vi, ri) + return polyhedral_complex(f, polyhedra, vr, far_vertices, L; non_redundant = non_redundant) end # TODO: Only works for this specific case; implement generalization using `iter.Acc` diff --git a/src/PolyhedralGeometry/PolyhedralComplex/properties.jl b/src/PolyhedralGeometry/PolyhedralComplex/properties.jl index 5d0c39d2ea3b..96122284a7c5 100644 --- a/src/PolyhedralGeometry/PolyhedralComplex/properties.jl +++ b/src/PolyhedralGeometry/PolyhedralComplex/properties.jl @@ -60,15 +60,13 @@ julia> matrix(QQ, vertices(PC)) """ vertices(as::Type{PointVector{T}}, PC::PolyhedralComplex{T}) where {T<:scalar_types} = lineality_dim(PC) == 0 ? _vertices(as, PC) : _empty_subobjectiterator(as, PC) -_vertices(as::Type{PointVector{T}}, PC::PolyhedralComplex) where {T<:scalar_types} = +_vertices(as::Type{PointVector{T}}, PC::PolyhedralComplex{T}) where {T<:scalar_types} = SubObjectIterator{as}(PC, _vertex_complex, length(_vertex_indices(pm_object(PC)))) -_vertex_complex( - ::Type{PointVector{T}}, PC::PolyhedralComplex, i::Base.Integer -) where {T<:scalar_types} = PointVector{T}( +_vertex_complex(U::Type{PointVector{T}}, PC::PolyhedralComplex{T}, i::Base.Integer) where {T<:scalar_types} = point_vector( coefficient_field(PC), @view pm_object(PC).VERTICES[_vertex_indices(pm_object(PC))[i], 2:end] -) +)::U _point_matrix(::Val{_vertex_complex}, PC::PolyhedralComplex; homogenized=false) = @view pm_object(PC).VERTICES[_vertex_indices(pm_object(PC)), (homogenized ? 1 : 2):end] @@ -85,28 +83,25 @@ function _all_vertex_indices(P::Polymake.BigObject) return vi end -function _vertex_or_ray_complex( - ::Type{Union{PointVector{T},RayVector{T}}}, PC::PolyhedralComplex, i::Base.Integer -) where {T<:scalar_types} +function _vertex_or_ray_complex(::Type{Union{PointVector{T}, RayVector{T}}}, PC::PolyhedralComplex{T}, i::Base.Integer) where T<:scalar_types A = pm_object(PC).VERTICES if iszero(A[_all_vertex_indices(pm_object(PC))[i], 1]) - return RayVector{T}( + return ray_vector( coefficient_field(PC), @view pm_object(PC).VERTICES[_all_vertex_indices(pm_object(PC))[i], 2:end] - ) + )::RayVector{T} else - return PointVector{T}( + return point_vector( coefficient_field(PC), @view pm_object(PC).VERTICES[_all_vertex_indices(pm_object(PC))[i], 2:end] - ) + )::PointVector{T} end end -_vertices( - as::Type{Union{RayVector{T},PointVector{T}}}, PC::PolyhedralComplex -) where {T<:scalar_types} = SubObjectIterator{as}( - PC, _vertex_or_ray_complex, length(_all_vertex_indices(pm_object(PC))) -) +_vertices(as::Type{Union{RayVector{T}, PointVector{T}}}, PC::PolyhedralComplex{T}) where {T<:scalar_types} = + SubObjectIterator{as}( + PC, _vertex_or_ray_complex, length(_all_vertex_indices(pm_object(PC))) + ) vertices_and_rays(PC::PolyhedralComplex{T}) where {T<:scalar_types} = _vertices(Union{PointVector{T},RayVector{T}}, PC) @@ -120,7 +115,7 @@ _vertices(::Type{PointVector}, PC::PolyhedralComplex{T}) where {T<:scalar_types} vertices(PC::PolyhedralComplex{T}) where {T<:scalar_types} = vertices(PointVector{T}, PC) _vertices(PC::PolyhedralComplex{T}) where {T<:scalar_types} = _vertices(PointVector{T}, PC) -_rays(as::Type{RayVector{T}}, PC::PolyhedralComplex) where {T<:scalar_types} = +_rays(as::Type{RayVector{T}}, PC::PolyhedralComplex{T}) where {T<:scalar_types} = SubObjectIterator{as}( PC, _ray_polyhedral_complex, length(_ray_indices_polyhedral_complex(pm_object(PC))) ) @@ -165,28 +160,18 @@ julia> RML.lineality_basis """ rays_modulo_lineality(PC::PolyhedralComplex{T}) where {T<:scalar_types} = rays_modulo_lineality( - NamedTuple{ - (:rays_modulo_lineality, :lineality_basis), - Tuple{SubObjectIterator{RayVector{T}},SubObjectIterator{RayVector{T}}}, - }, - PC, - ) -function rays_modulo_lineality( - as::Type{ - NamedTuple{ - (:rays_modulo_lineality, :lineality_basis), - Tuple{SubObjectIterator{RayVector{T}},SubObjectIterator{RayVector{T}}}, - }, - }, - PC::PolyhedralComplex, -) where {T<:scalar_types} - return ( - rays_modulo_lineality=_rays(RayVector{T}, PC), lineality_basis=lineality_space(PC) + NamedTuple{(:rays_modulo_lineality, :lineality_basis), Tuple{SubObjectIterator{RayVector{T}}, SubObjectIterator{RayVector{T}}}}, PC ) + +function rays_modulo_lineality(::Type{NamedTuple{(:rays_modulo_lineality, :lineality_basis), Tuple{SubObjectIterator{RayVector{T}}, SubObjectIterator{RayVector{T}}}}}, PC::PolyhedralComplex{T}) where {T<:scalar_types} + return ( + rays_modulo_lineality = _rays(RayVector{T}, PC), + lineality_basis = lineality_space(PC) + ) end -rays_modulo_lineality( - as::Type{RayVector{T}}, PC::PolyhedralComplex{T} -) where {T<:scalar_types} = _rays(RayVector{T}, PC) + +rays_modulo_lineality(U::Type{RayVector{T}}, PC::PolyhedralComplex{T}) where {T<:scalar_types} = + _rays(U, PC) @doc raw""" minimal_faces(as, PC::PolyhedralComplex) @@ -223,24 +208,19 @@ julia> MFPC.lineality_basis [0, 0, 1] ``` """ -minimal_faces(PC::PolyhedralComplex{T}) where {T<:scalar_types} = minimal_faces( - NamedTuple{ - (:base_points, :lineality_basis), - Tuple{SubObjectIterator{PointVector{T}},SubObjectIterator{RayVector{T}}}, - }, - PC, -) -function minimal_faces( - as::Type{ - NamedTuple{ - (:base_points, :lineality_basis), - Tuple{SubObjectIterator{PointVector{T}},SubObjectIterator{RayVector{T}}}, - }, - }, - PC::PolyhedralComplex{T}, -) where {T<:scalar_types} - return (base_points=_vertices(PointVector{T}, PC), lineality_basis=lineality_space(PC)) +minimal_faces(PC::PolyhedralComplex{T}) where {T<:scalar_types} = + minimal_faces( + NamedTuple{(:base_points, :lineality_basis), Tuple{SubObjectIterator{PointVector{T}}, SubObjectIterator{RayVector{T}}}}, + PC + ) + +function minimal_faces(::Type{NamedTuple{(:base_points, :lineality_basis), Tuple{SubObjectIterator{PointVector{T}}, SubObjectIterator{RayVector{T}}}}}, PC::PolyhedralComplex{T}) where {T<:scalar_types} + return ( + base_points = _vertices(PointVector{T}, PC), + lineality_basis = lineality_space(PC) + ) end + minimal_faces(as::Type{PointVector{T}}, PC::PolyhedralComplex{T}) where {T<:scalar_types} = _vertices(PointVector{T}, PC) @@ -272,12 +252,10 @@ _rays(PC::PolyhedralComplex{T}) where {T<:scalar_types} = _rays(RayVector{T}, PC _ray_indices_polyhedral_complex(PC::Polymake.BigObject) = collect(Polymake.to_one_based_indexing(PC.FAR_VERTICES)) -_ray_polyhedral_complex( - ::Type{RayVector{T}}, PC::PolyhedralComplex, i::Base.Integer -) where {T<:scalar_types} = RayVector{T}( +_ray_polyhedral_complex(U::Type{RayVector{T}}, PC::PolyhedralComplex{T}, i::Base.Integer) where {T<:scalar_types} = ray_vector( coefficient_field(PC), @view pm_object(PC).VERTICES[_ray_indices_polyhedral_complex(pm_object(PC))[i], 2:end] -) +)::U _matrix_for_polymake(::Val{_ray_polyhedral_complex}) = _vector_matrix @@ -286,10 +264,11 @@ _vector_matrix(::Val{_ray_polyhedral_complex}, PC::PolyhedralComplex; homogenize _ray_indices_polyhedral_complex(pm_object(PC)), (homogenized ? 1 : 2):end ] -_maximal_polyhedron( - ::Type{Polyhedron{T}}, PC::PolyhedralComplex, i::Base.Integer -) where {T<:scalar_types} = - Polyhedron{T}(Polymake.fan.polytope(pm_object(PC), i - 1), coefficient_field(PC)) +_maximal_polyhedron(::Type{Polyhedron{T}}, PC::PolyhedralComplex{T}, i::Base.Integer) where {T<:scalar_types} = + Polyhedron{T}( + Polymake.fan.polytope(pm_object(PC), i - 1), + coefficient_field(PC) + ) _vertex_indices(::Val{_maximal_polyhedron}, PC::PolyhedralComplex) = pm_object(PC).MAXIMAL_POLYTOPES[:, _vertex_indices(pm_object(PC))] @@ -479,28 +458,28 @@ end function _ith_polyhedron( ::Type{Polyhedron{T}}, - PC::PolyhedralComplex, + PC::PolyhedralComplex{T}, i::Base.Integer; - f_dim::Int=-1, - f_ind::Vector{Int64}=Vector{Int64}(), -) where {T<:scalar_types} + f_dim::Int = -1, + f_ind::Vector{Int64} = Vector{Int64}() + ) where {T<:scalar_types} pface = Polymake.row(Polymake.fan.cones_of_dim(pm_object(PC), f_dim), f_ind[i]) + V = pm_object(PC).VERTICES[collect(pface), :] + L = pm_object(PC).LINEALITY_SPACE + PT = _scalar_type_to_polymake(T) return Polyhedron{T}( - Polymake.polytope.Polytope{_scalar_type_to_polymake(T)}(; - VERTICES=pm_object(PC).VERTICES[collect(pface), :], - LINEALITY_SPACE=pm_object(PC).LINEALITY_SPACE, - ), - coefficient_field(PC), + Polymake.polytope.Polytope{PT}(VERTICES = V, LINEALITY_SPACE = L), + coefficient_field(PC) ) end lineality_space(PC::PolyhedralComplex{T}) where {T<:scalar_types} = SubObjectIterator{RayVector{T}}(PC, _lineality_complex, lineality_dim(PC)) -_lineality_complex( - ::Type{RayVector{T}}, PC::PolyhedralComplex, i::Base.Integer -) where {T<:scalar_types} = - RayVector{T}(coefficient_field(PC), @view pm_object(PC).LINEALITY_SPACE[i, 2:end]) +_lineality_complex(U::Type{RayVector{T}}, PC::PolyhedralComplex{T}, i::Base.Integer) where {T<:scalar_types} = ray_vector( + coefficient_field(PC), + @view pm_object(PC).LINEALITY_SPACE[i, 2:end] +)::U _generator_matrix(::Val{_lineality_complex}, PC::PolyhedralComplex; homogenized=false) = if homogenized diff --git a/src/PolyhedralGeometry/PolyhedralFan/constructors.jl b/src/PolyhedralGeometry/PolyhedralFan/constructors.jl index eaa018608ea9..0b6a16df06d5 100644 --- a/src/PolyhedralGeometry/PolyhedralFan/constructors.jl +++ b/src/PolyhedralGeometry/PolyhedralFan/constructors.jl @@ -64,33 +64,33 @@ julia> lineality_dim(PF) 1 ``` """ -function polyhedral_fan(f::Union{Type{T}, Field}, +function polyhedral_fan(f::scalar_type_or_field, Rays::AbstractCollection[RayVector], LS::Union{AbstractCollection[RayVector], Nothing}, Incidence::IncidenceMatrix; - non_redundant::Bool = false) where T<:scalar_types - parent_field, scalar_type = _determine_parent_and_scalar(f, Rays, LS) - RM = unhomogenized_matrix(Rays) - if isnothing(LS) - LM = Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(RM, 2)) - else - LM = unhomogenized_matrix(LS) - end - if non_redundant - return PolyhedralFan{scalar_type}(Polymake.fan.PolyhedralFan{_scalar_type_to_polymake(scalar_type)}( - RAYS = RM, - LINEALITY_SPACE = LM, - MAXIMAL_CONES = Incidence, - ), parent_field) - else - return PolyhedralFan{scalar_type}(Polymake.fan.PolyhedralFan{_scalar_type_to_polymake(scalar_type)}( - INPUT_RAYS = RM, - INPUT_LINEALITY = LM, - INPUT_CONES = Incidence, - ), parent_field) - end + non_redundant::Bool = false) + parent_field, scalar_type = _determine_parent_and_scalar(f, Rays, LS) + RM = unhomogenized_matrix(Rays) + if isnothing(LS) + LM = Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(RM, 2)) + else + LM = unhomogenized_matrix(LS) + end + if non_redundant + return PolyhedralFan{scalar_type}(Polymake.fan.PolyhedralFan{_scalar_type_to_polymake(scalar_type)}( + RAYS = RM, + LINEALITY_SPACE = LM, + MAXIMAL_CONES = Incidence, + ), parent_field) + else + return PolyhedralFan{scalar_type}(Polymake.fan.PolyhedralFan{_scalar_type_to_polymake(scalar_type)}( + INPUT_RAYS = RM, + INPUT_LINEALITY = LM, + INPUT_CONES = Incidence, + ), parent_field) + end end -polyhedral_fan(f::Union{Type{T}, Field}, Rays::AbstractCollection[RayVector], Incidence::IncidenceMatrix; non_redundant::Bool = false) where T<:scalar_types = polyhedral_fan(f, Rays, nothing, Incidence; non_redundant=non_redundant) +polyhedral_fan(f::scalar_type_or_field, Rays::AbstractCollection[RayVector], Incidence::IncidenceMatrix; non_redundant::Bool = false) = polyhedral_fan(f, Rays, nothing, Incidence; non_redundant=non_redundant) polyhedral_fan(Rays::AbstractCollection[RayVector], LS::Union{AbstractCollection[RayVector], Nothing}, Incidence::IncidenceMatrix; non_redundant::Bool = false) = polyhedral_fan(QQFieldElem, Rays, LS, Incidence; non_redundant = non_redundant) polyhedral_fan(Rays::AbstractCollection[RayVector], Incidence::IncidenceMatrix; non_redundant::Bool = false) = polyhedral_fan(QQFieldElem, Rays, Incidence; non_redundant = non_redundant) @@ -104,9 +104,9 @@ pm_object(PF::PolyhedralFan) = PF.pm_fan polyhedral_fan(itr::AbstractVector{Cone{T}}) where T<:scalar_types = PolyhedralFan{T}(Polymake.fan.check_fan_objects(pm_object.(itr)...), coefficient_field(iterate(itr)[1])) #Same construction for when the user gives Matrix{Bool} as incidence matrix -polyhedral_fan(f::Union{Type{T}, Field}, Rays::AbstractCollection[RayVector], LS::AbstractCollection[RayVector], Incidence::Matrix{Bool}) where T<:scalar_types = +polyhedral_fan(f::scalar_type_or_field, Rays::AbstractCollection[RayVector], LS::AbstractCollection[RayVector], Incidence::Matrix{Bool}) = polyhedral_fan(f, Rays, LS, IncidenceMatrix(Polymake.IncidenceMatrix(Incidence))) -polyhedral_fan(f::Union{Type{T}, Field}, Rays::AbstractCollection[RayVector], Incidence::Matrix{Bool}) where T<:scalar_types = +polyhedral_fan(f::scalar_type_or_field, Rays::AbstractCollection[RayVector], Incidence::Matrix{Bool}) = polyhedral_fan(f, Rays, IncidenceMatrix(Polymake.IncidenceMatrix(Incidence))) @@ -135,15 +135,15 @@ parent `Field`. group acting on the rays of the fan. """ -function polyhedral_fan_from_rays_action(f::Union{Type{T}, Field}, Rays::AbstractCollection[RayVector], MC_reps::IncidenceMatrix, perms::AbstractVector{PermGroupElem}) where T<:scalar_types - parent_field, scalar_type = _determine_parent_and_scalar(f, Rays) - pf = Polymake.fan.PolyhedralFan{_scalar_type_to_polymake(scalar_type)}() - Polymake.take(pf, "RAYS", Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(unhomogenized_matrix(Rays))) - d = length(Rays) - gp = _group_generators_to_pm_arr_arr(perms, d) - Polymake.take(pf, "GROUP.RAYS_ACTION.GENERATORS", gp) - Polymake.take(pf, "GROUP.MAXIMAL_CONES_ACTION.MAXIMAL_CONES_GENERATORS", MC_reps) - return PolyhedralFan{scalar_type}(pf, parent_field) +function polyhedral_fan_from_rays_action(f::scalar_type_or_field, Rays::AbstractCollection[RayVector], MC_reps::IncidenceMatrix, perms::AbstractVector{PermGroupElem}) + parent_field, scalar_type = _determine_parent_and_scalar(f, Rays) + pf = Polymake.fan.PolyhedralFan{_scalar_type_to_polymake(scalar_type)}() + Polymake.take(pf, "RAYS", Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(unhomogenized_matrix(Rays))) + d = length(Rays) + gp = _group_generators_to_pm_arr_arr(perms, d) + Polymake.take(pf, "GROUP.RAYS_ACTION.GENERATORS", gp) + Polymake.take(pf, "GROUP.MAXIMAL_CONES_ACTION.MAXIMAL_CONES_GENERATORS", MC_reps) + return PolyhedralFan{scalar_type}(pf, parent_field) end polyhedral_fan_from_rays_action(Rays::AbstractCollection[RayVector], MC_reps::IncidenceMatrix, perms::AbstractVector{PermGroupElem}) = polyhedral_fan_from_rays_action(QQFieldElem, Rays, MC_reps, perms) diff --git a/src/PolyhedralGeometry/PolyhedralFan/properties.jl b/src/PolyhedralGeometry/PolyhedralFan/properties.jl index 05bf4beb8e19..191c1a884a8e 100644 --- a/src/PolyhedralGeometry/PolyhedralFan/properties.jl +++ b/src/PolyhedralGeometry/PolyhedralFan/properties.jl @@ -46,13 +46,15 @@ julia> matrix(QQ, rays(NF)) rays(PF::_FanLikeType) = lineality_dim(PF) == 0 ? _rays(PF) : _empty_subobjectiterator(RayVector{_get_scalar_type(PF)}, PF) _rays(PF::_FanLikeType) = SubObjectIterator{RayVector{_get_scalar_type(PF)}}(PF, _ray_fan, _nrays(PF)) -_ray_fan(::Type{RayVector{T}}, PF::_FanLikeType, i::Base.Integer) where T<:scalar_types = RayVector{T}(coefficient_field(PF), view(pm_object(PF).RAYS, i, :)) +_ray_fan(U::Type{RayVector{T}}, PF::_FanLikeType, i::Base.Integer) where {T<:scalar_types} = + ray_vector(coefficient_field(PF), view(pm_object(PF).RAYS, i, :))::U _vector_matrix(::Val{_ray_fan}, PF::_FanLikeType; homogenized=false) = homogenized ? homogenize(pm_object(PF).RAYS, 0) : pm_object(PF).RAYS _matrix_for_polymake(::Val{_ray_fan}) = _vector_matrix -_maximal_cone(::Type{Cone{T}}, PF::_FanLikeType, i::Base.Integer) where T<:scalar_types = Cone{T}(Polymake.fan.cone(pm_object(PF), i - 1), coefficient_field(PF)) +_maximal_cone(::Type{Cone{T}}, PF::_FanLikeType, i::Base.Integer) where {T<:scalar_types} = + Cone{T}(Polymake.fan.cone(pm_object(PF), i - 1), coefficient_field(PF)) @doc raw""" @@ -87,14 +89,20 @@ julia> rays(NF) 0-element SubObjectIterator{RayVector{QQFieldElem}} ``` """ -rays_modulo_lineality(PF::_FanLikeType) = rays_modulo_lineality(NamedTuple{(:rays_modulo_lineality, :lineality_basis), Tuple{SubObjectIterator{RayVector{_get_scalar_type(PF)}}, SubObjectIterator{RayVector{_get_scalar_type(PF)}}}}, PF) -function rays_modulo_lineality(as::Type{NamedTuple{(:rays_modulo_lineality, :lineality_basis), Tuple{SubObjectIterator{RayVector{T}}, SubObjectIterator{RayVector{T}}}}}, F::_FanLikeType) where T<:scalar_types +rays_modulo_lineality(F::_FanLikeType) = + rays_modulo_lineality( + NamedTuple{(:rays_modulo_lineality, :lineality_basis), Tuple{SubObjectIterator{RayVector{_get_scalar_type(F)}}, SubObjectIterator{RayVector{_get_scalar_type(F)}}}}, + F + ) + +function rays_modulo_lineality(::Type{NamedTuple{(:rays_modulo_lineality, :lineality_basis), Tuple{SubObjectIterator{RayVector{T}}, SubObjectIterator{RayVector{T}}}}}, F::_FanLikeType) where {T<:scalar_types} return ( rays_modulo_lineality = _rays(F), lineality_basis = lineality_space(F) ) end -rays_modulo_lineality(as::Type{RayVector}, F::_FanLikeType) = _rays(F) + +rays_modulo_lineality(::Type{<:RayVector}, F::_FanLikeType) = _rays(F) @@ -159,7 +167,10 @@ function cones(PF::_FanLikeType, cone_dim::Int) end function _cone_of_dim(::Type{Cone{T}}, PF::_FanLikeType, i::Base.Integer; c_dim::Int = 0) where T<:scalar_types - return Cone{T}(Polymake.polytope.Cone{_scalar_type_to_polymake(T)}(RAYS = pm_object(PF).RAYS[collect(Polymake.row(Polymake.fan.cones_of_dim(pm_object(PF), c_dim), i)),:], LINEALITY_SPACE = pm_object(PF).LINEALITY_SPACE), coefficient_field(PF)) + R = pm_object(PF).RAYS[collect(Polymake.row(Polymake.fan.cones_of_dim(pm_object(PF), c_dim), i)), :] + L = pm_object(PF).LINEALITY_SPACE + PT = _scalar_type_to_polymake(T) + return Cone{T}(Polymake.polytope.Cone{PT}(RAYS = R, LINEALITY_SPACE = L), coefficient_field(PF)) end _ray_indices(::Val{_cone_of_dim}, PF::_FanLikeType; c_dim::Int = 0) = Polymake.fan.cones_of_dim(pm_object(PF), c_dim) @@ -385,7 +396,7 @@ julia> lineality_space(PF) """ lineality_space(PF::_FanLikeType) = SubObjectIterator{RayVector{_get_scalar_type(PF)}}(PF, _lineality_fan, lineality_dim(PF)) -_lineality_fan(::Type{RayVector{T}}, PF::_FanLikeType, i::Base.Integer) where T<:scalar_types = RayVector{T}(coefficient_field(PF), view(pm_object(PF).LINEALITY_SPACE, i, :)) +_lineality_fan(U::Type{RayVector{T}}, PF::_FanLikeType, i::Base.Integer) where {T<:scalar_types} = ray_vector(coefficient_field(PF), view(pm_object(PF).LINEALITY_SPACE, i, :))::U _generator_matrix(::Val{_lineality_fan}, PF::_FanLikeType; homogenized=false) = homogenized ? homogenize(pm_object(PF).LINEALITY_SPACE, 0) : pm_object(PF).LINEALITY_SPACE diff --git a/src/PolyhedralGeometry/Polyhedron/constructors.jl b/src/PolyhedralGeometry/Polyhedron/constructors.jl index 45250c090245..d5347eeb6780 100644 --- a/src/PolyhedralGeometry/Polyhedron/constructors.jl +++ b/src/PolyhedralGeometry/Polyhedron/constructors.jl @@ -30,8 +30,11 @@ polyhedron(A) = polyhedron(QQFieldElem, A) Construct a `Polyhedron` corresponding to a `Polymake.BigObject` of type `Polytope`. Scalar type and parent field will be detected automatically. To improve type stability and performance, please use [`Polyhedron{T}(p::Polymake.BigObject, f::Field) where T<:scalar_types`](@ref) instead, where possible. """ function polyhedron(p::Polymake.BigObject) - T, f = _detect_scalar_and_field(Polyhedron, p) - return Polyhedron{T}(p, f) + T, f = _detect_scalar_and_field(Polyhedron, p) + if T == EmbeddedElem{nf_elem} && Hecke.isquadratic_type(number_field(f))[1] && Polymake.bigobject_eltype(p) == "QuadraticExtension" + p = _polyhedron_qe_to_on(p, f) + end + return Polyhedron{T}(p, f) end @doc raw""" @@ -57,17 +60,17 @@ julia> polyhedron(A,b) Polyhedron in ambient dimension 2 ``` """ -polyhedron(f::Union{Type{T}, Field}, A::AnyVecOrMat, b::AbstractVector) where T<:scalar_types = polyhedron(f, (A, b)) +polyhedron(f::scalar_type_or_field, A::AnyVecOrMat, b::AbstractVector; non_redundant::Bool = false) = polyhedron(f, (A, b); non_redundant = non_redundant) -polyhedron(f::Union{Type{T}, Field}, A::AbstractVector, b::Any) where T<:scalar_types = polyhedron(f, ([A], [b])) +polyhedron(f::scalar_type_or_field, A::AbstractVector, b::Any; non_redundant::Bool = false) = polyhedron(f, ([A], [b]); non_redundant = non_redundant) -polyhedron(f::Union{Type{T}, Field}, A::AbstractVector, b::AbstractVector) where T<:scalar_types = polyhedron(f, ([A], b)) +polyhedron(f::scalar_type_or_field, A::AbstractVector, b::AbstractVector; non_redundant::Bool = false) = polyhedron(f, ([A], b); non_redundant = non_redundant) -polyhedron(f::Union{Type{T}, Field}, A::AbstractVector{<:AbstractVector}, b::Any) where T<:scalar_types = polyhedron(f, (A, [b])) +polyhedron(f::scalar_type_or_field, A::AbstractVector{<:AbstractVector}, b::Any; non_redundant::Bool = false) = polyhedron(f, (A, [b]); non_redundant = non_redundant) -polyhedron(f::Union{Type{T}, Field}, A::AbstractVector{<:AbstractVector}, b::AbstractVector) where T<:scalar_types = polyhedron(f, (A, b)) +polyhedron(f::scalar_type_or_field, A::AbstractVector{<:AbstractVector}, b::AbstractVector; non_redundant::Bool = false) = polyhedron(f, (A, b); non_redundant = non_redundant) -polyhedron(f::Union{Type{T}, Field}, A::AnyVecOrMat, b::Any) where T<:scalar_types = polyhedron(f, A, [b]) +polyhedron(f::scalar_type_or_field, A::AnyVecOrMat, b::Any; non_redundant::Bool = false) = polyhedron(f, A, [b]; non_redundant = non_redundant) @doc raw""" polyhedron(::Union{Type{T}, Field}, I::Union{Nothing, AbstractCollection[AffineHalfspace]}, E::Union{Nothing, AbstractCollection[AffineHyperplane]} = nothing) where T<:scalar_types @@ -106,17 +109,21 @@ julia> vertices(P) [0, 0] ``` """ -function polyhedron(f::Union{Type{T}, Field}, I::Union{Nothing, AbstractCollection[AffineHalfspace]}, E::Union{Nothing, AbstractCollection[AffineHyperplane]} = nothing; parent_field::Union{Nothing, Field} = nothing) where T<:scalar_types - parent_field, scalar_type = _determine_parent_and_scalar(f, I, E) - if isnothing(I) || _isempty_halfspace(I) - EM = affine_matrix_for_polymake(E) - IM = Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(EM, 2)) - else - IM = -affine_matrix_for_polymake(I) - EM = isnothing(E) || _isempty_halfspace(E) ? Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(IM, 2)) : affine_matrix_for_polymake(E) - end - +function polyhedron(f::scalar_type_or_field, I::Union{Nothing, AbstractCollection[AffineHalfspace]}, E::Union{Nothing, AbstractCollection[AffineHyperplane]} = nothing; non_redundant::Bool = false) + parent_field, scalar_type = _determine_parent_and_scalar(f, I, E) + if isnothing(I) || _isempty_halfspace(I) + EM = affine_matrix_for_polymake(E) + IM = Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(EM, 2)) + else + IM = -affine_matrix_for_polymake(I) + EM = isnothing(E) || _isempty_halfspace(E) ? Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(IM, 2)) : affine_matrix_for_polymake(E) + end + + if non_redundant + return Polyhedron{scalar_type}(Polymake.polytope.Polytope{_scalar_type_to_polymake(scalar_type)}(FACETS = remove_zero_rows(IM), AFFINE_HULL = remove_zero_rows(EM)), parent_field) + else return Polyhedron{scalar_type}(Polymake.polytope.Polytope{_scalar_type_to_polymake(scalar_type)}(INEQUALITIES = remove_zero_rows(IM), EQUATIONS = remove_zero_rows(EM)), parent_field) + end end """ @@ -126,15 +133,10 @@ Get the underlying polymake `Polytope`. """ pm_object(P::Polyhedron) = P.pm_polytope -function ==(P0::Polyhedron, P1::Polyhedron) - # TODO: Remove the following 3 lines, see #758 - for pair in Iterators.product([P0, P1], ["RAYS", "FACETS"]) - Polymake.give(pm_object(pair[1]),pair[2]) - end +function ==(P0::Polyhedron{T}, P1::Polyhedron{T}) where T<:scalar_types Polymake.polytope.equal_polyhedra(pm_object(P0), pm_object(P1)) end - ### Construct polyhedron from V-data, as the convex hull of points, rays and lineality. @doc raw""" convex_hull([::Union{Type{T}, Field} = QQFieldElem,] V [, R [, L]]; non_redundant::Bool = false) @@ -198,20 +200,20 @@ julia> XA = convex_hull(V, R, L) Polyhedron in ambient dimension 2 ``` """ -function convex_hull(f::Union{Type{T}, Field}, V::AbstractCollection[PointVector], R::Union{AbstractCollection[RayVector], Nothing} = nothing, L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) where T<:scalar_types - parent_field, scalar_type = _determine_parent_and_scalar(f, V, R, L) - # Rays and Points are homogenized and combined and - # Lineality is homogenized - points = stack(homogenized_matrix(V, 1), homogenized_matrix(R, 0)) - lineality = isnothing(L) || isempty(L) ? zero_matrix(QQ, 0, size(points,2)) : homogenized_matrix(L, 0) - - # These matrices are in the right format for polymake. - # given non_redundant can avoid unnecessary redundancy checks - if non_redundant - return Polyhedron{scalar_type}(Polymake.polytope.Polytope{_scalar_type_to_polymake(scalar_type)}(VERTICES = points, LINEALITY_SPACE = lineality), parent_field) - else - return Polyhedron{scalar_type}(Polymake.polytope.Polytope{_scalar_type_to_polymake(scalar_type)}(POINTS = remove_zero_rows(points), INPUT_LINEALITY = remove_zero_rows(lineality)), parent_field) - end +function convex_hull(f::scalar_type_or_field, V::AbstractCollection[PointVector], R::Union{AbstractCollection[RayVector], Nothing} = nothing, L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) + parent_field, scalar_type = _determine_parent_and_scalar(f, V, R, L) + # Rays and Points are homogenized and combined and + # Lineality is homogenized + points = stack(homogenized_matrix(V, 1), homogenized_matrix(R, 0)) + lineality = isnothing(L) || isempty(L) ? zero_matrix(QQ, 0, size(points,2)) : homogenized_matrix(L, 0) + + # These matrices are in the right format for polymake. + # given non_redundant can avoid unnecessary redundancy checks + if non_redundant + return Polyhedron{scalar_type}(Polymake.polytope.Polytope{_scalar_type_to_polymake(scalar_type)}(VERTICES = points, LINEALITY_SPACE = lineality), parent_field) + else + return Polyhedron{scalar_type}(Polymake.polytope.Polytope{_scalar_type_to_polymake(scalar_type)}(POINTS = remove_zero_rows(points), INPUT_LINEALITY = remove_zero_rows(lineality)), parent_field) + end end convex_hull(V::AbstractCollection[PointVector], R::Union{AbstractCollection[RayVector], Nothing} = nothing, L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) = convex_hull(QQFieldElem, V, R, L; non_redundant=non_redundant) diff --git a/src/PolyhedralGeometry/Polyhedron/properties.jl b/src/PolyhedralGeometry/Polyhedron/properties.jl index 33330dc1ddca..e2b074d9dc69 100644 --- a/src/PolyhedralGeometry/Polyhedron/properties.jl +++ b/src/PolyhedralGeometry/Polyhedron/properties.jl @@ -35,9 +35,12 @@ function faces(P::Polyhedron{T}, face_dim::Int) where T<:scalar_types return SubObjectIterator{Polyhedron{T}}(P, _face_polyhedron, length(rfaces), (f_dim = n, f_ind = rfaces)) end -function _face_polyhedron(::Type{Polyhedron{T}}, P::Polyhedron, i::Base.Integer; f_dim::Int = -1, f_ind::Vector{Int64} = Vector{Int64}()) where T<:scalar_types - pface = Polymake.to_one_based_indexing(Polymake.polytope.faces_of_dim(pm_object(P), f_dim))[f_ind[i]] - return Polyhedron{T}(Polymake.polytope.Polytope{_scalar_type_to_polymake(T)}(VERTICES = pm_object(P).VERTICES[collect(pface),:], LINEALITY_SPACE = pm_object(P).LINEALITY_SPACE), coefficient_field(P)) +function _face_polyhedron(::Type{Polyhedron{T}}, P::Polyhedron{T}, i::Base.Integer; f_dim::Int = -1, f_ind::Vector{Int64} = Vector{Int64}()) where T<:scalar_types + pface = Polymake.to_one_based_indexing(Polymake.polytope.faces_of_dim(pm_object(P), f_dim))[f_ind[i]] + V = pm_object(P).VERTICES[collect(pface), :] + L = pm_object(P).LINEALITY_SPACE + PT = _scalar_type_to_polymake(T) + return Polyhedron{T}(Polymake.polytope.Polytope{PT}(VERTICES = V, LINEALITY_SPACE = L), coefficient_field(P)) end function _vertex_indices(::Val{_face_polyhedron}, P::Polyhedron; f_dim = -1, f_ind::Vector{Int64} = Vector{Int64}()) @@ -54,9 +57,12 @@ end _incidencematrix(::Val{_face_polyhedron}) = _vertex_and_ray_indices -function _face_polyhedron_facet(::Type{Polyhedron{T}}, P::Polyhedron, i::Base.Integer) where T<:scalar_types - pface = pm_object(P).VERTICES_IN_FACETS[_facet_index(pm_object(P), i), :] - return Polyhedron{T}(Polymake.polytope.Polytope{_scalar_type_to_polymake(T)}(VERTICES = pm_object(P).VERTICES[collect(pface),:], LINEALITY_SPACE = pm_object(P).LINEALITY_SPACE), coefficient_field(P)) +function _face_polyhedron_facet(::Type{Polyhedron{T}}, P::Polyhedron{T}, i::Base.Integer) where T<:scalar_types + pface = pm_object(P).VERTICES_IN_FACETS[_facet_index(pm_object(P), i), :] + V = pm_object(P).VERTICES[collect(pface), :] + L = pm_object(P).LINEALITY_SPACE + PT = _scalar_type_to_polymake(T) + return Polyhedron{T}(Polymake.polytope.Polytope{PT}(VERTICES = V, LINEALITY_SPACE = L), coefficient_field(P)) end _vertex_indices(::Val{_face_polyhedron_facet}, P::Polyhedron) = vcat(pm_object(P).VERTICES_IN_FACETS[1:(_facet_at_infinity(pm_object(P)) - 1), _vertex_indices(pm_object(P))], pm_object(P).VERTICES_IN_FACETS[(_facet_at_infinity(pm_object(P)) + 1):end, _vertex_indices(pm_object(P))]) @@ -127,13 +133,13 @@ julia> minimal_faces(P) ``` """ minimal_faces(P::Polyhedron{T}) where T<:scalar_types = minimal_faces(NamedTuple{(:base_points, :lineality_basis), Tuple{SubObjectIterator{PointVector{T}}, SubObjectIterator{RayVector{T}}}}, P) -function minimal_faces(as::Type{NamedTuple{(:base_points, :lineality_basis), Tuple{SubObjectIterator{PointVector{T}}, SubObjectIterator{RayVector{T}}}}}, P::Polyhedron{T}) where T<:scalar_types - return ( - base_points = _vertices(PointVector{T}, P), - lineality_basis = lineality_space(P) - ) +function minimal_faces(::Type{NamedTuple{(:base_points, :lineality_basis), Tuple{SubObjectIterator{PointVector{T}}, SubObjectIterator{RayVector{T}}}}}, P::Polyhedron{T}) where T<:scalar_types + return ( + base_points = _vertices(PointVector{T}, P), + lineality_basis = lineality_space(P) + ) end -minimal_faces(as::Type{PointVector{T}}, P::Polyhedron{T}) where T<:scalar_types = _vertices(PointVector{T}, P) +minimal_faces(::Type{<:PointVector}, P::Polyhedron) = _vertices(P) @@ -163,13 +169,13 @@ julia> rmlP.lineality_basis ``` """ rays_modulo_lineality(P::Polyhedron{T}) where T<:scalar_types = rays_modulo_lineality(NamedTuple{(:rays_modulo_lineality, :lineality_basis), Tuple{SubObjectIterator{RayVector{T}}, SubObjectIterator{RayVector{T}}}}, P) -function rays_modulo_lineality(as::Type{NamedTuple{(:rays_modulo_lineality, :lineality_basis), Tuple{SubObjectIterator{RayVector{T}}, SubObjectIterator{RayVector{T}}}}}, P::Polyhedron) where T<:scalar_types - return ( - rays_modulo_lineality = _rays(P), - lineality_basis = lineality_space(P) - ) +function rays_modulo_lineality(::Type{NamedTuple{(:rays_modulo_lineality, :lineality_basis), Tuple{SubObjectIterator{RayVector{T}}, SubObjectIterator{RayVector{T}}}}}, P::Polyhedron{T}) where T<:scalar_types + return ( + rays_modulo_lineality = _rays(P), + lineality_basis = lineality_space(P) + ) end -rays_modulo_lineality(as::Type{RayVector}, P::Polyhedron) = _rays(P) +rays_modulo_lineality(::Type{<:RayVector}, P::Polyhedron) = _rays(P) @doc raw""" @@ -197,17 +203,17 @@ julia> vertices(PointVector, P) [1, 2] ``` """ -vertices(as::Type{PointVector{T}}, P::Polyhedron) where T<:scalar_types = lineality_dim(P) == 0 ? _vertices(as, P) : _empty_subobjectiterator(as, P) -_vertices(as::Type{PointVector{T}}, P::Polyhedron) where T<:scalar_types = SubObjectIterator{as}(P, _vertex_polyhedron, length(_vertex_indices(pm_object(P)))) +vertices(as::Type{PointVector{T}}, P::Polyhedron{T}) where T<:scalar_types = lineality_dim(P) == 0 ? _vertices(as, P) : _empty_subobjectiterator(as, P) +_vertices(as::Type{PointVector{T}}, P::Polyhedron{T}) where T<:scalar_types = SubObjectIterator{as}(P, _vertex_polyhedron, length(_vertex_indices(pm_object(P)))) -_vertex_polyhedron(::Type{PointVector{T}}, P::Polyhedron, i::Base.Integer) where T<:scalar_types = PointVector{T}(coefficient_field(P).(@view pm_object(P).VERTICES[_vertex_indices(pm_object(P))[i], 2:end])) +_vertex_polyhedron(U::Type{PointVector{T}}, P::Polyhedron{T}, i::Base.Integer) where T<:scalar_types = point_vector(coefficient_field(P), @view pm_object(P).VERTICES[_vertex_indices(pm_object(P))[i], 2:end])::U _point_matrix(::Val{_vertex_polyhedron}, P::Polyhedron; homogenized=false) = @view pm_object(P).VERTICES[_vertex_indices(pm_object(P)), (homogenized ? 1 : 2):end] _matrix_for_polymake(::Val{_vertex_polyhedron}) = _point_matrix -vertices(::Type{PointVector}, P::Polyhedron{T}) where T<:scalar_types = vertices(PointVector{T}, P) -_vertices(::Type{PointVector}, P::Polyhedron{T}) where T<:scalar_types = _vertices(PointVector{T}, P) +vertices(::Type{<:PointVector}, P::Polyhedron{T}) where T<:scalar_types = vertices(PointVector{T}, P) +_vertices(::Type{<:PointVector}, P::Polyhedron{T}) where T<:scalar_types = _vertices(PointVector{T}, P) _facet_indices(::Val{_vertex_polyhedron}, P::Polyhedron) = pm_object(P).FACETS_THRU_VERTICES[_vertex_indices(pm_object(P)),_facet_indices(pm_object(P))] @@ -311,10 +317,10 @@ julia> rays(RayVector, PO) [0, 1] ``` """ -rays(as::Type{RayVector{T}}, P::Polyhedron) where T<:scalar_types = lineality_dim(P) == 0 ? _rays(as, P) : _empty_subobjectiterator(as, P) -_rays(as::Type{RayVector{T}}, P::Polyhedron) where T<:scalar_types = SubObjectIterator{as}(P, _ray_polyhedron, length(_ray_indices(pm_object(P)))) +rays(as::Type{RayVector{T}}, P::Polyhedron{T}) where T<:scalar_types = lineality_dim(P) == 0 ? _rays(as, P) : _empty_subobjectiterator(as, P) +_rays(as::Type{RayVector{T}}, P::Polyhedron{T}) where T<:scalar_types = SubObjectIterator{as}(P, _ray_polyhedron, length(_ray_indices(pm_object(P)))) -_ray_polyhedron(::Type{RayVector{T}}, P::Polyhedron, i::Base.Integer) where T<:scalar_types = RayVector{T}(coefficient_field(P).(@view pm_object(P).VERTICES[_ray_indices(pm_object(P))[i], 2:end])) +_ray_polyhedron(U::Type{RayVector{T}}, P::Polyhedron{T}, i::Base.Integer) where T<:scalar_types = ray_vector(coefficient_field(P), @view pm_object(P).VERTICES[_ray_indices(pm_object(P))[i], 2:end])::U _facet_indices(::Val{_ray_polyhedron}, P::Polyhedron) = pm_object(P).FACETS_THRU_RAYS[_ray_indices(pm_object(P)), _facet_indices(pm_object(P))] @@ -324,8 +330,8 @@ _matrix_for_polymake(::Val{_ray_polyhedron}) = _vector_matrix _incidencematrix(::Val{_ray_polyhedron}) = _facet_indices -rays(::Type{RayVector}, P::Polyhedron{T}) where T<:scalar_types = rays(RayVector{T}, P) -_rays(::Type{RayVector}, P::Polyhedron{T}) where T<:scalar_types = _rays(RayVector{T}, P) +rays(::Type{<:RayVector}, P::Polyhedron{T}) where T<:scalar_types = rays(RayVector{T}, P) +_rays(::Type{<:RayVector}, P::Polyhedron{T}) where T<:scalar_types = _rays(RayVector{T}, P) @doc raw""" rays(P::Polyhedron) @@ -410,17 +416,17 @@ x₃ ≦ 1 """ facets(as::Type{T}, P::Polyhedron{S}) where {R, S<:scalar_types, T<:Union{AffineHalfspace{S}, Pair{R, S}, Polyhedron{S}}} = SubObjectIterator{as}(P, _facet_polyhedron, nfacets(P)) -function _facet_polyhedron(::Type{AffineHalfspace{S}}, P::Polyhedron, i::Base.Integer) where S<:scalar_types - h = decompose_hdata(view(pm_object(P).FACETS, [_facet_index(pm_object(P), i)], :)) - return affine_halfspace(coefficient_field(P), h[1], h[2][]) +function _facet_polyhedron(U::Type{AffineHalfspace{S}}, P::Polyhedron{S}, i::Base.Integer) where S<:scalar_types + h = decompose_hdata(view(pm_object(P).FACETS, [_facet_index(pm_object(P), i)], :)) + return affine_halfspace(coefficient_field(P), h[1], h[2][])::U end -function _facet_polyhedron(::Type{Pair{R, S}}, P::Polyhedron, i::Base.Integer) where {R, S<:scalar_types} - f = coefficient_field(P) - h = decompose_hdata(view(pm_object(P).FACETS, [_facet_index(pm_object(P), i)], :)) - return Pair{R, S}(f.(view(h[1], :, :)), f(h[2][])) # view_broadcast +function _facet_polyhedron(U::Type{Pair{R, S}}, P::Polyhedron{S}, i::Base.Integer) where {R, S<:scalar_types} + f = coefficient_field(P) + h = decompose_hdata(view(pm_object(P).FACETS, [_facet_index(pm_object(P), i)], :)) + return U(f.(view(h[1], :, :)), f(h[2][])) # view_broadcast end -function _facet_polyhedron(::Type{Polyhedron{T}}, P::Polyhedron, i::Base.Integer) where T<:scalar_types - return Polyhedron{T}(Polymake.polytope.facet(pm_object(P), _facet_index(pm_object(P), i)-1), coefficient_field(P)) +function _facet_polyhedron(::Type{Polyhedron{T}}, P::Polyhedron{T}, i::Base.Integer) where T<:scalar_types + return Polyhedron{T}(Polymake.polytope.facet(pm_object(P), _facet_index(pm_object(P), i) - 1), coefficient_field(P)) end _affine_inequality_matrix(::Val{_facet_polyhedron}, P::Polyhedron) = -_remove_facet_at_infinity(pm_object(P)) @@ -435,7 +441,7 @@ _vertex_and_ray_indices(::Val{_facet_polyhedron}, P::Polyhedron) = vcat(pm_objec _incidencematrix(::Val{_facet_polyhedron}) = _vertex_and_ray_indices -facets(::Type{Pair}, P::Polyhedron{T}) where T<:scalar_types = facets(Pair{Matrix{T}, T}, P) +facets(::Type{<:Pair}, P::Polyhedron{T}) where T<:scalar_types = facets(Pair{Matrix{T}, T}, P) facets(::Type{Polyhedron}, P::Polyhedron{T}) where T<:scalar_types = facets(Polyhedron{T}, P) @@ -461,9 +467,7 @@ x₃ ≦ 1 """ facets(P::Polyhedron{T}) where T<:scalar_types = facets(AffineHalfspace{T}, P) -facets(::Type{Halfspace}, P::Polyhedron{T}) where T<:scalar_types = facets(AffineHalfspace{T}, P) - -facets(::Type{AffineHalfspace}, P::Polyhedron{T}) where T<:scalar_types = facets(AffineHalfspace{T}, P) +facets(::Type{<:Halfspace}, P::Polyhedron{T}) where T<:scalar_types = facets(AffineHalfspace{T}, P) function _facet_index(P::Polymake.BigObject, i::Base.Integer) i < _facet_at_infinity(P) && return i @@ -562,7 +566,7 @@ julia> normalized_volume(C) 8 ``` """ -normalized_volume(P::Polyhedron{T}) where T<:scalar_types = coefficient_field(P)(factorial(dim(P))*(pm_object(P)).VOLUME) +normalized_volume(P::Polyhedron) = coefficient_field(P)(factorial(dim(P))*(pm_object(P)).VOLUME) @doc raw""" @@ -615,7 +619,7 @@ function lattice_points(P::Polyhedron{QQFieldElem}) return SubObjectIterator{PointVector{ZZRingElem}}(P, _lattice_point, size(pm_object(P).LATTICE_POINTS_GENERATORS[1], 1)) end -_lattice_point(::Type{PointVector{ZZRingElem}}, P::Polyhedron, i::Base.Integer) = PointVector{ZZRingElem}(@view pm_object(P).LATTICE_POINTS_GENERATORS[1][i, 2:end]) +_lattice_point(T::Type{PointVector{ZZRingElem}}, P::Polyhedron{QQFieldElem}, i::Base.Integer) = point_vector(ZZ, @view pm_object(P).LATTICE_POINTS_GENERATORS[1][i, 2:end])::T _point_matrix(::Val{_lattice_point}, P::Polyhedron; homogenized=false) = @view pm_object(P).LATTICE_POINTS_GENERATORS[1][:, (homogenized ? 1 : 2):end] @@ -646,7 +650,7 @@ function interior_lattice_points(P::Polyhedron{QQFieldElem}) return SubObjectIterator{PointVector{ZZRingElem}}(P, _interior_lattice_point, size(pm_object(P).INTERIOR_LATTICE_POINTS, 1)) end -_interior_lattice_point(::Type{PointVector{ZZRingElem}}, P::Polyhedron, i::Base.Integer) = PointVector{ZZRingElem}(@view pm_object(P).INTERIOR_LATTICE_POINTS[i, 2:end]) +_interior_lattice_point(T::Type{PointVector{ZZRingElem}}, P::Polyhedron{QQFieldElem}, i::Base.Integer) = point_vector(ZZ, @view pm_object(P).INTERIOR_LATTICE_POINTS[i, 2:end])::T _point_matrix(::Val{_interior_lattice_point}, P::Polyhedron; homogenized=false) = homogenized ? pm_object(P).INTERIOR_LATTICE_POINTS : @view pm_object(P).INTERIOR_LATTICE_POINTS[:, 2:end] @@ -686,7 +690,7 @@ function boundary_lattice_points(P::Polyhedron{QQFieldElem}) return SubObjectIterator{PointVector{ZZRingElem}}(P, _boundary_lattice_point, size(pm_object(P).BOUNDARY_LATTICE_POINTS, 1)) end -_boundary_lattice_point(::Type{PointVector{ZZRingElem}}, P::Polyhedron, i::Base.Integer) = PointVector{ZZRingElem}(@view pm_object(P).BOUNDARY_LATTICE_POINTS[i, 2:end]) +_boundary_lattice_point(T::Type{PointVector{ZZRingElem}}, P::Polyhedron{QQFieldElem}, i::Base.Integer) = point_vector(ZZ, @view pm_object(P).BOUNDARY_LATTICE_POINTS[i, 2:end])::T _point_matrix(::Val{_boundary_lattice_point}, P::Polyhedron; homogenized=false) = homogenized ? pm_object(P).BOUNDARY_LATTICE_POINTS : @view pm_object(P).BOUNDARY_LATTICE_POINTS[:, 2:end] @@ -753,7 +757,7 @@ julia> lineality_space(UH) """ lineality_space(P::Polyhedron{T}) where T<:scalar_types = SubObjectIterator{RayVector{T}}(P, _lineality_polyhedron, lineality_dim(P)) -_lineality_polyhedron(::Type{RayVector{T}}, P::Polyhedron, i::Base.Integer) where T<:scalar_types = RayVector{T}(coefficient_field(P).(@view pm_object(P).LINEALITY_SPACE[i, 2:end])) +_lineality_polyhedron(U::Type{RayVector{T}}, P::Polyhedron{T}, i::Base.Integer) where T<:scalar_types = ray_vector(coefficient_field(P), @view pm_object(P).LINEALITY_SPACE[i, 2:end])::U _generator_matrix(::Val{_lineality_polyhedron}, P::Polyhedron; homogenized=false) = homogenized ? pm_object(P).LINEALITY_SPACE : @view pm_object(P).LINEALITY_SPACE[:, 2:end] @@ -779,9 +783,9 @@ x₄ = 5 """ affine_hull(P::Polyhedron{T}) where T<:scalar_types = SubObjectIterator{AffineHyperplane{T}}(P, _affine_hull, size(pm_object(P).AFFINE_HULL, 1)) -function _affine_hull(::Type{AffineHyperplane{T}}, P::Polyhedron, i::Base.Integer) where T<:scalar_types - h = decompose_hdata(-view(pm_object(P).AFFINE_HULL, [i], :)) - return affine_hyperplane(coefficient_field(P), h[1], h[2][]) +function _affine_hull(U::Type{AffineHyperplane{T}}, P::Polyhedron{T}, i::Base.Integer) where T<:scalar_types + h = decompose_hdata(-view(pm_object(P).AFFINE_HULL, [i], :)) + return affine_hyperplane(coefficient_field(P), h[1], h[2][])::U end _affine_equation_matrix(::Val{_affine_hull}, P::Polyhedron) = pm_object(P).AFFINE_HULL @@ -1010,7 +1014,7 @@ julia> [1, -2] in PO false ``` """ -Base.in(v::AbstractVector, P::Polyhedron) = Polymake.polytope.contains(pm_object(P), [1; v])::Bool +Base.in(v::AbstractVector, P::Polyhedron) = Polymake.polytope.contains(pm_object(P), coefficient_field(P).([1; v]))::Bool @doc raw""" @@ -1213,7 +1217,7 @@ julia> matrix(QQ, vertices(square)) [ 1 1] ``` """ -relative_interior_point(P::Polyhedron{T}) where T<:scalar_types = PointVector{T}(coefficient_field(P).(view(dehomogenize(Polymake.common.dense(pm_object(P).REL_INT_POINT)), :))) #view_broadcast +relative_interior_point(P::Polyhedron{T}) where T<:scalar_types = point_vector(coefficient_field(P), view(dehomogenize(Polymake.common.dense(pm_object(P).REL_INT_POINT)), :))::PointVector{T} #view_broadcast @doc raw""" diff --git a/src/PolyhedralGeometry/Polyhedron/standard_constructions.jl b/src/PolyhedralGeometry/Polyhedron/standard_constructions.jl index 3dbab25fb135..e5edf65c7279 100644 --- a/src/PolyhedralGeometry/Polyhedron/standard_constructions.jl +++ b/src/PolyhedralGeometry/Polyhedron/standard_constructions.jl @@ -206,14 +206,14 @@ julia> normalized_volume(C) 120 ``` """ -function cube(f::Union{Type{T},Field}, d::Int) where {T<:scalar_types} +function cube(f::scalar_type_or_field, d::Int) parent_field, scalar_type = _determine_parent_and_scalar(f) return Polyhedron{scalar_type}( Polymake.polytope.cube{_scalar_type_to_polymake(scalar_type)}(d), parent_field ) end cube(d::Int) = cube(QQFieldElem, d) -function cube(f::Union{Type{T},Field}, d::Int, l, u) where {T<:scalar_types} +function cube(f::scalar_type_or_field, d::Int, l, u) parent_field, scalar_type = _determine_parent_and_scalar(f, l, u) return Polyhedron{scalar_type}( Polymake.polytope.cube{_scalar_type_to_polymake(scalar_type)}(d, u, l), parent_field @@ -622,14 +622,14 @@ julia> facets(t) x₁ + x₂ + x₃ + x₄ + x₅ + x₆ + x₇ ≦ 5 ``` """ -function simplex(f::Union{Type{T},Field}, d::Int, n) where {T<:scalar_types} +function simplex(f::scalar_type_or_field, d::Int, n) parent_field, scalar_type = _determine_parent_and_scalar(f, n) return Polyhedron{scalar_type}( Polymake.polytope.simplex{_scalar_type_to_polymake(scalar_type)}(d, n), parent_field ) end simplex(d::Int, n) = simplex(QQFieldElem, d, n) -function simplex(f::Union{Type{T},Field}, d::Int) where {T<:scalar_types} +function simplex(f::scalar_type_or_field, d::Int) parent_field, scalar_type = _determine_parent_and_scalar(f) return Polyhedron{scalar_type}( Polymake.polytope.simplex{_scalar_type_to_polymake(scalar_type)}(d), parent_field @@ -678,14 +678,14 @@ x₁ - x₂ - x₃ ≦ 2 -x₁ - x₂ - x₃ ≦ 2 ``` """ -function cross_polytope(f::Union{Type{T},Field}, d::Int64, n) where {T<:scalar_types} +function cross_polytope(f::scalar_type_or_field, d::Int64, n) parent_field, scalar_type = _determine_parent_and_scalar(f, n) return Polyhedron{scalar_type}( Polymake.polytope.cross{_scalar_type_to_polymake(scalar_type)}(d, n), parent_field ) end cross_polytope(d::Int64, n) = cross_polytope(QQFieldElem, d, n) -function cross_polytope(f::Union{Type{T},Field}, d::Int64) where {T<:scalar_types} +function cross_polytope(f::scalar_type_or_field, d::Int64) parent_field, scalar_type = _determine_parent_and_scalar(f) return Polyhedron{scalar_type}( Polymake.polytope.cross{_scalar_type_to_polymake(scalar_type)}(d), parent_field diff --git a/src/PolyhedralGeometry/SubdivisionOfPoints/constructors.jl b/src/PolyhedralGeometry/SubdivisionOfPoints/constructors.jl index 9a314dd42c03..ab6e22d7c9db 100644 --- a/src/PolyhedralGeometry/SubdivisionOfPoints/constructors.jl +++ b/src/PolyhedralGeometry/SubdivisionOfPoints/constructors.jl @@ -53,14 +53,14 @@ julia> MOAE = subdivision_of_points(moaepts, moaeimnonreg0) Subdivision of points in ambient dimension 3 ``` """ -function subdivision_of_points(f::Union{Type{T}, Field}, points::AbstractCollection[PointVector], cells::IncidenceMatrix) where T<:scalar_types - @req size(points)[1] == ncols(cells) "Number of points must be the same as columns of IncidenceMatrix" - parent_field, scalar_type = _determine_parent_and_scalar(f, points) - arr = @Polymake.convert_to Array{Set{Int}} Polymake.common.rows(cells) - SubdivisionOfPoints{scalar_type}(Polymake.fan.SubdivisionOfPoints{_scalar_type_to_polymake(scalar_type)}( - POINTS = homogenize(points,1), - MAXIMAL_CELLS = arr, - ), parent_field) +function subdivision_of_points(f::scalar_type_or_field, points::AbstractCollection[PointVector], cells::IncidenceMatrix) + @req size(points)[1] == ncols(cells) "Number of points must be the same as columns of IncidenceMatrix" + parent_field, scalar_type = _determine_parent_and_scalar(f, points) + arr = @Polymake.convert_to Array{Set{Int}} Polymake.common.rows(cells) + SubdivisionOfPoints{scalar_type}(Polymake.fan.SubdivisionOfPoints{_scalar_type_to_polymake(scalar_type)}( + POINTS = homogenize(points,1), + MAXIMAL_CELLS = arr, + ), parent_field) end @@ -91,13 +91,13 @@ julia> n_maximal_cells(SOP) 1 ``` """ -function subdivision_of_points(f::Union{Type{T}, Field}, points::AbstractCollection[PointVector], weights::AbstractVector) where T<:scalar_types - @req size(points)[1] == length(weights) "Number of points must equal number of weights" - parent_field, scalar_type = _determine_parent_and_scalar(f, points, weights) - SubdivisionOfPoints{scalar_type}(Polymake.fan.SubdivisionOfPoints{_scalar_type_to_polymake(scalar_type)}( - POINTS = homogenize(points,1), - WEIGHTS = weights, - ), parent_field) +function subdivision_of_points(f::scalar_type_or_field, points::AbstractCollection[PointVector], weights::AbstractVector) + @req size(points)[1] == length(weights) "Number of points must equal number of weights" + parent_field, scalar_type = _determine_parent_and_scalar(f, points, weights) + SubdivisionOfPoints{scalar_type}(Polymake.fan.SubdivisionOfPoints{_scalar_type_to_polymake(scalar_type)}( + POINTS = homogenize(points,1), + WEIGHTS = weights, + ), parent_field) end """ @@ -109,13 +109,13 @@ pm_object(SOP::SubdivisionOfPoints) = SOP.pm_subdivision #Same construction for when the user provides maximal cells -function subdivision_of_points(::Union{Type{T}, Field}, points::AbstractCollection[PointVector], cells::Vector{Vector{Int64}}) where T<:scalar_types - subdivision_of_points(T, points, IncidenceMatrix(cells)) +function subdivision_of_points(f::scalar_type_or_field, points::AbstractCollection[PointVector], cells::Vector{Vector{Int64}}) + subdivision_of_points(f, points, IncidenceMatrix(cells)) end #Same construction for when the user gives Matrix{Bool} as incidence matrix -function subdivision_of_points(::Union{Type{T}, Field}, Points::Union{Oscar.MatElem,AbstractMatrix}, cells::Matrix{Bool}) where T<:scalar_types - subdivision_of_points(T, points, IncidenceMatrix(Polymake.IncidenceMatrix(cells))) +function subdivision_of_points(f::scalar_type_or_field, Points::Union{Oscar.MatElem,AbstractMatrix}, cells::Matrix{Bool}) + subdivision_of_points(f, points, IncidenceMatrix(Polymake.IncidenceMatrix(cells))) end ############################################################################### diff --git a/src/PolyhedralGeometry/SubdivisionOfPoints/properties.jl b/src/PolyhedralGeometry/SubdivisionOfPoints/properties.jl index be7bc2d2e111..21205271f14f 100644 --- a/src/PolyhedralGeometry/SubdivisionOfPoints/properties.jl +++ b/src/PolyhedralGeometry/SubdivisionOfPoints/properties.jl @@ -29,11 +29,11 @@ julia> points(MOAE) [1, 1, 2] ``` """ -function points(SOP::SubdivisionOfPoints) - return SubObjectIterator{PointVector{QQFieldElem}}(SOP, _point, size(pm_object(SOP).POINTS, 1)) +function points(SOP::SubdivisionOfPoints{T}) where T<:scalar_types + return SubObjectIterator{PointVector{T}}(SOP, _point, size(pm_object(SOP).POINTS, 1)) end -_point(::Type{PointVector{QQFieldElem}}, SOP::SubdivisionOfPoints, i::Base.Integer) = PointVector{QQFieldElem}(pm_object(SOP).POINTS[i, 2:end]) +_point(U::Type{PointVector{T}}, SOP::SubdivisionOfPoints{T}, i::Base.Integer) where T<:scalar_types = point_vector(coefficient_field(SOP), pm_object(SOP).POINTS[i, 2:end])::U _point_matrix(::Val{_point}, SOP::SubdivisionOfPoints; homogenized=false) = pm_object(SOP).POINTS[:, (homogenized ? 1 : 2):end] @@ -185,7 +185,7 @@ julia> min_weights(SOP) 0 ``` """ -min_weights(SOP::SubdivisionOfPoints{T}) where T<:scalar_types = Vector{Int}(pm_object(SOP).MIN_WEIGHTS) +min_weights(SOP::SubdivisionOfPoints) = Vector{Int}(pm_object(SOP).MIN_WEIGHTS) @doc raw""" diff --git a/src/PolyhedralGeometry/helpers.jl b/src/PolyhedralGeometry/helpers.jl index 2ffedf4c83ac..9b6ec0955080 100644 --- a/src/PolyhedralGeometry/helpers.jl +++ b/src/PolyhedralGeometry/helpers.jl @@ -381,11 +381,14 @@ _get_scalar_type(::NormalToricVarietyType) = QQFieldElem const scalar_types = Union{FieldElem, Float64} const scalar_type_to_oscar = Dict{String, Type}([("Rational", QQFieldElem), - ("QuadraticExtension", Hecke.EmbeddedNumFieldElem{nf_elem}), - ("Float", Float64)]) + ("QuadraticExtension", Hecke.EmbeddedNumFieldElem{nf_elem}), + ("QuadraticExtension", Hecke.EmbeddedNumFieldElem{nf_elem}), + ("Float", Float64)]) const scalar_types_extended = Union{scalar_types, ZZRingElem} +const scalar_type_or_field = Union{Type{<:scalar_types}, Field} + _scalar_type_to_polymake(::Type{QQFieldElem}) = Polymake.Rational _scalar_type_to_polymake(::Type{<:FieldElem}) = Polymake.OscarNumber _scalar_type_to_polymake(::Type{Float64}) = Float64 @@ -433,7 +436,7 @@ function _find_parent_field(::Type{T}, x::AbstractArray) where T <: scalar_types return QQ end -_determine_parent_and_scalar(f::Field, x...) = (f, elem_type(f)) +_determine_parent_and_scalar(f::Union{Field, ZZRing}, x...) = (f, elem_type(f)) # isempty(x) => standard/trivial field? function _determine_parent_and_scalar(::Type{T}, x...) where T <: scalar_types if T == QQFieldElem @@ -478,37 +481,56 @@ _detect_default_field(::Type{QQFieldElem}, p::Polymake.BigObject) = QQ _detect_default_field(::Type{Float64}, p::Polymake.BigObject) = AbstractAlgebra.Floats{Float64}() function _detect_default_field(::Type{T}, p::Polymake.BigObject) where T<:FieldElem - # we only want to check existing properties - f = x -> Polymake.exists(p, string(x)) - propnames = intersect(propertynames(p), [:INPUT_RAYS, :POINTS, :RAYS, :VERTICES, :VECTORS, :INPUT_LINEALITY, :LINEALITY_SPACE, :FACETS, :INEQUALITIES, :EQUATIONS, :LINEAR_SPAN, :AFFINE_HULL]) - i = findfirst(f, propnames) - # find first OscarNumber wrapping a FieldElem - while !isnothing(i) - prop = getproperty(p, propnames[i]) - for el in prop - on = Polymake.unwrap(el) - if on isa T - return parent(on) - end - end - i = findnext(f, propnames, i + 1) + # we only want to check existing properties + propnames = intersect(Polymake.list_properties(p), ["INPUT_RAYS", "POINTS", "RAYS", "VERTICES", "VECTORS", "INPUT_LINEALITY", "LINEALITY_SPACE", "FACETS", "INEQUALITIES", "EQUATIONS", "LINEAR_SPAN", "AFFINE_HULL"]) + # find first OscarNumber wrapping a FieldElem + for pn in propnames + prop = getproperty(p, convert(String, pn)) + for el in prop + on = Polymake.unwrap(el) + if on isa T + return parent(on) + end end - throw(ArgumentError("BigObject does not contain information about a parent Field")) + end + throw(ArgumentError("BigObject does not contain information about a parent Field")) +end + +function _detect_wrapped_type_and_field(p::Polymake.BigObject) + # we only want to check existing properties + propnames = intersect(Polymake.list_properties(p), ["INPUT_RAYS", "POINTS", "RAYS", "VERTICES", "VECTORS", "INPUT_LINEALITY", "LINEALITY_SPACE", "FACETS", "INEQUALITIES", "EQUATIONS", "LINEAR_SPAN", "AFFINE_HULL"]) + # find first OscarNumber wrapping a FieldElem + for pn in propnames + prop = getproperty(p, convert(String, pn)) + for el in prop + on = Polymake.unwrap(el) + if on isa FieldElem + f = parent(on) + T = elem_type(f) + return (T, f) + end + end + end + throw(ArgumentError("BigObject does not contain information about a parent Field")) end function _detect_scalar_and_field(::Type{U}, p::Polymake.BigObject) where U<:PolyhedralObject - T = detect_scalar_type(U, p) + T = detect_scalar_type(U, p) + if isnothing(T) + return _detect_wrapped_type_and_field(p) + else return (T, _detect_default_field(T, p)) + end end # promotion helpers -function _promote_scalar_field(f::Field...) - try - x = sum([g(0) for g in f]) - p = parent(x) - return (elem_type(p), p) - catch e - throw(ArgumentError("Can not find a mutual parent field for $f.")) +function _promote_scalar_field(f::Union{Field, ZZRing}...) + try + x = sum([g(0) for g in f]) + p = parent(x) + return (elem_type(p), p) + catch e + throw(ArgumentError("Can not find a mutual parent field for $f.")) end end @@ -517,6 +539,10 @@ function _promote_scalar_field(a::AbstractArray{<:FieldElem}) return _promote_scalar_field(parent.(a)...) end +_parent_or_coefficient_field(r::Base.RefValue{<:Union{FieldElem, ZZRingElem}}) = parent(r.x) + +_parent_or_coefficient_field(v::AbstractVector{T}) where T<:Union{FieldElem, ZZRingElem} = _determine_parent_and_scalar(T, v)[1] + function _promoted_bigobject(::Type{T}, obj::PolyhedralObject{U}) where {T <: scalar_types, U <: scalar_types} T == U ? pm_object(obj) : Polymake.common.convert_to{_scalar_type_to_polymake(T)}(pm_object(obj)) end @@ -530,3 +556,28 @@ end function Polymake._fieldelem_is_rational(e::Hecke.EmbeddedNumFieldElem) return is_rational(e) end + +# convert a Polymake.BigObject's scalar from QuadraticExtension to OscarNumber (Polytope only) + +function _polyhedron_qe_to_on(x::Polymake.BigObject, f::Field) + res = Polymake.polytope.Polytope{Polymake.OscarNumber}() + for pn in Polymake.list_properties(x) + prop = Polymake.give(x, pn) + Polymake.take(res, string(pn), _property_qe_to_on(prop, f)) + end + return res +end + +_property_qe_to_on(x::Polymake.BigObject, f::Field) = Polymake.BigObject(Polymake.bigobject_type(x), x) + +_property_qe_to_on(x::Polymake.PropertyValue, f::Field) = x + +_property_qe_to_on(x::Polymake.QuadraticExtension{Polymake.Rational}, f::Field) = f(x) + +function _property_qe_to_on(x, f::Field) + if hasmethod(length, (typeof(x),)) && eltype(x) <: Polymake.QuadraticExtension{Polymake.Rational} + return f.(x) + else + return x + end +end diff --git a/src/PolyhedralGeometry/iterators.jl b/src/PolyhedralGeometry/iterators.jl index 5779e161ee60..fe0f1154b60c 100644 --- a/src/PolyhedralGeometry/iterators.jl +++ b/src/PolyhedralGeometry/iterators.jl @@ -2,221 +2,190 @@ ######## Vector types ################################################################################ -struct PointVector{U} <: AbstractVector{U} - p::Vector{U} - - PointVector{U}(p::AbstractVector) where U<:scalar_types_extended = new{U}(p) - PointVector(p::AbstractVector) = new{QQFieldElem}(p) -end - -Base.IndexStyle(::Type{<:PointVector}) = IndexLinear() - -Base.getindex(po::PointVector{T}, i::Base.Integer) where T<:scalar_types_extended = po.p[i] - -function Base.setindex!(po::PointVector, val, i::Base.Integer) - @boundscheck checkbounds(po.p, i) - po.p[i] = val - return val -end - -Base.firstindex(::PointVector) = 1 -Base.lastindex(iter::PointVector) = length(iter) -Base.size(po::PointVector) = size(po.p) - -# Forward multiplication with oscar matrices. -Base.:*(A::MatElem, v::PointVector) = A*v.p - -PointVector{U}(p::Union{Field, ZZRing}, v::AbstractVector) where U<:scalar_types_extended = PointVector{U}(p.(v)) +for (T, _t) in ((:PointVector, :point_vector), (:RayVector, :ray_vector)) + + @eval begin -PointVector{U}(p::Union{Field, ZZRing}, n::Base.Integer) where U<:scalar_types_extended = PointVector{U}(p.(zeros(Int, n))) + struct $T{U} <: AbstractVector{U} + p::MatElem{U} -PointVector{U}(::UndefInitializer, n::Base.Integer) where U<:scalar_types_extended = PointVector{U}(Vector{U}(undef, n)) - -PointVector{U}(n::UR) where {U<:scalar_types_extended, UR<:Base.AbstractUnitRange} = PointVector{U}(undef, length(n)) - -function Base.similar(X::PointVector, ::Type{S}, dims::Dims{1}) where S <: scalar_types_extended - return PointVector{S}(undef, dims...) -end + $T{U}(p::MatElem{U}) where U<:scalar_types_extended = new{U}(p) + end -Base.BroadcastStyle(::Type{<:PointVector}) = Broadcast.ArrayStyle{PointVector}() + Base.IndexStyle(::Type{<:$T}) = IndexLinear() + Base.getindex(po::$T, i::Base.Integer) = po.p[1, i] -function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{PointVector}}, ::Type{ElType}) where ElType - return PointVector{ElType}(axes(bc)...) -end + function Base.setindex!(po::$T, val, i::Base.Integer) + @boundscheck 1 <= length(po) <= i + po.p[1, i] = val + return val + end -################################################################################ + Base.firstindex(::$T) = 1 + Base.lastindex(iter::$T) = length(iter) + Base.size(po::$T) = (size(po.p, 2),) -struct RayVector{U} <: AbstractVector{U} - p::Vector{U} - - RayVector{U}(p::AbstractVector) where U<:scalar_types_extended = new{U}(p) - RayVector(p::AbstractVector) = new{QQFieldElem}(p) -end + coefficient_field(po::$T) = base_ring(po.p) -Base.IndexStyle(::Type{<:RayVector}) = IndexLinear() - -Base.getindex(po::RayVector{T}, i::Base.Integer) where T<:scalar_types_extended = po.p[i] + function $_t(p::Union{scalar_type_or_field, ZZRing}, v::AbstractVector) + parent_field, scalar_type = _determine_parent_and_scalar(p, v) + n = length(v) + mat = matrix(parent_field, 1, n, collect(v)) # collect: workaround for constructor typing + return $T{scalar_type}(mat) + end -function Base.setindex!(po::RayVector, val, i::Base.Integer) - @boundscheck checkbounds(po.p, i) - po.p[i] = val - return val -end + function $_t(p::Union{scalar_type_or_field, ZZRing}, n::Base.Integer) + parent_field, scalar_type = _determine_parent_and_scalar(p) + mat = zero_matrix(parent_field, 1, n) + return $T{scalar_type}(mat) + end -Base.firstindex(::RayVector) = 1 -Base.lastindex(iter::RayVector) = length(iter) -Base.size(po::RayVector) = size(po.p) + $_t(x::Union{AbstractVector, Base.Integer}) = $_t(QQ, x) -# Forward multiplication with oscar matrices. -Base.:*(A::MatElem, v::RayVector) = A*v.p + function Base.similar(X::$T, ::Type{S}, dims::Dims{1}) where S <: scalar_types_extended + return $_t(coefficient_field(X), dims...) + end -RayVector{U}(p::Union{Field, ZZRing}, v::AbstractVector) where U<:scalar_types_extended = RayVector{U}(p.(v)) + Base.BroadcastStyle(::Type{<:$T}) = Broadcast.ArrayStyle{$T}() -RayVector{U}(p::Union{Field, ZZRing}, n::Base.Integer) where U<:scalar_types_extended = RayVector{U}(p.(zeros(Int, n))) + _parent_or_coefficient_field(po::$T) = coefficient_field(po) -RayVector{U}(::UndefInitializer, n::Base.Integer) where U<:scalar_types_extended = RayVector{U}(Vector{U}(undef, n)) + function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{$T}}, ::Type{ElType}) where ElType + U, f = _promote_scalar_field(_parent_or_coefficient_field.(bc.args)...) + return $_t(f, axes(bc)...) + end -RayVector{U}(n::UR) where {U<:scalar_types_extended, UR<:Base.AbstractUnitRange} = RayVector{U}(undef, length(n)) + Base.:*(k::scalar_types_extended, po::$T) = k .* po -function Base.similar(X::RayVector, ::Type{S}, dims::Dims{1}) where S <: scalar_types_extended - return RayVector{S}(undef, dims...) + Base.:*(A::MatElem, v::$T) = A * transpose(v.p) + + end end -Base.BroadcastStyle(::Type{<:RayVector}) = Broadcast.ArrayStyle{RayVector}() +@doc """ + point_vector(p = QQ, v::AbstractVector) -function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{RayVector}}, ::Type{ElType}) where ElType - return RayVector{ElType}(axes(bc)...) -end +Return a `PointVector` resembling a point whose coordinates equal the entries of `v`. +`p` specifies the `Field` or `Type` of its coefficients. +""" +point_vector -################################################################################ +@doc """ + ray_vector(p = QQ, v::AbstractVector) -Base.:*(k::scalar_types_extended, po::Union{PointVector, RayVector}) = k .* po +Return a `RayVector` resembling a ray from the origin through the point whose coordinates equal the entries of `v`. +`p` specifies the `Field` or `Type` of its coefficients. +""" +ray_vector ################################################################################ ######## Halfspaces and Hyperplanes ################################################################################ -abstract type Halfspace{T} end - -################################################################################ - -@doc raw""" - Halfspace(a, b) +for (h, comp) in (("halfspace", "≤"), ("hyperplane", "=")) -One halfspace `H(a,b)` is given by a vector `a` and a value `b` such that -$$H(a,b) = \{ x | ax ≤ b \}.$$ -""" -struct AffineHalfspace{T} <: Halfspace{T} - a::Vector{T} - b::T - parent_field::Field + H = uppercasefirst(h) + Habs = Symbol(H) + Haff = Symbol("Affine", H) + Hlin = Symbol("Linear", H) + Fabs = Symbol(h) + Faff = Symbol("affine_", h) + Flin = Symbol("linear_", h) - AffineHalfspace{T}(p::Field, a::Union{MatElem, AbstractMatrix, AbstractVector}, b) where T<:scalar_types = new{T}(p.(vec(a)), p(b), p) -end + @eval begin -halfspace(a::Union{MatElem, AbstractMatrix, AbstractVector}, b) = affine_halfspace(a, b) -halfspace(f::Union{Type{T}, Field}, a::Union{MatElem, AbstractMatrix, AbstractVector}, b) where T<:scalar_types = affine_halfspace(f, a, b) + abstract type $Habs{T} end -invert(H::AffineHalfspace{T}) where T<:scalar_types = AffineHalfspace{T}(coefficient_field(H), -normal_vector(H), -negbias(H)) + # Affine types + struct $Haff{T} <: $Habs{T} + a::MatElem{T} + b::T -function affine_halfspace(f::Union{Type{T}, Field}, a::Union{MatElem, AbstractMatrix, AbstractVector}, b = 0) where T<:scalar_types - parent_field, scalar_type = _determine_parent_and_scalar(f, a, b) - return AffineHalfspace{scalar_type}(parent_field, a, b) -end - -affine_halfspace(a::Union{MatElem, AbstractMatrix, AbstractVector}, b = 0) = affine_halfspace(QQ, a, b) - -################################################################################ - -struct LinearHalfspace{T} <: Halfspace{T} - a::Vector{T} - parent_field::Field - - LinearHalfspace{T}(p::Field, a::Union{MatElem, AbstractMatrix, AbstractVector}) where T<:scalar_types = new{T}(p.(vec(a)), p) -end - -halfspace(a::Union{MatElem, AbstractMatrix, AbstractVector}) = linear_halfspace(a) -halfspace(f::Union{Type{T}, Field}, a::Union{MatElem, AbstractMatrix, AbstractVector}) where T<:scalar_types = linear_halfspace(f, a) - -invert(H::LinearHalfspace{T}) where T<:scalar_types = LinearHalfspace{T}(coefficient_field(H), -normal_vector(H)) + $Haff{T}(a::MatElem{T}, b::T) where T<:scalar_types = new{T}(a, b) + end -function linear_halfspace(f::Union{Type{T}, Field}, a::Union{MatElem, AbstractMatrix, AbstractVector}) where T<:scalar_types - parent_field, scalar_type = _determine_parent_and_scalar(f, a) - return LinearHalfspace{scalar_type}(parent_field, a) -end + $Fabs(a::Union{MatElem, AbstractMatrix, AbstractVector}, b) = $Faff(a, b) + $Fabs(f::scalar_type_or_field, a::Union{MatElem, AbstractMatrix, AbstractVector}, b) = $Faff(f, a, b) -linear_halfspace(a::Union{MatElem, AbstractMatrix, AbstractVector}) = linear_halfspace(QQ, a) + invert(H::$Haff{T}) where T<:scalar_types = $Haff{T}(-H.a, -negbias(H)) -coefficient_field(h::Halfspace) = h.parent_field + @doc """ + $($Faff)(p = QQ, a, b) -################################################################################ +Return the `$($Haff)` `H(a,b)`, which is given by a vector `a` and a value `b` such that +\$\$H(a,b) = \\{ x | ax $($comp) b \\}.\$\$ +`p` specifies the `Field` or `Type` of its coefficients. + """ + function $Faff(f::scalar_type_or_field, a::Union{MatElem, AbstractMatrix, AbstractVector}, b = 0) + parent_field, scalar_type = _determine_parent_and_scalar(f, a, b) + mat = matrix(parent_field, 1, length(a), collect(a)) + return $Haff{scalar_type}(mat, parent_field(b)) + end -abstract type Hyperplane{T} end + $Faff(a::Union{MatElem, AbstractMatrix, AbstractVector}, b = 0) = $Faff(QQ, a, b) -################################################################################ -@doc raw""" - AffineHyperplane(a, b) + # Linear types -One hyperplane `H(a,b)` is given by a vector `a` and a value `b` such that -$$H(a,b) = \{ x | ax = b \}.$$ -""" -struct AffineHyperplane{T} <: Hyperplane{T} - a::Vector{T} - b::T - parent_field::Field + struct $Hlin{T} <: $Habs{T} + a::MatElem{T} - AffineHyperplane{T}(p::Field, a::Union{MatElem, AbstractMatrix, AbstractVector}, b=0) where T<:scalar_types = new{T}(p.(vec(a)), p(b), p) -end + $Hlin{T}(a::MatElem{T}) where T<:scalar_types = new{T}(a) + end -hyperplane(a::Union{MatElem, AbstractMatrix, AbstractVector}, b) = affine_hyperplane(a, b) -hyperplane(f::Union{Type{T}, Field}, a::Union{MatElem, AbstractMatrix, AbstractVector}, b) where T<:scalar_types = affine_hyperplane(f, a, b) + $Fabs(a::Union{MatElem, AbstractMatrix, AbstractVector}) = $Flin(a) + $Fabs(f::scalar_type_or_field, a::Union{MatElem, AbstractMatrix, AbstractVector}) = $Flin(f, a) -function affine_hyperplane(f::Union{Type{T}, Field}, a::Union{MatElem, AbstractMatrix, AbstractVector}, b = 0) where T<:scalar_types - parent_field, scalar_type = _determine_parent_and_scalar(f, a, b) - return AffineHyperplane{scalar_type}(parent_field, a, b) -end + invert(H::$Hlin{T}) where T<:scalar_types = $Hlin{T}(-H.a) -affine_hyperplane(a::Union{MatElem, AbstractMatrix, AbstractVector}, b = 0) = affine_hyperplane(QQ, a, b) + @doc """ + $($Flin)(p = QQ, a, b) -################################################################################ - -struct LinearHyperplane{T} <: Hyperplane{T} - a::Vector{T} - parent_field::Field - - LinearHyperplane{T}(p::Field, a::Union{MatElem, AbstractMatrix, AbstractVector}) where T<:scalar_types = new{T}(p.(vec(a)), p) -end - -hyperplane(a::Union{MatElem, AbstractMatrix, AbstractVector}) = linear_hyperplane(a) -hyperplane(f::Union{Type{T}, Field}, a::Union{MatElem, AbstractMatrix, AbstractVector}) where T<:scalar_types = linear_hyperplane(f, a) +Return the `$($Hlin)` `H(a)`, which is given by a vector `a` such that +\$\$H(a,b) = \\{ x | ax $($comp) 0 \\}.\$\$ +`p` specifies the `Field` or `Type` of its coefficients. + """ + function $Flin(f::scalar_type_or_field, a::Union{MatElem, AbstractMatrix, AbstractVector}) + parent_field, scalar_type = _determine_parent_and_scalar(f, a) + mat = matrix(parent_field, 1, length(a), collect(a)) + return $Hlin{scalar_type}(mat) + end -function linear_hyperplane(f::Union{Type{T}, Field}, a::Union{MatElem, AbstractMatrix, AbstractVector}) where T<:scalar_types - parent_field, scalar_type = _determine_parent_and_scalar(f, a) - return LinearHyperplane{scalar_type}(parent_field, a) -end + $Flin(a::Union{MatElem, AbstractMatrix, AbstractVector}) = $Flin(QQ, a) -linear_hyperplane(a::Union{MatElem, AbstractMatrix, AbstractVector}) = linear_hyperplane(QQ, a) + coefficient_field(h::$Habs) = base_ring(h.a) -coefficient_field(h::Hyperplane) = h.parent_field + end -################################################################################ +end # Field access -negbias(H::Union{AffineHalfspace{T}, AffineHyperplane{T}}) where T<:scalar_types = H.b -negbias(H::Union{LinearHalfspace{T}, LinearHyperplane{T}}) where T<:scalar_types = coefficient_field(H)(0) -normal_vector(H::Union{Halfspace{T}, Hyperplane{T}}) where T <: scalar_types = Vector{T}(H.a) +negbias(H::Union{AffineHalfspace, AffineHyperplane}) = H.b +negbias(H::Union{LinearHalfspace, LinearHyperplane}) = coefficient_field(H)(0) +normal_vector(H::Union{Halfspace, Hyperplane}) = [H.a[1, i] for i in 1:length(H.a)] _ambient_dim(x::Union{Halfspace, Hyperplane}) = length(x.a) -# TODO: abstract notion of equality -Base.:(==)(x::AffineHalfspace, y::AffineHalfspace) = x.a == y.a && x.b == y.b - -Base.:(==)(x::LinearHalfspace, y::LinearHalfspace) = x.a == y.a - -Base.:(==)(x::AffineHyperplane, y::AffineHyperplane) = x.a == y.a && x.b == y.b +function Base.:(==)(x::Halfspace, y::Halfspace) + ax = normal_vector(x) + ay = normal_vector(y) + ix = findfirst(a -> !iszero(a), ax) + iy = findfirst(a -> !iszero(a), ay) + ix == iy || return false + r = y.a[iy]//x.a[ix] + r > 0 || return false + return (r .* ax == ay) && (r * negbias(x) == negbias(y)) +end -Base.:(==)(x::LinearHyperplane, y::LinearHyperplane) = x.a == y.a +function Base.:(==)(x::Hyperplane, y::Hyperplane) + ax = normal_vector(x) + ay = normal_vector(y) + ix = findfirst(a -> !iszero(a), ax) + iy = findfirst(a -> !iszero(a), ay) + ix == iy || return false + r = y.a[iy]//x.a[ix] + return (r .* ax == ay) && (r * negbias(x) == negbias(y)) +end Base.hash(x::T, h::UInt) where {T<:Union{AffineHalfspace,AffineHyperplane}} = hash((x.a, x.b), hash(T, h)) @@ -383,9 +352,8 @@ for f in ("_point_matrix", "_vector_matrix", "_generator_matrix") M = Symbol(f) @eval begin function $M(::Val{_empty_access}, P::PolyhedralObjectUnion; homogenized=false) - scalar_regexp = match(r"[^<]*<(.*)>[^>]*", String(Polymake.type_name(pm_object(P)))) - typename = scalar_regexp[1] - T = _scalar_type_to_polymake(scalar_type_to_oscar[typename]) + typename = Polymake.bigobject_eltype(pm_object(P)) + T = typename == "OscarNumber" ? Polymake.OscarNumber : _scalar_type_to_polymake(scalar_type_to_oscar[typename]) return Polymake.Matrix{T}(undef, 0, Polymake.polytope.ambient_dim(pm_object(P)) + homogenized) end end diff --git a/src/PolyhedralGeometry/linear_program.jl b/src/PolyhedralGeometry/linear_program.jl index dc0ee8882d4d..8f825528e149 100644 --- a/src/PolyhedralGeometry/linear_program.jl +++ b/src/PolyhedralGeometry/linear_program.jl @@ -44,13 +44,13 @@ function linear_program( end linear_program( - f::Union{Type{T},Field}, + f::scalar_type_or_field, A::AbstractCollection[AffineHalfspace], b, c::AbstractVector; k=0, convention=:max, -) where {T<:scalar_types} = +) = linear_program(polyhedron(f, A, b), c; k=k, convention=convention) pm_object(lp::LinearProgram) = lp.polymake_lp @@ -153,7 +153,7 @@ function optimal_vertex(lp::LinearProgram{T}) where {T<:scalar_types} opt_vert = lp.polymake_lp.MINIMAL_VERTEX end if !isnothing(opt_vert) - return PointVector{T}(coefficient_field(lp), view(dehomogenize(opt_vert), :)) + return point_vector(coefficient_field(lp), view(dehomogenize(opt_vert), :))::PointVector{T} else return nothing end diff --git a/src/PolyhedralGeometry/mixed_integer_linear_program.jl b/src/PolyhedralGeometry/mixed_integer_linear_program.jl index 614dfd14c34e..08809e82c9e6 100644 --- a/src/PolyhedralGeometry/mixed_integer_linear_program.jl +++ b/src/PolyhedralGeometry/mixed_integer_linear_program.jl @@ -176,7 +176,7 @@ function optimal_solution(milp::MixedIntegerLinearProgram{T}) where {T<:scalar_t opt_vert = milp.polymake_milp.MINIMAL_SOLUTION end if !isnothing(opt_vert) - return PointVector{T}(dehomogenize(opt_vert)) + return point_vector(coefficient_field(milp), dehomogenize(opt_vert))::PointVector{T} else return nothing end diff --git a/src/PolyhedralGeometry/solving_integrally.jl b/src/PolyhedralGeometry/solving_integrally.jl index 7a497471b753..85c0ec4d2f56 100644 --- a/src/PolyhedralGeometry/solving_integrally.jl +++ b/src/PolyhedralGeometry/solving_integrally.jl @@ -57,7 +57,7 @@ SubObjectIterator{PointVector{ZZRingElem}} julia> for x in it print(A*x," ") end -ZZRingElem[7] ZZRingElem[7] ZZRingElem[7] ZZRingElem[7] ZZRingElem[7] ZZRingElem[7] ZZRingElem[7] ZZRingElem[7] +[7] [7] [7] [7] [7] [7] [7] [7] ``` """ solve_mixed(as::Type{T}, A::ZZMatrix, b::ZZMatrix, C::ZZMatrix, d::ZZMatrix) where {T} = solve_mixed(T, A, b, C, d) @@ -106,7 +106,7 @@ SubObjectIterator{PointVector{ZZRingElem}} julia> for x in it print(A*x," ") end -ZZRingElem[3] ZZRingElem[3] ZZRingElem[3] ZZRingElem[3] +[3] [3] [3] [3] ``` """ solve_mixed(as::Type{T}, A::ZZMatrix, b::ZZMatrix, C::ZZMatrix) where {T} = solve_mixed(T, A, b, C, zero_matrix(FlintZZ, nrows(C), 1)) diff --git a/src/PolyhedralGeometry/type_functions.jl b/src/PolyhedralGeometry/type_functions.jl index ac6af70df478..73bfb92b6375 100644 --- a/src/PolyhedralGeometry/type_functions.jl +++ b/src/PolyhedralGeometry/type_functions.jl @@ -4,9 +4,8 @@ function detect_scalar_type(n::Type{T}, p::Polymake.BigObject) where T<:Union{Polyhedron, Cone, PolyhedralFan, SubdivisionOfPoints, PolyhedralComplex} - scalar_regexp = match(r"[^<]*<(.*)>[^>]*", String(Polymake.type_name(p))) - typename = scalar_regexp[1] - return scalar_type_to_oscar[typename] + typename = Polymake.bigobject_eltype(p) + return typename == "OscarNumber" ? nothing : scalar_type_to_oscar[typename] end scalar_type(::Union{Polyhedron{T}, Cone{T}, Hyperplane{T}, Halfspace{T}}) where T<:scalar_types = T @@ -28,13 +27,13 @@ end function halfspace_matrix_pair(iter::SubObjectIterator{<:Union{Halfspace{T}, Hyperplane{T}, Polyhedron{T}, Cone{T}, Pair{Matrix{T}, T}}}) where T<:scalar_types - try - f = coefficient_field(iter.Obj) - h = affine_matrix_for_polymake(iter) - return (A = matrix(f, h[:, 2:end]), b = Vector{T}(f.(-h[:, 1]))) - catch e - throw(ArgumentError("Halfspace-Matrix-Pair not defined in this context.")) - end + try + f = coefficient_field(iter.Obj) + h = affine_matrix_for_polymake(iter) + return (A = matrix(f, h[:, 2:end]), b = [f(x) for x in -h[:, 1]]) + catch e + throw(ArgumentError("Halfspace-Matrix-Pair not defined in this context.")) + end end for fun in ( diff --git a/src/exports.jl b/src/exports.jl index dd3e6a52d937..cb7bd8045768 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -1120,6 +1120,7 @@ export pitman_stanley_polytope export platonic_solid export point_coordinates export point_matrix +export point_vector export points export pol_elementary_divisors export polarize @@ -1195,6 +1196,7 @@ export rational_equivalence_class export rational_solutions export rational_to_continued_fraction_hirzebruch_jung export ray_indices +export ray_vector export rays export rays_modulo_lineality export real_projective_plane diff --git a/test/PolyhedralGeometry/cone.jl b/test/PolyhedralGeometry/cone.jl index 24cc9af84c45..2098289e1a6c 100644 --- a/test/PolyhedralGeometry/cone.jl +++ b/test/PolyhedralGeometry/cone.jl @@ -4,7 +4,7 @@ NF, sr2 = quadratic_field(2) Qx, x = QQ["x"] K, (a1, a2) = embedded_number_field([x^2 - 2, x^3 - 5], [(0, 2), (0, 2)]) -for f in (QQ, NF, K) +for f in (QQ, K) T = elem_type(f) diff --git a/test/PolyhedralGeometry/polyhedral_complex.jl b/test/PolyhedralGeometry/polyhedral_complex.jl index e9544c058f92..734852859a14 100644 --- a/test/PolyhedralGeometry/polyhedral_complex.jl +++ b/test/PolyhedralGeometry/polyhedral_complex.jl @@ -1,8 +1,8 @@ NF, sr2 = quadratic_field(2) Qx, x = QQ["x"] -K, (a1, a2) = embedded_number_field([x^2 - 2, x^3 - 5], [(0, 2), (0, 2)]) +K, a = Hecke.embedded_field(NF, real_embeddings(NF)[2]) -for f in (QQ, NF, K) +for f in (QQ, K) T = elem_type(f) @testset "PolyhedralComplex{$T}" begin diff --git a/test/PolyhedralGeometry/polyhedral_fan.jl b/test/PolyhedralGeometry/polyhedral_fan.jl index 4b5e9065c58c..f08d3b5e0b17 100644 --- a/test/PolyhedralGeometry/polyhedral_fan.jl +++ b/test/PolyhedralGeometry/polyhedral_fan.jl @@ -1,9 +1,7 @@ -NF, sr2 = quadratic_field(2) -ENF, sre2 = Hecke.embedded_field(NF, real_embeddings(NF)[2]) Qx, x = QQ["x"] K, (a1, a2) = embedded_number_field([x^2 - 2, x^3 - 5], [(0, 2), (0, 2)]) -for f in (QQ, ENF, K) +for f in (QQ, K) T = elem_type(f) @testset "PolyhedralFan{$T}" begin diff --git a/test/PolyhedralGeometry/polyhedron.jl b/test/PolyhedralGeometry/polyhedron.jl index 594edf872176..05d189d3e60a 100644 --- a/test/PolyhedralGeometry/polyhedron.jl +++ b/test/PolyhedralGeometry/polyhedron.jl @@ -483,11 +483,18 @@ for f in (QQ, ENF) end - if T == nf_elem + if T == EmbeddedElem{nf_elem} @testset "Dodecahedron" begin - R, a = quadratic_field(5) + D = polyhedron(Polymake.polytope.dodecahedron()) + R = coefficient_field(D) + NF = number_field(R) + let isq = Hecke.isquadratic_type(NF) + @test isq[1] + @test isq[2] == 5 + end + a = R(gens(NF)[]) V = [[1//2, a//4 + 3//4, 0], [-1//2, a//4 + 3//4, 0], @@ -510,9 +517,7 @@ for f in (QQ, ENF) [1//2, -a//4 - 3//4, 0], [-1//2, -a//4 - 3//4, 0]] - D = Polyhedron{T}(Polymake.polytope.dodecahedron(), R) @test D isa Polyhedron{T} - @test polyhedron(Polymake.polytope.dodecahedron()) == D @test nvertices(D) == 20 @test vertices(D) == V @@ -527,13 +532,14 @@ for f in (QQ, ENF) if S == Pair{Matrix{T}, T} @test facets(S, D) == [Pair(A[i], b[i]) for i in 1:12] elseif S == Polyhedron{T} - @test facets(S, D) == [polyhedron(R, A[i], b[i]) for i in 1:12] + @test nvertices.(facets(S, D)) == repeat([5], 12) else @test facets(S, D) == [affine_halfspace(R, A[i], b[i]) for i in 1:12] end @test length(facets(S, D)) == 12 @test affine_inequality_matrix(facets(S, D)) == matrix(R, hcat(-b, vcat(A...))) - @test halfspace_matrix_pair(facets(S, D)).A == vcat(A...) && halfspace_matrix_pair(facets(S, D)).b == b + @test halfspace_matrix_pair(facets(S, D)) == (A = matrix(R, vcat(A...)), b = b) + @test ray_indices(facets(S, D)) == IncidenceMatrix(12, 0) @test vertex_indices(facets(S, D)) == IncidenceMatrix([[1, 3, 5, 9, 10], [1, 2, 3, 4, 6], [1, 2, 5, 7, 8], [11, 12, 16, 18, 20], [5, 8, 10, 15, 17], [2, 4, 7, 11, 12], [3, 6, 9, 13, 14], [4, 6, 11, 13, 16], [13, 14, 16, 19, 20], [7, 8, 12, 15, 18], [9, 10, 14, 17, 19], [15, 17, 18, 19, 20]]) @test vertex_and_ray_indices(facets(S, D)) == IncidenceMatrix([[1, 3, 5, 9, 10], [1, 2, 3, 4, 6], [1, 2, 5, 7, 8], [11, 12, 16, 18, 20], [5, 8, 10, 15, 17], [2, 4, 7, 11, 12], [3, 6, 9, 13, 14], [4, 6, 11, 13, 16], [13, 14, 16, 19, 20], [7, 8, 12, 15, 18], [9, 10, 14, 17, 19], [15, 17, 18, 19, 20]]) @@ -562,9 +568,6 @@ for f in (QQ, ENF) @test faces(D, 0) == convex_hull.(T, V) @test isempty(affine_hull(D)) @test relative_interior_point(D) == [0, 0, 0] - - @test platonic_solid("dodecahedron") == D - end end diff --git a/test/PolyhedralGeometry/scalar_types.jl b/test/PolyhedralGeometry/scalar_types.jl index 2134693e2489..bca6a0799252 100644 --- a/test/PolyhedralGeometry/scalar_types.jl +++ b/test/PolyhedralGeometry/scalar_types.jl @@ -70,7 +70,7 @@ @test isq[2] == 2 end let j = johnson_solid(1) - jj = polyhedron(Polymake.polytope.Polytope{Polymake.QuadraticExtension{Polymake.Rational}}(POINTS=Oscar.pm_object(j).VERTICES)) + jj = polyhedron(Polymake.polytope.Polytope{Polymake.OscarNumber}(POINTS=Oscar.pm_object(j).VERTICES)) @test number_field(coefficient_field(j)) == number_field(coefficient_field(jj)) end end @@ -186,4 +186,36 @@ end end + @testset "QuadraticExtension-templated sub-objects" begin + + j = johnson_solid(2) + + let f = facets(Polyhedron, j) + fj = normal_vector.(facets(j)) + fj = [fj; -fj] + for i in 1:nfacets(j) + @test normal_vector(affine_hull(f[i])[]) in fj + end + @test halfspace_matrix_pair(f) isa NamedTuple{(:A, :b), Tuple{AbstractAlgebra.Generic.MatSpaceElem{EmbeddedElem{nf_elem}}, Vector{EmbeddedElem{nf_elem}}}} + g = halfspace_matrix_pair(f) + @test affine_halfspace(coefficient_field(j), g.A[1, :], g.b[1]) in facets(j) + end + for n in (1, 2) # faces which are facets have a different access function + let f = faces(j, n) + for i in 1:Int(f_vector(j)[n + 1]) + @test issubset(vertices(f[i]), vertices(j)) + end + end + end + + k = face_fan(j) + let f = cones(k, 3) + for i in 1:Int(f_vector(k)[3]) + @test issubset(rays(f[i]), rays(k)) + end + l = f[1] + end + + end + end diff --git a/test/PolyhedralGeometry/types.jl b/test/PolyhedralGeometry/types.jl index 756d43e26bca..3256888f9dfb 100644 --- a/test/PolyhedralGeometry/types.jl +++ b/test/PolyhedralGeometry/types.jl @@ -1,157 +1,155 @@ @testset "types" begin - @testset "IncidenceMatrix" begin + @testset "IncidenceMatrix" begin - im = IncidenceMatrix([[1,2,3],[4,5,6]]) - @test nrows(im) == 2 - @test ncols(im) == 6 - @test row(im, 1) isa Set{Int} - @test row(im, 1) == Set{Int}([1, 2, 3]) - @test column(im, 2) isa Set{Int} - @test column(im, 2) == Set{Int}([1]) + im = IncidenceMatrix([[1,2,3],[4,5,6]]) + @test nrows(im) == 2 + @test ncols(im) == 6 + @test row(im, 1) isa Set{Int} + @test row(im, 1) == Set{Int}([1, 2, 3]) + @test column(im, 2) isa Set{Int} + @test column(im, 2) == Set{Int}([1]) - end - - a = [1, 2, 3] - b = [8, 6, 4] + end - NF, sr2 = quadratic_field(2) + a = [1, 2, 3] + b = [8, 6, 4] - #TODO: RayVector - @testset "$T" for T in (PointVector, RayVector) + NF, sr2 = quadratic_field(2) + ENF, sre2 = Hecke.embedded_field(NF, real_embeddings(NF)[2]) - @test T(a) isa T{QQFieldElem} + @testset "$T" for (T, fun) in ((PointVector, point_vector), (RayVector, ray_vector)) - for f in (ZZ, QQ, NF) + @test fun(a) isa T{QQFieldElem} - U = elem_type(f) - @testset "$T{$U}" begin + for f in (ZZ, QQ, ENF) - @test T{U} <: AbstractVector - @test T{U} <: AbstractVector{U} + U = elem_type(f) + @testset "$T{$U}" begin - @test T{U}(f.(a)) isa T - @test T{U}(f.(a)) isa T{U} + @test T{U} <: AbstractVector + @test T{U} <: AbstractVector{U} - @test T{U}(f, 7) == f.(zeros(Int, 7)) + @test fun(f, a) isa T + @test fun(f, a) isa T{U} - A = T{U}(f.(a)) + @test fun(f, 7) == f.(zeros(Int, 7)) - @test A[1] isa U - @test A[2] == 2 - @test A[begin] == 1 - @test A[end] == 3 - A[3] = f(7) - @test A[end] == 7 - A[3] = f(3) + A = fun(f, a) - @test size(A) == (3,) + @test A[1] isa U + @test A[2] == 2 + @test A[begin] == 1 + @test A[end] == 3 + A[3] = f(7) + @test A[end] == 7 + A[3] = f(3) - @test_throws BoundsError A[0] - @test_throws BoundsError A[4] + @test size(A) == (3,) - for g in (ZZ, QQ, NF) - V = elem_type(g) + @test_throws BoundsError A[0] + @test_throws BoundsError A[4] - B = T{V}(g.(b)) + for g in (ZZ, QQ, ENF) + V = elem_type(g) - for op in [+, -] - @test op(A, B) isa T + B = fun(g, b) - @test op(A, B) == op(a, b) - end + for op in [+, -] + @test op(A, B) isa T - @test *(g(3), A) isa T + @test op(A, B) == op(a, b) + end - @test *(g(3), A) == 3 * a + @test *(g(3), A) isa T - @test [A; B] == [a; b] + @test *(g(3), A) == 3 * a - end - - end end + end end + end + - @testset "$T" for (T, f) in ((AffineHalfspace, affine_halfspace), (AffineHyperplane, affine_hyperplane)) + @testset "$T" for (T, f) in ((AffineHalfspace, affine_halfspace), (AffineHyperplane, affine_hyperplane)) - for p in [QQ, NF] - U = elem_type(p) - @test f(p, a, 0) isa T{U} - @test f(p, permutedims(a), 0) isa T{U} + for p in [QQ, ENF] + U = elem_type(p) + @test f(p, a, 0) isa T{U} + @test f(p, permutedims(a), 0) isa T{U} - @test f(p, a, 0) == f(p, permutedims(a), 0) == f(p, a) == f(p, permutedims(a)) + @test f(p, a, 0) == f(p, permutedims(a), 0) == f(p, a) == f(p, permutedims(a)) - A = f(p, a, 0) - B = f(p, b, 2) + A = f(p, a, 0) + B = f(p, b, 2) - @test A != B + @test A != B - @test normal_vector(A) isa Vector{U} - @test normal_vector(A) == a - @test negbias(A) isa U - @test negbias(A) == 0 + @test normal_vector(A) isa Vector{U} + @test normal_vector(A) == a + @test negbias(A) isa U + @test negbias(A) == 0 - @test normal_vector(B) == b - @test negbias(B) == 2 - end + @test normal_vector(B) == b + @test negbias(B) == 2 + end - @test f(a, 0) isa T{QQFieldElem} + @test f(a, 0) isa T{QQFieldElem} - end + end - @testset "$T" for (T, f) in ((LinearHalfspace, linear_halfspace), (LinearHyperplane, linear_hyperplane)) + @testset "$T" for (T, f) in ((LinearHalfspace, linear_halfspace), (LinearHyperplane, linear_hyperplane)) - for p in [QQ, NF] - U = elem_type(p) - @test f(p, a) isa T{U} - @test f(p, permutedims(a)) isa T{U} + for p in [QQ, ENF] + U = elem_type(p) + @test f(p, a) isa T{U} + @test f(p, permutedims(a)) isa T{U} - @test f(p, a) == f(p, permutedims(a)) + @test f(p, a) == f(p, permutedims(a)) - A = f(p, a) - B = f(p, b) + A = f(p, a) + B = f(p, b) - @test A != B + @test A != B - @test normal_vector(A) isa Vector{U} - @test normal_vector(A) == a - @test negbias(A) isa U - @test negbias(A) == 0 + @test normal_vector(A) isa Vector{U} + @test normal_vector(A) == a + @test negbias(A) isa U + @test negbias(A) == 0 - @test normal_vector(B) == b - @test negbias(B) == 0 - end + @test normal_vector(B) == b + @test negbias(B) == 0 + end - @test f(a) isa T{QQFieldElem} + @test f(a) isa T{QQFieldElem} - end + end - for p in [QQ, NF] - U = elem_type(p) - let A = linear_halfspace(p, a) - @test invert(A) isa LinearHalfspace{U} - Ai = invert(A) - @test normal_vector(Ai) == -a - end - let A = affine_halfspace(p, a, 8) - @test invert(A) isa AffineHalfspace{U} - Ai = invert(A) - @test normal_vector(Ai) == -a - @test negbias(Ai) == -8 - end + for p in [QQ, NF] + U = elem_type(p) + let A = linear_halfspace(p, a) + @test invert(A) isa LinearHalfspace{U} + Ai = invert(A) + @test normal_vector(Ai) == -a + end + let A = affine_halfspace(p, a, 8) + @test invert(A) isa AffineHalfspace{U} + Ai = invert(A) + @test normal_vector(Ai) == -a + @test negbias(Ai) == -8 end + end - @test halfspace(a) isa LinearHalfspace{QQFieldElem} - @test hyperplane(a) isa LinearHyperplane{QQFieldElem} - @test halfspace(a, 0) isa AffineHalfspace{QQFieldElem} - @test hyperplane(a, 0) isa AffineHyperplane{QQFieldElem} + @test halfspace(a) isa LinearHalfspace{QQFieldElem} + @test hyperplane(a) isa LinearHyperplane{QQFieldElem} + @test halfspace(a, 0) isa AffineHalfspace{QQFieldElem} + @test hyperplane(a, 0) isa AffineHyperplane{QQFieldElem} end