最小圆(球)覆盖问题

如果早知道,只要贪心就能Y……

原题

43rd ACM-ICPC Nanjing Problem D

题意

在三维地图中选一个点,使得这个点到所有点的距离的最大值最小。点的个数$N<100$,时间限制3s。

思路

很明显就是最小球覆盖问题。 一开始以为是二分题,穷举所有可能的半径,判断球是否相交,相交就在$[l, mid]$区间内搜索,否则在$[mid, r]$区间内搜索。想想还挺对的(误)结果样例都过不了。

然后徐臣开始了他的骚操作:每次选出3个点找到球面的中垂线,然后覆盖球的球心肯定在上面,重复判断更新答案,$O(n^4)$。最后没写出来(写得出来有鬼了)

然后是我的烂操作:二分半径,然后将球心放在第一个点上,列举所有点,如果点不在球内,就二分找到最近的能够将它覆盖的球心位置,然后把球心移过去。最后判断所有点是否在球心内,然后继续二分,$O(n^2 \lg n)$。然后也是样例都没过。

最后剩下10分钟反正做不出就开始扫雷+数独了。

随机增量法

算法原理

首先随机排序,打乱输入的顺序。然后按顺序处理规模逐渐增大的子问题。

具体实现

(2维的比3维简单多了,类比一下就可以了)

随机排序后首先选择前2个点,求出他们的最小外接球(就是两点连线中点),然后依次判断后面的点。对于点$p_i$:

  • $p_i$在球上或球内,什么都不做;
  • $p_i$在球外,则一定在新的最小包围球的边界上。以$p_1 p_i$为直径建立新的球,判断其他点是否在球上或球内。
    • $p_j$在球上/球内,什么都不做;
    • $p_j$在球外,则一定与$p_1$和$p_i$共球面。找到这个最小的球面(三点平面过圆心),然后判断其他点是否在球上或球内。
      • $p_k$在球上/球内,什么都不做;
      • $p_k$在球外,则一定与$p_1$、$p_i$、$p_j$共球面。由于4个点确定一个球,此处找到的最大半径的球是最终答案。

昏睡数学

淦啊我要疯了这是人能做的吗???

已知直径

球心就是直径的中点。

已知3个点

要求过3个点的球最小,说明圆心与三点共面,只要求出三角形外心即可。求三角形外心首先要求出平面方程:设三点为$A$、 $B$、 $C$,外心坐标为$O(x_0, y_0, z_0)$,平面$ABC$的法向量为$$\vec{n} = \vec{AB} \times \vec{AC}$$

则对于任意的点$P \in ABC$,有$$ \vec{AP} \cdot \vec{n} = 0 $$

将$O$点带入点面式方程,化简得平面方程$$Ax+By+Cz+D=0$$

由$O$是三角形外心,得方程组

$$ \begin{cases}
2(x_a – x_b) x + 2(y_a – y_b) y + 2(z_a – z_b) z = (x_a^2+y_a^2+z_a^2) + (x_b^2+y_b^2+z_b^2) \\
2(x_a – x_c) x + 2(y_a – y_c) y + 2(z_a – z_c) z = (x_a^2+y_a^2+z_a^2) + (x_c^2+y_c^2+z_c^2) \\
Ax + By + Cz + D = 0
\end{cases} $$

用Gauss消元法即可求出球心坐标。

已知4个点

吐了!但只要排除共线、共面的情况,直接求解线性方程组就可以了!

$$ \begin{cases}
2(x_a – x_b) x + 2(y_a – y_b) y + 2(z_a – z_b) z = (x_a^2+y_a^2+z_a^2) + (x_b^2+y_b^2+z_b^2) \\
2(x_a – x_c) x + 2(y_a – y_c) y + 2(z_a – z_c) z = (x_a^2+y_a^2+z_a^2) + (x_c^2+y_c^2+z_c^2) \\
2(x_a – x_d) x + 2(y_a – y_d) y + 2(z_a – z_d) z = (x_a^2+y_a^2+z_a^2) + (x_d^2+y_d^2+z_d^2)
\end{cases} $$

模拟退火法

物理原理

样卷学习法学的物理,不懂!

算法原理

