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