CoCalc Public FilesrelativeSimplicialComplex.sagews
Views : 355
Compute Environment: Ubuntu 18.04 (Deprecated)


#########################################
# A (N INCIPIENT) CLASS TO DEAL WITH
# RELATIVE SIMPLICIAL COMPLEXES
#########################################
class RelativeSimplicialComplex:
"""This is a class for relative simplicial complexes."""
###
# FUNCTION: constructor
#
# Delta is a SC, and Gamma is a SC in Delta
def __init__(self, _Delta_, _Gamma_=None):
"""
INPUT: Delta, a SC and Gamma a SC

EXAMPLE: Delta = [ [1,3,5],[1,3,6],[1,4,5],[1,4,6],[2,3,5],[2,3,6],[2,4,5] ]

"""
if _Delta_ == _Gamma_:
return None
#raise RuntimeError('Relative simplicial complex is void: It doesn\'t have faces')

self.Delta = copy(_Delta_)
self.Gamma = copy(_Gamma_)

if self.Gamma:
Phi = self.Delta.faces(subcomplex=self.Gamma)
Phi = {dim : Phi[dim] for dim in Phi if Phi[dim]}
else:
Phi = self.Delta.faces()

#
# generates a dictionary from the couple
# of simplicial complexes (Delta, Gamma)
# the key is the dimension and the value
# is a set with the elements in this dimension
self.RSC_dic = Phi

# generates a tuple from the couple
# of simplicial complexes (Delta, Gamma)
self.RSC_tup = (self.Delta, self.Gamma)

self._faces = self.faces()
self._faces_bunch = self.faces_bunch()
self._facets = self.facets()

self._h_vector = self.h_vector()
self._f_vector = self.f_vector()
self._is_pure = self.is_pure()

# used for improve the partitionability algorithm
# in the recursive version
self._h_counter = self._h_vector

###
# FUNCTION: dimension
#
# Dimension of a relative simplical complex
def dimension(self):
return max(self.RSC_dic, key=int)

###
# FUNCTION: f_vector
#
# f-vector of a relative simplical complex
def f_vector(self):
d = self.dimension()
return [0 if face_dim not in self.RSC_dic else len(self.RSC_dic[face_dim]) for face_dim in range(-1,d+1)]

###
# FUNCTION: faces
#
# Faces of the relative simplicial complex
def faces(self):
return(self.RSC_dic)

###
# FUNCTION: facets
#
# Facets of the relative simplicial complex
def facets(self):
return Set(s for s in self._faces_bunch if not filter(lambda x: s.is_face(x), self._faces_bunch - {s}))

###
# FUNCTION: max_faces
#
# Facets of the relative simplicial complex
def max_faces(self):
return self._facets

###
# FUNCTION: min_faces
#
# Minimal faces of the relative simplicial complex
def min_faces(self):
return Set(s for s in self._faces_bunch if not filter(lambda x: x.is_face(s), self._faces_bunch - {s}))

###
# FUNCTION: faces_bunch
#
# Faces of the relative simplicial complex
# disregarding the dimension
def faces_bunch(self):
return Set(reduce(lambda S, T: S | T, [self._faces[c_dim] for c_dim in self._faces]))

###
# FUNCTION: is_pure
#
# Decides purity of a relative simplicial complex
def is_pure(self):
if self.dimension() < 0:
return True
return max(self._facets, key=dimension) == min(self._facets, key=dimension)

###
# FUNCTION: h_vector
#
# h-vector of a relative simplical complex
def h_vector(self):
from sage.arith.all import binomial

d = self.dimension()
f = self.f_vector()
h = [0]*(d+2)

for j in range(0, d+2):
for i in range(-1, j):
h[j] += (-1)**(j-i-1) * binomial(d-i, j-i-1) * f[i+1]

return h

###
#
# Link of a relative simplical complex
bunch = {}
if self.Gamma:
bunch = Set(reduce(lambda S, T: S | T, [self.Gamma.faces()[c_dim] for c_dim in self.Gamma.faces()]))

if face in bunch:
else:

###
# FUNCTION: boundary
#
# Take the boundary of the relative simplicial complex
###
# CAUTION:
# * I'll use the definition I found in "Computer Graphics and Geometric Modelling for
#   Mathematics By Max K. Agoston page 330" for absolute simplicial complexes (see
#   https://goo.gl/FtKzFo) and I extended it for the relative case straightforwardly:
#
#     The boundary of a simplicial complex K, denoted dK is defined by:
#        dK = {t | t is a face of a simplex s^k in K that belongs to a
#                  unique (k+1)-simplex}
#        this definition is applied in a non-pure case
#
#  sage:
#    ->: RelativeSimplicialComplex(SimplicialComplex([[1,2,3],[2,3,4], [1,5]])).facets()
#    <-: {(2, 4), (3, 4), (5,), (1, 3), (1, 2)}
#
#  Prof. Russ Woodroofe suggests
#    <-: {(2, 4), (3, 4), (1, 5), (1, 3), (1, 2)}
#
# * Magma and Macaulay2 get a different approach but equivalent between them
#   Macaulay2:
#    ->: needsPackage "SimplicialComplexes"
#    ->: R = ZZ[a..e];
#    ->: boundary simplicialComplex{a*b*c, b*c*d, a*e}
#    <-: | e cd bd bc ac ab |
#   Magma:
#    ->: dDelta := Boundary(SimplicialComplex([{1,2,3},{1,2,4},{1,5}]));
#    ->: dDelta;
#    <-: Simplicial complex
#           [
#               { 1, 3 },
#               { 1, 4 },
#               { 2, 4 },
#               { 2, 3 },
#               { 1, 2 },
#               { 5 }
#           ]
###
def boundary(self):
from itertools import product

