#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>

using namespace lemon;

using namespace std;

int main(int argc, const char *argv[]) {

  srand(time(0));

  typedef SmartGraph Graph;

  GRAPH_TYPEDEFS(Graph);

  typedef Graph::EdgeMap<double> CostMap;

  const int n = argc > 1 ? atoi(argv[1]) : 100;

  const int e = argc > 2 ? atoi(argv[2]) : (int)(n * log(n));

  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) {

    int s = (int)(n * (double)rand() / (RAND_MAX + 1.0));

    int t = (int)(n * (double)rand() / (RAND_MAX + 1.0));

    double c = rand() / (1.0 + RAND_MAX) * 100.0 + 20.0;

    Edge edge = graph.addEdge(nodes[s], nodes[t]);

    cost[edge] = c;

  }

  Node source = nodes[(int)(n * (double)rand() / (RAND_MAX + 1.0))];

  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)) {

        LEMON_ASSERT(!tolerance.less(sum, cost[it]), "INVALID DUAL");

      }

      LEMON_ASSERT(!tolerance.less(cost[it], sum), "INVALID DUAL");

    }

  }

  LEMON_ASSERT(!tolerance.different(mca.dualValue(), mca.arborescenceValue()),

               "INVALID DUAL");

  LEMON_ASSERT(mca.reached(source), "INVALID ARBORESCENCE");

  for (EdgeIt it(graph); it != INVALID; ++it) {

    LEMON_ASSERT(!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;

      }

    }

    LEMON_ASSERT((it == source ? cnt == 0 : cnt == 1), "INVALID ARBORESCENCE");

  }

  return 0;

}