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

@

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;
        }
    }
}
posted @ 2019-04-02 21:09 Cwolf9 阅读(...) 评论(...) 编辑 收藏

Contact with me