专题3之模拟退火(蜜汁玄学)

这个算法真心玄学,网上还一堆假代码,假博客QAQ,还有并软用的样例。。。
还是总结总结吧
这次有四道题用了模拟退火Ellipsoid(也可以三分套三分)、Run Away、Groundhog Build Home、A Chocolate Manufacturer's Problem

常用

计算单个点到已知所有点的最值等问题

注意点

1、所有类型用double !!!!!。。。有三题分别WA在函数形参用int,中途变量用int,定义数组用int。。这种错误,在自己编译的时候编译器帮你强制转化了,放到OJ上就GG
2、根据题目要求的精度,设计T,delta和随机次数,太小了容易精度不够,太大了就容易TLE
3、GetVal的计算一定要仔细,笔者基本都是这里出错
4、HDU的。。。注意是多case!多case!多case!
5、不要漏要求输出的任何字符!!

接下来就发发题解

HDU 5017Ellipsoid

题意大致是求从原点到椭圆面的最短距离
解法大概有三分套三分,模拟退火。。。正所谓学了假高数,不会拉格朗日定理的我只能退火了

/*坑点大概就是一开始看了一份假题解
   后来WA就是没考虑方程的负解
   还有就是参数类型错误,一定要用double
*/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <algorithm>
#include <iostream>
using namespace std;
#define LL long long
#define CLR(x) memset(x, 0, sizeof(x))
const double INF = 1e10;
const double eps = 1e-9;

int dx[8] = {-1, -1, -1, 0, 0, 1, 1, 1}; //因为是三维坐标所有有八个方向走
int dy[8] = {-1, 0, 1, -1, 1, -1, 0, 1};

double a, b, c, d, e, f;

double dist(double x, double y, double z)
{
    return sqrt(x*x + y*y + z*z);
}

double GetVal(double x, double y) 
{ //啦啦解方程
    double A = c;
    double B = e*x + d*y;
    double C = a*x*x + b*y*y + f*x*y-1;
    double D  = B*B - 4*A*C; 
    if (D < 0) return INF;

    double det = sqrt(D);

    double res1 = (-B + det) / (2*A);
    double res2 = (-B - det) / (2*A);

    return dist(x,y,res1) < dist(x,y,res2) ? res1 : res2;
}

double Search()
{
    double T = 1; //初始温度
    double delta = 0.99;

    double x = 0, y = 0, z = sqrt(1/c);

    while (T > eps)
    {
        for (int i = 0; i < 8; i++)
        {
            double r = x + dx[i]*T;
            double c = y + dy[i]*T;
            double k = GetVal(r, c);
            if (k > INF) continue;
            if (dist(r,c,k) < dist(x,y,z))
            {
                x = r; 
                y = c;
                z = k;
            }
        }
        T *= delta;
    }
    return dist(x, y, z);

}

int main()
{
    while(scanf("%lf%lf%lf%lf%lf%lf", &a, &b, &c, &d, &e, &f) != EOF)
    {
        printf("%.7lf\n", Search());
    }    
}

HDU 1109Run Away

大致题意,给n个点,在平面内选一个点A,使A到所有点的距离和最小,并求最小距离。

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <algorithm>
#include <time.h>
#include <iostream>
using namespace std;
#define LL long long
#define CLR(x) memset(x, 0, sizeof(x))
const double INF = 1e8;
const double eps = 1e-3;

int dx[5] = {-1, 0, 1, 0}; //平面内
int dy[5] = {0, 1, 0, -1};

double X, Y;
int M;
double U[1005], V[1005];

struct node //数组和结构体都可用
{
    double x, y;
    double val;
} p[55];

double dist(double x, double y, double a, double b)
{
    return (x-a)*(x-a) + (y-b)*(y-b);
}

double GetSum(node h)
{
    double d = INF;
    for (int i = 0; i < M; i++)
        d = min(d, dist(h.x, h.y, U[i], V[i])); 
    return d;
}


