COIN-OR::LEMON - Graph Library

source: lemon/tools/lgf-gen.cc

tip
Last change on this file was 1308:89e1877e335f, checked in by Alpar Juttner <alpar@…>, 4 years ago

Clang compatibility fix in lgf-gen.cc (#480)

File size: 22.6 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        delete bit->second->second;
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, new 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, new 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())
404        {
405          delete bit->second->second;
406          spikeheap.erase(bit->second);
407        }
408      if (pbit->second != spikeheap.end())
409        {
410          delete pbit->second->second;
411          spikeheap.erase(pbit->second);
412        }
413      if (nbit->second != spikeheap.end())
414        {
415          delete nbit->second->second;
416          spikeheap.erase(nbit->second);
417        }
418     
419      beach.erase(nbit);
420      beach.erase(bit);
421      beach.erase(pbit);
422
423      SpikeHeap::iterator pit = spikeheap.end();
424      if (ppv != -1 && ppv != next &&
425          circle_form(points[ppv], points[prev], points[next])) {
426        double x = circle_point(points[ppv], points[prev], points[next]);
427        if (x < sweep) x = sweep;
428        pit = spikeheap.insert(std::make_pair(x, new BeachIt(beach.end())));
429        pit->second->it =
430          beach.insert(std::make_pair(Part(ppv, prev, next), pit));
431      } else {
432        beach.insert(std::make_pair(Part(ppv, prev, next), pit));
433      }
434
435      SpikeHeap::iterator nit = spikeheap.end();
436      if (nnt != -1 && prev != nnt &&
437          circle_form(points[prev], points[next], points[nnt])) {
438        double x = circle_point(points[prev], points[next], points[nnt]);
439        if (x < sweep) x = sweep;
440        nit = spikeheap.insert(std::make_pair(x, new BeachIt(beach.end())));
441        nit->second->it =
442          beach.insert(std::make_pair(Part(prev, next, nnt), nit));
443      } else {
444        beach.insert(std::make_pair(Part(prev, next, nnt), nit));
445      }
446
447    }
448  }
449
450  for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {
451    int curr = it->first.curr;
452    int next = it->first.next;
453
454    if (next == -1) continue;
455
456    std::pair<int, int> arc;
457
458    arc = curr < next ?
459      std::make_pair(curr, next) : std::make_pair(next, curr);
460
461    if (arcs.find(arc) == arcs.end()) {
462      arcs.insert(arc);
463      g.addEdge(nodes[curr], nodes[next]);
464      ++cnt;
465    }
466  }
467}
468
469void sparse(int d)
470{
471  Counter cnt("Number of arcs removed: ");
472  Bfs<ListGraph> bfs(g);
473  for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
474      ei!=arcs.rend();++ei)
475    {
476      Node a=g.u(*ei);
477      Node b=g.v(*ei);
478      g.erase(*ei);
479      bfs.run(a,b);
480      if(bfs.predArc(b)==INVALID || bfs.dist(b)>d)
481        g.addEdge(a,b);
482      else cnt++;
483    }
484}
485
486void sparse2(int d)
487{
488  Counter cnt("Number of arcs removed: ");
489  for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
490      ei!=arcs.rend();++ei)
491    {
492      Node a=g.u(*ei);
493      Node b=g.v(*ei);
494      g.erase(*ei);
495      ConstMap<Arc,int> cegy(1);
496      Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy);
497      int k=sur.run(a,b,2);
498      if(k<2 || sur.totalLength()>d)
499        g.addEdge(a,b);
500      else cnt++;
501//       else std::cout << "Remove arc " << g.id(a) << "-" << g.id(b) << '\n';
502    }
503}
504
505void sparseTriangle(int d)
506{
507  Counter cnt("Number of arcs added: ");
508  std::vector<Parc> pedges;
509  for(NodeIt n(g);n!=INVALID;++n)
510    for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
511      {
512        Parc p;
513        p.a=n;
514        p.b=m;
515        p.len=(coords[m]-coords[n]).normSquare();
516        pedges.push_back(p);
517      }
518  std::sort(pedges.begin(),pedges.end(),pedgeLess);
519  for(std::vector<Parc>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
520    {
521      Line li(pi->a,pi->b);
522      EdgeIt e(g);
523      for(;e!=INVALID && !cross(e,li);++e) ;
524      Edge ne;
525      if(e==INVALID) {
526        ConstMap<Arc,int> cegy(1);
527        Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy);
528        int k=sur.run(pi->a,pi->b,2);
529        if(k<2 || sur.totalLength()>d)
530          {
531            ne=g.addEdge(pi->a,pi->b);
532            arcs.push_back(ne);
533            cnt++;
534          }
535      }
536    }
537}
538
539template <typename Graph, typename CoordMap>
540class LengthSquareMap {
541public:
542  typedef typename Graph::Edge Key;
543  typedef typename CoordMap::Value::Value Value;
544
545  LengthSquareMap(const Graph& graph, const CoordMap& coords)
546    : _graph(graph), _coords(coords) {}
547
548  Value operator[](const Key& key) const {
549    return (_coords[_graph.v(key)] -
550            _coords[_graph.u(key)]).normSquare();
551  }
552
553private:
554
555  const Graph& _graph;
556  const CoordMap& _coords;
557};
558
559void minTree() {
560  std::vector<Parc> pedges;
561  Timer T;
562  std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
563  delaunay();
564  std::cout << T.realTime() << "s: Calculating spanning tree...\n";
565  LengthSquareMap<ListGraph, ListGraph::NodeMap<Point> > ls(g, coords);
566  ListGraph::EdgeMap<bool> tree(g);
567  kruskal(g, ls, tree);
568  std::cout << T.realTime() << "s: Removing non tree arcs...\n";
569  std::vector<Edge> remove;
570  for (EdgeIt e(g); e != INVALID; ++e) {
571    if (!tree[e]) remove.push_back(e);
572  }
573  for(int i = 0; i < int(remove.size()); ++i) {
574    g.erase(remove[i]);
575  }
576  std::cout << T.realTime() << "s: Done\n";
577}
578
579void tsp2()
580{
581  std::cout << "Find a tree..." << std::endl;
582
583  minTree();
584
585  std::cout << "Total arc length (tree) : " << totalLen() << std::endl;
586
587  std::cout << "Make it Euler..." << std::endl;
588
589  {
590    std::vector<Node> leafs;
591    for(NodeIt n(g);n!=INVALID;++n)
592      if(countIncEdges(g,n)%2==1) leafs.push_back(n);
593
594//    for(unsigned int i=0;i<leafs.size();i+=2)
595//       g.addArc(leafs[i],leafs[i+1]);
596
597    std::vector<Parc> pedges;
598    for(unsigned int i=0;i<leafs.size()-1;i++)
599      for(unsigned int j=i+1;j<leafs.size();j++)
600        {
601          Node n=leafs[i];
602          Node m=leafs[j];
603          Parc p;
604          p.a=n;
605          p.b=m;
606          p.len=(coords[m]-coords[n]).normSquare();
607          pedges.push_back(p);
608        }
609    std::sort(pedges.begin(),pedges.end(),pedgeLess);
610    for(unsigned int i=0;i<pedges.size();i++)
611      if(countIncEdges(g,pedges[i].a)%2 &&
612         countIncEdges(g,pedges[i].b)%2)
613        g.addEdge(pedges[i].a,pedges[i].b);
614  }
615
616  for(NodeIt n(g);n!=INVALID;++n)
617    if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 )
618      std::cout << "GEBASZ!!!" << std::endl;
619
620  for(EdgeIt e(g);e!=INVALID;++e)
621    if(g.u(e)==g.v(e))
622      std::cout << "LOOP GEBASZ!!!" << std::endl;
623
624  std::cout << "Number of arcs : " << countEdges(g) << std::endl;
625
626  std::cout << "Total arc length (euler) : " << totalLen() << std::endl;
627
628  ListGraph::EdgeMap<Arc> enext(g);
629  {
630    EulerIt<ListGraph> e(g);
631    Arc eo=e;
632    Arc ef=e;
633//     std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
634    for(++e;e!=INVALID;++e)
635      {
636//         std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
637        enext[eo]=e;
638        eo=e;
639      }
640    enext[eo]=ef;
641  }
642
643  std::cout << "Creating a tour from that..." << std::endl;
644
645  int nnum = countNodes(g);
646  int ednum = countEdges(g);
647
648  for(Arc p=enext[EdgeIt(g)];ednum>nnum;p=enext[p])
649    {
650//       std::cout << "Checking arc " << g.id(p) << std::endl;
651      Arc e=enext[p];
652      Arc f=enext[e];
653      Node n2=g.source(f);
654      Node n1=g.oppositeNode(n2,e);
655      Node n3=g.oppositeNode(n2,f);
656      if(countIncEdges(g,n2)>2)
657        {
658//           std::cout << "Remove an Arc" << std::endl;
659          Arc ff=enext[f];
660          g.erase(e);
661          g.erase(f);
662          if(n1!=n3)
663            {
664              Arc ne=g.direct(g.addEdge(n1,n3),n1);
665              enext[p]=ne;
666              enext[ne]=ff;
667              ednum--;
668            }
669          else {
670            enext[p]=ff;
671            ednum-=2;
672          }
673        }
674    }
675
676  std::cout << "Total arc length (tour) : " << totalLen() << std::endl;
677
678  std::cout << "2-opt the tour..." << std::endl;
679
680  tsp_improve();
681
682  std::cout << "Total arc length (2-opt tour) : " << totalLen() << std::endl;
683}
684
685
686int main(int argc,const char **argv)
687{
688  ArgParser ap(argc,argv);
689
690//   bool eps;
691  bool disc_d, square_d, gauss_d;
692//   bool tsp_a,two_a,tree_a;
693  int num_of_cities=1;
694  double area=1;
695  N=100;
696//   girth=10;
697  std::string ndist("disc");
698  ap.refOption("n", "Number of nodes (default is 100)", N)
699    .intOption("g", "Girth parameter (default is 10)", 10)
700    .refOption("cities", "Number of cities (default is 1)", num_of_cities)
701    .refOption("area", "Full relative area of the cities (default is 1)", area)
702    .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",
703               disc_d)
704    .optionGroup("dist", "disc")
705    .refOption("square", "Nodes are evenly distributed on a unit square",
706               square_d)
707    .optionGroup("dist", "square")
708    .refOption("gauss", "Nodes are located according to a two-dim Gauss "
709               "distribution", gauss_d)
710    .optionGroup("dist", "gauss")
711    .onlyOneGroup("dist")
712    .boolOption("eps", "Also generate .eps output (<prefix>.eps)")
713    .boolOption("nonodes", "Draw only the edges in the generated .eps output")
714    .boolOption("dir", "Directed graph is generated (each edge is replaced by "
715                "two directed arcs)")
716    .boolOption("2con", "Create a two connected planar graph")
717    .optionGroup("alg","2con")
718    .boolOption("tree", "Create a min. cost spanning tree")
719    .optionGroup("alg","tree")
720    .boolOption("tsp", "Create a TSP tour")
721    .optionGroup("alg","tsp")
722    .boolOption("tsp2", "Create a TSP tour (tree based)")
723    .optionGroup("alg","tsp2")
724    .boolOption("dela", "Delaunay triangulation graph")
725    .optionGroup("alg","dela")
726    .onlyOneGroup("alg")
727    .boolOption("rand", "Use time seed for random number generator")
728    .optionGroup("rand", "rand")
729    .intOption("seed", "Random seed", -1)
730    .optionGroup("rand", "seed")
731    .onlyOneGroup("rand")
732    .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
733    .run();
734
735  if (ap["rand"]) {
736    int seed = int(time(0));
737    std::cout << "Random number seed: " << seed << std::endl;
738    rnd = Random(seed);
739  }
740  if (ap.given("seed")) {
741    int seed = ap["seed"];
742    std::cout << "Random number seed: " << seed << std::endl;
743    rnd = Random(seed);
744  }
745
746  std::string prefix;
747  switch(ap.files().size())
748    {
749    case 0:
750      prefix="lgf-gen-out";
751      break;
752    case 1:
753      prefix=ap.files()[0];
754      break;
755    default:
756      std::cerr << "\nAt most one prefix can be given\n\n";
757      exit(1);
758    }
759
760  double sum_sizes=0;
761  std::vector<double> sizes;
762  std::vector<double> cum_sizes;
763  for(int s=0;s<num_of_cities;s++)
764    {
765      //         sum_sizes+=rnd.exponential();
766      double d=rnd();
767      sum_sizes+=d;
768      sizes.push_back(d);
769      cum_sizes.push_back(sum_sizes);
770    }
771  int i=0;
772  for(int s=0;s<num_of_cities;s++)
773    {
774      Point center=(num_of_cities==1?Point(0,0):rnd.disc());
775      if(gauss_d)
776        for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
777          Node n=g.addNode();
778          nodes.push_back(n);
779          coords[n]=center+rnd.gauss2()*area*
780            std::sqrt(sizes[s]/sum_sizes);
781        }
782      else if(square_d)
783        for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
784          Node n=g.addNode();
785          nodes.push_back(n);
786          coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*
787            std::sqrt(sizes[s]/sum_sizes);
788        }
789      else if(disc_d || true)
790        for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
791          Node n=g.addNode();
792          nodes.push_back(n);
793          coords[n]=center+rnd.disc()*area*
794            std::sqrt(sizes[s]/sum_sizes);
795        }
796    }
797
798//   for (ListGraph::NodeIt n(g); n != INVALID; ++n) {
799//     std::cerr << coords[n] << std::endl;
800//   }
801
802  if(ap["tsp"]) {
803    tsp();
804    std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
805  }
806  if(ap["tsp2"]) {
807    tsp2();
808    std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
809  }
810  else if(ap["2con"]) {
811    std::cout << "Make triangles\n";
812    //   triangle();
813    sparseTriangle(ap["g"]);
814    std::cout << "Make it sparser\n";
815    sparse2(ap["g"]);
816  }
817  else if(ap["tree"]) {
818    minTree();
819  }
820  else if(ap["dela"]) {
821    delaunay();
822  }
823
824
825  std::cout << "Number of nodes    : " << countNodes(g) << std::endl;
826  std::cout << "Number of arcs    : " << countEdges(g) << std::endl;
827  double tlen=0;
828  for(EdgeIt e(g);e!=INVALID;++e)
829    tlen+=std::sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
830  std::cout << "Total arc length  : " << tlen << std::endl;
831
832  if(ap["eps"])
833    graphToEps(g,prefix+".eps").scaleToA4().
834      scale(600).nodeScale(.005).arcWidthScale(.001).preScale(false).
835      coords(coords).hideNodes(ap.given("nonodes")).run();
836
837  if(ap["dir"])
838    DigraphWriter<ListGraph>(g,prefix+".lgf").
839      nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
840      nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
841      run();
842  else GraphWriter<ListGraph>(g,prefix+".lgf").
843         nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
844         nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
845         run();
846}
847
Note: See TracBrowser for help on using the repository browser.