out_faces_by_simplex = Set(Set(Set(F).subsets(F.dimension())) for F in self._facets)
int_faces_for_delete = Set(reduce(lambda A, B: A | B, (S[0] & S[1] for
S in product(out_faces_by_simplex,
out_faces_by_simplex)
if S[0] != S[1])))
dPhi_raw = Set(reduce(lambda A, B: A|B, out_faces_by_simplex)) - int_faces_for_delete
dPhi_cln = Set(S for S in dPhi_raw if not filter(lambda T: S.issubset(T), dPhi_raw - {S}))

return SimplicialComplex(list(dPhi_cln))

###
# FUNCTION: n_disjoint_copies
#
# Given the self Phi = (Delta, Gamma), returns the (absolute) simplicial
# complex nPhi constructed from n disjoint copies of Delta identified
# along the subcomplex Gamma
def n_copies_of_Delta_over_Gamma(self, n):
if n <= 0:
raise RuntimeError('n must be a positive integer')

if n is 1: # by convention
return self.Delta

vert_Gamma = self.Gamma.vertices()

#if not self.Gamma == self.Delta.generated_subcomplex(vert_Gamma):
#    raise RuntimeError('the relative subcomplex doesn\'t fulfill the requirements for this operation...')

list_Phi   = list(self._facets)#list(self._faces_bunch)
list_Gamma = list(self.Gamma.facets())#list(self.Gamma.face_iterator())
list_Rho_n = list_Gamma

for k in range(0, n):
list_Phi_mod = []
for face in list_Phi:
list_Phi_mod.append(map(lambda v: v if (v in vert_Gamma) else str(v)+'_'+str(k+1), face))
list_Rho_n += list_Phi_mod

return(SimplicialComplex(list_Rho_n))

###
# FUNCTION: is_cohen_macaulay
#
# Decides CM-ness of a relative simplical complex.
# This is strongly based on the official sages' function
# is_cohen_macaulay with some improvements I think
# comes from J.L. Martin
###
# CAUTION:
#  * I suggest in (***) to check plain old homology for Delta \ Gamma
#    and relative homology in Delta \cap Gamma
###
def is_cohen_macaulay(self, base_ring=QQ, ncpus=0):
if not self.is_pure():
raise RuntimeError('The SC is non-pure. I should check the behaviour in the non-pure case')

from sage.parallel.decorate import parallel

if not ncpus:
from sage.parallel.ncpus import ncpus as get_ncpus
ncpus = get_ncpus()

facs = [ face for face in self.Delta.face_iterator() ]
n = len(facs)

facs_divided = [ [] for i in range(ncpus) ]
for i in range(n):
facs_divided[i%ncpus].append(facs[i])

def all_homologies_vanish(s):
(lk_s_inDelta, lk_s_inGamma) = self.link(s)

#(***)
if lk_s_inGamma:
H = lk_s_inDelta.homology(subcomplex = lk_s_inGamma, base_ring = base_ring)
else:
H = lk_s_inDelta.homology(base_ring = base_ring)

if base_ring.is_field():
return all( H[i].dimension() == 0 for i in xrange(lk_s_inDelta.dimension()) )
else:
# I'm almost sure, J.L. Martin suggested this line
return any( H[i].invariants() for i in xrange(lk_s_inDelta.dimension()) )

@parallel(ncpus=ncpus)
def all_homologies_in_list_vanish(F):
return all( all_homologies_vanish(s) for s in F )

return all( ans[1] for ans in all_homologies_in_list_vanish(facs_divided) )

###
# FUNCTION: is_partitioned_by
#
def is_partitioned_by(self, interval_partition):
if not interval_partition:
if self:
return (False, True)
else: # remember, if Delta == Gamma is None
return (True, True)

from itertools import combinations

span_by_part = Set(Set(filter(lambda face:
Set(interval[0]).issubset(face),
Set(interval[1]).subsets()))
for interval in interval_partition)

faces_set_n1 = Set(map(lambda face: Set(face), self._faces_bunch))
faces_set_n2 = Set(reduce(lambda face1, face2: face1.union(face2), span_by_part))

b1 = True if faces_set_n1 == faces_set_n2 else False

# if it is covered but the set span_by_part has just one element
# there is not overlapping
if len(span_by_part) == 1:
return (True, True) if b1 else (False, True)  # ???

b2 = reduce(lambda a, b: a and b,
(f_i.intersection(f_j).is_empty() for f_i, f_j in combinations(span_by_part, 2)))

# new condition
b3 = Set(F[1] for F in interval_partition) == self._facets

return (b1, b2, b3)

###
# FUNCTION: possible_matchings
#
# some matchings are forced after some selection
# the function returns a list with the (partial and forced) partition,
# possibly empty, along with the possible matchings
def possible_matchings(self, A, B, h_test=False):  # A and B are sets of faces of a simplicial complex

###
# FUNCTION: match
#
# Tries to find a matching between the facets and the
# minimal faces the simplicial complex
def match(A, B, bottom_up=True):
from itertools import product
if bottom_up:
matching = {}; {matching.setdefault(s, []).append(t) for s, t in product(A, B) if s.is_face(t)}
else:  # top-down
matching = {}; {matching.setdefault(t, []).append(s) for s, t in product(A, B) if s.is_face(t)}
return matching

go_on = True
M, p = {}, Set()

if len(A) > len(B):
return ({}, Set())

while go_on:
M = match(A, B)
#print("    **A: {0} B : {1} M : {2}, h_counter : {3}".format(A, B, M, self._h_counter))

q = Set((key, M[key][0], M[key][0].dimension() - key.dimension()) for key in M if len(M[key]) == 1)

if not q:
go_on = False
else:
(X, Y) = ({x[0] for x in q}, {y[1] for y in q})

if len(X) > len(Y):  # (**)
return ({}, Set())
#if len(X) == len(Y) == 1 and X == Y: # (**) Interval with the form of [I_u, I_l] where I_u == I_l
#    return (match(X,Y), q)

