# HG changeset patch
# User deba
# Date 1181044159 0
# Node ID 260ce674cc65787aaeb3c4374587c9f2b9fec557
# Parent  dd20d76eed133615ed3293d38103a145f33b8361
Delaunay triangulation
Faster geometric minimum spanning tree

diff -r dd20d76eed13 -r 260ce674cc65 tools/lgf-gen.cc
--- a/tools/lgf-gen.cc	Tue Jun 05 10:59:16 2007 +0000
+++ b/tools/lgf-gen.cc	Tue Jun 05 11:49:19 2007 +0000
@@ -29,7 +29,7 @@
 #include <lemon/euler.h>
 #include <cmath>
 #include <algorithm>
-#include <lemon/unionfind.h>
+#include <lemon/kruskal.h>
 #include <lemon/time_measure.h>
 
 using namespace lemon;
@@ -146,32 +146,279 @@
 
 std::vector<UEdge> edges;
 
-void triangle()
-{
+namespace _delaunay_bits {
+
+  struct Part {
+    int prev, curr, next;
+
+    Part(int p, int c, int n) : prev(p), curr(c), next(n) {} 
+  };
+
+  inline std::ostream& operator<<(std::ostream& os, const Part& part) {
+    os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')';
+    return os;
+  }
+
+  inline double circle_point(const Point& p, const Point& q, const Point& r) {
+    double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y);
+    if (a == 0) return std::numeric_limits<double>::quiet_NaN();
+
+    double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) +
+      (q.x * q.x + q.y * q.y) * (r.y - p.y) +
+      (r.x * r.x + r.y * r.y) * (p.y - q.y);
+
+    double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) +
+      (q.x * q.x + q.y * q.y) * (r.x - p.x) +
+      (r.x * r.x + r.y * r.y) * (p.x - q.x);
+
+    double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) +
+      (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) +
+      (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y);
+
+    return d / (2 * a) + sqrt((d * d + e * e) / (4 * a * a) + f / a);
+  }
+
+  inline bool circle_form(const Point& p, const Point& q, const Point& r) {
+    return rot90(q - p) * (r - q) < 0.0;
+  }
+
+  inline double intersection(const Point& p, const Point& q, double sx) {
+    const double epsilon = 1e-8;
+
+    if (p.x == q.x) return (p.y + q.y) / 2.0;
+
+    if (sx < p.x + epsilon) return p.y;
+    if (sx < q.x + epsilon) return q.y;
+    
+    double a = q.x - p.x;
+    double b = (q.x - sx) * p.y - (p.x - sx) * q.y;    
+    double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare();
+    return (b - sqrt(d)) / a;
+  }
+
+  struct YLess {
+
+
+    YLess(const std::vector<Point>& points, double& sweep) 
+      : _points(points), _sweep(sweep) {}
+
+    bool operator()(const Part& l, const Part& r) const {
+      const double epsilon = 1e-8;
+
+      //      std::cerr << l << " vs " << r << std::endl;
+      double lbx = l.prev != -1 ?
+	intersection(_points[l.prev], _points[l.curr], _sweep) :
+	- std::numeric_limits<double>::infinity();
+      double rbx = r.prev != -1 ?
+	intersection(_points[r.prev], _points[r.curr], _sweep) :
+	- std::numeric_limits<double>::infinity();
+      double lex = l.next != -1 ?
+	intersection(_points[l.curr], _points[l.next], _sweep) :
+	std::numeric_limits<double>::infinity();
+      double rex = r.next != -1 ?
+	intersection(_points[r.curr], _points[r.next], _sweep) :
+	std::numeric_limits<double>::infinity();
+
+      if (lbx > lex) std::swap(lbx, lex);
+      if (rbx > rex) std::swap(rbx, rex);
+
+      if (lex < epsilon + rex && lbx + epsilon < rex) return true;
+      if (rex < epsilon + lex && rbx + epsilon < lex) return false;
+      return lex < rex;
+    }
+    
+    const std::vector<Point>& _points;
+    double& _sweep;
+  };
+  
+  struct BeachIt;
+  
+  typedef std::multimap<double, BeachIt> SpikeHeap;
+
+  typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach;
+
+  struct BeachIt {
+    Beach::iterator it;
+
+    BeachIt(Beach::iterator iter) : it(iter) {}
+  };
+
+}
+
+inline void delaunay() {
   Counter cnt("Number of edges added: ");
-  std::vector<Pedge> pedges;
-  for(NodeIt n(g);n!=INVALID;++n) 
-    for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
+  
+  using namespace _delaunay_bits;
+
+  typedef _delaunay_bits::Part Part;
+  typedef std::vector<std::pair<double, int> > SiteHeap;
+
+
+  std::vector<Point> points;
+  std::vector<Node> nodes;
+
+  for (NodeIt it(g); it != INVALID; ++it) {
+    nodes.push_back(it);
+    points.push_back(coords[it]);
+  }
+
+  SiteHeap siteheap(points.size());
+
+  double sweep;
+
+
+  for (int i = 0; i < int(siteheap.size()); ++i) {
+    siteheap[i] = std::make_pair(points[i].x, i);
+  }
+  
+  std::sort(siteheap.begin(), siteheap.end());
+  sweep = siteheap.front().first;
+  
+  YLess yless(points, sweep);
+  Beach beach(yless);
+
+  SpikeHeap spikeheap;
+
+  std::set<std::pair<int, int> > edges;
+
+  beach.insert(std::make_pair(Part(-1, siteheap[0].second, -1), 
+			      spikeheap.end()));
+  int siteindex = 1;
+
+  while (siteindex < int(points.size()) || !spikeheap.empty()) {
+
+    SpikeHeap::iterator spit = spikeheap.begin();
+
+    if (siteindex < int(points.size()) && 
+	(spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) {
+      int site = siteheap[siteindex].second;
+      sweep = siteheap[siteindex].first;
+          
+      Beach::iterator bit = beach.upper_bound(Part(site, site, site));
+      
+      if (bit->second != spikeheap.end()) {
+	spikeheap.erase(bit->second);	
+      }
+
+      int prev = bit->first.prev;
+      int curr = bit->first.curr;
+      int next = bit->first.next;
+
+      beach.erase(bit);
+      
+      SpikeHeap::iterator pit = spikeheap.end();
+      if (prev != -1 && 
+	  circle_form(points[prev], points[curr], points[site])) {
+	double x = circle_point(points[prev], points[curr], points[site]);
+	pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
+	pit->second.it =
+	  beach.insert(std::make_pair(Part(prev, curr, site), pit));
+      } else {
+	beach.insert(std::make_pair(Part(prev, curr, site), pit));
+      }
+
+      beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end()));
+      
+      SpikeHeap::iterator nit = spikeheap.end();
+      if (next != -1 && 
+	  circle_form(points[site], points[curr],points[next])) {
+	double x = circle_point(points[site], points[curr], points[next]);
+	nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
+	nit->second.it =
+	  beach.insert(std::make_pair(Part(site, curr, next), nit));
+      } else {
+	beach.insert(std::make_pair(Part(site, curr, next), nit));
+      }
+      
+      ++siteindex;
+    } else {
+      sweep = spit->first;      
+
+      Beach::iterator bit = spit->second.it;
+
+      int prev = bit->first.prev;
+      int curr = bit->first.curr;
+      int next = bit->first.next;
+
       {
-	Pedge p;
-	p.a=n;
-	p.b=m;
-	p.len=(coords[m]-coords[n]).normSquare();
-	pedges.push_back(p);
+	std::pair<int, int> edge;
+
+	edge = prev < curr ? 
+	  std::make_pair(prev, curr) : std::make_pair(curr, prev);
+	
+	if (edges.find(edge) == edges.end()) {
+	  edges.insert(edge);
+	  g.addEdge(nodes[prev], nodes[curr]);
+	  ++cnt;
+	}
+
+	edge = curr < next ? 
+	  std::make_pair(curr, next) : std::make_pair(next, curr);
+	
+	if (edges.find(edge) == edges.end()) {
+	  edges.insert(edge);
+	  g.addEdge(nodes[curr], nodes[next]);
+	  ++cnt;
+	}
       }
-  std::sort(pedges.begin(),pedges.end(),pedgeLess);
-  for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
-    {
-      Line li(pi->a,pi->b);
-      UEdgeIt e(g);
-      for(;e!=INVALID && !cross(e,li);++e) ;
-      UEdge ne;
-      if(e==INVALID) {
-	ne=g.addEdge(pi->a,pi->b);
-	edges.push_back(ne);
-	cnt++;
+      
+      Beach::iterator pbit = bit; --pbit;
+      int ppv = pbit->first.prev;
+      Beach::iterator nbit = bit; ++nbit;
+      int nnt = nbit->first.next;
+
+      if (bit->second != spikeheap.end()) spikeheap.erase(bit->second);
+      if (pbit->second != spikeheap.end()) spikeheap.erase(pbit->second);
+      if (nbit->second != spikeheap.end()) spikeheap.erase(nbit->second);
+
+      beach.erase(nbit);
+      beach.erase(bit);
+      beach.erase(pbit);
+
+      SpikeHeap::iterator pit = spikeheap.end();
+      if (ppv != -1 && ppv != next && 
+	  circle_form(points[ppv], points[prev], points[next])) {
+	double x = circle_point(points[ppv], points[prev], points[next]);
+	if (x < sweep) x = sweep;
+	pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
+	pit->second.it =
+	  beach.insert(std::make_pair(Part(ppv, prev, next), pit));
+      } else {
+	beach.insert(std::make_pair(Part(ppv, prev, next), pit));
       }
+
+      SpikeHeap::iterator nit = spikeheap.end();
+      if (nnt != -1 && prev != nnt &&
+	  circle_form(points[prev], points[next], points[nnt])) {
+	double x = circle_point(points[prev], points[next], points[nnt]);
+	if (x < sweep) x = sweep;
+	nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
+	nit->second.it =
+	  beach.insert(std::make_pair(Part(prev, next, nnt), nit));
+      } else {
+	beach.insert(std::make_pair(Part(prev, next, nnt), nit));
+      }
+      
     }
+  }
+
+  for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {
+    int curr = it->first.curr;
+    int next = it->first.next;
+
+    if (next == -1) continue;
+
+    std::pair<int, int> edge;
+
+    edge = curr < next ? 
+      std::make_pair(curr, next) : std::make_pair(next, curr);
+    
+    if (edges.find(edge) == edges.end()) {
+      edges.insert(edge);
+      g.addEdge(nodes[curr], nodes[next]);
+      ++cnt;
+    }
+  }
 }
 
 void sparse(int d) 
@@ -245,44 +492,43 @@
     }
 }
 
