JTS's profile - activity

I am studying a family of affine monoids. The programs Normaliz and libnormaliz are great for this, and some of their functionality is available in SageMath.

For instance, I can do

C = Cone([[1,1,0,0],[1,0,1,0],[1,0,0,1],[0,1,1,0],[0,1,0,1],[0,0,1,1]])
C.Hilbert_basis()

and from the output

N(1, 1, 0, 0),
N(1, 0, 1, 0),
N(1, 0, 0, 1),
N(0, 1, 1, 0),
N(0, 1, 0, 1),
N(0, 0, 1, 1),
N(1, 1, 0, 1),
N(1, 1, 1, 0),
N(1, 0, 1, 1),
N(0, 1, 1, 1)

one notes that the affine monoid generated by the extreme rays of the cone is not integrally closed (in the lattice Z^4).

However, Normaliz also allows one to find the sublattice that the monoid generates inside Z^4, as well as finding the Hilbert series of the (appropriately graded) affine monoid. Is this functionality exposed in SageMath? I could not find it.

I can certainly drag in libnormaliz and do it, but as you see below, it is cumbersome:

import PyNormaliz
from PyNormaliz import *
n = 4
C = Cone(cone = [[1,1,0,0],[1,0,1,0],[1,0,0,1],[0,1,1,0],[0,1,0,1],[0,0,1,1]], grading=[[int(1) for _ in range(n)]])
C.InternalIndex()
C.Grading()
print_series(C.HilbertSeries())

gives

2
[[1, 1, 1, 1], 1]
(1 -  t + 3t^2 +  t^3)
 ----------------------
 (1 - t) (1 - t^2)^3

This shows that the sublattice has index 2; reasonable, since all generators have even total sum. I got the Hilbert series of the affine monoid, not that of its integral closure --- nice. But I had to randomly convert some integers to python integers, switch from snake_case to CamelCase, overload Cone with the Normaliz version, which has a different syntax. This is irksome.

Questions:

1. Is the functionality for affine monoids already in Sagemath, if so where? I am running sage 106, but I can certainly upgrade.
2. If I need to import/call libnormaliz, is there some way of automagically turn sage integers to the type that libnormaliz prefers, but only when appropriate? In the example above, when inputting a cone, the extremal rays could be given by sage integers, but the grading vector needed python integers.

Edit: Max Alekseyev answered the Hilbert Series question. Guided by this, I found the "The Normaliz backend for polyhedral computations" section of the sagemath docs. It seems the most of the PyNormaliz functionality is implemented, but not the lattice algorithms of libnormaliz? For instance, in plain normaliz, using the input file

amb_space 4
cone_and_lattice 6
1 1 0 0
1 0 1 0
1 0 0 1
0 1 1 0
0 1 0 1 
0 0 1 1

total_degree
NoGradingDenom

HilbertBasis
HilbertSeries
Sublattice

I get at the end of the output file:

1 congruences:
1 1 1 1 2

4 basis elements of generated  lattice:
1 0 0 1
0 1 0 1
0 0 1 1
0 0 0 2

This tells me that all lattice points in the affine monoid lies in the sublattice consisting of those lattice points (a,b,c,d) with a+b+c+d congruent 0 mod 2, and that there are no further restrictions.

How can I use sagemath to do this calculation?

Edit 2: The functionality for finding the generated sublattice is available in PyNormaliz (somewhat poorly documented, I had to experiment):

import PyNormaliz
from PyNormaliz import *
n = 4
C = Cone(cone_and_lattice = [[1,1,0,0],[1,0,1,0],[1,0,0,1],[0,1,1,0],[0,1,0,1],[0,0,1,1]], grading=[[int(1) for _ in range(n)]])
print_matrix(C.Congruences())
print(C.Sublattice())

prints

 1 1 1 1 2
 [[[1, 0, 0, 1], [0, -1, -1, 0], [0, 1, 0, 1], [0, 0, 1, 1]], [[2, -1, -1, -1], [0, -1, 1, -1], [0, -1, -1, 1], [0, 1, 1, 1]], 2]

Yes, since it seems the desired functionality is already in PyNormaliz, see my second edit, it should be possible to add

Analyzing affine monoids I am studying a family of affine monoids. The programs Normaliz and libnormaliz are great for t

Thank you, that answers one of my questions. I would also like to find the lattice gp(M) inside Z^4 generated by the aff

Analyzing affine monoids I am studying a family of affine monoids. The programs Normaliz and libnormaliz are great for t

Analyzing affine monoids I am studying a family of affine monoids. The programs Normaliz and libnormaliz are great for t

Dear all, for polytopes there is a method P.is_combinatorially_isomorphic(Q) which tests if two polytopes are combinatorially equivalent, i.e. if their face lattices are isomorphic.

