gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
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.
0 2 0
default
2 files changed with 184 insertions and 318 deletions:
↑ Collapse diff ↑
Ignore white space 6 line context
... ...
@@ -71,23 +71,12 @@
71 71
  class NetworkSimplex
72 72
  {
73 73
  public:
74 74

	
75
    /// The flow type of the algorithm
75
    /// The type of the flow amounts, capacity bounds and supply values
76 76
    typedef V Value;
77
    /// The cost type of the algorithm
77
    /// The type of the arc costs
78 78
    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
90 79

	
91 80
  public:
92 81

	
93 82
    /// \brief Problem type constants for the \c run() function.
... ...
@@ -205,17 +194,13 @@
205 194
  private:
206 195

	
207 196
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
208 197

	
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

	
213 198
    typedef std::vector<Arc> ArcVector;
214 199
    typedef std::vector<Node> NodeVector;
215 200
    typedef std::vector<int> IntVector;
216 201
    typedef std::vector<bool> BoolVector;
217
    typedef std::vector<Value> FlowVector;
202
    typedef std::vector<Value> ValueVector;
218 203
    typedef std::vector<Cost> CostVector;
219 204

	
220 205
    // State constants for arcs
221 206
    enum ArcStateEnum {
... ...
@@ -231,36 +216,25 @@
231 216
    int _node_num;
232 217
    int _arc_num;
233 218

	
234 219
    // 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;
242 221
    SupplyType _stype;
243
    
244 222
    Value _sum_supply;
245 223

	
246
    // Result maps
247
    FlowMap *_flow_map;
248
    PotentialMap *_potential_map;
249
    bool _local_flow;
250
    bool _local_potential;
251

	
252 224
    // Data structures for storing the digraph
253 225
    IntNodeMap _node_id;
254
    ArcVector _arc_ref;
226
    IntArcMap _arc_id;
255 227
    IntVector _source;
256 228
    IntVector _target;
257 229

	
258 230
    // Node and arc data
259
    FlowVector _cap;
231
    ValueVector _lower;
232
    ValueVector _upper;
233
    ValueVector _cap;
260 234
    CostVector _cost;
261
    FlowVector _supply;
262
    FlowVector _flow;
235
    ValueVector _supply;
236
    ValueVector _flow;
263 237
    CostVector _pi;
264 238

	
265 239
    // Data for storing the spanning tree structure
266 240
    IntVector _parent;
... ...
@@ -688,14 +662,9 @@
688 662
    /// The constructor of the class.
689 663
    ///
690 664
    /// \param graph The digraph the algorithm runs on.
691 665
    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),
698 667
      INF(std::numeric_limits<Value>::has_infinity ?
699 668
          std::numeric_limits<Value>::infinity() :
700 669
          std::numeric_limits<Value>::max())
701 670
    {
... ...
@@ -703,14 +672,60 @@
703 672
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
704 673
        "The flow type of NetworkSimplex must be signed");
705 674
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
706 675
        "The cost type of NetworkSimplex must be signed");
707
    }
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;
708 682

	
709
    /// Destructor.
710
    ~NetworkSimplex() {
711
      if (_local_flow) delete _flow_map;
712
      if (_local_potential) delete _potential_map;
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;
713 728
    }
714 729

	
715 730
    /// \name Parameters
716 731
    /// The parameters of the algorithm can be specified using these
... ...
@@ -730,12 +745,11 @@
730 745
    ///
731 746
    /// \return <tt>(*this)</tt>
732 747
    template <typename LowerMap>
733 748
    NetworkSimplex& lowerMap(const LowerMap& map) {
734
      delete _plower;
735
      _plower = new ValueArcMap(_graph);
749
      _have_lower = true;
736 750
      for (ArcIt a(_graph); a != INVALID; ++a) {
737
        (*_plower)[a] = map[a];
751
        _lower[_arc_id[a]] = map[a];
738 752
      }
739 753
      return *this;
740 754
    }
741 755

	
... ...
@@ -752,12 +766,10 @@
752 766
    ///
753 767
    /// \return <tt>(*this)</tt>
754 768
    template<typename UpperMap>
755 769
    NetworkSimplex& upperMap(const UpperMap& map) {
756
      delete _pupper;
757
      _pupper = new ValueArcMap(_graph);
758 770
      for (ArcIt a(_graph); a != INVALID; ++a) {
759
        (*_pupper)[a] = map[a];
771
        _upper[_arc_id[a]] = map[a];
760 772
      }
761 773
      return *this;
762 774
    }
