COIN-OR::LEMON - Graph Library

Ticket #45: 88da4852da9e.patch

File 88da4852da9e.patch, 25.2 KB (added by Alpar Juttner, 10 years ago)
  • tools/Makefile.am

    # HG changeset patch
    # User Alpar Juttner <alpar@cs.elte.hu>
    # Date 1228241331 0
    # Node ID 88da4852da9e3a01312688edbfc4d5af35139796
    # Parent  f5eaf04b41c51005da9462da2f5c015741bcbce9
    Port lgf-gen from SVN -r3512 (#45)
    
     - apply the migrate script
     - apply the source unifyer
     - fix the compilation
    
    diff --git a/tools/Makefile.am b/tools/Makefile.am
    a b  
    11if WANT_TOOLS
    22
    33bin_PROGRAMS += \
    4         tools/dimacs-to-lgf
     4        tools/dimacs-to-lgf \
     5        tools/lgf-gen
    56
    67dist_bin_SCRIPTS += tools/lemon-0.x-to-1.x.sh
    78
    89endif WANT_TOOLS
    910
    1011tools_dimacs_to_lgf_SOURCES = tools/dimacs-to-lgf.cc
     12tools_lgf_gen_SOURCES = tools/lgf-gen.cc
  • new file tools/lgf-gen.cc

    diff --git a/tools/lgf-gen.cc b/tools/lgf-gen.cc
    new file mode 100644
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19///\ingroup tools
     20///\file
     21///\brief Special plane digraph generator.
     22///
     23///Digraph generator application for various types of plane digraphs.
     24///
     25///\verbatim
     26/// Usage:
     27///   ./tools/lgf-gen [-2con|-tree|-tsp|-tsp2|-dela] [-disc|-square|-gauss]
     28///      [-rand|-seed int] [--help|-h|-help] [-area num] [-cities int] [-dir]
     29///      [-eps] [-g int] [-n int] [prefix]
     30/// Where:
     31///   [prefix]
     32///      Prefix of the output files. Default is 'lgf-gen-out'
     33///   --help|-h|-help
     34///      Print a short help message
     35///   -2con
     36///      Create a two connected planar digraph
     37///   -area num
     38///      Full relative area of the cities (default is 1)
     39///   -cities int
     40///      Number of cities (default is 1)
     41///   -dela
     42///      Delaunay triangulation digraph
     43///   -dir
     44///      Directed digraph is generated (each arcs are replaced by two directed ones)
     45///   -disc
     46///      Nodes are evenly distributed on a unit disc (default)
     47///   -eps
     48///      Also generate .eps output (prefix.eps)
     49///   -g int
     50///      Girth parameter (default is 10)
     51///   -gauss
     52///      Nodes are located according to a two-dim gauss distribution
     53///   -n int
     54///      Number of nodes (default is 100)
     55///   -rand
     56///      Use time seed for random number generator
     57///   -seed int
     58///      Random seed
     59///   -square
     60///      Nodes are evenly distributed on a unit square
     61///   -tree
     62///      Create a min. cost spanning tree
     63///   -tsp
     64///      Create a TSP tour
     65///   -tsp2
     66///      Create a TSP tour (tree based)
     67///\endverbatim
     68/// \image html plane_tree.png
     69/// \image latex plane_tree.eps "Eucledian spanning tree" width=\textwidth
     70///
     71
     72
     73#include <algorithm>
     74#include <set>
     75#include <lemon/list_graph.h>
     76#include <lemon/random.h>
     77#include <lemon/dim2.h>
     78#include <lemon/bfs.h>
     79#include <lemon/counter.h>
     80#include <lemon/suurballe.h>
     81#include <lemon/graph_to_eps.h>
     82#include <lemon/lgf_writer.h>
     83#include <lemon/arg_parser.h>
     84#include <lemon/euler.h>
     85#include <lemon/math.h>
     86#include <lemon/kruskal.h>
     87#include <lemon/time_measure.h>
     88
     89using namespace lemon;
     90
     91typedef dim2::Point<double> Point;
     92
     93GRAPH_TYPEDEFS(ListGraph);
     94
     95bool progress=true;
     96
     97int N;
     98// int girth;
     99
     100ListGraph g;
     101
     102std::vector<Node> nodes;
     103ListGraph::NodeMap<Point> coords(g);
     104
     105
     106double totalLen(){
     107  double tlen=0;
     108  for(EdgeIt e(g);e!=INVALID;++e)
     109    tlen+=sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
     110  return tlen;
     111}
     112
     113int tsp_impr_num=0;
     114
     115const double EPSILON=1e-8;
     116bool tsp_improve(Node u, Node v)
     117{
     118  double luv=std::sqrt((coords[v]-coords[u]).normSquare());
     119  Node u2=u;
     120  Node v2=v;
     121  do {
     122    Node n;
     123    for(IncEdgeIt e(g,v2);(n=g.runningNode(e))==u2;++e) { }
     124    u2=v2;
     125    v2=n;
     126    if(luv+std::sqrt((coords[v2]-coords[u2]).normSquare())-EPSILON>
     127       std::sqrt((coords[u]-coords[u2]).normSquare())+
     128       std::sqrt((coords[v]-coords[v2]).normSquare()))
     129      {
     130         g.erase(findEdge(g,u,v));
     131         g.erase(findEdge(g,u2,v2));
     132        g.addEdge(u2,u);
     133        g.addEdge(v,v2);
     134        tsp_impr_num++;
     135        return true;
     136      }
     137  } while(v2!=u);
     138  return false;
     139}
     140
     141bool tsp_improve(Node u)
     142{
     143  for(IncEdgeIt e(g,u);e!=INVALID;++e)
     144    if(tsp_improve(u,g.runningNode(e))) return true;
     145  return false;
     146}
     147
     148void tsp_improve()
     149{
     150  bool b;
     151  do {
     152    b=false;
     153    for(NodeIt n(g);n!=INVALID;++n)
     154      if(tsp_improve(n)) b=true;
     155  } while(b);
     156}
     157
     158void tsp()
     159{
     160  for(int i=0;i<N;i++) g.addEdge(nodes[i],nodes[(i+1)%N]);
     161  tsp_improve();
     162}
     163
     164class Line
     165{
     166public:
     167  Point a;
     168  Point b;
     169  Line(Point _a,Point _b) :a(_a),b(_b) {}
     170  Line(Node _a,Node _b) : a(coords[_a]),b(coords[_b]) {}
     171  Line(const Arc &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
     172  Line(const Edge &e) : a(coords[g.u(e)]),b(coords[g.v(e)]) {}
     173};
     174
     175inline std::ostream& operator<<(std::ostream &os, const Line &l)
     176{
     177  os << l.a << "->" << l.b;
     178  return os;
     179}
     180
     181bool cross(Line a, Line b)
     182{
     183  Point ao=rot90(a.b-a.a);
     184  Point bo=rot90(b.b-b.a);
     185  return (ao*(b.a-a.a))*(ao*(b.b-a.a))<0 &&
     186    (bo*(a.a-b.a))*(bo*(a.b-b.a))<0;
     187}
     188
     189struct Parc
     190{
     191  Node a;
     192  Node b;
     193  double len;
     194};
     195
     196bool pedgeLess(Parc a,Parc b)
     197{
     198  return a.len<b.len;
     199}
     200
     201std::vector<Edge> arcs;
     202
     203namespace _delaunay_bits {
     204
     205  struct Part {
     206    int prev, curr, next;
     207
     208    Part(int p, int c, int n) : prev(p), curr(c), next(n) {}
     209  };
     210
     211  inline std::ostream& operator<<(std::ostream& os, const Part& part) {
     212    os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')';
     213    return os;
     214  }
     215
     216  inline double circle_point(const Point& p, const Point& q, const Point& r) {
     217    double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y);
     218    if (a == 0) return std::numeric_limits<double>::quiet_NaN();
     219
     220    double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) +
     221      (q.x * q.x + q.y * q.y) * (r.y - p.y) +
     222      (r.x * r.x + r.y * r.y) * (p.y - q.y);
     223
     224    double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) +
     225      (q.x * q.x + q.y * q.y) * (r.x - p.x) +
     226      (r.x * r.x + r.y * r.y) * (p.x - q.x);
     227
     228    double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) +
     229      (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) +
     230      (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y);
     231
     232    return d / (2 * a) + sqrt((d * d + e * e) / (4 * a * a) + f / a);
     233  }
     234
     235  inline bool circle_form(const Point& p, const Point& q, const Point& r) {
     236    return rot90(q - p) * (r - q) < 0.0;
     237  }
     238
     239  inline double intersection(const Point& p, const Point& q, double sx) {
     240    const double epsilon = 1e-8;
     241
     242    if (p.x == q.x) return (p.y + q.y) / 2.0;
     243
     244    if (sx < p.x + epsilon) return p.y;
     245    if (sx < q.x + epsilon) return q.y;
     246
     247    double a = q.x - p.x;
     248    double b = (q.x - sx) * p.y - (p.x - sx) * q.y;
     249    double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare();
     250    return (b - sqrt(d)) / a;
     251  }
     252
     253  struct YLess {
     254
     255
     256    YLess(const std::vector<Point>& points, double& sweep)
     257      : _points(points), _sweep(sweep) {}
     258
     259    bool operator()(const Part& l, const Part& r) const {
     260      const double epsilon = 1e-8;
     261
     262      //      std::cerr << l << " vs " << r << std::endl;
     263      double lbx = l.prev != -1 ?
     264        intersection(_points[l.prev], _points[l.curr], _sweep) :
     265        - std::numeric_limits<double>::infinity();
     266      double rbx = r.prev != -1 ?
     267        intersection(_points[r.prev], _points[r.curr], _sweep) :
     268        - std::numeric_limits<double>::infinity();
     269      double lex = l.next != -1 ?
     270        intersection(_points[l.curr], _points[l.next], _sweep) :
     271        std::numeric_limits<double>::infinity();
     272      double rex = r.next != -1 ?
     273        intersection(_points[r.curr], _points[r.next], _sweep) :
     274        std::numeric_limits<double>::infinity();
     275
     276      if (lbx > lex) std::swap(lbx, lex);
     277      if (rbx > rex) std::swap(rbx, rex);
     278
     279      if (lex < epsilon + rex && lbx + epsilon < rex) return true;
     280      if (rex < epsilon + lex && rbx + epsilon < lex) return false;
     281      return lex < rex;
     282    }
     283
     284    const std::vector<Point>& _points;
     285    double& _sweep;
     286  };
     287
     288  struct BeachIt;
     289
     290  typedef std::multimap<double, BeachIt> SpikeHeap;
     291
     292  typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach;
     293
     294  struct BeachIt {
     295    Beach::iterator it;
     296
     297    BeachIt(Beach::iterator iter) : it(iter) {}
     298  };
     299
     300}
     301
     302inline void delaunay() {
     303  Counter cnt("Number of arcs added: ");
     304
     305  using namespace _delaunay_bits;
     306
     307  typedef _delaunay_bits::Part Part;
     308  typedef std::vector<std::pair<double, int> > SiteHeap;
     309
     310
     311  std::vector<Point> points;
     312  std::vector<Node> nodes;
     313
     314  for (NodeIt it(g); it != INVALID; ++it) {
     315    nodes.push_back(it);
     316    points.push_back(coords[it]);
     317  }
     318
     319  SiteHeap siteheap(points.size());
     320
     321  double sweep;
     322
     323
     324  for (int i = 0; i < int(siteheap.size()); ++i) {
     325    siteheap[i] = std::make_pair(points[i].x, i);
     326  }
     327
     328  std::sort(siteheap.begin(), siteheap.end());
     329  sweep = siteheap.front().first;
     330
     331  YLess yless(points, sweep);
     332  Beach beach(yless);
     333
     334  SpikeHeap spikeheap;
     335
     336  std::set<std::pair<int, int> > arcs;
     337
     338  int siteindex = 0;
     339  {
     340    SiteHeap front;
     341
     342    while (siteindex < int(siteheap.size()) &&
     343           siteheap[0].first == siteheap[siteindex].first) {
     344      front.push_back(std::make_pair(points[siteheap[siteindex].second].y,
     345                                     siteheap[siteindex].second));
     346      ++siteindex;
     347    }
     348
     349    std::sort(front.begin(), front.end());
     350
     351    for (int i = 0; i < int(front.size()); ++i) {
     352      int prev = (i == 0 ? -1 : front[i - 1].second);
     353      int curr = front[i].second;
     354      int next = (i + 1 == int(front.size()) ? -1 : front[i + 1].second);
     355
     356      beach.insert(std::make_pair(Part(prev, curr, next),
     357                                  spikeheap.end()));
     358    }
     359  }
     360
     361  while (siteindex < int(points.size()) || !spikeheap.empty()) {
     362
     363    SpikeHeap::iterator spit = spikeheap.begin();
     364
     365    if (siteindex < int(points.size()) &&
     366        (spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) {
     367      int site = siteheap[siteindex].second;
     368      sweep = siteheap[siteindex].first;
     369
     370      Beach::iterator bit = beach.upper_bound(Part(site, site, site));
     371
     372      if (bit->second != spikeheap.end()) {
     373        spikeheap.erase(bit->second);
     374      }
     375
     376      int prev = bit->first.prev;
     377      int curr = bit->first.curr;
     378      int next = bit->first.next;
     379
     380      beach.erase(bit);
     381
     382      SpikeHeap::iterator pit = spikeheap.end();
     383      if (prev != -1 &&
     384          circle_form(points[prev], points[curr], points[site])) {
     385        double x = circle_point(points[prev], points[curr], points[site]);
     386        pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
     387        pit->second.it =
     388          beach.insert(std::make_pair(Part(prev, curr, site), pit));
     389      } else {
     390        beach.insert(std::make_pair(Part(prev, curr, site), pit));
     391      }
     392
     393      beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end()));
     394
     395      SpikeHeap::iterator nit = spikeheap.end();
     396      if (next != -1 &&
     397          circle_form(points[site], points[curr],points[next])) {
     398        double x = circle_point(points[site], points[curr], points[next]);
     399        nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
     400        nit->second.it =
     401          beach.insert(std::make_pair(Part(site, curr, next), nit));
     402      } else {
     403        beach.insert(std::make_pair(Part(site, curr, next), nit));
     404      }
     405
     406      ++siteindex;
     407    } else {
     408      sweep = spit->first;
     409
     410      Beach::iterator bit = spit->second.it;
     411
     412      int prev = bit->first.prev;
     413      int curr = bit->first.curr;
     414      int next = bit->first.next;
     415
     416      {
     417        std::pair<int, int> arc;
     418
     419        arc = prev < curr ?
     420          std::make_pair(prev, curr) : std::make_pair(curr, prev);
     421
     422        if (arcs.find(arc) == arcs.end()) {
     423          arcs.insert(arc);
     424          g.addEdge(nodes[prev], nodes[curr]);
     425          ++cnt;
     426        }
     427
     428        arc = curr < next ?
     429          std::make_pair(curr, next) : std::make_pair(next, curr);
     430
     431        if (arcs.find(arc) == arcs.end()) {
     432          arcs.insert(arc);
     433          g.addEdge(nodes[curr], nodes[next]);
     434          ++cnt;
     435        }
     436      }
     437
     438      Beach::iterator pbit = bit; --pbit;
     439      int ppv = pbit->first.prev;
     440      Beach::iterator nbit = bit; ++nbit;
     441      int nnt = nbit->first.next;
     442
     443      if (bit->second != spikeheap.end()) spikeheap.erase(bit->second);
     444      if (pbit->second != spikeheap.end()) spikeheap.erase(pbit->second);
     445      if (nbit->second != spikeheap.end()) spikeheap.erase(nbit->second);
     446
     447      beach.erase(nbit);
     448      beach.erase(bit);
     449      beach.erase(pbit);
     450
     451      SpikeHeap::iterator pit = spikeheap.end();
     452      if (ppv != -1 && ppv != next &&
     453          circle_form(points[ppv], points[prev], points[next])) {
     454        double x = circle_point(points[ppv], points[prev], points[next]);
     455        if (x < sweep) x = sweep;
     456        pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
     457        pit->second.it =
     458          beach.insert(std::make_pair(Part(ppv, prev, next), pit));
     459      } else {
     460        beach.insert(std::make_pair(Part(ppv, prev, next), pit));
     461      }
     462
     463      SpikeHeap::iterator nit = spikeheap.end();
     464      if (nnt != -1 && prev != nnt &&
     465          circle_form(points[prev], points[next], points[nnt])) {
     466        double x = circle_point(points[prev], points[next], points[nnt]);
     467        if (x < sweep) x = sweep;
     468        nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
     469        nit->second.it =
     470          beach.insert(std::make_pair(Part(prev, next, nnt), nit));
     471      } else {
     472        beach.insert(std::make_pair(Part(prev, next, nnt), nit));
     473      }
     474
     475    }
     476  }
     477
     478  for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {
     479    int curr = it->first.curr;
     480    int next = it->first.next;
     481
     482    if (next == -1) continue;
     483
     484    std::pair<int, int> arc;
     485
     486    arc = curr < next ?
     487      std::make_pair(curr, next) : std::make_pair(next, curr);
     488
     489    if (arcs.find(arc) == arcs.end()) {
     490      arcs.insert(arc);
     491      g.addEdge(nodes[curr], nodes[next]);
     492      ++cnt;
     493    }
     494  }
     495}
     496
     497void sparse(int d)
     498{
     499  Counter cnt("Number of arcs removed: ");
     500  Bfs<ListGraph> bfs(g);
     501  for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
     502      ei!=arcs.rend();++ei)
     503    {
     504      Node a=g.u(*ei);
     505      Node b=g.v(*ei);
     506      g.erase(*ei);
     507      bfs.run(a,b);
     508      if(bfs.predArc(b)==INVALID || bfs.dist(b)>d)
     509        g.addEdge(a,b);
     510      else cnt++;
     511    }
     512}
     513
     514void sparse2(int d)
     515{
     516  Counter cnt("Number of arcs removed: ");
     517  for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
     518      ei!=arcs.rend();++ei)
     519    {
     520      Node a=g.u(*ei);
     521      Node b=g.v(*ei);
     522      g.erase(*ei);
     523      ConstMap<Arc,int> cegy(1);
     524      Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy,a,b);
     525      int k=sur.run(2);
     526      if(k<2 || sur.totalLength()>d)
     527        g.addEdge(a,b);
     528      else cnt++;
     529//       else std::cout << "Remove arc " << g.id(a) << "-" << g.id(b) << '\n';
     530    }
     531}
     532
     533void sparseTriangle(int d)
     534{
     535  Counter cnt("Number of arcs added: ");
     536  std::vector<Parc> pedges;
     537  for(NodeIt n(g);n!=INVALID;++n)
     538    for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
     539      {
     540        Parc p;
     541        p.a=n;
     542        p.b=m;
     543        p.len=(coords[m]-coords[n]).normSquare();
     544        pedges.push_back(p);
     545      }
     546  std::sort(pedges.begin(),pedges.end(),pedgeLess);
     547  for(std::vector<Parc>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
     548    {
     549      Line li(pi->a,pi->b);
     550      EdgeIt e(g);
     551      for(;e!=INVALID && !cross(e,li);++e) ;
     552      Edge ne;
     553      if(e==INVALID) {
     554        ConstMap<Arc,int> cegy(1);
     555        Suurballe<ListGraph,ConstMap<Arc,int> >
     556          sur(g,cegy,pi->a,pi->b);
     557        int k=sur.run(2);
     558        if(k<2 || sur.totalLength()>d)
     559          {
     560            ne=g.addEdge(pi->a,pi->b);
     561            arcs.push_back(ne);
     562            cnt++;
     563          }
     564      }
     565    }
     566}
     567
     568template <typename Graph, typename CoordMap>
     569class LengthSquareMap {
     570public:
     571  typedef typename Graph::Edge Key;
     572  typedef typename CoordMap::Value::Value Value;
     573
     574  LengthSquareMap(const Graph& graph, const CoordMap& coords)
     575    : _graph(graph), _coords(coords) {}
     576
     577  Value operator[](const Key& key) const {
     578    return (_coords[_graph.v(key)] -
     579            _coords[_graph.u(key)]).normSquare();
     580  }
     581
     582private:
     583
     584  const Graph& _graph;
     585  const CoordMap& _coords;
     586};
     587
     588void minTree() {
     589  std::vector<Parc> pedges;
     590  Timer T;
     591  std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
     592  delaunay();
     593  std::cout << T.realTime() << "s: Calculating spanning tree...\n";
     594  LengthSquareMap<ListGraph, ListGraph::NodeMap<Point> > ls(g, coords);
     595  ListGraph::EdgeMap<bool> tree(g);
     596  kruskal(g, ls, tree);
     597  std::cout << T.realTime() << "s: Removing non tree arcs...\n";
     598  std::vector<Edge> remove;
     599  for (EdgeIt e(g); e != INVALID; ++e) {
     600    if (!tree[e]) remove.push_back(e);
     601  }
     602  for(int i = 0; i < int(remove.size()); ++i) {
     603    g.erase(remove[i]);
     604  }
     605  std::cout << T.realTime() << "s: Done\n";
     606}
     607
     608void tsp2()
     609{
     610  std::cout << "Find a tree..." << std::endl;
     611
     612  minTree();
     613
     614  std::cout << "Total arc length (tree) : " << totalLen() << std::endl;
     615
     616  std::cout << "Make it Euler..." << std::endl;
     617
     618  {
     619    std::vector<Node> leafs;
     620    for(NodeIt n(g);n!=INVALID;++n)
     621      if(countIncEdges(g,n)%2==1) leafs.push_back(n);
     622
     623//    for(unsigned int i=0;i<leafs.size();i+=2)
     624//       g.addArc(leafs[i],leafs[i+1]);
     625
     626    std::vector<Parc> pedges;
     627    for(unsigned int i=0;i<leafs.size()-1;i++)
     628      for(unsigned int j=i+1;j<leafs.size();j++)
     629        {
     630          Node n=leafs[i];
     631          Node m=leafs[j];
     632          Parc p;
     633          p.a=n;
     634          p.b=m;
     635          p.len=(coords[m]-coords[n]).normSquare();
     636          pedges.push_back(p);
     637        }
     638    std::sort(pedges.begin(),pedges.end(),pedgeLess);
     639    for(unsigned int i=0;i<pedges.size();i++)
     640      if(countIncEdges(g,pedges[i].a)%2 &&
     641         countIncEdges(g,pedges[i].b)%2)
     642        g.addEdge(pedges[i].a,pedges[i].b);
     643  }
     644
     645  for(NodeIt n(g);n!=INVALID;++n)
     646    if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 )
     647      std::cout << "GEBASZ!!!" << std::endl;
     648
     649  for(EdgeIt e(g);e!=INVALID;++e)
     650    if(g.u(e)==g.v(e))
     651      std::cout << "LOOP GEBASZ!!!" << std::endl;
     652
     653  std::cout << "Number of arcs : " << countEdges(g) << std::endl;
     654
     655  std::cout << "Total arc length (euler) : " << totalLen() << std::endl;
     656
     657  ListGraph::EdgeMap<Arc> enext(g);
     658  {
     659    EulerIt<ListGraph> e(g);
     660    Arc eo=e;
     661    Arc ef=e;
     662//     std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
     663    for(++e;e!=INVALID;++e)
     664      {
     665//         std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
     666        enext[eo]=e;
     667        eo=e;
     668      }
     669    enext[eo]=ef;
     670  }
     671
     672  std::cout << "Creating a tour from that..." << std::endl;
     673
     674  int nnum = countNodes(g);
     675  int ednum = countEdges(g);
     676
     677  for(Arc p=enext[EdgeIt(g)];ednum>nnum;p=enext[p])
     678    {
     679//       std::cout << "Checking arc " << g.id(p) << std::endl;
     680      Arc e=enext[p];
     681      Arc f=enext[e];
     682      Node n2=g.source(f);
     683      Node n1=g.oppositeNode(n2,e);
     684      Node n3=g.oppositeNode(n2,f);
     685      if(countIncEdges(g,n2)>2)
     686        {
     687//           std::cout << "Remove an Arc" << std::endl;
     688          Arc ff=enext[f];
     689          g.erase(e);
     690          g.erase(f);
     691          if(n1!=n3)
     692            {
     693              Arc ne=g.direct(g.addEdge(n1,n3),n1);
     694              enext[p]=ne;
     695              enext[ne]=ff;
     696              ednum--;
     697            }
     698          else {
     699            enext[p]=ff;
     700            ednum-=2;
     701          }
     702        }
     703    }
     704
     705  std::cout << "Total arc length (tour) : " << totalLen() << std::endl;
     706
     707  std::cout << "2-opt the tour..." << std::endl;
     708
     709  tsp_improve();
     710
     711  std::cout << "Total arc length (2-opt tour) : " << totalLen() << std::endl;
     712}
     713
     714
     715int main(int argc,const char **argv)
     716{
     717  ArgParser ap(argc,argv);
     718
     719//   bool eps;
     720  bool disc_d, square_d, gauss_d;
     721//   bool tsp_a,two_a,tree_a;
     722  int num_of_cities=1;
     723  double area=1;
     724  N=100;
     725//   girth=10;
     726  std::string ndist("disc");
     727  ap.refOption("n", "Number of nodes (default is 100)", N)
     728    .intOption("g", "Girth parameter (default is 10)", 10)
     729    .refOption("cities", "Number of cities (default is 1)", num_of_cities)
     730    .refOption("area", "Full relative area of the cities (default is 1)", area)
     731    .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",disc_d)
     732    .optionGroup("dist", "disc")
     733    .refOption("square", "Nodes are evenly distributed on a unit square", square_d)
     734    .optionGroup("dist", "square")
     735    .refOption("gauss",
     736            "Nodes are located according to a two-dim gauss distribution",
     737            gauss_d)
     738    .optionGroup("dist", "gauss")
     739//     .mandatoryGroup("dist")
     740    .onlyOneGroup("dist")
     741    .boolOption("eps", "Also generate .eps output (prefix.eps)")
     742    .boolOption("dir", "Directed digraph is generated (each arcs are replaced by two directed ones)")
     743    .boolOption("2con", "Create a two connected planar digraph")
     744    .optionGroup("alg","2con")
     745    .boolOption("tree", "Create a min. cost spanning tree")
     746    .optionGroup("alg","tree")
     747    .boolOption("tsp", "Create a TSP tour")
     748    .optionGroup("alg","tsp")
     749    .boolOption("tsp2", "Create a TSP tour (tree based)")
     750    .optionGroup("alg","tsp2")
     751    .boolOption("dela", "Delaunay triangulation digraph")
     752    .optionGroup("alg","dela")
     753    .onlyOneGroup("alg")
     754    .boolOption("rand", "Use time seed for random number generator")
     755    .optionGroup("rand", "rand")
     756    .intOption("seed", "Random seed", -1)
     757    .optionGroup("rand", "seed")
     758    .onlyOneGroup("rand")
     759    .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
     760    .run();
     761
     762  if (ap["rand"]) {
     763    int seed = time(0);
     764    std::cout << "Random number seed: " << seed << std::endl;
     765    rnd = Random(seed);
     766  }
     767  if (ap.given("seed")) {
     768    int seed = ap["seed"];
     769    std::cout << "Random number seed: " << seed << std::endl;
     770    rnd = Random(seed);
     771  }
     772
     773  std::string prefix;
     774  switch(ap.files().size())
     775    {
     776    case 0:
     777      prefix="lgf-gen-out";
     778      break;
     779    case 1:
     780      prefix=ap.files()[0];
     781      break;
     782    default:
     783      std::cerr << "\nAt most one prefix can be given\n\n";
     784      exit(1);
     785    }
     786
     787  double sum_sizes=0;
     788  std::vector<double> sizes;
     789  std::vector<double> cum_sizes;
     790  for(int s=0;s<num_of_cities;s++)
     791    {
     792      //         sum_sizes+=rnd.exponential();
     793      double d=rnd();
     794      sum_sizes+=d;
     795      sizes.push_back(d);
     796      cum_sizes.push_back(sum_sizes);
     797    }
     798  int i=0;
     799  for(int s=0;s<num_of_cities;s++)
     800    {
     801      Point center=(num_of_cities==1?Point(0,0):rnd.disc());
     802      if(gauss_d)
     803        for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
     804          Node n=g.addNode();
     805          nodes.push_back(n);
     806          coords[n]=center+rnd.gauss2()*area*
     807            std::sqrt(sizes[s]/sum_sizes);
     808        }
     809      else if(square_d)
     810        for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
     811          Node n=g.addNode();
     812          nodes.push_back(n);
     813          coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*
     814            std::sqrt(sizes[s]/sum_sizes);
     815        }
     816      else if(disc_d || true)
     817        for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
     818          Node n=g.addNode();
     819          nodes.push_back(n);
     820          coords[n]=center+rnd.disc()*area*
     821            std::sqrt(sizes[s]/sum_sizes);
     822        }
     823    }
     824
     825//   for (ListGraph::NodeIt n(g); n != INVALID; ++n) {
     826//     std::cerr << coords[n] << std::endl;
     827//   }
     828
     829  if(ap["tsp"]) {
     830    tsp();
     831    std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
     832  }
     833  if(ap["tsp2"]) {
     834    tsp2();
     835    std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
     836  }
     837  else if(ap["2con"]) {
     838    std::cout << "Make triangles\n";
     839    //   triangle();
     840    sparseTriangle(ap["g"]);
     841    std::cout << "Make it sparser\n";
     842    sparse2(ap["g"]);
     843  }
     844  else if(ap["tree"]) {
     845    minTree();
     846  }
     847  else if(ap["dela"]) {
     848    delaunay();
     849  }
     850
     851
     852  std::cout << "Number of nodes    : " << countNodes(g) << std::endl;
     853  std::cout << "Number of arcs    : " << countEdges(g) << std::endl;
     854  double tlen=0;
     855  for(EdgeIt e(g);e!=INVALID;++e)
     856    tlen+=sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
     857  std::cout << "Total arc length  : " << tlen << std::endl;
     858
     859  if(ap["eps"])
     860    graphToEps(g,prefix+".eps").scaleToA4().
     861      scale(600).nodeScale(.005).arcWidthScale(.001).preScale(false).
     862      coords(coords).run();
     863
     864  if(ap["dir"])
     865    DigraphWriter<ListGraph>(g,prefix+".lgf").
     866      nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
     867      nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
     868      run();
     869  else GraphWriter<ListGraph>(g,prefix+".lgf").
     870         nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
     871         nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
     872         run();
     873}
     874