← Back to branch summary

~ubuntu-branches/ubuntu/saucy/gfan/saucy-proposed

« back to all changes in this revision

Viewing changes to lp_soplexcdd.cpp

  • Committer: Package Import Robot
  • Author(s): Cédric Boutillier
  • Date: 2013-07-09 10:44:01 UTC
  • mfrom: (2.1.2 experimental)
  • Revision ID: package-import@ubuntu.com-20130709104401-5q66ozz5j5af0dak
Tags: 0.5+dfsg-3
* Upload to unstable.
* modify remove_failing_tests_on_32bits.patch to replace command of
  0009RenderStairCase test with an empty one instead of deleting it.
* remove lintian override about spelling error

Show diffs side-by-side

added added

removed removed

Lines of Context:
lp_soplexcdd.cpp
 
1
#include "lp_soplexcdd.h"
 
2
 
 
3
#include "printer.h"
 
4
 
 
5
#include "spxdefines.h"
 
6
#include "spxsolver.h"
 
7
 
 
8
#include "timer.h"
 
9
#include "spxpricer.h"
 
10
//#include "spxdefaultpr.h"
 
11
#include "spxparmultpr.h"
 
12
#include "spxdevexpr.h"
 
13
#include "spxhybridpr.h"
 
14
#include "spxsteeppr.h"
 
15
#include "spxweightpr.h"
 
16
#include "spxratiotester.h"
 
17
#include "spxharrisrt.h"
 
18
#include "spxdefaultrt.h"
 
19
#include "spxfastrt.h"
 
20
#include "spxsimplifier.h"
 
21
//#include "spxaggregatesm.h"
 
22
//#include "spxredundantsm.h"
 
23
//#include "spxrem1sm.h"
 
24
//#include "spxgeneralsm.h"
 
25
#include "spxscaler.h"
 
26
#include "spxequilisc.h"
 
27
#include "spxsumst.h"
 
28
#include "spxweightst.h"
 
29
#include "spxvectorst.h"
 
30
#include "slufactor.h"
 
31
#include "soplex.h"
 
32
#include "continuedfractions.h"
 
33
#include "matrix.h"
 
34
 
 
35
#include "log.h"
 
36
using namespace soplex;
 
37
 
 
38
/** Here comes a simple derived class from #SoPlex, which uses #terminate() as
 
39
 *  callback method for outputting statistics.
 
40
 */
 
41
class MySoPlex : public SoPlex
 