double Search()
{
    double T = X+Y;
    double delta = 0.9;
    for (int i = 0; i <= 50; i++)
    {
        p[i].x = (rand()%1000+1) * 1.0/1000.00 * X; //队内dalao说可以暴力,不用随机
        p[i].y = (rand()%1000+1) * 1.0/1000.00 * Y;
        p[i].val = GetSum(p[i]);
    }
    node tmp;
    while (T > eps)
    {
        for (int i = 0; i <= 50; i++)
        {
            for (int j = 0; j <= 50; j++)
            {
                double r = (double)(rand()%1000+1)*1.0/1000.00*T;
                double c = sqrt(T*T - r*r);

                if (rand()&1) r *= -1;
                if (rand()&1) c *= -1;

                tmp.x = p[i].x + r;
                tmp.y = p[i].y + c;

                if (tmp.x<0 || tmp.x>X || tmp.y<0 || tmp.y>Y) continue;

                if (GetSum(tmp) > p[i].val)
                {
                    p[i] = tmp;
                    p[i].val = GetSum(tmp);
                }
            }
        }
        T *= delta;
    }
    int cur = 1;
    for (int i = 0; i <= 50; i++)
        if (p[i].val > p[cur].val)
            cur = i;
    printf("The safest point is (%.1lf, %.1lf).\n", p[cur].x, p[cur].y);
   //不要漏了句号!!!!血一般的教训
}

int main()
{
    srand(unsigned(time(0)));
    int t;
    scanf("%d", &t);
    while(t--)
    {
        CLR(U);
        CLR(V);
        scanf("%lf%lf%d", &X, &Y, &M);
        for (int i = 0; i < M; i++)
            scanf("%lf%lf", &U[i], &V[i]);
        Search();
    }
}

HDU 3932Groundhog Build Home

大致题意是一直土拨鼠要找一个离他所有朋友家都最近的点,并输出到最近朋友家的最大距离。
其实就是平面上的最小圆覆盖问题,这种单点到多点的用退火很很合适,然后就是,为了避免局域最大值印象答案,在这题和下一题,笔者都用分块随机数

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <algorithm>
#include <iostream>
using namespace std;
#define LL long long
#define CLR(x) memset(x, 0, sizeof(x))
#define MAX 100000000*1LL
#define PI acos(-1)
const double eps = 1e-2;
const double INF = 1e8;

double X, Y;
int N;
double px[1005], py[1005];

struct node
{
    double x, y, val;
} p[1005];

double dist(double x, double y, double a, double b)
{
    return sqrt((x-a)*(x-a) + (y-b)*(y-b));
}

double GetMax(node t)
{
    double d = 0;
    for (int i = 0; i < N; i++)
        d = max(d, dist(t.x, t.y, px[i], py[i])); //注意,这里是最大值
    return d;
}

void Search()
{
    double delta = 0.6; //!!!太接近1会TLE
    double T = max(X, Y);
    node tmp;

    while(T > eps)
    {
        for(int i = 0; i < N; i++)
        {
            for(int j = 0; j < 20; j++) //这里也是,过大会TLE
            {
                double d = (rand()%1000+1)/1000.0*2*PI;
                tmp.x = p[i].x + cos(d)*T;
                tmp.y = p[i].y + sin(d)*T;

                if (tmp.x<=0 || tmp.x>=X || tmp.y<=0 || tmp.y>=Y) continue;

                if (GetMax(tmp) < p[i].val)
                {
                    p[i] = tmp;
                    p[i].val = GetMax(tmp);
                }
            }
        }
        T *= delta;
    }

    node res;
    res.val = INF;;
    for (int i = 0; i < N; i++) // 求出全局最优解
        if (p[i].val < res.val) res = p[i];

    printf("(%.1lf,%.1lf).\n", res.x, res.y);
    printf("%.1lf\n", res.val);
}

int main()
{
    srand(unsigned(time(0)));
    while(scanf("%lf%lf%d", &X, &Y, &N) != EOF)
    {
        CLR(px);
        CLR(py);
        for (int i = 0; i < N; i++)
            scanf("%lf%lf", &px[i], &py[i]);

        for (int i = 0; i < N; i++) /*随机取N个点。。。笔者原本是最直接改上一题的代码。。。
                                      结果出事情了。。。事实告诉我们随机数要谨慎*/
        {
            p[i].x = (rand()%1000+1)/1000.0*X;
            p[i].y = (rand()%1000+1)/1000.0*Y;
            p[i].val = GetMax(p[i]);
        }
        Search();
    }
}

HDU 3644A Chocolate Manufacturer's Problem

大致题意求一个多边形内最大的内接圆

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <algorithm>
#include <iostream>
using namespace std;
#define LL long long
#define CLR(x) memset(x, 0, sizeof(x))
#define MAX 100000000*1LL
#define PI 3.14159265358
//#define zero(x) (((x)>0?(x):(-x))<eps)

