classdef (ConstructOnLoad) Pipeline < bioinfo.pipeline.internal.SemanticModelObject
    %


    properties (Dependent)
        % BLOCKS Vector of block objects contained in the pipeline.
        %
        % See also BlockNames, addBlock, removeBlock, findBlock
        Blocks

        % BlockNames Vector of names for the blocks in the pipeline.
        %
        % See also renameBlock, addBlock, removeBlock, findBlock
        BlockNames
    end

    properties (Hidden)
        %BLOCKGRAPH Digraph storing Blocks (Nodes) and Connections (Edges).
        BlockGraph(1,1) digraph
    end

    properties (Hidden, Transient, Access = {?bioinfo.pipeline.internal.designer.PipelineController, ...
            ?bioinfo.pipeline.internal.Logger, ?BiopipelineEventLogger}, NonCopyable)
        MessageBroker(1,1) bioinfo.pipeline.internal.MessageBroker
    end

    properties (Access = private, NonCopyable)
        % RunContext the options used the last time the pipeline was run. Make-like
        % behavior is subject to reusing the same results folder.
        RunContext bioinfo.pipeline.internal.RunContext {mustBeScalarOrEmpty} = bioinfo.pipeline.internal.RunContext.empty
    end

    properties (Access = {?trunPipeline, ?tPipeline}, NonCopyable, Transient)
        Logger bioinfo.pipeline.internal.Logger
    end

    properties (Access = private, NonCopyable, Transient)
        Dispatcher  bioinfo.pipeline.internal.Dispatcher
        BlockSemanticModelChangeListeners(:,1) event.listener
        % RunData Data and metadata associated with a run of a pipeline. 
        RunData bioinfo.pipeline.internal.RunData {mustBeScalarOrEmpty} = bioinfo.pipeline.internal.RunData.empty
    end

    methods
        function obj = Pipeline()
            % PIPELINE Create Pipeline object.
            %
            % P = PIPELINE creates an empty pipeline P.
            %
            % Example:
            %
            % import bioinfo.pipeline.blocks.*
            % import bioinfo.pipeline.Pipeline
            %
            % P = Pipeline();
            % FCB = FileChooser(which("ex1.sam"));
            % SSB = SamSort;
            %
            % addBlock(P, [FCB, SSB]);
            % connect(P, FCB, SSB, ["Files", "SAMFile"]);
            %
            % run(P);
            % results(P, SSB)
            %
            % See also addBlock, removeBlock, connect, disconnect, run

            obj.MessageBroker   = bioinfo.pipeline.internal.MessageBroker;
            obj.Dispatcher      = bioinfo.pipeline.internal.Dispatcher(obj.MessageBroker.MessageQueue);
            obj.Logger          = bioinfo.pipeline.internal.Logger(obj.MessageBroker);

            import bioinfo.pipeline.Block
            import bioinfo.pipeline.Connection

            nodeTable = table(cell(0,1), Block.empty(0,1), 'VariableNames', ["Name", "Block"]);
            edgeTable = table(double.empty(0,2), Connection.empty(0,1), 'VariableNames', ["EndNodes", "Connection"]);
            obj.BlockGraph = digraph(edgeTable,nodeTable);
        end

        function blocks = addBlock(obj, blocks, names)
            % addBlock Add blocks to pipeline.
            %
            % Bout = addBlock(P, B) add the block B to the pipeline P. B is
            % a block object or a vector of block objects. Each unnamed
            % block in B is given a unique name, and the blocks are
            % returned in Bout.
            %
            % addBlock(P, B, N) adds each block in B to the pipeline with a
            % corresponding name N. If B is a single block, N is a
            % character vector or string scalar. If B contains multiple
            % blocks, N is a string vector or cell array of character
            % vectors with the same number of elements as B. Each block in
            % B will be assigned the corresponding name in N.
            %
            % Example:
            %
            % import bioinfo.pipeline.blocks.FileChooser
            %
            % P = bioinfo.pipeline.Pipeline;
            %
            % addBlock(P, FileChooser);
            % addBlock(P, FileChooser, "FileChooser2");
            % P.Blocks
            %
            % See also findBlock, removeBlock, renameBlock, connect, disconnect
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                blocks(:,1) bioinfo.pipeline.Block {mustBeNonempty}
                names(:,1) string = [];
            end

            % First ensure input blocks are unique
            if numel(blocks) ~= numel(unique(blocks))
                error(message("bioinfo:Pipeline:NonUniqueBlocks"));
            end

            % Ensure no block is the same as a block already in the
            % pipeline
            if size(obj.BlockGraph.Nodes, 1) > 0 && ~isempty(intersect(obj.BlockGraph.Nodes.Block(:), blocks))
                error(message("bioinfo:Pipeline:BlockAlreadyExists"));
            end

            % Assign default names if none are given
            if nargin < 3
                names = obj.assignBlockName(blocks);
            elseif ~all(arrayfun(@isvarname, names))
                error(message("bioinfo:Pipeline:InvalidBlockNames"));
            end

            % Enforce uniqueness of names
            if numel(names) ~= numel(unique(names))
                error(message("bioinfo:Pipeline:NonUniqueBlockNames"));
            end

            existingNames = string.empty;
            if ~isempty(obj.BlockGraph.Nodes)
                existingNames = convertStringsToChars(obj.BlockGraph.Nodes.Name);
            end
            conflicts = intersect(existingNames, names);

            if ~isempty(conflicts)
                error(message("bioinfo:Pipeline:BlockNameAlreadyExists", strjoin(conflicts, ', ')));
            end

            if numel(names) ~= numel(blocks)
                error(message("bioinfo:Pipeline:InvalidNumNames"));
            end

            % Listen to changes in the block's semantic model (i.e. ports)
            semanticModelListeners = arrayfun(@(blk) listener(blk, 'SemanticModelChanged', @(src, evt) obj.handleBlockModelChanged(evt)), blocks);

            blockProps = table(names, blocks, 'VariableNames', ["Name", "Block"]);
            obj.BlockGraph = addnode(obj.BlockGraph, blockProps);
            obj.BlockSemanticModelChangeListeners = [obj.BlockSemanticModelChangeListeners; semanticModelListeners];

            % Notify that new blocks were added
            modelChangedData = bioinfo.pipeline.internal.SemanticModelChangedEventData;
            for b = 1:numel(blocks)
                modelChangedData.add(blocks(b), obj, names(b))
            end

            if ~isempty(modelChangedData.Delta)
                notify(obj, 'SemanticModelChanged', modelChangedData);
            end

            obj.restoreRunData(blocks); % Restore run data for these blocks, if they are valid in this pipeline
        end

        function removeBlock(obj, blockIDs)
            % removeBlock Remove blocks from pipeline.
            %
            % removeBlock(P, B) removes the blocks specified by B from the
            % pipeline P. B is a block object, block name, or vector of
            % block objects or names.
            %
            % Any connections to or from B are also removed from P.
            %
            % Example:
            % import bioinfo.pipeline.blocks.FileChooser
            %
            % FCB1 = FileChooser(which("ex1.bam"));
            % FCB2 = FileChooser(which("ex1.sam"));
            % FCB3 = FileChooser(which("Myco_1_1.sam"));
            %
            %
            % P = bioinfo.pipeline.Pipeline;
            %
            % addBlock(P, [FCB1, FCB2, FCB3], ["F1", "F2", "F3"]);
            %
            % removeBlock(P, "F2");
            % P.Blocks
            %
            % removeBlock(P, [FCB1, FCB3]);
            % P.Blocks
            %
            % See also findBlock, addBlock, renameBlock, connect,
            % disconnect
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                blockIDs {mustBeNonempty, bioinfo.pipeline.internal.Validation.mustBeBlockID(blockIDs)}
            end

            toRemove = findBlockIndex(obj, blockIDs);

            % Find all connections that will also be deleted, to update
            % semantic model
            deletedConnections = arrayfun(@(bidx) [inedges(obj.BlockGraph, bidx); outedges(obj.BlockGraph, bidx)], toRemove, 'UniformOutput',false);
            deletedConnections = unique(vertcat(deletedConnections{:}));

            modelChangedData = bioinfo.pipeline.internal.SemanticModelChangedEventData;
            blocksToRemove = obj.BlockGraph.Nodes.Block(toRemove);

            for b = toRemove(:)'
                blkName = string(obj.BlockGraph.Nodes.Name{b});
                modelChangedData.remove(obj.BlockGraph.Nodes.Block(b), obj, blkName);
            end

            for c = deletedConnections(:)'
                modelChangedData.remove(obj.BlockGraph.Edges.Connection(c), obj, "");
            end
            
            % Remove and notify
            delete(obj.BlockSemanticModelChangeListeners(toRemove));
            obj.BlockSemanticModelChangeListeners(toRemove) = [];
            obj.BlockGraph = rmnode(obj.BlockGraph, toRemove);
            
            % Remove ProcessQueues
            if ~isempty(obj.RunData)
                pqsToRemove = blocksToRemove(isKey(obj.RunData.ProcessQueues, blocksToRemove));
                obj.RunData.ProcessQueues(pqsToRemove) = [];
            end
           
            if ~isempty(modelChangedData.Delta)
                notify(obj, 'SemanticModelChanged', modelChangedData);
            end
        end

        function postConnectMap = connect(obj, sourceBlock, targetBlock, newPortMap)
            % CONNECT Connect ports in pipeline.
            %
            % PORTMAP = CONNECT(P, SB, TB, NEWPORTMAP) connects the ports
            % in NEWPORTMAP between a source and a target block, SB and TB,
            % in P. SB and TB are block objects or block names. NEWPORTMAP
            % is an N-by-2 string vector. Its first column contains the
            % names of the output ports of SB to connect, and the second
            % column contains the names of the corresponding input ports of
            % TB. For example, PORTMAP = ["X", "M"] would indicate that the
            % output "X" of SB is mapped to the input "M" of TB. PORTMAP is
            % an array showing all connections between SB and TB after
            % adding the elements in NEWPORTMAP.
            %
            %
            % Example:
            %
            % import bioinfo.pipeline.blocks.FileChooser
            % import bioinfo.pipeline.blocks.SamSort
            %
            % P = bioinfo.pipeline.Pipeline();
            % FCB = FileChooser(which("ex1.sam"));
            % SSB = SamSort;
            %
            % addBlock(P, [FCB, SSB]);
            % connect(P, FCB, SSB, ["Files", "SAMFile"])
            %
            % See also disconnect, portMap, addBlock, removeBlock
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                sourceBlock {bioinfo.pipeline.internal.Validation.mustBeBlockID(sourceBlock, 1)}
                targetBlock {bioinfo.pipeline.internal.Validation.mustBeBlockID(targetBlock, 1)}
                newPortMap(:,2) string {mustBeNonempty}
            end

            srcBlockIdx = findBlockIndex(obj, sourceBlock);
            srcBlockObj = obj.Blocks(srcBlockIdx);

            tgtBlockIdx = findBlockIndex(obj, targetBlock);
            tgtBlockObj =  obj.Blocks(tgtBlockIdx);

            if srcBlockObj == tgtBlockObj
                error(message("bioinfo:Pipeline:NoSelfLoops"));
            end

            % Validate that all target port names are unique
            if numel(unique(newPortMap(:,2))) ~= size(newPortMap, 1)
                error(message("bioinfo:Pipeline:NonUniqueTargetPortNames"));
            end

            % Check that there are no existing connections to the requested port maps
            connInputs = obj.connectedInputs(tgtBlockIdx);
            conflictingInputs = intersect(newPortMap(:,2), connInputs);
            if numel(conflictingInputs) > 0
                error(message("bioinfo:Pipeline:InputAlreadyConnected", strjoin(conflictingInputs, ', ')));
            end

            modelChangedData = bioinfo.pipeline.internal.SemanticModelChangedEventData();
            for i=size(newPortMap,1):-1:1
                conn(i) = bioinfo.pipeline.Connection(srcBlockObj, tgtBlockObj, newPortMap(i,:));
                modelChangedData.add(conn(i), obj, "");
            end

            edgeProps = table(conn(:), 'VariableNames', "Connection");

            % expand srcBlockIdx & tgtBlockIdx
            srcBlockIdx = repmat(srcBlockIdx, size(conn));
            tgtBlockIdx = repmat(tgtBlockIdx, size(conn));

            % Add new connections to the graph, verify no cycles are
            % introduced
            newGraph = addedge(obj.BlockGraph, srcBlockIdx, tgtBlockIdx, edgeProps);

            if ~isdag(newGraph)
                error(message("bioinfo:Pipeline:CyclicPipeline"));
            end
            obj.BlockGraph = newGraph;

            postConnectMap = portMap(obj, srcBlockObj, tgtBlockObj);

            notify(obj, 'SemanticModelChanged', modelChangedData);
        end

        function postDisconnectMap = disconnect(obj, srcBlockID, targetBlockID, portMapToRemove)
            % DISCONNECT Remove connection between ports in pipeline.
            %
            % DISCONNECT(P, SB, TB) removes all connections between the
            % blocks specified by SB and TB in the pipeline P. SB and TB
            % are block objects or block names.
            %
            % PM = DISCONNECT(P, SB, TB, RC) removes the connections
            % between SB and TB indicated by the port map RC. PM contains
            % the port map for any remaining connections between SB and TB
            % after disconnecting the connections in RC.
            %
            % Example:
            % import bioinfo.pipeline.blocks.FileChooser
            % import bioinfo.pipeline.blocks.SamSort
            %
            % P = bioinfo.pipeline.Pipeline();
            % FCB = FileChooser(which("ex1.sam"));
            % SSB = SamSort;
            %
            % addBlock(P, [FCB, SSB]);
            % connect(P, FCB, SSB, ["Files", "SAMFile"])
            %
            % disconnect(P, FCB, SSB, ["Files", "SAMFile"])
            %
            % See also connect, portMap, addBlock, removeBlock
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                srcBlockID {bioinfo.pipeline.internal.Validation.mustBeBlockID(srcBlockID, 1)}
                targetBlockID {bioinfo.pipeline.internal.Validation.mustBeBlockID(targetBlockID, 1)}
                portMapToRemove(:,2) string {mustBeNonempty} = "all";
            end

            if nargin == 3
                toRemove = findConnectionIndex(obj, srcBlockID, targetBlockID);
                if toRemove == 0
                    error(message("bioinfo:Pipeline:NoConnectionsBetweenBlocks"));
                end
            else
                toRemove = findConnectionIndex(obj, srcBlockID, targetBlockID, portMapToRemove);
                if any(toRemove == 0)
                    error(message("bioinfo:Pipeline:ConnectionsNotFound"));
                end
            end


            removedConnections = obj.BlockGraph.Edges.Connection(toRemove);
            modelChangedData = bioinfo.pipeline.internal.SemanticModelChangedEventData();
            for c = 1:numel(toRemove)
                modelChangedData.remove(removedConnections(c), obj, "");
            end

            obj.BlockGraph = rmedge(obj.BlockGraph, toRemove);
            postDisconnectMap = portMap(obj, srcBlockID, targetBlockID);

            notify(obj, 'SemanticModelChanged', modelChangedData);

        end

        function block = renameBlock(obj, blockID, name)
            % renameBlock Rename block in pipeline.
            %
            % Bout = renameBlock(P, B, N) sets the name of block B in the
            % pipeline P to N. B is a block object or block name and N is
            % the new name as a character vector or string scalar. Bout is
            % the block object that was renamed.
            %
            % Example:
            %
            % P = bioinfo.pipeline.Pipeline
            % FCB = bioinfo.pipeline.blocks.FileChooser(which("ex1.bam"));
            % addBlock(P, FCB, "File");
            % renameBlock(P, "File", "BamFile")
            % renameBlock(P, FCB, "Ex1Bam")
            %
            % See also findBlock, addBlock, removeBlock, connect,
            % disconnect
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                blockID {bioinfo.pipeline.internal.Validation.mustBeBlockID(blockID, 1)}
                name(1,1) string
            end

            if ~isvarname(name)
                error(message("bioinfo:Pipeline:InvalidBlockNames"));
            end

            k = findBlockIndex(obj, blockID);
            oldName = string(obj.BlockGraph.Nodes.Name{k});

            if oldName == name
                return
            end

            if findnode(obj.BlockGraph, char(name)) ~= 0
                error(message('bioinfo:Pipeline:NonUniqueBlockNames'));
            end

            obj.BlockGraph.Nodes.Name{k} = char(name);
            block = obj.BlockGraph.Nodes.Block(k);

            if ~isempty(obj.RunData) && isKey(obj.RunData.ProcessQueues, block)
                pq = obj.RunData.ProcessQueues(block);
                pq.BlockName = string(name);
            end

            modelChangedData = bioinfo.pipeline.internal.SemanticModelChangedEventData;
            modelChangedData.rename(block, obj, [oldName, name]);

            notify(obj, 'SemanticModelChanged', modelChangedData);

        end

        function blockNames = get.BlockNames(obj)
            blockNames = string.empty(0,1);

            if size(obj.BlockGraph.Nodes, 1) > 0
                blockNames = string(obj.BlockGraph.Nodes.Name);
            end
        end

        function name = blockName(obj, block)
            % blockName Return name of specified block in pipeline.
            %
            % NAME = blockName(P, B) returns the name of the block B in the
            % pipeline P. B is a block object or array of block objects.
            % If B is an array, then the name of each block is returned,
            % and the output array has the same size as the input.
            %
            % P = bioinfo.pipeline.Pipeline;
            % FCB = bioinfo.pipeline.blocks.FileChooser(which("ex1.bam"));
            % addBlock(P, FCB, "MyFile");
            %
            % fcbName = blockName(P, FCB)
            %
            % See also BlockNames, addBlock, renameBlock, removeBlock
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                block bioinfo.pipeline.Block
            end
            k = findBlockIndex(obj, block);

            name = string(obj.BlockGraph.Nodes.Name);
            name = reshape(name(k), size(block));
        end

        function outStruct = results(obj, blockID)
            % RESULTS Get pipeline results.
            %
            % R = RESULTS(P, B) returns the results R for all the processes
            % associated with block B in the pipeline P. B is a block
            % object or block name. R is a struct whose fieldnames are the
            % output names of B and values are the output values from the
            % RUN method of the block. If the pipeline has not been run, or
            % the results are not available, R is a
            % bioinfo.pipeline.datatypes.Incomplete object.
            %
            %
            % Example:
            %
            % import bioinfo.pipeline.blocks.FileChooser
            % import bioinfo.pipeline.blocks.SamSort
            %
            % P = bioinfo.pipeline.Pipeline();
            % FCB = FileChooser(which("ex1.sam"));
            % SSB = SamSort;
            %
            % addBlock(P, [FCB, SSB]);
            % connect(P, FCB, SSB, ["Files", "SAMFile"]);
            %
            % run(P);
            % results(P, FCB)
            % results(P, SSB)
            %
            % See also inputs, fetchResults, run, processTable
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                blockID {bioinfo.pipeline.internal.Validation.mustBeBlockID(blockID, 1)}
            end
            
            k = findBlockIndex(obj, blockID);
            blk = obj.BlockGraph.Nodes.Block(k);

            obj.restoreRunData(k); % Restore data if it exists 

            % ProcessQueue for this block doesn't exist, return Incomplete
            if ~isKey(obj.RunData.ProcessQueues, blk)
                outStruct = bioinfo.pipeline.datatypes.Incomplete;
                return;
            end

            pq = obj.RunData.ProcessQueues(blk);
            [outStruct, ~] = pq.results();
        end

        function outStruct = fetchResults(obj, blockID)
            % fetchResults Wait for and retrieve pipeline results.
            %
            % fetchResults is equivalent to RESULTS, except that if
            % parallelization is enabled, the function waits for the block
            % to complete before returning.
            %
            % Example:
            %
            % import bioinfo.pipeline.blocks.FileChooser
            % import bioinfo.pipeline.blocks.SamSort
            %
            % P = bioinfo.pipeline.Pipeline();
            % FCB = FileChooser(which("ex1.sam"));
            % SSB = SamSort;
            %
            % addBlock(P, [FCB, SSB]);
            % connect(P, FCB, SSB, ["Files", "SAMFile"]);
            %
            % run(P, 'UseParallel', true);
            % fetchResults(P, FCB)
            %
            %
            % See also wait, results, run, processTable
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                blockID {bioinfo.pipeline.internal.Validation.mustBeBlockID(blockID, 1)}
            end
            wait(obj, blockID);
            outStruct = obj.results(blockID);
        end

        function procTable = processTable(obj, blockID, NVPArgs)
            % processTable Returns information about all processes in pipeline.
            %
            % T = processTable(P) returns a table T of information about
            % all processes in a pipeline.
            %
            % T = processTable(P, B) returns process information only for
            % the blocks specified by B in P. B is a block object, block name,
            % or vector of block objects or names.
            %
            % processTable(..., 'Expanded', TF) specifies whether to expand
            % the parallel executions for each block into separate rows.
            % The default is false.
            %
            % Example:
            % import bioinfo.pipeline.blocks.FileChooser
            % import bioinfo.pipeline.blocks.SamSort
            %
            % P = bioinfo.pipeline.Pipeline();
            % FCB = FileChooser(which("ex1.sam"));
            % SSB = SamSort;
            %
            % addBlock(P, [FCB, SSB]);
            % connect(P, FCB, SSB, ["Files", "SAMFile"]);
            %
            % run(P);
            %
            % T  = processTable(P)
            % T2 = processTable(P, SSB)
            %
            % See also run, results, cancel
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                blockID {mustBeNonempty, bioinfo.pipeline.internal.Validation.mustBeBlockID(blockID)} = "unspecified";
                NVPArgs.Expanded(1,1) logical = false;
            end

            % Restore run data, then populate missing process queues to ensure each
            % block has a processqueue
            obj.restoreRunData();

            procTable = table('Size', [0, 5], ...
                'VariableTypes', ["string", repmat("cell", 1, 4)], ...
                'VariableNames', {'Block', 'Status', 'RunStart', 'RunEnd', 'RunErrors'});

            if isempty(obj.RunData)
                return
            end

            if nargin == 1
                blocksToRetrieve = 1:size(obj.BlockGraph.Nodes,1);
                if isempty(blocksToRetrieve)
                    return
                end
            else
                blocksToRetrieve = findBlockIndex(obj, blockID);
            end

            blockTables = cell(1, numel(blocksToRetrieve));

            for b = 1:numel(blocksToRetrieve)
                blk = obj.BlockGraph.Nodes.Block(blocksToRetrieve(b));

                Block = string(obj.BlockGraph.Nodes.Name(blocksToRetrieve(b)));

                if isKey(obj.RunData.ProcessQueues, blk)
                    pq = obj.RunData.ProcessQueues(blk);
                    Status = {pq.RunStatus};
                    RunStart = {pq.RunStart};
                    RunEnd = {pq.RunEnd};
                    RunErrors = {pq.RunErrors};
                else
                    localEmptyDatetime = datetime.empty;
                    localEmptyDatetime.TimeZone = 'local';
                    Status = {bioinfo.pipeline.RunStatus.empty};
                    RunStart = {localEmptyDatetime};
                    RunEnd = {localEmptyDatetime};
                    RunErrors = {MException.empty};
                end

                if NVPArgs.Expanded
                    Status = Status{1}(:);
                    RunStart = RunStart{1}(:);
                    RunEnd = RunEnd{1}(:);
                    RunErrors = RunErrors{1}(:);

                    if isempty(Status)
                        continue % Expanded tables have no rows for blocks that are not run
                    end
                    Block = repmat(Block, numel(Status), 1);
                end

                blockTables{b} = table(Block, Status, RunStart, RunEnd, RunErrors);

            end

            % Filter empty tables, then concatenate
            blockTables(cellfun(@isempty, blockTables)) = [];
            procTable = vertcat(procTable, blockTables{:});
        end

        function blocks = get.Blocks(obj)
            blocks = findBlock(obj);
        end

        function blocks = findBlock(obj, blockID)
            % findBlock Access blocks in pipeline.
            %
            % B = findBlock(P) returns an array B of the blocks in the
            % pipeline P.
            %
            % B = findBlock(P, BN) returns only the blocks in P identified
            % by BN. BN is a block name or vector of block names.
            %
            % Example:
            % import bioinfo.pipeline.blocks.FileChooser
            %
            % FCB1 = FileChooser(which("ex1.bam"));
            % FCB2 = FileChooser(which("ex1.sam"));
            % FCB3 = FileChooser(which("Myco_1_1.sam"));
            %
            % P = bioinfo.pipeline.Pipeline;
            % addBlock(P, [FCB1, FCB2, FCB3], ["F1", "F2", "F3"]);
            %
            % B1 = findBlock(P)
            % B2 = findBlock(P, ["F2", "F3"])
            %
            % See also addBlock, removeBlock, Blocks, renameBlock, portMap
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                blockID {mustBeNonempty, bioinfo.pipeline.internal.Validation.mustBeBlockID(blockID)} = "-all";
            end

            blocks = obj.BlockGraph.Nodes;

            if size(blocks,1) > 0
                blocks = blocks.Block;
            else
                blocks = bioinfo.pipeline.Block.empty(0,1);
            end

            k = findBlockIndex(obj, blockID);
            blocks = blocks(k,:);

        end

        function pm = portMap(obj, srcID, tgtID)
            %PORTMAP Show connected ports between two blocks.
            %
            % PM = PORTMAP(P, SB, TB) shows the connected ports between the
            % blocks in P identified by SB and TB. SB and TB are block
            % objects or block names. PM is an N-by-2 string array where N
            % is the number of connected ports. Rows of PM indicate
            % connections, where PM(i,1) is the name of an output port of
            % SB connected to the input port PM(i,2) of TB.
            %
            % Example:
            % import bioinfo.pipeline.blocks.FileChooser
            % import bioinfo.pipeline.blocks.SamSort
            %
            % P = bioinfo.pipeline.Pipeline();
            % FCB = FileChooser(which("ex1.sam"));
            % SSB = SamSort;
            % addBlock(P, [FCB, SSB]);
            % connect(P, FCB, SSB, ["Files", "SAMFile"]);
            %
            % portMap(P, FCB, SSB)
            %
            % See also connect, disconnect, findBlock, addBlock,
            % removeBlock
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                srcID {bioinfo.pipeline.internal.Validation.mustBeBlockID(srcID,1)}
                tgtID {bioinfo.pipeline.internal.Validation.mustBeBlockID(tgtID,1)}
            end

            k = findConnectionIndex(obj, srcID, tgtID);
            pm = string.empty(0,2);
            if k ~= 0
                conns = obj.BlockGraph.Edges.Connection;
                pm = vertcat(conns(k).PortMap);
            end

        end

        function compile(obj, inputStruct)
            % COMPILE Verify pipeline structure and prepare for running.
            %
            % COMPILE(P) analyzes the blocks and connections in P and
            % reports any errors or warnings.
            %
            % COMPILE(P, S) analyzes the blocks and connections in P with
            % input structure S. The field names of S must match unconnected
            % ports in P.
            %
            % Blocks may implement additional compile time checks with
            % their own COMPILE method.
            %
            % Example:
            % import bioinfo.pipeline.blocks.FileChooser
            % import bioinfo.pipeline.blocks.SamSort
            %
            % P = bioinfo.pipeline.Pipeline();
            % FCB = FileChooser(which("ex1.sam"));
            % SSB = SamSort;
            %
            % addBlock(P, [FCB, SSB]);
            % connect(P, FCB, SSB, ["Files", "SAMFile"]);
            %
            % compile(P)
            %
            % See also run, addBlock, connect,
            % bioinfo.pipeline.Block.compile
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                inputStruct(1,1) struct = struct();
            end

            obj.restoreRunData();
            createMissingProcessQueues(obj);

            obj.Dispatcher.info(message("bioinfo:Pipeline:CompilingPipeline"), "PipelineCompileStart", obj.UUID);

            % Visit graph in topological order, and propagate data and process queue sizes
            %   downstream. Source Blocks have a process queue size of 1, and
            %   output DataQueues of size 1. All downstream blocks calculate queue
            %   sizes from the input connections, subject to any transformations.

            blockOrder = toposort(obj.BlockGraph, 'Order', 'stable');

            errors = MException.empty;
            for bidx = blockOrder
                errors = [errors,  obj.compileBlock(bidx, inputStruct)];
            end

            if ~isempty(errors)
                compileException = MException(message("bioinfo:Pipeline:CompilationError"));

                for e = 1:numel(errors)
                    compileException = addCause(compileException, errors(e));
                end

                obj.Dispatcher.error(compileException, "PipelineCompileFinishWithError", obj.UUID);
            end

            obj.Dispatcher.info(message("bioinfo:Pipeline:CompilePipelineSucceeded"), "PipelineCompileFinish", obj.UUID);
        end

        function run(obj, inputStruct, NameValueArgs)
            % RUN Run pipeline.
            %
            % RUN(P) runs the pipeline P. RUN compiles the pipeline to ensure all
            % blocks are properly configured and input ports are satisfied, and then
            % executes with the default settings.
            %
            % RUN(P, S) runs a pipeline with input struct S. The field
            % names of S must match unconnected ports in P. Use
            % inputNames(P) to retrieve a list of appropriate variable
            % names.
            %
            % RUN(..., 'ResultsDirectory', DIR) specifies where runtime
            % results are stored for each block. If DIR does
            % not exist, the function creates this folder before executing. By 
            % default, DIR is a folder with the name "PipelineResults" in the current
            % directory (PWD).
            %
            % RUN(..., 'DisplayLevel', DL) specifies how much information
            % to print to the command window while a pipeline is running.
            % Valid choices are:
            %
            %  "Off"   (0) - No messages are printed.
            %  "Error" (1) - Only print errors.
            %  "Warn"  (2) - Print warnings and errors (default).
            %  "Info"  (3) - Also print pipeline progress information.
            %  "Debug" (4) - Print more detailed debugging information.
            %
            % DL is a non-negative integer from 0-4, default is 2.
            %
            % RUN(..., 'UseParallel', TF) specifies whether to run in
            % parallel. TF is a logical scalar, default is false.
            %
            % RUN(..., 'RunMode', MODE) specifies the run mode of a
            % Pipeline, either 'Minimal' (default) or 'Full'. When MODE is
            % 'Full', the Pipeline runs all blocks even if they have
            % previously computed results in the results directory. When MODE is
            % 'Minimal', the pipeline runs only the blocks for which one of the following
            % statements is true:
            %
            %   - The block has not been run before
            %   - The block has no saved results in the results directory
            %   - The block object has been modified since the last time it
            %     was run.
            %   - Any input data to the block has changed since the last
            %     run, including new runtime inputs.
            %   - A block upstream of the block has been run more recently
            %     than the last time the block was run.
            %
            % If you specify a subset of blocks to run using "To", "From",
            % or "Only" name-value arguments, the rules above are applied
            % only to those selected blocks. It is recommended that you use
            % the default value ('Minimal') for this name-value argument
            % because skipping up-to-date blocks saves significant time
            % running the pipeline, especially when the pipeline has
            % long-running blocks that do not need to be rerun.
            %
            % RUN(..., 'From', FB) specifies the blocks from which to start
            % the pipeline. If FB is provided, only blocks downstream of FB
            % will be executed, including those in FB. If not provided
            % (default), the pipeline runs every block. This option is
            % incompatible with 'Only'.
            %
            % RUN(..., 'To', TB) specifies the blocks to end the pipeline.
            % If TB is provided, only blocks upstream of TB will be
            % executed, including those in TB. If not provided (default),
            % the pipeline will attempty to run every block. This option is
            % incompatible with 'Only'.
            %
            % RUN(..., 'Only', B) runs only the blocks specified by B. This
            % option is incompatible with 'To' and 'From'.
            %
            % RUN(..., 'SaveResults', B) saves the results of each
            % block specified in the vector B to MAT-files in the results folder. B is a
            % vector of block objects or names, or "-all", which saves results for all
            % blocks in the pipeline. Default is "-all".
            %
            % Example:
            % import bioinfo.pipeline.blocks.FileChooser
            % import bioinfo.pipeline.blocks.SamSort
            %
            % P = bioinfo.pipeline.Pipeline();
            % FCB = FileChooser(which("ex1.sam"));
            % SSB = SamSort;
            %
            % addBlock(P, [FCB, SSB]);
            % connect(P, FCB, SSB, ["Files", "SAMFile"]);
            %
            % run(P);
            % processTable(P)
            %
            % See also inputNames, inputs, processTable, results

            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                inputStruct(1,1) struct = struct()
                NameValueArgs.UseParallel(1,1) logical = false
                NameValueArgs.DisplayLevel(1,1) bioinfo.pipeline.LogLevel = 2
                NameValueArgs.RunMode(1,1) bioinfo.pipeline.RunMode = "Minimal"
                NameValueArgs.ResultsDirectory(1,1) string {mustBeNonzeroLengthText} = fullfile(pwd, 'PipelineResults')
                NameValueArgs.From {bioinfo.pipeline.internal.Validation.mustBeBlockID(NameValueArgs.From)}    = string.empty(1,0)
                NameValueArgs.To {bioinfo.pipeline.internal.Validation.mustBeBlockID(NameValueArgs.To)}        = string.empty(1,0)
                NameValueArgs.Only {bioinfo.pipeline.internal.Validation.mustBeBlockID(NameValueArgs.Only)}    = string.empty(1,0)
                NameValueArgs.SaveResults {bioinfo.pipeline.internal.Validation.mustBeBlockID(NameValueArgs.SaveResults)} = "-all"
            end

            if isRunning(obj)
                error(message("bioinfo:Pipeline:PipelineRunning"));
            end
            % Create new RunContext based on input args
            newRunContext = bioinfo.pipeline.internal.RunContext(NameValueArgs);            

            newRunContext.InputStruct      = inputStruct;

            % If we are switching run context (using a new results directory), or there
            % is no previous run context, attempt to restore run data from the new
            % results directory, or if none is available, create new empty process queues.

            if isempty(obj.RunContext) || ~hasSameResultsDirectory(obj, newRunContext)
               obj.RunContext = newRunContext;
               obj.validateAndCreateResultsDirectory();

               % Switched Run Context, blow away run data and restore from disk if it
               % exists
                obj.RunData = bioinfo.pipeline.internal.RunData;

               obj.restoreRunData(); % creates and populates runData, overwriting existing data
            else
                % If the results directory is the same as previous RunContext,
                % then simply copy over the existing ProcessQueues into new run data.
                obj.RunContext = newRunContext;
                obj.validateAndCreateResultsDirectory();

                newRunData = bioinfo.pipeline.internal.RunData;
                newRunData.ProcessQueues = obj.RunData.ProcessQueues;
                obj.RunData = newRunData;
            end

            % set Logger properties (ProjectDirectory must be created first)

            obj.Logger.ResultsDirectory = obj.RunContext.AbsoluteResultsDirectory;
            obj.Logger.DisplayLevel     = obj.RunContext.DisplayLevel;
            
            logFileName = fullfile(obj.RunContext.AbsoluteResultsDirectory, obj.Logger.LogFileName);
            obj.Dispatcher.info(message("bioinfo:Pipeline:StartingPipelineRun", logFileName), "PipelineRunStart", obj.UUID);

            % Begin compile (adds missing ProcessQueues)
            obj.compile(inputStruct);

            % Select & open parallel pool based on user preferences
            obj.Dispatcher.debug(message("bioinfo:Pipeline:EnsuringParPool"), "PipelineRun", obj.UUID);
            obj.selectParallelPool();

            blocksToRun = obj.getBlocksToRun(NameValueArgs.From, NameValueArgs.To, NameValueArgs.Only);
            
            % Blocks that are stale and are planned to run should have
            % ProcessQueues removed
            if obj.RunContext.RunMode == "Minimal"
                staleBlocks = find(obj.getBlocksStaleness());
                blocksToClear = intersect(blocksToRun, staleBlocks);
            else
                blocksToClear = blocksToRun;
            end

            for bid = blocksToClear(:)'
                blk = obj.BlockGraph.Nodes.Block(bid);
                obj.RunData.ProcessQueues(blk).clear();
            end
            obj.Dispatcher.info(message("bioinfo:Pipeline:BlocksToBeRunResultsCleared"), "PipelineRunResultsCleared", obj.UUID);

            blocksToSaveIdx = obj.findBlockIndex(obj.RunContext.SaveResults);
            isSaveResultsBlockInRun = ismember(blocksToSaveIdx, blocksToRun);
            if ~all(isSaveResultsBlockInRun) && ~isequal(obj.RunContext.SaveResults, "-all")
                blkNamesNotInRun = obj.BlockNames(blocksToSaveIdx(~isSaveResultsBlockInRun));
                msg = message("bioinfo:Pipeline:SavedBlockNotInRun", strjoin(blkNamesNotInRun, '\n\t'));
                obj.Dispatcher.warn(msg, "PipelineRunWarn", obj.UUID);
                warning(msg);
            end

            obj.Dispatcher.info(message("bioinfo:Pipeline:FoundBlocks", numel(blocksToRun)), "PipelineRun", obj.UUID);
            obj.Dispatcher.debug(message("bioinfo:Pipeline:InputStruct", evalc("disp(obj.RunContext.InputStruct)")), "PipelineRun", obj.UUID);
            obj.Dispatcher.debug(message("bioinfo:Pipeline:RunOptions",  strrep(evalc("disp(NameValueArgs)"), "\", "\\")), "PipelineRun", obj.UUID);

            obj.execute(blocksToRun);
        end

        function deleteResults(obj, blockIds, NVPArgs)
            % deleteResults Delete all computed results from pipeline.
            %
            % deleteResults(P) removes all computed results and run
            % metadata from the Pipeline P. P must not be running. To
            % delete results from a running pipeline you must first cancel
            % it with cancel(P).
            %
            % deleteResults(P, B) removes results of the specified blocks
            % in B only. B is a block object, block name, or vector of
            % block objects or names.
            %
            % deleteResults(..., 'INCLUDEFILES', TF) specifies whether to
            % also delete the runtime directories for the blocks in the
            % pipeline. TF is a logical scalar. The default is false.
            %
            % Example:
            % import bioinfo.pipeline.blocks.FileChooser
            % import bioinfo.pipeline.blocks.SamSort
            %
            % P = bioinfo.pipeline.Pipeline();
            % FCB = FileChooser(which("ex1.sam"));
            % SSB = SamSort;
            %
            % addBlock(P, [FCB, SSB]);
            % connect(P, FCB, SSB, ["Files", "SAMFile"]);
            %
            % run(P);
            % deleteResults(P)
            %
            % See also cancel, results, processTable
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                blockIds {mustBeNonempty, bioinfo.pipeline.internal.Validation.mustBeBlockID(blockIds)} = "-all";
                NVPArgs.IncludeFiles(1,1) logical = false;
            end

            % Check if the pipeline is running first. We still need to call
            % cancel() afterward in case parallel processes are started
            % after making this check
            if isRunning(obj)
                error(message("bioinfo:Pipeline:BlocksRunning"));
            end


            blocksToClear = findBlockIndex(obj, blockIds);

            if isempty(obj.RunContext) || isempty(obj.RunData)
                return; % No results to delete, pipeline has never been run
            end

            for bid = blocksToClear(:)'
                blk = obj.BlockGraph.Nodes.Block(bid);

                if obj.RunData.ProcessQueues.isKey(blk)
                    obj.RunData.ProcessQueues(blk).clear();
                end

                runtimeDir = obj.getRuntimeDirectory(bid);

                % Always delete runData.mat, even if IncludeFiles is false
                runDataMatfile = fullfile(runtimeDir, 'runData.mat');
                if isfile(runDataMatfile)
                    delete(runDataMatfile);
                end
                
                if NVPArgs.IncludeFiles && isfolder(runtimeDir)
                    % Warn if directories cannot be removed
                    [status, msg] = rmdir(runtimeDir, 's');
                    if ~status
                        warning(message("bioinfo:Pipeline:RMDIRFailed", runtimeDir, msg));
                    end
                end
            end
        end

        function cancel(obj)
            % CANCEL Cancel running blocks in pipeline.
            %
            % CANCEL(P) Cancels all blocks in the pipeline P that are currently running
            % in a parallel environment and prevents new blocks from being
            % queued.
            %
            % Example:
            % P = bioinfo.pipeline.Pipeline
            % PB = bioinfo.pipeline.blocks.UserFunction(@() pause(10));
            % addBlock(P, PB);
            %
            % run(P, 'UseParallel', true');
            % cancel(P)
            % processTable(P)
            %
            % See also run, processTable, results, fetchResults
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
            end

            % CANCEL works by traversing the graph in topological order,
            % cancelling process queues and propagation futures. This order
            % guarantees that all a block's predecessors will be cancelled
            % (or completed) before it, thereby preventing it from being
            % re-schehuled in the parallel case after it is cancelled.
            %
            % 
            % When the pipeline is run with UseParallel=false, CANCEL is not accessible
            % because the matlab thread is occupied, so a user must use Ctrl-C to interrupt. The
            % biopipelineDesigner does something similar by sending the interrupt
            % signal to the MVM running the pipeline in the serial case.

            if isempty(obj.RunData)
                obj.Dispatcher.info(message("bioinfo:Pipeline:RunCancelled"), "PipelineRunFinish", obj.UUID);
                return
            end

            try
                % In case any processes finish after the cancel is initiated
                % but before the process queue is shut down, set a flag so no
                % new processes will be started.
                obj.RunData.CancelActive = true;

                % Iterate over process queues, calling cancel().
                for n = 1:numel(obj.Blocks)
                    blk = obj.BlockGraph.Nodes.Block(n);
                    if isKey(obj.RunData.ProcessQueues, blk)
                        obj.RunData.ProcessQueues(blk).cancel();
                    end
                end

                % Now cancel any afterAlls that may still be active
                % Use a loop because thread-based futures and pool-based
                % futures cannot be cancelled together, and we may have both
                for pf = 1:numel(obj.RunData.PropagationFutures)
                    cancel(obj.RunData.PropagationFutures(pf));
                end
            catch CE
                if ~strcmp(CE.identifier, 'MATLAB:parallel:future:UnknownWaitErrors')
                    rethrow(CE);
                end
            end
            
            % No futures should be running now, clear flag
            obj.RunData.CancelActive = false;

            obj.RunData.PipelineCompleted = true;
            obj.Dispatcher.info(message("bioinfo:Pipeline:RunCancelled"), "PipelineRunFinish", obj.UUID);

        end

        function wait(obj, blockIds)
            % WAIT Wait for running blocks to complete in pipeline.
            %
            % WAIT(P) blocks execution until all Blocks in the pipeline P
            % have completed, or are guaranteed not to complete due to
            % upstream errors.
            %
            % WAIT(P, B) blocks execution until all Blocks specified by B
            % have completed, or are guaranteed not to complete due to
            % upstream errors. B is a block object, block name, or vector
            % of block objects or names.
            %
            % Example:
            % P = bioinfo.pipeline.Pipeline;
            % PB = bioinfo.pipeline.blocks.UserFunction(@() pause(5));
            % PB2 = bioinfo.pipeline.blocks.UserFunction(@() pause(1));
            % addBlock(P, [PB, PB2]);
            %
            % run(P, 'UseParallel', true');
            % wait(P, PB2);
            % processTable(P)
            %
            % wait(P);
            % processTable(P)
            %
            % See also run, processTable, results
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                blockIds {mustBeNonempty, bioinfo.pipeline.internal.Validation.mustBeBlockID(blockIds)} = "unspecified";
            end

            blockOrder = toposort(obj.BlockGraph, 'Order', 'stable');

            if nargin > 1
                % Filter blockOrder to only the given blockIds, and blocks
                % upstream of blockIds
                k = findBlocksBetween(obj, [], blockIds);
                blockOrder = intersect(blockOrder, k, 'stable');
            end

            if isempty(obj.RunData)
                return; % Pipeline has never been run (and can't currently be running), return early.
            end

            for bidx = blockOrder(:)'
                blk = obj.BlockGraph.Nodes.Block(bidx);
                if isKey(obj.RunData.ProcessQueues, blk)
                    obj.RunData.ProcessQueues(blk).wait();
                    if bidx <= numel(obj.RunData.PropagationFutures)
                        propagationFuture = obj.RunData.PropagationFutures(bidx);
    
                        if ~strcmpi(propagationFuture.State, 'unavailable')
                            wait(propagationFuture);
                        end
                    end
                end
            end

        end

        function inputVars = inputNames(obj, NVPArgs)
            % inputNames Get input port names for pipeline.
            %
            % NAMES = inputNames(P) returns the input port names for all the
            % blocks in the Pipeline P that are not connected. These names
            % can be used as fields of the input struct passed to run(P,
            % S). NAMES is a string vector.
            %
            % inputNames(P, 'IncludeOptional', TF) specifies whether to
            % include optional ports in the returned input port names. The default
            % is false.
            %
            % Examples:
            %
            %   P = bioinfo.pipeline.Pipeline;
            %   P.addBlock(bioinfo.pipeline.blocks.Bowtie2);
            %   inputNames(P)
            %   inputNames(P, 'IncludeOptional', true)
            %
            % See also run, inputs
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                NVPArgs.IncludeOptional(1,1) logical = false
            end

            inputVars = {};
            inputSrcBlk = {};

            % Loop over blocks, creating a struct
            for bidx = 1:size(obj.BlockGraph.Nodes, 1)
                blkName = obj.BlockGraph.Nodes.Name{bidx};
                blk = obj.BlockGraph.Nodes.Block(bidx);

                % Get Input struct for each block
                blkInputs = blk.emptyInputs('IncludeOptional', NVPArgs.IncludeOptional);

                % Filter out connected ports
                blkInEdges = inedges(obj.BlockGraph, bidx);
                connectedPorts = arrayfun(@(eid) obj.BlockGraph.Edges.Connection(eid).PortMap(:,2), blkInEdges, 'UniformOutput', false);
                connectedPorts = vertcat(connectedPorts{:});


                if ~isempty(connectedPorts)
                    fieldsToRemove = connectedPorts(isfield(blkInputs, connectedPorts));
                    blkInputs = rmfield(blkInputs, fieldsToRemove);
                end

                remainingNames = fieldnames(blkInputs)';

                if ~isempty(remainingNames)
                    inputVars = [inputVars, remainingNames];
                    inputSrcBlk = [inputSrcBlk, repmat({blkName}, 1, numel(remainingNames))];
                end
            end

            inputVars = unique(inputVars);

            inputVars = convertCharsToStrings(inputVars);
        end
    end

    methods (Hidden)
        function tf = isRunning(obj)
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
            end
            
            if isempty(obj.RunData)
                tf = false;
                return; % A pipeline with no run data can't be running.
            end

            tf = ~obj.RunData.PipelineCompleted;
        end

        % Hidden Connection APIs
        function connections = findConnection(obj, srcID, tgtID, portMap)
            % FINDCONNECTION is an internal method and is not recommended.

            % C = FINDCONNECTION(P) returns an array C of the connections
            % in the pipeline P.
            %
            % C = FINDCONNECTION(P, SRCID, TGTID) returns only the
            % connections between the blocks identified by SRCID and TGTID.
            % SRCID & TGTID are names of blocks or block objects in P.
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                srcID {bioinfo.pipeline.internal.Validation.mustBeBlockID(srcID,1)} = "unspecified"
                tgtID {bioinfo.pipeline.internal.Validation.mustBeBlockID(tgtID,1)} = "unspecified"
                portMap(:,2) string {mustBeNonempty} = ["", ""];
            end

            connections = bioinfo.pipeline.Connection.empty(0,1);

            if ~isempty(obj.BlockGraph.Edges)
                connections = obj.BlockGraph.Edges.Connection;
            end

            if nargin < 4
                portMap = [];
            end

            if nargin > 2
                connectionIdx = findConnectionIndex(obj, srcID, tgtID, portMap);

                if any(connectionIdx == 0)
                    if isempty(portMap)
                        error(message("bioinfo:Pipeline:NoConnectionsBetweenBlocks"));
                    else
                        error(message("bioinfo:Pipeline:ConnectionsNotFound"));
                    end
                end
                connections = connections(connectionIdx);
            end
        end

        function addConnection(obj, connections)
            % ADDCONNECTION is an internal method and not recommended.
            
            % Needed by the biopipelineDesigner in order to re-use
            % Connection objects.
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                connections(:,1) bioinfo.pipeline.Connection
            end

            try
                % Validate that these are not default-constructed
                % connections, and have source and target blocks, which
                % also exist in the pipeline.
                srcBlocks = [connections.SourceBlock];
                tgtBlocks = [connections.TargetBlock];

                assert(numel(srcBlocks) == numel(connections));
                assert(numel(tgtBlocks) == numel(connections));

                srcIdx = findBlockIndex(obj, srcBlocks);
                tgtIdx = findBlockIndex(obj, tgtBlocks);
            catch
                error(message("bioinfo:Pipeline:InvalidConnections"));
            end

            edgeProps = table(connections, 'VariableNames', "Connection");

            % Check for Cyclic Pipeline
            newGraph = addedge(obj.BlockGraph, srcIdx, tgtIdx, edgeProps);
            if ~isdag(newGraph)
                error(message("bioinfo:Pipeline:CyclicPipeline"));
            end
            obj.BlockGraph = newGraph;

            modelChangedData = bioinfo.pipeline.internal.SemanticModelChangedEventData;
            for c = 1:numel(connections)
                modelChangedData.add(connections(c), obj, "");
            end
            if ~isempty(modelChangedData.Delta)
                notify(obj, 'SemanticModelChanged', modelChangedData);
            end
        end

        function removeConnection(obj, connections)
            % REMOVECONNECTION is an internal method and not recommended.

            % Needed by the biopipelineDesigner in order to re-use
            % Connection objects.
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                connections(:,1) bioinfo.pipeline.Connection
            end

            existingConnections = findConnection(obj);

            [~, toRemove] = intersect(existingConnections, connections);

            if numel(toRemove) ~= numel(connections)
                error(message("bioinfo:Pipeline:ConnectionToRemoveNotFound"));
            end

            % remove connection
            modelChangedData = bioinfo.pipeline.internal.SemanticModelChangedEventData;
            for c = toRemove(:)'
                modelChangedData.remove(obj.BlockGraph.Edges.Connection(c), obj, "");
            end

            obj.BlockGraph = rmedge(obj.BlockGraph, toRemove);

            if ~isempty(modelChangedData.Delta)
                notify(obj, 'SemanticModelChanged', modelChangedData);
            end
        end

        function runStatus = blockRunStatus(obj, blockID)
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                blockID {bioinfo.pipeline.internal.Validation.mustBeBlockID(blockID, 1)}
            end

            k           = findBlockIndex(obj, blockID);
            blk         = obj.BlockGraph.Nodes.Block(k);
            if ~isempty(obj.RunData) && isKey(obj.RunData.ProcessQueues, blk)
                runStatus   = obj.RunData.ProcessQueues(blk).RunStatus;
            else
                runStatus = bioinfo.pipeline.RunStatus.empty;
            end
        end

        function tf = isRunnable(obj, blockIDs)
            % isRunnable helper method for the app to determine if specified
            % blocks can be run with "Run Selected". For this to be true,
            % predecessors blocks must:
            % - have results (previously ran) 
            % - not be stale
            %
            % (NOTE: blocks that have not run are NOT stale)
            arguments
                obj
                blockIDs {mustBeNonempty, bioinfo.pipeline.internal.Validation.mustBeBlockID(blockIDs)}
            end

            blkIdx = findBlockIndex(obj, blockIDs);

            isStale = obj.getBlocksStaleness();
            if isempty(obj.RunData)
                isRunStatusEmpty = true(size(obj.Blocks));
            else
                isRunStatusEmpty = arrayfun(@(b) ~isKey(obj.RunData.ProcessQueues, b) || isempty(obj.RunData.ProcessQueues(b).RunStatus), obj.Blocks);
            end
            isStaleOrRunStatusEmpty = isStale | isRunStatusEmpty;

            % Get unique predecessors, but remove any blocks in the query,
            % as these would eventually be run and therefore not stale
            allPredecessors = arrayfun(@(bidx) predecessors(obj.BlockGraph, bidx), blkIdx, 'UniformOutput', false);
            uniquePredecessorsIdx = setdiff(vertcat(allPredecessors{:}), blkIdx); 
            tf = ~any(isStaleOrRunStatusEmpty(uniquePredecessorsIdx));
        end

        function blockStaleTF = getBlocksStaleness(obj, blockIDs)
            arguments
                obj
                blockIDs(:,1) {bioinfo.pipeline.internal.Validation.mustBeBlockID} = 1:numel(obj.Blocks)
            end
            % Checks which blocks need to be evaluated in the pipeline (for block "stale flag")
            %
            % "Stale" denotes that re-running a block may produce new Results
            %
            % A Block is stale if ALL of the following is true:
            % - pipeline does not have ANY running/queued/waiting blocks
            % - block has results (new blocks are NOT stale)
            % - results folder deleted | upstream results empty | new properties | upstream newer | upstream stale

            blockIDs = obj.findBlockIndex(blockIDs);

            % return false early if pipeline is running
            numBlocks = numel(obj.Blocks);
            blockStaleTF = false(numBlocks, 1);
            if obj.isRunning() || isempty(obj.RunData)
                return
            end

            obj.restoreRunData;

            % Determine which blocks are downstream/upstream from each
            % other. For a given block index in "blkIsConnectedMatrix":
            % - Rows returns a logical vector of downstream blocks
            % - Columns returns a logical vector of upstream blocks
            %
            % upstream does not include the current block
            blkUnweightedDistMatrix = distances(obj.BlockGraph, 'Method', 'unweighted');
            blkIsConnectedMatrix = ~isinf(blkUnweightedDistMatrix) & blkUnweightedDistMatrix ~= 0;


            hasRunStatus = arrayfun(@(b) isKey(obj.RunData.ProcessQueues, b) && ~isempty(obj.RunData.ProcessQueues(b).RunStatus), obj.Blocks);

            blockOrder = toposort(obj.BlockGraph, "Order", "stable");

            for iB = 1:numBlocks
                upstreamBlkStaleTF = blkIsConnectedMatrix(:, blockOrder(iB));
                isUpstreamMissingRunStatus = any(~hasRunStatus(blkIsConnectedMatrix(:, blockOrder(iB))));
                blk = obj.BlockGraph.Nodes.Block(blockOrder(iB));

                if isKey(obj.RunData.ProcessQueues, blk)
                    pq = obj.RunData.ProcessQueues(blk);
                    hasError = any(pq.RunStatus == bioinfo.pipeline.RunStatus.Error, "all");
                    hasMergeFailure = pq.hasMergeFailure();
                    hasRunFolder = ~isempty(obj.RunContext) && isfolder(obj.getRuntimeDirectory(blockOrder(iB)));
                else
                    hasError = false;
                    hasMergeFailure = false;
                    hasRunFolder = false;
                end

                
                % NOTE: upstreamRunsNewer will error if the blockRunFolder does not exist
                blockStaleTF(blockOrder(iB)) = hasRunStatus(blockOrder(iB)) &&...
                    (any(blockStaleTF(upstreamBlkStaleTF)) || isUpstreamMissingRunStatus || hasError || hasMergeFailure || obj.isLastRunPropertiesStale(blockOrder(iB)) || ~hasRunFolder || obj.upstreamRunsNewer(blockOrder(iB)));

                blockStaleTF = blockStaleTF(blockIDs);
            end
        end

        function tf = isLastRunPropertiesStale(obj, bidx)
            % Check if the block, the input data, or the connections have
            % changed for a block since the last run.
            blk = obj.BlockGraph.Nodes.Block(bidx);
            pq = obj.RunData.ProcessQueues(blk);

            tf = pq.isStale() || ~isequaln(obj.getBlockConnectionDict(bidx), pq.LastRunProperties.ConnectionDict);
        end

        function connnDict = getBlockConnectionDict(obj, bidx)
            % return a dictionary where the block's Input UUIDs are keys
            % and upstream Output UUIDs are values for each connection. If
            % the Input is not connected, it is not included in the
            % dictionary.

            connnDict = dictionary(string.empty, string.empty);

            inConnIdxs = inedges(obj.BlockGraph, bidx);
            for connIdx = inConnIdxs'
                connection = obj.BlockGraph.Edges.Connection(connIdx);

                % Source
                srcPortName = connection.PortMap(1);
                srcBlk = connection.SourceBlock;
                srcPortUUID = srcBlk.Outputs.(srcPortName).UUID;

                % Target
                tgtPortName = connection.PortMap(2);
                tgtBlk = connection.TargetBlock;
                tgtPortUUID = tgtBlk.Inputs.(tgtPortName).UUID;

                % update dictionary
                connnDict(tgtPortUUID) = srcPortUUID;
            end
        end

        function ps = getParallelSize(obj, blockID)
            % Get the size of the blocks parallelization on the most recent
            % run
            arguments
                obj
                blockID {mustBeNonempty, bioinfo.pipeline.internal.Validation.mustBeBlockID(blockID, 1)}
            end

            k = obj.findBlockIndex(blockID);
            blk = obj.BlockGraph.Nodes.Block(k);
            pq = obj.RunData.ProcessQueues(blk);
            ps = pq.getParallelSize();
        end
    end

    methods(Access = private)
        function p = selectParallelPool(obj)
            % Select appropriate parallel pool.
            p = parallel.Pool.empty;
            useParallel = obj.RunContext.UseParallel;
            if useParallel
                if matlab.internal.parallel.isPCTInstalled && matlab.internal.parallel.isPCTLicensed
                    % PCT is available so make sure pool is open
                    p = gcp();
                    if isempty(p)
                        p = parpool; % force open pool
                    elseif isa(p, 'parallel.ThreadPool')
                        obj.Dispatcher.error(message("bioinfo:Pipeline:ThreadsNotSupported"), "PipelineRunError", obj.UUID);
                    end
                else
                    obj.Dispatcher.warn(message("bioinfo:Pipeline:PCTUnavailable"), "PipelineRunWarn", obj.UUID);
                end
            end

            obj.RunData.ParallelPool = p;
        end

    end

    methods(Access = ?bioinfo.pipeline.internal.designer.PipelineController)
        function runtimeDir = getRuntimeDirectory(obj, blockID)
            % gets run-time directory of block from last run (RunContext).
            % If pipeline has not run yet, returns empty. If no input args
            % are given, returns RunContext.AbsoluteResultsDirectory
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                blockID = []
            end

            if nargin > 1
                bioinfo.pipeline.internal.Validation.mustBeBlockID(blockID, 1);
            end

            runDir = obj.RunContext.AbsoluteResultsDirectory;

            if isempty(runDir)
                runtimeDir = [];
            else
                if nargin > 1
                    k           = findBlockIndex(obj, blockID);
                    blkName     = obj.BlockGraph.Nodes.Name{k};
                    runtimeDir  = fullfile(runDir, blkName);
                else
                    runtimeDir = runDir;
                end
            end
        end
    end

    methods (Access = private)
        function createMissingProcessQueues(obj)
            % createMissingProcessQueues create a new process queue for every block in
            % the pipeline that isn't currently in RunData
            if ~isempty(obj.RunData)
                blks = obj.Blocks;
                blkNames = obj.BlockNames;
                for bidx = 1:numel(blks)
                    if ~isKey(obj.RunData.ProcessQueues, blks(bidx))
                        obj.RunData.ProcessQueues(blks(bidx)) = bioinfo.pipeline.internal.ProcessQueue(blks(bidx), blkNames(bidx));
                    end
                end
            end
        end

        function saveRunData(obj, bidx)
            % Save RunData for the given block into the block's runtime directory. This
            % runData file is used to restore run data if the Pipeline object is
            % reloaded from a RunContext
            if isequal(obj.RunContext.SaveResults, "-all") || ismember(bidx, obj.findBlockIndex(obj.RunContext.SaveResults))
                blk = obj.BlockGraph.Nodes.Block(bidx);
                blockRunFolder = obj.getRuntimeDirectory(bidx);
                runDataFile = fullfile(blockRunFolder, 'runData.mat');
                ProcessQueue = obj.RunData.ProcessQueues(blk);
    
                % Force ProcessQueue to update before saving
                refresh(ProcessQueue);
                warnState = warning("off", "MATLAB:Figure:FigureSavedToMATFile");
                warnCleanup = onCleanup(@() warning(warnState));
                save(runDataFile,"ProcessQueue"); % Store ProcessQueue
            end
        end

        function restoreRunData(obj, blockIDs)
            % Restore run data from the filesystem. This is needed whenever the
            % pipeline is run, or results are requested for a RunContext whose run data
            % might no longer be in memory (for example, if you load a pipeline from a
            % MAT-file).
            %
            % If a pipeline has no run context, this is a no-op.
            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                blockIDs {bioinfo.pipeline.internal.Validation.mustBeBlockID} = 1:numel(obj.Blocks)
            end

            if isempty(obj.RunData)
                obj.RunData = bioinfo.pipeline.internal.RunData;
            end

            if isempty(obj.RunContext)
                % No way to restore from missing run context, but also skip if RunData
                % already exists unless forceRestore is true.
                return; % Can't restore any data without a valid run context
            end

            blockIDs = obj.findBlockIndex(blockIDs);

            if isfolder(obj.RunContext.AbsoluteResultsDirectory)

                for bidx = blockIDs(:)'
                    % For each block, check if runData exists. If it does, load the saved
                    % Process Queue
                    blockRunFolder = obj.getRuntimeDirectory(bidx);
                    runDataFile = fullfile(blockRunFolder, 'runData.mat');
                    blk = obj.Blocks(bidx);

                    % Only restore run data if the process queue for this block does not exist
                    % yet
                    if ~isKey(obj.RunData.ProcessQueues, blk) && isfile(runDataFile)
                        % Restore from existing, but only if the saved ProcessQueue's Block has the
                        % same UUID as the current block.
                        try
                            runData = load(runDataFile);
                            if runData.ProcessQueue.Block.UUID == obj.Blocks(bidx).UUID
                                runData.ProcessQueue.Block = obj.Blocks(bidx); % Overwrite so ProcessQueue uses the current block in the pipeline.
                                obj.RunData.ProcessQueues(obj.Blocks(bidx)) = runData.ProcessQueue;
                            end
                        catch
                            continue % Ignore failures
                        end
                    end
                end
            end
        end

        function tf = hasSameResultsDirectory(obj, newRunContext)
            % returns true if the current run context points to the same location as
            % the newRunContext

            tf = ~isempty(obj.RunContext) || ~isfolder(newRunContext.ResultsDirectory);

            if tf
                currentResultsDirectory = obj.RunContext.AbsoluteResultsDirectory;

                newResultsDirectory = what(newRunContext.ResultsDirectory);
                tf = ~isempty(newResultsDirectory) && isequal(newResultsDirectory(1).path, currentResultsDirectory);
            end
        end

        function handleBlockModelChanged(obj, blkChangedData)
            % Handler function for ModelChanged events coming from blocks.
            % When a block's ports are changed, the Pipeline must update
            % the connections, and then bundle all semantic updates into a
            % new event to notify the PipelineController
            
            blkIdx = findBlockIndex(obj, blkChangedData.Source);

            % Pipeline delta appends any connection changes to the Block
            % model delta
            blkDelta = blkChangedData.Delta;
            pipelineModelChangedData = bioinfo.pipeline.internal.SemanticModelChangedEventData();
            pipelineModelChangedData.Delta = blkDelta;

            inConnIdxs = inedges(obj.BlockGraph, blkIdx);
            outConnIdxs = outedges(obj.BlockGraph, blkIdx);

            % First apply all deletions, in case there is overlap between
            % deleted and renamed port names
            isInput = arrayfun(@(port) isa(port, 'bioinfo.pipeline.Input'), blkDelta.Model);
            
            inputDeletions = blkDelta(isInput & blkDelta.Action == "Removed", :);
            outputDeletions = blkDelta(~isInput & blkDelta.Action == "Removed", :);

            inputConnections = obj.BlockGraph.Edges.Connection(inConnIdxs);
            outputConnections = obj.BlockGraph.Edges.Connection(outConnIdxs);

            % Concatenate portmaps with empty string array to ensure
            % correct comparison in the empty case.
            inputPortMaps = vertcat(string.empty(0,2), inputConnections.PortMap);
            outputPortMaps = vertcat(string.empty(0,2), outputConnections.PortMap);

            % Compare target column of inputPort maps against any deleted
            % port names to find whether a connection should be removed.
            % Apply the same check to the source column of output
            % connections
            inConnsToRemove = any(inputPortMaps(:,2) == inputDeletions.Name(:,1)', 2);
            outConnsToRemove = any(outputPortMaps(:,1) == outputDeletions.Name(:,1)', 2);

            allConnectionsToRemove = [inConnIdxs(inConnsToRemove); outConnIdxs(outConnsToRemove)];

            for rmConnIdx = allConnectionsToRemove(:)'
                 % Add to ModelChangedDelta
                 conn = obj.BlockGraph.Edges.Connection(rmConnIdx);
                 pipelineModelChangedData.remove(conn, obj,""); 
            end

            obj.BlockGraph = rmedge(obj.BlockGraph, allConnectionsToRemove);

            % Connection indices may have changed, so regenerate them.

            inConnIdxs = inedges(obj.BlockGraph, blkIdx);
            outConnIdxs = outedges(obj.BlockGraph, blkIdx);

            % For ports that have been renamed, update the portMaps
%             renameDelta = blkDelta(blkDelta.Action == "Renamed",:);
            renamedInputs = blkDelta(isInput & blkDelta.Action == "Renamed", :);
            renamedOutputs = blkDelta(~isInput & blkDelta.Action == "Renamed", :);


            % Update PortMap for connections to inputs that have been
            % renamed
            for c = 1:numel(inConnIdxs)
                conn = obj.BlockGraph.Edges.Connection(inConnIdxs(c));
                renameIdx = find(conn.PortMap(1,2) == renamedInputs.Name(:,1));
                if ~isempty(renameIdx)
                    conn.PortMap(1,2) = renamedInputs.Name(renameIdx, 2);
                end
            end

            % Update PortMap for connections from outputs that have been
            % renamed
            for c = 1:numel(outConnIdxs)
                conn = obj.BlockGraph.Edges.Connection(outConnIdxs(c));
                renameIdx = find(conn.PortMap(1,1) == renamedOutputs.Name(:,1));
                if ~isempty(renameIdx)
                    conn.PortMap(1,1) = renamedOutputs.Name(renameIdx, 2);
                end
            end

            % Notify both block changes and connection changes for any
            % listening apps

            notify(obj, 'SemanticModelChanged', pipelineModelChangedData);
        end

        function errors = compileBlock(obj, bidx, inputStruct)
            %COMPILEBLOCK compile a single block, returning any errors
            %
            % This ensures that:
            %   - All required inputs are connected, or are passed as
            %   pipeline inputs
            %   - There are no conflicting connections (two input tables
            %   that share values used by this block)
            %   - This block's process queue is initialized

            blkName = obj.BlockGraph.Nodes.Name{bidx};
            blkUUID = obj.BlockGraph.Nodes.Block(bidx).UUID;

            obj.Dispatcher.info(message("bioinfo:Pipeline:CompilingBlock", blkName), "BlockCompileStart", blkUUID);

            %Init runtime directory, but only if RunData exists
            if ~isempty(obj.RunContext) && ~isempty(obj.RunData)
                runtimeDir = obj.getRuntimeDirectory(bidx);
                if exist(runtimeDir, 'dir') ~= 7
                    mkdir(runtimeDir);
                end
            end

            errors = MException.empty;

            % Block-specific compilation checks
            try
                compile(obj.findBlock(bidx));
            catch compileEx
                obj.Dispatcher.error(compileEx, "BlockCompileError", blkUUID, "Throw", false);
                errors = compileEx;
            end

            [satisfiedPorts, inputErrors] = obj.checkBlockInputs(bidx, inputStruct);

            errors = [errors, inputErrors];

            if ~isempty(obj.RunContext) && ~isempty(obj.RunData)
                % Only compile ProcessQueue if we have run data.
                blk = obj.BlockGraph.Nodes.Block(bidx);
                pq = obj.RunData.ProcessQueues(blk);
                pq.compile(satisfiedPorts);
            end

            if numel(errors) == 0
                obj.Dispatcher.info(message("bioinfo:Pipeline:CompileBlockSucceeded", blkName), "BlockCompileFinish", blkUUID);
            end
        end

        function [satisfiedPorts, errors] = checkBlockInputs(obj, bidx, inputStruct)
            % CHECKBLOCKINPUTS check that the inputs to this block are
            % valid. Also return the names of inputs that are satisfied by
            % incoming connections or by having an explicit value set on
            % the block
            %
            % Errors:
            % - Required block inputs are not connected
            % - Incoming connections have conflicting output names
            errors = MException.empty;

            blk = obj.BlockGraph.Nodes.Block(bidx);
            blkName = obj.BlockGraph.Nodes.Name{bidx};

            requiredPorts = blk.getRequiredInputs();

            % Ports are either satisfied by incoming connections, a set
            % value of the Input object, or the input struct.

            inPortNames = connectedInputs(obj, bidx);

            % Also check for blocks inputs with a set value
            blkInputNames = string(fieldnames(blk.Inputs));
            hasValue = arrayfun(@(inpName)  ~isequal(blk.Inputs.(inpName).Value, bioinfo.pipeline.datatypes.Unset), blkInputNames);
            blkInputsWithValue = blkInputNames(hasValue)';

            % Check the input struct

            % Get Inputs from pipeline, but filter out any
            pipelineInputVars = convertCharsToStrings(fieldnames(inputStruct));
            % Filter pipeline inputs to those valid for this block
            pipelineInputVars = intersect(pipelineInputVars, blkInputNames);

            satisfiedPorts = unique([inPortNames(:); blkInputsWithValue(:); pipelineInputVars(:)]');

            allSatisfied = all(ismember(requiredPorts, satisfiedPorts));

            if ~allSatisfied
                unSatisfiedPorts = setdiff(requiredPorts, satisfiedPorts);
                unSatisfiedPortUUIDs = arrayfun(@(x) blk.Inputs.(x).UUID, unSatisfiedPorts);
                newError = obj.Dispatcher.error(message("bioinfo:Pipeline:RequiredPortNotConnected", blkName), "InputPortCompileError", unSatisfiedPortUUIDs, "Throw", false);
                errors = [errors, newError];
            end
        end

        function obj = execute(obj, blocksToRun)

            % When running a pipeline, keep track of whether a block has
            % completed, but also whether the block has been "evaluated",
            % meaning the processes have been queued but may have not been
            % completed. This allows the Pipeline to avoid scheduling the
            % same block twice, and when to skip blocks whose
            % predecessors will never run.
            %
            % Set Evaluated & Completed to false only for blocks that will
            % be run.
            blockOrder = toposort(obj.BlockGraph, "Order", "stable");

            blocksToRun = intersect(blockOrder, blocksToRun, "stable");
            obj.RunData.BlocksToRun = blocksToRun;

            obj.RunData.BlocksEvaluated = true(size(obj.Blocks));
            obj.RunData.BlocksCompleted = true(size(obj.Blocks));
            obj.RunData.BlocksEvaluated(blocksToRun) = false;
            obj.RunData.BlocksCompleted(blocksToRun) = false;

            obj.insertRuntimeInputs();

            % Filter to start from only the source blocks in the subgraph
            % of blocks to run
            sourceBlocks = [];
            if ~isempty(blocksToRun)
                runSubgraph = subgraph(obj.BlockGraph, blocksToRun);
                blocksToRunNames = obj.BlockGraph.Nodes.Name(blocksToRun); % Node indices are different in subgraph, but Names are identical
                sourceBlocks = blocksToRun(cellfun(@(b) isempty(predecessors(runSubgraph, b)), blocksToRunNames));
            end

            % Need a DataQueue to check for pipeline completion events to
            % ensure dispatch only happens once
            obj.RunData.CompletionQueue = parallel.pool.DataQueue;
            afterEach(obj.RunData.CompletionQueue, @(~) obj.checkPipelineCompletion());

            % Do this in onCleanup in case run is terminated by Ctrl-C or another error
            interruptCleanup = bioinfo.pipeline.internal.Cleanup(@() obj.handleRunInterrupt());
            obj.RunData.PipelineCompleted = false;

            for bidx = sourceBlocks(:)'
                obj.checkExecuteBlock(bidx);
            end

            % Ensure we check pipeline completion at least once
            send(obj.RunData.CompletionQueue, [])

            % Because of the implementation of in-process parfeval, in
            % order for series-mode to work correctly we need to query the
            % state of one of the futures. This will happen automatically
            % from a matlab prompt when it becomes idle, but in test code,
            % immediatley requesting data from the pipeline like the
            % processTable does not look at the futures directly, so we
            % should wait() here if the ParallelPool is empty.
            if isempty(obj.RunData.ParallelPool)
                wait(obj);
            end
            
            % If we made it here, disable the interruptHandler
            removeTask(interruptCleanup);
        end

        function handleRunInterrupt(obj)
            obj.RunData.PipelineCompleted = true;
            obj.Dispatcher.info(message("bioinfo:Pipeline:RunCancelled"), "PipelineRunFinish", obj.UUID);
        end

        function checkPipelineCompletion(obj)
            % Check if the pipeline has completed all the blocksToRun, or
            % those blocks have been guaranteed to not run.
            if obj.RunData.PipelineCompleted
                return; % Don't notify twice
            end

            % Pipeline completion tracked on a block level
            if all(obj.RunData.BlocksCompleted)
                obj.RunData.PipelineCompleted = true;
                obj.Dispatcher.info(message("bioinfo:Pipeline:RunFinished"), "PipelineRunFinish", obj.UUID);
            end
        end

        function checkExecuteBlock(obj, bidx)
            % Checks whether the given block is ready to be evaluated in a
            % pipeline. If all upstream blocks are completed, and this
            % block has not been evaluated yet, it will be evaluated here
            % once.
            %
            % If the block has all it's inputs satisfied (propagated into
            % it's process queue), it is ready to run, and will run if none
            % of the Make-like behavior causes it to be skipped.
            %
            % If the inputs are not satisfied, the block is guaranteed
            % never to be ready to run, so mark the block as "completed",
            % and skip the block.
            %
            % If the block is evaluated in this call, then also propagate
            % the outputs (if any) to downstream blocks. And call
            % checkExecute on this block's successors. This way, block
            % completion drives progress through the graph.

            blk = obj.BlockGraph.Nodes.Block(bidx);
            pq = obj.RunData.ProcessQueues(blk);
            blockEvaluated = obj.RunData.BlocksEvaluated(bidx);
            upstreamBlockIds = predecessors(obj.BlockGraph, bidx);

            % Block can only run if all upstream blocks are complete and it
            % hasn't been evaluated
            if ~blockEvaluated && all(obj.RunData.BlocksCompleted(upstreamBlockIds))

                % Inputs not satisfied, but not waiting for anymore inputs.
                % Indicate that this block is skipped, and mark completed
                if ~obj.upstreamBlocksComplete(bidx)
                    msg = message("bioinfo:Pipeline:SkippingBlockWontRun", obj.BlockGraph.Nodes.Name{bidx});
                    obj.Dispatcher.info(msg, "BlockRunSkip", obj.BlockGraph.Nodes.Block(bidx).UUID);
                else
                    % Run the block if it hasn't been run yet (Status is empty)
                    % or runtime settings say to re-run blocks
                    shouldRun = true;

                    if obj.RunContext.RunMode == bioinfo.pipeline.RunMode.Full
                        runRationale = getString(message("bioinfo:Pipeline:RunRationaleFullRun"));
                    elseif obj.isLastRunPropertiesStale(bidx)
                        runRationale = getString(message("bioinfo:Pipeline:RunRationaleInputsChanged"));
                    elseif pq.hasMergeFailure()
                        runRationale = getString(message("bioinfo:Pipeline:RunRationaleMergeFailed"));
                    elseif obj.upstreamRunsNewer(bidx)
                        runRationale = getString(message("bioinfo:Pipeline:RunRationaleUpstreamNewer"));
                    elseif any(pq.RunStatus == bioinfo.pipeline.RunStatus.Error, "all")
                        runRationale = getString(message("bioinfo:Pipeline:RunRationaleBlockFailed"));
                    else
                        shouldRun = false;
                        obj.Dispatcher.info(message("bioinfo:Pipeline:SkippingBlockRun", obj.BlockGraph.Nodes.Name{bidx}), "BlockRunSkip", obj.BlockGraph.Nodes.Block(bidx).UUID);
                    end

                    if shouldRun
                        obj.Dispatcher.info(message("bioinfo:Pipeline:BlockQueuedToRun", obj.BlockGraph.Nodes.Name{bidx}, runRationale), "BlockRun", obj.BlockGraph.Nodes.Block(bidx).UUID);

                        obj.collectInputs(bidx);
                        blockDir = obj.getRuntimeDirectory(bidx);
                        saveBlkResults = ismember(obj.findBlockIndex(blk), obj.findBlockIndex(obj.RunContext.SaveResults));
                        connDict = obj.getBlockConnectionDict(bidx);

                        pq.run(obj.RunData.ParallelPool, obj.Dispatcher, blockDir, saveBlkResults, connDict);
                    end
                end

                % Don't evaluate the block again
                obj.RunData.BlocksEvaluated(bidx) = true;

                % Propagate results
                obj.RunData.PropagationFutures(bidx) = afterAll(pq, @(~) runDownstream(obj, bidx), obj.Dispatcher);
                if isempty(obj.RunData.ParallelPool)
                    wait(obj.RunData.PropagationFutures(bidx));
                end
            end
        end

        function tf = upstreamBlocksComplete(obj, bidx)
            % Check if upstream blocks are complete

            upstreamBids = predecessors(obj.BlockGraph, bidx);
            if isempty(upstreamBids)
                tf = true;
            else
                upstreamBlks = obj.BlockGraph.Nodes.Block(upstreamBids);
                upstreamPQs = obj.RunData.ProcessQueues(upstreamBlks);

                tf = all(arrayfun(@(x) x.isComplete(), upstreamPQs));
            end
        end

        function collectInputs(obj, referenceBidx)
            % Collect inputs for block (results from upstream blocks)
            % Save runData for the block (ProcessQueue) to runData.mat

            referenceBlk = obj.BlockGraph.Nodes.Block(referenceBidx);
            referencePQ = obj.RunData.ProcessQueues(referenceBlk);
            [inConnIds, upstreamBids] = inedges(obj.BlockGraph, referenceBidx);
            referenceBlockName = obj.BlockGraph.Nodes.Name{referenceBidx};

            msg = message("bioinfo:Pipeline:CollectingInputs", referenceBlockName);
            obj.Dispatcher.info(msg, "BlockRun", referenceBlk.UUID);

            for c = 1:numel(inConnIds)
                % Move data from propagateBidx to downstream bidx
                connIdx = inConnIds(c);
                upstreamBidx = upstreamBids(c);

                connPortMap = obj.BlockGraph.Edges.Connection(connIdx).PortMap;
                upstreamBlk = obj.BlockGraph.Nodes.Block(upstreamBidx);
                upstreamPQ = obj.RunData.ProcessQueues(upstreamBlk);

                for v = 1:size(connPortMap, 1)
                    outputValue = upstreamPQ.OutputData.value(connPortMap(v,1));
                    referencePQ.insertInput(connPortMap(v,2), outputValue);
                end
            end
        end

        function runDownstream(obj, propagateBidx)
            % Called by execute when a block completes to evaluate
            % downstream blocks. Reports to the pipeline CompletionQueue to
            % check for completion.

            obj.RunData.BlocksCompleted(propagateBidx) = true;

            obj.saveRunData(propagateBidx);

            if obj.RunData.CancelActive
                return;
            end

            % Only do this once per downstream block
            downstreamBids = successors(obj.BlockGraph, propagateBidx);

            for downstreamBidx = downstreamBids(:)'
                % Check if downstream blocks can run and propagate
                % their results.
                % Skip downstream blocks that shouldn't be run
                if ~ismember(downstreamBidx, obj.RunData.BlocksToRun)
                    continue
                end

                obj.checkExecuteBlock(downstreamBidx);
            end

            send(obj.RunData.CompletionQueue, []); % Send to completion queue to check if pipeline is completed
        end

        function tf = upstreamRunsNewer(obj, bidx)
            % Returns true if any blocks upstream of the given block have
            % run more recently, as determined by the presence of timestamp
            % files in their block run directory. See
            % ProcessQueue.writeTimestamp, ProcessQueue.readTimestamp.
            %
            % If there are no upstream blocks, this function returns false.

            import bioinfo.pipeline.internal.ProcessQueue.readTimestamp

            blockRunFolder = obj.getRuntimeDirectory(bidx);

            blockLastRunTime = readTimestamp(blockRunFolder);

            upstreamBids = predecessors(obj.BlockGraph, bidx);

            tf = isnat(blockLastRunTime);

            if ~isempty(upstreamBids)
                upstreamBlockLastRunTimes = arrayfun(@(bid) readTimestamp(getRuntimeDirectory(obj, bid)), upstreamBids);
                tf = tf ||  any(isnat(upstreamBlockLastRunTimes)) || any(upstreamBlockLastRunTimes > blockLastRunTime);
            end

        end

        function inputStruct = getRuntimeInputs(obj, bidx)
            % Retrieve the relevant pipeline inputs for the given block
            % from the RunContext struct. The input table will include
            % globals (unqualified inputs), and inputs that are scoped to
            % this block name

            runtimeInputNames = fieldnames(obj.RunContext.InputStruct);

            % Filter table to valid inputs for this block.
            validVarNames = fieldnames(obj.BlockGraph.Nodes.Block(bidx).Inputs);
            fieldsToRemove = setdiff(runtimeInputNames, validVarNames);

            inputStruct = rmfield(obj.RunContext.InputStruct, fieldsToRemove);
        end

        function insertRuntimeInputs(obj)
            % Before any blocks are run, this function inserts runtime
            % inputs into the relevant process queues. This negates the
            % need for a collectInputs step in addition to a propagation
            % step, as source blocks will not need to collectInputs and
            % downstream blocks will be notified when to check and run

            runtimeInputNames = fieldnames(obj.RunContext.InputStruct);

            unusedInputNames = runtimeInputNames;

            for bidx = 1:size(obj.BlockGraph.Nodes, 1)
                blk = obj.BlockGraph.Nodes.Block(bidx);
                pqueue = obj.RunData.ProcessQueues(blk);

                % Insert any runtime inputs that match the names of this
                % block's inputs, unless they are connected within the
                % pipeline.
                matchingVarNames = intersect(fieldnames(blk.Inputs), runtimeInputNames);
                connInputs = obj.connectedInputs(bidx);
                matchingVarNames = setdiff(matchingVarNames, connInputs);

                for mv = 1:numel(matchingVarNames)
                    pqueue.insertInput(matchingVarNames(mv), obj.RunContext.InputStruct.(matchingVarNames(mv)));
                end

                unusedInputNames = setdiff(unusedInputNames, matchingVarNames);

                % Then insert any inputs that are set via Port value,
                % unless they are connected within the pipeline.
                block = obj.BlockGraph.Nodes.Block(bidx);
                blkInputs = string(fieldnames(block.Inputs));
                inputsWithValues = arrayfun(@(inpName) ~isequal(block.Inputs.(inpName).Value, bioinfo.pipeline.datatypes.Unset), blkInputs);
                inputsWithValues = blkInputs(inputsWithValues);
                inputsWithValues = setdiff(inputsWithValues, connInputs);

                for inp = 1:numel(inputsWithValues)
                    pqueue.insertInput(inputsWithValues(inp), block.Inputs.(inputsWithValues(inp)).Value);
                end

            end

            if ~isempty(unusedInputNames)
                msg = message('bioinfo:Pipeline:UnusedInputFields', strjoin(unusedInputNames, ", "));
                obj.Dispatcher.warn(msg, "PipelineRunWarn", obj.UUID);
                warning(msg);
            end

        end

        function names = assignBlockName(obj, blocks)
            %ASSIGNBLOCKNAME assigns a unique block name to BLOCKS in the
            %context of this Pipeline.
            %
            % N = ASSIGNBLOCKNAME(P, B) assigns default names N to the
            % BLOCKs B in the pipeline P. The default is "<BlockType>_<N>"
            % where <BlockType> is class(B) and <N> is one more than the
            % largest integer that is already in the block graph.

            arguments
                obj(1,1) bioinfo.pipeline.Pipeline
                blocks(:,1) bioinfo.pipeline.Block % Vector of Blocks
            end

            % Build dictionary mapping block class name to largest _<N>
            blkNameCount = dictionary(string.empty, double.empty);
            blkNames = obj.BlockNames;

            for b = 1:numel(blkNames)
                match = regexp(blkNames(b), '^(?<type>.+)_(?<N>[0-9]+)$', 'names');
                % For each block name that matches the pattern store the
                % largest N

                if ~isempty(match)
                    if ~isKey(blkNameCount, match.type)
                        blkNameCount(match.type) = str2double(match.N);
                    else
                        % Keep the largest value
                        blkNameCount(match.type) = max(blkNameCount(match.type), str2double(match.N));
                    end
                end
            end

            % Now assign block names to each new block
            names = strings(numel(blocks), 1);

            % Now for each new block, use the map to determine the next N
            for nb = 1:numel(blocks)

                blkType = class(blocks(nb));

                % Strip packages from class name
                periods = strfind(blkType, '.');
                if ~isempty(periods)
                    blkType = blkType(periods(end)+1:end);
                end

                if isKey(blkNameCount, blkType)
                    n = blkNameCount(blkType) + 1;
                else
                    n = 1;
                end
                names(nb,1) = sprintf("%s_%d", blkType, n);
                blkNameCount(blkType) = n;
            end
        end

        function k = findBlockIndex(obj, blockID)
            % Locate block in BlockGraph
            %
            % K = FINDBLOCK(P, blockID) returns the numeric indices K for
            % each block specified by blockID in P.  blockID is either a
            % numeric index, cell array of character vectors or string
            % array of block names.
            arguments
                obj
                blockID
            end

            % This check should ideally belong in the customization of the load
            % process, but for V1 it can live here. If this pipeline was loaded from a
            % MAT file, the BlockSemanticModelChangeListeners will have been lost.
            % If a Block is modified without the model change listener, then the
            % pipeline's connections can be corrupted. Now, the only way to get a
            % handle to the Block object to modify it is to use findBlock, which passes
            % through this method and restores the model change listeners.
            %
            % NOTE: If the user saves a reference to the same block object elsewhere in
            % the MAT file, and uses that reference to modify the block without
            % triggering this check, then the Pipeline can be corrupted, although the
            % likelihood of a user doing that seems low.
            if isempty(obj.BlockSemanticModelChangeListeners) && ~isempty(obj.BlockGraph.Nodes)
                obj.BlockSemanticModelChangeListeners = arrayfun(@(blk) listener(blk, 'SemanticModelChanged', @(src, evt) obj.handleBlockModelChanged(evt)), obj.BlockGraph.Nodes.Block);
            end

            if isempty(blockID)
                k = [];
                return;
            end

            if isempty(obj.BlockGraph.Nodes) && ~isequal(blockID, "-all")
                throwAsCaller(MException(message("bioinfo:Pipeline:BlockNotFound")));
            end

            if isa(blockID, 'bioinfo.pipeline.Block')
                k = zeros(size(blockID));

                if size(obj.BlockGraph.Nodes, 1) > 0
                    for bid = 1:numel(blockID)
                        kb = find(blockID(bid) == obj.BlockGraph.Nodes.Block);

                        if ~isempty(kb)
                            k(bid) = kb;
                        end
                    end
                end
            elseif isequal(blockID, "-all")
                k = 1:numnodes(obj.BlockGraph);
            else
                % Pass-through for normal graph
                k = findnode(obj.BlockGraph, blockID);
            end

            % Error if block is not found
            if any(k == 0, 'all')
                throwAsCaller(MException(message("bioinfo:Pipeline:BlockNotFound")));
            end
        end

        function k = findConnectionIndex(obj, sourceBlockID, targetBlockID, portMap)
            % Internal method to locate connections in BlockGraph
            %
            % K = FINDCONNECTIONINDEX(P, S, T) returns the numeric index K
            % for each connection between the blockID pairs specified by S
            % and T. S & T are names of blocks in P, or arrays of block
            % objects.
            %
            % K = FINDCONNECTIONINDEX(P, S, T, PM) restricts the
            % connections to those matching the rows in the port map PM. If
            % PM is provided, K will have the same number of rows as PM,
            % and a 0 indicates rows that do not have a corresponding
            % connection.

            % A connection index K == 0 indicates edges that are not found
            % in the graph

            sourceK = findBlockIndex(obj, sourceBlockID);
            targetK = findBlockIndex(obj, targetBlockID);
            k = findedge(obj.BlockGraph, sourceK, targetK);

            % Select edges that match portMap (exactly), if given
            if nargin > 3 && ~isempty(portMap)

                connections = obj.BlockGraph.Edges.Connection;
                connections = connections(k);

                % vertcat the found connections
                connectionsPortMaps = vertcat(connections.PortMap);

                % Search for the input portmap in the connections' port
                % maps
                [foundK, connectionIdx] = ismember(portMap, connectionsPortMaps, 'rows');

                candidateK = k;

                % k now matches the size of the input portmap, and has
                % zeros in the correct location, but the indices
                % are into candidateK, not the original Connections vector.
                k = connectionIdx;

                % Update to match the real connection indices
                k(foundK) = candidateK(k(foundK));
            end
        end

        function ports = connectedInputs(obj, bidx)
            % CONNECTEDINPUTS return vector of connected input port names
            % for the given block index
            % Check incoming ports

            inEid = inedges(obj.BlockGraph, bidx);
            ports = string.empty(0,1);

            if ~isempty(inEid)
                inConnections = obj.BlockGraph.Edges.Connection(inEid);
                pMap = vertcat(inConnections.PortMap);
                ports = pMap(:,2);
            end
        end

        function k = findBlocksBetween(obj, sourceBlocks, targetBlocks)
            % Locate blocks between source and target blocks
            %
            % K = FINDBLOCKSBETWEEN(P, S, T) returns the numeric index K of
            % each block that is reachable by at least one in S, and
            % reaches at least one block in T. If either is empty, the
            % reachability is ignored.
            arguments
                obj
                sourceBlocks
                targetBlocks
            end

            % FINDBLOCKSBETWEEN(P, [], []) returns all blocks.
            sourceBlockIdx = findBlockIndex(obj, sourceBlocks);
            targetBlockIdx = findBlockIndex(obj, targetBlocks);

            % Algorithm: Compute distances for each pair of nodes
            dist = distances(obj.BlockGraph);

            % Assume all are reachable, then filter by reachability to/from
            % target & source blocks
            k = (1:size(obj.BlockGraph.Nodes, 1))';

            % If sourceBlocksIdx is not empty, filter to only blocks
            % reachable from any of the sourceBlocks
            if ~isempty(sourceBlockIdx)
                sourceDists = dist(sourceBlockIdx, :);
                reachableFromSource = find(any(sourceDists ~= Inf, 1));

                k = intersect(k, reachableFromSource);
            end

            % If targetBlocksIdx is not empty, filter to only blocks
            % that can reach any of the targetBlocks
            if ~isempty(targetBlockIdx)
                targetDists = dist(:, targetBlockIdx);
                canReachTarget = find(any(targetDists ~= Inf, 2));

                k = intersect(k, canReachTarget);
            end

        end

        function blocksToRun = getBlocksToRun(obj, fromBlock, toBlock, onlyBlocks)
            % Calculates which blocks to run given the To/From/Only args in
            % runNVPs
            if ~isempty(onlyBlocks)
                if ~isempty(fromBlock) || ~isempty(toBlock)
                    obj.Dispatcher.error(message("bioinfo:Pipeline:ConflictingRunModes"), "PipelineRunError", obj.UUID);
                end

                blocksToRun = findBlockIndex(obj, onlyBlocks);
            else
                blocksToRun = findBlocksBetween(obj, fromBlock, toBlock);

                if isempty(blocksToRun)
                    obj.Dispatcher.warn(message("bioinfo:Pipeline:NoBlocksToRun"), "PipelineRunWarn", obj.UUID);
                    warning(message("bioinfo:Pipeline:NoBlocksToRun"));
                end
            end
        end

        function validateAndCreateResultsDirectory(obj)

            if ~isfolder(obj.RunContext.ResultsDirectory)
                [status, msg] = mkdir(obj.RunContext.ResultsDirectory);

                if ~status
                    error(message("bioinfo:Pipeline:ResultsDirectoryNotValid", msg));
                end
            end

            s = what(obj.RunContext.ResultsDirectory);
            obj.RunContext.AbsoluteResultsDirectory  = s(1).path;
        end
    end

    methods (Access = protected)
        function tf = isEqualExcludeUUID(A, B, treatNaNsEqual)
            % ISEQUALEXCLUDEUUID True if pipelines are structurally equal.
            %
            % TF = isEqualExcludeUUID(A,B) returns TF=true if the two
            % Pipeline objects A and B are structurally equivalent. TNE
            % determines if NaNs are treated as equal. TF returns true if A
            % and B have the same:
            %
            % - metaclass information
            % - Blocks
            % - BlockNames
            % - Connections
            %
            % Blocks, Blocks, and Connections can be in any order. This
            % comparison does not include run artifacts such as results.
            arguments
                A(1,1)
                B(1,1)
                treatNaNsEqual(1,1) logical
            end

            mcA = metaclass(A);
            mcB = metaclass(B);

            if ~isequal(mcA, mcB)
                tf = false;
                return
            end

            [sortedBlockNamesA, ia] = sort(A.BlockNames);
            [sortedBlockNamesB, ib] = sort(B.BlockNames);

            connTblA = A.getConnectionTable();
            connTblB = B.getConnectionTable();

            % Pipelines are equal if the Blocks, BlockNames, and Connections
            if treatNaNsEqual
                tf = isequaln(sortedBlockNamesA, sortedBlockNamesB) &&...
                    isequaln(A.Blocks(ia), B.Blocks(ib)) &&...
                    isequaln(connTblA, connTblB);
            else
                tf = isequal(sortedBlockNamesA, sortedBlockNamesB) &&...
                    isequal(A.Blocks(ia), B.Blocks(ib)) &&...
                    isequal(connTblA, connTblB);
            end
        end

        function connTbl = getConnectionTable(obj)
            % GETCONNECTIONTABLE returns a table with all connection
            % information
            edgesTable = obj.BlockGraph.Edges;

            connTbl = array2table(strings(height(edgesTable), 4), 'VariableNames', {'SrcBlockName', 'OutputPortName', 'DstBlockName', 'InputPortName'});

            for i = 1:numel(edgesTable.Connection)
                newRow = {edgesTable.EndNodes(i,1),...
                    edgesTable.Connection(i).PortMap(1,1),...
                    edgesTable.EndNodes(i,2),...
                    edgesTable.Connection(i).PortMap(1,2)};
                connTbl(i, :) = newRow;
            end
            connTbl = sortrows(connTbl);
        end

        function cpObj = copyElement(obj)
            cpObj = copyElement@bioinfo.pipeline.internal.SemanticModelObject(obj);

            % NOTE: Copying Pipeline does not copy data processed in
            % previous runs

            % Create new MessageBroker, Dispatcher, and Logger
            cpObj.MessageBroker = bioinfo.pipeline.internal.MessageBroker;
            cpObj.Dispatcher    = bioinfo.pipeline.internal.Dispatcher(cpObj.MessageBroker.MessageQueue);
            cpObj.Logger        = bioinfo.pipeline.internal.Logger(cpObj.MessageBroker);

            if ~isempty(obj.Blocks)
                % Create new BlockGraph. Copy Blocks and create fresh
                % ProcessQueues
                cpObj.BlockGraph.Nodes.Block = copy(obj.Blocks);

                % Create new block listeners
                cpObj.BlockSemanticModelChangeListeners = arrayfun(@(blk) listener(blk, 'SemanticModelChanged', @(src, evt) obj.handleBlockModelChanged(evt)), cpObj.BlockGraph.Nodes.Block);

                connections = findConnection(obj);

                if ~isempty(connections)
                    % Create new connections (connections not copyable)
                    srcBlockIdx = findBlockIndex(obj, [connections.SourceBlock]);
                    tgtBlockIdx = findBlockIndex(obj, [connections.TargetBlock]);

                    blocksArray = cpObj.BlockGraph.Nodes.Block;

                    cpObj.BlockGraph.Edges.Connection = arrayfun(...
                        @(c) bioinfo.pipeline.Connection(blocksArray(srcBlockIdx(c)), blocksArray(tgtBlockIdx(c)), connections(c).PortMap),...
                        1:numel(connections))';
                end
            end
        end
    end
end