42
{
 
43
private:
 
44
   SLUFactor m_slu;
 
45
 
 
46
public:
 
47
   /// default constructor
 
48
   MySoPlex(SPxSolver::Type p_type = SPxSolver::LEAVE, SPxSolver::Representation p_rep = SPxSolver::COLUMN)
 
49
      : SoPlex(p_type, p_rep)
 
50
   {
 
51
 
 
52
 
 
53
 
 
54
   bool                   print_solution = false;
 
55
   bool                   print_quality  = false;
 
56
   NameSet                rownames;
 
57
   NameSet                colnames;
 
58
   SPxStarter*            starter        = 0;
 
59
   SPxSolver::Type           type           = SPxSolver::LEAVE;
 
60
   SPxSolver::Representation representation = SPxSolver::COLUMN;
 
61
   int                    precision;
 
62
   Real                   delta          = DEFAULT_BND_VIOL;
 
63
   Real                   epsilon        = DEFAULT_EPS_ZERO;
 
64
   int                    verbose        = 0;
 
65
   SLUFactor::UpdateType  update         = SLUFactor::FOREST_TOMLIN;
 
66
   Real                   timelimit      = 1.0;-1.0;
 
67
   SPxPricer*             pricer         = 0;
 
68
   SPxRatioTester*        ratiotester    = 0;
 
69
   SPxScaler*             scaler         = 0;
 
70
   SPxSimplifier*         simplifier     = 0;
 
71
 
 
72
   precision = int(-log10(delta)) + 1;
 
73
 
 
74
   Param::setEpsilon(epsilon);
 
75
   Param::setVerbose(verbose);
 
76
 
 
77
 
 
78
 
 
79
   //options -p4 -t2 -g1 -s0 -c0
 
80
   setUtype(update);
 
81
   setTerminationTime(timelimit);
 
82
   setDelta(delta);
 
83
 
 
84
   assert(isConsistent());
 
85
 
 
86
   pricer = new SPxSteepPR;
 
87
   setPricer(pricer);
 
88
   assert(isConsistent());
 
89
 
 
90
   ratiotester = new SPxFastRT;
 
91
   setTester(ratiotester);
 
92
   assert(isConsistent());
 
93
 
 
94
   /*   scaler = new SPxEquili(representation == SoPlex::COLUMN, true);
 
95
   setScaler(scaler);
 
96
   assert(isConsistent());
 
97
   */
 
98
   setSimplifier(simplifier);
 
99
   assert(isConsistent());
 
100
 
 
101
   setStarter(starter);
 
102
   assert(isConsistent());
 
103
   }
 
104
 
 
105
   virtual bool terminate()
 
106
   {
 
107
      /*      if (iteration() % 100 == 0)
 
108
         std::cout << iteration() << ":\t" << value() << std::endl;
 
109
      */
 
110
     return SoPlex::terminate();
 
111
   }
 
112
 
 
113
   void setUtype(SLUFactor::UpdateType tp)
 
114
   {
 
115
      m_slu.setUtype(tp);
 
116
   }
 
117
 
 
118
  void build(const IntegerVectorList &g, IntegerVectorList::const_iterator i)
 
119
   {
 
120
     int width=g.size()-1;
 
121
     int height=i->v.size();
 
122
     
 
123
     LPColSet cols(width,width*height);
 
124
     DSVector c1(height);
 
125
     
 
126
     //Build Cols
 
127
     for(IntegerVectorList::const_iterator j=g.begin();j!=g.end();j++)
 
128
       {
 
129
         c1.clear();
 
130
         for(int t=0;t<height;t++)
 
131
           if(j->v[t]!=0)
 
132
             c1.add(t,(double)(j->v[t]));
 
133
         
 
134
         if(j!=i)
 
135
           {
 
136
             Real obj=0;
 
137
             Real upper=infinity;
 
138
             Real lower=0;
 
139
             
 
140
             cols.add(obj,lower,c1,upper);
 
141
           }
 
142
       }
 
143
     
 
144
     LPRowSet rows(height,width*height);
 
145
     DSVector r1(width);
 
146
     
 
147
     //Change Rows
 
148
     for(int t=0;t<height;t++)
 
149
       rows.add(i->v[t],r1,i->v[t]);
 
150
     
 
151
     addRows(rows);
 
152
     addCols(cols);
 
153
     
 
154
     changeSense(SPxLP::MINIMIZE);
 
155
     
 
156
     assert(isConsistent());
 
157
   }
 
158
  void buildType2(int n, const IntegerVectorList &inequalities, const IntegerVectorList &equations)
 
159
   {
 
160
     int width=n;
 
161
     int nInequalities=inequalities.size();
 
162
     int height=inequalities.size()+equations.size();
 
163
 
 
164
     IntegerMatrix m=rowsToIntegerMatrix(inequalities,n);
 
165
     m.append(rowsToIntegerMatrix(equations,n));
 
166
     
 
167
     LPColSet cols(width,width*height);
 
168
     DSVector c1(height);
 
169
     
 
170
     //Build Cols
 
171
     for(int i=0;i<width;i++)
 
172
       {
 
173
         c1.clear();
 
174
         for(int t=0;t<height;t++)
 
175
           if(m[t][i]!=0)
 
176
             c1.add(t,(double)(m[t][i]));
 
177
 
 
178
         Real obj=0;
 
179
         Real upper=infinity;
 
180
         //Real lower=0;
 
181
         Real lower=-infinity;//is this correct?
 
182
 
 
183
         if(i==0)
 
184
           {
 
185
             upper=1;
 
186
             lower=1;
 
187
           }
 
188
 
 
189
         cols.add(obj,lower,c1,upper);
 
190
       }
 
191
     
 
192
     LPRowSet rows(height,width*height);
 
193
     DSVector r1(width);
 
194
     
 
195
     //Change Rows
 
196
     for(int t=0;t<nInequalities;t++)
 
197
       rows.add(0,r1,infinity);
 
198
     for(int t=nInequalities;t<height;t++)
 
199
       rows.add(0,r1,0);
 
200
     
 
201
     addRows(rows);
 
202
     addCols(cols);
 
203
     
 
204
     changeSense(SPxLP::MINIMIZE);
 
205
     
 
206
     assert(isConsistent());
 
207
   }
 
208
};
 
