COIN-OR::LEMON - Graph Library

Changeset 689:111698359429 in lemon for lemon


Ignore:
Timestamp:
04/29/09 16:54:27 (15 years ago)
Author:
Peter Kovacs <kpeter@…>
Branch:
default
Phase:
public
Message:

Less map copying in NetworkSimplex? (#234)

  • The graph is copied in the constructor instead of the init() function. It must not be modified after the class is constructed.
  • The maps are copied once (instead of twice).
  • Remove FlowMap?, PotentialMap? typedefs and flowMap(), pontentialMap() setter functions.
  • flowMap() and potentialMap() query functions copy the values into the given map (reference) instead of returning a const reference to a previously constructed map.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/network_simplex.h

    r688 r689  
    7373  public:
    7474
    75     /// The flow type of the algorithm
     75    /// The type of the flow amounts, capacity bounds and supply values
    7676    typedef V Value;
    77     /// The cost type of the algorithm
     77    /// The type of the arc costs
    7878    typedef C Cost;
    79 #ifdef DOXYGEN
    80     /// The type of the flow map
    81     typedef GR::ArcMap<Value> FlowMap;
    82     /// The type of the potential map
    83     typedef GR::NodeMap<Cost> PotentialMap;
    84 #else
    85     /// The type of the flow map
    86     typedef typename GR::template ArcMap<Value> FlowMap;
    87     /// The type of the potential map
    88     typedef typename GR::template NodeMap<Cost> PotentialMap;
    89 #endif
    9079
    9180  public:
     
    207196    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
    208197
    209     typedef typename GR::template ArcMap<Value> ValueArcMap;
    210     typedef typename GR::template ArcMap<Cost> CostArcMap;
    211     typedef typename GR::template NodeMap<Value> ValueNodeMap;
    212 
    213198    typedef std::vector<Arc> ArcVector;
    214199    typedef std::vector<Node> NodeVector;
    215200    typedef std::vector<int> IntVector;
    216201    typedef std::vector<bool> BoolVector;
    217     typedef std::vector<Value> FlowVector;
     202    typedef std::vector<Value> ValueVector;
    218203    typedef std::vector<Cost> CostVector;
    219204
     
    233218
    234219    // Parameters of the problem
    235     ValueArcMap *_plower;
    236     ValueArcMap *_pupper;
    237     CostArcMap *_pcost;
    238     ValueNodeMap *_psupply;
    239     bool _pstsup;
    240     Node _psource, _ptarget;
    241     Value _pstflow;
     220    bool _have_lower;
    242221    SupplyType _stype;
    243    
    244222    Value _sum_supply;
    245 
    246     // Result maps
    247     FlowMap *_flow_map;
    248     PotentialMap *_potential_map;
    249     bool _local_flow;
    250     bool _local_potential;
    251223
    252224    // Data structures for storing the digraph
    253225    IntNodeMap _node_id;
    254     ArcVector _arc_ref;
     226    IntArcMap _arc_id;
    255227    IntVector _source;
    256228    IntVector _target;
    257229
    258230    // Node and arc data
    259     FlowVector _cap;
     231    ValueVector _lower;
     232    ValueVector _upper;
     233    ValueVector _cap;
    260234    CostVector _cost;
    261     FlowVector _supply;
    262     FlowVector _flow;
     235    ValueVector _supply;
     236    ValueVector _flow;
    263237    CostVector _pi;
    264238
     
    690664    /// \param graph The digraph the algorithm runs on.
    691665    NetworkSimplex(const GR& graph) :
    692       _graph(graph),
    693       _plower(NULL), _pupper(NULL), _pcost(NULL),
    694       _psupply(NULL), _pstsup(false), _stype(GEQ),
    695       _flow_map(NULL), _potential_map(NULL),
    696       _local_flow(false), _local_potential(false),
    697       _node_id(graph),
     666      _graph(graph), _node_id(graph), _arc_id(graph),
    698667      INF(std::numeric_limits<Value>::has_infinity ?
    699668          std::numeric_limits<Value>::infinity() :
     
    705674      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
    706675        "The cost type of NetworkSimplex must be signed");
    707     }
    708 
    709     /// Destructor.
    710     ~NetworkSimplex() {
    711       if (_local_flow) delete _flow_map;
    712       if (_local_potential) delete _potential_map;
     676       
     677      // Resize vectors
     678      _node_num = countNodes(_graph);
     679      _arc_num = countArcs(_graph);
     680      int all_node_num = _node_num + 1;
     681      int all_arc_num = _arc_num + _node_num;
     682
     683      _source.resize(all_arc_num);
     684      _target.resize(all_arc_num);
     685
     686      _lower.resize(all_arc_num);
     687      _upper.resize(all_arc_num);
     688      _cap.resize(all_arc_num);
     689      _cost.resize(all_arc_num);
     690      _supply.resize(all_node_num);
     691      _flow.resize(all_arc_num);
     692      _pi.resize(all_node_num);
     693
     694      _parent.resize(all_node_num);
     695      _pred.resize(all_node_num);
     696      _forward.resize(all_node_num);
     697      _thread.resize(all_node_num);
     698      _rev_thread.resize(all_node_num);
     699      _succ_num.resize(all_node_num);
     700      _last_succ.resize(all_node_num);
     701      _state.resize(all_arc_num);
     702
     703      // Copy the graph (store the arcs in a mixed order)
     704      int i = 0;
     705      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
     706        _node_id[n] = i;
     707      }
     708      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
     709      i = 0;
     710      for (ArcIt a(_graph); a != INVALID; ++a) {
     711        _arc_id[a] = i;
     712        _source[i] = _node_id[_graph.source(a)];
     713        _target[i] = _node_id[_graph.target(a)];
     714        if ((i += k) >= _arc_num) i = (i % k) + 1;
     715      }
     716     
     717      // Initialize maps
     718      for (int i = 0; i != _node_num; ++i) {
     719        _supply[i] = 0;
     720      }
     721      for (int i = 0; i != _arc_num; ++i) {
     722        _lower[i] = 0;
     723        _upper[i] = INF;
     724        _cost[i] = 1;
     725      }
     726      _have_lower = false;
     727      _stype = GEQ;
    713728    }
    714729
     
    732747    template <typename LowerMap>
    733748    NetworkSimplex& lowerMap(const LowerMap& map) {
    734       delete _plower;
    735       _plower = new ValueArcMap(_graph);
     749      _have_lower = true;
    736750      for (ArcIt a(_graph); a != INVALID; ++a) {
    737         (*_plower)[a] = map[a];
     751        _lower[_arc_id[a]] = map[a];
    738752      }
    739753      return *this;
     
    754768    template<typename UpperMap>
    755769    NetworkSimplex& upperMap(const UpperMap& map) {
    756       delete _pupper;
    757       _pupper = new ValueArcMap(_graph);
    758770      for (ArcIt a(_graph); a != INVALID; ++a) {
    759         (*_pupper)[a] = map[a];
     771        _upper[_arc_id[a]] = map[a];
    760772      }
    761773      return *this;
     
    775787    template<typename CostMap>
    776788    NetworkSimplex& costMap(const CostMap& map) {
    777       delete _pcost;
    778       _pcost = new CostArcMap(_graph);
    779789      for (ArcIt a(_graph); a != INVALID; ++a) {
    780         (*_pcost)[a] = map[a];
     790        _cost[_arc_id[a]] = map[a];
    781791      }
    782792      return *this;
     
    797807    template<typename SupplyMap>
    798808    NetworkSimplex& supplyMap(const SupplyMap& map) {
    799       delete _psupply;
    800       _pstsup = false;
    801       _psupply = new ValueNodeMap(_graph);
    802809      for (NodeIt n(_graph); n != INVALID; ++n) {
    803         (*_psupply)[n] = map[n];
     810        _supply[_node_id[n]] = map[n];
    804811      }
    805812      return *this;
     
    825832    /// \return <tt>(*this)</tt>
    826833    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
    827       delete _psupply;
    828       _psupply = NULL;
    829       _pstsup = true;
    830       _psource = s;
    831       _ptarget = t;
    832       _pstflow = k;
     834      for (int i = 0; i != _node_num; ++i) {
     835        _supply[i] = 0;
     836      }
     837      _supply[_node_id[s]] =  k;
     838      _supply[_node_id[t]] = -k;
    833839      return *this;
    834840    }
     
    848854    }
    849855
    850     /// \brief Set the flow map.
    851     ///
    852     /// This function sets the flow map.
    853     /// If it is not used before calling \ref run(), an instance will
    854     /// be allocated automatically. The destructor deallocates this
    855     /// automatically allocated map, of course.
    856     ///
    857     /// \return <tt>(*this)</tt>
    858     NetworkSimplex& flowMap(FlowMap& map) {
    859       if (_local_flow) {
    860         delete _flow_map;
    861         _local_flow = false;
    862       }
    863       _flow_map = &map;
    864       return *this;
    865     }
    866 
    867     /// \brief Set the potential map.
    868     ///
    869     /// This function sets the potential map, which is used for storing
    870     /// the dual solution.
    871     /// If it is not used before calling \ref run(), an instance will
    872     /// be allocated automatically. The destructor deallocates this
    873     /// automatically allocated map, of course.
    874     ///
    875     /// \return <tt>(*this)</tt>
    876     NetworkSimplex& potentialMap(PotentialMap& map) {
    877       if (_local_potential) {
    878         delete _potential_map;
    879         _local_potential = false;
    880       }
    881       _potential_map = &map;
    882       return *this;
    883     }
    884    
    885856    /// @}
    886857
     
    895866    /// The paramters can be specified using functions \ref lowerMap(),
    896867    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
    897     /// \ref supplyType(), \ref flowMap() and \ref potentialMap().
     868    /// \ref supplyType().
    898869    /// For example,
    899870    /// \code
     
    907878    /// \ref reset() is called, thus only the modified parameters
    908879    /// have to be set again. See \ref reset() for examples.
     880    /// However the underlying digraph must not be modified after this
     881    /// class have been constructed, since it copies and extends the graph.
    909882    ///
    910883    /// \param pivot_rule The pivot rule that will be used during the
     
    929902    /// This function resets all the paramaters that have been given
    930903    /// before using functions \ref lowerMap(), \ref upperMap(),
    931     /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(),
    932     /// \ref flowMap() and \ref potentialMap().
     904    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
    933905    ///
    934906    /// It is useful for multiple run() calls. If this function is not
    935907    /// used, all the parameters given before are kept for the next
    936908    /// \ref run() call.
     909    /// However the underlying digraph must not be modified after this
     910    /// class have been constructed, since it copies and extends the graph.
    937911    ///
    938912    /// For example,
     
    958932    /// \return <tt>(*this)</tt>
    959933    NetworkSimplex& reset() {
    960       delete _plower;
    961       delete _pupper;
    962       delete _pcost;
    963       delete _psupply;
    964       _plower = NULL;
    965       _pupper = NULL;
    966       _pcost = NULL;
    967       _psupply = NULL;
    968       _pstsup = false;
     934      for (int i = 0; i != _node_num; ++i) {
     935        _supply[i] = 0;
     936      }
     937      for (int i = 0; i != _arc_num; ++i) {
     938        _lower[i] = 0;
     939        _upper[i] = INF;
     940        _cost[i] = 1;
     941      }
     942      _have_lower = false;
    969943      _stype = GEQ;
    970       if (_local_flow) delete _flow_map;
    971       if (_local_potential) delete _potential_map;
    972       _flow_map = NULL;
    973       _potential_map = NULL;
    974       _local_flow = false;
    975       _local_potential = false;
    976 
    977944      return *this;
    978945    }
     
    1002969    ///
    1003970    /// \pre \ref run() must be called before using this function.
    1004     template <typename Value>
    1005     Value totalCost() const {
    1006       Value c = 0;
    1007       if (_pcost) {
    1008         for (ArcIt e(_graph); e != INVALID; ++e)
    1009           c += (*_flow_map)[e] * (*_pcost)[e];
    1010       } else {
    1011         for (ArcIt e(_graph); e != INVALID; ++e)
    1012           c += (*_flow_map)[e];
     971    template <typename Number>
     972    Number totalCost() const {
     973      Number c = 0;
     974      for (ArcIt a(_graph); a != INVALID; ++a) {
     975        int i = _arc_id[a];
     976        c += Number(_flow[i]) * Number(_cost[i]);
    1013977      }
    1014978      return c;
     
    1027991    /// \pre \ref run() must be called before using this function.
    1028992    Value flow(const Arc& a) const {
    1029       return (*_flow_map)[a];
    1030     }
    1031 
    1032     /// \brief Return a const reference to the flow map.
    1033     ///
    1034     /// This function returns a const reference to an arc map storing
    1035     /// the found flow.
     993      return _flow[_arc_id[a]];
     994    }
     995
     996    /// \brief Return the flow map (the primal solution).
     997    ///
     998    /// This function copies the flow value on each arc into the given
     999    /// map. The \c Value type of the algorithm must be convertible to
     1000    /// the \c Value type of the map.
    10361001    ///
    10371002    /// \pre \ref run() must be called before using this function.
    1038     const FlowMap& flowMap() const {
    1039       return *_flow_map;
     1003    template <typename FlowMap>
     1004    void flowMap(FlowMap &map) const {
     1005      for (ArcIt a(_graph); a != INVALID; ++a) {
     1006        map.set(a, _flow[_arc_id[a]]);
     1007      }
    10401008    }
    10411009
     
    10471015    /// \pre \ref run() must be called before using this function.
    10481016    Cost potential(const Node& n) const {
    1049       return (*_potential_map)[n];
    1050     }
    1051 
    1052     /// \brief Return a const reference to the potential map
    1053     /// (the dual solution).
    1054     ///
    1055     /// This function returns a const reference to a node map storing
    1056     /// the found potentials, which form the dual solution of the
    1057     /// \ref min_cost_flow "minimum cost flow problem".
     1017      return _pi[_node_id[n]];
     1018    }
     1019
     1020    /// \brief Return the potential map (the dual solution).
     1021    ///
     1022    /// This function copies the potential (dual value) of each node
     1023    /// into the given map.
     1024    /// The \c Cost type of the algorithm must be convertible to the
     1025    /// \c Value type of the map.
    10581026    ///
    10591027    /// \pre \ref run() must be called before using this function.
    1060     const PotentialMap& potentialMap() const {
    1061       return *_potential_map;
     1028    template <typename PotentialMap>
     1029    void potentialMap(PotentialMap &map) const {
     1030      for (NodeIt n(_graph); n != INVALID; ++n) {
     1031        map.set(n, _pi[_node_id[n]]);
     1032      }
    10621033    }
    10631034
     
    10681039    // Initialize internal data structures
    10691040    bool init() {
    1070       // Initialize result maps
    1071       if (!_flow_map) {
    1072         _flow_map = new FlowMap(_graph);
    1073         _local_flow = true;
    1074       }
    1075       if (!_potential_map) {
    1076         _potential_map = new PotentialMap(_graph);
    1077         _local_potential = true;
    1078       }
    1079 
    1080       // Initialize vectors
    1081       _node_num = countNodes(_graph);
    1082       _arc_num = countArcs(_graph);
    1083       int all_node_num = _node_num + 1;
    1084       int all_arc_num = _arc_num + _node_num;
    10851041      if (_node_num == 0) return false;
    10861042
    1087       _arc_ref.resize(_arc_num);
    1088       _source.resize(all_arc_num);
    1089       _target.resize(all_arc_num);
    1090 
    1091       _cap.resize(all_arc_num);
    1092       _cost.resize(all_arc_num);
    1093       _supply.resize(all_node_num);
    1094       _flow.resize(all_arc_num);
    1095       _pi.resize(all_node_num);
    1096 
    1097       _parent.resize(all_node_num);
    1098       _pred.resize(all_node_num);
    1099       _forward.resize(all_node_num);
    1100       _thread.resize(all_node_num);
    1101       _rev_thread.resize(all_node_num);
    1102       _succ_num.resize(all_node_num);
    1103       _last_succ.resize(all_node_num);
    1104       _state.resize(all_arc_num);
    1105 
    1106       // Initialize node related data
    1107       bool valid_supply = true;
     1043      // Check the sum of supply values
    11081044      _sum_supply = 0;
    1109       if (!_pstsup && !_psupply) {
    1110         _pstsup = true;
    1111         _psource = _ptarget = NodeIt(_graph);
    1112         _pstflow = 0;
    1113       }
    1114       if (_psupply) {
    1115         int i = 0;
    1116         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
    1117           _node_id[n] = i;
    1118           _supply[i] = (*_psupply)[n];
    1119           _sum_supply += _supply[i];
    1120         }
    1121         valid_supply = (_stype == GEQ && _sum_supply <= 0) ||
    1122                        (_stype == LEQ && _sum_supply >= 0);
     1045      for (int i = 0; i != _node_num; ++i) {
     1046        _sum_supply += _supply[i];
     1047      }
     1048      if ( !(_stype == GEQ && _sum_supply <= 0 ||
     1049             _stype == LEQ && _sum_supply >= 0) ) return false;
     1050
     1051      // Remove non-zero lower bounds
     1052      if (_have_lower) {
     1053        for (int i = 0; i != _arc_num; ++i) {
     1054          Value c = _lower[i];
     1055          if (c >= 0) {
     1056            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
     1057          } else {
     1058            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
     1059          }
     1060          _supply[_source[i]] -= c;
     1061          _supply[_target[i]] += c;
     1062        }
    11231063      } else {
    1124         int i = 0;
    1125         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
    1126           _node_id[n] = i;
    1127           _supply[i] = 0;
    1128         }
    1129         _supply[_node_id[_psource]] =  _pstflow;
    1130         _supply[_node_id[_ptarget]] = -_pstflow;
    1131       }
    1132       if (!valid_supply) return false;
     1064        for (int i = 0; i != _arc_num; ++i) {
     1065          _cap[i] = _upper[i];
     1066        }
     1067      }
    11331068
    11341069      // Initialize artifical cost
     
    11441079      }
    11451080
     1081      // Initialize arc maps
     1082      for (int i = 0; i != _arc_num; ++i) {
     1083        _flow[i] = 0;
     1084        _state[i] = STATE_LOWER;
     1085      }
     1086     
    11461087      // Set data for the artificial root node
    11471088      _root = _node_num;
     
    11501091      _thread[_root] = 0;
    11511092      _rev_thread[0] = _root;
    1152       _succ_num[_root] = all_node_num;
     1093      _succ_num[_root] = _node_num + 1;
    11531094      _last_succ[_root] = _root - 1;
    11541095      _supply[_root] = -_sum_supply;
    1155       if (_sum_supply < 0) {
    1156         _pi[_root] = -ART_COST;
    1157       } else {
    1158         _pi[_root] = ART_COST;
    1159       }
    1160 
    1161       // Store the arcs in a mixed order
    1162       int k = std::max(int(std::sqrt(double(_arc_num))), 10);
    1163       int i = 0;
    1164       for (ArcIt e(_graph); e != INVALID; ++e) {
    1165         _arc_ref[i] = e;
    1166         if ((i += k) >= _arc_num) i = (i % k) + 1;
    1167       }
    1168 
    1169       // Initialize arc maps
    1170       if (_pupper && _pcost) {
    1171         for (int i = 0; i != _arc_num; ++i) {
    1172           Arc e = _arc_ref[i];
    1173           _source[i] = _node_id[_graph.source(e)];
    1174           _target[i] = _node_id[_graph.target(e)];
    1175           _cap[i] = (*_pupper)[e];
    1176           _cost[i] = (*_pcost)[e];
    1177           _flow[i] = 0;
    1178           _state[i] = STATE_LOWER;
    1179         }
    1180       } else {
    1181         for (int i = 0; i != _arc_num; ++i) {
    1182           Arc e = _arc_ref[i];
    1183           _source[i] = _node_id[_graph.source(e)];
    1184           _target[i] = _node_id[_graph.target(e)];
    1185           _flow[i] = 0;
    1186           _state[i] = STATE_LOWER;
    1187         }
    1188         if (_pupper) {
    1189           for (int i = 0; i != _arc_num; ++i)
    1190             _cap[i] = (*_pupper)[_arc_ref[i]];
    1191         } else {
    1192           for (int i = 0; i != _arc_num; ++i)
    1193             _cap[i] = INF;
    1194         }
    1195         if (_pcost) {
    1196           for (int i = 0; i != _arc_num; ++i)
    1197             _cost[i] = (*_pcost)[_arc_ref[i]];
    1198         } else {
    1199           for (int i = 0; i != _arc_num; ++i)
    1200             _cost[i] = 1;
    1201         }
    1202       }
    1203      
    1204       // Remove non-zero lower bounds
    1205       if (_plower) {
    1206         for (int i = 0; i != _arc_num; ++i) {
    1207           Value c = (*_plower)[_arc_ref[i]];
    1208           if (c > 0) {
    1209             if (_cap[i] < INF) _cap[i] -= c;
    1210             _supply[_source[i]] -= c;
    1211             _supply[_target[i]] += c;
    1212           }
    1213           else if (c < 0) {
    1214             if (_cap[i] < INF + c) {
    1215               _cap[i] -= c;
    1216             } else {
    1217               _cap[i] = INF;
    1218             }
    1219             _supply[_source[i]] -= c;
    1220             _supply[_target[i]] += c;
    1221           }
    1222         }
    1223       }
     1096      _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
    12241097
    12251098      // Add artificial arcs and initialize the spanning tree data structure
    12261099      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
     1100        _parent[u] = _root;
     1101        _pred[u] = e;
    12271102        _thread[u] = u + 1;
    12281103        _rev_thread[u + 1] = u;
    12291104        _succ_num[u] = 1;
    12301105        _last_succ[u] = u;
    1231         _parent[u] = _root;
    1232         _pred[u] = e;
    12331106        _cost[e] = ART_COST;
    12341107        _cap[e] = INF;
     
    15181391      }
    15191392
    1520       // Copy flow values to _flow_map
    1521       if (_plower) {
     1393      // Transform the solution and the supply map to the original form
     1394      if (_have_lower) {
    15221395        for (int i = 0; i != _arc_num; ++i) {
    1523           Arc e = _arc_ref[i];
    1524           _flow_map->set(e, (*_plower)[e] + _flow[i]);
    1525         }
    1526       } else {
    1527         for (int i = 0; i != _arc_num; ++i) {
    1528           _flow_map->set(_arc_ref[i], _flow[i]);
    1529         }
    1530       }
    1531       // Copy potential values to _potential_map
    1532       for (NodeIt n(_graph); n != INVALID; ++n) {
    1533         _potential_map->set(n, _pi[_node_id[n]]);
     1396          Value c = _lower[i];
     1397          if (c != 0) {
     1398            _flow[i] += c;
     1399            _supply[_source[i]] += c;
     1400            _supply[_target[i]] -= c;
     1401          }
     1402        }
    15341403      }
    15351404
Note: See TracChangeset for help on using the changeset viewer.