Gym101158 J 三分 or 模拟退火 Cover the Polygon with Your Disk

Gym101158 J: 求圆与给定凸多边形最大面积交

Gym101158 J: 求圆与给定凸多边形最大面积交

求 $10 $ 个点组成的凸多边形 $(convex\quad polygon) $ ，坐标范围 $[-100,100] $ ，与一个给定半径的圆的最大面积交。圆心的位置由你确定！

模拟退火坐标。

三分套三分求最优坐标，坑点是第二个三分的上下限需要动态求。

Steepest descent

Downhill simplex (Nelder–Mead)

Quasi-Newton (BFGS)

Evolution strategy (CMA-ES)

模拟退火

抄了个这场$Gym\quad rank1 $的模拟退火方法。

#include<bits/stdc++.h>

#define fi first

#define se second

#define iis std::ios::sync_with_stdio(false)

#define pb push_back

#define o2(x) (x)*(x)

#define db double

using namespace std;

typedef long long LL;

typedef pair<int, LL> pii;

#define cross(p1,p2,p3) ((p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y))

#define crossOp(p1,p2,p3) sign(cross(p1,p2,p3))

const int  MXN = 1e5 + 6;

const int INF = 0x3f3f3f3f;

const double eps = 1e-9;

const double PI = acos(-1.0);

int n, m;

struct P {

    db x, y;

    P(){};

    P(db _x, db _y):x(_x),y(_y){}

    P operator+(P p){return {x+p.x,y+p.y};}

    P operator-(P p){return {x-p.x,y-p.y};}

    db operator *(P s) {return x * s.x + y * s.y;}

    db operator ^(P s) {return x * s.y - y * s.x;}

    db dot(P p) {return x*p.x+y*p.y;}

    db det(P p) {return x*p.y-y*p.x;}

};

inline int sign(db a) {return a < -eps?-1:a>eps;}

inline int cmp(db a, db b) {return sign(a-b);}

db dis_pp(P a, P b) {return sqrt((b - a)*(b - a));}//两点距离

db dot_p(P a, P b, P c) {return (b - a) * (c - a);}//点乘

double disToLine(P c,P A,P B){return fabs(cross(A,B,c))/dis_pp(A,B);}

//

int judge(P a, P b, P c) { //判断c是否在ab线段上（前提是c在直线ab上）

    if(c.x >= min(a.x, b.x) && c.x <= max(a.x, b.x)

     && c.y >= min(a.y, b.y) && c.y <= max(a.y, b.y))return 1;

    return 0;

}

db area(P b, P c, db r) {

    P a(0.0, 0.0);

    if (dis_pp(b, c) < eps) return 0.0;

    db h = fabs(cross(a, b, c)) / dis_pp(b, c);

    if(dis_pp(a, b) > r - eps && dis_pp(a, c) > r - eps) { //两个端点都在圆的外面则分为两种情况

        db angle = acos(dot_p(a, b, c) / dis_pp(a, b) / dis_pp(a, c));

        if(h > r - eps) {

            return 0.5 * r * r * angle;

        } else if(dot_p(b, a, c) > 0 && dot_p(c, a, b) > 0) {

            db angle1 = 2 * acos(h / r);

            return 0.5 * r * r * fabs(angle - angle1) + 0.5 * r * r * sin(angle1);

        } else {

            return 0.5 * r * r * angle;

        }

    } else if(dis_pp(a, b) < r + eps && dis_pp(a, c) < r + eps) { //两个端点都在圆内的情况

        return 0.5 * fabs(cross(a, b, c));

    } else { //一个端点在圆上一个端点在圆内的情况

        if(dis_pp(a, b) > dis_pp(a, c)) { //默认b在圆内

            swap(b, c);

        }

        if(fabs(dis_pp(a, b)) < eps) { //ab距离为0直接返回0

            return 0.0;

        }

        if(dot_p(b, a, c) < eps) {

            db angle1 = acos(h / dis_pp(a, b));

            db angle2 = acos(h / r) - angle1;

            db angle3 = acos(h / dis_pp(a, c)) - acos(h / r);

            return 0.5 * dis_pp(a, b) * r * sin(angle2) + 0.5 * r * r * angle3;

        } else {

            db angle1 = acos(h / dis_pp(a, b));

            db angle2 = acos(h / r);

            db angle3 = acos(h / dis_pp(a, c)) - angle2;

            return 0.5 * r * dis_pp(a, b) * sin(angle1 + angle2) + 0.5 * r * r * angle3;

        }

    }

}

