# Ticket #12418: 12418_delsart_bounds.patch

File 12418_delsart_bounds.patch, 20.5 KB (added by dimpase, 8 years ago)

rebased for Sage 5.10 and fixed some outstanding issues

• ## doc/en/reference/coding/index.rst

# HG changeset patch
# User Dmitrii Pasechnik <dimpase@gmail.com>
# Date 1351442118 -28800
# Node ID 71260c4f7840a7527415e947a76a2e9901360206
# Parent  5d8c13e7b81a1e152855073f97a4d48ee97255d1
initial implementation of Delsarte bounds; added Krawchuk speedup by ppurka; more docstrings updated; trailing whitespaces removed; refactored

diff --git a/doc/en/reference/coding/index.rst b/doc/en/reference/coding/index.rst
 a sage/coding/code_constructions sage/coding/sd_codes sage/coding/code_bounds sage/coding/delsarte_bounds sage/coding/source_coding/huffman .. include:: ../footer.txt
• ## sage/coding/all.py

diff --git a/sage/coding/all.py b/sage/coding/all.py
 a from sd_codes import self_dual_codes_binary from delsarte_bounds import (Krawtchouk, delsarte_bound_hamming_space, delsarte_bound_additive_hamming_space)
• ## sage/coding/code_bounds.py

