mcf_new.t

A Library of Efficient Data Types and Algorithms,封装了常用的ADT及其相关算法的软件包

      if (v == s)
      {
        rc = cap[e] - flow[e];
        cp = cost[e] + pot[s] - pot[t];
      }
      else
      {
        rc = flow[e];
        cp = -cost[e] + pot[t] - pot[s];
      }

      return (rc > 0 && cp < 0) ? true : false;
    }




    void refine()
    {
      edge e;
      forall_edges(e,G)
      {
        if (fix[e] == true) continue;

        node v = source(e);
        node w = target(e);

        if (cost[e] + pot[v] - pot[w] < -epsilon)
        {
          T df = cap[e] - flow[e];

          exc[v]  -= df;
          exc[w]  += df;
          flow[e] += df;

          push_cnt[refine_cnt]++;
        }
        else
        {
          T df = flow[e];

          exc[v] += df;
          exc[w] -= df;
          flow[e] = 0;

          push_cnt[refine_cnt]++;
        }
      }

      node v;
      forall_nodes(v,G)
        if (exc[v] > 0) active.append(v);

      while (!active.empty())
      {
        node v = active.pop();
        if (exc[v] > 0) discharge(v);
      }

      edge_fixing();

      refine_cnt++;
    }




    void discharge(node v)
    {
      edge e = cur[v];
      if (!is_admissible(e,v)) relabel(v);

      while (exc[v] > 0)
      {
        edge e = cur[v];
        node w = G.opposite(v,e);

        if (exc[w] >= 0)
        {
          if (is_admissible(cur[w],w) || !relabel(w))
          {
            if (v == source(e))
            {
              T df = cap[e] - flow[e];
              if (df > exc[v]) df = exc[v];

              exc[v]  -= df;
              exc[w]  += df;
              flow[e] += df;
            }
            else
            {
              T df = flow[e];
              if (df > exc[v]) df = exc[v];

              exc[v]  -= df;
              exc[w]  += df;
              flow[e] -= df;
            }

            if (!active.is_member(w)) active.append(w);

            push_cnt[refine_cnt]++;
          }
          else
          {
            if (exc[w] > 0)
            {
              if (v == source(e))
              {
                T df = flow[e];
                if (df > exc[w]) df = exc[w];

                exc[w]  -= df;
                exc[v]  += df;
                flow[e] -= df;

                push_cnt[refine_cnt]++;
              }
              else
              {
                T df = cap[e] - flow[e];
                if (df > exc[w]) df = exc[w];

                exc[w]  -= df;
                exc[v]  += df;
                flow[e] += df;

                push_cnt[refine_cnt]++;
              }
            }
          }
        }
        else
        {
          if (v == source(e))
          {
            T df = cap[e] - flow[e];
            if (df > exc[v]) df = exc[v];

            exc[v]  -= df;
            exc[w]  += df;
            flow[e] += df;
          }
          else
          {
            T df = flow[e];
            if (df > exc[v]) df = exc[v];

            exc[v]  -= df;
            exc[w]  += df;
            flow[e] -= df;
          }

          if (!active.is_member(w)) active.append(w);

          push_cnt[refine_cnt]++;

          if (exc[w] > 0 && !is_admissible(cur[w],w)) relabel(w);
//          if (exc[w] > 0) relabel(w);
        }

        if (exc[v] == 0)
        {
          discharge_cnt[refine_cnt]++;
          return;
        }

        relabel(v);
      }
    }



    bool relabel(node v)
    {
      PT max_pot = mcf_detail::limits<PT>::min();
      PT   cur_p = max_pot;
      edge cur_e = cur[v];

      fix[cur[v]] = true;

      edge e;
      forall_out_edges(e,v)
      {
        if (fix[e] == true) continue;

        if (flow[e] < cap[e])          /* residual cap(e) > 0 */
        {
          node w = target(e);
          PT  dp = pot[w];

          if (cost[e] + pot[v] < dp)   /* reduced cost(e,v) < 0? */
          {
            fix[cur[v]] = false;

            cur[v] = e;
            return false;
          }

          dp -= cost[e];

          if (dp > cur_p)
          {
            cur_p = dp;
            cur_e = e;
          }
        }
      }


      forall_in_edges(e,v)
      {
        if (fix[e] == true) continue;

        if (flow[e] > 0)               /* residual cap(e) > 0 */
        {
          node w = G.source(e);
          PT  dp = pot[w];

          if (-cost[e] - pot[v] < dp)  /* reduced cost(e,v) < 0? */
          {
            fix[cur[v]] = false;

            cur[v] = e;
            return false;
          }

          dp += cost[e];

          if (dp > cur_p)
          {
            cur_p = dp;
            cur_e = e;
          }
        }
      }


      e = cur[v];
      fix[e] = false;

      if (v == source(e))
      {
        if (flow[e] < cap[e])
        {
          node w = target(e);
          PT  dp = pot[w];

          if (cost[e] + pot[v] < dp) return false;

          dp -= cost[e];

          if (dp > cur_p)
          {
            cur_p = dp;
            cur_e = e;
          }
        }
      }
      else
      {
        if (flow[e] > 0)
        {
          node w = source(e);
          PT  dp = pot[w];

          if (-cost[e] - pot[v] < dp) return false;

          dp += cost[e];

          if (dp > cur_p)
          {
            cur_p = dp;
            cur_e = e;
          }
        }
      }


      if (cur_p > max_pot)
      {
        cur[v] = cur_e;
        pot[v] = cur_p - epsilon;
      }
      else
      {
        if (exc[v] == 0)
          pot[v] = max_pot;
        else
          throw(mcf<T,GT>::error("mcf<T,GT>::relabel(): problem isn't feasible"));
      }

      relabel_cnt[refine_cnt]++;

      return true;
    }



    void edge_fixing()
    {
      if (refine_cnt < 2) return;

      double fixing_bound = 2 * G.number_of_nodes() * epsilon;

      edge e;
      forall_edges(e,G)
      {
        if (fix[e] == true) continue;

        PT ps = pot[source(e)];
        PT pt = pot[target(e)];

        if (abs(cost[e] + ps - pt) > fixing_bound)
        {
          if (abs(-cost[e] + pt - ps) > fixing_bound)
          {
            fix[e] = true;
            fixed_cnt[refine_cnt]++;
          }
        }
      }
    }



    void init_network()
    {
      refine_cnt = 0;

      for(int i=0; i<16; i++)
      {
        discharge_cnt[i] = 0;
        relabel_cnt[i]   = 0;
        push_cnt[i]      = 0;
        fixed_cnt[i]     = 0;
      }

      node v;
      forall_nodes(v,G)
      {
        pot.init(v,0);
        exc.init(v,bp[v]);
        edge e = G.first_adj_edge(v);
        cur.init(v, (e) ? e : G.first_in_edge(v));
      }

      T C = 0;

      edge e;
      forall_edges(e,G)
      {
        cap.init(e,up[e] - lp[e]);
        flow.init(e,0);
        cost.init(e,cp[e]);
        fix.init(e,false);

        if (cost[e] > C) C = cost[e];
      }

      epsilon = 1;
      while (epsilon < C) epsilon *= 2;
    }



    bool check_feasibility()
    {
/*      node s = G.new_node();
      node t = G.new_node();

      list<edge> s_edges;
      list<edge> t_edges;

      node v;
      forall_nodes(v,G)
      {
        if (v == s || v == t || exc[v] == 0) continue;

        if (exc[v] > 0)
          s_edges.append(G.new_edge(s,v));
        else
          t_edges.append(G.new_edge(v,t));
      }

      edge_array<T,GT> maxflow(G,0);
      edge_array<T,GT> maxcap(G,0);

      edge e;
      forall_edges(e,G) maxcap[e] =  cap[e];
      forall(e,s_edges) maxcap[e] =  exc[G.target(e)];
      forall(e,t_edges) maxcap[e] = -exc[G.source(e)];

      MAX_FLOW(G,s,t,maxcap,maxflow);

      forall(e,s_edges)
        if (maxcap[e] != maxflow[e])
        {
          G.del_node(s);
          G.del_node(t);

          return false;
        }

      G.del_node(s);
      G.del_node(t);
*/
      return true;
    }

  };

#ifdef MCF_EVENTS
  template <class T, class GT, int nslot, int eslot>
  EVENT1<mcf_detail::cs_data< mcf<T,GT,nslot,eslot> >&>
mcf<T,GT,nslot,eslot>::start_event;

  template <class T, class GT, int nslot, int eslot>
  EVENT0 mcf<T,GT,nslot,eslot>::end_event;

  template <class T, class GT, int nslot, int eslot>
  EVENT1<double> mcf<T,GT,nslot,eslot>::reduce_epsilon_event;

  template <class T, class GT, int nslot, int eslot>
  EVENT3<typename GT::edge,typename GT::node,double>
mcf<T,GT,nslot,eslot>::check_reduced_cost_event;
#endif

}; /* end namespace mincostflow_new */


LEDA_END_NAMESPACE

#endif