P myp[MXN], ar[MXN];

db get_s(int n, db R, P O) { //求圆与多边形面积并

    for (int i = 1; i <= n + 1; i++) myp[i] = ar[i];//顺或逆时针顺序

    for(int i = 1; i <= n + 1; i++) myp[i] = myp[i] - O;

    O = P(0, 0);

    db sum = 0;

    for(int i = 1; i <= n; i++) {

        int j = i + 1;

        db s = area(myp[i], myp[j], R);

        if(cross(O, myp[i], myp[j]) > 0) sum += s;

        else sum -= s;

    }

    return fabs(sum);

}

//

bool cmp1(P a, P b) { //3 4 1 2

    double p = atan2(a.y, a.x), q = atan2(b.y, b.x);

    if(p != q) return p < q;

    return a.x < b.x;

}

bool cmp2(P a, P b) { //2 3 4 1 2

    P c(0, 0);//尽量别单独使用

    double tmp = cross(a, b, c);

    if(tmp == 0) return a.x < b.x;

    return tmp > 0;

}

int Quadrant1(P a) { //象限排序,注意x负半轴

    if(a.x>0&&a.y>=0)  return 3;

    if(a.x<=0&&a.y>0)  return 4;

    if(a.x<0&&a.y<=0)  return 1;

    if(a.x>=0&&a.y<0)  return 2;

    return -1;

}

bool cmp5(P &a, P &b) {

    int qa = Quadrant1(a), qb = Quadrant1(b);

    if(qa == qb) return cmp2(a, b);

    return qa < qb;

}

bool intersect(db l1, db r1, db l2, db r2) {

    if(l1 > r1) swap(l1, r1);if(l2 > r2) swap(l2, r2);

    return !(cmp(r1,l2)==-1||cmp(r2,l1)==-1);

}

bool isSS(P p1, P p2, P q1, P q2) {//判断线段相交

    return intersect(p1.x,p2.x,q1.x,q2.x)&&intersect(p1.y,p2.y,q1.y,q2.y)&&

    crossOp(p1,p2,q1)*crossOp(p1,p2,q2)<=0&&crossOp(q1,q2,p1)*crossOp(q1,q2,p2)<=0;

}

bool isSS_strict(P p1, P p2, P q1, P q2) {//严格相交

    return crossOp(p1,p2,q1)*crossOp(p1,p2,q2)<0&&crossOp(q1,q2,p1)

    *crossOp(q1,q2,p2) < 0;

}

bool isMiddle(db a, db m, db b) {

    return sign(a-m)==0||sign(b-m)==0||(a<m!=b<m);

}

bool isMiddle(P a, P m, P b) {

    return isMiddle(a.x,m.x,b.x)&&isMiddle(a.y,m.y,b.y);

}

bool onSeg(P p1, P p2, P q) {

    return crossOp(p1,p2,q) == 0 && isMiddle(p1,q,p2);

}

bool ojbk(P p1, P p2, P q1, P q2) {

    P p = p2 - p1;

    P q = q2 - q1;

    if(q.det(p) == 0) return 1;

    return 0;

}

double rad(P p1, P p2) {

    return atan2(p1.det(p2),p1.dot(p2));

}

bool xielv(P p1, P p2, P q1, P q2) {

    P p = p2 - p1;

    P q = q2 - q1;

    if(q.det(p) == 0 && cmp(rad(p, q), 0) == 0) return 1;

    return 0;

}

double _R, ans;

P cur(0, 0);

mt19937 lh(std::clock());

double Rand() {return 1.0*(lh()%INF)/INF;}

