source: lemon-0.x/tools/lgf-gen.cc @ 2390:8450951a8e2d

#include <lemon/list_graph.h>
#include <lemon/graph_utils.h>
#include <lemon/random.h>
#include <lemon/dim2.h>
#include <lemon/bfs.h>
#include <lemon/counter.h>
#include <lemon/suurballe.h>
#include <lemon/graph_to_eps.h>
#include <lemon/graph_writer.h>
#include <lemon/arg_parser.h>
#include <cmath>
#include <algorithm>
#include <lemon/unionfind.h>

using namespace lemon;

typedef dim2::Point<double> Point;

UGRAPH_TYPEDEFS(ListUGraph);

int N;
int girth;

ListUGraph g;

std::vector<Node> nodes;
ListUGraph::NodeMap<Point> coords(g);

int tsp_impr_num=0;

const double EPSILON=1e-8; 
bool tsp_improve(Node u, Node v)
{
  double luv=std::sqrt((coords[v]-coords[u]).normSquare());
  Node u2=u;
  Node v2=v;
  do {
    Node n;
    for(IncEdgeIt e(g,v2);(n=g.runningNode(e))==u2;++e);
    u2=v2;
    v2=n;
    if(luv+std::sqrt((coords[v2]-coords[u2]).normSquare())-EPSILON>
       std::sqrt((coords[u]-coords[u2]).normSquare())+
       std::sqrt((coords[v]-coords[v2]).normSquare()))
      {
        g.erase(findUEdge(g,u,v));
        g.erase(findUEdge(g,u2,v2));
        g.addEdge(u2,u);
        g.addEdge(v,v2);
        tsp_impr_num++;
        return true;
      }
  } while(v2!=u);
  return false;
}

bool tsp_improve(Node u)
{
  for(IncEdgeIt e(g,u);e!=INVALID;++e)
    if(tsp_improve(u,g.runningNode(e))) return true;
  return false;
}

void tsp_improve()
{
  bool b;
  do {
    b=false;
    for(NodeIt n(g);n!=INVALID;++n)
      if(tsp_improve(n)) b=true;
  } while(b);
}

void tsp()
{
  for(int i=0;i<N;i++) g.addEdge(nodes[i],nodes[(i+1)%N]);
  tsp_improve();
}

