stc/exp/python/stc.py

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