function Vout = interp1(varargin)
%INTERP1 1-D interpolation (table lookup) for codistributed
%   Vq = interp1(V,Xq)
%   Vq = interp1(X,V,Xq)
%   Vq = interp1(V,Xq,METHOD)
%   Vq = interp1(X,V,Xq,METHOD)
%   Vq = interp1(V,Xq,METHOD,'extrap')
%   Vq = interp1(X,V,Xq,METHOD,'extrap')
%   Vq = interp1(V,Xq,METHOD,EXTRAPVAL)
%   Vq = interp1(X,V,Xq,METHOD,EXTRAPVAL)
%   
%   See also INTERP1, CODISTRIBUTED.
%   


%   Copyright 2012-2022 The MathWorks, Inc.

narginchk(2, 5);
distributedutil.CodistParser.verifyNonCodistributedInputs(varargin);
ndataarg = nargin;  % Number of X, V, Xq args.  Init to nargin and reduce.

% Default values
noX = false;
X = [];  % To satisfy assertion - it won't be used.
method = 'linear';
extrapval = 'none';

if nargin == 2 && isfloat(varargin{2})
    % INTERP1(V, Xq)
    V = varargin{1};
    Xq = varargin{2};
    noX = true;
    if matlab.internal.datatypes.isScalarText(Xq)
        error(message('MATLAB:interp1:nargin'));
    end
elseif nargin == 3 && isfloat(varargin{3})
    % INTERP1(X, V, Xq);
    X = varargin{1};
    V = varargin{2};
    Xq = varargin{3};
    if matlab.internal.datatypes.isScalarText(V)
        error(message('MATLAB:interp1:nargin'));
    end
else
    % INTERP1(X, V, Xq, ___)
    % INTERP1(V, Xq, ___)
    X = varargin{1};
    V = varargin{2};
    last = distributedutil.CodistParser.gatherIfCodistributed(varargin{nargin});
    penultimate = distributedutil.CodistParser.gatherIfCodistributed(varargin{nargin-1});
    [X, V, method, extrapval, ndataarg] = ...
        iParseInputs(X, V, penultimate, last, nargin, ndataarg);
    if ndataarg == 2
        Xq = varargin{2};
        noX = true;
    elseif ndataarg == 3
        Xq = varargin{3};
    else
        error(message('MATLAB:interp1:nargin'));
    end
end

% Interpolation must be done on the first dimension.  If either
% X or V is a vector, make it a column vector.
if ~isscalar(V) && isrow(V)
    V = V.';
end
if isrow(X)
    X = X.';
end

% Check variables before proceeding.
try
    [X, V, Xq] = iValidateData(X, V, Xq, noX);
    method = iValidateMethod(method);
    if strcmp(method, 'pp')
        error(message('parallel:array:Interp1UnsupportedppOption', 'codistributed'))
    end
catch err
    throw(err)
end

if any(strcmp(method, {'spline', 'pchip', 'makima'}))
    % Check V for NaNs and remove the rows containing them
    ndimV = ndims(V);
    nanV = find(sum(isnan(V), 2:ndimV)); % rows containing NaN
    if any(nanV)
        warning(message('MATLAB:interp1:NaNstrip'));
        V = deleteSlices(V, nanV, 1);
        X = deleteSlices(X, nanV, 1);
    end
end

szV = size(V);
szX = [szV(1) 1];
distGrid = iscodistributed(X) || iscodistributed(V);
distQuery = iscodistributed(Xq);

if distGrid
    % Create local GriddedInterpolant with padding to evaluate queries
    F = interpLocalGriddedInterpolant({X}, V, method, method, noX, szX(1), 1);
end

if distGrid && distQuery
    % Both grid and query are distributed
    Vout = iInterpDistGridDistQuery(F, X, V, Xq, method, extrapval, noX, szX, 1);

elseif distGrid
    % Only grid is distributed
    Vout = iInterpDistGridReplQuery(F, X, V, Xq, method, extrapval, noX, szX, 1);
    % If Vout is replicated then distribute it now.
    if ~iscodistributed(Vout)
        Vout = codistributed.pConstructFromReplicated(Vout);
    end

elseif distQuery
    % Only query is distributed
    Vout = iInterpReplGridDistQuery(X, V, Xq, method, extrapval, noX);

else
    % No data input was distributed. Gather and call in-memory.
    [varargin{:}] = gather(varargin{:});
    Vout = interp1(varargin{:});
end
end

