Changes between Version 9 and Version 10 of Ticket #12418


Ignore:
Timestamp:
06/08/12 14:44:48 (9 years ago)
Author:
dimpase
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • Ticket #12418

    • Property Dependencies changed from to #12533
  • Ticket #12418 – Description

    v9 v10  
    11Delsarte bound for codes, aka Linear Programming bound, is easy to implement in Sage.
    2 Here is a quick and dirty code that does it; so the problem would be to integrate this properly.
    3 {{{
    4 def Kra(n,q,l,i): # K^{n,q}_l(i), i.e Krawtchouk polynomial
    5    return sum([((-1)**j)*((q-1)**(l-j))*binomial(i,j)*binomial(n-i,l-j)
    6                               for j in range(l+1)])
    72
    8 def roundres(x): # this is a quick and unsafe way to round the result...
    9    import math
    10    tol = 0.0001
    11    if math.ceil(x)-x < tol:
    12       return int(math.ceil(x))
    13    if x-math.floor(x) < tol:
    14       return int(math.floor(x))
    15    return x
     3See the attached prototype implementation for details.
    164
    17 # @cached_function
    18 def delsarte_bound(n, q, d, d_star=1, q_base=0, return_log=True,\
    19                     isinteger=False, return_data=False):
    20    p = MixedIntegerLinearProgram(maximization=True)
    21    A = p.new_variable(integer=isinteger) # A>=0 is assumed
    22    p.set_objective(sum([A[r] for r in range(n+1)]))
    23    p.add_constraint(A[0]==1)
    24    for i in range(1,d):
    25       p.add_constraint(A[i]==0)
    26    for j in range(1,n+1):
    27       rhs = sum([Kra(n,q,j,r)*A[r] for r in range(n+1)])
    28       if j >= d_star:
    29         p.add_constraint(0*A[0] <= rhs)
    30       else: # rhs is proportional to j-th weight of the dual code
    31         p.add_constraint(0*A[0] == rhs)
    32    try:
    33       bd=p.solve()
    34    except sage.numerical.mip.MIPSolverException, exc:
    35       print "Solver exception: ", exc, exc.args
    36       if return_data:
    37          return A,p,False
    38       return False
    39    if q_base > 0: # rounding the bound down to the nearest power of q_base,
    40                   # for q=q_base^m
    41 #      qb = factor(q).radical()
    42 #      if len(qb) == 1:
    43 #         base = qb.expand()
    44 #         bd = base^int(log(bd, base=base))
    45          if return_log:
    46 #            bd = int(log(roundres(bd), base=q_base)) # unsafe:
    47                                                       # loss of precision
    48             bd = roundres(log(bd, base=q_base))
    49          else:
    50 #            bd = q_base^int(log(roundres(bd), base=q_base))
    51             bd = q_base^roundres(log(bd, base=q_base))
    52    if return_data:
    53       return A,p,bd
    54    return bd
    55 }}}
     5One 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 is a dependence for this ticket.
    566
    57 d_star is the minimal distance  of the dual code (if it exists at all)
    58 If q_base is 0, just compute the upper bound.
    59 If q_base is >0, it is assumed to be a prime power, and the code assumed
    60 to be additive over this field (i.e. the dual code exists,and its weight enumerator is
    61 obtained by applying {{{MacWilliams transform}}} --- the matrix A of the LP times the
    62 weight enumerator of our code), then the output is the corr. dimension, i.e.
    63 {{{floor(log(bound, q_base))}}}.
    64 
    65 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.
    66