function Vout = interpn(varargin)
%INTERPN N-D interpolation (table lookup) for codistributed
%   Vq = interpn(V,X1q,X2q,X3q,...)
%   Vq = interpn(X1,X2,X3,...,V,X1q,X2q,X3q,...)
%   Vq = interpn(V,X1q,X2q,X3q,...,METHOD)
%   Vq = interpn(X1,X2,X3,...,V,X1q,X2q,X3q,...,METHOD)
%   Vq = interpn(V,X1q,X2q,X3q,...,METHOD,EXTRAPVAL)
%   Vq = interpn(X1,X2,X3,...,V,X1q,X2q,X3q,...,,METHOD,EXTRAPVAL)
%   Vq = interpn(V,NTIMES)
%   Vq = interpn(V) % same as interpn(V, 1)
%   
%   The value argument V must be a double or single N-D array, real or complex. 
%   
%   The grid arguments X1, X2, X3, ... must have the following characteristics:
%     * Same type, double or single.
%     * Finite vectors or N-D arrays with increasing and non-repeating elements in corresponding dimension. 
%     * When X1, X2, X3, ... are N-D arrays, they must align with cartesian axes, as if they were 
%       produced by ndgrid.
%     * Dimensions consistent with those of value V. 
%   
%   The query arguments X1q, X2q, X3q, ... must be vectors or arrays of the same type (double or single). 
%   If X1q, X2q, X3q, ... are arrays, they must be the same size. If they are vectors with different lengths, 
%   one of them must have different orientation.
%   
%   Example:
%   spmd
%      f = @(x,y,z,t) t.*exp(-x.^2 - y.^2 - z.^2);
%      [x,y,z,t] = ndgrid(-1:0.2:1,-1:0.2:1,-1:0.2:1,0:2:10);
%      v = f(x,y,z,t);
%      [xq,yq,zq,tq] = ndgrid(-1:0.05:1,-1:0.08:1,-1:0.05:1,0:0.5:10);
%      vq = interpn(x,y,z,t,codistributed(v),xq,yq,zq,tq);        
%   end
%   
%      % Visualize the function:
%      nframes = size(tq, 4);
%      for j = 1:nframes
%         slice(yq(:,:,:,j), xq(:,:,:,j), zq(:,:,:,j), vq(:,:,:,j),0,0,0);
%         caxis([0 10]);
%         M(j) = getframe;
%      end
%      movie(M); 
%   
%   See also INTERP1, INTERP2, INTERP3, NDGRID, SCATTEREDINTERPOLANT, CODISTRIBUTED.


%   Copyright 2013-2023 The MathWorks, Inc.

narginchk(1, inf);
distributedutil.CodistParser.verifyNonCodistributedInputs(varargin);
if (nargin >= 3 && matlab.internal.datatypes.isScalarText(varargin{end-1})) && ...
   (~isnumeric(varargin{end}) || ~isscalar(varargin{end}))
    error(message('MATLAB:interpn:InvalidExtrapval'))
end

% Parse the method and extrap val
[narg, method, extrapval] = methodAndExtrapval(varargin{:});
stripNaNsForCubics = method == "spline" || method == "makima";
if stripNaNsForCubics
    extrap = method;
else
    extrap = 'none';
end

% 'pchip' is not supported in N-D and we switch to 'linear'
if method == "pchip" && ~iIsDimensionOne(narg, varargin)
    warning(message("MATLAB:griddedInterpolant:Pchip1DonlyWarnId"))
    method = 'linear';
end

