#include #include #include #include #include #include template struct mccf_graph{ int n; // input edges struct edge_with_function{ int from, to; std::function F; // must be discrete convex Flow cap, flow; }; // edges on the residual graph struct edge{ int to, rev, id; bool available; Cost cost; bool is_rev; }; std::vector> g; std::vector E; std::vector dual; mccf_graph() {} mccf_graph(int n) : n(n), g(n){ dual.resize(n, 0); } int add_edge(int from, int to, std::function F, Flow cap, Flow flow){ assert(0 <= from and from < n); assert(0 <= to and to < n); // assert(cap > 0); int m = int(E.size()); E.push_back(edge_with_function{from, to, F, cap, flow}); return m; } Cost calc_objective(){ Cost obj = 0; for(auto &e: E){ obj += e.F(e.flow); } return obj; } /* for debug void print_g(){ for(int i=0;i b){ // b is copied, not referenced assert((int)b.size() == n); std::vector dist(n); std::vector pv(n), pe(n); std::vector vis(n); Flow cap_max = 0; for(const auto &e: E) cap_max = std::max(cap_max, e.cap); Flow delta = 1; while(1){ if(delta * 2 > cap_max) break; delta *= 2; } Flow init_delta = delta; // construct a delta-residual graph int m = E.size(); for(int i=0;i e.cap){ g[e.from].push_back(edge{e.to, (int)g[e.to].size(), i, false, 0, false}); }else{ g[e.from].push_back(edge{e.to, (int)g[e.to].size(), i, true, e.F(delta) - e.F(0), false}); } g[e.to].push_back(edge{e.from, (int)g[e.from].size()-1, i, false, 0, true}); } auto check_RCO = [&](int e_id, Flow x, Flow delta)->bool{ // check edge-wise reduced cost optimality auto &e = E[e_id]; if(x < 0 or x > e.cap) return false; bool ok = true; if(x + delta <= e.cap){ Cost reduced_cost = (init_delta/delta)*(e.F(x+delta) - e.F(x)) - dual[e.to] + dual[e.from]; if(reduced_cost < 0) ok = false; } if(x - delta >= 0){ Cost reduced_cost = (init_delta/delta)*(e.F(x-delta) - e.F(x)) - dual[e.from] + dual[e.to]; if(reduced_cost < 0) ok = false; } return ok; }; auto delta_update_edge = [&](edge &e, Flow delta){ auto &e_func = E[e.id]; if(e.is_rev){ if(e_func.flow - delta < 0){ e.available = false; e.cost = 0; }else{ e.available = true; e.cost = e_func.F(e_func.flow-delta) - e_func.F(e_func.flow); e.cost *= (init_delta/delta); } }else{ if(e_func.flow + delta > e_func.cap){ e.available = false; e.cost = 0; }else{ e.available = true; e.cost = e_func.F(e_func.flow+delta) - e_func.F(e_func.flow); e.cost *= (init_delta/delta); } } return; }; auto saturate = [&](Flow delta){ // modify edge_with_functions for(int i=0;i::max()); std::fill(pv.begin(), pv.end(), -1); std::fill(pe.begin(), pe.end(), -1); std::fill(vis.begin(), vis.end(), false); struct Q { Cost key; int to; bool operator<(Q r) const { return key > r.key; } }; std::priority_queue que; for(int i=0;i=delta){ que.push(Q{0, i}); dist[i] = 0; } } int t = -1; while(!que.empty()){ int v = que.top().to; que.pop(); if(vis[v]) continue; vis[v] = true; if(b[v] <= -delta){ t = v; break; } for(int i=0;i<(int)g[v].size();i++){ auto &e = g[v][i]; if(vis[e.to] or !e.available) continue; Cost cost = e.cost - dual[e.to] + dual[v]; if(dist[e.to] - dist[v] > cost){ dist[e.to] = dist[v] + cost; pv[e.to] = v; pe[e.to] = i; que.push(Q{dist[e.to], e.to}); } } } if(t == -1) return -1; for(int v=0;v= delta) s = pv[v]; } assert(s != -1); b[s] -= delta; b[t] += delta; }; while(delta>0){ saturate(delta); while(true){ int t = delta_scaling_dual(delta); if(t == -1) break; delta_scaling_primal(delta, t); } delta /= 2; } for(int i=0;i> N >> M; std::vector A(N), B(N); for(int i=0;i> A[i]; for(int i=0;i> B[i]; std::vector> T(N, std::vector(N, 0)); for(int i=0;i> T[i][j]; } } // assert int A_sum = 0, B_sum = 0; for(int i=0;i= 0); A_sum += A[i]; } for(int i=0;i= 0); B_sum += B[i]; } assert(A_sum == M); assert(B_sum == M); // graph construction mccf_graph G(2*N); for(int i=0;i C(2*N); for(int i=0;i> res(N, std::vector(N, 0)); for(auto &e: G.E){ res[e.from][e.to-N] = e.flow; } // for(int i=0;i