new file mode 100644

@@ -0,0 +1,476 @@

+/*

+  max_weight_branch.c

+

+  Version : 1.1

+  Author  : Niko Beerenwinkel

+  

+  (c) Max-Planck-Institut fuer Informatik, 2004

+*/

+

+

+#include "include/max_weight_branch.h"

+

+

+// global variables:

+

+static map<node,string> Mut;  // node labels

+static map<edge,double> c;    // edge weights

+

+static array< array<node> > cycle_contr;

+// cycle_contr[i][k] == v_C <==> the k-th cycle C in B[i] has been contracted to v_C in G[i+1]

+

+static array< array<edge> > cheapest_edge;

+// cheapest_edge[i][k] is some cheapest edge on the k-th cycle in B[i]

+

+static map<edge,edge> Phi;

+// Phi[e'] == e  <==> e in G[i] has been contracted to e' in G[i+1]

+

+static array< array< map<edge,edge> > > alpha;  

+// e = (z, y), z not in C, y in C  <==>  alpha[e] == pred of y in C 

+

+

+

+// functions:

+

+int compare_weights(const edge& e1, const edge& e2)

+{

+  if (c[e1] < c[e2])  return -1;

+  if (c[e1] > c[e2])  return  1;

+  return 0;

+}

+

+

+

+double BRANCHING_WEIGHT(list<edge>& B, edge_array<double>& weight)

+{

+  // total weight of a branching

+  

+  double total_weight = 0;

+  edge e;

+

+  forall(e, B)

+    total_weight += weight[e];

+

+  return total_weight;

+}

+

+

+

+void UNCOVER_BRANCHING(graph& G, list<edge>& B)

+{

+  // turn G into the branching defined by B (B \subset E(G))

+

+  edge e;

+

+  edge_set E_B(G);  // edge set of B

+

+  forall(e, B)

+    E_B.insert(e);

+

+  edge_set T_D(G); //edges to delete

+

+  forall_edges(e, G)

+    if (! E_B.member(e))

+	T_D.insert(e);

+      //G.del_edge(e);

+

+  for(edge_set::iterator it = T_D.begin(); it != T_D.end(); ++it)

+      G.del_edge(*it);

+

+}

+

+

+

+void DOT(graph& G, map<node,string>& mut, map<edge,double>& cond_prob, char *filestem)

+{

+  // produces output in dot format (graphviz package)

+  // see http://www.research.att.com/sw/tools/graphviz/

+

+  node v;  edge e;

+

+  char filename[128];

+  sprintf(filename, "%s.dot", filestem);

+  

+  std::ofstream dotout(filename);

+  dotout << "digraph MWB {" << std::endl << std::endl;

+  

+  forall_nodes(v, G)

+    {

+      dotout << "\t \"" << v << "\"";

+      dotout << " [ label=\"" << mut[v] << "\", shape=\"plaintext\", height=\"0.3\", fontsize=\"12\", style=\"filled\", fillcolor=\"white\" ];" << std::endl;

+    }

+  dotout << std::endl;

+

+  forall_edges(e, G)

+    {

+      node s = G.source(e);

+      node t = G.target(e);

+      dotout.precision(2);

+      dotout << std::showpoint;

+      dotout << "\t \"" << s << "\" -> \"" << t << "\"";

+      dotout << " [ fontsize=\"10\", label=\"" << cond_prob[e] << "\" ];" << std::endl;

+    }

+

+  dotout << "}" << std::endl;

+  dotout.close();

+}

+

+

+double branching_weight_intern(list<edge>& B)

+{

+  // total weight of a branching

+  

+  double w = 0;

+

+  edge e;

+  forall(e, B)

+    w += c[e];

+  

+  return w;

+}

+

+

+

+list<edge> max_weight_subgraph_indeg_le_1(graph& G)

+{

+  // find maximum weight subgraph of G with indeg <= 1 for all v in V(G)

+  // a list of edges in G of this subgraph is returned

+  

+  node v;

+  list<edge> B;

+

+  forall_nodes(v, G)  // choose heaviest inedge

+    {

+      list<edge> L = G.in_edges(v);

+      if (! L.empty())

+	B.append(L.contents(L.max(*compare_weights)));

+    }

+

+  return B;

+}

+

+

+array< list<edge> > all_cycles(graph& G, list<edge>& B)

