function varargout = parallel_function(range, F, consume, supply, ...
                                       reduce, identity, concat, empty, ...
                                       M, divide, next_divide)
%PARALLEL_FUNCTION Basis of parfor
%   PARALLEL_FUNCTION is the basis of parfor, but is not officially supported,
%   either with this name or API.
%
%   See also parfor, spmd, parfeval, parpool, parcluster.

%   Copyright 1984-2019 The MathWorks, Inc.

%   The following is a help-style description of parallel_function.
%
%   PARALLEL_FUNCTION runs a function in parallel and consumes/reduces the results.
%   PARALLEL_FUNCTION(RANGE, F) takes RANGE and a function handle F.  RANGE must
%   be a two-element row vector of integers which we'll call BASE and LIMIT.
%   The range specifies a semi-open interval (BASE-1,LIMIT].  If BASE >= LIMIT, nothing
%   happens because the interval is empty.  Otherwise, let N denote LIMIT-(BASE-1),
%   i.e., the interval includes N integers.  PARALLEL_FUNCTION chooses k-1
%   intermediate points BASE < N_1 < ... < N_{k-1} < LIMIT, effectively dividing
%   the interval into k segments. Control over choosing k and the N_j will be
%   discussed below. We will use the convention that N_0 = BASE-1 and
%   N_k = LIMIT.  Independently of k and the intermediate N_j, (N_0,N_k] is the
%   union of the subintervals (N_0,N_1],...,(N_{k-1},N_k].
%
%   PARALLEL_FUNCTION's essence:
%
%     * In parallel, evaluate F(N_{j-1},N_j) in "worker" processes, j = 1,...,k.
%
%   The process on which PARALLEL_FUNCTION is called is the "client".
%
%   Because the choice of intermediate points is non-deterministic,
%   PARALLEL_FUNCTION will be non-deterministic unless F obeys the rule that the
%   following statement sequences are equivalent:
%
%     * F(N_{j-1}, N_j) % process a subinterval
%       F(N_j, N_{j+1}) % process the next subinterval
%
%     * F(N_{j-1}, N_{j+1}) % process the combined subintervals
%
%   Because the order of processing is non-deterministic, PARALLEL_FUNCTION will
%   be non-deterministic unless F obeys the rule that for any j1 and j2, these
%   are equivalent:
%
%     * F(N_{j1-1}, N_j1); F(N_{j2-1}, N_j2); % process j1 first
%
%     * F(N_{j2-1}, N_j2); F(N_{j1-1}, N_j1); % process j2 first
%
%   Actually, the requirement is stronger than this: the two statements here
%   must not interfere with one another's execution.  For example, if they are
%   both competing to write the same file, there will be trouble.  We will call
%   this the "non-interference" rule.
%
%   Not only N_{j-1} and N_j are sent from the client to the workers, but so is
%   the function handle F.  Handles to nested functions are particularly useful
%   as values of F, because the values in their workspace are sent along with F.
%   This is the standard way to communicate loop-invariant values to workers.
%
%   Note: PARALLEL_FUNCTION is the basis for PARFOR execution, and if PARFOR is
%   sufficient for your purposes, it is a more convenient interface to this
%   functionality.  PARFOR uses the nested function technique to transmit its
%   unsliced input variables.
%
%
%   PARALLEL_FUNCTION(RANGE, F, CONSUME) provides a way for workers to send
%   values back to the client.  If CONSUME is [], the behavior is as above, i.e,
%   F is evaluated with no expected results.  Otherwise, CONSUME must be a
%   function handle, and F must return a result.  Each worker does this:
%
%     * Evaluates O = F(N_{j-1}, N_j).
%
%     * Sends O to the client.
%
%   As each O is received, the client calls CONSUME(N_{j-1},N_j, O).
%
%   Because the choice of intermediate points is non-deterministic,
%   PARALLEL_FUNCTION will be non-deterministic unless F and consume obey the
%   rule that the following statement sequences are equivalent:
%
%     * O1 = F(N_{j-1}, N_j);      % process a subinterval
%       O2 = F(N_j, N_{j+1});      % process the next subinterval
%       CONSUME(N_{j-1}, N_j, O1)  % consume the first output
%       CONSUME(N_j, N_{j+1}, O2)  % consume the second output
%
%     * O = F(N_{j-1}, N_{j+1});     % process the combined subintervals
%       CONSUME(N_{j-1}, N_{j+1}, O) % consume the entire output
%
%   Because the order of processing is non-deterministic, PARALLEL_FUNCTION will
%   be non-deterministic unless F obeys the non-interference rule, extended for
%   outputs: for any distinct j1 and j2, the following may be executed
%   concurrently:
%
%     * O1 = F(N_{j1-1}, N_j1);
%
%     * O2 = F(N_{j2-1}, N_j2)
%
%   Calls on CONSUME occur in the client, and thus are serially executed.  Since
%   the order in which this happens is not guaranteed, PARALLEL_FUNCTION will be
%   non-deterministic unless CONSUME obeys the rule that for any j1 and j2, where
%   O1 and O2 arise as above, these are equivalent:
%
%     * CONSUME(N_{j1-1},N_j1, O1); CONSUME(N_{j2-1},N_j2, O2)
%
%     * CONSUME(N_{j2-1},N_j2, O2); CONSUME(N_{j1-1},N_j1, O1)
%
%   We say that CONSUME is "order-insensitive".
%
%   The final constraint comes from the observation that CONSUME runs on the
%   client and F on a worker, so there is an other form of the non-interference
%   requirement: for any distinct j1 and j2, the following may be executed
%   concurrently:
%
%     * O1 = F(N_{j1-1}, N_j1);
%
%     * CONSUME(N_{j2-1},N_j2, O2)
%
%   Note:  Even more so in this case, PARFOR is a convenient interface to this
%   functionality.  It uses CONSUME to implement its output sliced variables.
%
%
%   PARALLEL_FUNCTION(RANGE, F, CONSUME, SUPPLY) provides a way for the client
%   to send a worker a piece of data relevant only to its subinterval.  If
%   SUPPLY is [], the behavior is as above.  Otherwise, SUPPLY must be a
%   function handle (in practice, almost invariably to a nested function).
%   The client does the following:
%
%     * Evaluates I = SUPPLY(N_{j-1},N_j)
%
%     * Sends I to the worker along with N_{j-1} and N_j.
%
%   A worker evaluates F(N_{j-1}, N_j, I), with a result requested or not,
%   depending upon whether CONSUME is a function handle (as above).
%
%   Note: PARFOR uses SUPPLY for its sliced input variables.
%
%
%   R = PARALLEL_FUNCTION(RANGE, F, CONSUME, SUPPLY, REDUCE, IDENTITY)
%   A value of [] for REDUCE is ignored (same reason as for S), and will cause
%   IDENTITY to be ignored.  Otherwise, REDUCE must be a function handle and
%   IDENTITY must be present.  In this case, PARALLEL_FUNCTION returns a result,
%   computed as follows.  F is called in a worker with a "reduction" argument
%   and result (O is absent if CONSUME is [], and I is absent if SUPPLY is []):
%
%     * [O,Rout] = F(N_{j-1}, N_j, I, Rin);
%
%   For each call to PARALLEL_FUNCTION, the every time F is called in a worker,
%   Rin is IDENTITY.  Workers return their value of Rout to the client after each
%   evaluation of F, and do not retain any partially-reduced state in between
%   invocations of F.  The client combines all the partially-reduced states
%   computed on the workers.  The client combines these value by
%   REDUCE, pairwise as the results come in, obtaining a final value R that
%   becomes the result of PARALLEL_FUNCTION.
%
%   Evidently, REDUCE must take two arguments and produce one result.  Because
%   the order in which it is called is non-deterministic, REDUCE and IDENTITY
%   should obey these laws:
%
%     * identity:    REDUCE(r, IDENTITY) = REDUCE(IDENTITY, r) = r
%
%     * commutative: REDUCE(r1,r2) = REDUCE(r2,r1)
%
%     * associative: REDUCE(r1,REDUCE(r2,r3)) = REDUCE(REDUCE(r1,r2),r3)
%
%   Further, because of the unpredictability of the particular N_j and how the
%   subintervals are sent to workers, the following statement sequences must
%   produce the same value for r:
%
%     * [O1, r1] = F(N_{j-1}, N_j, I, IDENTITY);  % done in one worker
%       [O2, r2] = F(N_j, N_{j+1}, I2, IDENTITY); % done in another worker
%       r = REDUCE(r1, r2);                       % combined in the client
%
%     * [O, r] = F(N_{j1-1}, N_{j+1}, I1, IDENTITY); % done in one worker
%
%   REDUCE and F must also obey these laws:
%
%     * [O, r2] = F(N_{j-1}, N_j, I, r0); % start from arbitrary r0
%
%     * [O, r1] = F(N_{j-1}, N_j, I, IDENTITY); % start from IDENTITY
%       r2 = REDUCE(r0, r1);                    % combine with r0
%
%   A useful value of REDUCE is @plus, but there are others; see the PARFOR
%   documentation.  Also, note that for floating point numbers, @plus is only
%   approximately associative, so the errors introduced by it might vary from
%   run to run of the identical calls to PARALLEL_FUNCTION.
%
%   The server transmits the value of O and R to the client, as in the above
%   cases.
%
%   There is nothing you can do with REDUCE that cannot be done with CONSUME
%   alone.  However, the use of REDUCE can dramatically decrease communication
%   costs, because the size of its result may half the size of the sum of the
%   sizes of the inputs.
%
%
%   S = PARALLEL_FUNCTION(RANGE, F, CONSUME, SUPPLY, [], [], CONCAT, EMPTY)
%   [R, S] = PARALLEL_FUNCTION(RANGE, F, CONSUME, SUPPLY, REDUCE, IDENTITY,
%                              CONCAT, EMPTY)
%   CONCAT and EMPTY work in a way that is similar to REDUCE and IDENTITY,
%   except that CONCAT is not required to be commutative.  (The argument takes
%   its name from the fact that ordinary concatenation is associative but not
%   commutative.  Matrix multiply is another such function.)  CONCAT is given
%   arguments that come only from contiguous subintervals (N_{j-1},N_j] and
%   (N_j,N_{j+1}].  If CONCAT is non-empty, F is called in a worker with an
%   expected result, the we'll denote S, but unlike REDUCE, without a
%   corresponding input argument.
%
%     * [O, R, S] = F(N_{j-1}, N_j, I, R);
%
%   O is absent if CONSUME is [], I is absent if SUPPLY is [], and R is absent
%   (in both result and argument lists) if REDUCE is [].
%
%
%   ... = PARALLEL_FUNCTION(RANGE, F, CONSUME, SUPPLY, REDUCE, IDENTITY,
%                           CONCAT, EMPTY, M)
%   (Results are supplied by the whether REDUCE and CONCAT are empty, as above.)
%   The next argument is either a number M specifying the desired number of
%   workers.  The default value of M is inf, meaning "use as many workers as are
%   available, up to N."  If you specify M, it must be a positive integer or
%   inf, or some other valid means of specifying where to run the loop.
%   PARALLEL_FUNCTION will try to reserve min(N,M) workers.  In the
%   subsequent discussion W will denote the number of workers actually reserved,
%   so W <= min(N,M) <= N.
%
%   Valid values for M are:
%   - 0 - which means "run the loop in the client as quickly as possible"
%   - strictly positive numeric scalars - which means "run the loop on a parallel
%     pool of M (or fewer) workers (if one is already available or can be auto-created)"
%   - 'debug' - which means "run the loop in the client in a 'debug' mode"
%   - an instance of parallel.Pool, which explicitly runs the loop on the specified
%     pool
%   - an instance of parallel.Cluster, which explicitly runs the loop on the
%     specified cluster object, with no parallel pool required
%   - an instance of parallel.parfor.Options, which fully specifies where
%     to run the body of the loop.
%
%   In all cases, this argument, in combination with N, determines a set of
%   workers for this invocation of PARALLEL_FUNCTION.  At the end of the
%   computation, all reserved workers are released.
%
%
%   ... = PARALLEL_FUNCTION(RANGE, F, CONSUME, SUPPLY, REDUCE, IDENTITY,
%                           CONCAT, EMPTY, M, DIVIDE)
%   (Results are supplied by the whether REDUCE and CONCAT are empty, as above.)
%   The role of DIVIDE is to choose N_1, N_2,... .  If DIVIDE is numeric, it
%   must be scalar and at least 1; call this value D.  In this case, the N_i
%   are chosen by dividing RANGE into k = min(N,D*W) pieces as uniformly as
%   possible.  (The requirement that D is at least one guarantees that
%   W <= k <= N.)  If the computation time for each call on F is closely
%   proportional to the size of the interval for the call, then the best value
%   of D is 1.  But if the times are uneven, then larger values are better---
%   they keep all the workers busy, at the cost of more communication.
%
%   If DIVIDE is not number, it must be a function handle.  If it takes one
%   argument (as determined by NARGIN), then PARALLEL_FUNCTION calls it as
%   follows:
%
%     * k = min(N, max(W, DIVIDE(N))); % again, note W <= k <= N.
%
%   As before, the N_i are chosen to divide (BASE-1,LIMIT] as evenly as possible
%   into k subintervals.
%
%   If DIVIDE takes two arguments, it is called as follows:
%
%     * k = min(N, max(W, DIVIDE(N, W))); % again, note W <= k <= N.
%
%   The result k is used as before.
%
%   If DIVIDE does not take one or two arguments, it must take three, in which
%   case PARALLEL_FUNCTION calls it as follows:
%
%     * NN = DIVIDE(BASE-1, LIMIT, W)
%
%       - NN is a row vector, where NN(j) plays the role of N_j above.  Thus,
%         elements of NN must be positive, strictly increasing, and between
%         BASE and LIMIT.  The length of NN must be at least W.  If
%         NN(end) == LIMIT, then NN completely partitions (BASE,LIMIT], or
%         equivalently, the length of NN is what we have denoted by k.  In this
%         case, DIVIDE is allowed to be the last argument.  Otherwise, the next
%         argument is required; see below.
%
%
%   ... = PARALLEL_FUNCTION(N, F, CONSUME, SUPPLY, REDUCE, IDENTITY,
%                           CONCAT, EMPTY, M, DIVIDE, NEXT_DIVIDE)
%   (Results are supplied by the whether REDUCE and CONCAT are empty, as above.)
%   The purpose of NEXT_DIVIDE is to allow you to write an adaptive scheduler.  It
%   must be a function handle, and is called only when DIVIDE (which is called
%   only once per invocation of PARALLEL_DIVIDE) and all previous calls of
%   NEXT_DIVIDE have not specified the partition all the way to N.  NEXT_DIVIDE
%   takes an argument that supplies a summary of the computation so far, which is
%   what allows it to be adaptive.
%
%   The details are that PARALLEL_FUNCTION calls this argument as follows:
%
%     * NN_EXT = NEXT_DIVIDE(PREV_N, N, L_EXT)
%
%       - PREV_N is the endpoint of the last specified subinterval.
%
%       - N is the same as usual, so the task is to divide up (PREV_N,N].
%
%       - L_EXT specifies the minimum length of the result.  It is guaranteed that
%         N - PREV_N >= L_EXT.
%
%       - NN_EXT must be an extension to the partition specified so far.  It must
%         be a non-empty row vector of length at least L_EXT, where
%         PREV_N < NN_EXT(1), elements are strictly increasing, and
%         NN_EXT(end) <= N.  There will be no more calls on NEXT_DIVIDE if and
%         only if NN_EXT(end) == N.
%

