deba@2025: #include <iostream>
deba@2025: #include <set>
deba@2025: #include <vector>
deba@2025: #include <iterator>
deba@2025: 
deba@2025: #include <cmath>
deba@2025: #include <cstdlib>
deba@2025: 
deba@2025: #include <lemon/smart_graph.h>
deba@2025: #include <lemon/min_cost_arborescence.h>
deba@2025: 
deba@2025: #include <lemon/graph_utils.h>
deba@2025: #include <lemon/time_measure.h>
deba@2025: 
deba@2025: #include <lemon/tolerance.h>
deba@2025: 
deba@2180: #include "test_tools.h"
deba@2180: 
deba@2025: using namespace lemon;
deba@2025: using namespace std;
deba@2025: 
deba@2180: const int n = 10;
deba@2180: const int e = 22;
deba@2180: 
deba@2180: int sourceNode = 0;
deba@2180: 
deba@2180: int sources[e] = {
deba@2180:   1, 0, 2, 4, 4, 3, 9, 8, 9, 8,
deba@2180:   4, 2, 0, 6, 4, 1, 7, 2, 8, 6,
deba@2180:   1, 0
deba@2180: };
deba@2180: 
deba@2180: int targets[e] = {
deba@2180:   8, 3, 1, 1, 4, 9, 8, 1, 8, 0,
deba@2180:   3, 2, 1, 3, 1, 1, 2, 6, 3, 9,
deba@2180:   1, 3
deba@2180: };
deba@2180: 
deba@2180: double costs[e] = {
deba@2180:   107.444, 70.3069, 46.0496, 28.3962, 91.4325,
deba@2180:   76.9443, 61.986, 39.3754, 74.9575, 39.3153,
deba@2180:   45.7094, 34.6184, 100.156, 95.726, 22.3429,
deba@2180:   31.587, 51.6972, 29.6773, 115.038, 32.4137,
deba@2180:   60.0038, 40.1237
deba@2180: };
deba@2180: 
deba@2180: 
deba@2180: 
deba@2180: int main() {
deba@2025:   srand(time(0));
deba@2025:   typedef SmartGraph Graph;
deba@2025:   GRAPH_TYPEDEFS(Graph);
deba@2025: 
deba@2025:   typedef Graph::EdgeMap<double> CostMap;
deba@2025: 
deba@2025:   Graph graph;
deba@2025:   CostMap cost(graph);
deba@2025:   vector<Node> nodes;
deba@2025:   
deba@2025:   for (int i = 0; i < n; ++i) {
deba@2025:     nodes.push_back(graph.addNode());
deba@2025:   }
deba@2025: 
deba@2025:   for (int i = 0; i < e; ++i) {
deba@2180:     Edge edge = graph.addEdge(nodes[sources[i]], nodes[targets[i]]);
deba@2180:     cost[edge] = costs[i];
deba@2025:   }
deba@2025: 
deba@2180:   Node source = nodes[sourceNode];
deba@2025: 
deba@2025:   MinCostArborescence<Graph, CostMap> mca(graph, cost);
deba@2025:   mca.run(source);
deba@2025: 
deba@2025:   vector<pair<double, set<Node> > > dualSolution(mca.dualSize());
deba@2025: 
deba@2025:   for (int i = 0; i < mca.dualSize(); ++i) {
deba@2025:     dualSolution[i].first = mca.dualValue(i);
deba@2025:     for (MinCostArborescence<Graph, CostMap>::DualIt it(mca, i); 
deba@2025:          it != INVALID; ++it) {
deba@2025:       dualSolution[i].second.insert(it);
deba@2025:     }
deba@2025:   }
deba@2025: 
deba@2025:   Tolerance<double> tolerance;
deba@2025: 
deba@2025:   for (EdgeIt it(graph); it != INVALID; ++it) {
deba@2025:     if (mca.reached(graph.source(it))) {
deba@2025:       double sum = 0.0;
deba@2025:       for (int i = 0; i < (int)dualSolution.size(); ++i) {
deba@2025:         if (dualSolution[i].second.find(graph.target(it)) 
deba@2025:             != dualSolution[i].second.end() &&
deba@2025:             dualSolution[i].second.find(graph.source(it)) 
deba@2025:             == dualSolution[i].second.end()) {
deba@2025:           sum += dualSolution[i].first;
deba@2025:         }
deba@2025:       }
deba@2025:       if (mca.arborescence(it)) {
deba@2180:         check(!tolerance.less(sum, cost[it]), "INVALID DUAL");
deba@2025:       }
deba@2180:       check(!tolerance.less(cost[it], sum), "INVALID DUAL");
deba@2025:     }
deba@2025:   }
deba@2025: 
deba@2025: 
deba@2180:   check(!tolerance.different(mca.dualValue(), mca.arborescenceValue()),
deba@2025:                "INVALID DUAL");
deba@2025: 
deba@2025: 
deba@2180:   check(mca.reached(source), "INVALID ARBORESCENCE");
deba@2025:   for (EdgeIt it(graph); it != INVALID; ++it) {
deba@2180:     check(!mca.reached(graph.source(it)) || 
deba@2025:                  mca.reached(graph.target(it)), "INVALID ARBORESCENCE");
deba@2025:   }
deba@2025: 
deba@2025:   for (NodeIt it(graph); it != INVALID; ++it) {
deba@2025:     if (!mca.reached(it)) continue;
deba@2025:     int cnt = 0;
deba@2025:     for (InEdgeIt jt(graph, it); jt != INVALID; ++jt) {
deba@2025:       if (mca.arborescence(jt)) {
deba@2025:         ++cnt;
deba@2025:       }
deba@2025:     }
deba@2180:     check((it == source ? cnt == 0 : cnt == 1), "INVALID ARBORESCENCE");
deba@2025:   }
deba@2025:   
deba@2025:   return 0;
deba@2025: }