763 775

	
... ...
@@ -773,12 +785,10 @@
773 785
    ///
774 786
    /// \return <tt>(*this)</tt>
775 787
    template<typename CostMap>
776 788
    NetworkSimplex& costMap(const CostMap& map) {
777
      delete _pcost;
778
      _pcost = new CostArcMap(_graph);
779 789
      for (ArcIt a(_graph); a != INVALID; ++a) {
780
        (*_pcost)[a] = map[a];
790
        _cost[_arc_id[a]] = map[a];
781 791
      }
782 792
      return *this;
783 793
    }
784 794

	
... ...
@@ -795,13 +805,10 @@
795 805
    ///
796 806
    /// \return <tt>(*this)</tt>
797 807
    template<typename SupplyMap>
798 808
    NetworkSimplex& supplyMap(const SupplyMap& map) {
799
      delete _psupply;
800
      _pstsup = false;
801
      _psupply = new ValueNodeMap(_graph);
802 809
      for (NodeIt n(_graph); n != INVALID; ++n) {
803
        (*_psupply)[n] = map[n];
810
        _supply[_node_id[n]] = map[n];
804 811
      }
805 812
      return *this;
806 813
    }
807 814

	
... ...
@@ -823,14 +830,13 @@
823 830
    /// (i.e. the supply of \c s and the demand of \c t).
824 831
    ///
825 832
    /// \return <tt>(*this)</tt>
826 833
    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;
833 839
      return *this;
834 840
    }
835 841
    
836 842
    /// \brief Set the type of the supply constraints.
... ...
@@ -846,43 +852,8 @@
846 852
      _stype = supply_type;
847 853
      return *this;
848 854
    }
849 855

	
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
    
885 856
    /// @}
886 857

	
887 858
    /// \name Execution Control
888 859
    /// The algorithm can be executed using \ref run().
... ...
@@ -893,9 +864,9 @@
893 864
    ///
894 865
    /// This function runs the algorithm.
895 866
    /// The paramters can be specified using functions \ref lowerMap(),
896 867
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 
897
    /// \ref supplyType(), \ref flowMap() and \ref potentialMap().
868
    /// \ref supplyType().
898 869
    /// For example,
899 870
    /// \code
900 871
    ///   NetworkSimplex<ListDigraph> ns(graph);
901 872
    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
... ...
@@ -905,8 +876,10 @@
905 876
    /// This function can be called more than once. All the parameters
906 877
    /// that have been given are kept for the next call, unless
907 878
    /// \ref reset() is called, thus only the modified parameters
908 879
    /// 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.
909 882
    ///
910 883
    /// \param pivot_rule The pivot rule that will be used during the
911 884
    /// algorithm. For more information see \ref PivotRule.
912 885
    ///
... ...
@@ -927,14 +900,15 @@
927 900
    /// \brief Reset all the parameters that have been given before.
928 901
    ///
929 902
    /// This function resets all the paramaters that have been given
930 903
    /// 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().
933 905
    ///
934 906
    /// It is useful for multiple run() calls. If this function is not
935 907
    /// used, all the parameters given before are kept for the next
936 908
    /// \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.
937 911
    ///
938 912
    /// For example,
939 913
    /// \code
940 914
    ///   NetworkSimplex<ListDigraph> ns(graph);
... ...
@@ -956,25 +930,18 @@
956 930
    /// \endcode
957 931
    ///
958 932
    /// \return <tt>(*this)</tt>
959 933
    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;
969 943
      _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

	
977 944
      return *this;
978 945
    }
979 946

	
980 947
    /// @}
... ...
@@ -1000,17 +967,14 @@
1000 967
    /// type of the algorithm, which is the default return type of the
1001 968
    /// function.
1002 969
    ///
1003 970
    /// \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]);
1013 977
      }
1014 978
      return c;
1015 979
    }
1016 980

	
... ...
@@ -1025,19 +989,23 @@
1025 989
    /// This function returns the flow on the given arc.
1026 990
    ///
1027 991
    /// \pre \ref run() must be called before using this function.
1028 992
    Value flow(const Arc& a) const {
1029
      return (*_flow_map)[a];
993
      return _flow[_arc_id[a]];
1030 994
    }
1031 995

	
1032
    /// \brief Return a const reference to the flow map.
996
    /// \brief Return the flow map (the primal solution).
1033 997
    ///
1034
    /// This function returns a const reference to an arc map storing
1035
    /// the found flow.
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.
1036 1001
    ///