void simulateAnneal() {

    ans = get_s(n, _R, cur);

    double step = 30;

    for(int iter = 0; iter < 2048; ++ iter) {

        P best = cur;

        db bestArea = ans;

        bool improve = 0;

        db dir = Rand();

        const int K = 50;

        db stepAngle = PI * 2/ K;

        for(int i = 0; i < K; ++i) {

            P nex{best.x+cos(dir)*step, best.y+sin(dir)*step};

            db area = get_s(n, _R, nex);

            if(area > bestArea) {

                bestArea = area;

                best = nex;

                improve = 1;

            }

            dir += stepAngle;

        }

        if(improve) {

            ans = max(ans, bestArea);

            cur = best;

            step *= 1.2;

        }else {

            step *= 0.9;

        }

    }

}

int main(){

    scanf("%d%lf", &n, &_R);

    double L = 1000, R = -1000;

    for(int i = 1; i <= n; i++) {

        scanf("%lf%lf", &ar[i].x, &ar[i].y);

        L = min(L, ar[i].x); R = max(R, ar[i].x);

        cur = cur + ar[i];

    }

    cur.x /= n, cur.y /= n;;

    reverse(ar+1, ar+1+n);

    ar[n + 1] = ar[1];

    simulateAnneal();

    printf("%.8f\n", ans);

    return 0;

}

三分套三分

先三分 $x $ ，再三分 $y $ 。

第一个三分的上下限是最小横坐标和最大横坐标。

第二个三分的上下限是由三分$x$出的这条直线与凸多边形交点的最小纵坐标和最大纵坐标。

#include<bits/stdc++.h>

#define fi first

#define se second

#define iis std::ios::sync_with_stdio(false)

#define pb push_back

#define o2(x) (x)*(x)

#define db double

using namespace std;

typedef long long LL;

typedef pair<int, LL> pii;

#define cross(p1,p2,p3) ((p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y))

#define crossOp(p1,p2,p3) sign(cross(p1,p2,p3))

const int  MXN = 1e5 + 6;

const int INF = 0x3f3f3f3f;

const double eps = 1e-9;

const double PI = acos(-1.0);

int n, m;

struct P {

    db x, y;

    P(){};

    P(db _x, db _y):x(_x),y(_y){}

    P operator+(P p){return {x+p.x,y+p.y};}

    P operator-(P p){return {x-p.x,y-p.y};}

    db operator *(P s) {return x * s.x + y * s.y;}

    db operator ^(P s) {return x * s.y - y * s.x;}

    db dot(P p) {return x*p.x+y*p.y;}

    db det(P p) {return x*p.y-y*p.x;}

};

inline int sign(db a) {return a < -eps?-1:a>eps;}

inline int cmp(db a, db b) {return sign(a-b);}

db dis_pp(P a, P b) {return sqrt((b - a)*(b - a));}//两点距离

db dot_p(P a, P b, P c) {return (b - a) * (c - a);}//点乘

double disToLine(P c,P A,P B){return fabs(cross(A,B,c))/dis_pp(A,B);}

//

int judge(P a, P b, P c) { //判断c是否在ab线段上（前提是c在直线ab上）

    if(c.x >= min(a.x, b.x) && c.x <= max(a.x, b.x)

     && c.y >= min(a.y, b.y) && c.y <= max(a.y, b.y))return 1;

    return 0;

}

