码迷，mamicode.com
首页 > 其他好文 > 详细

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

时间：2015-01-30 09:06:51      阅读：338      评论：0      收藏：0      [点我收藏+]

标签：算法   算法导论   poj   dinkelbach   图论   

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

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

  • 转化：设<V‘,E‘>为全图的一个环路子图，其中点权和为U，边权和为W，令除积U/W = y，则可以得到一个关于y的函数Z(y) = W*y-U（不同于0-1规划可能有对子集的不同选取，在特定的环路中此函数是完全线性的）。
  • 再设所有环路中最大的除积为y_max，则有二分：
  1. 当y < y_max，对于最大除积的环路，Z(y_max) < 0，即我们利用新生成的边权w_new = W[u][v]*y-U[v]，可以推出此环为一负环，对于其他环路可能存在负环，也可能不存在。
  2. 当y >= y_max，根据函数式可以推出所有环路都是不存在负环的。
注意：负环可以用spfa和bellman_ford算法，可以首先初始化所有点距离为0，即把所有点都当成初始点在全图中寻找负环，加入简化的并查集优化可以加快速度。
  • Dinkelbach算法：
  1. 假设图中有k个环，其中的除积按升序排列为：y1 < y2 < ... < yk。
  2. 则我们初始化一个可知的下界y0 < y1，显然对于所有环路都是能得到负环的，而在寻找负环的过程中我们是可以利用并查集加上点权和数组找出这个负环的U和W的，也就可以找到这个环路所对应的yi = U/W。
  3. 注意到我们每次找到负环求得的yi是逐渐递增直到yk，但是每次出现哪个负环是不可知，这样就能在平均意义上加快搜索速度。
  • 负环搜索次数简单分析：假设剩下n个负环，接下来我们搜索到其中每个负环的概率都是相等的。令x[i]为遇到yi后还需要搜索的期望次数，i = 0,1,2,...,k，其中x[0]就是我们要求的平均次数。
  1. 易知x[k] = 0，x[k-1] = 1。
  2. 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)。
  3. 得到递推公式：x[i]-x[i+1] = 1/(k-i)。
  4. 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;
}


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

标签：算法   算法导论   poj   dinkelbach   图论   

原文地址：http://blog.csdn.net/dingzuoer/article/details/43282291
踩
(0)
赞
(0)
   
举报
评论 一句话评论（0）
登录后才能评论！
分享档案
更多>
2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)
周排行
mamicode.com排行更多图片
更多
友情链接
兰亭集智   国之画   百度统计   站长统计   阿里云   chrome插件   新版天听网
关于我们 - 联系我们 - 留言反馈
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码！