I would like to test two 0/1-polytopes, i.e. polytopes which are the convex hull of a subset of the vertices of a (canonically embedded) hypercube $Q_n \subset [0,1]^n$, for 0/1-equivalence. This means that there is an automorphism of $Q_n$ that maps $P$ to $Q$. Explicitly, the automorphism should be a combination of permuting the variables and flipping $x_i \to 1 - x_i$.

A reference is Ziegler, Lectures on 0/1-polytopes. He also mentions:

The full-dimensional 0/1-polytopes of dimension d = 4 were first enumerated by Alexx Below: There are 349 different 0/1-equivalence classes. In dimension 5 there are exactly 1226525 different 0/1-equivalence classes of 5-dimensional 0/1-polytopes. This classification was done by Oswin Aichholzer [2]: a considerable achievement, which was possible only by systematic use of all the symmetry that is in- herent in the problem.

It would be great to have access to a database of 0/1-polytopes, and be able to compare all facets of dimension at most 4 of some polytope P to candidates in the database, and classify them, something similar to

G.structure_description()

but for 0/1-polytopes.

Of secondary importance (to me) would be checking two polytopes for for affine equivalence and congruence.

Is this implemented i Sagemath?

Edit: the answer by Max Alekseyev points me in the right direction. I will investigate the polydb and polymake docs. Thanks for the help!

Dear all, I want to write a very simple class that inherits from existing polytopes functionality, creates one particular type of polytopes, and adds some functionality which is only relevant for this type of polytopes.

However, I am a python beginner, and am having a hard time navigating the sagemath src directory to find exactly from which class (Polyhedron or Polyhedra) that I should inherit, and how to do this.

Here is a toy example: suppose I want to investigate "depleted cubes", ie polytopes combinatorially equivalent with hypercubes with one corner neatly sliced off. I also want to easily retrieve the cut-off facet.

This works:

def depleted_cube(n: Integer) -> Polyhedron:
       """Returns the n-dimensional hypercube minus one vertex"""
       H = polytopes.hypercube(n)
       Verts = [v for v in H.vertices_list() if add(v) < n]
       C = Polyhedron(vertices=Verts)
       return C

def _strange_facet(P: Polyhedron):
      """Returns the strange facet of P, but only of P depleted cube"""
      n = P.dimension()
      L = [F for F in P.facets() if len(F.vertices()) == n]
      return L[0]

I can now go

sage: _strange_facet(depleted_cube(6)).vertices()
(A vertex at (-1, 1, 1, 1, 1, 1),
A vertex at (1, -1, 1, 1, 1, 1),
A vertex at (1, 1, -1, 1, 1, 1),
A vertex at (1, 1, 1, -1, 1, 1),
A vertex at (1, 1, 1, 1, -1, 1),
A vertex at (1, 1, 1, 1, 1, -1))

but I would like to be able to do

polytopes.depleted_cube(6).strange_facet().vertices()

so as not to dilute the precious namespace! I.e., the function should be a "method" to some "class" that inherits from, either Polyhedron or Polyhedra or polyhedron_base_with_backend_cpp or ... Can someone provide an example, skeleton implementation of a polytope class that inherits from somewhere? As I mentioned, I tried looking at the src directory, but the relevant files are extremely long and littered with docstrings and doctests and cruft; I can make neither head nor tails of them.

Edit: As per the answer by dan_fulea, it is not worthwhile to try to turn the functions into a method. I will undogmatically refrain from doing so, then. Thank you for the tour through the internals of sagemath; I was not aware of the f.?? method of getting the code of a function/object.

How to inherit from polytopes Dear all, I want to write a very simple class that inherits from existing polytopes functi

0/1-equivalence for polytopes Dear all, for polytopes there is a method P.is_combinatorially_isomorphic(Q) which tests

This looks promising, thanks! Many useful refeernces indeed. Poking around, I found clear instructions for accessing po

0/1-equivalence for polytopes Dear all, for polytopes there is a method P.is_combinatorially_isomorphic(Q) which tests

Here is my attempt: def smallerB(a,B,P): if len([b for b in B if P.is_lequal(a,b)]) > 0: return true else: r

How to inherit from polytopes Dear all, I want to write a very simple class that inherits from existing polytopes functi

Multivariate Pade approximation, recognizing rational function (Question edited) Suppose that I have a taylor polynomia

Let me show my work so far. S.<y> = QQ[] R.<x> = S[] T.<y,x> = PowerSeriesRing(QQ,default_prec=8) U.&

Multivariate Pade approximation, recognizing rational function Suppose that I have a taylor polynomial approximation of

Multivariate Pade approximation, recognizing rational function Suppose that I have a taylor polynomial approximation of

Multivariate Pade approximation, recognizing rational function Suppose that I have a taylor polynomial approximation of