% Indicate that we are entering this function
matlab.internal.incrementParallelFunctionDepth(1);
cleanupObject = onCleanup( @() matlab.internal.incrementParallelFunctionDepth(-1) );

if nargin < 2
    error(message('MATLAB:parfor:ArgumentMissing'))
end
% Validate range
if ~isequal(size(range), [1 2]) || ~isnumeric(range)
    error(message('MATLAB:parfor:InvalidArgumentRANGE'))
end
base = range(1)-1;
limit = range(2);
% Validate base and limit
if ~isreal(base) || ~isreal(limit)
    error(message('MATLAB:parfor:InvalidArgumentReal'))
end
if base ~= round(base) || limit ~= round(limit)
    error(message('MATLAB:parfor:InvalidArgumentInteger'))
end
N = double(limit - base); % Allows arbitrary arithmetic functions of N
% Validate F
if ~isa(F, 'function_handle')
    error(message('MATLAB:parfor:InvalidArgumentF'));
end
% The empty/function-handle status of consume, supply, and reduce is encoded in a single
% integer consume_supply_reduce as bit vector with mask bits:
consume_bit = 1;
supply_bit  = 2;
reduce_bit  = 4;
concat_bit  = 8;

consume_supply_reduce_concat = 0;
% Validate consume
if nargin >= 3 && ~isempty(consume)
    if ~isa(consume, 'function_handle')
        error(message('MATLAB:parfor:InvalidArgument', 'CONSUME'));
    end
    consume_supply_reduce_concat = consume_supply_reduce_concat + consume_bit;
