算法与程序设计(实验2）----分治法求最近点对问题

一．实验目的

掌握分治法思想。
学会最近点对问题求解方法。

二、实验内容

1. 对于平面上给定的N个点，给出具有最短距离的两点。

2. 要求随机生成N个点的平面坐标，应用蛮力法编程计算出所有点对的最短距离。

3. 要求随机生成N个点的平面坐标，应用分治法编程计算出所有点对的最短距离。

4. 分别对N=100000—1000000，统计算法运行时间，比较理论效率与实测效率的差异，同时对蛮力法和分治法的算法效率进行分析和比较。

5. 算法执行过程可视化。

三、实验步骤与结果

1.实验

（1）问题描述

对于平面上给定的N个点，给出所有点对的最短距离。即输入是平面上的N个点，输出是N点中具有最短距离的两点。

生成随机点

Creat_new_point //点的坐标(结构体)

struct Point

double x, y;

要保证两点不重合且最终有序，经过比较分析，本实验采用Set容器（集合）。

时间复杂度上，若采用数组，查重的时间复杂度最优情况为O（nlogn），排序的最优情况时间复杂度为O（nlogn）。而set 中的insert(x)函数可将 x 插入 set 容器中，并自动递增排序和去重，时间复杂度为 O（logn），其中 n 为集合内的元素个数。set集合存储点的信息，创建点集的时间复杂度为O（log1+log2+log3+...+nlogn）。则当n趋近于无穷时，nlogn与log（n!）为同阶无穷大。则该过程时间复杂度为O（logn!）

Creat_Random_Point（n） //n个点

while(set.size()!=n)

new_point.x=rand();

new_point.y=rand();

set.insert(new_point);

小规模测试生成随机点坐标，并验证：有序且不重复

图1 生成随机点坐标测试

重载比较函数：

bool operator == (const Point &p2)

bool operator < (const Point &p2)

return x == p2.x && y == p2.y

if x == p2.x

return y < p2.y;

return x < p2.x;

2.算法

（1）暴力算法

设置变量dis存放最短距离，初始化为INF。变量P1,P2存储相应的两点，初始化为-1。

循环遍历每一个点，每当遍历到一个点p时，循环遍历其他的点r到p的距离，如果距离小于dis，则更新dis和P1、P2的值。

图2 暴力穷举法流程图

伪代码：

Violent_exhaustive(*P)

q1 = 0 and q2 = 0 and min = INF

//初始化：当前最小距离为INF无穷大，其两端点下标为0

for i = 0 to n-1 //遍历所有的点

for j = i+1 to n

if dis(S[i]，S[j])<min

min = dis(S[i]，S[j]) //更新最短距离点信息

q1 = i

q2 = j

时间复杂度：遍历两遍循环，则时间复杂度为O（n^2）

表一暴力算法不同n运行时间

运行时间（s）	实际值	理论值	误差%
100000	28.042	28.042	0
200000	116.524	112.168	-3.883460524
300000	259.923	252.378	-2.989563274
400000	457.987	448.672	-2.076126881
500000	731.232	701.05	-4.305256401
600000	1028.925	1009.512	-1.923008345
700000	1456.279	1374.058	-5.983808544
800000	1889.198	1794.688	-5.266096391
900000	2415.253	2271.402	-6.333136979
1000000	2958.634	2804.2	-5.507239141

图3 暴力运行时间理论值与实际值比较

对于暴力穷举，由于算法未进行任何简化且算法执行次数与时间复杂度不随数据分布与数量级有任何变化，故整体拟合效果非常好，证明暴力法的时间复杂度为O(n^2)的结论正确。

（2）分治法

①思想

分治法可以将问题缩小化，将庞大的问题逐渐缩小成单个小的问题进行解决。

②步骤

a.预处理：根据输入点集S中的x轴和y轴坐标进行排序。

b. 分子集：当点数|S|>3时，将平面点集S分割成为大小大致相等的两个子集SL和SR，选取一个垂直线L作为分割直线，如图下图所示。

图4 分治法图示

c.求区间最值：两个递归调用，分别求出SL和SR中的最短距离为dl和dr。如图2所示，取d=min(dl, dr)，在直线L两边分别扩展d，得到边界区域Y，Y’是区域Y中的点按照y坐标值排序后得到的点集，Y'又可分为左右两个集合Y’L和Y’R。

递归调用的临界判断条件为：①该区域仅剩一个点时，返回INF；②该区域剩两个点时，返回两点距离③该区域有三个点时，返回三对点对之间距离最小值。④否则，将区间分成两份，分别执行c操作（递归求区间最值）。

图5 分治法划分约束区域

d.对于Y’L中的每一点，检查Y’R中的点与它的距离，更新所获得的最近距离。而更新最近距离时，不需要使用暴力穷举法，对[mid-d,mid+d]中的排序后，不妨设最终情况为pl,pr两点分别位于点集PL,PR，且该点对间距离小于d，则pl,pr两点一定位于下图中2d*d的矩形内。这是因为根据分析，当且仅当点对位于该蓝色区域内时，两点间纵坐标之差小于d，任何其他不在该范围内的点横坐标或纵坐标之差必定大于d，距离必定大于d，此时，该点对间距离一定不为最小，故无需进行比较。所以每个点最多比较6次，当与往后第i个点的距离大于当前最短距离变量d时（i<6）,无需比较剩下6-i个点。

图6 约束区间

图7 分治法流程图

④伪代码

#结构体点的构造：

struct Point

int x,y;

bool operator <(const Point &p2)

if(x==p2.x)

return y<p2.y;

return x<p2.x;

bool operator ==(const Point &p2)

return x==p2.x&&y==p2.y;

#重载比较函数：

bool cmpy(const Point &p1,const Point &p2) 重载比较函数

if(p1.y==p2.y)

return p1.x<p2.x;

return p1.y<p2.y;

#返回距离值函数：

double dis(Point p1,Point p2) 返回距离的值

d=sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));

