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