模拟退火 学习笔记
模拟退火挺没有正确性的。真的是个玄学。
基本思想：如果新的状态的解更优则修改答案，否则以一定的概率接受新状态，这个概率为 $\exp(\frac{-\Delta E}{t})$ 。
考前学学这个 ++rp 。
求 $n(n \le 10^3)$ 个点的带权类费马点。
三角形的不带权费马点就是 $PA + PB + PC$ 取到最小值的点 $P$ 。
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
#define N ((int)(1e3+15))
int n,x[N],y[N],w[N],st=clock();
double ansx, ansy, dis; int _;
mt19937_64 _rd(time(0) + (unsigned int)(&_));
double rd() { return (double)_rd() / UINT64_MAX; }
double Time() { return ((double)clock()-st) / CLOCKS_PER_SEC; }
double calc(double nx, double ny)
{
double res = 0;
for(int i=1; i<=n; i++)
{
double dx = x[i] - nx, dy = y[i] - ny;
res += sqrt(dx * dx + dy * dy) * w[i]; // 带权类费马点
}
if(res < dis) { dis = res, ansx = nx, ansy = ny; }
return res;
}
void sa() // simulateAnneal
{
double t = 1e5, x = ansx, y = ansy; // 初始值
for(; t > 1e-3; t *= 0.997)
{
double tx = x + t * (rd() * 2 - 1); // [-1,1]
double ty = y + t * (rd() * 2 - 1);
double delta = calc(tx,ty) - calc(x,y);
if(exp( -delta / t ) > rd()) tie(x,y) = tie(tx,ty);
}
for(int i=1; i<=1000; i++)
{
double tx = ansx + t * (rd() * 2 - 1);
double ty = ansy + t * (rd() * 2 - 1);
calc(tx,ty);
}
}
signed main()
{
ios::sync_with_stdio(0);
cin.tie(0);cout.tie(0);
// freopen("check.in","r",stdin);
// freopen("check.out","w",stdout);
cout << fixed << setprecision(3);
cin >> n;
for(int i=1; i<=n; i++)
{
cin >> x[i] >> y[i] >> w[i];
ansx += x[i]; ansy += y[i];
}
ansx /= n; ansy /= n; dis = calc(ansx, ansy); // 初值
while(Time() < 0.85) sa();
cout << ansx << ' ' << ansy << '\n';
return 0;
}这题求的是最大值，所以要改成（根据图像性质可知）
if(exp(-delta / t) < rd()) tie(nx,ny) = tie(tx,ty); // 最大值， > 改成 < 号如果要拿 AC ，建议直接面向数据调 seed 。
否则考试的时候就随机好了，亲测期望得分不低。（非暴力部分就搞这个，多香）
代码：（纯随机版，期望得分 60~80pts ）
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define double long double
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
const double eps = 1e-13;
#define N ((int)(1e3+15))
double st = clock(), ansx, ansy;
int n,m,R,ans,mx,x[N],y[N],r[N],p[N],q[N],_;
mt19937 _rd(time(0) + (unsigned int)&_);
double rd() { return (double)_rd() / UINT_MAX; }
double dis(double x, double y) { return sqrt(x * x + y * y);}
double Time() { return ((double)clock()-st) / CLOCKS_PER_SEC; }
int dcmp(double x) { if(fabs(x) <= eps) return 0; return x > eps ? 1 : -1; }
double calc(double nx, double ny)
{
double now = R; int res = 0;
for(int i=1; i<=n; i++)
{
double dx = nx - x[i], dy = ny - y[i];
now = min(now, dis(dx, dy) - r[i]); if(dcmp(now) < 0) return 0;
}
for(int i=1; i<=m; i++) if(dcmp(dis(nx - p[i],ny - q[i]) - now) <= 0) ++res;
if(res > ans) { ans = res; ansx = nx; ansy = ny; }
return res;
}
void sa()
{
double nx = ansx, ny = ansy, t = mx;
const double t0 = (double)mx / 1e10;
int cnt = 0;
for(; dcmp(t-t0) >= 0; ++cnt, t *= 0.99)
{
double tx = nx + t * (rd() * 2 - 1);
double ty = ny + t * (rd() * 2 - 1);
double delta = calc(tx,ty) - calc(nx,ny);
if(exp(-delta / t) < rd()) tie(nx,ny) = tie(tx,ty);
}
for(int i=1; i<=cnt; i++)
{
double tx = ansx + t * (rd() * 2 - 1);
double ty = ansy + t * (rd() * 2 - 1);
calc(tx,ty);
}
}
signed main()
{
ios::sync_with_stdio(0);
cin.tie(0);cout.tie(0);
// freopen("check.in","r",stdin);
// freopen("check.out","w",stdout);
cin >> n >> m >> R;
for(int i=1; i<=n; i++)
{
cin >> x[i] >> y[i] >> r[i];
mx = max({mx,x[i], y[i], r[i]});
}
for(int i=1; i<=m; i++)
{
cin >> p[i] >> q[i];
ansx += p[i]; ansy += q[i]; mx = max({mx,p[i],q[i]});
}
ansx /= m; ansy /= m; ans = calc(ansx,ansy);
while(Time() < 0.85) sa(); cout << ans << '\n';
return 0;
}
