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