Open in CoCalc
#########################################
# 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):

#(***)
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.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 combinatorics (com) or using 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="com", certificate=True)
timer.stop()

print("partitionable_com = {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