rooted-forests.sws

def apply(L,fn):
    res = []
    for elt in L:
        res = res + [fn(elt)]
    return res


# Input: Simplicial complex K (d-dimensional)
# Output: All spanning forests (in the sense of [BK15])


def SpanningForests(K): #gives all the SPANNING forests of a simplicial complex K
    zero_hom_grp = SimplicialComplex([0]).homology()[0]
    d = dim(K)
    faceList = list(K.faces()[d])
    co_K = K.chain_complex()
    third_cond = len(list(K.faces()[d])) - co_K.betti(d)
    possible_forests = [x for x in list(powerset(faceList)) if len(x) == third_cond]
    res = []
    for j in range(len(possible_forests)):
        facet_list = possible_forests[j]
        F_test = SimplicialComplex(facet_list)
        if F_test.homology()[dim(K)] == zero_hom_grp:
            res = res + [facet_list]
    return res

# Input: Triple (F,K,d) of a simplicial complex K dimension n, subset F of d-cells of K, dimension 0<=d<=n
# Output: Boolean True/False whether F is a forest of K

# Description: Checks whether the corresponding columns of \partial_d^F are independent

def is_forest(F,K,d):
    B = K.chain_complex().differential(d)
    kept_rows = apply(F,list(K.faces()[d]).index)
    corresponding_columns = [vector(QQ,list(B.columns()[n])) for n in kept_rows]
    V = QQ^(len(B.rows()))
    if len(V.linear_dependence(corresponding_columns)) == 0:
        return True
    else:
        return False

# Input: Pair (R,G) where G is a d-dimensional simplicial complex (forest -- not checked?) and R subset of (d-1)-cells
# Output: Boolean True/False on whether R is indeed a root of G

# Note: Does not check whether G is indeed a forest to save comnputational resources. Can be done by adding an is_forest check.

def is_root(R,G):
    d = dim(G)
    R_bar = list(G.faces()[d-1])
    for cell in R:
        R_bar.remove(cell)
    B = G.chain_complex().differential(d)
    kept_rows = apply(R_bar,list(G.faces()[d-1]).index)
    keyB_rows = []
    for index in kept_rows:
        keyB_rows = keyB_rows + [B.rows()[index]]
    keyB = Matrix(keyB_rows)
    if len(keyB.rows()) != len(keyB.columns()):
        return False
    elif keyB.is_invertible():  #techically transpose of keyB, but this is an equivalent condition
        return True

# Input: Pair (K,d) of simplicial complex K (of dimension n) and dimension 0<=d<=n
# Output: List of all d-dimensional forests F of K

# Note: Could use some optimization if intended to be used for large complexes.

def all_forests(K,d):
    faceList = K.faces()[d]
    possible_forests = list(powerset(faceList))[1:] 
    res =[]
    for F in possible_forests:
        if is_forest(F,K,d):
            res = res + [F]
    return res
    


#TODO: Fix below. Proper roots aren't getting returned. (Error with is_root function?)
#returns all roots R (list of (d-1)-faces) of a forest F (subset of d-1 faces of complex G)
def all_roots(F):
    F = SimplicialComplex(F)
    d = dim(F)
    faceList = F.faces()[d-1]
    possible_roots = list(powerset(faceList))[1:-1]
    res = []
    for R in possible_roots:
        if is_root(R,F):
            res = res + [R]
    return res


#K.homology(1,subcomplex=S) #relative homology!!















def rootedForests(K,d):
    forestList = K.faces()[d]
    rootList = K.faces()[d-1]
    possible_forests = list(powerset(forestList))[1:]
    possible_roots = list(powerset(rootList))[1:]
    D = K.chain_complex().differential(d)
    Ld = list(K.faces()[d])
    Ld.sort()
    L_d = list(K.faces()[d-1])
    L_d.sort()
    res = []
    for F in possible_forests:
        for R in possible_roots:
            R_bar = list(K.faces()[d-1])
            for cell in R:
                R_bar.remove(cell)
            if len(R_bar) == len(F): #|R_bar| = |F|
                cols = []
                for i in range(len(F)):
                    cols = cols + [Ld.index(F[i])]
                rows = []
                for i in range(len(R_bar)):
                    rows = rows + [L_d.index(R_bar[i])]
                rows.sort()
                cols.sort()
                if D[rows,cols].det() != 0: #relative homology H_d(F,R) = 0
                    F_res = apply(F,tuple)
                    R_res = apply(R,tuple)
                    F_res.sort()
                    R_res.sort()
                    res = res + [(F_res,R_res)] #apply function turns these from simplex cells to tuples in F and R, types in Sage
    return res


