lemon-0.x: comparison src/work/athos/mincostflow.h

equal deleted inserted replaced

-:97c30c04b2d4
+:336fbb1bba78
 typedef typename ResGraphType::Edge ResGraphEdge;
 class ModLengthMap {
 //typedef typename ResGraphType::template NodeMap<Length> NodeMap;
 typedef typename Graph::template NodeMap<Length> NodeMap;
-const ResGraphType& G;
+const ResGraphType& res_graph;
 //      const EdgeIntMap& rev;
 const LengthMap &ol;
 const NodeMap &pot;
 public :
 typedef typename LengthMap::KeyType KeyType;
 typedef typename LengthMap::ValueType ValueType;
 ValueType operator[](typename ResGraphType::Edge e) const {
-	if (G.forward(e))
+	if (res_graph.forward(e))
-	  return  ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
+	  return  ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);
 	else
-	  return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
+	  return -ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);
 }
-ModLengthMap(const ResGraphType& _G,
+ModLengthMap(const ResGraphType& _res_graph,
 		   const LengthMap &o,  const NodeMap &p) :
-	G(_G), /*rev(_rev),*/ ol(o), pot(p){};
+	res_graph(_res_graph), /*rev(_rev),*/ ol(o), pot(p){};
 };//ModLengthMap
 protected:
 //Input
-const Graph& G;
+const Graph& graph;
 const LengthMap& length;
 const SupplyDemandMap& supply_demand;//supply or demand of nodes
 //auxiliary variables
 public :
-MinCostFlow(Graph& _G, LengthMap& _length, SupplyDemandMap& _supply_demand) : G(_G),
+MinCostFlow(Graph& _graph, LengthMap& _length, SupplyDemandMap& _supply_demand) : graph(_graph),
-length(_length), supply_demand(_supply_demand), flow(_G), potential(_G){ }
+length(_length), supply_demand(_supply_demand), flow(_graph), potential(_graph){ }
 ///Runs the algorithm.
 ///Runs the algorithm.
 ///\todo May be it does make sense to be able to start with a nonzero
 /// feasible primal-dual solution pair as well.
-int run() {
+void run() {
 //Resetting variables from previous runs
 //total_length = 0;
 typedef typename Graph::template NodeMap<int> HeapMap;
 typedef Heap< Node, SupplyDemand, typename Graph::template NodeMap<int>,
 	std::greater<SupplyDemand> > 	HeapType;
 //A heap for the excess nodes
-HeapMap excess_nodes_map(G,-1);
+HeapMap excess_nodes_map(graph,-1);
 HeapType excess_nodes(excess_nodes_map);
 //A heap for the deficit nodes
-HeapMap deficit_nodes_map(G,-1);
+HeapMap deficit_nodes_map(graph,-1);
 HeapType deficit_nodes(deficit_nodes_map);
 //A container to store nonabundant arcs
 list<Edge> nonabundant_arcs;
-FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
+FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
 	flow.set(e,0);
 	nonabundant_arcs.push_back(e);
 }
 //Initial value for delta
 SupplyDemand delta = 0;
 typedef UnionFindEnum<Node, Graph::template NodeMap> UFE;
 //A union-find structure to store the abundant components
-UFE::MapType abund_comp_map(G);
+UFE::MapType abund_comp_map(graph);
 UFE abundant_components(abund_comp_map);
-FOR_EACH_LOC(typename Graph::NodeIt, n, G){
+FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
 	excess_deficit.set(n,supply_demand[n]);
 	//A supply node
 	if (excess_deficit[n] > 0){
 	  excess_nodes.push(n,excess_deficit[n]);
 	}
 //It'll be allright as an initial value, though this value
 //can be the maximum deficit here
 SupplyDemand max_excess = delta;
+//ConstMap<Edge,SupplyDemand> ConstEdgeMap;
 //We need a residual graph which is uncapacitated
-ResGraphType res_graph(G, flow);
+ResGraphType res_graph(graph, flow);
+//An EdgeMap to tell which arcs are abundant
+template typename Graph::EdgeMap<bool> abundant_arcs(graph);
+//Let's construct the sugraph consisting only of the abundant edges
+typedef ConstMap< typename Graph::Node, bool > ConstNodeMap;
+ConstNodeMap const_true_map(true);
+typedef SubGraphWrapper< Graph, ConstNodeMap,
+	 template typename Graph::EdgeMap<bool> >
+	AbundantGraph;
+AbundantGraph abundant_graph(graph, const_true_map, abundant_arcs );
+//Let's construct the residual graph for the abundant graph
+typedef ResGraphWrapper<const AbundantGraph,int,CapacityMap,EdgeIntMap>
+	ResAbGraph;
+//Again uncapacitated
+ResAbGraph res_ab_graph(abundant_graph, flow);
+//We need things for the bfs
+typename ResAbGraph::NodeMap<bool> bfs_reached(res_ab_graph);
+typename ResAbGraph::NodeMap<typename ResAbGraph::Edge>
+	bfs_pred(res_ab_graph);
+NullMap<typename ResAbGraph::Node, int> bfs_dist_dummy(res_ab_graph);
+//We want to run bfs-es (more) on this graph 'res_ab_graph'
+Bfs < ResAbGraph ,
+	typename ResAbGraph::NodeMap<bool>,
+	typename ResAbGraph::NodeMap<typename ResAbGraph::Edge>,
+	NullMap<typename ResAbGraph::Node, int> >
+	bfs(res_ab_graph, bfs_reached, bfs_pred, bfs_dist_dummy);
 ModLengthMap mod_length(res_graph, length, potential);
 Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length);
 	    if (flow[i]>=buf){
 	      Node a = abundant_components.find(res_graph.head(i));
 	      Node b = abundant_components.find(res_graph.tail(i));
 	      //Merge
 	      if (a != b){
-		//Find path and augment
-		//!!!
-		//Find path and augment
 		abundant_components.join(a,b);
+		//We want to push the smaller
+		//Which has greater absolut value excess/deficit
+		Node root=(abs(excess_deficit[a])>abs(excess_deficit[b]))?a:b;
+		//Which is the other
+		Node non_root = ( a == root ) ? b : a ;
+		abundant_components.makeRep(root);
+		SupplyDemand qty_to_augment = abs(excess_deficit[non_root]);
+		//Push the positive value
+		if (excess_deficit[non_root] < 0)
+		  swap(root, non_root);
+		//If the non_root node has excess/deficit at all
+		if (qty_to_augment>0){
+		  //Find path and augment
+		  bfs.run(non_root);
+		  //root should be reached
+		  //Augmenting on the found path
+		  Node n=root;
+		  ResGraphEdge e;
+		  while (n!=non_root){
+		    e = bfs_pred(n);
+		    n = res_graph.tail(e);
+		    res_graph.augment(e,qty_to_augment);
+		  }
+		  //We know that non_root had positive excess
+		  excess_nodes[non_root] -= qty_to_augment;
+		  //But what about root node
+		  //It might have been positive and so became larger
+		  if (excess_deficit[root]>0){
+		    excess_nodes[root] += qty_to_augment;
+		  }
+		  else{
+		    //Or negative but not turned into positive
+		    deficit_nodes[root] -= qty_to_augment;
+		  }
+		  //Update the excess_deficit map
+		  excess_deficit[non_root] -= qty_to_augment;
+		  excess_deficit[root] += qty_to_augment;
+		}
 	      }
 	      //What happens to i?
-	      nonabundant_arcs.erase(i);
+	      //Marci and Zsolt says I shouldn't do such things
+	      nonabundant_arcs.erase(i++);
+	      abundant_arcs[i] = true;
 	    }
 	    else
 	      ++i;
 	  }
 	}
 	Node s = excess_nodes.top();
 	SupplyDemand max_excess = excess_nodes[s];
 	Node t = deficit_nodes.top();
