分治法-最近点对问题

一．实验目的
（1）掌握分治法思想。
（2）学会最近点对问题求解方法。
二．实验步骤与结果
实验总体思路：
本实验求解问题为：随机生成点对，通过几种算法，获得最近点对距离，同时输出算法所花费得时间，进行算法之间的比较.
利用同一份数据进行多次拷贝，使用vector des(vector src)的方式，实现在同一种数据下，进行各种算法之间的比较。通过绘制各种算法消耗时间随数据量变化的曲线图，直观感受不同算法在不同时间复杂度下消耗的时间差异。

通过以下四种算法进行最近点对问题的求解：
[if !supportLists]1. [endif]蛮力法求解，通过计算所有点两点间距离来获得最近点对。
[if !supportLists]2. [endif]在蛮力法求解的基础上，使用仿函数重载运算符提高计算速度和代码的简洁性。
[if !supportLists]3. [endif]使用分治法的思想进行最近点对问题的求解
[if !supportLists]4. [endif]使用剪枝优化的方法进行最近点对问题的求解
实验中，使用pair<double,double>进行点的x，y数值存储，使用vector<pair<double,double>>进行点集存储。
[if !supportLists]l [endif]蛮力法求解，通过计算所有点两点间距离来获得最近点对。
def Original_find (s, left, right)
for (int i = left;i < right;i++)
for (int j = i + 1;j < right;j++)
double cur_distance = s[i] - s[j];
if (min_distance > cur_distance)
min_distance = cur_distance;
Return min_distance;

但蛮力法需要时间复杂度O(n^2),需要进行每两个点之间进行一次distance距离的求解。时间复杂度较高，故而我们需要对原始算法进行优化。
首先在需要不断调用的distance函数上进行优化，使用operator()-,进行仿函数的设计可以加快算法的时间。
[if !supportLists]l [endif]在蛮力法求解的基础上，使用仿函数重载运算符提高计算速度和代码的简洁性。
通过代码：
double operator
-(const pair<double, double>&
a, const pair<double, double>&
b) {
return sqrt((a.first

• b.first) * (a.first
• b.first) + (a.second
• b.second) * (a.second
• b.second));
  }

得到数据差别：
[if !vml]
[endif]
在数据量较大时，算法消耗的时间才开始有差异，但是在代码的结构上，更方便我们进行编码和操作，也符合其物理含义，通过减号获得两者之间的差距。
但我们发现，算法的时间复杂度还是非常高，故而考虑其他算法思想进行问题求解。

[if !supportLists]l [endif]使用分治法的思想进行最近点对问题的求解
分治法解最近点对问题算法流程图：

分治法求解问题的根本思想在于保证求解思路和求解过程相同的情况下，将大问题通过区域划分的方式分解成小问题进行求解，再将小问题的解返回回来筛选出我们需要的目标解。而分治法思想在最近点对问题的求解中则是使用了将大区域分割为小区域，再从小区域之中进行求解问题。
而考虑到需要进行区域的划分，我们使用sort对点集进行以first即(x)的升序排序，再进行划分。而sort的原因在于，经过sort后的数据，点集分布的下标与现实中的位置相符合。而且在后续的操作中，需要进行点集某个区域内的搜索，故而需要对原始数据进行以first的sort排序。
之后进行区域划分，分治求解，而在进行merge时，由于我们没有搜索中线部分左右两端的数据，故而需要求解他们其中最近的距离。
在一开始使用了中线左方对中线右方的每一个点都进行距离求解，算法复杂度为O(n^{2),故而总的算法复杂度为O(n}2)不符合我们的原意。所以我们需要对数据进行分析，我们可以发现，数据点之间具有互斥性，故而我们不需要对右方所有点进行测试，图（1.4）。同时我们发现，我们的数据在X轴上只分布在离中线distance的距离下，故而我们可以对这组数据进行重新关于Y轴排序，根据Y轴进行附近四个点的搜寻(鸽舍原理，点阵稀疏性)。但是算法效果较差
在一开始使用了sort函数以Y进行排序，则merge操作的时间复杂度为O(logn),所以总算法的时间复杂度为O(nloglogn)图(1.5)为merge操作所花时间。