%% Implementation for distributed grid, replicated query
function Vout = iInterpDistGridReplQuery(F, X, V, Xq, method, extrapval, noX, szX, dimInterp)
% Distributed grid, replicated query. Result will also be replicated since
% Vout must have same size and distribution as Xq.
if iscodistributed(V)
    codistributed.pVerifyUsing1d('interp1', V);
    codV = getCodistributor(V);
    if codV.Dimension ~= dimInterp
        localV = getLocalPart(V);
        if isempty(localV)
            % The worker has no values so return a compatibly sized empty
            localVout = reshape(localV, [numel(Xq) size(localV,2:ndims(localV))]);
        else
            localVout = iEvalGriddedInterpolant(F, Xq, extrapval);
        end
        szV = size(V);
        if isvector(Xq) && isvector(V)
            Vout = spmdBroadcast(find(codV.Partition), localVout);
            return
        elseif isvector(Xq)  % V is not a vector
            szLocalV = size(getLocalPart(V));
            localVout = reshape(localVout, [length(Xq), szLocalV(2:end)]);
            szVout = [length(Xq) szV(2:end)];
            distDim = codV.Dimension;
        else
            % Need to re-insert extra dimensions of Xq
            newLocalSz = [size(Xq) size(localVout, 2:ndims(localVout))];
            localVout = reshape(localVout, newLocalSz);
            szVout = [size(Xq) szV(2:end)];
            distDim = codV.Dimension + ndims(Xq) - 1;
        end
        VoutCodist = codistributor1d(distDim, codV.Partition, szVout);
        Vout = codistributed.pDoBuildFromLocalPart(localVout, VoutCodist);
        return

    elseif nnz(codV.Partition) == 1
        % V is distributed in the interpolation dimension, but all of V is on one worker.
        thisWorkerHasV = codV.Partition(spmdIndex);
        if thisWorkerHasV
            localVout = iEvalGriddedInterpolant(F, Xq, extrapval);
            if size(Xq, 2) > 1 && ~isrow(Xq) && ~ismatrix(localVout)
                % Need to re-insert extra dimensions of Xq
                newLocalSz = [size(Xq) size(localVout, 2:ndims(localVout))];
                localVout = reshape(localVout, newLocalSz);
            elseif ~iscolumn(Xq) && isvector(V)
                localVout = reshape(localVout, size(Xq));
            end
        else
            localVout = [];
        end

        % Broadcast the result to create replicated array
        Vout = spmdBroadcast(find(codV.Partition), localVout);
        return

    else % V is distributed in the interpolation dimension with some values
        % on different workers.
        codX = codistributor1d(dimInterp, codV.Partition, szX);
        if noX
            X = 1:max(szX);
            if dimInterp == 1
                X = X';
            end
            X = codistributed(X, codX);
        elseif ~iscodistributed(X)
            X = codistributed(X, codX);
        end
    end
end

if iscodistributed(X) && ~iscodistributed(V)
    % Check whether X is distributed in a scalar dimension
    codX = getCodistributor(X);
    if nnz(codX.Partition) == 1
        % All of X is on one worker and V is not distributed. That worker
        % can do all the work.
        if codX.Partition(spmdIndex)>0
            localVout = iEvalGriddedInterpolant(F, Xq, extrapval);
            if ~isvector(Xq)
                % Need to re-insert extra dimensions of Xq
                newLocalSz = [size(Xq) size(localVout, 2:ndims(localVout))];
                localVout = reshape(localVout, newLocalSz);
            end
        else
            localVout = [];
        end
        Vout = spmdBroadcast(find(codX.Partition), localVout);
        return;
    end

    % Make sure X and V have correct orientation and matching distribution.
    % Note that if V is not a vector it is always distributed in dimension
    % 1 since its first dimenion must match length(X).
    if isvector(V)
        V = reshape(V,size(X));
        codV = codistributor1d(codX.Dimension, codX.Partition, size(V));
    else
        codV = codistributor1d(1, codX.Partition, size(V));
    end
    V = codistributed(V, codV);
end

% By this point both X and V must be codistributed, and will have matching
% distribution schemes.
assert(iscodistributed(X) && iscodistributed(V), ...
    "Should only get here if both X and V are distributed.");

if any(strcmpi(method, {'pchip','cubic','v5cubic','makima'}))
    % Cubic-based methods may need two points on either side.
    padSz = 2;
else
    % Linear or nearest-neighbor methods need only one point on either side.
    padSz = 1;
end

% Figure out what Xq's to evaluate here
if isempty(Xq)
    Vout = [];
    return
end
partitionX = codX.Partition;
lastIndex = cumsum(partitionX);  % Index of last X on each worker
lengthX = lastIndex(end);
if ~isempty(F)
    localX = F.GridVectors{1};
else
    localX = [];
