tools/lgf-gen.cc
author kpeter
Fri, 29 Feb 2008 15:55:13 +0000
changeset 2586 37fb2c384c78
parent 2553 bfced05fa852
child 2618 6aa6fcaeaea5
permissions -rw-r--r--
Reimplemented Suurballe class.

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