function Vq = interpn(varargin)
%MATLAB Code Generation Library Function

%   Limitations and Notes:
%   1. For best performance supply the input grid, X1,X2,..., as vectors.
%   2. 'cubic' and 'spline' methods are not supported for 2-D and higher
%      interpolation.
%   3. If supplied, the method input must be a constant.
%   4. The *METHOD option is not supported.

%   Copyright 1984-2017 The MathWorks, Inc.
%#codegen

narginchk(1,inf);
coder.internal.prefer_const(varargin);
[NARG,METHOD,EXTRAPp] = coder.const(@methodAndExtrapP,varargin{:});
coder.internal.assert(NARG >= 1,'MATLAB:minrhs');
[Vidx,ND,AUTOGRID] = coder.const(@findVArg,varargin{1:NARG});
V = varargin{Vidx};
if ND == 1
    % Dispatch to INTERP1. This is required since the algorithms below
    % assume min(size(V)) >= 2.
    if NARG < 3
        % AUTOGRID is not supported in MATLAB via INTERPN for 1-D
        % interpolation, so this is necessarily to be interpreted as
        % N-times syntax. Since INTERP1 doesn't have an N-times syntax, we
        % will have to implement it here.
        if NARG == 1
            Xq = colon(1,0.5,length(V));
        else
            coder.internal.assert(isscalar(varargin{2}) && ...
                (isnumeric(varargin{2}) || islogical(varargin{2})), ...
                'Coder:toolbox:NMustBeANumericScalar');
            delta = pow2(-floor(double(varargin{2}(1))));
            Xq = colon(1,delta,length(V));
        end
    else
        Xq = varargin{NARG};
    end
    % Call INTERP1, but translate 'cubic' to 'v5cubic'.
    if nargin > NARG && METHOD == CUBIC
        Vq = interp1(varargin{1:Vidx},Xq,'v5cubic',varargin{NARG+2:end});
    else
        Vq = interp1(varargin{1:Vidx},Xq,varargin{NARG+1:end});
    end
elseif NARG <= 2
    % N-times syntax.
    if NARG == 2
        coder.internal.assert(isscalar(varargin{2}) && ...
            (isnumeric(varargin{2}) || islogical(varargin{2})), ...
            'Coder:toolbox:NMustBeANumericScalar');
        N = floor(double(varargin{2}(1)));
    else
        N = 1;
    end
    Vq = interpNtimes(ND,V,N,METHOD);
else
    vtype = coder.internal.scalarEg(V);
    if nargin == NARG + 2
        coder.internal.assert(METHOD == NEAREST || ...
            coder.internal.isBuiltInNumeric(varargin{end}), ...
            'MATLAB:interp2:extrapvalNotNumeric');
        coder.internal.assert(isreal(varargin{end}) || ~isreal(V), ...
            'Coder:toolbox:interpn_extrapvalNotReal');
        extrapval = cast(varargin{end},'like',vtype);
    else
        extrapval = cast(coder.internal.nan,'like',vtype);
    end
    if AUTOGRID
        Vq = interpnDispatch(ND,METHOD,EXTRAPp,extrapval, ...
            varargin{Vidx:(ND + 1)});
    else
        % Validate grid types.
        coder.internal.assert(NARG == 2*ND + 1, ...
            'MATLAB:interpn:nargin');
        for d = coder.unroll(1:ND)
            coder.internal.assert(isreal(varargin{d}), ...
                'Coder:toolbox:GridPointsMustBeReal');
            coder.internal.assert(isfloat(varargin{d}), ...
                'MATLAB:mathcgeo_catalog:NonFloatInpPtsErrId');
            coder.internal.assert(isa(varargin{d},class(varargin{1})),...
                'MATLAB:mathcgeo_catalog:GridOfMixedDataTypesErrId');
        end
        if coder.const(allVectors(varargin{1:ND}))
            % Validate grid vector lengths.
            for d = coder.unroll(1:ND)
                coder.internal.assert( ...
                    numel(varargin{d}) == size(V,d), ...
                    'MATLAB:mathcgeo_catalog:MixedGridCoordSizeErrId');
                coder.internal.assert(isincreasing(varargin{d}), ...
                    'Coder:toolbox:NonIncreasing');
            end
            % Interpolate.
            Vq = interpnDispatch(ND,METHOD,EXTRAPp,extrapval, ...
                varargin{Vidx:(2*ND + 1)},varargin{1:ND});
        else
            % Validate grid array sizes.
            for d = coder.unroll(1:ND)
                coder.internal.assert( ...
                    size(varargin{d},d) == size(V,d),...
                    'MATLAB:mathcgeo_catalog:MixedGridCoordSizeErrId');
            end
            % Verify that the grid arrays are plaid, as from NDGRID.
            coder.internal.assert(isplaid('ndgrid',varargin{1:ND}), ...
                'MATLAB:mathcgeo_catalog:NotNdGridErrId');
            gridv = cell(1,ND);
            for d = coder.unroll(1:ND)
                gridv{d} = unmeshgrid(varargin{d},d);
            end
            % Interpolate.
            Vq = interpnDispatch(ND,METHOD,EXTRAPp,extrapval, ...
                varargin{Vidx:(2*ND + 1)},gridv{:});
        end
    end
