Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

One should be able to compute primitive facets #4682

Open
lkastner opened this issue Mar 1, 2025 · 2 comments
Open

One should be able to compute primitive facets #4682

lkastner opened this issue Mar 1, 2025 · 2 comments
Assignees
Labels
enhancement New feature or request topic: polyhedral geometry Issue concerns polyhedral geometry code

Comments

@lkastner
Copy link
Member

lkastner commented Mar 1, 2025

Is your feature request related to a problem? Please describe.
For rays one has:

julia> c = positive_hull([3 1; 0 1])
Polyhedral cone in ambient dimension 2

julia> rays(c)
2-element SubObjectIterator{RayVector{QQFieldElem}}:
 [1, 1//3]
 [0, 1]

julia> matrix(ZZ,rays(c))
[3   1]
[0   1]

I would like something similar for facets. I know there is halfspace_matrix_pair:

julia> halfspace_matrix_pair(facets(c))
(A = [-1 0; 1 -3], b = QQFieldElem[0, 0])

julia> (A,b) = halfspace_matrix_pair(facets(c))
(A = [-1 0; 1 -3], b = QQFieldElem[0, 0])

julia> typeof(A)
QQMatrix

Describe the solution you'd like
Would something like halfspace_matrix_pair(ZZ, facets(c)) be feasible?

Describe alternatives you've considered
I considered spending time with my kids, but filed this issue instead.

Apologies if this already exists in some form, then please let me know.

@lkastner lkastner added the enhancement New feature or request label Mar 1, 2025
@thofma thofma added the topic: polyhedral geometry Issue concerns polyhedral geometry code label Mar 1, 2025
@lkastner
Copy link
Member Author

lkastner commented Mar 3, 2025

@benlorenz made me aware that there are two possibilities in the case of polytopes. Let the facets be given as $Ax\le b$, then one can make $A$ have primitive rows, or even the larger matrix $(A|b)$. Both cases are relevant.

@alexej-jordan
Copy link
Collaborator

Here is the definition of matrix(ZZ, ...) in the context of rays (in iterators.jl):

matrix(R::ZZRing, iter::SubObjectIterator{RayVector{QQFieldElem}}) =
  matrix(R, Polymake.common.primitive(matrix_for_polymake(iter)))

This relies just on polymake's primitive, which takes a polymake compatible matrix, in this case retrieved by matrix_for_polymake. In the case of halfspaces/hyperplanes we can also retrieve such a matrix with affine_matrix_for_polymake (or, depending on the context, linear_matrix_for_polymake). This especially happens within halfspace_matrix_pair (in type_functions.jl). It appears viable to me to extend this method or introduce additional specific methods or functions. This could look something like

...
    if ... //default case
      f = coefficient_field(iter.Obj)
      h = affine_matrix_for_polymake(iter)
      hm = matrix(f, h[:, 2:end])
      bv = [f(x) for x in -h[:, 1]]
    if ... //case (A|b)
      h = Polymake.common.primitive(affine_matrix_for_polymake(iter))
      hm = matrix(ZZ, h[:, 2:end])
      bv = [ZZ(x) for x in -h[:, 1]]
    elseif ... //case A
      f = coefficient_field(iter.Obj)
      h = affine_matrix_for_polymake(iter)
      hm = Polymake.common.primitive(h[:, 2:end])
      bv = Vector{elem_type(f)}(undef, size(hm, 1))
      for i in 1:size(hm, 1)
        j = findfirst(!is_zero, hm[i, :]
        r = hm[i, j]//h[i, j + 1]
        bv[i] = f(h[i, 0] * r)
    end
    return (A=hm, b=bv)
...

This still means we have to decide how the desired case can be set. matrix(ZZ, ...) usually is a special case of matrix(f, ...) with a specific parent f. This cannot be transferred to halfspace_matrix_pair. Nevertheless I think we should let the dispatch handle this, and halfspace_matrix_pair(ZZ, ...) is one way to do this, but here is no way to differentiate between case 2 and 3. We could also introduce new functions named like integral_halfspace_matrix_pair. And the last option i deem possible is using a keyword argument. As this does not utilize the dispatch, it is my least favorite.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request topic: polyhedral geometry Issue concerns polyhedral geometry code
Projects
None yet
Development

No branches or pull requests

3 participants