mcf_new.t

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

/*******************************************************************************
+
+  LEDA 4.5  
+
+
+  mcf_new.t
+
+
+  Copyright (c) 1995-2004
+  by Algorithmic Solutions Software GmbH
+  All rights reserved.
+ 
*******************************************************************************/

// $Revision: 1.4 $  $Date: 2004/02/06 11:20:18 $

#ifndef MINCOSTFLOW_NEW_H
#define MINCOSTFLOW_NEW_H

#include <LEDA/graph.h>
#include <LEDA/graph_alg.h>

LEDA_BEGIN_NAMESPACE

namespace mincostflow_new
{
  namespace mcf_detail
  {

    template <class E> struct limits {
      static E min() throw() { return E(); }
    };

    template <> struct limits <double> {
      static double min() { return -MAXDOUBLE; }
    };

    template <> struct limits <int> {
      static int min() { return -MAXINT; }
    };


    template <class E, class GraphObj, int slot>
    class cell
    {
      E*       p;
      GraphObj obj;

      public:

      static EVENT2<GraphObj,E> read_event;
      static EVENT2<GraphObj,E> write_event;

      cell(E& x, GraphObj o) : p(&x), obj(o) {}

      operator E& () const
      {
        read_event(obj,*p);
        return *p;
      }

      cell& operator=(const E& x)
      {
        *p = x;
        write_event(obj,*p);
        return *this;
      }
    };

    template <class E, class GraphObj, int slot>
    EVENT2<GraphObj,E> cell<E,GraphObj,slot>::read_event;

    template <class E,class GraphObj, int slot>
    EVENT2<GraphObj,E> cell<E,GraphObj,slot>::write_event;


    template <class E, class GraphObj, int slot, int base_sz>
    struct advance_cell_access
    {
      const cell<E,GraphObj,slot> operator[](GraphObj o) const
      { return cell<E,GraphObj,slot>(LEDA_CONST_ACCESS(E,
o->data(base_sz + slot)),o); }

      cell<E,GraphObj,slot>       operator[](GraphObj o)
      { return cell<E,GraphObj,slot>(LEDA_ACCESS(E, o->data(base_sz +
slot)),o); }
    };



    template <class E, class GraphObj, int slot, int base_sz>
    struct std_cell_access
    {
      const E& operator[](GraphObj o) const
      { return LEDA_CONST_ACCESS(E, o->data(base_sz + slot)); }

      E&       operator[](GraphObj o)
      { return LEDA_ACCESS(E, o->data(base_sz + slot)); }
    };


    template <class E, int slot, class GT = graph,
              class ACC = std_cell_access<E, typename GT::node, slot,
GT::n_base_sz> >
    struct node_array : public ACC
    {
      typedef typename GT::node node;

      void init(node v, const E& x)
      { v->data(GT::n_base_sz + slot) = leda_copy(x); }
    };


    template <class E, int slot, class GT = graph,
              class ACC = std_cell_access<E, typename GT::edge, slot,
GT::e_base_sz> >
    struct edge_array : public ACC
    {
      typedef typename GT::edge edge;

      void init(edge e, const E& x)
      { e->data(GT::e_base_sz + slot) = leda_copy(x); }
    };


    template <int slot, class GT = graph>
    class node_queue
    {
      typedef typename GT::node node;

      node_array<node,slot,GT> NEXT;

      node head;
      node tail;

      public:

#ifdef MCF_EVENTS
      static EVENT1<bool>      is_empty_event;
      static EVENT1<node>      append_event;
      static EVENT1<node>      pop_event;
      static EVENT2<node,bool> is_member_event;
#endif

      node_queue() : head(0), tail(0) {}

      bool empty() const
      {

#ifdef MCF_EVENTS
        is_empty_event(head == 0);
#endif
        return head == 0;
      }

      node pop()
      {
        node v = head;
        head = NEXT[v];
        if (head == v) head = tail = 0;
        NEXT[v] = 0;

#ifdef MCF_EVENTS
        pop_event(v);
#endif
        return v;
      }

      void append(node v)
      {

#ifdef MCF_EVENTS
        append_event(v);
#endif

        if (tail == 0)
          head = tail = NEXT[v] = v;
        else
        {
          NEXT[tail] = NEXT[v] = v;
          tail = v;
        }
      }

      bool is_member(node v) const
      {
#ifdef MCF_EVENTS
        is_member_event(v, NEXT[v]);
#endif
        return NEXT[v];
      }
    };


#ifdef MCF_EVENTS
    template <int slot, class GT>
    EVENT1<bool> node_queue<slot,GT>::is_empty_event;

    template <int slot, class GT>
    EVENT1<typename GT::node> node_queue<slot,GT>::pop_event;

    template <int slot, class GT>
    EVENT1<typename GT::node> node_queue<slot,GT>::append_event;

    template <int slot, class GT>
    EVENT2<typename GT::node, bool>
node_queue<slot,GT>::is_member_event;
#endif

  }; /* end of namespace mcf_detail */


  template <class T, class GT = graph, int nslot = 0, int eslot = 0>
  class mcf
  {
    typedef double PT;

    public:

    typedef typename GT::node node;
    typedef typename GT::edge edge;
    typedef GT graph_t;
    typedef T  value_t;
    typedef PT pot_t;

    private:

#ifndef MCF_EVENTS

    typedef mcf_detail::node_array<PT,  nslot,  GT> pot_array;
    typedef mcf_detail::node_array<T,   nslot+1,GT> exc_array;
    typedef mcf_detail::node_array<edge,nslot+2,GT> cur_array;

    typedef mcf_detail::node_queue<nslot+3,GT> node_queue;

