import networkx as nx
from networkx.algorithms import isomorphism
from math import *
import numpy as np 
from optparse import OptionParser


n= 5
nbmaxiter = 6
profM=n*pow(2,n)
prof = profM
 
def bin(elem,n):
    """Convertit un nombre en binaire"""
    q = -1
    res =[0 for i in range(n)]
    i = 1
    while q != 0:
        q = elem // 2
        r = elem % 2
        res[n-i] =  r
        elem = q
        i+=1
    return res

def dec(ch,n):
    l = len(ch)
    acc = 0
    for i in range(l):
        if ch[i]==1:
            acc = acc + 2**(n-i-1)        
    return acc


def calculetmix(epsilon):
    resTMIX=[]
    global treatedlist
    global n 
    if rf != False:
        treatedlist =get_treatedlist(rf)
    #print "taillle de la liste ",len(treatedlist)
    n= int(log(len(treatedlist[0]))/log(2))
    print "n",n
    l2n = range(int(pow(2,n)))
    pi0 = [pow(2,-n) for _ in range(int(pow(2,n)))]
    #treatedlist = [treatedlist[11]]
    cpt = 0
    for perm in treatedlist:
        M = np.zeros((pow(2,n),pow(2,n)))
        for conf in l2n:
            r = bin(perm[conf],n)
            for cp in range(n):
                confbp = bin(conf,n)
                confbp[cp] =  r[cp]
                M[conf-1][dec(confbp,n)-1] = 1/float(n)
        mx=0
        #print "Markov chain",M
        eigs = np.linalg.eigvals(M)
        #print "valp",eigs
        error = 1E-10
        modulusofegs = [abs(l) for l in eigs]
        modulusofegs = [l if (1-l)>error else 0 for l in modulusofegs]
        #print "modulusofegs",modulusofegs
        lambdastar = max(modulusofegs)
        #print "lambdastar",lambdastar
        trel = 1.0/(1-lambdastar)
        def tmixSup(epsilon) :
            pimin = 1.0/int(pow(2,n))
            #print "pimin",pimin 
            return np.log(1.0/(epsilon*pimin))*trel
        def tmixInf(epsilon) :
            pimin = 1.0/int(pow(2,n))
            #print "pimin",pimin 
            return (trel-1)*np.log(1.0/(2*epsilon))


        print "F_"+str(cpt+1), perm, int(tmixInf(epsilon)), int(tmixSup(epsilon))+1
        resTMIX +=[(perm, int(tmixInf(epsilon)), int(tmixSup(epsilon))+1)]
        cpt+=1

    def compare(x,y):
        return cmp(x[1],y[1])
    resTMIX.sort(cmp=compare)
    print resTMIX
    


def is_double_stoc(M):
    (nbl,nbc) =M.shape
    return all( [abs (sum([M[i,j] for i in range(nbl)])-1.0)<1E-6 for j in range(nbc)])
    


def matriceAdjacenceDe(G):
    l2n = range(int(pow(2,n)))
    P = np.zeros((2**n,2**n))
    L = G.edges()
    for i in l2n:
        k = 0
        for j in l2n:
            if i != j : 
                if (i,j) in L:
                    P[i][j] = 1/float(n)
                    k+=1
        P[i][i]=(n-k)/float(n)
    return P
        

def graphe_OK(G):
    M = matriceAdjacenceDe(G)
    r = True 
    r = r and is_double_stoc(M)
    r = r and nx.is_strongly_connected(G)
    return r


def graphe_de_liste(l):
    G = nx.DiGraph()
    l2n = range(int(pow(2,n)))
    for conf in l2n:
        G.add_node(conf)
        # r est l'image binaire de conf en parallele
        r = bin(l[conf],n)
        # on construit poentiellement les n image en unitaire
        for cp in range(n):
            # confbp est conf sauf en cp ou c'est r[cp] 
            confbp = bin(conf,n)
            confbp[cp] =  r[cp]
            # on ajoute l'arc
            G.add_edge(conf,dec(confbp,n))
    return G


