text:
function varargout = minimax(f, varargin)
%MINIMAX   Best polynomial or rational approximation for real valued
%          continuous functions. This code supersedes REMEZ.
%
%   P = MINIMAX(F, M) computes the minimax polynomial approximation of
%   degree M to the real function F using the Remez algorithm. F can
%   be either a CHEBFUN, a function handle or a string representation
%   of the function to approximate.  P is a CHEBFUN.
%
%   [P, Q] = MINIMAX(F, M, N) computes the minimax rational approximation
%   P/Q of type (M, N).  P and Q are CHEBFUNs, but in difficult cases
%   working with P/Q is numerically unstable.
%
%   [P, Q, R_HANDLE] = MINIMAX(F, M, N) additionally returns a numerically
%   stable function handle for evaluating P/Q (based on a barycentric
%   representation).
%
%   [...] = MINIMAX(..., [A, B]) takes the approximation domain to be
%   [A, B]. If a domain is not specified and F is a CHEBFUN, then the
%   domain of F is used. In all other cases, [-1, 1] is used.
%
%   [...] = MINIMAX(..., 'tol', TOL) uses the value TOL as the termination
%   tolerance on the relative equioscillation error.  The default is
%   approximately 1e-12 for polynomial approximation and 1e-4 for
%   rational approximation.
%
%   [...] = MINIMAX(..., 'display', 'iter') displays output at each
%   iteration.
%
%   [...] = MINIMAX(..., 'maxiter', MAXITER) sets the maximum number of
%   allowable iterations to MAXITER.
%
%   [...] = MINIMAX(..., 'init', XK) allows the user to specify the vector
%   XK as the starting reference.
%
%   [...] = MINIMAX(..., 'plot', 'on'), or equivalently
%   [...] = MINIMAX(..., 'plotfcns', 'error') plots the error after each
%   iteration while the algorithm executes.
%
%   [...] = MINIMAX(..., 'silent') turns off all messages regarding the
%   execution of the algorithm.
%
%   [P, ERR] = MINIMAX(...) and [P, Q, R_HANDLE, ERR] = MINIMAX(...)
%   return the maximum error estimate ERR.
%
%   [P, ERR, STATUS] = MINIMAX(...) and [P, Q, R_HANDLE, ERR, STATUS] =
%   MINIMAX(...) return a structure array STATUS with the following fields:
%       STATUS.DELTA - Tolerance obtained.
%       STATUS.ITER  - Number of iterations performed.
%       STATUS.DIFFX - Maximum correction in last trial reference.
%       STATUS.XK    - Last trial reference on which the error
%                      equioscillates.
%   In case we are doing rational approximation (denominator degree >=1),
%   two extra fields are computed:
%       STATUS.POL   - Poles of the minimax approximation.
%       STATUS.ZER   - Zeros of the minimax approximation.
%
%   This code is highly reliable for polynomial approximation but may
%   sometimes have difficulties in the rational case, though we believe
%   it is the most powerful rational minimax code available.
%
% Examples:
%   x = chebfun('x'); f = abs(x);
%   p = minimax(f, 20); plot(f-p)
%   [p, q, rh] = minimax(f, 10, 10);
%   xx = linspace(-1,1,10000); plot(xx, f(xx)-rh(xx))
%
% References:
%
%   [1] B. Beckermann, S. Filip, Y. Nakatsukasa and L. N. Trefethen,
%   "Rational minimax approximation via adaptive barycentric
%   representations", arXiv:1705.10132.
%
%   [2] R. Pachon and L. N. Trefethen, "Barycentric-Remez algorithms for
%   best polynomial approximation in the chebfun system", BIT Numerical
%   Mathematics, 49:721-742, 2009.
%
% See also AAA, CF, CHEBPADE, PADEAPPROX, RATINTERP, POLYFIT, POLYFITL1.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

if ( ~isa(f,'chebfun') ) % check if input is chebfun; if not, look for
                         % splitting points
    dom = [];
    domIndex = 0;

    for k = 1:length(varargin) % look for domain
        if isfloat(varargin{k})
            if length(varargin{k}) == 2
                dom = varargin{k};
                domIndex = k;
            end
        end
    end

    if isempty(dom) % domain not provided, default to [-1,1]
        dom = [-1 1];
    else
        varargin(domIndex) = [];
    end
    if ( ischar(f) )
        fHandle = str2op(vectorize(f));
    else
        fHandle = f;
    end
    f = chebfun(f, dom, 'splitting', 'on');
else % f is a chebfun input
    fHandle = @(x) feval(f, x);
end

if ( ~isreal(f) )
    error('CHEBFUN:CHEBFUN:minimax:real', ...
        'MINIMAX only supports real valued functions.');
end

if ( numColumns(f) > 1 )
    error('CHEBFUN:CHEBFUN:minimax:quasi', ...
        'MINIMAX does not currently support quasimatrices.');
end

% Parse the inputs.
[m, n, N, rationalMode, polyOutput, symFlag, xk, opts] = ...
                                   parseInputs(f, varargin{:});

% If m = -1, this means f = odd and input (m,n) = (0,n); return constant 0.
if ( m == -1 )
    q = chebfun(1, f.domain([1, end]));
    p = chebfun(0, f.domain([1, end]));
    varargout = {p, q, @(x) feval(p, x)./feval(q, x), norm(f,'inf'), []};
    return
end

if ( isempty(xk) ) % no initial reference is given by the user
    % Several initialization attempts are made
    if ( n == 0 ) % polynomial case
        % Try Chebyshev points
        xk = chebpts(N + 2, f.domain([1, end]));
        [p,err,status] = minimaxKernel(f, fHandle,m, n, N, rationalMode,...
                                       xk, opts, 1);
        q = chebfun('1', f.domain([1, end]));
        if ( polyOutput )
            varargout = {p, err, status};
        else
            varargout = {p, q, p, err, status};
        end
    else % rational case
        % A first attempt is using CF as an initial guess
        try
            xk = cfInit(f, fHandle, m, n);
            [p,q,rh,err,status] = minimaxKernel(f, fHandle,m, n, N, ...
                                                rationalMode, xk, opts, 1);
            varargout = {p, q, rh, err, status};
        catch
            if ~opts.silentFlag
                disp(['CF-based initialization failed,' ...
                    ' turning to AAA-Lawson...']);
            end
            status.success = 0; % CF didn't work
        end

        % If CF doesn't give a satisfactory answer, we try AAA-Lawson
        if ~status.success
            if ~opts.silentFlag
                disp('Trying AAA-Lawson-based initialization...');
            end
            xk = AAALawsonInit(f, fHandle, m, n);
            [p,q,rh,err,status] = minimaxKernel(f, fHandle, m, n, N, ...
                                                rationalMode, xk, opts, 1);
            varargout = {p, q, rh, err, status};
        end

        % A final attempt using cumulative distribution functions
        if ~status.success
            xk = cdfInit(f, fHandle, m, n, symFlag, opts, 1);
            [p,q,rh,err,status] = minimaxKernel(f, fHandle, m, n, N, ...
                                                rationalMode, xk, opts, 1);
            varargout = {p, q, rh, err, status};
        end

        if ~status.success
            xk = cdfInit(f, fHandle, m, n, symFlag, opts, 2);
            [p,q,rh,err,status] = minimaxKernel(f, fHandle, m, n, N, ...
                                                rationalMode, xk, opts, 1);
            varargout = {p, q, rh, err, status};
        end

        if ~status.success % all attempts failed
        	error('CHEBFUN:CHEBFUN:minimax:failure', ...
               ['MINIMAX failed to produce the best approximant. ' ...
                'If the accuracy is close to machine precision, ' ...
                'it may be that what you''ve asked for is unachievable ' ...
                'in floating-point arithmetic. In such a case, try ' ...
                'reducing the degree to get a clean best approximant.'])
        end
    end
else  % the user has also given a starting reference
    if ( n == 0 )
        [p,err,status] = minimaxKernel(f, fHandle, m, n, N, ...
                                       rationalMode, xk, opts, 1);
        q = chebfun('1', f.domain([1,end]));
        if ( polyOutput )
            varargout = {p, err, status};
        else
            varargout = {p, q, p, err, status};
        end
    else
        [p,q,rh,err,status] = minimaxKernel(f, fHandle, m, n, N, ...
                                            rationalMode, xk, opts, 1);
        varargout = {p, q, rh, err, status};
    end
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Functions implementing the different initialization strategies

% CF-based initialization
function xk = cfInit(f, fHandle, m, n)
    warning off
    if ( numel(f.funs) == 1 )
        [p, q] = cf(f, m, n);
        pqh = @(x) feval(p, x)./feval(q, x);
        [xk, ~, ~, flag] = exchange([], 0, 2, f, fHandle, p, pqh, m+n+2, n);
    else
        try
            [p, q] = cf(f, m, n, 50*(m+n) );
            pqh = @(x) feval(p, x)./feval(q, x);
            [xk, ~, ~, flag] = exchange([], 0, 2, f, fHandle, p, pqh, ...
                                        m+n+2, n);
        catch ME  % an error occured when calling cf (ignore it)
            flag = 0;
        end
    end
    warning on

    % If the above procedure failed to produce a reference
    % with enough oscillation points, use polynomial Remez.
    if ( flag == 0 )
        [~,~,status] = minimax(f, m+n); xk = status.xk;
    end
end

% Now turn to initialization via AAA-Lawson, this is more expensive than CF
% but less so than CDF (which follows if this fails).
%
function xk = AAALawsonInit(f,fHandle, m,n) % AAA-Lawson initialization for
                                            % functions with breakpoints
    NN = max(10*max(m,n),round(1e5/max(m,n)));
    dom = domain(f);
    Z = linspace(dom(1), dom(end), NN);
    F = fHandle(Z);
    [r,~,~,~,xk] = aaamn_lawson(F, Z, m, n);    % 1st AAA-Lawson
    xk = findReference(f, fHandle, r, m, n, xk);

    % Iterate twice during which sample points are refined.
    % This is done to find the nonsmooth parts of the function
    % and sample more densely there.
    for it = 1:2 % (maybe helps to run more)
    num = round(NN/length(xk));
    Z = [];
    for ii = 1:length(xk)-1
        Z = [Z linspace(xk(ii), xk(ii+1),num)]; % equispaced sampling
                                                % between each pair of
                                                % reference pts
    end
    Z = unique(Z); Z = Z(:); F = feval(f, Z);
    [r,~,~,~,xk] = aaamn_lawson(F, Z, m, n); % Do AAA-Lawson with updated
                                             % sample pts
    xk = findReference(f, fHandle, r, m, n, xk);
    end
end

