/**
* Closest pair
* Input: N points (x,y)
* Output: A pair of points (x1,y1),(x2,y2) which are closest.
* Algorithm: Line sweep (see http://w...content-available-to-author-only...r.com/tc?module=Static&d1=tutorials&d2=lineSweep)
* 
* Coder: Erel Segal http://t...content-available-to-author-only...s.fm/rentabrain
* Created: 2010-12-15
*/

#include <iostream>
#include <algorithm>
#include <set>
#include <map>
#include <vector>
#include <time.h>
#include <math.h>
#include <assert.h>
#include <limits>
using namespace std;


template <class Iterator> ostream& print (ostream& out, Iterator iFrom, Iterator iTo) {
	for (; iFrom!=iTo; ++iFrom) out << *iFrom << " ";
	return (out << endl);
}

template <class T> ostream& operator<< (ostream& out, const vector<T>& v) {
	return print (out, v.begin(), v.end());
}

template <class T> ostream& operator<< (ostream& out, const set<T>& v) {
	return print (out, v.begin(), v.end());
}

const int INF = 999999999;

struct Point {
	int x, y;
	Point() {x=y=INF;}
	Point(int _x, int _y) { x=_x; y=_y; }
	friend istream& operator>> (istream& in, Point& r) { in >> r.x >> r.y; return in; }
	friend ostream& operator<< (ostream& out, const Point& r) { out << r.x << "," << r.y; return out; }
};

bool lefter(const Point& A, const Point& B) { return A.x < B.x; }
bool righter(const Point& A, const Point& B) {  return A.x > B.x; }
bool lower(const Point& A, const Point& B) {  return A.y < B.y; }
bool upper(const Point& A, const Point& B) {  return A.y > B.y; }
float euclidean_distance(const Point& A, const Point& B) {  return sqrt( (A.x-B.x)*(A.x-B.x) + (A.y-B.y)*(A.y-B.y) ); }


class PointCollection {
	vector<Point> myPoints;
	typedef vector<Point>::const_iterator MyPointsIterator;
	
public:
	PointCollection(int N=0): myPoints(N) {}
	
	void clear() { myPoints.clear(); }
	
	void add(const Point& p) {
		myPoints.push_back(p);
	}

	void closestPairBruteForce() const {
		Point p1, p2;
		float shortestDistance = INF;
		for (MyPointsIterator i=myPoints.begin(); i!=myPoints.end(); ++i) {
			for (MyPointsIterator j=myPoints.begin(); j!=i; ++j) {
				float currentDistance = euclidean_distance(*i, *j);
				if (currentDistance<shortestDistance) {
					shortestDistance = currentDistance;
					p1 = *i;
					p2 = *j;
				}
			}
		}
		cout << "Brute Force Closest pair: " << p1 << " - " << p2 << " (distance: " << shortestDistance << ")" << endl;
	}

	/// @see http://w...content-available-to-author-only...r.com/tc?module=Static&d1=tutorials&d2=lineSweep
	void closestPairLineSweep() const {
		Point p1, p2;
		float shortestDistance = INF;
		//cout << "update shortest distance to " << shortestDistance << endl;

		// A. Order myPoints by the x coordinate:
		vector<Point> myPointsOrderedByX(myPoints);
		::sort(myPointsOrderedByX.begin(), myPointsOrderedByX.end(), &lefter);
		//cout << "myPointsOrderedByX: " << myPointsOrderedByX;

		// B. Maintain an active set of points - candidates to be a closest-pair with the current point.
		//    The set is ordered by the y coordinate:
		typedef set<Point, bool(*)(const Point&,const Point&)> ActiveSetType;
		typedef ActiveSetType::iterator ActiveSetIterator;
		ActiveSetType myActiveSet (&lower);
		
		MyPointsIterator leftestPointInActiveSet = myPointsOrderedByX.begin();
		
		// C. Line sweep the points by X order:
		for (MyPointsIterator point=myPointsOrderedByX.begin(); point!=myPointsOrderedByX.end(); ++point) {

			// C-1. Scan the active set for points that may be closer to the current point than the shortest distance:
			Point lowestPoint (point->x, point->y - shortestDistance);
			Point highestPoint (point->x, point->y + shortestDistance);
			ActiveSetIterator lowest = myActiveSet.lower_bound(lowestPoint);
			if (lowest!=myActiveSet.end()) {
				ActiveSetIterator highest = myActiveSet.upper_bound(highestPoint);
				//cout << "scan active set between " << lowestPoint << "<=" << *lowest << " and " << *highest << "<" << highestPoint << endl;

				for (ActiveSetIterator other=lowest; other!=highest; ++other) {
					//cout << "  compare " << *point << " with " << *other << endl;
					float currentDistance = euclidean_distance(*point,*other);
					if (currentDistance<shortestDistance) {
						shortestDistance = currentDistance;
						//cout << "    update shortest distance to " << shortestDistance << endl;
						p2 = *other;
						p1 = *point;
					}
				}
			}

			// C-2. Update the active set - keep only the points whose distance with the current point is at most the shortest distance
			//cout << "insert " << *point << " to active set" << endl;
			myActiveSet.insert(*point);
			while (point->x - leftestPointInActiveSet->x > shortestDistance) {
				//cout << "  remove " << *leftestPointInActiveSet << " from active set" << endl;
				myActiveSet.erase(*leftestPointInActiveSet);
				++leftestPointInActiveSet;
			}
		}
		
		cout << "Line Sweep  Closest pair: " << p1 << " - " << p2 << " (distance: " << shortestDistance << ")" << endl;
	}
	
	friend ostream& operator<< (ostream& out, PointCollection points) {
		return (out << points.myPoints);
	}
	
	void fillRandomX(int N) {
		Point p;
		for (int i=0; i<N; ++i) {
			p.x = rand()%1000;
			p.y = 500;
			add(p);
		}
	}
	
	void fillRandomY(int N) {
		Point p;
		for (int i=0; i<N; ++i) {
			p.y = rand()%1000;
			p.x = 500;
			add(p);
		}
	}
	
	void fillRandomXY(int N) {
		Point p;
		for (int i=0; i<N; ++i) {
			p.y = rand()%10000;
			p.x = rand()%10000;
			add(p);
		}
	}
	
	void fillRandomXYWithCheckpoints(int N) {
		// In this testcase, the min distance should always be 0, as there
		//	are two identical points.
		Point checkpoint(rand()%1000, rand()%1000);
		int i1 = rand()%N;
		int i2 = rand()%N;
		Point p;
		for (int i=0; i<N; ++i) {
			if (i==i1 || i==i2) {
				add(checkpoint);
			} else {
				p.y = rand()%1000;
				p.x = rand()%1000;
				add(p);
			}
		}
	}
};


int main() {
	srand(time(NULL));
	PointCollection points;
	points.clear();
	points.fillRandomXY(1000);
	points.closestPairLineSweep();
	points.closestPairBruteForce();
}