【深大算法设计与分析】实验二 分治法求最近点对问题 实验报告 附代码

目录

一、实验目的与要求

二、实验内容与方法

三、实验步骤与过程

四、实验结论与体会

尾注：

---

一、实验目的与要求

实验目的：

        1. 掌握分治法思想。

        2. 学会最近点对问题求解方法。

实验要求：

        1. 在blackboard提交电子版实验报告。

        2. 详细给出算法思想与实现代码之间的关系解释，不可直接粘贴代码。

        3. 作为实验报告附件上传源代码。

        4. 在实验课现场运行代码验证。

二、实验内容与方法

        1. 在平面上随机生成N个点的坐标，分别应用蛮力法和分治法编程计算出所有点对的最短距离。

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

三、实验步骤与过程

1. 准备工作，生成数据集和并且去重：

       首先定义结构体：“点”，包含两个整形数据x和y，分别表示横纵坐标。接着创建长度为N的“点”的数组，并用rand()函数给各个点赋值。

       考虑实际意义，如果数据集中两个点完全重合（距离为0），说明直接找到了距离最小的点对，接下来的运算将失去意义。因此需要去重。

        我们使用的去重方法是：先定义一个unordered_set容器（此容器使用的是哈希表存储数据，因此无论是查找还是添加元素都可以在O(1)的时间复杂度下完成）保存点，接下来遍历所有的点，每遍历一个点，就判定这个点是否以及在unordered_set中，如果在，则说明重复，重新生成这个点，再次判断，直到不重复为止；如果不在，就直接将其加入unordered_set。

2. 算法描述

（1）暴力算法：

        暴力算法较为简单，使用两层循环，先遍历每一个点，对于遍历到的某一个点a，去计算它和其他所有点的距离（实际上由于a到b的距离等于b到a的距离，只需要计算编号在a之后的所有点到a的距离），统计最小值即可。

        由于需要两层循环，暴力算法的时间复杂度为O(\(n^2\))。

        本实验中的暴力算法经过了计算的优化：如果找到的新一个点对的横坐标之差或纵坐标之差（的绝对值）已经大于当前找到的最小距，就直接跳过，且计算距离的函数不再是计算根号（坐标之差）平方和再返回，而是只计算平方和，避免了开根的复杂运算，加快效率。

代码核心步骤如图，其中，length表示的是数据集的规模，ans是最终返回的结果，dis()函数返回的是两点之间的距离。

图：暴力算法核心代码

（2）分治算法：

        暴力算法计算了所有点对的距离，但是随着计算的运行，可以逐渐排除一些距离过远的点对，尽可能避免多余的运算。而避免多余运算的一种重要方法正是使用分治法：

        由于分治法会将原始问题递归划分为一个个规模更小的问题，所以首先我们可以对规模最小的问题（不可继续分解的问题）进行讨论：

        对于一个区域，如果此区域只有一个点，那么此区域就不存在点对，最小点对距离可以认为是无穷；如果此区域有两个点，那么最小点对距离就是这个点对的距离；而如果有三个点，就需要求出三个距离，取最小值。

图：最小规模问题示意（取自实验文档）

        接下来讨论可以继续分解的问题：如果某区域点数大于3，说明此区域可以继续划分为点数更少的区域，到这一步，会出现三个问题需要解决：

        ①如何尽可能快的将此区域划分为两块拥有相同点数（或相差1）的小区域。

        ②如何解决划分后得到的小规模问题。

        ③由于两个小规模问题只考虑了两边各自的最小距离，如何处理最小距离点对横跨两个区域的情况。

       对于问题①，可以在数据处理阶段对所有点的横坐标（或纵坐标）进行排序，接下来对于某块区域，只需要知道最左、最右端点的编号，就可以直接找到划分线。而这一步效率的关键就在于排序算法的选择，如果使用了效率为O(\(n^2\))的算法，就会直接导致此方法失去效率优势。本实验中我们使用的是快速排序，时间复杂度为O(nlogn)。

代码：简单划分排序后的点集

        对于问题②，由于已经提出了最小规模问题的解决方案，故可以由递归函数返回解决。

        而问题③就是本实验最为复杂、巧妙的部分：

        对于任意状态的分割，设分界线为mid、左边区域最小点对距离为dl，右边的为dr，取d=min(dl,dr)作为两边找到的最小距离。对于合并后区域的任意一个点a，只要a的横坐标不在mid-d到mid+d的范围内，则a一定是可以排除的点，因为不在这个范围内的任点到另一区域某点的距离一定是大于最小点对距离的！实际编写代码时，为了方便处理，使用的不是“排除”不合条件的点，而是“找到”符合条件的点的方法。核心代码实现如下：

