BASH - Editorial

PROBLEM LINK:

Practice
Div-2 Contest
Div-1 Contest

Author: Sayantan Jana
Tester: Felipe Mota
Editorialist: Sayantan Jana

DIFFICULTY:

MEDIUM-HARD

PREREQUISITES:

Directed Minimum Spanning Tree, Flows

PROBLEM:

A Bash Matrix is a N*N matrix containing characters ‘U’,‘D’,‘L’,‘R’ on each cell, which act as instructions which cell to move next to, from current cell. Suppose one starts at a cell, he keeps on moving as per the directions of the current cell he is in and finally stops at the cell he had already visited.

Given the information of the stopping cell for each of the N^2 cells one can start movement at, find a Bash Matrix with the smallest possible cost. Putting characters ‘U’, ‘L’, ‘D’ or ‘R’ in a Bash matrix costs P_U, P_L, P_D or P_R per character respectively.

QUICK EXPLANATION:

  • Cyclic Nodes : For all the cells (i,j) such that corresponding stopping cell (s_i,s_j) = (i,j), should be in some cycles .
  • Non-cyclic Nodes : For all the cells (i,j) such that corresponding stopping cell (s_i,s_j) \neq (i,j), should not be in any cycle.
  • The cyclic nodes need to be put in some cycles. Notice that there can only be even length cycles. Now these even length cycles can always be broken into multiple cycles of length 2. It can be shown that the most effective way (minimum cost) of putting these cells into cycles necessarily needs to have the same total cost as trying to effectively (minimum cost) putting these nodes into cycles of length 2.
  • This can be solved by bipartitie matching. Create a bi-partition, a partition of those cells (i,j) such that (i+j) mod 2 == 0 and another partition of those cells (i,j) such that (i+j) mod 2 == 1. Notice there cannot be edges within a partition and always across a partition. Hence, for putting the cyclic nodes into cycles of length 2, we need to select a set of bidirectional edges across the partitions, which covers each of the vertices exactly once (perfect matching), which can be solved by any minimum cost bipartite matching algorithm .
  • For the non cyclic nodes, we need to put directions on the cells in a way that lead up to their terminals (stopping cells) and the total cost of doing this needs to be minimised. The terminals necessarily need to be cyclic nodes and hence once a directed path from a cell to its terminal is built, it is guaranteed to end up on a cycle and have the terminal as its first previously visited cell. This can be solved by directed minimum spanning tree. Check Edmond’s algorithm.
  • If it is not possible to put the cyclic nodes into cycles or if it is not possible to put directions in a way to that lead each cell to its terminal, the answer is -1.

EXPLANATION:

First let us discuss about the existence of solution, then we will move onto how to solve each of the parts as described in the Quick Explanation section.

Conditions for solution to exist

  • Cyclic Nodes : For all the cells (i,j) such that corresponding stopping cell (s_i,s_j) = (i,j), should be in some cycles. If it is not possible to put directions on the cyclic nodes in a way that all of them are in cycles, a solution cannot exist. Note that this is independent of the directions on the non cyclic nodes, since every node that is visited while starting at a cyclic node, needs to be itself in a cycle, i.e., it should be a self terminal [stopping cell (s_i,s_j) = (i,j)].
  • Non-cyclic Nodes : For all the cells (i,j) such that corresponding stopping cell (s_i,s_j) \neq (i,j), should not be in any cycle. If it is not possible to put the directions on the non cyclic nodes such that it leads to its terminal, a solution cannot exist. Note that the terminal of a non cyclic node must be a cyclic node since it is the first previously visited node for the traversal starting at the non cylic node. Consider a division of the non cylic nodes into several connected components such that each component has the same terminal (i.e., stopping cell). It is necessary that the terminal of a component be a neighbour of one of the non cyclic nodes in the component, otherwise a solution doesn’t exist. Also, notice that there cannot be 2 disconnected components with the same terminal. Putting directions on the non cyclic nodes again is independent of the directions on cyclic nodes.

As we understood in the conditions, putting directions on the cyclic nodes and non cyclic nodes are completely independent. Hence, once we classify which nodes are cyclic and which are non cyclic, we can split the problem into 2 branches.

For subtask 1, it is only required to find a solution without trying to minimise cost since cost for each of the cells is same.

Orientation in cylic nodes

This part only deals with the cyclic nodes, further discussion is only concerned with them in this part.

Let’s define the parity of a cell, parity(i,j) = (i+j) mod 2 .

Form 2 partitions, a partition consisting of even parity cells and another partition consisting of odd parity cells. Notice that all the edges in the matrix are to be formed are supposed to be across the partitions, that is from an even parity node to an odd parity node or from an odd parity node to even parity node.

Try constructing a bipartite graph from the 2 partitions by adding all the possible edges for each of the nodes. Note that this would necessarily be a bipartite graph since there cannot be edge between same parity nodes.

Now since this is a bipartite graph, there cannot exist an odd cycle. All the cycles to be formed are supposed to be even length cycles. Each even length cycle can be broken into multiple chains of 2 length cycles. A 2-length cycle is basically 2 neighbouring nodes which point at each other.

If there exists a way to put the directions on the cyclic nodes, such that each node is in cycle, it’s possible to break all the cycles into cycles of length 2. Hence, there exists a pairing cycle (2-length cycle) solution if any solution is possible. If we want to solve just subtask 1, we can get a solution by running any bipartite matching algorithm on the bipartite graph and check if we can have a perfect matching. This will ensure that each node on each partition has been linked to a node on the other partition, which helps get the pairing solution if a perfect matching exists.