end
% Validate supply
if nargin >= 4 && ~isempty(supply)
    if ~isa(supply, 'function_handle')
        error(message('MATLAB:parfor:InvalidArgument', 'SUPPLY'));
    end
    consume_supply_reduce_concat = consume_supply_reduce_concat + supply_bit;
end
% Validate reduce and identity
if nargin >= 5 && ~isempty(reduce)
    if ~isa(reduce, 'function_handle')
        error(message('MATLAB:parfor:InvalidArgument', 'REDUCE'));
    end
    if nargin == 5
        error(message('MATLAB:parfor:ArgumentNotSupplied', 'REDUCE', 'IDENTITY'))
    end
    consume_supply_reduce_concat = consume_supply_reduce_concat + reduce_bit;
else
    reduce = [];   % solely for the ...
    identity =  0; % ... benefit of ... distributed_execution
end
% Validate concat and empty
if nargin >= 7 && ~isempty(concat)
    if ~isa(concat, 'function_handle')
        error(message('MATLAB:parfor:InvalidArgument', 'CONCAT'));
    end
    if nargin == 7
        error(message('MATLAB:parfor:ArgumentNotSupplied', 'CONCAT', 'EMPTY'))
    end
    S = empty;
    consume_supply_reduce_concat = consume_supply_reduce_concat + concat_bit;
    concatenator = internal.parallel.parfor.online_concatenator(concat, empty);
