Source

qcode / j.sage

Full commit
# Askey-Wilson polynomials

# (a;q)_n:=prod_{k=0}^n-1 (1-aq^k)
def br_n(q,n,a): return prod([1-a*q**k for k in range(n)])

# (a_1,a_2,...,a_k;q)_n:=prod_{j=1}^k(a_j;q)_n
def brlist_n(q,n,*a): 
   def br(x): return prod([1-x*q**k for k in range(n)])
   return prod([br(aj) for aj in a])

# a and b are sequences
def q_hyper(q,z,a,b,ulim=oo):
   s = len(a)
   if s-1 == len(b):
      if ulim==oo: # do a symbolic summation
         var('kkk') 
         return sum(z**kkk*brlist_n(q,kkk,*a)/brlist_n(q,kkk,*(b+[q])), 
                     kkk, 0, ulim)
      else:
         return \
           reduce(lambda x,y: x+y, \
              [z**k*brlist_n(q,k,*a)/brlist_n(q,k,*(b+[q])) for \
                     k in range(ulim+1)])

# q-binomial coefficient
def qbinomial(q,v,j):
  if j < 0:
     return 0
  return q^(v*j-j*(j-1)/2)*(-1)^j*br_n(q,j,q^(-v))/br_n(q,j,q)

def Eber(v,k,i,u):
    if 2*k<=v:
     return sum(
       [(-1)^t*binomial(u,t)*binomial(k-u,i-t)*binomial(v-k-u,i-t)
          for t in range(i+1)])
    return Eber(v,v-k,i,u)

# P and Q for J(v,k), the Johnson scheme
def JP(v,k):
    def ff(i,j): return Eber(v,k,i,j)
    return matrix(k+1,k+1,ff)

def JQ(v,k):
    return binomial(v,k)*JP(v,k)^-1

def delsarte_bound_J(n, k, d, \
                    isinteger=False, return_data=False, solver="PPL"):
   
    
   from sage.numerical.mip import MixedIntegerLinearProgram
   if not isinstance(d, list):
      d = range(1, d)
   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 range(k+1)])) 
   p.add_constraint(A[0]==1)
   for i in d:
      p.add_constraint(A[i]==0)
   Q = JQ(n,k)
   Q = diagonal_matrix([y.denominator() for y in Q])*Q # we can scale, 
                                                       # as the RHS==0
   for j in range(1,k+1): 
      rhs = sum([Q[j,r]*A[r] for r in range(k+1)])
      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 return_data:
      return A,p,bd
   else: 
      return bd