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}