1037 1002
    /// \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
      }
1040 1008
    }
1041 1009

	
1042 1010
    /// \brief Return the potential (dual value) of the given node.
1043 1011
    ///
... ...
@@ -1045,92 +1013,59 @@
1045 1013
    /// given node.
1046 1014
    ///
1047 1015
    /// \pre \ref run() must be called before using this function.
1048 1016
    Cost potential(const Node& n) const {
1049
      return (*_potential_map)[n];
1017
      return _pi[_node_id[n]];
1050 1018
    }
1051 1019

	
1052
    /// \brief Return a const reference to the potential map
1053
    /// (the dual solution).
1020
    /// \brief Return the potential map (the dual solution).
1054 1021
    ///
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".
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.
1058 1026
    ///
1059 1027
    /// \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
      }
1062 1033
    }
1063 1034

	
1064 1035
    /// @}
1065 1036

	
1066 1037
  private:
1067 1038

	
1068 1039
    // Initialize internal data structures
1069 1040
    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;
1085 1041
      if (_node_num == 0) return false;
1086 1042

	
1087
      _arc_ref.resize(_arc_num);
1088
      _source.resize(all_arc_num);
1089
      _target.resize(all_arc_num);
1043
      // Check the sum of supply values
1044
      _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;
1090 1050

	
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;
1108
      _sum_supply = 0;
1109
      if (!_pstsup && !_psupply) {
1110
        _pstsup = true;
1111
        _psource = _ptarget = NodeIt(_graph);
1112
        _pstflow = 0;
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
        }
1063
      } else {
1064
        for (int i = 0; i != _arc_num; ++i) {
1065
          _cap[i] = _upper[i];
1066
        }
1113 1067
      }
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);
1123
      } 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;
1133 1068

	
1134 1069
      // Initialize artifical cost
1135 1070
      Cost ART_COST;
1136 1071
      if (std::numeric_limits<Cost>::is_exact) {
... ...
@@ -1142,95 +1077,33 @@
1142 1077
        }
1143 1078
        ART_COST = (ART_COST + 1) * _node_num;
1144 1079
      }
1145 1080

	
1081
      // Initialize arc maps
1082
      for (int i = 0; i != _arc_num; ++i) {
1083
        _flow[i] = 0;
1084
        _state[i] = STATE_LOWER;
1085
      }
1086
      
1146 1087
      // Set data for the artificial root node
1147 1088
      _root = _node_num;
1148 1089
      _parent[_root] = -1;
1149 1090
      _pred[_root] = -1;
1150 1091
      _thread[_root] = 0;
1151 1092
      _rev_thread[0] = _root;
1152
      _succ_num[_root] = all_node_num;
1093
      _succ_num[_root] = _node_num + 1;
1153 1094
      _last_succ[_root] = _root - 1;
1154 1095
      _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;
1224 1097

	
1225 1098
      // Add artificial arcs and initialize the spanning tree data structure
1226 1099
      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1100
        _parent[u] = _root;
1101
        _pred[u] = e;
1227 1102
        _thread[u] = u + 1;
1228 1103
        _rev_thread[u + 1] = u;
1229 1104
        _succ_num[u] = 1;
1230 1105
        _last_succ[u] = u;
1231
        _parent[u] = _root;
1232
        _pred[u] = e;
1233 1106
        _cost[e] = ART_COST;
1234 1107
        _cap[e] = INF;
1235 1108
        _state[e] = STATE_TREE;
1236 1109
        if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
... ...
@@ -1516,22 +1389,18 @@
1516 1389
          if (_flow[e] != 0) return INFEASIBLE;
1517 1390
        }
1518 1391
      }
1519 1392

	
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) {
1522 1395
        for (int i = 0; i != _arc_num; ++i) {
1523
          Arc e = _arc_ref[i];
1524
          _flow_map->set(e, (*_plower)[e] + _flow[i]);
1396
          Value c = _lower[i];
1397
          if (c != 0) {
1398
            _flow[i] += c;
1399
            _supply[_source[i]] += c;
1400
            _supply[_target[i]] -= c;
1401
          }
1525 1402
        }
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]]);
1534 1403
      }
1535 1404

	
1536 1405
      return OPTIMAL;
1537 1406
    }
Ignore white space 8 line context
... ...
@@ -83,9 +83,9 @@
83 83
  LEQ