+template <typename UGraph, typename CoordMap>
+class LengthSquareMap {
+public:
+  typedef typename UGraph::UEdge Key;
+  typedef typename CoordMap::Value::Value Value;
+
+  LengthSquareMap(const UGraph& ugraph, const CoordMap& coords)
+    : _ugraph(ugraph), _coords(coords) {}
+
+  Value operator[](const Key& key) const {
+    return (_coords[_ugraph.target(key)] -
+	    _coords[_ugraph.source(key)]).normSquare();
+  }
+
+private:
+
+  const UGraph& _ugraph;
+  const CoordMap& _coords;
+};
+
 void minTree() {
-  int en=0;
-  int pr=0;
   std::vector<Pedge> pedges;
   Timer T;
-  std::cout << T.realTime() << "s: Setting up the edges...\n";
-  for(NodeIt n(g);n!=INVALID;++n) 
-    {
-      for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
-	{
-	  Pedge p;
-	  p.a=n;
-	  p.b=m;
-	  p.len=(coords[m]-coords[n]).normSquare();
-	  pedges.push_back(p);
-	}
-      if(progress && en>=pr*double(N)/100) 
-	{
-	  std::cout << pr << "%  \r" << std::flush;
-	  pr++;
-	}
-    }
-  std::cout << T.realTime() << "s: Sorting the edges...\n";
-  std::sort(pedges.begin(),pedges.end(),pedgeLess);
-  std::cout << T.realTime() << "s: Creating the tree...\n";
-  ListUGraph::NodeMap<int> comp(g);
-  UnionFind<ListUGraph::NodeMap<int> > uf(comp);
-  for (NodeIt it(g); it != INVALID; ++it)
-    uf.insert(it);  
-  for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
-    {
-      if ( uf.join(pi->a,pi->b) ) {
-	g.addEdge(pi->a,pi->b);
-	en++;
-	if(en>=N-1)
-	  break;
-      }
-    }
+  std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
+  delaunay();
+  std::cout << T.realTime() << "s: Calculating spanning tree...\n";
+  LengthSquareMap<ListUGraph, ListUGraph::NodeMap<Point> > ls(g, coords);
+  ListUGraph::UEdgeMap<bool> tree(g);
+  kruskal(g, ls, tree);
+  std::cout << T.realTime() << "s: Removing non tree edges...\n";
+  std::vector<UEdge> remove;
+  for (UEdgeIt e(g); e != INVALID; ++e) {
+    if (!tree[e]) remove.push_back(e);
+  }
+  for(int i = 0; i < int(remove.size()); ++i) {
+    g.erase(remove[i]);
+  }
   std::cout << T.realTime() << "s: Done\n";
 }
 
