experiments/genDoubleStoc.py

   1 import networkx as nx
   2 from networkx.algorithms import isomorphism
   3 from math import *
   4 import numpy as np
   5 from optparse import OptionParser
   6
   7
   8 n= 5
   9 nbmaxiter = 6
  10 profM=n*pow(2,n)
  11 prof = profM
  12
  13 def bin(elem,n):
  14     """Convertit un nombre en binaire"""
  15     q = -1
  16     res =[0 for i in range(n)]
  17     i = 1
  18     while q != 0:
  19         q = elem // 2
  20         r = elem % 2
  21         res[n-i] =  r
  22         elem = q
  23         i+=1
  24     return res
  25
  26 def dec(ch,n):
  27     l = len(ch)
  28     acc = 0
  29     for i in range(l):
  30         if ch[i]==1:
  31             acc = acc + 2**(n-i-1)
  32     return acc
  33
  34
  35 def calculetmix(epsilon):
  36     resTMIX=[]
  37     global treatedlist
  38     global n
  39     if rf != False:
  40         treatedlist =get_treatedlist(rf)
  41     #print "taillle de la liste ",len(treatedlist)
  42     n= int(log(len(treatedlist[0]))/log(2))
  43     print "n",n
  44     l2n = range(int(pow(2,n)))
  45     pi0 = [pow(2,-n) for _ in range(int(pow(2,n)))]
  46     #treatedlist = [treatedlist[11]]
  47     cpt = 0
  48     for perm in treatedlist:
  49         M = np.zeros((pow(2,n),pow(2,n)))
  50         for conf in l2n:
  51             r = bin(perm[conf],n)
  52             for cp in range(n):
  53                 confbp = bin(conf,n)
  54                 confbp[cp] =  r[cp]
  55                 M[conf-1][dec(confbp,n)-1] = 1/float(n)
  56         mx=0
  57         #print "Markov chain",M
  58         eigs = np.linalg.eigvals(M)
  59         #print "valp",eigs
  60         error = 1E-10
  61         modulusofegs = [abs(l) for l in eigs]
  62         modulusofegs = [l if (1-l)>error else 0 for l in modulusofegs]
  63         #print "modulusofegs",modulusofegs
  64         lambdastar = max(modulusofegs)
  65         #print "lambdastar",lambdastar
  66         trel = 1.0/(1-lambdastar)
  67         def tmixSup(epsilon) :
  68             pimin = 1.0/int(pow(2,n))
  69             #print "pimin",pimin
  70             return np.log(1.0/(epsilon*pimin))*trel
  71         def tmixInf(epsilon) :
  72             pimin = 1.0/int(pow(2,n))
  73             #print "pimin",pimin
  74             return (trel-1)*np.log(1.0/(2*epsilon))
  75
  76
  77         print "F_"+str(cpt+1), perm, int(tmixInf(epsilon)), int(tmixSup(epsilon))+1
  78         resTMIX +=[(perm, int(tmixInf(epsilon)), int(tmixSup(epsilon))+1)]
  79         cpt+=1
  80
  81     def compare(x,y):
  82         return cmp(x[1],y[1])
  83     resTMIX.sort(cmp=compare)
  84     print resTMIX
  85
  86
  87
  88 def is_double_stoc(M):
  89     (nbl,nbc) =M.shape
  90     return all( [abs (sum([M[i,j] for i in range(nbl)])-1.0)<1E-6 for j in range(nbc)])
  91
  92
  93
  94 def matriceAdjacenceDe(G):
  95     l2n = range(int(pow(2,n)))
  96     P = np.zeros((2**n,2**n))
  97     L = G.edges()
  98     for i in l2n:
  99         k = 0
 100         for j in l2n:
 101             if i != j :
 102                 if (i,j) in L:
 103                     P[i][j] = 1/float(n)
 104                     k+=1
 105         P[i][i]=(n-k)/float(n)
 106     return P
 107
 108
 109 def graphe_OK(G):
 110     M = matriceAdjacenceDe(G)
 111     r = True
 112     r = r and is_double_stoc(M)
 113     r = r and nx.is_strongly_connected(G)
 114     return r
 115
 116
 117 def graphe_de_liste(l):
 118     G = nx.DiGraph()
 119     l2n = range(int(pow(2,n)))
 120     for conf in l2n:
 121         G.add_node(conf)
 122         # r est l'image binaire de conf en parallele
 123         r = bin(l[conf],n)
 124         # on construit poentiellement les n image en unitaire
 125         for cp in range(n):
 126             # confbp est conf sauf en cp ou c'est r[cp]
 127             confbp = bin(conf,n)
 128             confbp[cp] =  r[cp]
 129             # on ajoute l'arc
 130             G.add_edge(conf,dec(confbp,n))
 131     return G
 132
 133
 134 acc={}
 135 def genereGrapheHamiltonien(G):
 136     global acc
 137     def largeur(profondeur, effectues):
 138         #print "profondeur",profondeur
 139         global acc
 140         effectues_p = []
 141         for (c,vi) in effectues:
 142             (_,d) = c[len(c)-1]
 143             for v in G.neighbors(d):
 144                 if v not in vi:
 145                     effectues_p += [([_ for _ in c]+[(d,v)],set([_ for _ in vi]+[v]))]
 146         acc[profondeur] = effectues_p
 147         if profondeur <= pow(2,n)  :
 148             largeur(profondeur+1,effectues_p)
 149     largeur(4,[([(0,1)],set([0,1]))])
 150     res = acc[pow(2,n)+1]
 151     ret=[]
 152     for (c,vertices) in res :
 153         assert len(vertices) == pow(2,n)
 154         (_,t) = c[len(c)-1]
 155         if G.has_edge(t,0) :
 156             ret +=[c +[(t,0)]]
 157
 158     return ret
 159
 160
 161 def graphSansListe(G,l):
 162     Gp = nx.DiGraph()
 163     Gp.add_edges_from(G.edges())
 164     Gp.remove_edges_from(l)
 165     return Gp
 166
 167 def genereTousLesGraphesSansCycleHamiltonien(G):
 168     r = [graphSansListe(G,l) for l in genereGrapheHamiltonien(G)]
 169     res = [r[0]]
 170     for g in r[1:] :
 171         if all([not(isomorphic(g,g2)) for g2 in res]):
 172             res.append(g)
 173     return [grapheToList(G) for G in res]
 174
 175 def graphe_de_edges(l):
 176     G = nx.DiGraph()
 177     l2n = range(int(pow(2,n)))
 178     for conf in l2n:
 179         G.add_node(conf)
 180     for (o,e) in l:
 181         G.add_edge(o,e)
 182     return G
 183
 184 # s'il y a des doublons ce n'est pas une perm
 185 def isperm(L):
 186     return len(set(L)) == len(L)
 187
 188
 189 def grapheToList(G):
 190     l2n = range(int(pow(2,n)))
 191     edges = G.edges()
 192     images = range(int(pow(2,n)))
 193     for o in l2n:
 194         e = o
 195         eb= bin(e,n)
 196         for j in range(n):
 197             epb = bin(e,n)
 198             epb[j] = 1-epb[j]
 199             if (o,dec(epb,n)) in edges:
 200                 eb[j] = 1-eb[j]
 201         images[o]=dec(eb,n)
 202     return images
 203
 204 def isomorphic(g1,g2):
 205     GM = isomorphism.DiGraphMatcher(g1,g2)
 206     return GM.is_isomorphic()
 207
 208 def parcours_p(visited,resG,l):
 209
 210     t = len(l)
 211     for j in range(t):
 212         lp = l[0:j]+l[j+1:]
 213         G = graphe_de_edges(lp)
 214         if nx.is_strongly_connected(G) :
 215             if all([not(isomorphic(G,g2)) for g2 in visited]):
 216                 visited.append(G)
 217                 M = matriceAdjacenceDe(G)
 218                 if is_double_stoc(M):
 219                     resG.append(G)
 220                     gtl = grapheToList(G)
 221                     print gtl,len(resG)
 222                 parcours_p(visited,resG,G.edges())
 223
 224
 225 def parcours_l_dual(visited,resG,GGedges):
 226     global prof
 227     print "\n en entrant",int(prof),"a visiter",len(visited)
 228     prof -=1
 229     nvisited =[]
 230     for i in range(len(visited)) :
 231         #if i%10 == 0 :
 232         print "ds calcul",i,len(nvisited)
 233         ld = visited[i].edges()
 234         lr = list(GGedges-set(ld))
 235
 236         t = len(lr)
 237         for j in range(t):
 238             lp = [x for x in ld]+[lr[j]]
 239             G= graphe_de_edges(list(GGedges-set(lp)))
 240             if nx.is_strongly_connected(G) :
 241                 Gd= graphe_de_edges(lp)
 242                 if all([not(isomorphic(Gd,g2)) for g2 in nvisited]):
 243                     nvisited.append(Gd)
 244                     M = matriceAdjacenceDe(G)
 245                     if is_double_stoc(M):
 246                         resG.append(G)
 247                         gtl = grapheToList(G)
 248                         print "\n",gtl
 249     if nvisited != [] :
 250         parcours_l_dual(nvisited,resG,GGedges)
 251
 252
 253 def parcours_l(visited,resG):
 254     global prof
 255     if prof < profM - nbmaxiter  :
 256         assert False
 257     print "\n en entrant",int(prof),"a vistier",len(visited)
 258     prof -=1
 259     nvisited =[]
 260     for i in range(len(visited)) :
 261         #if i%100 == 0 :
 262         #print "ds calcul",i,len(nvisited)
 263         l = visited[i].edges()
 264         t = len(l)
 265         for j in range(t):
 266             lp = l[0:j]+l[j+1:]
 267             G = graphe_de_edges(lp)
 268             if nx.is_strongly_connected(G) :
 269                 if all([not(isomorphic(G,g2)) for g2 in nvisited]):
 270                     nvisited.append(G)
 271                     M = matriceAdjacenceDe(G)
 272                     if is_double_stoc(M):
 273                         resG.append(G)
 274                         gtl = grapheToList(G)
 275                         print "\n",gtl
 276     if nvisited != [] :
 277         parcours_l(nvisited,resG)
 278
 279
 280
 281
 282
 283 def main():
 284     global treatedlist
 285     l = range(int(pow(2,n)))
 286     l.reverse()
 287     G = graphe_de_liste(l)
 288     treatedlist = genereTousLesGraphesSansCycleHamiltonien(G)
 289     print treatedlist
 290     #Ginit = nx.DiGraph()
 291     #Ginit.add_edge(0,2)
 292     #resG = [G]
 293     # parcours en profondeur : deux graphes sont voinsin si l'un contient
 294     # un arc de plus que l'autre
 295     #visited=[]
 296     #parcours_p(visited,resG,G.edges())
 297     #parcours_l([G],resG)
 298     #parcours_l_dual([Ginit],resG,set(G.edges()))
 299
 300     #print [grapheToList(X) for X in  resG]
 301
 302
 303
 304
 305
 306 def main_old():
 307     resl=[]
 308     global treatedlist
 309     if rf != False:
 310         treatedlist =get_treatedlist(rf)
 311     else :
 312         treatedlist = [ perm for perm in permute_in_place(range(8))]
 313     print "taillle de la liste ",len(treatedlist)
 314     res = []
 315     resG= []
 316     global n
 317     n= int(log(len(treatedlist[0]))/log(2))
 318     print "n",n
 319
 320     nbit = 40
 321
 322     l2n = range(int(pow(2,n)))
 323     # pour chaque permutation de l ensemble [0,...,2^[n-1}}
 324     nnb = 0
 325     #for perm in permute_in_place(range(4)):
 326     D=[0 for x in range(len(treatedlist))]
 327     count=-1
 328     for perm in treatedlist:
 329       count +=1
 330       nnb +=1
 331       #if nnb%1000 == 0 :
 332       #  print nnb
 333       # si elle contient un element stable on l enleve (GTPIC n est pas SCC)
 334       if all([ ident!=perm[ident] for ident in l2n]):
 335         # on construit le GTPIC
 336         #G = nx.DiGraph()
 337         G = nx.DiGraph()
 338         # pour chaque noeud
 339         for conf in l2n:
 340             G.add_node(conf)
 341             # r est l'image binaire de conf en parallele
 342             r = bin(perm[conf])
 343             # on construit poentiellement les n image en unitaire
 344             for cp in range(n):
 345                 # confbp est conf sauf en cp ou c'est r[cp]
 346                 confbp = bin(conf)
 347                 confbp[cp] =  r[cp]
 348                 # on ajoute l'arc
 349                 G.add_edge(conf,dec(confbp))
 350         flg = False
 351         if nx.is_strongly_connected(G):
 352             if len(resG)==0:
 353                 res.append(perm)
 354                 resG.append(G)
 355                 flg=True
 356             else :
 357                 def isomorphic(g1,g2):
 358                     GM = isomorphism.DiGraphMatcher(g1,g2)
 359                     return GM.is_isomorphic()
 360                 if all([not(isomorphic(G,g2)) for g2 in resG]):
 361                     resG.append(G)
 362                     res.append(perm)
 363                     flg=True
 364         #else:
 365         #    print "pas scc",perm
 366         if flg == True :
 367             T=[0 for i in l2n]
 368             T[0] = 1
 369             P = np.zeros((2**n,2**n))
 370             for i in l2n:
 371                 k = 0
 372                 for j in l2n:
 373                     L = G.edges()
 374                     if i != j :
 375                         if (i,j) in L:
 376                             P[i][j] = 1/float(n)
 377                             k+=1
 378                 P[i][i]=(n-k)/float(n)
 379
 380
 381             if is_double_stoc(P):
 382                 d=pow(2,n)
 383                 t=0
 384                 unif = 1/float(pow(2,n))
 385                 cpt=0
 386                 while (d >0.00001 and t <10000):
 387                     cpt+=1
 388                     #print cpt
 389                     T =[ sum ([T[j]*P[j][i] for j in l2n]) for i in l2n]
 390                     #print T
 391                     dp = sum([abs(x-unif) for x in T])*(pow(2,n))
 392                     if dp < d :
 393                         D[count] = t
 394                         d = dp
 395                     t +=1
 396                 #print dp
 397                 print perm, t, d, D[count]
 398                 resl +=[perm]
 399     print resl
 400     return resl
 401 rf=False
 402
 403 def options():
 404     global rf
 405     parser = OptionParser()
 406     parser.add_option("-i", "--input", dest="i",
 407                       help="file of sequences")
 408     (options, args) = parser.parse_args()
 409     if (options.i != None):
 410         rf = options.i
 411
 412
 413 if __name__ == "__main__":
 414     options()
 415     print "------"
 416     main()
 417     calculetmix(1E-4)
 418