From ab444bbc46b5493588b2c3217ed35d86396d1260 Mon Sep 17 00:00:00 2001 From: Lars Kastner Date: Thu, 11 Apr 2024 15:36:13 +0200 Subject: [PATCH 1/3] PolyhedralGeometry: Format most of the files that aren't being worked on right now --- src/PolyhedralGeometry/Cone/constructors.jl | 110 +++-- src/PolyhedralGeometry/Cone/properties.jl | 11 +- .../Cone/standard_constructions.jl | 82 ++-- src/PolyhedralGeometry/LP_file_format.jl | 49 +- .../PolyhedralComplex/constructors.jl | 120 +++-- .../PolyhedralComplex/properties.jl | 104 +++-- .../standard_constructions.jl | 29 +- .../PolyhedralFan/constructors.jl | 117 +++-- .../PolyhedralFan/properties.jl | 117 +++-- .../PolyhedralFan/standard_constructions.jl | 110 +++-- src/PolyhedralGeometry/PolyhedralGeometry.jl | 3 +- .../Polyhedron/constructors.jl | 164 ++++--- .../SubdivisionOfPoints/constructors.jl | 62 ++- .../SubdivisionOfPoints/functions.jl | 6 +- .../SubdivisionOfPoints/properties.jl | 31 +- src/PolyhedralGeometry/groups.jl | 205 ++++----- src/PolyhedralGeometry/helpers.jl | 430 +++++++++++------- src/PolyhedralGeometry/iterators.jl | 363 +++++++++------ src/PolyhedralGeometry/linear_program.jl | 7 +- .../mixed_integer_linear_program.jl | 10 +- src/PolyhedralGeometry/solving_integrally.jl | 88 ++-- src/PolyhedralGeometry/triangulations.jl | 263 ++++++----- src/PolyhedralGeometry/type_functions.jl | 73 +-- src/PolyhedralGeometry/visualization.jl | 12 +- 24 files changed, 1512 insertions(+), 1054 deletions(-) diff --git a/src/PolyhedralGeometry/Cone/constructors.jl b/src/PolyhedralGeometry/Cone/constructors.jl index bcaf0bd014c7..5f8e0a6704a8 100644 --- a/src/PolyhedralGeometry/Cone/constructors.jl +++ b/src/PolyhedralGeometry/Cone/constructors.jl @@ -8,13 +8,13 @@ #interior description? struct Cone{T} <: PolyhedralObject{T} #a real polymake polyhedron - pm_cone::Polymake.BigObject - parent_field::Field - - # only allowing scalar_types; - # can be improved by testing if the template type of the `BigObject` corresponds to `T` - Cone{T}(c::Polymake.BigObject, f::Field) where T<:scalar_types = new{T}(c, f) - Cone{QQFieldElem}(c::Polymake.BigObject) = new{QQFieldElem}(c, QQ) + pm_cone::Polymake.BigObject + parent_field::Field + + # only allowing scalar_types; + # can be improved by testing if the template type of the `BigObject` corresponds to `T` + Cone{T}(c::Polymake.BigObject, f::Field) where {T<:scalar_types} = new{T}(c, f) + Cone{QQFieldElem}(c::Polymake.BigObject) = new{QQFieldElem}(c, QQ) end # default scalar type: `QQFieldElem` @@ -23,8 +23,8 @@ cone(x...; kwargs...) = Cone{QQFieldElem}(x...; kwargs...) # Automatic detection of corresponding OSCAR scalar type; # Avoid, if possible, to increase type stability function cone(p::Polymake.BigObject) - T, f = _detect_scalar_and_field(Cone, p) - return Cone{T}(p, f) + T, f = _detect_scalar_and_field(Cone, p) + return Cone{T}(p, f) end @doc raw""" @@ -62,7 +62,12 @@ julia> HS = positive_hull(R, L) Polyhedral cone in ambient dimension 2 ``` """ -function positive_hull(f::scalar_type_or_field, R::AbstractCollection[RayVector], L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) +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) @@ -70,28 +75,50 @@ function positive_hull(f::scalar_type_or_field, R::AbstractCollection[RayVector] 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) + 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) + 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(_guess_fieldelem_type(R, L), R, L; non_redundant=non_redundant) -cone(R::AbstractCollection[RayVector], L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) = positive_hull(_guess_fieldelem_type(R, L), R, L; non_redundant=non_redundant) -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) +positive_hull( + R::AbstractCollection[RayVector], + L::Union{AbstractCollection[RayVector],Nothing}=nothing; + non_redundant::Bool=false, +) = positive_hull(_guess_fieldelem_type(R, L), R, L; non_redundant=non_redundant) +cone( + R::AbstractCollection[RayVector], + L::Union{AbstractCollection[RayVector],Nothing}=nothing; + non_redundant::Bool=false, +) = positive_hull(_guess_fieldelem_type(R, L), R, L; non_redundant=non_redundant) +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 - return Polymake.polytope.equal_polyhedra(pm_object(C0), pm_object(C1))::Bool +function ==(C0::Cone{T}, C1::Cone{T}) where {T<:scalar_types} + return Polymake.polytope.equal_polyhedra(pm_object(C0), pm_object(C1))::Bool end # For a proper hash function for cones we should use a "normal form", # which would require a potentially expensive convex hull computation # (and even that is not enough). But hash methods should be fast, so we # just consider the ambient dimension and the precise type of the cone. -function Base.hash(x::T, h::UInt) where {T <: Cone} +function Base.hash(x::T, h::UInt) where {T<:Cone} h = hash(ambient_dim(x), h) h = hash(T, h) return h @@ -120,15 +147,34 @@ julia> rays(C) [1, 1] ``` """ -function cone_from_inequalities(f::scalar_type_or_field, I::AbstractCollection[LinearHalfspace], E::Union{Nothing, AbstractCollection[LinearHyperplane]} = nothing; non_redundant::Bool = false) +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) + EM = if isnothing(E) || isempty(E) + Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(IM, 2)) + else + linear_matrix_for_polymake(E) + end if non_redundant - return Cone{scalar_type}(Polymake.polytope.Cone{_scalar_type_to_polymake(scalar_type)}(FACETS = IM, LINEAR_SPAN = EM), parent_field) + 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) + return Cone{scalar_type}( + Polymake.polytope.Cone{_scalar_type_to_polymake(scalar_type)}(; + INEQUALITIES=IM, EQUATIONS=EM + ), + parent_field, + ) end end @@ -157,16 +203,21 @@ julia> dim(C) 1 ``` """ -function cone_from_equations(f::scalar_type_or_field, E::AbstractCollection[LinearHyperplane]; non_redundant::Bool = false) +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) + return cone_from_inequalities(f, IM, EM; non_redundant=non_redundant) end cone_from_inequalities(x...) = cone_from_inequalities(QQFieldElem, x...) -cone_from_equations(E::AbstractCollection[LinearHyperplane]; non_redundant::Bool = false) = cone_from_equations(_guess_fieldelem_type(E), E; non_redundant = non_redundant) +cone_from_equations(E::AbstractCollection[LinearHyperplane]; non_redundant::Bool=false) = + cone_from_equations(_guess_fieldelem_type(E), E; non_redundant=non_redundant) """ pm_object(C::Cone) @@ -175,16 +226,15 @@ Get the underlying polymake `Cone`. """ pm_object(C::Cone) = C.pm_cone - ############################################################################### ############################################################################### ### Display ############################################################################### ############################################################################### -function Base.show(io::IO, C::Cone{T}) where T<:scalar_types - print(io, "Polyhedral cone in ambient dimension $(ambient_dim(C))") - T != QQFieldElem && print(io, " with $T type coefficients") +function Base.show(io::IO, C::Cone{T}) where {T<:scalar_types} + print(io, "Polyhedral cone in ambient dimension $(ambient_dim(C))") + T != QQFieldElem && print(io, " with $T type coefficients") end Polymake.visual(C::Cone; opts...) = Polymake.visual(pm_object(C); opts...) diff --git a/src/PolyhedralGeometry/Cone/properties.jl b/src/PolyhedralGeometry/Cone/properties.jl index 9b769a6fa69a..cca656845867 100644 --- a/src/PolyhedralGeometry/Cone/properties.jl +++ b/src/PolyhedralGeometry/Cone/properties.jl @@ -645,11 +645,12 @@ _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) = if homogenized - homogenize(pm_object(C).HILBERT_BASIS_GENERATORS[1], 0) -else - pm_object(C).HILBERT_BASIS_GENERATORS[1] -end +_generator_matrix(::Val{_hilbert_generator}, C::Cone; homogenized=false) = + if homogenized + homogenize(pm_object(C).HILBERT_BASIS_GENERATORS[1], 0) + else + pm_object(C).HILBERT_BASIS_GENERATORS[1] + end _matrix_for_polymake(::Val{_hilbert_generator}) = _generator_matrix diff --git a/src/PolyhedralGeometry/Cone/standard_constructions.jl b/src/PolyhedralGeometry/Cone/standard_constructions.jl index b24c07c6886c..ee9070baab61 100644 --- a/src/PolyhedralGeometry/Cone/standard_constructions.jl +++ b/src/PolyhedralGeometry/Cone/standard_constructions.jl @@ -27,13 +27,12 @@ julia> dim(C01) ``` """ function intersect(C::Cone...) - T, f = _promote_scalar_field((coefficient_field(c) for c in C)...) - pmo = [pm_object(c) for c in C] - return Cone{T}(Polymake.polytope.intersection(pmo...), f) + T, f = _promote_scalar_field((coefficient_field(c) for c in C)...) + pmo = [pm_object(c) for c in C] + return Cone{T}(Polymake.polytope.intersection(pmo...), f) end intersect(C::AbstractVector{<:Cone}) = intersect(C...) - @doc raw""" polarize(C::Cone) @@ -54,11 +53,10 @@ julia> rays(Cv) [0, 1] ``` """ -function polarize(C::Cone{T}) where T<:scalar_types - return Cone{T}(Polymake.polytope.polarize(pm_object(C)), coefficient_field(C)) +function polarize(C::Cone{T}) where {T<:scalar_types} + return Cone{T}(Polymake.polytope.polarize(pm_object(C)), coefficient_field(C)) end - @doc raw""" transform(C::Cone{T}, A::AbstractMatrix) where T<:scalar_types @@ -94,49 +92,51 @@ julia> c == ctt true ``` """ -function transform(C::Cone{T}, A::Union{AbstractMatrix{<:Union{Number, FieldElem}}, MatElem{<:FieldElem}}) where {T<:scalar_types} +function transform( + C::Cone{T}, A::Union{AbstractMatrix{<:Union{Number,FieldElem}},MatElem{<:FieldElem}} +) where {T<:scalar_types} @assert ambient_dim(C) == nrows(A) "Incompatible dimension of cone and transformation matrix" @assert nrows(A) == ncols(A) "Transformation matrix must be square" @assert Polymake.common.rank(A) == nrows(A) "Transformation matrix must have full rank." return _transform(C, A) end -function _transform(C::Cone{T}, A::AbstractMatrix{<:FieldElem}) where T<:scalar_types - U, f = _promote_scalar_field(A) - V, g = _promote_scalar_field(coefficient_field(C), f) - OT = _scalar_type_to_polymake(V) - raymod = Polymake.Matrix{OT}(permutedims(A)) - facetmod = Polymake.Matrix{OT}(Polymake.common.inv(permutedims(raymod))) - return _transform(C, raymod, facetmod, g) +function _transform(C::Cone{T}, A::AbstractMatrix{<:FieldElem}) where {T<:scalar_types} + U, f = _promote_scalar_field(A) + V, g = _promote_scalar_field(coefficient_field(C), f) + OT = _scalar_type_to_polymake(V) + raymod = Polymake.Matrix{OT}(permutedims(A)) + facetmod = Polymake.Matrix{OT}(Polymake.common.inv(permutedims(raymod))) + return _transform(C, raymod, facetmod, g) end -function _transform(C::Cone{T}, A::AbstractMatrix{<:Number}) where T<:scalar_types - OT = _scalar_type_to_polymake(T) - raymod = Polymake.Matrix{OT}(permutedims(A)) - facetmod = Polymake.Matrix{OT}(Polymake.common.inv(permutedims(raymod))) - return _transform(C, raymod, facetmod, coefficient_field(C)) +function _transform(C::Cone{T}, A::AbstractMatrix{<:Number}) where {T<:scalar_types} + OT = _scalar_type_to_polymake(T) + raymod = Polymake.Matrix{OT}(permutedims(A)) + facetmod = Polymake.Matrix{OT}(Polymake.common.inv(permutedims(raymod))) + return _transform(C, raymod, facetmod, coefficient_field(C)) end -function _transform(C::Cone{T}, A::MatElem{U}) where {T<:scalar_types, U<:FieldElem} - V, f = _promote_scalar_field(coefficient_field(C), base_ring(A)) - OT = _scalar_type_to_polymake(V) - raymod = Polymake.Matrix{OT}(transpose(A)) - facetmod = Polymake.Matrix{OT}(inv(A)) - return _transform(C, raymod, facetmod, f) +function _transform(C::Cone{T}, A::MatElem{U}) where {T<:scalar_types,U<:FieldElem} + V, f = _promote_scalar_field(coefficient_field(C), base_ring(A)) + OT = _scalar_type_to_polymake(V) + raymod = Polymake.Matrix{OT}(transpose(A)) + facetmod = Polymake.Matrix{OT}(inv(A)) + return _transform(C, raymod, facetmod, f) end -function _transform(C::Cone{T}, raymod, facetmod, f::Field) where T<:scalar_types - U = elem_type(f) - OT = _scalar_type_to_polymake(U) - result = Polymake.polytope.Cone{OT}() - for prop in ("RAYS", "INPUT_RAYS", "LINEALITY_SPACE", "INPUT_LINEALITY") - if Polymake.exists(pm_object(C), prop) - resultprop = Polymake.Matrix{OT}(Polymake.give(pm_object(C), prop) * raymod) - Polymake.take(result, prop, resultprop) - end +function _transform(C::Cone{T}, raymod, facetmod, f::Field) where {T<:scalar_types} + U = elem_type(f) + OT = _scalar_type_to_polymake(U) + result = Polymake.polytope.Cone{OT}() + for prop in ("RAYS", "INPUT_RAYS", "LINEALITY_SPACE", "INPUT_LINEALITY") + if Polymake.exists(pm_object(C), prop) + resultprop = Polymake.Matrix{OT}(Polymake.give(pm_object(C), prop) * raymod) + Polymake.take(result, prop, resultprop) end - for prop in ("INEQUALITIES", "EQUATIONS", "LINEAR_SPAN", "FACETS") - if Polymake.exists(pm_object(C), prop) - resultprop = Polymake.Matrix{OT}(Polymake.give(pm_object(C), prop) * facetmod) - Polymake.take(result, prop, resultprop) - end + end + for prop in ("INEQUALITIES", "EQUATIONS", "LINEAR_SPAN", "FACETS") + if Polymake.exists(pm_object(C), prop) + resultprop = Polymake.Matrix{OT}(Polymake.give(pm_object(C), prop) * facetmod) + Polymake.take(result, prop, resultprop) end - return Cone{U}(result, f) + end + return Cone{U}(result, f) end diff --git a/src/PolyhedralGeometry/LP_file_format.jl b/src/PolyhedralGeometry/LP_file_format.jl index e83a34a414ab..6c9c2ac5f94f 100644 --- a/src/PolyhedralGeometry/LP_file_format.jl +++ b/src/PolyhedralGeometry/LP_file_format.jl @@ -3,27 +3,31 @@ function _internal_save_lp(io::IO, p::Polymake.BigObject, lp::Polymake.BigObject # the string to the IO object # to pass the object to the polymake shell we use a randomly generated variable name rstr = randstring(['A':'Z'; 'a':'z'], 12) - Polymake.call_function(:User, :set_shell_scalar, rstr*"_poly", p) - Polymake.call_function(:User, :set_shell_scalar, rstr*"_lp", lp) - out, err = Polymake.shell_execute("""polytope::poly2lp(\$$(rstr)_poly, \$$(rstr)_lp, $(max ? 1 : 0));""") + Polymake.call_function(:User, :set_shell_scalar, rstr * "_poly", p) + Polymake.call_function(:User, :set_shell_scalar, rstr * "_lp", lp) + out, err = Polymake.shell_execute( + """polytope::poly2lp(\$$(rstr)_poly, \$$(rstr)_lp, $(max ? 1 : 0));""" + ) write(io, out) - Polymake.call_function(:User, :set_shell_scalar, rstr*"_poly", 0) - Polymake.call_function(:User, :set_shell_scalar, rstr*"_lp", 0) + Polymake.call_function(:User, :set_shell_scalar, rstr * "_poly", 0) + Polymake.call_function(:User, :set_shell_scalar, rstr * "_lp", 0) nothing end function _internal_save_mps(io::IO, p::Polymake.BigObject, lp::Polymake.BigObject) rstr = randstring(['A':'Z'; 'a':'z'], 12) - Polymake.call_function(:User, :set_shell_scalar, rstr*"_poly", p) - Polymake.call_function(:User, :set_shell_scalar, rstr*"_lp", lp) + Polymake.call_function(:User, :set_shell_scalar, rstr * "_poly", p) + Polymake.call_function(:User, :set_shell_scalar, rstr * "_lp", lp) out, err = Polymake.shell_execute("""polytope::poly2mps(\$$(rstr)_poly, \$$(rstr)_lp);""") write(io, out) - Polymake.call_function(:User, :set_shell_scalar, rstr*"_poly", 0) - Polymake.call_function(:User, :set_shell_scalar, rstr*"_lp", 0) + Polymake.call_function(:User, :set_shell_scalar, rstr * "_poly", 0) + Polymake.call_function(:User, :set_shell_scalar, rstr * "_lp", 0) nothing end -function _internal_save_lp(filename::String, p::Polymake.BigObject, lp::Polymake.BigObject, max::Bool) +function _internal_save_lp( + filename::String, p::Polymake.BigObject, lp::Polymake.BigObject, max::Bool +) Polymake.polytope.poly2lp(p, lp, max, filename) nothing end @@ -33,8 +37,6 @@ function _internal_save_mps(filename::String, p::Polymake.BigObject, lp::Polymak nothing end - - @doc raw""" save_lp(target::Union{String, IO}, lp::Union{MixedIntegerLinearProgram,LinearProgram}) @@ -67,11 +69,13 @@ GENERAL END ``` """ -function save_lp(target::Union{String,IO}, lp::Union{MixedIntegerLinearProgram{QQFieldElem},LinearProgram{QQFieldElem}}) - _internal_save_lp(target, - pm_object(feasible_region(lp)), - pm_object(lp), - lp.convention == :max) +function save_lp( + target::Union{String,IO}, + lp::Union{MixedIntegerLinearProgram{QQFieldElem},LinearProgram{QQFieldElem}}, +) + _internal_save_lp( + target, pm_object(feasible_region(lp)), pm_object(lp), lp.convention == :max + ) end @doc raw""" @@ -119,10 +123,11 @@ BOUNDS ENDATA ``` """ -function save_mps(target::Union{String,IO}, lp::Union{MixedIntegerLinearProgram{QQFieldElem},LinearProgram{QQFieldElem}}) - _internal_save_mps(target, - pm_object(feasible_region(lp)), - pm_object(lp)) +function save_mps( + target::Union{String,IO}, + lp::Union{MixedIntegerLinearProgram{QQFieldElem},LinearProgram{QQFieldElem}}, +) + _internal_save_mps(target, pm_object(feasible_region(lp)), pm_object(lp)) end @doc raw""" @@ -145,7 +150,6 @@ function load_mps(file::String) end end - @doc raw""" load_lp(file::String) @@ -167,4 +171,3 @@ function load_lp(file::String) error("load_lp: cannot find LP or MILP subobject in polymake object") end end - diff --git a/src/PolyhedralGeometry/PolyhedralComplex/constructors.jl b/src/PolyhedralGeometry/PolyhedralComplex/constructors.jl index 233209cd3abf..74ad0d1382a6 100644 --- a/src/PolyhedralGeometry/PolyhedralComplex/constructors.jl +++ b/src/PolyhedralGeometry/PolyhedralComplex/constructors.jl @@ -5,22 +5,21 @@ ############################################################################### struct PolyhedralComplex{T} <: PolyhedralObject{T} - pm_complex::Polymake.BigObject - parent_field::Field - - PolyhedralComplex{T}(pm::Polymake.BigObject, p::Field) where T<:scalar_types = new{T}(pm, p) - PolyhedralComplex{QQFieldElem}(pm::Polymake.BigObject) = new{QQFieldElem}(pm, QQ) -end + pm_complex::Polymake.BigObject + parent_field::Field + PolyhedralComplex{T}(pm::Polymake.BigObject, p::Field) where {T<:scalar_types} = + new{T}(pm, p) + PolyhedralComplex{QQFieldElem}(pm::Polymake.BigObject) = new{QQFieldElem}(pm, QQ) +end function polyhedral_complex(p::Polymake.BigObject) - T, f = _detect_scalar_and_field(PolyhedralComplex, p) - return PolyhedralComplex{T}(p, f) + T, f = _detect_scalar_and_field(PolyhedralComplex, p) + return PolyhedralComplex{T}(p, f) end pm_object(pc::PolyhedralComplex) = pc.pm_complex - @doc raw""" polyhedral_complex(::T, polyhedra, vr, far_vertices, L) where T<:scalar_types @@ -75,50 +74,73 @@ julia> lineality_dim(PC) 1 ``` """ -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 - ) +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, +) 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) + LM = if isnothing(L) || isempty(L) + zero_matrix(QQ, 0, size(points, 2)) + else + homogenized_matrix(L, 0) + end # 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 + 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) + 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) + 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: guess from input, fallback QQ -polyhedral_complex(polyhedra::IncidenceMatrix, - vr::AbstractCollection[PointVector], - far_vertices::Union{Vector{Int}, Nothing} = nothing, - L::Union{AbstractCollection[RayVector], Nothing} = nothing; - non_redundant::Bool = false) = - polyhedral_complex(_guess_fieldelem_type(vr, L), 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) +polyhedral_complex( + polyhedra::IncidenceMatrix, + vr::AbstractCollection[PointVector], + far_vertices::Union{Vector{Int},Nothing}=nothing, + L::Union{AbstractCollection[RayVector],Nothing}=nothing; + non_redundant::Bool=false, +) = polyhedral_complex( + _guess_fieldelem_type(vr, L), + 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) + return polyhedral_complex(f, polyhedra, vr, far_vertices, L; non_redundant=non_redundant) end @doc raw""" @@ -126,7 +148,9 @@ end Assemble a polyhedral complex from a non-empty list of polyhedra. """ -function polyhedral_complex(polytopes::AbstractVector{Polyhedron{T}}; non_redundant::Bool=false) where T<:scalar_types +function polyhedral_complex( + polytopes::AbstractVector{Polyhedron{T}}; non_redundant::Bool=false +) where {T<:scalar_types} @req length(polytopes) > 0 "list of polytopes must be non-empty" if non_redundant pmcplx = Polymake.fan.complex_from_polytopes(pm_object.(polytopes)...) @@ -137,26 +161,26 @@ function polyhedral_complex(polytopes::AbstractVector{Polyhedron{T}}; non_redund end # shortcut for maximal polytopes of an existing complex -function polyhedral_complex(iter::SubObjectIterator{Polyhedron{T}}) where T<:scalar_types +function polyhedral_complex(iter::SubObjectIterator{Polyhedron{T}}) where {T<:scalar_types} if iter.Acc == _maximal_polyhedron && iter.Obj isa PolyhedralComplex return deepcopy(iter.Obj) end return polyhedral_complex(collect(iter)) end -polyhedral_complex(P::Polyhedron{T}) where T<:scalar_types = polyhedral_complex([P]) +polyhedral_complex(P::Polyhedron{T}) where {T<:scalar_types} = polyhedral_complex([P]) ############################################################################### ############################################################################### ### Display ############################################################################### ############################################################################### -function Base.show(io::IO, PC::PolyhedralComplex{T}) where T<:scalar_types - try - ad = ambient_dim(PC) - print(io, "Polyhedral complex in ambient dimension $(ad)") - T != QQFieldElem && print(io, " with $T type coefficients") - catch e - print(io, "Polyhedral complex without ambient dimension") - end +function Base.show(io::IO, PC::PolyhedralComplex{T}) where {T<:scalar_types} + try + ad = ambient_dim(PC) + print(io, "Polyhedral complex in ambient dimension $(ad)") + T != QQFieldElem && print(io, " with $T type coefficients") + catch e + print(io, "Polyhedral complex without ambient dimension") + end end diff --git a/src/PolyhedralGeometry/PolyhedralComplex/properties.jl b/src/PolyhedralGeometry/PolyhedralComplex/properties.jl index f434cfcad89a..f126316a3d78 100644 --- a/src/PolyhedralGeometry/PolyhedralComplex/properties.jl +++ b/src/PolyhedralGeometry/PolyhedralComplex/properties.jl @@ -63,7 +63,9 @@ vertices(as::Type{PointVector{T}}, PC::PolyhedralComplex{T}) where {T<:scalar_ty _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(U::Type{PointVector{T}}, PC::PolyhedralComplex{T}, i::Base.Integer) where {T<:scalar_types} = point_vector( +_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 @@ -83,7 +85,9 @@ function _all_vertex_indices(P::Polymake.BigObject) return vi end -function _vertex_or_ray_complex(::Type{Union{PointVector{T}, RayVector{T}}}, PC::PolyhedralComplex{T}, 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 ray_vector( @@ -98,10 +102,11 @@ function _vertex_or_ray_complex(::Type{Union{PointVector{T}, RayVector{T}}}, PC: end end -_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( + 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) @@ -160,18 +165,30 @@ 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 + 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) - ) +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(U::Type{RayVector{T}}, PC::PolyhedralComplex{T}) where {T<:scalar_types} = - _rays(U, 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) @@ -208,17 +225,24 @@ 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(::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} = @@ -252,7 +276,9 @@ _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(U::Type{RayVector{T}}, PC::PolyhedralComplex{T}, i::Base.Integer) where {T<:scalar_types} = ray_vector( +_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 @@ -264,11 +290,10 @@ _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{T}, 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))] @@ -460,26 +485,25 @@ function _ith_polyhedron( ::Type{Polyhedron{T}}, 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{PT}(VERTICES = V, LINEALITY_SPACE = L), - 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(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 +_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/PolyhedralComplex/standard_constructions.jl b/src/PolyhedralGeometry/PolyhedralComplex/standard_constructions.jl index 5d9f067170fa..04fc34151c04 100644 --- a/src/PolyhedralGeometry/PolyhedralComplex/standard_constructions.jl +++ b/src/PolyhedralGeometry/PolyhedralComplex/standard_constructions.jl @@ -32,15 +32,18 @@ julia> common_refinement(PC1,PC2) Polyhedral complex in ambient dimension 2 ``` """ -function common_refinement(PC1::PolyhedralComplex{T}, PC2::PolyhedralComplex{T}) where T<:scalar_types - U, f = _promote_scalar_field(coefficient_field(PC1), coefficient_field(PC2)) - pm_PC1 = pm_object(PC1) - pm_PC2 = pm_object(PC2) - result = Polymake.fan.PolyhedralComplex{_scalar_type_to_polymake(T)}(Polymake.fan.common_refinement(pm_PC1,pm_PC2)) - return PolyhedralComplex{T}(result, f) +function common_refinement( + PC1::PolyhedralComplex{T}, PC2::PolyhedralComplex{T} +) where {T<:scalar_types} + U, f = _promote_scalar_field(coefficient_field(PC1), coefficient_field(PC2)) + pm_PC1 = pm_object(PC1) + pm_PC2 = pm_object(PC2) + result = Polymake.fan.PolyhedralComplex{_scalar_type_to_polymake(T)}( + Polymake.fan.common_refinement(pm_PC1, pm_PC2) + ) + return PolyhedralComplex{T}(result, f) end - @doc raw""" k_skeleton(PC::PolyhedralComplex,k::Int) @@ -65,10 +68,10 @@ julia> k_skeleton(PC1,1) Polyhedral complex in ambient dimension 2 ``` """ -function k_skeleton(PC::PolyhedralComplex{T},k::Int) where T<:scalar_types - pm_PC = pm_object(PC) - ksk = Polymake.fan.PolyhedralComplex{_scalar_type_to_polymake(T)}(Polymake.fan.k_skeleton(pm_PC,k)) - return PolyhedralComplex{T}(ksk, coefficient_field(PC)) +function k_skeleton(PC::PolyhedralComplex{T}, k::Int) where {T<:scalar_types} + pm_PC = pm_object(PC) + ksk = Polymake.fan.PolyhedralComplex{_scalar_type_to_polymake(T)}( + Polymake.fan.k_skeleton(pm_PC, k) + ) + return PolyhedralComplex{T}(ksk, coefficient_field(PC)) end - - diff --git a/src/PolyhedralGeometry/PolyhedralFan/constructors.jl b/src/PolyhedralGeometry/PolyhedralFan/constructors.jl index 23e261232898..f54219634937 100644 --- a/src/PolyhedralGeometry/PolyhedralFan/constructors.jl +++ b/src/PolyhedralGeometry/PolyhedralFan/constructors.jl @@ -5,21 +5,21 @@ ############################################################################### struct PolyhedralFan{T} <: PolyhedralObject{T} - pm_fan::Polymake.BigObject - parent_field::Field + pm_fan::Polymake.BigObject + parent_field::Field - PolyhedralFan{T}(pm::Polymake.BigObject, f::Field) where T<:scalar_types = new{T}(pm, f) - PolyhedralFan{QQFieldElem}(pm::Polymake.BigObject) = new{QQFieldElem}(pm, QQ) + PolyhedralFan{T}(pm::Polymake.BigObject, f::Field) where {T<:scalar_types} = new{T}(pm, f) + PolyhedralFan{QQFieldElem}(pm::Polymake.BigObject) = new{QQFieldElem}(pm, QQ) end -const _FanLikeType = Union{NormalToricVarietyType, PolyhedralFan} -const _FanLikeTypeQQ = Union{NormalToricVarietyType, PolyhedralFan{QQFieldElem}} +const _FanLikeType = Union{NormalToricVarietyType,PolyhedralFan} +const _FanLikeTypeQQ = Union{NormalToricVarietyType,PolyhedralFan{QQFieldElem}} # Automatic detection of corresponding OSCAR scalar type; # Avoid, if possible, to increase type stability function polyhedral_fan(p::Polymake.BigObject) - T, f = _detect_scalar_and_field(PolyhedralFan, p) - return PolyhedralFan{T}(p, f) + T, f = _detect_scalar_and_field(PolyhedralFan, p) + return PolyhedralFan{T}(p, f) end @doc raw""" @@ -64,11 +64,13 @@ julia> lineality_dim(PF) 1 ``` """ -function polyhedral_fan(f::scalar_type_or_field, - Incidence::IncidenceMatrix, - Rays::AbstractCollection[RayVector], - LS::Union{AbstractCollection[RayVector], Nothing}; - non_redundant::Bool = false) +function polyhedral_fan( + f::scalar_type_or_field, + Incidence::IncidenceMatrix, + Rays::AbstractCollection[RayVector], + LS::Union{AbstractCollection[RayVector],Nothing}; + non_redundant::Bool=false, +) parent_field, scalar_type = _determine_parent_and_scalar(f, Rays, LS) RM = unhomogenized_matrix(Rays) if isnothing(LS) @@ -77,22 +79,36 @@ function polyhedral_fan(f::scalar_type_or_field, 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) + 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) + 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::scalar_type_or_field, Incidence::IncidenceMatrix, Rays::AbstractCollection[RayVector]; non_redundant::Bool = false) = polyhedral_fan(f, Incidence, Rays, nothing; non_redundant) -polyhedral_fan(Incidence::IncidenceMatrix, Rays::AbstractCollection[RayVector], LS::Union{AbstractCollection[RayVector], Nothing}; non_redundant::Bool = false) = polyhedral_fan(_guess_fieldelem_type(Rays, LS), Incidence, Rays, LS; non_redundant) -polyhedral_fan(Incidence::IncidenceMatrix, Rays::AbstractCollection[RayVector]; non_redundant::Bool = false) = polyhedral_fan(_guess_fieldelem_type(Rays), Incidence, Rays; non_redundant) +polyhedral_fan( + f::scalar_type_or_field, + Incidence::IncidenceMatrix, + Rays::AbstractCollection[RayVector]; + non_redundant::Bool=false, +) = polyhedral_fan(f, Incidence, Rays, nothing; non_redundant) +polyhedral_fan( + Incidence::IncidenceMatrix, + Rays::AbstractCollection[RayVector], + LS::Union{AbstractCollection[RayVector],Nothing}; + non_redundant::Bool=false, +) = polyhedral_fan(_guess_fieldelem_type(Rays, LS), Incidence, Rays, LS; non_redundant) +polyhedral_fan( + Incidence::IncidenceMatrix, Rays::AbstractCollection[RayVector]; non_redundant::Bool=false +) = polyhedral_fan(_guess_fieldelem_type(Rays), Incidence, Rays; non_redundant) """ pm_object(PF::PolyhedralFan) @@ -106,7 +122,9 @@ pm_object(PF::PolyhedralFan) = PF.pm_fan Assemble a polyhedral fan from a non-empty list of cones. """ -function polyhedral_fan(cones::AbstractVector{Cone{T}}; non_redundant::Bool=false) where T<:scalar_types +function polyhedral_fan( + cones::AbstractVector{Cone{T}}; non_redundant::Bool=false +) where {T<:scalar_types} @req length(cones) > 0 "list of cones must be non-empty" if non_redundant pmfan = Polymake.fan.fan_from_cones(pm_object.(cones)...) @@ -117,15 +135,20 @@ function polyhedral_fan(cones::AbstractVector{Cone{T}}; non_redundant::Bool=fals end #Same construction for when the user gives Matrix{Bool} as incidence matrix -polyhedral_fan(f::scalar_type_or_field, Incidence::Matrix{Bool}, Rays::AbstractCollection[RayVector], LS::AbstractCollection[RayVector]) = - polyhedral_fan(f, IncidenceMatrix(Polymake.IncidenceMatrix(Incidence)), Rays, LS) -polyhedral_fan(f::scalar_type_or_field, Incidence::Matrix{Bool}, Rays::AbstractCollection[RayVector]) = - polyhedral_fan(f, IncidenceMatrix(Polymake.IncidenceMatrix(Incidence)), Rays) - -polyhedral_fan(C::Cone{T}) where T<:scalar_types = polyhedral_fan([C]) +polyhedral_fan( + f::scalar_type_or_field, + Incidence::Matrix{Bool}, + Rays::AbstractCollection[RayVector], + LS::AbstractCollection[RayVector], +) = polyhedral_fan(f, IncidenceMatrix(Polymake.IncidenceMatrix(Incidence)), Rays, LS) +polyhedral_fan( + f::scalar_type_or_field, Incidence::Matrix{Bool}, Rays::AbstractCollection[RayVector] +) = polyhedral_fan(f, IncidenceMatrix(Polymake.IncidenceMatrix(Incidence)), Rays) + +polyhedral_fan(C::Cone{T}) where {T<:scalar_types} = polyhedral_fan([C]) # shortcut for maximal cones of an existing fan -function polyhedral_fan(iter::SubObjectIterator{Cone{T}}) where T<:scalar_types +function polyhedral_fan(iter::SubObjectIterator{Cone{T}}) where {T<:scalar_types} if iter.Acc == _maximal_cone && iter.Obj isa PolyhedralFan return deepcopy(iter.Obj) end @@ -152,25 +175,37 @@ parent `Field`. group acting on the rays of the fan. """ -function polyhedral_fan_from_rays_action(f::scalar_type_or_field, Rays::AbstractCollection[RayVector], MC_reps::IncidenceMatrix, perms::AbstractVector{PermGroupElem}) +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))) + 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) - +polyhedral_fan_from_rays_action( + Rays::AbstractCollection[RayVector], + MC_reps::IncidenceMatrix, + perms::AbstractVector{PermGroupElem}, +) = polyhedral_fan_from_rays_action(QQFieldElem, Rays, MC_reps, perms) ############################################################################### ############################################################################### ### Display ############################################################################### ############################################################################### -function Base.show(io::IO, PF::PolyhedralFan{T}) where T<:scalar_types - print(io, "Polyhedral fan in ambient dimension $(ambient_dim(PF))") - T != QQFieldElem && print(io, " with $T type coefficients") +function Base.show(io::IO, PF::PolyhedralFan{T}) where {T<:scalar_types} + print(io, "Polyhedral fan in ambient dimension $(ambient_dim(PF))") + T != QQFieldElem && print(io, " with $T type coefficients") end diff --git a/src/PolyhedralGeometry/PolyhedralFan/properties.jl b/src/PolyhedralGeometry/PolyhedralFan/properties.jl index 13da49581018..d853df94823c 100644 --- a/src/PolyhedralGeometry/PolyhedralFan/properties.jl +++ b/src/PolyhedralGeometry/PolyhedralFan/properties.jl @@ -43,20 +43,26 @@ julia> matrix(QQ, rays(NF)) [ 0 0 -1] ``` """ -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, _n_rays(PF)) +rays(PF::_FanLikeType) = + if lineality_dim(PF) == 0 + _rays(PF) + else + _empty_subobjectiterator(RayVector{_get_scalar_type(PF)}, PF) + end +_rays(PF::_FanLikeType) = + SubObjectIterator{RayVector{_get_scalar_type(PF)}}(PF, _ray_fan, _n_rays(PF)) _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 +_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)) - @doc raw""" rays_modulo_lineality(as, F::PolyhedralFan) @@ -89,22 +95,30 @@ julia> rays(NF) 0-element SubObjectIterator{RayVector{QQFieldElem}} ``` """ -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) - ) +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(::Type{<:RayVector}, F::_FanLikeType) = _rays(F) - - @doc raw""" maximal_cones(PF::PolyhedralFan) @@ -128,7 +142,8 @@ julia> for c in maximal_cones(PF) 4 ``` """ -maximal_cones(PF::_FanLikeType) = SubObjectIterator{Cone{_get_scalar_type(PF)}}(PF, _maximal_cone, n_maximal_cones(PF)) +maximal_cones(PF::_FanLikeType) = + SubObjectIterator{Cone{_get_scalar_type(PF)}}(PF, _maximal_cone, n_maximal_cones(PF)) _ray_indices(::Val{_maximal_cone}, PF::_FanLikeType) = pm_object(PF).MAXIMAL_CONES @@ -161,23 +176,31 @@ julia> cones(PF, 2) ``` """ function cones(PF::_FanLikeType, cone_dim::Int) - l = cone_dim - length(lineality_space(PF)) - l < 1 && return nothing - return SubObjectIterator{Cone{_get_scalar_type(PF)}}(PF, _cone_of_dim, size(Polymake.fan.cones_of_dim(pm_object(PF), l), 1), (c_dim = l,)) + l = cone_dim - length(lineality_space(PF)) + l < 1 && return nothing + return SubObjectIterator{Cone{_get_scalar_type(PF)}}( + PF, _cone_of_dim, size(Polymake.fan.cones_of_dim(pm_object(PF), l), 1), (c_dim=l,) + ) end -function _cone_of_dim(::Type{Cone{T}}, PF::_FanLikeType, i::Base.Integer; c_dim::Int = 0) where T<:scalar_types - R = pm_object(PF).RAYS[collect(Polymake.row(Polymake.fan.cones_of_dim(pm_object(PF), c_dim), i)), :] +function _cone_of_dim( + ::Type{Cone{T}}, PF::_FanLikeType, i::Base.Integer; c_dim::Int=0 +) where {T<:scalar_types} + 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)) + 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) +_ray_indices(::Val{_cone_of_dim}, PF::_FanLikeType; c_dim::Int=0) = + Polymake.fan.cones_of_dim(pm_object(PF), c_dim) _incidencematrix(::Val{_cone_of_dim}) = _ray_indices - @doc raw""" cones(PF::PolyhedralFan) @@ -204,13 +227,12 @@ julia> cones(PF) function cones(PF::_FanLikeType) pmo = pm_object(PF) ncones = pmo.HASSE_DIAGRAM.N_NODES - cones = [Polymake._get_entry(pmo.HASSE_DIAGRAM.FACES,i) for i in 0:(ncones-1)]; - cones = filter(x->!(-1 in x) && length(x)>0, cones); - cones = [Polymake.to_one_based_indexing(x) for x in cones]; + cones = [Polymake._get_entry(pmo.HASSE_DIAGRAM.FACES, i) for i in 0:(ncones - 1)] + cones = filter(x -> !(-1 in x) && length(x) > 0, cones) + cones = [Polymake.to_one_based_indexing(x) for x in cones] return IncidenceMatrix([Vector{Int}(x) for x in cones]) end - ############################################################################### ############################################################################### ### Access properties @@ -306,7 +328,6 @@ julia> n_rays(face_fan(cube(3))) n_rays(PF::_FanLikeType) = lineality_dim(PF) == 0 ? _n_rays(PF) : 0 _n_rays(PF::_FanLikeType) = pm_object(PF).N_RAYS::Int - @doc raw""" f_vector(PF::PolyhedralFan) @@ -338,12 +359,11 @@ julia> f_vector(nfc) ``` """ function f_vector(PF::_FanLikeType) - pmf = pm_object(PF) - ldim = pmf.LINEALITY_DIM - return Vector{ZZRingElem}(vcat(fill(0,ldim),pmf.F_VECTOR)) + pmf = pm_object(PF) + ldim = pmf.LINEALITY_DIM + return Vector{ZZRingElem}(vcat(fill(0, ldim), pmf.F_VECTOR)) end - @doc raw""" lineality_dim(PF::PolyhedralFan) @@ -394,11 +414,16 @@ julia> lineality_space(PF) [1, 0] ``` """ -lineality_space(PF::_FanLikeType) = SubObjectIterator{RayVector{_get_scalar_type(PF)}}(PF, _lineality_fan, lineality_dim(PF)) +lineality_space(PF::_FanLikeType) = + SubObjectIterator{RayVector{_get_scalar_type(PF)}}(PF, _lineality_fan, lineality_dim(PF)) -_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 +_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 +_generator_matrix(::Val{_lineality_fan}, PF::_FanLikeType; homogenized=false) = + homogenized ? homogenize(pm_object(PF).LINEALITY_SPACE, 0) : pm_object(PF).LINEALITY_SPACE _matrix_for_polymake(::Val{_lineality_fan}) = _generator_matrix @@ -431,7 +456,6 @@ julia> lineality_dim(nf) """ is_pointed(PF::_FanLikeType) = pm_object(PF).POINTED::Bool - @doc raw""" is_smooth(PF::PolyhedralFan{QQFieldElem}) @@ -465,7 +489,6 @@ false """ is_regular(PF::_FanLikeType) = pm_object(PF).REGULAR::Bool - @doc raw""" is_pure(PF::PolyhedralFan) @@ -481,7 +504,6 @@ false """ is_pure(PF::_FanLikeType) = pm_object(PF).PURE::Bool - @doc raw""" is_fulldimensional(PF::PolyhedralFan) @@ -498,7 +520,6 @@ true """ is_fulldimensional(PF::_FanLikeType) = pm_object(PF).FULL_DIM::Bool - @doc raw""" is_complete(PF::PolyhedralFan) @@ -514,7 +535,6 @@ true """ is_complete(PF::_FanLikeType) = pm_object(PF).COMPLETE::Bool - @doc raw""" is_simplicial(PF::PolyhedralFan) @@ -533,7 +553,6 @@ false """ is_simplicial(PF::_FanLikeType) = pm_object(PF).SIMPLICIAL::Bool - ############################################################################### ## Primitive collections ############################################################################### @@ -551,10 +570,8 @@ julia> primitive_collections(normal_fan(simplex(3))) ``` """ function primitive_collections(PF::_FanLikeType) - @req is_simplicial(PF) "PolyhedralFan must be simplicial." - I = ray_indices(maximal_cones(PF)) - K = simplicial_complex(I) - return minimal_nonfaces(K) + @req is_simplicial(PF) "PolyhedralFan must be simplicial." + I = ray_indices(maximal_cones(PF)) + K = simplicial_complex(I) + return minimal_nonfaces(K) end - - diff --git a/src/PolyhedralGeometry/PolyhedralFan/standard_constructions.jl b/src/PolyhedralGeometry/PolyhedralFan/standard_constructions.jl index 398f4d313249..37d4abca7e26 100644 --- a/src/PolyhedralGeometry/PolyhedralFan/standard_constructions.jl +++ b/src/PolyhedralGeometry/PolyhedralFan/standard_constructions.jl @@ -25,10 +25,10 @@ julia> rays(NF) [0, 0, -1] ``` """ -function normal_fan(P::Polyhedron{T}) where T<:scalar_types - pmp = pm_object(P) - pmnf = Polymake.fan.normal_fan(pmp) - return PolyhedralFan{T}(pmnf, coefficient_field(P)) +function normal_fan(P::Polyhedron{T}) where {T<:scalar_types} + pmp = pm_object(P) + pmnf = Polymake.fan.normal_fan(pmp) + return PolyhedralFan{T}(pmnf, coefficient_field(P)) end """ @@ -50,18 +50,16 @@ julia> n_maximal_cones(FF) == n_facets(C) true ``` """ -function face_fan(P::Polyhedron{T}) where T<:scalar_types - pmp = pm_object(P) - pmff = Polymake.fan.face_fan(pmp) - return PolyhedralFan{T}(pmff, coefficient_field(P)) +function face_fan(P::Polyhedron{T}) where {T<:scalar_types} + pmp = pm_object(P) + pmff = Polymake.fan.face_fan(pmp) + return PolyhedralFan{T}(pmff, coefficient_field(P)) end - ############################################################################### ## Star subdivision ############################################################################### - @doc raw""" star_subdivision(PF::PolyhedralFan, new_ray::AbstractVector{<:IntegerUnion}) @@ -98,25 +96,25 @@ julia> ray_indices(maximal_cones(star)) ``` """ function star_subdivision(Sigma::_FanLikeType, new_ray::AbstractVector{<:IntegerUnion}) - + # Check if new_ray is primitive in ZZ^d, i.e. gcd(new_ray)==1. @req ambient_dim(Sigma) == length(new_ray) "New ray is not $(ambient_dim(Sigma))-dimensional" @req gcd(new_ray) == 1 "The new ray r is not a primitive element of the lattice Z^d with d = length(r)" @req lineality_dim(Sigma) == 0 "star_subdivision does not work for polyhedral fans with lineality." - + old_rays = matrix(ZZ, rays(Sigma)) # In case the new ray is an old ray. - new_ray_index = findfirst(i->vec(old_rays[i,:])==new_ray, 1:nrows(old_rays)) + new_ray_index = findfirst(i -> vec(old_rays[i, :]) == new_ray, 1:nrows(old_rays)) new_rays = old_rays if isnothing(new_ray_index) new_rays = vcat(old_rays, matrix(ZZ, [new_ray])) - new_ray_index = n_rays(Sigma)+1 + new_ray_index = n_rays(Sigma) + 1 end mc_old = maximal_cones(IncidenceMatrix, Sigma) - + facet_normals = pm_object(Sigma).FACET_NORMALS refinable_cones = _get_maximal_cones_containing_vector(Sigma, new_ray) - @req length(refinable_cones)>0 "$new_ray not contained in support of fan." + @req length(refinable_cones) > 0 "$new_ray not contained in support of fan." new_cones = _get_refinable_facets(Sigma, new_ray, refinable_cones, facet_normals, mc_old) for nc in new_cones push!(nc, new_ray_index) @@ -126,23 +124,40 @@ function star_subdivision(Sigma::_FanLikeType, new_ray::AbstractVector{<:Integer push!(new_cones, Vector{Int}(Polymake.row(mc_old, i))) end end - - return polyhedral_fan(coefficient_field(Sigma), IncidenceMatrix([nc for nc in new_cones]), new_rays; non_redundant=true) + + return polyhedral_fan( + coefficient_field(Sigma), + IncidenceMatrix([nc for nc in new_cones]), + new_rays; + non_redundant=true, + ) end -function _get_refinable_facets(Sigma::_FanLikeType, new_ray::AbstractVector{<:IntegerUnion}, refinable_cones::Vector{Int}, facet_normals::AbstractMatrix, mc_old::IncidenceMatrix) +function _get_refinable_facets( + Sigma::_FanLikeType, + new_ray::AbstractVector{<:IntegerUnion}, + refinable_cones::Vector{Int}, + facet_normals::AbstractMatrix, + mc_old::IncidenceMatrix, +) new_cones = Vector{Int}[] v_facet_signs = _facet_signs(facet_normals, new_ray) R = pm_object(Sigma).RAYS hd = pm_object(Sigma).HASSE_DIAGRAM hd_graph = Graph{Directed}(hd.ADJACENCY) - hd_maximal_cones = inneighbors(hd_graph, hd.TOP_NODE+1) + hd_maximal_cones = inneighbors(hd_graph, hd.TOP_NODE + 1) Sfacets = pm_object(Sigma).MAXIMAL_CONES_FACETS for mc_index in refinable_cones mc_indices = Polymake.row(mc_old, mc_index) - mc_facet_indices = Polymake.findnz(Sfacets[mc_index,:])[1] - mc_hd_index = hd_maximal_cones[findfirst(i -> Polymake.to_one_based_indexing(Polymake._get_entry(hd.FACES, i-1)) == mc_indices, hd_maximal_cones)] - refinable_facets = _get_refinable_facets_of_cone(mc_hd_index, facet_normals, hd, hd_graph, v_facet_signs, R, mc_facet_indices) + mc_facet_indices = Polymake.findnz(Sfacets[mc_index, :])[1] + mc_hd_index = hd_maximal_cones[findfirst( + i -> + Polymake.to_one_based_indexing(Polymake._get_entry(hd.FACES, i - 1)) == mc_indices, + hd_maximal_cones, + )] + refinable_facets = _get_refinable_facets_of_cone( + mc_hd_index, facet_normals, hd, hd_graph, v_facet_signs, R, mc_facet_indices + ) append!(new_cones, refinable_facets) # If all facets contain new_ray, then the current maximal cone is just # new_ray. @@ -151,15 +166,25 @@ function _get_refinable_facets(Sigma::_FanLikeType, new_ray::AbstractVector{<:In return unique(new_cones) end -function _get_refinable_facets_of_cone(mc_hd_index::Int, facet_normals::AbstractMatrix, hd, hd_graph::Graph{Directed}, v_facet_signs::AbstractVector, R::AbstractMatrix, mcfi::Vector{Int}) +function _get_refinable_facets_of_cone( + mc_hd_index::Int, + facet_normals::AbstractMatrix, + hd, + hd_graph::Graph{Directed}, + v_facet_signs::AbstractVector, + R::AbstractMatrix, + mcfi::Vector{Int}, +) refinable_facets = Vector{Int}[] - mc_indices = Vector{Int}(Polymake.to_one_based_indexing(Polymake._get_entry(hd.FACES, mc_hd_index-1))) + mc_indices = Vector{Int}( + Polymake.to_one_based_indexing(Polymake._get_entry(hd.FACES, mc_hd_index - 1)) + ) for fc_index in inneighbors(hd_graph, mc_hd_index) - fc_indices = Polymake.to_one_based_indexing(Polymake._get_entry(hd.FACES, fc_index-1)) + fc_indices = Polymake.to_one_based_indexing(Polymake._get_entry(hd.FACES, fc_index - 1)) length(fc_indices) > 0 || return refinable_facets # The only facet was 0 - inner_ray = sum([R[i,:] for i in fc_indices]) + inner_ray = sum([R[i, :] for i in fc_indices]) fc_facet_signs = _facet_signs(facet_normals, inner_ray) - if(!_check_containment_via_facet_signs(v_facet_signs[mcfi], fc_facet_signs[mcfi])) + if (!_check_containment_via_facet_signs(v_facet_signs[mcfi], fc_facet_signs[mcfi])) push!(refinable_facets, Vector{Int}(fc_indices)) end end @@ -167,9 +192,10 @@ function _get_refinable_facets_of_cone(mc_hd_index::Int, facet_normals::Abstract end # FIXME: Small workaround, since sign does not work for polymake types. -_int_sign(e) = e>0 ? 1 : (e<0 ? -1 : 0) -_facet_signs(F::AbstractMatrix, v::AbstractVector{<:IntegerUnion}) = [_int_sign(e) for e in (F*Polymake.Vector{Polymake.Integer}(v))] -_facet_signs(F::AbstractMatrix, v::AbstractVector) = [_int_sign(e) for e in (F*v)] +_int_sign(e) = e > 0 ? 1 : (e < 0 ? -1 : 0) +_facet_signs(F::AbstractMatrix, v::AbstractVector{<:IntegerUnion}) = + [_int_sign(e) for e in (F * Polymake.Vector{Polymake.Integer}(v))] +_facet_signs(F::AbstractMatrix, v::AbstractVector) = [_int_sign(e) for e in (F * v)] function _check_containment_via_facet_signs(smaller::Vector{Int}, bigger::Vector{Int}) for a in zip(smaller, bigger) p = prod(a) @@ -180,14 +206,15 @@ function _check_containment_via_facet_signs(smaller::Vector{Int}, bigger::Vector end return true end -function _get_maximal_cones_containing_vector(Sigma::_FanLikeType, v::AbstractVector{<:IntegerUnion}) +function _get_maximal_cones_containing_vector( + Sigma::_FanLikeType, v::AbstractVector{<:IntegerUnion} +) # Make sure these are computed as otherwise performance degrades. pm_object(Sigma).FACET_NORMALS pm_object(Sigma).MAXIMAL_CONES_FACETS return findall(mc -> v in mc, maximal_cones(Sigma)) end - @doc raw""" star_subdivision(PF::PolyhedralFan, n::Int) @@ -223,20 +250,22 @@ function star_subdivision(Sigma::_FanLikeType, n::Int) tau = Polymake.row(cones_Sigma, n) @req length(tau) > 1 "Cannot subdivide cone $n as it is generated by a single ray" R = matrix(ZZ, rays(Sigma)) - newray = vec(sum([R[i,:] for i in tau])) + newray = vec(sum([R[i, :] for i in tau])) newray = newray ./ gcd(newray) return star_subdivision(Sigma, newray) end - - @doc raw""" transform(F::_FanLikeType, A::AbstractMatrix; check=true) where Get the image of a fan `F` under the matrix `A`. The default is to check whether the images of the maximal cones really form a fan. """ -function transform(F::_FanLikeType, A::Union{AbstractMatrix{<:Union{Number, FieldElem}}, MatElem{<:FieldElem}}; check::Bool=true) +function transform( + F::_FanLikeType, + A::Union{AbstractMatrix{<:Union{Number,FieldElem}},MatElem{<:FieldElem}}; + check::Bool=true, +) @req ncols(A) == ambient_dim(F) "Incompatible dimension of fan and transformation matrix" OT = _scalar_type_to_polymake(_get_scalar_type(F)) return _transform(F, Polymake.Matrix{OT}(A); check) @@ -253,12 +282,11 @@ function _transform(F::_FanLikeType, A::Polymake.Matrix; check::Bool) result = Polymake.fan.check_fan(R, MC, opt) return FT(result, coefficient_field(F)) else - result = Polymake.fan.PolyhedralFan{OT}(RAYS=R, LINEALITY_SPACE=L, MAXIMAL_CONES=MC) + result = Polymake.fan.PolyhedralFan{OT}(; RAYS=R, LINEALITY_SPACE=L, MAXIMAL_CONES=MC) return FT(result, coefficient_field(F)) end end - ############################################################################### ## Cartesian/Direct product ############################################################################### @@ -275,6 +303,6 @@ Polyhedral fan in ambient dimension 5 ``` """ function Base.:*(PF1::PolyhedralFan{QQFieldElem}, PF2::PolyhedralFan{QQFieldElem}) - prod = Polymake.fan.product(pm_object(PF1), pm_object(PF2)) - return PolyhedralFan{QQFieldElem}(prod, QQ) + prod = Polymake.fan.product(pm_object(PF1), pm_object(PF2)) + return PolyhedralFan{QQFieldElem}(prod, QQ) end diff --git a/src/PolyhedralGeometry/PolyhedralGeometry.jl b/src/PolyhedralGeometry/PolyhedralGeometry.jl index b08e782035ba..d506337241e9 100644 --- a/src/PolyhedralGeometry/PolyhedralGeometry.jl +++ b/src/PolyhedralGeometry/PolyhedralGeometry.jl @@ -1,4 +1,4 @@ -const AnyVecOrMat = Union{MatElem, AbstractVecOrMat} +const AnyVecOrMat = Union{MatElem,AbstractVecOrMat} include("helpers.jl") include("iterators.jl") @@ -26,7 +26,6 @@ include("visualization.jl") include("solving_integrally.jl") include("triangulations.jl") - # Some temporary aliases to avoid breaking all current PRs pm_cone(C::Cone) = pm_object(C) pm_fan(PF::PolyhedralFan) = pm_object(PF) diff --git a/src/PolyhedralGeometry/Polyhedron/constructors.jl b/src/PolyhedralGeometry/Polyhedron/constructors.jl index 25c7effb3a90..ab53214de77e 100644 --- a/src/PolyhedralGeometry/Polyhedron/constructors.jl +++ b/src/PolyhedralGeometry/Polyhedron/constructors.jl @@ -5,23 +5,23 @@ ############################################################################### struct Polyhedron{T<:scalar_types} <: PolyhedralObject{T} #a real polymake polyhedron - pm_polytope::Polymake.BigObject - parent_field::Field + pm_polytope::Polymake.BigObject + parent_field::Field - # only allowing scalar_types; - # can be improved by testing if the template type of the `BigObject` corresponds to `T` + # only allowing scalar_types; + # can be improved by testing if the template type of the `BigObject` corresponds to `T` - @doc raw""" - Polyhedron{T}(P::Polymake.BigObject, F::Field) where T<:scalar_types + @doc raw""" + Polyhedron{T}(P::Polymake.BigObject, F::Field) where T<:scalar_types - Construct a `Polyhedron` corresponding to a `Polymake.BigObject` of type `Polytope` with scalars from `Field` `F`. - """ - Polyhedron{T}(p::Polymake.BigObject, f::Field) where T<:scalar_types = new{T}(p, f) - Polyhedron{QQFieldElem}(p::Polymake.BigObject) = new{QQFieldElem}(p, QQ) + Construct a `Polyhedron` corresponding to a `Polymake.BigObject` of type `Polytope` with scalars from `Field` `F`. + """ + Polyhedron{T}(p::Polymake.BigObject, f::Field) where {T<:scalar_types} = new{T}(p, f) + Polyhedron{QQFieldElem}(p::Polymake.BigObject) = new{QQFieldElem}(p, QQ) end # default scalar type: guess from input and fall back to `QQFieldElem` -polyhedron(A, b) = polyhedron(_guess_fieldelem_type(A,b), A, b) +polyhedron(A, b) = polyhedron(_guess_fieldelem_type(A, b), A, b) polyhedron(A) = polyhedron(_guess_fieldelem_type(A), A) @doc raw""" @@ -31,9 +31,11 @@ Construct a `Polyhedron` corresponding to a `Polymake.BigObject` of type `Polyto """ function polyhedron(p::Polymake.BigObject) T, f = _detect_scalar_and_field(Polyhedron, p) - if T == EmbeddedNumFieldElem{AbsSimpleNumFieldElem} && Hecke.is_quadratic_type(number_field(f))[1] && Polymake.bigobject_eltype(p) == "QuadraticExtension" + if T == EmbeddedNumFieldElem{AbsSimpleNumFieldElem} && + Hecke.is_quadratic_type(number_field(f))[1] && + Polymake.bigobject_eltype(p) == "QuadraticExtension" p = _polyhedron_qe_to_on(p, f) - end + end return Polyhedron{T}(p, f) end @@ -60,17 +62,33 @@ julia> polyhedron(A,b) Polyhedron in ambient dimension 2 ``` """ -polyhedron(f::scalar_type_or_field, A::AnyVecOrMat, b::AbstractVector; non_redundant::Bool = false) = polyhedron(f, (A, b); non_redundant = non_redundant) - -polyhedron(f::scalar_type_or_field, A::AbstractVector, b::Any; non_redundant::Bool = false) = polyhedron(f, ([A], [b]); non_redundant = non_redundant) - -polyhedron(f::scalar_type_or_field, A::AbstractVector, b::AbstractVector; non_redundant::Bool = false) = polyhedron(f, ([A], b); non_redundant = non_redundant) - -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::scalar_type_or_field, A::AbstractVector{<:AbstractVector}, b::AbstractVector; non_redundant::Bool = false) = polyhedron(f, (A, b); non_redundant = non_redundant) - -polyhedron(f::scalar_type_or_field, A::AnyVecOrMat, b::Any; non_redundant::Bool = false) = polyhedron(f, A, [b]; non_redundant = non_redundant) +polyhedron( + f::scalar_type_or_field, A::AnyVecOrMat, b::AbstractVector; non_redundant::Bool=false +) = polyhedron(f, (A, b); non_redundant=non_redundant) + +polyhedron(f::scalar_type_or_field, A::AbstractVector, b::Any; non_redundant::Bool=false) = + polyhedron(f, ([A], [b]); non_redundant=non_redundant) + +polyhedron( + f::scalar_type_or_field, A::AbstractVector, b::AbstractVector; non_redundant::Bool=false +) = polyhedron(f, ([A], b); non_redundant=non_redundant) + +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::scalar_type_or_field, + A::AbstractVector{<:AbstractVector}, + b::AbstractVector; + non_redundant::Bool=false, +) = polyhedron(f, (A, b); non_redundant=non_redundant) + +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 @@ -109,20 +127,39 @@ julia> vertices(P) [0, 0] ``` """ -function polyhedron(f::scalar_type_or_field, I::Union{Nothing, AbstractCollection[AffineHalfspace]}, E::Union{Nothing, AbstractCollection[AffineHyperplane]} = nothing; non_redundant::Bool = false) +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) + EM = if isnothing(E) || _isempty_halfspace(E) + Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(IM, 2)) + else + affine_matrix_for_polymake(E) + end 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) + 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) + 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 @@ -133,15 +170,15 @@ Get the underlying polymake `Polytope`. """ pm_object(P::Polyhedron) = P.pm_polytope -function ==(P0::Polyhedron{T}, P1::Polyhedron{T}) where T<:scalar_types - Polymake.polytope.equal_polyhedra(pm_object(P0), pm_object(P1)) +function ==(P0::Polyhedron{T}, P1::Polyhedron{T}) where {T<:scalar_types} + Polymake.polytope.equal_polyhedra(pm_object(P0), pm_object(P1)) end # For a proper hash function for cones we should use a "normal form", # which would require a potentially expensive convex hull computation # (and even that is not enough). But hash methods should be fast, so we # just consider the ambient dimension and the precise type of the polyhedron. -function Base.hash(x::T, h::UInt) where {T <: Polyhedron} +function Base.hash(x::T, h::UInt) where {T<:Polyhedron} h = hash(ambient_dim(x), h) h = hash(T, h) return h @@ -210,44 +247,69 @@ julia> XA = convex_hull(V, R, L) Polyhedron in ambient dimension 2 ``` """ -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) +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) + lineality = if isnothing(L) || isempty(L) + zero_matrix(QQ, 0, size(points, 2)) + else + homogenized_matrix(L, 0) + end # 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) + 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) + 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(_guess_fieldelem_type(V,R,L), V, R, L; non_redundant=non_redundant) +convex_hull( + V::AbstractCollection[PointVector], + R::Union{AbstractCollection[RayVector],Nothing}=nothing, + L::Union{AbstractCollection[RayVector],Nothing}=nothing; + non_redundant::Bool=false, +) = convex_hull(_guess_fieldelem_type(V, R, L), V, R, L; non_redundant=non_redundant) ############################################################################### ############################################################################### ### Display ############################################################################### ############################################################################### -function Base.show(io::IO, P::Polyhedron{T}) where T<:scalar_types - pm_P = pm_object(P) - known_to_be_bounded = false - # if the vertices and rays are known, then it is easy to check if the polyhedron is bounded - if Polymake.exists(pm_P, "VERTICES") || Polymake.exists(pm_P, "BOUNDED") - known_to_be_bounded = pm_P.BOUNDED - end - poly_word = known_to_be_bounded ? "Polytope" : "Polyhedron" - try - ad = ambient_dim(P) - print(io, "$(poly_word) in ambient dimension $(ad)") - T != QQFieldElem && print(io, " with $T type coefficients") - catch e - print(io, "$(poly_word) without ambient dimension") - end +function Base.show(io::IO, P::Polyhedron{T}) where {T<:scalar_types} + pm_P = pm_object(P) + known_to_be_bounded = false + # if the vertices and rays are known, then it is easy to check if the polyhedron is bounded + if Polymake.exists(pm_P, "VERTICES") || Polymake.exists(pm_P, "BOUNDED") + known_to_be_bounded = pm_P.BOUNDED + end + poly_word = known_to_be_bounded ? "Polytope" : "Polyhedron" + try + ad = ambient_dim(P) + print(io, "$(poly_word) in ambient dimension $(ad)") + T != QQFieldElem && print(io, " with $T type coefficients") + catch e + print(io, "$(poly_word) without ambient dimension") + end end Polymake.visual(P::Polyhedron; opts...) = Polymake.visual(pm_object(P); opts...) diff --git a/src/PolyhedralGeometry/SubdivisionOfPoints/constructors.jl b/src/PolyhedralGeometry/SubdivisionOfPoints/constructors.jl index 8f4d630f9573..ce24c51966bc 100644 --- a/src/PolyhedralGeometry/SubdivisionOfPoints/constructors.jl +++ b/src/PolyhedralGeometry/SubdivisionOfPoints/constructors.jl @@ -5,14 +5,14 @@ ############################################################################### struct SubdivisionOfPoints{T} <: PolyhedralObject{T} - pm_subdivision::Polymake.BigObject - parent_field::Field + pm_subdivision::Polymake.BigObject + parent_field::Field - SubdivisionOfPoints{T}(pm::Polymake.BigObject, p::Field) where T<:scalar_types = new{T}(pm, p) - SubdivisionOfPoints{QQFieldElem}(pm::Polymake.BigObject) = new{QQFieldElem}(pm, QQ) + SubdivisionOfPoints{T}(pm::Polymake.BigObject, p::Field) where {T<:scalar_types} = + new{T}(pm, p) + SubdivisionOfPoints{QQFieldElem}(pm::Polymake.BigObject) = new{QQFieldElem}(pm, QQ) end - # default scalar type: guess from input, fallback to QQ subdivision_of_points(points::AbstractCollection[PointVector], cells) = subdivision_of_points(_guess_fieldelem_type(points), points, cells) @@ -22,8 +22,8 @@ subdivision_of_points(points::AbstractCollection[PointVector], weights::Abstract # Automatic detection of corresponding OSCAR scalar type; # Avoid, if possible, to increase type stability function subdivision_of_points(p::Polymake.BigObject) - T, f = _detect_scalar_and_field(SubdivisionOfPoints, p) - return SubdivisionOfPoints{T}(p, f) + T, f = _detect_scalar_and_field(SubdivisionOfPoints, p) + return SubdivisionOfPoints{T}(p, f) end @doc raw""" @@ -53,17 +53,20 @@ julia> MOAE = subdivision_of_points(moaepts, moaeimnonreg0) Subdivision of points in ambient dimension 3 ``` """ -function subdivision_of_points(f::scalar_type_or_field, points::AbstractCollection[PointVector], cells::IncidenceMatrix) +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 = homogenized_matrix(points,1), - MAXIMAL_CELLS = arr, - ), parent_field) + arr = Polymake.@convert_to Array{Set{Int}} Polymake.common.rows(cells) + SubdivisionOfPoints{scalar_type}( + Polymake.fan.SubdivisionOfPoints{_scalar_type_to_polymake(scalar_type)}(; + POINTS=homogenized_matrix(points, 1), MAXIMAL_CELLS=arr + ), + parent_field, + ) end - @doc raw""" subdivision_of_points([f = QQFieldElem,] points, weights) @@ -91,13 +94,17 @@ julia> n_maximal_cells(SOP) 1 ``` """ -function subdivision_of_points(f::scalar_type_or_field, points::AbstractCollection[PointVector], weights::AbstractVector) +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 = homogenized_matrix(points,1), - WEIGHTS = weights, - ), parent_field) + SubdivisionOfPoints{scalar_type}( + Polymake.fan.SubdivisionOfPoints{_scalar_type_to_polymake(scalar_type)}(; + POINTS=homogenized_matrix(points, 1), WEIGHTS=weights + ), + parent_field, + ) end """ @@ -107,14 +114,19 @@ Get the underlying polymake object, which can be used via Polymake.jl. """ pm_object(SOP::SubdivisionOfPoints) = SOP.pm_subdivision - #Same construction for when the user provides maximal cells -function subdivision_of_points(f::scalar_type_or_field, points::AbstractCollection[PointVector], cells::Vector{Vector{Int64}}) +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(f::scalar_type_or_field, Points::Union{Oscar.MatElem,AbstractMatrix}, cells::Matrix{Bool}) +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 @@ -123,7 +135,7 @@ end ### Display ############################################################################### ############################################################################### -function Base.show(io::IO, SOP::SubdivisionOfPoints{T}) where T<:scalar_types - print(io, "Subdivision of points in ambient dimension $(ambient_dim(SOP))") - T != QQFieldElem && print(io, " with $T type coefficients") +function Base.show(io::IO, SOP::SubdivisionOfPoints{T}) where {T<:scalar_types} + print(io, "Subdivision of points in ambient dimension $(ambient_dim(SOP))") + T != QQFieldElem && print(io, " with $T type coefficients") end diff --git a/src/PolyhedralGeometry/SubdivisionOfPoints/functions.jl b/src/PolyhedralGeometry/SubdivisionOfPoints/functions.jl index c7ad3843f0a7..50775eeb099f 100644 --- a/src/PolyhedralGeometry/SubdivisionOfPoints/functions.jl +++ b/src/PolyhedralGeometry/SubdivisionOfPoints/functions.jl @@ -28,8 +28,6 @@ julia> dim(C) 4 ``` """ -function secondary_cone(SOP::SubdivisionOfPoints{T}) where T<:scalar_types - Cone{T}(Polymake.fan.secondary_cone(pm_object(SOP)), coefficient_field(SOP)) +function secondary_cone(SOP::SubdivisionOfPoints{T}) where {T<:scalar_types} + Cone{T}(Polymake.fan.secondary_cone(pm_object(SOP)), coefficient_field(SOP)) end - - diff --git a/src/PolyhedralGeometry/SubdivisionOfPoints/properties.jl b/src/PolyhedralGeometry/SubdivisionOfPoints/properties.jl index a8090ad25998..133fac77face 100644 --- a/src/PolyhedralGeometry/SubdivisionOfPoints/properties.jl +++ b/src/PolyhedralGeometry/SubdivisionOfPoints/properties.jl @@ -4,7 +4,6 @@ ############################################################################### ############################################################################### - @doc raw""" points(SOP::SubdivisionOfPoints) @@ -29,20 +28,20 @@ julia> points(MOAE) [1, 1, 2] ``` """ -function points(SOP::SubdivisionOfPoints{T}) where T<:scalar_types +function points(SOP::SubdivisionOfPoints{T}) where {T<:scalar_types} return SubObjectIterator{PointVector{T}}(SOP, _point, size(pm_object(SOP).POINTS, 1)) 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( + 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] +_point_matrix(::Val{_point}, SOP::SubdivisionOfPoints; homogenized=false) = + pm_object(SOP).POINTS[:, (homogenized ? 1 : 2):end] _matrix_for_polymake(::Val{_point}) = _point_matrix - - - - @doc raw""" maximal_cells(SOP::SubdivisionOfPoints) @@ -86,11 +85,13 @@ julia> maximal_cells(MOAE) """ maximal_cells(SOP::SubdivisionOfPoints) = maximal_cells(Vector{Int}, SOP) function maximal_cells(::Type{Vector{Int}}, SOP::SubdivisionOfPoints) - return SubObjectIterator{Vector{Int}}(SOP, _maximal_cell, size(pm_object(SOP).MAXIMAL_CELLS, 1)) + return SubObjectIterator{Vector{Int}}( + SOP, _maximal_cell, size(pm_object(SOP).MAXIMAL_CELLS, 1) + ) end -_maximal_cell(::Type{Vector{Int}}, SOP::SubdivisionOfPoints, i::Base.Integer) = Vector{Int}(Polymake.row(pm_object(SOP).MAXIMAL_CELLS, i)) - +_maximal_cell(::Type{Vector{Int}}, SOP::SubdivisionOfPoints, i::Base.Integer) = + Vector{Int}(Polymake.row(pm_object(SOP).MAXIMAL_CELLS, i)) ############################################################################### ############################################################################### @@ -138,7 +139,6 @@ julia> ambient_dim(SOP) """ ambient_dim(SOP::SubdivisionOfPoints) = pm_object(SOP).VECTOR_AMBIENT_DIM::Int - 1 - @doc raw""" n_points(SOP::SubdivisionOfPoints) @@ -156,8 +156,6 @@ julia> n_points(SOP) """ n_points(SOP::SubdivisionOfPoints) = pm_object(SOP).N_POINTS::Int - - ############################################################################### ## Points properties ############################################################################### @@ -187,7 +185,6 @@ julia> min_weights(SOP) """ min_weights(SOP::SubdivisionOfPoints) = Vector{Int}(pm_object(SOP).MIN_WEIGHTS) - @doc raw""" maximal_cells(IncidenceMatrix, SOP::SubdivisionOfPoints) @@ -217,8 +214,8 @@ julia> maximal_cells(IncidenceMatrix, SOP) [1, 2, 3, 4, 5, 6] ``` """ -maximal_cells(::Type{IncidenceMatrix}, SOP::SubdivisionOfPoints) = pm_object(SOP).MAXIMAL_CELLS - +maximal_cells(::Type{IncidenceMatrix}, SOP::SubdivisionOfPoints) = + pm_object(SOP).MAXIMAL_CELLS ############################################################################### ## Boolean properties diff --git a/src/PolyhedralGeometry/groups.jl b/src/PolyhedralGeometry/groups.jl index be4c6597fdba..2725f0e061b3 100644 --- a/src/PolyhedralGeometry/groups.jl +++ b/src/PolyhedralGeometry/groups.jl @@ -1,51 +1,47 @@ function PermGroup_to_polymake_array(G::PermGroup) - generators = gens(G) - d = degree(G) - return _group_generators_to_pm_arr_arr(generators, d) + generators = gens(G) + d = degree(G) + return _group_generators_to_pm_arr_arr(generators, d) end - function _group_generators_to_pm_arr_arr(generators::AbstractVector{PermGroupElem}, d::Int) - if length(generators) == 0 - generators = elements(trivial_subgroup(symmetric_group(d))[1]) - end - result = Polymake.Array{Polymake.Array{Polymake.to_cxx_type(Int)}}(length(generators)) - i = 1 - for g in generators - array = Polymake.Array{Polymake.to_cxx_type(Int)}(d) - for j in 1:d - array[j] = g(j)-1 - end - result[i] = array - i = i+1 + if length(generators) == 0 + generators = elements(trivial_subgroup(symmetric_group(d))[1]) + end + result = Polymake.Array{Polymake.Array{Polymake.to_cxx_type(Int)}}(length(generators)) + i = 1 + for g in generators + array = Polymake.Array{Polymake.to_cxx_type(Int)}(d) + for j in 1:d + array[j] = g(j) - 1 end - return result + result[i] = array + i = i + 1 + end + return result end - function _gens_to_group(gens::Vector{PermGroupElem}) - @req length(gens) > 0 "List of generators empty, could not deduce degree" - S = parent(gens[1]) - return sub(S, gens)[1] + @req length(gens) > 0 "List of generators empty, could not deduce degree" + S = parent(gens[1]) + return sub(S, gens)[1] end -function _gens_to_group(gens::Dict{Symbol, Vector{PermGroupElem}}) - return Dict{Symbol, PermGroup}([k => _gens_to_group(v) for (k,v) in gens]) +function _gens_to_group(gens::Dict{Symbol,Vector{PermGroupElem}}) + return Dict{Symbol,PermGroup}([k => _gens_to_group(v) for (k, v) in gens]) end - function _pm_arr_arr_to_group_generators(M, n) - S=symmetric_group(n) - perm_bucket = Vector{Oscar.BasicGAPGroupElem{PermGroup}}() - for g in M - push!(perm_bucket, perm(S, Polymake.to_one_based_indexing(g))) - end - if length(M) == 0 - push!(perm_bucket, perm(S, Int[])) - end - return perm_bucket + S = symmetric_group(n) + perm_bucket = Vector{Oscar.BasicGAPGroupElem{PermGroup}}() + for g in M + push!(perm_bucket, perm(S, Polymake.to_one_based_indexing(g))) + end + if length(M) == 0 + push!(perm_bucket, perm(S, Int[])) + end + return perm_bucket end - @doc raw""" combinatorial_symmetries(P::Polyhedron) @@ -74,8 +70,8 @@ julia> order(G) 2 ``` """ -combinatorial_symmetries(P::Polyhedron) = automorphism_group(P; type=:combinatorial, action=:on_vertices) - +combinatorial_symmetries(P::Polyhedron) = + automorphism_group(P; type=:combinatorial, action=:on_vertices) @doc raw""" linear_symmetries(P::Polyhedron) @@ -113,7 +109,6 @@ julia> order(G) """ linear_symmetries(P::Polyhedron) = automorphism_group(P; type=:linear, action=:on_vertices) - @doc raw""" automorphism_group_generators(P::Polyhedron; type = :combinatorial, action = :all) @@ -183,32 +178,32 @@ Dict{Symbol, Vector{PermGroupElem}} with 2 entries: :on_facets => [(1,2)(3,4)] ``` """ -function automorphism_group_generators(P::Polyhedron; type = :combinatorial, action = :all) - @req is_bounded(P) "Automorphism groups not supported for unbounded polyhedra." - if type == :combinatorial - IM = vertex_indices(facets(P)) - if action == :all - result = automorphism_group_generators(IM; action = action) - return Dict{Symbol, Vector{PermGroupElem}}( - :on_vertices => result[:on_cols], - :on_facets => result[:on_rows]) - elseif action == :on_vertices - return automorphism_group_generators(IM; action = :on_cols) - elseif action == :on_facets - return automorphism_group_generators(IM; action = :on_rows) - else - throw(ArgumentError("Action $(action) not supported.")) - end - elseif type == :linear - return _linear_symmetries_generators(P; action = action) +function automorphism_group_generators(P::Polyhedron; type=:combinatorial, action=:all) + @req is_bounded(P) "Automorphism groups not supported for unbounded polyhedra." + if type == :combinatorial + IM = vertex_indices(facets(P)) + if action == :all + result = automorphism_group_generators(IM; action=action) + return Dict{Symbol,Vector{PermGroupElem}}( + :on_vertices => result[:on_cols], :on_facets => result[:on_rows] + ) + elseif action == :on_vertices + return automorphism_group_generators(IM; action=:on_cols) + elseif action == :on_facets + return automorphism_group_generators(IM; action=:on_rows) else - throw(ArgumentError("Action type $(type) not supported.")) + throw(ArgumentError("Action $(action) not supported.")) end - return Dict{Symbol, Vector{PermGroupElem}}(:on_vertices => vertex_action, - :on_facets => facet_action) + elseif type == :linear + return _linear_symmetries_generators(P; action=action) + else + throw(ArgumentError("Action type $(type) not supported.")) + end + return Dict{Symbol,Vector{PermGroupElem}}( + :on_vertices => vertex_action, :on_facets => facet_action + ) end - @doc raw""" automorphism_group_generators(IM::IncidenceMatrix; action = :all) @@ -256,51 +251,48 @@ julia> automorphism_group_generators(IM; action = :on_cols) (1,2)(3,4)(5,6)(7,8) ``` """ -function automorphism_group_generators(IM::IncidenceMatrix; action = :all) - gens = Polymake.graph.automorphisms(IM) - rows_action = _pm_arr_arr_to_group_generators([first(g) for g in gens], Polymake.nrows(IM)) - if action == :on_rows - return rows_action - end - cols_action = _pm_arr_arr_to_group_generators([last(g) for g in gens], Polymake.ncols(IM)) - if action == :on_cols - return cols_action - elseif action == :all - return Dict{Symbol, Vector{PermGroupElem}}( - :on_rows => rows_action, - :on_cols => cols_action - ) - else - throw(ArgumentError("Action $(action) not supported.")) - end +function automorphism_group_generators(IM::IncidenceMatrix; action=:all) + gens = Polymake.graph.automorphisms(IM) + rows_action = _pm_arr_arr_to_group_generators( + [first(g) for g in gens], Polymake.nrows(IM) + ) + if action == :on_rows + return rows_action + end + cols_action = _pm_arr_arr_to_group_generators([last(g) for g in gens], Polymake.ncols(IM)) + if action == :on_cols + return cols_action + elseif action == :all + return Dict{Symbol,Vector{PermGroupElem}}( + :on_rows => rows_action, :on_cols => cols_action + ) + else + throw(ArgumentError("Action $(action) not supported.")) + end end - - - -function _linear_symmetries_generators(P::Polyhedron; action = :all) - if is_bounded(P) - gp = Polymake.polytope.linear_symmetries(vertices(P)) - vaction = gp.PERMUTATION_ACTION - if action == :on_vertices - return _pm_arr_arr_to_group_generators(vaction.GENERATORS, vaction.DEGREE) - end - hgens = Polymake.group.induced_permutations(vaction.GENERATORS, vertex_indices(facets(P))) - if action == :on_facets - return _pm_arr_arr_to_group_generators(hgens, n_facets(P)) - end - return Dict{Symbol, Vector{PermGroupElem}}( - :on_vertices => _pm_arr_arr_to_group_generators(vaction.GENERATORS, vaction.DEGREE), - :on_facets => _pm_arr_arr_to_group_generators(hgens, n_facets(P)) - ) - else - throw(ArgumentError("Linear symmetries supported for bounded polyhedra only")) + +function _linear_symmetries_generators(P::Polyhedron; action=:all) + if is_bounded(P) + gp = Polymake.polytope.linear_symmetries(vertices(P)) + vaction = gp.PERMUTATION_ACTION + if action == :on_vertices + return _pm_arr_arr_to_group_generators(vaction.GENERATORS, vaction.DEGREE) + end + hgens = Polymake.group.induced_permutations( + vaction.GENERATORS, vertex_indices(facets(P)) + ) + if action == :on_facets + return _pm_arr_arr_to_group_generators(hgens, n_facets(P)) end + return Dict{Symbol,Vector{PermGroupElem}}( + :on_vertices => _pm_arr_arr_to_group_generators(vaction.GENERATORS, vaction.DEGREE), + :on_facets => _pm_arr_arr_to_group_generators(hgens, n_facets(P)), + ) + else + throw(ArgumentError("Linear symmetries supported for bounded polyhedra only")) + end end - - - - @doc raw""" automorphism_group(P::Polyhedron; type = :combinatorial, action = :all) @@ -309,12 +301,11 @@ values are the same as for [`automorphism_group_generators(P::Polyhedron; type = :combinatorial, action = :all)`](@ref) except that groups are returned instead of generators of groups. """ -function automorphism_group(P::Polyhedron; type = :combinatorial, action = :all) - result = automorphism_group_generators(P; type = type, action = action) - return _gens_to_group(result) +function automorphism_group(P::Polyhedron; type=:combinatorial, action=:all) + result = automorphism_group_generators(P; type=type, action=action) + return _gens_to_group(result) end - @doc raw""" automorphism_group(IM::IncidenceMatrix; action = :all) @@ -323,7 +314,7 @@ return values are the same as for [`automorphism_group_generators(IM::IncidenceMatrix; action = :all)`](@ref) except that groups are returned instead of generators of groups. """ -function automorphism_group(IM::IncidenceMatrix; action = :all) - result = automorphism_group_generators(IM; action = action) - return _gens_to_group(result) +function automorphism_group(IM::IncidenceMatrix; action=:all) + result = automorphism_group_generators(IM; action=action) + return _gens_to_group(result) end diff --git a/src/PolyhedralGeometry/helpers.jl b/src/PolyhedralGeometry/helpers.jl index 3b3e773c92b6..1733a098ce42 100644 --- a/src/PolyhedralGeometry/helpers.jl +++ b/src/PolyhedralGeometry/helpers.jl @@ -75,117 +75,142 @@ Set{Int64} with 1 element: """ column(i::IncidenceMatrix, n::Int) = convert(Set{Int}, Polymake.col(i, n)) -const _polymake_scalars = Union{Polymake.Integer, Polymake.Rational, Polymake.QuadraticExtension, Polymake.OscarNumber, Float64, Polymake.TropicalNumber} -const _polymake_compatible_scalars = Union{QQFieldElem, ZZRingElem, Base.Integer, Base.Rational, _polymake_scalars} - -function assure_matrix_polymake(m::Union{AbstractMatrix{Any}, AbstractMatrix{FieldElem}}) - a, b = size(m) - if a > 0 - i = findfirst(_cannot_convert_to_fmpq, m) - t = i === nothing ? QQFieldElem : typeof(m[i]) - if t <: _polymake_scalars - m = Polymake.Matrix{Polymake.convert_to_pm_type(t)}(m) - else - m = Polymake.Matrix{_scalar_type_to_polymake(t)}(m) - end +const _polymake_scalars = Union{ + Polymake.Integer, + Polymake.Rational, + Polymake.QuadraticExtension, + Polymake.OscarNumber, + Float64, + Polymake.TropicalNumber, +} +const _polymake_compatible_scalars = Union{ + QQFieldElem,ZZRingElem,Base.Integer,Base.Rational,_polymake_scalars +} + +function assure_matrix_polymake(m::Union{AbstractMatrix{Any},AbstractMatrix{FieldElem}}) + a, b = size(m) + if a > 0 + i = findfirst(_cannot_convert_to_fmpq, m) + t = i === nothing ? QQFieldElem : typeof(m[i]) + if t <: _polymake_scalars + m = Polymake.Matrix{Polymake.convert_to_pm_type(t)}(m) else - m = Polymake.Matrix{Polymake.Rational}(undef, a, b) + m = Polymake.Matrix{_scalar_type_to_polymake(t)}(m) end - return m + else + m = Polymake.Matrix{Polymake.Rational}(undef, a, b) + end + return m end -assure_matrix_polymake(m::AbstractMatrix{<:FieldElem}) = Polymake.Matrix{Polymake.OscarNumber}(m) +assure_matrix_polymake(m::AbstractMatrix{<:FieldElem}) = + Polymake.Matrix{Polymake.OscarNumber}(m) assure_matrix_polymake(m::MatElem) = Polymake.Matrix{_scalar_type_to_polymake(eltype(m))}(m) -assure_matrix_polymake(m::Union{Oscar.ZZMatrix, Oscar.QQMatrix, AbstractMatrix{<:_polymake_compatible_scalars}}) = m +assure_matrix_polymake( + m::Union{Oscar.ZZMatrix,Oscar.QQMatrix,AbstractMatrix{<:_polymake_compatible_scalars}} +) = m -assure_matrix_polymake(m::SubArray{T, 2, U, V, W}) where {T<:Union{_polymake_scalars}, U, V, W} = Polymake.Matrix{T}(m) +assure_matrix_polymake(m::SubArray{T,2,U,V,W}) where {T<:Union{_polymake_scalars},U,V,W} = + Polymake.Matrix{T}(m) -function assure_vector_polymake(v::Union{AbstractVector{Any}, AbstractVector{FieldElem}}) - i = findfirst(_cannot_convert_to_fmpq, v) - T = i === nothing ? QQFieldElem : typeof(v[i]) - return Polymake.Vector{_scalar_type_to_polymake(T)}(v) +function assure_vector_polymake(v::Union{AbstractVector{Any},AbstractVector{FieldElem}}) + i = findfirst(_cannot_convert_to_fmpq, v) + T = i === nothing ? QQFieldElem : typeof(v[i]) + return Polymake.Vector{_scalar_type_to_polymake(T)}(v) end -assure_vector_polymake(v::AbstractVector{<:FieldElem}) = Polymake.Vector{Polymake.OscarNumber}(v) +assure_vector_polymake(v::AbstractVector{<:FieldElem}) = + Polymake.Vector{Polymake.OscarNumber}(v) assure_vector_polymake(v::AbstractVector{<:_polymake_compatible_scalars}) = v -affine_matrix_for_polymake(x::Tuple{<:AnyVecOrMat, <:AbstractVector}) = augment(unhomogenized_matrix(x[1]), -Vector(assure_vector_polymake(x[2]))) -affine_matrix_for_polymake(x::Tuple{<:AnyVecOrMat, <:Any}) = homogenized_matrix(x[1], -x[2]) +affine_matrix_for_polymake(x::Tuple{<:AnyVecOrMat,<:AbstractVector}) = + augment(unhomogenized_matrix(x[1]), -Vector(assure_vector_polymake(x[2]))) +affine_matrix_for_polymake(x::Tuple{<:AnyVecOrMat,<:Any}) = homogenized_matrix(x[1], -x[2]) -_cannot_convert_to_fmpq(x::Any) = !hasmethod(convert, Tuple{Type{QQFieldElem}, typeof(x)}) +_cannot_convert_to_fmpq(x::Any) = !hasmethod(convert, Tuple{Type{QQFieldElem},typeof(x)}) -linear_matrix_for_polymake(x::Union{Oscar.ZZMatrix, Oscar.QQMatrix, AbstractMatrix}) = assure_matrix_polymake(x) +linear_matrix_for_polymake(x::Union{Oscar.ZZMatrix,Oscar.QQMatrix,AbstractMatrix}) = + assure_matrix_polymake(x) -linear_matrix_for_polymake(x::AbstractVector{<:AbstractVector}) = assure_matrix_polymake(stack(x...)) +linear_matrix_for_polymake(x::AbstractVector{<:AbstractVector}) = + assure_matrix_polymake(stack(x...)) -matrix_for_polymake(x::Union{Oscar.ZZMatrix, Oscar.QQMatrix, AbstractMatrix}) = assure_matrix_polymake(x) +matrix_for_polymake(x::Union{Oscar.ZZMatrix,Oscar.QQMatrix,AbstractMatrix}) = + assure_matrix_polymake(x) -number_of_rows(x::SubArray{T, 2, U, V, W}) where {T, U, V, W} = size(x, 1) +number_of_rows(x::SubArray{T,2,U,V,W}) where {T,U,V,W} = size(x, 1) -_isempty_halfspace(x::Pair{<:Union{Oscar.MatElem, AbstractMatrix}, Any}) = isempty(x[1]) +_isempty_halfspace(x::Pair{<:Union{Oscar.MatElem,AbstractMatrix},Any}) = isempty(x[1]) _isempty_halfspace(x) = isempty(x) -function Polymake.Matrix{T}(x::Union{MatElem,AbstractMatrix{<:FieldElem}}) where T<:_polymake_scalars - res = Polymake.Matrix{T}(size(x)...) - return res .= x +function Polymake.Matrix{T}( + x::Union{MatElem,AbstractMatrix{<:FieldElem}} +) where {T<:_polymake_scalars} + res = Polymake.Matrix{T}(size(x)...) + return res .= x end -Base.convert(::Type{Polymake.Matrix{T}}, x::MatElem) where T = Polymake.Matrix{T}(x) +Base.convert(::Type{Polymake.Matrix{T}}, x::MatElem) where {T} = Polymake.Matrix{T}(x) -Base.convert(::Type{Polymake.QuadraticExtension{Polymake.Rational}}, x::QQFieldElem) = Polymake.QuadraticExtension(convert(Polymake.Rational, x)) +Base.convert(::Type{Polymake.QuadraticExtension{Polymake.Rational}}, x::QQFieldElem) = + Polymake.QuadraticExtension(convert(Polymake.Rational, x)) -Base.convert(::Type{<:Polymake.Integer}, x::ZZRingElem) = GC.@preserve x return Polymake.new_integer_from_fmpz(x) +Base.convert(::Type{<:Polymake.Integer}, x::ZZRingElem) = + GC.@preserve x return Polymake.new_integer_from_fmpz(x) -Base.convert(::Type{<:Polymake.Rational}, x::QQFieldElem) = GC.@preserve x return Polymake.new_rational_from_fmpq(x) +Base.convert(::Type{<:Polymake.Rational}, x::QQFieldElem) = + GC.@preserve x return Polymake.new_rational_from_fmpq(x) -Base.convert(::Type{<:Polymake.Rational}, x::ZZRingElem) = GC.@preserve x return Polymake.new_rational_from_fmpz(x) +Base.convert(::Type{<:Polymake.Rational}, x::ZZRingElem) = + GC.@preserve x return Polymake.new_rational_from_fmpz(x) -Base.convert(::Type{<:Polymake.Integer}, x::QQFieldElem) = GC.@preserve x return Polymake.new_integer_from_fmpq(x) +Base.convert(::Type{<:Polymake.Integer}, x::QQFieldElem) = + GC.@preserve x return Polymake.new_integer_from_fmpq(x) function Base.convert(::Type{ZZRingElem}, x::Polymake.Integer) - res = ZZRingElem() - GC.@preserve x Polymake.new_fmpz_from_integer(x, pointer_from_objref(res)) - return res + res = ZZRingElem() + GC.@preserve x Polymake.new_fmpz_from_integer(x, pointer_from_objref(res)) + return res end function Base.convert(::Type{QQFieldElem}, x::Polymake.Rational) - res = QQFieldElem() - GC.@preserve x Polymake.new_fmpq_from_rational(x, pointer_from_objref(res)) - return res + res = QQFieldElem() + GC.@preserve x Polymake.new_fmpq_from_rational(x, pointer_from_objref(res)) + return res end function Base.convert(::Type{QQFieldElem}, x::Polymake.Integer) - res = QQFieldElem() - GC.@preserve x Polymake.new_fmpq_from_integer(x, pointer_from_objref(res)) - return res + res = QQFieldElem() + GC.@preserve x Polymake.new_fmpq_from_integer(x, pointer_from_objref(res)) + return res end function Base.convert(::Type{ZZRingElem}, x::Polymake.Rational) - res = ZZRingElem() - GC.@preserve x Polymake.new_fmpz_from_rational(x, pointer_from_objref(res)) - return res + res = ZZRingElem() + GC.@preserve x Polymake.new_fmpz_from_rational(x, pointer_from_objref(res)) + return res end Base.convert(::Type{Polymake.OscarNumber}, x::FieldElem) = Polymake.OscarNumber(x) -(::Type{T})(x::Polymake.OscarNumber) where T<:FieldElem = convert(T, Polymake.unwrap(x)) +(::Type{T})(x::Polymake.OscarNumber) where {T<:FieldElem} = convert(T, Polymake.unwrap(x)) (R::QQField)(x::Polymake.Rational) = convert(QQFieldElem, x) (Z::ZZRing)(x::Polymake.Rational) = convert(ZZRingElem, x) function (NF::Hecke.EmbeddedNumField)(x::Polymake.QuadraticExtension{Polymake.Rational}) - g = Polymake.generating_field_elements(x) - if g.r == 0 || g.b == 0 - return NF(convert(QQFieldElem, g.a)) - end - isq = Hecke.is_quadratic_type(number_field(NF)) - @req isq[1] "Can not construct non-trivial QuadraticExtension in non-quadratic number field." - @req isq[2] == base_field(number_field(NF))(g.r) "Source and target fields do not match." - a = NF(basis(number_field(NF))[2]) - return convert(QQFieldElem, g.a) + convert(QQFieldElem, g.b) * a + g = Polymake.generating_field_elements(x) + if g.r == 0 || g.b == 0 + return NF(convert(QQFieldElem, g.a)) + end + isq = Hecke.is_quadratic_type(number_field(NF)) + @req isq[1] "Can not construct non-trivial QuadraticExtension in non-quadratic number field." + @req isq[2] == base_field(number_field(NF))(g.r) "Source and target fields do not match." + a = NF(basis(number_field(NF))[2]) + return convert(QQFieldElem, g.a) + convert(QQFieldElem, g.b) * a end (F::Field)(x::Polymake.Rational) = F(QQ(x)) @@ -198,20 +223,29 @@ Polymake.convert_to_pm_type(::Type{ZZMatrix}) = Polymake.Matrix{Polymake.Integer Polymake.convert_to_pm_type(::Type{QQMatrix}) = Polymake.Matrix{Polymake.Rational} Polymake.convert_to_pm_type(::Type{ZZRingElem}) = Polymake.Integer Polymake.convert_to_pm_type(::Type{QQFieldElem}) = Polymake.Rational -Polymake.convert_to_pm_type(::Type{T}) where T<:FieldElem = Polymake.OscarNumber -Polymake.convert_to_pm_type(::Type{<:MatElem{T}}) where T = Polymake.Matrix{Polymake.convert_to_pm_type(T)} -Polymake.convert_to_pm_type(::Type{<:Graph{T}}) where T<:Union{Directed,Undirected} = Polymake.Graph{T} +Polymake.convert_to_pm_type(::Type{T}) where {T<:FieldElem} = Polymake.OscarNumber +Polymake.convert_to_pm_type(::Type{<:MatElem{T}}) where {T} = + Polymake.Matrix{Polymake.convert_to_pm_type(T)} +Polymake.convert_to_pm_type(::Type{<:Graph{T}}) where {T<:Union{Directed,Undirected}} = + Polymake.Graph{T} -Base.convert(::Type{<:Polymake.Graph{T}}, g::Graph{T}) where T<:Union{Directed,Undirected} = Oscar.pm_object(g) +Base.convert( + ::Type{<:Polymake.Graph{T}}, g::Graph{T} +) where {T<:Union{Directed,Undirected}} = Oscar.pm_object(g) function remove_zero_rows(A::AbstractMatrix) - A[findall(x->!iszero(x),collect(eachrow(A))),:] + A[findall(x -> !iszero(x), collect(eachrow(A))), :] end function remove_zero_rows(A::AbstractMatrix{Float64}) - A[findall(x->!isapprox(x,zero(x),atol=Polymake._get_global_epsilon()),collect(eachrow(A))),:] + A[ + findall( + x -> !isapprox(x, zero(x); atol=Polymake._get_global_epsilon()), collect(eachrow(A)) + ), + :, + ] end function remove_zero_rows(A::Oscar.MatElem) - remove_zero_rows(Matrix(A)) + remove_zero_rows(Matrix(A)) end # function remove_redundant_rows(A::Union{Oscar.MatElem,AbstractMatrix}) @@ -230,29 +264,31 @@ end # end function augment(vec::AbstractVector, val) - s = size(vec) - @req s[1] > 0 "cannot homogenize empty vector" - res = similar(vec, (s[1] + 1,)) - res[1] = val + zero(first(vec)) - res[2:end] = vec - return assure_vector_polymake(res) + s = size(vec) + @req s[1] > 0 "cannot homogenize empty vector" + res = similar(vec, (s[1] + 1,)) + res[1] = val + zero(first(vec)) + res[2:end] = vec + return assure_vector_polymake(res) end function augment(mat::AbstractMatrix, vec::AbstractVector) - s = size(mat) - res = similar(mat, promote_type(eltype(mat), eltype(vec)),(s[1], s[2] + 1)) - res[:, 1] = vec - res[:, 2:end] = mat - return assure_matrix_polymake(res) -end - -homogenize(vec::AbstractVector, val::Number = 0) = augment(vec, val) -homogenize(mat::AbstractMatrix, val::Number = 1) = augment(mat, fill(val, size(mat, 1))) -homogenize(mat::MatElem, val::Number = 1) = homogenize(Matrix(mat), val) -homogenize(nothing,val::Number)=nothing -homogenized_matrix(x::Union{AbstractVecOrMat,MatElem,Nothing}, val::Number) = homogenize(x, val) + s = size(mat) + res = similar(mat, promote_type(eltype(mat), eltype(vec)), (s[1], s[2] + 1)) + res[:, 1] = vec + res[:, 2:end] = mat + return assure_matrix_polymake(res) +end + +homogenize(vec::AbstractVector, val::Number=0) = augment(vec, val) +homogenize(mat::AbstractMatrix, val::Number=1) = augment(mat, fill(val, size(mat, 1))) +homogenize(mat::MatElem, val::Number=1) = homogenize(Matrix(mat), val) +homogenize(nothing, val::Number) = nothing +homogenized_matrix(x::Union{AbstractVecOrMat,MatElem,Nothing}, val::Number) = + homogenize(x, val) homogenized_matrix(x::AbstractVector, val::Number) = permutedims(homogenize(x, val)) -homogenized_matrix(x::AbstractVector{<:AbstractVector}, val::Number) = stack((homogenize(x[i], val) for i in 1:length(x))...) +homogenized_matrix(x::AbstractVector{<:AbstractVector}, val::Number) = + stack((homogenize(x[i], val) for i in 1:length(x))...) dehomogenize(vec::AbstractVector) = vec[2:end] dehomogenize(mat::AbstractMatrix) = mat[:, 2:end] @@ -260,7 +296,8 @@ dehomogenize(mat::AbstractMatrix) = mat[:, 2:end] unhomogenized_matrix(x::AbstractVector) = assure_matrix_polymake(stack(x)) unhomogenized_matrix(x::AbstractMatrix) = assure_matrix_polymake(x) unhomogenized_matrix(x::MatElem) = Matrix(assure_matrix_polymake(x)) -unhomogenized_matrix(x::AbstractVector{<:AbstractVector}) = unhomogenized_matrix(stack(x...)) +unhomogenized_matrix(x::AbstractVector{<:AbstractVector}) = + unhomogenized_matrix(stack(x...)) """ stack(A::AbstractVecOrMat, B::AbstractVecOrMat) @@ -299,9 +336,10 @@ julia> stack([1 2], []) stack(A::AbstractMatrix, ::Nothing) = A stack(::Nothing, B::AbstractMatrix) = B stack(A::AbstractMatrix, B::AbstractMatrix) = [A; B] -stack(A::AbstractMatrix, B::AbstractVector) = isempty(B) ? A : [A; permutedims(B)] +stack(A::AbstractMatrix, B::AbstractVector) = isempty(B) ? A : [A; permutedims(B)] stack(A::AbstractVector, B::AbstractMatrix) = isempty(A) ? B : [permutedims(A); B] -stack(A::AbstractVector, B::AbstractVector) = isempty(A) ? B : [permutedims(A); permutedims(B)] +stack(A::AbstractVector, B::AbstractVector) = + isempty(A) ? B : [permutedims(A); permutedims(B)] stack(A::AbstractVector, ::Nothing) = permutedims(A) stack(::Nothing, B::AbstractVector) = permutedims(B) stack(x, y, z...) = stack(stack(x, y), z...) @@ -331,13 +369,13 @@ _ambient_dim(x::MatElem) = ncols(x) Given a (homogeneous) polymake matrix split into vertices and rays and dehomogenize. """ function decompose_vdata(A::AbstractMatrix) - vertex_indices = findall(!iszero, view(A, :, 1)) - ray_indices = findall(iszero, view(A, :, 1)) - return (vertices = A[vertex_indices, 2:end], rays = A[ray_indices, 2:end]) + vertex_indices = findall(!iszero, view(A, :, 1)) + ray_indices = findall(iszero, view(A, :, 1)) + return (vertices=A[vertex_indices, 2:end], rays=A[ray_indices, 2:end]) end function decompose_hdata(A) - (A = -A[:, 2:end], b = A[:, 1]) + (A=-A[:, 2:end], b=A[:, 1]) end # TODO: different printing within oscar? if yes, implement the following method @@ -352,7 +390,7 @@ abstract type PolyhedralObject{T} end # several toric types need to inherit from abstract schemes and thus cannot # inherit from PolyhedralObject, but they still need to behave like polyhedral objects # for some operations. -const PolyhedralObjectUnion = Union{PolyhedralObject, NormalToricVarietyType} +const PolyhedralObjectUnion = Union{PolyhedralObject,NormalToricVarietyType} @doc raw""" coefficient_field(P::Union{Polyhedron{T}, Cone{T}, PolyhedralFan{T}, PolyhedralComplex{T}) where T<:scalar_types @@ -368,25 +406,27 @@ julia> coefficient_field(c) Rational field ``` """ -coefficient_field(x::PolyhedralObject) = x.parent_field +coefficient_field(x::PolyhedralObject) = x.parent_field coefficient_field(x::PolyhedralObject{QQFieldElem}) = QQ -_get_scalar_type(::PolyhedralObject{T}) where T = T +_get_scalar_type(::PolyhedralObject{T}) where {T} = T _get_scalar_type(::NormalToricVarietyType) = QQFieldElem ################################################################################ ######## Scalar types ################################################################################ -const scalar_types = Union{FieldElem, Float64} +const scalar_types = Union{FieldElem,Float64} -const scalar_type_to_oscar = Dict{String, Type}([("Rational", QQFieldElem), - ("QuadraticExtension", Hecke.EmbeddedNumFieldElem{AbsSimpleNumFieldElem}), - ("QuadraticExtension", Hecke.EmbeddedNumFieldElem{AbsSimpleNumFieldElem}), - ("Float", Float64)]) +const scalar_type_to_oscar = Dict{String,Type}([ + ("Rational", QQFieldElem), + ("QuadraticExtension", Hecke.EmbeddedNumFieldElem{AbsSimpleNumFieldElem}), + ("QuadraticExtension", Hecke.EmbeddedNumFieldElem{AbsSimpleNumFieldElem}), + ("Float", Float64), +]) -const scalar_types_extended = Union{scalar_types, ZZRingElem} +const scalar_types_extended = Union{scalar_types,ZZRingElem} -const scalar_type_or_field = Union{Type{<:scalar_types}, Field} +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 @@ -397,16 +437,16 @@ _scalar_type_to_polymake(::Type{Float64}) = Float64 #################################################################### function _embedded_quadratic_field(r::ZZRingElem) - if iszero(r) - R, = rationals_as_number_field() - return Hecke.embedded_field(R, real_embeddings(R)[]) - else - R, = quadratic_field(r) - return Hecke.embedded_field(R, real_embeddings(R)[2]) - end + if iszero(r) + R, = rationals_as_number_field() + return Hecke.embedded_field(R, real_embeddings(R)[]) + else + R, = quadratic_field(r) + return Hecke.embedded_field(R, real_embeddings(R)[2]) + end end -function _check_field_polyhedral(::Type{T}) where T +function _check_field_polyhedral(::Type{T}) where {T} @req !(T <: NumFieldElem) "Number fields must be embedded, e.g. via `embedded_number_field`." @req hasmethod(isless, (T, T)) "Field must be ordered and have `isless` method." end @@ -437,25 +477,25 @@ end _parent_or_coefficient_field(::Type{Float64}, x...) = AbstractAlgebra.Floats{Float64}() _parent_or_coefficient_field(::Type{ZZRingElem}, x...) = ZZ -_parent_or_coefficient_field(r::Base.RefValue{<:Union{FieldElem, ZZRingElem}}, x...) = parent(r.x) -_parent_or_coefficient_field(v::AbstractArray{T}) where T<:Union{FieldElem, ZZRingElem} = _parent_or_coefficient_field(T, v) +_parent_or_coefficient_field(r::Base.RefValue{<:Union{FieldElem,ZZRingElem}}, x...) = + parent(r.x) +_parent_or_coefficient_field(v::AbstractArray{T}) where {T<:Union{FieldElem,ZZRingElem}} = + _parent_or_coefficient_field(T, v) # QQ is done as a special case here to avoid ambiguities -function _parent_or_coefficient_field(::Type{T}, e::Any) where T<:FieldElem - hasmethod(parent, (typeof(e),)) && elem_type(parent(e)) <: T ? - parent(e) : - missing +function _parent_or_coefficient_field(::Type{T}, e::Any) where {T<:FieldElem} + hasmethod(parent, (typeof(e),)) && elem_type(parent(e)) <: T ? parent(e) : missing end -function _parent_or_coefficient_field(::Type{T}, c::MatElem) where T<:FieldElem +function _parent_or_coefficient_field(::Type{T}, c::MatElem) where {T<:FieldElem} elem_type(base_ring(c)) <: T ? base_ring(c) : missing end -function _parent_or_coefficient_field(::Type{T}, c::AbstractArray) where T<:FieldElem +function _parent_or_coefficient_field(::Type{T}, c::AbstractArray) where {T<:FieldElem} first([collect(skipmissing(_parent_or_coefficient_field.(Ref(T), c))); missing]) end -function _parent_or_coefficient_field(::Type{T}, x, y...) where T <: scalar_types +function _parent_or_coefficient_field(::Type{T}, x, y...) where {T<:scalar_types} for c in (x, y...) p = _parent_or_coefficient_field(T, c) if p !== missing && elem_type(p) <: T @@ -465,12 +505,12 @@ function _parent_or_coefficient_field(::Type{T}, x, y...) where T <: scalar_type missing end -function _determine_parent_and_scalar(f::Union{Field, ZZRing}, x...) +function _determine_parent_and_scalar(f::Union{Field,ZZRing}, x...) _check_field_polyhedral(elem_type(f)) return (f, elem_type(f)) end -function _determine_parent_and_scalar(::Type{T}, x...) where T <: scalar_types +function _determine_parent_and_scalar(::Type{T}, x...) where {T<:scalar_types} T == QQFieldElem && return (QQ, QQFieldElem) p = _parent_or_coefficient_field(T, x...) @req p !== missing "Scalars of type $T require specification of a parent field. Please pass the desired Field instead of the type or have a $T contained in your input data." @@ -478,39 +518,74 @@ function _determine_parent_and_scalar(::Type{T}, x...) where T <: scalar_types return (p, elem_type(p)) end -function _detect_default_field(::Type{Hecke.EmbeddedNumFieldElem{AbsSimpleNumFieldElem}}, p::Polymake.BigObject) - # 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 QuadraticExtension with root != 0 - # or first OscarNumber wrapping an embedded number field element - while !isnothing(i) - prop = getproperty(p, propnames[i]) - if eltype(prop) <: Polymake.QuadraticExtension - for el in prop - r = Polymake.generating_field_elements(el).r - iszero(r) || return _embedded_quadratic_field(ZZ(r))[1] - end - elseif eltype(prop) <: Polymake.OscarNumber - for el in prop - on = Polymake.unwrap(el) - if on isa Hecke.EmbeddedNumFieldElem{AbsSimpleNumFieldElem} - return parent(on) - end - end +function _detect_default_field( + ::Type{Hecke.EmbeddedNumFieldElem{AbsSimpleNumFieldElem}}, p::Polymake.BigObject +) + # 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 QuadraticExtension with root != 0 + # or first OscarNumber wrapping an embedded number field element + while !isnothing(i) + prop = getproperty(p, propnames[i]) + if eltype(prop) <: Polymake.QuadraticExtension + for el in prop + r = Polymake.generating_field_elements(el).r + iszero(r) || return _embedded_quadratic_field(ZZ(r))[1] + end + elseif eltype(prop) <: Polymake.OscarNumber + for el in prop + on = Polymake.unwrap(el) + if on isa Hecke.EmbeddedNumFieldElem{AbsSimpleNumFieldElem} + return parent(on) end - i = findnext(f, propnames, i + 1) + end end - return _embedded_quadratic_field(ZZ(0))[1] + i = findnext(f, propnames, i + 1) + end + return _embedded_quadratic_field(ZZ(0))[1] end _detect_default_field(::Type{QQFieldElem}, p::Polymake.BigObject) = QQ -_detect_default_field(::Type{Float64}, p::Polymake.BigObject) = AbstractAlgebra.Floats{Float64}() +_detect_default_field(::Type{Float64}, p::Polymake.BigObject) = + AbstractAlgebra.Floats{Float64}() -function _detect_default_field(::Type{T}, p::Polymake.BigObject) where T<:FieldElem +function _detect_default_field(::Type{T}, p::Polymake.BigObject) where {T<:FieldElem} # 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"]) + 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)) @@ -526,7 +601,23 @@ 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"]) + 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)) @@ -542,7 +633,9 @@ function _detect_wrapped_type_and_field(p::Polymake.BigObject) 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 +function _detect_scalar_and_field( + ::Type{U}, p::Polymake.BigObject +) where {U<:PolyhedralObject} T = detect_scalar_type(U, p) if isnothing(T) return _detect_wrapped_type_and_field(p) @@ -552,48 +645,53 @@ function _detect_scalar_and_field(::Type{U}, p::Polymake.BigObject) where U<:Pol end # promotion helpers -function _promote_scalar_field(f::Union{Field, ZZRing}...) +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 end function _promote_scalar_field(a::AbstractArray{<:FieldElem}) - isempty(a) && return (QQFieldElem, QQ) - return _promote_scalar_field(parent.(a)...) + isempty(a) && return (QQFieldElem, QQ) + return _promote_scalar_field(parent.(a)...) end -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)) +function _promoted_bigobject( + ::Type{T}, obj::PolyhedralObject{U} +) where {T<:scalar_types,U<:scalar_types} + if T == U + pm_object(obj) + else + Polymake.common.convert_to{_scalar_type_to_polymake(T)}(pm_object(obj)) + end end # oscarnumber helpers function Polymake._fieldelem_to_rational(e::EmbeddedNumFieldElem) - return Rational{BigInt}(QQ(e)) + return Rational{BigInt}(QQ(e)) end function Polymake._fieldelem_is_rational(e::EmbeddedNumFieldElem) - return is_rational(e) + return is_rational(e) end function Polymake._fieldelem_to_float(e::EmbeddedNumFieldElem) - return Float64(real(embedding(parent(e))(data(e), 32))) + return Float64(real(embedding(parent(e))(data(e), 32))) end function Polymake._fieldelem_to_float(e::QQBarFieldElem) - return Float64(ArbField(64)(e)) + return Float64(ArbField(64)(e)) end function Polymake._fieldelem_from_rational(::QQBarField, r::Rational{BigInt}) return QQBarFieldElem(QQFieldElem(r)) end - # convert a Polymake.BigObject's scalar from QuadraticExtension to OscarNumber (Polytope only) function _polyhedron_qe_to_on(x::Polymake.BigObject, f::Field) @@ -605,14 +703,16 @@ function _polyhedron_qe_to_on(x::Polymake.BigObject, f::Field) 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.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} + if hasmethod(length, (typeof(x),)) && + eltype(x) <: Polymake.QuadraticExtension{Polymake.Rational} return f.(x) else return x diff --git a/src/PolyhedralGeometry/iterators.jl b/src/PolyhedralGeometry/iterators.jl index 007a58e90c3c..4afaec777051 100644 --- a/src/PolyhedralGeometry/iterators.jl +++ b/src/PolyhedralGeometry/iterators.jl @@ -3,17 +3,15 @@ ################################################################################ for (T, _t) in ((:PointVector, :point_vector), (:RayVector, :ray_vector)) - @eval begin - struct $T{U} <: AbstractVector{U} p::MatElem{U} - $T{U}(p::MatElem{U}) where U<:scalar_types_extended = new{U}(p) + $T{U}(p::MatElem{U}) where {U<:scalar_types_extended} = new{U}(p) end Base.IndexStyle(::Type{<:$T}) = IndexLinear() - Base.getindex(po::$T{U}, i::Base.Integer) where U = po.p[1, i]::U + Base.getindex(po::$T{U}, i::Base.Integer) where {U} = po.p[1, i]::U function Base.setindex!(po::$T, val, i::Base.Integer) @boundscheck 1 <= length(po) <= i @@ -27,43 +25,47 @@ for (T, _t) in ((:PointVector, :point_vector), (:RayVector, :ray_vector)) coefficient_field(po::$T) = base_ring(po.p) - function $_t(p::Union{scalar_type_or_field, ZZRing}, v::AbstractVector) + 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 $_t(p::Union{scalar_type_or_field, ZZRing}, n::Base.Integer) + 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 - $_t(x::Union{AbstractVector, Base.Integer}) = $_t(QQ, x) + $_t(x::Union{AbstractVector,Base.Integer}) = $_t(QQ, x) - function Base.similar(X::$T, ::Type{S}, dims::Dims{1}) where S <: scalar_types_extended + function Base.similar(X::$T, ::Type{S}, dims::Dims{1}) where {S<:scalar_types_extended} return $_t(coefficient_field(X), dims...) end Base.BroadcastStyle(::Type{<:$T}) = Broadcast.ArrayStyle{$T}() - _parent_or_coefficient_field(::Type{TT}, po::$T{<:TT}) where TT <: FieldElem = coefficient_field(po) + _parent_or_coefficient_field(::Type{TT}, po::$T{<:TT}) where {TT<:FieldElem} = + coefficient_field(po) _find_elem_type(po::$T) = elem_type(coefficient_field(po)) - function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{$T}}, ::Type{ElType}) where ElType<:scalar_types_extended + function Base.similar( + bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{$T}}, ::Type{ElType} + ) where {ElType<:scalar_types_extended} e = bc.f(first.(bc.args)...) return $_t(parent(e), axes(bc)...) end - function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{$T}}, ::Type{ElType}) where ElType + function Base.similar( + bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{$T}}, ::Type{ElType} + ) where {ElType} return Vector{ElType}(undef, length(axes(bc)...)) end Base.:*(k::scalar_types_extended, po::$T) = k .* po Base.:*(A::MatElem, v::$T) = A * transpose(v.p) - end end @@ -88,7 +90,6 @@ ray_vector ################################################################################ for (h, comp) in (("halfspace", "≤"), ("hyperplane", "=")) - H = uppercasefirst(h) Habs = Symbol(H) Haff = Symbol("Affine", H) @@ -98,7 +99,6 @@ for (h, comp) in (("halfspace", "≤"), ("hyperplane", "=")) Flin = Symbol("linear_", h) @eval begin - abstract type $Habs{T} end # Affine types @@ -106,13 +106,14 @@ for (h, comp) in (("halfspace", "≤"), ("hyperplane", "=")) a::MatElem{T} b::T - $Haff{T}(a::MatElem{T}, b::T) where T<:scalar_types = new{T}(a, b) + $Haff{T}(a::MatElem{T}, b::T) where {T<:scalar_types} = new{T}(a, b) 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) + $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) - invert(H::$Haff{T}) where T<:scalar_types = $Haff{T}(-H.a, -negbias(H)) + invert(H::$Haff{T}) where {T<:scalar_types} = $Haff{T}(-H.a, -negbias(H)) @doc """ $($Faff)(p = QQ, a, b) @@ -121,27 +122,29 @@ Return the `$($Haff)` `H(a,b)`, which is given by a vector `a` and a value `b` s \$\$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) + 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 - $Faff(a::Union{MatElem, AbstractMatrix, AbstractVector}, b = 0) = $Faff(QQ, a, b) - + $Faff(a::Union{MatElem,AbstractMatrix,AbstractVector}, b=0) = $Faff(QQ, a, b) # Linear types struct $Hlin{T} <: $Habs{T} a::MatElem{T} - - $Hlin{T}(a::MatElem{T}) where T<:scalar_types = new{T}(a) + + $Hlin{T}(a::MatElem{T}) where {T<:scalar_types} = new{T}(a) end - $Fabs(a::Union{MatElem, AbstractMatrix, AbstractVector}) = $Flin(a) - $Fabs(f::scalar_type_or_field, a::Union{MatElem, AbstractMatrix, AbstractVector}) = $Flin(f, a) + $Fabs(a::Union{MatElem,AbstractMatrix,AbstractVector}) = $Flin(a) + $Fabs(f::scalar_type_or_field, a::Union{MatElem,AbstractMatrix,AbstractVector}) = + $Flin(f, a) - invert(H::$Hlin{T}) where T<:scalar_types = $Hlin{T}(-H.a) + invert(H::$Hlin{T}) where {T<:scalar_types} = $Hlin{T}(-H.a) @doc """ $($Flin)(p = QQ, a, b) @@ -150,29 +153,28 @@ 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}) + 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 - $Flin(a::Union{MatElem, AbstractMatrix, AbstractVector}) = $Flin(QQ, a) + $Flin(a::Union{MatElem,AbstractMatrix,AbstractVector}) = $Flin(QQ, a) coefficient_field(h::$Habs) = base_ring(h.a) _find_elem_type(h::$Habs) = elem_type(coefficient_field(h)) - _parent_or_coefficient_field(::Type{T}, h::$Habs{<:T}) where T <: FieldElem = coefficient_field(h) - + _parent_or_coefficient_field(::Type{T}, h::$Habs{<:T}) where {T<:FieldElem} = + coefficient_field(h) end - end # Field access -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)] +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) +_ambient_dim(x::Union{Halfspace,Hyperplane}) = length(x.a) function Base.:(==)(x::Halfspace, y::Halfspace) ax = normal_vector(x) @@ -224,20 +226,21 @@ Additional data required for specifying the property can be given using keyword arguments. """ struct SubObjectIterator{T} <: AbstractVector{T} - Obj::PolyhedralObjectUnion - Acc::Function - n::Int - options::NamedTuple + Obj::PolyhedralObjectUnion + Acc::Function + n::Int + options::NamedTuple end # `options` is empty by default -SubObjectIterator{T}(Obj::PolyhedralObjectUnion, Acc::Function, n::Base.Integer) where T = SubObjectIterator{T}(Obj, Acc, n, NamedTuple()) +SubObjectIterator{T}(Obj::PolyhedralObjectUnion, Acc::Function, n::Base.Integer) where {T} = + SubObjectIterator{T}(Obj, Acc, n, NamedTuple()) Base.IndexStyle(::Type{<:SubObjectIterator}) = IndexLinear() -function Base.getindex(iter::SubObjectIterator{T}, i::Base.Integer) where T - @boundscheck 1 <= i && i <= iter.n - return iter.Acc(T, iter.Obj, i; iter.options...)::T +function Base.getindex(iter::SubObjectIterator{T}, i::Base.Integer) where {T} + @boundscheck 1 <= i && i <= iter.n + return iter.Acc(T, iter.Obj, i; iter.options...)::T end Base.firstindex(::SubObjectIterator) = 1 @@ -247,101 +250,134 @@ Base.size(iter::SubObjectIterator) = (iter.n,) ################################################################################ # Incidence matrices -for (sym, name) in (("facet_indices", "Incidence matrix resp. facets"), ("ray_indices", "Incidence Matrix resp. rays"), ("vertex_indices", "Incidence Matrix resp. vertices"), ("vertex_and_ray_indices", "Incidence Matrix resp. vertices and rays")) - M = Symbol(sym) - _M = Symbol("_", sym) - @eval begin - $M(iter::SubObjectIterator) = $_M(Val(iter.Acc), iter.Obj; iter.options...) - $_M(::Any, ::PolyhedralObjectUnion) = throw(ArgumentError(string($name, " not defined in this context."))) - end +for (sym, name) in ( + ("facet_indices", "Incidence matrix resp. facets"), + ("ray_indices", "Incidence Matrix resp. rays"), + ("vertex_indices", "Incidence Matrix resp. vertices"), + ("vertex_and_ray_indices", "Incidence Matrix resp. vertices and rays"), +) + M = Symbol(sym) + _M = Symbol("_", sym) + @eval begin + $M(iter::SubObjectIterator) = $_M(Val(iter.Acc), iter.Obj; iter.options...) + $_M(::Any, ::PolyhedralObjectUnion) = + throw(ArgumentError(string($name, " not defined in this context."))) + end end # Matrices with rational or integer elements -for (sym, name) in (("point_matrix", "Point Matrix"), ("vector_matrix", "Vector Matrix"), ("generator_matrix", "Generator Matrix")) - M = Symbol(sym) - _M = Symbol("_", sym) - @eval begin - $M(iter::SubObjectIterator{<:AbstractVector{QQFieldElem}}) = matrix(QQ, $_M(Val(iter.Acc), iter.Obj; iter.options...)) - $M(iter::SubObjectIterator{<:AbstractVector{ZZRingElem}}) = matrix(ZZ, $_M(Val(iter.Acc), iter.Obj; iter.options...)) - $M(iter::SubObjectIterator{<:AbstractVector{<:FieldElem}}) = matrix(coefficient_field(iter.Obj), $_M(Val(iter.Acc), iter.Obj; iter.options...)) - $_M(::Any, ::PolyhedralObjectUnion) = throw(ArgumentError(string($name, " not defined in this context."))) - end +for (sym, name) in ( + ("point_matrix", "Point Matrix"), + ("vector_matrix", "Vector Matrix"), + ("generator_matrix", "Generator Matrix"), +) + M = Symbol(sym) + _M = Symbol("_", sym) + @eval begin + $M(iter::SubObjectIterator{<:AbstractVector{QQFieldElem}}) = + matrix(QQ, $_M(Val(iter.Acc), iter.Obj; iter.options...)) + $M(iter::SubObjectIterator{<:AbstractVector{ZZRingElem}}) = + matrix(ZZ, $_M(Val(iter.Acc), iter.Obj; iter.options...)) + $M(iter::SubObjectIterator{<:AbstractVector{<:FieldElem}}) = + matrix(coefficient_field(iter.Obj), $_M(Val(iter.Acc), iter.Obj; iter.options...)) + $_M(::Any, ::PolyhedralObjectUnion) = + throw(ArgumentError(string($name, " not defined in this context."))) + end end function matrix_for_polymake(iter::SubObjectIterator; homogenized=false) - if hasmethod(_matrix_for_polymake, Tuple{Val{iter.Acc}}) - return _matrix_for_polymake(Val(iter.Acc))(Val(iter.Acc), iter.Obj; homogenized=homogenized, iter.options...) - else - throw(ArgumentError("Matrix for Polymake not defined in this context.")) - end + if hasmethod(_matrix_for_polymake, Tuple{Val{iter.Acc}}) + return _matrix_for_polymake(Val(iter.Acc))( + Val(iter.Acc), iter.Obj; homogenized=homogenized, iter.options... + ) + else + throw(ArgumentError("Matrix for Polymake not defined in this context.")) + end end function IncidenceMatrix(iter::SubObjectIterator) - if hasmethod(_incidencematrix, Tuple{Val{iter.Acc}}) - return _incidencematrix(Val(iter.Acc))(Val(iter.Acc), iter.Obj; iter.options...) - else - throw(ArgumentError("IncidenceMatrix not defined in this context.")) - end + if hasmethod(_incidencematrix, Tuple{Val{iter.Acc}}) + return _incidencematrix(Val(iter.Acc))(Val(iter.Acc), iter.Obj; iter.options...) + else + throw(ArgumentError("IncidenceMatrix not defined in this context.")) + end end # primitive generators only for ray based iterators matrix(R::ZZRing, iter::SubObjectIterator{RayVector{QQFieldElem}}) = - matrix(R, Polymake.common.primitive(matrix_for_polymake(iter))) -matrix(R::ZZRing, iter::SubObjectIterator{<:Union{RayVector{ZZRingElem},PointVector{ZZRingElem}}}) = - matrix(R, matrix_for_polymake(iter)) -matrix(R::QQField, iter::SubObjectIterator{<:Union{RayVector{QQFieldElem}, PointVector{QQFieldElem}}}) = - matrix(R, matrix_for_polymake(iter)) -matrix(K, iter::SubObjectIterator{<:Union{RayVector{<:FieldElem}, PointVector{<:FieldElem}}}) = - matrix(K, matrix_for_polymake(iter)) - + matrix(R, Polymake.common.primitive(matrix_for_polymake(iter))) +matrix( + R::ZZRing, iter::SubObjectIterator{<:Union{RayVector{ZZRingElem},PointVector{ZZRingElem}}} +) = matrix(R, matrix_for_polymake(iter)) +matrix( + R::QQField, + iter::SubObjectIterator{<:Union{RayVector{QQFieldElem},PointVector{QQFieldElem}}}, +) = matrix(R, matrix_for_polymake(iter)) +matrix( + K, iter::SubObjectIterator{<:Union{RayVector{<:FieldElem},PointVector{<:FieldElem}}} +) = matrix(K, matrix_for_polymake(iter)) function linear_matrix_for_polymake(iter::SubObjectIterator) - if hasmethod(_linear_matrix_for_polymake, Tuple{Val{iter.Acc}}) - return _linear_matrix_for_polymake(Val(iter.Acc))(Val(iter.Acc), iter.Obj; iter.options...) - elseif hasmethod(_affine_matrix_for_polymake, Tuple{Val{iter.Acc}}) - res = _affine_matrix_for_polymake(Val(iter.Acc))(Val(iter.Acc), iter.Obj; iter.options...) - @req iszero(res[:, 1]) "Input not linear." - return res[:, 2:end] - end - throw(ArgumentError("Linear Matrix for Polymake not defined in this context.")) + if hasmethod(_linear_matrix_for_polymake, Tuple{Val{iter.Acc}}) + return _linear_matrix_for_polymake(Val(iter.Acc))( + Val(iter.Acc), iter.Obj; iter.options... + ) + elseif hasmethod(_affine_matrix_for_polymake, Tuple{Val{iter.Acc}}) + res = _affine_matrix_for_polymake(Val(iter.Acc))( + Val(iter.Acc), iter.Obj; iter.options... + ) + @req iszero(res[:, 1]) "Input not linear." + return res[:, 2:end] + end + throw(ArgumentError("Linear Matrix for Polymake not defined in this context.")) end function affine_matrix_for_polymake(iter::SubObjectIterator) - if hasmethod(_affine_matrix_for_polymake, Tuple{Val{iter.Acc}}) - return _affine_matrix_for_polymake(Val(iter.Acc))(Val(iter.Acc), iter.Obj; iter.options...) - elseif hasmethod(_linear_matrix_for_polymake, Tuple{Val{iter.Acc}}) - return homogenize(_linear_matrix_for_polymake(Val(iter.Acc))(Val(iter.Acc), iter.Obj; iter.options...), 0) - end - throw(ArgumentError("Affine Matrix for Polymake not defined in this context.")) + if hasmethod(_affine_matrix_for_polymake, Tuple{Val{iter.Acc}}) + return _affine_matrix_for_polymake(Val(iter.Acc))( + Val(iter.Acc), iter.Obj; iter.options... + ) + elseif hasmethod(_linear_matrix_for_polymake, Tuple{Val{iter.Acc}}) + return homogenize( + _linear_matrix_for_polymake(Val(iter.Acc))(Val(iter.Acc), iter.Obj; iter.options...), + 0, + ) + end + throw(ArgumentError("Affine Matrix for Polymake not defined in this context.")) end -Polymake.convert_to_pm_type(::Type{SubObjectIterator{RayVector{T}}}) where T = Polymake.Matrix{T} -Polymake.convert_to_pm_type(::Type{SubObjectIterator{PointVector{T}}}) where T = Polymake.Matrix{T} -Base.convert(::Type{<:Polymake.Matrix}, iter::SubObjectIterator) = assure_matrix_polymake(matrix_for_polymake(iter; homogenized=true)) +Polymake.convert_to_pm_type(::Type{SubObjectIterator{RayVector{T}}}) where {T} = + Polymake.Matrix{T} +Polymake.convert_to_pm_type(::Type{SubObjectIterator{PointVector{T}}}) where {T} = + Polymake.Matrix{T} +Base.convert(::Type{<:Polymake.Matrix}, iter::SubObjectIterator) = + assure_matrix_polymake(matrix_for_polymake(iter; homogenized=true)) -function homogenized_matrix(x::SubObjectIterator{<:PointVector}, v::Number = 1) - @req v == 1 "PointVectors can only be (re-)homogenized with parameter 1, please convert to a matrix first" - return matrix_for_polymake(x; homogenized=true) +function homogenized_matrix(x::SubObjectIterator{<:PointVector}, v::Number=1) + @req v == 1 "PointVectors can only be (re-)homogenized with parameter 1, please convert to a matrix first" + return matrix_for_polymake(x; homogenized=true) end -function homogenized_matrix(x::SubObjectIterator{<:RayVector}, v::Number = 0) - @req v == 0 "RayVectors can only be (re-)homogenized with parameter 0, please convert to a matrix first" - return matrix_for_polymake(x; homogenized=true) +function homogenized_matrix(x::SubObjectIterator{<:RayVector}, v::Number=0) + @req v == 0 "RayVectors can only be (re-)homogenized with parameter 0, please convert to a matrix first" + return matrix_for_polymake(x; homogenized=true) end -function homogenized_matrix(x::AbstractVector{<:PointVector}, v::Number = 1) - @req v == 1 "PointVectors can only be (re-)homogenized with parameter 1, please convert to a matrix first" - return stack((homogenize(x[i], v) for i in 1:length(x))...) +function homogenized_matrix(x::AbstractVector{<:PointVector}, v::Number=1) + @req v == 1 "PointVectors can only be (re-)homogenized with parameter 1, please convert to a matrix first" + return stack((homogenize(x[i], v) for i in 1:length(x))...) end -function homogenized_matrix(x::AbstractVector{<:RayVector}, v::Number = 0) - @req v == 0 "RayVectors can only be (re-)homogenized with parameter 0, please convert to a matrix first" - return stack((homogenize(x[i], v) for i in 1:length(x))...) +function homogenized_matrix(x::AbstractVector{<:RayVector}, v::Number=0) + @req v == 0 "RayVectors can only be (re-)homogenized with parameter 0, please convert to a matrix first" + return stack((homogenize(x[i], v) for i in 1:length(x))...) end -homogenized_matrix(::SubObjectIterator, v::Number) = throw(ArgumentError("Content of SubObjectIterator not suitable for homogenized_matrix.")) +homogenized_matrix(::SubObjectIterator, v::Number) = + throw(ArgumentError("Content of SubObjectIterator not suitable for homogenized_matrix.")) unhomogenized_matrix(x::SubObjectIterator{<:RayVector}) = matrix_for_polymake(x) -unhomogenized_matrix(x::AbstractVector{<:PointVector}) = throw(ArgumentError("unhomogenized_matrix only meaningful for RayVectors")) +unhomogenized_matrix(x::AbstractVector{<:PointVector}) = + throw(ArgumentError("unhomogenized_matrix only meaningful for RayVectors")) _ambient_dim(x::SubObjectIterator) = Polymake.polytope.ambient_dim(pm_object(x.Obj)) @@ -352,50 +388,57 @@ _ambient_dim(x::SubObjectIterator) = Polymake.polytope.ambient_dim(pm_object(x.O _empty_access() = nothing -function _empty_subobjectiterator(::Type{T}, Obj::PolyhedralObjectUnion) where T - return SubObjectIterator{T}(Obj, _empty_access, 0, NamedTuple()) +function _empty_subobjectiterator(::Type{T}, Obj::PolyhedralObjectUnion) where {T} + return SubObjectIterator{T}(Obj, _empty_access, 0, NamedTuple()) end for f in ("_point_matrix", "_vector_matrix", "_generator_matrix") - M = Symbol(f) - @eval begin - function $M(::Val{_empty_access}, P::PolyhedralObjectUnion; homogenized=false) - 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 + M = Symbol(f) + @eval begin + function $M(::Val{_empty_access}, P::PolyhedralObjectUnion; homogenized=false) + typename = Polymake.bigobject_eltype(pm_object(P)) + T = if typename == "OscarNumber" + Polymake.OscarNumber + else + _scalar_type_to_polymake(scalar_type_to_oscar[typename]) + end + return Polymake.Matrix{T}( + undef, 0, Polymake.polytope.ambient_dim(pm_object(P)) + homogenized + ) end + end end for f in ("_facet_indices", "_ray_indices", "_vertex_indices", "_vertex_and_ray_indices") - M = Symbol(f) - @eval begin - $M(::Val{_empty_access}, P::PolyhedralObjectUnion) = return Polymake.IncidenceMatrix(0, Polymake.polytope.ambient_dim(P)) - end + M = Symbol(f) + @eval begin + $M(::Val{_empty_access}, P::PolyhedralObjectUnion) = + return Polymake.IncidenceMatrix(0, Polymake.polytope.ambient_dim(P)) + end end for f in ("_linear_inequality_matrix", "_linear_equation_matrix") - M = Symbol(f) - @eval begin - function $M(::Val{_empty_access}, P::PolyhedralObjectUnion) - scalar_regexp = match(r"[^<]*<(.*)>[^>]*", String(Polymake.type_name(P))) - typename = scalar_regexp[1] - T = _scalar_type_to_polymake(scalar_type_to_oscar[typename]) - return Polymake.Matrix{T}(undef, 0, Polymake.polytope.ambient_dim(P)) - end + M = Symbol(f) + @eval begin + function $M(::Val{_empty_access}, P::PolyhedralObjectUnion) + scalar_regexp = match(r"[^<]*<(.*)>[^>]*", String(Polymake.type_name(P))) + typename = scalar_regexp[1] + T = _scalar_type_to_polymake(scalar_type_to_oscar[typename]) + return Polymake.Matrix{T}(undef, 0, Polymake.polytope.ambient_dim(P)) end + end end for f in ("_affine_inequality_matrix", "_affine_equation_matrix") - M = Symbol(f) - @eval begin - function $M(::Val{_empty_access}, P::PolyhedralObjectUnion) - scalar_regexp = match(r"[^<]*<(.*)>[^>]*", String(Polymake.type_name(P))) - typename = scalar_regexp[1] - T = _scalar_type_to_polymake(scalar_type_to_oscar[typename]) - return Polymake.Matrix{T}(undef, 0, Polymake.polytope.ambient_dim(P) + 1) - end + M = Symbol(f) + @eval begin + function $M(::Val{_empty_access}, P::PolyhedralObjectUnion) + scalar_regexp = match(r"[^<]*<(.*)>[^>]*", String(Polymake.type_name(P))) + typename = scalar_regexp[1] + T = _scalar_type_to_polymake(scalar_type_to_oscar[typename]) + return Polymake.Matrix{T}(undef, 0, Polymake.polytope.ambient_dim(P) + 1) end + end end _matrix_for_polymake(::Val{_empty_access}) = _point_matrix @@ -407,14 +450,38 @@ _matrix_for_polymake(::Val{_empty_access}) = _point_matrix # vector-like types: (un-)homogenized_matrix -> matrix_for_polymake # linear types: linear_matrix_for_polymake # affine types: affine_matrix_for_polymake -const AbstractCollection = Dict{UnionAll, Union}([(PointVector, AnyVecOrMat), - (RayVector, AnyVecOrMat), - (LinearHalfspace, Union{AbstractVector{<:Halfspace}, SubObjectIterator{<:Halfspace}, AnyVecOrMat}), - (LinearHyperplane, Union{AbstractVector{<:Hyperplane}, SubObjectIterator{<:Hyperplane}, AnyVecOrMat}), - (AffineHalfspace, Union{AbstractVector{<:Halfspace}, SubObjectIterator{<:Halfspace}, Tuple{<:AnyVecOrMat, <:Any}}), - (AffineHyperplane, Union{AbstractVector{<:Hyperplane}, SubObjectIterator{<:Hyperplane}, Tuple{<:AnyVecOrMat, <:Any}})]) - -affine_matrix_for_polymake(x::Union{Halfspace, Hyperplane}) = stack(augment(normal_vector(x), -negbias(x))) -affine_matrix_for_polymake(x::AbstractVector{<:Union{Halfspace, Hyperplane}}) = stack((affine_matrix_for_polymake(x[i]) for i in 1:length(x))...) -linear_matrix_for_polymake(x::Union{Halfspace, Hyperplane}) = negbias(x) == 0 ? stack(normal_vector(x)) : throw(ArgumentError("Input not linear.")) -linear_matrix_for_polymake(x::AbstractVector{<:Union{Halfspace, Hyperplane}}) = stack((linear_matrix_for_polymake(x[i]) for i in 1:length(x))...) +const AbstractCollection = Dict{UnionAll,Union}([ + (PointVector, AnyVecOrMat), + (RayVector, AnyVecOrMat), + ( + LinearHalfspace, + Union{AbstractVector{<:Halfspace},SubObjectIterator{<:Halfspace},AnyVecOrMat}, + ), + ( + LinearHyperplane, + Union{AbstractVector{<:Hyperplane},SubObjectIterator{<:Hyperplane},AnyVecOrMat}, + ), + ( + AffineHalfspace, + Union{ + AbstractVector{<:Halfspace},SubObjectIterator{<:Halfspace},Tuple{<:AnyVecOrMat,<:Any} + }, + ), + ( + AffineHyperplane, + Union{ + AbstractVector{<:Hyperplane}, + SubObjectIterator{<:Hyperplane}, + Tuple{<:AnyVecOrMat,<:Any}, + }, + ), +]) + +affine_matrix_for_polymake(x::Union{Halfspace,Hyperplane}) = + stack(augment(normal_vector(x), -negbias(x))) +affine_matrix_for_polymake(x::AbstractVector{<:Union{Halfspace,Hyperplane}}) = + stack((affine_matrix_for_polymake(x[i]) for i in 1:length(x))...) +linear_matrix_for_polymake(x::Union{Halfspace,Hyperplane}) = + negbias(x) == 0 ? stack(normal_vector(x)) : throw(ArgumentError("Input not linear.")) +linear_matrix_for_polymake(x::AbstractVector{<:Union{Halfspace,Hyperplane}}) = + stack((linear_matrix_for_polymake(x[i]) for i in 1:length(x))...) diff --git a/src/PolyhedralGeometry/linear_program.jl b/src/PolyhedralGeometry/linear_program.jl index 0831dd5b4d3f..33f8f7265f70 100644 --- a/src/PolyhedralGeometry/linear_program.jl +++ b/src/PolyhedralGeometry/linear_program.jl @@ -50,8 +50,7 @@ linear_program( c::AbstractVector; k=0, convention=:max, -) = - linear_program(polyhedron(f, A, b), c; k=k, convention=convention) +) = linear_program(polyhedron(f, A, b), c; k=k, convention=convention) pm_object(lp::LinearProgram) = lp.polymake_lp @@ -157,7 +156,9 @@ function optimal_vertex(lp::LinearProgram{T}) where {T<:scalar_types} opt_vert = lp.polymake_lp.MINIMAL_VERTEX end if !isnothing(opt_vert) - return point_vector(coefficient_field(lp), view(dehomogenize(opt_vert), :))::PointVector{T} + 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 a4d0fcaa2287..b3038e88484c 100644 --- a/src/PolyhedralGeometry/mixed_integer_linear_program.jl +++ b/src/PolyhedralGeometry/mixed_integer_linear_program.jl @@ -57,11 +57,15 @@ function mixed_integer_linear_program( integer_variables=Vector{Int64}([]), k=0, convention=:max, - ) where {T<:scalar_types} +) where {T<:scalar_types} P = polyhedron(T, A, b) return mixed_integer_linear_program( - P, c, coefficient_field(P); - integer_variables=integer_variables, k=k, convention=convention + P, + c, + coefficient_field(P); + integer_variables=integer_variables, + k=k, + convention=convention, ) end diff --git a/src/PolyhedralGeometry/solving_integrally.jl b/src/PolyhedralGeometry/solving_integrally.jl index 50230fdf74a6..7818e280a53e 100644 --- a/src/PolyhedralGeometry/solving_integrally.jl +++ b/src/PolyhedralGeometry/solving_integrally.jl @@ -1,24 +1,37 @@ -function solve_mixed(as::Type{SubObjectIterator{PointVector{ZZRingElem}}}, A::ZZMatrix, b::ZZMatrix, C::ZZMatrix, d::ZZMatrix; permit_unbounded=false) - @req ncols(A) == ncols(C) "solve_mixed(A,b,C,d): A and C must have the same number of columns." - @req nrows(A) == nrows(b) "solve_mixed(A,b,C,d): A and b must have the same number of rows." - @req nrows(C) == nrows(d) "solve_mixed(A,b,C,d): C and d must have the same number of rows." - @req ncols(b) == 1 "solve_mixed(A,b,C,d): b must be a matrix with a single column." - @req ncols(d) == 1 "solve_mixed(A,b,C,d): d must be a matrix with a single column." - P = polyhedron((-C, _vec(-d)), (A, _vec(b))) - if !permit_unbounded - return lattice_points(P) - else - sol = pm_object(P).LATTICE_POINTS_GENERATORS - return sol[1][:, 2:end] - end +function solve_mixed( + as::Type{SubObjectIterator{PointVector{ZZRingElem}}}, + A::ZZMatrix, + b::ZZMatrix, + C::ZZMatrix, + d::ZZMatrix; + permit_unbounded=false, +) + @req ncols(A) == ncols(C) "solve_mixed(A,b,C,d): A and C must have the same number of columns." + @req nrows(A) == nrows(b) "solve_mixed(A,b,C,d): A and b must have the same number of rows." + @req nrows(C) == nrows(d) "solve_mixed(A,b,C,d): C and d must have the same number of rows." + @req ncols(b) == 1 "solve_mixed(A,b,C,d): b must be a matrix with a single column." + @req ncols(d) == 1 "solve_mixed(A,b,C,d): d must be a matrix with a single column." + P = polyhedron((-C, _vec(-d)), (A, _vec(b))) + if !permit_unbounded + return lattice_points(P) + else + sol = pm_object(P).LATTICE_POINTS_GENERATORS + return sol[1][:, 2:end] + end end -function solve_mixed(as::Type{ZZMatrix}, A::ZZMatrix, b::ZZMatrix, C::ZZMatrix, d::ZZMatrix; permit_unbounded=false) - LP = solve_mixed(SubObjectIterator{PointVector{ZZRingElem}}, A, b, C, d; permit_unbounded) - return matrix(ZZ, LP) +function solve_mixed( + as::Type{ZZMatrix}, + A::ZZMatrix, + b::ZZMatrix, + C::ZZMatrix, + d::ZZMatrix; + permit_unbounded=false, +) + LP = solve_mixed(SubObjectIterator{PointVector{ZZRingElem}}, A, b, C, d; permit_unbounded) + return matrix(ZZ, LP) end - @doc raw""" solve_mixed(as::Type{T}, A::ZZMatrix, b::ZZMatrix, C::ZZMatrix, d::ZZMatrix) where {T} @@ -65,9 +78,11 @@ julia> for x in it [7] [7] [7] [7] [7] [7] [7] [7] ``` """ -solve_mixed(as::Type{T}, A::ZZMatrix, b::ZZMatrix, C::ZZMatrix, d::ZZMatrix; permit_unbounded=false) where {T} = solve_mixed(T, A, b, C, d; permit_unbounded) -solve_mixed(A::ZZMatrix, b::ZZMatrix, C::ZZMatrix, d::ZZMatrix; permit_unbounded=false) = solve_mixed(ZZMatrix, A, b, C, d; permit_unbounded) - +solve_mixed( + as::Type{T}, A::ZZMatrix, b::ZZMatrix, C::ZZMatrix, d::ZZMatrix; permit_unbounded=false +) where {T} = solve_mixed(T, A, b, C, d; permit_unbounded) +solve_mixed(A::ZZMatrix, b::ZZMatrix, C::ZZMatrix, d::ZZMatrix; permit_unbounded=false) = + solve_mixed(ZZMatrix, A, b, C, d; permit_unbounded) @doc raw""" solve_mixed(as::Type{T}, A::ZZMatrix, b::ZZMatrix, C::ZZMatrix) where {T} @@ -114,10 +129,11 @@ julia> for x in it [3] [3] [3] [3] ``` """ -solve_mixed(as::Type{T}, A::ZZMatrix, b::ZZMatrix, C::ZZMatrix; permit_unbounded=false) where {T} = solve_mixed(T, A, b, C, zero_matrix(FlintZZ, nrows(C), 1); permit_unbounded) -solve_mixed(A::ZZMatrix, b::ZZMatrix, C::ZZMatrix; permit_unbounded=false) = solve_mixed(ZZMatrix, A, b, C, zero_matrix(FlintZZ, nrows(C), 1); permit_unbounded) - - +solve_mixed( + as::Type{T}, A::ZZMatrix, b::ZZMatrix, C::ZZMatrix; permit_unbounded=false +) where {T} = solve_mixed(T, A, b, C, zero_matrix(FlintZZ, nrows(C), 1); permit_unbounded) +solve_mixed(A::ZZMatrix, b::ZZMatrix, C::ZZMatrix; permit_unbounded=false) = + solve_mixed(ZZMatrix, A, b, C, zero_matrix(FlintZZ, nrows(C), 1); permit_unbounded) @doc raw""" solve_ineq(as::Type{T}, A::ZZMatrix, b::ZZMatrix) where {T} @@ -154,10 +170,17 @@ julia> typeof(solve_ineq(SubObjectIterator{PointVector{ZZRingElem}}, A,b)) SubObjectIterator{PointVector{ZZRingElem}} ``` """ -solve_ineq(as::Type{T}, A::ZZMatrix, b::ZZMatrix; permit_unbounded=false) where {T} = solve_mixed(T, zero_matrix(FlintZZ, 0, ncols(A)), zero_matrix(FlintZZ,0,1), -A, -b; permit_unbounded) -solve_ineq(A::ZZMatrix, b::ZZMatrix; permit_unbounded=false) = solve_ineq(ZZMatrix, A, b; permit_unbounded) - - +solve_ineq(as::Type{T}, A::ZZMatrix, b::ZZMatrix; permit_unbounded=false) where {T} = + solve_mixed( + T, + zero_matrix(FlintZZ, 0, ncols(A)), + zero_matrix(FlintZZ, 0, 1), + -A, + -b; + permit_unbounded, + ) +solve_ineq(A::ZZMatrix, b::ZZMatrix; permit_unbounded=false) = + solve_ineq(ZZMatrix, A, b; permit_unbounded) @doc raw""" solve_non_negative(as::Type{T}, A::ZZMatrix, b::ZZMatrix) where {T} @@ -194,7 +217,8 @@ julia> typeof(solve_non_negative(SubObjectIterator{PointVector{ZZRingElem}}, A,b SubObjectIterator{PointVector{ZZRingElem}} ``` """ -solve_non_negative(as::Type{T}, A::ZZMatrix, b::ZZMatrix; permit_unbounded=false) where {T} = solve_mixed(T, A, b, identity_matrix(FlintZZ, ncols(A)); permit_unbounded) -solve_non_negative(A::ZZMatrix, b::ZZMatrix; permit_unbounded=false) = solve_non_negative(ZZMatrix, A, b; permit_unbounded) - - +solve_non_negative( + as::Type{T}, A::ZZMatrix, b::ZZMatrix; permit_unbounded=false +) where {T} = solve_mixed(T, A, b, identity_matrix(FlintZZ, ncols(A)); permit_unbounded) +solve_non_negative(A::ZZMatrix, b::ZZMatrix; permit_unbounded=false) = + solve_non_negative(ZZMatrix, A, b; permit_unbounded) diff --git a/src/PolyhedralGeometry/triangulations.jl b/src/PolyhedralGeometry/triangulations.jl index ad633b1b60b9..774f29a27e00 100644 --- a/src/PolyhedralGeometry/triangulations.jl +++ b/src/PolyhedralGeometry/triangulations.jl @@ -1,7 +1,9 @@ using JSON using TOPCOM_jll -function _postprocess_polymake_triangs(triangs::Polymake.Array{Polymake.Set{Polymake.Set{Polymake.to_cxx_type(Int)}}}) +function _postprocess_polymake_triangs( + triangs::Polymake.Array{Polymake.Set{Polymake.Set{Polymake.to_cxx_type(Int)}}} +) result = Vector{Vector{Int}}[] for triang in triangs push!(result, [Polymake.to_one_based_indexing(Vector{Int}(t)) for t in triang]) @@ -9,72 +11,76 @@ function _postprocess_polymake_triangs(triangs::Polymake.Array{Polymake.Set{Poly return result end -function topcom_regular_triangulations(pts::AbstractCollection[PointVector]; full::Bool=false) +function topcom_regular_triangulations( + pts::AbstractCollection[PointVector]; full::Bool=false +) input = homogenized_matrix(pts, 1) - PC = Polymake.polytope.PointConfiguration(POINTS=input) - result = full ? Polymake.polytope.topcom_fine_and_regular_triangulations(PC) : Polymake.polytope.topcom_regular_triangulations(PC) + PC = Polymake.polytope.PointConfiguration(; POINTS=input) + result = if full + Polymake.polytope.topcom_fine_and_regular_triangulations(PC) + else + Polymake.polytope.topcom_regular_triangulations(PC) + end return _postprocess_polymake_triangs(result) end -function topcom_regular_triangulation(pts::AbstractCollection[PointVector]; full::Bool=false) - input = homogenized_matrix(pts, 1) - inputstr = join(["["*join(input[i,:], ",")*"]" for i in 1:nrows(input)],",\n") - in = Pipe() - out = Pipe() - err = Pipe() - Base.link_pipe!(in, writer_supports_async=true) - Base.link_pipe!(out, reader_supports_async=true) - Base.link_pipe!(err, reader_supports_async=true) - cmd = Oscar.TOPCOM_jll.points2placingtriang() - if full - cmd = Oscar.TOPCOM_jll.points2finetriang() - end - proc = run(pipeline(`$(cmd) --regular`, stdin=in, stdout=out, stderr=err), wait=false) - task = @async begin - write(in, "[\n$inputstr\n]\n") - close(in) - end - close(in.out) - close(out.in) - close(err.in) - result = Vector{Vector{Vector{Int}}}() - for line in eachline(out) - m = match(r"{{.*}}", line) - triang = replace(m.match, "{"=>"[") - triang = replace(triang, "}"=>"]") - triang = convert(Vector{Vector{Int}},JSON.parse(triang)) - push!(result, Polymake.to_one_based_indexing(triang)) - end - wait(task) - if !success(proc) - msg = eof(err) ? "unknown error" : readchomp(err) - error("Failed to run TOPCOM: $msg") - end - return result +function topcom_regular_triangulation( + pts::AbstractCollection[PointVector]; full::Bool=false +) + input = homogenized_matrix(pts, 1) + inputstr = join(["[" * join(input[i, :], ",") * "]" for i in 1:nrows(input)], ",\n") + in = Pipe() + out = Pipe() + err = Pipe() + Base.link_pipe!(in; writer_supports_async=true) + Base.link_pipe!(out; reader_supports_async=true) + Base.link_pipe!(err; reader_supports_async=true) + cmd = Oscar.TOPCOM_jll.points2placingtriang() + if full + cmd = Oscar.TOPCOM_jll.points2finetriang() + end + proc = run(pipeline(`$(cmd) --regular`; stdin=in, stdout=out, stderr=err); wait=false) + task = @async begin + write(in, "[\n$inputstr\n]\n") + close(in) + end + close(in.out) + close(out.in) + close(err.in) + result = Vector{Vector{Vector{Int}}}() + for line in eachline(out) + m = match(r"{{.*}}", line) + triang = replace(m.match, "{" => "[") + triang = replace(triang, "}" => "]") + triang = convert(Vector{Vector{Int}}, JSON.parse(triang)) + push!(result, Polymake.to_one_based_indexing(triang)) + end + wait(task) + if !success(proc) + msg = eof(err) ? "unknown error" : readchomp(err) + error("Failed to run TOPCOM: $msg") + end + return result end - ################################################################################ function _is_triangulation(sop::SubdivisionOfPoints{QQFieldElem}) ad = ambient_dim(sop) for mc in maximal_cells(sop) - ad == length(mc)-1 || return false + ad == length(mc) - 1 || return false end return true end - - function _is_full_triangulation(sop::SubdivisionOfPoints{QQFieldElem}) _is_triangulation(sop) || return false length(Base.union(maximal_cells(sop)...)) == n_points(sop) || return false return true end -_is_star_triangulation(sop::SubdivisionOfPoints{QQFieldElem}) = all(cell -> 1 in cell, maximal_cells(sop)) - - +_is_star_triangulation(sop::SubdivisionOfPoints{QQFieldElem}) = + all(cell -> 1 in cell, maximal_cells(sop)) function _find_full_star_triangulation(pts::ZZMatrix; seed::Int=-1) seed == -1 || Random.seed!(seed) @@ -91,16 +97,14 @@ function _find_full_star_triangulation(pts::ZZMatrix; seed::Int=-1) return collect(maximal_cells(sop)) end - - function _is_star_triangulation(triang::Vector{Vector{Int}}) - u = Set{Int}() - for v in triang - if !(1 in v) - return false - end + u = Set{Int}() + for v in triang + if !(1 in v) + return false end - return true + end + return true end @doc raw""" @@ -135,14 +139,17 @@ julia> all_triangulations(V) ``` """ function all_triangulations(pts::AbstractCollection[PointVector]; full::Bool=false) - input = homogenized_matrix(pts, 1) - PC = Polymake.polytope.PointConfiguration(POINTS=input) - PC.FULL_DIM::Bool || error("Input points must have full rank.") - result = full ? Polymake.polytope.topcom_fine_triangulations(PC) : Polymake.polytope.topcom_all_triangulations(PC) - return _postprocess_polymake_triangs(result) + input = homogenized_matrix(pts, 1) + PC = Polymake.polytope.PointConfiguration(; POINTS=input) + PC.FULL_DIM::Bool || error("Input points must have full rank.") + result = if full + Polymake.polytope.topcom_fine_triangulations(PC) + else + Polymake.polytope.topcom_all_triangulations(PC) + end + return _postprocess_polymake_triangs(result) end - @doc raw""" all_triangulations(P::Polyhedron) @@ -167,12 +174,11 @@ julia> all_triangulations(c) ``` """ function all_triangulations(P::Polyhedron) - is_fulldimensional(P) || error("Input polytope must be full-dimensional.") - is_bounded(P) || error("Input polytope must be bounded.") - return all_triangulations(vertices(P); full=false) + is_fulldimensional(P) || error("Input polytope must be full-dimensional.") + is_bounded(P) || error("Input polytope must be bounded.") + return all_triangulations(vertices(P); full=false) end - @doc raw""" star_triangulations(pts::AbstractCollection[PointVector]; full::Bool=false, regular::Bool=false) @@ -186,13 +192,15 @@ a simplex as the set of indices of the vertices of the simplex. I.e. the `Vector{Int}` `[1,2,4]` corresponds to the simplex that is the convex hull of the first, second, and fourth input point. """ -function star_triangulations(pts::AbstractCollection[PointVector]; full::Bool=false, regular::Bool=false) - if regular - result = regular_triangulations(pts; full=full) - else - result = all_triangulations(pts; full=full) - end - return [t for t in result if _is_star_triangulation(t)] +function star_triangulations( + pts::AbstractCollection[PointVector]; full::Bool=false, regular::Bool=false +) + if regular + result = regular_triangulations(pts; full=full) + else + result = all_triangulations(pts; full=full) + end + return [t for t in result if _is_star_triangulation(t)] end @doc raw""" @@ -238,19 +246,16 @@ julia> star_triangulations(P) ``` """ function star_triangulations(P::Polyhedron; full::Bool=false, regular::Bool=false) - is_fulldimensional(P) || error("Input polytope must be full-dimensional.") - is_bounded(P) || error("Input polytope must be bounded.") - zero = [0 for i in 1:ambient_dim(P)] - @req zero in P "Input polyhedron must contain origin." - V = vertices(P) - V = [Vector{QQFieldElem}(v) for v in V if !iszero(v)] - pts = vcat(matrix(QQ, transpose(zero)), matrix(QQ, transpose(hcat(V...)))) - return pts, star_triangulations(pts; full=full, regular=regular) + is_fulldimensional(P) || error("Input polytope must be full-dimensional.") + is_bounded(P) || error("Input polytope must be bounded.") + zero = [0 for i in 1:ambient_dim(P)] + @req zero in P "Input polyhedron must contain origin." + V = vertices(P) + V = [Vector{QQFieldElem}(v) for v in V if !iszero(v)] + pts = vcat(matrix(QQ, transpose(zero)), matrix(QQ, transpose(hcat(V...)))) + return pts, star_triangulations(pts; full=full, regular=regular) end - - - @doc raw""" regular_triangulations(pts::AbstractCollection[PointVector]; full=false) @@ -287,13 +292,12 @@ julia> regular_triangulations(V) ``` """ function regular_triangulations(pts::AbstractCollection[PointVector]; full::Bool=false) - input = homogenized_matrix(pts, 1) - PC = Polymake.polytope.PointConfiguration(POINTS=input) - PC.FULL_DIM::Bool || error("Input points must have full rank.") - return topcom_regular_triangulations(pts; full=full) + input = homogenized_matrix(pts, 1) + PC = Polymake.polytope.PointConfiguration(; POINTS=input) + PC.FULL_DIM::Bool || error("Input points must have full rank.") + return topcom_regular_triangulations(pts; full=full) end - @doc raw""" regular_triangulations(P::Polyhedron) @@ -323,15 +327,11 @@ julia> regular_triangulations(c) ``` """ function regular_triangulations(P::Polyhedron) - is_fulldimensional(P) || error("Input polytope must be full-dimensional.") - is_bounded(P) || error("Input polytope must be bounded.") - return regular_triangulations(vertices(P); full=false) + is_fulldimensional(P) || error("Input polytope must be full-dimensional.") + is_bounded(P) || error("Input polytope must be bounded.") + return regular_triangulations(vertices(P); full=false) end - - - - @doc raw""" regular_triangulation(pts::AbstractCollection[PointVector]; full=false) @@ -369,13 +369,12 @@ julia> regular_triangulation(V) ``` """ function regular_triangulation(pts::AbstractCollection[PointVector]; full::Bool=false) - input = homogenized_matrix(pts, 1) - PC = Polymake.polytope.PointConfiguration(POINTS=input) - PC.FULL_DIM::Bool || error("Input points must have full rank.") - return topcom_regular_triangulation(pts; full=full) + input = homogenized_matrix(pts, 1) + PC = Polymake.polytope.PointConfiguration(; POINTS=input) + PC.FULL_DIM::Bool || error("Input points must have full rank.") + return topcom_regular_triangulation(pts; full=full) end - @doc raw""" regular_triangulation(P::Polyhedron) @@ -406,14 +405,11 @@ julia> regular_triangulation(c) ``` """ function regular_triangulation(P::Polyhedron) - is_fulldimensional(P) || error("Input polytope must be full-dimensional.") - is_bounded(P) || error("Input polytope must be bounded.") - return regular_triangulation(vertices(P); full=false) + is_fulldimensional(P) || error("Input polytope must be full-dimensional.") + is_bounded(P) || error("Input polytope must be bounded.") + return regular_triangulation(vertices(P); full=false) end - - - @doc raw""" secondary_polytope(P::Polyhedron) @@ -431,8 +427,10 @@ julia> sc = secondary_polytope(c) Polytope in ambient dimension 8 ``` """ -function secondary_polytope(P::Polyhedron{T}) where T<:scalar_types - return Polyhedron{T}(Polymake.polytope.secondary_polytope(pm_object(P)), coefficient_field(P)) +function secondary_polytope(P::Polyhedron{T}) where {T<:scalar_types} + return Polyhedron{T}( + Polymake.polytope.secondary_polytope(pm_object(P)), coefficient_field(P) + ) end @doc raw""" @@ -453,14 +451,10 @@ true ``` """ function is_regular(pts::AbstractCollection[PointVector], cells::Vector{Vector{Int64}}) - as_sop = subdivision_of_points(pts,cells) - is_regular(as_sop) + as_sop = subdivision_of_points(pts, cells) + is_regular(as_sop) end - - - - @doc raw""" subdivision_of_points(P::Polyhdron, cells::IncidenceMatrix) @@ -483,8 +477,8 @@ julia> S = subdivision_of_points(C, cells) Subdivision of points in ambient dimension 2 ``` """ -subdivision_of_points(P::Polyhedron, cells::IncidenceMatrix) = subdivision_of_points(vertices(P), cells) - +subdivision_of_points(P::Polyhedron, cells::IncidenceMatrix) = + subdivision_of_points(vertices(P), cells) @doc raw""" subdivision_of_points(P::Polyhdron, weights::AbstractVector) @@ -508,14 +502,17 @@ julia> S = subdivision_of_points(C, weights) Subdivision of points in ambient dimension 2 ``` """ -subdivision_of_points(P::Polyhedron, weights::AbstractVector) = subdivision_of_points(vertices(P), weights) -subdivision_of_points(P::Polyhedron, cells::Vector{Vector{Int64}}) = subdivision_of_points(vertices(P), IncidenceMatrix(cells)) -subdivision_of_points(Iter::SubObjectIterator{<:PointVector}, cells::IncidenceMatrix) = subdivision_of_points(point_matrix(Iter), cells) -subdivision_of_points(Iter::SubObjectIterator{<:PointVector}, weights::AbstractVector) = subdivision_of_points(point_matrix(Iter), weights) -subdivision_of_points(Iter::SubObjectIterator{<:PointVector}, cells::Vector{Vector{Int64}}) = subdivision_of_points(point_matrix(Iter), IncidenceMatrix(cells)) - - - +subdivision_of_points(P::Polyhedron, weights::AbstractVector) = + subdivision_of_points(vertices(P), weights) +subdivision_of_points(P::Polyhedron, cells::Vector{Vector{Int64}}) = + subdivision_of_points(vertices(P), IncidenceMatrix(cells)) +subdivision_of_points(Iter::SubObjectIterator{<:PointVector}, cells::IncidenceMatrix) = + subdivision_of_points(point_matrix(Iter), cells) +subdivision_of_points(Iter::SubObjectIterator{<:PointVector}, weights::AbstractVector) = + subdivision_of_points(point_matrix(Iter), weights) +subdivision_of_points( + Iter::SubObjectIterator{<:PointVector}, cells::Vector{Vector{Int64}} +) = subdivision_of_points(point_matrix(Iter), IncidenceMatrix(cells)) @doc raw""" gkz_vector(SOP::SubdivisionOfPoints) @@ -539,14 +536,14 @@ julia> gkz_vector(Triang) ``` """ function gkz_vector(SOP::SubdivisionOfPoints) - V = SOP.pm_subdivision.POINTS - T = SOP.pm_subdivision.MAXIMAL_CELLS - n = ambient_dim(SOP) - for i in 1:size(T,1) - @assert sum(T[i,:]) == n+1 #poor check that subdivision is triangulation - end - TT = [Polymake.to_zero_based_indexing(Polymake.row(T,i)) for i in 1:Polymake.nrows(T)] - F = coefficient_field(SOP) - tt = elem_type(F) - tt[F(e) for e in Polymake.call_function(:polytope, :gkz_vector, V, TT)] + V = SOP.pm_subdivision.POINTS + T = SOP.pm_subdivision.MAXIMAL_CELLS + n = ambient_dim(SOP) + for i in 1:size(T, 1) + @assert sum(T[i, :]) == n + 1 #poor check that subdivision is triangulation + end + TT = [Polymake.to_zero_based_indexing(Polymake.row(T, i)) for i in 1:Polymake.nrows(T)] + F = coefficient_field(SOP) + tt = elem_type(F) + tt[F(e) for e in Polymake.call_function(:polytope, :gkz_vector, V, TT)] end diff --git a/src/PolyhedralGeometry/type_functions.jl b/src/PolyhedralGeometry/type_functions.jl index 73bfb92b6375..01ab8641de95 100644 --- a/src/PolyhedralGeometry/type_functions.jl +++ b/src/PolyhedralGeometry/type_functions.jl @@ -2,49 +2,68 @@ ######## Scalar types ################################################################################ - -function detect_scalar_type(n::Type{T}, p::Polymake.BigObject) where T<:Union{Polyhedron, Cone, PolyhedralFan, SubdivisionOfPoints, PolyhedralComplex} - typename = Polymake.bigobject_eltype(p) - return typename == "OscarNumber" ? nothing : scalar_type_to_oscar[typename] +function detect_scalar_type( + n::Type{T}, p::Polymake.BigObject +) where {T<:Union{Polyhedron,Cone,PolyhedralFan,SubdivisionOfPoints,PolyhedralComplex}} + 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 +scalar_type( + ::Union{Polyhedron{T},Cone{T},Hyperplane{T},Halfspace{T}} +) where {T<:scalar_types} = T ################################################################################ ######## SubObjectIterator ################################################################################ # Matrices with rational only elements -for (sym, name) in (("linear_inequality_matrix", "Linear Inequality Matrix"), ("affine_inequality_matrix", "Affine Inequality Matrix"), ("linear_equation_matrix", "Linear Equation Matrix"), ("affine_equation_matrix", "Affine Equation Matrix")) - M = Symbol(sym) - _M = Symbol("_", sym) - @eval begin - $M(iter::SubObjectIterator{<:Union{Halfspace{QQFieldElem}, Hyperplane{QQFieldElem}, Polyhedron{QQFieldElem}, Cone{QQFieldElem}, Pair{Matrix{QQFieldElem}, QQFieldElem}}}) = matrix(QQ, Matrix{QQFieldElem}($_M(Val(iter.Acc), iter.Obj; iter.options...))) - $M(iter::SubObjectIterator{<:Union{Halfspace{T}, Hyperplane{T}, Polyhedron{T}, Cone{T}, Pair{Matrix{T}, T}}}) where T<:scalar_types = matrix(coefficient_field(iter.Obj), $_M(Val(iter.Acc), iter.Obj; iter.options...)) - $_M(::Any, ::PolyhedralObject) = throw(ArgumentError(string($name, " not defined in this context."))) - end +for (sym, name) in ( + ("linear_inequality_matrix", "Linear Inequality Matrix"), + ("affine_inequality_matrix", "Affine Inequality Matrix"), + ("linear_equation_matrix", "Linear Equation Matrix"), + ("affine_equation_matrix", "Affine Equation Matrix"), +) + M = Symbol(sym) + _M = Symbol("_", sym) + @eval begin + $M( + iter::SubObjectIterator{ + <:Union{ + Halfspace{QQFieldElem}, + Hyperplane{QQFieldElem}, + Polyhedron{QQFieldElem}, + Cone{QQFieldElem}, + Pair{Matrix{QQFieldElem},QQFieldElem}, + }, + }, + ) = matrix(QQ, Matrix{QQFieldElem}($_M(Val(iter.Acc), iter.Obj; iter.options...))) + $M( + iter::SubObjectIterator{ + <:Union{Halfspace{T},Hyperplane{T},Polyhedron{T},Cone{T},Pair{Matrix{T},T}} + }, + ) where {T<:scalar_types} = + matrix(coefficient_field(iter.Obj), $_M(Val(iter.Acc), iter.Obj; iter.options...)) + $_M(::Any, ::PolyhedralObject) = + throw(ArgumentError(string($name, " not defined in this context."))) + end end - -function halfspace_matrix_pair(iter::SubObjectIterator{<:Union{Halfspace{T}, Hyperplane{T}, Polyhedron{T}, Cone{T}, Pair{Matrix{T}, T}}}) where T<:scalar_types +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 = [f(x) for x in -h[:, 1]]) + 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 ( - cones, - faces, - facets, - maximal_cones, - maximal_polyhedra, - rays, - vertices, - ) - F = Symbol(fun) - @eval $F(::Type{IncidenceMatrix}, x...) = IncidenceMatrix($F(x...)) +for fun in (cones, faces, facets, maximal_cones, maximal_polyhedra, rays, vertices) + F = Symbol(fun) + @eval $F(::Type{IncidenceMatrix}, x...) = IncidenceMatrix($F(x...)) end diff --git a/src/PolyhedralGeometry/visualization.jl b/src/PolyhedralGeometry/visualization.jl index 4afa44ff29fe..e7c9d64175ee 100644 --- a/src/PolyhedralGeometry/visualization.jl +++ b/src/PolyhedralGeometry/visualization.jl @@ -7,32 +7,34 @@ Four-dimensional polytopes are visualized as a Schlegel diagram, which is a proj In higher dimensions there is no standard method; use projections to lower dimensions or try ideas from [GJRW10](@cite). """ -function visualize(P::Union{Polyhedron{T}, Cone{T}, PolyhedralFan{T}, PolyhedralComplex{T}}) where T<:Union{Float64, FieldElem} +function visualize( + P::Union{Polyhedron{T},Cone{T},PolyhedralFan{T},PolyhedralComplex{T}} +) where {T<:Union{Float64,FieldElem}} d = ambient_dim(P) b = P isa Polyhedron if b && d == 4 @req is_fulldimensional(P) "Can only visualize full-dimensional $(typeof(P)) of ambient dimension $d" else - @req d < 4 "Can not visualize $(typeof(P)) of ambient dimension $d. Supported range: 1 <= d <= $(3 + b)" + @req d < 4 "Can not visualize $(typeof(P)) of ambient dimension $d. Supported range: 1 <= d <= $(3 + b)" end # polymake will by default use 0:n-1 as ray labels so we assign labels # starting from 1 here if there are no labels yet # (note: labels are mutable, i.e. they can be changed again later) if !Polymake.exists(pm_object(P), "RAY_LABELS") - pm_object(P).RAY_LABELS = string.(1:Oscar.pm_object(P).N_RAYS) + pm_object(P).RAY_LABELS = string.(1:(Oscar.pm_object(P).N_RAYS)) end pmo = pm_object(P) Polymake.visual(pmo) end -function visualize(P::SubdivisionOfPoints{T}) where T<:Union{FieldElem, Float64} +function visualize(P::SubdivisionOfPoints{T}) where {T<:Union{FieldElem,Float64}} d = ambient_dim(P) @req d <= 3 "Can not visualize $(typeof(P)) of ambient dimension $d. Supported range: 1 <= d <= 3" # polymake will by default use 0:n-1 as labels so we assign labels # starting from 1 here if there are no labels yet # (note: labels are mutable, i.e. they can be changed again later) if !Polymake.exists(pm_object(P), "POINT_LABELS") - pm_object(P).POINT_LABELS = string.(1:Oscar.pm_object(P).N_POINTS) + pm_object(P).POINT_LABELS = string.(1:(Oscar.pm_object(P).N_POINTS)) end pmo = pm_object(P) Polymake.visual(pmo) From f1196ba58591f3ecc357161eb918fa93bb41b711 Mon Sep 17 00:00:00 2001 From: Lars Kastner Date: Thu, 11 Apr 2024 19:09:32 +0200 Subject: [PATCH 2/3] Add workflow for monitoring the formatting --- .github/workflows/JuliaFormatterCI.yml | 30 +++++++ etc/test_formatting.jl | 104 +++++++++++++++++++++++++ 2 files changed, 134 insertions(+) create mode 100644 .github/workflows/JuliaFormatterCI.yml create mode 100644 etc/test_formatting.jl diff --git a/.github/workflows/JuliaFormatterCI.yml b/.github/workflows/JuliaFormatterCI.yml new file mode 100644 index 000000000000..4268285b55d9 --- /dev/null +++ b/.github/workflows/JuliaFormatterCI.yml @@ -0,0 +1,30 @@ +name: JuliaFormatter test + +on: + push: + branches: + - master + - 'release-*' + pull_request: + workflow_dispatch: + +concurrency: + # group by workflow and ref; the last slightly strange component ensures that for pull + # requests, we limit to 1 concurrent job, but for the master branch we don't + group: ${{ github.workflow }}-${{ github.ref }}-${{ github.ref != 'refs/heads/master' || github.run_number }} + # Cancel intermediate builds, but only if it is a pull request build. + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} + +jobs: + check-consistent-formatting: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + - uses: julia-actions/cache@v1 + - name: 'Check format with JuliaFormatter' + run: | + julia etc/test_formatting.jl + shell: bash + env: + jf-version: ${{ inputs.version }} diff --git a/etc/test_formatting.jl b/etc/test_formatting.jl new file mode 100644 index 000000000000..bc9d4805632a --- /dev/null +++ b/etc/test_formatting.jl @@ -0,0 +1,104 @@ +using Pkg +Pkg.activate(temp=true) +Pkg.add("Test") +Pkg.add("JuliaFormatter") +using Test +using JuliaFormatter + +# Disclaimer: +# The goal should be to work with +# https://github.com/julia-actions/julia-format +# . However, currently too many files have broken formatting and since there +# are many ongoing pull requests, we will need to extend proper formatting to +# the entire codebase in a step by step fashion. + + +file = @__FILE__ +oscardir = replace(file, "etc/test_formatting.jl"=>"") + +result = 0 + +@testset "Formatting" begin + + function _gather_source_files(path::String) + queue = [path] + result = String[] + while length(queue) > 0 + current_folder = pop!(queue) + next_batch = readdir(current_folder; join=true) + for elem in next_batch + if isfile(elem) && endswith(elem, ".jl") + push!(result, elem) + elseif isdir(elem) + push!(queue, elem) + end + end + end + return result + end + + # Since right now only very few files are formatted, we also use a whitelist + # approach. + enabled = [ + "src/PolyhedralGeometry", + "src/aliases.jl", + "experimental/LieAlgebras", + "experimental/BasisLieHighestWeight", + ] + skip = [ + "src/PolyhedralGeometry/Polyhedron/standard_constructions.jl", + "src/PolyhedralGeometry/Polyhedron/properties.jl", + "experimental/LieAlgebras/test/AbstractLieAlgebra-test.jl", + "experimental/LieAlgebras/test/LieAlgebraModule-test.jl", + "experimental/LieAlgebras/src/LieAlgebra.jl", + "experimental/LieAlgebras/src/LieAlgebraHom.jl", + "experimental/LieAlgebras/src/LieSubalgebra.jl", + "experimental/LieAlgebras/src/LinearLieAlgebra.jl", + "experimental/LieAlgebras/src/RootSystem.jl", + "experimental/LieAlgebras/src/Util.jl", + "experimental/BasisLieHighestWeight/src/MainAlgorithm.jl", + ] + + # Collect all code files. + entire = [ + _gather_source_files(joinpath(oscardir, "src/")); + _gather_source_files(joinpath(oscardir, "experimental/")) + ] + + failed = String[] + + for file in entire + is_enabled = !isnothing(findfirst(e->occursin(e, file), enabled)) + should_skip = !isnothing(findfirst(e->occursin(e, file), skip)) + if is_enabled && !should_skip + # Do not actually format file, only check whether format is ok. + res = @test format(file; overwrite=false) + if res isa Test.Fail + result = 1 + filename = replace(file, oscardir=>"") + push!(failed, filename) + else + filename = replace(file, oscardir=>"") + end + else + @test format(file; overwrite=false) skip=true + end + end + + # Produce some nice output at the end so users do not have to scroll through + # entire output. + if length(failed) > 0 + println(""" +The files +$failed +were not properly formatted. You can fix the formatting by running +``` +using JuliaFormatter""") + for fn in failed + println("format(\"$fn\")") + end + println("```\n") + end +end + +exit(result) From 1be61d8b63850916c14c1f81337d4c5bcd7f8d2c Mon Sep 17 00:00:00 2001 From: Lars Kastner Date: Thu, 11 Apr 2024 19:19:32 +0200 Subject: [PATCH 3/3] Add .git-blame-ignore-revs so `git blame` does not get diluted from formatting commits --- .git-blame-ignore-revs | 2 ++ etc/test_formatting.jl | 3 +++ 2 files changed, 5 insertions(+) create mode 100644 .git-blame-ignore-revs diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs new file mode 100644 index 000000000000..7e3565b6c6de --- /dev/null +++ b/.git-blame-ignore-revs @@ -0,0 +1,2 @@ +# Ran JuliaFormatter on most of src/PolyhedralGeometry +ab444bbc46b5493588b2c3217ed35d86396d1260 diff --git a/etc/test_formatting.jl b/etc/test_formatting.jl index bc9d4805632a..686cee2a2bc6 100644 --- a/etc/test_formatting.jl +++ b/etc/test_formatting.jl @@ -11,6 +11,9 @@ using JuliaFormatter # . However, currently too many files have broken formatting and since there # are many ongoing pull requests, we will need to extend proper formatting to # the entire codebase in a step by step fashion. +# +# In case you format some code, also add the commit hash to +# .git-blame-ignore-revs for `git blame` to ignore these commits. file = @__FILE__