For minimum cost, for each of the edges associate a cost with it. For each vertical edge in the biparite graph , put cost = cost(U)+cost(D) and for each horizontal edge put cost = cost(L)+cost(R). We claim that running a minimum cost perfect matching gives us a solution with minimum cost.

The solution can be proved correct by contradiction, i.e., we assume a possible way of putting directions on cyclic nodes with lower minimum cost. Since in any possible cycle set construction, the number of upward and downward edges are supposed to be equal , similarly the number of left and right edges are supposed to be equal . Also, the individual direction cost doesn’t matter, but cost(U)+cost(D) vs cost(L)+cost(R). Since count(U) = count(D) and count(L) = count(R), the total no. of nodes having an outward vertical edge is even , the total no. of nodes having an outward horizontal edge is also even . Hence there is a possible pairing of the nodes in the alternative solution . Therefore, a contradiction arise that our algorithm that was designated for pairing found a construction with higher cost . Hence our assumption was wrong .

Thus running a minimum cost perfect matching gives us a solution with minimum cost. We can use general minimum cost maximum flow algorithm here, though we could have tried for an optimised solution since our graph at hand has certain special properties. A O(F*(V+E)*log(V)) solution works for us, where F is the amount of the flow, V is the number of vertices (O(N^2)) and E is the number of edges (O(N^2)), and hence overall O(N^4*log(N)) complexity in our case.

Orientation in non-cylic nodes

Recall the division of non-cyclic nodes into components stated in point 2 of the conditions for existence of solution. In each component we need to assign the directions so that each cell lead upto the terminal, that we already mentioned needs to be a neighbour of one of the nodes in the component and hence a solution always exists.

If we want to solve just subtask 1, we can start at the terminal and keep visiting the neighbouring nodes with the same terminal (these are the nodes forming the component), hence pulling an edge from the neighbouring cell to the current cell.

For assigning directions and ensure minimum cost, the above strategy doesn’t work. Our problem at hand is choose a minimum cost directed spanning tree of edges from the component (including the terminal here) such that there is a path from each of the nodes to the terminal. Consider all the edges reversed, then this problem has a change in condition, there should be a path from the terminal to each of the nodes. This is actually a known problem, directed MST. There exists Edmond’s algorithm to solve the problem. Check the link in prerequisites to know more on the algorithm. It is also quite interesting to note how running a Prim’s algorithm starting at the terminal node doesn’t guarantee a minimum cost solution for directed MST. There exists a better complexity for the algorithm but O(V*E) solution works for us, where V is the number of vertices(O(N^2)) and E is the number of edges (O(N^2)), and hence overall O(N^4) complexity in our case.

Overall Complexity : O(N^4*log(N))

SOLUTIONS:

Setter's Solution
#include <bits/stdc++.h>
using namespace std;
 
// minCost maxFlow
 
