COIN-OR::LEMON - Graph Library

source: lemon/tools/lgf-gen.cc @ 670:7c1324b35d89

Last change on this file since 670:7c1324b35d89 was 670:7c1324b35d89, checked in by Peter Kovacs <kpeter@…>, 10 years ago

Modify the interface of Suurballe (#266, #181)

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