计蒜客 A2141 - Country Meow

题目链接：ICPC 2018 南京现场赛 D 题 - Country Meow

思路

最小球覆盖问题，可以用模拟退火做。

昨天补 LC 第 49 场双周赛的时候看学会了模拟退火，于是翻出这道题。

模拟退火的框架是这样的：

int ans;

void score() {
    // ...
    // 计算分数

    // [!IMPORTANT] 在每次计算分数的时候更新最优解
    ans = max(ans, current_score);
}

void simulated_anneling() {
    // 随机取一个状态作为初始状态
    // 这个 Status 可能是 int，可能是坐标，etc.
    Status s = rand_status();

    // t 表示「温度」，应当介于 t_init 和 t_term 之间
    // rate 是一个接近 1 且小于 1 的数
    // rate 越接近 1，计算答案的过程就越平缓，但是相对也越慢
    double t_init = 1e6, t_term = 1e-5, rate = 0.98;
    for (double t = t_init; t > t_term; t *= rate) {
        // 随机一个增量
        // [!IMPORTANT] 这个增量是和 t 相关的
        // 例如 int d = t * ((double)rand() / RAND_MAX * 2 - 1)
        // 它意味着如果温度高，变化就快一些，温度低变化就慢一些
        // 这个有点类似于梯度下降
        Status d = rand_delta(t);

        // 计算增量前和增量后的分数
        int pre_score = score(s);
        int aft_score = score(s + d);

        // [!IMPORTANT]
        // 假设 aft_score > pre_score 表示增量后更优
        // 1. 如果 score_delta > 0，表示 s + d 更优，那么 s + d 直接被 s 接受
        //    即 s <- s + d
        // 2. 否则，让 s + d 以一定的概率被 s 接受
        int score_delta = aft_score - pre_score;
        if (exp(-score_delta / t) > (double)rand() / RAND_MAX) {
            s = s + d;
        }
    }
}

int main() {
    // ...
    // 输入数据

    // 单次模拟退火找到的很可能是局部最优解
    // 所以要多做几次
    int repeat_time = 50;
    for (int i = 0; i < repeat_time; ++i) {
        simulated_anneling();
    }

    // ...
    // 输出答案

    return 0;
}

• WA 时候可以尝试增加 rate 或者 repeate_time，或者缩小 t 的范围
• TLE 时可以尝试缩小 rate 或者 repeate_time，或者增大 t 的范围

在这一题中，我们把点的三维坐标作为状态，score 计算的就是当前点距离最远点的距离。

代码

#include <bits/stdc++.h>

using namespace std;

const int MAX_N = 100 + 5;

int n;
double ans;

struct Coordinate {
    double x, y, z;
} c[MAX_N], mn, mx;

double dist(Coordinate a, Coordinate b) {
    return sqrt(
        (a.x - b.x) * (a.x - b.x) + 
        (a.y - b.y) * (a.y - b.y) + 
        (a.z - b.z) * (a.z - b.z));
}

double calculate(Coordinate r) {
    double far_dis = 0;
    for (int i = 0; i < n; ++i) {
        far_dis = max(far_dis, dist(r, c[i]));
    }

    ans = min(ans, far_dis);
    return far_dis;
}

void silumated_anneling() {
    Coordinate cur = { 
        rand() % (int(2e5) + 1) - 1e5,
        rand() % (int(2e5) + 1) - 1e5,
        rand() % (int(2e5) + 1) - 1e5,
    };

    for (double t = 1e6; t > 1e-6; t *= 0.99) {
        Coordinate nex = { 
            cur.x + t * ((double)rand() / RAND_MAX * 2 - 1),
            cur.y + t * ((double)rand() / RAND_MAX * 2 - 1),
            cur.z + t * ((double)rand() / RAND_MAX * 2 - 1),
        };
        
        double cur_dis = calculate(cur);
        double nex_dis = calculate(nex);
        double delta = nex_dis - cur_dis;
        if (exp(-delta / t) > (double)rand() / RAND_MAX) {
            cur = nex;
        }
    }
}

int main() {
    scanf("%d", &n);

    for (int i = 0; i < n; i++) {
        scanf("%lf%lf%lf", &c[i].x, &c[i].y, &c[i].z);
    }

    srand(time(0));
    ans = INT_MAX;

    for (int i = 0; i < 100; ++i) {
        silumated_anneling();
    }

    printf("%.10f\n", ans);

    return 0;
}