1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/tools/lgf-gen.cc Thu Dec 10 17:05:35 2009 +0100
1.3 @@ -0,0 +1,834 @@
1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library.
1.7 + *
1.8 + * Copyright (C) 2003-2009
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +/// \ingroup tools
1.23 +/// \file
1.24 +/// \brief Special plane graph generator.
1.25 +///
1.26 +/// Graph generator application for various types of plane graphs.
1.27 +///
1.28 +/// See
1.29 +/// \code
1.30 +/// lgf-gen --help
1.31 +/// \endcode
1.32 +/// for more information on the usage.
1.33 +
1.34 +#include <algorithm>
1.35 +#include <set>
1.36 +#include <ctime>
1.37 +#include <lemon/list_graph.h>
1.38 +#include <lemon/random.h>
1.39 +#include <lemon/dim2.h>
1.40 +#include <lemon/bfs.h>
1.41 +#include <lemon/counter.h>
1.42 +#include <lemon/suurballe.h>
1.43 +#include <lemon/graph_to_eps.h>
1.44 +#include <lemon/lgf_writer.h>
1.45 +#include <lemon/arg_parser.h>
1.46 +#include <lemon/euler.h>
1.47 +#include <lemon/math.h>
1.48 +#include <lemon/kruskal.h>
1.49 +#include <lemon/time_measure.h>
1.50 +
1.51 +using namespace lemon;
1.52 +
1.53 +typedef dim2::Point<double> Point;
1.54 +
1.55 +GRAPH_TYPEDEFS(ListGraph);
1.56 +
1.57 +bool progress=true;
1.58 +
1.59 +int N;
1.60 +// int girth;
1.61 +
1.62 +ListGraph g;
1.63 +
1.64 +std::vector<Node> nodes;
1.65 +ListGraph::NodeMap<Point> coords(g);
1.66 +
1.67 +
1.68 +double totalLen(){
1.69 + double tlen=0;
1.70 + for(EdgeIt e(g);e!=INVALID;++e)
1.71 + tlen+=std::sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
1.72 + return tlen;
1.73 +}
1.74 +
1.75 +int tsp_impr_num=0;
1.76 +
1.77 +const double EPSILON=1e-8;
1.78 +bool tsp_improve(Node u, Node v)
1.79 +{
1.80 + double luv=std::sqrt((coords[v]-coords[u]).normSquare());
1.81 + Node u2=u;
1.82 + Node v2=v;
1.83 + do {
1.84 + Node n;
1.85 + for(IncEdgeIt e(g,v2);(n=g.runningNode(e))==u2;++e) { }
1.86 + u2=v2;
1.87 + v2=n;
1.88 + if(luv+std::sqrt((coords[v2]-coords[u2]).normSquare())-EPSILON>
1.89 + std::sqrt((coords[u]-coords[u2]).normSquare())+
1.90 + std::sqrt((coords[v]-coords[v2]).normSquare()))
1.91 + {
1.92 + g.erase(findEdge(g,u,v));
1.93 + g.erase(findEdge(g,u2,v2));
1.94 + g.addEdge(u2,u);
1.95 + g.addEdge(v,v2);
1.96 + tsp_impr_num++;
1.97 + return true;
1.98 + }
1.99 + } while(v2!=u);
1.100 + return false;
1.101 +}
1.102 +
1.103 +bool tsp_improve(Node u)
1.104 +{
1.105 + for(IncEdgeIt e(g,u);e!=INVALID;++e)
1.106 + if(tsp_improve(u,g.runningNode(e))) return true;
1.107 + return false;
1.108 +}
1.109 +
1.110 +void tsp_improve()
1.111 +{
1.112 + bool b;
1.113 + do {
1.114 + b=false;
1.115 + for(NodeIt n(g);n!=INVALID;++n)
1.116 + if(tsp_improve(n)) b=true;
1.117 + } while(b);
1.118 +}
1.119 +
1.120 +void tsp()
1.121 +{
1.122 + for(int i=0;i<N;i++) g.addEdge(nodes[i],nodes[(i+1)%N]);
1.123 + tsp_improve();
1.124 +}
1.125 +
1.126 +class Line
1.127 +{
1.128 +public:
1.129 + Point a;
1.130 + Point b;
1.131 + Line(Point _a,Point _b) :a(_a),b(_b) {}
1.132 + Line(Node _a,Node _b) : a(coords[_a]),b(coords[_b]) {}
1.133 + Line(const Arc &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
1.134 + Line(const Edge &e) : a(coords[g.u(e)]),b(coords[g.v(e)]) {}
1.135 +};
1.136 +
1.137 +inline std::ostream& operator<<(std::ostream &os, const Line &l)
1.138 +{
1.139 + os << l.a << "->" << l.b;
1.140 + return os;
1.141 +}
1.142 +
1.143 +bool cross(Line a, Line b)
1.144 +{
1.145 + Point ao=rot90(a.b-a.a);
1.146 + Point bo=rot90(b.b-b.a);
1.147 + return (ao*(b.a-a.a))*(ao*(b.b-a.a))<0 &&
1.148 + (bo*(a.a-b.a))*(bo*(a.b-b.a))<0;
1.149 +}
1.150 +
1.151 +struct Parc
1.152 +{
1.153 + Node a;
1.154 + Node b;
1.155 + double len;
1.156 +};
1.157 +
1.158 +bool pedgeLess(Parc a,Parc b)
1.159 +{
1.160 + return a.len<b.len;
1.161 +}
1.162 +
1.163 +std::vector<Edge> arcs;
1.164 +
1.165 +namespace _delaunay_bits {
1.166 +
1.167 + struct Part {
1.168 + int prev, curr, next;
1.169 +
1.170 + Part(int p, int c, int n) : prev(p), curr(c), next(n) {}
1.171 + };
1.172 +
1.173 + inline std::ostream& operator<<(std::ostream& os, const Part& part) {
1.174 + os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')';
1.175 + return os;
1.176 + }
1.177 +
1.178 + inline double circle_point(const Point& p, const Point& q, const Point& r) {
1.179 + double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y);
1.180 + if (a == 0) return std::numeric_limits<double>::quiet_NaN();
1.181 +
1.182 + double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) +
1.183 + (q.x * q.x + q.y * q.y) * (r.y - p.y) +
1.184 + (r.x * r.x + r.y * r.y) * (p.y - q.y);
1.185 +
1.186 + double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) +
1.187 + (q.x * q.x + q.y * q.y) * (r.x - p.x) +
1.188 + (r.x * r.x + r.y * r.y) * (p.x - q.x);
1.189 +
1.190 + double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) +
1.191 + (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) +
1.192 + (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y);
1.193 +
1.194 + return d / (2 * a) + std::sqrt((d * d + e * e) / (4 * a * a) + f / a);
1.195 + }
1.196 +
1.197 + inline bool circle_form(const Point& p, const Point& q, const Point& r) {
1.198 + return rot90(q - p) * (r - q) < 0.0;
1.199 + }
1.200 +
1.201 + inline double intersection(const Point& p, const Point& q, double sx) {
1.202 + const double epsilon = 1e-8;
1.203 +
1.204 + if (p.x == q.x) return (p.y + q.y) / 2.0;
1.205 +
1.206 + if (sx < p.x + epsilon) return p.y;
1.207 + if (sx < q.x + epsilon) return q.y;
1.208 +
1.209 + double a = q.x - p.x;
1.210 + double b = (q.x - sx) * p.y - (p.x - sx) * q.y;
1.211 + double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare();
1.212 + return (b - std::sqrt(d)) / a;
1.213 + }
1.214 +
1.215 + struct YLess {
1.216 +
1.217 +
1.218 + YLess(const std::vector<Point>& points, double& sweep)
1.219 + : _points(points), _sweep(sweep) {}
1.220 +
1.221 + bool operator()(const Part& l, const Part& r) const {
1.222 + const double epsilon = 1e-8;
1.223 +
1.224 + // std::cerr << l << " vs " << r << std::endl;
1.225 + double lbx = l.prev != -1 ?
1.226 + intersection(_points[l.prev], _points[l.curr], _sweep) :
1.227 + - std::numeric_limits<double>::infinity();
1.228 + double rbx = r.prev != -1 ?
1.229 + intersection(_points[r.prev], _points[r.curr], _sweep) :
1.230 + - std::numeric_limits<double>::infinity();
1.231 + double lex = l.next != -1 ?
1.232 + intersection(_points[l.curr], _points[l.next], _sweep) :
1.233 + std::numeric_limits<double>::infinity();
1.234 + double rex = r.next != -1 ?
1.235 + intersection(_points[r.curr], _points[r.next], _sweep) :
1.236 + std::numeric_limits<double>::infinity();
1.237 +
1.238 + if (lbx > lex) std::swap(lbx, lex);
1.239 + if (rbx > rex) std::swap(rbx, rex);
1.240 +
1.241 + if (lex < epsilon + rex && lbx + epsilon < rex) return true;
1.242 + if (rex < epsilon + lex && rbx + epsilon < lex) return false;
1.243 + return lex < rex;
1.244 + }
1.245 +
1.246 + const std::vector<Point>& _points;
1.247 + double& _sweep;
1.248 + };
1.249 +
1.250 + struct BeachIt;
1.251 +
1.252 + typedef std::multimap<double, BeachIt> SpikeHeap;
1.253 +
1.254 + typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach;
1.255 +
1.256 + struct BeachIt {
1.257 + Beach::iterator it;
1.258 +
1.259 + BeachIt(Beach::iterator iter) : it(iter) {}
1.260 + };
1.261 +
1.262 +}
1.263 +
1.264 +inline void delaunay() {
1.265 + Counter cnt("Number of arcs added: ");
1.266 +
1.267 + using namespace _delaunay_bits;
1.268 +
1.269 + typedef _delaunay_bits::Part Part;
1.270 + typedef std::vector<std::pair<double, int> > SiteHeap;
1.271 +
1.272 +
1.273 + std::vector<Point> points;
1.274 + std::vector<Node> nodes;
1.275 +
1.276 + for (NodeIt it(g); it != INVALID; ++it) {
1.277 + nodes.push_back(it);
1.278 + points.push_back(coords[it]);
1.279 + }
1.280 +
1.281 + SiteHeap siteheap(points.size());
1.282 +
1.283 + double sweep;
1.284 +
1.285 +
1.286 + for (int i = 0; i < int(siteheap.size()); ++i) {
1.287 + siteheap[i] = std::make_pair(points[i].x, i);
1.288 + }
1.289 +
1.290 + std::sort(siteheap.begin(), siteheap.end());
1.291 + sweep = siteheap.front().first;
1.292 +
1.293 + YLess yless(points, sweep);
1.294 + Beach beach(yless);
1.295 +
1.296 + SpikeHeap spikeheap;
1.297 +
1.298 + std::set<std::pair<int, int> > arcs;
1.299 +
1.300 + int siteindex = 0;
1.301 + {
1.302 + SiteHeap front;
1.303 +
1.304 + while (siteindex < int(siteheap.size()) &&
1.305 + siteheap[0].first == siteheap[siteindex].first) {
1.306 + front.push_back(std::make_pair(points[siteheap[siteindex].second].y,
1.307 + siteheap[siteindex].second));
1.308 + ++siteindex;
1.309 + }
1.310 +
1.311 + std::sort(front.begin(), front.end());
1.312 +
1.313 + for (int i = 0; i < int(front.size()); ++i) {
1.314 + int prev = (i == 0 ? -1 : front[i - 1].second);
1.315 + int curr = front[i].second;
1.316 + int next = (i + 1 == int(front.size()) ? -1 : front[i + 1].second);
1.317 +
1.318 + beach.insert(std::make_pair(Part(prev, curr, next),
1.319 + spikeheap.end()));
1.320 + }
1.321 + }
1.322 +
1.323 + while (siteindex < int(points.size()) || !spikeheap.empty()) {
1.324 +
1.325 + SpikeHeap::iterator spit = spikeheap.begin();
1.326 +
1.327 + if (siteindex < int(points.size()) &&
1.328 + (spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) {
1.329 + int site = siteheap[siteindex].second;
1.330 + sweep = siteheap[siteindex].first;
1.331 +
1.332 + Beach::iterator bit = beach.upper_bound(Part(site, site, site));
1.333 +
1.334 + if (bit->second != spikeheap.end()) {
1.335 + spikeheap.erase(bit->second);
1.336 + }
1.337 +
1.338 + int prev = bit->first.prev;
1.339 + int curr = bit->first.curr;
1.340 + int next = bit->first.next;
1.341 +
1.342 + beach.erase(bit);
1.343 +
1.344 + SpikeHeap::iterator pit = spikeheap.end();
1.345 + if (prev != -1 &&
1.346 + circle_form(points[prev], points[curr], points[site])) {
1.347 + double x = circle_point(points[prev], points[curr], points[site]);
1.348 + pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
1.349 + pit->second.it =
1.350 + beach.insert(std::make_pair(Part(prev, curr, site), pit));
1.351 + } else {
1.352 + beach.insert(std::make_pair(Part(prev, curr, site), pit));
1.353 + }
1.354 +
1.355 + beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end()));
1.356 +
1.357 + SpikeHeap::iterator nit = spikeheap.end();
1.358 + if (next != -1 &&
1.359 + circle_form(points[site], points[curr],points[next])) {
1.360 + double x = circle_point(points[site], points[curr], points[next]);
1.361 + nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
1.362 + nit->second.it =
1.363 + beach.insert(std::make_pair(Part(site, curr, next), nit));
1.364 + } else {
1.365 + beach.insert(std::make_pair(Part(site, curr, next), nit));
1.366 + }
1.367 +
1.368 + ++siteindex;
1.369 + } else {
1.370 + sweep = spit->first;
1.371 +
1.372 + Beach::iterator bit = spit->second.it;
1.373 +
1.374 + int prev = bit->first.prev;
1.375 + int curr = bit->first.curr;
1.376 + int next = bit->first.next;
1.377 +
1.378 + {
1.379 + std::pair<int, int> arc;
1.380 +
1.381 + arc = prev < curr ?
1.382 + std::make_pair(prev, curr) : std::make_pair(curr, prev);
1.383 +
1.384 + if (arcs.find(arc) == arcs.end()) {
1.385 + arcs.insert(arc);
1.386 + g.addEdge(nodes[prev], nodes[curr]);
1.387 + ++cnt;
1.388 + }
1.389 +
1.390 + arc = curr < next ?
1.391 + std::make_pair(curr, next) : std::make_pair(next, curr);
1.392 +
1.393 + if (arcs.find(arc) == arcs.end()) {
1.394 + arcs.insert(arc);
1.395 + g.addEdge(nodes[curr], nodes[next]);
1.396 + ++cnt;
1.397 + }
1.398 + }
1.399 +
1.400 + Beach::iterator pbit = bit; --pbit;
1.401 + int ppv = pbit->first.prev;
1.402 + Beach::iterator nbit = bit; ++nbit;
1.403 + int nnt = nbit->first.next;
1.404 +
1.405 + if (bit->second != spikeheap.end()) spikeheap.erase(bit->second);
1.406 + if (pbit->second != spikeheap.end()) spikeheap.erase(pbit->second);
1.407 + if (nbit->second != spikeheap.end()) spikeheap.erase(nbit->second);
1.408 +
1.409 + beach.erase(nbit);
1.410 + beach.erase(bit);
1.411 + beach.erase(pbit);
1.412 +
1.413 + SpikeHeap::iterator pit = spikeheap.end();
1.414 + if (ppv != -1 && ppv != next &&
1.415 + circle_form(points[ppv], points[prev], points[next])) {
1.416 + double x = circle_point(points[ppv], points[prev], points[next]);
1.417 + if (x < sweep) x = sweep;
1.418 + pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
1.419 + pit->second.it =
1.420 + beach.insert(std::make_pair(Part(ppv, prev, next), pit));
1.421 + } else {
1.422 + beach.insert(std::make_pair(Part(ppv, prev, next), pit));
1.423 + }
1.424 +
1.425 + SpikeHeap::iterator nit = spikeheap.end();
1.426 + if (nnt != -1 && prev != nnt &&
1.427 + circle_form(points[prev], points[next], points[nnt])) {
1.428 + double x = circle_point(points[prev], points[next], points[nnt]);
1.429 + if (x < sweep) x = sweep;
1.430 + nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
1.431 + nit->second.it =
1.432 + beach.insert(std::make_pair(Part(prev, next, nnt), nit));
1.433 + } else {
1.434 + beach.insert(std::make_pair(Part(prev, next, nnt), nit));
1.435 + }
1.436 +
1.437 + }
1.438 + }
1.439 +
1.440 + for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {
1.441 + int curr = it->first.curr;
1.442 + int next = it->first.next;
1.443 +
1.444 + if (next == -1) continue;
1.445 +
1.446 + std::pair<int, int> arc;
1.447 +
1.448 + arc = curr < next ?
1.449 + std::make_pair(curr, next) : std::make_pair(next, curr);
1.450 +
1.451 + if (arcs.find(arc) == arcs.end()) {
1.452 + arcs.insert(arc);
1.453 + g.addEdge(nodes[curr], nodes[next]);
1.454 + ++cnt;
1.455 + }
1.456 + }
1.457 +}
1.458 +
1.459 +void sparse(int d)
1.460 +{
1.461 + Counter cnt("Number of arcs removed: ");
1.462 + Bfs<ListGraph> bfs(g);
1.463 + for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
1.464 + ei!=arcs.rend();++ei)
1.465 + {
1.466 + Node a=g.u(*ei);
1.467 + Node b=g.v(*ei);
1.468 + g.erase(*ei);
1.469 + bfs.run(a,b);
1.470 + if(bfs.predArc(b)==INVALID || bfs.dist(b)>d)
1.471 + g.addEdge(a,b);
1.472 + else cnt++;
1.473 + }
1.474 +}
1.475 +
1.476 +void sparse2(int d)
1.477 +{
1.478 + Counter cnt("Number of arcs removed: ");
1.479 + for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
1.480 + ei!=arcs.rend();++ei)
1.481 + {
1.482 + Node a=g.u(*ei);
1.483 + Node b=g.v(*ei);
1.484 + g.erase(*ei);
1.485 + ConstMap<Arc,int> cegy(1);
1.486 + Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy);
1.487 + int k=sur.run(a,b,2);
1.488 + if(k<2 || sur.totalLength()>d)
1.489 + g.addEdge(a,b);
1.490 + else cnt++;
1.491 +// else std::cout << "Remove arc " << g.id(a) << "-" << g.id(b) << '\n';
1.492 + }
1.493 +}
1.494 +
1.495 +void sparseTriangle(int d)
1.496 +{
1.497 + Counter cnt("Number of arcs added: ");
1.498 + std::vector<Parc> pedges;
1.499 + for(NodeIt n(g);n!=INVALID;++n)
1.500 + for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
1.501 + {
1.502 + Parc p;
1.503 + p.a=n;
1.504 + p.b=m;
1.505 + p.len=(coords[m]-coords[n]).normSquare();
1.506 + pedges.push_back(p);
1.507 + }
1.508 + std::sort(pedges.begin(),pedges.end(),pedgeLess);
1.509 + for(std::vector<Parc>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
1.510 + {
1.511 + Line li(pi->a,pi->b);
1.512 + EdgeIt e(g);
1.513 + for(;e!=INVALID && !cross(e,li);++e) ;
1.514 + Edge ne;
1.515 + if(e==INVALID) {
1.516 + ConstMap<Arc,int> cegy(1);
1.517 + Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy);
1.518 + int k=sur.run(pi->a,pi->b,2);
1.519 + if(k<2 || sur.totalLength()>d)
1.520 + {
1.521 + ne=g.addEdge(pi->a,pi->b);
1.522 + arcs.push_back(ne);
1.523 + cnt++;
1.524 + }
1.525 + }
1.526 + }
1.527 +}
1.528 +
1.529 +template <typename Graph, typename CoordMap>
1.530 +class LengthSquareMap {
1.531 +public:
1.532 + typedef typename Graph::Edge Key;
1.533 + typedef typename CoordMap::Value::Value Value;
1.534 +
1.535 + LengthSquareMap(const Graph& graph, const CoordMap& coords)
1.536 + : _graph(graph), _coords(coords) {}
1.537 +
1.538 + Value operator[](const Key& key) const {
1.539 + return (_coords[_graph.v(key)] -
1.540 + _coords[_graph.u(key)]).normSquare();
1.541 + }
1.542 +
1.543 +private:
1.544 +
1.545 + const Graph& _graph;
1.546 + const CoordMap& _coords;
1.547 +};
1.548 +
1.549 +void minTree() {
1.550 + std::vector<Parc> pedges;
1.551 + Timer T;
1.552 + std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
1.553 + delaunay();
1.554 + std::cout << T.realTime() << "s: Calculating spanning tree...\n";
1.555 + LengthSquareMap<ListGraph, ListGraph::NodeMap<Point> > ls(g, coords);
1.556 + ListGraph::EdgeMap<bool> tree(g);
1.557 + kruskal(g, ls, tree);
1.558 + std::cout << T.realTime() << "s: Removing non tree arcs...\n";
1.559 + std::vector<Edge> remove;
1.560 + for (EdgeIt e(g); e != INVALID; ++e) {
1.561 + if (!tree[e]) remove.push_back(e);
1.562 + }
1.563 + for(int i = 0; i < int(remove.size()); ++i) {
1.564 + g.erase(remove[i]);
1.565 + }
1.566 + std::cout << T.realTime() << "s: Done\n";
1.567 +}
1.568 +
1.569 +void tsp2()
1.570 +{
1.571 + std::cout << "Find a tree..." << std::endl;
1.572 +
1.573 + minTree();
1.574 +
1.575 + std::cout << "Total arc length (tree) : " << totalLen() << std::endl;
1.576 +
1.577 + std::cout << "Make it Euler..." << std::endl;
1.578 +
1.579 + {
1.580 + std::vector<Node> leafs;
1.581 + for(NodeIt n(g);n!=INVALID;++n)
1.582 + if(countIncEdges(g,n)%2==1) leafs.push_back(n);
1.583 +
1.584 +// for(unsigned int i=0;i<leafs.size();i+=2)
1.585 +// g.addArc(leafs[i],leafs[i+1]);
1.586 +
1.587 + std::vector<Parc> pedges;
1.588 + for(unsigned int i=0;i<leafs.size()-1;i++)
1.589 + for(unsigned int j=i+1;j<leafs.size();j++)
1.590 + {
1.591 + Node n=leafs[i];
1.592 + Node m=leafs[j];
1.593 + Parc p;
1.594 + p.a=n;
1.595 + p.b=m;
1.596 + p.len=(coords[m]-coords[n]).normSquare();
1.597 + pedges.push_back(p);
1.598 + }
1.599 + std::sort(pedges.begin(),pedges.end(),pedgeLess);
1.600 + for(unsigned int i=0;i<pedges.size();i++)
1.601 + if(countIncEdges(g,pedges[i].a)%2 &&
1.602 + countIncEdges(g,pedges[i].b)%2)
1.603 + g.addEdge(pedges[i].a,pedges[i].b);
1.604 + }
1.605 +
1.606 + for(NodeIt n(g);n!=INVALID;++n)
1.607 + if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 )
1.608 + std::cout << "GEBASZ!!!" << std::endl;
1.609 +
1.610 + for(EdgeIt e(g);e!=INVALID;++e)
1.611 + if(g.u(e)==g.v(e))
1.612 + std::cout << "LOOP GEBASZ!!!" << std::endl;
1.613 +
1.614 + std::cout << "Number of arcs : " << countEdges(g) << std::endl;
1.615 +
1.616 + std::cout << "Total arc length (euler) : " << totalLen() << std::endl;
1.617 +
1.618 + ListGraph::EdgeMap<Arc> enext(g);
1.619 + {
1.620 + EulerIt<ListGraph> e(g);
1.621 + Arc eo=e;
1.622 + Arc ef=e;
1.623 +// std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
1.624 + for(++e;e!=INVALID;++e)
1.625 + {
1.626 +// std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
1.627 + enext[eo]=e;
1.628 + eo=e;
1.629 + }
1.630 + enext[eo]=ef;
1.631 + }
1.632 +
1.633 + std::cout << "Creating a tour from that..." << std::endl;
1.634 +
1.635 + int nnum = countNodes(g);
1.636 + int ednum = countEdges(g);
1.637 +
1.638 + for(Arc p=enext[EdgeIt(g)];ednum>nnum;p=enext[p])
1.639 + {
1.640 +// std::cout << "Checking arc " << g.id(p) << std::endl;
1.641 + Arc e=enext[p];
1.642 + Arc f=enext[e];
1.643 + Node n2=g.source(f);
1.644 + Node n1=g.oppositeNode(n2,e);
1.645 + Node n3=g.oppositeNode(n2,f);
1.646 + if(countIncEdges(g,n2)>2)
1.647 + {
1.648 +// std::cout << "Remove an Arc" << std::endl;
1.649 + Arc ff=enext[f];
1.650 + g.erase(e);
1.651 + g.erase(f);
1.652 + if(n1!=n3)
1.653 + {
1.654 + Arc ne=g.direct(g.addEdge(n1,n3),n1);
1.655 + enext[p]=ne;
1.656 + enext[ne]=ff;
1.657 + ednum--;
1.658 + }
1.659 + else {
1.660 + enext[p]=ff;
1.661 + ednum-=2;
1.662 + }
1.663 + }
1.664 + }
1.665 +
1.666 + std::cout << "Total arc length (tour) : " << totalLen() << std::endl;
1.667 +
1.668 + std::cout << "2-opt the tour..." << std::endl;
1.669 +
1.670 + tsp_improve();
1.671 +
1.672 + std::cout << "Total arc length (2-opt tour) : " << totalLen() << std::endl;
1.673 +}
1.674 +
1.675 +
1.676 +int main(int argc,const char **argv)
1.677 +{
1.678 + ArgParser ap(argc,argv);
1.679 +
1.680 +// bool eps;
1.681 + bool disc_d, square_d, gauss_d;
1.682 +// bool tsp_a,two_a,tree_a;
1.683 + int num_of_cities=1;
1.684 + double area=1;
1.685 + N=100;
1.686 +// girth=10;
1.687 + std::string ndist("disc");
1.688 + ap.refOption("n", "Number of nodes (default is 100)", N)
1.689 + .intOption("g", "Girth parameter (default is 10)", 10)
1.690 + .refOption("cities", "Number of cities (default is 1)", num_of_cities)
1.691 + .refOption("area", "Full relative area of the cities (default is 1)", area)
1.692 + .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",
1.693 + disc_d)
1.694 + .optionGroup("dist", "disc")
1.695 + .refOption("square", "Nodes are evenly distributed on a unit square",
1.696 + square_d)
1.697 + .optionGroup("dist", "square")
1.698 + .refOption("gauss", "Nodes are located according to a two-dim Gauss "
1.699 + "distribution", gauss_d)
1.700 + .optionGroup("dist", "gauss")
1.701 + .onlyOneGroup("dist")
1.702 + .boolOption("eps", "Also generate .eps output (<prefix>.eps)")
1.703 + .boolOption("nonodes", "Draw only the edges in the generated .eps output")
1.704 + .boolOption("dir", "Directed graph is generated (each edge is replaced by "
1.705 + "two directed arcs)")
1.706 + .boolOption("2con", "Create a two connected planar graph")
1.707 + .optionGroup("alg","2con")
1.708 + .boolOption("tree", "Create a min. cost spanning tree")
1.709 + .optionGroup("alg","tree")
1.710 + .boolOption("tsp", "Create a TSP tour")
1.711 + .optionGroup("alg","tsp")
1.712 + .boolOption("tsp2", "Create a TSP tour (tree based)")
1.713 + .optionGroup("alg","tsp2")
1.714 + .boolOption("dela", "Delaunay triangulation graph")
1.715 + .optionGroup("alg","dela")
1.716 + .onlyOneGroup("alg")
1.717 + .boolOption("rand", "Use time seed for random number generator")
1.718 + .optionGroup("rand", "rand")
1.719 + .intOption("seed", "Random seed", -1)
1.720 + .optionGroup("rand", "seed")
1.721 + .onlyOneGroup("rand")
1.722 + .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
1.723 + .run();
1.724 +
1.725 + if (ap["rand"]) {
1.726 + int seed = int(time(0));
1.727 + std::cout << "Random number seed: " << seed << std::endl;
1.728 + rnd = Random(seed);
1.729 + }
1.730 + if (ap.given("seed")) {
1.731 + int seed = ap["seed"];
1.732 + std::cout << "Random number seed: " << seed << std::endl;
1.733 + rnd = Random(seed);
1.734 + }
1.735 +
1.736 + std::string prefix;
1.737 + switch(ap.files().size())
1.738 + {
1.739 + case 0:
1.740 + prefix="lgf-gen-out";
1.741 + break;
1.742 + case 1:
1.743 + prefix=ap.files()[0];
1.744 + break;
1.745 + default:
1.746 + std::cerr << "\nAt most one prefix can be given\n\n";
1.747 + exit(1);
1.748 + }
1.749 +
1.750 + double sum_sizes=0;
1.751 + std::vector<double> sizes;
1.752 + std::vector<double> cum_sizes;
1.753 + for(int s=0;s<num_of_cities;s++)
1.754 + {
1.755 + // sum_sizes+=rnd.exponential();
1.756 + double d=rnd();
1.757 + sum_sizes+=d;
1.758 + sizes.push_back(d);
1.759 + cum_sizes.push_back(sum_sizes);
1.760 + }
1.761 + int i=0;
1.762 + for(int s=0;s<num_of_cities;s++)
1.763 + {
1.764 + Point center=(num_of_cities==1?Point(0,0):rnd.disc());
1.765 + if(gauss_d)
1.766 + for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
1.767 + Node n=g.addNode();
1.768 + nodes.push_back(n);
1.769 + coords[n]=center+rnd.gauss2()*area*
1.770 + std::sqrt(sizes[s]/sum_sizes);
1.771 + }
1.772 + else if(square_d)
1.773 + for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
1.774 + Node n=g.addNode();
1.775 + nodes.push_back(n);
1.776 + coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*
1.777 + std::sqrt(sizes[s]/sum_sizes);
1.778 + }
1.779 + else if(disc_d || true)
1.780 + for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
1.781 + Node n=g.addNode();
1.782 + nodes.push_back(n);
1.783 + coords[n]=center+rnd.disc()*area*
1.784 + std::sqrt(sizes[s]/sum_sizes);
1.785 + }
1.786 + }
1.787 +
1.788 +// for (ListGraph::NodeIt n(g); n != INVALID; ++n) {
1.789 +// std::cerr << coords[n] << std::endl;
1.790 +// }
1.791 +
1.792 + if(ap["tsp"]) {
1.793 + tsp();
1.794 + std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
1.795 + }
1.796 + if(ap["tsp2"]) {
1.797 + tsp2();
1.798 + std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
1.799 + }
1.800 + else if(ap["2con"]) {
1.801 + std::cout << "Make triangles\n";
1.802 + // triangle();
1.803 + sparseTriangle(ap["g"]);
1.804 + std::cout << "Make it sparser\n";
1.805 + sparse2(ap["g"]);
1.806 + }
1.807 + else if(ap["tree"]) {
1.808 + minTree();
1.809 + }
1.810 + else if(ap["dela"]) {
1.811 + delaunay();
1.812 + }
1.813 +
1.814 +
1.815 + std::cout << "Number of nodes : " << countNodes(g) << std::endl;
1.816 + std::cout << "Number of arcs : " << countEdges(g) << std::endl;
1.817 + double tlen=0;
1.818 + for(EdgeIt e(g);e!=INVALID;++e)
1.819 + tlen+=std::sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
1.820 + std::cout << "Total arc length : " << tlen << std::endl;
1.821 +
1.822 + if(ap["eps"])
1.823 + graphToEps(g,prefix+".eps").scaleToA4().
1.824 + scale(600).nodeScale(.005).arcWidthScale(.001).preScale(false).
1.825 + coords(coords).hideNodes(ap.given("nonodes")).run();
1.826 +
1.827 + if(ap["dir"])
1.828 + DigraphWriter<ListGraph>(g,prefix+".lgf").
1.829 + nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
1.830 + nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
1.831 + run();
1.832 + else GraphWriter<ListGraph>(g,prefix+".lgf").
1.833 + nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
1.834 + nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
1.835 + run();
1.836 +}
1.837 +