Combinatorial Criteria for affine sections

This notebook contains code to compute

• combinatorial types of affine sections,

• the affine slice which maximized the number of k-dimensional faces.

Please see the input section below.

Defining some functions

_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)

def remove_duplicates(lis):
    res = []
    for el in lis:
        if not el in res:
            res.append(el)
    return res

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

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

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)]

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]]

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]))

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)

def sweeping_polytope(P):
    Z = sum([Polyhedron(vertices = [vector(pair[0]),vector(pair[1])]) for pair in Subsets(P.vertices(),2)])
    return Z

def is_new_comb_type(Qsample,comb_types):
    for combtype in comb_types:
        if Qsample.is_combinatorially_isomorphic(combtype):
            return False
    return True

Input

Specify your polytope here.

P = polytopes.cross_polytope(3)

Choose the problem you would like to compute.

max-number-k-faces computes the central slice whith the maximum number of k-dimensional faces. k=0 computes the maximum number of vertices.

comb-types computes representatives of all combinatorial types of polytopes that arise as affine slices.

comb_problems = ["max-number-k-faces","comb-types"]
#comb_problem = "comb-types"
comb_problem = "max-number-k-faces"
k = 0 #dimension of k-faces

Computation

variables = [var("x"+str(i)) for i in range(1,P.dimension()+1)]
direction = variables
beta = var("beta")

if not P.dimension() == P.ambient_dim():
    print('the polytope is not full dimensional, the output is not meaningful')

d = P.ambient_dim()
SFan = sweeping_polytope(P).normal_fan()
print("computed the sweeping arrangement")

computed the sweeping arrangement

#get all lower-dimensional cones from the fan
cones = [cone for cones_of_dim in SFan.cones() for cone in cones_of_dim if cone.dim()>0]

samples = []
j = 0
comb_types = []
for region in cones:
    direction_sample = sum([r.sparse_vector() for r in region.rays()])/(sum([r.sparse_vector() for r in region.rays()]).norm())
    Fan = get_intervals_beta(P,direction_sample)

    i = _sage_const_0
    res = []
    for interval in Fan:
        if not j % _sage_const_10 :
            if not i % _sage_const_10 :
                print(j, "of", len(cones), i, 'of', len(Fan))
        i+=_sage_const_1


        #sample on first lower-dimensional face
        bsample = interval[0]
        Qsample = Polyhedron(eqns = [[-bsample]+list(direction_sample)], base_ring = AA).intersection(P)

        if comb_problem == "max-number-k-faces":
            samples.append([len(Qsample.vertices()),region,[bsample],Qsample])
        elif comb_problem == "comb-types" and Qsample.dim()==d-1:
            if is_new_comb_type(Qsample,comb_types):
                comb_types.append(Qsample)
                new_type = Qsample


        #sample on interior of interval
        bsample = (interval[0]+interval[1])/2
        Qsample = Polyhedron(eqns = [[-bsample]+list(direction_sample)], base_ring = AA).intersection(P)
        if comb_problem == "max-number-k-faces":
            samples.append([len(Qsample.vertices()),region,interval,Qsample])
        elif comb_problem == "comb-types" and Qsample.dim()==d-1:
            if is_new_comb_type(Qsample,comb_types):
                comb_types.append(Qsample)
                new_type = Qsample

    #sample on last lower-dimensional face
    bsample = Fan[-1][-1]
    Qsample = Polyhedron(eqns = [[-bsample]+list(direction_sample)], base_ring = AA).intersection(P)
    if comb_problem == "max-number-k-faces":
        samples.append([len(Qsample.vertices()),region,[bsample],Qsample])
    elif comb_problem == "comb-types" and Qsample.dim()==d-1:
        if is_new_comb_type(Qsample,comb_types):
            comb_types.append(Qsample)
    j+=1


print("done")

0 of 146 0 of 2
10 of 146 0 of 1
20 of 146 0 of 1
30 of 146 0 of 4
40 of 146 0 of 3
50 of 146 0 of 3
60 of 146 0 of 4
70 of 146 0 of 3
80 of 146 0 of 3
90 of 146 0 of 4
100 of 146 0 of 5
110 of 146 0 of 5
120 of 146 0 of 5
130 of 146 0 of 5

Output

if comb_problem == "max-number-k-faces":
    output = max(samples)
    print("maximum number of k-faces:", output[0])
    print("region of normal directions:", output[1])
    print("interval of affine values:", output[2])
    print("maximizer:", output[3])
elif comb_problem == "comb-types":
    print(comb_types)