const double eps = 1e-3;
const double INF = 1e8;

int n;
double R;
double MidX, MidY;

struct node
{
    double x, y;
    double R;
} p[105];

int zero(double x)
{ //不知道什么缘故,宏定义出事故了。
    if (fabs(x) < eps) return 0;
    else return x > 0 ? 1 : -1;
}

double dist(node a, node b)
{
    return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));
}

double cross(node k1, node k2, node k0)
{
    return (k1.x-k0.x)*(k2.y-k0.y) - (k2.x-k0.x)*(k1.y-k0.y);
}

double Disptosege(node& l1, node& l2, node& k ) //求解点到线段的距离!!
{
    double a, b, c;
    a = dist(l1, k);
    b = dist(l2, k);
    c = dist(l1, l2);

    if(zero(a*a+c*c-b*b) < 0) return a;
    else if(zero(b*b+c*c-a*a) < 0) return b;
    else return fabs(cross(k, l1, l2)/c);
}

bool inside_convex(node a)  //判断点是不是在多边形内
{
    int cnt = 0;
    double tmp;

    for(int i = 0; i < n; i++)
    {
        if( (p[i].y<=a.y && p[i+1].y>a.y) || (p[i+1].y<=a.y && p[i].y>a.y))
        {
            if (zero(p[i].y-p[i+1].y))
                tmp = p[i+1].x - (p[i+1].x-p[i].x)*(p[i+1].y-a.y)/(p[i+1].y-p[i].y);
            else
            {
                if(!zero(p[i].y-a.y)) cnt++;
                tmp = -INF;
            }
            if(zero(tmp - a.x) >= 0) cnt++;
        }
    }
    return cnt&1;
}


void cal(node& a)
{
    double t;
    a.R = INF;
    for( int i = 0; i < n; ++i )
    {
        t = Disptosege(p[i], p[i+1], a);
        a.R = min(a.R, t);
    }
}

node mid[105];

void Search()
{
    bool flag = false;
    double delta = 0.90;
    double T = sqrt(MidX*MidX + MidY*MidY)/2.0;
    int k = 0;
    while(!flag && T > eps)
    {
        for (int i = 0; i < n; i++)
        {
            k++;
            if (flag) break;
            for (int j = 0; j < 5; j++)
            {
                if (flag) break;
                node tmp;
                double angle = rand();
                tmp.x = mid[i].x + cos(angle)*T;
                tmp.y = mid[i].y + sin(angle)*T;

                if (inside_convex(tmp))
                {
                    cal(tmp);
                    if (tmp.R > mid[i].R)
                    {
                        mid[i] = tmp;
                        //mid[i].R = GetVal(tmp);
                        if (zero(tmp.R - R) >= 0) flag = true;
                    }
                }
            }
        }
        T *= delta;
    }
    //printf("%d\n", k);
    if (flag) printf("Yes\n");
    else printf("No\n");
}

int main()
{
    //srand(unsigned(time(0)));
    while(scanf("%d", &n)==1 && n)
    {
        double MaxX=0, MaxY = 0;
        double MinX=INF, MinY=INF;
        for (int i = 0; i < n; i++)
        {
            scanf("%lf%lf", &p[i].x, &p[i].y);
            MaxX = max(MaxX, p[i].x);
            MaxY = max(MaxY, p[i].y);
            MinX = min(MinX, p[i].x);
            MinY = min(MinY, p[i].y);
        }
        MidX = MaxX - MinX;
        MidY = MaxY - MinY;

        scanf("%lf", &R);
        p[n] = p[0];
        for (int i = 0; i < n; i++)
        {
            mid[i].x = (p[i].x + p[i+1].x) / 2; //取中点
                                               //队内一位dalao取右上点也过了,但是左下点一定WA
            mid[i].y = (p[i].y + p[i+1].y) / 2;
            mid[i].R = 0;
        }
        Search();
    }
}

总结
数据只能判断大体没问题,概率,温度,如何随机去点还是要根据情况定,如何逼近最优解

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,080评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,422评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,630评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,554评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,662评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,856评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,014评论 3 408
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,752评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,212评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,541评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,687评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,347评论 4 331
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,973评论 3 315
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,777评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,006评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,406评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,576评论 2 349

推荐阅读更多精彩内容