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

function [bound, ret] = pfenderlp(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 (w.r.t. Delsarte, we also need c0 here)
  f = zeros(1,Hmax+1);
  f(1) = 1; % obj fun is max c0
  % constr RHS (one more ineq w.r.t. Delsarte)
  b = zeros(Tmax+1,1); % all ineqs but last are <= 0
  b(Tmax+1,1) = 1; % only last ineq is <= 1
  % constraint matrix (one more ineq w.r.t. Delsarte)
  A = zeros(Tmax+1, Hmax+1);
  A(:,1) = 1; % add +c0 to each inequality 
  for j = 1:Tmax
    for h = 1:Hmax
      A(j,h+1) = gegenbauer(K,h,t(j));
    end
  end
  % last constraint in Pfender (absent in Delsarte)
  for h = 1:Hmax
    A(Tmax+1,h+1) = gegenbauer(K,h,1); % =1 if scaled Gegenbauer
  end
  % var bounds: [0,infty]
  lb = zeros(Hmax+1,1); % >= 0
  ub = Inf*ones(Hmax+1,1); % no upper bound
  % constraint, objective sense and variable types
  consense = char("U"*ones(1,Tmax+1)); % <= ineq
  objsense = -1; % max c0
  vartype = char("C"*ones(1,Hmax+1)); % all continuous variables
  param.msglev = 1; % 
  % call GLPK
  [c,objval,flag,glpkret] = glpk(f,A,b,lb,ub,consense,vartype,objsense,param);
  bound = floor(1/objval); % different w.r.t. Delsarte bound
  
  ret.glpk = glpkret;
  ret.c = c;
  ret.objval = objval;
  ret.bound = bound;
  ret.cpu = toc;
end