209
 
 
210
static int toint(float r)
 
211
{
 
212
   return *((int*)&r);
 
213
}
 
214
 
 
215
static void printLP(SPxLP &w)
 
216
{
 
217
      std::cout << "LP has " 
 
218
             << w.nRows() 
 
219
             << "\trows and\n       "
 
220
             << w.nCols() 
 
221
             << "\tcolumns" 
 
222
             << std::endl;
 
223
      int nr=w.nRows();
 
224
      int nc=w.nCols();
 
225
 
 
226
      for(int i=0;i<nr;i++)
 
227
         {      
 
228
            for(int j=0;j<nc;j++)
 
229
               {
 
230
                  LPRow R;
 
231
                  w.getRow(i,R);
 
232
                  //  if(j<R.rowVector().size())
 
233
                     std::cout<<R.rowVector()[j]<<" ";
 
234
                  //  else
 
235
                     //                     std::cout<<(0.0)<<" ";
 
236
               }
 
237
            std::cout<<std::endl;
 
238
         }
 
239
      std::cout<<std::endl;
 
240
 
 
241
      for(int i=0;i<nr;i++)
 
242
         {      
 
243
            for(int j=0;j<nc;j++)
 
244
               {
 
245
                  LPCol C;
 
246
                  w.getCol(j,C);
 
247
                  //                  if(i<C.colVector().size())
 
248
                     std::cout<<C.colVector()[i]<<" ";
 
249
                     // else
 
250
                     // std::cout<<(0.0)<<" ";
 
251
               }
 
252
            std::cout<<std::endl;
 
253
         }
 
254
 
 
255
      std::cout<<"cols:"<<std::endl;
 
256
 
 
257
      for(int j=0;j<nc;j++)
 
258
         {
 
259
            LPCol C;
 
260
            w.getCol(j,C);
 
261
            std::cout<<C.lower()<<" "<<C.upper()<<" "<<C.obj()<<std::endl;
 
262
            std::cout<<toint(C.lower())<<" "<<toint(C.upper())<<" "<<toint(C.obj())<<std::endl;
 
263
         }
 
264
 
 
265
      std::cout<<"rows:"<<std::endl;
 
266
 
 
267
      for(int i=0;i<nr;i++)
 
268
         {
 
269
            LPRow R;
 
270
            w.getRow(i,R);
 
271
            std::cout<<toint(R.lhs())<<" "<<toint(R.rhs())<<" "<<R.type()<<std::endl;
 
272
         }
 
273
}
 
274
 
 
275
 
 
276
MySoPlex work(SPxSolver::LEAVE, SPxSolver::COLUMN);
 
277
 
 
278
 
 
279
static bool isFeasibleSolution(IntegerVector const &solution, int denominator, IntegerVectorList const &g, IntegerVectorList::const_iterator i)
 
280
{
 
281
  if(denominator<=0)return false;
 
282
  // To do: Truncate  
 
283
  IntegerVector sum=denominator*(*i);
 
284
 
 
285
  int t=0;
 
286
  for(IntegerVectorList::const_iterator j=g.begin();j!=g.end();j++)
 
287
    if(j!=i)
 
288
      {
 
289
        if(solution[t]<0)return false;
 
290
        sum-=solution[t]*(*j);
 
291
        t++;
 
292
      }
 
293
  return sum.isZero();
 
294
}
 