end

%--------------------------------------------------------------------------

function Vq = interpNtimes(ND,V,N,METHOD)
% N-times interpolation. N must be an integer-valued, scalar double.
coder.internal.prefer_const(ND,N,METHOD);
delta = pow2(-N);
gridv = cell(1,ND);
for d = coder.unroll(1:ND)
    gridv{d} = colon(1,delta,size(V,d));
end
EXTRAPp = false;
extrapval = cast(0,'like',V);
grid = cell(1,ND);
[grid{:}] = ndgrid(gridv{:});
Vq = interpnDispatch(ND,METHOD,EXTRAPp,extrapval,V,grid{:});

%--------------------------------------------------------------------------

function Vq = interpnDispatch(ND,METHOD,EXTRAPp,extrapval,V,varargin)
% This layer will become useful when we add support for CUBIC.
coder.inline('always');
coder.internal.prefer_const(ND,METHOD,EXTRAPp);
coder.internal.assert(METHOD == NEAREST || METHOD == LINEAR, ...
    'Coder:toolbox:InterpnUnsupportedMethod');
coder.internal.assert(ndims(V) == ND,'MATLAB:interpn:nargin');
coder.internal.assert(METHOD == NEAREST || isfloat(V), ...
    'MATLAB:mathcgeo_catalog:NonFloatValuesErrId');
for d = coder.unroll(1:ND)
    coder.internal.assert(size(V,d) >= 2, ...
        'Coder:toolbox:NotEnoughPoints');
    coder.internal.assert( ...
        isequal(size(varargin{d}),size(varargin{1})), ...
        'MATLAB:mathcgeo_catalog:InputMixSizeErrId');
    coder.internal.assert(isreal(varargin{d}), ...
        'Coder:toolbox:QueryPointsMustBeReal');
    coder.internal.assert(isfloat(varargin{d}) && ...
        isa(varargin{d},class(varargin{1})), ...
        'MATLAB:mathcgeo_catalog:QueryOfMixedDataTypesErrId');
end
Vq = interpnLocal(ND,METHOD,EXTRAPp,extrapval,V,varargin{:});

%--------------------------------------------------------------------------

function Vq = interpnLocal(ND,METHOD,EXTRAPp,extrapval,V,varargin)
% The local interpolation methods, linear and nearest.
% This function assumes that all grid inputs have lengths at least 2 and
% the sizes are consistent with the size of V. The inputs varargin{1:ND}
% are the values at which to interpolate. Inputs varargin{ND+1:2*ND}, if
% supplied, are the grid vectors.
coder.internal.prefer_const(ND,METHOD,EXTRAPp);
AUTOGRID = length(varargin) == ND;
if METHOD == LINEAR
    interpf = @scalarMultilinearInterp;
else
    interpf = @scalarNearestInterp;
end
Vq = coder.nullcopy(zeros(size(varargin{1}),'like',extrapval));
n = coder.internal.indexInt(numel(Vq));
if coder.internal.isConst(size(V))
    [vbbidx,stride] = coder.const(@BoundingBoxIndexArrays, ...
        false(size(V)));
else
    [vbbidx,stride] = BoundingBoxIndexArrays(V);
end
if AUTOGRID
    xmin = ones(1,ND,coder.internal.indexIntClass);
    xmax = coder.nullcopy(zeros(1,ND,coder.internal.indexIntClass));
    for d = 1:ND
        xmax(d) = coder.internal.indexInt(size(V,d));
    end
else
    gtype = coder.internal.scalarEg(varargin{ND+1:2*ND});
    xmin = coder.nullcopy(zeros(1,ND,'like',gtype));
    xmax = coder.nullcopy(zeros(1,ND,'like',gtype));
    for d = coder.unroll(1:ND)
        xmin(d) = varargin{ND + d}(1);
        xmax(d) = varargin{ND + d}(end);
    end
