POJ 3621--Sightseeing Cows(0-1规划求最大密度)

最新推荐文章于 2019-12-02 17:10:24 发布

置顶 fu_ding1991

最新推荐文章于 2019-12-02 17:10:24 发布

阅读量682

点赞数

分类专栏： BFS/DFS POJ 图论数学文章标签：算法算法导论 poj Dinkelbach 图论

本文链接：https://blog.csdn.net/dingzuoer/article/details/43282291

版权

POJ 同时被 3 个专栏收录

66 篇文章 0 订阅

订阅专栏

BFS/DFS

32 篇文章 0 订阅

订阅专栏

图论

24 篇文章 0 订阅

订阅专栏

题意：给定一个<V,E>规模分别为L和P的有向图，每个点有点权，每个有向边有边权，请找出一个环路，使此环的点权和与边权和的除积最大。

题解：此类求两个集合和的最大除积，类似于0-1规划问题可以先转化为线性式子，然后利用二分或者Dinkelbach算法求得最大或者最小值。

转化：设<V',E'>为全图的一个环路子图，其中点权和为U，边权和为W，令除积U/W = y，则可以得到一个关于y的函数Z(y) = W*y-U（不同于0-1规划可能有对子集的不同选取，在特定的环路中此函数是完全线性的）。
再设所有环路中最大的除积为y_max，则有二分：

当y < y_max，对于最大除积的环路，Z(y_max) < 0，即我们利用新生成的边权w_new = W[u][v]*y-U[v]，可以推出此环为一负环，对于其他环路可能存在负环，也可能不存在。
当y >= y_max，根据函数式可以推出所有环路都是不存在负环的。

注意：负环可以用spfa和bellman_ford算法，可以首先初始化所有点距离为0，即把所有点都当成初始点在全图中寻找负环，加入简化的并查集优化可以加快速度。

Dinkelbach算法：

假设图中有k个环，其中的除积按升序排列为：y1 < y2 < ... < yk。
则我们初始化一个可知的下界y0 < y1，显然对于所有环路都是能得到负环的，而在寻找负环的过程中我们是可以利用并查集加上点权和数组找出这个负环的U和W的，也就可以找到这个环路所对应的yi = U/W。
注意到我们每次找到负环求得的yi是逐渐递增直到yk，但是每次出现哪个负环是不可知，这样就能在平均意义上加快搜索速度。

负环搜索次数简单分析：假设剩下n个负环，接下来我们搜索到其中每个负环的概率都是相等的。令x[i]为遇到yi后还需要搜索的期望次数，i = 0,1,2,...,k，其中x[0]就是我们要求的平均次数。

易知x[k] = 0，x[k-1] = 1。
x[i] = {x[i+1]+1+x[i+2]+1+...+x[k]+1}/(k-i) = 1+sum{x[j],j = i+1,i+2,...,k}/(k-i)。
得到递推公式：x[i]-x[i+1] = 1/(k-i)。
x[0] = 1+1/2+1/3+...+1/k。前k个倒数和近似于k的对数复杂度。

注意：因为精度问题，若直接用yi作为下一个求负环的输入，则仍然有可能得到此负环的除积yi，而导致错误的结果，可以加上一个小于输出精度的数，使上一个负环不成立并且不影响最终结果。

bellman_ford&二分：

#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;
#define maxN 1005
#define maxP 5005
#define eps (1e-4)

struct node
{
    short x,y;
    short val;
};

class solve
{
private:
    short funVal[maxN];
    node edges[maxP];
    int L,P;
    double maxFunUnit;
public:
    solve(int l,int p):L(l),P(p)
    {
        processIn();
        maxFunUnit = two_partion();
        if(maxFunUnit != -1)
        {
            printf("%.2lf\n",maxFunUnit);
        }
        else
        {
            printf("0.00\n");
        }
    }
    int processIn();
    double two_partion();
    int bellman_ford(double funUnit);
};