def rootedForestsCount(K,d):
    S = rootedForests(K,d)
    L = [0..len(K.faces()[d-1])]
    for j in range(len(L)):
        L[j] = len([x for x in S if len(x[1]) == j])
    return L



#[x for x in S[5] if not(x in G[5])]






def rootedForestsList(K,d):
    S = rootedForests(K,d)
    L = [0..len(K.faces()[d-1])]
    for j in range(len(L)):
        L[j] = [x for x in S if len(x[1]) == j]
    L.reverse()
    return L


#using grad_list_raw(K,d)[n] gives all the gradients of K corresponding to particular coeff. of Laplacian (explicitly note which eventually, TODO)
def grad_list_raw(K,d):
    vec_list = morse_set_vectors(K,d)
    #pool = limited_acyclicmatchings(K)
    grads = []
    pool = limited_acyclicmatchings2(K,d)
    for k in range(len(vec_list)):
        vec = vec_list[k]
        res = strainer(K,d,pool,vec)
        grads = grads + [res]
    return grads



#TODO.
#returns all the rooted d-forests coming from gradients, in the format of (F,R,G), where G is the gradient -- use fn[n] for induced forests with |F| = n
def inducedRootedForests(K,d,remove_grads = False):
    rawList = grad_list_raw(K,d)
    rootCells = apply(list(K.faces()[d-1]),tuple)
    res = []
    for n in range(len(rawList)):
        res_temp = []
        gradient_list = rawList[n]
        for grad in gradient_list:
            induced_forest = []
            for m in range(len(grad)):
                induced_forest = induced_forest + [grad[m][0]]
            induced_root = list(rootCells) #this is where it gets complicated, need to look at elts not in gradient
            for m in range(len(grad)):
                induced_root.remove(grad[m][1])
            induced_forest.sort()
            induced_root.sort()
            if remove_grads == True:
                res_temp = res_temp + [(induced_forest,induced_root)]
            elif remove_grads == False:
                res_temp = res_temp + [(induced_forest,induced_root,grad)] #where to decide whether or not to put grad as third argument
        res = res + [res_temp]
    return res


#input rooted d-forest (F,R), returns \mathcal{H}_{(F,R)}, optional third argument when True provides image
def rootedForest_hasse(F,R,im = False):
    d = len(F[0])-1
    K = SimplicialComplex(F)
    H = hasse_levels(K,[d-1,d])
    H_elts = []
    for n in range(d-1,d+1):
        H_elts = H_elts + apply(list(K.faces()[n]),list)
    R = apply(R,list)
    for cell in R:
        if cell in H_elts:
            H_elts.remove(cell)
    H_elts = apply(H_elts,tuple)
    res = H.subposet(H_elts)
    if im == False:
        return res
    if im == True:
        return plot(res.hasse_diagram().reverse())

M = SimplicialComplex([[0,1,4],[0,4,3],[4,3,2],[3,2,1],[2,1,0]])




#use checker(K,d) first, see if "False", hence having a problematic forest
def problems(K,d):
    L = inducedRootedForests(K,d,True)
    induced_list = []
    for n in range(len(L)):
        induced_list = induced_list + L[n]
    rootedforest_list = rootedForests(K,d)
    res = [x for x in rootedforest_list if x not in induced_list]
    return res


#K = SimplicialComplex([[1,2,3,4]])
#d = 2

#S = all_forests(K,d)
#S = [F for F in S if len(F) == 2]
#res = 0
#for forest in S:
#    res = res + len(all_roots(forest))
#res


#code for list of forests after setting L = proj_nice list
#sage: for n in range(len(L)):
#....:     print(str(n))
#....:     rootedForest_hasse(L[n][0],L[n][1],True)
#....:


#for n in range(len(L)): rootedForest_hasse(L[n][0],L[n][1],True).save('/Users/atawfeek/Desktop/rooted_forests/' + str(n) + '.png')


#to save outputted images in sage:
# p = G.plot() #use this instead of show()
# p.save('/Users/atawfeek/Desktop/test.png')