noX = false;
if narg == 1 || narg == 2
    %  interpn(V,NTIMES)
    V = varargin{1};
    szV = size(V);

    if isvector(V)
        ndimsV = 1;
        szV = length(V);
        if isrow(V)
            V = V.';
        end
    else
        ndimsV = ndims(V);
    end
    if narg == 1
        ntimes = 1;
    else
        ntimes = varargin{2}(1);
        if ~isscalar(ntimes) || ~isreal(ntimes)
            error(message('MATLAB:interpn:NtimesInvalid'))
        end
        ntimes = distributedutil.CodistParser.gatherIfCodistributed(ntimes);
        ntimes = floor(ntimes);
    end
    
    Xq = cell(1, ndimsV);
    if isinf(ntimes)
        Xq = cellfun(@(x) zeros(1, 0), Xq, 'UniformOutput', false);
    elseif isnan(ntimes)
        Xq = cellfun(@(x) NaN, Xq, 'UniformOutput', false);
    elseif ntimes < 0
        Xq = arrayfun(@(x) 1:1/(2^ntimes):szV(x), 1:ndimsV, 'UniformOutput', false);
    else
        Xq = arrayfun(@(x) linspace(1, szV(x), (szV(x)-1)*2^ntimes+1), 1:ndimsV,...
            'UniformOutput', false);
    end
    if ndimsV > 1
        Xq{1} = Xq{1}.';  % Flip orientation to denote compact grid.
    end

    szVout = prod(cellfun(@numel, Xq));
    % By default, the query is replicated.  If the query is really large,
    % this uses a lot of memory on all workers, so distributing the query
    % cuts the memory pressure.
    if szVout > 100000
        [~, longestVec] = max(cellfun(@numel, Xq));
        if iscodistributed(V)
            codV = getCodistributor(V);
            distDim = codV.Dimension;
        end
        if (~iscodistributed(V) || distDim > ndimsV)
            distDim = longestVec;
        end
        distDimXq = min(distDim, 2);
        Xq{distDim} = codistributed(Xq{distDim}, codistributor1d(distDimXq));
    end
    noX = true;
else
    arg1 = varargin{1};
    ndimsArg1 = ndims(arg1);
    if ~isvector(arg1) && narg == ndimsArg1 + 1
        %  interpn(V,X1q,X2q,X3q,...)
        V = arg1;
        Xq = varargin(2:narg);
        noX = true;
    elseif rem(narg, 2) == 1
        %  interpn(X1,X2,X3,...,V,X1q,X2q,X3q,...)
        mididx = ceil(narg/2);
        X = varargin(1:(mididx-1));
        V = varargin{mididx};
        Xq = varargin((mididx+1):narg);
        iCheckX(X);
        if ~isscalar(V) && isrow(V)
            V = V.';
        end
        if ~all(cellfun(@isvector, X))
            % At least one X is not a vector; make sure they define an ndgrid and
            % extract the grid vectors.
            [isMeshGrid, isNdGrid, exception] = checkGrid(X{:}, V);
            if exception == "ArgMismatch"
                sizeX = cell2mat(cellfun(@size, X, 'UniformOutput', false));
                if any(diff(sizeX), 'all')
                    % Grid arrays aren't the same size.
                    error(message('MATLAB:griddedInterpolant:BadGridErrId'))
                else  % One or more grid arrays doesn't match V.
                    error(message('MATLAB:griddedInterpolant:NumCoordsGridNdimsMismatchErrId'))
                end
            elseif exception == "DimensionMismatch"
                error(message('MATLAB:griddedInterpolant:DegenerateGridErrId'))
            elseif exception == "TypeMismatch"
                error(message('MATLAB:griddedInterpolant:NumCoordsGridNdimsMismatchErrId'));   
            elseif isMeshGrid
                if numel(X) == 2
                    error(message('MATLAB:griddedInterpolant:NdgridNotMeshgrid2DErrId'))
                elseif numel(X) == 3
                    error(message('MATLAB:griddedInterpolant:NdgridNotMeshgrid3DErrId'))
                else
                    error(message('MATLAB:griddedInterpolant:BadGridErrId'))
                end
            elseif ~isNdGrid % neither meshgrid nor ndgrid.
                error(message('MATLAB:griddedInterpolant:BadGridErrId'))
            end
            [X{:}, ~] = getGridVectors(false, isNdGrid, X{:}, V);
        else
            X = iValidateGridVectors(X, V);
        end
    else
        error(message('MATLAB:interpn:nargin'))
    end
end

% Every call passes through this section.
if isvector(V)
    nInterp = 1;
    if isrow(V)
        V = V.';
    end
else
    nInterp = ndims(V);
end

if noX
    [X{1:nInterp}] = ndgridVectors(V);
