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 ### # FUNCTION: link # # Link of a relative simplical complex 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) ### # 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 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) 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"] #,["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"] ]) print(Gadget.is_shellable(True), Gadget.f_vector(), Gadget.h_vector()) ######################################### # 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