COIN-OR::LEMON - Graph Library

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

Last change on this file was 2618:6aa6fcaeaea5, checked in by Balazs Dezso, 11 years ago

G++-4.3 compatibility changes

File size: 22.9 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2008
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 */
18
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
72
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>
83#include <lemon/euler.h>
84#include <lemon/math.h>
85#include <algorithm>
86#include <lemon/kruskal.h>
87#include <lemon/time_measure.h>
88
89using namespace lemon;
90
91typedef dim2::Point<double> Point;
92
93UGRAPH_TYPEDEFS(ListUGraph);
94
95bool progress=true;
96
97int N;
98// int girth;
99
100ListUGraph g;
101
102std::vector<Node> nodes;
103ListUGraph::NodeMap<Point> coords(g);
104
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
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
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() {
303  Counter cnt("Number of edges added: ");
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
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  }
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
416      {
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        }
436      }
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));
461      }
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     
475    }
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  }
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
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
588void minTree() {
589  std::vector<Pedge> pedges;
590  Timer T;
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  }
605  std::cout << T.realTime() << "s: Done\n";
606}
607
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);
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);
643  }
644
645  for(NodeIt n(g);n!=INVALID;++n)
646    if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 )
647      std::cout << "GEBASZ!!!" << std::endl;
648 
649  for(UEdgeIt e(g);e!=INVALID;++e)
650    if(g.source(e)==g.target(e))
651      std::cout << "LOOP GEBASZ!!!" << std::endl;
652 
653  std::cout << "Number of edges : " << countUEdges(g) << std::endl;
654 
655  std::cout << "Total edge length (euler) : " << totalLen() << std::endl;
656
657  ListUGraph::UEdgeMap<Edge> enext(g);
658  {
659    UEulerIt<ListUGraph> e(g);
660    Edge eo=e;
661    Edge ef=e;
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  }
671   
672  std::cout << "Creating a tour from that..." << std::endl;
673 
674  int nnum = countNodes(g);
675  int ednum = countUEdges(g);
676 
677  for(Edge p=enext[UEdgeIt(g)];ednum>nnum;p=enext[p])
678    {
679//       std::cout << "Checking edge " << g.id(p) << std::endl;     
680      Edge e=enext[p];
681      Edge f=enext[e];
682      Node n2=g.source(f);
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;
688          Edge ff=enext[f];
689          g.erase(e);
690          g.erase(f);
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          }
702        }
703    }
704
705  std::cout << "Total edge length (tour) : " << totalLen() << std::endl;
706
707  std::cout << "2-opt the tour..." << std::endl;
708 
709  tsp_improve();
710 
711  std::cout << "Total edge length (2-opt tour) : " << totalLen() << std::endl;
712}
713
714
715int main(int argc,const char **argv)
716{
717  ArgParser ap(argc,argv);
718
719//   bool eps;
720  bool disc_d, square_d, gauss_d;
721//   bool tsp_a,two_a,tree_a;
722  int num_of_cities=1;
723  double area=1;
724  N=100;
725//   girth=10;
726  std::string ndist("disc");
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)
732    .optionGroup("dist", "disc")
733    .refOption("square", "Nodes are evenly distributed on a unit square", square_d)
734    .optionGroup("dist", "square")
735    .refOption("gauss",
736            "Nodes are located according to a two-dim gauss distribution",
737            gauss_d)
738    .optionGroup("dist", "gauss")
739//     .mandatoryGroup("dist")
740    .onlyOneGroup("dist")
741    .boolOption("eps", "Also generate .eps output (prefix.eps)")
742    .boolOption("dir", "Directed graph is generated (each edges are replaced by two directed ones)")
743    .boolOption("2con", "Create a two connected planar graph")
744    .optionGroup("alg","2con")
745    .boolOption("tree", "Create a min. cost spanning tree")
746    .optionGroup("alg","tree")
747    .boolOption("tsp", "Create a TSP tour")
748    .optionGroup("alg","tsp")
749    .boolOption("tsp2", "Create a TSP tour (tree based)")
750    .optionGroup("alg","tsp2")
751    .boolOption("dela", "Delaunay triangulation graph")
752    .optionGroup("alg","dela")
753    .onlyOneGroup("alg")
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")
759    .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
760    .run();
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  }
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    }
824
825//   for (ListUGraph::NodeIt n(g); n != INVALID; ++n) {
826//     std::cerr << coords[n] << std::endl;
827//   }
828 
829  if(ap["tsp"]) {
830    tsp();
831    std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
832  }
833  if(ap["tsp2"]) {
834    tsp2();
835    std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
836  }
837  else if(ap["2con"]) {
838    std::cout << "Make triangles\n";
839    //   triangle();
840    sparseTriangle(ap["g"]);
841    std::cout << "Make it sparser\n";
842    sparse2(ap["g"]);
843  }
844  else if(ap["tree"]) {
845    minTree();
846  }
847  else if(ap["dela"]) {
848    delaunay();
849  }
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;
858
859  if(ap["eps"])
860    graphToEps(g,prefix+".eps").scaleToA4().
861      scale(600).nodeScale(.2).edgeWidthScale(.001).preScale(false).
862      coords(coords).run();
863 
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();
873}
874
Note: See TracBrowser for help on using the repository browser.