

Thank you jbernadas, your comment helped me a lot.
I still have question: in tutorial in procedure Reduce Cost (pi) there writes that
c_rev[i, j]
should always, for each edge (I suppose that E is whole set), be 0. I understand that c[i, j] will be 0 if edge (i, j) is part of the shortest part (pi[j] = pi[i] + c[i, j]), but I don't see why c_rev should be 0 for each edge.
So, why I am asking/writing that  in second part, figure 4, graph (g, h, i), edge (3, 2) has cost 6 and it is c_rev[2, 3] (please correct me if I am wrong). In graph (a) there is no edge (3, 2), so (3, 2) could be only c_rev[2, 3] (as result of residual network) and according to Reduce Cost, its cost should be 0, but as I see it's not always after Reduce Cost. Probably I don't understand something, but please show me that part. Am I missing something? Thank you. 

Hello, I will try to explain this trick.
First of all in Reduce_Cost (pi) there is not E but E_x (the set of reduced edges). That is why edge (3,2) is in E_x. Now about my trick: if the edge (i,j) is the part of the shortest path then the reduced costs c(i,j) and c_rev(i,j) will be zero (it is clear and you have understood it). But if there are no negativecost cycles, reduced costs will always be nonnegative. That is why I assign c_rev(i,j)=0 in both cases (it could not be less than zero and if reduced c(i,j)>0 then the edge rev(i,j) is not in E_x). Actually I would have to write c_rev(i,j)=c(i,j). But if there are no multiple edges in the input, this trick will always work. If there are multiple edges or both (i,j) and (j,i) are given in the input, you will have to make 4 edges: (i,j), rev(i,j), (j,i) and rev(j,i) because all of them could have different costs. But I used to solving problems with no multiple edges :)
Now about that program you have shown me. In that program authors use adjacency list. So, the code
for( int i = 0; i < n; i++ )
for( int j = 0; j < n; j++ )
if( cap[i][j]  cap[j][i] ) adj[i][deg[i]++] = j;
creates 4 edges as I said (if both (i,j) and (j,i) are given). Authors don't use something like c_rev(i,j)=0 because they don't keep the reduced costs. They reduce costs "in process" with help of Pot(u,v) macro. The implementation you have shown me is better than mine, but it seems to me little more complex.
As for my mistype with t instead of i  you are absolutely right! Thank you. 

This problem makes me crazy :
First of all in Reduce_Cost (pi) there is not E but E_x (the set of reduced edges)
Sorry, my (printer) mistake  on my printed version there wasn't _x. As I see, it is still whole set of edges (c and c_rev), but not only edges on the shortest path  according to "What is Ex?..." and Figure 4 (g) ((3, 2) is not part of the shortest path, but its cost is changed). Am I correct?
I tried to change code according to my understanding of all yours explanations, but I am getting WA on http://icpcres.ecs.baylor.edu/onlinejudge/ 10594  Data Flow after 0.76 sec (with before written code it AC after 0.8 sec).
Please, can you provide code for Successive Shortest Path with potentials method that uses ReduceCost and really change costs (I suppose you have it written)? I think it also will help to others (as already posted codes help) to understand. I was talking with 3 my friends who are also doing TC and they also have trouble to code it and make to work.
As for my mistype with t instead of i  you are absolutely right! Thank you. You are welcome. Thank you for all your patience to help me. 