else
    empty = []; % solely for the benefit of empty loops.
    concatenator = []; % solely for the benefit of distributed_execution
end

% M is already validated by parallel_M_check, simply set the value if none specified.
if nargin < 9
    M = Inf;
end

% Validate nargout
if nargout > 3
    error(message('MATLAB:parfor:InvalidOutputArgument'))
elseif nargout > 2 ...
        && (bitand(consume_supply_reduce_concat, reduce_bit) == 0 ...
            || bitand(consume_supply_reduce_concat, concat_bit) == 0)
    error(message('MATLAB:parfor:InvalidArgsForThreeOutputs'));
elseif nargout > 1 ...
        && bitand(consume_supply_reduce_concat, reduce_bit) == 0 ...
        && bitand(consume_supply_reduce_concat, concat_bit) == 0
    error(message('MATLAB:parfor:InvalidArgsForTwoOutputs'))
end

if N <= 0
    S = empty;
    R = identity;
    varargout = results(nargout);
    return
end

% Now that we know that the range is not empty, choose where to run based on M.

% First, resolve any raw Cluster or Pool objects to the corresponding Options object.
if isa(M, 'parallel.Cluster') || isa(M, 'parallel.Pool')
    M = parforOptions(M);
elseif isnumeric(M) && M > 0 && PCTInstalled
    % Resolve positive M into parforOptions.
    [isPoolRunning, pool] = pctTryCreatePoolIfNecessary();
    if isPoolRunning
        M = parforOptions(pool, 'MaxNumWorkers', M);
    else
        % No pool, and auto-create didn't create one - fall
        % back to the "0" case.
        M = 0;
    end