end
if coder.internal.useParforConst('interpn_orders_0_1',n)
    if AUTOGRID
        parfor k = 1:n
            Xq = ndslice(k,varargin{1:ND}); %#ok<PFBNS>
            Vq(k) = localLoopBody(ND,interpf,EXTRAPp,extrapval, ...
                vbbidx,stride,xmin,xmax,V,Xq);
        end
    else
        parfor k = 1:n
            Xq = ndslice(k,varargin{1:ND}); %#ok<PFBNS>
            Vq(k) = localLoopBody(ND,interpf,EXTRAPp,extrapval, ...
                vbbidx,stride,xmin,xmax,V,Xq,varargin{ND+1:2*ND});
        end
    end
else
    if AUTOGRID
        for k = 1:n
            Xq = ndslice(k,varargin{1:ND});
            Vq(k) = localLoopBody(ND,interpf,EXTRAPp,extrapval, ...
                vbbidx,stride,xmin,xmax,V,Xq);
        end
    else
        for k = 1:n
            Xq = ndslice(k,varargin{1:ND});
            Vq(k) = localLoopBody(ND,interpf,EXTRAPp,extrapval, ...
                vbbidx,stride,xmin,xmax,V,Xq,varargin{ND+1:2*ND});
        end
    end
end

%--------------------------------------------------------------------------

function Vqk = localLoopBody(ND,interpf,EXTRAPp,extrapval, ...
    vbbidx,stride,xmin,xmax,V,Xq,varargin)
% For AUTOGRID, varargin is not supplied, otherwise varargin is N elements
% long: X1,X2,...,XN. XK is a vector of grid points in the K-th dimension.
coder.inline('always');
coder.internal.prefer_const(ND,interpf,EXTRAPp,vbbidx,stride);
AUTOGRID = coder.const(isempty(varargin));
if ~EXTRAPp && isOutOfRange(ND,Xq,xmin,xmax)
    Vqk = extrapval;
else
    voffset = coder.internal.indexInt(0);
    % Construct the bounding box of the interpolation point.
    if AUTOGRID
        xbox = coder.nullcopy(zeros(2,ND));
        for d = 1:ND
            k = autogridLookup(Xq(d),xmax(d) - 1);
            voffset = voffset + (k - 1)*stride(d);
            xbox(1,d) = k;
            xbox(2,d) = k + 1;
        end
    else
        xbox = coder.nullcopy(zeros(2,ND,'like',varargin{1}));
        for d = coder.unroll(1:ND)
            k = coder.internal.bsearch(varargin{d},Xq(d));
            voffset = voffset + (k - 1)*stride(d);
            xbox(1,d) = varargin{d}(k);
            xbox(2,d) = varargin{d}(k + 1);
        end
    end
    % Evaluate the "local" interpolation rule over the bounding box, with the
    % reduced array of values V at each vertex.
    vbox = V(vbbidx + voffset);
    Vqk = interpf(ND,AUTOGRID,xbox,vbox,Xq);
end

%--------------------------------------------------------------------------

function vq = scalarMultilinearInterp(ND,AUTOGRID,x,v,xq)
% Let N be eml_ndims(v) and suppose that 1-D interpolation is excluded. x
% is an 2-by-N matrix of nearest grid neighbors bracketing the point xq,
% i.e. x defines an N-D bounding box for xq. v is an array with pow2(N)
% elements, laid out like a 2-by-2-by-2-... array, values at the grid
% points defined by x. To illustrate, if N = 4, v(1) is the value at x(1,:)
% and v(11) = v(1,2,1,2) is the value at (x(1,1),x(2,2),x(1,3),x(2,4)).
% xq is the point to interpolate, where, unless extrapolating,
% x(1,k) <= xq(k) <= x(2,k).  If AUTOGRID is true, x(1,k) == k and x(2,k)
% == k + 1.
coder.inline('always');
coder.internal.prefer_const(ND,AUTOGRID);
ONE = coder.internal.indexInt(1);
if isa(x,'single') && isa(v,'single') && isa(xq,'single')
    % Use single precision only.
    % Note that AUTOGRID effectively implies a double precision grid since
    % x belongs to coder.internal.indexIntClass in that case.
    RATIOCLASS = 'single';
else
    RATIOCLASS = 'double';
end
for d = 1:ND
    m = eml_lshift(ONE,ND - d);
    if xq(d) == x(1,d)
        for k = 1:m
            v(k) = v(2*k - 1);
        end
    elseif xq(d) == x(2,d)
        for k = 1:m
            v(k) = v(2*k);
        end
    else
        xqd = cast(xq(d),RATIOCLASS);
        x1d = cast(x(1,d),RATIOCLASS);
        if AUTOGRID
            r = xqd - x1d;
        else
            x2d = cast(x(2,d),RATIOCLASS);
            r = (xqd - x1d)/(x2d - x1d);
        end
        onemr = 1 - r;
        for k = 1:m
            v(k) = onemr*v(2*k - 1) + r*v(2*k);
        end
    end
