Implement three pivot rules

src/smt/network_flow.h

@@ -23,8 +23,6 @@ Notes:
 
     It remains unclear how to convert DL assignment to a basic feasible solution of Network Simplex.
     A naive approach is to run an algorithm on max flow in order to get a spanning tree.
-
-    The network_simplex class hasn't had multiple pivoting strategies yet.
    
 --*/
 #ifndef _NETWORK_FLOW_H_
@@ -36,20 +34,207 @@ Notes:
 
 namespace smt {
 
+    enum pivot_rule {
+        // First eligible edge pivot rule
+        // Edges are traversed in a wraparound fashion
+        FIRST_ELIGIBLE,
+        // Best eligible edge pivot rule
+        // The best edge is selected in every iteration
+        BEST_ELIGIBLE,
+        // Candidate list pivot rule
+        // Major iterations: candidate list is built from eligible edges (in a wraparound way)
+        // Minor iterations: the best edge is selected from the list
+        CANDIDATE_LIST
+    };
+
     // Solve minimum cost flow problem using Network Simplex algorithm
     template<typename Ext>
     class network_flow : private Ext {
+    private:
         enum edge_state {
             LOWER = 1,
             BASIS = 0,
             UPPER = -1
         };
+
         typedef dl_var node;
         typedef dl_edge<Ext> edge;
         typedef dl_graph<Ext> graph;     
         typedef typename Ext::numeral numeral;
         typedef typename Ext::fin_numeral fin_numeral;
 
+        class pivot_rule_impl {
+        protected:
+            graph & m_graph;
+            svector<edge_state> & m_states;
+            vector<numeral> & m_potentials;
+            edge_id & m_enter_id;
+
+        public: 
+            pivot_rule_impl() {}
+            pivot_rule_impl(graph & g, vector<numeral> & potentials, 
+                            svector<edge_state> & states, edge_id & enter_id) 
+                : m_graph(g),
+                  m_potentials(potentials),
+                  m_states(states),
+                  m_enter_id(enter_id) {
+            }
+            bool choose_entering_edge() {return false;};
+        };
+        
+        class first_eligible_pivot : pivot_rule_impl {
+        private:
+            edge_id m_next_edge;
+
+        public:
+            first_eligible_pivot(graph & g, vector<numeral> & potentials, 
+                                 svector<edge_state> & states, edge_id & enter_id) : 
+                pivot_rule_impl(g, potentials, states, enter_id),
+                m_next_edge(0) {
+            }
+
+            bool choose_entering_edge() {
+                TRACE("network_flow", tout << "choose_entering_edge...\n";);        
+                unsigned num_edges = m_graph.get_num_edges();
+                for (unsigned i = m_next_edge; i < m_next_edge + num_edges; ++i) {
+                    edge_id id = (i >= num_edges) ? (i - num_edges) : i;
+                    node src = m_graph.get_source(id);
+                    node tgt = m_graph.get_target(id);
+                    if (m_states[id] != BASIS) {
+                        numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
+                        if (cost.is_pos()) {
+                            m_enter_id = id;
+                            TRACE("network_flow", {
+                                tout << "Found entering edge " << id << " between node ";
+                                tout << src << " and node " << tgt << " with reduced cost = " << cost << "...\n";
+                            });
+                            return true;
+                        }
+                    }
+                }
+                TRACE("network_flow", tout << "Found no entering edge...\n";);
+                return false;
+            };
+        };
+
+        class best_eligible_pivot : pivot_rule_impl {
+        public:
+            best_eligible_pivot(graph & g, vector<numeral> & potentials, 
+                                 svector<edge_state> & states, edge_id & enter_id) : 
+                pivot_rule_impl(g, potentials, states, enter_id) {
+            }
+
+            bool choose_entering_edge() {
+                TRACE("network_flow", tout << "choose_entering_edge...\n";);        
+                unsigned num_edges = m_graph.get_num_edges();
+                numeral max = numeral::zero();
+                for (unsigned i = 0; i < num_edges; ++i) {
+                    node src = m_graph.get_source(i);
+                    node tgt = m_graph.get_target(i);
+                    if (m_states[i] != BASIS) {
+                        numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(i);
+                        if (cost > max) {
+                            max = cost;
+                            m_enter_id = i;
+                        }
+                    }
+                }
+                if (max.is_pos()) {
+                    TRACE("network_flow", {
+                        tout << "Found entering edge " << m_enter_id << " between node ";
+                        tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
+                        tout << " with reduced cost = " << max << "...\n";
+                    });
+                    return true;
+                }
+                TRACE("network_flow", tout << "Found no entering edge...\n";);
+                return false;
+            };
+        };
+
+        class candidate_list_pivot : pivot_rule_impl {
+        private:
+            edge_id m_next_edge;
+            svector<edge_id> m_candidates;
+            unsigned num_candidates;
+            static const unsigned NUM_CANDIDATES = 10;
+
+        public:
+            candidate_list_pivot(graph & g, vector<numeral> & potentials, 
+                                 svector<edge_state> & states, edge_id & enter_id) : 
+                pivot_rule_impl(g, potentials, states, enter_id),
+                m_next_edge(0),
+                num_candidates(NUM_CANDIDATES),
+                m_candidates(num_candidates) {
+            }
+
+            bool choose_entering_edge() {
+                if (m_candidates.empty()) {
+                    // Build the candidate list
+                    unsigned num_edges = m_graph.get_num_edges();
+                    numeral max = numeral::zero();
+                    unsigned count = 0;
+                    for (unsigned i = m_next_edge; i < m_next_edge + num_edges; ++i) {
+                        edge_id id = (i >= num_edges) ? i - num_edges : i;
+                        node src = m_graph.get_source(id);
+                        node tgt = m_graph.get_target(id);
+                        if (m_states[id] != BASIS) {
+                            numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
+                            if (cost.is_pos()) {
+                                m_candidates[count++] = id;
+                                if (cost > max) {
+                                    max = cost;
+                                    m_enter_id = id;
+                                }
+                            }
+                            if (count >= num_candidates) break;
+                        }
+                    }
+                    m_next_edge = m_enter_id;
+                    if (max.is_pos()) {
+                        TRACE("network_flow", {
+                            tout << "Found entering edge " << m_enter_id << " between node ";
+                            tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
+                            tout << " with reduced cost = " << max << "...\n";
+                        });
+                        return true;
+                    }
+                    TRACE("network_flow", tout << "Found no entering edge...\n";);
+                    return false;
+                }
+                else {
+                    numeral max = numeral::zero();
+                    unsigned last = m_candidates.size();
+                    for (unsigned i = 0; i < last; ++i) {
+                        edge_id id = m_candidates[i];
+                        node src = m_graph.get_source(id);
+                        node tgt = m_graph.get_target(id);
+                        if (m_states[id] != BASIS) {
+                            numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
+                            if (cost > max) {
+                                max = cost;
+                                m_enter_id = id;
+                            }
+                            // Remove stale candidates
+                            if (!cost.is_pos()) {
+                                m_candidates[i] = m_candidates[--last];
+                            }
+                        }
+                    }
+                    if (max.is_pos()) {
+                        TRACE("network_flow", {
+                            tout << "Found entering edge " << m_enter_id << " between node ";
+                            tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
+                            tout << " with reduced cost = " << max << "...\n";
+                        });
+                        return true;
+                    }
+                    TRACE("network_flow", tout << "Found no entering edge...\n";);
+                    return false;
+                }
+            };
+        };
+        
         graph m_graph;
         thread_spanning_tree<Ext> m_tree;
 
@@ -76,9 +261,7 @@ namespace smt {
 
         void update_flows();
             
-        // If all reduced costs are non-negative, return false since the current spanning tree is optimal
+        bool choose_entering_edge(pivot_rule pr);
-        // Otherwise return true and update m_entering_edge
-        bool choose_entering_edge();
 
         // Send as much flow as possible around the cycle, the first basic edge with flow 0 will leave
         // Return false if the problem is unbounded
@@ -99,7 +282,7 @@ namespace smt {
         
         // Minimize cost flows
         // Return true if found an optimal solution, and return false if unbounded
-        bool min_cost();
+        bool min_cost(pivot_rule pr = FIRST_ELIGIBLE);
 
         // Compute the optimal solution
         numeral get_optimal_solution(vector<numeral> & result, bool is_dual);

src/smt/network_flow_def.h

@@ -129,30 +129,6 @@ namespace smt {
         TRACE("network_flow", tout << pp_vector("Flows", m_flows, true););
     }
     
-    template<typename Ext>
-    bool network_flow<Ext>::choose_entering_edge() {
-        TRACE("network_flow", tout << "choose_entering_edge...\n";);        
-        unsigned num_edges = m_graph.get_num_edges();
-        for (unsigned i = 0; i < num_edges; ++i) {
-            node src = m_graph.get_source(i);
-            node tgt = m_graph.get_target(i);
-            if (m_states[i] != BASIS) {
-                numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(i);
-                // TODO: add multiple pivoting strategies
-                if (cost.is_pos()) {
-                    m_enter_id = i;
-                    TRACE("network_flow", {
-                        tout << "Found entering edge " << i << " between node ";
-                        tout << src << " and node " << tgt << " with reduced cost = " << cost << "...\n";
-                    });
-                    return true;
-                }
-            }
-        }
-        TRACE("network_flow", tout << "Found no entering edge...\n";);
-        return false;
-    }
-
     template<typename Ext>
     bool network_flow<Ext>::choose_leaving_edge() {
         TRACE("network_flow", tout << "choose_leaving_edge...\n";);
@@ -191,12 +167,29 @@ namespace smt {
         m_tree.update(m_enter_id, m_leave_id);
     }
 
+    // FIXME: should declare pivot as a pivot_rule_impl and refactor
+    template<typename Ext>
+    bool network_flow<Ext>::choose_entering_edge(pivot_rule pr) {
+        if (pr == FIRST_ELIGIBLE) {
+            first_eligible_pivot pivot(m_graph, m_potentials, m_states, m_enter_id);
+            return pivot.choose_entering_edge();
+        }
+        else if (pr == BEST_ELIGIBLE) {
+            best_eligible_pivot pivot(m_graph, m_potentials, m_states, m_enter_id);
+            return pivot.choose_entering_edge();
+        }
+        else {
+            candidate_list_pivot pivot(m_graph, m_potentials, m_states, m_enter_id);
+            return pivot.choose_entering_edge();
+        } 
+    }
+
     // Minimize cost flows
     // Return true if found an optimal solution, and return false if unbounded
     template<typename Ext>
-    bool network_flow<Ext>::min_cost() {
+    bool network_flow<Ext>::min_cost(pivot_rule pr) {
         initialize();
-        while (choose_entering_edge()) { 
+        while (choose_entering_edge(pr)) { 
             bool bounded = choose_leaving_edge();
             if (!bounded) return false;
             update_flows();

src/smt/spanning_tree_def.h

@@ -162,7 +162,7 @@ namespace smt {
             tout << u << ", " << v << ") leaves\n";
         });
 
-        node old_pred = m_pred[q];        
+        node old_pred = m_pred[q]; 
         // Update stem nodes from q to v
         if (q != v) {
             for (node n = q; n != u; ) {
@@ -175,18 +175,23 @@ namespace smt {
                 old_pred = next_old_pred;
             }     
         }
-        else {
-            node x = get_final(p);
-            node y = m_thread[x];
-            node z = get_final(q);
-            node t = m_thread[get_final(v)];
-            node r = find_rev_thread(v);
-            m_thread[z] = y;
-            m_thread[x] = q;
-            m_thread[r] = t;
-        }
 
-        m_pred[q] = p;        
+        // Old threads: alpha -> q -*-> f(q) -> beta | p -*-> f(p) -> gamma
+        // New threads: alpha -> beta                | p -*-> f(p) -> q -*-> f(q) -> gamma
+
+        node f_p = get_final(p);
+        node f_q = get_final(q);
+        node alpha = find_rev_thread(q);
+        node beta = m_thread[f_q];
+        node gamma = m_thread[f_p];
+
+        if (q != gamma) {
+            m_thread[alpha] = beta;
+            m_thread[f_p] = q;
+            m_thread[f_q] = gamma;
+        }
+       
+        m_pred[q] = p; 
         m_tree[q] = enter_id;
         m_root_t2 = q;
 
@@ -211,7 +216,6 @@ namespace smt {
         Spanning tree of m_graph + root is represented using:
         
         svector<edge_state> m_states;      edge_id |-> edge_state
-        svector<bool> m_upwards;           node |-> bool
         svector<node> m_pred;              node |-> node
         svector<int>  m_depth;             node |-> int
         svector<node> m_thread;            node |-> node
@@ -224,9 +228,7 @@ namespace smt {
         m_thread is a linked list traversing all nodes.
         Furthermore, the nodes linked in m_thread follows a 
         depth-first traversal order.
-
+        
-        m_upwards direction of edge from i to m_pred[i] m_graph
-                
     */
     template<typename Ext>
     bool thread_spanning_tree<Ext>::check_well_formed() {