K Subsequence【优先队列优化Dijkstra费用流】【2019HDU多校】【费用流势函数理解】

HDU-6611

---

  被卡了好多东西，T了好几次，然后竟然有MLE了，最后发现是reserve()好像会消耗很多内存的说，改成clear就过了。这道题卡了链式前向星的建边，因为这里有可能会变成建立了N*N的边，而其中N又特别大，所以这里用vector<>来存就好了。

  这里用到了优先队列的Dijkstra来优化费用流，因为SPFA的复杂度是O(N * M)，而优化的Dijkstra的复杂度只有O(N * N * logN)，似乎快了好多啊！

  然后这里理了好久的势函数，就是理不懂，Dijkstra是不能跑负边权的，但是加了这样的一个势函数就可以了，有点迷，我再去学一下，然后后期补上（已经补上）。

  我们对网络G中的每个点i设置一个势函数h[i]，在算法运行的每一时刻，对于每条残余网络中的边(x, y)，都要满足h[x] + w[x][y] – h[y] >= 0（三角不等式）。维护h[i]的方法下文再讨论。

  构建新网络G’，V’ = V，w'[x][y] = h[x] + w[x][y] – h[y]，这样在新网络中所有边的权值均非负，可以使用Dijkstra求解最短路。

G’中的一条路径p[1], p[2], …, p[m]的权值为

w’[path] =

w’[p[1]][p[2]] + w’[p[2]][p[3]] + … + w’[p[m – 1]][p[m]] =

h[p[1]] – h[p[2]] + w[p[1]][p[2]] +

h[p[2]] – h[p[3]] + w[p[2]][p[3]] +

…+

h[p[m – 1]] – h[p[m]] + w[p[m – 1]][p[m]] =

w[p[1]][p[2]] + w[p[2]][p[3]] + … + w[p[m – 1]][p[m]] + h[p[1]] – h[p[m]] =

w[path] + h[p[1]] – h[p[m]]

而当p[1]和p[m]确定时，最优化w[path]与最优化w’[path]是等价的，所以我们可以在G’中求解S到T的最短路w’[path]，G中S到T的最短路w[path] = w’[path] – h[S] + h[T]。

下面我们来讨论如何维护势函数h[]。

初始情况下，如果所有的权值都为正，那么可以简单地将所有h[i]设置为0；如果有负权值，那么我们做一遍spfa，让h[]等于距离函数。如果初始网络为DAG，也可以采用递推h[i] = min(h[j] + w[j][i], g[j][i] = true)。

让d’[i]表示G’中S到点i的距离，当某次增广结束后，G’中会新加入某些边(j, i)，而这些(j, i)必定满足d’[i] + w’[i][j] = d’[j]（否则(i, j)就不会在增广路中）。对上式进行一些代数变换：

d’[i] + w’[i][j] = d’[j]

d’[i] + w[i][j] + h[i] – h[j] = d’[j]

(d’[j] + h[j]) - (d’[i] + h[i]) – w[i][j] = 0

(d’[j] + h[j]) - (d’[i] + h[i]) + w[j][i] = 0（因为是费用流，所以有w[i][j] = –w[j][i]）

因此让所有h[i] += d’[i]后，新加入的边(j, i)也会满足势函数的性质。

同时我们有：

d’[i] + w’[i][j] – d’[j] >= 0

d’[i] + h[i] – h[j] + w[i][j] – d’[j] >= 0

(d’[i] + h[i]) - (d’[j] + h[j]) + w[i][j] >= 0

因此修改h[i]后，(i, j)依然会满足势函数的性质。

算法过程如下：

S1 初始化h[]

S2 在残留网络中做Dijkstra

S3 若S到T有可行路径，则修改增广路上的边的容量并所有h[i] += d’[i]，转S2，否则退出

算法时间复杂度为：O(spfa() + K * Dijkstra())或O(K * Dijkstra())，这取决于初始化h[]时是否调用spfa，K表示增广次数。