% Cumulative distribution function-based initialization of the rational
% version of the exchange algorithm.
function xk = cdfInit(f, fHandle, m, n, symFlag, opts, step)
    newlineCounter = 0;
    stepSize = step; % Increase in degree of numerator and/or denominator
                     % at each step.
    if ( symFlag > 0 ) % Dealing with symmetry (even or odd).
        stepSize = 2*step;
    end
    if ~opts.silentFlag
        text = ['Trying CDF-based initialization with step size ', ...
                                                num2str(stepSize),'...'];
        disp(text);
    end

    % The approximation is close to being diagonal; start from an
    % approximation where both the numerator and denominator degrees
    % are decreased by the same value.
    if ( abs(m-n) <= 2 )
        minValue = min(m, n);
        minValue = minValue - rem(minValue, stepSize);
        k = minValue / stepSize;
        k = k - (3 - step);
        % Starting small degree problem.
        startM = m - stepSize * k;
        startN = n - stepSize * k;

        % We need an initialization strategy that has a high chance
        % of working without problem for the small degree case;
        % CF is used for now.
        xk = cfInit(f, fHandle, startM, startN);
        [~,~,~,~,status] = minimaxKernel(f, fHandle, startM, startN, ...
                                        startM+startN, true, xk, opts, 0);
        % If Remez worked on the small degree problem, start increasing
        % the degrees in both the numerator and denominator.
        if status.success
            while(startM < m - stepSize && (status.success == 1))
                startM = startM + stepSize;
                startN = startN + stepSize;
                newlineCounter = newlineCounter + 1;
                if(newlineCounter == 10)
                    newlineCounter = 0;
                    if ~opts.silentFlag
                        fprintf('\n');
                    end
                end
                if ~opts.silentFlag
                    fprintf('(%d,%d) ',startM, startN);
                end
                % Use the distribution information from the previous
                % reference to construct a starting reference for the new,
                % larger degree problem.
                xk = refGen(f, status.xk, startM + startN + 2, symFlag);
                [~,~,~,~,status] = minimaxKernel(f, fHandle, startM, ...
                                    startN, startM+startN, true, xk, ...
                                    opts, 0);
            end
        end

        if status.success
            xk = refGen(f, status.xk, m + n + 2, symFlag);
            if ~opts.silentFlag
                fprintf('(%d,%d)\n',m, n);
            end
        else
            % There was a failure somewhere, use reference from polynomial
            % Remez (might work sometimes).
            if ~opts.silentFlag
                fprintf('\n');
                text = ['Initialization failed using CDF with step', ...
                                    ' size ', num2str(stepSize) '...'];
                disp(text);
            end
            [~,~,status] = minimax(f, m+n); xk = status.xk;
        end
    else
        % Similar strategy to the diagonal case.
        if ( m < n )
            % Construct the 'corner' instance.
            % (m - stepSize*k, n - stepSize*k), where m - stepSize * k will
            % usually be 0 or 1.
            minValue = m - rem(m, stepSize);
            k = minValue / stepSize;
            hM = m - stepSize * k;
            hN = n - stepSize * k;

            % Now decrease the degree only in the denominator.
            hminValue = hN - rem(hN, stepSize);
            hk = hminValue / stepSize;
            startN = hN - stepSize * hk;
            xk = cfInit(f, fHandle, hM, startN);
            [~,~,~,~,status] = minimaxKernel(f, fHandle, hM, startN, ...
                                        hM+startN, true, xk, opts, 0);
            % Construct approximations by successively increasing the
            % denominator degree.
            if status.success
                while ( startN < hN - stepSize && status.success )
                    startN = startN + stepSize;
                    xk = refGen(f, status.xk, hM + startN + 2, 0);
                    newlineCounter = newlineCounter + 1;
                    if ( newlineCounter == 10 )
                        newlineCounter = 0;
                        if ~opts.silentFlag
                            fprintf('\n');
                        end
                    end
                    if ~opts.silentFlag
                        fprintf('(%d,%d) ',hM, startN);
                    end
                    [~,~,~,~,status] = minimaxKernel(f, fHandle, hM, ...
                                        startN, hM+startN, true, xk, ...
                                        opts, 0);
                end
            end

            if ~status.success
                [~,~,status] = minimax(f, hM + hN);
            end

            % Go to the initial degree by now simultaneously increasing
            % both numerator and denominator degree.
            status.success = 1;
            while ( hM < m - stepSize && status.success )
                    hM = hM + stepSize;
                    hN = hN + stepSize;
                    xk = refGen(f, status.xk, hM + hN + 2, symFlag);
                    newlineCounter = newlineCounter + 1;
                    if(newlineCounter == 10)
                        newlineCounter = 0;
                        if ~opts.silentFlag
                            fprintf('\n');
                        end
                    end
                    if ~opts.silentFlag
                        fprintf('(%d,%d) ', hM, hN);
                    end
                    [~,~,~,~,status] = minimaxKernel(f, fHandle, hM, ...
                                                    hN, hM+hN, true, ...
                                                    xk, opts, 0);
            end

            if status.success
                xk = refGen(f, status.xk, m + n + 2, 0);
                if ~opts.silentFlag
                    fprintf('(%d,%d)\n', m, n);
                end
            else
                if ~opts.silentFlag
                    fprintf('\n');
                    text = ['Initialization failed using CDF with ', ...
                                'step size ', num2str(stepSize)];
                    disp(text);
                end
                [~,~,status] = minimax(f, m+n); xk = status.xk;
            end

        else % m > n
            % Construction of the 'corner' instance by decreasing the
            % degree in both numerator and denominator.
            minValue = n - rem(n, stepSize);
            k = minValue / stepSize;
            startM = m - stepSize * k;
            startN = n - stepSize * k;
            xk = cfInit(f, fHandle, startM, startN);
            [~,~,~,~,status] = minimaxKernel(f, fHandle, startM, ...
                                    startN, startM + startN, true, ...
                                    xk, opts, 0);

            if status.success
                while ( startM < m - stepSize && status.success )
                    startM = startM + stepSize;
                    startN = startN + stepSize;
                    newlineCounter = newlineCounter + 1;
                    if ( newlineCounter == 10 )
                        newlineCounter = 0;
                        if ~opts.silentFlag
                            fprintf('\n');
                        end
                    end
                    if ~opts.silentFlag
                        fprintf('(%d,%d) ', startM, startN);
                    end
                    xk = refGen(f, status.xk, startM + startN + 2, ...
                                    symFlag);
                    [~,~,~,~,status] = minimaxKernel(f, fHandle, ...
                                        startM, startN, startM+startN, ...
                                        true, xk, opts, 0);
                end
            end

            if status.success
                xk = refGen(f, status.xk, m + n + 2, symFlag);
                if ~opts.silentFlag
                    fprintf('(%d,%d)\n', m, n);
                end
            else
                if ~opts.silentFlag
                    fprintf('\n');
                    text = ['Initialization failed using CDF with ', ...
                            'step size ', num2str(stepSize)];
                    disp(text);
                end
                [~,~,status] = minimax(f, m+n); xk = status.xk;
            end
        end
    end
    if ~opts.silentFlag
        fprintf('\n');
    end
end

function varargout = minimaxKernel(f, fHandle, m, n, N, rationalMode, ...
                                   xk, opts, dialogFlag)

% This core function should only ever be called with a nonempty initial set
% of xk reference values
normf = opts.normf;
dom = opts.dom;

% If m = -1, this means f = odd and input (m,n) = (0,n); return constant 0.
if ( m == -1 && dialogFlag)
    q = chebfun(1, dom);
    p = chebfun(0, dom);
    varargout = {p, q, @(x) feval(p, x)./feval(q, x), norm(f,inf), []};
    return
elseif ( m == -1 )
    q = [];
    p = [];
    varargout = {p, q, [], [], []};
    return
end

% With zero denominator degree, the denominator polynomial is trivial.
if ( n == 0 )
    q = chebfun(1, dom);
end

% Initial values for some parameters.
iter = 0;                 % Iteration count.
delta = max(normf, eps);  % Value for stopping criterion.
deltamin = inf;           % Minimum error encountered.
diffx = 1;                % Maximum correction to trial reference

xo = xk;

% Print header for text output display if requested.
if ( opts.displayIter && dialogFlag)
    disp('It.   Max(|Error|)     |ErrorRef|    Delta ErrorRef    Delta Ref     m  n')
end

err = normf;
% Initialise the levelled error such that one iteration always executes
h = 2*err + 1;
interpSuccess = 1;
% Run the main algorithm.
while ( (abs(abs(h)-abs(err))/abs(err) > opts.tol) && ...
    (iter < opts.maxIter) && (diffx > 0) && interpSuccess )
    hpre = h;
    % Approximation error is at the level of machine precision, stop.
    if ( abs(abs(h)-abs(err))/normf < 1e-14 )
        break
    end
    % Compute trial function and levelled reference error.
    if ( n == 0 )
        fk = fHandle(xk);      % Evaluate on the exchange set.
        w = baryWeights(xk);   % Barycentric weights for exchange set.
        [p, h] = computeTrialFunctionPolynomial(fk, xk, w, m, N, dom);

        % Perturb exactly-zero values of the levelled error.
        if ( h == 0 )
            h = 1e-19;
        end

        rh = @(x) 0;
        % Update the exchange set using the Remez algorithm
        % with full exchange rule.
        [xk, err, err_handle, ~] = exchange(xk, h, 2, f, fHandle, p, ...
                                            rh, N + 2, n);

        % If overshoot, recompute with one-point exchange rule.
        if ( err/normf > 1e5 )
            [xk, err, err_handle, ~] = exchange(xo, h, 1, f, fHandle, ...
                                            p, rh, N + 2, n);
        end

        % Update max. correction to trial reference and stopping criterion.
        diffx = max(abs(xo - xk));
        delta = err - abs(h);

        % Store approximation with minimum norm.
        if ( delta < deltamin )
            pmin = p;
            errmin = err;
            xkmin = xk;
            deltamin = delta;
        end

    else
        err = inf;
        [p, q, rh, h, interpSuccess, tk, alpha, beta] = ...
            computeTrialFunctionRational(f, fHandle, xk, m, n, hpre, ...
                                         dialogFlag, opts.silentFlag);

        % Perturb exactly-zero values of the levelled error.
        if ( h == 0 )
            h = 1e-19;
        end

        if(interpSuccess == 1)
            [xk, err, err_handle, ~] = exchange(xk, h, 2, f, fHandle, ...
                                            p, rh, N+2, n);
            diffx = max(abs(xo - xk));
            delta = err - abs(h);

            if opts.tol*norm(err_handle(xk),inf) < normf*1e-14
                % Relative tolerance below machine precision, make it
                % reasonable.
                opts.tol = normf*1e-13/norm(err_handle(xk),inf);
                opts.tol = min(opts.tol, 0.1);
            end
        end
    end

    % Display diagnostic information as requested.
    if ( opts.plotIter && interpSuccess && dialogFlag )
        doPlotIter(xo, xk, err_handle, h, dom);
    end

    if ( opts.displayIter && dialogFlag )
        doDisplayIter(iter, err, h, delta, normf, diffx, m, n);
    end

    xo = xk;
    iter = iter + 1;
end

if ( n == 0 )
    % Take best results of all the iterations we ran.
    p = pmin;
    err = errmin;
    xk = xkmin;
    delta = deltamin;
end

% Warn the user if we failed to converge.
if ( abs(abs(h)-abs(err))/abs(err) > opts.tol && ...
     abs(abs(h)-abs(err))/normf >= 1e-14 && dialogFlag && interpSuccess )
    warning('CHEBFUN:CHEBFUN:minimax:convergence', ...
        ['minimax algorithm did not converge after ', num2str(iter), ...
         ' iterations to the tolerance ', num2str(opts.tol), '.']);
end

% Form the outputs.
status.delta = delta/normf;
status.iter = iter;
status.diffx = diffx;
status.xk = xk;
status.success = interpSuccess;
% Compute the poles and zeros in case of a rational approximation
if status.success && dialogFlag && rationalMode
    [status.zer, status.pol] = pzeros(tk, alpha, beta, rh, m, n, dom);
else
    status.zer = []; status.pol = [];
end

if( ~isempty(p))
    p = simplify(p);
end
if ( rationalMode )
    varargout = {p, q, rh, err, status};
else
    varargout = {p, err, status};
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input parsing.

function [m, n, N, rationalMode, polyOutput, symFlag, xk, opts] = ...
                                             parseInputs(f, varargin)

opts.silentFlag = 0;
isSilent = 0;
for k = 1:length(varargin)
    if ( ischar(varargin{k}) && strcmpi('silent', varargin{k}) )
        opts.silentFlag = k;
        isSilent = 1;
    end
end

if opts.silentFlag
    varargin(opts.silentFlag) = [];
end

% Detect polynomial / rational approximation type and parse degrees.
polyOutput = true;
if ( ~mod(nargin - isSilent, 2) ) % Even number of inputs --> polynomial
                                  % case.
    m = varargin{1};
    n = 0;
    rationalMode = false;
    symFlag = 0;
    varargin = varargin(2:end);
else                              % Odd number of inputs --> rational case.
    polyOutput = false;
    [m, n, symFlag] = adjustDegreesForSymmetries(f, varargin{1}, ...
                                                 varargin{2});
    if ( n == 0 )
        rationalMode = false;
    else
        rationalMode = true;
    end
    varargin = varargin(3:end);
end

N = m + n;

% Parse name-value option pairs.
if rationalMode
    opts.tol = 1e-4;                     % Relative tolerance for deciding
                                         % convergence.
    opts.maxIter = 10+round(max(m,n)/2); % Maximum number of allowable
                                         % iterations.
else
    opts.tol = 1e-14*(N^2 + 10); % Polynomial case is much more robust.
    opts.maxIter = 30;           % Maximum number of allowable iterations.
end

opts.displayIter = false;    % Print output after each iteration.
opts.plotIter = false;       % Plot approximation at each iteration.
opts.dom = f.domain([1, end]);
opts.normf = norm(f);
xk = [];

for k = 1:2:length(varargin)
    if ( strcmpi('tol', varargin{k}) )
        opts.tol = varargin{k+1};
    elseif ( strcmpi('maxiter', varargin{k}) )
        opts.maxIter = varargin{k+1};
    elseif ( strcmpi('display', varargin{k}) )
        opts.displayIter = true;
    elseif ( strcmpi('plotfcns', varargin{k}) )
        opts.plotIter = true;
    elseif ( strcmpi('plot', varargin{k}) )
        opts.plotIter = true;
    elseif ( strcmpi('init', varargin{k}) )
        xk = varargin{k+1};
    else
        error('CHEBFUN:CHEBFUN:minimax:badInput', ...
            'Unrecognized sequence of input parameters.')
    end
end

end

