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
   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 = lnm
  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     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         l = range(w)
 103         l.reverse()
 104         for j in l:   # pour chaque colonne de H_hat
 105             y[indx] = lit(path,(indx,state))
 106             state = xor(state,y[indx]*H_hat[j],h)
 107             indx -=1
 108         state = 2*state + message[indm]
 109         indm -=1
 110         i +=1
 111     return [int(t) for t in y]
 112
 113
 114
 115
 116
 117
 118 def trouve_H_hat(n,m,h):
 119     assert h ==7
 120     alpha = float(n)/m
 121     assert alpha >= 2
 122     index = min(int(alpha),9)
 123     mat = {
 124         2 : [71,109],
 125         3 : [95, 101, 121],
 126         4 : [81, 95, 107, 121],
 127         5 : [75, 95, 97, 105, 117],
 128         6 : [73, 83, 95, 103, 109, 123],
 129         7 : [69, 77, 93, 107, 111, 115, 121],
 130         8 : [69, 79, 81, 89, 93, 99, 107, 119],
 131         9 : [69, 79, 81, 89, 93, 99, 107, 119, 125]
 132         }[index]
 133     return(mat, index*m)
 134
 135
 136 def stc(x,rho,message):
 137     lnm = len(message)
 138     (mat,taille_suff) = trouve_H_hat(len(x),len(message),7)
 139     x_b = x[:taille_suff]
 140     (start,path) = forward(mat,x_b,message,lnm,rho)
 141     return (x_b,backward(start,mat,x_b,message,lnm,path),mat)
 142
 143
 144
 145
 146
 147 def nbdif(x,y):
 148     r,it = 0,0
 149     l = len(y)
 150     while it < l :
 151         if x[it] != y[it] :
 152             r +=1
 153         it += 1
 154     return float(r)/l
 155
 156
 157
 158
 159
 160
 161 def prod(H_hat,lnm,y):
 162     (h,w) = int(log(max(H_hat),2))+1, len(H_hat)
 163     i=0
 164     H =[]
 165     V=[0 for _ in range(len(y))]
 166     sol=[]
 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 xrange(min([i+1,h])): #nbre de blocks presents sur la ligne i
 172             for l in xrange(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         sol.append(np.dot(np.array(V),np.array(y)))
 177         i+=1
 178         #H += [V]
 179     #H = np.array(H)
 180     #y = np.array(y)
 181     #print "dot",np.dot(H,y),H.shape
 182     #print "sol",sol
 183     return sol#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 x = [randint(0,1) for _ in xrange(65000)]
 199 rho = [randint(1,9) for _ in xrange(65000)]
 200 message = [randint(0,1) for _ in xrange(26000)]
 201
 202
 203
 204 #(x_b,y,H_hat) = stc(x,rho,message)
 205 # x_b est la sous partie de x qui va etre modifiee
 206 # y est le vecteur des bits modifies
 207 # H_hat est la sous-matrice retenue qui est embarquee dans H
 208
 209 #print nbdif(x_b,y)
 210
 211 #print message
 212 #print x
 213 #print rho
 214 #print y
 215 #message2 = [x%2 for x in prod(H_hat,len(message),y)]
 216 #print message2
 217 print max([message[l]-message2[l] for l in xrange(len(message))])
 218 print equiv(message,message2)
 219
 220 """
 221 count = 0
 222
 223 while count < 100 and eval :
 224     lx = randint(500,1000)
 225     x = [randint(0,1) for _ in xrange(lx)]
 226     rho = [randint(1,9) for _ in xrange(lx)]
 227
 228     lm = randint(lx/10,lx/2)
 229     message = [randint(0,1) for _ in xrange(lm)]
 230     (x_b,y,H_hat) = stc(x,rho,message)
 231     eval = equiv(message, prod(H_hat,len(message),y))
 232     if not (eval):
 233         print x
 234         print message
 235
 236     count +=1
 237     if count % 10 == 0 :
 238         print count
 239
 240 print count
 241