#include <iostream>
#include <set>
#include <vector>
#include <iterator>

#include <cmath>
#include <cstdlib>

#include <lemon/smart_graph.h>
#include <lemon/min_cost_arborescence.h>

#include <lemon/graph_utils.h>
#include <lemon/time_measure.h>

#include <lemon/tolerance.h>

#include "test_tools.h"

using namespace lemon;
using namespace std;

const int n = 10;
const int e = 22;

int sourceNode = 0;

int sources[e] = {
  1, 0, 2, 4, 4, 3, 9, 8, 9, 8,
  4, 2, 0, 6, 4, 1, 7, 2, 8, 6,
  1, 0
};

int targets[e] = {
  8, 3, 1, 1, 4, 9, 8, 1, 8, 0,
  3, 2, 1, 3, 1, 1, 2, 6, 3, 9,
  1, 3
};

double costs[e] = {
  107.444, 70.3069, 46.0496, 28.3962, 91.4325,
  76.9443, 61.986, 39.3754, 74.9575, 39.3153,
  45.7094, 34.6184, 100.156, 95.726, 22.3429,
  31.587, 51.6972, 29.6773, 115.038, 32.4137,
  60.0038, 40.1237
};



int main() {
  typedef SmartGraph Graph;
  GRAPH_TYPEDEFS(Graph);

  typedef Graph::EdgeMap<double> CostMap;

  Graph graph;
  CostMap cost(graph);
  vector<Node> nodes;
  
  for (int i = 0; i < n; ++i) {
    nodes.push_back(graph.addNode());
  }

  for (int i = 0; i < e; ++i) {
    Edge edge = graph.addEdge(nodes[sources[i]], nodes[targets[i]]);
    cost[edge] = costs[i];
  }

  Node source = nodes[sourceNode];

  MinCostArborescence<Graph, CostMap> mca(graph, cost);
  mca.run(source);

  vector<pair<double, set<Node> > > dualSolution(mca.dualSize());

  for (int i = 0; i < mca.dualSize(); ++i) {
    dualSolution[i].first = mca.dualValue(i);
    for (MinCostArborescence<Graph, CostMap>::DualIt it(mca, i); 
         it != INVALID; ++it) {
      dualSolution[i].second.insert(it);
    }
  }

  Tolerance<double> tolerance;

  for (EdgeIt it(graph); it != INVALID; ++it) {
    if (mca.reached(graph.source(it))) {
      double sum = 0.0;
      for (int i = 0; i < (int)dualSolution.size(); ++i) {
        if (dualSolution[i].second.find(graph.target(it)) 
            != dualSolution[i].second.end() &&
            dualSolution[i].second.find(graph.source(it)) 
            == dualSolution[i].second.end()) {
          sum += dualSolution[i].first;
        }
      }
      if (mca.arborescence(it)) {
        check(!tolerance.less(sum, cost[it]), "INVALID DUAL");
      }
      check(!tolerance.less(cost[it], sum), "INVALID DUAL");
    }
  }


  check(!tolerance.different(mca.dualValue(), mca.arborescenceValue()),
               "INVALID DUAL");


  check(mca.reached(source), "INVALID ARBORESCENCE");
  for (EdgeIt it(graph); it != INVALID; ++it) {
    check(!mca.reached(graph.source(it)) || 
                 mca.reached(graph.target(it)), "INVALID ARBORESCENCE");
  }

  for (NodeIt it(graph); it != INVALID; ++it) {
    if (!mca.reached(it)) continue;
    int cnt = 0;
    for (InEdgeIt jt(graph, it); jt != INVALID; ++jt) {
      if (mca.arborescence(jt)) {
        ++cnt;
      }
    }
    check((it == source ? cnt == 0 : cnt == 1), "INVALID ARBORESCENCE");
  }
  
  return 0;
}