end
vq = v(1);

%--------------------------------------------------------------------------

function vq = scalarNearestInterp(ND,AUTOGRID,x,v,xq)
% Let N be eml_ndims(v) and suppose that 1-D interpolation is excluded. x
% is an 2-by-N matrix of nearest grid neighbors bracketing the point xq,
% i.e. x defines an N-D bounding box for xq. v is an array with pow2(N)
% elements, laid out like a 2-by-2-by-2-... array, values at the grid
% points defined by x. To illustrate, if N = 4, v(1) is the value at x(1,:)
% and v(11) = v(1,2,1,2) is the value at (x(1,1),x(2,2),x(1,3),x(2,4)).
% xq is the point to interpolate, where, unless extrapolating,
% x(1,k) <= xq(k) <= x(2,k).
coder.inline('always');
coder.internal.prefer_const(ND,AUTOGRID);
ONE = coder.internal.indexInt(1);
% Construct the linear index. For example, if we settled on v(1,2,1,2)
% being the correct value to return, we will have arrived at it by
% constructing the linear index into v as follows:
% idx = 1 + (1 - 1)*0 + (2 - 1)*2 + (1 - 1)*4 + (2 - 1)*8 = 11.
% Consequently, we return v(11).
idx = ONE;
for d = 1:ND
    if abs(xq(d) - x(2,d)) <= abs(xq(d) - x(1,d))
        idx = idx + eml_lshift(ONE,d - ONE);
    end
end
vq = v(idx);

%--------------------------------------------------------------------------

function m = StringToMethodID(method)
% Convert a constant string input to a method ID.
coder.internal.prefer_const(method);
coder.internal.assert(coder.internal.isConst(method), ...
    'Coder:toolbox:MethodMustBeConstant');
n = coder.const(eml_max(1,strlength(method)));
islinear = coder.const(strncmpi(method,'linear',n));
isnearest = coder.const(strncmpi(method,'nearest',n));
iscubic = coder.const(strncmpi(method,'cubic',n));
isspline = coder.const(strncmpi(method,'spline',n));
if isnearest
    m = NEAREST;
elseif iscubic
    m = CUBIC;
elseif isspline
    m = SPLINE;
elseif islinear
    m = LINEAR;
else
    coder.internal.assert(false, ...
        'MATLAB:mathcgeo_catalog:BadInterpTypeErrId');
end

%--------------------------------------------------------------------------

function p = defaultEXTRAPp(METHOD)
if METHOD == SPLINE
    p = true;
else
    p = false;
end

%--------------------------------------------------------------------------

function [NARG,METHOD,EXTRAPp] = methodAndExtrapP(varargin)
% Return the number of grid/data arguments, the interpolation method, and
% whether or not to extrapolate.
coder.inline('always');
coder.internal.prefer_const(varargin);
if nargin >= 2 && coder.internal.isCharOrScalarString(varargin{end - 1})
    coder.internal.assert(coder.internal.isConst(varargin{end - 1}), ...
        'Coder:toolbox:MethodMustBeConstant');
    NARG = coder.internal.indexInt(nargin - 2);
    METHOD = coder.const(StringToMethodID(varargin{end - 1}));
    EXTRAPp = false;
elseif coder.internal.isCharOrScalarString(varargin{end})
    coder.internal.assert(coder.internal.isConst(varargin{end}), ...
        'Coder:toolbox:MethodMustBeConstant');
    NARG = coder.internal.indexInt(nargin - 1);
    METHOD = coder.const(StringToMethodID(varargin{end}));
    EXTRAPp = defaultEXTRAPp(METHOD);
else
    NARG = coder.internal.indexInt(nargin);
    METHOD = LINEAR;
    EXTRAPp = defaultEXTRAPp(METHOD);
end

%--------------------------------------------------------------------------

