COIN-OR::LEMON - Graph Library

source: lemon/tools/lgf-gen.cc @ 570:d9e43511d11c

Last change on this file since 570:d9e43511d11c was 570:d9e43511d11c, checked in by Alpar Juttner <alpar@…>, 15 years ago

Port lgf-gen from SVN -r3512 (#45)

  • apply the migrate script
  • apply the source unifyer
  • fix the compilation
File size: 22.2 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/// \verbatim
27///  lgf-gen --help
28/// \endverbatim
29/// for more info on the usage.
30///
31
32
33#include <algorithm>
34#include <set>
35#include <lemon/list_graph.h>
36#include <lemon/random.h>
37#include <lemon/dim2.h>
38#include <lemon/bfs.h>
39#include <lemon/counter.h>
40#include <lemon/suurballe.h>
41#include <lemon/graph_to_eps.h>
42#include <lemon/lgf_writer.h>
43#include <lemon/arg_parser.h>
44#include <lemon/euler.h>
45#include <lemon/math.h>
46#include <lemon/kruskal.h>
47#include <lemon/time_measure.h>
48
49using namespace lemon;
50
51typedef dim2::Point<double> Point;
52
53GRAPH_TYPEDEFS(ListGraph);
54
55bool progress=true;
56
57int N;
58// int girth;
59
60ListGraph g;
61
62std::vector<Node> nodes;
63ListGraph::NodeMap<Point> coords(g);
64
65
66double totalLen(){
67  double tlen=0;
68  for(EdgeIt e(g);e!=INVALID;++e)
69    tlen+=sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
70  return tlen;
71}
72
73int tsp_impr_num=0;
74
75const double EPSILON=1e-8;
76bool tsp_improve(Node u, Node v)
77{
78  double luv=std::sqrt((coords[v]-coords[u]).normSquare());
79  Node u2=u;
80  Node v2=v;
81  do {
82    Node n;
83    for(IncEdgeIt e(g,v2);(n=g.runningNode(e))==u2;++e) { }
84    u2=v2;
85    v2=n;
86    if(luv+std::sqrt((coords[v2]-coords[u2]).normSquare())-EPSILON>
87       std::sqrt((coords[u]-coords[u2]).normSquare())+
88       std::sqrt((coords[v]-coords[v2]).normSquare()))
89      {
90         g.erase(findEdge(g,u,v));
91         g.erase(findEdge(g,u2,v2));
92        g.addEdge(u2,u);
93        g.addEdge(v,v2);
94        tsp_impr_num++;
95        return true;
96      }
97  } while(v2!=u);
98  return false;
99}
100
101bool tsp_improve(Node u)
102{
103  for(IncEdgeIt e(g,u);e!=INVALID;++e)
104    if(tsp_improve(u,g.runningNode(e))) return true;
105  return false;
106}
107
108void tsp_improve()
109{
110  bool b;
111  do {
112    b=false;
113    for(NodeIt n(g);n!=INVALID;++n)
114      if(tsp_improve(n)) b=true;
115  } while(b);
116}
117
118void tsp()
119{
120  for(int i=0;i<N;i++) g.addEdge(nodes[i],nodes[(i+1)%N]);
121  tsp_improve();
122}
123
124class Line
125{
126public:
127  Point a;
128  Point b;
129  Line(Point _a,Point _b) :a(_a),b(_b) {}
130  Line(Node _a,Node _b) : a(coords[_a]),b(coords[_b]) {}
131  Line(const Arc &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
132  Line(const Edge &e) : a(coords[g.u(e)]),b(coords[g.v(e)]) {}
133};
134
135inline std::ostream& operator<<(std::ostream &os, const Line &l)
136{
137  os << l.a << "->" << l.b;
138  return os;
139}
140
141bool cross(Line a, Line b)
142{
143  Point ao=rot90(a.b-a.a);
144  Point bo=rot90(b.b-b.a);
145  return (ao*(b.a-a.a))*(ao*(b.b-a.a))<0 &&
146    (bo*(a.a-b.a))*(bo*(a.b-b.a))<0;
147}
148
149struct Parc
150{
151  Node a;
152  Node b;
153  double len;
154};
155
156bool pedgeLess(Parc a,Parc b)
157{
158  return a.len<b.len;
159}
160
161std::vector<Edge> arcs;
162
163namespace _delaunay_bits {
164
165  struct Part {
166    int prev, curr, next;
167
168    Part(int p, int c, int n) : prev(p), curr(c), next(n) {}
169  };
170
171  inline std::ostream& operator<<(std::ostream& os, const Part& part) {
172    os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')';
173    return os;
174  }
175
176  inline double circle_point(const Point& p, const Point& q, const Point& r) {
177    double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y);
178    if (a == 0) return std::numeric_limits<double>::quiet_NaN();
179
180    double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) +
181      (q.x * q.x + q.y * q.y) * (r.y - p.y) +
182      (r.x * r.x + r.y * r.y) * (p.y - q.y);
183
184    double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) +
185      (q.x * q.x + q.y * q.y) * (r.x - p.x) +
186      (r.x * r.x + r.y * r.y) * (p.x - q.x);
187
188    double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) +
189      (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) +
190      (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y);
191
192    return d / (2 * a) + sqrt((d * d + e * e) / (4 * a * a) + f / a);
193  }
194
195  inline bool circle_form(const Point& p, const Point& q, const Point& r) {
196    return rot90(q - p) * (r - q) < 0.0;
197  }
198
199  inline double intersection(const Point& p, const Point& q, double sx) {
200    const double epsilon = 1e-8;
201
202    if (p.x == q.x) return (p.y + q.y) / 2.0;
203
204    if (sx < p.x + epsilon) return p.y;
205    if (sx < q.x + epsilon) return q.y;
206
207    double a = q.x - p.x;
208    double b = (q.x - sx) * p.y - (p.x - sx) * q.y;
209    double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare();
210    return (b - sqrt(d)) / a;
211  }
212
213  struct YLess {
214
215
216    YLess(const std::vector<Point>& points, double& sweep)
217      : _points(points), _sweep(sweep) {}
218
219    bool operator()(const Part& l, const Part& r) const {
220      const double epsilon = 1e-8;
221
222      //      std::cerr << l << " vs " << r << std::endl;
223      double lbx = l.prev != -1 ?
224        intersection(_points[l.prev], _points[l.curr], _sweep) :
225        - std::numeric_limits<double>::infinity();
226      double rbx = r.prev != -1 ?
227        intersection(_points[r.prev], _points[r.curr], _sweep) :
228        - std::numeric_limits<double>::infinity();
229      double lex = l.next != -1 ?
230        intersection(_points[l.curr], _points[l.next], _sweep) :
231        std::numeric_limits<double>::infinity();
232      double rex = r.next != -1 ?
233        intersection(_points[r.curr], _points[r.next], _sweep) :
234        std::numeric_limits<double>::infinity();
235
236      if (lbx > lex) std::swap(lbx, lex);
237      if (rbx > rex) std::swap(rbx, rex);
238
239      if (lex < epsilon + rex && lbx + epsilon < rex) return true;
240      if (rex < epsilon + lex && rbx + epsilon < lex) return false;
241      return lex < rex;
242    }
243
244    const std::vector<Point>& _points;
245    double& _sweep;
246  };
247
248  struct BeachIt;
249
250  typedef std::multimap<double, BeachIt> SpikeHeap;
251
252  typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach;
253
254  struct BeachIt {
255    Beach::iterator it;
256
257    BeachIt(Beach::iterator iter) : it(iter) {}
258  };
259
260}
261
262inline void delaunay() {
263  Counter cnt("Number of arcs added: ");
264
265  using namespace _delaunay_bits;
266
267  typedef _delaunay_bits::Part Part;
268  typedef std::vector<std::pair<double, int> > SiteHeap;
269
270
271  std::vector<Point> points;
272  std::vector<Node> nodes;
273
274  for (NodeIt it(g); it != INVALID; ++it) {
275    nodes.push_back(it);
276    points.push_back(coords[it]);
277  }
278
279  SiteHeap siteheap(points.size());
280
281  double sweep;
282
283
284  for (int i = 0; i < int(siteheap.size()); ++i) {
285    siteheap[i] = std::make_pair(points[i].x, i);
286  }
287
288  std::sort(siteheap.begin(), siteheap.end());
289  sweep = siteheap.front().first;
290
291  YLess yless(points, sweep);
292  Beach beach(yless);
293
294  SpikeHeap spikeheap;
295
296  std::set<std::pair<int, int> > arcs;
297
298  int siteindex = 0;
299  {
300    SiteHeap front;
301
302    while (siteindex < int(siteheap.size()) &&
303           siteheap[0].first == siteheap[siteindex].first) {
304      front.push_back(std::make_pair(points[siteheap[siteindex].second].y,
305                                     siteheap[siteindex].second));
306      ++siteindex;
307    }
308
309    std::sort(front.begin(), front.end());
310
311    for (int i = 0; i < int(front.size()); ++i) {
312      int prev = (i == 0 ? -1 : front[i - 1].second);
313      int curr = front[i].second;
314      int next = (i + 1 == int(front.size()) ? -1 : front[i + 1].second);
315
316      beach.insert(std::make_pair(Part(prev, curr, next),
317                                  spikeheap.end()));
318    }
319  }
320
321  while (siteindex < int(points.size()) || !spikeheap.empty()) {
322
323    SpikeHeap::iterator spit = spikeheap.begin();
324
325    if (siteindex < int(points.size()) &&
326        (spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) {
327      int site = siteheap[siteindex].second;
328      sweep = siteheap[siteindex].first;
329
330      Beach::iterator bit = beach.upper_bound(Part(site, site, site));
331
332      if (bit->second != spikeheap.end()) {
333        spikeheap.erase(bit->second);
334      }
335
336      int prev = bit->first.prev;
337      int curr = bit->first.curr;
338      int next = bit->first.next;
339
340      beach.erase(bit);
341
342      SpikeHeap::iterator pit = spikeheap.end();
343      if (prev != -1 &&
344          circle_form(points[prev], points[curr], points[site])) {
345        double x = circle_point(points[prev], points[curr], points[site]);
346        pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
347        pit->second.it =
348          beach.insert(std::make_pair(Part(prev, curr, site), pit));
349      } else {
350        beach.insert(std::make_pair(Part(prev, curr, site), pit));
351      }
352
353      beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end()));
354
355      SpikeHeap::iterator nit = spikeheap.end();
356      if (next != -1 &&
357          circle_form(points[site], points[curr],points[next])) {
358        double x = circle_point(points[site], points[curr], points[next]);
359        nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
360        nit->second.it =
361          beach.insert(std::make_pair(Part(site, curr, next), nit));
362      } else {
363        beach.insert(std::make_pair(Part(site, curr, next), nit));
364      }
365
366      ++siteindex;
367    } else {
368      sweep = spit->first;
369
370      Beach::iterator bit = spit->second.it;
371
372      int prev = bit->first.prev;
373      int curr = bit->first.curr;
374      int next = bit->first.next;
375
376      {
377        std::pair<int, int> arc;
378
379        arc = prev < curr ?
380          std::make_pair(prev, curr) : std::make_pair(curr, prev);
381
382        if (arcs.find(arc) == arcs.end()) {
383          arcs.insert(arc);
384          g.addEdge(nodes[prev], nodes[curr]);
385          ++cnt;
386        }
387
388        arc = curr < next ?
389          std::make_pair(curr, next) : std::make_pair(next, curr);
390
391        if (arcs.find(arc) == arcs.end()) {
392          arcs.insert(arc);
393          g.addEdge(nodes[curr], nodes[next]);
394          ++cnt;
395        }
396      }
397
398      Beach::iterator pbit = bit; --pbit;
399      int ppv = pbit->first.prev;
400      Beach::iterator nbit = bit; ++nbit;
401      int nnt = nbit->first.next;
402
403      if (bit->second != spikeheap.end()) spikeheap.erase(bit->second);
404      if (pbit->second != spikeheap.end()) spikeheap.erase(pbit->second);
405      if (nbit->second != spikeheap.end()) spikeheap.erase(nbit->second);
406
407      beach.erase(nbit);
408      beach.erase(bit);
409      beach.erase(pbit);
410
411      SpikeHeap::iterator pit = spikeheap.end();
412      if (ppv != -1 && ppv != next &&
413          circle_form(points[ppv], points[prev], points[next])) {
414        double x = circle_point(points[ppv], points[prev], points[next]);
415        if (x < sweep) x = sweep;
416        pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
417        pit->second.it =
418          beach.insert(std::make_pair(Part(ppv, prev, next), pit));
419      } else {
420        beach.insert(std::make_pair(Part(ppv, prev, next), pit));
421      }
422
423      SpikeHeap::iterator nit = spikeheap.end();
424      if (nnt != -1 && prev != nnt &&
425          circle_form(points[prev], points[next], points[nnt])) {
426        double x = circle_point(points[prev], points[next], points[nnt]);
427        if (x < sweep) x = sweep;
428        nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
429        nit->second.it =
430          beach.insert(std::make_pair(Part(prev, next, nnt), nit));
431      } else {
432        beach.insert(std::make_pair(Part(prev, next, nnt), nit));
433      }
434
435    }
436  }
437
438  for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {
439    int curr = it->first.curr;
440    int next = it->first.next;
441
442    if (next == -1) continue;
443
444    std::pair<int, int> arc;
445
446    arc = curr < next ?
447      std::make_pair(curr, next) : std::make_pair(next, curr);
448
449    if (arcs.find(arc) == arcs.end()) {
450      arcs.insert(arc);
451      g.addEdge(nodes[curr], nodes[next]);
452      ++cnt;
453    }
454  }
455}
456
457void sparse(int d)
458{
459  Counter cnt("Number of arcs removed: ");
460  Bfs<ListGraph> bfs(g);
461  for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
462      ei!=arcs.rend();++ei)
463    {
464      Node a=g.u(*ei);
465      Node b=g.v(*ei);
466      g.erase(*ei);
467      bfs.run(a,b);
468      if(bfs.predArc(b)==INVALID || bfs.dist(b)>d)
469        g.addEdge(a,b);
470      else cnt++;
471    }
472}
473
474void sparse2(int d)
475{
476  Counter cnt("Number of arcs removed: ");
477  for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
478      ei!=arcs.rend();++ei)
479    {
480      Node a=g.u(*ei);
481      Node b=g.v(*ei);
482      g.erase(*ei);
483      ConstMap<Arc,int> cegy(1);
484      Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy,a,b);
485      int k=sur.run(2);
486      if(k<2 || sur.totalLength()>d)
487        g.addEdge(a,b);
488      else cnt++;
489//       else std::cout << "Remove arc " << g.id(a) << "-" << g.id(b) << '\n';
490    }
491}
492
493void sparseTriangle(int d)
494{
495  Counter cnt("Number of arcs added: ");
496  std::vector<Parc> pedges;
497  for(NodeIt n(g);n!=INVALID;++n)
498    for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
499      {
500        Parc p;
501        p.a=n;
502        p.b=m;
503        p.len=(coords[m]-coords[n]).normSquare();
504        pedges.push_back(p);
505      }
506  std::sort(pedges.begin(),pedges.end(),pedgeLess);
507  for(std::vector<Parc>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
508    {
509      Line li(pi->a,pi->b);
510      EdgeIt e(g);
511      for(;e!=INVALID && !cross(e,li);++e) ;
512      Edge ne;
513      if(e==INVALID) {
514        ConstMap<Arc,int> cegy(1);
515        Suurballe<ListGraph,ConstMap<Arc,int> >
516          sur(g,cegy,pi->a,pi->b);
517        int k=sur.run(2);
518        if(k<2 || sur.totalLength()>d)
519          {
520            ne=g.addEdge(pi->a,pi->b);
521            arcs.push_back(ne);
522            cnt++;
523          }
524      }
525    }
526}
527
528template <typename Graph, typename CoordMap>
529class LengthSquareMap {
530public:
531  typedef typename Graph::Edge Key;
532  typedef typename CoordMap::Value::Value Value;
533
534  LengthSquareMap(const Graph& graph, const CoordMap& coords)
535    : _graph(graph), _coords(coords) {}
536
537  Value operator[](const Key& key) const {
538    return (_coords[_graph.v(key)] -
539            _coords[_graph.u(key)]).normSquare();
540  }
541
542private:
543
544  const Graph& _graph;
545  const CoordMap& _coords;
546};
547
548void minTree() {
549  std::vector<Parc> pedges;
550  Timer T;
551  std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
552  delaunay();
553  std::cout << T.realTime() << "s: Calculating spanning tree...\n";
554  LengthSquareMap<ListGraph, ListGraph::NodeMap<Point> > ls(g, coords);
555  ListGraph::EdgeMap<bool> tree(g);
556  kruskal(g, ls, tree);
557  std::cout << T.realTime() << "s: Removing non tree arcs...\n";
558  std::vector<Edge> remove;
559  for (EdgeIt e(g); e != INVALID; ++e) {
560    if (!tree[e]) remove.push_back(e);
561  }
562  for(int i = 0; i < int(remove.size()); ++i) {
563    g.erase(remove[i]);
564  }
565  std::cout << T.realTime() << "s: Done\n";
566}
567
568void tsp2()
569{
570  std::cout << "Find a tree..." << std::endl;
571
572  minTree();
573
574  std::cout << "Total arc length (tree) : " << totalLen() << std::endl;
575
576  std::cout << "Make it Euler..." << std::endl;
577
578  {
579    std::vector<Node> leafs;
580    for(NodeIt n(g);n!=INVALID;++n)
581      if(countIncEdges(g,n)%2==1) leafs.push_back(n);
582
583//    for(unsigned int i=0;i<leafs.size();i+=2)
584//       g.addArc(leafs[i],leafs[i+1]);
585
586    std::vector<Parc> pedges;
587    for(unsigned int i=0;i<leafs.size()-1;i++)
588      for(unsigned int j=i+1;j<leafs.size();j++)
589        {
590          Node n=leafs[i];
591          Node m=leafs[j];
592          Parc p;
593          p.a=n;
594          p.b=m;
595          p.len=(coords[m]-coords[n]).normSquare();
596          pedges.push_back(p);
597        }
598    std::sort(pedges.begin(),pedges.end(),pedgeLess);
599    for(unsigned int i=0;i<pedges.size();i++)
600      if(countIncEdges(g,pedges[i].a)%2 &&
601         countIncEdges(g,pedges[i].b)%2)
602        g.addEdge(pedges[i].a,pedges[i].b);
603  }
604
605  for(NodeIt n(g);n!=INVALID;++n)
606    if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 )
607      std::cout << "GEBASZ!!!" << std::endl;
608
609  for(EdgeIt e(g);e!=INVALID;++e)
610    if(g.u(e)==g.v(e))
611      std::cout << "LOOP GEBASZ!!!" << std::endl;
612
613  std::cout << "Number of arcs : " << countEdges(g) << std::endl;
614
615  std::cout << "Total arc length (euler) : " << totalLen() << std::endl;
616
617  ListGraph::EdgeMap<Arc> enext(g);
618  {
619    EulerIt<ListGraph> e(g);
620    Arc eo=e;
621    Arc ef=e;
622//     std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
623    for(++e;e!=INVALID;++e)
624      {
625//         std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
626        enext[eo]=e;
627        eo=e;
628      }
629    enext[eo]=ef;
630  }
631
632  std::cout << "Creating a tour from that..." << std::endl;
633
634  int nnum = countNodes(g);
635  int ednum = countEdges(g);
636
637  for(Arc p=enext[EdgeIt(g)];ednum>nnum;p=enext[p])
638    {
639//       std::cout << "Checking arc " << g.id(p) << std::endl;
640      Arc e=enext[p];
641      Arc f=enext[e];
642      Node n2=g.source(f);
643      Node n1=g.oppositeNode(n2,e);
644      Node n3=g.oppositeNode(n2,f);
645      if(countIncEdges(g,n2)>2)
646        {
647//           std::cout << "Remove an Arc" << std::endl;
648          Arc ff=enext[f];
649          g.erase(e);
650          g.erase(f);
651          if(n1!=n3)
652            {
653              Arc ne=g.direct(g.addEdge(n1,n3),n1);
654              enext[p]=ne;
655              enext[ne]=ff;
656              ednum--;
657            }
658          else {
659            enext[p]=ff;
660            ednum-=2;
661          }
662        }
663    }
664
665  std::cout << "Total arc length (tour) : " << totalLen() << std::endl;
666
667  std::cout << "2-opt the tour..." << std::endl;
668
669  tsp_improve();
670
671  std::cout << "Total arc length (2-opt tour) : " << totalLen() << std::endl;
672}
673
674
675int main(int argc,const char **argv)
676{
677  ArgParser ap(argc,argv);
678
679//   bool eps;
680  bool disc_d, square_d, gauss_d;
681//   bool tsp_a,two_a,tree_a;
682  int num_of_cities=1;
683  double area=1;
684  N=100;
685//   girth=10;
686  std::string ndist("disc");
687  ap.refOption("n", "Number of nodes (default is 100)", N)
688    .intOption("g", "Girth parameter (default is 10)", 10)
689    .refOption("cities", "Number of cities (default is 1)", num_of_cities)
690    .refOption("area", "Full relative area of the cities (default is 1)", area)
691    .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",disc_d)
692    .optionGroup("dist", "disc")
693    .refOption("square", "Nodes are evenly distributed on a unit square", square_d)
694    .optionGroup("dist", "square")
695    .refOption("gauss",
696            "Nodes are located according to a two-dim gauss distribution",
697            gauss_d)
698    .optionGroup("dist", "gauss")
699//     .mandatoryGroup("dist")
700    .onlyOneGroup("dist")
701    .boolOption("eps", "Also generate .eps output (prefix.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).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.