tools/lgf-gen.cc
author Peter Kovacs <kpeter@inf.elte.hu>
Tue, 24 Mar 2009 00:18:25 +0100
changeset 604 8c3112a66878
parent 523 d9e43511d11c
child 570 ab6da8cf5ab2
permissions -rw-r--r--
Use XTI implementation instead of ATI in NetworkSimplex (#234)

XTI (eXtended Threaded Index) is an imporved version of the widely
known ATI (Augmented Threaded Index) method for storing and updating
the spanning tree structure in Network Simplex algorithms.

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