function xleja = lejapoints(n,varargin)
% LEJAPOINTS Computes (symmetrized) Leja points.
%
% LEJAPOINTS(N) computes the first N Leja points in [-2,2]. 
% Given the first Leja point X(1) = 2, the Leja points are 
% defined recursively in such a way that the function
%
% F(Y,X(1:N-1)) = ABS(PROD(Y-X(1:N-1)))
%
% is maximized over [-2,2] by X(N).
% If N is not provided, the first 10 Leja poins are computed.
% The maximum of the function above is computed by FMINBND 
% (the algorithm is based on golden section search and parabolic 
% interpolation).
%
% LEJAPOINTS(N,TYPE) computes the first N Leja points in [-2,2] if
% TYPE = 1 and the first N-(1-MOD(N,2)) symmetrized Leja points in
% [-2,2] if TYPE = 2. Symmetrized Leja points are defined as
%
% X(0) = 0
% X(1) = 2
% X(2) = -2
% ...
% X(N) = argmax(F(Y,X(1:N-1))) if N is odd
% X(N+1) = -X(N)               if N is even
%
% LEJAPOINTS(N,X) sorts the vector X in the Leja order, if LENGTH(X) > 1.
% It means that the function above is maximized over X.
if (nargin == 0)
  n = 10;
end
if (nargin > 1)
  x = varargin{1};
else
  x = 1;
end
if (length(x) > 1)
  xleja = lejaorder(x,n);
  return
end
options = [];
xleja = zeros(n,1);
xmaxloc = zeros(n-2,1);
fmaxloc = xmaxloc;
if (x == 1)
  xleja(1) = 2;
  xleja(2) = -2;
  xleja(3) = 0;
else
  xleja(1) = 0;
  xleja(2) = 2;
  xleja(3) = -2;
  n = n-(1-mod(n,2));
end
xsort = xleja;
for k = 4:x:n
  xsort = sort(xsort(1:k-1));
  for i = 1:(k-2)/x
    [xmaxloc(i),fmaxloc(i)] = fminbnd(@fleja,xsort(i),xsort(i+1),options,xleja,k-1);
  end
  [fmax index] = max(-fmaxloc(1:k-2));
  xleja(k) = xmaxloc(index);
  xleja(k+1) = -xmaxloc(index);
  xsort(k) = xleja(k);
  xsort(k+1) = xleja(k+1);
end
xleja = xleja(1:n);


function z = fleja(x,xleja,n)
% the opposite of the Leja function to maximize
z = -abs(prod(x-xleja(1:n)));


function z = lejaorder(x,n)
% Leja order of the first n components of the vector x
m = length(x);
n = min(n,m);
z = zeros(size(x));
[v index] = max(abs(x));
z(1) = x(index);
temp = abs(x-z(1));
[v index] = max(temp);
z(2) = x(index);
for k = 2:n-1
  for i = 1:m
    temp(i) = temp(i)*abs(x(i)-z(k));
  end
  [v index] = max(temp);
  z(k+1) = x(index);
end
z = z(1:n);

function [min minval] = fminbnd(func,lb,ub, options, varargin)
  delta = 1e-17;
  gr = (sqrt(5)-1)/2;
  width = (ub-lb);
  out = [ lb:(width/3):ub ];
  out(2) = out(4)-gr*width;
  out(3) = out(1)+gr*width;
  upper = feval(func,out(3), varargin{:});
  lower = feval(func,out(2), varargin{:});
  while((out(3)-out(2)) > delta) 
    if (upper > lower)
      out(4) = out(3);
      out(3) = out(2);
      width = out(4)-out(1);
      out(2) = out(4)-gr*width;
      upper = lower;
      lower = feval(func,out(2), varargin{:});
    else
      out(1) = out(2);
      out(2) = out(3);
      width = out(4)-out(1);
      out(3) = out(1)+width*gr;
      lower = upper;
      upper = feval(func,out(3), varargin{:});
    end
  end
  min = out(2);
  minval = feval(func,out(2), varargin{:});

