COIN-OR::LEMON - Graph Library

source: lemon/tools/lgf-gen.cc @ 1030:a80381c43760

Last change on this file since 1030:a80381c43760 was 701:a312f84d86c6, checked in by Peter Kovacs <kpeter@…>, 15 years ago

Doc fixes for lgf-gen (#282)

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 graph 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 information 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)",
690               disc_d)
691    .optionGroup("dist", "disc")
692    .refOption("square", "Nodes are evenly distributed on a unit square",
693               square_d)
694    .optionGroup("dist", "square")
695    .refOption("gauss", "Nodes are located according to a two-dim Gauss "
696               "distribution", gauss_d)
697    .optionGroup("dist", "gauss")
698    .onlyOneGroup("dist")
699    .boolOption("eps", "Also generate .eps output (<prefix>.eps)")
700    .boolOption("nonodes", "Draw only the edges in the generated .eps output")
701    .boolOption("dir", "Directed graph is generated (each edge is replaced by "
702                "two directed arcs)")
703    .boolOption("2con", "Create a two connected planar graph")
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 graph")
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 = int(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+=std::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.