Don't worry about this problem, You will understand it, I believe in you! :)
Don't blame your printer. It is my mistake, I wrote 'E_x', but this 'x' became invisible :( It is not whole set of edges; it is the set of edges with positive residual capacities. On figure 4(g) cost of edge (3,2) has been changed because it was 0, but pi[3]=6 and pi[2]=0. So, c(3,2)=0+60=6. This picture is correct.
I have one realization of this method. But I have written it many years ago and this implementation is not very efficient. It uses adjacency matrix and doesnâ€™t allow using multiedges or undirected edges (as you need in 10594 Data Flow task). But this program is rather clear to understand. It really works because I have passed a lot of tasks using this code:
/**
* 
*  MinCostFlow implementation 
*  Successive Shortest Path 
* 
*
* WARNINGS:
* (1) No multiedges are allowed.
* (2) Only one edge (i,j) or (j,i) could be
* between nodes i and j.
* (3) This program uses adjacency matrix
* (could be slow).
*
*/
#include <fstream>
//
#define N 256 // Number of nodes
#define oo 1000000000 // Infinity
//
using namespace std;
ifstream cin("file.in"); // Read from
ofstream cout("file.out"); // Write to
ofstream cerr("file.err"); // Write errors to
int n, m, source, sink;
int G[N][N], // Graph matrix
C[N][N], // Reduced costs matrix
C_UsedToBe[N][N], // Given costs (not reduced)
F[N][N]; // Flow matrix
// For Dijkstra's algorithm:
int inq[N], // (Flag) Is node in the queue?
d[N], // Distance label
pi[N]; // Previous vertex at the shortest path
// It extracts the nearest node
// from the Queue
int extr_min() {
int i, j(0);
for(i=1; i<=n;i++)
if(inq[i]>0&&d[i]<d[j])
j=i;
// Extract the node from the queue
inq[j]=0;
return j;
}
// It reduces costs
void reduce_cost() {
int i, j;
for(i=1; i<=n; i++)
for(j=1; j<=n; j++)
// If (i,j) is the residual edge
// (It's residual capacity is positive)
if(G[i][j]>0) {
C[i][j]+=d[i]d[j];
C[j][i]=0;
}
}
// It finds the shortest path
// using Dijkstra's algorithm
// and augments a flow along it.
int augment() {
int i, j;
for(i=0; i<=n; i++) {
d[i]=oo; // Nodes are unreachable
inq[i]=1; // All nodes are in the queue
}
d[source]=0;
// While queue is not empty
for( ; (i=extr_min())>0; ) {
for(j=1; j<=n; j++) {
// If there is no residual edge (i,j) or
// node j is not in the queue
if(G[i][j] == 0  inq[j] == 0)
continue;
// If the distance label of j
// bigger than the distance label
// of i plus (residual) distance from i to j
if(d[j]>d[i]+C[i][j]) {
pi[j]=i; // Make i previous for j
d[j]=d[i]+C[i][j]; // Relabel
}
}
}
// Reduce costs
reduce_cost();
// If there are no shortest (and any)
// augmenting paths
if(inq[sink]==1)
return 0;
// Calculating the "bottleneck" of the path
int width(oo);
for(j=sink; j!=source; j=i) {
i=pi[j];
if(width>G[i][j])
width=G[i][j];
}
// Make the augmentation
for(j=sink; j!=source; j=i) {
i=pi[j];
G[i][j]=width; F[i][j]+=width;
G[j][i]+=width; F[j][i]=width;
}
return width;
}
int find_max_flow_min_cost() {
int ans(0), i, j;
// While there is any shortest
// augmenting path, augment flow
// along it.
while(augment());
// No calculate the total cost
// of the maximum flow
for(i=1; i<=n; i++) {
for(j=1; j<=n; j++)
if(F[i][j]>0)
ans+=F[i][j]*C_UsedToBe[i][j];
}
return ans;
}
int main() {
int i, p, q, u, c;
// Read #Nodes, #Edges, Source node and Sink node.
cin >> n >> m >> source >> sink;
// Read edges
for(i=1; i<=m; i++) {
// "from" "to" "capacity" "cost"
cin >> p >> q >> u >> c;
// If there is a multiedge
if(G[p][q]>0) {
cerr << "Multiedge!" << endl;
return 0;
}
G[p][q]=u;
C[p][q]=C_UsedToBe[p][q]=c;
}
// Output the cost of the minimumcost flow
cout << find_max_flow_min_cost() << endl;
return 0;
}
Input format is very simple: #Nodes #Edges Source Sink from_1 to_1 capacity_1 cost_1 from_2 to_2 capacity_2 cost_2 ... ... ... ... from_m to_m capacity_m cost_m
I wish you good luck! 

