partitioning.c

#include<LEDA/graph.h>
#include<LEDA/graphwin.h>
#include<LEDA/graph_misc.h>
#include<LEDA/graph_alg.h>
#include<LEDA/delaunay.h>

#include<scil/scil.h>

using namespace LEDA;
using namespace SCIL;

#define RHO 0.4

class mysep : public sym_constraint {
  graph& G;
  var_map<node>& VM;
  var_map<edge>& EM;

  public:
  mysep(graph& G_, var_map<node>& VM_, var_map<edge>& EM_)
    : G(G_), VM(VM_), EM(EM_) {
  };

  status separate(subproblem& s) {
    status sta=no_constraint_found;
/*
    edge_array<double> E(G);
    node_array<double> C(G);
    node_array<edge> P(G);
    edge e,f;
    forall_edges(e, G) {
      E[e]=s.value(EM[e]);
      if(E[e]<=0) E[e]=0;
    };

    G.sort_edges(E);

    list<edge> L=MIN_SPANNING_TREE(G,E);

    forall(e, L) {
      G.hide_edge(e);
    }

    list<edge> R=G.all_edges();

    forall(e, R) G.hide_edge(e);
    forall(e, L) G.restore_edge(e);

    edge_array<int> Cp(G);
    node_array<bool> Co(G);
    int sz;
    double d;
    node v;

    forall(e,R) {
      G.restore_edge(e);
 
      BICONNECTED_COMPONENTS(G, Cp);

      forall_nodes(v, G) Co[v]=false;
      sz=0; d=0;
      row r;
      forall_edges(f,G) if(Cp[f]==Cp[e]) {
        if(Co[source(f)]==false) { sz++; Co[source(f)]=true; }
        if(Co[target(f)]==false) { sz++; Co[target(f)]=true; }
        r+=EM[f]; d+=E[f];
      }

      if ((sz>RHO*G.number_of_nodes()) && (d<1.999)) {
        s.add_basic_constraint(r>=2);
        sta=constraint_found;
      }
    }
*/
    return sta;
  };
};

void MY_DFS(graph& G, node s, node_array<bool>& R, edge_array<int>& C, edge_array<int>& F) {
  R[s]=true;
  edge e;
  forall_out_edges(e,s) if(!R[target(e)]) {
    if(C[e]!=F[e]) MY_DFS(G,target(e),R,C,F);
  }
  forall_in_edges(e,s) if(!R[source(e)]) {
    if(F[e]>0) MY_DFS(G,source(e),R,C,F);
  }
};

/*
class myheuristic : public primal_heuristic {
private:
  graph& G;
  node r;
  var_map<edge>& EM;
  var_map<node>& NM;
  edge_array<double>& Cost;

public:
  myheuristic(graph& G_, var_map<edge>& EM_, var_map<node>& NM_, edge_array<double>& Cost_, node r_) 
  : G(G_), EM(EM_), NM(NM_), Cost(Cost_), r(r_) {
  };

  int heuristic(subproblem& Sub) {
    node_array<double> D(G);

    node v;
    forall_nodes(v, G) D[v]=Sub.value(NM[v]);

    node_array<int> D1(G);
    BFS(G,r,D1);
    forall_nodes(v, G) D[v]!=((double) D1[v])/1000.0;

    G.sort_nodes(D);

    int i=0;
    graph H;
    node s=H.new_node();
    node t=H.new_node();
    node_array<node> GtoH(G);
    forall_nodes(v, G) {
      if(i<RHO*G.number_of_nodes()) GtoH[v]=s;
      else {
        if(i>(1-RHO)*G.number_of_nodes()) GtoH[v]=t;
        else GtoH[v]=H.new_node();
      }
      i++;
    }

    edge_array<int> C(H,2*G.number_of_edges(), 0);
    edge e;
    forall_edges(e, G) { 
      C[H.new_edge(GtoH[source(e)], GtoH[target(e)])]=(int) 10000.0*Cost[e];
      C[H.new_edge(GtoH[target(e)], GtoH[source(e)])]=(int) 10000.0*Cost[e];
    }

    edge_array<int> F(H);
    MAX_FLOW(H,s,t,C,F);

    node_array<bool> R(H, false);
    MY_DFS(H,s,R,C,F);

    solution S;
//    forall_nodes(v, G) if(!R[GtoH[v]]) S.set_value(NM[v],1);
 //   forall_edges(e, G) if(R[GtoH[source(e)]]!=R[GtoH[target(e)]]) S.set_value(EM[e],1);
    Sub.save_solution(S);
  };
};
*/

