]> AND Private Git Repository - canny.git/blob - stc/exp/python/stc_lght.py
Logo AND Algorithmique Numérique Distribuée

Private GIT Repository
Relecture
[canny.git] / stc / exp / python / stc_lght.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(5000)]
192 message = [randint(0,1) for _ in xrange(2500)]
193
194
195
196 (x_b,y,H_hat) = stc(x,message)
197 # x_b est la sous partie de x qui va etre modifiée
198 # y est le vecteur des bits modifiés
199 # H_hat est la sous-matrice retenue qui est embarquee dans H  
200
201 print nbdif(x_b,y)
202
203
204
205
206