(A, B) = (A - X, B - Y)
p |= q

##### < CAUTION
# ** Note that this routine doesn't modifies the attribute _h_counter
# - test if this matching fits with the h_vector
tmp_count = [self._h_counter[k + 1] - map(lambda x: dimension(x), M).count(k) for k in
range(-1, self.dimension() + 1)] if h_test else None
# - Don't take in account the matching M if it generates lower bounds outside the h_vector
if h_test and any([I_l < 0 for I_l in tmp_count]):
return ({}, Set())
##### > CAUTION

return (match(A, B, bottom_up=False), p)

###
# FUNCTION: is_partitionable_rec
#
# This is a recursive version based on matchings. Recursion is explicit
def is_partitionable_rec(self, faces, h_test=False, random_choice=False):

from random import sample

def get_min_faces(my_set):
return Set(s for s in my_set if not filter(lambda x: x.is_face(s), my_set - {s})) if my_set else None

def get_max_faces(my_set):
return Set(s for s in my_set if not filter(lambda x: s.is_face(x), my_set - {s}) and s in self._facets) if my_set else None

def p_is_a_partition_of(interval_partition, checking_set):
if not interval_partition:
if checking_set:
return (False, True)
else:
return (True, True)

from itertools import combinations

span_by_part = Set(Set(filter(lambda face:
Set(interval[0]).issubset(face),
Set(interval[1]).subsets()))
for interval in interval_partition)

faces_set_n1 = Set(map(lambda face: Set(face), checking_set))
faces_set_n2 = Set(reduce(lambda face1, face2: face1.union(face2), span_by_part))

b1 = True if faces_set_n1 == faces_set_n2 else False

# if it is covered but the set span_by_part has just one element
# there is not overlapping
if len(span_by_part) == 1:
return (True, True) if b1 else (False, True)  # ???

b2 = reduce(lambda a, b: a and b,
(f_i.intersection(f_j).is_empty() for f_i, f_j in combinations(span_by_part, 2)))

return (b1, b2)

# 0. In case of empty complex
if faces.is_empty():
return (True, Set())

# 1. the minimal faces in the relative simplicial complex can
#    force some selections, now we need to decide partitionability
#    in the remaining complex
A = get_min_faces(faces)
B = get_max_faces(faces)
(M, p) = self.possible_matchings(A, B, h_test)

# testing before proceed
if not M and not p:
return (False, Set())

##### < CAUTION
# ** Note that this routine can modify the _h_counter
# - Copy for backtraking
tmp_forced_cntr = self._h_counter[0:] if h_test else None
# - test if this matching fits with the h_vector and modify _h_counter
self._h_counter = [self._h_counter[k + 1] - map(lambda x: dimension(x), [I[0] for I in p]).count(k)
for k in range(-1, self.dimension() + 1)] if h_test else None
##### > CAUTION

# 1. 2. If we are lucky enough, the partition is complete
b = (covered, non_overlapping) = p_is_a_partition_of(p, faces)
if all(b):
return (True, p)
elif not non_overlapping:
##### < CAUTION
self._h_counter = tmp_forced_cntr if h_test else None
##### > CAUTION
return (False, Set())

# 2. Iterate the recursive calling:
for I_u in sample(M, k=len(M)) if random_choice else [I_u for I_u in M]:
for I_l in sample(M[I_u], k=len(M[I_u])) if random_choice else M[I_u]:
# Generate a new interval and included in the partition
new_interval = ((I_l, I_u, I_u.dimension() - I_l.dimension()),)

##### > CAUTION
# ** Note that this routine can modify the _h_counter
# Use h_counter for checking if the selection fits with the h_vector restrictions and update _h_counter
# - Copies for backtracking
q, tmp_count = Set(p), self._h_counter[0:] if h_test else None
# - Test
if h_test:
if self._h_counter[I_l.dimension() + 1] - 1 < 0:
return (False, Set())
# update _h_counter
self._h_counter[I_l.dimension() + 1] -= 1
##### > CAUTION

# The effect of choosing the new_interval
F_c = faces.difference(Set(Simplex(f) for f in Set(I_u).subsets() if Set(I_l).issubset(f)))
(b, r) = self.is_partitionable_rec(F_c, h_test)

# put the interval in the partition along with the return and p
p = p | r | Set(new_interval)
if all(p_is_a_partition_of(p, faces)):
return (True, p)

##### < CAUTION
# - Reset for backtracking
p, self._h_counter = Set(q), tmp_count if h_test else None
##### > CAUTION

##### < CAUTION
self._h_counter = tmp_forced_cntr if h_test else None
##### > CAUTION
return (False, Set())

###
# FUNCTION: is_partitionable_ilp
#
# This is strongly based on the official sages' function
# recently developed for J.L. Martin. It uses an Integer
# Linear program.
def is_partitionable_ilp(self, certificate=False):
from sage.numerical.mip import MixedIntegerLinearProgram

facets = self._facets
faces = self._faces_bunch

RFPairs = [(r, f, len(Set(f)) - len(Set(r)))
for f in facets for r in faces if Set(r).issubset(Set(f))]

n = len(RFPairs)
IP = MixedIntegerLinearProgram()
y = IP.new_variable(binary=True)

for i0, pair0 in enumerate(RFPairs):
for i1, pair1 in enumerate(RFPairs):
if (i0 < i1 and pair0[0].is_face(pair1[1]) and
pair1[0].is_face(pair0[1])):
IP.add_constraint(y[i0] + y[i1] <= 1)

IP.set_objective(sum(2 ** RFPairs[i][2] * y[i] for i in range(n)))
sol = round(IP.solve())