class Line
{
public:
  Point a;
  Point b;
  Line(Point _a,Point _b) :a(_a),b(_b) {}
  Line(Node _a,Node _b) : a(coords[_a]),b(coords[_b]) {}
  Line(const Edge &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
  Line(const UEdge &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
};
  
inline std::ostream& operator<<(std::ostream &os, const Line &l)
{
  os << l.a << "->" << l.b;
  return os;
}

bool cross(Line a, Line b) 
{
  Point ao=rot90(a.b-a.a);
  Point bo=rot90(b.b-b.a);
  return (ao*(b.a-a.a))*(ao*(b.b-a.a))<0 &&
    (bo*(a.a-b.a))*(bo*(a.b-b.a))<0;
}

struct Pedge
{
  Node a;
  Node b;
  double len;
};

bool pedgeLess(Pedge a,Pedge b)
{
  return a.len<b.len;
}

std::vector<UEdge> edges;

void triangle()
{
  Counter cnt("Number of edges added: ");
  std::vector<Pedge> pedges;
  for(NodeIt n(g);n!=INVALID;++n) 
    for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
      {
        Pedge p;
        p.a=n;
        p.b=m;
        p.len=(coords[m]-coords[n]).normSquare();
        pedges.push_back(p);
      }
  std::sort(pedges.begin(),pedges.end(),pedgeLess);
  for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
    {
      Line li(pi->a,pi->b);
      UEdgeIt e(g);
      for(;e!=INVALID && !cross(e,li);++e) ;
      UEdge ne;
      if(e==INVALID) {
        ne=g.addEdge(pi->a,pi->b);
        edges.push_back(ne);
        cnt++;
      }
    }
}

void sparse(int d) 
{
  Counter cnt("Number of edges removed: ");
  Bfs<ListUGraph> bfs(g);
  for(std::vector<UEdge>::reverse_iterator ei=edges.rbegin();
      ei!=edges.rend();++ei)
    {
      Node a=g.source(*ei);
      Node b=g.target(*ei);
      g.erase(*ei);
      bfs.run(a,b);
      if(bfs.predEdge(b)==INVALID || bfs.dist(b)>d)
        g.addEdge(a,b);
      else cnt++;
    }
}

void sparse2(int d) 
{
  Counter cnt("Number of edges removed: ");
  for(std::vector<UEdge>::reverse_iterator ei=edges.rbegin();
      ei!=edges.rend();++ei)
    {
      Node a=g.source(*ei);
      Node b=g.target(*ei);
      g.erase(*ei);
      ConstMap<Edge,int> cegy(1);
      Suurballe<ListUGraph,ConstMap<Edge,int> > sur(g,cegy,a,b);
      int k=sur.run(2);
      if(k<2 || sur.totalLength()>d)
        g.addEdge(a,b);
      else cnt++;
//       else std::cout << "Remove edge " << g.id(a) << "-" << g.id(b) << '\n';
    }
}

void sparseTriangle(int d)
{
  Counter cnt("Number of edges added: ");
  std::vector<Pedge> pedges;
  for(NodeIt n(g);n!=INVALID;++n) 
    for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
      {
        Pedge p;
        p.a=n;
        p.b=m;
        p.len=(coords[m]-coords[n]).normSquare();
        pedges.push_back(p);
      }
  std::sort(pedges.begin(),pedges.end(),pedgeLess);
  for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
    {
      Line li(pi->a,pi->b);
      UEdgeIt e(g);
      for(;e!=INVALID && !cross(e,li);++e) ;
      UEdge ne;
      if(e==INVALID) {
        ConstMap<Edge,int> cegy(1);
        Suurballe<ListUGraph,ConstMap<Edge,int> >
          sur(g,cegy,pi->a,pi->b);
        int k=sur.run(2);
        if(k<2 || sur.totalLength()>d)
          {
            ne=g.addEdge(pi->a,pi->b);
            edges.push_back(ne);
            cnt++;
          }
      }
    }
}

void minTree() {
  std::vector<Pedge> pedges;
  for(NodeIt n(g);n!=INVALID;++n) 
    for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
      {
        Pedge p;
        p.a=n;
        p.b=m;
        p.len=(coords[m]-coords[n]).normSquare();
        pedges.push_back(p);
      }
  std::sort(pedges.begin(),pedges.end(),pedgeLess);
  ListUGraph::NodeMap<int> comp(g);
  UnionFind<ListUGraph::NodeMap<int> > uf(comp);
  for (NodeIt it(g); it != INVALID; ++it)
    uf.insert(it);

  int en=0;
  for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
    {
      if ( uf.join(pi->a,pi->b) ) {
        g.addEdge(pi->a,pi->b);
        en++;
        if(en>=N-1) return;
      }
    }
}



int main(int argc,char **argv) 
{
  ArgParser ap(argc,argv);

  bool eps;
  bool disc_d, square_d, gauss_d;
  bool tsp_a,two_a,tree_a;
  int num_of_cities=1;
  double area=1;
  N=100;
  girth=10;
  std::string ndist("disc");
  ap.option("n", "Number of nodes (default is 100)", N)
    .option("g", "Girth parameter (default is 10)", girth)
    .option("cities", "Number of cities (default is 1)", num_of_cities)
    .option("area", "Full relative area of the cities (default is 1)", area)
    .option("disc", "Nodes are evenly distributed on a unit disc (default)",disc_d)
    .optionGroup("dist", "disc")
    .option("square", "Nodes are evenly distributed on a unit square", square_d)
    .optionGroup("dist", "square")
    .option("gauss",
            "Nodes are located according to a two-dim gauss distribution",
            gauss_d)
    .optionGroup("dist", "gauss")
//     .mandatoryGroup("dist")
    .onlyOneGroup("dist")
    .option("eps", "Also generate .eps output (prefix.eps)",eps)
    .option("2con", "Create a two connected planar graph",two_a)
    .optionGroup("alg","2con")
    .option("tree", "Create a min. cost spanning tree",tree_a)
    .optionGroup("alg","tree")
    .option("tsp", "Create a TSP tour",tsp_a)
    .optionGroup("alg","tsp")
    .onlyOneGroup("alg")
    .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
    .run();
  
  std::string prefix;
  switch(ap.files().size()) 
    {
    case 0:
      prefix="lgf-gen-out";
      break;
    case 1:
      prefix=ap.files()[0];
      break;
    default:
      std::cerr << "\nAt most one prefix can be given\n\n";
      exit(1);
    }
  
  double sum_sizes=0;
  std::vector<double> sizes;
  std::vector<double> cum_sizes;
  for(int s=0;s<num_of_cities;s++) 
    {
      //        sum_sizes+=rnd.exponential();
      double d=rnd();
      sum_sizes+=d;
      sizes.push_back(d);
      cum_sizes.push_back(sum_sizes);
    }
  int i=0;
  for(int s=0;s<num_of_cities;s++) 
    {
      Point center=(num_of_cities==1?Point(0,0):rnd.disc());
      if(gauss_d)
        for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
          Node n=g.addNode();
          nodes.push_back(n);
          coords[n]=center+rnd.gauss2()*area*
            std::sqrt(sizes[s]/sum_sizes);
        }
      else if(square_d)
        for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
          Node n=g.addNode();
          nodes.push_back(n);
          coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*
            std::sqrt(sizes[s]/sum_sizes);
        }
      else if(disc_d || true)
        for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
          Node n=g.addNode();
          nodes.push_back(n);
          coords[n]=center+rnd.disc()*area*
            std::sqrt(sizes[s]/sum_sizes);
        }
    }
  
  if(tsp_a) {
    tsp();
    std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
  }
  else if(two_a) {
    std::cout << "Make triangles\n";
    //   triangle();
    sparseTriangle(girth);
    std::cout << "Make it sparser\n";
    sparse2(girth);
  }
  else if(tree_a) {
    minTree();
  }
  

  std::cout << "Number of nodes    : " << countNodes(g) << std::endl;
  std::cout << "Number of edges    : " << countUEdges(g) << std::endl;
  double tlen=0;
  for(UEdgeIt e(g);e!=INVALID;++e)
    tlen+=sqrt((coords[g.source(e)]-coords[g.target(e)]).normSquare());
  std::cout << "Total edge length  : " << tlen << std::endl;
  if(eps)
    graphToEps(g,prefix+".eps").
      scale(600).nodeScale(.2).edgeWidthScale(.001).preScale(false).
      coords(coords).run();

  UGraphWriter<ListUGraph>(prefix+".lgf",g).
    writeNodeMap("coordinates_x",scaleMap(xMap(coords),600)).
    writeNodeMap("coordinates_y",scaleMap(yMap(coords),600)).
    run();
}