    typedef mcf_detail::edge_array<T,eslot,     GT> cap_array;
    typedef mcf_detail::edge_array<T,eslot+1,   GT> flow_array;
    typedef mcf_detail::edge_array<T,eslot+2,   GT> cost_array;
    typedef mcf_detail::edge_array<bool,eslot+3,GT> fix_array;

#else

    typedef mcf_detail::node_array<PT,  nslot,
GT,mcf_detail::advance_cell_access<PT,  node,nslot,  GT::n_base_sz> >
pot_array;
    typedef mcf_detail::node_array<T,
nslot+1,GT,mcf_detail::advance_cell_access<T,
node,nslot+1,GT::n_base_sz> > exc_array;
    typedef
mcf_detail::node_array<edge,nslot+2,GT,mcf_detail::advance_cell_access<edge,node,nslot+2,GT::n_base_sz>
> cur_array;

    typedef mcf_detail::node_queue<nslot+3,GT> node_queue;

    typedef mcf_detail::edge_array<T,eslot,
GT,mcf_detail::advance_cell_access<T,edge,eslot,  GT::e_base_sz> >
cap_array;
    typedef
mcf_detail::edge_array<T,eslot+1,GT,mcf_detail::advance_cell_access<T,edge,eslot+1,GT::e_base_sz>
> flow_array;
    typedef
mcf_detail::edge_array<T,eslot+2,GT,mcf_detail::advance_cell_access<T,edge,eslot+2,GT::e_base_sz>
> cost_array;
    typedef
mcf_detail::edge_array<T,eslot+3,GT,mcf_detail::advance_cell_access<T,edge,eslot+3,GT::e_base_sz>
> fix_array;

#endif

    const node_array<T,GT>& bp;
    const edge_array<T,GT>& lp;
    const edge_array<T,GT>& up;
    const edge_array<T,GT>& cp;

    GT& G;

    pot_array  pot;
    exc_array  exc;
    cur_array  cur;

    cap_array  cap;
    flow_array flow;
    cost_array cost;
    fix_array  fix;

    node_queue active;

    double epsilon;
    double scale_factor;

    int discharge_cnt[16];
    int relabel_cnt[16];
    int fixed_cnt[16];
    int push_cnt[16];
    int refine_cnt;

    float init_time;

    public:

    enum heuristics
    {
      lookahead  = (1 << 0),
      edgefixing = (1 << 1)
    };


#ifdef MCF_EVENTS
    static EVENT1<mcf_detail::cs_data< mcf<T,GT,nslot,eslot> >&>
start_event;
    static EVENT0
end_event;
    static EVENT1<double>
reduce_epsilon_event;
    static EVENT3<edge,node,double>
check_reduced_cost_event;
#endif


    class error
    {
      char* m;

      public:

      error(char* m0) : m(m0) {}
      char* message() const { return m; }
    };


    mcf(GT& G0,
        const node_array<T,GT>& b,
        const edge_array<T,GT>& c, const edge_array<T,GT>& l,
        const edge_array<T,GT>& u) : G(G0), bp(b), cp(c), lp(l), up(u)
{}


    void cost_scaling(edge_array<T,GT>& f, double sc, bool check =
false)
    {
      float t = used_time();

      scale_factor = sc;
      init_network();

      if (check)
        if (!check_feasibility())
          throw(mcf<T,GT>::error("mcf<T,GT>::cost_scaling(): problem isn't feasible!"));

      init_time = used_time(t);

      double bound = 1 / double(G.number_of_nodes());
      while (epsilon > bound)
      {
        //reduce_eps();
      epsilon /= scale_factor;
        refine();
      }

      f.init(G);

      edge e;
      forall_edges(e,G)
        f[e] = flow[e] + lp[e];
    }


/*
    void cost_scaling_infos(ostream& os = cout)
    {
      os << "c cost scaling version 1.0\n";
      os << "c\n";
      os << "c heuristics: push lookahead, edge fixing\n";
      os << "c features:   advance data access, own excess queue, event
support\n";
      os << "c\n";
      os << "c init.time         " << string("%10.2f\n",init_time);
      os << "c scale-factor      " << string("%10.0f\n",scale_factor);
      os << "c refines           " << string("%10d\n",refine_cnt);
      os << "c\n";
      os << "c ";
      os.width(15); os << "discharges";
      os.width(15); os << "relabels";
      os.width(15); os << "pushes";
      os.width(20); os << "fixed edges\n";
      os << "c\n";

      int sum[4], i;
      for (i=0; i<4; i++) sum[i] = 0;

      for (i=0; i<refine_cnt; i++)
      {
        os << "c ";
        os.width(15); os << string("%d",discharge_cnt[i]);
        os.width(15); os << string("%d",relabel_cnt[i]);
        os.width(15); os << string("%d",push_cnt[i]);
        os.width(20); os << string("%d\n",fixed_cnt[i]);

        sum[0] += discharge_cnt[i];
        sum[1] += relabel_cnt[i];
        sum[2] += push_cnt[i];
        sum[3] += fixed_cnt[i];
      }

      os << "c\n";
      os << "c ";
      os.width(15); os << string("%d",sum[0]);
      os.width(15); os << string("%d",sum[1]);
      os.width(15); os << string("%d",sum[2]);
      os.width(20); os << string("%d\n",sum[3]);

      os << "\n";
    }
*/


    private:


    void reduce_eps()
    {
      epsilon /= scale_factor;

#ifdef MCF_EVENTS
      reduce_epsilon_event(epsilon);
#endif

    }



    bool is_admissible(edge e, node v)
    {
      node s = source(e);
      node t = target(e);

      T  rc;
      PT cp;