Thank you very much. I think I've got it completely now  thanks to code and your explanations. I will try to prove my understanding solving some tasks. 

I too I'm trying to solve the 10594 problem in UVa. I still haven't figured out how negative costs came in to play in this problem... can someone give me some pointers.
My current solution is using a implementation of FordFulkerson with Dijkstra for path finding and using a cost of (time*minflowinpath). This is maybe too greedy, because sometimes choosing the shorter paths dont maximize the flow. 

Btw, here is nothing told about special meaning of these potentials. I've met the following fact: these potentials are the solution of a linear programming problem, dual to a usual mincostmaxflow problem. (the article of Sokkalingam, Ahuja and Orlin "INVERSE SPANNING TREE PROBLEMS: FORMULATIONS AND ALGORITHMS"; they showed that inverseMST problem can be expressed as a dual problem to an assignment problem).
But in the article there is no proof of the fact (which is, I think, very interesting and useful). They refer to their book "Network Flows: Theory, Algorithms, and Applications", but I haven't found it in the Internet. And, what is very strange, I didn't see the fact in any other book... 

Here is my implementation:
#include <algorithm>
#include <iostream>
#include <sstream>
#include <string>
#include <utility>
#include <vector>
#include <cmath>
#include <queue>
#include <map>
using namespace std;
bool bellman_ford(const std::vector<std::vector<int> >& graph, const std::vector<std::vector<int> >& enabled, int start, std::vector<int>& distances)
{
int N = graph.size();
distances.assign(N, INT_MAX>>1);
distances[start] = 0;
for(int s = 0; s < N1; ++s)
{
for(int i = 0; i < N; ++i)
for(int j = 0; j < N; ++j)
if (enabled[i][j])
distances[j] = std::min<int>(distances[j], distances[i] + graph[i][j]);
}
for(int i = 0; i < N; ++i)
{
for(int j = 0; j < N; ++j)
{
if (enabled[i][j] && distances[j] > distances[i] + graph[i][j])
return false;
}
}
return true;
}
std::pair<std::vector<int>, std::vector<int> > dijkstra(const std::vector<std::vector<int> >& cost, const std::vector<std::vector<int> >& enabled, const std::vector<int>& potentials, int start)
{
int N = cost.size();
std::vector<int> distances(N, INT_MAX>>1);
std::vector<int> previous(N, 1);
std::priority_queue<std::pair<int,int>, std::vector<std::pair<int,int> >, std::greater<std::pair<int,int> > > q;
q.push(std::make_pair(0, start) );
while(!q.empty())
{
int v = q.top().second, d = q.top().first; q.pop();
distances[v] = d;
for (int i = 0; i < N; ++i)
{
if (!enabled[v][i])
continue;
int new_d = d+potentials[v]potentials[i] + cost[v][i];
if (distances[i] > new_d)
{
q.push(std::make_pair(new_d, i));
previous[i] = v;
}
}
}
return std::make_pair(distances, previous);
}
std::pair<int,int> min_cost_max_flow(const std::vector<std::vector<int> >& capacities, const std::vector<std::vector<int> >& flow_cost, int source, int sink)
{
std::vector<std::vector<int> > residual(capacities);
std::vector<std::vector<int> > cost(flow_cost);
int N = residual.size();
for(int i = 0; i < N; ++i)
for(int j = 0; j < N; ++j)
if (capacities[i][j] > 0)
cost[j][i] = cost[i][j];
std::vector<int> potentials;
while(!bellman_ford(cost, residual, source, potentials))
{
//fill up the negative cycles
}
while(true)
{
std::pair<std::vector<int>, std::vector<int> > dijkstra_result = dijkstra(cost, residual, potentials, source);
const std::vector<int>& distances = dijkstra_result.first;
const std::vector<int>& previous = dijkstra_result.second;
potentials = distances;
if (previous[sink] == 1)
break;
std::vector<int> path;
int v=sink;
while(true)
{
if (v==1)
break;
path.push_back(v);
v = previous[v];
}
std::reverse(path.begin(), path.end());
int flow = INT_MAX>>1;
for(int i = 0; i < path.size()1; ++i)
flow = std::min<int>(flow, residual[path[i]][path[i+1]]);
for(int i = 0; i < path.size()1; ++i)
{
residual[path[i]][path[i+1]] = flow;
residual[path[i+1]][path[i]] += flow;
}
}
int flow = 0;
for(int i = 0; i < N; ++i)
flow += capacities[source][i]  residual[source][i];
int cost_flow = 0;
for(int i = 0; i < N; ++i)
{
for(int j = 0; j < N; ++j)
{
if (capacities[i][j])
cost_flow += (capacities[i][j]  residual[i][j])*cost[i][j];
}
}
return std::make_pair(flow, cost_flow);
}
int main(int argc, char* argv[])
{
freopen("in.txt","r", stdin);
int n,m,source,sink;
cin >> n >> m >> source >> sink;
vector<vector<int> > capacities(n, vector<int>(n));
vector<vector<int> > costs(n, vector<int>(n));
for(int i=0; i<m; i++)
{
int from, to, capacity, cost;
cin >> from >> to >> capacity >> cost;
capacities[from][to] = capacity;
costs[from][to] = cost;
}
pair<int,int> res = min_cost_max_flow(capacities, costs, source, sink);
cout << "flow = " << res.first << ", cost = " << res.second << "\n";
}


