Integrals over affine sections
This notebook contains code to compute rational function expressions for
the volume of all affine sections,
the integral of a polynomial over all affine sections.
Please see the input section below.
Defining some functions
[1]:
_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)
[2]:
def remove_duplicates(lis):
res = []
for el in lis:
if not el in res:
res.append(el)
return res
[3]:
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
'''
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
[4]:
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
[5]:
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)]
[6]:
def get_parametric_intervals_beta(P,direction_sample):
values_vcs = sorted([[vector(direction_sample).dot_product(vector(v)),vector(v)] for v in P.vertices()])
return [[vector(variables).dot_product(values_vcs[i][1]),vector(variables).dot_product(values_vcs[i+1][1])] for i in range(len(values_vcs)-1) if not values_vcs[i][0]==values_vcs[i+1][0]]
[7]:
def get_intervals_beta(P,direction):
values = sorted([vector(direction).dot_product(vector(v)) for v in P.vertices()])
return [[values[i],values[i+1]] for i in range(len(values)-1) if not values[i]==values[i+1]]
[8]:
def parametric_interval_to_string(interval):
return str(interval[0])+" <= beta, beta <= "+str(interval[1])
[9]:
def intersecting_edges(polytope,cell):
#get all edges of P that have nonempty intersection with cell
intersecting_edges = []
for edge in polytope.faces(_sage_const_1 ):
if edge.as_polyhedron().intersection(cell).dimension()>-_sage_const_1 :
intersecting_edges.append(edge)
return intersecting_edges
def get_parametrized_vertex(edge):
[a,b] = list(edge.vertices())
a = vector(a)
b = vector(b)
u = vector(variables)
intersection_point = (beta*(b-a) + u.dot_product(b)*a - u.dot_product(a)*b) / u.inner_product(b-a)
return intersection_point
def get_parametrized_vertices(simplex):
intersecting_edges_cell = intersecting_edges(P,simplex)
vertices_of_cell = [get_parametrized_vertex(edge) for edge in intersecting_edges_cell]
return vertices_of_cell
def get_detM(vertices_of_cell,bsample,direction_sample):
M = matrix([v-vertices_of_cell[0] for v in vertices_of_cell[1:]]+[list(vector(direction)/vector(direction).norm())])
det_M = M.determinant()
if det_M:
dic = {}
for j in range(d):
dic[direction[j]] = direction_sample[j]
dic[beta] = bsample
sgn = sign(det_M(dic))
return sgn*det_M
def integrate_over_simplex(h,simplex,bsample,direction_sample):
simplex_full = get_parametrized_vertices(simplex)
integral_on_simplex = 0
for index in range(len(h.exponents())):
alpha = h.exponents()[index]
c_alpha = h.coefficients()[index]
PP = list(cartesian_product([range(a+1) for a in alpha]))
t = sum(alpha)
K = [list(k) for k in LatticePolytope([list(v) for v in (t*polytopes.simplex(d-1)).vertices()]).points()]
factor = c_alpha / factorial(t + d-1)
intermediate_res = 0
for p in PP:
inner_sum = sum([prod([vector(p).dot_product(vector(simplex_full[j]))^k[j] for j in range(len(simplex_full))]) for k in K])
intermediate_res += (-1)^(t-sum(p))*prod([binomial(alpha[j],p[j]) for j in range(len(alpha))])*inner_sum
integral_on_simplex += factor * intermediate_res
return get_detM(simplex_full,bsample,direction_sample)*integral_on_simplex*1/sqrt(sum([x^2 for x in direction]))
[10]:
def sweeping_arrangement(P):
arrangement = []
for pair in Subsets(P.vertices(),2):
H = Polyhedron(eqns = [[0]+list(vector(pair[0])-vector(pair[1]))])
if not H in arrangement:
arrangement.append(H)
return arrangement
def sweeping_fan(P):
AA = sweeping_arrangement(P)
return hyperplane_arrangement_to_fan(AA)
[11]:
def get_inequalities_as_str(region,variables):
return " ,".join([str(ieq.A()*vector(variables) + ieq.b())+" >= 0" for ieq in region.inequalities()])
Input
Specify your polytope here.
[12]:
P = Polyhedron(vertices=[ [-1,-1], [1,-1], [1,1], [0,2], [-1,1]])
[13]:
variables = [var("x"+str(i)) for i in range(1,P.dimension()+1)]
R = PolynomialRing(QQ,variables)
Specify your polynomial here. If you want to compute the volume, choose h=R(1).
[14]:
h = R(1)
Computation
[15]:
direction = variables
beta = var("beta")
d = P.ambient_dim()
[16]:
if not P.dimension() == P.ambient_dim():
print('the polytope is not full dimensional, the output is not meaningful')
SFan = sweeping_fan(P)
print("computed the sweeping arrangement")
computed the sweeping arrangement
[17]:
final_res = []
j = 0
for region in SFan:
direction_sample = sum([vector(r) for r in region.rays()])/sum([vector(r) for r in region.rays()]).norm()
Fan = get_intervals_beta(P,direction_sample)
Fan_parametric = get_parametric_intervals_beta(P,direction_sample)
i = _sage_const_0
#res = []
for k in range(len(Fan)):
if not i % _sage_const_10 :
print(j, "of", len(SFan), i, 'of', len(Fan))
i+=_sage_const_1
interval = Fan[k]
interval_parametric = Fan_parametric[k]
bsample = (interval[0]+interval[1])/2
Qsample = Polyhedron(eqns = [[-bsample]+list(direction_sample)], base_ring = AA).intersection(P)
#triangulation time!
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(Qsample.vertices()))
triangulation = regular_subdivision(Qsample.vertices(), weight, plot=False)
#check that this is an actual triangulation
for cell in triangulation:
if not len(cell.vertices())==d:
t += _sage_const_1
print('the regular subdivision is not a triangulation! Trying t =',t)
regular = False
integral = sum([integrate_over_simplex(h,cell,bsample,direction_sample) for cell in triangulation])
final_res.append(str(integral.simplify_full().factor())+", "+parametric_interval_to_string(interval_parametric)+", "+get_inequalities_as_str(region,variables))
j+=1
print("done")
0 of 12 0 of 4
1 of 12 0 of 4
2 of 12 0 of 4
3 of 12 0 of 4
4 of 12 0 of 4
5 of 12 0 of 4
6 of 12 0 of 4
7 of 12 0 of 4
8 of 12 0 of 4
9 of 12 0 of 4
10 of 12 0 of 4
11 of 12 0 of 4
done
Output
The output is stored in the list final_res. Each entry of final_res is again a list: The first entry of this list is a rational function, and the remaining entries are inequalities on the variables on which the rational function is valid. These inequalities define the corresponding region and chamber. We assume that x is a unit vector, i.e. additionally satisfies x1^2 + ... + xd^2 == 1.
[18]:
final_res[0]
[18]:
'(beta - x1 - x2)/((x1 - x2)*x2), x1 + x2 <= beta, beta <= x1 - x2, -x2 >= 0 ,-x1 + 3*x2 >= 0'
[19]:
final_res[1]
[19]:
'(beta - 3*x1 + x2)/((x1 - x2)*x1), x1 - x2 <= beta, beta <= 2*x2, -x2 >= 0 ,-x1 + 3*x2 >= 0'