295
 
 
296
static bool isInfeasibilityCertificate(IntegerVector const &certificate, IntegerVectorList const &g, IntegerVectorList::const_iterator i)
 
297
{
 
298
  IntegerVector c=certificate;
 
299
  // To do: add truncation on c
 
300
  if(dotLong(c,*i)<=0)return false;
 
301
 
 
302
  for(IntegerVectorList::const_iterator j=g.begin();j!=g.end();j++)
 
303
    if(i!=j)
 
304
      if(dotLong(c,*j)>0)return false;
 
305
  return true;
 
306
}
 
307
/*           for(IntegerVectorList::const_iterator j=g.begin();j!=g.end();j++)
 
308
               {
 
309
                 /*              double prod=0;
 
310
                 for(int i=0;i<work.nRows();i++)
 
311
                   {prod+=(*j)[i]*certificate[i];
 
312
                     //          fprintf(stderr,"%f \n",prod);
 
313
                   }
 
314
                 int num,den;doubleToFraction(prod,num,den);
 
315
                 */
 
316
                 //              fprintf(stderr,":%f: %i/%i\n",prod,num,den);
 
317
/*               fprintf(stderr,":%i\n",dotLong(c,*j));
 
318
               }
 
319
        */       
 
320
 
 
321
 
 
322
static bool isFeasibleSolutionType2(IntegerVector const &s, IntegerVectorList const &inequalities, IntegerVectorList const &equations)
 
323
{
 
324
  // To do: do truncation
 
325
  if(s[0]<=0)return false;
 
326
  for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
 
327
    if(dotLong(s,*i)<0)return false;
 
328
  for(IntegerVectorList::const_iterator i=equations.begin();i!=equations.end();i++)
 
329
    if(dotLong(s,*i)!=0)return false;
 
330
  return true;
 
331
}
 
332
 
 
333
 
 
334
static bool isInfeasibilityCertificateType2(IntegerVector const &c, int n, IntegerVectorList const &inequalities, IntegerVectorList const &equations)
 
335
{
 
336
  // To do: truncation
 
337
  int nInequalities=inequalities.size();
 
338
  if(!c.subvector(0,nInequalities).isNonNegative())return false;
 
339
 
 
340
  IntegerVector sum(n);
 
341
  int j=0;
 
342
  for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++,j++)
 
343
    sum+=c[j]* *i;
 
344
  for(IntegerVectorList::const_iterator i=equations.begin();i!=equations.end();i++,j++)
 
345
    sum+=c[j]* *i;
 
346
 
 
347
  if(sum[0]>=0)return false;
 
348
 
 
349
  for(int i=1;i<n;i++)
 
350
    if(sum[i]<0)
 
351
      return false;
 
352
 
 
353
  return true;
 
354
}
 
355
 
 
356
 
 
357
bool LpSolverSoPlexCddGmp::isFacet(const IntegerVectorList &g, IntegerVectorList::const_iterator I)
 