-	if (max_excess < dificit_nodes[t]){
+	if (max_excess < deficit_nodes[t]){
-	  max_excess = dificit_nodes[t];
+	  max_excess = deficit_nodes[t];
 	}
-	while(max_excess > 0){
+	while(max_excess > (n-1)*delta/n){
 	  //s es t valasztasa
 	  //Dijkstra part
 	  dijkstra.run(s);
 	  /*We know from theory that t can be reached
 	  if (!dijkstra.reached(t)){
 	    //There are no k paths from s to t
 	    break;
 	  };
 	      total_length += length[e];
 	    else
 	      total_length -= length[e];
 	    */
 	  }
+	  //Update the excess_deficit map
+	  excess_deficit[s] -= delta;
+	  excess_deficit[t] += delta;
 	  //Update the excess_nodes heap
 	  if (delta >= excess_nodes[s]){
 	    if (delta > excess_nodes[s])
 	      deficit_nodes.push(s,delta - excess_nodes[s]);
 	  }
 	  else{
 	    deficit_nodes[t] -= delta;
 	  }
 	  //Dijkstra part ends here
+	  //Choose s and t again
+	  s = excess_nodes.top();
+	  max_excess = excess_nodes[s];
+	  t = deficit_nodes.top();
+	  if (max_excess < deficit_nodes[t]){
+	    max_excess = deficit_nodes[t];
+	  }
 	}
 	/*
 	 * End of the delta scaling phase
 	*/
 	delta = delta / 2;
 	/*This is not necessary here
 	//Update the max_excess
 	max_excess = 0;
-	FOR_EACH_LOC(typename Graph::NodeIt, n, G){
+	FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
 	  if (max_excess < excess_deficit[n]){
 	    max_excess = excess_deficit[n];
 	  }
 	}
 	*/
 ///
 ///\todo Is this OK here?
 bool checkComplementarySlackness(){
 Length mod_pot;
 Length fl_e;
-FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
+FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
 	//C^{\Pi}_{i,j}
-	mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)];
+	mod_pot = length[e]-potential[graph.head(e)]+potential[graph.tail(e)];
 	fl_e = flow[e];
 	//	std::cout << fl_e << std::endl;
 	if (0<fl_e && fl_e<capacity[e]){
 	  if (mod_pot != 0)
 	    return false;

changeset 660	edb42cb9d352
parent 657	531fc5f575ef
child 661	d306e777117e