template <class Cap, class Cost> struct mcf_graph {
  public:
    mcf_graph() {}
    mcf_graph(int n) : _n(n), g(n) {}
 
    int add_edge(int from, int to, Cap cap, Cost cost) {
        assert(0 <= from && from < _n);
        assert(0 <= to && to < _n);
        int m = int(pos.size());
        pos.push_back({from, int(g[from].size())});
        int from_id = int(g[from].size());
        int to_id = int(g[to].size());
        if (from == to) to_id++;
        g[from].push_back(_edge{to, to_id, cap, cost});
        g[to].push_back(_edge{from, from_id, 0, -cost});
        return m;
    }
 
    struct edge {
        int from, to;
        Cap cap, flow;
        Cost cost;
    };
 
    edge get_edge(int i) {
        int m = int(pos.size());
        assert(0 <= i && i < m);
        auto _e = g[pos[i].first][pos[i].second];
        auto _re = g[_e.to][_e.rev];
        return edge{
            pos[i].first, _e.to, _e.cap + _re.cap, _re.cap, _e.cost,
        };
    }
    std::vector<edge> edges() {
        int m = int(pos.size());
        std::vector<edge> result(m);
        for (int i = 0; i < m; i++) {
            result[i] = get_edge(i);
        }
        return result;
    }
 
    std::pair<Cap, Cost> flow(int s, int t) {
        return flow(s, t, std::numeric_limits<Cap>::max());
    }
    std::pair<Cap, Cost> flow(int s, int t, Cap flow_limit) {
        return slope(s, t, flow_limit).back();
    }
    std::vector<std::pair<Cap, Cost>> slope(int s, int t) {
        return slope(s, t, std::numeric_limits<Cap>::max());
    }
    std::vector<std::pair<Cap, Cost>> slope(int s, int t, Cap flow_limit) {
        assert(0 <= s && s < _n);
        assert(0 <= t && t < _n);
        assert(s != t);
        std::vector<Cost> dual(_n, 0), dist(_n);
        std::vector<int> pv(_n), pe(_n);
        std::vector<bool> vis(_n);
        auto dual_ref = [&]() {
            std::fill(dist.begin(), dist.end(),
                      std::numeric_limits<Cost>::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<Q> que;
            dist[s] = 0;
            que.push(Q{0, s});
            while (!que.empty()) {
                int v = que.top().to;
                que.pop();
                if (vis[v]) continue;
                vis[v] = true;
                if (v == t) break;
                for (int i = 0; i < int(g[v].size()); i++) {
                    auto e = g[v][i];
                    if (vis[e.to] || !e.cap) 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 (!vis[t]) {
                return false;
            }
 
            for (int v = 0; v < _n; v++) {
                if (!vis[v]) continue;
                dual[v] -= dist[t] - dist[v];
            }
            return true;
        };
        Cap flow = 0;
        Cost cost = 0, prev_cost_per_flow = -1;
        std::vector<std::pair<Cap, Cost>> result;
        result.push_back({flow, cost});
        while (flow < flow_limit) {
            if (!dual_ref()) break;
            Cap c = flow_limit - flow;
            for (int v = t; v != s; v = pv[v]) {
                c = std::min(c, g[pv[v]][pe[v]].cap);
            }
            for (int v = t; v != s; v = pv[v]) {
                auto& e = g[pv[v]][pe[v]];
                e.cap -= c;
                g[v][e.rev].cap += c;
            }
            Cost d = -dual[s];
            flow += c;
            cost += c * d;
            if (prev_cost_per_flow == d) {
                result.pop_back();
            }
            result.push_back({flow, cost});
            prev_cost_per_flow = d;
        }
        return result;
    }
 
  private:
    int _n;
 
    struct _edge {
        int to, rev;
        Cap cap;
        Cost cost;
    };
 
    std::vector<std::pair<int, int>> pos;
    std::vector<std::vector<_edge>> g;
};
 
// var declarations and helper functions 
 
int n;
int cost;
int cost_up,cost_down,cost_left,cost_right;
vector<int> parity_nodes[2];
int node_rep[51][51];
int x_rep[51*51],y_rep[51*51];
int stop_x[51][51],stop_y[51][51];
vector<string> direction;
bool visited[51][51],in_cycle[51][51];
map<char,int> price;
 
bool is_part_of_cycle(int i, int j)
{
    if(stop_x[i][j] == i && stop_y[i][j] == j)
        return true;
    else
        return false;
}
 
bool has_same_stop(int i, int j, int x, int y)
{
    if(stop_x[i][j] == stop_x[x][y] && stop_y[i][j] == stop_y[x][y])
        return true;
    else
        return false;
}
 
// Directed MST
 
const int N = 100005;
const int E = 1000000;
const int inf = 2139062143;
 
struct edge {
  int u, v, len, id;
}e[N];
int chosen[E];
int pre[N], inwei[N], innum[N], id[N], vis[N];
int inc[E], decr[E], pos;
 
inline int mst(int n, int m, int root) {
  int ans = 0, eg = m;
  pos = m;
  while (true) {
    memset(inwei, 127, sizeof(int) * (n + 2));
    memset(id, 0, sizeof(int) * (n + 2));
    memset(vis, 0, sizeof(int) * (n + 2));
    for (int i = 1; i <= m; i++)
      if (e[i].len < inwei[e[i].v]) {
    inwei[e[i].v] = e[i].len;
    pre[e[i].v] = e[i].u;
    innum[e[i].v] = e[i].id;
      }
    inwei[root] = 0;
    for (int i = 1; i <= n; i++)
      if (inwei[i] == inf) return -1;
    int tot = 0;
    for (int i = 1; i <= n; i++) {
      ans += inwei[i];
      if (i != root) chosen[innum[i]]++;
      int j = i;
      while (vis[j] != i && j != root && !id[j])
    vis[j] = i, j = pre[j];
      if (j != root && !id[j]) {
    id[j] = ++tot;
    for (int k = pre[j]; k != j; k = pre[k]) id[k] = tot;
      }
    }
    if (!tot) break;
    for (int i = 1; i <= n; i++)
      if (!id[i]) id[i] = ++tot;
    n = tot;
    for (int i = 1; i <= m; i++) {
      int v = e[i].v;
      e[i].v = id[e[i].v];
      e[i].u = id[e[i].u];
      if (e[i].v != e[i].u) {
    e[i].len -= inwei[v];
    pos ++;
    inc[pos] = e[i].id;
    decr[pos] = innum[v];
    e[i].id = pos;
      }
      else swap(e[i--], e[m--]);
    }
    root = id[root];
  }
  while (pos > eg) {
    if (chosen[pos] > 0) {
      chosen[inc[pos]]++;
      chosen[decr[pos]]--;
    }
    pos--;
  }
  return ans;
}
bool bad[N];
char orientation[N];
int x_coord[N],y_coord[N];
int edge_cnt = 0;
int vertices = 0;
map<pair<int,int>,int> mp;
 
void put_edge(int s,int t,int i,int j,char d)
{
    edge_cnt++;
    bad[edge_cnt] = 1;
    if(!mp[make_pair(s,t)])
        mp[make_pair(s,t)] = ++vertices;
    if(!mp[make_pair(i,j)])
        mp[make_pair(i,j)] = ++vertices;
    e[edge_cnt].u = mp[make_pair(s,t)];
    e[edge_cnt].v = mp[make_pair(i,j)];
    e[edge_cnt].len = price[d];
    e[edge_cnt].id = edge_cnt;
    x_coord[edge_cnt] = i;
    y_coord[edge_cnt] = j;
    orientation[edge_cnt] = d;
}
 
void form_graph(int i,int j)
{
    if(visited[i][j]) return;
    visited[i][j] = true;
    if(i >= 2 && has_same_stop(i,j,i-1,j))
    {
        put_edge(i,j,i-1,j,'D');
        form_graph(i-1,j);
    }
    if(i <= n-1 && has_same_stop(i,j,i+1,j))
    {
        put_edge(i,j,i+1,j,'U');
        form_graph(i+1,j);
    }
    if(j >= 2 && has_same_stop(i,j,i,j-1))
    {
        put_edge(i,j,i,j-1,'R');
        form_graph(i,j-1);
    }
    if(j <= n-1 && has_same_stop(i,j,i,j+1))
    {
        put_edge(i,j,i,j+1,'L');
        form_graph(i,j+1);
    }
}
 
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout.tie(0);
 
    cin >> n;
    cin >> cost_up >> cost_left >> cost_down >> cost_right;
    price['U'] = cost_up;
    price['L'] = cost_left;
    price['D'] = cost_down;
    price['R'] = cost_right;
    for(int i = 1; i <= n; ++i)
        for(int j = 1; j <= n; ++j)
            cin >> stop_x[i][j] >> stop_y[i][j];
 
    direction.resize(n+1);
    for(int i = 1; i <= n; ++i)
        for(int j =1; j <= n; ++j)
            direction[i].push_back('N');
 
    int node_cnt = 0;
 
    for(int i = 1; i <= n; ++i)
    {
        for(int j = 1; j <= n; ++j)
        {
            if(is_part_of_cycle(i,j))
            {
                node_rep[i][j] == ++node_cnt;
                x_rep[node_cnt] = i;
                y_rep[node_cnt] = j;
                parity_nodes[(i+j)%2].push_back(node_cnt);
                form_graph(i,j);
                int tmp = edge_cnt;
                cost += mst(vertices,edge_cnt,mp[make_pair(i,j)]);
                for(int k = 1; k <= tmp; ++k)
                    if (chosen[k] && bad[k]) direction[x_coord[k]][y_coord[k]] = orientation[k]; 
 
                edge_cnt = pos = vertices = 0;
                mp.clear();
                for(int k = 0; k <= 10000; ++k)
                {
                    x_coord[k] = y_coord[k] = 0;
                    pre[k] = inwei[k] = innum[k] = id[k] = vis[k] = 0;
                    chosen[k] = bad[k] = 0;
                    orientation[k] = 'N';
                    e[k].u = e[k].v = e[k].len = e[k].id = 0;
                    inc[k] = decr[k] = 0;
                }
            }
        }
    }
 
    for(int i = 1; i <= n; ++i)
    {
        for(int j = 1; j <= n; ++j)
        {
            if(!visited[i][j])
            {
                // cout << "didn't get allocated a path leading to cycle" << endl;
                cout << "-1" << "\n";
                return 0;
            }
        }
    }
 
    if((int)parity_nodes[0].size() != (int)parity_nodes[1].size())
    {
        // cout << "parity mismatch" << endl;
        cout << "-1" << "\n";
        return 0;
    }
 
    mcf_graph<int, long long> g(node_cnt + 2);
    int s = 0, t = node_cnt + 1;
 
    for(auto u:parity_nodes[0])
        g.add_edge(s,u,1,0);
    for(auto v:parity_nodes[1])
        g.add_edge(v,t,1,0);
 
    for(auto u:parity_nodes[0])
    {
        for(auto v:parity_nodes[1])
        {
            if(x_rep[u] == x_rep[v])
            {
                if(abs(y_rep[u]-y_rep[v]) == 1)
                    g.add_edge(u,v,1,cost_left+cost_right);
            }
            if(y_rep[u] == y_rep[v])
            {
                if(abs(x_rep[u]-x_rep[v]) == 1)
                    g.add_edge(u,v,1,cost_up+cost_down);
            }
        }
    }
 
    auto result = g.flow(s, t);
    // cout << "flow done : " << result.first << " " << result.second << "\n";
    if(result.first != (int)parity_nodes[0].size())
        cout << "-1" << "\n";
    else
    {
        auto edges = g.edges();
        for (auto e : edges) {
            if (e.from == s || e.to == t || e.flow == 0) continue;
 
            auto u = e.from;
            auto v = e.to;
            if(x_rep[u] == x_rep[v])
            {
                if(y_rep[u] == y_rep[v]+1)
                    swap(u,v);
                direction[x_rep[u]][y_rep[u]] = 'R';
                direction[x_rep[v]][y_rep[v]] = 'L';
            }
            else
            {
                if(x_rep[u] == x_rep[v]+1)
                    swap(u,v);
                direction[x_rep[u]][y_rep[u]] = 'D';
                direction[x_rep[v]][y_rep[v]] = 'U';
            }
        }
 
        cout << result.second+cost << endl;
        for(int i = 1; i <= n; ++i)
        {
            for(int j = 1; j <= n; ++j)
                cout << direction[i][j];
            cout << "\n";
        }
 
    }
    return 0;
}
Tester's Solution
#include <bits/stdc++.h>
using namespace std;
template<typename T = int> vector<T> create(size_t n){ return vector<T>(n); }
template<typename T, typename... Args> auto create(size_t n, Args... args){ return vector<decltype(create<T>(args...))>(n, create<T>(args...)); }
long long readInt(long long l,long long r,char endd){
    long long x=0;
    int cnt=0;
    int fi=-1;
    bool is_neg=false;
    while(true){
        char g=getchar();
        if(g=='-'){
            assert(fi==-1);
            is_neg=true;
            continue;
        }
        if('0'<=g && g<='9'){
            x*=10;
            x+=g-'0';
            if(cnt==0){
                fi=g-'0';
            }
            cnt++;
            assert(fi!=0 || cnt==1);
            assert(fi!=0 || is_neg==false);
 
            assert(!(cnt>19 || ( cnt==19 && fi>1) ));
        } else if(g==endd){
            if(is_neg){
                x= -x;
            }
            assert(l<=x && x<=r);
            return x;
        } else {
            assert(false);
        }
    }
}
string readString(int l,int r,char endd){
    string ret="";
    int cnt=0;
    while(true){
        char g=getchar();
        assert(g!=-1);
        if(g==endd){
            break;
        }
        cnt++;
        ret+=g;
    }
    assert(l<=cnt && cnt<=r);
    return ret;
}
long long readIntSp(long long l,long long r){
    return readInt(l,r,' ');
}
long long readIntLn(long long l,long long r){
    return readInt(l,r,'\n');
}
string readStringLn(int l,int r){
    return readString(l,r,'\n');
}
string readStringSp(int l,int r){
    return readString(l,r,' ');
}
// source: https://github.com/kth-competitive-programming/kactl/blob/master/content/graph/DirectedMST.h
#define rep(i, from, to) for (int i = from; i < (to); ++i)
#define all(x) x.begin(), x.end()
#define sz(x) (int)(x).size()
typedef long long ll;
typedef pair<int, int> pii;
typedef vector<int> vi;
struct RollbackUF {
    vi e; vector<pii> st;
    RollbackUF(int n) : e(n, -1) {}
    int size(int x) { return -e[find(x)]; }
    int find(int x) { return e[x] < 0 ? x : find(e[x]); }
    int time() { return sz(st); }
    void rollback(int t) {
        for (int i = time(); i --> t;)
            e[st[i].first] = st[i].second;
        st.resize(t);
    }
    bool join(int a, int b) {
        a = find(a), b = find(b);
        if (a == b) return false;
        if (e[a] > e[b]) swap(a, b);
        st.push_back({a, e[a]});
        st.push_back({b, e[b]});
        e[a] += e[b]; e[b] = a;
        return true;
    }
};
struct Edge { int a, b; ll w; };
struct Node { /// lazy skew heap node
    Edge key;
    Node *l, *r;
    ll delta;
    void prop() {
        key.w += delta;
        if (l) l->delta += delta;
        if (r) r->delta += delta;
        delta = 0;
    }
    Edge top() { prop(); return key; }
};
Node *merge(Node *a, Node *b) {
    if (!a || !b) return a ?: b;
    a->prop(), b->prop();
    if (a->key.w > b->key.w) swap(a, b);
    swap(a->l, (a->r = merge(b, a->r)));
    return a;
}
void pop(Node*& a) { a->prop(); a = merge(a->l, a->r); }
 
pair<ll, vi> dmst(int n, int r, vector<Edge>& g) {
    RollbackUF uf(n);
    vector<Node*> heap(n);
    for (Edge e : g) heap[e.b] = merge(heap[e.b], new Node{e});
    ll res = 0;
    vi seen(n, -1), path(n), par(n);
    seen[r] = r;
    vector<Edge> Q(n), in(n, {-1,-1}), comp;
    deque<tuple<int, int, vector<Edge>>> cycs;
    rep(s,0,n) {
        int u = s, qi = 0, w;
        while (seen[u] < 0) {
            if (!heap[u]) return {-1,{}};
            Edge e = heap[u]->top();
            heap[u]->delta -= e.w, pop(heap[u]);
            Q[qi] = e, path[qi++] = u, seen[u] = s;
            res += e.w, u = uf.find(e.a);
            if (seen[u] == s) { /// found cycle, contract
                Node* cyc = 0;
                int end = qi, time = uf.time();
                do cyc = merge(cyc, heap[w = path[--qi]]);
                while (uf.join(u, w));
                u = uf.find(u), heap[u] = cyc, seen[u] = -1;
                cycs.push_front({u, time, {&Q[qi], &Q[end]}});
            }
        }
        rep(i,0,qi) in[uf.find(Q[i].b)] = Q[i];
    }
 
    for (auto& [u,t,comp] : cycs) { // restore sol (optional)
        uf.rollback(t);
        Edge inEdge = in[u];
        for (auto& e : comp) in[uf.find(e.b)] = e;
        in[uf.find(inEdge.b)] = inEdge;
    }
    rep(i,0,n) par[i] = in[i].a;
    return {res, par};
}
struct SSP{
    int n, source, sink;
    vector<vector<int> > adj;
    struct edge {
        int u, flow, cap, cost, origcost;
    };
    vector<edge> edges;
    SSP(int n, int source, int sink): n(n),source(source),sink(sink),adj(n) {}
    int add_edge(int a, int b, int cap, int cost){
        int result = edges.size();
        adj[a].push_back(edges.size());
        edges.push_back({b, 0, cap, cost, cost});
        adj[b].push_back(edges.size());
        edges.push_back({a, 0, 0, -cost, -cost});
        return result;
    }
    pair<int, int> djk(){
        struct state_t{
            int dist, prev, maxf;
        };
        priority_queue<pair<int, int>,vector<pair<int, int> >, greater<pair<int, int> > > pq;
        vector<state_t> states(n, {INT_MAX/2, -1, 0});
        states[source] = {0, -2, INT_MAX};
        pq.push(make_pair(states[source].dist, source));
        while(!pq.empty()){
            int d, v; tie(d, v) = pq.top(); pq.pop();
            if(d != states[v].dist) continue;
            for(int edid : adj[v]){
                edge & ed = edges[edid];
                if(ed.flow < ed.cap && states[ed.u].dist > states[v].dist + ed.cost){
                    states[ed.u].dist = states[v].dist + ed.cost;
                    states[ed.u].prev = edid;
                    states[ed.u].maxf = min(states[v].maxf, ed.cap - ed.flow);
                    pq.push(make_pair(states[ed.u].dist, ed.u));
                }
            }
        }
        int lastid = states[sink].prev, extraflow = states[sink].maxf, atcost = 0;
        if(lastid == -1) return {0, 0};
        for(int id = lastid; id != -2; id = states[edges[id^1].u].prev){
            edges[id].flow += extraflow;
            edges[id^1].flow -= extraflow;
            atcost += edges[id].origcost;
        }
        for(int id = 0; id < (int) edges.size(); id++){
            int v = edges[id^1].u, u = edges[id].u;
            if(states[v].prev != -1 && states[u].prev != -1)
                edges[id].cost += states[v].dist - states[u].dist;
        }
        return {extraflow, atcost};
    }
    pair<int, int> mfmc(){
        int flow = 0, cost = 0;
        if(any_of(edges.begin(), edges.end(), [](const edge & ed){
            return ed.flow > 0 && ed.cost < 0;
        })) bellamanford();
        int extraflow, atcost;
        while(tie(extraflow, atcost) = djk(), extraflow > 0)
            flow += extraflow, cost += extraflow * atcost;
        return {flow, cost};
    }
    void bellamanford(){
        vector<int> dist(n, INT_MAX/2);
        dist[source] = 0;
        for(int rep = 0; rep < n; rep++){
            bool changed = false;
            for(int id = 0; id < (int)edges.size(); id++){
                int v = edges[id^1].u, u = edges[id].u;
                if(edges[id].flow < edges[id].cap && dist[v] != INT_MAX/2 &&
                    dist[u] > dist[v] + edges[id].cost){
                    dist[u] = dist[v] + edges[id].cost;
                    changed = true;
                    assert(rep < n - 1);
                }
            }
            if(!changed) break;
        }
        for(int id = 0; id < (int) edges.size(); id++){
            int v = edges[id^1].u, u = edges[id].u;
            if(dist[v] < INT_MAX/2 && dist[u] < INT_MAX/2)
                edges[id].cost += dist[v] - dist[u];
        }
    }
};
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);
 
    int n = readIntLn(1, 50);
    vector<int> p(4);
    for(int i = 0; i < 4; i++){
        if(i == 3) p[i] = readIntLn(1, 100000);
        else p[i] = readIntSp(1, 100000);
    }
 
    auto grid = create<pair<int,int>>(n, n);
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++){
            grid[i][j].first = readIntSp(1, n);
            if(j == n - 1) grid[i][j].second = readIntLn(1, n);
            else grid[i][j].second = readIntSp(1, n);
            grid[i][j].first -= 1;
            grid[i][j].second -= 1;
        }
 
    auto res = create<char>(n, n);
    auto dist = create<int>(n, n);
    auto seen = create<int>(n, n);
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++){
            dist[i][j] = (1<<30);
            res[i][j] = '*';
        }
 
    vector<int> cost, chr;
    vector<pair<int,int>> dir;
    dir.push_back({0, 1}); chr.push_back('R'); cost.push_back(p[3]);
    dir.push_back({1, 0}); chr.push_back('D'); cost.push_back(p[2]);
    dir.push_back({0, -1}); chr.push_back('L'); cost.push_back(p[1]);
    dir.push_back({-1, 0}); chr.push_back('U'); cost.push_back(p[0]);
 
    vector<pair<int,int>> roots;
    auto id = create<int>(n, n);
    for(int x = 0; x < n; x++) for(int y = 0; y < n; y++) {
        if(grid[x][y] != make_pair(x, y))
            continue;
        roots.push_back({x, y});
        vector<pair<int,int>> nodes;
        vector<Edge> edges;
        vector<char> fchars;
        for(int i = 0; i < n; i++) for(int j = 0; j < n; j++) {
            if(grid[i][j] == make_pair(x, y)){
                id[i][j] = nodes.size();
                nodes.push_back({i, j});
            }
        }
        map<pair<int,int>, int> rid;
        for(int i = 0; i < n; i++) for(int j = 0; j < n; j++) {
            if(grid[i][j] != make_pair(x, y)) continue;
            for(int r = 0; r < 4; r++){
                int ni = i + dir[r].first;
                int nj = j + dir[r].second;
                if(ni >= 0 && ni < n && nj >= 0 && nj < n){
                    if(grid[ni][nj] == make_pair(x, y)){
                        rid[{id[i][j], id[ni][nj]}] = edges.size();
                        edges.push_back((Edge){id[i][j], id[ni][nj], cost[r ^ 2]});
                        fchars.push_back(chr[r ^ 2]);
                    }
                }
            }
        }
        int nds = nodes.size();
        auto ans = dmst(nds, id[x][y], edges);
        if(ans.second.size() != nds){
            cout << -1 << '\n';
            return 0;
        }
        for(int i = 0; i < nds; i++){
            if(i == id[x][y]) continue;
            int u = ans.second[i], v = i;
            int eid = rid[{u, v}];
            res[nodes[v].first][nodes[v].second] = fchars[eid];
        }
    }
 
    auto color = [&](int i){
        return (roots[i].first + roots[i].second)%2;
    };
 
    auto mandist = [&](int i, int j){
        return abs(roots[i].first - roots[j].first) + abs(roots[i].second - roots[j].second);
    };
 
    int len = roots.size();
    SSP ssp(len + 2, len, len + 1);
    vector<int> to_check;
    for(int i = 0; i < len; i++){
        if(!color(i)) ssp.add_edge(len, i, 1, 0);
        else ssp.add_edge(i, len + 1, 1, 0);
        for(int j = i + 1; j < len; j++){
            if(mandist(i, j) == 1){
                if(roots[i].first == roots[j].first){
                    if(!color(i)) ssp.add_edge(i, j, 1, p[1] + p[3]);
                    else ssp.add_edge(j, i, 1, p[1] + p[3]);
                } else {
                    if(!color(i)) ssp.add_edge(i, j, 1, p[0] + p[2]);
                    else ssp.add_edge(j, i, 1, p[0] + p[2]);
                }
                to_check.push_back(ssp.edges.size() - 2);
            }
        }
    }
 
    ssp.mfmc();
 
    for(int edid : to_check){
        int i = ssp.edges[edid ^ 1].u, j = ssp.edges[edid].u;
        if(ssp.edges[edid].flow == 1){
            if(roots[i].first == roots[j].first){
                if(roots[i].second > roots[j].second) swap(i, j);
                res[roots[i].first][roots[i].second] = 'R';
                res[roots[j].first][roots[j].second] = 'L';
            } else {
                if(roots[i].first > roots[j].first) swap(i, j);
                res[roots[i].first][roots[i].second] = 'D';
                res[roots[j].first][roots[j].second] = 'U';
            }
        }
    }
 
    int ans = 0;
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++){
            if(res[i][j] == '*'){
                cout << -1 << '\n';
                return 0;
            }
            int at = find(chr.begin(), chr.end(), res[i][j]) - chr.begin();
            ans += cost[at];
        }
    
    cout << ans << '\n';
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            cout << res[i][j];
        }
        cout << '\n';
    }
 
    auto stopped = create<pair<int,int>>(n, n);
    for(int x = 0; x < n; x++){
        for(int y = 0; y < n; y++){
            auto visited = create<int>(n, n);
            visited[x][y] = 1;
            queue<pair<int,int>> q;
            q.push({x, y});
            while(!q.empty()){
                auto [i, j] = q.front(); q.pop();
                int at = find(chr.begin(), chr.end(), res[i][j]) - chr.begin();
                int ni = i + dir[at].first, nj = j + dir[at].second;
                if(ni >= 0 && ni < n && nj >= 0 && nj < n){
                    if(visited[ni][nj]){
                        assert(grid[x][y] == make_pair(ni, nj));
                        break;
                    } else {
                        visited[ni][nj] = 1;
                        q.push({ni, nj});
                    }
                } else {
                    assert(false);
                }
            }
        }
    }
    return 0;
} 

