# This file was *autogenerated* from the file intersection_body_functions.sage
from sage.all_cmdline import * # import sage library
_sage_const_0 = Integer(0); _sage_const_1 = Integer(1); _sage_const_10 = Integer(10); _sage_const_2 = Integer(2); _sage_const_3 = Integer(3)
# functions for computing the radial function and the algebraic boundary
# of the intersection body of a polytope P
# written in SageMath-9.2
## preliminary functions
def hyperplane_arrangement_to_fan(hyperplane_arrangement):
'''
Input: Linear hyperplane arrangement
Output: Maximal cones of the induced fan. Each maximal cone is the closure of a connected component
of R^d after removal of the hyperplane arrangement (i.e. the realization of the oriented matroid )
'''
n = len(hyperplane_arrangement)
p = list(powerset(range(n)))
d = hyperplane_arrangement[_sage_const_0 ].ambient_dim()
cones = []
for entry in p:
signs = [_sage_const_1 if i in entry else -_sage_const_1 for i in range(n)]
inequalities = []
for i in range(n):
hyperplane = hyperplane_arrangement[i]
A = hyperplane.Hrepresentation()[_sage_const_0 ].A()
ieq = [_sage_const_0 ]+[signs[i]*a for a in A]
inequalities.append(ieq)
C = Polyhedron(ieqs = inequalities)
if C.dimension() == d:
cones.append(C)
return cones
def intersecting_edges(polytope, u):
'''
Input: A polytope and a normal vector u
Output: Set of edges of polytope that has nonempty intersection with the hyperplane orthogonal to u
'''
u_perp = Polyhedron(eqns= [ [_sage_const_0 ]+list(u) ])
intersecting_edges = []
for edge in polytope.faces(_sage_const_1 ):
if edge.as_polyhedron().intersection(u_perp).dimension()>-_sage_const_1 :
intersecting_edges.append(edge)
return intersecting_edges
## main functions
def radial_function_IP(polytope,variables):
'''
computes the radial function of the intersection body of a polytope in any dimension
'''
def regular_subdivision(points, weights, lower_hull = None, upper_hull=None, plot=True):
'''
computes the regular subdidision of a list of points induced by weights
default: lower hull computation
'''
def lower_hull_subdivision(points, weights):
if not len(points)==len(weights):
print('number of points must equal number of weights')
lifted_points = [points[i]+[weights[i]] for i in range(len(points))]
lifted_polytope = Polyhedron(vertices = lifted_points)
if lifted_polytope.dimension() == Polyhedron(vertices = points).dimension():
# no proper subdivision, nothing to do here
return [Polyhedron(vertices = points)]
else:
lower_hull = []
for F in lifted_polytope.facets():
inequality = get_inequality(F.ambient_Hrepresentation())
A = inequality.A()
if A[-_sage_const_1 ]>_sage_const_0 :
lower_hull.append(proj(F))
return lower_hull
def proj(F):
vcs = [v[:-_sage_const_1 ] for v in F.vertices()]
return Polyhedron(vertices = vcs)
def get_inequality(ambient_H_representation_of_facet):
for representation in ambient_H_representation_of_facet:
if representation.is_inequality():
return representation
if (lower_hull and upper_hull) or (lower_hull==False and upper_hull==False):
print('can only compute lower or upper hull. Computing lower hull (default)')
if lower_hull or upper_hull==False:
subdiv_faces = lower_hull_subdivision(points, weights)
elif upper_hull or lower_hull==False:
neg_weghts = [-w for w in weights]
subdiv_faces = lower_hull_subdivision(points, neg_weghts)
else:
subdiv_faces = lower_hull_subdivision(points, weights)
return subdiv_faces
def moment_curve(t,n):
return [t**(i+_sage_const_1 ) for i in range(n)]
P = polytope
if not P.dimension() == P.ambient_dim():
print('the polytope is not full dimensional, the output is not meaningful')
d = P.ambient_dim()
origin = [_sage_const_0 ]*d
x = vector(variables)
print('computing hyperplane arrangement...')
H = list(set([Polyhedron(eqns = [[_sage_const_0 ]+list(vertex)]) for vertex in P.vertices() if not list(vertex)==origin]))
print('completing it to a fan...')
Fan = hyperplane_arrangement_to_fan(H)
alg_plt = None
res = []
print('checking if the origin is inside the polytope...')
if P.contains(origin):
origin_outside_P = False
else:
origin_outside_P = True
i = _sage_const_0
for C in Fan:
if not i % _sage_const_10 :
print(i, 'of', len(Fan))
i+=_sage_const_1
u_sample = sum([vector(ray) for ray in C.rays()])
u_perp = Polyhedron(eqns= [ [_sage_const_0 ]+list(u_sample) ])
intersecting_facets = []
if P.intersection(u_perp).dimension()==d-_sage_const_1 :
for facet in P.facets():
if not facet.as_polyhedron().contains(origin):
if facet.as_polyhedron().intersection(u_perp).dimension() == d-_sage_const_2 :
intersecting_facets.append(facet)
if d<=_sage_const_3 :
intersecting_edges_of_facets = []
for facet in intersecting_facets:
edges = intersecting_edges(facet.as_polyhedron(),u_sample)
intersecting_edges_of_facets.append(edges)
intersecting_edges_of_cells = intersecting_edges_of_facets
else:
#triangulation time
cells_in_triangulation = []
for facet in intersecting_facets:
facet_Q_sample = u_perp.intersection(facet.as_polyhedron())
points = facet_Q_sample.vertices()
regular = False
t = _sage_const_3
# sometimes the moment curve for t=3 does not yield a regular triangulation.
while not regular:
regular = True
weight = moment_curve(t,len(points))
triangulation = regular_subdivision(points, weight, plot=False)
#check that this is an actual triangulation
for cell in triangulation:
if not len(cell.vertices())==(d-_sage_const_1 ):
#print('the regular subdivision is not a triangulation!')
t += _sage_const_1
print('the regular subdivision is not a triangulation! Trying t =',t)
regular = False
cells_in_triangulation+=triangulation
intersecting_edges_of_cells = []
for cell in cells_in_triangulation:
cell_vertices = [list(vertex) for vertex in cell.vertices()]
edges = []
for edge in P.faces(_sage_const_1 ):
inter = u_perp.intersection(edge.as_polyhedron())
if inter.dimension() == _sage_const_0 :
intersection_point = list(inter.vertices()[_sage_const_0 ])
if intersection_point in cell_vertices:
edges.append(edge)
intersecting_edges_of_cells.append(edges)
v = []
for edges_of_cell in intersecting_edges_of_cells:
intersecting_points_of_cells = []
for edge in edges_of_cell:
[a,b] = list(edge.vertices())
a = vector(a)
b = vector(b)
intersection_point = (b.inner_product(x) / (b-a).inner_product(x)) * (a - b) + b
intersecting_points_of_cells.append(intersection_point)
v.append(intersecting_points_of_cells)
determinants = []
for simplex in v:
M = matrix(
simplex + [x]
)
det_M = M.determinant()
if det:
dic = {}
for j in range(d):
dic[variables[j]] = u_sample[j]
sgn = sign(det_M(dic))
if origin_outside_P:
#if the origin is not contained in P, then we have to correct the sign a second time
if P.dimension()<=_sage_const_3 :
intersec = intersecting_facets[v.index(simplex)].as_polyhedron().intersection(u_perp)
simplex_sample_vertices = [list(vert) for vert in intersec.vertices()]
else:
simplex_sample_vertices = [list(vert) for vert in cells_in_triangulation[v.index(simplex)].vertices()]
simplex_sample_vertices.append([_sage_const_0 ]*d)
simplex_sample = Polyhedron(vertices = simplex_sample_vertices)
if not P.intersection(simplex_sample).dimension() == d-_sage_const_1 :
sgn = -sgn
determinants.append([sgn,M.determinant()])
correcting_factor = _sage_const_1 /(factorial(d-1) * sum([x_i**_sage_const_2 for x_i in variables]))
rho = correcting_factor*sum([sng*det for [sng,det] in determinants])
res.append([rho.factor(),C])
else:
res.append([_sage_const_0 ,C])
print('done')
return res
def alg_boundary(radial_fct):
res = []
for rho, C in radial_fct:
if rho:
p, q = rho.numerator_denominator()
alg = q-p
res.append([alg, C])
else:
res.append([rho,C])
return res