int solve::bellman_ford(double funUnit)
{
    short fa[maxN];
    double dis[maxN];
    double tmpDis;
    int i,j;
    char IsRelax;
    short x,y;
    for(i = 1;i <= L;i++)
    {
        dis[i] = 1e9;
    }
    for(i = 1;i <= L;i++)
    {
        fa[i] = i;
    }
    for(i = 0;i < L-1;i++)
    {
        IsRelax = false;
        for(j = 0;j < P;j++)
        {
            x = edges[j].x;
            y = edges[j].y;
            tmpDis = dis[x]+((double)edges[j].val)*funUnit-(double)funVal[y];
            if(dis[y] > tmpDis)
            {
                if(fa[x] == y)
                    return 1;
                dis[y] = tmpDis;
                fa[y] = fa[x];
                IsRelax = true;
            }
        }
        if(!IsRelax)
            return 0;
    }
    for(j = 0;j < P;j++)
    {
        x = edges[j].x;
        y = edges[j].y;
        tmpDis = dis[x]+((double)edges[j].val)*funUnit-(double)funVal[y];
        if(dis[y] > tmpDis)
        {
            return 1;
        }
    }
    return 0;
}

double solve::two_partion()
{
    if(!bellman_ford(-1.0))
        return -1;
    double left,right,mid;
    left = 1e-3;
    right = 1000;
    while(right-left > eps)
    {
        mid = (left+right)/2;
        if(bellman_ford(mid))      //有负环
        {
            left = mid;
        }
        else
        {
            right = mid;
        }
    }
    return right;
}

int solve::processIn()
{
    int i;
    for(i = 1;i <= L;i++)
    {
        scanf("%d",funVal+i);
    }
    for(i = 0;i < P;i++)
    {
        scanf("%d%d%d",&(edges[i].x),&(edges[i].y),&(edges[i].val));
    }
    return 0;
}

int main()
{
    int l,p;
    while(~scanf("%d%d",&l,&p))
    {
        solve sightseeing(l,p);
    }
    return 0;
}

spfa&二分：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<vector>
#include<queue>
using namespace std;
#define maxN 1005
#define eps (1e-4)

struct node
{
    int adjNode;
    int val;
};

class solve
{
private:
    short funVal[maxN];
    int L,P;
    double maxFunUnit;
    vector<vector<node> > adj;
public:
    solve(int l,int p):L(l),P(p)
    {
        processIn();
        maxFunUnit = two_partion();
        if(maxFunUnit != -1)
        {
            printf("%.2lf\n",maxFunUnit);
        }
        else
        {
            printf("0.00\n");
        }
    }
    int processIn();
    double two_partion();
    int spfa(double funUnit);
};

int solve::spfa(double funUnit)
{
    char vis[maxN];
    short fa[maxN];
    short depth[maxN];
    double dis[maxN];
    double tmpDis;
    int adjSize;
    int i;
    short now,next;
    queue<short> q;
    for(i = 1;i <= L;i++)
    {
        dis[i] = 1e9;
        q.push(i);
        fa[i] = i;
    }
    memset(vis,1,sizeof(vis));
    memset(depth,0,sizeof(depth));
    while(!q.empty())
    {
        now = q.front();
        q.pop();
        vis[now] = 0;
        adjSize = adj[now].size();
        for(i = 0;i < adjSize;i++)
        {
            next = adj[now][i].adjNode;
            tmpDis = dis[now]+((double)(adj[now][i].val))*funUnit-(double)funVal[next];
            if(dis[next] > tmpDis)
            {
                if(fa[now] == next)
                    return 1;
                dis[next] = tmpDis;
                fa[next] = fa[now];
                if(!vis[next])
                {
                    vis[next] = 1;
                    q.push(next);
                    depth[next]++;
                    if(depth[next] > L)
                    {
                        return 1;
                    }
                }
            }
        }
    }
    return 0;
}

double solve::two_partion()
{
    if(!spfa(-1.0))
        return -1;
    double left,right,mid;
    left = 1e-3;
    right = 1000;
    while(right-left > eps)
    {
        mid = (left+right)/2;
        if(spfa(mid))      //有负环
        {
            left = mid;
        }
        else
        {
            right = mid;
        }
    }
    return right;
}