而根据分治法的思想，我们考虑到merge排序算法也是通过分治法实现，故而可以将merge排序算法嵌入其中，将复杂度降低为O(nlogn)。
伪代码：
def Binary_divide_find(s,left,right):
if
right - left == 1 : val <-
INT_MAX; ret val
if
right - left == 2: val <- distance(left,left+1),merge(left,right),ret val;
mid = (right-left)/2 + left;
left_min_distance <- Binary_divide_find(s,left,mid)
right_min_distance <-Binary_divide_find(s, mid, right)
mid_point = s[mid];
cur_min_distance <- min(left_min_distance, right_min_distance);
merge(s, left, right);
search_arr(right - left) // distance(mid,Elem) < cur_min_distance
for (int i = 0;i < index;i++) :
for (int j = i + 1;j < i + 7 && j < index && abs(search_arr[i].second - search_arr[j].second) < cur_min_distance;j++) :
tmp_min_distance = min(tmp_min_distance, search_arr[i] - search_arr[j]);
Return
min(cur_min_distance,tmp_min_distance)
而根据鸽舍原理与实验结果(图1.6-1.8)我们可以发现，问题求解结果有较大概率在临近点中出现，故而我们考虑是否在前部分所说的暴力法中也不需要进行在X轴上与该点距离大于distance的点集呢(三角形定理).
]使用剪枝优化的方法进行最近点对问题的求解、
剪枝,就是通过某种判断,避免一些不必要的遍历过程,剪枝优化有三个原则: 正确、准确、高效。 而在本问题求解中，我们之所以能进行剪枝优化，是因为根据鸽舍原理，我们只需要进行其附近常数个点的搜索。故而通过该定理我们可以提前结束每一个点与后续远于distance的点的搜索。

