COIN-OR::LEMON - Graph Library

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

Last change on this file since 570:ab6da8cf5ab2 was 570:ab6da8cf5ab2, checked in by Akos Ladanyi <ladanyi@…>, 10 years ago

Fix compilation with MSVC (#259)

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