if sol < sum(self._f_vector):
return False
elif not certificate:
return True
else:
x = IP.get_values(y)
return [RFPairs[i] for i in range(n) if x[i] == 1]

###
# FUNCTION: is_partitionable
#
# Decides partitionability of a relative simplical complex
# using either recursion (rec) or an Integer Linear program (ilp)
def is_partitionable(self, algorithm="rec", certificate=False, randomize=False):

if algorithm is "ilp":
return self.is_partitionable_ilp(certificate)
elif algorithm is "rec":
# If it is pure, check conditions on h-vector and f-vector
# TODO. Use the h-triangle and the f-triangle for the non-pure cases
pure = self._is_pure

if pure:
# first check the h-vector
h = self._h_vector
f = self._f_vector

# every entry must be non-negative
if {h_i for h_i in h if h_i < 0}:
return False

if sum(h) != f[self.dimension() + 1]:
return False

# decide partitionability
(b, p) = self.is_partitionable_rec(faces = self._faces_bunch,
random_choice = randomize,
h_test = True if pure else False,
)
return p if b else False

from sage.doctest.util import Timer

#########################################
# Formating prints
format_print = lambda x,y: str(x) + ",\n\t\t" + str(y)

[1, 2, "b"], ["b", "x_1", 2], ["x_1", 2, 3], ["x_1", 3, 1], ["x_1", 1, "a"], ["a", 3, 1], ["a", 3, 2],
["a", 2, "x_0"],[2, "x_0", 3], [3, "x_0", 1], ["x_0", 1, "b"] #,["x_0", "a", "x_1"], ["x_0", "b", "x_1"]
,
[11, 22, "b"], ["b", "x_2", 22], ["x_2", 22, 33], ["x_2", 33, 11], ["x_2", 11, "a"], ["a", 33, 11], ["a", 33, 22],
["a", 22, "x_1"],[22, "x_1", 33], [33, "x_1", 11], ["x_1", 11, "b"]
,
[111, 222, "b"], ["b", "x_0", 222], ["x_0", 222, 333], ["x_0", 333, 111], ["x_0", 111, "a"], ["a", 333, 111], ["a", 333, 222],
["a", 222, "x_2"], [222, "x_2", 333], [333, "x_2", 111], ["x_2", 111, "b"]
])

#########################################
# An example of a non-partitionable simplicial complex
X = RelativeSimplicialComplex(SimplicialComplex([ [1,3,5],[1,3,6],[1,4,5],[1,4,6],[2,3,5],[2,3,6],[2,4,5] ]))
timer = Timer().start()
b = X.is_partitionable(algorithm="ilp", certificate=True)
timer.stop()

print("partitionable_ilp = {0}, time = {2}\n  -- certificate: \n  {1}".format(True, b, timer))

timer = Timer().start()
b = X.is_partitionable(algorithm="rec", certificate=True)
timer.stop()

print("partitionable_rec = {0}, time = {2}\n  -- certificate: \n  {1}".format(True, b, timer))

#(False, [1, 14, 48, 33], [1, 11, 23, -2])
#{'cputime': 0.7386550000000005, 'walltime': 0.752079963684082}
#partitionable_ilp = True, time = {'cputime': 0.7386550000000005, 'walltime': 0.752079963684082}
#  -- certificate:
#  [((4, 5), (1, 4, 5), 1), ((2, 3), (2, 3, 5), 1), ((6,), (2, 3, 6), 2), ((1, 6), (1, 3, 6), 1), ((), (1, 3, 5), 3), ((4,), (1, 4, 6), 2), ((2,), (2, 4, 5), 2)]
#{'cputime': 0.00027099999999968816, 'walltime': 0.00026679039001464844}
#partitionable_com = True, time = {'cputime': 0.00027099999999968816, 'walltime': 0.00026679039001464844}
#  -- certificate:
#  None

#########################################
# Rudin's ball
#########################################

print("\n-------------------\nINSIDE RUDIN'S BALL\n")

bigChunk = SimplicialComplex([
[11, 12, 13, 14], [10, 11, 12, 14], [ 9, 11, 13, 14], [ 8, 12, 13, 14],
[ 7, 11, 12, 13], [ 5,  8, 12, 13], [ 3, 10, 11, 14], [ 3,  6, 10, 14],
[ 3,  6, 10, 11], [ 6, 10, 11, 12], [ 2,  6, 10, 12], [ 2, 10, 12, 14],
[ 2,  4, 10, 14], [ 2,  8, 12, 14], [ 2,  4,  8, 14], [ 2,  6,  8, 12],
[ 4,  8, 13, 14], [ 4, 10, 13, 14], [ 4,  5,  8, 13], [ 4,  5,  8, 12],
[ 6,  8, 11, 12], [ 6, 10, 13, 14], [ 6,  9, 13, 14], [ 4,  8, 11, 12],
[ 4,  7, 11, 12], [ 1,  7, 11, 13], [ 1,  9, 11, 13]]);

noBall = SimplicialComplex([
[ 1,  5,  7, 11], [ 1,  5,  9, 11], [ 5,  9, 11, 14], [ 5,  7, 11, 14],
[ 5,  6,  9, 14],

[ 1,  3,  7, 13], [ 1,  3,  9, 13], [ 3,  9, 12, 13], [ 5,  9, 12, 13],
[ 3,  7, 12, 13], [ 3,  4,  7, 12], [ 5,  6,  9, 13], [ 3,  4,  7, 11],
[ 3,  7, 11, 14]]);

Disk_list = [[ 1,  7, 11], [ 1,  7, 13], [ 1,  9, 11], [ 1,  9, 13], [ 3, 11, 14],
[ 4,  7, 11], [ 4,  7, 12], [ 5, 12, 13], [ 6,  9, 13], [ 6,  9, 14],
[ 7, 12, 13], [ 9, 11, 14]]
Disk = SimplicialComplex(Disk_list);