Behind the scenes :
When I had created the problem, I didn’t know that Prim’s won’t provide us with a minimum cost solution in directed graphs. While testing, Felipe found that there exists a better solution in one of the testcases with his solution, it wasn’t still the best though, then he introuduced to me the directed MST problem and Edmond’s algorithm to solve it. I used cescmentation_folch’s submission on Road Repairs problem as a template for my solution for Bash Matrix. Do give a try to the Road Repairs problem, if you are new to the concept of directed MST as well.

VIDEO EDITORIAL:

4 Likes

I really enjoyed the problem! The directed MST twist was a huge bonus, so congratulations to both setter and tester!

6 Likes

Hi,

I am relatively new to competitive programming and I program only in Python for now. May I ask where I can find an efficient implementation of max flow min cost in Python? I tried to look online, but I couldn’t find any except inside the package Nertworkx, and it seems complicated to extract the code from there and paste it inside my code.

After a few searches, C++ seems to have some of these implementations. Is it usually the case that it is much easier to find self-contained implementations of more advanced algorithms in C++ than in Python? I am assuming that the answer is yes since C++ is much more popular (and fast) for competitive programming, but I didn’t expect it to be that hard to do it in Python.

Very nice problem!

1 Like

Here is one submission I found by zkou on the problem E - MinCostFlow.
I sorted the submissions on the problem by execution times and it was one of the fastest solutions that didn’t use atcoder library. Hopefully it helps!