I have go through the explanation of the trick:
1. c(i,j) is on the shortest path  then c(i,j) = c(j,i) = 0
2. c(i,j) is not on the shortest path  then c(j,i) can be assigned to 0 because there will be no positive residual capacity along the arc c(j,i) and it flows from the fact that there are no negative cost cycles Am I correct? And what is the purpose of this trick? Why cannot we set c[j][i] = c[i][j] at the beginning? 

You're right, also I suggest you to read this post. But I also have some trouble with understanding potentials. There is an illustration from the Ahuja/Magnanti/Orlin book. At the figure d the potential of the node "3" changes to 3. Why is the cost of the edge 31 is still 2? Why is it not 2+(3)0? This edge is present in the residual network, and has nonzero residual capacity.
Potentials are negative because in Ahuja's notation cpi(i, j) = c(i, j) + pi(j)  pi(i), so pi = d. 

A good explanation of potentials is in wiki: Johnson's algorithm. 

# include <iostream> # include <stdio.h> # include <string.h> # include <string> # include <utility> # include <vector> # include <algorithm> # include <queue> # include <map> # include <bitset> # include <math.h> # include <limits.h> using namespace std;
// macros #define fr(i, st, ed) for (int i = st ; i < ed; ++i ) #define ipInt(x) scanf("%d", &x) #define MAX_NODES 110 #define inf LLONG_MAX
//typedefinitions typedef long long int lli; typedef pair<int, int> pii; typedef vector<int> vi; typedef vector<lli> vlli; typedef vector<pii> vii; typedef vector<vii> matrix; typedef map<string, int> mSI; typedef pair<int , pii> edge; //typedefinitions matrix adjMat(MAX_NODES); int n, m;
// MAxFlow variables lli minCost, cost[MAX_NODES][MAX_NODES], costToAdd, flowToIncrease, maxFlow, flow[MAX_NODES][MAX_NODES], cap[MAX_NODES][MAX_NODES]; vi parentBellmen; vlli dist; // MAxFlow variables
// problem specific variables lli flowToGet, k; // problem specific variables
lli resCapacity (int from, int to) { // cap[from][to]  flow[from][to]; return k  flow[from][to]; }
void makeEdge (int from, int to, lli weight) { // because here apacity is constant in this question if not constant // uncomment these lines and give the necessary input to the function // cap[from][to] = weight; // cap[to][from] = weight; cost[from][to] = cost[to][from] = weight; adjMat[from].push_back(pii(to, 0)); adjMat[to].push_back(pii(from, 0)); } void augmentFlow (int source, int sink) { if (parentBellmen[sink] == 1) return; lli minResCap = inf; int node = sink; while (parentBellmen[node] != 1) { int parent = parentBellmen[node]; minResCap = min(minResCap, resCapacity(parent, node)); node = parent; } node = sink; while (parentBellmen[node] != 1) { int parent = parentBellmen[node]; flow[parent][node] = flow[parent][node] + minResCap; flow[node][parent] = flow[node][parent]  minResCap; node = parent; } flowToIncrease = minResCap; costToAdd = dist[sink]; }
//djakstra void SPFA (int source) { //initialize dist queue<int> q; dist.assign(MAX_NODES, inf); vi visited(MAX_NODES, 0); parentBellmen.assign(MAX_NODES, 1); //initialization end //start dist[source] = 0; visited[source] = 1; q.push(source); while (!q.empty()) { int u = q.front(); q.pop(); visited[u] = 0; int neighbours = adjMat[u].size(); fr (i, 0, neighbours) { int v = adjMat[u][i].first; lli costVal = (flow[v][u] > 0) ? cost[u][v] : cost[u][v]; if (resCapacity(u, v) && dist[u] + costVal < dist[v]) { parentBellmen[v] = u; dist[v] = dist[u] + costVal; if (!visited[v]) { visited[v] = 1; q.push(v); } } }
} }
void bellmenFord (int source) { //initialize dist queue<int> q; dist.assign(MAX_NODES, inf); parentBellmen.assign(MAX_NODES, 1); //initialization ends
//start dist[source] = 0; fr (turn, 0, n  1){ fr (v, 0, n) { if (dist[v] != inf) { int neighbours = adjMat[v].size(); // going on every edge of v fr (i, 0, neighbours) { int edgeToNode = adjMat[v][i].first; lli edgeWeight = (flow[edgeToNode][v] > 0) ? cost[v][edgeToNode] : cost[v][edgeToNode]; if ((resCapacity(v, edgeToNode) != 0) && (dist[v] + edgeWeight < dist[edgeToNode])) { dist[edgeToNode] = dist[v] + edgeWeight; parentBellmen[edgeToNode] = v; } } } } } }
// this now becomes minCost Max Flow void runMaxFlow (int s, int t) { // while i can increase / augment flow i will increase the maxFlow value while (true) { flowToIncrease = 0; costToAdd = 0; SPFA(s); augmentFlow(s, t); if (flowToIncrease == 0) { break; } if (maxFlow + flowToIncrease >= flowToGet) { lli num = costToAdd * (flowToGet  maxFlow); minCost += num; maxFlow = flowToGet; break; } maxFlow = maxFlow + flowToIncrease; minCost += costToAdd * flowToIncrease; } } void initialize () { memset(cap, 0, sizeof cap); memset(flow, 0, sizeof flow); fr(i, 0, MAX_NODES) fr(j, 0, MAX_NODES) cost[i][j] = inf; fr(i, 0, adjMat.size()) { adjMat[i].clear(); } maxFlow = 0; minCost = 0; }
int main () { int s, t, c; while (scanf("%d %d", &n, &m) != EOF) { initialize(); while(m) { int from, to; lli timeVal; ipInt(from);ipInt(to);scanf("%lld", &timeVal); makeEdge(from  1 , to  1 , timeVal); } scanf("%lld %lld", &flowToGet, &k); s = 0; t = n  1; runMaxFlow(s, t); if (maxFlow == flowToGet) { printf("%lld\n", minCost); } else { printf("Impossible.\n"); } } return 0; }
// In this question capacity of each node was same i.e why to get the cost we just needed to do dist[sink] * flowToIncrease; because if i empty the flow // in a running stream in this question I only empty the stream because every edge has same capacity and when i emptied the flow, the flow also had // the same capacity, so whenever i empty the flow i make the flow empty.. // But if capacities are different then What ? // in that case costToIncrease != dist[sink] * flowToIncrease.... // i have to increase the constToIncrease while augmenting like this // costToIncrease += cost[parent[node]][node] * minResCap;
// if reverting the flow then canceling the cost which was added earlier // lli costVal = (flow[node][parent] > 0) ? cost[parent][node] : cost[parent][node]; // costToIncrease += flowToIncrease * costVal; 
