%  [bound, ret] = delsartelp(K, Tmax, Hmax)
%
% Find the Delsarte bound in K dimensions with constraint discretization
% T = 1..Tmax and Gegenbauer polynomials up to degree Hmax
%
% by Liberti 180209

function [bound, ret] = delsartelp(K, Tmax, Hmax)
  tic;
  
  % discretization indices
  t = zeros(1,Tmax);
  % discretization
  discreps = 1.5 / Tmax;
  for j=1:Tmax
    t(j) = -1.0 + (j-1)*discreps;
  end
  t(Tmax) = 0.5;
  % initial value for SILP feasibility verification epsilon
  epsilon = 0;  

  %% construction of the LP: Hmax vars and Tmax constrs  
  % obj fun vector
  f = ones(1,Hmax);
  % constr RHS
  b = -ones(Tmax,1);
  % constraint matrix
  A = zeros(Tmax, Hmax);
  for j = 1:Tmax
    for h = 1:Hmax
      A(j,h) = gegenbauer(K,h,t(j));
    end
  end
  % var bounds: [0,infty]
  lb = zeros(Hmax,1); % >= 0
  ub = Inf*ones(Hmax,1); % no upper bound
  % constraint, objective sense and variable types
  consense = char("U"*ones(1,Tmax)); % <= ineq
  objsense = 1; % minimization
  vartype = char("C"*ones(1,Hmax)); % all continuous variables
  param.msglev = 1; % 
  % call GLPK
  [c,objval,flag,glpkret] = glpk(f,A,b,lb,ub,consense,vartype,objsense,param);
  bound = floor(objval + 1);
  
  ret.glpk = glpkret;
  ret.c = c;
  ret.objval = objval;
  ret.bound = bound;
  ret.cpu = toc;
end