elseif isa(M, 'parallel.parfor.Options')
    % No action required
else
    % No PCT
    M = 0;
end

% There should be three remaining cases: numeric, 'debug', or Options object.
engineBuilder = []; % Implies use fall-back executor.
if isequal(M, 'debug')
    % Use the debug engine. Fabricate an effective number of workers.
    localW = max(1, floor(log2(N)));
    engineBuilder = @(varargin) localparfor(localW, varargin{:});
elseif isa(M, 'parallel.parfor.Options')
    % Build the engine from the options
    engineBuilder = @(initData, F) createEngine(M, initData, F);
else
    % Must be left with the M == 0 case.
    assert(isequal(M, 0), ...
           'Expected to find M==0.');
end

% Get the engine
if isempty(engineBuilder)
    engine = [];
else
    if logical(bitand(consume_supply_reduce_concat, reduce_bit))
        parfor_C = {consume_supply_reduce_concat, F, identity};
    else
        parfor_C = {consume_supply_reduce_concat, F};
    end
    try
        engine = engineBuilder({@make_channel, parfor_C}, F);
    catch E
        if strcmp( E.identifier, 'parallel:lang:parfor:IllegalComposite' )
            % In this case, we wish to abort parfor execution
            E2 = CatalogException(message('MATLAB:parfor:InvalidComposite'));
            E2 = addCause( E2, E );
            throwAsCaller( E2 );
        else
            % There are circumstances where we try to make a remoteparfor
            % controller and it throws an error (because it fails to acquire
            % the right resources for example). In that case we wish to fall
            % through to running locally
            warning(message('MATLAB:remoteparfor:ParforRunningLocally', E.getReport));
            engine = [];
        end
    end
end

% The code below is the short-circuit code to run the parfor on the local
% machine as fast as possible. The range is not divided up, it is simply
% run as a single interval.
if isempty(engine)
    try
        switch consume_supply_reduce_concat
          case 0
            F(base, limit);
          case 1
            consume(base, limit, F(base, limit));
          case 2
            F(base, limit, supply(base, limit));
          case 3
            consume(base, limit, F(base, limit, supply(base, limit)));
          case 4
            R = F(base, limit, identity);
          case 5
            [O, R] = F(base, limit, identity);
            consume(base, limit, O);
          case 6
            R = F(base, limit, supply(base, limit), identity);
          case 7
            [O, R] = F(base, limit, supply(base, limit), identity);
            consume(base, limit, O);
          case 8
            S = F(base, limit);
          case 9
            [O, S] = F(base, limit);
            consume(base, limit, O);
          case 10
            S = F(base, limit, supply(base, limit));
          case 11
            [O, S] = F(base, limit, supply(base, limit));
            consume(base, limit, O);
          case 12
            [R, S] = F(base, limit, identity);
          case 13
            [O, R, S] = F(base, limit, identity);
            consume(base, limit, O);
          case 14
            [R, S] = F(base, limit, supply(base, limit), identity);
          case 15
            [O, R, S] = F(base, limit, supply(base, limit), identity);
            consume(base, limit, O);
        end
    catch E % work around bug in cleaning up this workspace
        throwAsCaller(iBuildLocalParallelException(E));
    end
    varargout = results(nargout);
    return