end

if nInterp == 1 && isrow(X{1})
    X{1} = X{1}.';
end

if all(method ~= ["linear", "nearest", "cubic", "makima", "spline", "pchip"])
    error(message('MATLAB:griddedInterpolant:BadInterpTypeErrId'))
end

[V, origvtype] = iConvertV(V, method, Xq);
[X{:}, V] = checkMonotonicGridVectors(X{:}, V);
if stripNaNsForCubics
    [X, V] = iStripNaNWrapper(X, V);
end

iscodistX = cellfun(@iscodistributed, X);
iscodistXq = cellfun(@iscodistributed, Xq);
isDistGrid = any(iscodistX) || iscodistributed(V);
isDistQuery = any(iscodistXq);

if isDistGrid
    % Create local griddedInterpolant with padding to evaluate queries
    szX = cellfun(@length, X);
    F = iMakeGriddedInterpolant(X, V, method, extrap, false, szX, 1:numel(X));
    if method == "cubic"
        if ~isempty(F)
            method = F.Method;
        end
        methods = spmdCat({method}, 1);
        if ~all(methods == "cubic")  % at least one was changed.
            [X, V] = iStripNaNWrapper(X, V);
            method = "spline";
            F = iMakeGriddedInterpolant(X, V, method, extrap, false, szX, 1:numel(X));
        end
    end
end

if isDistGrid && isDistQuery
    % Both grid and query are distributed
    Vout = iInterpDistGridDistQuery(F, X, V, Xq, method, extrapval);

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

elseif isDistQuery
    % 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 = interpn(varargin{:});
end

if origvtype ~= "double"
    Vout = distributedutil.syncOnError(@cast, Vout, origvtype);
end
end

%% Implementation for distributed grid, replicated query
function Vout = iInterpDistGridReplQuery(F, X, V, Xq, method, extrapval)
% Distributed grid, replicated query. Result will also be replicated since
% Vout must have same size and distribution as Xq and Yq.
isCodX = cellfun(@iscodistributed, X);
distDimX = find(isCodX, 1);
if iscodistributed(V)
    codistributed.pVerifyUsing1d('interpn', V);
    codV = getCodistributor(V);
    distDimV = codV.Dimension;
    if nnz(codV.Partition) == 1
        % V is distributed in an interpolation dimension, but all of V is on one worker.
        thisWorkerHasV = codV.Partition(spmdIndex);
        if thisWorkerHasV
            localVout = iEvalGriddedInterpolant(F, Xq, extrapval);
            if numel(Xq) == 1 && ~iscolumn(Xq{1}) && isvector(V)
                localVout = reshape(localVout, size(Xq{1}));
            end
        else
            localVout = [];
        end

        % Broadcast the result to create replicated array
        Vout = gather(localVout);
        return

    else % V is distributed in an interpolation dimension with some values
        % on different workers.
        codX = codistributor1d(1, codV.Partition, size(X{distDimV}));
        if ~iscodistributed(X(distDimV))
            X{distDimV} = codistributed(X{distDimV}, codistributor1d(1, codV.Partition));
            distDimX = distDimV;
        end
    end
end

if any(isCodX) && ~iscodistributed(V)
    % Check whether the X is distributed in a scalar dimension
    codX = getCodistributor(X{distDimX});
    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);
        else
            localVout = [];
        end
        Vout = gather(localVout);
        return
    end

    % Make sure X and V have correct orientation and matching distribution.
    V = codistributed(V, codistributor1d(distDimX, codX.Partition, size(V)));
elseif ~any(isCodX) && ~iscodistributed(V) 
    % Case where 'cubic' gets converted to 'spline' because of NaNs.
    localVout = iEvalGriddedInterpolant(F, Xq, extrapval);
    Vout = gather(localVout);
    return
end

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


if method == "makima"
    % makima in interpn needs three points for full accuracy.
    padSz = 3;
elseif method == "cubic" || method == "spline" || method == "pchip"
    % 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
partitionX = codX.Partition;
lastIndex = cumsum(partitionX);  % Index of last X on each worker
lengthX = lastIndex(end);
if ~isempty(F)
    localX = F.GridVectors{distDimX};
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