dbigChunk  = RelativeSimplicialComplex(bigChunk).boundary();
print("dbigChunk = {0}".format(sorted(dbigChunk.facets())))

dnoBall  = RelativeSimplicialComplex(noBall).boundary();
print("dnoBall = {0}".format(sorted(dnoBall.facets())))

ddnoBall = RelativeSimplicialComplex(SimplicialComplex(dnoBall.facets())).boundary();
print("ddnoBall = {0}\n".format(sorted(ddnoBall.facets())))

print("bigChunk:  CM? = {0},  pseudo_MF?= {4}, f_vtr = {1}, h_vtr = {2}, hom = {3}".format
( bigChunk.is_cohen_macaulay(), bigChunk.f_vector(), bigChunk.h_vector(), bigChunk.homology(), bigChunk.is_pseudomanifold()))
print("dbigChunk: CM? = {0},  pseudo_MF?= {4},  f_vtr = {1},     h_vtr = {2},    hom = {3}".format
(dbigChunk.is_cohen_macaulay(),dbigChunk.f_vector(),dbigChunk.h_vector(),dbigChunk.homology(),dbigChunk.is_pseudomanifold()))
print("noBall:    CM? = {0}, pseudo_MF?= {4}, f_vtr = {1}, h_vtr = {2},  hom = {3}".format
(   noBall.is_cohen_macaulay(),   noBall.f_vector(),   noBall.h_vector(),   noBall.homology(),   noBall.is_pseudomanifold()))
print("dnoBall:   CM? = {0}, pseudo_MF?= {4}, f_vtr = {1},     h_vtr = {2},     hom = {3}".format
(  dnoBall.is_cohen_macaulay(),  dnoBall.f_vector(),  dnoBall.h_vector(),  dnoBall.homology(),dnoBall.is_pseudomanifold()))
print("Disk:      CM? = {0},  pseudo_MF?= {4}, f_vtr = {1},     h_vtr = {2},      hom = {3}, shelling_order = {5}".format
(     Disk.is_cohen_macaulay(),     Disk.f_vector(),     Disk.h_vector(),     Disk.homology(),   Disk.is_pseudomanifold(), Disk.is_shellable(certificate=True)))
print("")

X = RelativeSimplicialComplex(noBall, Disk);
timer = Timer().start()
b = X.is_partitionable(algorithm="rec", certificate=True)
timer.stop()

if b:
print("partitionable_rec = {0}, time = {1}\n  -- certificate:".format(True, timer))
for i in sorted(b):
print(i)
print("CM? = {0}, f_vector = {1}, h_vector = {2}, ".format(X.is_cohen_macaulay(), X.f_vector(), X.h_vector()))
print("")

sorted(dbigChunk.facets().intersection(dnoBall.facets()))
sorted(Disk.facets())

#FP = noBall.face_poset()
#FP.plot(figsize=100, border=True, element_shape='s',
#        element_size=100, element_color='white',
#        #element_colors={'pink': F, 'cyan': L.double_irreducibles()},
#        cover_color='lightgray',
#        #cover_colors={'black': F_internal},
#        title='poset space X', fontsize=10)

#D12 = posets.DivisorLattice(12)
#d = {(a, b): b/a for a, b in D12.cover_relations()}
#D12.plot(cover_labels=d, cover_color='gray', cover_style='dotted')

#L = LatticePoset({0: [1, 2, 3, 4], 1: [12], 2: [6, 7],
#                   3: [5, 9], 4: [5, 6, 10, 11], 5: [13],
#                   6: [12], 7: [12, 8, 9], 8: [13], 9: [13],
#                   10: [12], 11: [12], 12: [13]})
#F = L.frattini_sublattice()
#F_internal = [c for c in F.cover_relations() if c in L.cover_relations()]
#L.plot(figsize=12, border=True, element_shape='s',
#        element_size=400, element_color='white',
#        element_colors={'pink': F, 'cyan': L.double_irreducibles()},
#        cover_color='lightgray', cover_colors={'black': F_internal},
#        title='The Frattini\nsublattice in blue', fontsize=10)


