Pre knowledge

• Maximum flow (refer to blog: https://www.cnblogs.com/konjakhzx/p/15561745.html)

It is to find the minimum cost in the case of maximum flow (and the maximum cost maximum flow can take the cost as a negative number to run the minimum cost maximum flow).
So we thought of an idea. When looking for the expansion path, we should change the shortest expansion path to the one with the lowest cost.
The correct way of thinking is like this.
Then there are two ideas

ZKW expense flow

Invented by Zhang kunwei, gold medalist of noi 2008

solution

First, let's set a dis value, which is a bit like the d value in sap, but here dis[i] represents the shortest path from the end point to point I.
Then we can directly find the shortest path according to the dis value. For example, i want to know whether it is the shortest path from point i to point j. when dis[j] + the distance from point i to point j = = dis[j], the shortest path must be taken here, otherwise it is not.
When you can't find the augmented path, there is a problem with your dis value, so we consider how to change the dis value.

We first divide the node into two parts, one is the S part and the other is the T part.
Part S is the point that can be reached by dis value from the source point, and the remaining points are part T.
First, we know that if we add all the dis values of the points of the S part to a value, the properties of the S part are satisfied.
Therefore, we only need to find a value so that the S part can be connected with the T part after adding this value.
Because we seek the shortest path, we find a point from part S to part T, which is the minimum value that dis needs to increase.

Then we run the maximum flow until we can'T get from the source point to the sink point, that is, the point of part S can never be connected with the point of part T.

How to record the answer is to multiply the amount of your expansion by the sum of the expenses you have passed along the way, that is, the shortest path from the source point to the sink point by the expansion.

code

bool label()
{
    int mindis=INF;
    for(int i=1;i<=n;i++)
        if(vis[i])
            for(int j=last[i];j;j=nxt[j])
                if(!vis[to[j]]&&a[j])
                    mindis=min(w[j]+dis[to[j]]-dis[i],mindis);//Find the minimum difference
    if(mindis==INF)
        return false;//If no edge can flow from part S to part T, it means there is no widening path
    for(int i=1;i<=n;i++)
        if(vis[i])
            dis[i]+=mindis;//Add a value to all dis of points in part S
    return true;
}
int aug(int s,int augcc)
{
    if(s==n)
    {
        ans+=dis[1]*augcc;//Multiply the expansion by the sum of the costs along the way is the price you spend
        return augcc;
    }
    int augc=augcc,delta,son;
    vis[s]=1;
    for(int i=last[s];i&&augc;i=nxt[i])
    {
        son=to[i];
        if(a[i]>0&&!vis[son]&&dis[s]==dis[son]+w[i])
        {
            delta=aug(son,min(a[i],augc));
            a[i]-=delta,a[i^1]+=delta;//i^1 is the reverse side, and the side is stored from side 2
            augc-=delta;
        }
    }
    return augcc-augc;
}
void augca()
{
    do
    {
        do
        {
            memset(vis,0,sizeof(vis));
        }while(aug(1,INF));
    }while(label());//Update the dis value when the augmentation path cannot be found
}

Find the shortest path with SPFA

solution

Before finding the augmentation path each time, use SPFA to find the shortest path, and then use the shortest path for augmentation.
We use dis value to record the shortest path, and then expand it according to zkw cost flow expansion method.
In other words, every time it is illegal, we do not change the dis value based on the original dis value, but recalculate the dis value.

code

bool spfa()
{
    queue<int> myq;
    myq.push(target);
    memset(dis,0x3f,sizeof dis);
    dis[target]=0,inq[target]=1;
    while(!myq.empty())
    {
        int t=myq.front();
        myq.pop(),inq[t]=0;
        for(int i=fir[t];i;i=e[i].next)
        {
            int v=e[i].v;
            if(e[i^1].cap>0)
                if(dis[v]>dis[t]+e[i^1].cost)
                {
                    dis[v]=dis[t]+e[i^1].cost;
                    if(inq[v]==0) myq.push(v),inq[v]=1;
                }
        }
    }
    if(dis[src]==0x3f3f3f3f)return 0;
    else return 1;
}
int aug(int s,int augco)
{
    if(s==target)
    {
        ans+=dis[src]*augco;
        return augco;
    }
    int tempcap=augco,delta;
    for(int i=fir[s];i&&tempcap>0;i=e[i].next)
    {
        if(!vis[e[i].v]&&e[i].cap>0)
        {
            if(dis[e[i].v]==dis[s]-e[i].cost)
            {
                vis[e[i].v]=1;
                delta=aug(e[i].v,min(e[i].cap,tempcap));
                vis[e[i].v]=0;
                tempcap-=delta;
                e[i].cap-=delta;
                e[i^1].cap+=delta;
            }
        }
    }
    return augco-tempcap;
}