[1153] | 1 | // -*- c++ -*- |
---|
| 2 | #include <iostream> |
---|
| 3 | #include <fstream> |
---|
| 4 | |
---|
| 5 | #include <lemon/time_measure.h> |
---|
| 6 | #include <lp_solver_base.h> |
---|
| 7 | |
---|
| 8 | using std::cout; |
---|
| 9 | using std::endl; |
---|
| 10 | using namespace lemon; |
---|
| 11 | |
---|
| 12 | /* |
---|
| 13 | for n=3,4 , the program is very fast |
---|
| 14 | for n=5, with glpk, the run takes hours |
---|
| 15 | */ |
---|
| 16 | |
---|
| 17 | int main(int, char **) { |
---|
| 18 | const int n=4; |
---|
| 19 | const double row_sum=(1.0+n*n)*n/2; |
---|
| 20 | Timer ts; |
---|
| 21 | ts.reset(); |
---|
| 22 | typedef LPGLPK LPSolver; |
---|
| 23 | typedef LPSolver::Col Col; |
---|
| 24 | LPSolver lp; |
---|
| 25 | typedef std::map<std::pair<int, int>, Col> Coords; |
---|
| 26 | Coords x; |
---|
| 27 | // we create a new variable for each entry |
---|
| 28 | // of the magic square |
---|
| 29 | for (int i=1; i<=n; ++i) { |
---|
| 30 | for (int j=1; j<=n; ++j) { |
---|
| 31 | Col col=lp.addCol(); |
---|
| 32 | x[std::make_pair(i,j)]=col; |
---|
| 33 | lp.setColLowerBound(col, 1.0); |
---|
| 34 | lp.setColUpperBound(col, double(n*n)); |
---|
| 35 | } |
---|
| 36 | } |
---|
| 37 | LPSolver::Expression expr3, expr4; |
---|
| 38 | for (int i=1; i<=n; ++i) { |
---|
| 39 | LPSolver::Expression expr1, expr2; |
---|
| 40 | for (int j=1; j<=n; ++j) { |
---|
| 41 | expr1+=x[std::make_pair(i, j)]; |
---|
| 42 | expr2+=x[std::make_pair(j, i)]; |
---|
| 43 | } |
---|
| 44 | // sum of rows and columns |
---|
| 45 | lp.addRow(expr1==row_sum); |
---|
| 46 | lp.addRow(expr2==row_sum); |
---|
| 47 | expr3+=x[std::make_pair(i, i)]; |
---|
| 48 | expr4+=x[std::make_pair(i, (n+1)-i)]; |
---|
| 49 | } |
---|
| 50 | // sum of the diagonal entries |
---|
| 51 | lp.addRow(expr3==row_sum); |
---|
| 52 | lp.addRow(expr4==row_sum); |
---|
| 53 | lp.solveSimplex(); |
---|
| 54 | cout << "elapsed time: " << ts << endl; |
---|
| 55 | for (int i=1; i<=n; ++i) { |
---|
| 56 | for (int j=1; j<=n; ++j) { |
---|
| 57 | cout << "x("<<i<<","<<j<<")="<<lp.getPrimal(x[std::make_pair(i,j)]) |
---|
| 58 | << endl; |
---|
| 59 | } |
---|
| 60 | } |
---|
| 61 | // we make new binary variables for each pair of |
---|
| 62 | // entries of the square to achieve that |
---|
| 63 | // the values of different entries are different |
---|
| 64 | lp.setMIP(); |
---|
| 65 | for (Coords::const_iterator it=x.begin(); it!=x.end(); ++it) { |
---|
| 66 | Coords::const_iterator jt=it; ++jt; |
---|
| 67 | for(; jt!=x.end(); ++jt) { |
---|
| 68 | Col col1=(*it).second; |
---|
| 69 | Col col2=(*jt).second; |
---|
| 70 | Col col=lp.addCol(); |
---|
| 71 | lp.setColLowerBound(col, 0.0); |
---|
| 72 | lp.setColUpperBound(col, 1.0); |
---|
| 73 | lp.addRow(double(-n*n+1.0)<=1.0*col2-1.0*col1-double(n*n)*col<=-1.0); |
---|
| 74 | lp.setColInt(col); |
---|
| 75 | } |
---|
| 76 | } |
---|
| 77 | cout << "elapsed time: " << ts << endl; |
---|
| 78 | lp.solveSimplex(); |
---|
| 79 | // let's solve the integer problem |
---|
| 80 | lp.solveBandB(); |
---|
| 81 | cout << "elapsed time: " << ts << endl; |
---|
| 82 | for (int i=1; i<=n; ++i) { |
---|
| 83 | for (int j=1; j<=n; ++j) { |
---|
| 84 | cout << "x("<<i<<","<<j<<")="<<lp.getMIPPrimal(x[std::make_pair(i,j)]) |
---|
| 85 | << endl; |
---|
| 86 | } |
---|
| 87 | } |
---|
| 88 | } |
---|