358
{
 
359
   SPxSolver::Type           type           = SPxSolver::LEAVE;
 
360
   SPxSolver::Representation representation = SPxSolver::COLUMN;
 
361
 
 
362
   int lp_status=0;
 
363
 
 
364
   work.clear();
 
365
 
 
366
   work.build(g,I);
 
367
 
 
368
 retry:
 
369
 
 
370
   //   std::cerr<< work;
 
371
   work.solve();
 
372
 
 
373
   SPxSolver::Status stat = work.status();
 
374
 
 
375
   switch (stat)
 
376
     {
 
377
     case SPxSolver::OPTIMAL:
 
378
       {
 
379
         DVector objx(work.nCols());
 
380
         
 
381
         if( work.getPrimal(objx) != SPxSolver::ERROR )
 
382
           {
 
383
             vector<double> solution(work.nCols());
 
384
             for(int i=0;i<work.nCols();i++)
 
385
               solution[i]=objx[i];
 
386
             
 
387
             vector<int> solutionNum;
 
388
             int denominator;
 
389
             doubleVectorToFractions(solution,solutionNum,denominator);
 
390
             IntegerVector s(solution.size());
 
391
             for(int i=0;i<s.size();i++)s[i]=solutionNum[i];
 
392
             
 
393
             if(isFeasibleSolution(s,denominator,g,I))
 
394
               {
 
395
                 log3 fprintf(Stderr,"Solution OK.\n");
 
396
                 return false;
 
397
               }
 
398
             log3 fprintf(Stderr,"Solution failed .\n");
 
399
             goto fallBack;
 
400
           }
 
401
       }
 
402
       break;
 
403
     case SPxSolver::UNBOUNDED:
 
404
       std::cerr << "LP is unbounded";
 
405
       lp_status=1;
 
406
       break;
 
407
     case SPxSolver::INFEASIBLE:
 
408
       {
 
409
         DVector farkasx(work.nRows());
 
410
         
 
411
         if( work.getDualfarkas(farkasx) != SPxSolver::ERROR )
 
412
           {
 
413
             vector<double> certificate(work.nRows());
 
414
             for(int i=0;i<work.nRows();i++)
 
415
               certificate[i]=farkasx[i];
 
416
             
 
417
             vector<int> certificateNum;
 
418
             int denominator;
 
419
             doubleVectorToFractions(certificate,certificateNum,denominator);
 
420
             IntegerVector c(certificate.size());
 
421
             for(int i=0;i<c.size();i++)c[i]=certificateNum[i];
 
422
             
 
423
             if(isInfeasibilityCertificate(c,g,I))
 
424
               {
 
425
                 log3 fprintf(Stderr,"Certificate for infeasibility OK.\n");
 
426
                 return true;
 
427
               }
 
428
             log3 fprintf(Stderr,"Certificate failed.\n");
 
429
           }
 
430
         else
 
431
           {
 
432
             log3 fprintf(Stderr,"Error while producing certificate\n");
 
433
           }
 
434
         goto fallBack;
 
435
       }
 
436
       break;
 
437
     case SPxSolver::ABORT_TIME:
 
438
       std::cout << "aborted due to time limit";
 
439
       lp_status=1;
 
440
       break;
 
441
     case SPxSolver::ABORT_ITER:
 
442
       std::cout << "aborted due to iteration limit";
 
443
       lp_status=1;
 
444
       break;
 
445
     case SPxSolver::ABORT_VALUE:
 
446
       std::cout << "aborted due to objective value limit";
 
447
       lp_status=1;
 
448
       break;
 
449
     default:
 
450
       std::cout << "An error occurred during the solution process";
 
451
       lp_status=1;
 
452
       break;
 
453
     }
 
454
 
 
455
 fallBack:
 
456
   log0 fprintf(Stderr,"Falling back on CddLib\n");
 
457
   return LpSolverCddGmp::isFacet(g,I);
 
458
}
 
459
 
 
460
 
 
461
 
 
462
bool LpSolverSoPlexCddGmp::hasHomogeneousSolution(int n, const IntegerVectorList &inequalities, const IntegerVectorList &equations)
 