function [m, n, symFlag] = adjustDegreesForSymmetries(f, m, n)
%ADJUSTDEGREESFORSYMMETRIES   Adjust rational approximation degrees to
%   account for function symmetries.
%
%   [M, N] = ADJUSTDEGREESFORSYMMETRIES(F, M, N) returns new degrees M and
%   N to correct the defect of the rational approximation if the target
%   function is even or odd.  In either case, the Walsh table is covered
%   with blocks of size 2x2, e.g.  for even function the best rational
%   approximant is the same for types [m/n], [m+1/n], [m/n+1] and
%   [m+1/n+1], with m and n even. This strategy is similar to the one
%   proposed by van Deun and Trefethen for CF approximation in Chebfun
%   (see @chebfun/cf.m).

% Sample piecewise-smooth CHEBFUNs.
if ( (numel(f.funs) > 1) || (length(f) > 128) )
  f = chebfun(f, f.domain([1, end]), 128);
end

% Assume no symmetry at the outset.
symFlag = 0;
% Compute the Chebyshev coefficients.
c = chebcoeffs(f, length(f));
c(1) = 2*c(1);

% Check for symmetries and reduce degrees accordingly.
if ( max(abs(c(2:2:end)))/vscale(f) < eps )     % f is even.
    symFlag = 1;
    if ( mod(m, 2) == 1 )
        m = max(m - 1, 0);
    end
    if ( mod(n, 2) == 1 )
        n = max(n - 1, 0);
    end
elseif ( max(abs(c(1:2:end)))/vscale(f) < eps ) % f is odd.
    symFlag = 2;
    if ( ~mod(m, 2) )
        m = m - 1;
    end
    if ( mod(n, 2) )
        n = n - 1;
    end
end

end

function [p, h] = computeTrialFunctionPolynomial(fk, xk, w, m, N, dom)

% Vector of alternating signs.
sigma = ones(N + 2, 1);
sigma(2:2:end) = -1;

h = (w'*fk) / (w'*sigma);           % Levelled reference error.
pk = (fk - h*sigma);                % Vals. of r*q in reference.

% Trial polynomial.
p = chebfun(@(x) bary(x, pk, xk, w), dom, m + 1);

end

function [p, q, rh, h, interpSuccess,xsupport, wN, wD] = ...
    computeTrialFunctionRational(f, fHandle, xk, m, n, hpre, ...
                                 dialogFlag, silentFlag)
% computeTrialFunctionRational finds a rational approximation to f at an
% iteration of the Remez algorithm. It uses a barycentric representation
% for improved numerical stability.
% f:          function
% xk:         approximate reference points
% m,n:        type
% hpre:       levelled error at the previous iteration

% The function values at the current reference points
fk = fHandle(xk);
% Take barycentric support points to be alternating values of two
% reference points
    xsupport = xk(2:2:end);
    xsuppind = 2:2:length(xk);
    xadd = xk(1:2:end);
    xother = xadd;
    xotherind = 1:2:length(xk);

if m~=n % need to add more support points
    xadd = xother;
    % take Leja points from the remaining ref pts and add
    [xadd, ~] = leja(xadd, 1, length(xadd));
    xsupport = [xsupport;xadd(1:max(m,n)+1-length(xsupport))];
    xother = zeros(m+n+2-max(m,n)-1,1);
    xotherind = zeros(m+n+2-max(m,n)-1,1);
    xsuppind = zeros(max(m,n)+1,1);

    iother = 1; isupp = 1;
    for ii = 1:length(xk) % keep indices for later use
        if ~ismember(xk(ii),xsupport)
            xother(iother) = xk(ii);
            xotherind(iother) = ii;
            iother = iother+1;
        else
            xsuppind(isupp) = ii;
            isupp = isupp+1;
        end
    end

end
    xsupport = sort(xsupport,'ascend');

if m~=n
    % projection matrices that force coefficients to lie in null space
    % of Vandermonde matrix
    % projection subspace
    Qmn = orthspace(xsupport,abs(m-n),ones(length(xsupport),1));
    [Qmnall,~] = qr(Qmn);
