看到題面蒟蒻忽然心塞...

題目大意:

求二維平面上質點組重心.

模擬退火是一種奇妙的搜索算法(似乎比之前的A*更惡心???)

引用一下wiki(不存在的):

模擬退火是一種通用概率算法，用來在固定時間內尋求在一個大的搜尋空間內找到的最優解。
模擬退火來自冶金學的專有名詞退火。退火是將材料加熱後再經特定速率冷卻，目的是增大晶粒的體積，並且減少晶格中的缺陷。材料中的原子原來會停留在使內能有局部最小值的位置，加熱使能量變大，原子會離開原來位置，而隨機在其他位置中移動。退火冷卻時速度較慢，使得原子有較多可能可以找到內能比原先更低的位置。
模擬退火的原理也和金屬退火的原理近似：我們將熱力學的理論套用到統計學上，將搜尋空間內每一點想像成空氣內的分子；分子的能量，就是它本身的動能；而搜尋空間內的每一點，也像空氣分子一樣帶有“能量”，以表示該點對命題的合適程度。算法先以搜尋空間內一個任意點作起始：每一步先選擇一個“鄰居”，然後再計算從現有位置到達“鄰居”的概率。

讓我們小小的舉個例子:

求這組數中最大的數,如果我們用貪心法,從左開始我們依次搜索1->2->3.發現3的值小於2,取得局部最優解.

然而最高點是4,mmp,看來貪心算法在此題是行不通的.

這時我們就要掏出神器:模擬退火!

事實上,模擬退火是一個腦抽+玄學的算法.比如此題中,如果我們在搜到3後,又腦抽的搜索了4,就能發現4的值比局部最優解更優.

模擬退火算法允許在合適範圍內嘗試選取較差解以達到整體解最優的情況.

然而,所謂的合適範圍該如何選取呢?翻開萬能wiki,我們找到如下說明:

1. 由一個產生函數從當前解產生一個位於解空間的新解；為便於後續的計算和接受，減少算法耗時，通常選擇由當前新解經過簡單地變換即可產生新解的方法，如對構成新解的全部或部分元素進行置換、互換等，註意到產生新解的變換方法決定了當前新解的鄰域結構，因而對冷卻進度表的選取有一定的影響。
2. 計算與新解所對應的目標函數差。因為目標函數差僅由變換部分產生，所以目標函數差的計算最好按增量計算。事實表明，對大多數應用而言，這是計算目標函數差的最快方法。
3. 判斷新解是否被接受，判斷的依據是一個接受準則，最常用的接受準則是Metropolis準則：若Δt > 0則接受S′作為新的當前解S，否則以概率exp（Δt/T）接受S′作為新的當前解S。
4. 當新解被確定接受時，用新解代替當前解，這只需將當前解中對應於產生新解時的變換部分予以實現，同時修正目標函數值即可。此時，當前解實現了一次叠代。可在此基礎上開始下一輪試驗。而當新解被判定為舍棄時，則在原當前解的基礎上繼續下一輪試驗。

(話說Metropolis是哪位大神)

(還有對Metropolis準則我做了一些小改動,原文中Δt是前值-後值,似乎與我們(好像只是我)的認知不符)

接下來讓我們解釋(kou hu)一下第三點的Metropolis準則:

Δt>0時後解優於前解,此時顯然優先選擇後解.

Δt<0時,聯系起模擬退火的內涵,我們默默掏出熱力學公式:

\({P_{dE}} = {e^{\frac{{dE}}{{kT}}}}\)

P表示出現玄學解的概率.

於是核心代碼如下:

1         sub = stat(now) - stat(tmp);
2         if(sub >= 0 || exp(sub / T) >= grand()) now = tmp;
3         T *= c_kappa;

然後,本題中,求解重心位置等同於找到點A,使圖上的點\({P_1},{P_2},...,{P_n}\)滿足\(\sum\limits_{i = 1}^n {d(a,{p_i})*{w_{{P_i}}}} \).

代碼如下:

 1 #include <cstdio>
 2 #include <cmath>
 3 #include <cstring>
 4 #include <algorithm>
 5 using namespace std;
 6 const int MAXS = 60 * 1024 * 32,N = 10005;
 7 const double inf = 1e17;double total = inf;
 8 const double c_epsilon = 1e-3;
 9 const double c_kappa = 0.99;
10 const double c_pi = acos(-1);
11 char buf[MAXS];
12 int n;
13 inline char __getchar()
14 {
15     static char buf[MAXS];static int len = 0,pos = 0;
16     if (pos == len)pos = 0,len = fread(buf,1,MAXS,stdin);if(pos == len)return EOF;return buf[pos++];
17 }
18 inline double getdouble()
19 {
20     int num = 0,sgn = 1;char ch = __getchar();for(;ch < ‘0‘ || ch > ‘9‘;ch = __getchar())if(ch = ‘-‘)sgn = 0;
21     for(;ch >= ‘0‘ && ch <= ‘9‘;ch = __getchar())num = (num << 3) + (num << 1) + (ch ^ ‘0‘);
22     return (double) (sgn?num:(-num));
23 }
24 struct point
25 {
26     double x,y,w;
27     point(double x_,double y_,double w_){x = x_;y = y_;w = w_;}
28     point(double x_,double y_){x = x_;y = y_;}
29     point(){}
30     void operator += (const point &a){x += a.x;y += a.y;}
31     void operator /= (const int &n){x /= n;y /= n;}
32     void operator = (const point &a){x = a.x;y = a.y;w = a.w;}
33 }p[N],ans,now;
34 inline double grand(){return (rand() % 1000 + 1) / 1000.0;}
35 inline double dist(const point &a,const point &b){return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));}
36 inline double stat(point a)
37 {
38     double res = 0.0;
39     for(int i = 1;i <= n;i++) res += dist(a,p[i]) * p[i].w;
40     if(res < total) total = res,ans = a;
41     return res;
42 }
43 void simanl()
44 {
45     double T = 100000.0,alpha,sub;
46     while(T > c_epsilon)
47     {
48         alpha = 2.0 * c_pi * grand();
49         point tmp(now.x + T * cos(alpha),now.y + T * sin(alpha));
50         sub = stat(now) - stat(tmp);
51         if(sub >= 0 || exp(sub / T) >= grand()) now = tmp;
52         T *= c_kappa;
53     }
54     T = c_epsilon;
55     for(int i = 1;i <= 1000;i++)
56     {
57         alpha = 2.0 * c_pi * grand();
58         point tmp(ans.x + T * cos(alpha) * grand(),ans.y + T * sin(alpha) * grand());
59         stat(tmp);
60     }
61 }
62 int main()
63 {
64     srand(19260817);
65     n = getdouble();for(int i = 1;i <= n;i++){double x = getdouble();double y = getdouble();double w = getdouble();p[i] = {x,y,w};now += p[i];}now /= n;
66     simanl();printf("%.3f %.3f\n",ans.x,ans.y);
67     return 0;
68 }

[bzoj 3680]吊打XXX