COIN-OR::LEMON - Graph Library

source: lemon-0.x/tools/lgf-gen.cc @ 2552:5f711e4668f5

Last change on this file since 2552:5f711e4668f5 was 2493:6231d9d3957b, checked in by Balazs Dezso, 12 years ago

Bad documentation

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