void compute_partitioning(graph& G, edge_array<double>& C, list<node>& L) {

  G.make_undirected();

  ILP_Problem IP(Optsense_Min);

  var_map<node> VM;
  node u;
  row r;
  forall_nodes(u, G) {
    VM[u]=IP.add_variable(0,0,1,Vartype_Integer);
    r+=VM[u];
  };

  u=G.choose_node();
  IP.add_basic_constraint((row) VM[u]==0);
  IP.add_basic_constraint(r>=(int) (RHO*G.number_of_nodes()));
  IP.add_basic_constraint(r<=(int) ((1-RHO)*G.number_of_nodes()));

  var_map<edge> EM;
  edge e;
  forall_edges(e,G) {
    EM[e]=IP.add_variable(1/C[e],0,1,Vartype_Integer);
    IP.add_basic_constraint((row) EM[e]>=(row) VM[source(e)]-VM[target(e)]);
    IP.add_basic_constraint((row) EM[e]>=(row) VM[target(e)]-VM[source(e)]);
  }

  //IP.add_sym_constraint(new mysep(G, VM, EM));

  //IP.set_primal_heuristic(new myheuristic(G, EM, VM, C, u));

  IP.optimize();

  forall_nodes(u, G) {
    if(IP.get_solution(VM[u])>0.5) L.append(u);
  }

};

void compute_partitioning_2(graph& G, edge_array<double>& C, list<node>& L) {
  G.make_undirected();

  ILP_Problem IP(Optsense_Min);

  var_map<node> SM;
  var_map<node> TM;
  node u;
  row rs;
  row rt;
  forall_nodes(u, G) {
    SM[u]=IP.add_variable(-1,0,1,Vartype_Integer);
    rs+=SM[u];
    TM[u]=IP.add_variable(-1,0,1,Vartype_Integer);
    rt+=TM[u];
    IP.add_basic_constraint(SM[u]+TM[u]<=1);
  };

cout<<"l\n";
  u=G.choose_node();
  IP.add_basic_constraint((row) SM[u]-TM[u]<=0);
cout<<"d\n";
  IP.add_basic_constraint(rs<=(int) ((1-RHO)*G.number_of_nodes()));
  IP.add_basic_constraint(rt<=(int) ((1-RHO)*G.number_of_nodes()));

cout<<"Hier\n";
  edge e;
  forall_edges(e,G) {
    IP.add_basic_constraint(SM[source(e)]+TM[target(e)]<=1);
    IP.add_basic_constraint(TM[source(e)]+SM[target(e)]<=1);
  }

cout<<"done\n";

  IP.optimize();

  forall_nodes(u, G) {
    if(IP.get_solution(SM[u])+IP.get_solution(TM[u])<0.5) L.append(u);
  }

};


void triang_points(GraphWin& GW) 
{
  graph& G=GW.get_graph();
  int n=G.number_of_nodes();
  if (n!=0) {
  
    GW.set_flush(false);

    point_set T;  
    node_array<node> tr(T,n,nil);
    node_array<int> num(T,n,0);
    node v,v1;
    int i=0;

    double x,y;
    forall_nodes(v1,G) 
    {
      x=GW.get_position(v1).xcoord();
      y=GW.get_position(v1).ycoord();
      point p(x,y);
      v=T.insert(p);
      tr[v]=v1;
      num[v]=i++;
     }

    G.del_all_edges();
    GW.update_graph();
    edge e;
    forall_edges(e,T) if (num[T.source(e)]<num[T.target(e)])
                        GW.new_edge(tr[T.source(e)],tr[T.target(e)]);
    GW.redraw();
    GW.set_flush(true);
    GW.redraw();
  }
}

void generate_del(GraphWin& GW) 
{
  GW.set_flush(false);

  graph& G=GW.get_graph();
  if (!G.empty()) GW.clear_graph();

  int n=10;

  panel p;
  p.int_item("Number of nodes : ",n);
  p.open(GW.get_window());
 
  random_source ran;

  double x,y;
  for(int i=0;i<n;i++) 
  { ran >> x >> y;
    x=(x+0.08)*(GW.get_xmax()-GW.get_xmin())*0.9+GW.get_xmin();
    y=(y+0.08)*(GW.get_ymax()-GW.get_ymin())*0.9+GW.get_ymin();
    point p1(x,y);
    GW.new_node(p1);
   }
  triang_points(GW);

  GW.set_flush(true);
  GW.redraw();
}

int main()
{
  GraphWin GW;

// Setup of the graph window
  GW.set_node_width(8);
  GW.set_node_height(8);
  GW.set_flush(true);
  GW.set_edge_direction(undirected_edge);

// Additional functionality of the graph window; the implementation is hidden
  int gn=gw_add_menu(GW,"Generation");
  gw_add_simple_call(GW,triang_points,"Triangulate defined points",gn);
  gw_add_simple_call(GW,generate_del, "Generate random Delaunay-Graph",gn);

  GW.display();

  while(GW.edit())
  {
    graph& G=GW.get_graph();

    edge e;
    edge_array<double> C(G);
    forall_edges(e, G) 
      C[e]=GW.get_position(source(e)).distance(GW.get_position(target(e)));

    list<node> T;
    compute_partitioning_2(G, C, T);

    GW.set_flush(false);
    node u;
    forall_nodes(u,G) {
      GW.set_color(u, black);
    }
    forall(u,T) { 
      GW.set_color(u,red);
    }
    GW.set_flush(true);
    GW.redraw();
  }
  return 0;
}

---