COIN-OR::LEMON - Graph Library

source: lemon-1.2/tools/lgf-gen.cc @ 614:19b6f20e0ea2

Last change on this file since 614:19b6f20e0ea2 was 584:33c6b6e755cd, checked in by Peter Kovacs <kpeter@…>, 11 years ago

Small doc improvements (#263)

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+=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) + 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 - 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,a,b);
484      int k=sur.run(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> >
515          sur(g,cegy,pi->a,pi->b);
516        int k=sur.run(2);
517        if(k<2 || sur.totalLength()>d)
518          {
519            ne=g.addEdge(pi->a,pi->b);
520            arcs.push_back(ne);
521            cnt++;
522          }
523      }
524    }
525}
526
527template <typename Graph, typename CoordMap>
528class LengthSquareMap {
529public:
530  typedef typename Graph::Edge Key;
531  typedef typename CoordMap::Value::Value Value;
532
533  LengthSquareMap(const Graph& graph, const CoordMap& coords)
534    : _graph(graph), _coords(coords) {}
535
536  Value operator[](const Key& key) const {
537    return (_coords[_graph.v(key)] -
538            _coords[_graph.u(key)]).normSquare();
539  }
540
541private:
542
543  const Graph& _graph;
544  const CoordMap& _coords;
545};
546
547void minTree() {
548  std::vector<Parc> pedges;
549  Timer T;
550  std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
551  delaunay();
552  std::cout << T.realTime() << "s: Calculating spanning tree...\n";
553  LengthSquareMap<ListGraph, ListGraph::NodeMap<Point> > ls(g, coords);
554  ListGraph::EdgeMap<bool> tree(g);
555  kruskal(g, ls, tree);
556  std::cout << T.realTime() << "s: Removing non tree arcs...\n";
557  std::vector<Edge> remove;
558  for (EdgeIt e(g); e != INVALID; ++e) {
559    if (!tree[e]) remove.push_back(e);
560  }
561  for(int i = 0; i < int(remove.size()); ++i) {
562    g.erase(remove[i]);
563  }
564  std::cout << T.realTime() << "s: Done\n";
565}
566
567void tsp2()
568{
569  std::cout << "Find a tree..." << std::endl;
570
571  minTree();
572
573  std::cout << "Total arc length (tree) : " << totalLen() << std::endl;
574
575  std::cout << "Make it Euler..." << std::endl;
576
577  {
578    std::vector<Node> leafs;
579    for(NodeIt n(g);n!=INVALID;++n)
580      if(countIncEdges(g,n)%2==1) leafs.push_back(n);
581
582//    for(unsigned int i=0;i<leafs.size();i+=2)
583//       g.addArc(leafs[i],leafs[i+1]);
584
585    std::vector<Parc> pedges;
586    for(unsigned int i=0;i<leafs.size()-1;i++)
587      for(unsigned int j=i+1;j<leafs.size();j++)
588        {
589          Node n=leafs[i];
590          Node m=leafs[j];
591          Parc p;
592          p.a=n;
593          p.b=m;
594          p.len=(coords[m]-coords[n]).normSquare();
595          pedges.push_back(p);
596        }
597    std::sort(pedges.begin(),pedges.end(),pedgeLess);
598    for(unsigned int i=0;i<pedges.size();i++)
599      if(countIncEdges(g,pedges[i].a)%2 &&
600         countIncEdges(g,pedges[i].b)%2)
601        g.addEdge(pedges[i].a,pedges[i].b);
602  }
603
604  for(NodeIt n(g);n!=INVALID;++n)
605    if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 )
606      std::cout << "GEBASZ!!!" << std::endl;
607
608  for(EdgeIt e(g);e!=INVALID;++e)
609    if(g.u(e)==g.v(e))
610      std::cout << "LOOP GEBASZ!!!" << std::endl;
611
612  std::cout << "Number of arcs : " << countEdges(g) << std::endl;
613
614  std::cout << "Total arc length (euler) : " << totalLen() << std::endl;
615
616  ListGraph::EdgeMap<Arc> enext(g);
617  {
618    EulerIt<ListGraph> e(g);
619    Arc eo=e;
620    Arc ef=e;
621//     std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
622    for(++e;e!=INVALID;++e)
623      {
624//         std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
625        enext[eo]=e;
626        eo=e;
627      }
628    enext[eo]=ef;
629  }
630
631  std::cout << "Creating a tour from that..." << std::endl;
632
633  int nnum = countNodes(g);
634  int ednum = countEdges(g);
635
636  for(Arc p=enext[EdgeIt(g)];ednum>nnum;p=enext[p])
637    {
638//       std::cout << "Checking arc " << g.id(p) << std::endl;
639      Arc e=enext[p];
640      Arc f=enext[e];
641      Node n2=g.source(f);
642      Node n1=g.oppositeNode(n2,e);
643      Node n3=g.oppositeNode(n2,f);
644      if(countIncEdges(g,n2)>2)
645        {
646//           std::cout << "Remove an Arc" << std::endl;
647          Arc ff=enext[f];
648          g.erase(e);
649          g.erase(f);
650          if(n1!=n3)
651            {
652              Arc ne=g.direct(g.addEdge(n1,n3),n1);
653              enext[p]=ne;
654              enext[ne]=ff;
655              ednum--;
656            }
657          else {
658            enext[p]=ff;
659            ednum-=2;
660          }
661        }
662    }
663
664  std::cout << "Total arc length (tour) : " << totalLen() << std::endl;
665
666  std::cout << "2-opt the tour..." << std::endl;
667
668  tsp_improve();
669
670  std::cout << "Total arc length (2-opt tour) : " << totalLen() << std::endl;
671}
672
673
674int main(int argc,const char **argv)
675{
676  ArgParser ap(argc,argv);
677
678//   bool eps;
679  bool disc_d, square_d, gauss_d;
680//   bool tsp_a,two_a,tree_a;
681  int num_of_cities=1;
682  double area=1;
683  N=100;
684//   girth=10;
685  std::string ndist("disc");
686  ap.refOption("n", "Number of nodes (default is 100)", N)
687    .intOption("g", "Girth parameter (default is 10)", 10)
688    .refOption("cities", "Number of cities (default is 1)", num_of_cities)
689    .refOption("area", "Full relative area of the cities (default is 1)", area)
690    .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",disc_d)
691    .optionGroup("dist", "disc")
692    .refOption("square", "Nodes are evenly distributed on a unit square", square_d)
693    .optionGroup("dist", "square")
694    .refOption("gauss",
695            "Nodes are located according to a two-dim gauss distribution",
696            gauss_d)
697    .optionGroup("dist", "gauss")
698//     .mandatoryGroup("dist")
699    .onlyOneGroup("dist")
700    .boolOption("eps", "Also generate .eps output (prefix.eps)")
701    .boolOption("nonodes", "Draw the edges only in the generated .eps")
702    .boolOption("dir", "Directed digraph is generated (each arcs are replaced by two directed ones)")
703    .boolOption("2con", "Create a two connected planar digraph")
704    .optionGroup("alg","2con")
705    .boolOption("tree", "Create a min. cost spanning tree")
706    .optionGroup("alg","tree")
707    .boolOption("tsp", "Create a TSP tour")
708    .optionGroup("alg","tsp")
709    .boolOption("tsp2", "Create a TSP tour (tree based)")
710    .optionGroup("alg","tsp2")
711    .boolOption("dela", "Delaunay triangulation digraph")
712    .optionGroup("alg","dela")
713    .onlyOneGroup("alg")
714    .boolOption("rand", "Use time seed for random number generator")
715    .optionGroup("rand", "rand")
716    .intOption("seed", "Random seed", -1)
717    .optionGroup("rand", "seed")
718    .onlyOneGroup("rand")
719    .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
720    .run();
721
722  if (ap["rand"]) {
723    int seed = time(0);
724    std::cout << "Random number seed: " << seed << std::endl;
725    rnd = Random(seed);
726  }
727  if (ap.given("seed")) {
728    int seed = ap["seed"];
729    std::cout << "Random number seed: " << seed << std::endl;
730    rnd = Random(seed);
731  }
732
733  std::string prefix;
734  switch(ap.files().size())
735    {
736    case 0:
737      prefix="lgf-gen-out";
738      break;
739    case 1:
740      prefix=ap.files()[0];
741      break;
742    default:
743      std::cerr << "\nAt most one prefix can be given\n\n";
744      exit(1);
745    }
746
747  double sum_sizes=0;
748  std::vector<double> sizes;
749  std::vector<double> cum_sizes;
750  for(int s=0;s<num_of_cities;s++)
751    {
752      //         sum_sizes+=rnd.exponential();
753      double d=rnd();
754      sum_sizes+=d;
755      sizes.push_back(d);
756      cum_sizes.push_back(sum_sizes);
757    }
758  int i=0;
759  for(int s=0;s<num_of_cities;s++)
760    {
761      Point center=(num_of_cities==1?Point(0,0):rnd.disc());
762      if(gauss_d)
763        for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
764          Node n=g.addNode();
765          nodes.push_back(n);
766          coords[n]=center+rnd.gauss2()*area*
767            std::sqrt(sizes[s]/sum_sizes);
768        }
769      else if(square_d)
770        for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
771          Node n=g.addNode();
772          nodes.push_back(n);
773          coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*
774            std::sqrt(sizes[s]/sum_sizes);
775        }
776      else if(disc_d || true)
777        for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
778          Node n=g.addNode();
779          nodes.push_back(n);
780          coords[n]=center+rnd.disc()*area*
781            std::sqrt(sizes[s]/sum_sizes);
782        }
783    }
784
785//   for (ListGraph::NodeIt n(g); n != INVALID; ++n) {
786//     std::cerr << coords[n] << std::endl;
787//   }
788
789  if(ap["tsp"]) {
790    tsp();
791    std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
792  }
793  if(ap["tsp2"]) {
794    tsp2();
795    std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
796  }
797  else if(ap["2con"]) {
798    std::cout << "Make triangles\n";
799    //   triangle();
800    sparseTriangle(ap["g"]);
801    std::cout << "Make it sparser\n";
802    sparse2(ap["g"]);
803  }
804  else if(ap["tree"]) {
805    minTree();
806  }
807  else if(ap["dela"]) {
808    delaunay();
809  }
810
811
812  std::cout << "Number of nodes    : " << countNodes(g) << std::endl;
813  std::cout << "Number of arcs    : " << countEdges(g) << std::endl;
814  double tlen=0;
815  for(EdgeIt e(g);e!=INVALID;++e)
816    tlen+=sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
817  std::cout << "Total arc length  : " << tlen << std::endl;
818
819  if(ap["eps"])
820    graphToEps(g,prefix+".eps").scaleToA4().
821      scale(600).nodeScale(.005).arcWidthScale(.001).preScale(false).
822      coords(coords).hideNodes(ap.given("nonodes")).run();
823
824  if(ap["dir"])
825    DigraphWriter<ListGraph>(g,prefix+".lgf").
826      nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
827      nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
828      run();
829  else GraphWriter<ListGraph>(g,prefix+".lgf").
830         nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
831         nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
832         run();
833}
834
Note: See TracBrowser for help on using the repository browser.