test/lp_test.cc
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 17 Apr 2009 18:04:36 +0200
changeset 609 e6927fe719e6
parent 461 08d495d48089
child 537 0fec6a017ead
child 540 9db62975c32b
permissions -rw-r--r--
Support >= and <= constraints in NetworkSimplex (#219, #234)

By default the same inequality constraints are supported as by
Circulation (the GEQ form), but the LEQ form can also be selected
using the problemType() function.

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