function Vout = interp2(varargin)
%INTERP2 2-D interpolation (table lookup) for codistributed
%   Vq = interp2(V,Xq,Yq)
%   Vq = interp2(X,Y,V,Xq,Yq)
%   Vq = interp2(V,Xq,Yq,METHOD)
%   Vq = interp2(X,Y,V,Xq,Yq,METHOD)
%   Vq = interp2(V,Xq,Yq,METHOD,EXTRAPVAL)
%   Vq = interp2(X,Y,V,Xq,Yq,METHOD,EXTRAPVAL)
%   Vq = interp2(V,NTIMES)
%   Vq = interp2(V) % same as interp2(V, 1)
%   
%   The value argument V must be a double or single 2-D array, real or complex. V cannot be a vector.
%   
%   The grid arguments X and Y must have the following characteristics:
%     * same type, double or single.
%     * finite vectors or 2-D arrays with increasing and non-repeating elements in corresponding dimension. 
%     * When X and Y are non-vector 2-D arrays, they must align with cartesian axes, as if they were 
%       produced by meshgrid.
%     * dimensions consistent with those of value V. 
%   
%   The query arguments Xq and Yq must be vectors or arrays of the same type (double or single). 
%   If Xq and Yq are arrays, they must be the same size. If they are vectors with different lengths, 
%   they must have different orientations.
%   
%   NTIMES currently only supports non-negative real scalars.
%   
%   The extrapolation for the out-of-boundary input is not supported.
%   
%   Example:
%   spmd
%       [X,Y,V] = peaks(10);
%       xq = codistributed.colon(-3, 0.01, 3);
%       yq = codistributed.colon(-3, 0.01, 3);
%       Vq = interp2(X, Y, V, xq', yq);
%       [Xq,Yq] = meshgrid(xq, yq);  
%       mesh(Xq,Yq,Vq);
%   end
%   
%   creates an array Vq from bilinearly interpolating function values in V. 
%   The size of Vq equals to that of Xq or Yq. Or the number of rows in Vq equals to the length of yq. 
%   The number of columns in Vq equals to the length of xq.  
%   
%   See also INTERP2, INTERP3, INTERPN, MESHGRID, SCATTEREDINTERPOLANT, CODISTRIBUTED.


%   Copyright 2013-2023 The MathWorks, Inc.

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

% Parse the method and extrap val
[narg, method, extrapval] = methodAndExtrapval(varargin{:});
if narg ~= 1 && narg ~= 2 && narg ~= 3 && narg ~= 5
    error(message('MATLAB:interp2:nargin'));
end
stripNaNsForCubics = method == "spline" || method == "makima";
if stripNaNsForCubics
    extrap = method;
else
    extrap = 'none';
end
if method == "pchip"
    warning(message("MATLAB:griddedInterpolant:Pchip1DonlyWarnId"))
    method = 'linear';
end

noX = false;
if narg == 1 || narg == 2
    %  interp2(V,NTIMES)
    V = varargin{1};
    if isvector(V)
        error(message('MATLAB:griddedInterpolant:DegenerateGridErrId'))
    end
    [nrows,ncols] = size(V);
    if narg == 1
        ntimes = 1;
    else
        ntimes = varargin{2}(1);
        if ~isscalar(ntimes) || ~isreal(ntimes)
            error(message('MATLAB:interp2:NtimesInvalid'));
        end
        if iscodistributed(ntimes)
            gather(ntimes);
        end
        ntimes = floor(ntimes);
    end
    if isinf(ntimes)
        Xq = [];
        Yq = [];
    elseif isnan(ntimes)
        Xq = NaN;
        Yq = NaN;
    else
        Xq = linspace(1, ncols, (ncols-1)*2^ntimes+1)';
        Yq = linspace(1, nrows, (nrows-1)*2^ntimes+1);
    end
    noX = true;
elseif narg == 3
    %  interp2(V,Xq,Yq)
    [V, Xq, Yq] = deal(varargin{1:3});
    noX = true;
elseif narg == 5
    %  interp2(X,Y,V,Xq,Yq)
    [X, Y, V, Xq, Yq] = deal(varargin{1:5});
    iCheckXY(X, Y)
    if ~isvector(X) || ~isvector(Y)
        % X and Y are matrices, make sure they define a meshgrid and
        % extract the grid vectors.
        [isMeshGrid, isNdGrid, exception] = checkGrid(X, Y, V);
        if strcmpi(exception, 'SizeMismatch')
            if ~all(size(X) == size(Y))  % Grid vectors don't match each other.
                error(message('MATLAB:griddedInterpolant:BadGridErrId'))
            else  % One or both grid vectors doesn't match V.
                error(message('MATLAB:griddedInterpolant:MixedGridCoordSizeErrId'))
            end
        elseif strcmpi(exception, 'DimensionMismatch')
            if ~ismatrix(X) || ~ismatrix(Y)
                error(message('MATLAB:transpose:NDArray'))
            else
                error(message('MATLAB:griddedInterpolant:BadGridErrId'))
            end
        elseif isNdGrid
            error(message('MATLAB:interp2:InvalidMeshgrid'));
        elseif ~isMeshGrid
            error(message('MATLAB:griddedInterpolant:BadGridErrId'))
        end
        [X, Y, ~] = getGridVectors(isMeshGrid, false, X, Y, V);
    else
        [X, Y] = iValidateGridVectors(X, Y, V);
    end
end

% Every call passes through this section.
if noX
    [X, Y] = ndgridVectors(V);
end

if ~any(ismember(method, {'linear', 'nearest', 'cubic', 'makima', 'spline'}))
    error(message('MATLAB:griddedInterpolant:BadInterpTypeErrId'))
end

[V, origvtype] = iConvertV(V, method, Xq, Yq);

[X, Y, V] = checkMonotonicGridVectors(X, Y, V);
if stripNaNsForCubics
    [X, Y, V] = iStripNaNWrapper(X, Y, V);
end

isDistGrid = iscodistributed(X) || iscodistributed(Y) || iscodistributed(V);
isDistQuery = iscodistributed(Xq) || iscodistributed(Yq);

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

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

elseif isDistGrid
    % Only grid is distributed
    Vout = iInterpDistGridReplQuery(F, {X, Y}, V, {Xq, Yq}, 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({Y, X}, V, {Xq, Yq}, method,...
        extrapval, noX);  % grid vectors are reversed because this calls
                          % in-memory interp2, which wants a meshgrid.

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

if ~strcmp(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);
Xq = fliplr(Xq);   % Change from meshgrid to ndgrid.
if iscodistributed(V)
    codistributed.pVerifyUsing1d('interp2', 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);
        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 = spmdBroadcast(find(codX.Partition), 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 interp2 needs three points for full accuracy.
    padSz = 3;
elseif method == "cubic" || method == "spline"
    % 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{1}, Xq{2}, Xq{1});
    if exception == "SizeMismatch"
        error(message('MATLAB:griddedInterpolant:InputMixSizeErrId'))
    end
    if isNdGrid
        [Xq{1}, Xq{2}, ~] = getGridVectors(false, isNdGrid, Xq{1}, Xq{2}, Xq{1}); % getGridVectors returns ndgrid vectors.
        Xq{2} = Xq{2}.';   %  getGridVectors returns column vectors
        isMgVec = true;
    end
elseif isMgVec && isrow(Xq{1})  % if inputs are row and column, transpose both.
    Xq{1} = Xq{1}.';
    Xq{2} = Xq{2}.';
end
if isMgVec
    localSzVout = [numel(Xq{1}) numel(Xq{2})];
    distDimVout = distDimX;
else
    if ~all(size(Xq{1}) == size(Xq{2}))
        error(message('MATLAB:griddedInterpolant:InputMixSizeErrId'))
    end
    szVout = size(Xq{1});
    % make them column vectors for evaluation - reshape at end.
    Xq{1} = Xq{1}(:);
    Xq{2} = Xq{2}(:);
    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 any(isnan(Xq{1})) || any(isnan(Xq{2}))
        queriesOnThisWorker = isnan(Xq{1}) | isnan(Xq{2}) | queriesOnThisWorker;
    end
    localXq = cellfun(@(x) x(queriesOnThisWorker), Xq, 'UniformOutput', false);
end
if method == "spline"
    localXqInDistDim = localXq{distDimX};
    [uLocalXqInDistDim, ~, idxIn] = unique(localXqInDistDim);
    notDistDimX = mod(distDimX, 2) + 1;
    if distDimX == 2
        localXq{distDimX} = cast(localX.', 'like', uLocalXqInDistDim);
    else
        localXq{distDimX} = cast(localX, 'like', uLocalXqInDistDim);
    end
    if isMgVec
        localXq{notDistDimX} = Xq{notDistDimX};
        if ~compactGridFormat(localXq{:})
            localXq{notDistDimX} = localXq{notDistDimX}.';
        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});
        localXq{distDimX} = reshape(repmat(localXq{distDimX}, 1, nXq), nX*nXq, 1);
        localXq{notDistDimX} = reshape(repmat(Xq{notDistDimX}.', nX, 1), nX*nXq, 1);
        origExtrapval = extrapval;
        extrapval = [];
    end
end
localVout = iEvalGriddedInterpolant(F, localXq, extrapval);
if method == "spline" && ~isempty(localVout)
    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).';
            splineVout = splineVout(idxIn, :);
            localVout = ipermute(double(splineVout), Fsp.dims);
            localXq{distDimX} = localXqInDistDim;       % Restore localXq
        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 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

% 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, Yq)
%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, Xq, and Yq.

origvtype = underlyingType(V);
if ~isfloat(V) && (strcmpi(method, 'nearest') || (isscalar(Xq) && isscalar(Yq)))
    V = double(V);
end
if ~ismatrix(V)
    error(message('MATLAB:transpose:NDArray'));
elseif ~isfloat(V)
    error(message('MATLAB:griddedInterpolant:NonFloatInpPtsErrId'))
end
if ~isfloat(Xq) || ~isfloat(Yq)
    error(message('MATLAB:griddedInterpolant:NonFloatInpPtsErrId'))
elseif ~isreal(Xq) || ~isreal(Yq)
    error(message('MATLAB:griddedInterpolant:ComplexDataPointErrId'))
end
codistributed.pVerifyUsing1d('interp2', Xq);
codistributed.pVerifyUsing1d('interp2', Yq);
end

%-------------------------------------------------------------------------%
function [X, Y, V] = iStripNaNWrapper(X, Y, V)
%ISTRIPNANWRAPPER removes rows or columns of V containing NaNs.

[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
if ncolnan > nrownan
    [Y, X, V] = stripNaNsForSpline(Y, X, V.');
    V = V.';
else
    [X, Y, V] = stripNaNsForSpline(X, Y, V);
end
warning(message('MATLAB:interp2:NaNstrip'))
if isempty(V) || isvector(V)
    error(message('MATLAB:interp2:NotEnoughPointsNanStrip'))
end
if isempty(V) || isvector(V)
    error(message('MATLAB:interp2:NotEnoughPointsNanStrip'))
end
end

%-------------------------------------------------------------------------%
function [X, Y] = iValidateGridVectors(X, Y, V)
% Convert grid vectors from meshgrid to ndgrid format and make sure that
% only one 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', 'INTERP2', func2str(V)))
elseif isa(V, 'struct')
    error(message('MATLAB:invalidConversion', 'struct', underlyingType(V)))
elseif iscell(V)
    error(message('MATLAB:invalidConversion', 'cell', underlyingType(V)))
elseif isscalar(V) && ~iscell(V) || isvector(V)
    error(message('MATLAB:griddedInterpolant:DegenerateGridErrId'));
elseif numel(X) ~= szV(2)
    error(message('MATLAB:griddedInterpolant:CompVecValueMismatchErrId', 1, szV(2)));
elseif numel(Y) ~= szV(1)
    error(message('MATLAB:griddedInterpolant:CompVecValueMismatchErrId', 2, szV(1)));
end
if ~strcmpi(underlyingType(X), underlyingType(Y))
    error(message('MATLAB:griddedInterpolant:GridOfMixedDataTypesErrId'))
end
[Y, X] = deal(X, Y);  % Swap to make ndgrid
if isrow(X)           % Transpose rows to make columns
    X = X.';
end
if isrow(Y)
    Y = Y.';
end
isCodX = iscodistributed(X);
isCodY = iscodistributed(Y);
if iscodistributed(V)
    codV = getCodistributor(V);
    distDimV = codV.Dimension;
    if distDimV == 1
        newCodX = codistributor1d(1, codV.Partition, [numel(X), 1]);
        if isCodX
            codX = getCodistributor(X);
            if ~all(codX.Partition == codV.Partition)
                X = redistribute(X, newCodX);
            end
        else
            X = codistributed(X, newCodX);
        end
        % If grid vector X is codistributed, Y should not be.
        Y = distributedutil.CodistParser.gatherIfCodistributed(Y);
    else  % distDimV = 2
        newCodY = codistributor1d(1, codV.Partition, [numel(Y), 1]);
        if isCodY
            codY = getCodistributor(Y);
            if ~all(codY.Partition == codV.Partition)
                Y = redistribute(Y, newCodY);
            end
        else
            Y = codistributed(Y, newCodY);
        end
        % If grid vector Y is codistributed, X should not be.
        X = distributedutil.CodistParser.gatherIfCodistributed(X);
    end
elseif isCodX && isCodY  % Both grid vectors shouldn't be codistributed.
    Y = gather(Y);
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 isinteger(Xq{1}) || isinteger(Xq{2})
    error(message('MATLAB:norm:inputType'))
elseif ~isfloat(Xq{1}) || ~isfloat(Xq{2})
    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.
if ~isscalar(Xq{1}) && ((isrow(Xq{1}) && iscolumn(Xq{2})) || ...
    iscolumn(Xq{1}) && isrow(Xq{2}))  % Xq is a pair of meshgrid vectors.
    if isrow(Xq{1}) % make both column vectors
        Xq{1} = Xq{1}.';
    else
        Xq{2} = Xq{2}.';
    end
    % Because input is interpreted as meshgrid vectors, swap elements of Xq.
    Xq = fliplr(Xq);
    isXqDist = fliplr(isXqDist);
    firstDistrib = find(isXqDist, 1);
    isXqDist(firstDistrib) = false;
    distDimXq = firstDistrib;
    if any(isXqDist)  % Other vector is codistributed.
        Xq{isXqDist} = gather(Xq{isXqDist});
    end
    szVout = [numel(Xq{1}) numel(Xq{2})];
    codistributed.pVerifyUsing1d('interp2', 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.
    localXq = fliplr(localXq);
else
    if ~all(size(Xq{1}) == size(Xq{2}))
        error(message('MATLAB:griddedInterpolant:InputMixSizeErrId'))
    end
    if nnz(isXqDist)==1
        XqDist = find(isXqDist);
        codXq = getCodistributor(Xq{XqDist});
        XqNotDistDim = mod(XqDist, 2) + 1;
        Xq{XqNotDistDim} = codistributed(Xq{XqNotDistDim}, codXq);
    else
        codXq = getCodistributor(Xq{1});
        codYq = getCodistributor(Xq{2});
        if (codXq.Dimension ~= codYq.Dimension) ||...
                ~all(codXq.Partition == codYq.Partition)
            Xq{2} = redistribute(Xq{2}, codXq);
        end
    end
    szVout = size(Xq{1});
    localXq = {getLocalPart(Xq{1}), getLocalPart(Xq{2})};
    distDimXq = codXq.Dimension;
end
if ~isequal(underlyingType(Xq{1}), underlyingType(Xq{2}))
    error(message('MATLAB:griddedInterpolant:GridOfMixedDataTypesErrId'))
end
end

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

gime = 'noerror';
[prevWarn, prevWarnID] = lastwarn;
lastwarn('nowarn');
try
    F = interpLocalGriddedInterpolant(X, V, method, extrap, noX, szX, dimInterp);
catch gime
end
gimes = spmdCat({gime}, 1);  % get errors from all workers.
if any(~strcmpi(gimes, 'noerror'))
    firstError = find(~strcmpi(gimes, 'noerror'), 1);
    rethrow(gimes{firstError});
end
[warn, warnID] = lastwarn;
warns = spmdCat([{warn}, {warnID}], 1);  % get warnings from all workers.
if any(~strcmpi(warns(:, 1), 'nowarn'))
    firstWarn = find(~strcmpi(warns(:, 1), 'nowarn'), 1);
    notWarned = strcmpi(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.

if isempty(F) || isempty(Xq{1}) || isempty(Xq{2})
    localVout = [];
    return
end
isMgGrid = false;
isNdGrid = false;
iscompact = compactGridFormat(Xq{:});
if (~all(isfinite(Xq{1}), 'all') || ~all(isfinite(Xq{2}), 'all')) && ~isvector(Xq{1}) ...
        && ~isvector(Xq{2})
    Xq = fliplr(Xq);
elseif ~iscompact && ~isvector(Xq{1}) && ~isvector(Xq{2}) && ~isempty(Xq{1})...
        && ~isempty(Xq{2})
    [isMgGrid, isNdGrid] = checkGrid(Xq{1}, Xq{2}, Xq{1});
end
if iscompact || (isscalar(Xq{1}) && isscalar(Xq{2}) && F.Method == "spline")
    localVout = F(Xq);
elseif isNdGrid || isMgGrid
    [Xq{1}, Xq{2}, ~] = getGridVectors(isMgGrid, isNdGrid, Xq{1}, Xq{2}, Xq{1}); % getGridVectors returns ndgrid vectors.
    if isMgGrid  % This is equivalent to transposing 
        temp = flipud(Xq{2});
        Xq{2} = flipud(Xq{1});
        Xq{1}= temp;
    end
    Xq{2} = Xq{2}.';   %  getGridVectors returns column vectors
    iscompact = true;
    localVout = F(Xq);
else
    localVout = F(Xq{1}, Xq{2});
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 iCheckXY(X, Y)
% Throws an error if X and Y are not both floats with finite values.

if iscell(X) || iscell(Y)
    error(message('parallel:array:UnsupportedFunction', 'INTERP2', 'cell'))
end
if ~isfloat(X) || ~isfloat(Y)
    error(message('MATLAB:griddedInterpolant:NonFloatInpPtsErrId'))
end
if any(~isfinite(X), 'all') || any(~isfinite(Y), 'all')
    error(message('MATLAB:griddedInterpolant:BadGridErrId'))
end
codistributed.pVerifyUsing1d('interp2', X);
codistributed.pVerifyUsing1d('interp2', Y);
end