(False, [1, 14, 48, 33], [1, 11, 23, -2]) {'cputime': 0.8246079999999978, 'walltime': 0.8372659683227539} partitionable_ilp = True, time = {'cputime': 0.8246079999999978, 'walltime': 0.8372659683227539} -- certificate: [((4, 5), (1, 4, 5), 1), ((2, 3), (2, 3, 5), 1), ((6,), (2, 3, 6), 2), ((1, 6), (1, 3, 6), 1), ((), (1, 3, 5), 3), ((4,), (1, 4, 6), 2), ((2,), (2, 4, 5), 2)] {'cputime': 0.030622000000001037, 'walltime': 0.030647993087768555} partitionable_rec = True, time = {'cputime': 0.030622000000001037, 'walltime': 0.030647993087768555} -- certificate: {((4,), (1, 4, 5), 2), ((2, 4), (2, 4, 5), 1), ((1,), (1, 3, 6), 2), ((4, 6), (1, 4, 6), 1), ((2, 5), (2, 3, 5), 1), ((), (2, 3, 6), 3), ((5,), (1, 3, 5), 2)} ------------------- INSIDE RUDIN'S BALL dbigChunk = [(1, 7, 11), (1, 7, 13), (1, 9, 11), (1, 9, 13), (2, 4, 8), (2, 4, 10), (2, 6, 8), (2, 6, 10), (3, 6, 11), (3, 6, 14), (3, 11, 14), (4, 5, 12), (4, 5, 13), (4, 7, 11), (4, 7, 12), (4, 8, 11), (4, 10, 13), (5, 12, 13), (6, 8, 11), (6, 9, 13), (6, 9, 14), (6, 10, 13), (7, 12, 13), (9, 11, 14)] dnoBall = [(1, 3, 7), (1, 3, 9), (1, 5, 7), (1, 5, 9), (1, 7, 11), (1, 7, 13), (1, 9, 11), (1, 9, 13), (3, 4, 11), (3, 4, 12), (3, 7, 14), (3, 9, 12), (3, 11, 14), (4, 7, 11), (4, 7, 12), (5, 6, 13), (5, 6, 14), (5, 7, 14), (5, 9, 12), (5, 12, 13), (6, 9, 13), (6, 9, 14), (7, 12, 13), (9, 11, 14)] ddnoBall = [()] bigChunk: CM? = True, pseudo_MF?= False, f_vtr = [1, 14, 52, 66, 27], h_vtr = [1, 10, 16, 0, 0], hom = {0: 0, 1: 0, 2: 0, 3: 0} dbigChunk: CM? = True, pseudo_MF?= True, f_vtr = [1, 14, 36, 24], h_vtr = [1, 11, 11, 1], hom = {0: 0, 1: 0, 2: Z} noBall: CM? = False, pseudo_MF?= False, f_vtr = [1, 11, 36, 40, 14], h_vtr = [1, 7, 9, -3, 0], hom = {0: 0, 1: 0, 2: 0, 3: 0} dnoBall: CM? = False, pseudo_MF?= False, f_vtr = [1, 11, 34, 24], h_vtr = [1, 8, 15, 0], hom = {0: 0, 1: Z, 2: Z} Disk: CM? = True, pseudo_MF?= False, f_vtr = [1, 11, 22, 12], h_vtr = [1, 8, 3, 0], hom = {0: 0, 1: 0, 2: 0}, shelling_order = ((1, 7, 11), (1, 7, 13), (1, 9, 11), (9, 11, 14), (7, 12, 13), (4, 7, 11), (4, 7, 12), (1, 9, 13), (6, 9, 13), (5, 12, 13), (6, 9, 14), (3, 11, 14)) {'cputime': 6.274062000000008, 'walltime': 6.304924964904785} partitionable_rec = True, time = {'cputime': 6.274062000000008, 'walltime': 6.304924964904785} -- certificate: ((1, 3), (1, 3, 7, 13), 2) ((1, 5), (1, 5, 9, 11), 2) ((3, 4), (3, 4, 7, 11), 2) ((3, 7), (3, 7, 11, 14), 2) ((3, 9), (1, 3, 9, 13), 2) ((3, 12), (3, 4, 7, 12), 2) ((3, 13), (3, 7, 12, 13), 2) ((5, 6), (5, 6, 9, 13), 2) ((5, 7), (1, 5, 7, 11), 2) ((5, 9), (5, 9, 12, 13), 2) ((5, 11), (5, 9, 11, 14), 2) ((5, 14), (5, 6, 9, 14), 2) ((7, 14), (5, 7, 11, 14), 2) ((9, 12), (3, 9, 12, 13), 2) CM? = True, f_vector = [0, 0, 14, 28, 14], h_vector = [0, 0, 14, 0, 0], [(1, 7, 11), (1, 7, 13), (1, 9, 11), (1, 9, 13), (3, 11, 14), (4, 7, 11), (4, 7, 12), (5, 12, 13), (6, 9, 13), (6, 9, 14), (7, 12, 13), (9, 11, 14)] [(1, 7, 11), (1, 7, 13), (1, 9, 11), (1, 9, 13), (3, 11, 14), (4, 7, 11), (4, 7, 12), (5, 12, 13), (6, 9, 13), (6, 9, 14), (7, 12, 13), (9, 11, 14)]


#########################################
# Ziegler's ball
#########################################

print("\n-------------------\nINSIDE ZIEGLER'S BALL\n")

bigChunk = SimplicialComplex([
[1, 2, 4, 9], [1, 2, 6, 9], [1, 5, 6, 9], [1, 5, 8, 9], [1, 4, 8, 9], [1, 4, 5, 8], [1, 4, 5, 7],
[4, 5, 7, 8], [1, 2, 5, 6], [0, 1, 2, 5], [0, 2, 5, 6], [0, 1, 2, 3], [1, 2, 3, 4], [1, 3, 4, 7]]);
print("bigChunk = {0}".format(sorted(bigChunk.facets())))
shllBall = SimplicialComplex([
[0, 2, 3, 7], [0, 2, 6, 7], [2, 3, 6, 7], [2, 3, 6, 8], [2, 3, 4, 8], [3, 6, 7, 8], [3, 4, 7, 8]]);
print("shllBall = {0}".format(sorted(shllBall.facets())))
dbigChunk = RelativeSimplicialComplex(bigChunk).boundary();
print("dbigChunk = {0}".format(sorted(dbigChunk.facets())))
dshllBall = RelativeSimplicialComplex(shllBall).boundary();
print("dshllBall = {0}".format(sorted(dshllBall.facets())))
patch = shllBall.intersection(bigChunk);
print("patch = bigChunk n shllBall = dbigChunk n dshllBall = {0}".format(sorted(patch.facets())))
semiPatch = SimplicialComplex([
[0, 2, 3], (0, 2, 6), (2, 3, 4), (3, 4, 7)])
print("semiPatch = {0}".format(sorted(semiPatch.facets())))
print("")

