experiments/genHamiltonian.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 import random as rd
   7 from copy import *
   8 from combinaisons import *
   9
  10
  11 # Implemantation of the balanced gray code algorithm,construction B
  12 # detailled in "Totally Balanced and exponentially Balanced Gray codes",
  13 # A J. van Zanten, I.N. Suparta 2004
  14
  15
  16
  17 def bin(elem,n):
  18     """Convertit un nombre en binaire"""
  19     q = -1
  20     res =[0 for i in range(n)]
  21     i = 1
  22     while q != 0:
  23         q = elem // 2
  24         r = elem % 2
  25         res[n-i] =  r
  26         elem = q
  27         i+=1
  28     return res
  29
  30 def dec(ch,n):
  31     l = len(ch)
  32     acc = 0
  33     for i in range(l):
  34         if ch[i]==1:
  35             acc = acc + 2**(n-i-1)
  36     return acc
  37
  38
  39 def flatten(x):
  40     result = []
  41     for el in x:
  42         if hasattr(el, "__iter__") and not isinstance(el, basestring):
  43             result.extend(flatten(el))
  44         else:
  45             result.append(el)
  46     return result
  47
  48
  49
  50 def genereToutesLesSousListes(n,su):
  51     # specificique puisque la taille mini est 2 et commence par [0,1]
  52     print "n",n,"su",su
  53     r= [[0,1]]
  54     for j in xrange(su-2):
  55         print j
  56         rp = []
  57         for li in r :
  58             for el in range(li[-1]+1,n):
  59                 rp += [ li+[el]]
  60         r= [x for x in rp]
  61     return r
  62
  63
  64
  65
  66
  67 def extraitSousListe(Snm2,sublisIndex,n):
  68     #l is sublistdindices
  69     nm2 = len(Snm2)
  70     indexOfSnm2 = xrange(nm2)
  71     # l is an even positive number
  72     #assert l >= 0 and l%2 == 0 and l<= nm2
  73     l = len(sublisIndex)
  74     #print "l",l,"n",n,"len(Snm2)",nm2
  75     # we have to select l elements, but the two first are always selected
  76     #Silindex = [0,1]+rd.sample(xrange(2,nm2),l-2)
  77     Silindex = sublisIndex
  78     Silindex.sort()
  79     Sil = [Snm2[x] for x in Silindex]
  80
  81     dif = list(set(indexOfSnm2)-set(Silindex))
  82     dif.sort()
  83
  84     Uinitial=[[]]
  85     UinitialIndex = [[]]
  86     cur = 2
  87
  88     for k in xrange(2,l):
  89         Uinitial += [[Snm2[x] for x in range(cur,Silindex[k])]]
  90         UinitialIndex = [[x for x in range(cur,Silindex[k])]]
  91         cur = Silindex[k]+1
  92     v = [Snm2[x] for x in range(cur,nm2)]
  93
  94     #print "Snm2",Snm2,"Silindex",Silindex,"Sil",Sil,"dif",dif,"Uinitial",Uinitial,"UinitialIndex",UinitialIndex,"v",v
  95     return (Silindex,Sil,Uinitial,v)
  96
  97
  98
  99
 100 def replace(Silindex,Sil,u,n,l,v):
 101     def uExtended(u,j,n):
 102         # If j is odd, we are in u(n-1,n), and u(n,n-1) otherwise
 103         (X,Y) = (n-1,n) if i%2 == 1 else (n,n-1)
 104
 105         uir = [_ for _ in u[i]]
 106         uir.reverse()
 107         R = [_ for _ in u[i]]+ [X] + uir + [Y] + [_ for _ in u[i]]
 108         return R
 109
 110     lu = len(u)
 111
 112     Up = [n-1]+[uExtended(u,i,n) for i in range(1,lu)]
 113
 114     U=[]
 115     for j in xrange(l-2+1):
 116         U += [Sil[j]]
 117         U += [Up[j]]
 118     U += [Sil[-1]]
 119     U += v
 120     return flatten(U)
 121
 122
 123 def generateGrayCode(Snm2,sublistdindices):
 124
 125     #Snm2 = [3,2,3,1,3,2,3,1]
 126     n= int(log(len(Snm2))/log(2))+2
 127     l=int(pow(2,n)/n)
 128     #step 1
 129     (Silindex,Sil,Uinitial,v) = extraitSousListe(Snm2,sublistdindices,n)
 130
 131
 132     #step 2 : replace in Snm2 and genereate U
 133     U = replace(Silindex,Sil,Uinitial,n,l,v)
 134
 135     #step 3 : compute V and W
 136     V = [_ for _ in v]
 137     V.reverse()
 138     V += [n] + [_ for _ in v]
 139     W = [n-1] + [_ for _ in Snm2] +[n]
 140
 141     #step 4 : swap value n-1 and s1 in W
 142     Wp = [W[1]] + [n-1] + [_ for _ in W[2:]]
 143
 144     #step 5 : generate Sn
 145     Ur = [_ for _ in U]
 146     Ur.reverse()
 147
 148     Sn = Ur + V + Wp
 149     return Sn
 150
 151
 152 def is_totallyBalanced(l):
 153     res = [0 for _ in range(max(l))]
 154     for j in l :
 155         res[j-1] +=1
 156     return (max(res) == min(res))
 157
 158 def is_onlyBalanced(l):
 159     res = [0 for _ in range(max(l))]
 160     for j in l :
 161         res[j-1] +=1
 162     return max(res)- min(res) < 3
 163
 164
 165 def grapheToList(G,n):
 166     l2n = range(int(pow(2,n)))
 167     edges = G.edges()
 168     images = range(int(pow(2,n)))
 169     for o in l2n:
 170         e = o
 171         eb= bin(e,n)
 172         for j in range(n):
 173             epb = bin(e,n)
 174             epb[j] = 1-epb[j]
 175             if (o,dec(epb,n)) in edges:
 176                 eb[j] = 1-eb[j]
 177         images[o]=dec(eb,n)
 178     return images
 179
 180
 181 def Sn2map(l):
 182     n= int(log(len(l))/log(2))
 183     init = bin(0,n)
 184     cheminb=[]
 185     a = copy(init)
 186     for i in l:
 187         a[n-i] = 1 - a[n-i]
 188         cheminb += [copy(a)]
 189     chemin=[dec(x,n) for x in cheminb]
 190     lc = len(chemin)
 191     enl=[(chemin[i%lc],chemin[(i+1)%lc]) for i in range(lc)]
 192     assert len(set(range(int(pow(2,n))))-set(chemin)) == 0
 193     # cheminp est le chemin hamiltonien correspondant a l
 194     # on construit la fonction sans ce chemin hamiltonien
 195     G = nx.DiGraph()
 196     for i in range(2**n):
 197         ib = bin(i,n)
 198         for k in range(n):
 199             ibp = [_ for _ in ib]
 200             ibp[k]=1-ibp[k]
 201             # si l'arrete ne doit pas etre enlevee
 202             if (i,dec(ibp,n)) not in enl :
 203                 G.add_edge(i,dec(ibp,n))
 204
 205     #resl = grapheToList(G,n)
 206     return (G,grapheToList(G,n))
 207
 208
 209
 210
 211
 212
 213 def main():
 214     snm2 = [[1,2,1,2]]
 215     #snm2 = [[4, 5, 6, 2, 6, 5, 1, 5, 6, 2, 6, 5, 3, 5, 6, 1, 6, 5, 4, 5, 6, 1, 3, 2, 3, 4, 1, 4, 6, 4, 1, 4, 3, 2, 3, 5, 3, 2, 3, 4, 1, 4, 3, 5, 2, 6, 2, 5, 3, 4, 1, 4, 3, 2, 3, 1, 4, 1, 3, 2, 1, 2, 4, 6]]
 216     resG=[]
 217     p = 3
 218     for k in range(p):
 219         print "taille de snm2", len(snm2)
 220         snm2p = []
 221         for sn in snm2 :
 222             n= int(log(len(sn))/log(2))+2
 223             ln= int(pow(2,n)/n)
 224             #subls =  genereToutesLesSousListes(int(pow(2,n-2)),ln)
 225             #print "subls",subls
 226             ll = range(int(pow(2,n-2)))
 227             #for l in subls:
 228             compteur = 0
 229             for l in xuniqueCombinations(ll[2:], ln-2):
 230                 compteur +=1
 231                 l = [0,1]+l
 232                 nouveauCode =  generateGrayCode(sn,l)
 233
 234
 235                 flg = False
 236
 237                 (G,nc) = Sn2map(nouveauCode)
 238                 #on verifie que la fonction est nouvelle
 239                 if nx.is_strongly_connected(G):
 240                     if len(resG)==0:
 241                         resG.append(G)
 242                         flg=True
 243                     else :
 244                         def isomorphic(g1,g2):
 245                             GM = isomorphism.DiGraphMatcher(g1,g2)
 246                             return GM.is_isomorphic()
 247                         if all([not(isomorphic(G,g2)) for g2 in resG]):
 248                             resG.append(G)
 249                             flg=True
 250                     if is_onlyBalanced(nouveauCode) :
 251                         snm2p +=[nouveauCode]
 252                         if is_totallyBalanced(nouveauCode) :
 253                             print "TB", nc
 254                         else :
 255                             print "OB", nc
 256         snm2=deepcopy(snm2p)
 257
 258
 259
 260
 261
 262
 263
 264 def options():
 265     global rf
 266     parser = OptionParser()
 267     parser.add_option("-i", "--input", dest="i",
 268                       help="file of sequences")
 269     (options, args) = parser.parse_args()
 270     if (options.i != None):
 271         rf = options.i
 272
 273
 274 if __name__ == "__main__":
 275     options()
 276     print "------"
 277     main()
 278
 279