acc={}
def genereGrapheHamiltonien(G):
    global acc
    def largeur(profondeur, effectues):
        #print "profondeur",profondeur
        global acc
        effectues_p = []
        for (c,vi) in effectues: 
            (_,d) = c[len(c)-1]
            for v in G.neighbors(d):
                if v not in vi:
                    effectues_p += [([_ for _ in c]+[(d,v)],set([_ for _ in vi]+[v]))]
        acc[profondeur] = effectues_p 
        if profondeur <= pow(2,n)  : 
            largeur(profondeur+1,effectues_p)
    largeur(4,[([(0,1)],set([0,1]))])
    res = acc[pow(2,n)+1]
    ret=[]
    for (c,vertices) in res :
        assert len(vertices) == pow(2,n)
        (_,t) = c[len(c)-1]
        if G.has_edge(t,0) : 
            ret +=[c +[(t,0)]]
    
    return ret


def graphSansListe(G,l):
    Gp = nx.DiGraph()
    Gp.add_edges_from(G.edges())
    Gp.remove_edges_from(l)
    return Gp
    
def genereTousLesGraphesSansCycleHamiltonien(G):
    r = [graphSansListe(G,l) for l in genereGrapheHamiltonien(G)]
    res = [r[0]]
    for g in r[1:] :
        if all([not(isomorphic(g,g2)) for g2 in res]):
            res.append(g)
    return [grapheToList(G) for G in res]
        
def graphe_de_edges(l):
    G = nx.DiGraph()
    l2n = range(int(pow(2,n)))
    for conf in l2n:
        G.add_node(conf)
    for (o,e) in l:
        G.add_edge(o,e)
    return G

# s'il y a des doublons ce n'est pas une perm
def isperm(L):
    return len(set(L)) == len(L)


def grapheToList(G):
    l2n = range(int(pow(2,n)))
    edges = G.edges()
    images = range(int(pow(2,n)))
    for o in l2n:
        e = o
        eb= bin(e,n)
        for j in range(n):
            epb = bin(e,n)
            epb[j] = 1-epb[j]
            if (o,dec(epb,n)) in edges:
                eb[j] = 1-eb[j]
        images[o]=dec(eb,n)
    return images
    
def isomorphic(g1,g2):
    GM = isomorphism.DiGraphMatcher(g1,g2)
    return GM.is_isomorphic()

def parcours_p(visited,resG,l):
    
    t = len(l)
    for j in range(t):
        lp = l[0:j]+l[j+1:]
        G = graphe_de_edges(lp)
        if nx.is_strongly_connected(G) :
            if all([not(isomorphic(G,g2)) for g2 in visited]):
                visited.append(G)
                M = matriceAdjacenceDe(G)
                if is_double_stoc(M):
                    resG.append(G)
                    gtl = grapheToList(G)
                    print gtl,len(resG)
                parcours_p(visited,resG,G.edges())


def parcours_l_dual(visited,resG,GGedges):
    global prof
    print "\n en entrant",int(prof),"a visiter",len(visited)
    prof -=1
    nvisited =[]
    for i in range(len(visited)) :
        #if i%10 == 0 :
        print "ds calcul",i,len(nvisited)
        ld = visited[i].edges()
        lr = list(GGedges-set(ld))

        t = len(lr)
        for j in range(t):
            lp = [x for x in ld]+[lr[j]]
            G= graphe_de_edges(list(GGedges-set(lp)))             
            if nx.is_strongly_connected(G) :
                Gd= graphe_de_edges(lp)
                if all([not(isomorphic(Gd,g2)) for g2 in nvisited]):
                    nvisited.append(Gd)
                    M = matriceAdjacenceDe(G)
                    if is_double_stoc(M):
                        resG.append(G)
                        gtl = grapheToList(G)
                        print "\n",gtl
    if nvisited != [] :
        parcours_l_dual(nvisited,resG,GGedges)


def parcours_l(visited,resG):
    global prof
    if prof < profM - nbmaxiter  :
        assert False
    print "\n en entrant",int(prof),"a vistier",len(visited)
    prof -=1
    nvisited =[]
    for i in range(len(visited)) :
        #if i%100 == 0 :
        #print "ds calcul",i,len(nvisited)
        l = visited[i].edges()
        t = len(l)
        for j in range(t):
            lp = l[0:j]+l[j+1:]
            G = graphe_de_edges(lp)
            if nx.is_strongly_connected(G) :
                if all([not(isomorphic(G,g2)) for g2 in nvisited]):
                    nvisited.append(G)
                    M = matriceAdjacenceDe(G)
                    if is_double_stoc(M):
                        resG.append(G)
                        gtl = grapheToList(G)
                        print "\n",gtl
    if nvisited != [] :
        parcours_l(nvisited,resG)



            