模拟退火法就是纯贪心的修改版。在搜索到某个局部最优解之后,普通的贪心算法就结束搜索了,但是模拟退火法会以一定的概率进行下一次搜索,这个概率随着时间的变化而减小。

伪代码

我的理解:

while (T > T_min)
  delta = S(n+1) - S(n)
  if (delta > 0)
    n++
  else
    if (exp(delta / T) > rand(0, 1))
      n++
      T = r * T // annealing
    else
      break

Wiki的理解:

Let s = s0
For k = 0 through kmax (exclusive):
  T ← temperature(k ∕ kmax)
  Pick a random neighbour, snew ← neighbour(s)
  If P(E(s), E(snew), T) ≥ random(0, 1):
    s ← snew
Output: the final state s

具体实现

选择任意一个点作为初始解,然后将所有点按到当前圆/球心的距离排序,向距离最远的点靠近。每次靠近的距离都会减小(相当于模拟退火法的概率逐渐减小),最后圆/球心原来越靠近最优解,半径趋近稳定。

对于每次靠近,靠近的距离都是相差的距离乘以一个祖传比例系数$r/temp$,其中$r$是当前位置的最大距离,$temp$是当前温度。至于为什么是这个系数,淦啊全世界是不是只有中国人用这个方法解决问题,我连篇论文都找不到!

复杂度分析

时间复杂度$O(n \lg T)$,其中$T$是选定的初始温度。$T$越大,答案越精确,花费时间越长。

补题

POJ 2069 Super Star (模拟退火)

网上找了好几个题解对着抄都WA,直接贴题解跑也WA,哭了 CSDN全是复读机真是垃圾

POJ是真坑爹,不支持C++0x和1x,不能用时间做随机种子(会RE),输出用%lf就WA,搞了半天改成%f就对了……

#include <cstdio>
#include <cmath>
#include <ctime>
#include <algorithm>
using namespace std;
const double e = exp(1.0);
const double eps = 1e-6;

class point {
 public:
  double x, y, z;
};

int n = 0;
point s[105] = {};

double getDis(point &a, point &b) {
  double dis2 = 0;
  dis2 += pow(a.x - b.x, 2);
  dis2 += pow(a.y - b.y, 2);
  dis2 += pow(a.z - b.z, 2);
  return sqrt(dis2);
}

double getMaxDis(point &c) {
  double maxDis = 0, nexDis = 0;
  for (int i = 0; i < n; ++i) {
    nexDis = getDis(c, s[i]);
    if (nexDis > maxDis) {
      maxDis = nexDis;
    }
  }
  return maxDis;
}

point getNextC(point &c, double temp) {
  int chosen = 0;
  double maxR = getDis(c, s[0]), nexR = 0;
  for (int i = 1; i < n; ++i) {
    nexR = getDis(c, s[i]);
    if (nexR > maxR) {
      maxR = nexR;
      chosen = i;
    }
  }
  return (point) {
      c.x + (s[chosen].x - c.x) / maxR * temp,
      c.y + (s[chosen].y - c.y) / maxR * temp,
      c.z + (s[chosen].z - c.z) / maxR * temp
  };
}

double getRandPct() {
  return (double) (rand() % 10000) / 10000.0;
}

double solve() {
  point c = (point) {0, 0, 0};
  double ans = getMaxDis(c);
  double temp = 1e5, rate = 0.998;
  while (temp > eps) {
    point nc = getNextC(c, temp);
    double nAns = getMaxDis(nc);
    if (nAns < ans) {
      c = nc;
      ans = nAns;
    } else {
      if (exp((ans - nAns) / temp) > getRandPct()) {
        c = nc;
        ans = nAns;
      }
    }
    temp *= rate;
  }
  return ans;
}

int main() {
  srand(12450); // f*ck POJ!
  while (scanf("%d", &n) != EOF && n != 0) {
    for (int i = 0; i < n; ++i) {
      scanf("%lf%lf%lf", &s[i].x, &s[i].y, &s[i].z);
    }
    printf("%.5f\n", solve()); // f*ck POJ!
  }
  return 0;
}

1条评论

发表评论

电子邮件地址不会被公开。 必填项已用*标注