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

function [bnd, ret] = pfenderbnd(K, Tmax, Hmax, epstol)
  pkg load gsl;
  if nargin < 4
    epstol = 1e-4;
  end
  [b,r] = pfenderlp(K,Tmax,Hmax);
  gegenF = @ (t) (-pfenderF(K,r.c,t)); %% for maximization using sqp()
  ret.pfender = 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 = pfenderF(K,r.c,tp);
    %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("pfenderbnd: invalid bound (positivity exceeds %f)\n", epstol);
    bnd = inf;
  else
    if gstar > 0
      fprintf("pfenderbnd: bound positivity in [0,%f], accepting\n", epstol);
    end
    bnd = b;
  end
end

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