模拟退火合集

hdu 3644 A Chocolate Manufacturer's Problem

首次接触模拟退火，看来还是挺神奇的。

主要参考这篇博文：http://blog.sina.com.cn/s/blog_7da04dd30100wlqj.html

题意判断多边形内部能否容纳一个半径为R的圆，即在有限的平面内找最优范围。遗传算法的结果难以掌控，爬山算法又没法保证跳出局部最优，所以基于贪心原则的模拟退火算法还是值得考虑的。

然后就是设定每次变化的步长和演化方式。该题可以从每条边的中点开始，向一个随机方向延伸，使找出的点到各条边的最小距离最大。

#include <stdio.h>
#include <vector>
#include <string>
#include <algorithm>
#include <queue>
#include <cstring>
#include <map>
#include <set>
#include <iostream>
#include <cmath>
using namespace std;
#ifdef __GNUC__
typedef long long LL;
inline void opt64(LL a) {
printf("%lld", a);
}
#else
typedef __int64 LL;
inline void opt64(LL a) {
printf("%I64d", a);
}
#endif
const int MAXN = 55, maxstps = 6;
const double inf = 1e20, eps = 1e-3, aeps = 1e-3;
int n;
double radius;
inline int cmp(double x, double y) {
return ((x-y)>aeps) - ((x-y)<-aeps);
}
inline int cmp(double x) {
return (x>aeps) - (x<-aeps);
}
struct _point
{
double x, y, r;
_point(double xx=0.0, double yy=0.0):x(xx),y(yy) {}
void ipt() {
scanf("%lf%lf", &x, &y);
}
double dot(const _point &a) {
return x*a.x+y*a.y;
}
double cross(const _point &a) {
return x*a.y-y*a.x;
}
_point operator + (const _point &a) {
return _point(x+a.x, y+a.y);
}
_point operator - (const _point &a) {
return _point(x-a.x, y-a.y);
}
double dis() {
return sqrt(x*x + y*y);
}
} da[MAXN], db[MAXN];
bool point_in_poly(_point p)
{
_point q = p, a, b;
q.x = inf;
int res = 0;
double x;
for (int i = 0; i< n; ++i)
{
a = da[i];
b = da[1+i];
if (cmp(a.y, b.y) == 0 ||
cmp(p.y, min(a.y, b.y)) == -1 ||
cmp(p.y, max(a.y, b.y)) >= 0 ) continue;
x = (p.y-a.y)*(b.x - a.x) / (b.y - a.y) + a.x;
if (cmp(x, p.x) == 0) return 0;
if (cmp(x, p.x) == 1) ++res;
}
return res%2;
}
double dis_Point_to_seg(_point o, _point a, _point b)
{
_point l1 = o-a, l2 = o-b, l3 = a-b;
if (cmp(l1.dot(l3)) == 1) return l1.dis();
if (cmp(l2.dot(l3)) == -1) return l2.dis();
return fabs(l1.cross(l2)) / l3.dis();
}
void cal(_point * a)
{
a->r = inf;
for (int i = 0; i< n; ++i)
a->r = min( dis_Point_to_seg(*a, da[i], da[i+1]), a->r);
}
int solve(double r)
{
_point tmp;
double agl;
while (r > eps)
{
for (int i = 0; i< n; ++i)
{
for (int j = 0; j< maxstps; ++j)
{
agl = rand();
tmp.x = db[i].x + cos(agl)*r;
tmp.y = db[i].y + sin(agl)*r;
if (point_in_poly(tmp))
{
cal(&tmp);
if (tmp.r > db[i].r)
{
db[i] = tmp;
if (cmp(tmp.r, radius) >= 0)
return 1;
}
}

}
}
r *= 0.9;
}
return 0;
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt", "r", stdin);
#endif
while (scanf("%d", &n) && n)
{
double maxx = -1, maxy = -1, minx = inf, miny = inf;
for (int i = 0; i< n; ++i)
da[i].ipt(), maxx = max(maxx, da[i].x), maxy = max(maxy, da[i].y),
minx = min(minx, da[i].x), miny = min(miny, da[i].y);
da
= da[0];
for (int i = 0; i< n; ++i)
{
db[i].x = (da[i].x + da[i+1].x)/2.0;
db[i].y = (da[i].y + da[i+1].y)/2.0;
db[i].r = 0.0;
}
maxx = maxx - minx;
maxy = maxy - miny;
scanf("%lf", &radius);
if (solve(sqrt(maxx*maxx + maxy*maxy)/2.0)) puts("Yes");
else puts("No");
}
return 0;
}

hdu 5017 Ellipsoid