end
lowBreak = -Inf;
highBreak = Inf;
if spmdIndex>1 && lastIndex(spmdIndex) >= padSz && ...
        lastIndex(spmdIndex-1) <= lengthX - padSz && ...
        lastIndex(spmdIndex-1) >= padSz
    lowBreak = (localX(padSz)+localX(padSz+1))/2;
elseif spmdIndex>1 && lastIndex(spmdIndex-1) > lengthX - padSz
    lowBreak = Inf;
end
Breaks = spmdCat(lowBreak);
if spmdIndex < spmdSize
    highBreak = Breaks(spmdIndex+1);
end
origSzXq = size(Xq);
if ~isvector(Xq)  % Need to stretch out matrix for resolving query
    Xq = Xq(:).'; 
elseif ~isvector(V)
    origSzXq = numel(Xq);
end
queriesOnThisWorker = (Xq > lowBreak) & (Xq <= highBreak);

% Work out the output size
szV = size(V);
if isvector(V)
    szVout = size(Xq);
else
    szVout = [numel(Xq) szV(2:end)];
end

localXq = Xq(queriesOnThisWorker);
localVout = iEvalGriddedInterpolant(F, localXq, extrapval);
Vout = nan(szVout, "like", localVout);
dimCombine = 1;
if numel(queriesOnThisWorker) ~= szVout(1) && numel(queriesOnThisWorker) == szVout(2)
    dimCombine = 2;
end

Vout = assignSlices(Vout, queriesOnThisWorker, localVout, dimCombine);
% Each worker will have filled just the elements it had the grid for. We
% now need to combine results across all workers.
Vout = combineSlicesAcrossWorkers(Vout, queriesOnThisWorker, 1);

% Reshape because Xq was stretched into a column vector.
Vout = reshape(Vout, [origSzXq, szV(2:end)]);

% NB Vout is replicated to match Xq

end

%% Implementation for replicated grid, distributed query
function Vout = iInterpReplGridDistQuery(X, V, Xq, method, extrapval, noX)
codistributed.pVerifyUsing1d('interp1', Xq);
codXq = getCodistributor(Xq);
localXq = getLocalPart(Xq);
localVout = callLocalInterp({X}, V, {localXq}, method, extrapval, noX);
if codXq.Partition(spmdIndex) == 1 && codXq.Dimension == ndims(Xq)
    % localVout is missing a dimension.
    szV = size(V);
    if isvector(Xq)
        localVout = reshape(localVout, [1, szV(2:end)]);
    elseif ismatrix(Xq)
        localVout = reshape(localVout, [size(localXq), szV(2:end)]);
    else
        localVout = reshape(localVout, [size(localXq), 1, szV(2:end)]);
    end
end
partCodXq = codXq.Partition;

% Determine size of output
if isvector(V)
    VoutCodist = codXq;
else
    szV = size(V);
    if isvector(Xq)
        szVout = [length(Xq) szV(2:end)];
        XqDistDim = 1;
        if sum(partCodXq) == 1  % All of Xq is on worker 1
            partCodXq(1) = length(Xq);
        end
    else
        szVout = [size(Xq) szV(2:end)];
        XqDistDim = codXq.Dimension;
    end
    VoutCodist = codistributor1d(XqDistDim, partCodXq, szVout);
end
Vout = codistributed.pDoBuildFromLocalPart(localVout, VoutCodist);
end

%% Implementation for distributed grid, distributed query
function Vout = iInterpDistGridDistQuery(F, X, V, Xq, method, extrapval, noX, szX, dimInterp)
% For every worker that contains a part of Xq, the algorithm goes through it by
% turn and then replicates the data across all workers. Following this it
% uses the distributed grid, replicated query algorithm to generate the
% result for that query across all workers.
codXq = getCodistributor(Xq);
partXq = codXq.Partition;
localVout = [];
for ii = 1:spmdSize
    if partXq(ii) == 0
        % Local part empty so skip
        continue;
    end
    % Replicate one local part across all workers
    if ii == spmdIndex
        sharedXq = spmdBroadcast(ii,getLocalPart(Xq));
    else
        sharedXq = spmdBroadcast(ii);
    end
    % Have all workers process the replicated query with the distributed
    % grid.
    replVout = iInterpDistGridReplQuery(F, X, V, sharedXq, method, extrapval, noX, szX, dimInterp);
    % Might have come back distributed if V was distributed, so gather onto
    % the worker who's local part we are processing.
    if iscodistributed(replVout)
        replVout = gather(replVout, ii);
    end
    % If this was our local part, store it
    if ii == spmdIndex
        localVout = replVout;
        localXq = sharedXq;
    end