end

% Validate divide and next_divide after call to run locally as this isn't used in the local case
W = engine.getNumWorkers();
if nargin >= 10 && ~isempty(divide)
    if isequal(divide, 'debug')
        divide = divide_uniform(@(N, W) max(W, round(N/floor(sqrt(N)))));
    elseif isnumeric(divide)
        if ~isscalar(divide) || divide < 1 || isnan(divide)
            error(message('MATLAB:parfor:InvalidArgumentNumericDIVIDE'))
        end
        divide = divide_uniform(@(N, W) min(N, double(divide)*W));
    elseif isa(divide, 'function_handle')
        switch nargin(divide)
          case 1
            divide = divide_uniform_min_max(@(N, W) divide(N));
          case 2
            divide = divide_uniform_min_max(divide);
          case 3
            % Use divide as it came in.
          otherwise
            error(message('MATLAB:parfor:InvalidArgumentFHandleDIVIDE'))
        end
    else
        error(message('MATLAB:parfor:InvalidArgumentDIVIDE'));
    end
    NN = divide(base, limit, W);
else
    % Use the options structure to divide the range
    NN = engine.divideRangeIntoIntervals(base, limit, W);
end

% Validate next_divide
if nargin >= 11 && ~isempty(next_divide)
    if ~isa(next_divide, 'function_handle')
        error(message('MATLAB:parfor:InvalidArgument', 'NEXT_DIVIDE'));
    end
end

% The partition:
k = length(NN);
NN0 = [base, NN(1:k-1)];     % index by j to get the jth base
if ~all(NN0 < NN) || NN(k) > limit
    error(message('MATLAB:parfor:DivideNotIncreasing'))
end

try
    R = distributed_execution(...
        engine, base, limit, W, k, NN, NN0, consume, supply, ...
        reduce, identity, concatenator, ...
        consume_supply_reduce_concat, ...
        consume_bit, supply_bit, reduce_bit, concat_bit);
catch err
    if ~isa(err, 'ParallelException')
        err = iBuildLocalParallelException(err);
    end
    throwAsCaller(err);
end

if bitand(consume_supply_reduce_concat, concat_bit) ~= 0
    S = concatenator.final();
end
varargout = results(nargout);

    % Provides the varargout value for parallel_function
    function C = results(nresults)
    if bitand(consume_supply_reduce_concat, reduce_bit) ~= 0 % always result R in this case
        if bitand(consume_supply_reduce_concat, concat_bit) ~= 0
            assert(nresults <= 2, 'Only up to 2 outputs supported.');
            if nresults == 2
                C = {R, S};
            else
                C = {R};
            end
        else
            assert(nresults <= 1, 'Only 1 output supported for this case.');
            C = {R};
        end
    elseif bitand(consume_supply_reduce_concat, concat_bit) ~= 0 % always result S in this case
        assert(nresults <= 1, 'Only 1 output supported for this case.');
        C = {S};
    else
        assert(nresults == 0, 'No outputs supported for this case.');
        C = {};
    end
    end

end

%----------------------------------------------------------------------------
%
%----------------------------------------------------------------------------
function d = divide_uniform(f)
    d = @div;
    function NN = div(base, limit, W)
        N = limit - base;
        k = f(N, W);
        NN = base + round((1:k)*N/k);
    end
end

function d = divide_uniform_min_max(f)
    d = divide_uniform(@(N, W) min(N, max(W, f(N, W))));
end

function NN = divide_harmonic(base, limit, W)
N = limit - base;
% Iterate for the while loop
i = 0;
% Minimum chunk size such that there are no more than 10*W iterates.
minChunk = max(ceil(N/(10*W)), 1);
maxChunk = max(ceil(N/(1.5*W)), 1);
curr = 0;
% Output size guess - heuristic - assert that N>0 & W>0
outputSize = ceil(10*W);
% Allocate output array - ensure that it is the correct class
NN = zeros(1, outputSize, 'like', base);
while curr < N
    i = i+1;
    curr = curr + min(max(ceil((N - curr)/W), minChunk), maxChunk);
    NN(i) = curr;
end
% Force end point to be N as required - NOTE NN(i-1) is
% necessarily less than N by the loop condition so we are able to
% guarantee that NN is monotonically increasing.
NN(i) = N;
% And trim to expected output
NN = base + NN(1:i);
end

%----------------------------------------------------------------------------
%
%----------------------------------------------------------------------------
function R = distributed_execution(...
               P, base, limit, W, k, NN, NN0, consume, supply, ...
               reduce, R, concatenator, ...
               consume_supply_reduce_concat, ...
               consume_bit, supply_bit, reduce_bit, concat_bit)