463
{
 
464
   SPxSolver::Type           type           = SPxSolver::LEAVE;
 
465
   SPxSolver::Representation representation = SPxSolver::COLUMN;
 
466
 
 
467
   int lp_status=0;
 
468
 
 
469
   work.clear();
 
470
 
 
471
   work.buildType2(n,inequalities,equations);
 
472
 
 
473
 retry:
 
474
 
 
475
   //   std::cerr<< work;
 
476
 
 
477
   //   assert(0);
 
478
   work.solve();
 
479
 
 
480
   SPxSolver::Status stat = work.status();
 
481
 
 
482
   switch (stat)
 
483
     {
 
484
     case SPxSolver::OPTIMAL:
 
485
       {
 
486
         DVector objx(work.nCols());
 
487
         
 
488
         if( work.getPrimal(objx) != SPxSolver::ERROR )
 
489
           {
 
490
             vector<double> solution(work.nCols());
 
491
             for(int i=0;i<work.nCols();i++)
 
492
               solution[i]=objx[i];
 
493
             
 
494
             vector<int> solutionNum;
 
495
             int denominator;
 
496
             doubleVectorToFractions(solution,solutionNum,denominator);
 
497
             IntegerVector s(solution.size());
 
498
             for(int i=0;i<s.size();i++)s[i]=solutionNum[i];
 
499
             
 
500
             //  AsciiPrinter(Stderr).printVector(s);
 
501
             if(isFeasibleSolutionType2(s,inequalities,equations))
 
502
               {
 
503
                 log3 fprintf(Stderr,"Solution OK.\n");
 
504
                 return true;
 
505
               }
 
506
             log2 fprintf(Stderr,"Solution failed (Type2).\n");
 
507
 
 
508
             /*      for(int i=0;i<work.nCols();i++)
 
509
               {
 
510
                 std::cerr<<solution[i]<<',';
 
511
               }
 
512
             std::cerr<< work;
 
513
             AsciiPrinter(Stderr).printVector(s);
 
514
             AsciiPrinter(Stderr).printVectorList(inequalities);
 
515
             AsciiPrinter(Stderr).printVectorList(equations);
 
516
             assert(0);
 
517
             */
 
518
             goto fallBack;
 
519
           }
 
520
       }
 
521
       break;
 
522
     case SPxSolver::UNBOUNDED:
 
523
       std::cerr << "LP is unbounded";
 
524
       lp_status=1;
 
525
       break;
 
526
     case SPxSolver::INFEASIBLE:
 
527
       {
 
528
         DVector farkasx(work.nRows());
 
529
         
 
530
         if( work.getDualfarkas(farkasx) != SPxSolver::ERROR )
 
531
           {
 
532
             vector<double> certificate(work.nRows());
 
533
             for(int i=0;i<work.nRows();i++)
 
534
               certificate[i]=farkasx[i];
 
535
             
 
536
             vector<int> certificateNum;
 
537
             int denominator;
 
538
             doubleVectorToFractions(certificate,certificateNum,denominator);
 
539
             IntegerVector c(certificate.size());
 
540
             for(int i=0;i<c.size();i++)c[i]=certificateNum[i];
 
541
             
 
542
             if(isInfeasibilityCertificateType2(c,n,inequalities,equations))
 
543
               {
 
544
                 log3 fprintf(Stderr,"Certificate for infeasibility OK.\n");
 
545
                 return false;
 
546
               }
 
547
             
 
548
               log2 fprintf(Stderr,"Certificate failed (Type2).\n");
 
549
               /* std::cerr<< work;
 
550
              std::cerr<< farkasx;
 
551
             AsciiPrinter(Stderr).printVector(c);
 
552
             AsciiPrinter(Stderr).printVectorList(inequalities);
 
553
             AsciiPrinter(Stderr).printVectorList(equations);
 
554
             assert(0);
 
555
             */
 
556
           }
 
557
         else
 
558
           {
 
559
             log3 fprintf(Stderr,"Error while producing certificate\n");
 
560
           }
 
561
         goto fallBack;
 
562
         }
 
563
       goto fallBack;
 
564
       break;
 
565
     case SPxSolver::ABORT_TIME:
 
566
       std::cout << "aborted due to time limit";
 
567
       lp_status=1;
 
568
       break;
 
569
     case SPxSolver::ABORT_ITER:
 
570
       std::cout << "aborted due to iteration limit";
 
571
       lp_status=1;
 
572
       break;
 
573
     case SPxSolver::ABORT_VALUE:
 
574
       std::cout << "aborted due to objective value limit";
 
575
       lp_status=1;
 
576
       break;
 
577
     default:
 
578
       std::cout << "An error occurred during the solution process";
 
579
       lp_status=1;
 
580
       break;
 
581
     }
 
582
 
 
583
 fallBack:
 
584
   log0 fprintf(Stderr,"Falling back on CddLib\n");
 
585
   return LpSolverCddGmp::hasHomogeneousSolution(n, inequalities,equations);
 
586
}
 
587
 
 
588
 
 
589
static LpSolverSoPlexCddGmp theLpSolverSoPlexCdd;