Ticket #12418: delsarte_bound.py

File delsarte_bound.py, 6.7 KB (added by dimpase, 9 years ago)

a prototype implementation

Line 
1def Krawtchouk(n,q,l,i): 
2   """
3   Compute ``K^{n,q}_l(i)``, the Krawtchouk polynomial:
4   see ``http://en.wikipedia.org/wiki/Kravchuk_polynomials``
5
6   """
7
8   return sum([((-1)**j)*((q-1)**(l-j))*binomial(i,j)*binomial(n-i,l-j) 
9                              for j in range(l+1)])
10
11def delsarte_bound_hamming_space(n, q, d, \
12                    isinteger=False, return_data=False, solver="PPL"):
13   """
14   Find the classical Delsarte bound on codes in Hamming space
15   ``H_q^n`` of minimal distance ``d``
16
17
18   INPUT:
19
20   - ``n`` -- the code length
21
22   - ``q`` -- the length of the alphabet
23
24   - ``d`` -- the (lower bound on) minimal distance of the code
25
26   - ``isinteger`` -- if ``True``, uses an integer programming solver (ILP), rather
27     that an LP solver. Can be very slow if set to ``True``.
28
29   - ``return_data`` -- if ``True``, return a weights vector, which actually need not
30     be a proper weight enumerator, or even have integer entries, and the LP.
31
32   - ``solver`` -- the LP/ILP solved to be used. Defaults to ``PPL``. It is arbitrary
33     precision, thus there will be no rounding errors. With other solvers, you are on
34     your own!
35
36
37   EXAMPLES:
38
39   The bound on the size of the F_2-codes of length 11 and minimal distance 6::     
40 
41       sage: delsarte_bound_hamming_space(11, 2, 6)
42       12 
43       sage: a, p, val = delsarte_bound_hamming_space(11, 2, 6)
44       sage: [j for i,j in p.get_values(a).iteritems()]
45       [1, 0, 0, 0, 0, 0, 11, 0, 0, 0, 0, 0]
46       
47   The bound on the size of the F_2-codes of length 24 and minimal distance
48   8, i.e. parameters of the extened binary Golay code::     
49
50       sage: a,p,x=delsarte_bound_hamming_space(24,2,8,return_data=True)
51       sage: x
52       4096
53       sage: [j for i,j in p.get_values(a).iteritems()]
54       [1, 0, 0, 0, 0, 0, 0, 0, 759, 0, 0, 0, 2576, 0, 0, 0, 759, 0, 0, 0, 0, 0, 0, 0, 1]
55
56   The bound on the size of F_4-codes of length 11 and minimal distance 3::
57       
58       sage: delsarte_bound_hamming_space(11,4,3)
59       327680/3
60
61   """
62
63   p = MixedIntegerLinearProgram(maximization=True, solver=solver)
64   A = p.new_variable(integer=isinteger) # A>=0 is assumed
65   p.set_objective(sum([A[r] for r in range(n+1)])) 
66   p.add_constraint(A[0]==1)
67   for i in range(1,d):
68      p.add_constraint(A[i]==0)
69   for j in range(1,n+1): 
70      rhs = sum([Krawtchouk(n,q,j,r)*A[r] for r in range(n+1)])
71      p.add_constraint(0*A[0] <= rhs) 
72   try:
73      bd=p.solve()
74   except sage.numerical.mip.MIPSolverException, exc:
75      print "Solver exception: ", exc, exc.args
76      if return_data:
77         return A,p,False
78      return False
79   
80   if return_data:
81      return A,p,bd
82   else: 
83      return bd
84
85def delsarte_bound_additive_hamming_space(n, q, d, d_star=1, q_base=0, 
86                    isinteger=False, return_data=False, solver="PPL"):
87   """
88   Find the Delsarte LP bound on ``F_{q_base}``-dimension of additive codes in
89   Hamming space ``H_q^n`` of minimal distance ``d`` with minimal distance of the dual
90   code at least ``d_star``.  If ``q_base`` is set to
91   non-zero, then  ``q`` is a power of ``q_base``, and the code is, formally, linear over
92   ``F_{q_base}``. Otherwise it is assumed that ``q_base==q``.
93
94
95   INPUT:
96
97   - ``n`` -- the code length
98
99   - ``q`` -- the length of the alphabet
100
101   - ``d`` -- the (lower bound on) minimal distance of the code
102
103   - ``d_star`` -- the (lower bound on) minimal distance of the dual code;
104     only makes sense for additive codes.
105
106   - ``q_base`` -- if ``0``, the code is assumed to be nonlinear. Otherwise,
107     ``q=q_base^m`` and the code is linear over ``F_{q_base}``.
108
109   - ``isinteger`` -- if ``True``, uses an integer programming solver (ILP), rather
110     that an LP solver. Can be very slow if set to ``True``.
111
112   - ``return_data`` -- if ``True``, return a weights vector, which actually need not
113     be a proper weight enumerator, or even have integer entries, and the LP.
114
115   - ``solver`` -- the LP/ILP solved to be used. Defaults to ``PPL``. It is arbitrary
116     precision, thus there will be no rounding errors. With other solvers, you are on
117     your own!
118
119
120   EXAMPLES:
121
122   The bound on dimension of linear F_2-codes of length 11 and minimal distance 6::     
123 
124       sage: delsarte_bound_additive_hamming_space(11, 2, 6)
125       3
126       sage: a,p,val=delsarte_bound_additive_hamming_space(11, 2, 6,\
127                                      return_data=True)
128       sage: [j for i,j in p.get_values(a).iteritems()]
129       [1, 0, 0, 0, 0, 0, 5, 2, 0, 0, 0, 0]
130
131   The bound on the dimension of linear F_4-codes of length 11 and minimal distance 3::
132       
133       sage: delsarte_bound_additive_hamming_space(11,4,3)
134       8
135
136   The bound on the F_2-dimension of additive F_4-codes of length 11 and minimal
137   distance 3::
138     
139       sage: delsarte_bound_additive_hamming_space(11,4,3,q_base=2)
140       16
141         
142   """
143
144   if q_base == 0:
145      q_base = q
146
147   kk = 0
148   while q_base**kk < q:
149      kk += 1
150
151   if q_base**kk != q:
152      print "Wrong q_base=", q_base, " for q=", q, kk
153      return False
154
155   # this implementation assumes that our LP solver to be unable to do a hot
156   # restart with an adjusted constraint
157
158   m = kk*n # this is to emulate repeat/until block
159   bd = q**n+1 
160
161   while q_base**m < bd: # need to solve the LP repeatedly, as this is a new constraint!
162                         # we might become infeasible. More precisely, after rounding down
163                         # to the closest value of q_base^m, the LP, with the constraint that
164                         # the objective function is at most q_base^m,
165      p = MixedIntegerLinearProgram(maximization=True, solver=solver)
166      A = p.new_variable(integer=isinteger) # A>=0 is assumed
167      p.set_objective(sum([A[r] for r in xrange(n+1)])) 
168      # the following constraint comes from rounding
169      p.add_constraint(sum([A[r] for r in xrange(n+1)]) <= q_base**m) 
170      p.add_constraint(A[0]==1)
171      for i in xrange(1,d):
172        p.add_constraint(A[i]==0)
173      for j in xrange(1,n+1): 
174        rhs = sum([Krawtchouk(n,q,j,r)*A[r] for r in xrange(n+1)])
175        if j >= d_star: 
176          p.add_constraint(0*A[0] <= rhs) 
177        else: # rhs is proportional to j-th weight of the dual code
178          p.add_constraint(0*A[0] == rhs) 
179      try:
180        bd=p.solve()
181      except sage.numerical.mip.MIPSolverException, exc:
182        print "Solver exception: ", exc, exc.args
183        if return_data:
184           return A,p,False
185        return False
186    # rounding the bound down to the nearest power of q_base, for q=q_base^m
187#      bd_r = roundres(log(bd, base=q_base))
188      m = -1
189      while q_base**(m+1) < bd:
190        m += 1
191      if q_base**(m+1) == bd:
192        m += 1
193         
194   if return_data: 
195      return A, p, m
196   else:
197      return m
198