class RelativeSimplicialComplex:
"""This is a class for relative simplicial complexes."""
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
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()
self.RSC_dic = Phi
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()
self._h_counter = self._h_vector
def dimension(self):
return max(self.RSC_dic, key=int)
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)]
def faces(self):
return(self.RSC_dic)
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}))
def max_faces(self):
return self._facets
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}))
def faces_bunch(self):
return Set(reduce(lambda S, T: S | T, [self._faces[c_dim] for c_dim in self._faces]))
def is_pure(self):
if self.dimension() < 0:
return True
return max(self._facets, key=dimension) == min(self._facets, key=dimension)
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
def link(self, face):
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:
return (self.Delta.link(face), self.Gamma.link(face))
else:
return (self.Delta.link(face), None)
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))
def n_copies_of_Delta_over_Gamma(self, n):
if n <= 0:
raise RuntimeError('n must be a positive integer')
if n is 1:
return self.Delta
vert_Gamma = self.Gamma.vertices()
list_Phi = list(self._facets)
list_Gamma = list(self.Gamma.facets())
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))
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:
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) )
def is_partitioned_by(self, interval_partition):
if not interval_partition:
if self:
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), 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 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)))
b3 = Set(F[1] for F in interval_partition) == self._facets
return (b1, b2, b3)
def possible_matchings(self, A, B, h_test=False):
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:
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)
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())
(A, B) = (A - X, B - Y)
p |= q
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
if h_test and any([I_l < 0 for I_l in tmp_count]):
return ({}, Set())
return (match(A, B, bottom_up=False), p)
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 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)
if faces.is_empty():
return (True, Set())
A = get_min_faces(faces)
B = get_max_faces(faces)
(M, p) = self.possible_matchings(A, B, h_test)
if not M and not p:
return (False, Set())
tmp_forced_cntr = self._h_counter[0:] if h_test else None
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
b = (covered, non_overlapping) = p_is_a_partition_of(p, faces)
if all(b):
return (True, p)
elif not non_overlapping:
self._h_counter = tmp_forced_cntr if h_test else None
return (False, Set())
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]:
new_interval = ((I_l, I_u, I_u.dimension() - I_l.dimension()),)
q, tmp_count = Set(p), self._h_counter[0:] if h_test else None
if h_test:
if self._h_counter[I_l.dimension() + 1] - 1 < 0:
return (False, Set())
self._h_counter[I_l.dimension() + 1] -= 1
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)
p = p | r | Set(new_interval)
if all(p_is_a_partition_of(p, faces)):
return (True, p)
p, self._h_counter = Set(q), tmp_count if h_test else None
self._h_counter = tmp_forced_cntr if h_test else None
return (False, Set())
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]
def is_partitionable(self, algorithm="rec", certificate=False, randomize=False):
if algorithm is "ilp":
return self.is_partitionable_ilp(certificate)
elif algorithm is "rec":
pure = self._is_pure
if pure:
h = self._h_vector
f = self._f_vector
if {h_i for h_i in h if h_i < 0}:
return False
if sum(h) != f[self.dimension() + 1]:
return False
(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
format_print = lambda x,y: str(x) + ",\n\t\t" + str(y)
Gadget = SimplicialComplex([
[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"]
,
[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"]
])
print(Gadget.is_shellable(True), Gadget.f_vector(), Gadget.h_vector())
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))
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())
(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)]