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