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