print("bigChunk:  CM? = {0}, pseudo_MF?= {4}, shellable?= {5}, f_vtr = {1}, h_vtr = {2}, hom = {3}".format
(bigChunk.is_cohen_macaulay(),
bigChunk.f_vector(), bigChunk.h_vector(),
bigChunk.homology(),
bigChunk.is_pseudomanifold(),
bigChunk.is_shellable()))
print("dbigChunk: CM? = {0}, pseudo_MF?=  {4}, shellable?= {5}, f_vtr = {1},     h_vtr = {2},    hom = {3}".format
(dbigChunk.is_cohen_macaulay(),
dbigChunk.f_vector(), dbigChunk.h_vector(),
dbigChunk.homology(),
dbigChunk.is_pseudomanifold(),
dbigChunk.is_shellable()))
print("shllBall:  CM? = {0}, pseudo_MF?= {4}, shellable?= {5}, f_vtr = {1},   h_vtr = {2}, hom = {3}".format
(shllBall.is_cohen_macaulay(),
shllBall.f_vector(), shllBall.h_vector(),
shllBall.homology(),
shllBall.is_pseudomanifold(),
shllBall.is_shellable()))
print("dshllBall: CM? = {0}, pseudo_MF?=  {4}, shellable?= {5}, f_vtr = {1},      h_vtr = {2},    hom = {3}".format
(dshllBall.is_cohen_macaulay(),
dshllBall.f_vector(), dshllBall.h_vector(),
dshllBall.homology(),
dshllBall.is_pseudomanifold(),
dshllBall.is_shellable()))
print("patch:     CM? = {0}, pseudo_MF?= {4}, shellable?= {5}, f_vtr = {1},       h_vtr = {2},    hom = {3}".format
(patch.is_cohen_macaulay(),
patch.f_vector(), patch.h_vector(),
patch.homology(),
patch.is_pseudomanifold(),
patch.is_shellable()))
print("")

'''
# 1. (shllBall, semiPatch)
#####
print("Relative complex: (shllBall, semiPatch)")
X = RelativeSimplicialComplex(shllBall, semiPatch);

print("CM? = {0}, f_vector = {1}, h_vector = {2}, ".format(X.is_cohen_macaulay(), X.f_vector(), X.h_vector()))

timer = Timer().start()
b = X.is_partitionable(algorithm="rec", certificate=True)
timer.stop()

if b:
print("partitionable_rec = {0}, time = {1}\n  -- certificate:".format(True, timer))
for i in sorted(b):
print(i)
else:
print("partitionable_rec = {0}, time = {1}".format(False, timer))
print("---")

# 1.a. n copies (shllBall, semiPatch)... iterate
#####
for i in range(2,21):
print("{0} copies of (shllBall, semiPatch)".format(i))
Xi = X.n_copies_of_Delta_over_Gamma(i)

print("CM? = {0}, f_vector = {1}, h_vector = {2}, ".format(Xi.is_cohen_macaulay(), Xi.f_vector(), Xi.h_vector()))

timer = Timer().start()
b = RelativeSimplicialComplex(Xi).is_partitionable(algorithm="rec", certificate=True)
timer.stop()

if b:
print("partitionable_rec = {0}, time = {1}\n  -- certificate:".format(True, timer))
for i in sorted(b):
print(i)
else:
print("partitionable_rec = {0}, time = {1}".format(False, timer))
print("---")

# 2. (shllBall, patch)
#####
print("Relative complex: (shllBall, patch)")
Y = RelativeSimplicialComplex(shllBall, patch);

print("CM? = {0}, f_vector = {1}, h_vector = {2}, ".format(Y.is_cohen_macaulay(), Y.f_vector(), Y.h_vector()))

timer = Timer().start()
b = Y.is_partitionable(algorithm="rec", certificate=True)
timer.stop()

if b:
print("partitionable_rec = {0}, time = {1}\n  -- certificate:".format(True, timer))
for i in sorted(b):
print(i)
else:
print("partitionable_rec = {0}, time = {1}".format(False, timer))
print("---")

# 2.a. n copies (shllBall, patch)... iterate
#####
for i in range(2,3):
print("{0} copies of (shllBall, patch)".format(i))
Yi = Y.n_copies_of_Delta_over_Gamma(i)

print("CM? = {0}, f_vector = {1}, h_vector = {2}, ".format(Yi.is_cohen_macaulay(), Yi.f_vector(), Yi.h_vector()))

timer = Timer().start()
b = RelativeSimplicialComplex(Yi).is_partitionable(algorithm="rec", certificate=True)
timer.stop()

if b:
print("partitionable_rec = {0}, time = {1}\n  -- certificate:".format(True, timer))
for i in sorted(b):
print(i)
else:
print("partitionable_rec = {0}, time = {1}".format(False, timer))
print("---")
'''

# 3. (bigChunk, patch)
#####
print("Relative complex: (bigChunk, patch)")

Z = RelativeSimplicialComplex(bigChunk, patch);
print("CM? = {0}, f_vector = {1}, h_vector = {2}, ".format(Z.is_cohen_macaulay(), Z.f_vector(), Z.h_vector()))

timer = Timer().start()
b = Z.is_partitionable(algorithm="rec", certificate=True)
timer.stop()

if b:
print("partitionable_rec = {0}, time = {1}\n  -- certificate:".format(True, timer))
for i in sorted(b):
print(i)
else:
print("partitionable_rec = {0}, time = {1}".format(False, timer))
print("---")


