1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/athos/lp/magic_square.cc Tue Mar 22 11:45:47 2005 +0000
1.3 @@ -0,0 +1,90 @@
1.4 +// -*- c++ -*-
1.5 +#include <iostream>
1.6 +#include <fstream>
1.7 +
1.8 +#include <lemon/time_measure.h>
1.9 +#include <lp_solver_base.h>
1.10 +
1.11 +using std::cout;
1.12 +using std::endl;
1.13 +using namespace lemon;
1.14 +
1.15 +/*
1.16 + On an 1537Mhz PC, the run times with
1.17 + glpk are the following.
1.18 + for n=3,4, some secondes
1.19 + for n=5, 25 hours
1.20 + */
1.21 +
1.22 +int main(int, char **) {
1.23 + const int n=4;
1.24 + const double row_sum=(1.0+n*n)*n/2;
1.25 + Timer ts;
1.26 + ts.reset();
1.27 + typedef LPGLPK LPSolver;
1.28 + typedef LPSolver::Col Col;
1.29 + LPSolver lp;
1.30 + typedef std::map<std::pair<int, int>, Col> Coords;
1.31 + Coords x;
1.32 + // we create a new variable for each entry
1.33 + // of the magic square
1.34 + for (int i=1; i<=n; ++i) {
1.35 + for (int j=1; j<=n; ++j) {
1.36 + Col col=lp.addCol();
1.37 + x[std::make_pair(i,j)]=col;
1.38 + lp.setColLowerBound(col, 1.0);
1.39 + lp.setColUpperBound(col, double(n*n));
1.40 + }
1.41 + }
1.42 + LPSolver::Expression expr3, expr4;
1.43 + for (int i=1; i<=n; ++i) {
1.44 + LPSolver::Expression expr1, expr2;
1.45 + for (int j=1; j<=n; ++j) {
1.46 + expr1+=x[std::make_pair(i, j)];
1.47 + expr2+=x[std::make_pair(j, i)];
1.48 + }
1.49 + // sum of rows and columns
1.50 + lp.addRow(expr1==row_sum);
1.51 + lp.addRow(expr2==row_sum);
1.52 + expr3+=x[std::make_pair(i, i)];
1.53 + expr4+=x[std::make_pair(i, (n+1)-i)];
1.54 + }
1.55 + // sum of the diagonal entries
1.56 + lp.addRow(expr3==row_sum);
1.57 + lp.addRow(expr4==row_sum);
1.58 + lp.solveSimplex();
1.59 + cout << "elapsed time: " << ts << endl;
1.60 + for (int i=1; i<=n; ++i) {
1.61 + for (int j=1; j<=n; ++j) {
1.62 + cout << "x("<<i<<","<<j<<")="<<lp.getPrimal(x[std::make_pair(i,j)])
1.63 + << endl;
1.64 + }
1.65 + }
1.66 + // we make new binary variables for each pair of
1.67 + // entries of the square to achieve that
1.68 + // the values of different entries are different
1.69 + lp.setMIP();
1.70 + for (Coords::const_iterator it=x.begin(); it!=x.end(); ++it) {
1.71 + Coords::const_iterator jt=it; ++jt;
1.72 + for(; jt!=x.end(); ++jt) {
1.73 + Col col1=(*it).second;
1.74 + Col col2=(*jt).second;
1.75 + Col col=lp.addCol();
1.76 + lp.setColLowerBound(col, 0.0);
1.77 + lp.setColUpperBound(col, 1.0);
1.78 + lp.addRow(double(-n*n+1.0)<=1.0*col2-1.0*col1-double(n*n)*col<=-1.0);
1.79 + lp.setColInt(col);
1.80 + }
1.81 + }
1.82 + cout << "elapsed time: " << ts << endl;
1.83 + lp.solveSimplex();
1.84 + // let's solve the integer problem
1.85 + lp.solveBandB();
1.86 + cout << "elapsed time: " << ts << endl;
1.87 + for (int i=1; i<=n; ++i) {
1.88 + for (int j=1; j<=n; ++j) {
1.89 + cout << "x("<<i<<","<<j<<")="<<lp.getMIPPrimal(x[std::make_pair(i,j)])
1.90 + << endl;
1.91 + }
1.92 + }
1.93 +}