MCMF python template
import heapq

class mcf_graph:   
    def __init__(self, n):
        self.n = n
        self.pos = []
        self.g = [[] for _ in range(n)]
 
    def add_edge(self, from_, to, cap, cost):
        # assert 0 <= from_ < self.n
        # assert 0 <= to < self.n
        m = len(self.pos)
        self.pos.append((from_, len(self.g[from_])))
        self.g[from_].append(self.__class__._edge(to, len(self.g[to]), cap, cost))
        self.g[to].append(self.__class__._edge(from_, len(self.g[from_]) - 1, 0, -cost))
        return m 
 
    class edge:
        def __init__(self, from_, to, cap, flow, cost):
            self.from_ = from_
            self.to = to
            self.cap = cap
            self.flow = flow
            self.cost = cost
 
    def get_edge(self, i):
        _e = self.g[self.pos[i][0]][self.pos[i][1]]
        _re = self.g[_e.to][_e.rev]
        return self.__class__.edge(self.pos[i][0], _e.to, _e.cap + _re.cap, _re.cap, _e.cost)
 
    def edges(self):
        ret = []
        for i in range(len(self.pos)):
            _e = self.g[self.pos[i][0]][self.pos[i][1]]
            _re = self.g[_e.to][_e.rev]
            ret.append(self.__class__.edge(self.pos[i][0], _e.to, _e.cap + _re.cap, _re.cap, _e.cost))
        return ret
 
    def _dual_ref(self, s, t):
        self.dist = [float('inf')] * self.n
        self.pv = [-1] * self.n
        self.pe = [-1] * self.n
        self.vis = [False] * self.n
        que = [(0, s)]
        self.dist[s] = 0
        while que:
            _, v = heapq.heappop(que)
            if self.vis[v]:
                continue
            self.vis[v] = True
            if v == t:
                break
            for i in range(len(self.g[v])):
                e = self.g[v][i]
                if self.vis[e.to] or e.cap == 0:
                    continue
                cost = e.cost - self.dual[e.to] + self.dual[v]
                if self.dist[e.to] > self.dist[v] + cost:
                    self.dist[e.to] = self.dist[v] + cost
                    self.pv[e.to] = v
                    self.pe[e.to] = i
                    heapq.heappush(que, (self.dist[e.to], e.to))
        if not self.vis[t]:
            return False
        for v in range(self.n):
            if not self.vis[v]:
                continue
            self.dual[v] -= self.dist[t] - self.dist[v]
        
        return True
 
    def slope(self, s, t, flow_limit=float('inf')):
        self.dual = [0] * self.n
        self.dist = [float('inf')] * self.n
        self.pv = [-1] * self.n
        self.pe = [-1] * self.n
        self.vis = [False] * self.n
        flow = 0
        cost = 0
        prev_cost = -1
        result = [(flow, cost)]
        while flow < flow_limit:
            if not self._dual_ref(s, t):
                break
            c = flow_limit - flow
            v = t
            while v != s:
                c = min(c, self.g[self.pv[v]][self.pe[v]].cap)
                v = self.pv[v]
            v = t
            while v != s:
                e = self.g[self.pv[v]][self.pe[v]]
                e.cap -= c
                self.g[v][e.rev].cap += c
                v = self.pv[v]
            d = -self.dual[s]
            flow += c
            cost += c * d
            if prev_cost == d:
                result.pop()
            result.append((flow, cost))
            prev_cost = cost
        return result
 
    def flow(self, s, t, flow_limit=float('inf')):
        return self.slope(s, t, flow_limit)[-1]
    
    class _edge:
        def __init__(self, to, rev, cap, cost):
            self.to = to
            self.rev = rev
            self.cap = cap
            self.cost = cost