找出椭球面上距离球心最近的点

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;
const double eps = 1e-10, r = 0.99;
double a, b, c, d, e, f, res;
double dir[8][2] = {0,1, 0,-1, -1,0, 1,0, -1,1, -1,-1, 1,1, 1,-1};
inline double dis (double x, double y, double z)
{
return sqrt(x*x + y*y + z*z);
}
inline int cmp(double x)
{
return (x>eps) - (x<-eps);
}
int cal(double A, double B, double C, double &z)
{
double q = B*B - 4*A*C;
if (cmp(q) < 0) return 0;
z = (sqrt(q) - B)/A/2.0;
return 1;
}
void solve()
{
double x=0, y=0, z, k = 1.0, ax=0, ay=0, pp;
cal(c, 0, -1, z);
res = dis(0, 0, z);
while (k > eps)
{
x = ax; y = ay;
for (int i = 0; i< 8; ++i)
{
double xx = x+dir[i][0]*k, yy = y+dir[i][1]*k;
if (cal(c, d*yy+e*xx, a*xx*xx+b*yy*yy+f*xx*yy-1, z) == 0) continue;
if ( cmp((pp=dis(xx, yy, z)) - res) < 0)
ax = xx, ay = yy, res = pp;
}
k *= r;
}
printf("%.7lf\n", res);
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt", "r", stdin);
//freopen("out.txt", "w", stdout);
#endif // ONLINE_JUDGE
while (scanf("%lf%lf%lf%lf%lf%lf", &a, &b, &c, &d, &e, &f) != EOF)
{
solve();
}
return 0;
}

hdu 3932 Groundhog Build Home

最小覆盖圆

#include <stdio.h>
#include <vector>
#include <string>
#include <algorithm>
#include <queue>
#include <cstring>
#include <map>
#include <set>
#include <iostream>
#include <cmath>
using namespace std;
#ifdef __GNUC__
typedef long long LL;
inline void opt64(LL a) {
printf("%lld", a);
}
#else
typedef __int64 LL;
inline void opt64(LL a) {
printf("%I64d", a);
}
#endif
const int MAXN = 1005, maxstps = 5;
const double inf = 1e20, eps = 1e-4, aeps = 1e-6;
int n;
double radius;
double maxx = -1, maxy = -1, minx = inf, miny = inf;
inline int cmp(double x, double y) {
return ((x-y)>aeps) - ((x-y)<-aeps);
}
inline int cmp(double x) {
return (x>aeps) - (x<-aeps);
}
struct _point
{
double x, y;
_point(double xx=0.0, double yy=0.0):x(xx),y(yy) {}
void ipt() {
scanf("%lf%lf", &x, &y);
}
_point operator - (const _point &a) {
return _point(x-a.x, y-a.y);
}
double dis() {
return sqrt(x*x + y*y);
}
} da[MAXN];
inline bool point_in_rect(double x, double y)
{
if (cmp(x-maxx) == 1 ||
cmp(x-minx) == -1 ||
cmp(y-maxy) == 1 ||
cmp(y-miny) == -1) return 0;
return 1;
}
void solve(double x, double y, double r, double mxdis)
{
double angle, xx=x, yy=y, ax, ay, d, c;
while (r > eps)
{
x = xx; y = yy;
for (angle = 0.0; angle< 6.28; angle += 0.1)
{
ax = x+cos(angle)*r; ay=y+sin(angle)*r;
if (point_in_rect(ax,ay))
{
_point ap(ax,ay);
c = -1;
for (int j = 0; j< n; ++j)
{
d = (ap-da[j]).dis();
if (cmp(d-c)==1)
c = d;
}
if (cmp(c-mxdis)==-1)
{
mxdis = c;
xx = ax; yy = ay;
}
}
}
r *= 0.7;
}
printf("(%.1lf,%.1lf).\n", x, y);
printf("%.1lf\n", mxdis);
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt", "r", stdin);
#endif
while (scanf("%lf%lf%d", &maxx, &maxy, &n) != EOF)
{
minx = 0.0, miny = 0.0;
double x1=inf, x2=-1, y1=inf, y2=-1, c=-1;
for (int i = 0; i< n; ++i)
da[i].ipt(), x1=min(x1,da[i].x),
x2=max(x2, da[i].x), y1=min(y1,da[i].y), y2=max(y2,da[i].y);
_point pp((x1+x2)/2, (y1+y2)/2);
for (int i = 0; i< n; ++i)
c = max(c, (pp-da[i]).dis());
solve(pp.x, pp.y, sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))*0.5, c);
}
return 0;
}