% It is important that we make these logical as we are going to cast the
% values to double and assume that the values are zero or one (a bit of a hack that
% makes the code a little cleaner)
doConsume = logical(bitand(consume_supply_reduce_concat, consume_bit));
doSupply  = logical(bitand(consume_supply_reduce_concat, supply_bit));
doReduce  = logical(bitand(consume_supply_reduce_concat, reduce_bit));
doConcat  = logical(bitand(consume_supply_reduce_concat, concat_bit));

% As we will see, a channel here is a function handle, which is called
% as follows:
%   channel(distcompserialize({base, limit, I})) -- send the next subinterval
%     * I is omitted if supply is [].
%     * Unless consume is [], channel responds by yielding distcompserialize(O).
%   channel(distcompserialize({}))               -- tell the worker to close
%     * This yields no results.

totalReceived = 0; % total intervals received
% How many intervals to dispatch in one go - and receive in one go
[initialDispatchSize, normalDispatchSize] = P.getDispatchSizes(k, W);
assert(initialDispatchSize > 0 && initialDispatchSize <= k, ...
    'Incorrect initial dispatch size: %d', initialDispatchSize);

thisDispatchSize = min(k, initialDispatchSize);

noLoopErrorDetected = true;
j = 0; % this indexes the subinterval
prev = base;
try
    % Loop whilst there is work to be done
    while prev < limit
        if noLoopErrorDetected
            for i = 1:thisDispatchSize
                j = j+1; % index of the next subinterval.
                next = NN(j);
                if doSupply
                    subinterval_C = {prev, next, supply(prev, next)};
                else
                    subinterval_C = {prev, next};
                end
                noLoopErrorDetected = P.addInterval(j, subinterval_C);
                prev = next;
                % Did we successfully add the interval - an error
                % might have occurred on a remote machine. Don't add
                % any more as there is an error that will be
                % thrown during getCompleteIntervals
                if ~noLoopErrorDetected, break, end
            end
        end
        % Setup for the next loop before getting the results
        thisDispatchSize = min(k-j, normalDispatchSize);
        % Get the results from the computation - remember that tag is an index
        % into the allLimits array
        [tags, out] = P.getCompleteIntervals(min(thisDispatchSize, j));
        totalReceived = totalReceived + numel(tags);
        % Now consume and/or concat the latest results.
        for i = 1:numel(tags)
            if doConsume || doConcat
                consume_reduce(tags(i), out{i, 1});
            end
            if doReduce
                R = reduce(R, out{i, 2});
            end
        end
    end

    % Notify that all intervals have been added.
    P.allIntervalsAdded();

    % Number of intervals that need to be picked up for consumption
    numToConsume = (j - totalReceived);
    % Define a maximum number to consume or reduce
    maxSize = round(W/floor(sqrt(W)));
    while numToConsume > 0
        % What chunk size should we take this time
        chunkSize = min(maxSize, numToConsume);
        % Get some of the remaining intervals
        [tags, out] = P.getCompleteIntervals(chunkSize);
        % Loop over these intervals doing the right thing
        for i = 1:numel(tags)
            thisTag = tags(i);
            % If we are consuming and this tag is a normal dispatch tag
            % then consume this interval and decrement the number to
            % consume
            if thisTag <= j
                if doConsume || doConcat
                    consume_reduce(thisTag, out{i, 1});
                end
                if doReduce
                    R = reduce(R, out{i, 2});
                end
                numToConsume = numToConsume - 1;
            end
        end
    end
catch err
    P.complete(true);
    rethrow(err);
end
% Tell this remoteparfor that we have completed
P.complete(false);

    function consume_reduce(j, out)
        % Called only when doConsume || doConcat is true
        if doConsume
            if doConcat
                concatenator.concat(j, out{2})
                out = out{1};
            end
            consume(NN0(j), NN(j), out);
        else
            concatenator.concat(j, out)
        end
    end

end

%----------------------------------------------------------------------------
%
%----------------------------------------------------------------------------
function channel = make_channel(C)

consume_bit = 1;
reduce_bit  = 4;
concat_bit  = 8;

consume_supply_reduce_concat = C{1};
F = C{2};

% Need to know if we are doing consume or reduce to deal with the
% number of output arguments from the channel function.
doConsume = bitand(consume_supply_reduce_concat, consume_bit);
doReduce  = bitand(consume_supply_reduce_concat, reduce_bit);
doConcat  = bitand(consume_supply_reduce_concat, concat_bit);

if numel(C) < 3
    identity = [];
else
    identity = C{3};
end
channel = make_general_channel(F, identity, doConsume, doReduce, doConcat);

