【模板】模拟退火 费马点以及最小球覆盖
最近学了一波模拟退火。个人感觉，就是随机算法，然后我们的目标点，一开始温度T高的时候会移动的剧烈，T小的时候移动的缓和（所以这就是为什么每一次移动的距离都是乘T）。然后真正的模拟退火是如果当前的tem比ans优，那就毫不犹豫地接受，否则则以一定概率接受。也就是那个exp（dE/T）> rand 那个。
然后还有爬山算法，就是只会一直找更优解，不接受差解，具体就是在模拟退火基础上，一直找最优解，找不到就降温（所以会陷入局部最优解的坑）。在网上嫖了一份代码（https://blog.csdn.net/just_sort/article/details/51648958）
1 while(t>eps) 2 { 3 bool fuck = 1; 4 while(fuck) 5 { 6 fuck = 0; 7 for(int i=0; i<4; i++) 8 { 9 z.x = s.x+dir[i][0]*t; 10 z.y = s.y+dir[i][1]*t; 11 double tmp = getSum(p,n,z); 12 if(ans>tmp) 13 { 14 ans = tmp; 15 s = z; 16 fuck= 1; 17 } 18 } 19 } 20 t*=delta; 21 }
求费马点：poj2420
题目链接：https://vjudge.net/problem/POJ-2420
一开始就以坐标平均点作为起始点，然后移动的时候虽然可以无脑的向四周移动，但是我们可以这样看，求出每一个点和当前点的dx，dy，然后肯定是当前点朝dx和dy方向移动。这样会更好
（这个是带上一定概率接受更差解的）
1 #include<iostream> 2 #include<cmath> 3 #include<cstdio> 4 #include<ctime> 5 #include<algorithm> 6 using namespace std; 7 #define eps 1e-6 8 int n; 9 const int N = 109; 10 struct Point{ 11 double x,y; 12 }p[N],now; 13 double dis(Point a,Point b){ 14 return sqrt( (a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y - b.y)); 15 } 16 double getsum(Point u){ 17 double res = 0; 18 for(int i = 1;i<=n;++i){ 19 res += dis(p[i],u); 20 } 21 return res; 22 } 23 int main(){ 24 scanf("%d",&n); 25 for(int i = 1;i<=n;++i){ 26 scanf("%lf %lf",&p[i].x,&p[i].y); 27 now.x += p[i].x; 28 now.y += p[i].y; 29 } 30 now.x /= n; now.y /= n; 31 double ans = getsum(now); 32 double T = 100; 33 srand(23333); 34 while(T > eps){ 35 double dx = 0,dy = 0; 36 for(int i = 1;i<=n;++i){ 37 dx += (p[i].x - now.x)/dis(now,p[i]); 38 dy += (p[i].y - now.y)/dis(now,p[i]); 39 } 40 Point next = (Point){now.x + dx*T, now.y + dy*T}; 41 double tem = getsum(next); 42 if(tem < ans){ 43 ans = tem; 44 now = next; 45 } 46 else if(exp( (ans - tem) / T) > (rand()%10000)/10000.0){ 47 ans = tem; 48 now = next; 49 } 50 T *= 0.98; 51 } 52 printf("%.0f",ans); 53 return 0; 54 }
最小球覆盖：poj2069
题目链接：https://vjudge.net/problem/POJ-2069
和上一题差不多，先以平均点做起始点，然后每次移动肯定是往距离当前点的最远点移动。
（这个没有接受更差解了）
1 #include<iostream> 2 #include<cstdio> 3 #include<cmath> 4 #define eps 1e-6 5 using namespace std; 6 const int N = 40; 7 int n; 8 struct Point{ 9 double x,y,z; 10 }p[N],now; 11 double dis(Point u,Point v){ 12 return sqrt( (u.x-v.x)*(u.x-v.x) + (u.y-v.y)*(u.y-v.y) + (u.z - v.z)*(u.z - v.z)); 13 } 14 int main(){ 15 while(~scanf("%d",&n) && n){ 16 for(int i = 1 ; i <= n ;++i){ 17 scanf("%lf %lf %lf",&p[i].x,&p[i].y,&p[i].z); 18 now.x+=p[i].x; now.y+=p[i].y; now.z += p[i].z; 19 } 20 now.x/=n; now.y/=n; now.z/=n; 21 double T = 100; 22 double ans = 1e99; 23 while(T > eps){ 24 int k = 1; 25 for(int i = 1;i <= n ;++i){ 26 if(dis(now,p[i]) > dis(now,p[k]) ){ 27 k = i; 28 } 29 } 30 double tem = dis(p[k],now); 31 ans = min(ans,tem); 32 now.x += ( (p[k].x - now.x)/tem )*T; 33 now.y += ( (p[k].y - now.y)/tem )*T; 34 now.z += ( (p[k].z - now.z)/tem )*T; 35 T *= 0.98; 36 } 37 printf("%.6f\n",ans); 38 } 39 return 0; 40 }
至于最小园覆盖，模拟退火的话和最小求覆盖一样啦。不过最小圆覆盖有线性的点增量法（下一篇讲到），所以就不妨在这里了。
简单多边形塞入最大的圆：hdu3644
题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=3644
真的是模拟退火的奇淫怪招。这题区别就是，之前的模拟退火都是一个点在跳来跳去，这一题是多个点（初始化为所有边的中点）跳来跳去，并且各自记录这多个点之前跳到最优解的状态，并且一个T每个点跳多次。
1 /************************************************************************* 2 > File Name: hdu3644.cpp 3 # File Name: hdu3644.cpp 4 # Author : xiaobuxie 5 # QQ : 760427180 6 # Email:760427180@qq.com 7 # Created Time: 2019年10月31日 星期四 15时57分41秒 8 ************************************************************************/ 9 10 #include<bits/stdc++.h> 11 #include<iostream> 12 #include<cstdio> 13 #include<map> 14 #include<cmath> 15 #include<cstring> 16 #include<set> 17 #include<queue> 18 #include<vector> 19 #include<algorithm> 20 using namespace std; 21 typedef long long ll; 22 #define inf 0x3f3f3f3f 23 #define eps 1e-5 24 const double pi = acos(-1.0); 25 int sgn(double x){ 26 if(fabs(x) < eps) return 0; 27 if(x<0) return -1; 28 return 1; 29 } 30 const int N = 55; 31 struct Point{ 32 double x,y,r; 33 Point(){x=y=0;} 34 Point(double a,double b){x=a,y=b;} 35 Point operator - (const Point& b)const 36 {return Point(x-b.x,y-b.y);} 37 Point operator + (const Point& b)const 38 {return Point(x+b.x,y+b.y);} 39 Point operator * (const double& b)const 40 {return Point(x*b,y*b);} 41 double dot (const Point& b)const 42 {return x*b.x+y*b.y;} 43 double cross (const Point& b,const Point& c)const 44 {return (b.x - x)*(c.y - y) - (c.x - x)*(b.y - y);} 45 double dis (const Point& b)const 46 {return sqrt( (x-b.x)*(x-b.x) + (y-b.y)*(y-b.y) );} 47 bool OnLine (const Point& st,Point& ed)const 48 {return !sgn(cross(st,ed));} 49 bool OnSeg (const Point& st,Point& ed)const 50 {return OnLine(st,ed) && (*this - ed).dot(*this - st) < eps;} 51 }p[N],tp[N]; 52 struct Segment{ 53 Point s,e; 54 }; 55 int SegCross(Segment a,Segment b){ 56 double x1 = a.s.cross(a.e,b.s); 57 double x2 = a.s.cross(a.e,b.e); 58 double x3 = b.s.cross(b.e,a.s); 59 double x4 = b.s.cross(b.e,a.e); 60 if( b.e.OnLine(a.s,a.e) && b.s.OnSeg(a.s,a.e) || 61 b.s.OnLine(a.s,a.e) && b.e.OnSeg(a.s,a.e) || 62 a.e.OnLine(b.s,b.e) && a.s.OnSeg(b.s,b.e) || 63 a.s.OnLine(b.s,b.e) && a.e.OnSeg(b.s,b.e) ) 64 return 2; 65 else if( sgn(x1*x2) <= 0 && sgn(x3*x4) <= 0) return 1; 66 return 0; 67 } 68 int PointIn(Point p0,Point p[],int n){ 69 Segment L,S; 70 Point temh,teml; 71 L.s = p0; L.e = Point(inf,p0.y); 72 int cnt = 0; 73 for(int i = 1;i<=n;++i){ 74 S.s = p[i-1]; S.e = p[i]; 75 if(p0.OnSeg(S.s,S.e)) return 0; 76 if(S.s.y == S.e.y) continue; 77 if(S.s.y > S.e.y) temh = S.s,teml = S.e; 78 else temh = S.e,teml = S.s; 79 if(temh.OnSeg(L.s,L.e)) cnt++; 80 else if(SegCross(L,S) == 1 && teml.OnSeg(L.s,L.e) == 0) ++cnt; 81 } 82 if(cnt%2) return -1; 83 return 1; 84 } 85 double seg_point_dis(Point s,Point e,Point o){ 86 double a,b,c; 87 a = o.dis(s); 88 b = o.dis(e); 89 c = s.dis(e); 90 if(sgn(a*a + c*c - b*b) < 0) return a; 91 if(sgn(b*b + c*c - a*a) < 0) return b; 92 return fabs( o.cross(s,e)/c ); 93 94 } 95 double cal(Point p0,Point p[],int n){ 96 double res = inf; 97 for(int i = 0;i<n;++i){ 98 double d = seg_point_dis(p[i],p[i+1],p0); 99 res = min(res,d); 100 } 101 return res; 102 } 103 int main(){ 104 srand(time(0)); 105 int n; 106 double r; 107 while(~scanf("%d",&n) && n){ 108 Point now = (Point){0,0}; 109 double mxx = 0,mix = inf ,mxy = 0, miy = inf; 110 for(int i = 0;i<n;++i){ 111 scanf("%lf %lf",&p[i].x,&p[i].y); 112 now.x += p[i].x; 113 now.y += p[i].y; 114 mxx = max(mxx,now.x); mix = min(mix,now.x); 115 mxy = max(mxy,now.y); miy = min(miy,now.y); 116 } 117 p[n] = p[0]; 118 for(int i = 0;i<n;++i){ 119 tp[i].x = (p[i].x + p[i+1].x)/2; 120 tp[i].y = (p[i].y + p[i+1].y)/2; 121 tp[i].r = 0; 122 } 123 mxx -= mix; mxy -= miy; 124 scanf("%lf",&r); 125 double T = sqrt(mxx*mxx + mxy*mxy)/2; 126 bool ok = 0; 127 while(T > eps && (!ok) ){ 128 for(int i = 0;i<n && (!ok) ;++i){ 129 for(int u = 1;u<=30 && (!ok);++u){ 130 double ang = (rand()%1000)/1000.0 * 2*pi; 131 Point tem = Point(tp[i].x + T*cos(ang),tp[i].y + T*sin(ang)); 132 if(PointIn(tem,p,n) == -1 ){ 133 tem.r = cal(tem,p,n); 134 if(tem.r > tp[i].r ) tp[i] = tem; 135 if(sgn(tp[i].r - r) >= 0){ 136 ok = 1; 137 } 138 } 139 } 140 } 141 T*=0.6; 142 } 143 if(ok) puts("Yes"); 144 else puts("No"); 145 } 146 return 0; 147 }
所以模拟退火的姿势还是很重要的。。。
