1.1 --- a/tools/lgf-gen.cc Tue Jun 05 10:59:16 2007 +0000
1.2 +++ b/tools/lgf-gen.cc Tue Jun 05 11:49:19 2007 +0000
1.3 @@ -29,7 +29,7 @@
1.4 #include <lemon/euler.h>
1.5 #include <cmath>
1.6 #include <algorithm>
1.7 -#include <lemon/unionfind.h>
1.8 +#include <lemon/kruskal.h>
1.9 #include <lemon/time_measure.h>
1.10
1.11 using namespace lemon;
1.12 @@ -146,32 +146,279 @@
1.13
1.14 std::vector<UEdge> edges;
1.15
1.16 -void triangle()
1.17 -{
1.18 +namespace _delaunay_bits {
1.19 +
1.20 + struct Part {
1.21 + int prev, curr, next;
1.22 +
1.23 + Part(int p, int c, int n) : prev(p), curr(c), next(n) {}
1.24 + };
1.25 +
1.26 + inline std::ostream& operator<<(std::ostream& os, const Part& part) {
1.27 + os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')';
1.28 + return os;
1.29 + }
1.30 +
1.31 + inline double circle_point(const Point& p, const Point& q, const Point& r) {
1.32 + double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y);
1.33 + if (a == 0) return std::numeric_limits<double>::quiet_NaN();
1.34 +
1.35 + double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) +
1.36 + (q.x * q.x + q.y * q.y) * (r.y - p.y) +
1.37 + (r.x * r.x + r.y * r.y) * (p.y - q.y);
1.38 +
1.39 + double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) +
1.40 + (q.x * q.x + q.y * q.y) * (r.x - p.x) +
1.41 + (r.x * r.x + r.y * r.y) * (p.x - q.x);
1.42 +
1.43 + double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) +
1.44 + (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) +
1.45 + (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y);
1.46 +
1.47 + return d / (2 * a) + sqrt((d * d + e * e) / (4 * a * a) + f / a);
1.48 + }
1.49 +
1.50 + inline bool circle_form(const Point& p, const Point& q, const Point& r) {
1.51 + return rot90(q - p) * (r - q) < 0.0;
1.52 + }
1.53 +
1.54 + inline double intersection(const Point& p, const Point& q, double sx) {
1.55 + const double epsilon = 1e-8;
1.56 +
1.57 + if (p.x == q.x) return (p.y + q.y) / 2.0;
1.58 +
1.59 + if (sx < p.x + epsilon) return p.y;
1.60 + if (sx < q.x + epsilon) return q.y;
1.61 +
1.62 + double a = q.x - p.x;
1.63 + double b = (q.x - sx) * p.y - (p.x - sx) * q.y;
1.64 + double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare();
1.65 + return (b - sqrt(d)) / a;
1.66 + }
1.67 +
1.68 + struct YLess {
1.69 +
1.70 +
1.71 + YLess(const std::vector<Point>& points, double& sweep)
1.72 + : _points(points), _sweep(sweep) {}
1.73 +
1.74 + bool operator()(const Part& l, const Part& r) const {
1.75 + const double epsilon = 1e-8;
1.76 +
1.77 + // std::cerr << l << " vs " << r << std::endl;
1.78 + double lbx = l.prev != -1 ?
1.79 + intersection(_points[l.prev], _points[l.curr], _sweep) :
1.80 + - std::numeric_limits<double>::infinity();
1.81 + double rbx = r.prev != -1 ?
1.82 + intersection(_points[r.prev], _points[r.curr], _sweep) :
1.83 + - std::numeric_limits<double>::infinity();
1.84 + double lex = l.next != -1 ?
1.85 + intersection(_points[l.curr], _points[l.next], _sweep) :
1.86 + std::numeric_limits<double>::infinity();
1.87 + double rex = r.next != -1 ?
1.88 + intersection(_points[r.curr], _points[r.next], _sweep) :
1.89 + std::numeric_limits<double>::infinity();
1.90 +
1.91 + if (lbx > lex) std::swap(lbx, lex);
1.92 + if (rbx > rex) std::swap(rbx, rex);
1.93 +
1.94 + if (lex < epsilon + rex && lbx + epsilon < rex) return true;
1.95 + if (rex < epsilon + lex && rbx + epsilon < lex) return false;
1.96 + return lex < rex;
1.97 + }
1.98 +
1.99 + const std::vector<Point>& _points;
1.100 + double& _sweep;
1.101 + };
1.102 +
1.103 + struct BeachIt;
1.104 +
1.105 + typedef std::multimap<double, BeachIt> SpikeHeap;
1.106 +
1.107 + typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach;
1.108 +
1.109 + struct BeachIt {
1.110 + Beach::iterator it;
1.111 +
1.112 + BeachIt(Beach::iterator iter) : it(iter) {}
1.113 + };
1.114 +
1.115 +}
1.116 +
1.117 +inline void delaunay() {
1.118 Counter cnt("Number of edges added: ");
1.119 - std::vector<Pedge> pedges;
1.120 - for(NodeIt n(g);n!=INVALID;++n)
1.121 - for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
1.122 +
1.123 + using namespace _delaunay_bits;
1.124 +
1.125 + typedef _delaunay_bits::Part Part;
1.126 + typedef std::vector<std::pair<double, int> > SiteHeap;
1.127 +
1.128 +
1.129 + std::vector<Point> points;
1.130 + std::vector<Node> nodes;
1.131 +
1.132 + for (NodeIt it(g); it != INVALID; ++it) {
1.133 + nodes.push_back(it);
1.134 + points.push_back(coords[it]);
1.135 + }
1.136 +
1.137 + SiteHeap siteheap(points.size());
1.138 +
1.139 + double sweep;
1.140 +
1.141 +
1.142 + for (int i = 0; i < int(siteheap.size()); ++i) {
1.143 + siteheap[i] = std::make_pair(points[i].x, i);
1.144 + }
1.145 +
1.146 + std::sort(siteheap.begin(), siteheap.end());
1.147 + sweep = siteheap.front().first;
1.148 +
1.149 + YLess yless(points, sweep);
1.150 + Beach beach(yless);
1.151 +
1.152 + SpikeHeap spikeheap;
1.153 +
1.154 + std::set<std::pair<int, int> > edges;
1.155 +
1.156 + beach.insert(std::make_pair(Part(-1, siteheap[0].second, -1),
1.157 + spikeheap.end()));
1.158 + int siteindex = 1;
1.159 +
1.160 + while (siteindex < int(points.size()) || !spikeheap.empty()) {
1.161 +
1.162 + SpikeHeap::iterator spit = spikeheap.begin();
1.163 +
1.164 + if (siteindex < int(points.size()) &&
1.165 + (spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) {
1.166 + int site = siteheap[siteindex].second;
1.167 + sweep = siteheap[siteindex].first;
1.168 +
1.169 + Beach::iterator bit = beach.upper_bound(Part(site, site, site));
1.170 +
1.171 + if (bit->second != spikeheap.end()) {
1.172 + spikeheap.erase(bit->second);
1.173 + }
1.174 +
1.175 + int prev = bit->first.prev;
1.176 + int curr = bit->first.curr;
1.177 + int next = bit->first.next;
1.178 +
1.179 + beach.erase(bit);
1.180 +
1.181 + SpikeHeap::iterator pit = spikeheap.end();
1.182 + if (prev != -1 &&
1.183 + circle_form(points[prev], points[curr], points[site])) {
1.184 + double x = circle_point(points[prev], points[curr], points[site]);
1.185 + pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
1.186 + pit->second.it =
1.187 + beach.insert(std::make_pair(Part(prev, curr, site), pit));
1.188 + } else {
1.189 + beach.insert(std::make_pair(Part(prev, curr, site), pit));
1.190 + }
1.191 +
1.192 + beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end()));
1.193 +
1.194 + SpikeHeap::iterator nit = spikeheap.end();
1.195 + if (next != -1 &&
1.196 + circle_form(points[site], points[curr],points[next])) {
1.197 + double x = circle_point(points[site], points[curr], points[next]);
1.198 + nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
1.199 + nit->second.it =
1.200 + beach.insert(std::make_pair(Part(site, curr, next), nit));
1.201 + } else {
1.202 + beach.insert(std::make_pair(Part(site, curr, next), nit));
1.203 + }
1.204 +
1.205 + ++siteindex;
1.206 + } else {
1.207 + sweep = spit->first;
1.208 +
1.209 + Beach::iterator bit = spit->second.it;
1.210 +
1.211 + int prev = bit->first.prev;
1.212 + int curr = bit->first.curr;
1.213 + int next = bit->first.next;
1.214 +
1.215 {
1.216 - Pedge p;
1.217 - p.a=n;
1.218 - p.b=m;
1.219 - p.len=(coords[m]-coords[n]).normSquare();
1.220 - pedges.push_back(p);
1.221 + std::pair<int, int> edge;
1.222 +
1.223 + edge = prev < curr ?
1.224 + std::make_pair(prev, curr) : std::make_pair(curr, prev);
1.225 +
1.226 + if (edges.find(edge) == edges.end()) {
1.227 + edges.insert(edge);
1.228 + g.addEdge(nodes[prev], nodes[curr]);
1.229 + ++cnt;
1.230 + }
1.231 +
1.232 + edge = curr < next ?
1.233 + std::make_pair(curr, next) : std::make_pair(next, curr);
1.234 +
1.235 + if (edges.find(edge) == edges.end()) {
1.236 + edges.insert(edge);
1.237 + g.addEdge(nodes[curr], nodes[next]);
1.238 + ++cnt;
1.239 + }
1.240 }
1.241 - std::sort(pedges.begin(),pedges.end(),pedgeLess);
1.242 - for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
1.243 - {
1.244 - Line li(pi->a,pi->b);
1.245 - UEdgeIt e(g);
1.246 - for(;e!=INVALID && !cross(e,li);++e) ;
1.247 - UEdge ne;
1.248 - if(e==INVALID) {
1.249 - ne=g.addEdge(pi->a,pi->b);
1.250 - edges.push_back(ne);
1.251 - cnt++;
1.252 +
1.253 + Beach::iterator pbit = bit; --pbit;
1.254 + int ppv = pbit->first.prev;
1.255 + Beach::iterator nbit = bit; ++nbit;
1.256 + int nnt = nbit->first.next;
1.257 +
1.258 + if (bit->second != spikeheap.end()) spikeheap.erase(bit->second);
1.259 + if (pbit->second != spikeheap.end()) spikeheap.erase(pbit->second);
1.260 + if (nbit->second != spikeheap.end()) spikeheap.erase(nbit->second);
1.261 +
1.262 + beach.erase(nbit);
1.263 + beach.erase(bit);
1.264 + beach.erase(pbit);
1.265 +
1.266 + SpikeHeap::iterator pit = spikeheap.end();
1.267 + if (ppv != -1 && ppv != next &&
1.268 + circle_form(points[ppv], points[prev], points[next])) {
1.269 + double x = circle_point(points[ppv], points[prev], points[next]);
1.270 + if (x < sweep) x = sweep;
1.271 + pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
1.272 + pit->second.it =
1.273 + beach.insert(std::make_pair(Part(ppv, prev, next), pit));
1.274 + } else {
1.275 + beach.insert(std::make_pair(Part(ppv, prev, next), pit));
1.276 }
1.277 +
1.278 + SpikeHeap::iterator nit = spikeheap.end();
1.279 + if (nnt != -1 && prev != nnt &&
1.280 + circle_form(points[prev], points[next], points[nnt])) {
1.281 + double x = circle_point(points[prev], points[next], points[nnt]);
1.282 + if (x < sweep) x = sweep;
1.283 + nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
1.284 + nit->second.it =
1.285 + beach.insert(std::make_pair(Part(prev, next, nnt), nit));
1.286 + } else {
1.287 + beach.insert(std::make_pair(Part(prev, next, nnt), nit));
1.288 + }
1.289 +
1.290 }
1.291 + }
1.292 +
1.293 + for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {
1.294 + int curr = it->first.curr;
1.295 + int next = it->first.next;
1.296 +
1.297 + if (next == -1) continue;
1.298 +
1.299 + std::pair<int, int> edge;
1.300 +
1.301 + edge = curr < next ?
1.302 + std::make_pair(curr, next) : std::make_pair(next, curr);
1.303 +
1.304 + if (edges.find(edge) == edges.end()) {
1.305 + edges.insert(edge);
1.306 + g.addEdge(nodes[curr], nodes[next]);
1.307 + ++cnt;
1.308 + }
1.309 + }
1.310 }
1.311
1.312 void sparse(int d)
1.313 @@ -245,44 +492,43 @@
1.314 }
1.315 }
1.316
1.317 +template <typename UGraph, typename CoordMap>
1.318 +class LengthSquareMap {
1.319 +public:
1.320 + typedef typename UGraph::UEdge Key;
1.321 + typedef typename CoordMap::Value::Value Value;
1.322 +
1.323 + LengthSquareMap(const UGraph& ugraph, const CoordMap& coords)
1.324 + : _ugraph(ugraph), _coords(coords) {}
1.325 +
1.326 + Value operator[](const Key& key) const {
1.327 + return (_coords[_ugraph.target(key)] -
1.328 + _coords[_ugraph.source(key)]).normSquare();
1.329 + }
1.330 +
1.331 +private:
1.332 +
1.333 + const UGraph& _ugraph;
1.334 + const CoordMap& _coords;
1.335 +};
1.336 +
1.337 void minTree() {
1.338 - int en=0;
1.339 - int pr=0;
1.340 std::vector<Pedge> pedges;
1.341 Timer T;
1.342 - std::cout << T.realTime() << "s: Setting up the edges...\n";
1.343 - for(NodeIt n(g);n!=INVALID;++n)
1.344 - {
1.345 - for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
1.346 - {
1.347 - Pedge p;
1.348 - p.a=n;
1.349 - p.b=m;
1.350 - p.len=(coords[m]-coords[n]).normSquare();
1.351 - pedges.push_back(p);
1.352 - }
1.353 - if(progress && en>=pr*double(N)/100)
1.354 - {
1.355 - std::cout << pr << "% \r" << std::flush;
1.356 - pr++;
1.357 - }
1.358 - }
1.359 - std::cout << T.realTime() << "s: Sorting the edges...\n";
1.360 - std::sort(pedges.begin(),pedges.end(),pedgeLess);
1.361 - std::cout << T.realTime() << "s: Creating the tree...\n";
1.362 - ListUGraph::NodeMap<int> comp(g);
1.363 - UnionFind<ListUGraph::NodeMap<int> > uf(comp);
1.364 - for (NodeIt it(g); it != INVALID; ++it)
1.365 - uf.insert(it);
1.366 - for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
1.367 - {
1.368 - if ( uf.join(pi->a,pi->b) ) {
1.369 - g.addEdge(pi->a,pi->b);
1.370 - en++;
1.371 - if(en>=N-1)
1.372 - break;
1.373 - }
1.374 - }
1.375 + std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
1.376 + delaunay();
1.377 + std::cout << T.realTime() << "s: Calculating spanning tree...\n";
1.378 + LengthSquareMap<ListUGraph, ListUGraph::NodeMap<Point> > ls(g, coords);
1.379 + ListUGraph::UEdgeMap<bool> tree(g);
1.380 + kruskal(g, ls, tree);
1.381 + std::cout << T.realTime() << "s: Removing non tree edges...\n";
1.382 + std::vector<UEdge> remove;
1.383 + for (UEdgeIt e(g); e != INVALID; ++e) {
1.384 + if (!tree[e]) remove.push_back(e);
1.385 + }
1.386 + for(int i = 0; i < int(remove.size()); ++i) {
1.387 + g.erase(remove[i]);
1.388 + }
1.389 std::cout << T.realTime() << "s: Done\n";
1.390 }
1.391
1.392 @@ -408,9 +654,27 @@
1.393 .optionGroup("alg","tsp")
1.394 .boolOption("tsp2", "Create a TSP tour (tree based)")
1.395 .optionGroup("alg","tsp2")
1.396 + .boolOption("dela", "Delaunay triangulation graph")
1.397 + .optionGroup("alg","dela")
1.398 .onlyOneGroup("alg")
1.399 + .boolOption("rand", "Use time seed for random number generator")
1.400 + .optionGroup("rand", "rand")
1.401 + .intOption("seed", "Random seed", -1)
1.402 + .optionGroup("rand", "seed")
1.403 + .onlyOneGroup("rand")
1.404 .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
1.405 .run();
1.406 +
1.407 + if (ap["rand"]) {
1.408 + int seed = time(0);
1.409 + std::cout << "Random number seed: " << seed << std::endl;
1.410 + rnd = Random(seed);
1.411 + }
1.412 + if (ap.given("seed")) {
1.413 + int seed = ap["seed"];
1.414 + std::cout << "Random number seed: " << seed << std::endl;
1.415 + rnd = Random(seed);
1.416 + }
1.417
1.418 std::string prefix;
1.419 switch(ap.files().size())
1.420 @@ -463,6 +727,10 @@
1.421 std::sqrt(sizes[s]/sum_sizes);
1.422 }
1.423 }
1.424 +
1.425 +// for (ListUGraph::NodeIt n(g); n != INVALID; ++n) {
1.426 +// std::cerr << coords[n] << std::endl;
1.427 +// }
1.428
1.429 if(ap["tsp"]) {
1.430 tsp();
1.431 @@ -482,6 +750,9 @@
1.432 else if(ap["tree"]) {
1.433 minTree();
1.434 }
1.435 + else if(ap["dela"]) {
1.436 + delaunay();
1.437 + }
1.438
1.439
1.440 std::cout << "Number of nodes : " << countNodes(g) << std::endl;