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

first hg version after svn

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/core/tools/DiGraph.m	Tue Feb 10 15:05:51 2015 +0000
@@ -0,0 +1,853 @@
+
+% 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
\ No newline at end of file