def main():
    global treatedlist
    l = range(int(pow(2,n)))
    l.reverse()
    G = graphe_de_liste(l)
    treatedlist = genereTousLesGraphesSansCycleHamiltonien(G)
    print treatedlist
    #Ginit = nx.DiGraph()
    #Ginit.add_edge(0,2)
    #resG = [G]
    # parcours en profondeur : deux graphes sont voinsin si l'un contient
    # un arc de plus que l'autre
    #visited=[]
    #parcours_p(visited,resG,G.edges())
    #parcours_l([G],resG)
    #parcours_l_dual([Ginit],resG,set(G.edges()))

    #print [grapheToList(X) for X in  resG]
    




def main_old():
    resl=[]
    global treatedlist
    if rf != False:
        treatedlist =get_treatedlist(rf)
    else : 
        treatedlist = [ perm for perm in permute_in_place(range(8))]
    print "taillle de la liste ",len(treatedlist)
    res = []
    resG= []
    global n 
    n= int(log(len(treatedlist[0]))/log(2))
    print "n",n

    nbit = 40

    l2n = range(int(pow(2,n)))
    # pour chaque permutation de l ensemble [0,...,2^[n-1}}
    nnb = 0
    #for perm in permute_in_place(range(4)):
    D=[0 for x in range(len(treatedlist))]
    count=-1
    for perm in treatedlist:
      count +=1
      nnb +=1
      #if nnb%1000 == 0 :
      #  print nnb
      # si elle contient un element stable on l enleve (GTPIC n est pas SCC)
      if all([ ident!=perm[ident] for ident in l2n]):
        # on construit le GTPIC  
        #G = nx.DiGraph()
        G = nx.DiGraph()
        # pour chaque noeud 
        for conf in l2n:
            G.add_node(conf)
            # r est l'image binaire de conf en parallele
            r = bin(perm[conf])
            # on construit poentiellement les n image en unitaire
            for cp in range(n):
                # confbp est conf sauf en cp ou c'est r[cp] 
                confbp = bin(conf)
                confbp[cp] =  r[cp]
                # on ajoute l'arc
                G.add_edge(conf,dec(confbp))
        flg = False
        if nx.is_strongly_connected(G):
            if len(resG)==0:
                res.append(perm)
                resG.append(G)
                flg=True
            else :
                def isomorphic(g1,g2):
                    GM = isomorphism.DiGraphMatcher(g1,g2)
                    return GM.is_isomorphic()
                if all([not(isomorphic(G,g2)) for g2 in resG]):
                    resG.append(G)
                    res.append(perm)
                    flg=True
        #else:
        #    print "pas scc",perm
        if flg == True :
            T=[0 for i in l2n]
            T[0] = 1
            P = np.zeros((2**n,2**n))
            for i in l2n:
                k = 0
                for j in l2n:
                    L = G.edges()
                    if i != j : 
                        if (i,j) in L:
                            P[i][j] = 1/float(n)
                            k+=1
                P[i][i]=(n-k)/float(n)
                
        
            if is_double_stoc(P):    
                d=pow(2,n)
                t=0
                unif = 1/float(pow(2,n))
                cpt=0
                while (d >0.00001 and t <10000):
                    cpt+=1
                    #print cpt
                    T =[ sum ([T[j]*P[j][i] for j in l2n]) for i in l2n]
                    #print T
                    dp = sum([abs(x-unif) for x in T])*(pow(2,n))
                    if dp < d :
                        D[count] = t
                        d = dp
                    t +=1 
                #print dp
                print perm, t, d, D[count]
                resl +=[perm]
    print resl
    return resl
rf=False

def options():
    global rf
    parser = OptionParser()
    parser.add_option("-i", "--input", dest="i",
                      help="file of sequences")
    (options, args) = parser.parse_args()
    if (options.i != None):
        rf = options.i


if __name__ == "__main__":
    options()
    print "------"
    main()
    calculetmix(1E-4)