+{

+  // find all cycles in B

+

+  edge e;

+

+  array< list<edge> > CA;  // array of cycles

+  int k = -1;              // cycle index

+

+  // construct subgraph:

+  GRAPH<node,edge> G_B;  // the subgraph of G induced by B

+  map<node,node> proj;   // proj : V(G) -->> V(G_B)

+

+  forall(e, B)

+    {

+      node s = G.source(e);

+      if (! proj.defined(s))

+	{

+	  node s_G = G_B.new_node();

+	  G_B[s_G] = s;  // needed?

+	  proj[s] = s_G;

+	}

+      node t = G.target(e);

+      if (! proj.defined(t))

+	{

+	  node t_G = G_B.new_node();

+	  G_B[t_G] = t;  // needed?

+	  proj[t] = t_G;

+	}

+      G_B[G_B.new_edge(proj[s], proj[t])] = e;

+    }

+  

+  // identify cycles in G_B:

+  list<edge> CEL;  // cycle edge list, 

+  // will contain one edge (in G_B) from each cycle after Is_Acyclic call

+

+  if (! Is_Acyclic(G_B, CEL))

+    {

+      edge first_e;

+      forall(first_e, CEL)

+	{

+ 	  list<edge> C;  // a cycle is stored as the list of its edges in G

+	  

+	  // iterate over cycle following inedges:

+	  node s = G_B.target(first_e);

+	  node first_v = s;  // starting node;

+	  edge e = NULL; //edge e = nil;

+	  while (s != first_v || e == NULL)//nil)

+	    {

+	      e = G_B.in_edges(s).head();

+	      s = G_B.source(e);  // move on

+	      C.push(G_B[e]);  // store corresponding edge in G (!)

+	      // the order in C follows the edge directions

+	    }

+

+	  k++;  CA.resize(k+1);

+ 	  CA[k] = C;

+ 	}

+    }

+  

+  return CA;

+}

+

+

+

+edge predecessor_in_cycle(node y, list<edge>& C)

+{

+  // predecessor edge of y in C

+  

+  edge e;

+  forall(e, C)

+    {

+      if (target(e) == y)

+	return e;

+    }

+

+  return NULL; //nil;

+}

+

+

+void contract_cycle_edge(edge& e, double alpha_cost, double cheapest_cost, GRAPH<node,edge>& H, node& v_C, map<edge,edge>& edge_contr, node_set& VC)

+{

+  // note:  e \in E(H)

+

+  node s = H.source(e);  // \in G(H)

+  node t = H.target(e);  // \in G(H)

+

+  if (VC.member(s) && (! VC.member(t)))  // e leaves C

+    {

+      edge e_tilde = H.new_edge(v_C, t);  // new edge

+      c[e_tilde] = c[e];

+      edge_contr[H[e]] = e_tilde;

+

+      H[e_tilde] = H[e];  // update

+      Phi[e_tilde] = H[e];

+

+      H.del_edge(e);

+    }

+

+  if ((! VC.member(s)) && VC.member(t))  // e enters C

+    {

+      edge e_prime = H.new_edge(s, v_C);  // new edge

+      c[e_prime] = c[e] - alpha_cost + cheapest_cost;

+      edge_contr[H[e]] = e_prime;

+

+      H[e_prime] = H[e];  // update

+      Phi[e_prime] = H[e];

+

+      H.del_edge(e);

+    }

+

+  if (VC.member(s) && VC.member(t))  // e is in C

+    H.del_edge(e);

+  

+}

+

+

+

+string contract_cycle_node(node& v, GRAPH<node,edge>& H, node& v_C, map<node,node>& node_contr, node_set& VC)

+{

+  string label;

+

+  if (VC.member(v))  // v is on C

+    {

+      label += Mut[v];

+      H.del_node(v);

+      node_contr[H[v]] = v_C;

+    }

+  

+  return label;

+}

+

+

+

+void contract_cycle(int i, array< list<edge> >& CA, int k, GRAPH<node,edge>& H, map<node,node>& node_contr, map<edge,edge>& edge_contr)

+{

+  node v;  edge e;

+

+  list<edge> C = CA[k];  // cycle to be contracted (indexed by k)

+  

+  // node set in H indiced by cycle (in G):

+  node_set VC(H);  // V(C) \subset H

+  forall(e, C)

+    {

+      VC.insert(H.source(edge_contr[e]));

+      VC.insert(H.target(edge_contr[e]));

+    }

+  

+  // some cheapest edge on C

+  cheapest_edge.resize(i+1);  cheapest_edge[i].resize(k+1);

+  cheapest_edge[i][k] = C.contents(C.min(*compare_weights));

+  

+  // new node v_C for contracted cycle C

+  node v_C = H.new_node();

+  cycle_contr.resize(i+1);  cycle_contr[i].resize(k+1);

+  cycle_contr[i][k] = v_C;  // in G[i+1]

+

+  // contract edges in H

+  list<edge> EL = H.all_edges(); 

+  forall(e, EL)

+    {

+      alpha.resize(i+1);  alpha[i].resize(k+1);

+      alpha[i][k][H[e]] = predecessor_in_cycle(target(H[e]), C);  // predecessor edge of target(e) in C

+      contract_cycle_edge(e, c[alpha[i][k][H[e]]], c[cheapest_edge[i][k]], H, v_C, edge_contr, VC);

+    }

+  

+  // contract nodes in H

+  string new_label("(");

+  list<node> NL = H.all_nodes();

+  forall(v, NL)

+    new_label += contract_cycle_node(v, H, v_C, node_contr, VC);

+  Mut[v_C] = new_label + ")";

+

+}

+

+

+

+void contract_all_cycles(graph& G, int i, array< list<edge> > CL, GRAPH<node,edge>& H)