@@ -408,9 +654,27 @@
     .optionGroup("alg","tsp")
     .boolOption("tsp2", "Create a TSP tour (tree based)")
     .optionGroup("alg","tsp2")
+    .boolOption("dela", "Delaunay triangulation graph")
+    .optionGroup("alg","dela")
     .onlyOneGroup("alg")
+    .boolOption("rand", "Use time seed for random number generator")
+    .optionGroup("rand", "rand")
+    .intOption("seed", "Random seed", -1)
+    .optionGroup("rand", "seed")
+    .onlyOneGroup("rand")
     .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
     .run();
+
+  if (ap["rand"]) {
+    int seed = time(0);
+    std::cout << "Random number seed: " << seed << std::endl;
+    rnd = Random(seed);
+  }
+  if (ap.given("seed")) {
+    int seed = ap["seed"];
+    std::cout << "Random number seed: " << seed << std::endl;
+    rnd = Random(seed);
+  }
   
   std::string prefix;
   switch(ap.files().size()) 
@@ -463,6 +727,10 @@
 	    std::sqrt(sizes[s]/sum_sizes);
 	}
     }
+
+//   for (ListUGraph::NodeIt n(g); n != INVALID; ++n) {
+//     std::cerr << coords[n] << std::endl;
+//   }
   
   if(ap["tsp"]) {
     tsp();
@@ -482,6 +750,9 @@
   else if(ap["tree"]) {
     minTree();
   }
+  else if(ap["dela"]) {
+    delaunay();
+  }
   
 
   std::cout << "Number of nodes    : " << countNodes(g) << std::endl;