db area(P b, P c, db r) {

    P a(0.0, 0.0);

    if (dis_pp(b, c) < eps) return 0.0;

    db h = fabs(cross(a, b, c)) / dis_pp(b, c);

    if(dis_pp(a, b) > r - eps && dis_pp(a, c) > r - eps) { //两个端点都在圆的外面则分为两种情况

        db angle = acos(dot_p(a, b, c) / dis_pp(a, b) / dis_pp(a, c));

        if(h > r - eps) {

            return 0.5 * r * r * angle;

        } else if(dot_p(b, a, c) > 0 && dot_p(c, a, b) > 0) {

            db angle1 = 2 * acos(h / r);

            return 0.5 * r * r * fabs(angle - angle1) + 0.5 * r * r * sin(angle1);

        } else {

            return 0.5 * r * r * angle;

        }

    } else if(dis_pp(a, b) < r + eps && dis_pp(a, c) < r + eps) { //两个端点都在圆内的情况

        return 0.5 * fabs(cross(a, b, c));

    } else { //一个端点在圆上一个端点在圆内的情况

        if(dis_pp(a, b) > dis_pp(a, c)) { //默认b在圆内

            swap(b, c);

        }

        if(fabs(dis_pp(a, b)) < eps) { //ab距离为0直接返回0

            return 0.0;

        }

        if(dot_p(b, a, c) < eps) {

            db angle1 = acos(h / dis_pp(a, b));

            db angle2 = acos(h / r) - angle1;

            db angle3 = acos(h / dis_pp(a, c)) - acos(h / r);

            return 0.5 * dis_pp(a, b) * r * sin(angle2) + 0.5 * r * r * angle3;

        } else {

            db angle1 = acos(h / dis_pp(a, b));

            db angle2 = acos(h / r);

            db angle3 = acos(h / dis_pp(a, c)) - angle2;

            return 0.5 * r * dis_pp(a, b) * sin(angle1 + angle2) + 0.5 * r * r * angle3;

        }

    }

}

P myp[MXN], ar[MXN];

db get_s(int n, db R, P O) { //求圆与多边形面积并

    for (int i = 1; i <= n + 1; i++) myp[i] = ar[i];//顺或逆时针顺序

    for(int i = 1; i <= n + 1; i++) myp[i] = myp[i] - O;

    O = P(0, 0);

    db sum = 0;

    for(int i = 1; i <= n; i++) {

        int j = i + 1;

        db s = area(myp[i], myp[j], R);

        if(cross(O, myp[i], myp[j]) > 0) sum += s;

        else sum -= s;

    }

    return fabs(sum);

}

//

bool cmp1(P a, P b) { //3 4 1 2

    double p = atan2(a.y, a.x), q = atan2(b.y, b.x);

    if(p != q) return p < q;

    return a.x < b.x;

}

bool cmp2(P a, P b) { //2 3 4 1 2

    P c(0, 0);//尽量别单独使用

    double tmp = cross(a, b, c);

    if(tmp == 0) return a.x < b.x;

    return tmp > 0;

}

int Quadrant1(P a) { //象限排序,注意x负半轴

    if(a.x>0&&a.y>=0)  return 3;

    if(a.x<=0&&a.y>0)  return 4;

    if(a.x<0&&a.y<=0)  return 1;

    if(a.x>=0&&a.y<0)  return 2;

    return -1;

}

bool cmp5(P &a, P &b) {

    int qa = Quadrant1(a), qb = Quadrant1(b);

    if(qa == qb) return cmp2(a, b);

    return qa < qb;

}

bool intersect(db l1, db r1, db l2, db r2) {

    if(l1 > r1) swap(l1, r1);if(l2 > r2) swap(l2, r2);

    return !(cmp(r1,l2)==-1||cmp(r2,l1)==-1);

}

bool isSS(P p1, P p2, P q1, P q2) {//判断线段相交

    return intersect(p1.x,p2.x,q1.x,q2.x)&&intersect(p1.y,p2.y,q1.y,q2.y)&&

    crossOp(p1,p2,q1)*crossOp(p1,p2,q2)<=0&&crossOp(q1,q2,p1)*crossOp(q1,q2,p2)<=0;

}

bool isSS_strict(P p1, P p2, P q1, P q2) {//严格相交

    return crossOp(p1,p2,q1)*crossOp(p1,p2,q2)<0&&crossOp(q1,q2,p1)

    *crossOp(q1,q2,p2) < 0;

}

bool isMiddle(db a, db m, db b) {

    return sign(a-m)==0||sign(b-m)==0||(a<m!=b<m);

}

bool isMiddle(P a, P m, P b) {

    return isMiddle(a.x,m.x,b.x)&&isMiddle(a.y,m.y,b.y);

}

bool onSeg(P p1, P p2, P q) {

    return crossOp(p1,p2,q) == 0 && isMiddle(p1,q,p2);

}