我们只需要目标点上圈出一个小矩形进行搜索便可以的出该点与那个点距离最近：
算法伪代码：
def sort_find(s, left, right)
sort(s);
for (int i = left;i < right;i++) {
for (int j = i + 1;j < right;j++) {
double cur_distance = s[i] - s[j];
if (min_distance > cur_distance)
min_distance = cur_distance;
if (s[j].first - s[i].first > min_distance) break;
if (s[j].first - s[i].first == min_distance && s[j].second- s[i].second > min_distance) break;
ret min_distance;

四．实验心得

首先，通过本次实验，对最近点对距离有了比较好的思考过程。
其次，深刻意识到，好的算法的求解过程，需要一步一步进行优化，还需要有很多理论方面的论据证明。同时，着重感受到，同种算法思想下的两种算法可以相互嵌套，减少所需要的时间复杂度。
最后，感受到我们需要对问题的特殊性和数据的特殊性有充分的思考，而往往由于问题的特殊性和数据的特殊性我们可以从中提取到更多有用的信息来编写我们的算法。
下面为完整代码

#include <iostream>
#include <vector>
#include <algorithm>
#include <ctime>
#include <iomanip>
#include<windows.h>
#include<random>
using namespace std;
pair<double, double> myrand();
double get_distance(const pair<double, double>& a, const pair<double, double>& b);
double operator -(const pair<double, double>& a, const pair<double, double>& b);
void Create_s(vector<pair<double, double>>& s, int k);
void display(vector<pair<double, double>>& s);
double sort_find(vector<pair<double, double>>& s, int left, int right);
double Binary_divide_find(vector<pair<double, double>>& s, int left, int right);
double Intersection_search(vector<pair<double, double>>& search_arr, int len, double cur_min_distance);
double Original_find_plus(vector<pair<double, double>>& s, int left, int right);
double Original_find(vector<pair<double, double>>& s, int left, int right);
void merge(vector<pair<double, double>>& s, int left, int right);
int main() {
	int k = 10;
	int k1 = 1;
	int size = 10000;
	while (k1 < k) {
		vector<pair<double, double>> s;
		vector<pair<double, double>> s1(size);
		vector<pair<double, double>> s2(size);
		Create_s(s, size);
		copy(s.begin(), s.end(), s1.begin());
		copy(s.begin(), s.end(), s2.begin());
		clock_t  start, end;
	/*	cout << "Search min distance from point in size = " << size << endl;
		start = clock();
		cout << "Oringnal_find result " << Original_find(s, 0, size) << endl;
		end = clock();
		cout<< "Oringnal_find cost "<< end - start << "ms" << endl;
		start = clock();
		cout << "Oringnal_find_plus result " << Original_find_plus(s, 0, size) << endl;
		end = clock();
		cout << "Oringnal_find_plus cost " << end - start << "ms" << endl;*/
		start = clock();
		cout << "sort_find result " << sort_find(s1, 0, size) << endl;
		end = clock();
		cout << "sort_find cost " << end - start << "ms" << endl;
		start = clock();
		//sort(s2.begin(), s2.end());
		cout << "Binary_divide_find result " << Binary_divide_find(s2, 0, size) << endl;
		end = clock();
		cout << "Binary_divide_find cost " << end - start << "ms" << endl;
		k1++;			
	}

}
pair<double, double> myrand() {
	random_device rd;
	mt19937 gen(rd());
	uniform_real_distribution<> dis(0, 10000.0);
	return { dis(gen),dis(gen) };

}

double get_distance(const pair<double, double>& a, const pair<double, double>& b)
{
	return sqrt((a.first - b.first) * (a.first - b.first) + (a.second - b.second) * (a.second - b.second));
}

double operator -(const pair<double, double>& a, const pair<double, double>& b) {
	return sqrt((a.first - b.first) * (a.first - b.first) + (a.second - b.second) * (a.second - b.second));
}

void Create_s(vector<pair<double, double>>& s, int k)
{
	Sleep(1000);
	srand(time(0));
	s.resize(k);
	generate(s.begin(), s.end(), myrand);
}

void display(vector<pair<double, double>>& s) {
	for (const auto& point : s) {
		cout << "( " << point.first << " , " << point.second << ") " << endl;
	}
}
double sort_find(vector<pair<double, double>>& s, int left, int right)
{
	sort(s.begin(), s.end());
	double min_distance = INT_MAX;
	pair<int, int> min_index;
	for (int i = left;i < right;i++) {
		for (int j = i + 1;j < right;j++) {
			double cur_distance = s[i] - s[j];
			if (min_distance > cur_distance) {
				min_index = { i,j };
				min_distance = cur_distance;
			}
			if (s[j].first - s[i].first > sqrt(min_distance)) break;
		}
	}
	return min_distance;
}

double Binary_divide_find(vector<pair<double, double>>& s, int left, int right)
{
	if (right - left == 1) {
		return INT_MAX;
	}
	else if (right - left == 2) {
		if (s[left].second > s[left + 1].second) swap(s[left], s[left + 1]);
		return s[right - 1] - s[left];
	}
	int mid = (right - left) / 2 + left;
	pair<double, double>  mid_point = s[mid];
	double left_min_distance = Binary_divide_find(s, left, mid);
	double right_min_distance = Binary_divide_find(s, mid, right);
	double cur_min_distance = min(left_min_distance, right_min_distance);
	merge(s, left, right);
	vector<pair<double, double>> search_arr(right - left);
	int index = 0;
	for (int i = left;i < right;i++) {
		if (fabs(s[i].first-mid_point.first) < cur_min_distance) search_arr[index++] = s[i];
	}
	for (int i = 0;i < index;i++) {
		for (int j = i + 1;j < i + 5 && j < index && abs(search_arr[i].second - search_arr[j].second) < cur_min_distance;j++) {
			cur_min_distance = min(cur_min_distance, search_arr[i] - search_arr[j]);
		}
	}
	
	return cur_min_distance;
}

double Intersection_search(vector<pair<double, double>>& search_arr, int len, double cur_min_distance)
{
	double min_distance = INT_MAX;
	for (int i = 0;i < len;i++) {
		for (int j = i + 1;j < i + 5 && j < len && abs(search_arr[i].second - search_arr[j].second) < cur_min_distance;j++) {
			min_distance = min(min_distance, search_arr[i] - search_arr[j]);
		}
	}
	return min_distance;
}
double Original_find_plus(vector<pair<double, double>>& s, int left, int right)
{
	time_t start = time(0);
	double min_distance = INT_MAX;
	pair<int, int> min_index;
	for (int i = left;i < right;i++) {
		for (int j = i + 1;j < right;j++) {
			double cur_distance = s[i] - s[j];
			if (min_distance > cur_distance) {
				min_index = { i,j };
				min_distance = cur_distance;
			}
		}
	}
	time_t end = time(0);
	return min_distance;
}
double Original_find(vector<pair<double, double>>& s, int left, int right)
{
	time_t start = time(0);
	double min_distance = INT_MAX;
	pair<int, int> min_index;
	for (int i = left;i < right;i++) {
		for (int j = i + 1;j < right;j++) {
			double cur_distance = get_distance(s[i], s[j]);
			if (min_distance > cur_distance) {
				min_index = { i,j };
				min_distance = cur_distance;
			}
		}
	}
	time_t end = time(0);
	return min_distance;
}
void merge(vector<pair<double, double>>& s, int left, int right) {
	if (right - left == 1) return;
	if (right - left == 2) {
		if (s[left].second > s[left + 1].second) swap(s[left], s[left + 1]);
		return;
	}
	vector<pair<double, double>> res(right - left);
	int mid = (right - left) / 2 + left;
	int left_cur = left;
	int right_cur = mid;
	int main_cur = 0;
	while (left_cur < mid && right_cur < right) {
		while (left_cur < mid && right_cur < right&&s[left_cur].second <= s[right_cur].second) res[main_cur++] = s[left_cur++];
		while (left_cur < mid && right_cur < right&&s[left_cur].second > s[right_cur].second) res[main_cur++] = s[right_cur++];
	}
	while (left_cur < mid)res[main_cur++] = s[left_cur++];
	while (right_cur < right) res[main_cur++] = s[right_cur++];
	copy(res.begin(), res.end(), s.begin() + left);
	//cout << right - left << endl;
}