COIN-OR::LEMON - Graph Library

source: lemon-0.x/tools/lgf-gen.cc @ 2553:bfced05fa852

Last change on this file since 2553:bfced05fa852 was 2553:bfced05fa852, checked in by Alpar Juttner, 12 years ago

Happy New Year to LEMON (+ better update-copyright-header script)

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