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 | } |
---|