#12418: adding Delsarte bound for codes
---------------------------------+------------------------------------------
Reporter: dimpase | Owner: wdj
Type: enhancement | Status: new
Priority: major | Milestone: sage-5.1
Component: coding theory | Resolution:
Keywords: | Work issues:
Report Upstream: N/A | Reviewers:
Authors: | Merged in:
Dependencies: | Stopgaps:
---------------------------------+------------------------------------------
Description changed by dimpase:
Old description:
> 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. We are planning to add
> an [http://www.sagemath.org/doc/reference/sage/libs/ppl.html 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.
New description:
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.
--
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/12418#comment:9>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en.