bool ojbk(P p1, P p2, P q1, P q2) {

    P p = p2 - p1;

    P q = q2 - q1;

    if(q.det(p) == 0) return 1;

    return 0;

}

double rad(P p1, P p2) {

    return atan2(p1.det(p2),p1.dot(p2));

}

bool xielv(P p1, P p2, P q1, P q2) {

    P p = p2 - p1;

    P q = q2 - q1;

    if(q.det(p) == 0 && cmp(rad(p, q), 0) == 0) return 1;

    return 0;

}

double _R;

double check1(double x) {

    double L = 10000, R = -10000, midl, midr, tmpl, tmpr, ans = 0;

    for(int i = 1; i <= n; ++i) {

        if(sign((ar[i].x-x)*(ar[i+1].x-x)) <= 0) {

            double tmp = ar[i].y + (ar[i+1].y-ar[i].y)/(ar[i+1].x-ar[i].x)*(x-ar[i].x);

            L = min(L, tmp);

            R = max(R, tmp);

        }

    }

    for(int i = 0; i < 80; ++i) {

        midl = (L + L + R) / 3;

        midr = (L + R + R) / 3;

        tmpl = get_s(n, _R, P{x, midl}); tmpr = get_s(n, _R, P{x, midr});

        if(tmpl > tmpr) {

            ans = max(ans, tmpl);

            R = midr;

        }else {

            ans = max(ans, tmpr);

            L = midl;

        }

    }

    return ans;

}

int main(){

    scanf("%d%lf", &n, &_R);

    double L = 1000, R = -1000;

    for(int i = 1; i <= n; i++) {

        scanf("%lf%lf", &ar[i].x, &ar[i].y);

        L = min(L, ar[i].x); R = max(R, ar[i].x);

    }

    reverse(ar+1, ar+1+n);

    ar[n + 1] = ar[1];

    double midl, midr, tmpl, tmpr, ans = 0;

    for(int i = 0; i < 80; ++i) {

        midl = (L + L + R) / 3;

        midr = (L + R + R) / 3;

        tmpl = check1(midl); tmpr = check1(midr);

        if(tmpl > tmpr) {

            ans = max(ans, tmpl);

            R = midr;

        }else {

            ans = max(ans, tmpr);

            L = midl;

        }

    }

    printf("%.8f\n", ans);

    return 0;

}

模拟退火套路

for(int tc = 0; tc < 100; tc++) {//模拟退火次数

    pdd o;//随机起点

    double d = 400;//初始温度

    double area;//最优解

    while (d > eps) {//终止温度

        pdd nex = o + rand();//随机移动

        double tot;//当前答案

        if (tot > area) o = nex, area = tot;//更新最优解

        d *= 0.99;//降温系数0.998,0.975

    }

    ans = max(ans, area);//更新答案

}

double _R, ans;

P cur(0, 0);

mt19937 lh(std::clock());

double Rand() {return 1.0*(lh()%INF)/INF;}

void simulateAnneal() {

    ans = get_s(n, _R, cur);

    double step = 30;//看题目范围咯

    for(int iter = 0; iter < 2048; ++ iter) {//或者改成>eps

        P best = cur;

        db bestArea = ans;

        bool improve = 0;

        db dir = Rand();

        const int K = 50;//试探50个方向

        db stepAngle = PI * 2 / K;

        for(int i = 0; i < K; ++i) {

            P nex{best.x+cos(dir)*step, best.y+sin(dir)*step};

            db area = get_s(n, _R, nex);

            if(area > bestArea) {

                bestArea = area;

                best = nex;

                improve = 1;

            }

            dir += stepAngle;

        }

        if(improve) {

            ans = max(ans, bestArea);

            cur = best;

            step *= 1.2;

        }else {

            step *= 0.9;

        }

    }

}

巴特西

Gym101158 J 三分 or 模拟退火 Cover the Polygon with Your Disk

Gym101158 J: 求圆与给定凸多边形最大面积交

模拟退火

三分套三分

模拟退火套路

最新文章

热门文章