% If the query is two vectors with different orientation, they
% define a meshgrid.  Otherwise, Xq and Yq are used in pairs to
% define discrete sample points.
isXqScalar = all(cellfun(@isscalar, Xq));
isXqRow = cellfun(@isrow, Xq);
isXqCol = cellfun(@iscolumn, Xq);
isXqVec = all(isXqRow | isXqCol);
isMgVec = ~isXqScalar && isXqVec && (any(isXqRow) && any(isXqCol));
isXqEmpty = any(cellfun(@isempty, Xq));
hasXqNaN = any(cellfun(@(x) any(isnan(x), 'all'), Xq));
if ~isMgVec && ~isXqEmpty && ~hasXqNaN
    [~, isNdGrid, exception] = checkGrid(Xq{:}, Xq{1});
    if exception == "SizeMismatch"
        error(message('MATLAB:griddedInterpolant:InputMixSizeErrId'))
    end
    if isNdGrid
        [Xq{:}, ~] = getGridVectors(false, isNdGrid, Xq{:}, Xq{1}); % getGridVectors returns ndgrid vectors.
        Xq{2} = Xq{2}.';   %  getGridVectors returns column vectors
        isMgVec = true;
    end
end
nInterp = numel(Xq);
if isMgVec
    localSzVout = cellfun(@numel, Xq);
    distDimVout = distDimX;