代码：找出在此范围内的所有点

        接下来的处理即是本次实验的精髓。此时我们已经找出了符合范围的点，显然，暴力算法可以解决此问题：遍历某一边所有的点对于每一个点，找到其对另一边所有点的最小距离，最终求最小。但是这种方法并未达到线性效率，当此区域的点较多时，仍会占用大量运算时间。

        那么如何解决呢？首先我们要对范围内所有的点进行纵坐标排序，从下往上从小到大。从小到大遍历该区域的点，对于某一个点a，如果另一侧存在一个点b可以刷新最小点对距，那么b的位置必然出现在以y=ya,y=ya+d,x=mid-d,x=mid+d四条直线围成的矩形中（因为如果b出现在此矩形外，则ab间距显然大于d，矛盾）。

图：ab位置关系简要示意

       那么接下来只需要证明：此矩形内的点数存在一个常数上限（且远小于数据集大小），就可以保证此步骤的效率达到线性。

       而显然，此矩形内可放入的最多点数为6，分布如图。

图：矩形内点最多的情况

        这是因为矩形本身的大小由两个划分的最小点距决定，这导致了矩形的左边（或右边）的小正方形中每个点之间的距离均大于等于d，要装下最多的点只能是四个点分布于四个顶点的情况，将这两个小正方形拼起来，排除两个重合的点，最终只能装下六个点。故证明了此步骤的效率可以达到线性。

        核心代码如下：

代码：解决跨划分点对问题

        时间复杂度计算：

        横纵坐标的快速排序时间复杂度为O(nlogn)；

        对于找异边点的最小点对距离中，获取点的时间复杂度为O(n)，对点排序O(nlogn)，找点对距离是n循环嵌套一层常数循环，时间复杂度为O(n)，故总时间复杂度为O(nlogn)。

3. 运行效率分析：

       不同数据集规模下，两种算法的运行时间如表（单位为秒）：

表：运行时间统计

       以数据集规模N=100000 的运行时间为基础，计算两种算法的理论运行时间，并绘制图表，结果如图：

图：暴力算法理论效率与实际效率比较

 

图：分治算法理论效率与实际效率比较

       理论效率与实际效率相吻合，可认为实验算法编写无误。

       效率分析：比较暴力法与分治法，可见分治法极大地提升了算法的运行效率，这是由于分治法通过巧妙地划分原始问题，避免了相当多的无效运算。

        两种方法的代码如下：

#include<iostream>
#include<map>
#include<unordered_set>
#include<math.h>
#include<chrono>
using namespace std::chrono;
using namespace std;

double ans = DBL_MAX;
double average_force = 0, average_DivideandConquer = 0;
//先初始化数组
const int length = 1000000;
struct Point {
    double x, y;
};
Point points[length];
void init() {
    for (int i = 0; i < length; ++i) {
        points[i].x = rand() + rand() / 10000.0;
        points[i].y = rand() + rand() / 10000.0;
    }
    cout << "未去重的数组初始化成功->";
}

struct eqfunc {
    bool operator()(const Point& p1, const Point& p2) const {
        return ((p1.x == p2.x) && (p1.y == p2.y));
    }
};
struct hashfunc {
    size_t operator()(const Point& P) const {
        return size_t(P.x * 201928 + P.y);
    }
};
void Deduplication() {
    //利用哈希函数进行去重
    unordered_set<Point, hashfunc, eqfunc> hash;
    //如果没有找到就直接加入，如果有找到就重新生成
    int i = 0;
    while (i < length) {

        //find函数：寻找某个元素是否已经在容器中，如果不在，就返回hash.end()
        if (hash.find(points[i]) == hash.end()) {
            hash.insert(points[i]);
            ++i;
        }
        else {
            points[i].x = rand() + rand() / 10000.0;
            points[i].y = rand() + rand() / 10000.0;
            continue;
        }
    }
    cout << "去重并重新赋值的数组生成成功" << endl;
}
//定义计算距离函数
double dis(Point p1, Point p2)
{
    if (abs(p1.x - p2.x) > ans || abs(p1.y - p2.y) > ans)
        return ans;
    return pow((p1.x - p2.x), 2) + pow((p1.y - p2.y), 2);
}
double dist(Point p1, Point p2)
{
    return sqrt(pow((p1.x - p2.x), 2) + pow((p1.y - p2.y), 2));
}
//蛮力算法
double force() {
    ans = DBL_MAX;
    for (int i = 0; i <= length - 1; i++) {
        if (i % 10000 == 0)cout << "i=" << i << ' ';//进度条
        for (int j = i + 1; j <= length - 1; j++) {
            ans = min(ans, dis(points[i], points[j]));
        }
    }
    cout << endl;
    return sqrt(ans);
}

