1.1 --- a/tools/Makefile.am Mon Nov 03 11:59:54 2008 +0000
1.2 +++ b/tools/Makefile.am Mon Feb 23 11:30:15 2009 +0000
1.3 @@ -1,10 +1,12 @@
1.4 if WANT_TOOLS
1.5
1.6 bin_PROGRAMS += \
1.7 - tools/dimacs-to-lgf
1.8 + tools/dimacs-to-lgf \
1.9 + tools/lgf-gen
1.10
1.11 dist_bin_SCRIPTS += tools/lemon-0.x-to-1.x.sh
1.12
1.13 endif WANT_TOOLS
1.14
1.15 tools_dimacs_to_lgf_SOURCES = tools/dimacs-to-lgf.cc
1.16 +tools_lgf_gen_SOURCES = tools/lgf-gen.cc
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/tools/lgf-gen.cc Mon Feb 23 11:30:15 2009 +0000
2.3 @@ -0,0 +1,834 @@
2.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
2.5 + *
2.6 + * This file is a part of LEMON, a generic C++ optimization library.
2.7 + *
2.8 + * Copyright (C) 2003-2009
2.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
2.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
2.11 + *
2.12 + * Permission to use, modify and distribute this software is granted
2.13 + * provided that this copyright notice appears in all copies. For
2.14 + * precise terms see the accompanying LICENSE file.
2.15 + *
2.16 + * This software is provided "AS IS" with no warranty of any kind,
2.17 + * express or implied, and with no claim as to its suitability for any
2.18 + * purpose.
2.19 + *
2.20 + */
2.21 +
2.22 +/// \ingroup tools
2.23 +/// \file
2.24 +/// \brief Special plane digraph generator.
2.25 +///
2.26 +/// Graph generator application for various types of plane graphs.
2.27 +///
2.28 +/// See
2.29 +/// \verbatim
2.30 +/// lgf-gen --help
2.31 +/// \endverbatim
2.32 +/// for more info on the usage.
2.33 +///
2.34 +
2.35 +
2.36 +#include <algorithm>
2.37 +#include <set>
2.38 +#include <lemon/list_graph.h>
2.39 +#include <lemon/random.h>
2.40 +#include <lemon/dim2.h>
2.41 +#include <lemon/bfs.h>
2.42 +#include <lemon/counter.h>
2.43 +#include <lemon/suurballe.h>
2.44 +#include <lemon/graph_to_eps.h>
2.45 +#include <lemon/lgf_writer.h>
2.46 +#include <lemon/arg_parser.h>
2.47 +#include <lemon/euler.h>
2.48 +#include <lemon/math.h>
2.49 +#include <lemon/kruskal.h>
2.50 +#include <lemon/time_measure.h>
2.51 +
2.52 +using namespace lemon;
2.53 +
2.54 +typedef dim2::Point<double> Point;
2.55 +
2.56 +GRAPH_TYPEDEFS(ListGraph);
2.57 +
2.58 +bool progress=true;
2.59 +
2.60 +int N;
2.61 +// int girth;
2.62 +
2.63 +ListGraph g;
2.64 +
2.65 +std::vector<Node> nodes;
2.66 +ListGraph::NodeMap<Point> coords(g);
2.67 +
2.68 +
2.69 +double totalLen(){
2.70 + double tlen=0;
2.71 + for(EdgeIt e(g);e!=INVALID;++e)
2.72 + tlen+=sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
2.73 + return tlen;
2.74 +}
2.75 +
2.76 +int tsp_impr_num=0;
2.77 +
2.78 +const double EPSILON=1e-8;
2.79 +bool tsp_improve(Node u, Node v)
2.80 +{
2.81 + double luv=std::sqrt((coords[v]-coords[u]).normSquare());
2.82 + Node u2=u;
2.83 + Node v2=v;
2.84 + do {
2.85 + Node n;
2.86 + for(IncEdgeIt e(g,v2);(n=g.runningNode(e))==u2;++e) { }
2.87 + u2=v2;
2.88 + v2=n;
2.89 + if(luv+std::sqrt((coords[v2]-coords[u2]).normSquare())-EPSILON>
2.90 + std::sqrt((coords[u]-coords[u2]).normSquare())+
2.91 + std::sqrt((coords[v]-coords[v2]).normSquare()))
2.92 + {
2.93 + g.erase(findEdge(g,u,v));
2.94 + g.erase(findEdge(g,u2,v2));
2.95 + g.addEdge(u2,u);
2.96 + g.addEdge(v,v2);
2.97 + tsp_impr_num++;
2.98 + return true;
2.99 + }
2.100 + } while(v2!=u);
2.101 + return false;
2.102 +}
2.103 +
2.104 +bool tsp_improve(Node u)
2.105 +{
2.106 + for(IncEdgeIt e(g,u);e!=INVALID;++e)
2.107 + if(tsp_improve(u,g.runningNode(e))) return true;
2.108 + return false;
2.109 +}
2.110 +
2.111 +void tsp_improve()
2.112 +{
2.113 + bool b;
2.114 + do {
2.115 + b=false;
2.116 + for(NodeIt n(g);n!=INVALID;++n)
2.117 + if(tsp_improve(n)) b=true;
2.118 + } while(b);
2.119 +}
2.120 +
2.121 +void tsp()
2.122 +{
2.123 + for(int i=0;i<N;i++) g.addEdge(nodes[i],nodes[(i+1)%N]);
2.124 + tsp_improve();
2.125 +}
2.126 +
2.127 +class Line
2.128 +{
2.129 +public:
2.130 + Point a;
2.131 + Point b;
2.132 + Line(Point _a,Point _b) :a(_a),b(_b) {}
2.133 + Line(Node _a,Node _b) : a(coords[_a]),b(coords[_b]) {}
2.134 + Line(const Arc &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
2.135 + Line(const Edge &e) : a(coords[g.u(e)]),b(coords[g.v(e)]) {}
2.136 +};
2.137 +
2.138 +inline std::ostream& operator<<(std::ostream &os, const Line &l)
2.139 +{
2.140 + os << l.a << "->" << l.b;
2.141 + return os;
2.142 +}
2.143 +
2.144 +bool cross(Line a, Line b)
2.145 +{
2.146 + Point ao=rot90(a.b-a.a);
2.147 + Point bo=rot90(b.b-b.a);
2.148 + return (ao*(b.a-a.a))*(ao*(b.b-a.a))<0 &&
2.149 + (bo*(a.a-b.a))*(bo*(a.b-b.a))<0;
2.150 +}
2.151 +
2.152 +struct Parc
2.153 +{
2.154 + Node a;
2.155 + Node b;
2.156 + double len;
2.157 +};
2.158 +
2.159 +bool pedgeLess(Parc a,Parc b)
2.160 +{
2.161 + return a.len<b.len;
2.162 +}
2.163 +
2.164 +std::vector<Edge> arcs;
2.165 +
2.166 +namespace _delaunay_bits {
2.167 +
2.168 + struct Part {
2.169 + int prev, curr, next;
2.170 +
2.171 + Part(int p, int c, int n) : prev(p), curr(c), next(n) {}
2.172 + };
2.173 +
2.174 + inline std::ostream& operator<<(std::ostream& os, const Part& part) {
2.175 + os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')';
2.176 + return os;
2.177 + }
2.178 +
2.179 + inline double circle_point(const Point& p, const Point& q, const Point& r) {
2.180 + double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y);
2.181 + if (a == 0) return std::numeric_limits<double>::quiet_NaN();
2.182 +
2.183 + double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) +
2.184 + (q.x * q.x + q.y * q.y) * (r.y - p.y) +
2.185 + (r.x * r.x + r.y * r.y) * (p.y - q.y);
2.186 +
2.187 + double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) +
2.188 + (q.x * q.x + q.y * q.y) * (r.x - p.x) +
2.189 + (r.x * r.x + r.y * r.y) * (p.x - q.x);
2.190 +
2.191 + double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) +
2.192 + (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) +
2.193 + (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y);
2.194 +
2.195 + return d / (2 * a) + sqrt((d * d + e * e) / (4 * a * a) + f / a);
2.196 + }
2.197 +
2.198 + inline bool circle_form(const Point& p, const Point& q, const Point& r) {
2.199 + return rot90(q - p) * (r - q) < 0.0;
2.200 + }
2.201 +
2.202 + inline double intersection(const Point& p, const Point& q, double sx) {
2.203 + const double epsilon = 1e-8;
2.204 +
2.205 + if (p.x == q.x) return (p.y + q.y) / 2.0;
2.206 +
2.207 + if (sx < p.x + epsilon) return p.y;
2.208 + if (sx < q.x + epsilon) return q.y;
2.209 +
2.210 + double a = q.x - p.x;
2.211 + double b = (q.x - sx) * p.y - (p.x - sx) * q.y;
2.212 + double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare();
2.213 + return (b - sqrt(d)) / a;
2.214 + }
2.215 +
2.216 + struct YLess {
2.217 +
2.218 +
2.219 + YLess(const std::vector<Point>& points, double& sweep)
2.220 + : _points(points), _sweep(sweep) {}
2.221 +
2.222 + bool operator()(const Part& l, const Part& r) const {
2.223 + const double epsilon = 1e-8;
2.224 +
2.225 + // std::cerr << l << " vs " << r << std::endl;
2.226 + double lbx = l.prev != -1 ?
2.227 + intersection(_points[l.prev], _points[l.curr], _sweep) :
2.228 + - std::numeric_limits<double>::infinity();
2.229 + double rbx = r.prev != -1 ?
2.230 + intersection(_points[r.prev], _points[r.curr], _sweep) :
2.231 + - std::numeric_limits<double>::infinity();
2.232 + double lex = l.next != -1 ?
2.233 + intersection(_points[l.curr], _points[l.next], _sweep) :
2.234 + std::numeric_limits<double>::infinity();
2.235 + double rex = r.next != -1 ?
2.236 + intersection(_points[r.curr], _points[r.next], _sweep) :
2.237 + std::numeric_limits<double>::infinity();
2.238 +
2.239 + if (lbx > lex) std::swap(lbx, lex);
2.240 + if (rbx > rex) std::swap(rbx, rex);
2.241 +
2.242 + if (lex < epsilon + rex && lbx + epsilon < rex) return true;
2.243 + if (rex < epsilon + lex && rbx + epsilon < lex) return false;
2.244 + return lex < rex;
2.245 + }
2.246 +
2.247 + const std::vector<Point>& _points;
2.248 + double& _sweep;
2.249 + };
2.250 +
2.251 + struct BeachIt;
2.252 +
2.253 + typedef std::multimap<double, BeachIt> SpikeHeap;
2.254 +
2.255 + typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach;
2.256 +
2.257 + struct BeachIt {
2.258 + Beach::iterator it;
2.259 +
2.260 + BeachIt(Beach::iterator iter) : it(iter) {}
2.261 + };
2.262 +
2.263 +}
2.264 +
2.265 +inline void delaunay() {
2.266 + Counter cnt("Number of arcs added: ");
2.267 +
2.268 + using namespace _delaunay_bits;
2.269 +
2.270 + typedef _delaunay_bits::Part Part;
2.271 + typedef std::vector<std::pair<double, int> > SiteHeap;
2.272 +
2.273 +
2.274 + std::vector<Point> points;
2.275 + std::vector<Node> nodes;
2.276 +
2.277 + for (NodeIt it(g); it != INVALID; ++it) {
2.278 + nodes.push_back(it);
2.279 + points.push_back(coords[it]);
2.280 + }
2.281 +
2.282 + SiteHeap siteheap(points.size());
2.283 +
2.284 + double sweep;
2.285 +
2.286 +
2.287 + for (int i = 0; i < int(siteheap.size()); ++i) {
2.288 + siteheap[i] = std::make_pair(points[i].x, i);
2.289 + }
2.290 +
2.291 + std::sort(siteheap.begin(), siteheap.end());
2.292 + sweep = siteheap.front().first;
2.293 +
2.294 + YLess yless(points, sweep);
2.295 + Beach beach(yless);
2.296 +
2.297 + SpikeHeap spikeheap;
2.298 +
2.299 + std::set<std::pair<int, int> > arcs;
2.300 +
2.301 + int siteindex = 0;
2.302 + {
2.303 + SiteHeap front;
2.304 +
2.305 + while (siteindex < int(siteheap.size()) &&
2.306 + siteheap[0].first == siteheap[siteindex].first) {
2.307 + front.push_back(std::make_pair(points[siteheap[siteindex].second].y,
2.308 + siteheap[siteindex].second));
2.309 + ++siteindex;
2.310 + }
2.311 +
2.312 + std::sort(front.begin(), front.end());
2.313 +
2.314 + for (int i = 0; i < int(front.size()); ++i) {
2.315 + int prev = (i == 0 ? -1 : front[i - 1].second);
2.316 + int curr = front[i].second;
2.317 + int next = (i + 1 == int(front.size()) ? -1 : front[i + 1].second);
2.318 +
2.319 + beach.insert(std::make_pair(Part(prev, curr, next),
2.320 + spikeheap.end()));
2.321 + }
2.322 + }
2.323 +
2.324 + while (siteindex < int(points.size()) || !spikeheap.empty()) {
2.325 +
2.326 + SpikeHeap::iterator spit = spikeheap.begin();
2.327 +
2.328 + if (siteindex < int(points.size()) &&
2.329 + (spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) {
2.330 + int site = siteheap[siteindex].second;
2.331 + sweep = siteheap[siteindex].first;
2.332 +
2.333 + Beach::iterator bit = beach.upper_bound(Part(site, site, site));
2.334 +
2.335 + if (bit->second != spikeheap.end()) {
2.336 + spikeheap.erase(bit->second);
2.337 + }
2.338 +
2.339 + int prev = bit->first.prev;
2.340 + int curr = bit->first.curr;
2.341 + int next = bit->first.next;
2.342 +
2.343 + beach.erase(bit);
2.344 +
2.345 + SpikeHeap::iterator pit = spikeheap.end();
2.346 + if (prev != -1 &&
2.347 + circle_form(points[prev], points[curr], points[site])) {
2.348 + double x = circle_point(points[prev], points[curr], points[site]);
2.349 + pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
2.350 + pit->second.it =
2.351 + beach.insert(std::make_pair(Part(prev, curr, site), pit));
2.352 + } else {
2.353 + beach.insert(std::make_pair(Part(prev, curr, site), pit));
2.354 + }
2.355 +
2.356 + beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end()));
2.357 +
2.358 + SpikeHeap::iterator nit = spikeheap.end();
2.359 + if (next != -1 &&
2.360 + circle_form(points[site], points[curr],points[next])) {
2.361 + double x = circle_point(points[site], points[curr], points[next]);
2.362 + nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
2.363 + nit->second.it =
2.364 + beach.insert(std::make_pair(Part(site, curr, next), nit));
2.365 + } else {
2.366 + beach.insert(std::make_pair(Part(site, curr, next), nit));
2.367 + }
2.368 +
2.369 + ++siteindex;
2.370 + } else {
2.371 + sweep = spit->first;
2.372 +
2.373 + Beach::iterator bit = spit->second.it;
2.374 +
2.375 + int prev = bit->first.prev;
2.376 + int curr = bit->first.curr;
2.377 + int next = bit->first.next;
2.378 +
2.379 + {
2.380 + std::pair<int, int> arc;
2.381 +
2.382 + arc = prev < curr ?
2.383 + std::make_pair(prev, curr) : std::make_pair(curr, prev);
2.384 +
2.385 + if (arcs.find(arc) == arcs.end()) {
2.386 + arcs.insert(arc);
2.387 + g.addEdge(nodes[prev], nodes[curr]);
2.388 + ++cnt;
2.389 + }
2.390 +
2.391 + arc = curr < next ?
2.392 + std::make_pair(curr, next) : std::make_pair(next, curr);
2.393 +
2.394 + if (arcs.find(arc) == arcs.end()) {
2.395 + arcs.insert(arc);
2.396 + g.addEdge(nodes[curr], nodes[next]);
2.397 + ++cnt;
2.398 + }
2.399 + }
2.400 +
2.401 + Beach::iterator pbit = bit; --pbit;
2.402 + int ppv = pbit->first.prev;
2.403 + Beach::iterator nbit = bit; ++nbit;
2.404 + int nnt = nbit->first.next;
2.405 +
2.406 + if (bit->second != spikeheap.end()) spikeheap.erase(bit->second);
2.407 + if (pbit->second != spikeheap.end()) spikeheap.erase(pbit->second);
2.408 + if (nbit->second != spikeheap.end()) spikeheap.erase(nbit->second);
2.409 +
2.410 + beach.erase(nbit);
2.411 + beach.erase(bit);
2.412 + beach.erase(pbit);
2.413 +
2.414 + SpikeHeap::iterator pit = spikeheap.end();
2.415 + if (ppv != -1 && ppv != next &&
2.416 + circle_form(points[ppv], points[prev], points[next])) {
2.417 + double x = circle_point(points[ppv], points[prev], points[next]);
2.418 + if (x < sweep) x = sweep;
2.419 + pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
2.420 + pit->second.it =
2.421 + beach.insert(std::make_pair(Part(ppv, prev, next), pit));
2.422 + } else {
2.423 + beach.insert(std::make_pair(Part(ppv, prev, next), pit));
2.424 + }
2.425 +
2.426 + SpikeHeap::iterator nit = spikeheap.end();
2.427 + if (nnt != -1 && prev != nnt &&
2.428 + circle_form(points[prev], points[next], points[nnt])) {
2.429 + double x = circle_point(points[prev], points[next], points[nnt]);
2.430 + if (x < sweep) x = sweep;
2.431 + nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
2.432 + nit->second.it =
2.433 + beach.insert(std::make_pair(Part(prev, next, nnt), nit));
2.434 + } else {
2.435 + beach.insert(std::make_pair(Part(prev, next, nnt), nit));
2.436 + }
2.437 +
2.438 + }
2.439 + }
2.440 +
2.441 + for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {
2.442 + int curr = it->first.curr;
2.443 + int next = it->first.next;
2.444 +
2.445 + if (next == -1) continue;
2.446 +
2.447 + std::pair<int, int> arc;
2.448 +
2.449 + arc = curr < next ?
2.450 + std::make_pair(curr, next) : std::make_pair(next, curr);
2.451 +
2.452 + if (arcs.find(arc) == arcs.end()) {
2.453 + arcs.insert(arc);
2.454 + g.addEdge(nodes[curr], nodes[next]);
2.455 + ++cnt;
2.456 + }
2.457 + }
2.458 +}
2.459 +
2.460 +void sparse(int d)
2.461 +{
2.462 + Counter cnt("Number of arcs removed: ");
2.463 + Bfs<ListGraph> bfs(g);
2.464 + for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
2.465 + ei!=arcs.rend();++ei)
2.466 + {
2.467 + Node a=g.u(*ei);
2.468 + Node b=g.v(*ei);
2.469 + g.erase(*ei);
2.470 + bfs.run(a,b);
2.471 + if(bfs.predArc(b)==INVALID || bfs.dist(b)>d)
2.472 + g.addEdge(a,b);
2.473 + else cnt++;
2.474 + }
2.475 +}
2.476 +
2.477 +void sparse2(int d)
2.478 +{
2.479 + Counter cnt("Number of arcs removed: ");
2.480 + for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
2.481 + ei!=arcs.rend();++ei)
2.482 + {
2.483 + Node a=g.u(*ei);
2.484 + Node b=g.v(*ei);
2.485 + g.erase(*ei);
2.486 + ConstMap<Arc,int> cegy(1);
2.487 + Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy,a,b);
2.488 + int k=sur.run(2);
2.489 + if(k<2 || sur.totalLength()>d)
2.490 + g.addEdge(a,b);
2.491 + else cnt++;
2.492 +// else std::cout << "Remove arc " << g.id(a) << "-" << g.id(b) << '\n';
2.493 + }
2.494 +}
2.495 +
2.496 +void sparseTriangle(int d)
2.497 +{
2.498 + Counter cnt("Number of arcs added: ");
2.499 + std::vector<Parc> pedges;
2.500 + for(NodeIt n(g);n!=INVALID;++n)
2.501 + for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
2.502 + {
2.503 + Parc p;
2.504 + p.a=n;
2.505 + p.b=m;
2.506 + p.len=(coords[m]-coords[n]).normSquare();
2.507 + pedges.push_back(p);
2.508 + }
2.509 + std::sort(pedges.begin(),pedges.end(),pedgeLess);
2.510 + for(std::vector<Parc>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
2.511 + {
2.512 + Line li(pi->a,pi->b);
2.513 + EdgeIt e(g);
2.514 + for(;e!=INVALID && !cross(e,li);++e) ;
2.515 + Edge ne;
2.516 + if(e==INVALID) {
2.517 + ConstMap<Arc,int> cegy(1);
2.518 + Suurballe<ListGraph,ConstMap<Arc,int> >
2.519 + sur(g,cegy,pi->a,pi->b);
2.520 + int k=sur.run(2);
2.521 + if(k<2 || sur.totalLength()>d)
2.522 + {
2.523 + ne=g.addEdge(pi->a,pi->b);
2.524 + arcs.push_back(ne);
2.525 + cnt++;
2.526 + }
2.527 + }
2.528 + }
2.529 +}
2.530 +
2.531 +template <typename Graph, typename CoordMap>
2.532 +class LengthSquareMap {
2.533 +public:
2.534 + typedef typename Graph::Edge Key;
2.535 + typedef typename CoordMap::Value::Value Value;
2.536 +
2.537 + LengthSquareMap(const Graph& graph, const CoordMap& coords)
2.538 + : _graph(graph), _coords(coords) {}
2.539 +
2.540 + Value operator[](const Key& key) const {
2.541 + return (_coords[_graph.v(key)] -
2.542 + _coords[_graph.u(key)]).normSquare();
2.543 + }
2.544 +
2.545 +private:
2.546 +
2.547 + const Graph& _graph;
2.548 + const CoordMap& _coords;
2.549 +};
2.550 +
2.551 +void minTree() {
2.552 + std::vector<Parc> pedges;
2.553 + Timer T;
2.554 + std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
2.555 + delaunay();
2.556 + std::cout << T.realTime() << "s: Calculating spanning tree...\n";
2.557 + LengthSquareMap<ListGraph, ListGraph::NodeMap<Point> > ls(g, coords);
2.558 + ListGraph::EdgeMap<bool> tree(g);
2.559 + kruskal(g, ls, tree);
2.560 + std::cout << T.realTime() << "s: Removing non tree arcs...\n";
2.561 + std::vector<Edge> remove;
2.562 + for (EdgeIt e(g); e != INVALID; ++e) {
2.563 + if (!tree[e]) remove.push_back(e);
2.564 + }
2.565 + for(int i = 0; i < int(remove.size()); ++i) {
2.566 + g.erase(remove[i]);
2.567 + }
2.568 + std::cout << T.realTime() << "s: Done\n";
2.569 +}
2.570 +
2.571 +void tsp2()
2.572 +{
2.573 + std::cout << "Find a tree..." << std::endl;
2.574 +
2.575 + minTree();
2.576 +
2.577 + std::cout << "Total arc length (tree) : " << totalLen() << std::endl;
2.578 +
2.579 + std::cout << "Make it Euler..." << std::endl;
2.580 +
2.581 + {
2.582 + std::vector<Node> leafs;
2.583 + for(NodeIt n(g);n!=INVALID;++n)
2.584 + if(countIncEdges(g,n)%2==1) leafs.push_back(n);
2.585 +
2.586 +// for(unsigned int i=0;i<leafs.size();i+=2)
2.587 +// g.addArc(leafs[i],leafs[i+1]);
2.588 +
2.589 + std::vector<Parc> pedges;
2.590 + for(unsigned int i=0;i<leafs.size()-1;i++)
2.591 + for(unsigned int j=i+1;j<leafs.size();j++)
2.592 + {
2.593 + Node n=leafs[i];
2.594 + Node m=leafs[j];
2.595 + Parc p;
2.596 + p.a=n;
2.597 + p.b=m;
2.598 + p.len=(coords[m]-coords[n]).normSquare();
2.599 + pedges.push_back(p);
2.600 + }
2.601 + std::sort(pedges.begin(),pedges.end(),pedgeLess);
2.602 + for(unsigned int i=0;i<pedges.size();i++)
2.603 + if(countIncEdges(g,pedges[i].a)%2 &&
2.604 + countIncEdges(g,pedges[i].b)%2)
2.605 + g.addEdge(pedges[i].a,pedges[i].b);
2.606 + }
2.607 +
2.608 + for(NodeIt n(g);n!=INVALID;++n)
2.609 + if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 )
2.610 + std::cout << "GEBASZ!!!" << std::endl;
2.611 +
2.612 + for(EdgeIt e(g);e!=INVALID;++e)
2.613 + if(g.u(e)==g.v(e))
2.614 + std::cout << "LOOP GEBASZ!!!" << std::endl;
2.615 +
2.616 + std::cout << "Number of arcs : " << countEdges(g) << std::endl;
2.617 +
2.618 + std::cout << "Total arc length (euler) : " << totalLen() << std::endl;
2.619 +
2.620 + ListGraph::EdgeMap<Arc> enext(g);
2.621 + {
2.622 + EulerIt<ListGraph> e(g);
2.623 + Arc eo=e;
2.624 + Arc ef=e;
2.625 +// std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
2.626 + for(++e;e!=INVALID;++e)
2.627 + {
2.628 +// std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
2.629 + enext[eo]=e;
2.630 + eo=e;
2.631 + }
2.632 + enext[eo]=ef;
2.633 + }
2.634 +
2.635 + std::cout << "Creating a tour from that..." << std::endl;
2.636 +
2.637 + int nnum = countNodes(g);
2.638 + int ednum = countEdges(g);
2.639 +
2.640 + for(Arc p=enext[EdgeIt(g)];ednum>nnum;p=enext[p])
2.641 + {
2.642 +// std::cout << "Checking arc " << g.id(p) << std::endl;
2.643 + Arc e=enext[p];
2.644 + Arc f=enext[e];
2.645 + Node n2=g.source(f);
2.646 + Node n1=g.oppositeNode(n2,e);
2.647 + Node n3=g.oppositeNode(n2,f);
2.648 + if(countIncEdges(g,n2)>2)
2.649 + {
2.650 +// std::cout << "Remove an Arc" << std::endl;
2.651 + Arc ff=enext[f];
2.652 + g.erase(e);
2.653 + g.erase(f);
2.654 + if(n1!=n3)
2.655 + {
2.656 + Arc ne=g.direct(g.addEdge(n1,n3),n1);
2.657 + enext[p]=ne;
2.658 + enext[ne]=ff;
2.659 + ednum--;
2.660 + }
2.661 + else {
2.662 + enext[p]=ff;
2.663 + ednum-=2;
2.664 + }
2.665 + }
2.666 + }
2.667 +
2.668 + std::cout << "Total arc length (tour) : " << totalLen() << std::endl;
2.669 +
2.670 + std::cout << "2-opt the tour..." << std::endl;
2.671 +
2.672 + tsp_improve();
2.673 +
2.674 + std::cout << "Total arc length (2-opt tour) : " << totalLen() << std::endl;
2.675 +}
2.676 +
2.677 +
2.678 +int main(int argc,const char **argv)
2.679 +{
2.680 + ArgParser ap(argc,argv);
2.681 +
2.682 +// bool eps;
2.683 + bool disc_d, square_d, gauss_d;
2.684 +// bool tsp_a,two_a,tree_a;
2.685 + int num_of_cities=1;
2.686 + double area=1;
2.687 + N=100;
2.688 +// girth=10;
2.689 + std::string ndist("disc");
2.690 + ap.refOption("n", "Number of nodes (default is 100)", N)
2.691 + .intOption("g", "Girth parameter (default is 10)", 10)
2.692 + .refOption("cities", "Number of cities (default is 1)", num_of_cities)
2.693 + .refOption("area", "Full relative area of the cities (default is 1)", area)
2.694 + .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",disc_d)
2.695 + .optionGroup("dist", "disc")
2.696 + .refOption("square", "Nodes are evenly distributed on a unit square", square_d)
2.697 + .optionGroup("dist", "square")
2.698 + .refOption("gauss",
2.699 + "Nodes are located according to a two-dim gauss distribution",
2.700 + gauss_d)
2.701 + .optionGroup("dist", "gauss")
2.702 +// .mandatoryGroup("dist")
2.703 + .onlyOneGroup("dist")
2.704 + .boolOption("eps", "Also generate .eps output (prefix.eps)")
2.705 + .boolOption("dir", "Directed digraph is generated (each arcs are replaced by two directed ones)")
2.706 + .boolOption("2con", "Create a two connected planar digraph")
2.707 + .optionGroup("alg","2con")
2.708 + .boolOption("tree", "Create a min. cost spanning tree")
2.709 + .optionGroup("alg","tree")
2.710 + .boolOption("tsp", "Create a TSP tour")
2.711 + .optionGroup("alg","tsp")
2.712 + .boolOption("tsp2", "Create a TSP tour (tree based)")
2.713 + .optionGroup("alg","tsp2")
2.714 + .boolOption("dela", "Delaunay triangulation digraph")
2.715 + .optionGroup("alg","dela")
2.716 + .onlyOneGroup("alg")
2.717 + .boolOption("rand", "Use time seed for random number generator")
2.718 + .optionGroup("rand", "rand")
2.719 + .intOption("seed", "Random seed", -1)
2.720 + .optionGroup("rand", "seed")
2.721 + .onlyOneGroup("rand")
2.722 + .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
2.723 + .run();
2.724 +
2.725 + if (ap["rand"]) {
2.726 + int seed = time(0);
2.727 + std::cout << "Random number seed: " << seed << std::endl;
2.728 + rnd = Random(seed);
2.729 + }
2.730 + if (ap.given("seed")) {
2.731 + int seed = ap["seed"];
2.732 + std::cout << "Random number seed: " << seed << std::endl;
2.733 + rnd = Random(seed);
2.734 + }
2.735 +
2.736 + std::string prefix;
2.737 + switch(ap.files().size())
2.738 + {
2.739 + case 0:
2.740 + prefix="lgf-gen-out";
2.741 + break;
2.742 + case 1:
2.743 + prefix=ap.files()[0];
2.744 + break;
2.745 + default:
2.746 + std::cerr << "\nAt most one prefix can be given\n\n";
2.747 + exit(1);
2.748 + }
2.749 +
2.750 + double sum_sizes=0;
2.751 + std::vector<double> sizes;
2.752 + std::vector<double> cum_sizes;
2.753 + for(int s=0;s<num_of_cities;s++)
2.754 + {
2.755 + // sum_sizes+=rnd.exponential();
2.756 + double d=rnd();
2.757 + sum_sizes+=d;
2.758 + sizes.push_back(d);
2.759 + cum_sizes.push_back(sum_sizes);
2.760 + }
2.761 + int i=0;
2.762 + for(int s=0;s<num_of_cities;s++)
2.763 + {
2.764 + Point center=(num_of_cities==1?Point(0,0):rnd.disc());
2.765 + if(gauss_d)
2.766 + for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
2.767 + Node n=g.addNode();
2.768 + nodes.push_back(n);
2.769 + coords[n]=center+rnd.gauss2()*area*
2.770 + std::sqrt(sizes[s]/sum_sizes);
2.771 + }
2.772 + else if(square_d)
2.773 + for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
2.774 + Node n=g.addNode();
2.775 + nodes.push_back(n);
2.776 + coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*
2.777 + std::sqrt(sizes[s]/sum_sizes);
2.778 + }
2.779 + else if(disc_d || true)
2.780 + for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
2.781 + Node n=g.addNode();
2.782 + nodes.push_back(n);
2.783 + coords[n]=center+rnd.disc()*area*
2.784 + std::sqrt(sizes[s]/sum_sizes);
2.785 + }
2.786 + }
2.787 +
2.788 +// for (ListGraph::NodeIt n(g); n != INVALID; ++n) {
2.789 +// std::cerr << coords[n] << std::endl;
2.790 +// }
2.791 +
2.792 + if(ap["tsp"]) {
2.793 + tsp();
2.794 + std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
2.795 + }
2.796 + if(ap["tsp2"]) {
2.797 + tsp2();
2.798 + std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
2.799 + }
2.800 + else if(ap["2con"]) {
2.801 + std::cout << "Make triangles\n";
2.802 + // triangle();
2.803 + sparseTriangle(ap["g"]);
2.804 + std::cout << "Make it sparser\n";
2.805 + sparse2(ap["g"]);
2.806 + }
2.807 + else if(ap["tree"]) {
2.808 + minTree();
2.809 + }
2.810 + else if(ap["dela"]) {
2.811 + delaunay();
2.812 + }
2.813 +
2.814 +
2.815 + std::cout << "Number of nodes : " << countNodes(g) << std::endl;
2.816 + std::cout << "Number of arcs : " << countEdges(g) << std::endl;
2.817 + double tlen=0;
2.818 + for(EdgeIt e(g);e!=INVALID;++e)
2.819 + tlen+=sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
2.820 + std::cout << "Total arc length : " << tlen << std::endl;
2.821 +
2.822 + if(ap["eps"])
2.823 + graphToEps(g,prefix+".eps").scaleToA4().
2.824 + scale(600).nodeScale(.005).arcWidthScale(.001).preScale(false).
2.825 + coords(coords).run();
2.826 +
2.827 + if(ap["dir"])
2.828 + DigraphWriter<ListGraph>(g,prefix+".lgf").
2.829 + nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
2.830 + nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
2.831 + run();
2.832 + else GraphWriter<ListGraph>(g,prefix+".lgf").
2.833 + nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
2.834 + nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
2.835 + run();
2.836 +}
2.837 +