load("Cubulation.sage")
load("Script_iso_sig_standard.py")

#####

# Script in Sage to find the 63 non-equivalent states for the 24-cell whose ascending and descending links all collapse to points

# Find all orbits of states. We can supppose that the first three digits are zero.
# May take roughly one hour

states = get_states(C24(), C24_colouring(), k=3)
print ("States found: ", len(states))

# It has found 2632 states

G = C24()
colouring = C24_colouring()

# Select only those states where all the 24 cusps fiber

good_states = []
for state in states:
    if all_cusps_fiber(G, 4, colouring, state):
        good_states.append(state[:])
print ("Sates with fibering cusps: ", len(good_states))        

# It finds 148 of them. Now select only those whose ascending and descending links collapse to points or circles.

very_good_states = []
for state in good_states:
    print (state)
    if links_collapse(G, colouring, state, 1, True, True) == True:
        very_good_states.append(state[:])
print ("States found: ", len(very_good_states))

# It finds the 63 states we are interested in
# Now for each state we determine the singular fiber.

import snappy, regina

manifolds = []
counter = 1
for state in very_good_states:
    print ("State", counter, " = ", state)
    cubes = get_cubulation(C24(), 4, C24_colouring(), True, state, extended = True)
    diag = get_diagonal_complex(cubes)
    print ("Peel")
    diag_peeled = peel(diag)
    top_peeled = []
    for simplex in diag_peeled:
        if len(simplex) == 4:
            top_peeled.append(simplex[:])
    print ("Calculate the iso_signature")
    iso_sigs = get_iso_signature(top_peeled)
    iso_sig = iso_sigs[0]
    print ("Send to Regina")
    M = regina.Triangulation3(iso_sig)
    print ("Convert finite vertices to ideal")
    M.finiteToIdeal()
    print ("Simplify")
    M.simplifyToLocalMinimum()
    print ("We get a triangulation with ", M.countTetrahedra(), "tetrahedra")
    new_iso_sig = M.isoSig()
    print ("Send it to SnapPy")
    N = snappy.Manifold(new_iso_sig)
    manifolds.append(N.copy())
    print ("Volume = ", N.volume())
    counter = counter + 1
    
# The list "manifolds" contains the 63 manifolds in the SnapPy format to be investigated

# For each manifold we write a string 
# [volume, symmetry group, num tetrahedra in canonical decomposition, state]

manifolds_data = []
counter = 0
for state in very_good_states:
    M = manifolds[counter]
    M.simplify()
    for tentative in range(3):
        M.randomize()
        M.simplify()
    M.canonize()
    new_data = [M.volume(), M.symmetry_group(), M.num_tetrahedra(), state]
    print (counter, new_data)
    counter = counter + 1
    manifolds_data.append(new_data[:])
    
##########

# Proof that the 120-cell has a state that fulfills all the requirement to furnish a perfect Morse function

G = C120()  # May take some minutes
colouring = C120_colouring()
state = C120_state()
links_collapse(G, colouring, state, 1, True, False)

# It returns "True", which means that all links collapse to connected 1-dimensional complexes