function [idx,stride] = BoundingBoxIndexArrays(V)
% Construct the index vectors needed to extract the elements of V for a
% bounding box. For example, if the bounding box starts at the V(2,5,1,3),
% then the value array for the bounding box defined by the points with
% coodinates x1(2:3),x2(5:6),x3(1:2),x4(3:4) is V(2:3,5:6,1:2,3:4). These
% values can be obtained in the form of a column vector as V(idx + offset),
% where offset = sub2ind(size(V),2,5,1,3) = 1 +
% (2-1)*stride(1) + (5-1)*stride(2) + (1-1)*stride(3) + (3-1)*stride(4).
ONE = coder.internal.indexInt(1);
ND = coder.internal.indexInt(eml_ndims(V));
n = coder.const(eml_lshift(ONE,ND));
idx = ones(n,1,coder.internal.indexIntClass);
stride = zeros(ND,1,coder.internal.indexIntClass);
s = ONE;
stride(1) = ONE;
for j = 2:2:n
    idx(j) = idx(j) + 1;
end
for k = 2:ND
    s = s*coder.internal.indexInt(size(V,k - 1));
    stride(k) = s;
    period = eml_lshift(ONE,k);
    halfperiod = eml_rshift(period,ONE);
    nperiods = eml_rshift(n,k);
    for j = 1:nperiods
        offset = (j - 1)*period + halfperiod;
        for i = 1:halfperiod
            idx(offset + i) = idx(offset + i) + s;
        end
    end
end

%--------------------------------------------------------------------------

function v = ndslice(k,varargin)
% Create an n-vector representing the k-th point in N-D space. TYPE is the
% an example of the return type. VARARGIN is a set of N vectors,
% VARARGIN{K} is the vector of abscissae in the K-th dimension.
coder.inline('always');
n = coder.const(numel(varargin));
v = coder.nullcopy(zeros(1,n,'like',varargin{1}));
for j = coder.unroll(1:n)
    v(j) = varargin{j}(k);
end

%--------------------------------------------------------------------------

function idx = autogridLookup(x,maxidx)
% Substitute for bsearch when the grid is automatic.
coder.inline('always')
if x <= 1
    idx = coder.internal.indexInt(1);
elseif x <= maxidx
    idx = coder.internal.indexInt(floor(x));
else
    idx = coder.internal.indexInt(maxidx);
end

%--------------------------------------------------------------------------

function p = allVectors(varargin)
% Returns true iff all inputs are known to be vectors when generating code.
p = true;
for k = coder.unroll(1:nargin)
    p = p && coder.const( ...
        coder.internal.isConst(isvector(varargin{k})) && ...
        isvector(varargin{k}));
end

%--------------------------------------------------------------------------

function p = isOutOfRange(ND,Xq,xmin,xmax)
% Determine whether the ND point Xq is out of range.
coder.internal.prefer_const(ND);
p = false;
for d = 1:ND
    p = p || ~(Xq(d) >= xmin(d) && Xq(d) <= xmax(d));
end

%--------------------------------------------------------------------------

function [Vidx,ND,AUTOGRID] = findVArg(varargin)
% Find the index of V in varargin, the dimensionality of the problem, and
% whether or not a grid is supplied.
narginchk(1,inf);
ONE = coder.internal.indexInt(1);
TWO = coder.internal.indexInt(2);
if nargin < 3
    Vidx = ONE;
    vectorV = coder.const( ...
        coder.internal.isConst(isvector(varargin{Vidx})) && ...
        isvector(varargin{Vidx}));
    if vectorV
        ND = ONE;
    else
        ND = coder.internal.indexInt(eml_ndims(varargin{Vidx}));
    end
    AUTOGRID = true;
elseif nargin == 3
    if coder.const(coder.internal.isConst(isvector(varargin{1})) && ...
            isvector(varargin{1}))
        % This is a 1-D problem.
        Vidx = TWO;
        ND = ONE;
        AUTOGRID = false;
    else
        % This is a 2-D problem.
        Vidx = ONE;
        ND = TWO;
        AUTOGRID = true;
    end
else
    NARG = coder.internal.indexInt(nargin);
    NNONV = coder.const(NARG - ONE);
    NDV1 = coder.internal.indexInt(eml_ndims(varargin{1}));
    % It might be possible to fool this algorithm if eml_ndims(varargin{1})
    % is greater than ndims(varargin{1}), but we need this function to
    % constant-fold, so we can't use ndims.
    if coder.const(NNONV > NDV1)
        ND = coder.const(eml_rshift(NNONV,ONE));
        Vidx = coder.const(ND + ONE);
        AUTOGRID = false;
    else
        Vidx = ONE;
        ND = NNONV;
        AUTOGRID = true;
    end
end

%--------------------------------------------------------------------------

function m = LINEAR
coder.inline('always');
m = uint8(0);

function m = NEAREST
coder.inline('always');
m = uint8(1);

function m = CUBIC
coder.inline('always');
m = uint8(3);

function m = SPLINE
coder.inline('always');
m = uint8(4);

%--------------------------------------------------------------------------