end

function channel = make_general_channel(F, identity, doConsume, doReduce, doConcat)
    function OR = channel_general(C)
        % Exit condition that sends back the reduce variable if requested
        if isempty(C)
            % 'final' intervals are no longer used.
            OR = [];
            return;
        end
        % Otherwise we are calling F with the appropriate input and output
        % arguments
        R = [];
        if doConsume
            if doReduce
                if doConcat
                    [O, R, S] = F(C{:}, identity);
                    O = {O, S}; % bundle S with O for shipping to client
                else
                    [O, R] = F(C{:}, identity);
                end
            else
                if doConcat
                    [O, S] = F(C{:});
                    O = {O, S}; % bundle S with O for shipping to client
                else
                    O = F(C{:});
                end
            end
        else
            O = [];
            if doReduce
                if doConcat
                    [R, O] = F(C{:}, identity); % ship "S" to client, so call it O
                else
                    R = F(C{:}, identity);
                end
            else
                if doConcat
                    O = F(C{:}); % ship "S" to client, so call it O
                else
                    F(C{:});
                end
            end
        end
        % Always return O and R - workers do not store partially-reduced values, and
        % always start again from the identity.
        OR = {O, R};
    end
    channel = @channel_general;
end

%----------------------------------------------------------------------------
% Convenience wrappers areound PCT installation/licence checks.
function OK = PCTInstalled
OK = matlab.internal.parallel.isPCTInstalled;
end
function OK = PCTLicensed
OK = matlab.internal.parallel.isPCTLicensed;
end
%----------------------------------------------------------------------------


%----------------------------------------------------------------------------
% Implementation of localparfor - this is intended to be a local
% approximation to the job of distcomp.remoteparfor, ParforControllerImpl
% and remoteParallelFunction. Most of the code here is stolen from the
% distcomp.remoteparfor and remoteParallelFunction.
%----------------------------------------------------------------------------
function P = localparfor(W, initData, ~)
P = struct('addInterval', @nAddInterval, ...
           'getNumWorkers', @() W, ...
           'allIntervalsAdded', @() [], ...
           'getCompleteIntervals', @nGetCompleteIntervals, ...
           'divideRangeIntoIntervals', @divide_harmonic, ...
           'getDispatchSizes', @(k, W) deal(k, 2*W), ...
           'complete', @(varargin)[]); % 'complete' does nothing for local execution.
% Hold the channels in a cell array
channels = cell(W, 1);
% Define the current channel index to use;
currentChannelIndex = 1;
% Do we want to try serialization
if PCTInstalled && PCTLicensed
    send = @serialize;
else
    send = @noserialize;
end
% This queue will hold the list of intervals to do - we will fill it up
% from the addInterval and addFinalInterval methods
pendingQueue = cell(0, 2);

    function OK = nAddInterval(tag, varargin)
        OK = true;
        % Add this interval to the queue
        pendingQueue(end+1, :) = {tag varargin};
    end

    function [tags, results] = nGetCompleteIntervals(numIntervals)
        % Pre-allocate the output from the function
        results = cell(numIntervals, 2);
        tags = ones(numIntervals, 1);
        % Loop over the requested intervals
        for i = 1:numIntervals
            % First make the channel if it doesn't exists - this ensures
            % that any error during channel construction is thrown in
            % getCompleteIntervals as it would do in a true distributed
            % parfor.
            if isempty(channels{currentChannelIndex})
                C = send(initData);
                channels{currentChannelIndex} = feval(C{1}, C{2:end});
            end
            % Get the first item from the stack of pending intervals
            [tags(i), data] = deal(pendingQueue{1, :});
            % Remove it
            pendingQueue(1, :) = [];
            % Send to the channel, evaluate and send back
            channelArg = send(data);
            try
                out = feval(channels{currentChannelIndex}, channelArg{:});
            catch err
                ex = iBuildLocalParallelException(send(err));
                throw(ex);
            end
            results(i, :) = send(out);
            % Make sure we use the next available channel for the next
            % interval
            currentChannelIndex = mod(currentChannelIndex, W) + 1;
        end
    end
end

function x = noserialize(x)
end

function x = serialize(x)
bufH = distcompMakeByteBufferHandle( distcompserialize( x ) );
x = distcompdeserialize(distcompByteBuffer2MxArray( bufH.get ));
bufH.free;
end

% ------------------------------------------------
function me = CatalogException(message)
% Helper function for MException

me = MException(message.Identifier, '%s', message.getString);

end

function newEx = iBuildLocalParallelException(originalErr)
% Create a new exception that stitches together the client and worker stack
% but that doesn't include this file in the stack.
newEx = ParallelException.hBuildLocalParallelException(originalErr, mfilename);
end