#include <iostream>
#include <cstdio>
#include <cmath>
#include <string>
#include <cstring>
#include <algorithm>
#include <limits>
#include <vector>
#include <stack>
#include <queue>
#include <set>
#include <map>
#define lowbit(x) ( x&(-x) )
#define pi 3.141592653589793
#define e 2.718281828459045
#define INF 0x3f3f3f3f
#define HalF (l + r)>>1
#define lsn rt<<1
#define rsn rt<<1|1
#define Lson lsn, l, mid
#define Rson rsn, mid+1, r
#define QL Lson, ql, qr
#define QR Rson, ql, qr
#define myself rt, l, r
#define MP(x, y) make_pair(x, y)
using namespace std;
typedef unsigned long long ull;
typedef long long ll;
const int maxN = 4e3 + 7, ss = 0;
int N, K, a[maxN], head[maxN], cnt, S, tt, T, dist[maxN], h[maxN];
struct Eddge
{
    int nex, u, v, flow, cost;
    Eddge(int a=-1, int b=0, int c=0, int d=0, int f=0):nex(a), u(b), v(c), flow(d), cost(f) {}
};
vector<Eddge> G[maxN];
inline void _add(int u, int v, int flow, int cost)
{
    G[u].push_back(Eddge((int)G[v].size(), u, v, flow, cost));
    G[v].push_back(Eddge((int)G[u].size() - 1, v, u, 0, -cost));
}
struct node
{
    int id, val;
    node(int a=0, int b=0):id(a), val(b) {}
    friend bool operator < (node e1, node e2) { return e1.val > e2.val; }
};
int preP[maxN], preE[maxN];
priority_queue<node> Q;
inline int MaxFlow_MinCost(int Flow)
{
    int ans = 0;
    for(int i=0; i<=T; i++) h[i] = 0;
    while(Flow)
    {
        while(!Q.empty()) Q.pop();
        for(int i=0; i<=T; i++) dist[i] = INF;
        dist[ss] = 0;
        Q.push(node(ss, 0));
        while(!Q.empty())
        {
            node now = Q.top(); Q.pop();
            int  u = now.id;
            if(dist[u] < now.val) continue;
            int len = (int)G[u].size();
            for(int i=0, v, f, c; i<len; i++)
            {
                v = G[u][i].v; f = G[u][i].flow; c = G[u][i].cost;
                if(G[u][i].flow && dist[v] > dist[u] + c + h[u] - h[v])
                {
                    dist[v] = dist[u] + c + h[u] - h[v];
                    preP[v] = u; preE[v] = i;
                    Q.push(node(v, dist[v]));
                }
            }
        }
        if(dist[T] == INF) break;
        for(int i=0; i<=T; i++) h[i] += dist[i];
        int Capa = Flow;
        for(int u = T; u; u = preP[u]) Capa = min(Capa, G[preP[u]][preE[u]].flow);
        Flow -= Capa;
        ans += Capa * h[T];
        for(int u = T; u; u = preP[u])
        {
            Eddge &E = G[preP[u]][preE[u]];
            E.flow -= Capa;
            G[E.v][E.nex].flow += Capa;
        }
    }
    return -ans;
}
inline void init()
{
    cnt = 0;
    S = (N << 1) + 1;
    tt = S + 1;
    T = tt + 1;
    for(int i=0; i<=T; i++)
    {
        head[i] = -1;
        h[i] = 0;
        G[i].clear();
    }
    _add(ss, S, K, 0);
    _add(tt, T, K, 0);
}
int main()
{
    int Cas; scanf("%d", &Cas);
    while(Cas--)
    {
        scanf("%d%d", &N, &K);
        init();
        for(int i=1; i<=N; i++) scanf("%d", &a[i]);
        for(int i=1; i<=N; i++)
        {
            _add(S, i, 1, 0);
            _add(i, i + N, 1, -a[i]);
            _add(i + N, tt, 1, 0);
            for(int j = i + 1; j<=N; j++)
            {
                if(a[j] >= a[i])
                {
                    _add(i + N, j, 1, 0);
                }
            }
        }
        printf("%d\n", MaxFlow_MinCost(K));
    }
    return 0;
}