else
    szXq = cell2mat(cellfun(@size, Xq, 'UniformOutput', false)');
    if nInterp > 1 && any(diff(szXq), 'all')
        error(message('MATLAB:griddedInterpolant:InputMixSizeErrId'))
    end
    szVout = size(Xq{1});
    % make them column vectors for evaluation - reshape at end.
    Xq = cellfun(@(x) x(:), Xq, 'UniformOutput', false);
    localSzVout = size(Xq{1});
    distDimVout = 1;
end
queriesOnThisWorker = (Xq{distDimX} >= lowBreak) & (Xq{distDimX} < highBreak);
if isMgVec
    localXq = Xq;
    localXq{distDimX} = Xq{distDimX}(queriesOnThisWorker);
else
    if hasXqNaN
        nanLoc = cell2mat(cellfun(@isnan, Xq, 'UniformOutput', false));
        nanLoc = reshape(nanLoc, nInterp, []);
        nanLoc = reshape(max(nanLoc), size(Xq{1}));
        queriesOnThisWorker = nanLoc | queriesOnThisWorker;
    end
    localXq = cellfun(@(x) x(queriesOnThisWorker), Xq, 'UniformOutput', false);
end
if method == "spline" && nInterp > 1
    localXqInDistDim = localXq{distDimX};
    [uLocalXqInDistDim, ~, idxIn] = unique(localXqInDistDim);
    notDistDimX = 1:nInterp;
    notDistDimX(distDimX) = [];
    if distDimX == 2
        localXq{distDimX} = cast(localX.', 'like', uLocalXqInDistDim);
    else
        localXq{distDimX} = cast(localX, 'like', uLocalXqInDistDim);
    end
    if isMgVec
        localXq(notDistDimX) = cellfun(@(x) x, Xq(notDistDimX), 'UniformOutput', false);
        if ~compactGridFormat(localXq{:})
            localXq{notDistDimX(1)} = localXq{notDistDimX(1)}.';
        end
    else % This block creates a set of parallel query vectors so that the griddedInterpolant
         % does a spline interpolation onto each point of the local part of the distributed
         % grid vector using the undistributed query points.
        origLocalXq = localXq;
        nX = numel(localX);
        nXq = numel(Xq{notDistDimX(1)});
        localXq{distDimX} = reshape(repmat(localXq{distDimX}, 1, nXq), nX*nXq, 1);
        [localXq(notDistDimX)] = arrayfun(@(x) reshape(repmat(Xq{x}.', nX, 1), nX*nXq, 1),...
            notDistDimX, 'UniformOutput', false);
        origExtrapval = extrapval;
        extrapval = [];
    end
end
localVout = iEvalGriddedInterpolant(F, localXq, extrapval);
if method == "spline" && ~isempty(localVout) && nInterp > 1
    if isMgVec
        Fsp = interpSpline(X{distDimX}, localXq{distDimX}, localVout, localVout, distDimX);
    else
        localVout = reshape(localVout, nX, nXq);
        Fsp = interpSpline(X{distDimX}, origLocalXq{distDimX}, localVout, localVout, 1);
    end        
    if ~isempty(uLocalXqInDistDim)
        if isMgVec % The query is a grid, so we calculate all combinations.
            if ~isrow(uLocalXqInDistDim)
                uLocalXqInDistDim = uLocalXqInDistDim.';
            end
            splineVout = ppval(Fsp, uLocalXqInDistDim);
            if ndims(splineVout) < nInterp
                splineVout = shiftdim(splineVout, -1);
            else
                splineVout = shiftdim(splineVout, nInterp-1);
            end
            splineVout = extractSlices(splineVout, idxIn, 1);
            localVout = ipermute(double(splineVout), Fsp.dims);
            localXq{distDimX} = localXqInDistDim;
        else  % The query is a set of points, so we calculate one result for each point.
            localVout = double(ppValDiag(Fsp, Xq{distDimX}));
            localXq = Xq;
            extrapval = origExtrapval;
        end
    else
        localVout = [];
    end            
    if ~isempty(extrapval)
        localVout = imposeExtrapval(localXq, F.GridVectors, localVout,...
            extrapval, isMgVec);
    end
end
Vout = nan(localSzVout, 'like', localVout);
if isempty(localVout)
    queriesOnThisWorker = false(size(queriesOnThisWorker));
end
if iscolumn(Vout) && method == "spline"
    Vout(queriesOnThisWorker) = localVout(queriesOnThisWorker);
else
    Vout = assignSlices(Vout, queriesOnThisWorker, localVout, distDimVout);
end

% 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, distDimVout);

% If Xq was reshaped, we need to put the result back into the right shape.
if ~isMgVec && ~isscalar(Xq{1})
    Vout = reshape(Vout, szVout);
end
if isempty(Vout)
    Vout = double.empty(zeros(1, nInterp));
end

% NB Vout is replicated to match Xq

end

%% Implementation for replicated grid, distributed query
function Vout = iInterpReplGridDistQuery(X, V, Xq, method, extrapval, noX)

[localXq, codXq, szVout, VoutDistDim] = iValidateDistQuery(Xq);

localVout = callLocalInterp(X, V, localXq, method, extrapval, noX);
% Build Vout from local parts.
VoutCodist = codistributor1d(VoutDistDim, codXq.Partition, szVout);
Vout = codistributed.pDoBuildFromLocalPart(localVout, VoutCodist);
end

%% Implementation for distributed grid, distributed query
function Vout = iInterpDistGridDistQuery(F, X, V, Xq, method, extrapval)
% 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.
[localXq, codXq, szVout, VoutDistDim] = iValidateDistQuery(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
    sharedXq = spmdBroadcast(ii, localXq);
    % Have all workers process the replicated query with the distributed
    % grid.
    replVout = iInterpDistGridReplQuery(F, X, V, sharedXq, method, extrapval);
    % Might have come back distributed if V was distributed, so gather onto
    % the worker whose 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;
    end
end

% If V is a vector the output size and distribution matches Xq. If a matrix
% then the dimensions need combining.
VoutCodist = codistributor1d(VoutDistDim, codXq.Partition, szVout);
Vout = codistributed.pDoBuildFromLocalPart(localVout, VoutCodist);
end

%% Helpers
%-------------------------------------------------------------------------%
function [V, origvtype] = iConvertV(V, method, Xq)
%ICONVERTV changes the type of V to double in most cases, but returns the
%original type in any case.  It also checks the validity of V and the Xq's.

origvtype = underlyingType(V);
isScalarXq = cellfun(@isscalar, Xq);
if ~isfloat(V) && (method == "nearest" || all(isScalarXq))
    V = double(V);
end
if ~isfloat(V)
    error(message('MATLAB:griddedInterpolant:NonFloatInpPtsErrId'))
end
isFloatXq = cellfun(@isfloat, Xq);
isRealXq = cellfun(@isreal, Xq);
if ~all(isFloatXq)
    error(message('MATLAB:griddedInterpolant:NonFloatInpPtsErrId'))
elseif ~all(isRealXq)
    error(message('MATLAB:griddedInterpolant:ComplexDataPointErrId'))
end
cellfun(@(x) codistributed.pVerifyUsing1d('interpn', x), Xq);
end

%-------------------------------------------------------------------------%
function [X, V] = iStripNaNWrapper(X, V)
%ISTRIPNANWRAPPER removes slices of V containing NaNs.

numvbefore = numel(V);
if isvector(V)
    ndimsbefore = 1;
else
    ndimsbefore = ndims(V);
end
[inan, jnan] = find(isnan(V));
inan = distributedutil.CodistParser.gatherIfCodistributed(inan);
jnan = distributedutil.CodistParser.gatherIfCodistributed(jnan);
ncolnan = length(unique(jnan));
nrownan = length(unique(inan));
if ncolnan == 0 && nrownan == 0
    return
end
[X{:}, V] = stripNaNsForSpline(X{:}, V);
numvafter = numel(V);
if numvbefore > numvafter
    warning(message('MATLAB:interpn:NaNstrip'))
end

if isempty(V)
    ndimsafter = 0;
elseif isvector(V)
    ndimsafter = 1;
else
    ndimsafter = ndims(V);
end

if ndimsbefore > ndimsafter
    error(message('MATLAB:interpn:NotEnoughPointsNanStrip'))
end
end

%-------------------------------------------------------------------------%
function gridVec = iValidateGridVectors(X, V)
% Make sure that only one grid vector is distributed and its distribution
% matches V if V is codistributed. Vectors are returned as column vectors.

szV = size(V);
if isa(V, 'function_handle')
    error(message('parallel:array:UnsupportedFunction', 'INTERPN', func2str(V)))
elseif isa(V, 'struct')
    error(message('MATLAB:invalidConversion', 'struct', underlyingType(V)))
elseif iscell(V)
    error(message('MATLAB:invalidConversion', 'cell', underlyingType(V)))
end
numelX = cellfun(@numel, X);
nInterp = numel(X);
if nInterp == 1  % If V is a vector, we are concerned only with its length.
    szV(2) = [];
end
if any(szV < 2)
    error(message('MATLAB:griddedInterpolant:DegenerateGridErrId'))
end
if any(numelX-szV)
    mismatchX = find(numelX-szV, 1);
    error(message('MATLAB:griddedInterpolant:CompVecValueMismatchErrId',...
        mismatchX, szV(mismatchX)));
end
gridVec = X;  % Grid vectors in a cell array.
clx = cellfun(@(x) string(underlyingType(x)), gridVec);
if any(clx ~= clx{1})
    error(message('MATLAB:griddedInterpolant:GridOfMixedDataTypesErrId'))
end
iscolX = cellfun(@iscolumn, gridVec);
gridVec(~iscolX) = cellfun(@transpose, gridVec(~iscolX), 'UniformOutput', false);
isCodist = cellfun(@iscodistributed, gridVec);
if iscodistributed(V)
    codV = getCodistributor(V);
    distDimV = codV.Dimension;
    newCodX = codistributor1d(1, codV.Partition, [numel(gridVec{distDimV}), 1]);
    if isCodist(distDimV)
        codX = getCodistributor(gridVec{distDimV});
        if ~all(codX.Partition == codV.Partition)
            gridVec{distDimV} = redistribute(gridVec{distDimV}, newCodX);
        end
    else
        gridVec{distDimV} = codistributed(gridVec{distDimV}, newCodX);
    end
    % Once one grid vector X is codistributed, the others should not be.
    dims = 1:nInterp;
    dims(distDimV) = [];
    gridVec(dims) = distributedutil.CodistParser.gatherElements(gridVec(dims));
elseif sum(isCodist) > 1  % More than one vector shouldn't be codistributed.
    szGridVec = cellfun(@numel, gridVec);
    [~, maxIsCodist] = max(szGridVec(isCodist)); % Find which of codistributed is largest.
    idxIsCodist = find(isCodist); % Find indices of codistributed gridvectors.
    toGather = isCodist;  % Potentially gather all that are codistributed.
    toGather(idxIsCodist(maxIsCodist)) = false;  % leave largest one codistributed.
    [gridVec{toGather}] = gather(gridVec{toGather});
end
end

%-------------------------------------------------------------------------%
function [localXq, codXq, szVout, distDimXq] = iValidateDistQuery(Xq)
% Find out which Xq's are distributed and gather all but one.
% If the output is grid vectors, Xq is a column, Yq is a row.

if any(cellfun(@isinteger, Xq))
    error(message('MATLAB:norm:inputType'))
elseif ~all(cellfun(@isfloat, Xq))
    error(message('MATLAB:griddedInterpolant:NonFloatInpPtsErrId'))
end
isXqDist = cellfun(@iscodistributed, Xq);
assert(any(isXqDist), 'iValidateDistQuery called with undistributed query')

% Determine size of output and get local part of query.
isXqScalar = all(cellfun(@isscalar, Xq));
isXqRow = cellfun(@isrow, Xq);
isXqCol = cellfun(@iscolumn, Xq);
isXqVec = all(isXqRow | isXqCol);
isMgVec = ~isXqScalar && isXqVec && (any(isXqRow) && any(isXqCol));
if isMgVec  % Xq is a set of meshgrid vectors.
    % make all Xq vectors column vectors
    Xq(isXqRow) = cellfun(@transpose, Xq(isXqRow), 'UniformOutput', false);
    firstDistrib = find(isXqDist, 1);
    isXqDist(firstDistrib) = false;
    distDimXq = firstDistrib;
    if any(isXqDist)  % Other vector is codistributed.
        [Xq{isXqDist}] = gather(Xq{isXqDist});
    end
    szVout = cellfun(@numel, Xq);
    codistributed.pVerifyUsing1d('interpn', Xq{distDimXq});
    codXq = getCodistributor(Xq{distDimXq});
    if sum(codXq.Partition) == 1   % Xq vector is distributed in wrong dimension,
        % so all of its values are on one worker.
        Xq{distDimXq} = redistribute(Xq{distDimXq}, codistributor1d(1));
        codXq = getCodistributor(Xq{distDimXq});
    end
    localXq = Xq;
    localXq{distDimXq} = getLocalPart(Xq{distDimXq});
    localXq{2} = localXq{2}.';  % Make Yq a row vector.
else  % Xq is an array of points to be evaluated
    szXq = cell2mat(cellfun(@size, Xq, 'UniformOutput', false)');
    if numel(Xq) > 1 && any(diff(szXq), 'all')
        error(message('MATLAB:griddedInterpolant:InputMixSizeErrId'))
    end
    if nnz(isXqDist)==1
        XqDist = find(isXqDist);
        codXq = getCodistributor(Xq{XqDist});
        XqNotDistDim = 1:numel(Xq);
        XqNotDistDim(XqDist) = [];
        Xq(XqNotDistDim) = cellfun(@(x) codistributed(x, codXq), ...
            Xq(XqNotDistDim), 'UniformOutput', false);
    else
        firstXqDist = find(isXqDist, 1);
        codXq = getCodistributor(Xq{firstXqDist});
        isXqDist(firstXqDist) = false;
        for i = find(isXqDist)
            codYq = getCodistributor(Xq{i});
            if (codXq.Dimension ~= codYq.Dimension) ||...
                  ~all(codXq.Partition == codYq.Partition)
                Xq{i} = redistribute(Xq{i}, codXq);
            end
        end
        for i = find(~isXqDist)
            Xq{i} = codistributed(Xq{i}, codXq);
        end
    end
    szVout = size(Xq{1});
    localXq = cellfun(@getLocalPart, Xq, 'UniformOutput', false);
    distDimXq = codXq.Dimension;
end
clx = cellfun(@(x) string(underlyingType(x)), Xq);
if any(clx ~= clx{1})
    error(message('MATLAB:griddedInterpolant:QueryOfMixedDataTypesErrId'))
end
end

%-------------------------------------------------------------------------%
function F = iMakeGriddedInterpolant(X, V, method, extrap, noX, szX, dimInterp)
%IMAKEGRIDDEDINTERPOLANT wrapper for interpLocalGriddedInterpolant'

[prevWarn, prevWarnID] = lastwarn;
lastwarn('nowarn');
if isscalar(X)
    szX = szX(1);
end
F = distributedutil.syncOnError(@interpLocalGriddedInterpolant, X, V, method,...
    extrap, noX, szX, dimInterp);

[warn, warnID] = lastwarn;
warns = spmdCat([{warn}, {warnID}], 1);  % get warnings from all workers.
if any(warns(:, 1) ~= "nowarn")
    firstWarn = find(warns(:, 1) ~= "nowarn", 1);
    notWarned = warns(:, 1) == "nowarn";
    notWarned(firstWarn) = false;  % already warned.
    if notWarned(spmdIndex)
        warning(warns{firstWarn, 2}, warns{firstWarn, 1});
    end
else
    lastwarn(prevWarn, prevWarnID);  % If no warnings this time, restore.
end
end

%-------------------------------------------------------------------------%
function localVout = iEvalGriddedInterpolant(F, Xq, extrapval)
% Use a gridded interpolant to calculate interpolated values on a worker.
% localVout has size determined by Xq.

isXqEmpty = cellfun(@isempty, Xq);
if isempty(F) || any(isXqEmpty)
    localVout = zeros(zeros([1 numel(Xq)]));
    return
end
isMgGrid = false;
isNdGrid = false;
isXqVector = cellfun(@isvector, Xq);
isXqScalar = cellfun(@isscalar, Xq);
iscompact = compactGridFormat(Xq{:});
if ~iscompact && ~any(isXqVector) && ~any(isXqEmpty)
    [isMgGrid, isNdGrid] = checkGrid(Xq, Xq{1});
end
if iscompact || (all(isXqScalar) && F.Method == "spline")
    if isa(F, 'griddedInterpolant')
        localVout = F(Xq);
    else  % V is a vector and method == "spline"
        localVout = ppval(F, Xq{:}).';
    end
elseif isNdGrid || isMgGrid
    [Xq{:}, ~] = getGridVectors(isMgGrid, isNdGrid, Xq{:}, Xq{1}); % getGridVectors returns ndgrid vectors.
    if isMgGrid  % This is equivalent to transposing 
        Xq = cellfun(@flipud, Xq, 'UniformOutput', false);
        [Xq(1), Xq(2)] = deal(Xq(2), Xq(1));
    end
    Xq{2} = Xq{2}.';   %  getGridVectors returns column vectors
    iscompact = true;
    localVout = F(Xq);
else
    localVout = F(Xq{:});
end

% Impose the extrapolation value to the queries that lie outside the
% domain.
if ~isempty(extrapval)
    localVout = imposeExtrapval(Xq, F.GridVectors, localVout,...
        extrapval, iscompact);
end
end

%-------------------------------------------------------------------------%
function iCheckX(X)
% Throws an error if not all of X are floats with finite values.

if ~all(cellfun(@(x) all(isfloat(x), 'all'), X))
    error(message('MATLAB:griddedInterpolant:NonFloatInpPtsErrId'))
end
if ~all(cellfun(@(x) distributedutil.CodistParser.gatherIfCodistributed...
        (all(isfinite(x), 'all')), X))
    error(message('MATLAB:griddedInterpolant:BadGridErrId'))
end
cellfun(@(x) codistributed.pVerifyUsing1d('interpn', x), X);
end

%-------------------------------------------------------------------------%
function tf = iIsDimensionOne(numNumericArgs, args)
% Returns true if args correspond to a 1d problem.

% Corresponds to interpn(x, v, xq, _)
tf = numNumericArgs == 3;
if numNumericArgs == 1
    % Corresponds to interpn(v, _)
    tf = isvector(args{1});
elseif numNumericArgs == 2
    % Corresponds to interpn(v, k, _)
    tf = isvector(args{2});
end
end
