view core/tools/DiGraph.m @ 0:e9a9cd732c1e tip

first hg version after svn

% GENERAL GRAPH THEORY CLASS: DiGraph
%
% CAVE: Any special application oriented stuff 
% should be outside here

classdef DiGraph < handle
    
properties
    
    V = sparse(0,1); % boolean for nodes 
    Vlbl = {};
    
    E = sparse(0,0); % sparse weighted edge graph 
end
     
methods

% ---
% Constructor
% ---
function G = DiGraph(V, E, Vlbl)
     
    if nargin == 1
        if  ~isnumeric(V)
            
            G = DiGraph.copy(V);
        else
            % build graph from Edge Adjacency Matrix
            G.V = sparse(find(sum(V + V')) > 0);
            G.E = V;
        end
    end
    
    if nargin >= 2 
        G = DiGraph();
        
        if numel(V) ~= max(size(E))
            error 'wrong dimensions';
        end
        
        % ---
        % We copy existent nodes and edges 
        % ---
        if numel(V) == size(V,2)
            G.V = V;
        else
            G.V = V';
        end
        G.E = E;
        
        % ---
        % We copy the labels for existent nodes, 
        % if supplied.
        % NOTE: the FIND call is neccessary, as logical indexing
        % does not work in this case 
        % TODO: WHY NOT?
        % ---
        if (nargin == 3) && ~isempty(Vlbl)
            G.Vlbl = cell(numel(G.V),1);
            G.Vlbl(find(G.V)) = Vlbl(find(G.V));
        end
        
    end
end
      
% get list of node ids
function out = nodes(G)
    out = find(G.V);
end

% get list of node ids
function out = last_node(G)
    out = find(G.V,1,'last');
end

function w = edge(G, V1, V2)
    
    w = G.E(V1, V2);
end

% get edge list representation
function [val, V1, V2] = edges(G, V1)
    
    if nargin == 1
        
        % get all edges of the graph
        [V1, V2] = find(G.E);
   
        val = zeros(numel(V1), 1);
        
        for i = 1:numel(V1)
            val(i) = G.E(V1(i), V2(i));
        end
        
    elseif nargin == 2
        
        % get all edges coming from or entering this node
        V2 = find(G.E(V1, :));
        val = G.E(V1, V2);
        
        V2 = [ V2 find(G.E(:,V1))'];
        val = [ val G.E(V2, V1)'];
    end
end

% compare two graphs
function out = eq(a,b)
    
    % ---
    % compare graph stats
    % necessary
    % ---
    if a.order ~= b.order
        out = false;
        cprint(3, 'eq: graph node numbers do not match'); 
        return
    end
    
    if a.num_edges ~= b.num_edges
        out = false;
        cprint(3, 'eq: graph edge numbers do not match'); 
        return
    end
    
    
    % ---
    % compare labels and reindex graph b if 
    % necessary
    % ---
    lbla = a.label();
    lblb = b.label();
    
    for i = 1:numel(lbla)
        if ~isempty(lbla{i})
            tmpidx = strcellfind(lblb, lbla{i});
            
            % check for substring problems
            for j = 1:numel(tmpidx)
                if strcmp(lblb{tmpidx(j)}, lbla{i})
                    idx(i) = tmpidx(j);
                    break;
                end
            end
        end
    end
    
    % test on found labels
    num_matching_lbl = sum(idx > 0);
    if ~isempty(lbla) && (num_matching_lbl < a.order)
        out = false;
        cprint(3, 'eq: graph labels do not match'); 
        return
    end
    
    % ---
    % reindex edges and nodes: only replace valid indexes
    % ---
    val_idx = idx > 0;
    idx = idx(val_idx);

    bE(val_idx,val_idx) = b.E(idx,idx);
    bV(val_idx) = b.V(idx);
    
    if sum(a.V ~= bV) > 0 
        out = false;
        cprint(3, 'eq: nodes do not match'); 
        return
    end
    
    if sum(sum(a.E ~= bE)) > 0 
        out = false;
        cprint(3, 'eq: edges do not match'); 
        return
    end
    
    % ---
    % OK, seems these Graphs are the same!!!
    % ---
    out = true;
    
end

% find node via its label
function Vid = node(G, label)
    lbl = G.label();
    
    Vid = strcellfind(lbl, label,1);
end

% add a node
function add_node(G,V)
    
    % set node indicator
    G.V(V) = 1;
    
    G.E(V,V) = 0;
end

% remove a node
function remove_node(G,V)
    
    % remove node 
    G.V(V) = 0;
    
    % remove all edges
    G.E(V,:) = 0;
    G.E(:,V) = 0;
end

% remove an edge
function remove_edge(G, V1, V2)

    G.E(V1, V2) = 0;
end

% add an edge
function add_edge(G, V1, V2, weight)

    G.add_node(V1);
	G.add_node(V2);
    
    if  G.E(V1,V2) == 0

        % save weight
        set_edge(G, V1, V2, weight);
    else
        
        join_edge(G, V1, V2, weight)
    end
end

% ---
% implementation of edge joining,
% to be overwritten for several
% purposes
% ---
function join_edge(G, V1, V2, weight)

    set_edge(G, V1, V2, G.edge(V1, V2) + weight);
end

% ---
% sets the edge without considering earlier weights
% ---
function set_edge(G, V1, V2, weight)
    
    G.E(V1, V2) = weight;
end

% ---
% Graph-specific functions
% ---
function c = children(G, Vid)

    c = find(G.E(Vid,:) > 0);
end

function p = parents(G, Vid)

    p = find(G.E(:,Vid) > 0)';
end

% ---
% Vertex Degree
% ---
function out = degree(G, V)
    
    out = sum(G.E(V,:) > 0) + sum(G.E(:,V) > 0);
end

% ---
% Vertex Degree
% ---
function out = degree_in(G, V)
    
    out = sum(G.E(V,:) > 0);
end

% ---
% Max Degree In
% ---
function out = max_degree_in(G)
    
    out = max(sum(G.E > 0, 2));
end

% ---
% Vertex Degree
% ---
function out = degree_out(G, V)
    
    out = sum(G.E(:,V) > 0);
end

% ---
% Max Degree In
% ---
function out = max_degree_out(G)
    
    out = max(sum(G.E > 0, 1));
end


% ---
% Max Degree
% ---
function out = max_degree(G)
    
    out = max(sum(G.E > 0, 1) + sum(G.E > 0, 2)');
end

% ---
% Max weight
% ---
function out = max_weight(G)
    
    out = max(max(G.E));
end

% ---
% Number of Vertices in Graph
% ---
function out = order(G)
    out = sum(G.V);
end

% ---
% Number of Edges in Graph
% ---
function out = num_edges(G)
    out = sum(sum(G.E ~= 0));
end

% ---
% Number of Vertices in Graph
% ---
function out = cardinality(G)
    out = order(G);
end

function Gs = unconnected_subgraph(G)
    Vtmp = (sum(G.E + G.E') == 0);
    Etmp = sparse(size(G.E, 1), size(G.E, 1));

    Gs = DiGraph(Vtmp, Etmp, G.label());
end

% ---
% return string labels for a (set of) node(S)
% ---
function out = label(G, Vid)
    
    out = {};
    if nargin == 1
        % maybe much faster for whole graph
        if (numel(G.Vlbl) < G.order() || isempty(G.Vlbl))
            
            out = cell(numel(G.V), 1);
            for i = 1:numel(Vid)
                out{i} = sprintf('%d', Vid(i));
            end
        else
            out = G.Vlbl;
        end
        
    elseif nargin == 2
        for i = 1:numel(Vid)

            if (numel(G.Vlbl) < Vid(i)) || isempty(G.Vlbl{Vid(i)})
                
                if numel(Vid) > 1
                    out{i} = sprintf('%d', Vid(i));
                else
                    out = sprintf('%d', Vid(i));
                end
            else
                if numel(Vid) > 1
                    out{i} = G.Vlbl{Vid(i)};
                else
                    out = G.Vlbl{Vid(i)};
                end
            end
        end
    end
end


% -----------------------------------------------------------------
% Graph theory algorithms
% ---


% ---
% sets all the main diagonal edges to zero
% ---
function remove_cycles_length1(G)
    
    for i = G.nodes();
        G.E(i,i) = 0;
    end  
end


% ---
% Returns whether a given graph has cycles or not,
% using the SCC search
% ---
function out = isAcyclic(G)
     % ---
     % We get all the sccs in the DiGraph, and 
     % return true if none of the sccs has more than 
     % one node
     % ---
     
     % check for self-loops
     if sum(abs(diag(G.E))) > 0
         out = 0;
         return;
     end
     
     [~, s, ~] = strongly_connected_components(G);
     
     if max(s) < 2
         out = 1;
     else
         out = 0;
     end
end

% ---
% All SCC's ordered by number of nodes in graph
% 
% this should also be able to return 
% the SCC of just one node
% ---
function [Gs, s, id] = strongly_connected_components(G, Vin)

    % marking
    marked = zeros(size(G.V));
    
    % ---
    % two stacks, 
    % one: has not been assigned to scc
    % two: unclear if in same scc
    % ---
    stack1 = CStack();
    stack2 = CStack();
    
    % initialise scc ids
    id = zeros(size(G.V)) - 1;  % assigned scc?
    idctr = 1;
        
    % initialise graph ordering
    preorder = zeros(G.order, 1);
    prectr = 0;
   
    for v = G.nodes();
        if ~marked(v)
          dfs(G, v);
        end
    end
    
    % ---
    % create subgraphs (DiGraph here) for sccs
    % ---
    if nargin == 1
        
        s = zeros(idctr-1,1);
        for idctr2 = 1:idctr-1
            Vtmp = (id == idctr2);
            Emask = sparse(size(G.E, 1), size(G.E, 1));
            Emask(Vtmp,Vtmp) = 1;
            Etmp = G.E .* Emask;
    
            Gs(idctr2) = DiGraph(Vtmp, Etmp, G.Vlbl);
            s(idctr2) = Gs(idctr2).order();
        end
         
        % ---
        % order by number of nodes
        % ---
        [s, idx] = sort(s,'descend');
        Gs = Gs(idx);
        Gmax = Gs(1);
        
    else
        % ---
        % get just the scc for the questioned node
        % ---
        Vtmp = (id == id(Vin));
        Emask = sparse(size(G.E, 1), size(G.E, 1));
        Emask(Vtmp,Vtmp) = 1;
        Etmp = G.E .* Emask;

        Gs = DiGraph(Vtmp, Etmp);
        s = Gs.order();
    end

    % ---
    % NOTE: THIS IS A NESTED DFS BASED GRAPH ORDERING
    % ---
    function dfs(G, v)
        
        % mark this node as visited
        marked(v) = 1;
        
        preorder(v) = prectr;
        prectr = prectr + 1;
        
        % push into both stacks
        stack1.push(v);
        stack2.push(v);
        
        % ---
        % dfs
        % ---
        for w = G.children(v)
            if ~marked(w)
                % ---
                % traverse into dfs if not yet visited
                % ---
                dfs(G, w);
                  
            elseif id(w) == -1
                
                % ---
                % w has not yet been assigned to a strongly connected
                % component
                % ---
                while ((preorder(stack2.top()) > preorder(w)))
                    stack2.pop();
                end
            end   
        end
        
        % ---
        % found scc containing v
        % ---
        if (stack2.top() == v)
            stack2.pop();
            
            w = -1;
            while (w ~= v)
                
                % ---
                % collect all the nodes of this scc
                % ---
                w = stack1.pop();
                id(w) = idctr;
            end
            idctr = idctr + 1;
        end
        
    end % function dfs 
end

function [Gs, s, id] = connected_components(G, varargin)
    % ---
    % get all connected subgraphs:
    % ---
    
    % make edge matrix undirected
    G2 = DiGraph(Graph(G));
    
    [GsGraph, s, id] = strongly_connected_components(G2, varargin{:});
    
    % get the actual subgraps
    
    for i =1:numel(GsGraph)
        Gs(i) = G.subgraph(GsGraph(i).nodes);
    end
end

% ---
% creates new graph just containing the 
% specified nodes and optionally specified edges
% nodes can be specified using 
% a. the binary G.V structure
% b. or node indices 
% ---
function G2 = subgraph(G, V, E)
    if nargin == 2
        E = [];
    end
    
    % ---
    % create new graph as copy ofthe old
    % ---
    
    G2 = feval(class(G), G);
    
    % ---
    % reset nodes and edges
    % NOTE: we determine the input format first
    % ---
    if (max(V) == 1 && numel(V) > 1) || max(V) == 0
        
        G2.remove_node(find(~V));
    else
        
        G2.remove_node(setdiff(1:numel(G.V), V));
    end
    if ~isempty(E)   
        G2.E = E;
    end
end

% ---
% joins the information of graph G2 with
% this GRaph, not duplicating nodes.
% ---
function add_graph(G, G2)
    
    % determine if symmetric  edges have to be added
    clas = cat(1, class(G2), superclasses(G2));
    if strcellfind(clas, 'Graph')
        add_symm = 1;
    else
        add_symm = 0;
    end
            
    % ---
    % add all nodes and edges in G2
    % ---
    for V = G2.nodes();
        
        G.add_node(V);
    end
    
    % ---
    % NOTE / TODO: 
    % this LABEL inheritance is far too expensive when 
    % creating many copiesof the same graph
    % Its also unnessecary for lots of them 
    % except for debugging :(. 
    % ---
    G.Vlbl = cell(1,numel(G2.V));
    G.Vlbl(G2.nodes()) = G2.label(G2.nodes());
    
    % ---
    % get all edges in G2
    % ---
    [val, V1, V2] = edges(G2);
    
    for i = 1:numel(val)
        
        % ---
        % add edges to graph
        % ---
        G.add_edge(V1(i), V2(i), val(i));
        
        if add_symm
            % ---
            % add symmetric edges to graph
            % ---
            G.add_edge(V2(i), V1(i), val(i));
        end
        
    end
end

% ---
% substracts the edges in G2 from
% this GRaph, removing not connected nodes.
% ---
function Gout = minus(a, b)
    Gout = feval(class(a),a);
    Gout.remove_graph(b);
end

function remove_graph(G, G2)
    
    % determine if symmetric  edges have to be added
    clas = cat(1, class(G2), superclasses(G2));
    if strcellfind(clas, 'Graph')
        symm = 1;
    else
        symm = 0;
    end
    
    % remove edges
    [val, V1, V2] = G2.edges();
    for j = 1:numel(val)
        
        % remove specified edges
        G.remove_edge(V1(j), V2(j));
        
        % remove symmetric edges if subtracting symm graph
        if symm

            G.remove_edge(V2(j), V1(j));
        end
    end

    % ---
    % Note : we only remove nodes with no 
    % remaining incoming edges
    % ---
    V = G2.nodes();
    for j = 1:numel(V)

        if G.degree(V(j)) == 0

            G.remove_node(V(j));
        end
    end
    
end

% ---
% compact graph representation
% ---
function [E, labels] = compact(G)
    
    % ---
    % get nodes and create a reverse index
    % ---
    Vidx = sparse(G.nodes,1,1:G.order());
    [w, V1, V2]  = G.edges();
    
    % create compact adjacency matrix
    E = sparse(Vidx(V1), Vidx(V2),w, G.order(), G.order());
    
    labels = G.label(G.nodes());
end

% ---
% determines if Edges in G2 are the same as in G
% ---
function [out] = isSubgraph(G, G2)
    
    [val, V1, V2] = G2.edges();
    validE = false(numel(V1), 1);
    
    i = 1;
    while i <= numel(V1)
        
        % ---
        % Test if edge exists in other graph
        % ---
        if G.edge(V1(i),V2(i)) == val(i) 
           out = 0;
           return;
        end
        i = i +1 ;
    end
    
    % ---
    % Test labels
    % ---
    if ~isempty(G.Vlbl) && ~isempty(G2.Vlbl)
        
        V = G2.nodes();
        i = 1;
        while i <= numel(V)
            if strcmp(G.label(V(i)), G2.label(V(i))) ~= 0
                out = 0;
                return;
            end
            i = i + 1;
        end 
    end
    
    out = 1;
end

% ---
% Visualise the Similarity Graph
% ---
function visualise(G)

   % ---
   % get colormap for edge weights
   % ---
   cmap = jet(100);
   
   % ---
   % NOTE: we now get the weight colors for all edges
   % get maximum weight and all edges
   % ---
   [colors, V1, V2]  = G.edges();
   w = G.max_weight();
   
   % normalise colors
   colors = max(1,round((colors / w) * 100));

   % get node labels
   V1lbl = G.label(V1);
   V2lbl = G.label(V2);

   % ---
   % compose edgecolor matrix
   % ---
   ec = cell(numel(V1), 3);
   for i = 1:numel(V1)
       ec(i,:) = {V1lbl{i}, V2lbl{i}, cmap(colors(i),:)};
   end
   
   % ---
   % For Performance Issues
   % We get the compact Graph and
   % !hope! for the labels to correspond to those above
   % ---
   [E, labels] = compact(G);
   
   % determine if symmetric Graph
    clas = cat(1, class(G), superclasses(G));
    if strcellfind(clas, 'Graph')
        symm = 1;
    else
        symm = 0;
    end
   
   graphViz4Matlab('-adjMat',E, ...
    '-nodeLabels',labels, ...
    '-edgeColors', ec, ...
    '-undirected', symm); 
end
end


methods (Static)
   function G = copy(Gin)
       
       if strcmp(class(Gin), 'DiGraph')
           
           G = DiGraph(Gin.V, Gin.E);
           G.Vlbl = Gin.Vlbl;
           
       else
           warning ('cannot copy graph, casting instead')
           G = DiGraph.cast(Gin);
       end
        
   end
   
   function G = cast(Gin)
       
        % ---
        % this uses an imput grpaph
        % to create a new digraph
        % ---
        G = DiGraph();
        
        % ---
        % Add the input graph to the empty one
        % ---
        G.add_graph(Gin);
        
   end
end
end