Opened 10 years ago

Last modified 8 years ago

#12418 closed enhancement

adding Delsarte bound for codes — at Version 9

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)])) 
   p.add_constraint(A[0]==1)
   for i in range(1,d):
      p.add_constraint(A[i]==0)
   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) 
   try:
      bd=p.solve()
   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))
         else:
#            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 in Sage. This is (almost - i.e. the corresponding ticket is still not 100% ready, as reviewers think) taken care of by #12533, which should be made a dependence for this ticket.

Change History (9)

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)

comment:5 Changed 10 years ago by ncohen

(GLPK can solve non-integer rational LP. It is not exposed, but may not be too hard either)

comment:6 Changed 10 years ago by dimpase

  • Cc ppurka added

comment:7 Changed 10 years ago by ppurka

The function named delsarte_bound should be renamed to something like delsarte_bound_hamming_space. This is so that in future other functions like delsarte_bound_johnson_space, delsarte_bound_permutation_space, etc can be added easily, without having inconsistencies in naming.

comment:8 Changed 10 years ago by dimpase

  • Cc ptrrsn_1 added

comment:9 Changed 10 years ago by dimpase

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