3 Likes

cool problem. Directed MST was a new concept for me. Couldn’t figure it out during contest.

1 Like

Thanks for the problem.
I spent 2 days just thinking that my code was not working because of my Bipartite matching code, and Prim not being optimal did not even enter my dreams. Then I found something called Arborescence accidentally that changed my perspective.
What I like about this problem is the easy reduction to a very non-trivial algorithm that gives great learning.

1 Like

Hi @shaanknight,

Thank you very much; you just opened my eyes :slight_smile: ! I know little about competitive programming. I am a former expert in combinatorics, and I am preparing for coding interviews. I always struggled to find efficient algos and data structures in Python for competitive programming. I was a bit reluctant to try to read tough papers and implement them. First, it takes a lot of time, and second, it might not be fast enough, especially with Python that is slow. I looked for a while for max flow min cost algos in Python for this problem, and I couldn’t find any. I also looked for the code of NTT for the “XOR Sums” problem of the same contest, but I didn’t find anything in Python except the code in Sympy. I extracted the code from Sympy but it was way too slow for the “XOR Sums” problem. I will try to see if they have a (fast) Python implementation of NTT (problem F - Convolution on AtCoder). Similarly, I couldn’t find an implementation of lazy segment tree, and there seems to be some there in Python also. I never participated at AtCoder, and I am still a bit confused by that library (in particular, there are some C++ algos here also: GitHub - atcoder/ac-library: AtCoder Library), but your post changes my view on CP a lot. I began to think that to do really well, you have to drop Python in favor of C++. But with all these algos in Python, there is hope.

Thank you enormously for your post @shaanknight, it is extremely helpful to me!!

1 Like

@shaanknight All the cycles will be disjoint ,right?

Actually all the algo’s and data structures you mentioned are available with python on Pyrival. But I agree that some questions just can’t work python even with several optimisations.

1 Like

Each node is supposed to have outdegree 1, so it’s supposedly not possible to have 2 cycles intersecting at a node.

Hi @anon10218360,
Thank you a lot for mentioning this repo! I will bookmark both for sure. And thanks for your comment about Python also. I am also new to this, and this confirms that at a certain (advanced) level, you need C++ to do well (or at least something else than Python). In particular, the monthly Long Challenge on CodeChef is very interesting, but Python seems to be too slow for some problems. But I have seen the same issues with bunch of other platforms. Time to learn C++ I guess!
Thanks a lot for your help guys :slight_smile:

Thanks for the problem. Learnt a lot regarding directed minimum spanning tree.

1 Like