stc/exp/python/stc_raphus.py

   1 from random import *
   2 import numpy as np
   3 from math import *
   4 import gc
   5
   6 infinity = 100000000
   7
   8
   9
  10 # forward part
  11
  12
  13 def dec(ch,n):
  14     l = len(ch)
  15     acc = 0
  16     for i in xrange(l):
  17         if ch[i]==1:
  18             acc = acc + 2**(n-i-1)
  19     return acc
  20
  21
  22 def bin(elem,n):
  23     """Convertit un nombre en binaire"""
  24     q = -1
  25     res = [0 for i in xrange(n)]
  26     i = 1
  27     while q != 0:
  28         q = elem // 2
  29         r = elem % 2
  30         res[n-i] =  r
  31         elem = q
  32         i+=1
  33     return res
  34
  35
  36
  37 def xorb(a,b):
  38     return 1 if a != b else 0
  39
  40 def xor(e1,e2,h):
  41     e1b,e2b  = bin(e1,h),bin(e2,h)
  42     d = dec([xorb(e1b[j],e2b[j]) for j in xrange(h)],h)
  43     return d
  44
  45 def lit(d,(indx,indy)):
  46     if (indx,indy) in d :
  47         return d[(indx,indy)]
  48     else :
  49         return 0
  50
  51
  52
  53
  54 def forward(H_hat,x,message,lnm,rho):
  55     (h,w) = int(log(max(H_hat),2))+1, len(H_hat)
  56     path = dict()
  57     nbblock = len(x)/w
  58     wght = [infinity for _ in xrange(int(2**h))]
  59     wght[0]=0
  60     newwght = [0 for _ in xrange(int(2**h))]
  61 #    rho = 1
  62 #    rho= [1 for _ in xrange(len(x))]
  63     indx,indm = 0,0
  64     i=0
  65     while i < nbblock: # pour chaque bit du message
  66         for j in xrange(w):   # pour chaque colonne de H_hat
  67             print indx, "en entrant",wght
  68             k = 0
  69             while k < int(2**h): # pour chaque ligne de H
  70                 w0 = wght[k] + x[indx]*rho[indx]
  71                 w1 = wght[xor(k,H_hat[j],h)] + (1-x[indx])*rho[indx]
  72                 if w1 < w0 :
  73                     path[(indx,k)] = 1
  74                 else :
  75                     if (indx,k) in path:
  76                         del path[(indx,k)]
  77                 newwght[k] = min(w0,w1)
  78                 k +=1
  79             indx +=1
  80             wght = [t for t in newwght]
  81             print " apres calcul",wght
  82
  83         for j in xrange(int(2**(h-1))):   # pour chaque colonne de H
  84             wght[j] = wght[2*j + message[indm]]
  85         wght = wght[:int(pow(2,h-1))] + [infinity for _ in xrange(int(pow(2,h)-pow(2,h-1)))]
  86         indm +=1
  87         i +=1
  88     # juste sur le modulo
  89     reste_a_faire = len(x) % w
  90     if reste_a_faire != 0 :
  91         for j in xrange(reste_a_faire):   # pour chaque colonne de H_hat qui reste
  92             #print indx, "en entrant",wght
  93             k = 0
  94             while k < int(2**h): # pour chaque ligne de H
  95                 w0 = wght[k] + x[indx]*rho[indx]
  96                 w1 = wght[xor(k,H_hat[j],h)] + (1-x[indx])*rho[indx]
  97                 if w1 < w0 :
  98                     path[(indx,k)] = 1
  99                 else :
 100                     if (indx,k) in path:
 101                         del path[(indx,k)]
 102                 newwght[k] = min(w0,w1)
 103                 k +=1
 104             indx +=1
 105             wght = [t for t in newwght]
 106             #print " apres calcul",wght
 107
 108         for j in xrange(int(2**(h-1))):   # pour chaque colonne de H
 109             wght[j] = wght[2*j + message[indm]]
 110         wght = wght[:int(pow(2,h-1))] + [infinity for _ in xrange(int(pow(2,h)-pow(2,h-1)))]
 111         indm +=1
 112     #fin du modulo
 113
 114     start = np.argmin(wght)
 115     return (start,path)
 116
 117
 118 def backward(start,H_hat,x,message,lnm,path):
 119     (h,w) = int(log(max(H_hat),2))+1, len(H_hat)
 120     indx,indm = len(x)-1,lnm-1
 121     state = 2*start + message[indm]
 122     indm -=1
 123     # l'initialisation de state n'est pas optimale...
 124     nbblock = len(x)/w
 125     y=np.zeros(len(x))
 126     i=0
 127     # sur le reste
 128     reste_a_faire = len(x) % w
 129     if reste_a_faire != 0 :
 130         l = range(reste_a_faire)
 131         l.reverse()
 132         for j in l:   # pour chaque colonne qui reste a faire
 133             y[indx] = lit(path,(indx,state))
 134             state = xor(state,y[indx]*H_hat[j],h)
 135             indx -=1
 136         state = 2*state + message[indm]
 137         indm -=1
 138     # fin du reste
 139     while i < nbblock:
 140         l = range(w)
 141         l.reverse()
 142         for j in l:   # pour chaque colonne de H_hat
 143             y[indx] = lit(path,(indx,state))
 144             state = xor(state,y[indx]*H_hat[j],h)
 145             indx -=1
 146         state = 2*state + message[indm]
 147         indm -=1
 148         i +=1
 149     return [int(t) for t in y]
 150
 151
 152
 153
 154
 155
 156 def trouve_H_hat(n,m,h):
 157     assert h ==7
 158     alpha = float(n)/m
 159     assert alpha >= 1
 160     index = min(int(alpha),9)
 161     matr = {
 162         2 : [71,109],
 163         3 : [95, 101, 121],
 164         4 : [81, 95, 107, 121],
 165         5 : [75, 95, 97, 105, 117],
 166         6 : [73, 83, 95, 103, 109, 123],
 167         7 : [69, 77, 93, 107, 111, 115, 121],
 168         8 : [69, 79, 81, 89, 93, 99, 107, 119],
 169         9 : [69, 79, 81, 89, 93, 99, 107, 119, 125]
 170         }
 171     mat = []
 172     if index not in matr:
 173          mat = matr[2]
 174     else :
 175         mat = matr[index]
 176     return mat
 177
 178
 179 def stc(x,rho,message):
 180     lnm = len(message)
 181     mat = trouve_H_hat(len(x),len(message),7)
 182     x_b = [i for i in x]
 183     mat=[3,1]
 184     (start,path) = forward(mat,x_b,message,lnm,rho)
 185     return (x_b,backward(start,mat,x_b,message,lnm,path),mat)
 186
 187
 188
 189
 190
 191 def nbdif(x,y):
 192     r,it = 0,0
 193     l = len(y)
 194     while it < l :
 195         if x[it] != y[it] :
 196             r +=1
 197         it += 1
 198     return float(r)/l
 199
 200
 201
 202
 203
 204
 205 def prod(H_hat,lnm,y):
 206     (h,w) = int(log(max(H_hat),2))+1, len(H_hat)
 207     i=0
 208     H =[]
 209     V=[0 for _ in range(len(y))]
 210     sol=[]
 211     Vp =[]
 212     while i < lnm: # pour chaque ligne
 213         V=[0 for _ in range(len(y))]
 214         k = max([(i-h+1)*w,0])
 215         dec = max([i-h+1,0])
 216         for j in xrange(min([i+1,h])): #nbre de blocks presents sur la ligne i
 217             for l in xrange(w): # pour chaque collone de H_hat
 218                 if k < len(y):
 219                     V[k] = bin(H_hat[l],h)[h-i-1+j+dec]
 220                     k+=1
 221         sol.append(np.dot(np.array(V),np.array(y)))
 222         i+=1
 223         Vp +=[V]
 224     Vp =  np.array(Vp)
 225     print Vp
 226     return np.dot(Vp,np.array(y))
 227
 228
 229 def equiv(x,y):
 230     lx = len(x)
 231     assert lx == len(y)
 232     i=0
 233     while i < lx :
 234         if x[i] % 2 != y[i]%2 :
 235             return False
 236         i += 1
 237     return True
 238
 239
 240 ################
 241 """
 242 x = [randint(0,1) for _ in xrange(65000)]
 243 rho = [randint(1,9) for _ in xrange(65000)]
 244 message = [randint(0,1) for _ in xrange(26000)]
 245 """
 246 x = [0, 0, 1, 0, 0, 1]
 247 rho = [1 for _ in xrange(100)]
 248 message = [0, 1, 0, 1]
 249
 250
 251
 252 (x_b,y,H_hat) = stc(x,rho,message)
 253 print "message", message
 254 print "y", y
 255 print "x", x
 256 print "H_hat", H_hat
 257
 258
 259
 260 # x_b est la sous partie de x qui va etre modifiee
 261 # y est le vecteur des bits modifies
 262 # H_hat est la sous-matrice retenue qui est embarquee dans H
 263
 264 print "avec stc :", nbdif(x_b,y)
 265 print "sans stc :", nbdif(message,x[:len(message)])
 266
 267 #print message
 268 #print x
 269 #print rho
 270 #print y
 271 message2 = [x%2 for x in prod(H_hat,len(message),y)]
 272 print "messag2", message2
 273
 274 #print message2
 275 print equiv(message,message2)
 276