test/lp_test.cc
author Balazs Dezso <deba@inf.elte.hu>
Tue, 02 Dec 2008 21:40:33 +0100
changeset 481 7afc121e0689
child 482 ed54c0d13df0
permissions -rw-r--r--
Port LP and MIP solvers from SVN -r3509 (#44)
deba@481
     1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
deba@481
     2
 *
deba@481
     3
 * This file is a part of LEMON, a generic C++ optimization library.
deba@481
     4
 *
deba@481
     5
 * Copyright (C) 2003-2008
deba@481
     6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
deba@481
     7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
deba@481
     8
 *
deba@481
     9
 * Permission to use, modify and distribute this software is granted
deba@481
    10
 * provided that this copyright notice appears in all copies. For
deba@481
    11
 * precise terms see the accompanying LICENSE file.
deba@481
    12
 *
deba@481
    13
 * This software is provided "AS IS" with no warranty of any kind,
deba@481
    14
 * express or implied, and with no claim as to its suitability for any
deba@481
    15
 * purpose.
deba@481
    16
 *
deba@481
    17
 */
deba@481
    18
deba@481
    19
#include <sstream>
deba@481
    20
#include <lemon/lp_skeleton.h>
deba@481
    21
#include "test_tools.h"
deba@481
    22
#include <lemon/tolerance.h>
deba@481
    23
deba@481
    24
#ifdef HAVE_CONFIG_H
deba@481
    25
#include <lemon/config.h>
deba@481
    26
#endif
deba@481
    27
deba@481
    28
#ifdef HAVE_GLPK
deba@481
    29
#include <lemon/lp_glpk.h>
deba@481
    30
#endif
deba@481
    31
deba@481
    32
#ifdef HAVE_CPLEX
deba@481
    33
#include <lemon/lp_cplex.h>
deba@481
    34
#endif
deba@481
    35
deba@481
    36
#ifdef HAVE_SOPLEX
deba@481
    37
#include <lemon/lp_soplex.h>
deba@481
    38
#endif
deba@481
    39
deba@481
    40
using namespace lemon;
deba@481
    41
deba@481
    42
void lpTest(LpSolverBase & lp)
deba@481
    43
{
deba@481
    44
deba@481
    45
deba@481
    46
deba@481
    47
  typedef LpSolverBase LP;
deba@481
    48
deba@481
    49
  std::vector<LP::Col> x(10);
deba@481
    50
  //  for(int i=0;i<10;i++) x.push_back(lp.addCol());
deba@481
    51
  lp.addColSet(x);
deba@481
    52
  lp.colLowerBound(x,1);
deba@481
    53
  lp.colUpperBound(x,1);
deba@481
    54
  lp.colBounds(x,1,2);
deba@481
    55
#ifndef GYORSITAS
deba@481
    56
deba@481
    57
  std::vector<LP::Col> y(10);
deba@481
    58
  lp.addColSet(y);
deba@481
    59
deba@481
    60
  lp.colLowerBound(y,1);
deba@481
    61
  lp.colUpperBound(y,1);
deba@481
    62
  lp.colBounds(y,1,2);
deba@481
    63
deba@481
    64
  std::map<int,LP::Col> z;
deba@481
    65
deba@481
    66
  z.insert(std::make_pair(12,INVALID));
deba@481
    67
  z.insert(std::make_pair(2,INVALID));
deba@481
    68
  z.insert(std::make_pair(7,INVALID));
deba@481
    69
  z.insert(std::make_pair(5,INVALID));
deba@481
    70
deba@481
    71
  lp.addColSet(z);
deba@481
    72
deba@481
    73
  lp.colLowerBound(z,1);
deba@481
    74
  lp.colUpperBound(z,1);
deba@481
    75
  lp.colBounds(z,1,2);
deba@481
    76
deba@481
    77
  {
deba@481
    78
    LP::Expr e,f,g;
deba@481
    79
    LP::Col p1,p2,p3,p4,p5;
deba@481
    80
    LP::Constr c;
deba@481
    81
deba@481
    82
    p1=lp.addCol();
deba@481
    83
    p2=lp.addCol();
deba@481
    84
    p3=lp.addCol();
deba@481
    85
    p4=lp.addCol();
deba@481
    86
    p5=lp.addCol();
deba@481
    87
deba@481
    88
    e[p1]=2;
deba@481
    89
    e.constComp()=12;
deba@481
    90
    e[p1]+=2;
deba@481
    91
    e.constComp()+=12;
deba@481
    92
    e[p1]-=2;
deba@481
    93
    e.constComp()-=12;
deba@481
    94
deba@481
    95
    e=2;
deba@481
    96
    e=2.2;
deba@481
    97
    e=p1;
deba@481
    98
    e=f;
deba@481
    99
deba@481
   100
    e+=2;
deba@481
   101
    e+=2.2;
deba@481
   102
    e+=p1;
deba@481
   103
    e+=f;
deba@481
   104
deba@481
   105
    e-=2;
deba@481
   106
    e-=2.2;
deba@481
   107
    e-=p1;
deba@481
   108
    e-=f;
deba@481
   109
deba@481
   110
    e*=2;
deba@481
   111
    e*=2.2;
deba@481
   112
    e/=2;
deba@481
   113
    e/=2.2;
deba@481
   114
deba@481
   115
    e=((p1+p2)+(p1-p2)+(p1+12)+(12+p1)+(p1-12)+(12-p1)+
deba@481
   116
       (f+12)+(12+f)+(p1+f)+(f+p1)+(f+g)+
deba@481
   117
       (f-12)+(12-f)+(p1-f)+(f-p1)+(f-g)+
deba@481
   118
       2.2*f+f*2.2+f/2.2+
deba@481
   119
       2*f+f*2+f/2+
deba@481
   120
       2.2*p1+p1*2.2+p1/2.2+
deba@481
   121
       2*p1+p1*2+p1/2
deba@481
   122
       );
deba@481
   123
deba@481
   124
deba@481
   125
    c = (e  <= f  );
deba@481
   126
    c = (e  <= 2.2);
deba@481
   127
    c = (e  <= 2  );
deba@481
   128
    c = (e  <= p1 );
deba@481
   129
    c = (2.2<= f  );
deba@481
   130
    c = (2  <= f  );
deba@481
   131
    c = (p1 <= f  );
deba@481
   132
    c = (p1 <= p2 );
deba@481
   133
    c = (p1 <= 2.2);
deba@481
   134
    c = (p1 <= 2  );
deba@481
   135
    c = (2.2<= p2 );
deba@481
   136
    c = (2  <= p2 );
deba@481
   137
deba@481
   138
    c = (e  >= f  );
deba@481
   139
    c = (e  >= 2.2);
deba@481
   140
    c = (e  >= 2  );
deba@481
   141
    c = (e  >= p1 );
deba@481
   142
    c = (2.2>= f  );
deba@481
   143
    c = (2  >= f  );
deba@481
   144
    c = (p1 >= f  );
deba@481
   145
    c = (p1 >= p2 );
deba@481
   146
    c = (p1 >= 2.2);
deba@481
   147
    c = (p1 >= 2  );
deba@481
   148
    c = (2.2>= p2 );
deba@481
   149
    c = (2  >= p2 );
deba@481
   150
deba@481
   151
    c = (e  == f  );
deba@481
   152
    c = (e  == 2.2);
deba@481
   153
    c = (e  == 2  );
deba@481
   154
    c = (e  == p1 );
deba@481
   155
    c = (2.2== f  );
deba@481
   156
    c = (2  == f  );
deba@481
   157
    c = (p1 == f  );
deba@481
   158
    //c = (p1 == p2 );
deba@481
   159
    c = (p1 == 2.2);
deba@481
   160
    c = (p1 == 2  );
deba@481
   161
    c = (2.2== p2 );
deba@481
   162
    c = (2  == p2 );
deba@481
   163
deba@481
   164
    c = (2 <= e <= 3);
deba@481
   165
    c = (2 <= p1<= 3);
deba@481
   166
deba@481
   167
    c = (2 >= e >= 3);
deba@481
   168
    c = (2 >= p1>= 3);
deba@481
   169
deba@481
   170
    e[x[3]]=2;
deba@481
   171
    e[x[3]]=4;
deba@481
   172
    e[x[3]]=1;
deba@481
   173
    e.constComp()=12;
deba@481
   174
deba@481
   175
    lp.addRow(LP::INF,e,23);
deba@481
   176
    lp.addRow(LP::INF,3.0*(x[1]+x[2]/2)-x[3],23);
deba@481
   177
    lp.addRow(LP::INF,3.0*(x[1]+x[2]*2-5*x[3]+12-x[4]/3)+2*x[4]-4,23);
deba@481
   178
deba@481
   179
    lp.addRow(x[1]+x[3]<=x[5]-3);
deba@481
   180
    lp.addRow(-7<=x[1]+x[3]-12<=3);
deba@481
   181
    lp.addRow(x[1]<=x[5]);
deba@481
   182
deba@481
   183
    std::ostringstream buf;
deba@481
   184
deba@481
   185
deba@481
   186
    //Checking the simplify function
deba@481
   187
deba@481
   188
//     //How to check the simplify function? A map gives no information
deba@481
   189
//     //on the question whether a given key is or is not stored in it, or
deba@481
   190
//     //it does?
deba@481
   191
//   Yes, it does, using the find() function.
deba@481
   192
    e=((p1+p2)+(p1-p2));
deba@481
   193
    e.simplify();
deba@481
   194
    buf << "Coeff. of p2 should be 0";
deba@481
   195
    //    std::cout<<e[p1]<<e[p2]<<e[p3]<<std::endl;
deba@481
   196
    check(e.find(p2)==e.end(), buf.str());
deba@481
   197
deba@481
   198
deba@481
   199
deba@481
   200
deba@481
   201
    e=((p1+p2)+(p1-0.99*p2));
deba@481
   202
    //e.prettyPrint(std::cout);
deba@481
   203
    //(e<=2).prettyPrint(std::cout);
deba@481
   204
    double tolerance=0.001;
deba@481
   205
    e.simplify(tolerance);
deba@481
   206
    buf << "Coeff. of p2 should be 0.01";
deba@481
   207
    check(e[p2]>0, buf.str());
deba@481
   208
deba@481
   209
    tolerance=0.02;
deba@481
   210
    e.simplify(tolerance);
deba@481
   211
    buf << "Coeff. of p2 should be 0";
deba@481
   212
    check(e.find(p2)==e.end(), buf.str());
deba@481
   213
deba@481
   214
deba@481
   215
  }
deba@481
   216
deba@481
   217
  {
deba@481
   218
    LP::DualExpr e,f,g;
deba@481
   219
    LP::Row p1 = INVALID, p2 = INVALID, p3 = INVALID,
deba@481
   220
      p4 = INVALID, p5 = INVALID;
deba@481
   221
deba@481
   222
    e[p1]=2;
deba@481
   223
    e[p1]+=2;
deba@481
   224
    e[p1]-=2;
deba@481
   225
deba@481
   226
    e=p1;
deba@481
   227
    e=f;
deba@481
   228
deba@481
   229
    e+=p1;
deba@481
   230
    e+=f;
deba@481
   231
deba@481
   232
    e-=p1;
deba@481
   233
    e-=f;
deba@481
   234
deba@481
   235
    e*=2;
deba@481
   236
    e*=2.2;
deba@481
   237
    e/=2;
deba@481
   238
    e/=2.2;
deba@481
   239
deba@481
   240
    e=((p1+p2)+(p1-p2)+
deba@481
   241
       (p1+f)+(f+p1)+(f+g)+
deba@481
   242
       (p1-f)+(f-p1)+(f-g)+
deba@481
   243
       2.2*f+f*2.2+f/2.2+
deba@481
   244
       2*f+f*2+f/2+
deba@481
   245
       2.2*p1+p1*2.2+p1/2.2+
deba@481
   246
       2*p1+p1*2+p1/2
deba@481
   247
       );
deba@481
   248
  }
deba@481
   249
deba@481
   250
#endif
deba@481
   251
}
deba@481
   252
deba@481
   253
void solveAndCheck(LpSolverBase& lp, LpSolverBase::SolutionStatus stat,
deba@481
   254
                   double exp_opt) {
deba@481
   255
  using std::string;
deba@481
   256
  lp.solve();
deba@481
   257
  //int decimal,sign;
deba@481
   258
  std::ostringstream buf;
deba@481
   259
  buf << "Primalstatus should be: " << int(stat);
deba@481
   260
deba@481
   261
  //  itoa(stat,buf1, 10);
deba@481
   262
  check(lp.primalStatus()==stat, buf.str());
deba@481
   263
deba@481
   264
  if (stat ==  LpSolverBase::OPTIMAL) {
deba@481
   265
    std::ostringstream sbuf;
deba@481
   266
    sbuf << "Wrong optimal value: the right optimum is " << exp_opt;
deba@481
   267
    check(std::abs(lp.primalValue()-exp_opt) < 1e-3, sbuf.str());
deba@481
   268
    //+ecvt(exp_opt,2)
deba@481
   269
  }
deba@481
   270
}
deba@481
   271
deba@481
   272
void aTest(LpSolverBase & lp)
deba@481
   273
{
deba@481
   274
  typedef LpSolverBase LP;
deba@481
   275
deba@481
   276
 //The following example is very simple
deba@481
   277
deba@481
   278
  typedef LpSolverBase::Row Row;
deba@481
   279
  typedef LpSolverBase::Col Col;
deba@481
   280
deba@481
   281
deba@481
   282
  Col x1 = lp.addCol();
deba@481
   283
  Col x2 = lp.addCol();
deba@481
   284
deba@481
   285
deba@481
   286
  //Constraints
deba@481
   287
  Row upright=lp.addRow(x1+x2 <=1);
deba@481
   288
  lp.addRow(x1+x2 >=-1);
deba@481
   289
  lp.addRow(x1-x2 <=1);
deba@481
   290
  lp.addRow(x1-x2 >=-1);
deba@481
   291
  //Nonnegativity of the variables
deba@481
   292
  lp.colLowerBound(x1, 0);
deba@481
   293
  lp.colLowerBound(x2, 0);
deba@481
   294
  //Objective function
deba@481
   295
  lp.obj(x1+x2);
deba@481
   296
deba@481
   297
  lp.max();
deba@481
   298
deba@481
   299
  //Testing the problem retrieving routines
deba@481
   300
  check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!");
deba@481
   301
  check(lp.isMax(),"This is a maximization!");
deba@481
   302
  check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!");
deba@481
   303
  //  std::cout<<lp.colLowerBound(x1)<<std::endl;
deba@481
   304
  check(  lp.colLowerBound(x1)==0,
deba@481
   305
          "The lower bound for variable x1 should be 0.");
deba@481
   306
  check(  lp.colUpperBound(x1)==LpSolverBase::INF,
deba@481
   307
          "The upper bound for variable x1 should be infty.");
deba@481
   308
  LpSolverBase::Value lb,ub;
deba@481
   309
  lp.getRowBounds(upright,lb,ub);
deba@481
   310
  check(  lb==-LpSolverBase::INF,
deba@481
   311
          "The lower bound for the first row should be -infty.");
deba@481
   312
  check(  ub==1,"The upper bound for the first row should be 1.");
deba@481
   313
  LpSolverBase::Expr e = lp.row(upright);
deba@481
   314
  check(  e.size() == 2, "The row retrieval gives back wrong expression.");
deba@481
   315
  check(  e[x1] == 1, "The first coefficient should 1.");
deba@481
   316
  check(  e[x2] == 1, "The second coefficient should 1.");
deba@481
   317
deba@481
   318
  LpSolverBase::DualExpr de = lp.col(x1);
deba@481
   319
  check(  de.size() == 4, "The col retrieval gives back wrong expression.");
deba@481
   320
  check(  de[upright] == 1, "The first coefficient should 1.");
deba@481
   321
deba@481
   322
  LpSolverBase* clp = lp.copyLp();
deba@481
   323
deba@481
   324
  //Testing the problem retrieving routines
deba@481
   325
  check(clp->objCoeff(x1)==1,"First term should be 1 in the obj function!");
deba@481
   326
  check(clp->isMax(),"This is a maximization!");
deba@481
   327
  check(clp->coeff(upright,x1)==1,"The coefficient in question is 1!");
deba@481
   328
  //  std::cout<<lp.colLowerBound(x1)<<std::endl;
deba@481
   329
  check(  clp->colLowerBound(x1)==0,
deba@481
   330
          "The lower bound for variable x1 should be 0.");
deba@481
   331
  check(  clp->colUpperBound(x1)==LpSolverBase::INF,
deba@481
   332
          "The upper bound for variable x1 should be infty.");
deba@481
   333
deba@481
   334
  clp->getRowBounds(upright,lb,ub);
deba@481
   335
  check(  lb==-LpSolverBase::INF,
deba@481
   336
          "The lower bound for the first row should be -infty.");
deba@481
   337
  check(  ub==1,"The upper bound for the first row should be 1.");
deba@481
   338
  e = clp->row(upright);
deba@481
   339
  check(  e.size() == 2, "The row retrieval gives back wrong expression.");
deba@481
   340
  check(  e[x1] == 1, "The first coefficient should 1.");
deba@481
   341
  check(  e[x2] == 1, "The second coefficient should 1.");
deba@481
   342
deba@481
   343
  de = clp->col(x1);
deba@481
   344
  check(  de.size() == 4, "The col retrieval gives back wrong expression.");
deba@481
   345
  check(  de[upright] == 1, "The first coefficient should 1.");
deba@481
   346
deba@481
   347
  delete clp;
deba@481
   348
deba@481
   349
  //Maximization of x1+x2
deba@481
   350
  //over the triangle with vertices (0,0) (0,1) (1,0)
deba@481
   351
  double expected_opt=1;
deba@481
   352
  solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
deba@481
   353
deba@481
   354
  //Minimization
deba@481
   355
  lp.min();
deba@481
   356
  expected_opt=0;
deba@481
   357
  solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
deba@481
   358
deba@481
   359
  //Vertex (-1,0) instead of (0,0)
deba@481
   360
  lp.colLowerBound(x1, -LpSolverBase::INF);
deba@481
   361
  expected_opt=-1;
deba@481
   362
  solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
deba@481
   363
deba@481
   364
  //Erase one constraint and return to maximization
deba@481
   365
  lp.eraseRow(upright);
deba@481
   366
  lp.max();
deba@481
   367
  expected_opt=LpSolverBase::INF;
deba@481
   368
  solveAndCheck(lp, LpSolverBase::INFINITE, expected_opt);
deba@481
   369
deba@481
   370
  //Infeasibilty
deba@481
   371
  lp.addRow(x1+x2 <=-2);
deba@481
   372
  solveAndCheck(lp, LpSolverBase::INFEASIBLE, expected_opt);
deba@481
   373
deba@481
   374
  //Change problem and forget to solve
deba@481
   375
  lp.min();
deba@481
   376
  check(lp.primalStatus()==LpSolverBase::UNDEFINED,
deba@481
   377
        "Primalstatus should be UNDEFINED");
deba@481
   378
deba@481
   379
deba@481
   380
//   lp.solve();
deba@481
   381
//   if (lp.primalStatus()==LpSolverBase::OPTIMAL){
deba@481
   382
//     std::cout<< "Z = "<<lp.primalValue()
deba@481
   383
//              << " (error = " << lp.primalValue()-expected_opt
deba@481
   384
//              << "); x1 = "<<lp.primal(x1)
deba@481
   385
//              << "; x2 = "<<lp.primal(x2)
deba@481
   386
//              <<std::endl;
deba@481
   387
deba@481
   388
//   }
deba@481
   389
//   else{
deba@481
   390
//     std::cout<<lp.primalStatus()<<std::endl;
deba@481
   391
//     std::cout<<"Optimal solution not found!"<<std::endl;
deba@481
   392
//   }
deba@481
   393
deba@481
   394
deba@481
   395
deba@481
   396
}
deba@481
   397
deba@481
   398
deba@481
   399
int main()
deba@481
   400
{
deba@481
   401
  LpSkeleton lp_skel;
deba@481
   402
  lpTest(lp_skel);
deba@481
   403
deba@481
   404
#ifdef HAVE_GLPK
deba@481
   405
  LpGlpk lp_glpk1,lp_glpk2;
deba@481
   406
  lpTest(lp_glpk1);
deba@481
   407
  aTest(lp_glpk2);
deba@481
   408
#endif
deba@481
   409
deba@481
   410
#ifdef HAVE_CPLEX
deba@481
   411
  LpCplex lp_cplex1,lp_cplex2;
deba@481
   412
  lpTest(lp_cplex1);
deba@481
   413
  aTest(lp_cplex2);
deba@481
   414
#endif
deba@481
   415
deba@481
   416
#ifdef HAVE_SOPLEX
deba@481
   417
  LpSoplex lp_soplex1,lp_soplex2;
deba@481
   418
  lpTest(lp_soplex1);
deba@481
   419
  aTest(lp_soplex2);
deba@481
   420
#endif
deba@481
   421
deba@481
   422
  return 0;
deba@481
   423
}