end

% If V is a vector the output size and distribution matches Xq. If a matrix
% then the dimensions need combining.
szV = size(V, 2:ndims(V));
if isvector(V)
    szVout = size(Xq);
    voutDistDim = codXq.Dimension;
elseif isvector(Xq)
    szVout = [numel(Xq), szV];
    voutDistDim = 1;
    if sum(partXq) == 1  % All of Xq is on one worker
        partXq(1) = length(Xq);
    end
else
    szVout = [size(Xq), szV];
    voutDistDim = codXq.Dimension;
    if partXq(spmdIndex) == 1 && voutDistDim == ndims(Xq)
        % localVout is missing a dimension.
        if isvector(localXq)
            localVout = reshape(localVout, [size(localXq), szV]);
        else
            localVout = reshape(localVout, [size(localXq), 1, szV]);
        end
    end
end
VoutCodist = codistributor1d(voutDistDim, partXq, szVout);
Vout = codistributed.pDoBuildFromLocalPart(localVout, VoutCodist);
end

%% Helpers
%-------------------------------------------------------------------------%
function [X,V,method, extrapval, ndataarg] = ...
    iParseInputs(X, V, penultimate, last, ninputs, ndataarg)
% Default values
method = 'linear';
extrapval = 'none';

if matlab.internal.datatypes.isScalarText(last)
    if strcmp(last,'extrap')
        if (ninputs ~= 4 && ninputs ~= 5)
            error(message('MATLAB:interp1:nargin'));
        end
        if ~isempty(penultimate) && ~matlab.internal.datatypes.isScalarText(penultimate)
            error(message('MATLAB:interp1:ExtrapNoMethod'));
        end
        method = iValidateMethod(penultimate);
        ndataarg = ninputs-2;
        if strcmp(method,'cubic')
            warning(message('MATLAB:interp1:NoExtrapForV5cubic'))
        else
            extrapval = 'extrap';
        end
    else
        if matlab.internal.datatypes.isScalarText(penultimate)
            if strcmp(last, 'pp')
                if ~strcmp(penultimate, 'linear')
                    error(message('parallel:array:Interp1UnsupportedppOption', 'codistributed'))
                end
            else
                error(message('MATLAB:interp1:InvalidSpecPPExtrap'))
            end
        end
        method = iValidateMethod(last);
        needextrapval = ~any(strcmp(method,{'spline','pchip','makima'}));
        if ~needextrapval
            extrapval = 'extrap';
        end
        ndataarg = ninputs-1;
    end
elseif isscalar(last)
    extrapval = last;
    if matlab.internal.datatypes.isScalarText(penultimate)
        ndataarg = ninputs-2;
        method = penultimate;
    elseif isempty(penultimate) && (ninputs == 4 || ninputs == 5)
        % default method via []
        ndataarg = ninputs-2;
    elseif isfloat(penultimate) && (ninputs == 4 || ninputs == 5)
        error(message('MATLAB:interp1:nargin'))
    end
elseif isempty(last)
    % This is potentially ambiguous, the assumed intent is case I
    % I)    X,  V, []  Empty query
    % II)   V, [], []  Empty query and empty method
    % III)  V, Xq, []  Empty method
    if ninputs ~= 3
        ndataarg = ninputs-1;
    end
end
if ndataarg < 3
    V = X;
    X = [];
end
end

function [X,V,Xq] = iValidateData(X, V, Xq, noX)
szV = size(V);
vecV = isvector(V);

if issparse(V)
    error(message('MATLAB:griddedInterpolant:SparseInterpValuesErrId'))
elseif isempty(V) && szV(2)>0
    error(message('MATLAB:griddedInterpolant:DegenerateGridErrId'))
elseif isempty(X) && isempty(V) && isempty(Xq) && ~noX
    error(message('MATLAB:interp1:Xvector'))
elseif ~isfloat(X)
    error(message('MATLAB:interp1:Xnumeric'))
elseif ~noX && ~isvector(X)
    error(message('MATLAB:interp1:Xvector'))
elseif any(isinf(X)) || any(isnan(X))
    error(message('MATLAB:griddedInterpolant:NonFiniteInputPtsErrId'))
elseif ~isreal(X)
    error(message('MATLAB:griddedInterpolant:ComplexDataPointErrId'))
elseif ~isfloat(V)
    error(message('MATLAB:interp1:NonFloatValues'))
elseif ~noX && vecV && length(V) ~= length(X)
    error(message('MATLAB:interp1:YVectorInvalidNumRows'))
elseif numel(V)<2
    error(message('MATLAB:griddedInterpolant:DegenerateGridErrId'))