return d;

#递归求区域内两点的最小距离函数：

double get_min(int begin,int end)

if (begin >= end) 区域内只有一个点时：返回无穷大

return INF;

if (end - begin == 1) 区域内有两个点：返回两点距离

return dis(P[begin], P[end]);

if (end - begin == 2) 区域内有三个点：返回区域三条连线中最短的线的距离

d = min(dis(P[begin], P[begin + 1]), dis(P[begin], P[begin + 1]));

d = min(dis(P[end], P[begin + 1]), d);

return d;

mid = (begin + end) / 2; 中点坐标作为划分区域的直线

d = min(get_min(begin, mid), get_min(mid, end)); 更新最短距离

#求中间区域[mid-d,mid+d]的最小值：

cnt = 0;

Point *py=new Point[end-begin+1]; 设立新数组存放中间区域的点

for (int i = begin; i <= end; i++)

if (P[i].x >= P[mid].x-d)

if(P[i].x>P[mid].x+d)

break;

py[cnt++] = P[i];

sort(py, py + cnt, cmpy); 按y坐标排序

for (int i = 0; i < cnt; i++)

for (int j = i + 1; j < cnt; j++)

if (py[j].y - py[i].y >=d)

break;

d= min(d, dis(py[i], py[j])); 更新最短距离

return d; 返回最短距离

表二分治算法不同n运行时间

运行时间（ms）	实际值	理论值	误差%
100000	53.7	53.7	0
200000	109.9	113.8661243	3.483146837
300000	164	176.4728468	7.067856076
400000	235	240.6644972	2.353690425
500000	271.1	306.0346892	11.41527103
600000	349.6	372.3440666	6.108346719
700000	420	439.4344706	4.422609501
800000	493.8	507.1934917	2.640706535
900000	545.3	575.537081	5.253715523
1000000	588.5	644.4	8.674736189

图8 分治算法不同n的运行时间

对于分治法，由于分治法需要开辟空间并实现递归操作，递归栈的迭代与开辟内存空间均需要消耗一部分时间。对于数量级较小的数据时，整个算法实际运行时间比较短，开辟内存空间与递归栈迭代消耗的时间的影响将被一定程度上放大。但当数量级较大时，这种影响被冲淡。因此会体现出当数量级较小时拟合效果较差，但当数量级较大时，拟合效果比较好。

(3)算法性能

图9 两种时间复杂度比较

通过对算法时间的分析知道，当对于小数量级的数据时（10^3以下），暴力穷举法要优于分治法，但当数量级较大时，分治法明显优于暴力穷举法。这是由于分治法需要开辟内存空间并通过递归栈实现递归。因此将消耗更多时间，并对小数量级下的最终时间消耗产生较大影响。

因此对于“最近点对”问题最优的方法是当数量级小于10^3时选择暴力穷举法进行运算，当数量级大于10^3时选择分治法进行运算。

表3 对小规模数据下测试两种算法的性能

穷举法			分治法
数量级	平均时间	理论时间	数量级	平均时间	理论时间
10	0.00025	0.00025	10	0.0043	0.001078
100	0.02805	0.02805	100	0.03371	0.0215687
1000	2.79846	2.79846	1000	0.38379	0.3235305
10000	278.374	278.374	10000	4.31374	4.31374
100000	28042	28042	100000	52.4283	53.92175

三．实验心得

实验成功！在本次实验中我了解并掌握了分治法的思想、学会了最近点对问题求解的方法。在实验过程中，也遇到了许多问题，经过努力都得以解决。值得一提的是，在测试的时候发现，当测试规模相同时，每次测试生成的随机数相同，查阅资料得知rand()方法生成的是伪随机数，是根据一个数按照某个公式推算出来的，这个数我们称之为“种子”，但是这个种子在系统启动之后就是一个定值，所以造成了这种问题。后来想到了解决办法：rand() 产生随机数时，用srand(seed) 用时间作种子 srand(time(NULL)) 播下后，因为每次运行程序的时间肯定是不相同的，所以产生的随机数必然不同。time()返回值是以秒为单位，故在每一次测试后休眠一秒钟

算法过程的可视化

为了对实验的算法过程有更清晰深入的了解，我使用Python中的matplotlib函数库，实现了小数据规模下算法过程的可视化，清晰的展示了排序的具体过程。请查看附件中的两个gif文件。