112    MinCostFlow(Graph& _graph, CostMap& _cost, SupplyDemandMap& _supply_demand):  | 
   110    MinCostFlow(Graph& _graph, CostMap& _cost, SupplyDemandMap& _supply_demand):  | 
   113      graph(_graph),   | 
   111      graph(_graph),   | 
   114      cost(_cost),   | 
   112      cost(_cost),   | 
   115      supply_demand(_supply_demand),   | 
   113      supply_demand(_supply_demand),   | 
   116      flow(_graph),   | 
   114      flow(_graph),   | 
   117      potential(_graph),  | 
   115      potential(_graph){ } | 
   118      excess_deficit(_graph){ } | 
         | 
   119   | 
   116   | 
   120       | 
   117       | 
   121     ///Runs the algorithm.  | 
   118     ///Runs the algorithm.  | 
   122   | 
   119   | 
   123     ///Runs the algorithm.  | 
   120     ///Runs the algorithm.  | 
   124   | 
   121   | 
   125     ///\todo May be it does make sense to be able to start with a nonzero   | 
   122     ///\todo May be it does make sense to be able to start with a nonzero   | 
   126     /// feasible primal-dual solution pair as well.  | 
   123     /// feasible primal-dual solution pair as well.  | 
   127     void run() { | 
   124     void run() { | 
   128   | 
   125   | 
         | 
   126       //To store excess-deficit values  | 
         | 
   127       SupplyDemandMap excess_deficit(graph);  | 
         | 
   128   | 
   129       //Resetting variables from previous runs  | 
   129       //Resetting variables from previous runs  | 
   130       //total_cost = 0;  | 
   130       //total_cost = 0;  | 
         | 
   131   | 
   131   | 
   132   | 
   132       typedef typename Graph::template NodeMap<int> HeapMap;  | 
   133       typedef typename Graph::template NodeMap<int> HeapMap;  | 
   133       typedef BinHeap< Node, SupplyDemand, typename Graph::template NodeMap<int>,  | 
   134       typedef BinHeap< Node, SupplyDemand, typename Graph::template NodeMap<int>,  | 
   134 	std::greater<SupplyDemand> > 	HeapType;  | 
   135 	std::greater<SupplyDemand> > 	HeapType;  | 
   135   | 
   136   | 
   194       //An EdgeMap to tell which arcs are abundant  | 
   195       //An EdgeMap to tell which arcs are abundant  | 
   195       typename Graph::template EdgeMap<bool> abundant_arcs(graph);  | 
   196       typename Graph::template EdgeMap<bool> abundant_arcs(graph);  | 
   196   | 
   197   | 
   197       //Let's construct the sugraph consisting only of the abundant edges  | 
   198       //Let's construct the sugraph consisting only of the abundant edges  | 
   198       typedef ConstMap< typename Graph::Node, bool > ConstNodeMap;  | 
   199       typedef ConstMap< typename Graph::Node, bool > ConstNodeMap;  | 
         | 
   200   | 
   199       ConstNodeMap const_true_map(true);  | 
   201       ConstNodeMap const_true_map(true);  | 
   200       typedef SubGraphWrapper< const Graph, ConstNodeMap,   | 
   202       typedef SubGraphWrapper< const Graph, ConstNodeMap,   | 
   201 	 typename Graph::template EdgeMap<bool> >   | 
   203 	 typename Graph::template EdgeMap<bool> >   | 
   202 	AbundantGraph;  | 
   204 	AbundantGraph;  | 
   203       AbundantGraph abundant_graph(graph, const_true_map, abundant_arcs );  | 
   205       AbundantGraph abundant_graph(graph, const_true_map, abundant_arcs );  | 
   368 	  excess_deficit[s] -= delta;  | 
   370 	  excess_deficit[s] -= delta;  | 
   369 	  excess_deficit[t] += delta;  | 
   371 	  excess_deficit[t] += delta;  | 
   370 	    | 
   372 	    | 
   371   | 
   373   | 
   372 	  //Update the excess_nodes heap  | 
   374 	  //Update the excess_nodes heap  | 
   373 	  if (delta >= excess_nodes[s]){ | 
   375 	  if (delta > excess_nodes[s]){ | 
   374 	    if (delta > excess_nodes[s])  | 
   376 	    if (delta > excess_nodes[s])  | 
   375 	      deficit_nodes.push(s,delta - excess_nodes[s]);  | 
   377 	      deficit_nodes.push(s,delta - excess_nodes[s]);  | 
   376 	    excess_nodes.pop();  | 
   378 	    excess_nodes.pop();  | 
   377 	      | 
   379 	      | 
   378 	  }   | 
   380 	  }   | 
   379 	  else{ | 
   381 	  else{ | 
   380 	    excess_nodes.set(s, excess_nodes[s] - delta);  | 
   382 	    excess_nodes.set(s, excess_nodes[s] - delta);  | 
   381 	  }  | 
   383 	  }  | 
   382 	  //Update the deficit_nodes heap  | 
   384 	  //Update the deficit_nodes heap  | 
   383 	  if (delta >= deficit_nodes[t]){ | 
   385 	  if (delta > deficit_nodes[t]){ | 
   384 	    if (delta > deficit_nodes[t])  | 
   386 	    if (delta > deficit_nodes[t])  | 
   385 	      excess_nodes.push(t,delta - deficit_nodes[t]);  | 
   387 	      excess_nodes.push(t,delta - deficit_nodes[t]);  | 
   386 	    deficit_nodes.pop();  | 
   388 	    deficit_nodes.pop();  | 
   387 	      | 
   389 	      | 
   388 	  }   | 
   390 	  }   | 
   442     /// \pre \ref run() must be called before using this function.  | 
   444     /// \pre \ref run() must be called before using this function.  | 
   443     const PotentialMap &getPotential() const { return potential;} | 
   445     const PotentialMap &getPotential() const { return potential;} | 
   444   | 
   446   | 
   445     ///This function checks, whether the given solution is optimal  | 
   447     ///This function checks, whether the given solution is optimal  | 
   446     ///Running after a \c run() should return with true  | 
   448     ///Running after a \c run() should return with true  | 
   447     ///In this "state of the art" this only check optimality, doesn't bother with feasibility  | 
   449     ///In this "state of the art" this only checks optimality, doesn't bother with feasibility  | 
   448     ///  | 
   450     ///  | 
   449     ///\todo Is this OK here?  | 
   451     ///\todo Is this OK here?  | 
   450     bool checkComplementarySlackness(){ | 
   452     bool checkComplementarySlackness(){ | 
   451       Cost mod_pot;  | 
   453       Cost mod_pot;  | 
   452       Cost fl_e;  | 
   454       Cost fl_e;  | 
   453       FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){ | 
   455       FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){ | 
   454 	//C^{\Pi}_{i,j} | 
   456 	//C^{\Pi}_{i,j} | 
   455 	mod_pot = cost[e]-potential[graph.head(e)]+potential[graph.tail(e)];  | 
   457 	mod_pot = cost[e]-potential[graph.head(e)]+potential[graph.tail(e)];  | 
   456 	fl_e = flow[e];  | 
   458 	fl_e = flow[e];  | 
   457 	//	std::cout << fl_e << std::endl;  | 
   459 	//	std::cout << fl_e << std::endl;  | 
   458 	if (0<fl_e && fl_e<capacity[e]){ | 
   460 	if (mod_pot > 0 && fl_e != 0)  | 
   459 	  if (mod_pot != 0)  | 
   461 	  return false;  | 
   460 	    return false;  | 
   462   | 
   461 	}  | 
         | 
   462 	else{ | 
         | 
   463 	  if (mod_pot > 0 && fl_e != 0)  | 
         | 
   464 	    return false;  | 
         | 
   465 	  if (mod_pot < 0 && fl_e != capacity[e])  | 
         | 
   466 	    return false;  | 
         | 
   467 	}  | 
         | 
   468       }  | 
   463       }  | 
   469       return true;  | 
   464       return true;  | 
   470     }  | 
   465     }  | 
   471       | 
   466   | 
         | 
   467     /*  | 
         | 
   468     //For testing purposes only  | 
         | 
   469     //Lists the node_properties  | 
         | 
   470     void write_property_vector(const SupplyDemandMap& a,  | 
         | 
   471 			       char* prop_name="property"){ | 
         | 
   472       FOR_EACH_LOC(typename Graph::NodeIt, i, graph){ | 
         | 
   473 	cout<<"Node id.: "<<graph.id(i)<<", "<<prop_name<<" value: "<<a[i]<<endl;  | 
         | 
   474       }  | 
         | 
   475       cout<<endl;  | 
         | 
   476     }  | 
         | 
   477     */  | 
         | 
   478     bool checkFeasibility(){ | 
         | 
   479       SupplyDemandMap supdem(graph);  | 
         | 
   480       FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){ | 
         | 
   481   | 
         | 
   482 	if ( flow[e] < 0){ | 
         | 
   483   | 
         | 
   484 	  return false;  | 
         | 
   485 	}  | 
         | 
   486 	supdem[graph.tail(e)] += flow[e];  | 
         | 
   487 	supdem[graph.head(e)] -= flow[e];  | 
         | 
   488       }  | 
         | 
   489       //write_property_vector(supdem, "supdem");  | 
         | 
   490       //write_property_vector(supply_demand, "supply_demand");  | 
         | 
   491   | 
         | 
   492       FOR_EACH_LOC(typename Graph::NodeIt, n, graph){ | 
         | 
   493   | 
         | 
   494 	if ( supdem[n] != supply_demand[n]){ | 
         | 
   495 	  //cout<<"Node id.: "<<graph.id(n)<<" : "<<supdem[n]<<", should be: "<<supply_demand[n]<<endl;  | 
         | 
   496 	  return false;  | 
         | 
   497 	}  | 
         | 
   498       }  | 
         | 
   499   | 
         | 
   500       return true;  | 
         | 
   501     }  | 
         | 
   502   | 
         | 
   503     bool checkOptimality(){ | 
         | 
   504       return checkFeasibility() && checkComplementarySlackness();  | 
         | 
   505     }  | 
   472   | 
   506   | 
   473   }; //class MinCostFlow  | 
   507   }; //class MinCostFlow  | 
   474   | 
   508   | 
   475   ///@}  | 
   509   ///@}  | 
   476   | 
   510   |