Multivariate Pade approximation, recognizing rational function Suppose that I have a taylor polynomial approximation of

Multivariate Pade approximation, recognizing rational function Suppose that I have a taylor polynomial approximation of

Your matrix is not symmetric: A = matrix(QQbar, [[2,-1,0],[-1,2,1],[0,-1,2]]) A.is_symmetric()

Adding dot2tex to the Sage installation, as suggested by FrédéricC and John Palmieri, did work for me and produces Hass

latex for Hasse diagram of poset not properly laid out The plot method on a poset displays the Hasse diagram in a satisf

latex for Hasse diagram of poset not properly laid out The plot method on a poset displays the Hasse diagram in a satisf

You could expand around t=0 instead, for t = w + 0.1 var('t') bmac = b.substitute(w=t-0.1).taylor(t,0,3) cmac = c.subst

You could expand around t=0 instead, for t = w + 0. var('t') bmac = b.substitute(w=t-0.1).taylor(t,0,3) cmac = c.substi

Displaying and typesetting partitions A partition can be viewed as a list or pretty-printed using ascii art. This works

I am trying to simplify some binomial expressions using Sage.

var('a,n,c,m')
H = binomial(2*m-c,m-c)-binomial(2*m-c,m+1)
Cn = binomial(2*n,n)/(n+1)
Uan = (H.substitute(m=n,c=a+1) + a*H.substitute(m=n-1,c=a))/Cn 
assume(a>0)
assume(n>0)
assume(n,'integer')
assume(a,'integer')
UUan=Uan.full_simplify()
UUan

Now I want the limit of this expression as n -> +Infininty. Maple can do it, it correctly returns

                                                             2                  (-a)
                                                          (a  + 3 a + 4) 2
                                                          --------------------
                                                                   4

I can get Sage to calculate the limit for specific values of a:

limit(UUan.substitute(a=5),n=+Infinity)

However, when I try

limit(UUan,n=+Infinity)

I get an error from the underlying Maxima engine:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-70-59aa81a6cc92> in <module>
----> 1 limit(UUan,n=+Infinity)

~/sage/local/var/lib/sage/venv-python3.10/lib/python3.10/site-packages/sage/calculus/calculus.py in limit(ex, dir, taylor, algorithm, **argv)
   1415     if algorithm == 'maxima':
   1416         if dir is None:
-> 1417             l = maxima.sr_limit(ex, v, a)
   1418         elif dir in dir_plus:
   1419             l = maxima.sr_limit(ex, v, a, 'plus')

~/sage/local/var/lib/sage/venv-python3.10/lib/python3.10/site-packages/sage/interfaces/maxima_lib.py in sr_limit(self, expr, v, a, dir)
    985             elif dir == "minus":
    986                 L.append(max_minus)
--> 987             return max_to_sr(maxima_eval(([max_limit], L)))
    988         except RuntimeError as error:
    989             s = str(error)

~/sage/local/var/lib/sage/venv-python3.10/lib/python3.10/site-packages/sage/libs/ecl.pyx in sage.libs.ecl.EclObject.__call__ (build/cythonized/sage/libs/ecl.c:8511)()
    836         """
    837         lispargs = EclObject(list(args))
--> 838         return ecl_wrap(ecl_safe_apply(self.obj,(<EclObject>lispargs).obj))
    839 
    840     def __richcmp__(left, right, int op):

~/sage/local/var/lib/sage/venv-python3.10/lib/python3.10/site-packages/sage/libs/ecl.pyx in sage.libs.ecl.ecl_safe_apply (build/cythonized/sage/libs/ecl.c:6053)()
    357             raise KeyboardInterrupt("ECL says: {}".format(message))
    358         else:
--> 359             raise RuntimeError("ECL says: {}".format(message))
    360     else:
    361         return ret

RuntimeError: ECL says: BINDING-STACK overflow at size 10240. Stack can probably be resized.
Proceed with caution.

---------------------------------------------------------------------------------------------------------------------------------------------

Is SAGE able to compute this limit? Should i provide other provisos? Is there another approach I should try?

Edit: The following works:

Uan.limit(n=Infinity,dir='+',taylor=True)

so I guess all is well. I do not know if the fact that Maxima throws an exception when the taylor keyword is ommitted is a bug or not. The documentation for limit says

• "taylor" - (default: False); if True, use Taylor series, which allows more limits to be computed (but may also crash in some obscure cases due to bugs in Maxima).

Maybe this should be amended to indicate that in some cases, taylor=False causes crashes.

The function limit() has a keyword taylor, which might help or hinder. In this case, Uan.limit(n=Infinity,dir='+',taylo

The function limit() has a keyword taylor, which might help or hinder. In this case, Uan.limit(n=Infinity,dir='+'