function Vout = interp3(varargin)
%INTERP3 3-D interpolation (table lookup) for codistributed
%   Vq = interp3(V,Xq,Yq,Zq)
%   Vq = interp3(X,Y,Z,V,Xq,Yq,Zq)
%   Vq = interp3(V,Xq,Yq,Zq,METHOD)
%   Vq = interp3(X,Y,Z,V,Xq,Yq,Zq,METHOD)
%   Vq = interp3(V,Xq,Yq,Zq,METHOD,EXTRAPVAL)
%   Vq = interp3(X,Y,Z,V,Xq,Yq,Zq,METHOD,EXTRAPVAL)
%   Vq = interp3(V,NTIMES)
%   Vq = interp3(V) % same as interp3(V, 1)
%   
%   The value argument V must be a double or single 3-D array, real or complex. 
%   
%   The grid arguments X, Y and Z must have the following characteristics:
%     * Same type, double or single.
%     * Finite vectors or 3-D arrays with increasing and non-repeating elements in corresponding dimension. 
%     * When X, Y, and Z are 3-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, Yq and Zq must be vectors or arrays of the same type (double or single). 
%   If Xq, Yq and Zq are arrays, they must be the same size. If they are vectors with different lengths, 
%   one of them must have different orientation.
%   
%   Example:
%   spmd
%       [X,Y,Z,V] = flow(10);
%       xq = codistributed.colon(0.1, 0.1, 10);
%       yq = codistributed.colon(-3, 0.1, 3);
%       zq = codistributed.colon(-3, 0.1, 3);
%       Vq = interp3(X, Y, Z, V, xq', yq, zq); % Vq is 61-by-100-by-61
%       [Xq,Yq,Zq] = meshgrid(xq, yq, zq);  
%       slice(Xq,Yq,Zq,Vq,[6 9.5],2,[-2 .2]); % shading flat
%   end
%   
%   creates an array Vq from linearly interpolating function values in V. 
%   The size of Vq equals that of Xq or Yq or Zq. The number of rows in Vq equals the length of yq. 
%   The number of columns in Vq equals the length of xq. The number of elements in the third dimension
%   equals the length of zq.  
%   
%   See also INTERP1, INTERP2, INTERPN, MESHGRID, SCATTEREDINTERPOLANT, CODISTRIBUTED.


%   Copyright 2013-2023 The MathWorks, Inc.

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

% Parse the method and extrap val
[narg, method, extrapval] = methodAndExtrapval(varargin{:});
if narg ~= 1 && narg ~= 2 && narg ~= 4 && narg ~= 7
    error(message('MATLAB:interp3: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
    %  interp3(V,NTIMES)
    V = varargin{1};
    [nrows, ncols, npages] = size(V);
    if narg == 1
        ntimes = 1;
    else
        ntimes = varargin{2}(1);
        if ~isscalar(ntimes) || ~isreal(ntimes)
            error(message('MATLAB:interp3:NtimesInvalid'))
        end
        ntimes = distributedutil.CodistParser.gatherIfCodistributed(ntimes);
        ntimes = floor(ntimes);
    end
    if isinf(ntimes)
        Xq = [];
        Yq = [];
        Zq = [];
    elseif isnan(ntimes)
        Xq = NaN;
        Yq = NaN;
        Zq = NaN;
    elseif ntimes < 0
        Xq = (1:1/(2^ntimes):ncols)';
        Yq = 1:1/(2^ntimes):nrows;
        Zq = 1:1/(2^ntimes):npages;
    else
        Xq = linspace(1, ncols, (ncols-1)*2^ntimes+1)';
        Yq = linspace(1, nrows, (nrows-1)*2^ntimes+1);
        Zq = linspace(1, npages, (npages-1)*2^ntimes+1);
    end
    szVout = numel(Xq)*numel(Yq)*numel(Zq);
    % 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([numel(Xq), numel(Yq), numel(Zq)]);
        if iscodistributed(V)
            codV = getCodistributor(V);
            distDim = codV.Dimension;
        end
        if (~iscodistributed(V) || distDim > 3)
            distDim = longestVec;
        end
        switch distDim
            case 1
                Xq = codistributed(Xq, codistributor1d(1));
            case 2
                Yq = codistributed(Yq, codistributor1d(2));
            case 3
                Zq = codistributed(Zq, codistributor1d(2));
        end
    end
    noX = true;
elseif narg == 4
    %  interp3(V,Xq,Yq,Zq)
    [V, Xq, Yq, Zq] = deal(varargin{1:4});
    noX = true;
elseif narg == 7
    %  interp3(X,Y,Z,V,Xq,Yq,Zq)
    [X, Y, Z, V, Xq, Yq, Zq] = deal(varargin{1:7});
    iCheckXYZ(X, Y, Z)
    if ~isvector(X) || ~isvector(Y) || ~isvector(Z)
        % X, Y, and Z are not vectors; make sure they define a meshgrid and
        % extract the grid vectors.
        [isMeshGrid, isNdGrid, exception] = checkGrid(X, Y, Z, V);
        if exception == "ArgMismatch"
            if ~all(size(X) == size(Y)) || ~all(size(X) == size(Z)) 
                % 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:BadGridErrId'))
        elseif exception == "TypeMismatch"
            error(message('MATLAB:griddedInterpolant:NumCoordsGridNdimsMismatchErrId'));   
        elseif isNdGrid
            error(message('MATLAB:interp3:InvalidMeshgrid'));
        elseif ~isMeshGrid
            error(message('MATLAB:griddedInterpolant:BadGridErrId'))
        end
        [X, Y, Z, ~] = getGridVectors(isMeshGrid, false, X, Y, Z, V);
    else
        [X, Y, Z] = iValidateGridVectors(X, Y, Z, V);
    end
end

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

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

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

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

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

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

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

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

else
    % No data input was distributed. Gather and call in-memory.
    [varargin{:}] = gather(varargin{:});
    Vout = interp3(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);
[Xq(1), Xq(2)] = deal(Xq(2), Xq(1));   % Change from meshgrid to ndgrid.
if iscodistributed(V)
    codistributed.pVerifyUsing1d('interp3', 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 = 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 interp3 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{3}, Xq{1});
    if exception == "SizeMismatch"
        error(message('MATLAB:griddedInterpolant:InputMixSizeErrId'))
    end
    if isNdGrid
        [Xq{1}, Xq{2}, Xq{3}, ~] = getGridVectors(false, isNdGrid, Xq{1}, Xq{2}, Xq{3}, Xq{1}); % getGridVectors returns ndgrid vectors.
        Xq{2} = Xq{2}.';   %  getGridVectors returns column vectors
        isMgVec = true;
    end
end
if isMgVec
    localSzVout = [numel(Xq{1}) numel(Xq{2}) numel(Xq{3})];
    distDimVout = distDimX;
else
    szXq = cell2mat(cellfun(@size, Xq, 'UniformOutput', false)');
    if any(diff(szXq), 'all')
        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}(:);
    Xq{3} = Xq{3}(:);
    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
        queriesOnThisWorker = isnan(Xq{1}) | isnan(Xq{2}) | isnan(Xq{3})| queriesOnThisWorker;
    end
    localXq = cellfun(@(x) x(queriesOnThisWorker), Xq, 'UniformOutput', false);
end
if method == "spline"
    localXqInDistDim = localXq{distDimX};
    [uLocalXqInDistDim, ~, idxIn] = unique(localXqInDistDim);
    notDistDimX = 1:3;
    notDistDimX(distDimX) = [];
    if distDimX == 2
        localXq{distDimX} = cast(localX.', 'like', localXqInDistDim);
    else
        localXq{distDimX} = cast(localX, 'like', localXqInDistDim);
    end
    if isMgVec
        localXq{notDistDimX(1)} = Xq{notDistDimX(1)};
        localXq{notDistDimX(2)} = Xq{notDistDimX(2)};
        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(1)} = reshape(repmat(Xq{notDistDimX(1)}.', nX, 1), nX*nXq, 1);
        localXq{notDistDimX(2)} = reshape(repmat(Xq{notDistDimX(2)}.', 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) || any(queriesOnThisWorker)
        if isMgVec % The query is a grid, so we calculate all combinations.
            if ~isrow(uLocalXqInDistDim)
                uLocalXqInDistDim = uLocalXqInDistDim.';
            end
            splineVout = permute(ppval(Fsp, uLocalXqInDistDim), [3 1 2]);
            splineVout = splineVout(idxIn, :, :);
            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 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 = zeros(0, 0, 0);
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, Zq)
%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, Yq and Zq.

origvtype = underlyingType(V);
if ~isfloat(V) && (method == "nearest" || (isscalar(Xq) && ...
        isscalar(Yq) && isscalar(Zq)))
    V = double(V);
end
if ndims(V) > 3
    error(message('MATLAB:permute:tooFewIndices'));
elseif ~isfloat(V)
    error(message('MATLAB:griddedInterpolant:NonFloatInpPtsErrId'))
end
if ~isfloat(Xq) || ~isfloat(Yq) || ~isfloat(Zq)
    error(message('MATLAB:griddedInterpolant:NonFloatInpPtsErrId'))
elseif ~isreal(Xq) || ~isreal(Yq) || ~isreal(Zq)
    error(message('MATLAB:griddedInterpolant:ComplexDataPointErrId'))
end
codistributed.pVerifyUsing1d('interp3', Xq);
codistributed.pVerifyUsing1d('interp3', Yq);
codistributed.pVerifyUsing1d('interp3', Zq);
end

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

[inan, jnan, knan] = find(isnan(V));
inan = distributedutil.CodistParser.gatherIfCodistributed(inan);
jnan = distributedutil.CodistParser.gatherIfCodistributed(jnan);
knan = distributedutil.CodistParser.gatherIfCodistributed(knan);
ncolnan = length(unique(jnan));
nrownan = length(unique(inan));
npagenan = length(unique(knan));
if ncolnan == 0 && nrownan == 0 && npagenan == 0
    return
end
if ncolnan > nrownan
    permutedV = permute(V, [2 1 3]);
    [Y, X, V] = stripNaNsForSpline(Y, X, Z, permutedV);
    V = ipermute(V, [2 1 3]);
else
    [X, Y, Z, V] = stripNaNsForSpline(X, Y, Z, V);
end
warning(message('MATLAB:interp3:NaNstrip'))
if isempty(V) || ismatrix(V)
    error(message('MATLAB:interp3:NotEnoughPointsNanStrip'))
end
end

%-------------------------------------------------------------------------%
function [X, Y, Z] = iValidateGridVectors(X, Y, Z, 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);
[Y, X] = deal(X, Y);  % Swap to make ndgrid
if isa(V, 'function_handle')
    error(message('parallel:array:UnsupportedFunction', 'INTERP3', 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) || ismatrix(V)
    error(message('MATLAB:griddedInterpolant:DegenerateGridErrId'));
elseif numel(X) ~= szV(1)
    error(message('MATLAB:griddedInterpolant:CompVecValueMismatchErrId', 1, szV(1)));
elseif numel(Y) ~= szV(2)
    error(message('MATLAB:griddedInterpolant:CompVecValueMismatchErrId', 2, szV(2)));
elseif numel(Z) ~= szV(3)
    error(message('MATLAB:griddedInterpolant:CompVecValueMismatchErrId', 3, szV(3)));
end
gridVec = {X, Y, Z};  % Grid vectors in a cell array.
clx = cellfun(@(x) string(underlyingType(x)), gridVec);
if any(clx ~= clx{1})
    error(message('MATLAB:griddedInterpolant:NumCoordsGridNdimsMismatchErrId'))
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:3;
    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
[X, Y, Z] = deal(gridVec{:});  % Retrieve grid vectors from cell array.
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);
    % Because input is interpreted as meshgrid vectors, swap first two elements of Xq.
    [Xq(2), Xq(1)] = deal(Xq(1), Xq(2));
    [isXqDist(2), isXqDist(1)] = deal(isXqDist(1), isXqDist(2));
    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}) numel(Xq{3})];
    codistributed.pVerifyUsing1d('interp3', 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(2), localXq(1)] = deal(localXq(1), localXq(2));  % put it back to meshgrid
else  % Xq is an array of points to be evaluated
    szXq = cell2mat(cellfun(@size, Xq, 'UniformOutput', false)');
    if any(diff(szXq), 'all')
        error(message('MATLAB:griddedInterpolant:InputMixSizeErrId'))
    end
    if nnz(isXqDist)==1
        XqDist = find(isXqDist);
        codXq = getCodistributor(Xq{XqDist});
        XqNotDistDim = 1:3;
        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 = 1:firstXqDist-1
            Xq{i} = codistributed(Xq{i}, codXq);
        end
    end
    szVout = size(Xq{1});
    localXq = {getLocalPart(Xq{1}), getLocalPart(Xq{2}), getLocalPart(Xq{3})};
    distDimXq = codXq.Dimension;
end
clx = cellfun(@(x) string(underlyingType(x)), Xq);
if any(clx ~= clx{1})
    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(gimes ~= "noerror")
    firstError = find(gimes ~= "noerror", 1);
    rethrow(gimes{firstError});
end
[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.

if isempty(F) || isempty(Xq{1}) || isempty(Xq{2}) || isempty(Xq{3})
    localVout = zeros(0, 0, 0);
    return
end
isMgGrid = false;
isNdGrid = false;
isXqFinite = cellfun(@(x) all(isfinite(x), 'all'), Xq);
isXqVector = cellfun(@isvector, Xq);
isXqEmpty = cellfun(@isempty, Xq);
isXqScalar = cellfun(@isscalar, Xq);
iscompact = compactGridFormat(Xq{:});
if (~all(isXqFinite)) && ~any(isXqVector)
    [Xq(2), Xq(1)] = deal(Xq(1), Xq(2));
elseif ~iscompact && ~any(isXqVector) && ~any(isXqEmpty)
    [isMgGrid, isNdGrid] = checkGrid(Xq{1}, Xq{2}, Xq{3}, Xq{1});
end
if iscompact || (all(isXqScalar) && F.Method == "spline")
    localVout = F(Xq);
elseif isNdGrid || isMgGrid
    [Xq{1}, Xq{2}, Xq{3}, ~] = getGridVectors(isMgGrid, isNdGrid, Xq{1}, Xq{2}, Xq{3}, 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;
        Xq{3} = flipud(Xq{3});
    end
    Xq{2} = Xq{2}.';   %  getGridVectors returns column vectors
    iscompact = true;
    localVout = F(Xq);
else
    localVout = F(Xq{1}, Xq{2}, Xq{3});
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 iCheckXYZ(X, Y, Z)
% Throws an error if X, Y, and Z are not all floats with finite values.

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