function [r,pol,res,zer,err,S] = aaax(f,deg,nl,tol,plt)
%AAAX   Rational approximation on [-1,1].
%   [r,pol,res,zer,err,S] = aaax(f,deg,nlawson,tol,plt)
%Examples:
%   aaax(@exp); aaax(@exp,5); aaax(@exp,5,20);
%   [r,pol] = aaax(@(x) tanh(20*pi*x)); pol
%   aaax(@abs,150,0,1e-10);

                                               % INITIALIZATIONS
if nargin<2, deg = 150; end                    % arg 2: (max) degree
if nargin<3, nl = 0; end                       % arg 3: Lawson steps, eg 20
if nargin<4, tol = 1e-13; end                  % arg 4: relative tolerance
if nargin<5, plt = 1; end                      % arg 5: plt=0 stops plotting
blue = [0 0 .8]; red = [.8 0 0];
grn = [.5 .8 .5];  MS = 'markersize';
CO = 'color'; LW = 'linewidth'; 
S = [-1 1];                                    % row vector of support pts
f0 = f([S'; XS(S,10)]); err = std(f0);
if (err==0) | (err/abs(mean(f0))<=tol) ...     % constant function or deg 0:
            | (deg==0)                         % return (no plot)
   r = @(x) mean(f0) + 0*x; return
end 
m0 = 2; S0 = S; w0 = [1; -1]; err0 = err;

while true                                     % MAIN AAA LOOP
   m = length(S);
   X = XS(S,max(3,16-m));                      % column vector of sample pts
   C = 1./(X-S);                               % Cauchy matrix
   A = (f(X)-f(S)).*C;                         % Loewner matrix
   [~,~,V] = svd(A,0);                         % SVD
   w = V(:,end);                               % barycentric wt vector
   R = (C*(w.*f(S')))./(C*w);                  % rat approx at sample pts
   err = norm(f(X)-R,inf);                     % abs error of this approx
   pol = prz(S',f(S'),w);                      % poles
   bad = any(imag(pol)==0 & abs(pol)<=1);      % check for bad poles
   co = blue;  
   fmax = norm(f([S';X]),inf);
   if ~bad & err < err0
      m0 = m; S0 = S; w0 = w; err0 = err;      % save latest success
   end
   if plt
      subplot(211), if bad, co = red; end
      semilogy(m-1,err,'.',MS,10,CO,co), hold on
      title convergence, ylim([1e-16 1e4]), grid on
      set(gca,'ytick',10.^(-16:8:0)), xlabel('degree (=m-1)')
   end
   if ~bad & (err/fmax<=tol) | ...             % stop if converged
      (m == deg+1) | ...                       % stop if max degree reached
      ((m-m0>=10) & (err0/fmax<1e-2))          % stop if stagnated
      break
   end 
   [~,j] = max(abs(f(X)-R)); S = [S X(j)];     % next support point
end
m = m0; S = S0; w = w0; err = err0; 
[pol,res,zer] = prz(S',f(S'),w); 
r = @(x) reval(x,S',f(S'),w);                  % function handle for r

if nl > 0                                      % LAWSON ITERATION
   X = XS(S,20); n = length(X);                % switch to finer grid
   wt = ones(m+n,1);                           % Lawson wt vector
   C = 1./(X-S);                               % Cauchy matrix
   R = (C*(w.*f(S')))./(C*w);                  % rat approx at sample pts
   F = [f(X); f(S')]; R = [R; f(S')];          % f and r values
   L = f(X).*C;
   A = [C; eye(m)]; B = [L; diag(f(S))];       % numer and f*denom rows
   for stepno = 1:nl                           % take nl Lawson steps
      [~,~,V] = svd(sqrt(wt).*[A B/fmax],0);
      c = V(:,2*m);
      a = c(1:m); w = -c(m+1:2*m)/fmax;        % numer wts, denom wts
      R = (A*a)./(A*w);
      wt = wt.*abs(F-R);                       % iterative reweighting
      wt = wt/norm(wt,inf);
   end
   r = @(x) reval(x,S',a./w,w);                % function handle for r
end

xx = sort([XS(S,30); S']);                     % COMPUTE ERROR AND PLOT
ee = f(xx)-r(xx); err = norm(ee,inf);
if plt
   semilogy(m-1,max(err,1e-16),'o',LW,2,CO,grn)
   hold off, subplot(212), co = [0 .45 .74];
   if nl > 0, co = [.7 0 .9]; end, plot(xx,ee,LW,.7,CO,co)
   s = sprintf('error: max =%8.1e',err); ylim(1.5*(err+realmin)*[-1 1])
   title(s), grid on, xlabel x
end
end

function XS = XS(S,p)                          % SAMPLE POINTS
S = sort(S); d = (1:p)/(p+1);
XS = S(1:end-1) + d'.*diff(S);
XS = XS(:);
end

function [pol,res,zer] = prz(xj,fj,wj)         % POLES, RESIDUES, ZEROS
m = length(wj);                                %    (from Chebfun)
if any(wj==0)
   ii = find(wj~=0); m = length(ii);
   xj = xj(ii); fj = fj(ii); wj = wj(ii);
end   
B = eye(m+1); B(1,1) = 0;
E = [0 wj.'; ones(m,1) diag(xj)];
pol = eig(E,B); pol = pol(~isinf(pol));
if nargout >= 2
   N = @(t)(1./(t-xj.'))*(fj.*wj);
   Ddiff = @(t) -((1./(t-xj.')).^2)*wj;
   res = N(pol)./Ddiff(pol);
   E = [0 (wj.*fj).'; ones(m,1) diag(xj)];
   zer = eig(E,B); zer = zer(~isinf(zer));
end
end

function r = reval(zz,xj,fj,wj)                % FUNCTION HANDLE FOR R
z = zz(:); C = 1./(z-xj.');                    %    (from Chebfun)    
r = (C*(wj.*fj))./(C*wj);
r(isinf(z)) = sum(wj.*fj)./sum(wj);
ii = find(isnan(r));
for j = 1:length(ii)
   r(ii(j)) = fj(z(ii(j))==xj);
end
r = reshape(r,size(zz));
end