void QuickSort(int begin, int end) {// 快速排序
    if (begin >= end)return;
    double key = points[begin].x;// 选取第一个数作为key
    int p1 = begin, p2 = end;
    while (p2 > p1) {
        while (points[p2].x >= key && p1 < p2) {
            p2--;
        }
        swap(points[p1], points[p2]);
        while (points[p1].x < key && p1 < p2) {
            p1++;
        }
        swap(points[p1], points[p2]);
    }
    QuickSort(begin, p1 - 1);
    QuickSort(p1 + 1, end);
}
void QuickSortY(int begin, int end, Point* temp) {// 快速排序
    if (begin >= end) return;
    double key = temp[begin].y;// 选取第一个数作为key
    int p1 = begin, p2 = end;

    // 循环：最终目的是使p1==p2，a[p1]==key，p1左边的数全部小于key，右边的数全部大于等于key
    while (p2 > p1) {
        while (temp[p2].y >= key && p1 < p2) {
            p2--;
        }
        swap(temp[p1], temp[p2]);
        while (temp[p1].y < key && p1 < p2) {
            p1++;
        }
        swap(temp[p1], temp[p2]);
    }
    QuickSortY(begin, p1 - 1, temp);
    QuickSortY(p1 + 1, end, temp);
}
double Divide(int left, int right) {
    if (left == right)
        return DBL_MAX;
    if (left + 1 == right)
        return dist(points[left], points[right]);
    int mid = (right + left) >> 1;
    double min_dist = min(Divide(left, mid), Divide(mid + 1, right));

    // 使用temp数组记录符合条件的点
    Point* temp = new Point[right - left + 1];
    int i_size = 0; // 符合条件的点的数量
    for (int i = left; i <= right; i++) {
        //如果点在最小距离内则记录入temp数组
        if (points[i].x > points[mid].x - min_dist && points[i].x < points[mid].x + min_dist)
            temp[i_size++] = points[i];
    }

    QuickSortY(0, i_size - 1, temp);// 对纵坐标进行排序

    for (int i = 0; i < i_size; i++) {
        for (int j = i + 1; j < i_size; j++) {
            // 第二个点出现在矩形外，说明不存在可更新最小距的点对
            if (temp[j].y - temp[i].y > min_dist)
                break;

            // 找到了新的最小距
            if (dist(temp[i], temp[j]) <= min_dist)
                min_dist = dist(temp[i], temp[j]);
        }
    }
    return min_dist;
}
//蛮力算法求解
double AnswerBy_Force() {
    auto begintime = system_clock::now();
    double ans = force();
    duration<double> diff = system_clock::now() - begintime;
    cout << "其中最小的点对距离是：" << ans << endl;
    cout << "用时是：" << diff.count() << "秒" << endl;
    average_force += diff.count();
    return ans;
}

double AnswerBy_DivideandConquer() {
    QuickSort(0, length - 1);
    auto begintime = system_clock::now();
    double ans = Divide(0, length - 1);
    duration<double> diff = system_clock::now() - begintime;
    cout << "其中最小的点对距离是：" << ans << endl;
    cout << "用时是：" << diff.count() << "秒" << endl;
    average_DivideandConquer += diff.count();
    return ans;
}

int main() {
    int testtimes = 10;
    const int i = testtimes;
    while (testtimes--) {
        srand(testtimes);//随机数种子，使实验可以完全重现
        cout << "第" << i - testtimes << "次重复测试，此时的数据库大小是" << length << ":" << endl;
        init();
        Deduplication();

        cout << "蛮力算法的结果如下：" << endl;
        double ans1 = AnswerBy_Force();
        cout << endl;

        cout << "分治算法的结果如下：" << endl;
        double ans2 = AnswerBy_DivideandConquer();
        if (abs(ans1 - ans2) > 0.00001) {// 浮点数可能存在溢出问题，判断不设置为(ans1!=ans2)
            cout << "----------------------------------------------------------------wrong!!" << endl;
        }

        cout << "测试完成" << endl;
        cout << "----------------------------------------------------------------------------------" << endl
            << "----------------------------------------------------------------------------------" << endl
            << endl;
    }
    cout << "在" << i << "次重复的情况下，蛮力算法平均用时为" << average_force / i << endl;
    cout << "在" << i << "次重复的情况下，分治算法平均用时为" << average_DivideandConquer / i << endl;
}

四、实验结论与体会

实验结论：

        本次实验通过编写代码，分别使用暴力法与分治法解决了不同数据集规模下平面空间中的最小点对距离问题，对两种算法的运行效率进行了统计、比较、绘图与分析，圆满完成了实验任务。

实验体会：

        本次实验的暴力算法如果不进行优化，将会导致运行时间极为漫长，为了避免这种情况，提高实验效率，我们将计算距离的函数与暴力算法函数本身进行了优化，取得了显著的效果。由此可知，当编写运行时间较长的算法时，要特别注意多次执行的函数或语句的效率，对这些部分进行优化也往往能取得较好的结果。

        除此之外，在实验初期，我们试图输出算法运行过程中的所有结果，但是后来查阅资料得知，过多的输出会导致程序运行效率降低，因此我们改变了输出逻辑：暴力算法每遍历一万个点就输出一次当前的遍历次数，保证过程可视的同时尽可能避免了输出带来的效率损失。

尾注：

        本实验是此课程的第二次实验，相比第一次实验有所进步，但仍未做到最佳优化，尚有许多不完善不成熟之处，仅供参考。后续的实验将会拥有更高的质量，敬请期待！

        如有疑问欢迎讨论，如有好的建议与意见欢迎提出，如有发现错误则恳请指正！