COIN-OR::LEMON - Graph Library

source: lemon-main/tools/lgf-gen.cc @ 903:3ffd46dc8e01

Last change on this file since 903:3ffd46dc8e01 was 654:a312f84d86c6, checked in by Peter Kovacs <kpeter@…>, 16 years ago

Doc fixes for lgf-gen (#282)

File size: 22.3 KB
RevLine 
[523]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
[654]21/// \brief Special plane graph generator.
[523]22///
23/// Graph generator application for various types of plane graphs.
24///
25/// See
[584]26/// \code
27///   lgf-gen --help
28/// \endcode
[654]29/// for more information on the usage.
[523]30
31#include <algorithm>
32#include <set>
[570]33#include <ctime>
[523]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)
[612]68    tlen+=std::sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
[523]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
[612]191    return d / (2 * a) + std::sqrt((d * d + e * e) / (4 * a * a) + f / a);
[523]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();
[612]209    return (b - std::sqrt(d)) / a;
[523]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);
[623]483      Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy);
484      int k=sur.run(a,b,2);
[523]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);
[623]514        Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy);
515        int k=sur.run(pi->a,pi->b,2);
[523]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)
[654]689    .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",
690               disc_d)
[523]691    .optionGroup("dist", "disc")
[654]692    .refOption("square", "Nodes are evenly distributed on a unit square",
693               square_d)
[523]694    .optionGroup("dist", "square")
[654]695    .refOption("gauss", "Nodes are located according to a two-dim Gauss "
696               "distribution", gauss_d)
[523]697    .optionGroup("dist", "gauss")
698    .onlyOneGroup("dist")
[654]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")
[523]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")
[654]711    .boolOption("dela", "Delaunay triangulation graph")
[523]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"]) {
[616]723    int seed = int(time(0));
[523]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)
[612]816    tlen+=std::sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
[523]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).
[524]822      coords(coords).hideNodes(ap.given("nonodes")).run();
[523]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.