wolffd@0: 
wolffd@0: % GENERAL GRAPH THEORY CLASS: DiGraph
wolffd@0: %
wolffd@0: % CAVE: Any special application oriented stuff 
wolffd@0: % should be outside here
wolffd@0: 
wolffd@0: classdef DiGraph < handle
wolffd@0:     
wolffd@0: properties
wolffd@0:     
wolffd@0:     V = sparse(0,1); % boolean for nodes 
wolffd@0:     Vlbl = {};
wolffd@0:     
wolffd@0:     E = sparse(0,0); % sparse weighted edge graph 
wolffd@0: end
wolffd@0:      
wolffd@0: methods
wolffd@0: 
wolffd@0: % ---
wolffd@0: % Constructor
wolffd@0: % ---
wolffd@0: function G = DiGraph(V, E, Vlbl)
wolffd@0:      
wolffd@0:     if nargin == 1
wolffd@0:         if  ~isnumeric(V)
wolffd@0:             
wolffd@0:             G = DiGraph.copy(V);
wolffd@0:         else
wolffd@0:             % build graph from Edge Adjacency Matrix
wolffd@0:             G.V = sparse(find(sum(V + V')) > 0);
wolffd@0:             G.E = V;
wolffd@0:         end
wolffd@0:     end
wolffd@0:     
wolffd@0:     if nargin >= 2 
wolffd@0:         G = DiGraph();
wolffd@0:         
wolffd@0:         if numel(V) ~= max(size(E))
wolffd@0:             error 'wrong dimensions';
wolffd@0:         end
wolffd@0:         
wolffd@0:         % ---
wolffd@0:         % We copy existent nodes and edges 
wolffd@0:         % ---
wolffd@0:         if numel(V) == size(V,2)
wolffd@0:             G.V = V;
wolffd@0:         else
wolffd@0:             G.V = V';
wolffd@0:         end
wolffd@0:         G.E = E;
wolffd@0:         
wolffd@0:         % ---
wolffd@0:         % We copy the labels for existent nodes, 
wolffd@0:         % if supplied.
wolffd@0:         % NOTE: the FIND call is neccessary, as logical indexing
wolffd@0:         % does not work in this case 
wolffd@0:         % TODO: WHY NOT?
wolffd@0:         % ---
wolffd@0:         if (nargin == 3) && ~isempty(Vlbl)
wolffd@0:             G.Vlbl = cell(numel(G.V),1);
wolffd@0:             G.Vlbl(find(G.V)) = Vlbl(find(G.V));
wolffd@0:         end
wolffd@0:         
wolffd@0:     end
wolffd@0: end
wolffd@0:       
wolffd@0: % get list of node ids
wolffd@0: function out = nodes(G)
wolffd@0:     out = find(G.V);
wolffd@0: end
wolffd@0: 
wolffd@0: % get list of node ids
wolffd@0: function out = last_node(G)
wolffd@0:     out = find(G.V,1,'last');
wolffd@0: end
wolffd@0: 
wolffd@0: function w = edge(G, V1, V2)
wolffd@0:     
wolffd@0:     w = G.E(V1, V2);
wolffd@0: end
wolffd@0: 
wolffd@0: % get edge list representation
wolffd@0: function [val, V1, V2] = edges(G, V1)
wolffd@0:     
wolffd@0:     if nargin == 1
wolffd@0:         
wolffd@0:         % get all edges of the graph
wolffd@0:         [V1, V2] = find(G.E);
wolffd@0:    
wolffd@0:         val = zeros(numel(V1), 1);
wolffd@0:         
wolffd@0:         for i = 1:numel(V1)
wolffd@0:             val(i) = G.E(V1(i), V2(i));
wolffd@0:         end
wolffd@0:         
wolffd@0:     elseif nargin == 2
wolffd@0:         
wolffd@0:         % get all edges coming from or entering this node
wolffd@0:         V2 = find(G.E(V1, :));
wolffd@0:         val = G.E(V1, V2);
wolffd@0:         
wolffd@0:         V2 = [ V2 find(G.E(:,V1))'];
wolffd@0:         val = [ val G.E(V2, V1)'];
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0: % compare two graphs
wolffd@0: function out = eq(a,b)
wolffd@0:     
wolffd@0:     % ---
wolffd@0:     % compare graph stats
wolffd@0:     % necessary
wolffd@0:     % ---
wolffd@0:     if a.order ~= b.order
wolffd@0:         out = false;
wolffd@0:         cprint(3, 'eq: graph node numbers do not match'); 
wolffd@0:         return
wolffd@0:     end
wolffd@0:     
wolffd@0:     if a.num_edges ~= b.num_edges
wolffd@0:         out = false;
wolffd@0:         cprint(3, 'eq: graph edge numbers do not match'); 
wolffd@0:         return
wolffd@0:     end
wolffd@0:     
wolffd@0:     
wolffd@0:     % ---
wolffd@0:     % compare labels and reindex graph b if 
wolffd@0:     % necessary
wolffd@0:     % ---
wolffd@0:     lbla = a.label();
wolffd@0:     lblb = b.label();
wolffd@0:     
wolffd@0:     for i = 1:numel(lbla)
wolffd@0:         if ~isempty(lbla{i})
wolffd@0:             tmpidx = strcellfind(lblb, lbla{i});
wolffd@0:             
wolffd@0:             % check for substring problems
wolffd@0:             for j = 1:numel(tmpidx)
wolffd@0:                 if strcmp(lblb{tmpidx(j)}, lbla{i})
wolffd@0:                     idx(i) = tmpidx(j);
wolffd@0:                     break;
wolffd@0:                 end
wolffd@0:             end
wolffd@0:         end
wolffd@0:     end
wolffd@0:     
wolffd@0:     % test on found labels
wolffd@0:     num_matching_lbl = sum(idx > 0);
wolffd@0:     if ~isempty(lbla) && (num_matching_lbl < a.order)
wolffd@0:         out = false;
wolffd@0:         cprint(3, 'eq: graph labels do not match'); 
wolffd@0:         return
wolffd@0:     end
wolffd@0:     
wolffd@0:     % ---
wolffd@0:     % reindex edges and nodes: only replace valid indexes
wolffd@0:     % ---
wolffd@0:     val_idx = idx > 0;
wolffd@0:     idx = idx(val_idx);
wolffd@0: 
wolffd@0:     bE(val_idx,val_idx) = b.E(idx,idx);
wolffd@0:     bV(val_idx) = b.V(idx);
wolffd@0:     
wolffd@0:     if sum(a.V ~= bV) > 0 
wolffd@0:         out = false;
wolffd@0:         cprint(3, 'eq: nodes do not match'); 
wolffd@0:         return
wolffd@0:     end
wolffd@0:     
wolffd@0:     if sum(sum(a.E ~= bE)) > 0 
wolffd@0:         out = false;
wolffd@0:         cprint(3, 'eq: edges do not match'); 
wolffd@0:         return
wolffd@0:     end
wolffd@0:     
wolffd@0:     % ---
wolffd@0:     % OK, seems these Graphs are the same!!!
wolffd@0:     % ---
wolffd@0:     out = true;
wolffd@0:     
wolffd@0: end
wolffd@0: 
wolffd@0: % find node via its label
wolffd@0: function Vid = node(G, label)
wolffd@0:     lbl = G.label();
wolffd@0:     
wolffd@0:     Vid = strcellfind(lbl, label,1);
wolffd@0: end
wolffd@0: 
wolffd@0: % add a node
wolffd@0: function add_node(G,V)
wolffd@0:     
wolffd@0:     % set node indicator
wolffd@0:     G.V(V) = 1;
wolffd@0:     
wolffd@0:     G.E(V,V) = 0;
wolffd@0: end
wolffd@0: 
wolffd@0: % remove a node
wolffd@0: function remove_node(G,V)
wolffd@0:     
wolffd@0:     % remove node 
wolffd@0:     G.V(V) = 0;
wolffd@0:     
wolffd@0:     % remove all edges
wolffd@0:     G.E(V,:) = 0;
wolffd@0:     G.E(:,V) = 0;
wolffd@0: end
wolffd@0: 
wolffd@0: % remove an edge
wolffd@0: function remove_edge(G, V1, V2)
wolffd@0: 
wolffd@0:     G.E(V1, V2) = 0;
wolffd@0: end
wolffd@0: 
wolffd@0: % add an edge
wolffd@0: function add_edge(G, V1, V2, weight)
wolffd@0: 
wolffd@0:     G.add_node(V1);
wolffd@0: 	G.add_node(V2);
wolffd@0:     
wolffd@0:     if  G.E(V1,V2) == 0
wolffd@0: 
wolffd@0:         % save weight
wolffd@0:         set_edge(G, V1, V2, weight);
wolffd@0:     else
wolffd@0:         
wolffd@0:         join_edge(G, V1, V2, weight)
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % implementation of edge joining,
wolffd@0: % to be overwritten for several
wolffd@0: % purposes
wolffd@0: % ---
wolffd@0: function join_edge(G, V1, V2, weight)
wolffd@0: 
wolffd@0:     set_edge(G, V1, V2, G.edge(V1, V2) + weight);
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % sets the edge without considering earlier weights
wolffd@0: % ---
wolffd@0: function set_edge(G, V1, V2, weight)
wolffd@0:     
wolffd@0:     G.E(V1, V2) = weight;
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % Graph-specific functions
wolffd@0: % ---
wolffd@0: function c = children(G, Vid)
wolffd@0: 
wolffd@0:     c = find(G.E(Vid,:) > 0);
wolffd@0: end
wolffd@0: 
wolffd@0: function p = parents(G, Vid)
wolffd@0: 
wolffd@0:     p = find(G.E(:,Vid) > 0)';
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % Vertex Degree
wolffd@0: % ---
wolffd@0: function out = degree(G, V)
wolffd@0:     
wolffd@0:     out = sum(G.E(V,:) > 0) + sum(G.E(:,V) > 0);
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % Vertex Degree
wolffd@0: % ---
wolffd@0: function out = degree_in(G, V)
wolffd@0:     
wolffd@0:     out = sum(G.E(V,:) > 0);
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % Max Degree In
wolffd@0: % ---
wolffd@0: function out = max_degree_in(G)
wolffd@0:     
wolffd@0:     out = max(sum(G.E > 0, 2));
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % Vertex Degree
wolffd@0: % ---
wolffd@0: function out = degree_out(G, V)
wolffd@0:     
wolffd@0:     out = sum(G.E(:,V) > 0);
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % Max Degree In
wolffd@0: % ---
wolffd@0: function out = max_degree_out(G)
wolffd@0:     
wolffd@0:     out = max(sum(G.E > 0, 1));
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: % ---
wolffd@0: % Max Degree
wolffd@0: % ---
wolffd@0: function out = max_degree(G)
wolffd@0:     
wolffd@0:     out = max(sum(G.E > 0, 1) + sum(G.E > 0, 2)');
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % Max weight
wolffd@0: % ---
wolffd@0: function out = max_weight(G)
wolffd@0:     
wolffd@0:     out = max(max(G.E));
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % Number of Vertices in Graph
wolffd@0: % ---
wolffd@0: function out = order(G)
wolffd@0:     out = sum(G.V);
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % Number of Edges in Graph
wolffd@0: % ---
wolffd@0: function out = num_edges(G)
wolffd@0:     out = sum(sum(G.E ~= 0));
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % Number of Vertices in Graph
wolffd@0: % ---
wolffd@0: function out = cardinality(G)
wolffd@0:     out = order(G);
wolffd@0: end
wolffd@0: 
wolffd@0: function Gs = unconnected_subgraph(G)
wolffd@0:     Vtmp = (sum(G.E + G.E') == 0);
wolffd@0:     Etmp = sparse(size(G.E, 1), size(G.E, 1));
wolffd@0: 
wolffd@0:     Gs = DiGraph(Vtmp, Etmp, G.label());
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % return string labels for a (set of) node(S)
wolffd@0: % ---
wolffd@0: function out = label(G, Vid)
wolffd@0:     
wolffd@0:     out = {};
wolffd@0:     if nargin == 1
wolffd@0:         % maybe much faster for whole graph
wolffd@0:         if (numel(G.Vlbl) < G.order() || isempty(G.Vlbl))
wolffd@0:             
wolffd@0:             out = cell(numel(G.V), 1);
wolffd@0:             for i = 1:numel(Vid)
wolffd@0:                 out{i} = sprintf('%d', Vid(i));
wolffd@0:             end
wolffd@0:         else
wolffd@0:             out = G.Vlbl;
wolffd@0:         end
wolffd@0:         
wolffd@0:     elseif nargin == 2
wolffd@0:         for i = 1:numel(Vid)
wolffd@0: 
wolffd@0:             if (numel(G.Vlbl) < Vid(i)) || isempty(G.Vlbl{Vid(i)})
wolffd@0:                 
wolffd@0:                 if numel(Vid) > 1
wolffd@0:                     out{i} = sprintf('%d', Vid(i));
wolffd@0:                 else
wolffd@0:                     out = sprintf('%d', Vid(i));
wolffd@0:                 end
wolffd@0:             else
wolffd@0:                 if numel(Vid) > 1
wolffd@0:                     out{i} = G.Vlbl{Vid(i)};
wolffd@0:                 else
wolffd@0:                     out = G.Vlbl{Vid(i)};
wolffd@0:                 end
wolffd@0:             end
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: % -----------------------------------------------------------------
wolffd@0: % Graph theory algorithms
wolffd@0: % ---
wolffd@0: 
wolffd@0: 
wolffd@0: % ---
wolffd@0: % sets all the main diagonal edges to zero
wolffd@0: % ---
wolffd@0: function remove_cycles_length1(G)
wolffd@0:     
wolffd@0:     for i = G.nodes();
wolffd@0:         G.E(i,i) = 0;
wolffd@0:     end  
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: % ---
wolffd@0: % Returns whether a given graph has cycles or not,
wolffd@0: % using the SCC search
wolffd@0: % ---
wolffd@0: function out = isAcyclic(G)
wolffd@0:      % ---
wolffd@0:      % We get all the sccs in the DiGraph, and 
wolffd@0:      % return true if none of the sccs has more than 
wolffd@0:      % one node
wolffd@0:      % ---
wolffd@0:      
wolffd@0:      % check for self-loops
wolffd@0:      if sum(abs(diag(G.E))) > 0
wolffd@0:          out = 0;
wolffd@0:          return;
wolffd@0:      end
wolffd@0:      
wolffd@0:      [~, s, ~] = strongly_connected_components(G);
wolffd@0:      
wolffd@0:      if max(s) < 2
wolffd@0:          out = 1;
wolffd@0:      else
wolffd@0:          out = 0;
wolffd@0:      end
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % All SCC's ordered by number of nodes in graph
wolffd@0: % 
wolffd@0: % this should also be able to return 
wolffd@0: % the SCC of just one node
wolffd@0: % ---
wolffd@0: function [Gs, s, id] = strongly_connected_components(G, Vin)
wolffd@0: 
wolffd@0:     % marking
wolffd@0:     marked = zeros(size(G.V));
wolffd@0:     
wolffd@0:     % ---
wolffd@0:     % two stacks, 
wolffd@0:     % one: has not been assigned to scc
wolffd@0:     % two: unclear if in same scc
wolffd@0:     % ---
wolffd@0:     stack1 = CStack();
wolffd@0:     stack2 = CStack();
wolffd@0:     
wolffd@0:     % initialise scc ids
wolffd@0:     id = zeros(size(G.V)) - 1;  % assigned scc?
wolffd@0:     idctr = 1;
wolffd@0:         
wolffd@0:     % initialise graph ordering
wolffd@0:     preorder = zeros(G.order, 1);
wolffd@0:     prectr = 0;
wolffd@0:    
wolffd@0:     for v = G.nodes();
wolffd@0:         if ~marked(v)
wolffd@0:           dfs(G, v);
wolffd@0:         end
wolffd@0:     end
wolffd@0:     
wolffd@0:     % ---
wolffd@0:     % create subgraphs (DiGraph here) for sccs
wolffd@0:     % ---
wolffd@0:     if nargin == 1
wolffd@0:         
wolffd@0:         s = zeros(idctr-1,1);
wolffd@0:         for idctr2 = 1:idctr-1
wolffd@0:             Vtmp = (id == idctr2);
wolffd@0:             Emask = sparse(size(G.E, 1), size(G.E, 1));
wolffd@0:             Emask(Vtmp,Vtmp) = 1;
wolffd@0:             Etmp = G.E .* Emask;
wolffd@0:     
wolffd@0:             Gs(idctr2) = DiGraph(Vtmp, Etmp, G.Vlbl);
wolffd@0:             s(idctr2) = Gs(idctr2).order();
wolffd@0:         end
wolffd@0:          
wolffd@0:         % ---
wolffd@0:         % order by number of nodes
wolffd@0:         % ---
wolffd@0:         [s, idx] = sort(s,'descend');
wolffd@0:         Gs = Gs(idx);
wolffd@0:         Gmax = Gs(1);
wolffd@0:         
wolffd@0:     else
wolffd@0:         % ---
wolffd@0:         % get just the scc for the questioned node
wolffd@0:         % ---
wolffd@0:         Vtmp = (id == id(Vin));
wolffd@0:         Emask = sparse(size(G.E, 1), size(G.E, 1));
wolffd@0:         Emask(Vtmp,Vtmp) = 1;
wolffd@0:         Etmp = G.E .* Emask;
wolffd@0: 
wolffd@0:         Gs = DiGraph(Vtmp, Etmp);
wolffd@0:         s = Gs.order();
wolffd@0:     end
wolffd@0: 
wolffd@0:     % ---
wolffd@0:     % NOTE: THIS IS A NESTED DFS BASED GRAPH ORDERING
wolffd@0:     % ---
wolffd@0:     function dfs(G, v)
wolffd@0:         
wolffd@0:         % mark this node as visited
wolffd@0:         marked(v) = 1;
wolffd@0:         
wolffd@0:         preorder(v) = prectr;
wolffd@0:         prectr = prectr + 1;
wolffd@0:         
wolffd@0:         % push into both stacks
wolffd@0:         stack1.push(v);
wolffd@0:         stack2.push(v);
wolffd@0:         
wolffd@0:         % ---
wolffd@0:         % dfs
wolffd@0:         % ---
wolffd@0:         for w = G.children(v)
wolffd@0:             if ~marked(w)
wolffd@0:                 % ---
wolffd@0:                 % traverse into dfs if not yet visited
wolffd@0:                 % ---
wolffd@0:                 dfs(G, w);
wolffd@0:                   
wolffd@0:             elseif id(w) == -1
wolffd@0:                 
wolffd@0:                 % ---
wolffd@0:                 % w has not yet been assigned to a strongly connected
wolffd@0:                 % component
wolffd@0:                 % ---
wolffd@0:                 while ((preorder(stack2.top()) > preorder(w)))
wolffd@0:                     stack2.pop();
wolffd@0:                 end
wolffd@0:             end   
wolffd@0:         end
wolffd@0:         
wolffd@0:         % ---
wolffd@0:         % found scc containing v
wolffd@0:         % ---
wolffd@0:         if (stack2.top() == v)
wolffd@0:             stack2.pop();
wolffd@0:             
wolffd@0:             w = -1;
wolffd@0:             while (w ~= v)
wolffd@0:                 
wolffd@0:                 % ---
wolffd@0:                 % collect all the nodes of this scc
wolffd@0:                 % ---
wolffd@0:                 w = stack1.pop();
wolffd@0:                 id(w) = idctr;
wolffd@0:             end
wolffd@0:             idctr = idctr + 1;
wolffd@0:         end
wolffd@0:         
wolffd@0:     end % function dfs 
wolffd@0: end
wolffd@0: 
wolffd@0: function [Gs, s, id] = connected_components(G, varargin)
wolffd@0:     % ---
wolffd@0:     % get all connected subgraphs:
wolffd@0:     % ---
wolffd@0:     
wolffd@0:     % make edge matrix undirected
wolffd@0:     G2 = DiGraph(Graph(G));
wolffd@0:     
wolffd@0:     [GsGraph, s, id] = strongly_connected_components(G2, varargin{:});
wolffd@0:     
wolffd@0:     % get the actual subgraps
wolffd@0:     
wolffd@0:     for i =1:numel(GsGraph)
wolffd@0:         Gs(i) = G.subgraph(GsGraph(i).nodes);
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % creates new graph just containing the 
wolffd@0: % specified nodes and optionally specified edges
wolffd@0: % nodes can be specified using 
wolffd@0: % a. the binary G.V structure
wolffd@0: % b. or node indices 
wolffd@0: % ---
wolffd@0: function G2 = subgraph(G, V, E)
wolffd@0:     if nargin == 2
wolffd@0:         E = [];
wolffd@0:     end
wolffd@0:     
wolffd@0:     % ---
wolffd@0:     % create new graph as copy ofthe old
wolffd@0:     % ---
wolffd@0:     
wolffd@0:     G2 = feval(class(G), G);
wolffd@0:     
wolffd@0:     % ---
wolffd@0:     % reset nodes and edges
wolffd@0:     % NOTE: we determine the input format first
wolffd@0:     % ---
wolffd@0:     if (max(V) == 1 && numel(V) > 1) || max(V) == 0
wolffd@0:         
wolffd@0:         G2.remove_node(find(~V));
wolffd@0:     else
wolffd@0:         
wolffd@0:         G2.remove_node(setdiff(1:numel(G.V), V));
wolffd@0:     end
wolffd@0:     if ~isempty(E)   
wolffd@0:         G2.E = E;
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % joins the information of graph G2 with
wolffd@0: % this GRaph, not duplicating nodes.
wolffd@0: % ---
wolffd@0: function add_graph(G, G2)
wolffd@0:     
wolffd@0:     % determine if symmetric  edges have to be added
wolffd@0:     clas = cat(1, class(G2), superclasses(G2));
wolffd@0:     if strcellfind(clas, 'Graph')
wolffd@0:         add_symm = 1;
wolffd@0:     else
wolffd@0:         add_symm = 0;
wolffd@0:     end
wolffd@0:             
wolffd@0:     % ---
wolffd@0:     % add all nodes and edges in G2
wolffd@0:     % ---
wolffd@0:     for V = G2.nodes();
wolffd@0:         
wolffd@0:         G.add_node(V);
wolffd@0:     end
wolffd@0:     
wolffd@0:     % ---
wolffd@0:     % NOTE / TODO: 
wolffd@0:     % this LABEL inheritance is far too expensive when 
wolffd@0:     % creating many copiesof the same graph
wolffd@0:     % Its also unnessecary for lots of them 
wolffd@0:     % except for debugging :(. 
wolffd@0:     % ---
wolffd@0:     G.Vlbl = cell(1,numel(G2.V));
wolffd@0:     G.Vlbl(G2.nodes()) = G2.label(G2.nodes());
wolffd@0:     
wolffd@0:     % ---
wolffd@0:     % get all edges in G2
wolffd@0:     % ---
wolffd@0:     [val, V1, V2] = edges(G2);
wolffd@0:     
wolffd@0:     for i = 1:numel(val)
wolffd@0:         
wolffd@0:         % ---
wolffd@0:         % add edges to graph
wolffd@0:         % ---
wolffd@0:         G.add_edge(V1(i), V2(i), val(i));
wolffd@0:         
wolffd@0:         if add_symm
wolffd@0:             % ---
wolffd@0:             % add symmetric edges to graph
wolffd@0:             % ---
wolffd@0:             G.add_edge(V2(i), V1(i), val(i));
wolffd@0:         end
wolffd@0:         
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % substracts the edges in G2 from
wolffd@0: % this GRaph, removing not connected nodes.
wolffd@0: % ---
wolffd@0: function Gout = minus(a, b)
wolffd@0:     Gout = feval(class(a),a);
wolffd@0:     Gout.remove_graph(b);
wolffd@0: end
wolffd@0: 
wolffd@0: function remove_graph(G, G2)
wolffd@0:     
wolffd@0:     % determine if symmetric  edges have to be added
wolffd@0:     clas = cat(1, class(G2), superclasses(G2));
wolffd@0:     if strcellfind(clas, 'Graph')
wolffd@0:         symm = 1;
wolffd@0:     else
wolffd@0:         symm = 0;
wolffd@0:     end
wolffd@0:     
wolffd@0:     % remove edges
wolffd@0:     [val, V1, V2] = G2.edges();
wolffd@0:     for j = 1:numel(val)
wolffd@0:         
wolffd@0:         % remove specified edges
wolffd@0:         G.remove_edge(V1(j), V2(j));
wolffd@0:         
wolffd@0:         % remove symmetric edges if subtracting symm graph
wolffd@0:         if symm
wolffd@0: 
wolffd@0:             G.remove_edge(V2(j), V1(j));
wolffd@0:         end
wolffd@0:     end
wolffd@0: 
wolffd@0:     % ---
wolffd@0:     % Note : we only remove nodes with no 
wolffd@0:     % remaining incoming edges
wolffd@0:     % ---
wolffd@0:     V = G2.nodes();
wolffd@0:     for j = 1:numel(V)
wolffd@0: 
wolffd@0:         if G.degree(V(j)) == 0
wolffd@0: 
wolffd@0:             G.remove_node(V(j));
wolffd@0:         end
wolffd@0:     end
wolffd@0:     
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % compact graph representation
wolffd@0: % ---
wolffd@0: function [E, labels] = compact(G)
wolffd@0:     
wolffd@0:     % ---
wolffd@0:     % get nodes and create a reverse index
wolffd@0:     % ---
wolffd@0:     Vidx = sparse(G.nodes,1,1:G.order());
wolffd@0:     [w, V1, V2]  = G.edges();
wolffd@0:     
wolffd@0:     % create compact adjacency matrix
wolffd@0:     E = sparse(Vidx(V1), Vidx(V2),w, G.order(), G.order());
wolffd@0:     
wolffd@0:     labels = G.label(G.nodes());
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % determines if Edges in G2 are the same as in G
wolffd@0: % ---
wolffd@0: function [out] = isSubgraph(G, G2)
wolffd@0:     
wolffd@0:     [val, V1, V2] = G2.edges();
wolffd@0:     validE = false(numel(V1), 1);
wolffd@0:     
wolffd@0:     i = 1;
wolffd@0:     while i <= numel(V1)
wolffd@0:         
wolffd@0:         % ---
wolffd@0:         % Test if edge exists in other graph
wolffd@0:         % ---
wolffd@0:         if G.edge(V1(i),V2(i)) == val(i) 
wolffd@0:            out = 0;
wolffd@0:            return;
wolffd@0:         end
wolffd@0:         i = i +1 ;
wolffd@0:     end
wolffd@0:     
wolffd@0:     % ---
wolffd@0:     % Test labels
wolffd@0:     % ---
wolffd@0:     if ~isempty(G.Vlbl) && ~isempty(G2.Vlbl)
wolffd@0:         
wolffd@0:         V = G2.nodes();
wolffd@0:         i = 1;
wolffd@0:         while i <= numel(V)
wolffd@0:             if strcmp(G.label(V(i)), G2.label(V(i))) ~= 0
wolffd@0:                 out = 0;
wolffd@0:                 return;
wolffd@0:             end
wolffd@0:             i = i + 1;
wolffd@0:         end 
wolffd@0:     end
wolffd@0:     
wolffd@0:     out = 1;
wolffd@0: end
wolffd@0: 
wolffd@0: % ---
wolffd@0: % Visualise the Similarity Graph
wolffd@0: % ---
wolffd@0: function visualise(G)
wolffd@0: 
wolffd@0:    % ---
wolffd@0:    % get colormap for edge weights
wolffd@0:    % ---
wolffd@0:    cmap = jet(100);
wolffd@0:    
wolffd@0:    % ---
wolffd@0:    % NOTE: we now get the weight colors for all edges
wolffd@0:    % get maximum weight and all edges
wolffd@0:    % ---
wolffd@0:    [colors, V1, V2]  = G.edges();
wolffd@0:    w = G.max_weight();
wolffd@0:    
wolffd@0:    % normalise colors
wolffd@0:    colors = max(1,round((colors / w) * 100));
wolffd@0: 
wolffd@0:    % get node labels
wolffd@0:    V1lbl = G.label(V1);
wolffd@0:    V2lbl = G.label(V2);
wolffd@0: 
wolffd@0:    % ---
wolffd@0:    % compose edgecolor matrix
wolffd@0:    % ---
wolffd@0:    ec = cell(numel(V1), 3);
wolffd@0:    for i = 1:numel(V1)
wolffd@0:        ec(i,:) = {V1lbl{i}, V2lbl{i}, cmap(colors(i),:)};
wolffd@0:    end
wolffd@0:    
wolffd@0:    % ---
wolffd@0:    % For Performance Issues
wolffd@0:    % We get the compact Graph and
wolffd@0:    % !hope! for the labels to correspond to those above
wolffd@0:    % ---
wolffd@0:    [E, labels] = compact(G);
wolffd@0:    
wolffd@0:    % determine if symmetric Graph
wolffd@0:     clas = cat(1, class(G), superclasses(G));
wolffd@0:     if strcellfind(clas, 'Graph')
wolffd@0:         symm = 1;
wolffd@0:     else
wolffd@0:         symm = 0;
wolffd@0:     end
wolffd@0:    
wolffd@0:    graphViz4Matlab('-adjMat',E, ...
wolffd@0:     '-nodeLabels',labels, ...
wolffd@0:     '-edgeColors', ec, ...
wolffd@0:     '-undirected', symm); 
wolffd@0: end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: methods (Static)
wolffd@0:    function G = copy(Gin)
wolffd@0:        
wolffd@0:        if strcmp(class(Gin), 'DiGraph')
wolffd@0:            
wolffd@0:            G = DiGraph(Gin.V, Gin.E);
wolffd@0:            G.Vlbl = Gin.Vlbl;
wolffd@0:            
wolffd@0:        else
wolffd@0:            warning ('cannot copy graph, casting instead')
wolffd@0:            G = DiGraph.cast(Gin);
wolffd@0:        end
wolffd@0:         
wolffd@0:    end
wolffd@0:    
wolffd@0:    function G = cast(Gin)
wolffd@0:        
wolffd@0:         % ---
wolffd@0:         % this uses an imput grpaph
wolffd@0:         % to create a new digraph
wolffd@0:         % ---
wolffd@0:         G = DiGraph();
wolffd@0:         
wolffd@0:         % ---
wolffd@0:         % Add the input graph to the empty one
wolffd@0:         % ---
wolffd@0:         G.add_graph(Gin);
wolffd@0:         
wolffd@0:    end
wolffd@0: end
wolffd@0: end