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