% [bnd,ret] = delsartebnd(K, Tmax, Hmax)
%
% Verify that the value returned by delsartelp() is a valid KNP bound
%  (epsilon method by Odlyzko & Sloane; epstol is a "positive zero")

function [bnd, ret] = delsartebnd(K, Tmax, Hmax, epstol)
  pkg load gsl;
  if nargin < 4
    epstol = 1e-4;
  end
  [b,r] = delsartelp(K,Tmax,Hmax);
  gegenF = @ (t) (-delsarteF(K,r.c,t)); %% for maximization using sqp()
  ret.delsarte = r;
  
  % discretization  
  t = zeros(1,Tmax);
  discreps = 1.5 / Tmax;
  for j=1:Tmax
    t(j) = -1.0 + (j-1)*discreps;
  end
  t(Tmax) = 0.5;
  
  % maximize LHS in constraints w.r.t. t
  gstar = -inf;
  tstar = 1;
  for j = 1:Tmax
    [tp,gval,info] = sqp(t(j),gegenF,[],[],[-1],[0.5],10,1e-3);
    gval = delsarteF(K,r.c,tp) + 1; % because test is <= -1
    %fprintf("j=%d, t=%f, gval=%f\n", j,tp,gval); 
    if gval > gstar
      gstar = gval;
      tstar = tp;
    end
  end
  
  ret.gstar = gstar;
  ret.tstar = tstar;
  
  if gstar > epstol
    % solution infeasible in SILP
    fprintf("delsartebnd: invalid bound (positivity exceeds %f)\n", epstol);
    bnd = inf;
  else
    if gstar > 0
      fprintf("delsartebnd: bound positivity in [0,%f], accepting\n", epstol);
    end
    bnd = b;
  end
end

% the 'F' function (conic comb of Gegenbauer poly)
function g = delsarteF(K,c,t)
  Hmax = length(c);
  g = 0;
  for i=1:Hmax
    g = g + c(i)*gegenbauer(K,i,t);   
  end
end

% derivative of F (see <http://mathworld.wolfram.com/GegenbauerPolynomial.html>)
function gd = delsarteFd(K,c,t)
  Hmax = length(c);
  gd = 0;
  Kder = K+2; % so (Kder-2)/2 = K/2 = (K-2)/2+1
  for i=1:Hmax
    gd = gd + c(i)*((K-2)/2)*gegenbauer(Kder,i-1,t);
  end
end