stc/exp/python/stc_r.py

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