elseif iscell(Xq)
    error(message('parallel:distributed:InterpInputsNotNumeric', 3))
elseif isa(Xq, 'function_handle')
    error(message('MATLAB:lang:colonInput'))
elseif ~isfloat(Xq)
    error(message('MATLAB:griddedInterpolant:NonFloatInpPtsErrId'))
elseif ~isreal(Xq)
    error(message('MATLAB:griddedInterpolant:ComplexDataPointErrId'))
end

% Enforce sorting and check unique
[unqX, idx] = unique(X);
if ~isempty(X) && length(X) > length(unqX)
    error(message('MATLAB:griddedInterpolant:NonUniqueCompVecsPtsErrId'))
end
if ~issorted(X)
    X = unqX;
    V = extractSlices(V, idx, 1);
end

end

%-------------------------------------------------------------------------%
%     'nearest'  - nearest neighbor interpolation
%     'next'     - next neighbor interpolation
%     'previous' - previous neighbor interpolation
%     'linear'   - linear interpolation
%     'spline'   - piecewise cubic spline interpolation (SPLINE)
%     'pchip'    - shape-preserving piecewise cubic interpolation
%     'cubic'    - same as 'pchip'
%     'v5cubic'  - the cubic interpolation from MATLAB 5
%     'makima'   - modified Akima cubic interpolation
function methodname = iValidateMethod(method_in)
method = char(method_in); % string support
if isempty(method)
    methodname = 'linear';
else
    if method(1) == '*'
        method(1) = [];
    end
    if matlab.internal.math.partialMatch(method, 'linear')
        methodname = 'linear';
    elseif matlab.internal.math.partialMatch(method, 'makima')
        methodname = 'makima';
    elseif matlab.internal.math.partialMatch(method, 'spline')
        methodname = 'spline';
    elseif matlab.internal.math.partialMatch(method, 'next', 3)
        methodname = 'next';
    elseif matlab.internal.math.partialMatch(method, 'nearest')
        methodname = 'nearest';
    elseif matlab.internal.math.partialMatch(method, 'pchip')
        methodname = 'pchip';
    elseif matlab.internal.math.partialMatch(method, 'previous', 2)
        methodname = 'previous';
    elseif matlab.internal.math.partialMatch(method, 'cubic')
        methodname = 'cubic';
    elseif matlab.internal.math.partialMatch(method, 'pp', 2)
        methodname = 'pp';
    elseif matlab.internal.math.partialMatch(method, 'v5cubic')
        methodname = 'cubic';
    else
        error(message('MATLAB:interp1:InvalidMethod'));
    end
end
end

function localVout = iEvalGriddedInterpolant(F, Xq, extrapval)
% Use a gridded interpolant to calculate interpolated values on a worker.
% localVout has size [numel(Xq) size(localV,2:ndims(localV))].  If
% Xq is not a vector, the additional dimensions must be inserted in
% localVout.

if isempty(F)
    localVout = [];
    return
end
X = F.GridVectors{1};
V = F.Values;
szV = size(F.Values);
Xqcol = Xq(:);
doExtrapPostprocessing = ~strcmp(extrapval, 'extrap');
if doExtrapPostprocessing
    if ~isempty(Xqcol) && isfloat(Xqcol) && isreal(Xqcol)
        extrapMask = Xqcol < X(1) | Xqcol > X(end);
        if any(extrapMask)
            Xqcol = Xq(~extrapMask);
        else
            doExtrapPostprocessing = false;
        end
    else
        doExtrapPostprocessing = false;
    end
end
if isa(F, 'griddedInterpolant')
    localVout = F(Xqcol);
else
    localVout = ppval(F, Xqcol);
    newOrder = [F.dims(end) F.dims(1:end-1)];
    localVout = permute(localVout, newOrder);
end
if doExtrapPostprocessing
    if strcmp(extrapval, 'none')
        if ~isreal(V)
            extrapval = NaN + 1i*NaN;
        else
            extrapval = NaN;
        end
    end
    % localVout is too small because elems of Xqcol were removed.
    newLocalVout = zeros([numel(Xqcol), szV(2:end)], 'like', Xq);
    newLocalVout = assignSlices(newLocalVout, ~extrapMask, localVout, 1);
    newLocalVout = assignSlices(newLocalVout, extrapMask, extrapval, 1);
    if ~isscalar(localVout) && isrow(localVout)
        localVout = newLocalVout.';
    else
        localVout = newLocalVout;
    end
end
if isrow(Xq) && iscolumn(localVout)
    localVout = localVout.';
end
localVout = cast(localVout, superiorfloat(X, V, Xq));
end
