classdef ProcessQueue < handle & bioinfo.pipeline.internal.VersionMixin
    %

    %   Copyright 2020-2021 The MathWorks, Inc.

    properties (Dependent)
        OutputData
        RunErrors
        RunStatus
        RunStart
        RunEnd
    end

    properties (SetAccess = private)
        InputData(1,1) bioinfo.pipeline.internal.PipelineData

        % Copy of the settings used during the last run block. When the
        % processes in this queue are run. These properties are used to
        % determine whether the ProcessQueue is stale, see ISSTALE().
        LastRunProperties(1,1) struct = struct('Block', [], 'InputData', [], 'ConnectionDict', [], 'BlockName', []);
    end

    properties (Constant, Access = private)
        TimeStampFileName = '.block_run.ts';
    end

    properties (Access = private)
        % When a process queue is (re)-compiled, store which inputs are
        % connected. Which inputs are connected determines whether a given
        % input row is complete or not, if there are optional inputs that
        % are not connected, they will be ignored. If an optional input is
        % connected, the row must wait until that data is available.
        ExpectedInputs(:,1) string = string.empty

        % Private versions of run outputs
        RunStructs cell = {};

        Merged(1,1) logical = false;

        % Private versions of public properties
        privateRunErrors cell = {}
        privateRunStatus bioinfo.pipeline.RunStatus = bioinfo.pipeline.RunStatus.empty;
        privateRunStart datetime
        privateRunEnd datetime
        privateOutputData(1,1) bioinfo.pipeline.internal.PipelineData
        privateHasMergeFailure
    end

    properties (Access = private, Transient)
        RunFutures = parallel.Future.empty;
        DispatchFutures = parallel.Future.empty;
    end

    properties
        % Handle to Block in pipeline
        Block bioinfo.pipeline.Block = bioinfo.pipeline.Block.empty;
        BlockName(1,1) string
    end

    methods
        function obj = ProcessQueue(blk, blockName)
            arguments
                blk(1,1) bioinfo.pipeline.Block
                blockName(1,1) string = ""
            end

            obj.Block = blk;
            obj.BlockName = blockName;

            obj.InputData = bioinfo.pipeline.internal.PipelineData();
            obj.privateOutputData = bioinfo.pipeline.internal.PipelineData();
        end

        function refresh(obj)
            obj.refreshPrivateProperties();
        end

        function value = get.OutputData(obj)
            obj.refreshPrivateProperties();
            value = obj.privateOutputData;
        end

        function value = get.RunErrors(obj)
            obj.refreshPrivateProperties();
            value = obj.privateRunErrors;
        end

        function value = get.RunStatus(obj)
            obj.refreshPrivateProperties();
            value = obj.privateRunStatus;
        end

        function value = get.RunStart(obj)
            obj.refreshPrivateProperties();
            value = obj.privateRunStart;
        end

        function value = get.RunEnd(obj)
            obj.refreshPrivateProperties();
            value = obj.privateRunEnd;
        end

        function obj = insertInput(obj, varName, varValue)
            obj.InputData.insert(varName, varValue);
        end

        function tf = isRunning(obj)
            % running/queued/waiting
            tf = any(ismember(obj.RunStatus, [bioinfo.pipeline.RunStatus.Waiting bioinfo.pipeline.RunStatus.Queued bioinfo.pipeline.RunStatus.Running]), "all");
        end

        function tf = isComplete(obj)
            tf  = ~isempty(obj.RunStatus) && all(obj.RunStatus == bioinfo.pipeline.RunStatus.Completed, "all") && obj.Merged;
        end

        function tf = hasMergeFailure(obj)
            tf = obj.privateHasMergeFailure;
        end

        function tf = isStale(obj)
            % check if block and Input data have changed since the last
            % run. Note that this does NOT capture if incoming connections
            % have changed, which would necessitate re-runing the block
            tf = obj.Block.isStaleComparedTo(obj.LastRunProperties.Block) || ...
                ~isequaln(obj.LastRunProperties.InputData, obj.InputData) || ...
                ~isequaln(obj.LastRunProperties.BlockName, obj.BlockName);
        end

        function ps = getParallelSize(obj)
            % Get the size of the processes queues parallelization on the
            % most recent run
            ps = size(obj.RunFutures);
        end

        function run(obj, pool, dispatcher, baseDir, saveBlkResults, connectionDict)
            arguments
                obj
                pool(:,1) parallel.Pool
                dispatcher(1,1) bioinfo.pipeline.internal.Dispatcher
                baseDir(1,1) string
                saveBlkResults(1,1) logical
                connectionDict(1,1) dictionary
            end

            % Clear output data
            obj.privateOutputData = bioinfo.pipeline.internal.PipelineData();

            % Split input data into block args, run, then rejoin
            try
                runArgs = obj.createArgs(dispatcher);
                parallelSize = size(runArgs);

                dispatcher.debug(message("bioinfo:ProcessQueue:BlockRunArgsCreated", mat2str(parallelSize)), "BlockRun", obj.Block.UUID);

                % Only clear errors if runArgs creation was successful
                obj.privateRunErrors = cell(parallelSize);
                obj.privateRunErrors(:) = {MException.empty};
                obj.privateRunStatus = repmat(bioinfo.pipeline.RunStatus.Waiting, parallelSize);

            catch createArgsError
                runArgs = [];
                parallelSize = [1 1];

                % If args creation failed, set that as the runtime error
                obj.privateRunErrors = {createArgsError};
                obj.privateRunStatus = bioinfo.pipeline.RunStatus.Error;

                msg = message("bioinfo:ProcessQueue:InvalidBlockRunArgs", obj.BlockName, createArgsError.message);
                dispatcher.error(msg, "BlockRunError", obj.Block.UUID, "Throw", false);
            end

            % Init properties even if createArgs failed, but don't run
            obj.RunFutures = repmat(parallel.FevalFuture(), parallelSize);
            obj.DispatchFutures = repmat(parallel.FevalFuture(), parallelSize);
            obj.privateRunStart = repmat(datetime(NaT, 'TimeZone', 'local'), parallelSize);
            obj.privateRunEnd = obj.RunStart;
            obj.RunStructs = cell(parallelSize);
            obj.RunStructs(:) = {struct.empty};
            obj.Merged = false;
            obj.privateHasMergeFailure = false;

            % Cache copy of runtime settings. This allows isStale() to
            % determine when this process queue must be rerun
            obj.LastRunProperties.Block = copy(obj.Block);
            obj.LastRunProperties.BlockName = obj.BlockName;
            obj.LastRunProperties.InputData = copy(obj.InputData);
            obj.LastRunProperties.ConnectionDict = connectionDict;

            % Also drop a timestamp file in the Block Dir
            bioinfo.pipeline.internal.ProcessQueue.writeTimestamp(baseDir);

            % If splitting results in 0 run structs, warn and insert empty
            % results
            expectedOutputFields = string(fieldnames(obj.Block.Outputs));

            if any(parallelSize == 0)
                for outputName = expectedOutputFields
                    if obj.Block.Outputs.(outputName).UniformOutput
                        obj.privateOutputData.insert(outputName, []);
                    else
                        obj.privateOutputData.insert(outputName, {});
                    end
                end
                obj.privateRunStatus = bioinfo.pipeline.RunStatus.Completed;
                obj.Merged = true; % No run results to merge
            end

            if isempty(runArgs)
                return;
            end

            obj.createBlockSubFolders(baseDir, numel(runArgs));

            % For Futures, private status fields will be updated each
            % time the public version is requested
            numRuns = numel(runArgs);
            blkName = obj.BlockName;
            for r = 1:numRuns
                inStruct = runArgs(r);
                blk = obj.Block;
                blkDir = fullfile(baseDir, num2str(r));


                if isempty(pool)
                    % Serial case. Avoid using in-process parfeval for the block running,
                    % because serial parfeval duplicates figures in fetchOutputs, and because
                    % UDD objects are not yet supported correctly (g2891157).
                    % buggy. Unlike the parallel case Status, Start, End, and Errors must be
                    % updated here since they are not dependent of the futures returned by
                    % parfeval. Similarly, ensure the struct returned by run() is inserted into
                    % RunStructs, as it will not need to be retrieved via fetchOutputs in
                    % mergeOutputStructs.
                    try
                        obj.privateRunStatus(r) = bioinfo.pipeline.RunStatus.Running;
                        obj.privateRunStart(r) = datetime('now', 'TimeZone', 'local');

                        % Run the following clean-up operations if user issues CTRL-C during block run:
                        % - Set RunStatus to Error
                        % - dispatch error
                        %
                        % otherwise block will think it is still running
                        interruptCleanup = bioinfo.pipeline.internal.Cleanup(@() obj.handleRunInterrupt(r, dispatcher, blkName));

                        obj.RunStructs{r} = bioinfo.pipeline.internal.ProcessQueue.blockRunWrapper(blk, inStruct, r, expectedOutputFields, dispatcher, blkName, numRuns, saveBlkResults, blkDir);

                        dispatcher.info(message("bioinfo:ProcessQueue:BlockRunFinished", blkName, r, numRuns), "BlockRunFinish", obj.Block.UUID);
                        obj.privateRunStatus(r) = bioinfo.pipeline.RunStatus.Completed;
                    catch BlkRunError
                        dispatcher.error(BlkRunError, "BlockRunError", obj.Block.UUID, Throw=false);
                        obj.privateRunErrors{r} = BlkRunError;
                        obj.privateRunStatus(r) = bioinfo.pipeline.RunStatus.Error;

                    end
                    removeTask(interruptCleanup);

                    obj.privateRunEnd(r) = datetime('now', 'TimeZone', 'local');

                else
                    obj.RunFutures(r) = parfeval(pool, @() bioinfo.pipeline.internal.ProcessQueue.blockRunWrapper(blk, inStruct, r, expectedOutputFields, dispatcher, blkName, numRuns, saveBlkResults, blkDir), 1);                    
                    obj.DispatchFutures(r) = afterEach(obj.RunFutures(r), @(runFevalFuture) dispatchAfterEachRun(runFevalFuture, dispatcher, blkName, r, numRuns, blk.UUID), 0, "PassFuture", true);
                    obj.privateRunStatus(r) = bioinfo.pipeline.RunStatus.Queued;
                end
            end

            obj.RunFutures = reshape(obj.RunFutures, size(runArgs));

            function dispatchAfterEachRun(runFevalFuture, dispatcher, blkName, runIdx, numRuns, blkUUID)
                % FevalFuture states: 'pending' | 'queued' | 'running' | 'finished' | 'failed' | 'unavailable'

                if ~isempty(runFevalFuture.Error)
                    % future can be 'finished', but with Error
                    dispatcher.error(runFevalFuture.Error, "BlockRunError", blkUUID);
                end
                dispatcher.info(message("bioinfo:ProcessQueue:BlockRunFinished", blkName, runIdx, numRuns), "BlockRunFinish", blkUUID);
            end
        end

        function handleRunInterrupt(obj, r, dispatcher, blkName)
            % called if run is terminated by Ctrl-C

            runInteruptedME = MException(message("bioinfo:ProcessQueue:BlockRunInterrupted", blkName));
            dispatcher.error(runInteruptedME, "BlockRunError", obj.Block.UUID, Throw=false);
            obj.privateRunErrors{r} = runInteruptedME;
            obj.privateRunStatus(r) = bioinfo.pipeline.RunStatus.Error;
        end

        function compile(obj, expectedInputs)
            obj.ExpectedInputs = expectedInputs;
        end

        function outFuture = afterAll(obj, postBlockFinishFcn, dispatcher)
            % Run blkRunAfterFinishFcn after all futures in the queue have
            % completed. Note that unlike regular afterAll, the given
            % function will always run but does not take any futures as
            % input
            arguments
                obj(1,1) bioinfo.pipeline.internal.ProcessQueue
                postBlockFinishFcn(1,1) function_handle % must not return anything
                dispatcher
            end

            futures = obj.DispatchFutures;

            isUnavailable = strcmp({futures.State}, "unavailable");
            futures(isUnavailable) = [];

            if isempty(futures)
                % Create a no-op future to run after, since afterAll does
                % not work with empty futures arrays. This is needed in
                % case createArgs fails or produces an empty array of input
                % arguments
                futures = parfeval(parallel.Pool.empty, @() [], 0);
            end

            outFuture = afterAll(futures, @(~) mergeThenPostBlockFinishFcn(), 0, 'PassFuture', true);

            function mergeThenPostBlockFinishFcn()

                futures; outFuture; %#ok<VUNUS> % avoid R2024b nested function incompatibility by referencing variable explicitly; see g3089026
                
                allCompleted = ~isempty(obj.RunStatus) && all(arrayfun(@(status) isequal(status, bioinfo.pipeline.RunStatus.Completed), obj.RunStatus(:)));

                if allCompleted
                    try
                        obj.mergeOutputStructs();
                        msg = message("bioinfo:ProcessQueue:BlockRunFinishedAll", obj.BlockName);
                        dispatcher.info(msg, "BlockRunFinishAll", obj.Block.UUID);

                    catch ME
                        obj.privateHasMergeFailure = true;
                        msg = message("bioinfo:ProcessQueue:OutputCantMerge");
                        dispatcher.error(msg, "BlockRunError", obj.Block.UUID, Throw=false);
                    end
                else
                    msg = message("bioinfo:ProcessQueue:BlockRunFinishedAll", obj.BlockName);
                    dispatcher.info(msg, "BlockRunFinishAll", obj.Block.UUID);
                end

                postBlockFinishFcn();
            end
        end

        function cancel(obj)
            cancel(obj.RunFutures);

            % Status/Errors will be updated by refreshPrivateProperties
        end

        function clear(obj)
            obj.RunFutures = parallel.Future.empty;
            obj.RunStructs = {};
            obj.Merged = false;


            obj.privateRunErrors = {};
            obj.privateRunStatus = bioinfo.pipeline.RunStatus.empty;
            obj.privateRunStart(:) = [];
            obj.privateRunEnd(:) = [];

            obj.privateRunStart.TimeZone = 'local';
            obj.privateRunEnd.TimeZone = 'local';

            % Clear input & output data
            obj.InputData = bioinfo.pipeline.internal.PipelineData();
            obj.privateOutputData = bioinfo.pipeline.internal.PipelineData();

            obj.LastRunProperties = struct('Block', [], 'InputData', [], 'ConnectionDict', []);
        end

        function wait(obj)
            futures = obj.DispatchFutures;

            isUnavailable = strcmp({futures.State}, "unavailable");
            futures(isUnavailable) = [];

            wait(futures);
        end

        function [resultsStruct, outputsModified] = results(obj)
              [outputsToRetrieve, outputsModified] = obj.lastRunPorts("Outputs");

              if isComplete(obj)
                  resultsStruct = struct(obj.OutputData, outputsToRetrieve);
              else
                  % Results Incomplete, return sentinel value instead of struct.
                  resultsStruct = bioinfo.pipeline.datatypes.Incomplete;
              end
        end

        function [inputsStruct, inputsModified] = inputs(obj)

            [inputsToRetrieve, inputsModified] = obj.lastRunPorts("Inputs");
            inputsStruct = struct(obj.InputData, inputsToRetrieve);
        end
    end


    methods (Access = private)

        function [lastRunPortNames, portsModified] = lastRunPorts(obj, type)
            % lastRunPorts helper to return the names of the ports of the
            % last block run. Also return whether the Block's current
            % definition differs, so that Pipeline can issue a warning
            arguments
                obj
                type(1,1) string {mustBeMember(type, ["Inputs", "Outputs"])}
            end

            lastRunBlock = obj.LastRunProperties.Block;

             if isempty(lastRunBlock)
                 lastRunPortNames = fieldnames(obj.Block.(type));
                 portsModified = false;
             else
                 lastRunPortNames = fieldnames(lastRunBlock.(type));

                 % Check if block definition has changed. Warn from
                 % Pipeline Object
                 portsModified = ~isequal(lastRunPortNames, fieldnames(obj.Block.(type)));
             end

        end

        function refreshPrivateProperties(obj)
            arguments
                obj
            end

            % Update results of any pending futures
            if ~isempty(obj.RunFutures)
                futuresState = {obj.RunFutures.State};
                futuresState = reshape(futuresState, size(obj.RunFutures));

                runningFutures = strcmp('running', futuresState);

                % Update status of running processes
                obj.privateRunStatus(runningFutures) = bioinfo.pipeline.RunStatus.Running;

                % Update the start time
                if any(runningFutures)
                    obj.privateRunStart(runningFutures) = [obj.RunFutures(runningFutures).StartDateTime];
                end

                assumedRunningProcesses = (obj.privateRunStatus == bioinfo.pipeline.RunStatus.Running | ...
                    obj.privateRunStatus == bioinfo.pipeline.RunStatus.Queued);
                completedFutures = (strcmp('finished', futuresState) | strcmp('failed', futuresState));

                % Note that any processes that ran in serial have already updated the
                % private properties, and their "futures" will remain 'unavailable'.
                toRefresh = find(assumedRunningProcesses & completedFutures);

                for i = toRefresh(:)'
                    future = obj.RunFutures(i);

                    % If it failed, store the errors, and set status to failed
                    if ~isempty(future.Error)
                        obj.privateRunErrors{i} = future.Error;
                        obj.privateRunStatus(i) = bioinfo.pipeline.RunStatus.Error;
                    else
                        obj.privateRunStatus(i) = bioinfo.pipeline.RunStatus.Completed;
                    end

                    % Finally, update completion time
                    obj.privateRunStart(i) = future.StartDateTime;
                    obj.privateRunEnd(i) = future.FinishDateTime;
                end
            end
        end

        function runArgs = createArgs(obj, dispatcher)
            % Create a struct array of eval args based on the parallelized
            % dimensions in Block and in PipelineData

            % First, split input data values based on any parallel
            % dimensions
            inputVarNames = obj.ExpectedInputs;

            % Simple case, no input variables, the run struct is a scalar
            % struct with no fields
            if isempty(inputVarNames)
                runArgs = struct();
                return;
            end

            inputVarValues = obj.InputData.values(inputVarNames);
            inputSplitDims = arrayfun(@(inputVarName) obj.Block.Inputs.(inputVarName).SplitDimension, inputVarNames, 'UniformOutput', false);

            % Calculate maxdim, to create a matrix for the split sizes for
            % each input
            maxdim = max(cellfun(@ndims, inputVarValues)); % Max dimension of the inputs
            explicitDims = inputSplitDims(cellfun(@(c) ~isequal(c, "all"), inputSplitDims)); % Filter out "all"
            maxdim = max([maxdim, [explicitDims{:}]]); % Take maximium of explicit dimensions

            inputSplitSizes = ones(numel(inputVarValues), maxdim); % SplitSizes for each input/dim pair up to maxdim, default is 1

            % Calculate split size for each input/dim pair, expanding
            % "all", while warning for scalars
            for iv = 1:numel(inputVarNames)
                if isequal(inputSplitDims{iv}, "all")

                    if isscalar(inputVarValues{iv})
                        % Warn that input is scalar
                        msg = message("bioinfo:ProcessQueue:ScalarInput", obj.BlockName, inputVarNames(iv));
                        dispatcher.warn(msg, "BlockRunWarn", obj.Block.UUID);
                        warning(msg);
                    end

                    % Split Size is the size of the input value across all
                    % dimensions
                    inputSplitSizes(iv, :) = size(inputVarValues{iv}, 1:maxdim);
                    inputSplitDims{iv} = 1:maxdim;
                elseif ~isempty(inputSplitDims{iv})
                    % Otherwise, split size is size of input in each split
                    % dimension
                    inputSplitSizes(iv, inputSplitDims{iv}) = size(inputVarValues{iv}, inputSplitDims{iv});

                    % Warn if size is scalar for any of these dims
                    warningScalarDims = inputSplitSizes(iv, inputSplitDims{iv}) == 1;
                    warningScalarDims = inputSplitDims{iv}(warningScalarDims);

                    if ~isempty(warningScalarDims)
                        for wd = warningScalarDims(:)'
                            msg = message("bioinfo:ProcessQueue:ScalarSplitDimension", obj.BlockName, inputVarNames(iv), wd);
                            dispatcher.warn(msg, "BlockRunWarn", obj.Block.UUID);
                            warning(msg);
                        end
                    end
                end

                % Warn if any split dimensions are size 0
                warningEmptyDim = find(inputSplitSizes(iv, :) == 0, 1);
                if ~isempty(warningEmptyDim)
                    msg = message("bioinfo:ProcessQueue:EmptySplitInput", obj.BlockName, inputVarNames(iv), warningEmptyDim);
                    dispatcher.warn(msg, "BlockRunWarn", obj.Block.UUID);
                    warning(msg);
                end
            end

            % Check Size compatibility.
            % Each dimension must have (at most) one unique non-scalar size
            outputSize = ones(1, maxdim);

            for d = 1:maxdim
                dSizes = inputSplitSizes(:,d);

                % There must be at most 1 non-scalar size in this dimension
                dSizes(dSizes == 1) = [];
                dSizes = unique(dSizes);

                if numel(dSizes) > 1
                    input1Name = inputVarNames(find(inputSplitSizes(:,d) == dSizes(1), 1));
                    input2Name = inputVarNames(find(inputSplitSizes(:,d) == dSizes(2), 1));

                    msg = message("bioinfo:ProcessQueue:IncompatibleSizes", obj.BlockName, d, input1Name, dSizes(1), input2Name, dSizes(2));
                    dispatcher.error(msg, "BlockRunError", obj.Block.UUID);
                elseif numel(dSizes) == 1
                    % If there is a non-scalar size, then this is the
                    % output size
                    outputSize(d) = dSizes(1);
                end
            end

            splitVars = cell(1, numel(inputVarNames));

            % Split each variable into a cell array of the output size
            for iv = 1:numel(inputVarNames)
                try
                    splitVars{iv} = bioinfo.pipeline.internal.ProcessQueue.splitVariable(inputVarValues{iv}, outputSize, inputSplitDims{iv});
                catch splitExceptionCause
                    splitException = MException(message("bioinfo:ProcessQueue:SplitVariableFailed", obj.BlockName, inputVarNames(iv), mat2str(inputSplitDims{iv})));
                    splitException = addCause(splitException, splitExceptionCause);
                    dispatcher.error(splitException, "BlockRunError", obj.Block.UUID);
                end
            end

            % Finally, create a struct using the variable cell arrays
            structArgs = cell(1, 2*numel(inputVarNames));

            structArgs(1:2:2*numel(inputVarNames)) = cellstr(inputVarNames);
            structArgs(2:2:1+2*numel(inputVarNames)) = splitVars;

            runArgs = struct(structArgs{:});
        end

        function mergeOutputStructs(obj)
            % If this was run in parallel, extract the returned structs
            % from RunFutures. Futures will be unavailable if they were placeholders
            % before running in series, and RunStructs has already been populated with
            % the outputs from the blockRunWrappers.
            if ~isempty(obj.RunFutures) && ~any(strcmp('unavailable', {obj.RunFutures.State}))
                obj.RunStructs = arrayfun(@fetchOutputs, obj.RunFutures, 'UniformOutput', false);
            end

            sz = size(obj.RunStructs);
            varNames = fieldnames(obj.RunStructs{1});

            pipelineData = bioinfo.pipeline.internal.PipelineData();
            for vn = 1:numel(varNames)

                varValue = cellfun(@(s) s.(varNames{vn}), obj.RunStructs, 'UniformOutput', false);
                varValue = reshape(varValue, sz);

                % For UniformOutputs, elements must be scalar
                if obj.LastRunProperties.Block.Outputs.(varNames{vn}).UniformOutput

                    % UniformOutput shouldn't be required if we are merging
                    % a scalar, so just extract the value. Otherwise ensure
                    % that the values are scalar and use cellfun to extract
                    if isequal(sz, [1 1])
                        varValue = varValue{1};
                    elseif any(cellfun(@(c) ~isscalar(c), varValue))
                        error(message("bioinfo:ProcessQueue:NonScalarOutput", varNames{vn}));
                    else
                        % Output already guaranteed to be uniform
                        varValue = cellfun(@(c) c, varValue, 'UniformOutput', true);
                    end
                end

                pipelineData.insert(varNames{vn}, varValue);
            end

            obj.privateOutputData = pipelineData;

            obj.Merged = true;
        end

        function createBlockSubFolders(obj, baseDir, nDirs)
            arguments
                obj
                baseDir
                nDirs(1,1) double {mustBePositive, mustBeInteger}
            end

            for d = 1:nDirs
                subdir = fullfile(baseDir, num2str(d));
                if exist(subdir, 'dir') ~= 7
                    mkdir(subdir);
                end
            end

        end
    end

    methods (Static, Access = private)
        function outStruct = blockRunWrapper(block, inStruct, runIdx, expectedOutputFields, dispatcher, blkName, numRuns, saveBlkResults, blkDir)
            % Wrapper around Block.run() that changes into the correct
            % directory, and in the event of an error in block.run, returns
            % the corresponding value from the error handler function.

            origDir = cd(blkDir);
            cdCleanup = onCleanup(@() cd(origDir));

            dispatcher.info(message("bioinfo:ProcessQueue:RunningBlock", blkName, runIdx, numRuns), "BlockRunStart", block.UUID); % add runIdx here in msg

            % blow away previous results, if they exist
            % RAM and FileSystem must always be in sync for a given run
            blockResultsFile = fullfile(blkDir, "blockResults.mat");
            if isfile(blockResultsFile)
                delete(blockResultsFile);
            end

            try
                outStruct = block.run(inStruct);

            catch runME
                if isempty(block.ErrorHandler)
                    rethrow(runME);

                else
                    errHandlerStruct.identifier = runME.identifier;
                    errHandlerStruct.message = runME.message;
                    errHandlerStruct.index = runIdx;
                    try
                        outStruct = block.ErrorHandler(errHandlerStruct, inStruct);

                        % Validate outStruct still has the expected fields
                        if ~isstruct(outStruct) || ~isscalar(outStruct) || numel(intersect(fieldnames(outStruct), expectedOutputFields)) ~= numel(expectedOutputFields)
                            msg = message("bioinfo:ProcessQueue:InvalidErrorHandler", strjoin(expectedOutputFields, ", "));
                            error(msg)
                        end

                    catch errHandlerME
                        rethrow(errHandlerME);

                    end
                end
            end

            % Save Results if specified
            if saveBlkResults
                % save all outputs as individual variables in a mat file
                % (including figures). This will allow users to load
                % specific outputs they may want
                warnState = warning("off", "MATLAB:Figure:FigureSavedToMATFile");
                warnCleanup = onCleanup(@() warning(warnState));

                save(blockResultsFile, "-struct", "outStruct");
            end
        end

        function writeTimestamp(inDir)
            % WRITETIMESTAMP write a timestamp file in the given directory.
            arguments
                inDir(1,1) string {mustBeFolder}
            end

            now = datetime('now', 'TimeZone', 'local');

            timeStampFile = fullfile(inDir, bioinfo.pipeline.internal.ProcessQueue.TimeStampFileName);
            tsFID = fopen(timeStampFile, 'w');
            if tsFID == -1
                error(message("bioinfo:ProcessQueue:CouldntWriteTimestamp"));
            end

            fprintf(tsFID, '%.15f', datenum(now));
            fclose(tsFID);
        end

        function splitData = splitVariable(varValue, outputSize, parallelDims)
            numDims = ndims(varValue);

            mat2cellArgs = cell(1,numDims);

            % For any dimension that is parallel, the mat2cell arg for that
            % dimension should be ones(1, size(varValue, d)), indicating
            % that that dimension is split. For any non-parallel dimension,
            % there should be a scalar of size(varValue, d), so that value
            % will not be split.

            for d = 1:numDims
                if ismember(d, parallelDims)
                    mat2cellArgs{d} = ones(1, size(varValue, d));
                else
                    mat2cellArgs{d} = size(varValue, d);
                end

            end

            % mat2cell requires paren-indexing, so for scalars just wrap
            % value as is in a cell. If this throws an exception, will be
            % caught and dispatched in createArgs
            if ~isequal(mat2cellArgs, {1, 1})
                splitData = mat2cell(varValue, mat2cellArgs{:});
            else
                splitData = {varValue};
            end

            % After parallelization is taken care of, scalar expand across
            % any other dimensions of output Size

            repMatArgs = cell(1,numDims);

            for d = 1:numDims
                if outputSize(d) == size(splitData, d)
                    repMatArgs{d} = 1;
                else
                    % Scalar expand in this dimension
                    repMatArgs{d} = outputSize(d);
                end
            end

            splitData = repmat(splitData, repMatArgs{:});
        end
    end

    methods (Static)
        function dt = readTimestamp(inDir)
            %READTIMESTAMP read the timestamp file in the given direcory.
            % Returns the datetime of the last modification date, or NaT if
            % the file does not exist.
            arguments
                inDir(1,1) string {mustBeFolder}
            end

            dt = NaT('TimeZone', 'local'); % Default if no timestamp file exists

            timeStampFile = fullfile(inDir, bioinfo.pipeline.internal.ProcessQueue.TimeStampFileName);

            tsInfo = dir(timeStampFile);
            if numel(tsInfo) == 1 && ~ tsInfo.isdir
                try
                    dtDatenum = fileread(timeStampFile);
                    dt = datetime(str2double(dtDatenum), 'ConvertFrom', 'datenum');
                catch ME
                    warning(ME.message, ME.identifier); % Warn when failing to read timestamp correctly, but return NaT
                end
            end
        end
    end
end