diff --git a/sage/coding/code_bounds.py b/sage/coding/code_bounds.py
 a integer, added example to elias_bound_asymp. - " (2009-05): removed all calls to Guava but left it as an option. - Dima Pasechnik (2012-10): added LP bounds. Let F be a finite field (we denote the finite field with q elements by \GF{q}). A subset C of V=F^n is called a code of .. math:: R={\frac{\log_q\vert C\vert}{n}}, R={\frac{\log_q\vert C\vert}{n}}, where \vert C\vert denotes the number of elements of .. math:: d({\bf v},{\bf w}) =\vert\{i\ \vert\ 1\leq i\leq n,\ v_i\not= w_i\}\vert d({\bf v},{\bf w}) =\vert\{i\ \vert\ 1\leq i\leq n,\ v_i\not= w_i\}\vert to be the Hamming distance between {\bf v} and .. math:: \delta = d/n. \delta = d/n. A linear code with length n, dimension k, and minimum distance -  mrrw1_bound_asymp(delta,q), "first" asymptotic McEliese-Rumsey-Rodemich-Welsh bound for the information rate. -  Delsarte (a.k.a. Linear Programming (LP)) upper bounds. PROBLEM: In this module we shall typically either (a) seek bounds on k, given n, d, q, (b) seek bounds on R, delta, q (assuming n is from sage.rings.arith import factorial from sage.functions.all import log, sqrt from sage.misc.decorators import rename_keyword from delsarte_bounds import delsarte_bound_hamming_space, \ delsarte_bound_additive_hamming_space @rename_keyword(deprecation=6094, method="algorithm") def codesize_upper_bound(n,d,q,algorithm=None): r""" This computes the minimum value of the upper bound using the algorithms of Singleton, Hamming, Plotkin and Elias. This computes the minimum value of the upper bound using the methods of Singleton, Hamming, Plotkin, and Elias. If algorithm="gap" then this returns the best known upper bound A(n,d)=A_q(n,d) for the size of a code of length n, minimum distance d over a field of size q. The function first checks for trivial cases (like d=1 or n=d), and if the value is in the built-in table. Then it calculates the minimum value of the upper bound using the algorithms of Singleton, Hamming, Johnson, Plotkin and Elias. If the code is binary, If algorithm="gap" then this returns the best known upper bound A(n,d)=A_q(n,d) for the size of a code of length n, minimum distance d over a field of size q. The function first checks for trivial cases (like d=1 or n=d), and if the value is in the built-in table. Then it calculates the minimum value of the upper bound using the algorithms of Singleton, Hamming, Johnson, Plotkin and Elias. If the code is binary, A(n, 2\ell-1) = A(n+1,2\ell), so the function takes the minimum of the values obtained from all algorithms for the parameters (n, 2\ell-1) and (n+1, 2\ell). This wraps GUAVA's UpperBound( n, d, q ). wraps GUAVA's (i.e. GAP's package Guava) UpperBound( n, d, q ). If algorithm="LP" then this returns the Delsarte (a.k.a. Linear Programming) upper bound. EXAMPLES:: sage: codesize_upper_bound(10,3,2) 93 sage: codesize_upper_bound(24,8,2,algorithm="LP") 4096 sage: codesize_upper_bound(10,3,2,algorithm="gap")  # optional - gap_packages (Guava package) 85 sage: codesize_upper_bound(11,3,4,algorithm=None) 123361 sage: codesize_upper_bound(11,3,4,algorithm="gap")  # optional - gap_packages (Guava package) 123361 sage: codesize_upper_bound(11,3,4,algorithm="LP") 109226 """ if algorithm=="gap": return int(gap.eval("UpperBound(%s,%s,%s)"%( n, d, q ))) if algorithm=="LP": return int(delsarte_bound_hamming_space(n,d,q)) else: eub = elias_upper_bound(n,q,d) gub = griesmer_upper_bound(n,q,d) sub = singleton_upper_bound(n,q,d) return min([eub,gub,hub,pub,sub]) def dimension_upper_bound(n,d,q): @rename_keyword(deprecation=6094, method="algorithm") def dimension_upper_bound(n,d,q,algorithm=None): r""" Returns an upper bound B(n,d) = B_q(n,d) for the dimension of a linear code of length n, minimum distance d over a field of size q. Parameter "algorithm" has the same meaning as in :func:codesize_upper_bound EXAMPLES:: sage: dimension_upper_bound(10,3,2) 6 sage: dimension_upper_bound(30,15,4) 13 sage: dimension_upper_bound(30,15,4,algorithm="LP") 12 """ q = ZZ(q) return int(log(codesize_upper_bound(n,d,q),q)) if algorithm=="LP": return delsarte_bound_additive_hamming_space(n,d,q) else:       # algorithm==None or algorithm="gap": return int(log(codesize_upper_bound(n,d,q,algorithm=algorithm),q)) def volume_hamming(n,q,r): r""" Returns number of elements in a Hamming ball of radius r in \GF{q}^n. Agrees with Guava's SphereContent(n,r,GF(q)). EXAMPLES:: sage: volume_hamming(10,2,3) 176 """ r""" Returns lower bound for number of elements in the largest code of minimum distance d in \GF{q}^n. EXAMPLES:: sage: gilbert_lower_bound(10,2,3) 128/7 """ def plotkin_upper_bound(n,q,d, algorithm=None): r""" Returns Plotkin upper bound for number of elements in the largest code of minimum distance d in \GF{q}^n. code of minimum distance d in \GF{q}^n. The algorithm="gap" option wraps Guava's UpperBoundPlotkin. EXAMPLES:: sage: plotkin_upper_bound(10,2,3) sage: plotkin_upper_bound(10,2,3) 192 sage: plotkin_upper_bound(10,2,3,algorithm="gap")  # optional - gap_packages (Guava package) 192 Returns the Griesmer upper bound for number of elements in the largest code of minimum distance d in \GF{q}^n. Wraps GAP's UpperBoundGriesmer. EXAMPLES:: sage: griesmer_upper_bound(10,2,3) 128 sage: griesmer_upper_bound(10,2,3,algorithm="gap")  # optional - gap_packages (Guava package) if not(add == 1): if d%den==0: add = int(d/den) else: else: add = int(d/den)+1 s = s + add den = den * q Returns the Elias upper bound for number of elements in the largest code of minimum distance d in \GF{q}^n. Wraps GAP's UpperBoundElias. EXAMPLES:: sage: elias_upper_bound(10,2,3) 232 sage: elias_upper_bound(10,2,3,algorithm="gap")  # optional - gap_packages (Guava package) Returns the Hamming upper bound for number of elements in the largest code of minimum distance d in \GF{q}^n. Wraps GAP's UpperBoundHamming. The Hamming bound (also known as the sphere packing bound) returns an upper bound on the size of a code of length n, minimum distance d, over a field of size q. The Hamming bound is obtained by \GF{q}^n by the contents of a ball with radius floor((d-1)/2). As all these balls are disjoint, they can never contain more than the whole vector space. .. math:: M \leq {q^n \over V(n,e)}, M \leq {q^n \over V(n,e)}, where M is the maximum number of codewords and V(n,e) is equal to the contents of a ball of radius e. This bound is useful for small values of d. Codes for which equality holds are called perfect. EXAMPLES:: sage: hamming_upper_bound(10,2,3) sage: hamming_upper_bound(10,2,3) 93 """ return int((q**n)/(volume_hamming(n, q, int((d-1)/2)))) Returns the Singleton upper bound for number of elements in the largest code of minimum distance d in \GF{q}^n. Wraps GAP's UpperBoundSingleton. This bound is based on the shortening of codes. By shortening an (n, M, d) code d-1 times, an (n-d+1,M,1) code results, with M \leq q^n-d+1. Thus .. math:: M \leq q^{n-d+1}. M \leq q^{n-d+1}. Codes that meet this bound are called maximum distance separable (MDS). EXAMPLES:: sage: singleton_upper_bound(10,2,3) sage: singleton_upper_bound(10,2,3) 256 """ return q**(n - d + 1) """ GV lower bound for information rate of a q-ary code of length n minimum distance delta\*n EXAMPLES:: sage: RDF(gv_info_rate(100,1/4,3)) 0.367049926083 """ def gv_bound_asymp(delta,q): """ Computes the asymptotic GV bound for the information rate, R. EXAMPLES:: sage: RDF(gv_bound_asymp(1/4,2)) 0.188721875541 sage: f = lambda x: gv_bound_asymp(x,2) def hamming_bound_asymp(delta,q): """ Computes the asymptotic Hamming bound for the information rate. EXAMPLES:: sage: RDF(hamming_bound_asymp(1/4,2)) 0.4564355568 sage: f = lambda x: hamming_bound_asymp(x,2) def singleton_bound_asymp(delta,q): """ Computes the asymptotic Singleton bound for the information rate. EXAMPLES:: sage: singleton_bound_asymp(1/4,2) 3/4 sage: f = lambda x: singleton_bound_asymp(x,2) """ Computes the asymptotic Plotkin bound for the information rate, provided 0 < \delta < 1-1/q. EXAMPLES:: sage: plotkin_bound_asymp(1/4,2) 1/2 """ def elias_bound_asymp(delta,q): """ Computes the asymptotic Elias bound for the information rate, provided 0 < \delta 1-1/q. provided 0 < \delta < 1-1/q. EXAMPLES:: sage: elias_bound_asymp(1/4,2) 0.39912396330... """ """ Computes the first asymptotic McEliese-Rumsey-Rodemich-Welsh bound for the information rate, provided 0 < \delta < 1-1/q. EXAMPLES:: sage: mrrw1_bound_asymp(1/4,2) 0.354578902665 """ return RDF(entropy((q-1-delta*(q-2)-2*sqrt((q-1)*delta*(1-delta)))/q,q))
• ## new file sage/coding/delsarte_bounds.py

diff --git a/sage/coding/delsarte_bounds.py b/sage/coding/delsarte_bounds.py
new file mode 100644
 - r""" Delsarte, a.k.a. Linear Programming (LP), upper bounds. This module provides  LP upper bounds for the parameters of codes. Exact LP solver, PPL, is used by defaut, ensuring that no rounding/overflow problems occur. AUTHORS: - Dmitrii V. (Dima) Pasechnik (2012-10): initial implementation. """ #***************************************************************************** #       Copyright (C) 2012 Dima Pasechnik # #  Distributed under the terms of the GNU General Public License (GPL) # #                  http://www.gnu.org/licenses/ #***************************************************************************** def Krawtchouk(n,q,l,i): """ Compute K^{n,q}_l(i), the Krawtchouk polynomial: see en.wikipedia.org/wiki/Kravchuk_polynomials _. It is given by .. math:: K^{n,q}_l(i)=\sum_{j=0}^l (-1)^j(q-1)^{(l-j)}{i \choose j}{n-i \choose l-j} EXAMPLES:: sage: Krawtchouk(24,2,5,4) 2224 sage: Krawtchouk(12300,4,5,6) 567785569973042442072 """ from sage.rings.arith import binomial # Use the expression in equation (55) of MacWilliams & Sloane, pg 151 # We write jth term = some_factor * (j-1)th term kraw = jth_term = (q-1)**l * binomial(n, l) # j=0 for j in range(1,l+1): jth_term *= -q*(l-j+1)*(i-j+1)/((q-1)*j*(n-j+1)) kraw += jth_term return kraw def _delsarte_LP_building(n, d, d_star, q, isinteger,  solver, maxc = 0): """ LP builder - common for the two functions; not exported. """ from sage.numerical.mip import MixedIntegerLinearProgram p = MixedIntegerLinearProgram(maximization=True, solver=solver) A = p.new_variable(integer=isinteger) # A>=0 is assumed p.set_objective(sum([A[r] for r in xrange(n+1)])) p.add_constraint(A[0]==1) for i in xrange(1,d): p.add_constraint(A[i]==0) for j in xrange(1,n+1): rhs = sum([Krawtchouk(n,q,j,r)*A[r] for r in xrange(n+1)]) p.add_constraint(0*A[0] <= rhs) 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) if maxc > 0: p.add_constraint(sum([A[r] for r in xrange(n+1)]), max=maxc) return A, p def delsarte_bound_hamming_space(n, d, q, isinteger=False, return_data=False, solver="PPL"): """ Find the classical Delsarte bound [1]_ on codes in Hamming space H_q^n of minimal distance d INPUT: - n -- the code length - d -- the (lower bound on) minimal distance of the code - q -- the size of the alphabet - isinteger -- if True, uses an integer programming solver (ILP), rather that an LP solver. Can be very slow if set to True. - return_data -- if True, return a triple (W,LP,bound), where W is a weights vector,  and LP the Delsarte bound LP; both of them are Sage LP data.  W need not be a weight distribution of a code, or, if isinteger==False, even have integer entries. - solver -- the LP/ILP solver to be used. Defaults to PPL. It is arbitrary precision, thus there will be no rounding errors. With other solvers (see :class:MixedIntegerLinearProgram for the list), you are on your own! EXAMPLES: The bound on the size of the F_2-codes of length 11 and minimal distance 6:: sage: delsarte_bound_hamming_space(11, 6, 2) 12 sage: a, p, val = delsarte_bound_hamming_space(11, 6, 2, return_data=True) sage: [j for i,j in p.get_values(a).iteritems()] [1, 0, 0, 0, 0, 0, 11, 0, 0, 0, 0, 0] The bound on the size of the F_2-codes of length 24 and minimal distance 8, i.e. parameters of the extened binary Golay code:: sage: a,p,x=delsarte_bound_hamming_space(24,8,2,return_data=True) sage: x 4096 sage: [j for i,j in p.get_values(a).iteritems()] [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] The bound on the size of F_4-codes of length 11 and minimal distance 3:: sage: delsarte_bound_hamming_space(11,3,4) 327680/3 REFERENCES: .. [1] P. Delsarte, An algebraic approach to the association schemes of coding theory, Philips Res. Rep., Suppl., vol. 10, 1973. """ A, p = _delsarte_LP_building(n, d, 0, q, isinteger,  solver) 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 return_data: return A,p,bd else: return bd def delsarte_bound_additive_hamming_space(n, d, q, d_star=1, q_base=0, isinteger=False, return_data=False, solver="PPL"): """ Find the Delsarte LP bound on F_{q_base}-dimension of additive codes in Hamming space H_q^n of minimal distance d with minimal distance of the dual code at least d_star.  If q_base is set to non-zero, then  q is a power of q_base, and the code is, formally, linear over F_{q_base}. Otherwise it is assumed that q_base==q. INPUT: - n -- the code length - d -- the (lower bound on) minimal distance of the code - q -- the size of the alphabet - d_star -- the (lower bound on) minimal distance of the dual code; only makes sense for additive codes. - q_base -- if 0, the code is assumed to be nonlinear. Otherwise, q=q_base^m and the code is linear over F_{q_base}. - isinteger -- if True, uses an integer programming solver (ILP), rather that an LP solver. Can be very slow if set to True. - return_data -- if True, return a triple (W,LP,bound), where W is a weights vector,  and LP the Delsarte bound LP; both of them are Sage LP data.  W need not be a weight distribution of a code, or, if isinteger==False, even have integer entries. - solver -- the LP/ILP solver to be used. Defaults to PPL. It is arbitrary precision, thus there will be no rounding errors. With other solvers (see :class:MixedIntegerLinearProgram for the list), you are on your own! EXAMPLES: The bound on dimension of linear F_2-codes of length 11 and minimal distance 6:: sage: delsarte_bound_additive_hamming_space(11, 6, 2) 3 sage: a,p,val=delsarte_bound_additive_hamming_space(11, 6, 2,\ return_data=True) sage: [j for i,j in p.get_values(a).iteritems()] [1, 0, 0, 0, 0, 0, 5, 2, 0, 0, 0, 0] The bound on the dimension of linear F_4-codes of length 11 and minimal distance 3:: sage: delsarte_bound_additive_hamming_space(11,3,4) 8 The bound on the F_2-dimension of additive F_4-codes of length 11 and minimal distance 3:: sage: delsarte_bound_additive_hamming_space(11,3,4,q_base=2) 16 """ if q_base == 0: q_base = q kk = 0 while q_base**kk < q: kk += 1 if q_base**kk != q: print "Wrong q_base=", q_base, " for q=", q, kk return False # this implementation assumes that our LP solver to be unable to do a hot # restart with an adjusted constraint m = kk*n # this is to emulate repeat/until block bd = q**n+1 while q_base**m < bd: # need to solve the LP repeatedly, as this is a new constraint! # we might become infeasible. More precisely, after rounding down # to the closest value of q_base^m, the LP, with the constraint that # the objective function is at most q_base^m, A, p = _delsarte_LP_building(n, d, d_star, q, isinteger,  solver, q_base**m) 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 # rounding the bound down to the nearest power of q_base, for q=q_base^m #      bd_r = roundres(log(bd, base=q_base)) m = -1 while q_base**(m+1) < bd: m += 1 if q_base**(m+1) == bd: m += 1 if return_data: return A, p, m else: return m