src/work/athos/mincostflow.h
changeset 672 6c7bd0edd1d7
parent 671 708df4dc6ab6
child 921 818510fa3d99
equal deleted inserted replaced
7:e4fcfda0d161 8:8428da634758
    97     //To store the flow
    97     //To store the flow
    98     FlowMap flow; 
    98     FlowMap flow; 
    99     //To store the potential (dual variables)
    99     //To store the potential (dual variables)
   100     typedef typename Graph::template NodeMap<Cost> PotentialMap;
   100     typedef typename Graph::template NodeMap<Cost> PotentialMap;
   101     PotentialMap potential;
   101     PotentialMap potential;
   102     //To store excess-deficit values
       
   103     SupplyDemandMap excess_deficit;
       
   104     
   102     
   105 
   103 
   106     Cost total_cost;
   104     Cost total_cost;
   107 
   105 
   108 
   106 
   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 );
   318 	  }
   320 	  }
   319 	}
   321 	}
   320 
   322 
   321 
   323 
   322 	Node s = excess_nodes.top(); 
   324 	Node s = excess_nodes.top(); 
   323 	SupplyDemand max_excess = excess_nodes[s];
   325 	max_excess = excess_nodes[s];
   324 	Node t = deficit_nodes.top(); 
   326 	Node t = deficit_nodes.top(); 
   325 	if (max_excess < deficit_nodes[t]){
   327 	if (max_excess < deficit_nodes[t]){
   326 	  max_excess = deficit_nodes[t];
   328 	  max_excess = deficit_nodes[t];
   327 	}
   329 	}
   328 
   330 
   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