int solve::processIn()
{
    int i;
    int a;
    node tmpNode;
    for(i = 1;i <= L;i++)
    {
        scanf("%d",funVal+i);
    }
    adj.resize(L+1);
    for(i = 0;i < P;i++)
    {
        scanf("%d%d%d",&a,&(tmpNode.adjNode),&(tmpNode.val));
        adj[a].push_back(tmpNode);
    }
    return 0;
}

int main()
{
    int l,p;
    while(~scanf("%d%d",&l,&p))
    {
        solve sightseeing(l,p);
    }
    return 0;
}

spfa&Dinkelbach：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<vector>
#include<queue>
#include<cmath>
using namespace std;
#define maxN 1005
#define eps (1e-3)

struct node
{
    int adjNode;
    int val;
};

class solve
{
private:
    short funVal[maxN];
    double sumVal[maxN];
    int L,P;
    double maxFunUnit;
    double nextFunUnit;
    vector<vector<node> > adj;
public:
    solve(int l,int p):L(l),P(p)
    {
        processIn();
        maxFunUnit = two_partion();
        printf("%.2lf\n",maxFunUnit);
    }
    int processIn();
    double two_partion();
    int spfa(double funUnit);
};

int solve::spfa(double funUnit)
{
    char vis[maxN];
    short fa[maxN];
    short depth[maxN];
    double dis[maxN];
    double tmpDis;
    int adjSize;
    int i;
    short now,next;
    queue<short> q;
    for(i = 1;i <= L;i++)
    {
        q.push(i);
        fa[i] = i;
    }
    memset(dis,0,sizeof(dis));
    memset(sumVal,0,sizeof(sumVal));
    memset(vis,1,sizeof(vis));
    memset(depth,0,sizeof(depth));
    while(!q.empty())
    {
        now = q.front();
        q.pop();
        vis[now] = 0;
        adjSize = adj[now].size();
        for(i = 0;i < adjSize;i++)
        {
            next = adj[now][i].adjNode;
            tmpDis = dis[now]+((double)(adj[now][i].val))*funUnit-(double)funVal[next];
            if(dis[next] > tmpDis)
            {
                if(fa[now] == next)
                {
                    nextFunUnit = funUnit/((tmpDis/(sumVal[now]+funVal[next]))+1.0);
                    //计算此负环的除积
                    return 1;
                }
                dis[next] = tmpDis;
                fa[next] = fa[now];
                sumVal[next] = sumVal[now]+funVal[next];
                //记录从祖先结点到此节点的所有新生成的边权和
                if(!vis[next])
                {
                    vis[next] = 1;
                    q.push(next);
                    depth[next]++;
                    if(depth[next] > L)
                    {
                        nextFunUnit = funUnit/((tmpDis/(sumVal[now]+funVal[next]))+1.0);
                        return 1;
                    }
                }
            }
        }
    }
    return 0;
}

double solve::two_partion()
{
    double left,pre_Left;
    int tmpRe;
    left = 1e-3;
    pre_Left = 100;
    while(fabs(left-pre_Left) > eps)
    {
        tmpRe = spfa(left);
        pre_Left = left;
        if(tmpRe)
        {
            left = nextFunUnit+0.00001;		
            //加上精度防止再找到此负环
        }
    }
    return left;
}

int solve::processIn()
{
    int i;
    int a;
    node tmpNode;
    for(i = 1;i <= L;i++)
    {
        scanf("%d",funVal+i);
    }
    adj.resize(L+1);
    for(i = 0;i < P;i++)
    {
        scanf("%d%d%d",&a,&(tmpNode.adjNode),&(tmpNode.val));
        adj[a].push_back(tmpNode);
    }
    return 0;
}

int main()
{
    int l,p;
    while(~scanf("%d%d",&l,&p))
    {
        solve sightseeing(l,p);
    }
    return 0;
}

fu_ding1991

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
POJ 3621--Sightseeing Cows(0-1规划求最大密度)

利用0-1规划中的二分以及Dinkelbach算法求解有向图中的最大密度问题，并且对此问题中的Dinkelbach算法的期望搜索迭代次数作了简要分析。
复制链接

扫一扫