end

    C = 1./bsxfun(@minus,xk,xsupport.');    % Cauchy matrix

    % form matrix Cstar = sqrt(|Delta|)*C
    % Cstar(ii,jj) = |wt(xi)/sqrt(wx'(xi))|/(xi-tj)
    Xkdiff = abs(bsxfun(@minus, xk, xk.'));
    Xkdiff(eye(size(Xkdiff))~=0) = 1;       % inf to 0
    Xtdiff = abs(bsxfun(@minus, xother, xsupport.'));
    ST = sum(log(Xtdiff.')); SX = sum(log(Xkdiff));
    VV = exp(ST.'-0.5*SX(xotherind).');
    VV = VV*ones(1,length(xsupport));
    Div = bsxfun(@minus,xother,xsupport.');
    C1 = VV./Div; % odd columns of Cstar

    % diag elements Cstar(ii,jj) = |wt(xi)/sqrt(wx'(xi))|/(xi-tj)
    Xtdiff = abs(bsxfun(@minus, xsupport, xsupport.'));
    Xtdiff(eye(size(Xtdiff))~=0) = 1;
    ST = sum(log(Xtdiff.'));SX = sum(log(Xkdiff));
    C2 = diag(exp(ST.'-0.5*SX(xsuppind).'));
    Cstar = [C1;C2];
    Cstar(xsuppind,:) = C2;
    Cstar(xotherind,:) = C1;

% prepare QR factorizations; these lead to a symmetric eigenproblem
% we need to be careful how to do QR as rows have
% large dynamical range (though orthogonal columns when m=n)
% do Householder QR with row sorting, better than [Q,R] = qr(Cstar,0);
if ( m == n )
    nrm = zeros(1, size(Cstar, 1));
    for ii = 1:length(nrm)
        nrm(ii) = norm(Cstar(ii,:));
    end
    [~,ix] = sort(nrm, 'descend');
    %[~,ix] = sort(norms(Cstar'), 'descend');
    [Q,R] = qr(Cstar(ix,:),0);
    ixx(ix) = 1:length(ix);    Q = Q(ixx,:);

    %{
    nrm = zeros(1,m+1); % Cholesky QR with col-scaling, this works too
    for ii = 1:m+1, nrm(ii) = norm(Cstar(:,ii));    end
    Cstar = Cstar/diag(nrm);    % this normalized Cstar is orthogonal
    CTC = Cstar'*Cstar;
    R = chol(CTC);    Q = Cstar/R;    R = R*diag(nrm);
    %}
elseif ( m > n )
    nrm = zeros(1, size(Cstar, 1));
    for ii = 1:length(nrm)
        nrm(ii) = norm(Cstar(ii,:));
    end
    [~,ix] = sort(nrm, 'descend');
    %[~,ix] = sort(norms(Cstar'),'descend');
    [Q,R] = qr(Cstar(ix,:)*Qmn,0);
    [Qall,Rall] = qr(Cstar(ix,:)*Qmnall,0);
    ixx(ix) = 1:length(ix);
    Q = Q(ixx,:);        Qall = Qall(ixx,:);

else % m<n
    nrm = zeros(1, size(Cstar, 1));
    for ii = 1:length(nrm)
        nrm(ii) = norm(Cstar(ii,:));
    end
    [~,ix] = sort(nrm, 'descend');
    %[~,ix] = sort(norms(Cstar'),'descend');
    [Q,R] = qr(Cstar(ix,:),0);
    [Qpart,Rpart] = qr(Cstar(ix,:)*Qmn,0);
    ixx(ix) = 1:length(ix);
    Q = Q(ixx,:);        Qpart = Qpart(ixx,:);
end

S = diag((-1).^(0:length(xk)-1));
%Q2 = S*Q; for sanity check svd(Q'*Q2) or svd(Qpart'*Q2) when m<n,
% should be O(eps)

QSQ = Q'*S*diag(fk)*Q;
QSQ = (QSQ+QSQ')/2; % force symmetry as it's supposed to be

% key operation; this forces (F+hsigma)N=D, where N/D is rational
% approximant. The eigenvector VR containing the coefficients for
% D(x)= sum_i VR_{i}/(x-xsupport_{i}).
[VR,d] = eig(-QSQ); % symmetric eigenproblem
beta = R\VR;        % Denominator coefficients in barycentric form

% obtain alpha (the Numerator coefficients in bary form) from beta
if ( m == n )
    alpha = R\(-Q'*diag(fk)*Q*VR);
elseif ( m > n )
    alpha = Qmnall*((Rall)\((Qall'*diag(-fk)*Q*VR)));
else % m<n
    alpha = (Rpart)\(Qpart'*diag(-fk)*Q*VR);
end
vt = [alpha;beta];

% Among the n+1 eigenvalues, only one can be the solution. The correct
% one needs to have no sign changes in the denominator polynomial
% D(x)*node(x), where node(x) = prod(x-xsupport).

if ( m <= n ) % values of D at xk
    % Dvals = C(:,1:n+1)*vt(m+1+1:end,:);
    bet = vt(m+1+1:end,:);
else
    bet = Qmn*vt(m+1+1:end,:);
end
    Dvals = C*bet;

node = @(z) prod(z-xsupport); % needed to check sign
nodevec = xother;
for ii = 1:length(xother)
    nodevec(ii) = node(xother(ii));   % values of node polynomial
end
% Find position without sign changes in D*node.
% Evaluate this separately for xsupport and xother.
% Ignore ones with too small Dvals.

checksign = zeros(length(xk),n+1);
% sign at supp pts
checksign(1:length(xsupport),:) = ...
    diag((-1).^(max(m,n):max(m,n)+length(xsupport)-1))*bet;
% sign at other pts
signs = sign(diag(nodevec)*Dvals(xotherind,:));
checksign(length(xsupport)+1:end,:) = signs;
pos = find(abs(sum(sign(checksign))) == m+n+2 & sum(abs(Dvals))>1e-7);

if isempty(pos)  % Unfortunately, no solution with same signs.
                 % Try old remez.

% Take barycentric support points to be alternating values of two
% reference points
xsupport = (xk(1:2:end-1)+xk(2:2:end))/2;
xadd = (xk(2:2:end-1)+xk(3:2:end))/2; % when m~=n, we need more support
                                      % points

if ismember(f.domain(1),xk) == 0      % if endpoints aren't included,
                                      % add them
    xadd = [(f.domain(1)+xk(1))/2;xadd];
end
if ismember(f.domain(end),xk) == 0
    xadd = [(f.domain(end)+xk(end))/2;xadd];
end
num = abs((max(m,n)+1-length(xsupport)));
[xadd, ~] = leja(xadd, 1, num);  % take Leja points from the
                                 % remaining ref pts

if m~=n
    % add any lacking supp pts
    xsupport = [xsupport;xadd(1:max(m,n)+1-length(xsupport))];
end
xsupport = sort(xsupport, 'ascend');

C = 1./bsxfun(@minus,xk,xsupport.');

% find Delta diag matrix
Delta = zeros( 1,length(xk) );
for ii = 1:length(xk)
% wt(ii) = prod(xk(ii)-xsupport);
% wxdiff(ii) = prod(xk(ii)-xk([1:ii-1 ii+1:end]));
% Delta = diag(-(wt.^2)./wxdiff); do in a way that avoids
%                                 underflow, overflow
    Delta(ii) = -exp(2*sum(log(abs(prod(xk(ii)-xsupport)))) ...
        - sum(log(abs(xk(ii)-xk([1:ii-1 ii+1:end])))));
end
Delta = diag(Delta);

%DD = diag(1./norms(sqrt(abs(Delta))*C)); % scaling, might help stability
DD = eye(size(C,2));

% prepare QR factorizations; these lead to symmetric eigenproblem
if ( m == n )
    [Q,R] = qr(sqrt(abs(Delta))*C,0);
elseif ( m > n )
    [Q,R] = qr(sqrt(abs(Delta))*C*Qmn,0);
    [Qall,Rall] = qr(sqrt(abs(Delta))*C*Qmnall,0);
else % m<n
    [Q,R] = qr(sqrt(abs(Delta))*C,0);
    [Qpart,Rpart] = qr(sqrt(abs(Delta))*C*Qmn,0);
end

S = diag((-1).^(0:length(xk)-1));
% Q2 = S*Q; for sanity check svd(Q'*Q2) or svd(Qpart'*Q2) when m<n,
% should be O(eps)

QSQ = Q'*S*diag(fk)*Q;
QSQ = (QSQ+QSQ')/2; % force symmetry as it's supposed to be

% key operation; this forces (F+hsigma)N=D, where N/D is rational
% approximant. The eigenvector VR for which beta=R\VR contains the
% coefficients for D(x)= sum_i beta_{i}/(x-xsupport_{i}).
[VR,d] = eig(-QSQ); % symmetric eigenproblem
beta = R\VR;        % Denominator coefficients in barycentric form

% obtain alpha (the Numerator coefficients in bary form) from beta
if ( m == n )
    alpha = R\(-Q'*diag(fk)*Q*VR);
elseif ( m > n )
    alpha = Qmnall*((Rall)\((Qall'*diag(-fk)*Q*VR)));
else % m<n
    alpha = (Rpart)\(Qpart'*diag(-fk)*Q*VR);
end
vt = [alpha;beta];

% conditioning check, might help
% disp([cond(C) cond(sqrt(abs(Delta))*C) cond(sqrt(abs(Delta))) ...
%     cond(C*DD) cond(sqrt(abs(Delta))*C*DD) m n])

if ( m <= n ) % values of D at xk
    Dvals = C(:,1:n+1)*(DD*vt(m+1+1:end,:));
else
    Dvals = C*(Qmn*vt(m+1+1:end,:));
end
node = @(z) prod(z-xsupport);     % needed to check sign

nodevec = xk;
for ii = 1:length(xk)
    nodevec(ii) = node(xk(ii));   % values of node polynomial
end
% Find position without sign changes in D*node.
pos = find(abs(sum(sign(diag(nodevec)*Dvals))) == m+n+2 & ...
                                        sum(abs(Dvals))>1e-4);

if isempty(pos) % still no solution, give up
    if ( dialogFlag && ~silentFlag )
        disp('Trial interpolant too far from optimal...')
    end
    interpSuccess = 0;
    p = []; q = []; rh = []; h = 1e-19; wD = []; wN = [];
    return
elseif ( length(pos) > 1 ) % more than one solution with no sign changes...
    [~,ix] = min(abs(hpre)-diag(abs(d(pos,pos))));
    pos = pos(ix);
end

end

h = -d(pos, pos);                 % levelled reference error.

% coefficients for barycentric representations
if ( m <= n )
    wD = vt(m+2:end,pos);
else
    wD = Qmn*vt(m+2:end,pos);
end
if ( m >= n )
    wN = vt(1:m+1,pos);
else
    wN = Qmn*vt(1:m+1,pos);
end

D = @(x) 0; N = @(x) 0;    % form function handle rh = N/D
for ii = 1:length(xsupport)
   D = @(x) D(x) + wD(ii)./(x-xsupport(ii));
   N = @(x) N(x) + wN(ii)./(x-xsupport(ii));
end
D = @(x)-D(x); % flip back sign

rh = @(zz) reval(zz, xsupport, N, D, wN, wD);

interpSuccess = 1; % declare success

% Form chebfuns of p and q (note: could be numerically unstable, but
% provided for convenience).
% Find values of node polynomial at Chebyshev points
if dialogFlag
    x = chebpts(m+n+1,f.domain([1,end]));
    nodex = zeros(length(x),1);
    for ii = 1:length(x)
        nodex(ii) = node(x(ii));
    end
    qvals = nodex.*feval(D,x);  % Values of p and q at Chebyshev points
    pvals = nodex.*feval(N,x);
    % If certain Chebyshev points map to support points, inf values will
    % get propagated, so we need to handle them separately
    for ii = 1:length(xsupport)
        for jj = 1:length(x)
            if x(jj) == xsupport(ii)
                nodei = 1.0;
                for kk = 1:length(xsupport)
                    if ~(kk == ii)
                        nodei = nodei * (x(jj) - xsupport(kk));
                    end
                end
                qvals(jj) = -nodei*wD(ii);
                pvals(jj) = nodei*wN(ii);
            end
        end
    end

    p = chebfun(pvals,f.domain([1,end]));
    q = chebfun(qvals,f.domain([1,end]));
    p = simplify(p); q = simplify(q);
else
    p = [];
    q = [];
end

end

function r = reval(zz, xsupport, N, D, wN, wD)
zv = zz(:);
r = N(zv)./D(zv);
ii = find(isnan(r));
for jj = 1:length(ii)
    if ( isnan(zv(ii(jj))) || ~any(zv(ii(jj)) == xsupport) )
        % r(NaN) = NaN is fine.
        % The second case may happen if r(zv(ii)) = 0/0 at some point.
    else
        % Clean up values NaN = inf/inf at support points.
        % Find the corresponding node and set entry to correct value:
        pos = zv(ii(jj)) == xsupport;
        r(ii(jj)) = -wN(pos)./wD(pos);
    end
end
r = reshape(r, size(zz));
end

function [xx, pos] = leja(x, startIndex, nPts)
% put NPTS from X in a Leja sequence
% starting from x(startIndex)
n = length(x);
p = zeros(n,1);
pos = zeros(nPts, 1);
xx = zeros(nPts, 1);
xx(1) = x(startIndex);
pos(1) = startIndex;

for j = 2:nPts
    % we want to pick the jth point now:
    for i = 1:n
        %p(i) = prod(abs(x(i) - xx(1:j-1)));
        p(i) = sum(log(abs(x(i) - xx(1:j-1)))); % no overflow
    end
    [~,pos(j)] = max(p);
    xx(j) = x(pos(j));
end

end

function [xk, norme, err_handle, flag] = exchange(xk, h, method, f, ...
                                                  fHandle, p, rh, Npts, n)
%EXCHANGE   Modify an equioscillation reference using the Remez algorithm.
%   EXCHANGE(XK, H, METHOD, F, P, RH, NPTS, N) performs one step of the
%   Remez algorithm for the best rational approximation of the CHEBFUN F
%   of the target function according to the first method (METHOD = 1),
%   i.e., exchanges only one point, or the second method (METHOD = 2),
%   i.e., exchanges all the reference points. XK is a column vector with
%   the reference, H is the levelled error, P is the numerator, and RH is a
%   function handle, NPTS is the required number of alternation points,
%   and N is the denominator degree.
%
%   [XK, NORME, E_HANDLE, FLAG] = EXCHANGE(...) returns the modified
%   reference XK, the supremum norm of the error NORME (included as an
%   output argument, since it is readily computed in EXCHANGE and is used
%   later in MINIMAX), a function handle E_HANDLE for the error, and a FLAG
%   indicating whether there were at least N+2 alternating extrema of the
%   error to form the next reference (FLAG = 1) or not (FLAG = 0).

% Compute extrema of the error.
if(n == 0) % polynomial case
    % Function handle output for evaluating the error.
    rh = @(x) feval(p,x);
    rr = findExtrema(f, fHandle, rh, xk);
    err_handle = @(x) fHandle(x) - feval(p, x);
else       % Rational case.
    rr = findExtrema(f, fHandle, rh, xk);
    err_handle = @(x) fHandle(x) - rh(x);
end

% Select exchange method.
if ( method == 1 )                             % One-point exchange.
    [~, pos] = max(abs(feval(err_handle, rr)));
    pos = pos(1);
else                                           % Full exchange.
    pos = find(abs(err_handle(rr)) >= abs(h)); % Values above levelled
                                               % error.
end

% Add extrema nearest to those which are candidates for exchange to the
% existing exchange set.
[r, m] = sort([rr(pos) ; xk]);
v = ones(Npts, 1);
v(2:2:end) = -1;
er = [feval(err_handle, rr(pos)) ; v*h];
er = er(m);

% Delete repeated points.
repeated = diff(r) == 0;
r(repeated) = [];
er(repeated) = [];

% Determine points and values to be kept for the reference set.
s = r(1);    % Points to be kept.
es = er(1);  % Values to be kept.
for i = 2:length(r)
    if ( (sign(er(i)) == sign(es(end))) && (abs(er(i)) > abs(es(end))) )
        % Given adjacent points with the same sign, keep one with largest
        % value.
        s(end) = r(i);
        es(end) = er(i);
    elseif ( sign(er(i)) ~= sign(es(end)) )
        % Keep points which alternate in sign.
        s = [s ; r(i)];    %#ok<AGROW>
        es = [es ; er(i)]; %#ok<AGROW>
    end
end

% Of the points we kept, choose n + 2 consecutive ones that include the
% maximum of the error.
[norme, index] = max(abs(es));
d = max(index - Npts + 1, 1);
if ( Npts <= length(s) )
    xk = s(d:d+Npts-1);
    flag = 1;
else
    xk = s;
    flag = 0;
end

end

function rts = findExtrema(f, fHandle, rh,xk)
% Finds all the local maxima and minima of f-rh.
% xk is the present reference
% rh is a handle to p/q

err_handle = @(x) fHandle(x) - rh(x);

doms = unique([f.domain'; xk]).';
doms = sort(doms,'ascend');

% Initial trial
if ( isempty(xk) )
sample_points = linspace(f.domain(1),f.domain(end),5000);
scale_of_error = norm(err_handle(sample_points),inf);
relTol =  1e-15 * (vscale(f)/scale_of_error);
    warning off
    ek = chebfun(@(x) err_handle(x), f.domain, 'eps', relTol, ...
                                               'splitting', 'on');
    warning on
    rts = roots(diff(ek), 'nobreaks');
else
    nn = 2^3; % sampling pts in each subinterval (try low number first)
    mid = (doms(1:end-1)+doms(2:end))/2; % midpoints
    rad = (doms(2:end)-doms(1:end-1))/2; % radius
    xx = ones(nn+1,1)*mid+cos(pi*((nn:-1:0).')/nn)*rad; % sample matrix
    valerr = feval(f,xx)-rh(xx);
    rts = zeros(5*length(xk),1);
    pos = 1;
    for k = 1:length(doms)-1
       rnow = rootsdiff(valerr(:,k),[doms(k) doms(k+1)],err_handle);
       rts(pos:pos+length(rnow)-1) = rnow; % update reference points
       pos = pos+length(rnow);
    end
    rts(pos:end) = [];
end

% Append end points of the domain.
rts = unique([f.domain' ; rts]);
rts = sort(rts,'ascend');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find extrema of error function in AAA-Lawson
function xk = findReference(f,fHandle,r,m,n,z)
    % f: function
    % r: rational approximant
    % m,n: (m,n) is the type of our rational approximant
    % z: barycentric support points
    % OUTPUT: reference points xk

    % Find extrema points as usual.
    xk = findExtrema(f,fHandle, r, sort(z,'ascend'));

    % Deal with length(xk) not equal to the desired m+n+2.
    if length(xk) > m+n+2 % Reduce reference pts becuse too many found.

        xkdiff = diff(xk);
        [~,ix] = sort(xkdiff,'descend'); % Take those with largest gaps.
        xk = [xk(1);xk(1+ix(1:m+n+1))];
        xk = sort(xk,'ascend');
    elseif length(xk) < m+n+2 % Increase # of reference pts
                              % if too few found.

        xkdiff = diff(xk);
        add = m+n+2-length(xk);  % We need to add this many reference
                                 % points.
        % Take those with largest gaps and fill midpoints.
        [~,ix] = sort(xkdiff,'descend');
        xk = [xk;(xk(ix(1:add))+xk(ix(1:add)+1))/2];
        xk = sort(xk,'ascend');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Functions for displaying diagnostic information.

% Function called when opts.plotIter is set.
function doPlotIter(xo, xk, err_handle, h, dom)

xxk = linspace(dom(1), dom(end), 10000);
plot(xo, err_handle(xo), 'or', 'MarkerSize', 4)   % Old reference.
holdState = ishold;
hold on
plot(xk, err_handle(xk), '*k', 'MarkerSize', 4)   % New reference.
plot(xxk, err_handle(xxk))               % Error function.
plot(xxk, ones(size(xxk))*h,'r');
plot(xxk, -ones(size(xxk))*h,'r');
if ( ~holdState )                        % Return to previous hold state.
    hold off
end
xlim(dom)
err = norm(err_handle(xk),'inf');
ylim(2*[-err,err])
legend('Current Ref.', 'Next Ref.', 'Error')
drawnow
end

% Function called when opts.displayIter is set.
function doDisplayIter(iter, err, h, delta, normf, diffx, m, n)

disp([num2str(iter,'%3g'), '      ', num2str(err, '%5.4e'), '      ', ...
    num2str(abs(h), '%5.4e'), '      ', ...
    num2str(delta/normf, '%5.4e'), '      ', num2str(diffx, '%5.4e'),...
    '    ', num2str(m, '%4g'),'  ', num2str(n, '%4g')])

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Functions used by the CDF-based initialization routine

% Constructs a piecewise linear function which interpolates the data
% (xd(i), yd(i)) when i goes from 1 to length(xd). It then computes
% the value of this piecewise linear approximation at the xi nodes
% (ni is the number of xi nodes)
function yi = pwiselin(xd, yd, ni, xi)

  nd = length(xd);
  xd = xd(:);
  yd = yd(:);
  xi = xi(:);

  if ( nd == 1 )
    yi(1:ni,1) = yd;
    return
  end

  [~, ~,k] = histcounts(xi, xd);

  k ( k == 0 ) = 1;
  k ( k == nd ) = nd - 1;

  t = ( xi - xd(k,1) ) ./ ( xd(k+1,1) - xd(k,1) );
  yi = ( 1 - t ) .* yd(k) + t .* yd(k+1);

  return
end

% Generate a set of n reference points which follow a distribution
% of the xk nodes. The symType flag tells us if we are dealing with
% an even or odd function, in which case the new reference nodes are
% taken to be symmetric with respect to the middle of the approximation
% domain
function nxk = refGen(f, xk, n, symType)

xx = linspace(-1,1,length(xk));

if(symType == 0)

    nxk = pwiselin(xx, xk, n, linspace(-1,1,n));

% handling of even symmetries

elseif(symType == 1)

    halfSize = length(xx)/2;
    halfn = n/2;

    if (xk(1) == f.domain(1))
        nxk = pwiselin(xx(halfSize+1:end), xk(halfSize+1:end),halfn, ...
            linspace(xx(halfSize+1),xx(end),halfn));
        nxk = [nxk; -nxk(2:end); f.domain(1)];
        nxk = sort(nxk,'ascend');
    elseif (xk(end) == f.domain(end))
        nxk = pwiselin(xx(1:halfSize), xk(1:halfSize),halfn, ...
            linspace(xx(1),xx(halfSize),halfn));
        nxk = [nxk; -nxk(1:end-1); f.domain(end)];
        nxk = sort(nxk,'ascend');
    else
        nxk = pwiselin(xx, xk, n, linspace(-1,1,n));
    end

% handling of odd symmetries
else

halfSize = (length(xx)-1) / 2;
halfn = (n-1)/2;

    if (xk(1) == f.domain(1))
        nxk = pwiselin(xx(halfSize+2:end), xk(halfSize+2:end),halfn, ...
            linspace(xx(halfSize+2),xx(end),halfn));
        nxk = [nxk; -nxk(1:end); f.domain(1)];
        nxk = sort(nxk,'ascend');
    elseif (xk(end) == f.domain(end))
        nxk = pwiselin(xx(1:halfSize+1), xk(1:halfSize+1),halfn, ...
            linspace(xx(1),xx(halfSize+1),halfn));
        nxk = [nxk; -nxk(1:end); f.domain(end)];
        nxk = sort(nxk,'ascend');
    else
        nxk = pwiselin(xx, xk, n, linspace(-1,1,n));
    end

end

end

function [r, pol, res, zer, z, Z, f, w, wf, errvec, p, q] = ...
                                                aaamn_lawson(F, varargin)
%AAAMN_Lawson   near-best rational approximation of F.
%
% R = aaamn_lawson(F) computes a rational aproximant of (default) type
% (10,10) on the default interval [-1,1]
%
% [R, POL, RES, ZER] = aaamn_lawson(...) outputs the poles, residues and
% zeros of R. The poles, zeros will approximate those of F (not well if
% R-F is not small)
%
% [R, POL, RES, ZER, z, Z, F, W,WF, ERRVEC, P, Q] = aaamn_lawson(...)
% outputs additionally the sample points Z, support points z, values
% F=f(Z), weights w and wf (see below) and AAA errvec, and p,q = chebfuns
% s.t. r= p/q (note this can be numerically unstable)
%
% [...] = aaamn_lawson(F,m,n) specifies the type to (m,n).
%
% [...] = aaamn_lawson(F,Z,m,n) also specifies the sample points Z
% (recommended).
%
% [...] = aaamn_lawson(F,Z,m,n,'plot','on') will plot the error functions
% as the Lawson iterations proceed.
%
% [...] = aaamn_lawson(F,m,n,'dom',[-1,2]) specifies the domain (this has
% no effect if Z is specified)
%
% [...] = aaamn_lawson(F,m,n,'tol',1e-5) specifies the Lawson iterate
% relative stopping criterion (here to 1e-5)
%
% [...] = aaamn_lawson(F,m,n,'iter',10) limits the Lawson maximum
% iterations to 10.
%
%
% The algorithm first finds a AAA rational approximant to F \approx r(Z),
% then attempts to refine the approximant by a Lawson process, i.e. an
% iterative reweighting.
%
% This code is designed for computing good reference points for the
% rational minimax code to follow, but can be used independently for
% constructing a rational approximation r that can be much closer than AAA
% to the best rational approximant.
%
% Input:  Z = vector of sample points
%         F = vector of data values, or a function handle
%         m, n: max type is (m,n), set to 10 if omitted
%         tol = relative tolerance tol, default: 1e-13
%         Lawsoniter: max. iteration number of Lawson updates (default 10)
%         doplot: 1 to plot error curve history (default 0)
%     R = AAAMN_Lawson(F, Z, m, n, NAME, VALUE) sets the following
%     parameters:
%           - 'tol', TOL: relative tolerance for Lawson
%           iteration (default 1e-5)
%           - 'iter', IT: maximal number of Lawson iterations
%           (default MMAX = max([5 min([20, m, n])])).
%
% Output: r = AAA-Lawson approximant to F (function handle)
%         pol,res,zer = vectors of poles, residues, zeros
%         errvec = vector of errors at each step
%         z,f,w,wf = vectors of support pts, function values, weights
%         s.t. r = N/D, N(x) = sum_i wf(i)/(x-z(i)) and
%         D(x) = sum_i w(i)/(x-z(i)).
%         p,q = chebfuns s.t. r= p/q (note this can be numerically unstable)
%
% Examples:
%    r = aaamn_lawson(@abs,10,10)
%    r = aaamn_lawson(@abs,10,10,'plot','on')
%    [r, pol, res, zer, z, f, w, wf, errvec, p, q] = aaamn_lawson(@abs,...
%                                        10,10,'plot','on','dom',[-1 2])
%    r = aaamn_lawson(@exp,4,2,'plot','on','dom',[-1 2])
%    r = aaamn_lawson(@(x)log(1.1-x),5,5,'plot','on')
%
%    f = chebfun(@(x)-1./(log(abs(x)).^2),[-.1,.1],'splitting','on');
%    [r,pol,res] = aaamn_lawson(f,linspace(-.1,.1,1e4),18,18,'plot','on')
%
%

% parse inputs
[F, Z, m, n, Lawsoniter, tolLawson, doplot, tol ] = ...
    parseInputslawson(F, varargin{:});

M = length(Z);                             % number of sample points
mmax = m+1; nmax = n+1;                    % for coding convenience
if ( (nargin < 6) || isempty(Lawsoniter) ) % number of Lawson updates
    Lawsoniter = max([5 min([20,mmax,nmax])]);
end
if ~isfloat(F), F = feval(F,Z); end        % convert function handle to
                                           % vector
 Z = Z(:); F = F(:);                       % work with column vectors
 SF = spdiags(F,0,M,M);                    % left scaling matrix
 J = 1:M;                                  % indices that are not support
                                           % pts
 z = []; f = []; C = [];                   % initializations
 errvec = []; R = mean(F);
for mn = 1:max(mmax,nmax)
  [~,j] = max(abs(F-R));                  % select next support point
  z = [z; Z(j)];                          % update set of support points
  f = [f; F(j)];                          % update set of data values
  J(J==j) = [];                           % update index vector
  C = [C 1./(Z-Z(j))];                    % next column of Cauchy matrix
  Sf = diag(f);                           % right scaling matrix
  A = SF*C - C*Sf;                        % Loewner matrix

        if ( mn > min(nmax,mmax) ) % nondiagonal case, find projection
                                   % subspace
             if mmax < nmax
             q = f(:);
             else
             q = ones(length(z),1);
             end
             Q = orthspace(z,mn-min(mmax,nmax),q);  % projection subspace
             [~,~,V] = svd(A(J,:)*Q,0);             % SVD on projected
                                                    % subspace
             w = Q*V(:,end);
        else
             [~,~,V] = svd(A(J,:),0);               % SVD, no projection
                                                    % needed
             w = V(:,mn);                           % weight vector
        end
  wf = w.*f;
  N = C*(w.*f); D = C*w;                  % numerator and denominator
  R = F; R(J) = N(J)./D(J);               % rational approximation
  err = norm(F-R,inf);
  errvec = [errvec; err];                 % max error at sample points
  if ( err < tol*norm(F,inf) ), break, end    % stop if converged
end
    r = @(zz) feval(@rrint,zz,z,w,f);            % AAA approximant as
                                                 % function handle
    Rori = R;

% now start Lawson, in this mode we leave interpolation and work with
% 'alpha-beta' mode.
          wei = ones(length(J),1);
          nrmbest = inf;
          if ( mn > min(nmax,mmax) )  % Deal with projection for m neq n
              if mn>nmax
                A =[SF*C*Q -C];
              else % need to redefine Q as not the same as AAA above
                q = ones(length(z),1);
                Q = orthspace(z,mn-min(mmax,nmax),q);     % projection
                                                          % subspace
                A =[SF*C -C*Q];
              end
          else
            A =[SF*C -C];   % diagonal case
          end

          rate = 1;         % default Lawson rate, will shrink if not
                            % converging
          nrmincreased = 0; % initialization
	      for it = 1:Lawsoniter
              weiold = wei;
              wei = wei .* power(abs(F(J)-R(J)),rate); % update Lawson
                                                       % weights
              wei = wei/sum(wei);                      % normalize
              if norm(weiold-wei)/norm(wei)< tolLawson % declare Lawson
                                                       % converged
                  break
              end
              % diagonal weight matrix
              D = spdiags(sqrt(wei),0,length(wei),length(wei));

              [~,~,V] = svd(D*A(J,:),0);  % weighted least-squares via SVD

          if ( mn > min(nmax,mmax) )      % deal with nondiagonal case
              if ( mn > nmax )
                w = Q*V(1:nmax,end); wf = V(nmax+1:end,end);
              else
                w = V(1:nmax,end); wf = Q*V(nmax+1:end,end);
              end
          else
            w = V(1:mn,end); wf = V(mn+1:2*mn,end);
          end
            f = wf./w;                         % for compatibility with
                                               % interpolatory-aaa
            N = C*wf; D = C*w;                 % numerator and denominator
            R = F; R(J) = N(J)./D(J);          % rational approximation
            err = norm(F-R,inf);
            errvec = [errvec; err];            % max error at sample points
            if ( err < nrmbest )    % adopt best so far
                nrmbest = norm(F-R,'inf');
                weibest = wei;                      % store best weight
                r = @(zz) feval(@rrab,zz,z,w,wf,f); % AAA approximant as
                                                    % function handle
            else
                nrmincreased = nrmincreased + 1;
            end
            if ( nrmincreased >= 3 )     % perhaps not converging,
            rate = max( rate/2,0.01 );   % make Lawson update conservative
            if doplot
                warning(['Lawson rate made conservative to ',num2str(rate)])
            end
              nrmincreased = 0;
            end

            if doplot  % plot error functions
                       % (hopefully near-equioscillating)
                subplot(2,1,1)
                plot(Z,F-Rori,'r.','markersize',8)
                title('AAA error')
                grid on, hold on
                if exist('hh','var')
                    set(hh,'color',(0.8)*[1 1 1]);
                end
                subplot(2,1,2)
                title('AAA-Lawson error')
                hh = plot(Z,F-R,'k.','markersize',8);
                grid on, hold on
                h2 = plot(z,0*z,'m.','markersize',12);
                ylim(err*[-1 1]); drawnow, shg
                if ( it == Lawsoniter ) % plot best function
                plot(Z,F-r(Z),'b.','markersize',10);
                end
            end
	      end

    % compute poles and roots
    B = eye(mn+1); B(1,1) = 0;
    E = [0 wf.'; ones(mn,1) diag(z)];
    zer = eig(E,B); zer = zer(~isinf(zer));   % zeros
    E = [0 w.'; ones(mn,1) diag(z)];
    pol = eig(E,B); pol = pol(~isinf(pol));   % poles
    dz = 1e-5*exp(2i*pi*(1:4)/4);
    res = r(bsxfun(@plus,pol,dz))*dz.'/4;     % residues

    if ( nargout > 8 )  % form p and q via N/D, NOTE: not always stable
    D = @(x) 0;    N = @(x) 0;     % form r=N/D in barycentric form
    for ii = 1:length(z)
       D = @(x) D(x) + w(ii)./(x-z(ii));
       N = @(x) N(x) + wf(ii)./(x-z(ii));
    end
    dom = [min(real(Z)) max(real(Z))];  % set domain
    x = chebpts(mmax+nmax+1,dom);
    node = @(x) prod(x-z);              % needed to check sign
    nodex = zeros(length(x),1);         % setup node values
    for ii = 1:length(x), nodex(ii) = node(x(ii)); end
    qvals = nodex.*feval(D,x);          % values of p,q
    pvals = nodex.*feval(N,x);
    p = chebfun(pvals,dom); q = chebfun(qvals,dom); % form chebfuns
    end
end

% parse Inputs for Lawson:

function [F, Z, m, n, Lawsoniter, tolLawson, doplot, tol ] = ...
    parseInputslawson(F, varargin)
% Input parsing for AAAmn_lawson.

% Check if F is empty:
if ( isempty(F) )
    error('CHEBFUN:aaamn_lawson:emptyF', 'No function given.')
elseif ( isa(F, 'chebfun') )
    if ( size(F, 2) ~= 1 )
        error('CHEBFUN:aaamn_lawson:nColF', ...
            'Input chebfun must have one column.')
    end
end

% Domain:
if ( isa(F, 'chebfun') )
    dom = F.domain([1, end]);
else
    dom = [-1, 1];
end

% Sample points:
if ( ~isempty(varargin) && isfloat(varargin{1}) )
    if length(varargin{1})>2   % sample points Z given.
    Z = varargin{1};
    varargin(1) = [];
    else                       % sample points not given.

    end
end

% m,n
if ( ~isempty(varargin) && isfloat(varargin{1}) )
    if length(varargin{1})>1 % input (f,Z,[m n])
    mn = varargin{1};
    m = mn(1); n = mn(2);
    varargin(1) = [];
    elseif isfloat(varargin{2}) % input (f,Z,m,n)
    m = varargin{1};
    n = varargin{2};
    varargin([1, 2]) = [];
    else
    m = varargin{1};
    n = m; % input (f,Z,m), default to diagonal type (m,m)
    varargin(1) = [];
    end
end

if ( ~exist('m', 'var') )
     warning(['CHEBFUN:aaamn_lawson: type (m,n) not specified,', ...
         ' default to (10,10)'])
     m = 10; n = 10;
end

% Set defaults for other parameters:
tolLawson = 1e-5;                       % Relative tolerance for
                                        % Lawson update.
tol = 1e-15;                            % AAA tolerance
Lawsoniter = max([5 min([20, m, n])]);  % Maximum number of terms.
doplot = 0;                             % Don't plot intermediate functions
                                        % unless specified

% Check if parameters have been provided:
while ( ~isempty(varargin) )
    if ( strncmpi(varargin{1}, 'tol', 3) ||  ...
            strncmpi(varargin{1}, 'tolLawson', 3) )
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 1]) )
            tolLawson = varargin{2};   % Lawson tolerance
        end
        varargin([1, 2]) = [];

    elseif ( strncmpi(varargin{1}, 'iter', 4) || ...
            strncmpi(varargin{1}, 'maxit', 5))
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 1]) )
            Lawsoniter = varargin{2};  % maximum Lawson iterations
        else
        warning(['CHEBFUN:aaamn_lawson:iter unspecified,', ...
            ' use default itermax ', num2str(Lawsoniter)])
        end
        varargin([1, 2]) = [];

    elseif ( strncmpi(varargin{1}, 'dom', 3) )
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 2]) )
            dom = varargin{2};
        end
        varargin([1, 2]) = [];
        if ( isa(F, 'chebfun') )
            if ( ~isequal(dom, F.domain([1, end])) )
                warning('CHEBFUN:aaamn_lawson:dom', ...
                    ['Given domain does not match the domain ', ...
                    'of the chebfun.\n', 'Results may be inaccurate.'])
            end
        end

    elseif strncmpi(varargin{1}, 'plot', 4)  % plot error functions
        if isfloat(varargin{2})
            doplot = varargin{2};
        elseif ( strncmpi(varargin{2}, 'true', 4) || ...
                strncmpi(varargin{2}, 'on', 2) )
            doplot = 1;
        end
        varargin([1, 2]) = [];
    else
        error('CHEBFUN:aaamn_lawson:UnknownArg', 'Argument unknown.')
    end
end

% Deal with Z and F:
if ( ~exist('Z', 'var') && isfloat(F) )
    % F is given as data values, pick same number of sample points:
    Z = linspace(dom(1), dom(2), length(F)).';
end

if ( exist('Z', 'var') )
    % Work with column vector:
    Z = Z(:);
    M = length(Z);

    % Function values:
    if ( isa(F, 'function_handle') || isa(F, 'chebfun') )
        % Sample F on Z:
        F = F(Z);
    elseif ( isnumeric(F) )
        % Work with column vector and check that it has correct length.
        F = F(:);
        if ( length(F) ~= M )
            error('CHEBFUN:aaamn_lawson:lengthFZ', ...
                'Inputs F and Z must have the same length.')
        end
    else
        error('CHEBFUN:aaamn_lawson:UnknownF', ...
            'Input for F not recognized.')
    end

else
    % Z was not given.  Set flag that Z needs to be determined.
    % Also set Z and M since they are needed as output.
    % in AAA this is done adaptively. This can be done with Lawson, but
    % probably safe to take as many points as reasonably possible here.
    Z = linspace(dom(1), dom(end), 4000).';
end

end % End of PARSEINPUT().

% generate function handle, interpolatory mode
function r = rrint(zz,z,w,f)                 % evaluate r at zz
zv = zz(:);                               % vectorize zz if necessary
CC = 1./bsxfun(@minus,zv,z.');            % Cauchy matrix
r = (CC*(w.*f))./(CC*w);                  % AAA approx as vector
ii = find(isnan(r));                      % Find values NaN = 0/0 if any
for j = 1:length(ii)
  r(ii(j)) = f(find(zv(ii(j))==z));       % Force interpolation there
end
r = reshape(r,size(zz));                  % AAA approx
end

% generate function handle, non-interpolatory mode
function r = rrab(zz,z,w,wf,f)                 % evaluate r at zz
zv = zz(:);                               % vectorize zz if necessary
CC = 1./bsxfun(@minus,zv,z.');            % Cauchy matrix
r = (CC*(wf))./(CC*w);                  % AAA approx as vector

ii = find(isnan(r));                      % Find values NaN = 0/0 if any
for j = 1:length(ii)
  r(ii(j)) = f(find(zv(ii(j))==z));       % Force interpolation there
end

r = reshape(r,size(zz));                  % AAA approx
end

% find null space for nondiagonal case m~=n
function Q = orthspace(z,dim,q)    % orthonormal projection space for (m,n)
if ( dim == 0 ), Q = eye(length(z)); end
if ( nargin < 3 ), q = ones(length(z),1); end
                Q = q/norm(q);
                for ii = 2:dim               % orthogonal complement via
                                             % Lanczos-type process
                Qtmp = diag(z)*Q(:,end);
                Qtmp = Qtmp - Q*(Q'*Qtmp);   % orthogonalize
                Qtmp = Qtmp - Q*(Q'*Qtmp);   % orthogonalize again (CGS2)
                Qtmp = Qtmp/norm(Qtmp);      % normalize
                Q = [Q Qtmp];
                end
            [Q,~] = qr(Q); Q = conj(Q(:,dim+1:end)); % desired null space
end

function r = rootsdiff(vals,dom,err_handle) % vals is either a function or
                                            % values at chebpts
% returns the roots of diff(vals) via ChebyshevU-colleague
if length(vals)<=1
if nargin<3, n = 2^5; end
xx = chebpts(n+1,dom);
vals = vals(xx);
else
n = size(vals,1)-1;
end

tol = 1e-3; % no need for high tolerance
c = 1;   % initialize
while ( (abs(c(end)/c(1))>tol) && (n<=2^6) )% sample until happy
cc = fft([vals(end:-1:1);vals(2:end-1)])/n;
cc(1) = cc(1)/2;
c = real(cc(1:n+1)); % coeffs of f=err_handle in T
cU = c(2:end).*(1:length(c)-1).'; % coeffs of df in U
% simplify; no need to get full accuracy.
% Then reorder to highest coeffs first.
len = max( find((abs(cU)/norm(cU)>1e-14)) ); cU = flipud(cU(1:len));
%if ( length(cU)<=1 || norm(cU./max(abs(cU)))<1e-14 ), r = []; return; end
% constant function
if ( length(cU)<=1 ), r = []; return; end
if abs(c(end)/c(1))>tol  % resample at finer grid
    n = 2*n;
    vals = feval(err_handle,(dom(1)+dom(end))/2 + ...
        cos(pi*((n:-1:0).')/n)*(dom(end)-dom(1))/2);
end
end

if length(cU)<=1, r = []; return; end % constant function

if (length(cU)==2)
    % degree 1 polynomial: just take the root
    ei = -cU(2)/(2*cU(1));
    % remove root if outside domain
    ei = ei(abs(ei)<=1+1e-7);
else
    % now construct colleague matrix for ChebyshevU
    oh = ones(len-2,1)/2;
    C = diag(oh,1) + diag(oh,-1);
    cU = -cU(2:end)/cU(1)/2;cU(2) = cU(2)+.5;
    C(1,:) = cU.';
    ei = eig(C);
    % remove irrelevant roots
    ei = real(ei(abs(imag(ei))<1e-5 & abs(ei)<=1+1e-7));
end
r = (dom(1)+dom(2))/2 + ei*(dom(2)-dom(1))/2; % map back to the subinterval
end

% Compute the poles zeros of the barycentric approximant with
% weights alpha and beta.
function [zer, pol] = pzeros(zj, alpha, beta, rh, m, n, dom)

if ( n == 0 || isempty(beta) )
    pol = [];
else
    l = length(beta);

    % Compute poles via generalized eigenvalue problem:
    B = eye(l+1);
    B(1,1) = 0;
    E = [0 beta.'; ones(l, 1) diag(zj)];
    pol = eig(E, B);
    % Remove zeros of denominator at infinity:
    pol = pol(~isinf(pol));

    rad = 1e-5; % radius for approximating residual

    if ( l - 1 > n )% superdiagonal case, remove irrelevant poles
    dz = rad*exp(2i*pi*(1:4)/4);
    res = rh(bsxfun(@plus, pol, dz))*dz.'/4; % residues
    ix = find( abs(res) > 1e-10 ); % pole with suff. residues
    pol = pol(ix);
    zerBern = abs(pol-dom(1)) + abs(pol-dom(2)); % Bernstein ellipse radius
    [~,ix] = sort( zerBern, 'ascend'); % sort wrt radius
    pol = pol( ix(1:min(n,end)) ); % choose <=n zeros with largest residues
    end
end

if ( m == 0 && isempty(alpha) )
    zer = [];
else
    l = length(alpha);

    % Compute zeros via generalized eigenvalue problem:
    B = eye(l+1);
    B(1,1) = 0;
    E = [0 alpha.'; ones(l, 1) diag(zj)];

    zer = eig(E,B);
    % Remove zeros of numerator at infinity:
    zer = zer(~isinf(zer));

    % subdiagonal case, remove irrelevant zeros:
    if ( l - 1 > m )
        rad = 1e-5; % radius for approximating derivative
        dz = rad*exp(2i*pi*(1:4)/4);
        deriv = sum(abs(bsxfun(@minus,rh(zer),rh(dz))./bsxfun(@minus,zer,dz)),2);
        deriv = deriv/4;
        ix = find( abs(deriv) > 1e-10 );

        zer = zer(ix);
        zerBern = abs(zer-dom(1)) + abs(zer-dom(2));
        [~,ix] = sort(zerBern,'ascend');
        zer = zer( ix(1:min(m,end)) );
    end

end

end % End of PZEROS().

function op = str2op(op)
    % Convert string inputs to either numeric format or function_handles.
    sop = str2num(op);
    if ( ~isempty(sop) )
        op = sop;
    else
        depVar = symvar(op);
        if ( numel(depVar) ~= 1 )
            error('CHEBFUN:CHEBFUN:str2op:indepvars', ...
             'Incorrect number of independent variables in string input.');
        end
        op = eval(['@(' depVar{:} ')', op]);
    end
end


text:
function fo=rdir(name,opt)
%RDIR Recursive list directory.
%   F = RDIR(NAME) returns a cell array of file names matching NAME.
%   NAME may be a cell array of strings, and may contain cascading
%   pathnames separated by '/' (or '\') followed by a filename. Wildcards
%   (*) may be used both in the final filename, as in DIR, but also in the
%   intermediate pathnames. For example, F=RDIR('mydir*/*.m') returns all
%   the M-Files contained in each directory that begins with 'mydir'.
%
%   File number indexing: Brackets [] may be used in the file or directory
%   names for convenient number indexing. If the input string is in the form
%   'PP[RANGE]SS', all files or directories matching 'PP0000nSS' will be
%   returned, where 'PP' and 'SS' are prefix and suffix substrings, and n
%   is an index lying in the range RANGE (paded with 5 zeros by default).
%   RANGE is a vector, that can be in the form [N1 N2 N3..], or START:END,
%   or START:STEP:END, or any other MATLAB valid syntax.
%
%   Use 'PP[RANGE,NZ]SS') to specify the number of zeros to pad the index
%   string (NZ=5 by default). For example, 'B[1:4,2].v*' gives
%   {'B01.v*','B02.v*','B03.v*','B04.v*'}.
%
%   Use 'PP[RANGE,NZ.NP]SS' to specify the number of digits NP after the
%   decimal point (NP=0 by default). For example, 'dt=[1:0.5:2,2.3]s' gives
%   {'dt=1.000s','dt=1.500s','dt=2.000s'}.
%
%   If the input string has more than one bracket pair [], RDIR is
%   called recursively for each pair. For example, 'B[1:4,2]_[1 2,1]'
%   gives {'B01_1','B01_2','B02_1','B02_2','B03_1','B03_2','B04_1','B04_2'}
%
%   RDIR file_name  or  RDIR('file_name') displays the result.
%
%   F = RDIR(NAME,'fileonly') only returns file names.
%   F = RDIR(NAME,'dironly') only returns directory names, but does not
%   list their content.
%   F = RDIR(NAME,'filedir') returns both files and directories (as DIR),
%   by default.
%   (Note that, when options 'dironly' or 'filedir' are used, the
%   fictitious directories '.' (current) and '..' (parent) are not
%   returned, contrarily to DIR).
%
%   Examples:
%
%     F = RDIR('set*/B[1:8,2].v*')  returns the file names matching
%     'B01.v*' to 'B08.*' in each directory matching 'set*'.
%
%     F = RDIR('set[1 2 3]*/*/B01.vec')  returns the file name 'B01.vec',
%     if present, in all the subdirectories of each directory matching
%     'set1*', 'set2*', 'set3*'.
%
%     F = RDIR('Myfiles_dt=[1:0.5:2,5.2]s.*') returns the file names
%     'Myfiles_dt=01.00s.*', 'Myfiles_dt=01.50s.*', 'Myfiles_dt=02.00s.*'
%
%   See also DIR, LSW, CDW.

%   F. Moisy, moisy_at_fast.u-psud.fr
%   Revision: 1.20,  Date: 2010/05/12

% History:
% 2004/10/04: v1.00, first version.
% 2006/11/29: v1.01, 'filedir' is now the new default parameter
% 2010/01/19: v1.10, use filesep ('' or '/')
% 2010/05/01: v1.11, replace all / and \ with filesep, and add a *.
% 2010/05/05: v1.12, bug fixed with absolute path names (eg, C:\*)
% 2010/05/07: v1.13, removes the prefix './' if any
% 2010/05/12: v1.12, accepts non-integer filenumber indexing (based on
% expandstr v1.10; actually, only the help text has been changed in rdir)

error(nargchk(1,2,nargin));

verb=0; % set verb=1 for debugging (verbose), or verb=0 for normal use.

if nargin==1, opt='filedir'; end    % changed V1.01

if iscell(name),
    f={};
    for i=1:length(name),
        if verb, disp(['im in ' pwd ' and im calling rdir for ' name{i}]); end
        ff=rdir(name{i},opt);
        f={f{:} ff{:}};
    end
else
    % new v1.11:
    name = strrep(name,'\',filesep);
    name = strrep(name,'/',filesep);

    % new v1.13: removes any './'
    name = strrep(name, ['.' filesep], '');

    p=findstr(name,filesep);

    % new v1.12: 'C:\' not considered has having a '\'
    if findstr(name,[':' filesep]) & length(p)==1
        p=[];
    end

    if isempty(p) % if the input string contains no filesep
        if verb, disp(['im in ' pwd ', with ' name ' that contains no ' filesep]); end

        name=expandstr(name); % tries to expand the file name
        if ~iscell(name), % if there were nothing to expand
            % then simply returns the filenames matching NAME:
            if verb, disp(['im in ' pwd ' and im going to call dir for ' name]); end
            f=dir(name);
            switch lower(opt)
                case 'fileonly', f=f(~[f.isdir]);
                case 'dironly', f=f([f.isdir]);
                case {'filedir','dirfile'}, %f=f;
                otherwise, error(['Unknown option ''' opt '''']);
            end
            f={f.name};
            ff={};
            for i=1:length(f),
                if ((~strcmp(f{i},'..')) && (~strcmp(f{i},'.'))), % excludes '.' and '..'
                    ff={ff{:} f{i}};
                end
            end
            f=ff;
        else
            % something has been expanded, so calls rdir for each file name
            f={};
            for i=1:length(name),
                if verb, disp(['im in ' pwd ' and im going to call rdir for ' name{i}]); end
                ff=rdir(name{i},opt);
                for j=1:length(ff),
                    %f={f{:} [name{i} '/' ff{j}]};
                    f={f{:} ff{j}};
                end
            end
        end
    else
        if verb, disp(['im in ' pwd ', with ' name ' that contains ' filesep]); end
        pname=name(1:(p-1)); % first pathname of the input string
        fname=name((p+1):end); % next pathnames and final filename
        pname=expandstr(pname); % tries to expand the first path name
        if verb, disp(['pname=' pname ',  fname=' fname]); end
        if ~iscell(pname), % if there were nothing to expand,
            if isempty(findstr(pname,'*'))
                if verb, disp([pname ' contains no *']); end
                if exist(pname,'dir')
                    d={pname};
                else
                    d={};
                end
            else
                if verb, disp([pname ' contains *(s)']); end
                d=dir(pname); % looks for pathnames matching the first pathname
                d=d([d.isdir]); % keeps only directories
                d={d.name}; % extracts the names
            end
            f={};
            for i=1:length(d),
                if ((~strcmp(d{i},'..')) && (~strcmp(d{i},'.'))), % excludes '.' and '..'
                    if verb, disp(['im in ' pwd ', and im entering in ' d{i}]); end

                    trep=d{i};  % new v1.12
                    if isequal(trep(end),':')
                        trep = [trep filesep];
                    end
                    cd(trep)

                    % calls (possibly recursively) rdir for the next pathnames:
                    if verb, disp(['im in ' pwd ' and im calling rdir for ' fname]); end
                    ff=rdir(fname,opt);
                    % builds the output cell array:
                    for j=1:length(ff),
                        f={f{:} [d{i} filesep ff{j}]};
                    end
                    cd ..
                end
            end
        else
            % if pname has been expanded, then calls rdir for each pathname
            f={};
            for i=1:length(pname),
                ff=rdir([pname{i} filesep fname],opt);
                f={f{:} ff{:}};
            end
        end
    end
end

if nargout==0
    disp(f');
else
    fo=f;
end

author: Sable

repo: mcbench-benchmarks


text:
function e = end(X,k,n)
%END Last index of indexing expression for ktensor.
%
%   The expression X(end,:,:) will call END(X,1,3) to determine
%   the value of the first index.
%
%   See also KTENSOR, KTENSOR/SUBSREF, END.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%TODO (after 2.0 release): Resolve ambiguity w.r.t X{end}and X(end,1,1)
%for 1st-order tensors.

if n > ndims(X)
  error('Subscript out of range.');
end

if (n ~= 1)
  e = size(X,k);
else
  e = ndims(X);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ THD, ph, amp ] = compute_THD( t,x, freq )
% function [ THD, ph, amp ] = compute_THD( t,x, freq )
%
% Written by Dr. Yoash Levron
% February 2013.
%
% computes the Total-Harmonic-Distortion (THD)
% of a signal x(t). The amplitude and phase of the
% basic harmonic are also computed. These values
% are typically useful in power systems, audio signal
% processing, and other related fields.
%
% DC offset does not affect THD.
%
% The function computes the basic harmonic
% of the signal, in the form:
% x(t) = amp*cos(w*t - ph) + (higher Harmonics)
% where :  w = 2*pi*freq
% so 'amp' and 'ph' are the phase and amplitude
% of the basic harmonic.
%
% inputs:
% t - [sec] time vector. (should be periodical with basic harmonic 'freq')
% x - signal vector.
% freq - [Hz] frequency of the basic harmonic.
%
% outputs:
% THD - total harmonic distortion (the scale is 1 = 100%).
% ph - [rad] phase of the basic harmonic.
% amp - Amplitude of the basic harmonic.

%%%%%%%%%%%% start function %%%%%%%%%%

%%% check that t,x are the same length
if (length(t) ~= length(x))
      'Error: t and x should be the same length'
      THD = NaN;   ph = NaN; amp = NaN;
      beep; return
end

if (size(t,2) == 1)
    t = t.';
end
if (size(x,2) == 1)
    x = x.';
end

%%% condition input time vector
input_error = 0;
%%% add two samples to complete the last cycle.
t = t - t(1);   % remove any time shift
dtt = t(end) - t(end-1);
t = [t (t(end)+dtt) (t(end)+2*dtt)];
x = [x x(end) x(end)];
T = t(end);
if (T < (1/freq) )
    input_error = 1;
end
if (input_error)
        'Error: Input time vector is illegal. Time samples should expand at least to 1/freq'
        THD = NaN;   ph = NaN; amp = NaN;
        beep; return
end

% truncate extra samples, to fit in an integer number of cycles of freq
T = floor(T * freq) /freq;

% resample on a linear grid:
% t1, x1 is the new input, not including the last sample
x = x - sum(x)/length(x);  % remove any DC offset
N = max(1e6, length(x));  % number of samples
dt = T/N;
t1 = 0:dt:(T-dt);
x1 = interp1(t,x,t1,'cubic');

%%% compute cos-sin fourier coefficients
w = 2*pi*freq;
acs = (2/T) * sum(x1.*cos(w*t1))*dt;  % basic frequency cos coefficient.
bsn = (2/T) * sum(x1.*sin(w*t1))*dt;  % basic frequency sin coefficient.
amp = (acs^2 + bsn^2)^0.5;
ph = pi/2 - sign(acs) * acos(bsn/amp);

rms22 = (2/T) * sum(x1.^2) * dt;
THD = (rms22/amp^2 - 1)^0.5;

% correct phase to be in the range [-pi : pi]
if (ph > pi)
    ph = ph - 2*pi;
end
if (ph < -pi)
    ph = ph + 2*pi;
end

end


text:
function [ meta, symmetry ] = meta2sicd_tiffrs( tiff_filename, gen )
%META2SICD_TIFFRS Compile SICD metadata for RS2 or RCM data package
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Initialize outputs
meta=struct();
symmetry=[0 0 0];

% RS2/RCM TIFF doesn't have much metadata in the TIFF itself, but it should
% have a co-located product.xml file
[pathname, fileonly, ext] = fileparts(tiff_filename);
jFile=java.io.File(tiff_filename);
tiff_filename = char(jFile.getCanonicalPath);
if strcmp(gen, 'RS2')
    productxmlfile = fullfile(pathname,'product.xml'); % Main metadata file
elseif strcmp(gen, 'RCM')
    productxmlfile = fullfile(fileparts(pathname),'metadata','product.xml'); % Main metadata file
end
files_found = dir(productxmlfile);
if length(files_found)==1
    xp=javax.xml.xpath.XPathFactory.newInstance.newXPath();
    rsxml_meta=xmlread(productxmlfile);
    meta=meta2sicd_rs_xml(rsxml_meta);

    % XML describes all images of a polarimetric collection
    % Select metadata associated with only this TIFF file
    % Must be a more elegant way to do this in a single call to XPath...
    for i=1:length(meta)
        if strcmp(gen, 'RS2')
            current_file_pol=char(xp.evaluate(['/*[local-name()=''product'']'...
                                       '/*[local-name()=''imageAttributes'']'...
                                       '/*[local-name()=''fullResolutionImageData'']'...
                                       '[' num2str(i) ']'],...
                                      rsxml_meta)); % Check name of one image
        elseif strcmp(gen, 'RCM')
            current_file_pol=char(xp.evaluate(['/*[local-name()=''product'']'...
                                       '/*[local-name()=''sceneAttributes'']'...
                                       '/*[local-name()=''imageAttributes'']'...
                                       '[' num2str(i) ']'...
                                       '/*[local-name()=''ipdf'']'],...
                                      rsxml_meta)); % Check name of one image
        end
        jFile=java.io.File(fullfile(fileparts(productxmlfile), current_file_pol));
        current_file_pol = char(jFile.getCanonicalPath);  % RCM IDPF are given as relative paths
        if strcmp(tiff_filename,current_file_pol) % Does it match this TIFF?
            meta=meta{i}; % If so, use the pole value from this image
            meta.native.product_xml=rsxml_meta;
            if strcmp(gen, 'RS2')
                pol=char(xp.evaluate(['/*[local-name()=''product'']'...
                                      '/*[local-name()=''imageAttributes'']'...
                                      '/*[local-name()=''fullResolutionImageData'']'...
                                      '[' num2str(i) ']/@pole'],...
                                      rsxml_meta));
                pol = {pol(1) pol(2)};
            elseif strcmp(gen, 'RCM')
                pol=char(xp.evaluate(['/*[local-name()=''product'']'...
                                      '/*[local-name()=''sceneAttributes'']'...
                                      '/*[local-name()=''imageAttributes'']'...
                                      '[' num2str(i) ']'...
                                      '/*[local-name()=''ipdf'']/@pole'],...
                                      rsxml_meta));
                M = struct('H','H','V','V','C','RHC');
                pol = {M.(pol(1)) M.(pol(2))};
            end
            meta.ImageFormation.RcvChanProc.ChanIndex = i;
            meta.ImageFormation.TxRcvPolarizationProc=[pol{1} ':' pol{2}];
            break;
        end
    end
    if ~iscell(meta) % File was one of the components pointed to by this XML
        % Compute symmetry
        lookdir=upper(char(xp.evaluate(['/*[local-name()=''product'']'...
                               '/*[local-name()=''sourceAttributes'']'...
                               '/*[local-name()=''radarParameters'']'...
                               '/*[local-name()=''antennaPointing'']'],...
                              rsxml_meta)));
        if strcmp(gen, 'RS2')
            ia_str = 'imageAttributes';
        elseif strcmp(gen, 'RCM')
            ia_str = 'imageReferenceAttributes';
        end
        lineOrder=char(xp.evaluate(['/*[local-name()=''product'']'...
                                    '/*[local-name()=''' ia_str ''']'...
                                    '/*[local-name()=''rasterAttributes'']'...
                                    '/*[local-name()=''lineTimeOrdering'']'],...
                                    rsxml_meta));
        sampleOrder=char(xp.evaluate(['/*[local-name()=''product'']'...
                                 '/*[local-name()=''' ia_str ''']'...
                                 '/*[local-name()=''rasterAttributes'']'...
                                 '/*[local-name()=''pixelTimeOrdering'']'],...
                                rsxml_meta));
        symmetry(1)=xor(strcmp(lineOrder,'Decreasing'),lookdir(1)=='L');
        symmetry(2)=strcmp(sampleOrder,'Decreasing');
    end
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function test_bug3297

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/bug3297/testMEG001_1200hz_20170517_05.ds');

event = ft_read_event(filename);

% it only returned 3 'trial' events upon reading with the verion up to 22 May 2017
% after inlcuding type 20 as trigger, it returned 11 events
assert(numel(event)>3, 'incorrect number of events read from file');

% I looked it up in "CTF MEG(TM) File Formats" pdf, Release 5.2.1
%
% eMEGReference     0 Reference magnetometer channel
% eMEGReference1    1 Reference 1st-order gradiometer channel
% eMEGReference2    2 Reference 2nd-order gradiometer channel
% eMEGReference3    3 Reference 3rd-order gradiometer channel
% eMEGSensor        4 Sensor magnetometer channel located in head shell
% eMEGSensor1       5 Sensor 1st-order gradiometer channel located in head shell
% eMEGSensor2       6 Sensor 2nd-order gradiometer channel located in head shell
% eMEGSensor3       7 Sensor 3rd-order gradiometer channel located in head shell
% eEEGRef           8 EEG unipolar sensors not on the scalp
% eEEGSensor        9 EEG unipolar sensors on the scalp
% eADCRef           10 (see eADCAmpRef below)
% eADCAmpRef        10 ADC amp channels from HLU or PIU (old electronics)
% eStimRef          11 Stimulus channel for MEG41
% eTimeRef          12 Time reference coming from video channel
% ePositionRef      13 Not used
% eDACRef           14 DAC channel from ECC or HLU
% eSAMSensor        15 SAM channel derived through data analysis
% eVirtualSensor    16 Virtual channel derived by combining two or more physical channels
% eSystemTimeRef    17 System time showing elapsed time since trial started
% eADCVoltRef       18 ADC volt channels from ECC
% eStimAnalog       19 Analog trigger channels
% eStimDigital      20 Digital trigger channels
% eEEGBipolar       21 EEG bipolar sensor not on the scalp
% eEEGAflg          22 EEG ADC over range flags
% eMEGReset         23 MEG resets (counts sensor jumps for crosstalk pur- poses)
% eDipSrc           24 Dipole source
% eSAMSensorNorm    25 Normalized SAM channel derived through data analy- sis
% eAngleRef         26 Orientation of head localization field
% eExtractionRef    27 Extracted signal from each sensor of field generated by each localization coil
% eFitErr           28 Fit error from each head localization coil
% eOtherRef         29 Any other type of sensor not mentioned but still valid
% eInvalidType      30 An invalid sensor

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function edges = meshEdges(faces, varargin)
%MESHEDGES Computes array of edge vertex indices from face array.
%
%   EDGES = meshEdges(FACES);
%
%   Example
%     meshEdges
%
%   See also
%     meshes3d, meshEdgeFaces, meshFaceEdges

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-06-28,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

%   HISTORY
%   2013-08-22 rename from computeMeshEdges to meshEdges, add more control
%       on inputs

%% Process input arguments

if isstruct(faces) && isfield(faces, 'faces')
    % if input is a mesh structure, extract the 'faces' field
    faces = faces.faces;
elseif nargin > 1
    % if two arguments are given, keep the second one
    faces = varargin{1};
end

if ~iscell(faces)
    %% Process faces given as numeric array
    % all faces have same number of vertices, stored in nVF variable

    % compute total number of edges
    nFaces  = size(faces, 1);
    nVF     = size(faces, 2);
    nEdges  = nFaces * nVF;

    % create all edges (with double ones)
    edges = zeros(nEdges, 2);
    for i = 1:nFaces
        f = faces(i, :);
        edges(((i-1)*nVF+1):i*nVF, :) = [f' f([2:end 1])'];
    end

else
    %% faces are given as a cell array
    % faces may have different number of vertices

    % number of faces
    nFaces  = length(faces);

    % compute the number of edges
    nEdges = 0;
    for i = nFaces
        nEdges = nEdges + length(faces{i});
    end

    % allocate memory
    edges = zeros(nEdges, 2);
    ind = 0;

    % fillup edge array
    for i = 1:nFaces
        % get vertex indices, ensuring horizontal array
        f = faces{i}(:)';
        nVF = length(f);
        edges(ind+1:ind+nVF, :) = [f' f([2:end 1])'];
        ind = ind + nVF;
    end

end

% keep only unique edges, and return sorted result
edges = sortrows(unique(sort(edges, 2), 'rows'));

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function title = p00_title ( problem )

%*****************************************************************************80
%
%% P00_TITLE returns a title for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the number of the problem.
%
%    Output, string TITLE, a title for the problem.
%
  if ( problem == 1 )
    title = p01_title ( );
  elseif ( problem == 2 )
    title = p02_title ( );
  elseif ( problem == 3 )
    title = p03_title ( );
  elseif ( problem == 4 )
    title = p04_title ( );
  elseif ( problem == 5 )
    title = p05_title ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_TITLE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of PROBLEM = %d', problem );
    error ( 'P00_TITLE - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
 function out = subsref(ob, arg)
%function out = subsref(ob, arg)
% evaluate reference of the form ob.ref or ob(ref,:)

%
% ob.?
% get from Gblock object if possible, otherwise from base object
%

if arg(1).type == '.'

	if isfield(struct(ob), arg(1).subs)
		out = struct(ob);
		out = out.(arg(1).subs);
		if length(arg) > 1
			out = subsref(out, arg(2:end));
		end
		return
	end

	if isfield(struct(ob.base), arg(1).subs)
		out = struct(ob.base);
		out = out.(arg(1).subs);
		if length(arg) > 1
			out = subsref(out, arg(2:end));
		end
		return
	end

	% otherwise, error
	disp 'field names of base system'
	disp(fieldnames(ob.base))
	printf('subscript "%s" not a field of object', arg.subs)
	error bug

elseif length(arg) > 1
	disp '-------------'
	disp(arg)
%	disp(arg(:).type)
%	disp(arg(:).subs)
	error 'not done'

%
% ob{i_block}
% this is primary method 'handled' by this object!
%
elseif arg.type == '{}' && length(arg.subs) == 1

	out = ob;
	out.i_block = arg.subs{1};
	if out.i_block < 1 || out.i_block > out.nblock
		error 'bad block index'
	end

%
% ob(i_block)
% included for backward compatibility
%
elseif arg.type == '()' && length(arg.subs) == 1

	warning 'please use {} instead of () now for block objects'
	out = ob;
	out.i_block = arg.subs{1};
	if out.i_block < 1 || out.i_block > out.nblock
		error 'bad block index'
	end

%
% anything else gets passed to base object
% and we hope that it can handle it!
%
else
	tmp = subsref(ob.base, arg);
	if isa(tmp, class(ob.base))
		out = ob;
		out.base = tmp;
	else
		out = tmp;
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [ V ] = T2V( T )
%T2V converts 4x4 transformation matrix into a 1x6 vector
%   Inputs -
%   T - a standard 4x4 transformation matrix
%
%   Outputs -
%   V - 1x6 vector of form [x,y,z,rx,ry,rz] where x,y,z is the translation
%   and rx,ry,rz is an angle-axis representation of the angle where the
%   unit vector representing the axis has been multipled by the angle of
%   rotation about it

validateattributes(T, {'numeric'},{'size',[4,4]});

T = double(T);

V(1:3) = T(1:3,4);
V(4:6) = R2V(T(1:3,1:3));

end

author: ZacharyTaylor

repo: Camera-to-Arm-Calibration


text:
function r8sp_print ( m, n, nz_num, row, col, a, title )

%*****************************************************************************80
%
%% R8SP_PRINT prints a R8SP matrix.
%
%  Discussion:
%
%    This version of R8SP_PRINT has been specifically modified to allow,
%    and correctly handle, the case in which a single matrix location
%    A(I,J) is referenced more than once by the sparse matrix structure.
%    In such cases, the routine prints out the sum of all the values.
%
%    The R8SP storage format stores the row, column and value of each nonzero
%    entry of a sparse matrix.
%
%    It is possible that a pair of indices (I,J) may occur more than
%    once.  Presumably, in this case, the intent is that the actual value
%    of A(I,J) is the sum of all such entries.  This is not a good thing
%    to do, but I seem to have come across this in MATLAB.
%
%    The R8SP format is used by CSPARSE ("sparse triplet"), DLAP/SLAP
%    ("nonsymmetric SLAP triad"), by MATLAB, and by SPARSEKIT ("COO" format).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, integer NZ_NUM, the number of nonzero elements in the matrix.
%
%    Input, integer ROW(NZ_NUM), COL(NZ_NUM), the row and column indices
%    of the nonzero elements.
%
%    Input, real A(NZ_NUM), the nonzero elements of the matrix.
%
%    Input, string TITLE, a title.
%
  r8sp_print_some ( m, n, nz_num, row, col, a, 1, 1, m, n, title );

  return
end


text:
function [WRA, WF, WR] = compute_SWTFiltResp(params)
%% Load parameters
rs = params.rs;
cs = params.cs;
nlev = params.Wavelet.nlev;

lod = params.Wavelet.lod; %% Normalize filter for easier SWT implementation
hid = params.Wavelet.hid; %% Normalize filter for easier SWT implementation

% Get the shift-invariant wavelet filters
indicator = zeros(rs, cs);
indicator(1,1) = 1;
SWC = myswt2(indicator, nlev, lod, hid);

WF = zeros(rs, cs, 3*nlev+1);
WR = WF;
WRA = WR;

for ib = 1:3
    for ilev = 1:nlev
        WF(:,:,(ib-1)*nlev+ilev) = SWC(:,:,(ib-1)*nlev+ilev);
        WR(:,:,(ib-1)*nlev+ilev) = fft2(WF(:,:,(ib-1)*nlev+ilev));
        WRA(:,:,(ib-1)*nlev+ilev) = abs(WR(:,:,(ib-1)*nlev+ilev)).^2;
    end
end
WF(:,:,3*nlev+1) = SWC(:,:,3*nlev+1);
WR(:,:,3*nlev+1) = fft2(WF(:,:,3*nlev+1));
WRA(:,:,3*nlev+1) = abs(WR(:,:,3*nlev+1)).^2;

% DCR = 0;
% for ib = 1:3*nlev
%     DCR = DCR + WRA(:,:,ib);
% end
% dif = DCR-ACR;
% max(abs(dif(:)))