------------------- INSIDE ZIEGLER'S BALL dbigChunk = [(0, 1, 3), (0, 1, 5), (0, 2, 3), (0, 2, 6), (0, 5, 6), (1, 3, 7), (1, 5, 7), (2, 3, 4), (2, 4, 9), (2, 6, 9), (3, 4, 7), (4, 7, 8), (4, 8, 9), (5, 6, 9), (5, 7, 8), (5, 8, 9)] dshllBall = [(0, 2, 3), (0, 2, 6), (0, 3, 7), (0, 6, 7), (2, 3, 4), (2, 4, 8), (2, 6, 8), (3, 4, 7), (4, 7, 8), (6, 7, 8)] patch = bigChunk n shllBall = dbigChunk n dshllBall = [(0, 2, 3), (0, 2, 6), (2, 3, 4), (3, 4, 7), (4, 7, 8)] semiPatch = [(0, 2, 3), (0, 2, 6), (2, 3, 4), (3, 4, 7)] bigChunk: CM? = True, pseudo_MF?= False, shellable?= True, f_vtr = [1, 10, 31, 36, 14], h_vtr = [1, 6, 7, 0, 0], hom = {0: 0, 1: 0, 2: 0, 3: 0} dbigChunk: CM? = True, pseudo_MF?= True, shellable?= True, f_vtr = [1, 10, 24, 16], h_vtr = [1, 7, 7, 1], hom = {0: 0, 1: 0, 2: Z} shllBall: CM? = True, pseudo_MF?= False, shellable?= True, f_vtr = [1, 7, 18, 19, 7], h_vtr = [1, 3, 3, 0, 0], hom = {0: 0, 1: 0, 2: 0, 3: 0} dshllBall: CM? = True, pseudo_MF?= True, shellable?= True, f_vtr = [1, 7, 15, 10], h_vtr = [1, 4, 4, 1], hom = {0: 0, 1: 0, 2: Z} patch: CM? = True, pseudo_MF?= False, shellable?= True, f_vtr = [1, 7, 11, 5], h_vtr = [1, 4, 0, 0], hom = {0: 0, 1: 0, 2: 0} Relative complex: (shllBall, semiPatch) CM? = True, f_vector = [0, 1, 9, 15, 7], h_vector = [0, 1, 6, 0, 0], {'cputime': 19.661591, 'walltime': 21.23753595352173} partitionable_rec = False, time = {'cputime': 19.661591, 'walltime': 21.23753595352173} Two copies of (shllBall, semiPatch) CM? = True, f_vector = [1, 8, 23, 27, 11], h_vector = [1, 4, 5, 1, 0], {'cputime': 0.8820270000000008, 'walltime': 0.995081901550293} partitionable_rec = True, time = {'cputime': 0.8820270000000008, 'walltime': 0.995081901550293} -- certificate: ((), (2, 3, 4, '8_1'), 4) ((0,), (0, 2, 3, 7), 3) ((0, 6), (0, 2, 6, 7), 2) ((2, '8_2'), (2, 3, 4, '8_2'), 2) ((4, 7), (3, 4, 7, '8_1'), 2) ((6,), (2, 3, 6, '8_2'), 3) ((6, 7, '8_2'), (3, 6, 7, '8_2'), 1) ((6, '8_1'), (2, 3, 6, '8_1'), 2) ((7,), (2, 3, 6, 7), 3) ((7, '8_1'), (3, 6, 7, '8_1'), 2) (('8_2',), (3, 4, 7, '8_2'), 3) Three copies of (shllBall, semiPatch) CM? = True, f_vector = [1, 9, 28, 35, 15], h_vector = [1, 5, 7, 2, 0], {'cputime': 8.501892000000002, 'walltime': 9.671139001846313} partitionable_rec = True, time = {'cputime': 8.501892000000002, 'walltime': 9.671139001846313} -- certificate: ((), (2, 3, 4, '8_1'), 4) ((0,), (0, 2, 3, 7), 3) ((0, 6), (0, 2, 6, 7), 2) ((2, 4, '8_2'), (2, 3, 4, '8_2'), 1) ((2, 6), (2, 3, 6, '8_1'), 2) ((2, 7), (2, 3, 6, 7), 2) ((2, '8_2'), (2, 3, 6, '8_2'), 2) ((4, '8_3'), (2, 3, 4, '8_3'), 2) ((6,), (3, 6, 7, '8_2'), 3) ((6, 7, '8_3'), (3, 6, 7, '8_3'), 1) ((6, '8_1'), (3, 6, 7, '8_1'), 2) ((7,), (3, 4, 7, '8_3'), 3) ((7, '8_1'), (3, 4, 7, '8_1'), 2) (('8_2',), (3, 4, 7, '8_2'), 3) (('8_3',), (2, 3, 6, '8_3'), 3) Relative complex: (shllBall, patch) CM? = True, f_vector = [0, 0, 7, 14, 7], h_vector = [0, 0, 7, 0, 0], {'cputime': 13.546595000000003, 'walltime': 14.97779893875122} partitionable_rec = False, time = {'cputime': 13.546595000000003, 'walltime': 14.97779893875122} Two copies of (shllBall, patch) CM? = True, f_vector = [1, 7, 18, 19, 7], h_vector = [1, 3, 3, 0, 0], {'cputime': 0.2323769999999996, 'walltime': 0.24488186836242676} partitionable_rec = True, time = {'cputime': 0.2323769999999996, 'walltime': 0.24488186836242676} -- certificate: ((), (2, 3, 4, 8), 4) ((0,), (0, 2, 3, 7), 3) ((0, 6), (0, 2, 6, 7), 2) ((4, 7), (3, 4, 7, 8), 2) ((6,), (2, 3, 6, 8), 3) ((7,), (2, 3, 6, 7), 3) ((7, 8), (3, 6, 7, 8), 2) Three copies of (shllBall, patch) CM? = True, f_vector = [1, 7, 18, 19, 7], h_vector = [1, 3, 3, 0, 0], {'cputime': 0.19152799999999814, 'walltime': 0.199692964553833} partitionable_rec = True, time = {'cputime': 0.19152799999999814, 'walltime': 0.199692964553833} -- certificate: ((), (2, 3, 4, 8), 4) ((0,), (0, 2, 3, 7), 3) ((0, 6), (0, 2, 6, 7), 2) ((4, 7), (3, 4, 7, 8), 2) ((6,), (2, 3, 6, 8), 3) ((7,), (2, 3, 6, 7), 3) ((7, 8), (3, 6, 7, 8), 2) Relative complex: (bigChunk, patch) CM? = True, f_vector = [0, 3, 20, 31, 14], h_vector = [0, 3, 11, 0, 0],