84 84
};
85 85

	
86 86
// Check the interface of an MCF algorithm
87
template <typename GR, typename Flow, typename Cost>
87
template <typename GR, typename Value, typename Cost>
88 88
class McfClassConcept
89 89
{
90 90
public:
91 91

	
... ...
@@ -94,56 +94,48 @@
94 94
    void constraints() {
95 95
      checkConcept<concepts::Digraph, GR>();
96 96

	
97 97
      MCF mcf(g);
98
      const MCF& const_mcf = mcf;
98 99

	
99 100
      b = mcf.reset()
100 101
             .lowerMap(lower)
101 102
             .upperMap(upper)
102 103
             .costMap(cost)
103 104
             .supplyMap(sup)
104 105
             .stSupply(n, n, k)
105
             .flowMap(flow)
106
             .potentialMap(pot)
107 106
             .run();
108
      
109
      const MCF& const_mcf = mcf;
110

	
111
      const typename MCF::FlowMap &fm = const_mcf.flowMap();
112
      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
113 107

	
114 108
      c = const_mcf.totalCost();
115
      double x = const_mcf.template totalCost<double>();
109
      x = const_mcf.template totalCost<double>();
116 110
      v = const_mcf.flow(a);
117 111
      c = const_mcf.potential(n);
118
      
119
      v = const_mcf.INF;
120

	
121
      ignore_unused_variable_warning(fm);
122
      ignore_unused_variable_warning(pm);
123
      ignore_unused_variable_warning(x);
112
      const_mcf.flowMap(fm);
113
      const_mcf.potentialMap(pm);
124 114
    }
125 115

	
126 116
    typedef typename GR::Node Node;
127 117
    typedef typename GR::Arc Arc;
128
    typedef concepts::ReadMap<Node, Flow> NM;
129
    typedef concepts::ReadMap<Arc, Flow> FAM;
118
    typedef concepts::ReadMap<Node, Value> NM;
119
    typedef concepts::ReadMap<Arc, Value> VAM;
130 120
    typedef concepts::ReadMap<Arc, Cost> CAM;
121
    typedef concepts::WriteMap<Arc, Value> FlowMap;
122
    typedef concepts::WriteMap<Node, Cost> PotMap;
131 123

	
132 124
    const GR &g;
133
    const FAM &lower;
134
    const FAM &upper;
125
    const VAM &lower;
126
    const VAM &upper;
135 127
    const CAM &cost;
136 128
    const NM &sup;
137 129
    const Node &n;
138 130
    const Arc &a;
139
    const Flow &k;
140
    Flow v;
141
    Cost c;
131
    const Value &k;
132
    FlowMap fm;
133
    PotMap pm;
142 134
    bool b;
143

	
144
    typename MCF::FlowMap &flow;
145
    typename MCF::PotentialMap &pot;
135
    double x;
136
    typename MCF::Value v;
137
    typename MCF::Cost c;
146 138
  };
147 139

	
148 140
};
149 141

	
... ...
@@ -220,26 +212,31 @@
220 212
               SupplyType type = EQ )
221 213
{
222 214
  check(mcf_result == result, "Wrong result " + test_id);
223 215
  if (optimal) {
224
    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
216
    typename GR::template ArcMap<typename SM::Value> flow(gr);
217
    typename GR::template NodeMap<typename CM::Value> pi(gr);
218
    mcf.flowMap(flow);
219
    mcf.potentialMap(pi);
220
    check(checkFlow(gr, lower, upper, supply, flow, type),
225 221
          "The flow is not feasible " + test_id);
226 222
    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
227
    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
228
                         mcf.potentialMap()),
223
    check(checkPotential(gr, lower, upper, cost, supply, flow, pi),
229 224
          "Wrong potentials " + test_id);
230 225
  }
231 226
}
232 227

	
233 228
int main()
234 229
{
235 230
  // Check the interfaces
236 231
  {
237
    typedef int Flow;
238
    typedef int Cost;
239 232
    typedef concepts::Digraph GR;
240
    checkConcept< McfClassConcept<GR, Flow, Cost>,
241
                  NetworkSimplex<GR, Flow, Cost> >();
233
    checkConcept< McfClassConcept<GR, int, int>,
234
                  NetworkSimplex<GR> >();
235
    checkConcept< McfClassConcept<GR, double, double>,
236
                  NetworkSimplex<GR, double> >();
237
    checkConcept< McfClassConcept<GR, int, double>,
238
                  NetworkSimplex<GR, int, double> >();
242 239
  }
243 240

	
244 241
  // Run various MCF tests
245 242
  typedef ListDigraph Digraph;
0 comments (0 inline)