################################################################## ################################################################## ################################################################## #### #### #### GAP-Program #### #### #### #### `DiscreteMorse' #### #### #### #### Version May/2015 #### #### by Frank H. Lutz, TU Berlin, Germany #### #### #### #### available at #### #### #### #### http://page.math.tu-berlin.de/~lutz/stellar/ #### #### #### ################################################################## ################################################################## ################################################################## ################################################################## ################################################################## ## ## ## The GAP-program `DiscreteMorse' provides three heuristics ## ## - random_random ## ## - random_lex_first ## ## - random_lex_last ## ## to search for discrete Morse vectors with few critical cells ## ## for a given simplicial complex. ## ## ## ## For an exposition of the program see the reference: ## ## ## ## B. Benedetti and F. H. Lutz. ## ## Random discrete Morse theory and a new library ## ## of triangulations. ## ## Exp. Math. 23, 66-94 (2014). ## ## ## ## Further references: ## ## ## ## K. A. Adiprasito, B. Benedetti, and F. H. Lutz. ## ## Extremal examples of collapsible complexes ## ## and random discrete Morse theory. ## ## Preprint, 20 pages, 2014; arXiv:1404.4239. ## ## ## ## B. Benedetti and F. H. Lutz. ## ## Library of Triangulations, 2013-2015, ## ## http://page.math.tu-berlin.de/~lutz/stellar/ ## ## library_of_triangulations/ ## ## ## ## M. Joswig, F. H. Lutz, and M. Tsuruga. ## ## Sphere Recognition: Heuristics and Examples. ## ## Preprint, 30 pages, 2015; arXiv:1302.6856v2. ## ## ## ## The GAP Group. ## ## GAP -- Groups, Algorithms, and Programming. ## ## Version 4.7.7, 2015, ## ## http://www.gap-system.org ## ## ## ################################################################## ################################################################## Print(StringTime(Runtime()),"\n"); ### fix random sequence ### R_N:=2; # <--- make your choice to alter random sequence R_X:=[66318732,86395905,22233618,21989103,237245480,264566285,240037038, 264902875,9274660,180361945,94688010,24032135,106293216,27264613, 126456102,243761907,80312412,2522186,59575208,70682510,228947516, 173992210,175178224,250788150,73030390,210575942,128491926,194508966, 201311350,63569414,185485910,62786150,213986102,88913350,94904086, 252860454,247700982,233113990,75685846,196780518,74570934,7958751, 130274620,247708693,183364378,82600777,28385464,184547675,20423483, 75041763,235736203,54265107,49075195,100648387,114539755]; # <--- make your choice to alter random sequence ### adjust parameters ### rounds:=1000; # <--- fix number of rounds flag_manifold:=1; # <--- value for closed manifold = 1 # value for general complex = 0 strategy:=1; # <--- value for random_random = 1 # value for random_lex_first = 2 # value for random_lex_last = 3 ### maximal faces ### # example: double_trefoil facets:=[ [ 1, 2, 5, 6 ], [ 1, 2, 5, 12 ], [ 1, 2, 6, 12 ], [ 1, 3, 7, 8 ], [ 1, 3, 7, 11 ], [ 1, 3, 8, 11 ], [ 1, 4, 5, 6 ], [ 1, 4, 5, 16 ], [ 1, 4, 6, 12 ], [ 1, 4, 10, 13 ], [ 1, 4, 10, 16 ], [ 1, 4, 12, 13 ], [ 1, 5, 12, 13 ], [ 1, 5, 13, 16 ], [ 1, 7, 8, 9 ], [ 1, 7, 9, 11 ], [ 1, 8, 9, 14 ], [ 1, 8, 10, 14 ], [ 1, 8, 10, 15 ], [ 1, 8, 11, 15 ], [ 1, 9, 11, 15 ], [ 1, 9, 14, 15 ], [ 1, 10, 13, 14 ], [ 1, 10, 15, 16 ], [ 1, 13, 14, 16 ], [ 1, 14, 15, 16 ], [ 2, 3, 4, 13 ], [ 2, 3, 4, 15 ], [ 2, 3, 13, 15 ], [ 2, 4, 7, 8 ], [ 2, 4, 7, 15 ], [ 2, 4, 8, 16 ], [ 2, 4, 10, 13 ], [ 2, 4, 10, 16 ], [ 2, 5, 6, 14 ], [ 2, 5, 12, 14 ], [ 2, 6, 8, 12 ], [ 2, 6, 8, 16 ], [ 2, 6, 9, 14 ], [ 2, 6, 9, 16 ], [ 2, 7, 8, 9 ], [ 2, 7, 9, 10 ], [ 2, 7, 10, 13 ], [ 2, 7, 13, 15 ], [ 2, 8, 9, 14 ], [ 2, 8, 12, 14 ], [ 2, 9, 10, 16 ], [ 3, 4, 12, 13 ], [ 3, 4, 12, 15 ], [ 3, 5, 6, 7 ], [ 3, 5, 6, 14 ], [ 3, 5, 7, 8 ], [ 3, 5, 8, 11 ], [ 3, 5, 11, 14 ], [ 3, 6, 7, 16 ], [ 3, 6, 9, 14 ], [ 3, 6, 9, 16 ], [ 3, 7, 11, 14 ], [ 3, 7, 14, 16 ], [ 3, 9, 12, 13 ], [ 3, 9, 12, 16 ], [ 3, 9, 13, 15 ], [ 3, 9, 14, 15 ], [ 3, 12, 15, 16 ], [ 3, 14, 15, 16 ], [ 4, 5, 6, 7 ], [ 4, 5, 7, 8 ], [ 4, 5, 8, 16 ], [ 4, 6, 7, 15 ], [ 4, 6, 12, 15 ], [ 5, 8, 11, 13 ], [ 5, 8, 13, 16 ], [ 5, 11, 12, 13 ], [ 5, 11, 12, 14 ], [ 6, 7, 13, 15 ], [ 6, 7, 13, 16 ], [ 6, 8, 12, 15 ], [ 6, 8, 13, 15 ], [ 6, 8, 13, 16 ], [ 7, 9, 10, 12 ], [ 7, 9, 11, 12 ], [ 7, 10, 12, 14 ], [ 7, 10, 13, 14 ], [ 7, 11, 12, 14 ], [ 7, 13, 14, 16 ], [ 8, 10, 12, 14 ], [ 8, 10, 12, 15 ], [ 8, 11, 13, 15 ], [ 9, 10, 12, 16 ], [ 9, 11, 12, 13 ], [ 9, 11, 13, 15 ], [ 10, 12, 15, 16 ] ]; ################################################################### ################################################################### ### Main Part ### ################################################################### ################################################################### morse_vector_collection:=[]; morse_vector_count:=[]; if strategy = 1 then rounds_A:=rounds; rounds_B:=1; elif strategy = 2 or strategy =3 then rounds_A:=1; rounds_B:=rounds; fi; for round_B in [1..rounds_B] do if strategy = 2 or strategy = 3 then ########################## # relabeling of vertices # ########################## vertices:=[]; for facet in facets do UniteSet(vertices,facet); od; new_vertices:=[1..Length(vertices)]; new_labels:=[]; for i in [1..Length(vertices)] do new_labels[i]:=RandomList(new_vertices); RemoveSet(new_vertices,new_labels[i]); od; new_facets:=[]; for facet in facets do new_facet:=[]; for k in facet do AddSet(new_facet,new_labels[Position(vertices,k)]); od; AddSet(new_facets,new_facet); od; facets:=ShallowCopy(new_facets); fi; ############################################################### # Computing the dimension dim=max-1 of the simplicial complex # ############################################################### max_start:=0; for facet in facets do if Length(facet) > max_start then max_start:=Length(facet); fi; od; ####################### # Computing all faces # ####################### faces:=[]; for i in [1..max_start] do faces[i]:=[]; od; for i in [1..max_start] do for facet in facets do if Length(facet) >= i then UniteSet(faces[i],Combinations(facet,i)); fi; od; od; if round_B = 1 then F:=[]; Print("f-vector: ("); for i in [1..(max_start-1)] do F[i]:=Length(faces[i]); Print(F[i],","); od; F[max_start]:=Length(faces[max_start]); Print(F[max_start],")\n"); fi; ########################### # Computing Hasse diagram # ########################### initial_hasse:=[]; downward_hasse:=[]; for i in [1..(max_start-1)] do initial_hasse[i]:=[]; downward_hasse[i+1]:=[]; for a in [1..Length(faces[i])] do initial_hasse[i][a]:=[]; od; for b in [1..Length(faces[i+1])] do downward_hasse[i+1][b]:=[]; sub_faces:=Combinations(faces[i+1][b],i); for face in sub_faces do AddSet(initial_hasse[i][Position(faces[i],face)],b); AddSet(downward_hasse[i+1][b],Position(faces[i],face)); od; od; od; ## ## ## ## ## ## ## ## ## ## ## ## ## ###################################### ### ## Computing critical cells ## ### ###################################### ## ## ## ## ## ## ## ## ## ## ## ## ## if round_B = 1 then Print(StringTime(Runtime()),"\n"); fi; for round_A in [1..rounds_A] do max:=max_start; hasse:=[]; for i in [1..(max_start-1)] do hasse[i]:=[]; for a in [1..Length(initial_hasse[i])] do hasse[i][a]:=ShallowCopy(initial_hasse[i][a]); od; od; current_faces:=[]; for i in [1..max_start] do current_faces[i]:=[1..F[i]]; od; removed_vertices:=[]; pairs:=[]; if flag_manifold = 0 then # initialize pairs for a in [1..Length(hasse[max-1])] do if Length(hasse[max-1][a]) = 1 then AddSet(pairs,[a,hasse[max-1][a][1]]); fi; od; fi; ######################### # Sequence of collapses # ######################### morse_vector:=[]; for k in [1..max_start] do morse_vector[k]:=0; od; while max > 1 do if Length(pairs) > 0 then if strategy = 1 then random_pair:=RandomList(pairs); elif strategy = 2 then random_pair:=pairs[1]; elif strategy = 3 then random_pair:=pairs[Length(pairs)]; fi; RemoveSet(current_faces[max-1],random_pair[1]); RemoveSet(current_faces[max],random_pair[2]); for t in downward_hasse[max][random_pair[2]] do RemoveSet(hasse[max-1][t],random_pair[2]); if Length(hasse[max-1][t]) = 0 then RemoveSet(pairs,[t,random_pair[2]]); elif Length(hasse[max-1][t]) = 1 then AddSet(pairs,[t,hasse[max-1][t][1]]); fi; od; if max > 2 then for s in downward_hasse[max-1][random_pair[1]] do RemoveSet(hasse[max-2][s],random_pair[1]); od; elif max = 2 then AddSet(removed_vertices,faces[1][random_pair[1]]); fi; RemoveSet(pairs,random_pair); else if current_faces[max] = [] then max:=max-1; if max > 1 then # initialize pairs for a in [1..Length(hasse[max-1])] do if Length(hasse[max-1][a]) = 1 then AddSet(pairs,[a,hasse[max-1][a][1]]); fi; od; fi; else if strategy = 1 then random_index:=RandomList(current_faces[max]); elif strategy = 2 then random_index:=current_faces[max][1]; elif strategy = 3 then random_index:=current_faces[max][Length(current_faces[max])]; fi; morse_vector[max]:=morse_vector[max]+1; RemoveSet(current_faces[max],random_index); if max > 1 then for t in downward_hasse[max][random_index] do RemoveSet(hasse[max-1][t],random_index); if Length(hasse[max-1][t]) = 1 then AddSet(pairs,[t,hasse[max-1][t][1]]); fi; od; fi; fi; fi; od; remaining_vertices:=ShallowCopy(faces[1]); SubtractSet(remaining_vertices,removed_vertices); morse_vector[1]:=Length(remaining_vertices); #Print(round," Critical cells: ",morse_vector,"\n"); if morse_vector in morse_vector_collection then morse_vector_count[Position(morse_vector_collection,morse_vector)]:=morse_vector_count[Position(morse_vector_collection,morse_vector)]+1; else Add(morse_vector_collection,morse_vector); morse_vector_count[Length(morse_vector_collection)]:=1; fi; od; od; for k in [1..Length(morse_vector_collection)] do Print(morse_vector_collection[k],": ",morse_vector_count[k],"\n"); od; plain:=0; normalized:=0; for k in [1..Length(morse_vector_collection)] do sum:=0; for j in [1..Length(morse_vector_collection[k])] do sum:=sum+morse_vector_collection[k][j]; od; plain:=plain+morse_vector_count[k]*sum; normalized:=normalized+morse_vector_count[k]*(sum-2*(morse_vector_collection[k][1]-1)); od; plain:=plain/rounds; normalized:=normalized/rounds; Print("average number of critical cells (plain) = ",plain,"\n"); Print("average number of critical cells (normalized) = ",normalized,"\n"); Print(StringTime(Runtime()),"\n"); #QUIT;