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         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
  71                 w1 = wght[xor(k,H_hat[j],h)] + (1-x[indx])*rho
  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     start = np.argmin(wght)
  89     return (start,path)
  90
  91
  92 def backward(start,H_hat,x,message,lnm,path):
  93     (h,w) = int(log(max(H_hat),2))+1, len(H_hat)
  94     indx,indm = len(x)-1,lnm-1
  95     state = 2*start + message[indm]
  96     indm -=1
  97     # l'initialisation de state n'est pas optimale...
  98     nbblock = lnm
  99     y=np.zeros(len(x))
 100     i=0
 101     while i < nbblock:
 102         if i%10000 == 0 :
 103             print i
 104         l = range(w)
 105         l.reverse()
 106         for j in l:   # pour chaque colonne de H_hat
 107             y[indx] = lit(path,(indx,state))
 108             state = xor(state,y[indx]*H_hat[j],h)
 109             indx -=1
 110         state = 2*state + message[indm]
 111         indm -=1
 112         i +=1
 113     return [int(t) for t in y]
 114
 115
 116
 117
 118
 119
 120 def trouve_H_hat(n,m,h):
 121     assert h ==7
 122     alpha = float(n)/m
 123     assert alpha >= 2
 124     index = min(int(alpha),9)
 125     mat = {
 126         2 : [71,109],
 127         3 : [95, 101, 121],
 128         4 : [81, 95, 107, 121],
 129         5 : [75, 95, 97, 105, 117],
 130         6 : [73, 83, 95, 103, 109, 123],
 131         7 : [69, 77, 93, 107, 111, 115, 121],
 132         8 : [69, 79, 81, 89, 93, 99, 107, 119],
 133         9 : [69, 79, 81, 89, 93, 99, 107, 119, 125]
 134         }[index]
 135     return(mat, index*m)
 136
 137
 138 def stc(x,message):
 139     lnm = len(message)
 140     (mat,taille_suff) = trouve_H_hat(len(x),len(message),7)
 141     x_b = x[:taille_suff]
 142     (start,path) = forward(mat,x_b,message,lnm)
 143     return (x_b,backward(start,mat,x_b,message,lnm,path),mat)
 144
 145
 146
 147
 148
 149 def nbdif(x,y):
 150     r,it = 0,0
 151     l = len(y)
 152     while it < l :
 153         if x[it] != y[it] :
 154             r +=1
 155         it += 1
 156     return float(r)/l
 157
 158
 159
 160
 161
 162
 163 def prod(H_hat,lnm,y):
 164     (h,w) = int(log(max(H_hat),2))+1, len(H_hat)
 165     i=0
 166     H =[]
 167     while i < lnm: # pour chaque ligne
 168         V=[0 for _ in range(len(y))]
 169         k = max([(i-h+1)*w,0])
 170         dec = max([i-h+1,0])
 171         for j in range(min([i+1,h])): #nbre de blocks presents sur la ligne i
 172             for l in range(w): # pour chaque collone de H_hat
 173                 V[k] = bin(H_hat[l],h)[h-i-1+j+dec]
 174                 k+=1
 175
 176         i+=1
 177         H += [V]
 178     H = np.array(H)
 179     y = np.array(y)
 180     return list(np.dot(H,y))
 181
 182 def equiv(x,y):
 183     lx = len(x)
 184     assert lx == len(y)
 185     i=0
 186     while i < lx :
 187         if x[i] % 2 != y[i]%2 :
 188             return False
 189         i += 1
 190     return True
 191
 192
 193 ################
 194
 195 x = [randint(0,1) for _ in xrange(50000)]
 196 message = [randint(0,1) for _ in xrange(10000)]
 197
 198
 199 (x_b,y,H_hat) = stc(x,message)
 200 # x_b est la sous partie de x qui va etre modifiee
 201 # y est le vecteur des bits modifies
 202 # H_hat est la sous-matrice retenue qui est embarquee dans H
 203
 204
 205 print nbdif(x_b,y)
 206
 207
 208
 209
 210