Opened 10 years ago

Last modified 8 years ago

#12418 closed enhancement

adding Delsarte bound for codes — at Version 4

Reported by: dimpase Owned by: wdj
Priority: major Milestone: sage-5.12
Component: coding theory Keywords:
Cc: jpang, kini, wdj, ncohen, ppurka, ptrrsn_1 Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Status badges

Description (last modified by dimpase)

Delsarte bound for codes, aka Linear Programming bound, is easy to implement in Sage. Here is a quick and dirty code that does it; so the problem would be to integrate this properly.

def Kra(n,q,l,i): # K^{n,q}_l(i), i.e Krawtchouk polynomial
   return sum([((-1)**j)*((q-1)**(l-j))*binomial(i,j)*binomial(n-i,l-j) 
                              for j in range(l+1)])

def roundres(x): # this is a quick and unsafe way to round the result...
   import math
   tol = 0.0001
   if math.ceil(x)-x < tol:
      return int(math.ceil(x))
   if x-math.floor(x) < tol:
      return int(math.floor(x))
   return x

# @cached_function
def delsarte_bound(n, q, d, d_star=1, q_base=0, return_log=True,\
                    isinteger=False, return_data=False):
   p = MixedIntegerLinearProgram(maximization=True)
   A = p.new_variable(integer=isinteger) # A>=0 is assumed
   p.set_objective(sum([A[r] for r in range(n+1)])) 
   for i in range(1,d):
   for j in range(1,n+1): 
      rhs = sum([Kra(n,q,j,r)*A[r] for r in range(n+1)])
      if j >= d_star: 
        p.add_constraint(0*A[0] <= rhs) 
      else: # rhs is proportional to j-th weight of the dual code
        p.add_constraint(0*A[0] == rhs) 
   except sage.numerical.mip.MIPSolverException, exc:
      print "Solver exception: ", exc, exc.args
      if return_data:
         return A,p,False
      return False
   if q_base > 0: # rounding the bound down to the nearest power of q_base,
                  # for q=q_base^m
#      qb = factor(q).radical()
#      if len(qb) == 1:
#         base = qb.expand()
#         bd = base^int(log(bd, base=base))
         if return_log: 
#            bd = int(log(roundres(bd), base=q_base)) # unsafe: 
                                                      # loss of precision
            bd = roundres(log(bd, base=q_base))
#            bd = q_base^int(log(roundres(bd), base=q_base))
            bd = q_base^roundres(log(bd, base=q_base))
   if return_data:
      return A,p,bd
   return bd

d_star is the minimal distance of the dual code (if it exists at all) If q_base is 0, just compute the upper bound. If q_base is >0, it is assumed to be a prime power, and the code assumed to be additive over this field (i.e. the dual code exists,and its weight enumerator is obtained by applying MacWilliams transform --- the matrix A of the LP times the weight enumerator of our code), then the output is the corr. dimension, i.e. floor(log(bound, q_base)).

One obstacle for this to work well in big dimensions is a lack of arbitrary precision LP solver backend available. We are planning to add an PPL backend (PPL is a polyhedral library with a Cython interface and arbitrary precision rational arithmetic, already a standard package in Sage). Once a ticked on this is opened, we'll cite it here.

Change History (4)

comment:1 Changed 10 years ago by dimpase

Last edited 9 years ago by dimpase (previous) (diff)

comment:2 Changed 10 years ago by dimpase

  • Cc kini added

comment:3 Changed 10 years ago by dimpase

  • Cc wdj added

comment:4 Changed 10 years ago by dimpase

  • Cc ncohen added
  • Description modified (diff)
Note: See TracTickets for help on using tickets.