001package algs44; // section 6.5
002import stdlib.*;
003
004public class XAssignmentProblemDense {
005        private static final int UNMATCHED = -1;
006
007        private final int N;              // number of rows and columns
008        private final double[][] weight;  // the N-by-N cost matrix
009        private final double[] px;        // px[i] = dual variable for row i
010        private final double[] py;        // py[j] = dual variable for col j
011        private final int[] xy;           // xy[i] = j means i-j is a match
012        private final int[] yx;           // yx[j] = i means i-j is a match
013
014
015        public XAssignmentProblemDense(double[][] weight) {
016                this.weight = weight;
017                N = weight.length;
018
019                // dual variables
020                px = new double[N];
021                py = new double[N];
022
023                // initial matching is empty
024                xy = new int[N];
025                yx = new int[N];
026                for (int i = 0; i < N; i++) xy[i] = UNMATCHED;
027                for (int j = 0; j < N; j++) yx[j] = UNMATCHED;
028
029                // add N edges to matching
030                for (int k = 0; k < N; k++) {
031                        StdOut.println(k);
032                        assert isDualFeasible();
033                        assert isComplementarySlack();
034                        augment();
035                }
036                assert check();
037        }
038
039        // find shortest augmenting path and upate
040        private void augment() {
041
042                // build residual graph
043                EdgeWeightedDigraph G = new EdgeWeightedDigraph(2*N+2);
044                int s = 2*N, t = 2*N+1;
045                for (int i = 0; i < N; i++) {
046                        if (xy[i] == UNMATCHED) G.addEdge(new DirectedEdge(s, i, 0.0));
047                }
048                for (int j = 0; j < N; j++) {
049                        if (yx[j] == UNMATCHED) G.addEdge(new DirectedEdge(N+j, t, py[j]));
050                }
051                for (int i = 0; i < N; i++) {
052                        for (int j = 0; j < N; j++) {
053                                if (xy[i] == j) G.addEdge(new DirectedEdge(N+j, i, 0.0));
054                                else            G.addEdge(new DirectedEdge(i, N+j, reduced(i, j)));
055                        }
056                }
057
058                // compute shortest path from s to every other vertex
059                // DenseDijkstraSP spt = new DenseDijkstraSP(G, s);
060                DijkstraSP spt = new DijkstraSP(G, s);
061
062                // augment along alternating path
063                for (DirectedEdge e : spt.pathTo(t)) {
064                        int v = e.from(), w = e.to() - N;
065                        if (v < N) {
066                                xy[v] = w;
067                                yx[w] = v;
068                        }
069                }
070
071                // update dual variables
072                for (int i = 0; i < N; i++) px[i] += spt.distTo(i);
073                for (int j = 0; j < N; j++) py[j] += spt.distTo(N+j);
074        }
075
076        // reduced cost of i-j
077        private double reduced(int i, int j) {
078                return weight[i][j] + px[i] - py[j];
079        }
080
081        // total weight of min weight perfect matching
082        public double weight() {
083                double total = 0.0;
084                for (int i = 0; i < N; i++) {
085                        if (xy[i] != UNMATCHED)
086                                total += weight[i][xy[i]];
087                }
088                return total;
089        }
090
091        public int sol(int i) {
092                return xy[i];
093        }
094
095        // check that dual variables are feasible
096        private boolean isDualFeasible() {
097                // check that all edges have >= 0 reduced cost
098                for (int i = 0; i < N; i++) {
099                        for (int j = 0; j < N; j++) {
100                                if (reduced(i, j) < 0) {
101                                        StdOut.println("Dual variables are not feasible");
102                                        return false;
103                                }
104                        }
105                }
106                return true;
107        }
108
109        // check that primal and dual variables are complementary slack
110        private boolean isComplementarySlack() {
111
112                // check that all matched edges have 0-reduced cost
113                for (int i = 0; i < N; i++) {
114                        if ((xy[i] != UNMATCHED) && (reduced(i, xy[i]) != 0)) {
115                                StdOut.println("Primal and dual variables are not complementary slack");
116                                return false;
117                        }
118                }
119                return true;
120        }
121
122        // check that primal variables are a perfect matching
123        private boolean isPerfectMatching() {
124
125                // check that xy[] is a perfect matching
126                boolean[] perm = new boolean[N];
127                for (int i = 0; i < N; i++) {
128                        if (perm[xy[i]]) {
129                                StdOut.println("Not a perfect matching");
130                                return false;
131                        }
132                        perm[xy[i]] = true;
133                }
134
135                // check that xy[] and yx[] are inverses
136                for (int j = 0; j < N; j++) {
137                        if (xy[yx[j]] != j) {
138                                StdOut.println("xy[] and yx[] are not inverses");
139                                return false;
140                        }
141                }
142                for (int i = 0; i < N; i++) {
143                        if (yx[xy[i]] != i) {
144                                StdOut.println("xy[] and yx[] are not inverses");
145                                return false;
146                        }
147                }
148
149                return true;
150        }
151
152
153        // check optimality conditions
154        private boolean check() {
155                return isPerfectMatching() && isDualFeasible() && isComplementarySlack();
156        }
157
158        public static void main(String[] args) {
159
160                int N = Integer.parseInt(args[0]);
161                double[][] weight = new double[N][N];
162                for (int i = 0; i < N; i++)
163                        for (int j = 0; j < N; j++)
164                                weight[i][j] = StdRandom.random();
165
166                XAssignmentProblemDense assignment = new XAssignmentProblemDense(weight);
167                StdOut.println("weight = " + assignment.weight());
168                for (int i = 0; i < N; i++)
169                        StdOut.println(i + "-" + assignment.sol(i));
170        }
171
172}