/*
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;
}