+{

+  // construct the contraction (G, c) --> (H, c)

+

+  node v;  edge e;

+  

+  // {node, edge}_contr : (G, c) -->> (H, c')

+  map<node,node> node_contr;  // node contraction map

+  map<edge,edge> edge_contr;  // edge contraction map

+

+  // initialize as identity map

+  forall_nodes(v, G)

+    {

+      node_contr[v] = H.new_node();

+      H[node_contr[v]] = v;

+      Mut[node_contr[v]] = Mut[v];  // labels

+    } 

+  forall_edges(e, G)

+    {

+      edge_contr[e] = H.new_edge(node_contr[G.source(e)], node_contr[G.target(e)]);

+      H[edge_contr[e]] = e;

+      c[edge_contr[e]] = c[e];  // weights

+    }

+

+  // contract cycles in H

+  for (int k=CL.low(); k<=CL.high(); k++)

+    {

+      contract_cycle(i, CL, k, H, node_contr, edge_contr);

+      list<edge> empty_list;

+    }

+}

+

+

+

+void reconstruct_branching(int i, GRAPH<node,edge>& H, list<edge>& B_H, list<edge>& B_G, array< list<edge> >& CA)

+{

+  //  reconstruct branching B_G from the contraction

+  //  (G, B_G==B[i-1]) --> (H==G[i], B_H==B[i]) 

+  //  CA is the array of cycles in G

+

+  edge e;

+  

+  B_G.clear();

+  

+  node_set CCS(H);  // contracted cycles in H

+  for (int k=CA.low(); k<=CA.high(); k++)

+    CCS.insert(cycle_contr[i-1][k]);

+  

+  // E(B_G) := E(B_H) \ {e_prime=(z, v_C) | z \in V(B_H), C cycle in B_G}

+  forall(e, B_H)

+    if (! CCS.member(H.target(e))) 

+      B_G.append(H[e]);

+ 

+  // check all cycles in B_G

+  for (int k=CA.high(); k>=CA.low(); k--)

+    {

+      list<edge> C = CA[k];  // cycle edges

+      

+      // look for e'==(z, v_C) in B_H

+      edge e_prime = NULL; //nil;

+      forall(e, B_H)

+	if (H.target(e) == cycle_contr[i-1][k])  

+	  e_prime = e;

+      

+      if (e_prime != NULL) //nil)  // e' == (z, v_C) in B_H

+	{

+	  B_G.append(Phi[e_prime]);

+	  forall(e, C)

+	    if (e != alpha[i-1][k][Phi[e_prime]])

+	      B_G.append(e);

+	}

+      else  // no edge (z, v_C) in B_H

+	{

+	  forall(e, C)

+	    if (e != cheapest_edge[i-1][k])

+	      B_G.append(e);

+	}

+      

+    }

+}

+

+

+

+list<edge> MAX_WEIGHT_BRANCHING(graph& G0, map<node,string>& node_label, edge_array<double>& weight)

+{

+  // Compute a maximum weight branching of the weighted digraph (G, weight)

+  // The list of edges of the branching is returned.

+  // The algorithm with all notation is taken from the book

+  // Korte, Vygen: Combinatorial Optimization, Springer 2000, pp. 121-125

+

+  node v;  edge e;

+

+  // step 0: data structures:

+  array< GRAPH<node,edge> > G(1);  // G[i+1] is subgraph of G[i]

+  array< list<edge> > B(1);  // B[i] \subset E(G[i])

+  array< array< list<edge> > > CAA(1);  // array of cycle arrays, each cycle is a subset of B[i]

+  

+  // step 1: initialization.

+  CopyGraph(G[0], G0);  

+  forall_edges(e, G[0])  // edge weights

+    c[e] = weight[G[0][e]];

+

+  forall_nodes(v, G[0])  // node labels

+    Mut[v] = node_label[G[0][v]];

+

+  // step 2: construct greedy subgraph B[i].

+  B[0] = max_weight_subgraph_indeg_le_1(G[0]);

+

+  // step 3: for all cycles in B[i].

+  CAA[0] = all_cycles(G[0], B[0]);

+

+  int i = 0;

+  while (CAA[i].size() > 0)

+    {

+      // step 4: construct (G[i+1], c[i+1]) from (G[i], c[i]).

+      G.resize(i+2);  B.resize(i+2);  CAA.resize(i+2);

+

+      contract_all_cycles(G[i], i, CAA[i], G[i+1]);

+

+      i++;

+      B[i] = max_weight_subgraph_indeg_le_1(G[i]);  // (step 2)

+

+      CAA[i] = all_cycles(G[i], B[i]);  // (step 3)

+    }

+  

+  // step 5+6: reconstruct the maximum branching B

+  while (i > 0)

+    {

+      reconstruct_branching(i, G[i], B[i], B[i-1], CAA[i-1]);

+      i--;

+    }

+

+  list<edge> B0;

+  forall(e, B[i]){

+    B0.append(G[i][e]);

+  }

+

+  //clear the global static variables

+  Mut.clear();

+  c.clear();

+  cycle_contr.clear();

+  cheapest_edge.clear();

+  Phi.clear();

+  alpha.clear();

+

+  return B0;

+}