function [ParentList DiscoveryList DistanceList] = DijkstraDigraph(EdgeListG, Weights, startV, n)
%[ParentList DiscoveryList DistanceList] = DijkstraDigraph(EdgeListG,Weights, startV, n)

%We initialize the tree growing edge lists:
EdgeListT = []; %initially T consists of just the starting vertex startV (no edges)
%The boundary edges are all edges directed from startV
EdgeListGHeadEndpoints = EdgeListG(:,2); %only the terminal vertices (second column) are listed.
EdgeListGTailEndpts = EdgeListG(:,1); %Head endpoints of edges
a = find(EdgeListGTailEndpts == startV); 
BdyEdges = EdgeListG(a,:);
%The boundary edges are the only edges that are initially "used"
UnusedEdges = setdiff(EdgeListG, BdyEdges, 'rows');


DiscoveryList = [startV];
%Initially, only the distance from the root to itself is known = 0.
DistanceList = [0];
ParentList = zeros(1,n);
%A zero entry in ParentList will mean the vertex has no parent "NULL"
%In the initial tree T = [startV], no vertex has a parent (and the root
%startV never will)
CurrDist = Inf*ones(1,n);
CurrDist(startV)= 0;

%Before entering into the while loop, we do the first iteration separately
%because initially EdgeListT is empty (so does not have startV)
%First Edge Selection Step (Dijkstra)
%Need to find a lightest boundary edge
[lightest ind] = min(Weights(BdyEdges(:,3)));
e = BdyEdges(ind(1),:); %third entry is the index of edge w.r.t. EdgeListG
nonTreeEndpt = e(2);

DiscoveryList = [DiscoveryList nonTreeEndpt];
ParentList(nonTreeEndpt) = startV;
CurrDist(nonTreeEndpt) = Weights(BdyEdges(ind(1),3));
DistanceList = [DistanceList CurrDist(nonTreeEndpt)];


%We need to delete any edge of BdyEdges that has nonTreeEndpt
%as one of its endpoints.

BdyEdgesEndpts = BdyEdges(:,1:2);
[a b] = find(BdyEdgesEndpts == nonTreeEndpt); 
BdyEdges(a,:) = []; BdyEdgesEndpts(a,:) = [];
%Next, we find the unused edges that have tail =  nonTreeEndpt.  
UnusedEdgesTailEndpts = UnusedEdges(:,1);
a = find(UnusedEdgesTailEndpts==nonTreeEndpt);
%We add these edges to BdyEdges if their head endpoints are not in the tree
%(i.e., are in DiscoveryList):
for i = 1:length(a)
    ai = a(i);
    if ~ismember(UnusedEdges(ai,2), DiscoveryList)
        BdyEdges = [BdyEdges; UnusedEdges(ai,:)];
    end
end
%Now we delete these same edges from Unused edges:
for i = 1:length(a)
    ai = a(i);
    if ~ismember(UnusedEdges(ai,2), DiscoveryList)
        UnusedEdges(ai,:)=[];
    end
end
BdyEdgesEndpts = BdyEdges(:,1:2);
UnusedEdgesEndpts = UnusedEdges(:,1:2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
EdgeListT = [EdgeListT; e];
EdgeListTEndpts= EdgeListT(:,1:2);
EdgeListTTailEndpts= EdgeListT(:,1);
EdgeListTHeadEndpts= EdgeListT(:,2);
BdyEdgesHeadEndpts = BdyEdges(:,2);
BdyEdgesTailEndpts = BdyEdges(:,1);

while length(DiscoveryList)<n
    %Next Edge Selection Step (Dijkstra)
    %Need to find a boundary edge such that when its weight is added to
    %CurrDist of its tree endpt (=the shortest distance to root) it is a
    %minimum
    CurrDistTemp = CurrDist;
    Best = [Inf 0 0 0]; %will store min non tree endpt new distance (first component), and 
    %non tree endpt (second component intialized to 0), tree endpt (third component), 
    %and index of corresponding bdy edge (forth component)  
    for i = 1:n
            if ismember(i,EdgeListTEndpts) %i is a tree vertex
            a = find(BdyEdgesTailEndpts == i); %a = indices of bdy edges starting at i
            for j = 1:length(a)
                
                    nonTreeEndpt = BdyEdges(a(j),2);
                
                if CurrDist(i) + Weights(BdyEdges(a(j),3)) < Best(1) %nonTreeEndpt is candidate 
                                                    %for next discovery vertex
                    Best = [CurrDist(i) + Weights(BdyEdges(a(j),3))  nonTreeEndpt i BdyEdges(a(j),3)];
                end
            end
            end
    end
    nonTreeEndpt = Best(2);
    DiscoveryList = [DiscoveryList Best(2)];
    ParentList(Best(2)) = Best(3);
    CurrDist(Best(2)) = Best(1);
    DistanceList = [DistanceList Best(1)];
    %%%%%%BdyEdges Update for Directed Program:
    %We need to delete any edge of BdyEdges that has nonTreeEndpt
    %as one of its endpoints.
    BdyEdgesEndpts = BdyEdges(:,1:2);
    [a b] = find(BdyEdgesEndpts == nonTreeEndpt);
    BdyEdges(a,:) = []; BdyEdgesEndpts(a,:) = [];
    %Next, we find the unused edges that are incident to nonTreeEndpt. 
    UnusedEdgesTailEndpts = UnusedEdges(:,1);
    [a b] = find(UnusedEdgesTailEndpts==nonTreeEndpt);
    %We add these edges to BdyEdges if their head endpoints are not in the tree
    %(i.e., are in DiscoveryList):
    for i = 1:length(a)
        ai = a(i);
       if ~ismember(UnusedEdges(ai,2), DiscoveryList)
          BdyEdges = [BdyEdges; UnusedEdges(ai,:)];
       end
    end
    %Now we delete these same edges from UnusedEdges
     for i = 1:length(a)
        ai = a(i);
       if ~ismember(UnusedEdges(ai,2), DiscoveryList)
          UnusedEdges(ai,:)=[];
       end
    end  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    BdyEdgesEndpts = BdyEdges(:,1:2);
    UnusedEdgesEndpts = UnusedEdges(:,1:2);
    EdgeListT = [EdgeListT; [Best([3 2 4])]];
    EdgeListTEndpts= EdgeListT(:,1:2);
    EdgeListTTailEndpts= EdgeListT(:,1);
    EdgeListTHeadEndpts= EdgeListT(:,2);
    BdyEdgesHeadEndpts = BdyEdges(:,2);
    BdyEdgesTailEndpts = BdyEdges(:,1);
end
%Finally we format the outputs
ParentList2 = zeros(2,n);
ParentList2(1,:) = 1:n;
ParentList2(2,:) = ParentList;
ParentList = ParentList2';

