1 // Copyright (C) 2005-2009, Pierre Bonami and others. All Rights Reserved.
2 // Author: Pierre Bonami
3 // Tepper School of Business
4 // Carnegie Mellon University, Pittsburgh, PA 15213
5 // Date: 21/07/05
6 //
7 // $Id$
8 //
9 // This code is licensed under the terms of the Eclipse Public License (EPL).
10 //---------------------------------------------------------------------------
11 #include "CglLandPSimplex.hpp"
12 #include "CoinTime.hpp"
13 #ifdef COIN_HAS_OSICLP
14 #include "OsiClpSolverInterface.hpp"
15 #endif
16 #include "CoinIndexedVector.hpp"
17 #include <cassert>
18 #include <iterator>
19
20 #include <list>
21 #include <algorithm>
22
23 #define REMOVE_LOG 0
24
25 #define RED_COST_CHECK 1e-6
26
27
28 //#define OUT_CGLP_PIVOTS
29 #ifdef OUT_CGLP_PIVOTS
30 #include "CglLandPOutput.hpp"
31 #endif
32 #ifdef DEBUG_LAP
33
34 /* The function is not used anywhere (LL)
35 static void MyAssertFunc(bool c, const std::string &s, const std::string& file, unsigned int line){
36 if (c != true){
37 fprintf(stderr, "Failed MyAssertion: %s in %s line %i.\n", s.c_str(), file.c_str(), line);
38 throw -1;
39 }
40 }
41 */
42
DblGtAssertFunc(const double & a,const std::string & a_s,const double & b,const std::string & b_s,const std::string & file,unsigned int line)43 static void DblGtAssertFunc(const double& a, const std::string &a_s, const double&b, const std::string& b_s,
44 const std::string& file, unsigned int line)
45 {
46 if (a<b)
47 {
48 fprintf(stderr, "Failed comparison: %s = %f < %s =%f in %s line %i.\n", a_s.c_str(),
49 a, b_s.c_str(), b, file.c_str(), line);
50 throw -1;
51 }
52 }
53
DblEqAssertFunc(const double & a,const std::string & a_s,const double & b,const std::string & b_s,const std::string & file,unsigned int line)54 static void DblEqAssertFunc(const double& a, const std::string &a_s, const double&b, const std::string& b_s,
55 const std::string& file, unsigned int line)
56 {
57 CoinRelFltEq eq(1e-7);
58 if (!eq(a,b))
59 {
60 fprintf(stderr, "Failed comparison: %s = %f != %s =%f in %s line %i.\n", a_s.c_str(),
61 a, b_s.c_str(), b, file.c_str(), line);
62 throw -1;
63 }
64 }
65
VecModEqAssertFunc(const CoinIndexedVector & a,const std::string a_s,const CoinIndexedVector & b,const std::string b_s,const std::string file,unsigned int line)66 static void VecModEqAssertFunc(const CoinIndexedVector& a, const std::string a_s,
67 const CoinIndexedVector& b, const std::string b_s,
68 const std::string file, unsigned int line)
69 {
70 CoinRelFltEq eq(1e-7);
71 assert(a.capacity()==b.capacity());
72 int n = a.capacity();
73 const double * a_v = a.denseVector();
74 const double * b_v = b.denseVector();
75 bool failed=false;
76 for (int i = 0 ; i < n ; i++)
77 {
78 double a = a_v[i] - floor(a_v[i]);
79 double b = b_v[i] - floor(b_v[i]);
80 if (!eq(a,b) && !eq(a_v[i],b_v[i]))
81 {
82 fprintf(stderr, "Failed comparison: %s[%i] = %.15f != %s[%i] =%.15f in %s line %i.\n", a_s.c_str(),i,
83 a_v[i], b_s.c_str(), i, b_v[i], file.c_str(), line);
84 failed = true;
85 }
86 }
87 if (failed) throw -1;
88 }
89
VecEqAssertFunc(const CoinIndexedVector & a,const std::string a_s,const CoinIndexedVector & b,const std::string b_s,const std::string file,unsigned int line)90 static void VecEqAssertFunc(const CoinIndexedVector& a, const std::string a_s,
91 const CoinIndexedVector& b, const std::string b_s,
92 const std::string file, unsigned int line)
93 {
94 CoinRelFltEq eq(1e-7);
95 assert(a.capacity()==b.capacity());
96 int n = a.capacity();
97 const double * a_v = a.denseVector();
98 const double * b_v = b.denseVector();
99 bool failed=false;
100 for (int i = 0 ; i < n ; i++)
101 {
102 if (!eq(a_v[i],b_v[i]))
103 {
104 fprintf(stderr, "Failed comparison: %s[%i] = %f != %s[%i] =%f in %s line %i.\n", a_s.c_str(),i,
105 a_v[i], b_s.c_str(), i, b_v[i], file.c_str(), line);
106 failed = true;
107 }
108 }
109 if (failed) throw -1;
110 }
111
112
113 #define MAKE_STRING(exp) std::string(#exp)
114 #define MyAssert(exp) MyAssertFunc(exp, MAKE_STRING(exp), __FILE__, __LINE__);
115 #define DblEqAssert(a,b) DblEqAssertFunc(a,MAKE_STRING(a),b,MAKE_STRING(b), __FILE__, __LINE__);
116 #define DblGtAssert(a,b) DblGtAssertFunc(a,MAKE_STRING(a),b,MAKE_STRING(b), __FILE__, __LINE__);
117
118 #define VecEqAssert(a,b) VecEqAssertFunc(a, MAKE_STRING(a), b, MAKE_STRING(b), __FILE__, __LINE__);
119 #define VecModEqAssert(a,b) VecModEqAssertFunc(a, MAKE_STRING(a), b, MAKE_STRING(b), __FILE__, __LINE__);
120
checkVecFunc(const CoinIndexedVector & v)121 void checkVecFunc(const CoinIndexedVector &v)
122 {
123 CoinIndexedVector v2 = v;
124 double * x = v2.denseVector();
125 const int * idx = v2.getIndices();
126 int n = v2.getNumElements();
127 for (int i = 0 ; i < n ; i++)
128 {
129 x[idx[i]] = 0.;
130 }
131 n = v2.capacity();
132 for (int i = 0 ; i < n ; i++)
133 {
134 DblEqAssert(x[i],0.);
135 }
136 }
137
138 #define checkVec(a) checkVecFunc(a)
139 #else
140 #define MyAssert(exp)
141 #define DblEqAssert(a,b)
142 #define DblGtAssert(a,b)
143 #define VecEqAssert(a,b)
144 #define VecModEqAssert(a,b)
145
146 #define checkVec(a)
147 #endif
148
149 #include <algorithm>
150 //#define TEST_M3
151 namespace LAP
152 {
153
printTableau(std::ostream & os)154 void CglLandPSimplex::printTableau(std::ostream & os)
155 {
156 int width = 9;
157 os<<"Tableau at current basis"<<std::endl;
158 os<<" ";
159 //Head with non basics indices
160 for (int i = 0 ; i < ncols_orig_ ; i++)
161 {
162 os.width(width);
163 os.setf(std::ios_base::right, std::ios_base::adjustfield);
164 std::cout<<nonBasics_[i]<<" ";
165 }
166
167 os.width(width);
168 os.setf(std::ios_base::right, std::ios_base::adjustfield);
169 std::cout<<'b';
170
171 os<<std::endl;
172
173 //print row by row
174 for (int i = 0 ; i < nrows_ ; i++)
175 {
176 //int ind = basics_[i];
177 row_i_.num = i;
178 pullTableauRow(row_i_);
179 row_i_.print(os, width, nonBasics_, ncols_orig_);
180
181 }
182
183 }
184
185 bool
checkBasis()186 CglLandPSimplex::checkBasis()
187 {
188 //Check that basics_ is correct
189 int * basic2 = new int [nrows_];
190 si_->getBasics(basic2);
191 for (int i = 0; i < nrows_ ; i++)
192 assert(basics_[i]==basic2[i]);
193 delete [] basic2;
194 return true;
195 }
196
197
CglLandPSimplex(const OsiSolverInterface & si,const CglLandP::CachedData & cached,const CglLandP::Parameters & params,Validator & validator)198 CglLandPSimplex::CglLandPSimplex(const OsiSolverInterface &si,
199 const CglLandP::CachedData &cached,
200 const CglLandP::Parameters ¶ms,
201 Validator& validator):
202 #ifdef COIN_HAS_OSICLP
203 clp_(NULL),
204 #endif
205 row_k_(this),
206 original_row_k_(this),
207 row_i_(this),
208 new_row_(this),
209 gammas_(false),
210 rowFlags_(NULL),
211 col_in_subspace(),
212 colCandidateToLeave_(NULL),
213 basics_(NULL), nonBasics_(NULL),
214 M1_(), M2_(), M3_(),
215 sigma_(0), basis_(NULL), colsolToCut_(NULL),
216 colsol_(NULL),
217 ncols_orig_(0),nrows_orig_(0),
218 inDegenerateSequence_(false),
219 chosenReducedCostVal_(1e100),
220 original_index_(),
221 si_(NULL),
222 validator_(validator),
223 numPivots_(0),
224 numSourceRowEntered_(0),
225 numIncreased_(0)
226 {
227 ncols_orig_ = si.getNumCols();
228 nrows_orig_ = si.getNumRows();
229 handler_ = new CoinMessageHandler();
230 handler_->setLogLevel(2);
231 messages_ = LandPMessages();
232
233 si_ = const_cast<OsiSolverInterface *>(&si);
234
235 #ifdef COIN_HAS_OSICLP
236 OsiClpSolverInterface * clpSi = dynamic_cast<OsiClpSolverInterface *>(si_);
237 if (clpSi)
238 {
239 clp_ = clpSi;
240 }
241 #endif
242
243 int rowsize = ncols_orig_ + nrows_orig_ + 1;
244 row_k_.reserve(rowsize);
245 #ifndef NDEBUG
246 new_row_.reserve(rowsize);
247 #endif
248
249 lo_bounds_.resize(ncols_orig_ + nrows_orig_);
250 up_bounds_.resize(ncols_orig_ + nrows_orig_);
251
252 CoinCopyN(si.getColLower(),ncols_orig_, &lo_bounds_[0]);
253 CoinCopyN(si.getColUpper(),ncols_orig_,&up_bounds_[0]);
254 const double * rowUpper = si.getRowUpper();
255 const double * rowLower = si.getRowLower();
256 double infty = si.getInfinity();
257 int i=ncols_orig_;
258 for (int iRow = 0; iRow < nrows_orig_ ; iRow++, i++)
259 {
260 if (rowUpper[iRow] < infty)
261 lo_bounds_[i]=0.;
262 else lo_bounds_[i]= - infty;
263 if (rowLower[iRow] <= - infty)
264 up_bounds_[i] = infty;
265 else if (rowUpper[iRow] < infty)
266 {
267 lo_bounds_[i] = rowLower[iRow] - rowUpper[iRow];
268 up_bounds_[i] = 0;
269 }
270 else
271 up_bounds_[i] = 0.;
272 }
273 cuts_.resize(ncols_orig_);
274 if (params.pivotLimit != 0)
275 {
276 own_ = true;
277 rWk1_.resize(nrows_orig_);
278 rWk2_.resize(nrows_orig_);
279 rWk3_.resize(nrows_orig_);
280 rWk4_.resize(nrows_orig_);
281 rIntWork_.resize(nrows_orig_);
282
283 row_i_.reserve(rowsize);
284 rowFlags_ = new bool[nrows_orig_];
285 col_in_subspace.resize(ncols_orig_ + nrows_orig_);
286 colCandidateToLeave_ = new bool[ncols_orig_];
287 basics_ = new int[nrows_orig_];
288 nonBasics_ = new int[ncols_orig_];
289
290 colsolToCut_ = new double[ncols_orig_ + nrows_orig_];
291 colsol_ = new double[ncols_orig_ + nrows_orig_];
292 original_index_.resize(ncols_orig_ + nrows_orig_);
293 CoinIotaN(&original_index_[0],ncols_orig_ + nrows_orig_, 0);
294
295
296
297 }
298 else
299 {
300 nrows_ = nrows_orig_;
301 ncols_ = ncols_orig_;
302 original_index_.resize(ncols_orig_ + nrows_orig_);
303 CoinIotaN(&original_index_[0],ncols_orig_ + nrows_orig_, 0);
304 own_ = false;
305 si_->enableSimplexInterface(0);
306 basis_ = new CoinWarmStartBasis(*cached.basis_);
307 }
308 cacheUpdate(cached,params.sepSpace != CglLandP::Full);
309 if (params.normalization)
310 {
311 computeWeights(params.lhs_norm, params.normalization, params.rhsWeightType);
312 }
313 else rhs_weight_ = 1;
314 }
315
~CglLandPSimplex()316 CglLandPSimplex::~CglLandPSimplex()
317 {
318 delete handler_;
319 handler_ = NULL;
320 delete basis_;
321 basis_ = NULL;
322 if (own_)
323 {
324 delete [] rowFlags_;
325 rowFlags_ = NULL;
326 delete [] colCandidateToLeave_;
327 colCandidateToLeave_ = NULL;
328 delete [] basics_;
329 basics_ = NULL;
330 delete [] nonBasics_;
331 nonBasics_ = NULL;
332 delete [] colsolToCut_;
333 colsolToCut_ = NULL;
334 delete [] colsol_;
335 colsol_ = NULL;
336 }
337 else
338 {
339 si_->disableSimplexInterface();
340 }
341 }
342
343 /** Compute normalization weights.*/
344 void
computeWeights(CglLandP::LHSnorm norm,CglLandP::Normalization type,CglLandP::RhsWeightType rhs)345 CglLandPSimplex::computeWeights(CglLandP::LHSnorm norm, CglLandP::Normalization type,
346 CglLandP::RhsWeightType rhs)
347 {
348 norm_weights_.clear();
349 norm_weights_.resize(ncols_orig_ ,1.);
350 norm_weights_.resize(ncols_orig_ + nrows_orig_, 0.);
351
352 double * rows_weights = &norm_weights_[ncols_orig_];
353 #ifndef INTEL_COMPILER
354 std::vector<int> nnz(nrows_orig_,0);
355 #else
356 std::vector<int> nnz(nrows_orig_);
357 #endif
358 const CoinPackedMatrix * m = si_->getMatrixByCol();
359 const double * val = m->getElements();
360 const int * ind = m->getIndices();
361 const int * length = m->getVectorLengths();
362 const CoinBigIndex * start = m->getVectorStarts();
363
364 rhs_weight_ = 1;
365
366 if (type== CglLandP::WeightRHS)
367 {
368 if (rhs == CglLandP::Fixed)
369 rhs_weight_ = (ncols_orig_ + 1);//params.rhsWeight;
370 else if (rhs == CglLandP::Dynamic)
371 {
372 throw -1;
373 }
374 }
375
376 if (norm == CglLandP::Infinity)
377 {
378 for (int i = 0 ; i < ncols_orig_ ; i++)
379 {
380 CoinBigIndex begin = start[i];
381 CoinBigIndex end = begin + length[i];
382 for (CoinBigIndex k = begin ; k < end ; k++)
383 {
384 rows_weights[ind[k]] = CoinMax(fabs(val[k]), rows_weights[ind[k]]);
385 rhs_weight_ += fabs(val[k]);
386 nnz[ind[k]] ++;
387 }
388 }
389 }
390 else if (norm == CglLandP::L1 ||
391 norm == CglLandP::Average)
392 {
393 for (int i = 0 ; i < ncols_orig_ ; i++)
394 {
395 CoinBigIndex begin = start[i];
396 CoinBigIndex end = begin + length[i];
397 for (CoinBigIndex k = begin ; k < end ; k++)
398 {
399 rows_weights[ind[k]] += fabs(val[k]);
400 nnz[ind[k]] ++;
401 }
402 }
403 if (norm == CglLandP::Average)
404 {
405 for (int i = 0 ; i < nrows_orig_ ; i++)
406 {
407 rows_weights[i] = static_cast<double>(nnz[i]);
408 }
409 }
410 if (type== CglLandP::WeightBoth)
411 {
412 rhs_weight_ += (ncols_orig_ + 1);
413 std::cout<<"rhs_weight : "<<rhs_weight_<<std::endl;
414 }
415 }
416 else if (norm == CglLandP::L2)
417 {
418 for (int i = 0 ; i < ncols_orig_ ; i++)
419 {
420 CoinBigIndex begin = start[i];
421 CoinBigIndex end = begin + length[i];
422 for (CoinBigIndex k = begin ; k < end ; k++)
423 {
424 rows_weights[ind[k]] += (val[k])*(val[k]);
425 nnz[ind[k]] ++;
426 rhs_weight_ += fabs(val[k]);
427 }
428 }
429 for (int i = 0 ; i < nrows_orig_ ; i++)
430 {
431 rows_weights[i] = sqrt(rows_weights[i]);;
432 }
433 if (type== CglLandP::WeightBoth)
434 {
435 rhs_weight_ = (ncols_orig_ + 1);
436 }
437 }
438 else if (norm == CglLandP::SupportSize)
439 {
440 for (int i = 0 ; i < ncols_orig_ ; i++)
441 {
442 CoinBigIndex begin = start[i];
443 CoinBigIndex end = begin + length[i];
444 for (CoinBigIndex k = begin ; k < end ; k++)
445 {
446 nnz[ind[k]] ++;
447 }
448 }
449
450 for (int i = 0 ; i < nrows_orig_ ; i++)
451 {
452 rows_weights[i] = 1./ static_cast<double> (nnz[i]);
453 }
454
455 if (type== CglLandP::WeightBoth)
456 {
457 rhs_weight_ = (ncols_orig_ + 1);
458 }
459
460 }
461 else if (norm ==CglLandP::Uniform)
462 {
463 for (int i = 0 ; i < nrows_orig_ ; i++)
464 {
465 rows_weights[i] = static_cast<double> (1);
466 }
467 if (type== CglLandP::WeightBoth)
468 {
469 rhs_weight_ = (ncols_orig_ + 1);
470 }
471
472 }
473
474 }
475 void
cacheUpdate(const CglLandP::CachedData & cached,bool reducedSpace)476 CglLandPSimplex::cacheUpdate(const CglLandP::CachedData &cached, bool reducedSpace)
477 {
478 integers_ = cached.integers_;
479 if (own_)
480 {
481 CoinCopyN(cached.basics_, nrows_orig_, basics_);
482 CoinCopyN(cached.nonBasics_, ncols_orig_, nonBasics_);
483 CoinCopyN(cached.colsol_, nrows_orig_ + ncols_orig_, colsol_);
484 for (int i = 0 ; i < ncols_orig_ ; i++)
485 {
486 colsol_[nonBasics_[i]] = 0;
487 }
488 CoinCopyN(cached.colsol_, nrows_orig_ + ncols_orig_, colsolToCut_);
489 //Zero all non basics in colsol setup the reduced space
490 col_in_subspace.resize(0);
491 col_in_subspace.resize(ncols_orig_+nrows_orig_,true);
492 for (int i = 0 ; i < ncols_orig_ ; i++)
493 {
494 setColsolToCut(nonBasics_[i],0.);
495 colsol_[nonBasics_[i]] = 0;
496 }
497 /** Mark the variables at zero in solution to cut so that we know that their contribution to reduced cost has to be computed*/
498 int n_removed =0;
499 if (reducedSpace)
500 {
501 for (int ii = 0; ii < ncols_orig_ ; ii++)
502 {
503 if (getColsolToCut(ii) - up_bounds_[ii] > 1e-08 || getColsolToCut(ii) - lo_bounds_[ii] < 1e-08)
504 {
505 col_in_subspace[ii]=false;
506 n_removed++;
507 }
508 }
509 }
510 }
511 else
512 {
513 basics_ = cached.basics_;
514 nonBasics_ = cached.nonBasics_;
515 }
516 }
517
resetSolver(const CoinWarmStartBasis *)518 bool CglLandPSimplex::resetSolver(const CoinWarmStartBasis * /*basis*/)
519 {
520 si_->disableSimplexInterface();
521 return 0;
522 }
523
524 int
generateExtraCuts(const CglLandP::CachedData & cached,const CglLandP::Parameters & params)525 CglLandPSimplex::generateExtraCuts(const CglLandP::CachedData & cached,
526 const CglLandP::Parameters& params)
527 {
528 int ret_val = 0;
529 for (int i = 0 ; i < nrows_ && cuts_.numberCuts() < params.extraCutsLimit; i++)
530 {
531 if (basics_[i] < ncols_)
532 ret_val += generateExtraCut(i, cached, params);
533 }
534 return ret_val;
535 }
536
537 int
generateExtraCut(int i,const CglLandP::CachedData & cached,const CglLandP::Parameters & params)538 CglLandPSimplex::generateExtraCut(int i, const CglLandP::CachedData & cached,
539 const CglLandP::Parameters& params)
540 {
541 const int & iCol = basics_[i];
542 if (!isInteger(iCol) || int_val(colsol_[iCol], params.away) ||
543 !int_val(getColsolToCut(iCol), params.away) ||
544 colsol_[iCol] < getLoBound(iCol) || colsol_[iCol] > getUpBound(iCol) ||
545 (cuts_.rowCut(iCol) != NULL) )
546 {
547 return false;
548 }
549
550 #ifdef DBG_OUT
551 printf("var: %i, basic in row %i. integer=%i, colsol_=%f, colsolToCut_=%f.\n",iCol, i,
552 isInteger(iCol), colsol_[iCol],
553 getColsolToCut(iCol));
554
555 printf("generating extra cut....\n");
556 #endif
557
558 OsiRowCut * cut = new OsiRowCut;
559 generateMig(i, *cut, params);
560 assert(fabs(row_k_.rhs - colsol_[iCol]) < 1e-10);
561
562 int code = validator_(*cut, cached.colsol_, *si_, params, &lo_bounds_[0], &up_bounds_[0]);
563 if (code)
564 {
565 delete cut;
566 return false;
567 }
568 else
569 {
570 cuts_.insert(iCol, cut);
571 return true;
572 }
573 }
574
575
576 void
genThisBasisMigs(const CglLandP::CachedData & cached,const CglLandP::Parameters & params)577 CglLandPSimplex::genThisBasisMigs(const CglLandP::CachedData &cached,
578 const CglLandP::Parameters & params)
579 {
580 for (int i = 0 ; i < cached.nBasics_ ; i++)
581 {
582 const int iCol = basics_[i];
583 if (iCol >= ncols_ ||
584 !cached.integers_[iCol] ||
585 int_val(colsol_[iCol], params.away))
586 continue;
587 OsiRowCut * cut = new OsiRowCut;
588 generateMig(i, *cut, params);
589 int code = validator_(*cut, cached.colsol_, *si_, params, &lo_bounds_[0], &up_bounds_[0]);
590 if (code)
591 {
592 delete cut;
593 continue;
594 }
595 cut->setEffectiveness(cut->violated(cached.colsol_));
596 if (cuts_.rowCut(iCol) == NULL || cut->effectiveness() > cuts_.rowCut(iCol)->effectiveness())
597 {
598 cuts_.insert(iCol,cut);
599 }
600 else
601 delete cut;
602 }
603 }
604
605
606 bool
generateMig(int row,OsiRowCut & cut,const CglLandP::Parameters & params)607 CglLandPSimplex::generateMig(int row, OsiRowCut & cut,
608 const CglLandP::Parameters & params)
609 {
610 row_k_.num = row;
611 pullTableauRow(row_k_);
612 row_k_.rhs = row_k_.rhs - floor(row_k_.rhs);
613 if (params.strengthen || params.modularize)
614 createMIG(row_k_, cut);
615 else
616 createIntersectionCut(row_k_, cut);
617
618 return 1;//At this point nothing failed, always generate a cut
619 }
620
621 bool
optimize(int row,OsiRowCut & cut,const CglLandP::CachedData & cached,const CglLandP::Parameters & params)622 CglLandPSimplex::optimize
623 (int row, OsiRowCut & cut,const CglLandP::CachedData &cached,const CglLandP::Parameters & params)
624 {
625 bool optimal = false;
626 int nRowFailed = 0;
627
628 double timeLimit = CoinMin(params.timeLimit, params.singleCutTimeLimit);
629 timeLimit += CoinCpuTime();
630 // double timeBegin = CoinCpuTime();
631 /** Copy the cached information */
632 nrows_ = nrows_orig_;
633 ncols_ = ncols_orig_;
634 CoinCopyN(cached.basics_, nrows_, basics_);
635 CoinCopyN(cached.nonBasics_, ncols_, nonBasics_);
636 CoinCopyN(cached.colsol_, nrows_+ ncols_, colsol_);
637 CoinCopyN(cached.colsol_, nrows_+ ncols_, colsolToCut_);
638
639 delete basis_;
640 basis_ = new CoinWarmStartBasis(*cached.basis_);
641 #define CACHED_SOLVER
642 #ifndef CACHED_SOLVER
643 si_->enableSimplexInterface(0);
644 #else
645 delete si_;
646 si_ = cached.solver_->clone();
647 #ifdef COIN_HAS_OSICLP
648 OsiClpSolverInterface * clpSi = dynamic_cast<OsiClpSolverInterface *>(si_);
649 OsiClpSolverInterface * clpSiRhs = dynamic_cast<OsiClpSolverInterface *>(cached.solver_);
650 if (clpSi)
651 {
652 clp_ = clpSi;
653 clpSi->getModelPtr()->copyEnabledStuff(clpSiRhs->getModelPtr());;
654 }
655 #endif
656 #endif
657 #ifdef APPEND_ROW
658 if (params.modularize)
659 {
660 append_row(row, params.modularize);
661 row_k_.modularize(integers_);
662 }
663 #endif
664 for (int i = 0 ; i < ncols_; i++)
665 {
666 setColsolToCut(nonBasics_[i], 0.);
667 colsol_[nonBasics_[i]] = 0;
668 }
669
670 #ifdef APPEND_ROW
671 if (params.modularize)
672 row_k_.num = nrows_ - 1;
673 else
674 #endif
675 row_k_.num = row;
676
677 pullTableauRow(row_k_);
678 row_k_.rhs = row_k_.rhs - floor(row_k_.rhs);
679
680 if (params.modularize)
681 row_k_.modularize(integers_);
682
683 updateM1_M2_M3(row_k_, 0., params.perturb);
684 sigma_ = computeCglpObjective(row_k_);
685
686 handler_->message(Separating,messages_)<<basics_[row]<<sigma_<<CoinMessageEol<<CoinMessageEol;
687 handler_->message(LogHead, messages_)<<CoinMessageEol<<CoinMessageEol;
688
689 //Save the variable basic in this row
690 // int var_k = cached.basics_[k_];
691
692 //Put a flag on each row to say if we want to continue trying to use it
693 CoinFillN(rowFlags_,nrows_,true);
694
695 int numberConsecutiveDegenerate = 0;
696 bool allowDegeneratePivot = numberConsecutiveDegenerate < params.degeneratePivotLimit;
697 bool beObstinate = 0;
698 int numPivots = 0;
699 int saveNumSourceEntered = numSourceRowEntered_;
700 int saveNumIncreased = numIncreased_;
701 int numCycle = 0;
702 int numFailedPivots = 0;
703 bool hasFlagedRow = false;
704 int maxTryRow = 5;
705 while ( !optimal && numPivots < params.pivotLimit)
706 {
707 if (timeLimit - CoinCpuTime() < 0.) break;
708
709 updateM1_M2_M3(row_k_, 0., params.perturb);
710 sigma_ = computeCglpObjective(row_k_);
711 int direction = 0;
712 int gammaSign = 0;
713 int leaving = -1;
714 int incoming = -1;
715 double bestSigma;
716 if (params.pivotSelection != CglLandP::initialReducedCosts || numPivots == 0)
717 {
718 leaving = fastFindCutImprovingPivotRow(direction, gammaSign, params.pivotTol,
719 params.pivotSelection == CglLandP::initialReducedCosts);
720 #if 0
721 plotCGLPobj(direction, params.pivotTol, params.pivotTol, true, true, false);
722 exit(1);
723 #endif
724 if (leaving >= 0)
725 {
726 if (params.pivotSelection == CglLandP::mostNegativeRc ||
727 (params.pivotSelection == CglLandP::initialReducedCosts && numPivots == 0))
728 {
729
730 if (params.pivotSelection == CglLandP::initialReducedCosts)
731 rowFlags_[leaving] = false;
732 incoming = fastFindBestPivotColumn(direction, gammaSign,
733 params.pivotTol, params.away,
734 (params.sepSpace==CglLandP::Fractional),
735 allowDegeneratePivot,
736 bestSigma, false//params.modularize
737 );
738 while (incoming < 0 && !optimal &&
739 nRowFailed < maxTryRow) // if no improving was found rescan the tables of reduced cost to find a good one
740 {
741 if (incoming == -1 || params.countMistakenRc) nRowFailed ++;
742 rowFlags_[leaving] = false;
743 hasFlagedRow = true;
744 leaving = rescanReducedCosts(direction, gammaSign, params.pivotTol);
745 if (leaving >= 0)
746 {
747 incoming = fastFindBestPivotColumn(direction, gammaSign,
748 params.pivotTol,
749 params.away,
750 (params.sepSpace==CglLandP::Fractional),
751 allowDegeneratePivot,
752 bestSigma, false//params.modularize
753 );
754 }
755 else optimal = true;
756 }
757 }
758 else if (params.pivotSelection == CglLandP::bestPivot)
759 {
760 incoming = findBestPivot(leaving, direction, params);
761 }
762 }
763 }
764 else if (params.pivotSelection == CglLandP::initialReducedCosts)
765 {
766 assert(numPivots > 0);
767 while (incoming < 0 && !optimal) // if no improving was found rescan the tables of reduced cost to find a good one
768 {
769 if (!hasFlagedRow)
770 hasFlagedRow = true;
771 leaving = rescanReducedCosts(direction, gammaSign, params.pivotTol);
772 rowFlags_[leaving] = false;
773 if (leaving >= 0)
774 {
775 incoming = fastFindBestPivotColumn(direction, gammaSign,
776 params.pivotTol,
777 params.away,
778 (params.sepSpace==CglLandP::Fractional),
779 allowDegeneratePivot,
780 bestSigma, false//params.modularize
781 );
782 }
783 else optimal = true;
784 }
785 }
786 if (leaving >= 0)
787 {
788 if ( incoming >= 0 && !optimal)
789 {
790 if (inDegenerateSequence_) //flag leaving row
791 {
792 numberConsecutiveDegenerate++;
793 allowDegeneratePivot = numberConsecutiveDegenerate < params.degeneratePivotLimit;
794 rowFlags_[leaving] = false;
795 }
796 else
797 {
798 beObstinate = 0;
799 numberConsecutiveDegenerate = 0;
800 allowDegeneratePivot = numberConsecutiveDegenerate < params.degeneratePivotLimit;
801 }
802 double gamma = - row_k_[nonBasics_[incoming]] / row_i_[nonBasics_[incoming]];
803 #ifdef COIN_HAS_OSICLP
804 if (numPivots && ( numPivots % 40 == 0 ) && clp_)
805 {
806 clp_->getModelPtr()->factorize();
807 }
808 #endif
809
810 #ifndef OLD_COMPUTATION
811 bool recompute_source_row = (numPivots && (numPivots % 10 == 0 ||
812 fabs(gamma) < 1e-05));
813 #endif
814
815 std::pair<int, int> cur_pivot(nonBasics_[incoming],basics_[leaving]);
816
817 bool pivoted = changeBasis(incoming,leaving,direction,
818 #ifndef OLD_COMPUTATION
819 recompute_source_row,
820 #endif
821 false);//params.modularize);
822
823 if (leaving == row)
824 {
825 numSourceRowEntered_++;
826 }
827
828 if (params.generateExtraCuts == CglLandP::WhenEnteringBasis &&
829 basics_[leaving] < ncols_ && cuts_.numberCuts() < params.extraCutsLimit)
830 generateExtraCut(leaving, cached, params);
831 if (pivoted)
832 {
833 numPivots++;
834
835 double lastSigma = sigma_;
836 if (params.modularize)
837 {
838 row_k_.modularize(integers_);
839 }
840 sigma_ = computeCglpObjective(row_k_);
841
842 if (sigma_ - lastSigma > -1e-4*(lastSigma))
843 {
844 if(sigma_ > 0) return 0;
845 #if 0
846 if (sigma_ > 0 || sigma_ - lastSigma > 1e1*(-lastSigma))
847 return 0;
848 #endif
849 }
850
851
852 handler_->message(PivotLog,messages_)<<numPivots<<sigma_<<
853 nonBasics_[incoming]<<basics_[leaving]<<direction<<gamma<<inDegenerateSequence_<<CoinMessageEol<<CoinMessageEol;
854 }
855 else //pivot failed
856 {
857 numFailedPivots++;
858 //check wether sigma has changed if it has exit cut generation if it has not continue
859 double lastSigma = sigma_;
860 sigma_ = computeCglpObjective(row_k_);
861 if ( sigma_-lastSigma>1e-8)
862 {
863 handler_->message(PivotFailedSigmaIncreased,messages_)<<CoinMessageEol<<CoinMessageEol;
864 //break;
865 return 0;
866 }
867 handler_->message(PivotFailedSigmaUnchanged,messages_)<<CoinMessageEol<<CoinMessageEol;
868 numFailedPivots = params.failedPivotLimit + 1;
869 return 0;
870 if (numFailedPivots > params.failedPivotLimit)
871 break;
872
873 }
874 }
875 else //attained max number of leaving vars tries with no improvement
876 {
877 handler_->message(WarnGiveUpRow,messages_)<<nRowFailed<<CoinMessageEol<<CoinMessageEol;
878 break;
879 }
880 }
881 else
882 {
883 if (hasFlagedRow && beObstinate)
884 {
885 //Reset row flags
886 CoinFillN(rowFlags_,nrows_,true);
887 hasFlagedRow = false;
888 if (inDegenerateSequence_)
889 {
890 allowDegeneratePivot = false;
891 beObstinate = false;
892 }
893 }
894 else
895 {
896 //could perturb but Ionut skipped that will see later
897 optimal = true;
898 handler_->message(FinishedOptimal, messages_)<<sigma_<<numPivots<<CoinMessageEol<<CoinMessageEol;
899 }
900 }
901 }
902
903 if (!optimal && numPivots >= params.pivotLimit)
904 {
905 std::string limit="pivots";
906 handler_->message(HitLimit, messages_)<<limit<<numPivots<<CoinMessageEol<<CoinMessageEol;
907 }
908 if (!optimal && numFailedPivots >= params.failedPivotLimit)
909 {
910 std::string limit="failed pivots";
911 handler_->message(HitLimit, messages_)<<limit<<numPivots<<CoinMessageEol<<CoinMessageEol;
912 }
913 //Create the cut
914
915 //pullTableauRow(row_k_);
916 //row_k_.rhs = row_k_.rhs - floor(row_k_.rhs);
917
918 // double normalization = 100*normCoef(row_k_);
919 {
920 if (params.strengthen || params.modularize)
921 createMIG(row_k_, cut);
922 else
923 createIntersectionCut(row_k_, cut);
924 }
925
926 if (params.generateExtraCuts == CglLandP::AtOptimalBasis)
927 {
928 generateExtraCuts(cached, params);
929 }
930 handler_->message(CutStat, messages_)<<row<<numPivots
931 <<numSourceRowEntered_ - saveNumSourceEntered
932 <<numIncreased_- saveNumIncreased
933 <<numCycle<<CoinMessageEol;
934 return 1;//At this point nothing failed, always generate a cut
935 }
936
937
938 bool
changeBasis(int incoming,int leaving,int leavingStatus,bool recompute_source_row,bool modularize)939 CglLandPSimplex::changeBasis(int incoming, int leaving, int leavingStatus,
940 #ifndef OLD_COMPUTATION
941 bool recompute_source_row,
942 #endif
943 bool modularize)
944 {
945 double infty = si_->getInfinity();
946 int clpLeavingStatus = leavingStatus;
947
948 #ifdef COIN_HAS_OSICLP
949 if (clp_)
950 {
951 if (basics_[leaving] >= ncols_)
952 clpLeavingStatus = - leavingStatus;
953 }
954 #endif
955
956 int code = 0;
957
958 code = si_->pivot(nonBasics_[incoming],basics_[leaving], clpLeavingStatus);
959 if (code)
960 {
961 #ifdef OLD_COMPUTATION
962 if (!modularize)
963 {
964 pullTableauRow(row_k_);
965 row_k_.rhs = row_k_.rhs - floor(row_k_.rhs);
966 }
967 else{
968 int & indexLeaving = basics_[leaving];
969 if (leavingStatus==1)
970 {
971 setColsolToCut(indexLeaving, getUpBound(indexLeaving) - getColsolToCut(indexLeaving));
972 }
973 else
974 {
975 setColsolToCut(indexLeaving, getColsolToCut(indexLeaving) + getLoBound(indexLeaving));
976 }
977 }
978 #endif
979 return 0;
980 }
981 numPivots_ ++;
982 //swap bounds
983 int & indexLeaving = basics_[leaving];
984 #ifdef OLD_COMPUTATION
985 if (!modularize)
986 {
987 if (leavingStatus==1)
988 {
989 setColsolToCut(indexLeaving, getUpBound(indexLeaving) - getColsolToCut(indexLeaving));
990 }
991 else
992 {
993 setColsolToCut(indexLeaving, getColsolToCut(indexLeaving) - getLoBound(indexLeaving));
994 }
995 }
996 #endif
997
998 if (indexLeaving < ncols_)
999 {
1000 basis_->setStructStatus(indexLeaving, leavingStatus==1 ? CoinWarmStartBasis::atUpperBound : CoinWarmStartBasis::atLowerBound);
1001 }
1002 else
1003 {
1004 int iRow = basics_[leaving] - ncols_;
1005 basis_->setArtifStatus(iRow, leavingStatus==1 ? CoinWarmStartBasis::atUpperBound : CoinWarmStartBasis::atLowerBound);
1006 // assert(leavingStatus==-1 || (rowLower_[iRow]>-1e50 && rowUpper_[iRow] < 1e50));
1007 }
1008
1009 if (nonBasics_[incoming] < ncols_)
1010 {
1011 int & indexIncoming = nonBasics_[incoming];
1012 CoinWarmStartBasis::Status status = basis_->getStructStatus(indexIncoming);
1013 if (status==CoinWarmStartBasis::atUpperBound)
1014 setColsolToCut(indexIncoming, getUpBound(indexIncoming) - getColsolToCut(indexIncoming));
1015 else
1016 setColsolToCut(indexIncoming, getColsolToCut(indexIncoming) + getLoBound(indexIncoming));
1017 basis_->setStructStatus(indexIncoming, CoinWarmStartBasis::basic);
1018 }
1019 else
1020 {
1021 int iRow = nonBasics_[incoming] - ncols_;
1022 int & indexIncoming = nonBasics_[incoming];
1023
1024 if (basis_->getArtifStatus(iRow)==CoinWarmStartBasis::atUpperBound)
1025 setColsolToCut(indexIncoming, getUpBound(indexIncoming) - getColsolToCut(indexIncoming));
1026 else
1027 setColsolToCut(indexIncoming, getColsolToCut(indexIncoming) + getLoBound(indexIncoming));
1028
1029 basis_->setArtifStatus(iRow, CoinWarmStartBasis::basic);
1030 }
1031
1032 int swap = basics_[leaving];
1033 basics_[leaving] = nonBasics_[incoming];
1034 nonBasics_[incoming] = swap;
1035 //update solution of leaving variable
1036 colsol_[nonBasics_[incoming]] = 0;
1037
1038 //update solution for basics
1039 const double * lpSol = si_->getColSolution();
1040 const double * rowAct = si_->getRowActivity();
1041 const double * rowLower = si_->getRowLower();
1042 const double * rowUpper = si_->getRowUpper();
1043
1044 for (int i = 0 ; i < nrows_ ; i++)
1045 {
1046 int& iCol = basics_[i];
1047 if (iCol<ncols_)
1048 colsol_[iCol] = lpSol[iCol];
1049 else // Osi does not give direct acces to the value of slacks
1050 {
1051 int iRow = iCol - ncols_;
1052 colsol_[iCol] = - rowAct[iRow];
1053 if (rowLower[iRow]> -infty)
1054 {
1055 colsol_[iCol] += rowLower[iRow];
1056 }
1057 else
1058 {
1059 colsol_[iCol] += rowUpper[iRow];
1060 }
1061 }
1062 }
1063
1064 // basics_ may unfortunately change reload
1065 int k = basics_[row_k_.num];
1066 si_->getBasics(basics_);
1067 if (basics_[row_k_.num] != k)
1068 {
1069 for (int ii = 0 ; ii < nrows_ ; ii++)
1070 {
1071 if (basics_[ii]==k)
1072 {
1073 row_k_.num= ii;
1074 break;
1075 }
1076 }
1077 }
1078
1079 #ifndef OLD_COMPUTATION
1080 if (!modularize && recompute_source_row)
1081 {
1082 #else
1083 if (!modularize)
1084 {
1085 #endif
1086 pullTableauRow(row_k_);
1087 row_k_.rhs = row_k_.rhs - floor(row_k_.rhs);
1088 }
1089 #if 0
1090 }
1091 #endif
1092 else //Update row k by hand
1093 {
1094 double gamma = - row_k_[basics_[leaving]] / row_i_[basics_[leaving]];
1095 row_k_[basics_[leaving]] = 0;
1096 row_k_.quickAdd(nonBasics_[incoming], gamma);
1097 if(1 || fabs(gamma) > 1e-9){
1098 int nnz = row_i_.getNumElements();
1099 const int * indices = row_i_.getIndices();
1100 for (int i = 0 ; i < nnz; i++)
1101 {
1102 if(row_k_.getNumElements() > row_k_.capacity() - 2) row_k_.scan();
1103 if (indices[i] != nonBasics_[incoming] && indices[i] != basics_[leaving])
1104 row_k_.quickAdd(indices[i], gamma * row_i_[indices[i]]);
1105 }
1106 row_k_.rhs += gamma * row_i_.rhs;
1107 }
1108 row_k_.scan();
1109 row_k_.clean(1e-10);
1110 checkVec(row_k_);
1111 #if 0
1112 TabRow test_row(this);
1113 int rowsize = ncols_orig_ + nrows_orig_ + 1;
1114 test_row.reserve(rowsize);
1115 test_row.num = row_k_.num;
1116 pullTableauRow(test_row);
1117 test_row.rhs = test_row.rhs - floor(test_row.rhs);
1118 VecModEqAssert(row_k_, test_row);
1119 #endif
1120 }
1121
1122 return true;
1123 }
1124
1125 /** Find a row which can be used to perform an improving pivot return index of the cut or -1 if none exists
1126 * (i.e., find the leaving variable).*/
1127 int
findCutImprovingPivotRow(int & direction,int & gammaSign,double tolerance)1128 CglLandPSimplex::findCutImprovingPivotRow( int &direction, int &gammaSign, double tolerance)
1129 {
1130 bool bestRed = 0;
1131 tolerance = -10*tolerance;
1132 int bestRow = -1;
1133 int bestDirection = 0;
1134 int bestGamma = 0;
1135 double infty = si_->getInfinity();
1136 for (row_i_.num = 0 ; row_i_.num < nrows_; row_i_.num++)
1137 {
1138 //if ( (row_k_.modularized_ || row_i_.num != row_k_.num)//obviously not necessary to combine row k with itself (unless modularized)
1139 if ( (row_i_.num != row_k_.num)//obviously not necessary to combine row k with itself (unless modularized)
1140 && rowFlags_[row_i_.num] //row has not been flaged
1141 // && fabs(getUpBound(basics_[row_i_.num]) - getLoBound(basics_[row_i_.num]))>1e-09 //variable is not fixed
1142 )
1143 {
1144 pullTableauRow(row_i_);
1145 double tau = computeRedCostConstantsInRow();
1146
1147 if (getLoBound(basics_[row_i_.num]) > -infty)
1148 // variable can leave at its lower bound
1149 //Compute reduced cost with basics_[i] value decreasing
1150 {
1151 direction = -1;
1152
1153 gammaSign = -1;
1154 double redCost = computeCglpRedCost(direction, gammaSign, tau);
1155 if (redCost<tolerance)
1156 {
1157 if (bestRed)
1158 {
1159 tolerance = redCost;
1160 bestRow = row_i_.num;
1161 bestDirection = direction;
1162 bestGamma = gammaSign;
1163 }
1164 else return row_i_.num;
1165 }
1166 gammaSign = 1;
1167 redCost = computeCglpRedCost(direction, gammaSign, tau);
1168 if (redCost<tolerance)
1169 {
1170 if (bestRed)
1171 {
1172 tolerance = redCost;
1173 bestRow = row_i_.num;
1174 bestDirection = direction;
1175 bestGamma = gammaSign;
1176 }
1177 else return row_i_.num;
1178 }
1179 }
1180 if ( getUpBound(basics_[row_i_.num])<infty) // variable can leave at its upper bound
1181 //Compute reduced cost with basics_[i] value decreasing
1182 {
1183 direction = 1;
1184 // adjustTableauRow(i_, row_i_, rhs_i_, direction);
1185 gammaSign = -1;
1186 double redCost = computeCglpRedCost(direction, gammaSign, tau);
1187 if (redCost<tolerance)
1188 {
1189 if (bestRed)
1190 {
1191 tolerance = redCost;
1192 bestRow = row_i_.num;
1193 bestDirection = direction;
1194 bestGamma = gammaSign;
1195 }
1196 else return row_i_.num;
1197 }
1198 gammaSign = 1;
1199 redCost = computeCglpRedCost(direction, gammaSign, tau);
1200 if (redCost<tolerance)
1201 {
1202 if (bestRed)
1203 {
1204 tolerance = redCost;
1205 bestRow = row_i_.num;
1206 bestDirection = direction;
1207 bestGamma = gammaSign;
1208 }
1209 else return row_i_.num;
1210 }
1211 }
1212 rowFlags_[row_i_.num]=false;
1213 }
1214 }
1215 direction = bestDirection;
1216 gammaSign = bestGamma;
1217 row_i_.num=bestRow;
1218 if (row_i_.num >=0 && row_i_.num != bestRow)
1219 {
1220 row_i_.num=bestRow;
1221 pullTableauRow(row_i_);
1222 }
1223 assert (bestRow<0||direction!=0);
1224 return bestRow;
1225 }
1226
1227
1228 /** Find a row which can be used to perform an improving pivot the fast way
1229 * (i.e., find the leaving variable).
1230 \return index of the cut or -1 if none exists. */
1231 int
fastFindCutImprovingPivotRow(int & direction,int & gammaSign,double tolerance,bool flagPositiveRows)1232 CglLandPSimplex::fastFindCutImprovingPivotRow( int &direction, int &gammaSign,
1233 double tolerance, bool flagPositiveRows)
1234 {
1235 bool modularize = false;
1236 double sigma = sigma_ /rhs_weight_;
1237 //Fill vector to compute contribution to reduced cost of variables in M1 and M2 (nz non-basic vars in row_k_.row).
1238 // 1. Put the values
1239 // 2. Post multiply by basis inverse
1240 double * rWk1bis_ =NULL;
1241 CoinFillN(&rWk1_[0],nrows_,static_cast<double> (0));
1242 if (modularize)
1243 CoinFillN(rWk1bis_, nrows_, static_cast<double> (0));
1244 int capacity = 0;
1245 const CoinPackedMatrix* mat = si_->getMatrixByCol();
1246
1247 const CoinBigIndex* starts = mat->getVectorStarts();
1248 const int * lengths = mat->getVectorLengths();
1249 const int * indices = mat->getIndices();
1250 const double * elements = mat->getElements();
1251
1252 for (unsigned int i = 0 ; i < M1_.size() ; i++)
1253 {
1254 const int& ii = M1_[i];
1255 if (ii < ncols_)
1256 {
1257 const CoinBigIndex& begin = starts[ii];
1258 const CoinBigIndex end = begin + lengths[ii];
1259 bool swap = false;
1260 if (basis_->getStructStatus(ii)==CoinWarmStartBasis::atUpperBound) swap = true;
1261 for (CoinBigIndex k = begin ; k < end ; k++)
1262 {
1263 if (swap)
1264 {
1265 rWk1_[indices[k]] += normedCoef(elements[k] * sigma, ii);
1266 if (modularize)
1267 {
1268 rWk1bis_[indices[k]] += elements[k] * (getColsolToCut(ii) - normedCoef(sigma, ii));
1269 }
1270
1271 }
1272 else
1273 {
1274 rWk1_[indices[k]] -= normedCoef(elements[k] * sigma, ii);
1275 if (modularize)
1276 {
1277 rWk1bis_[indices[k]] -= elements[k] * (getColsolToCut(ii) - normedCoef(sigma, ii));
1278 }
1279 }
1280 }
1281 }
1282 else
1283 {
1284 bool swap = false;
1285 if (basis_->getArtifStatus(ii - ncols_orig_)==CoinWarmStartBasis::atUpperBound) swap = true;
1286 if (swap)
1287 {
1288 rWk1_[ii - ncols_] += normedCoef(sigma, ii);
1289 if (modularize)
1290 {
1291 rWk1bis_[ii - ncols_] += (getColsolToCut(ii) - normedCoef(sigma, ii));
1292 }
1293 }
1294 else
1295 {
1296 rWk1_[ii - ncols_] -= normedCoef(sigma, ii);
1297 if (modularize)
1298 {
1299 rWk1bis_[ii - ncols_] -= (getColsolToCut(ii) - normedCoef(sigma, ii));
1300 }
1301 }
1302 }
1303 }
1304 for (unsigned int i = 0 ; i < M2_.size(); i++)
1305 {
1306 const int& ii = M2_[i];
1307 if (ii<ncols_)
1308 {
1309 const CoinBigIndex& begin = starts[ii];
1310 const CoinBigIndex end = begin + lengths[ii];
1311 bool swap = false;
1312 if (basis_->getStructStatus(ii)==CoinWarmStartBasis::atUpperBound) swap = true;
1313 for (CoinBigIndex k = begin ; k < end ; k++)
1314 {
1315 if (swap)
1316 {
1317 rWk1_[indices[k]] += elements[k] * (getColsolToCut(ii) - normedCoef(sigma, ii));
1318 if (modularize)
1319 rWk1bis_[indices[k]] += elements[k] * normedCoef(sigma, ii);
1320 }
1321 else
1322 {
1323 rWk1_[indices[k]] -= elements[k] * (getColsolToCut(ii) - normedCoef(sigma, ii));
1324 if (modularize)
1325 rWk1bis_[indices[k]] -= elements[k] * normedCoef(sigma, ii);
1326 }
1327 }
1328 }
1329 else
1330 {
1331 bool swap = false;
1332 if (basis_->getArtifStatus(M2_[i] - ncols_orig_)==CoinWarmStartBasis::atUpperBound) swap = true;
1333 if (swap)
1334 {
1335 rWk1_[ii - ncols_] += (getColsolToCut(ii) - normedCoef(sigma, ii));
1336 if (modularize)
1337 rWk1bis_[ii - ncols_] += normedCoef(sigma, ii);
1338 }
1339 else
1340 {
1341 rWk1_[ii - ncols_] -= (getColsolToCut(ii) - normedCoef(sigma, ii));
1342 if (modularize)
1343 rWk1bis_[ii - ncols_] -= normedCoef(sigma, ii);
1344 }
1345 }
1346 }
1347
1348 for (int i = 0 ; i < nrows_ ; i++)
1349 {
1350 if (rWk1_[i])
1351 rIntWork_[capacity++] = i;
1352 }
1353 CoinIndexedVector indexed;
1354 indexed.borrowVector(nrows_, capacity, &rIntWork_[0], &rWk1_[0]);
1355
1356 #ifdef COIN_HAS_OSICLP
1357 if (clp_)
1358 clp_->getBInvACol(&indexed);
1359 else
1360 #endif
1361 throw CoinError("Function not implemented in this OsiSolverInterface",
1362 "getBInvACol","CglLandpSimplex");
1363 indexed.returnVector();
1364 if (modularize)
1365 {
1366 capacity = 0;
1367 for (int i = 0 ; i < nrows_ ; i++)
1368 {
1369 if (rWk1bis_[i])
1370 rIntWork_[capacity++] = i;
1371 }
1372
1373 indexed.borrowVector(nrows_, capacity, &rIntWork_[0], rWk1bis_);
1374 #ifdef COIN_HAS_OSICLP
1375 if (clp_)
1376 clp_->getBInvACol(&indexed);
1377 else
1378 #endif
1379 indexed.returnVector();
1380 }
1381 //Now compute the contribution of the variables in M3_
1382 //Need to get the column of the tableau in rW3_ for each of these and
1383 //add up with correctly in storage for multiplier for negative gamma (named rW3_) and
1384 //for positive gamma (which is named rW4_)
1385 if (!M3_.empty())
1386 {
1387 CoinFillN(&rWk3_[0],nrows_,0.);
1388 CoinFillN(&rWk4_[0],nrows_,0.);
1389 if (modularize)
1390 {
1391 double * rWk3bis_ = NULL;
1392 double * rWk4bis_ = NULL;
1393 CoinFillN(rWk3bis_,nrows_,0.);
1394 CoinFillN(rWk4bis_,nrows_,0.);
1395 }
1396 }
1397 for (unsigned int i = 0 ; i < M3_.size() ; i++)
1398 {
1399 const int & ii = M3_[i];
1400 si_->getBInvACol(ii, &rWk2_[0]);
1401 bool swap = false;
1402 if (ii < ncols_orig_ && basis_->getStructStatus(ii)==CoinWarmStartBasis::atUpperBound) swap = true;
1403 if (ii >= ncols_orig_ && basis_->getArtifStatus(ii - ncols_orig_)==CoinWarmStartBasis::atUpperBound) swap = true;
1404
1405 for (int j = 0 ; j < nrows_ ; j++)
1406 {
1407 if (swap)
1408 rWk2_[j] = - rWk2_[j];
1409 if (rWk2_[j] > 0.)
1410 {
1411 //is in M1 for multiplier with negative gamma
1412 rWk3_[j] -= normedCoef(sigma*rWk2_[j], ii);
1413 //is in M2 for multiplier with positive gamma
1414 rWk4_[j] -= (getColsolToCut(ii) - normedCoef(sigma, ii))*rWk2_[j];
1415 }
1416 else if (rWk2_[j] < 0.)
1417 {
1418 //is in M2 for multiplier with negative gamma
1419 rWk3_[j] -=(getColsolToCut(ii) - normedCoef(sigma, ii))*rWk2_[j];
1420
1421 //is in M1 for multiplier with positive gamma
1422 rWk4_[j] -= normedCoef(sigma, ii)*rWk2_[j];
1423 }
1424 }
1425
1426 }
1427 //Now, we just need to add up everything correctly for each of the reduced
1428 //cost. Compute the Tau in rWk2_ which is not used anymore then compute the reduced cost in rWk1_ for u^l_j
1429 // rwk2_ in u^u_j rWk3_ for v^u_j rWk4_ for v^u_j
1430 // Let's rename not to get too much confused
1431 double * ul_i = &rWk1_[0];
1432 double * uu_i = &rWk2_[0];
1433 double * vl_i = &rWk3_[0];
1434 double * vu_i = &rWk4_[0];
1435 int bestRow = -1;
1436 int bestDirection = 0;
1437 int bestGammaSign = 0;
1438 // double infty = si_->getInfinity();
1439
1440
1441 nNegativeRcRows_ = 0;//counter
1442 int nZeroRc = 0;
1443 int nPositiveRc = 0;
1444
1445 double fzero = getColsolToCut(basics_[row_k_.num]) - floor(getColsolToCut(basics_[row_k_.num]));
1446 // fzero = row_k_.rhs;
1447 //for (int i = 0 ; i < ncols_orig_ ; i++) {
1448 // fzero -= getColsolToCut(nonBasics_[i]) * row_k_[nonBasics_[i]];
1449 //}
1450
1451 double bestReducedCost = -tolerance;
1452 for (int i = 0 ; i < nrows_ ; i++)
1453 {
1454 //if ((!row_k_.modularized_ && i == row_k_.num)//obviously not necessary to combine row k with itself
1455 if ((i == row_k_.num)//obviously not necessary to combine row k with itself
1456 // && fabs(getUpBound(basics_[row_i_.num]) - getLoBound(basics_[row_i_.num]))>1e-09 //variable is not fixed
1457 || col_in_subspace[basics_[i]] == false
1458 )
1459 {
1460 ul_i[i]=uu_i[i]=vl_i[i]=vu_i[i]=10.;
1461 rowFlags_[i] = false;
1462 continue;
1463 }
1464
1465 double tau1 = rWk1_[i];
1466 double tau2 = rWk1_[i];
1467 double tau3 = rWk1_[i];
1468 double tau4 = rWk1_[i];
1469 if (!M3_.empty())
1470 {
1471 tau1 += rWk3_[i];
1472 tau2 += rWk4_[i];
1473 tau3 += rWk4_[i];
1474 tau4 += rWk3_[i];
1475 }
1476 if (modularize)
1477 {
1478 tau1 = rWk1_[i] + rWk3_[i];
1479 tau2 = rWk1bis_[i] + rWk3_[i];
1480 tau3 = - rWk1_[i] - rWk4_[i];
1481 tau4 = - rWk1bis_[i] - rWk3_[i];
1482 }
1483
1484
1485 double redCost;
1486 bool hasNegativeRc = false;
1487
1488 double loBound = getLoBound(basics_[i]);
1489
1490 if (loBound > -1e50)
1491 {
1492 redCost = - normedCoef(sigma,basics_[i]) + (tau1)
1493 + (1 - fzero) * ( colsol_[basics_[i]]
1494 - loBound);
1495
1496
1497 if (redCost < -tolerance)
1498 {
1499 ul_i[i] = redCost;
1500 hasNegativeRc = true;
1501 }
1502 else
1503 {
1504 if (fabs(redCost) < tolerance) nZeroRc++;
1505 else nPositiveRc ++;
1506 ul_i[i] = 10.;
1507 }
1508 if (redCost < bestReducedCost
1509 && rowFlags_[i] ) //row has not been flaged
1510 {
1511 bestDirection = -1;
1512 bestGammaSign = -1;
1513 bestReducedCost = redCost;
1514 bestRow = i;
1515 }
1516
1517
1518 redCost = -normedCoef(sigma,basics_[i]) - (tau2)
1519 - (1 - fzero) * ( colsol_[basics_[i]]
1520 - loBound)
1521 - loBound + getColsolToCut(basics_[i]);
1522
1523 if (redCost < -tolerance)
1524 {
1525 vl_i[i] = redCost;
1526 hasNegativeRc = true;
1527 }
1528 else
1529 {
1530 if (fabs(redCost) < tolerance) nZeroRc++;
1531 else nPositiveRc ++;
1532 vl_i[i]=10.;
1533 }
1534
1535 if (redCost < bestReducedCost
1536 && rowFlags_[i]) //row has not been flaged
1537 {
1538 bestDirection = -1;
1539 bestGammaSign = 1;
1540 bestReducedCost = redCost;
1541 bestRow = i;
1542 }
1543
1544
1545
1546 }
1547 else
1548 {
1549 ul_i[i] = 10.;
1550 vl_i[i] = 10.;
1551 }
1552 double upBound = getUpBound(basics_[i]);
1553 if (getUpBound(basics_[i]) < 1e50)
1554 {
1555 redCost = - normedCoef(sigma,basics_[i]) - (tau3)
1556 + (1 - fzero) * ( - colsol_[basics_[i]]
1557 + upBound);
1558
1559 if (redCost < -tolerance)
1560 {
1561 uu_i[i] = redCost;
1562 hasNegativeRc = true;
1563 }
1564 else
1565 {
1566 if (fabs(redCost) < tolerance) nZeroRc++;
1567 else nPositiveRc ++;
1568 uu_i[i] = 10.;
1569 }
1570
1571 if (redCost < bestReducedCost
1572 && rowFlags_[i]) //row has not been flaged
1573 {
1574 bestDirection = 1;
1575 bestGammaSign = -1;
1576 bestReducedCost = redCost;
1577 bestRow = i;
1578 }
1579
1580 redCost = -normedCoef(sigma,basics_[i]) + (tau4)
1581 - (1 - fzero) * ( - colsol_[basics_[i]]
1582 + upBound)
1583 + upBound - getColsolToCut(basics_[i]);
1584
1585 if (redCost < -tolerance)
1586 {
1587 vu_i[i] = redCost;
1588 hasNegativeRc = true;
1589 }
1590
1591
1592 else
1593 {
1594 if (fabs(redCost) < tolerance) nZeroRc++;
1595 else nPositiveRc ++;
1596 vu_i[i] = 10.;
1597 }
1598
1599 if (redCost < bestReducedCost
1600 && rowFlags_[i]) //row has not been flaged
1601 {
1602 bestDirection = 1;
1603 bestGammaSign = 1;
1604 bestReducedCost = redCost;
1605 bestRow = i;
1606 }
1607 }
1608 else
1609 {
1610 uu_i[i] = 10.;
1611 vu_i[i] = 10.;
1612 }
1613 if (hasNegativeRc) nNegativeRcRows_ ++;
1614 else if (flagPositiveRows) rowFlags_[i] = false;
1615 }
1616 handler_->message(NumberNegRc, messages_)<<nNegativeRcRows_<<CoinMessageEol;
1617 handler_->message(NumberZeroRc, messages_)<<nZeroRc<<CoinMessageEol;
1618 handler_->message(NumberPosRc, messages_)<<nPositiveRc<<CoinMessageEol;
1619 // throw -1;
1620 direction = bestDirection;
1621 gammaSign = bestGammaSign;
1622 // row_i_.num=bestRow;
1623 if (bestRow != -1)
1624 {
1625 chosenReducedCostVal_ = bestReducedCost;
1626 row_i_.num=bestRow;
1627 pullTableauRow(row_i_);
1628 handler_->message(FoundImprovingRow, messages_)<<
1629 bestRow<<basics_[bestRow]<<direction<<gammaSign<<bestReducedCost
1630 <<CoinMessageEol;
1631 }
1632 assert (bestRow<0||direction!=0);
1633 return bestRow;
1634 }
1635
1636
1637 /** Find a row which can be used to perform an improving pivot tables are already filled.
1638 \return index of the cut or -1 if none exists. */
1639 int
rescanReducedCosts(int & direction,int & gammaSign,double tolerance)1640 CglLandPSimplex::rescanReducedCosts( int &direction, int &gammaSign, double tolerance)
1641 {
1642 // The reduced cost are already here in rWk1_ is u^l_j
1643 // rwk2_ is u^u_j rWk3_ is v^u_j rWk4_ is v^u_j
1644 // Let's rename not to get too much confused
1645 double * ul_i = &rWk1_[0];
1646 double * uu_i = &rWk2_[0];
1647 double * vl_i = &rWk3_[0];
1648 double * vu_i = &rWk4_[0];
1649 int bestRow = -1;
1650 int bestDirection = 0;
1651 int bestGammaSign = 0;
1652 // double infty = si_->getInfinity();
1653 double bestReducedCost = -tolerance;
1654 for (int i = 0 ; i < nrows_ ; i++)
1655 {
1656 if (i == row_k_.num//obviously not necessary to combine row k with itself
1657 || !rowFlags_[i] //row has not been flaged
1658 // && fabs(getUpBound(basics_[row_i_.num]) - getLoBound(basics_[row_i_.num]))>1e-09 //variable is not fixed
1659 )
1660 continue;
1661
1662 if (ul_i[i] < bestReducedCost
1663 && rowFlags_[i]) //row has not been flaged
1664 {
1665 bestDirection = -1;
1666 bestGammaSign = -1;
1667 bestReducedCost = ul_i[i];
1668 bestRow = i;
1669 }
1670
1671 if (vl_i[i] < bestReducedCost
1672 && rowFlags_[i]) //row has not been flaged
1673 {
1674 bestDirection = -1;
1675 bestGammaSign = 1;
1676 bestReducedCost = vl_i[i];
1677 bestRow = i;
1678 }
1679
1680
1681 if (uu_i[i] < bestReducedCost
1682 && rowFlags_[i]) //row has not been flaged
1683 {
1684 bestDirection = 1;
1685 bestGammaSign = -1;
1686 bestReducedCost = uu_i[i];
1687 bestRow = i;
1688 }
1689 if (vu_i[i] < bestReducedCost
1690 && rowFlags_[i]) //row has not been flaged
1691 {
1692 bestDirection = 1;
1693 bestGammaSign = 1;
1694 bestReducedCost = vu_i[i];
1695 bestRow = i;
1696 }
1697
1698 }
1699 direction = bestDirection;
1700 gammaSign = bestGammaSign;
1701 // row_i_.num=bestRow;
1702 if (bestRow != -1)
1703 {
1704 chosenReducedCostVal_ = bestReducedCost;
1705 row_i_.num=bestRow;
1706 pullTableauRow(row_i_);
1707 handler_->message(FoundImprovingRow, messages_)<<
1708 bestRow<<basics_[bestRow]<<direction<<gammaSign<<bestReducedCost
1709 <<CoinMessageEol;
1710 }
1711 assert (bestRow<0||direction!=0);
1712 return bestRow;
1713 }
1714
1715
1716 void
compute_p_q_r_s(double gamma,int gammaSign,double & p,double & q,double & r,double & s)1717 CglLandPSimplex::compute_p_q_r_s(double gamma, int gammaSign, double &p, double & q, double & r , double &s)
1718 {
1719 for (int i = 0 ; i < ncols_ ; i++)
1720 {
1721 const int &ii = nonBasics_[i];//True index
1722 if (colCandidateToLeave_[i]==false) continue;
1723 const double& val = getColsolToCut(ii); //value in solution to cut
1724 const double& row_k = row_k_[ii]; // coefficient in row k
1725 const double& row_i = row_i_[ii]; // coefficient in row i
1726 double coeff = row_k + gammaSign * gamma*row_i;
1727 if (coeff>0.)
1728 {
1729 if (gammaSign > 0)
1730 {
1731 p += row_k * val;
1732 }
1733 else
1734 {
1735 // if(fabs(getColsolToCut(nonBasics_[i])) > 0)
1736 {
1737 p += row_k * val;
1738 q += row_i * val;
1739 }
1740 }
1741 r += normedCoef(row_k, ii) ;
1742 s += normedCoef(row_i, ii);
1743 }
1744 else if (coeff< 0.)
1745 {
1746 if (gammaSign > 0)
1747 {
1748 q -= row_i * val;
1749 }
1750 r -= normedCoef(row_k,ii);
1751 s -= normedCoef(row_i,ii);
1752 }
1753 else
1754 {
1755 if (gammaSign > 0 && row_i < 0)
1756 {
1757 q -= row_i * val;
1758 }
1759 else if (gammaSign < 0 && row_i < 0)
1760 {
1761 q += row_i * val;
1762 }
1763 s += normedCoef(gammaSign*fabs(row_i),ii);
1764 }
1765 }
1766 }
1767
1768 // double f_plus(double gamma, int gammaSign){
1769 // double num = row_k_.
1770 // for(int i = 0 ; i < ncols_ ; i++){
1771 // }
1772 //}
1773
1774 /** Find the column which leads to the best cut (i.e., find incoming variable).*/
1775 int
fastFindBestPivotColumn(int direction,int gammaSign,double pivotTol,double rhsTol,bool reducedSpace,bool allowDegenerate,double & bestSigma,bool modularize)1776 CglLandPSimplex::fastFindBestPivotColumn(int direction, int gammaSign,
1777 double pivotTol, double rhsTol,
1778 bool reducedSpace, bool allowDegenerate,
1779 double & bestSigma, bool modularize)
1780 {
1781 gammas_.clear();
1782 pivotTol = 1e-05;
1783 adjustTableauRow(basics_[row_i_.num], row_i_, direction);
1784
1785 double fzero = getColsolToCut(basics_[row_k_.num]) - floor(getColsolToCut(basics_[row_k_.num]));
1786
1787
1788
1789 double p = 0;
1790 double q = 0;
1791 if(!modularize){//Take a shortcut
1792 p = -row_k_.rhs * (1 - fzero);
1793 q = row_i_.rhs * fzero;
1794
1795 if (gammaSign < 0){
1796 q -= row_i_.rhs;
1797 }
1798 }
1799 double r = 1.;
1800 double s = normedCoef( static_cast<double> (gammaSign), basics_[row_i_.num]);
1801
1802 bool haveSmallGammaPivot = false;
1803 double gammaTolerance = 0;
1804 if (allowDegenerate)
1805 gammaTolerance = 0;
1806 //fill the array with the gammas of correct sign
1807 for (int i = 0 ; i < ncols_ ; i++)
1808 {
1809 const int &ii = nonBasics_[i];//True index
1810 const double& val = getColsolToCut(ii); //value in solution to cut
1811 const double& row_k = row_k_[ii]; // coefficient in row k
1812 const double& row_i = row_i_[ii]; // coefficient in row i
1813 if(modularize){
1814 p-=row_k_.rhs*row_k*val;
1815 q-=row_i_.rhs*row_k*val;
1816 }
1817
1818 if (reducedSpace && colCandidateToLeave_[i]==false)
1819 {
1820 assert(col_in_subspace[ii]==false);
1821 continue;
1822 }
1823 double gamma = 1;
1824 if (fabs(row_i) > gammaTolerance && fabs(row_k) > gammaTolerance)
1825 {
1826 gamma = - row_k/row_i;
1827 if (gamma * gammaSign > gammaTolerance)
1828 {
1829 gammas_.insert(i,gamma*gammaSign);
1830 }
1831 }
1832 gamma = fabs(gamma); // we already know the sign of gamma, its absolute value is more usefull
1833 if (row_k>gammaTolerance)
1834 {
1835 if (gammaSign > 0)
1836 {
1837 p += row_k * val;
1838 }
1839 else
1840 {
1841 {
1842 p += row_k * val;
1843 q += row_i * val;
1844 }
1845 }
1846 r += normedCoef(row_k, ii) ;
1847 s += normedCoef(row_i, ii);
1848 }
1849 else if (row_k< gammaTolerance)
1850 {
1851 if (gammaSign > 0)
1852 {
1853 q -= row_i * val;
1854 }
1855 r -= normedCoef(row_k,ii);
1856 s -= normedCoef(row_i,ii);
1857 }
1858 else
1859 {
1860 haveSmallGammaPivot |= true;
1861 if (gammaSign > 0 && row_i < 0)
1862 {
1863 q -= row_i * val;
1864 }
1865 else if (gammaSign < 0 && row_i < 0)
1866 {
1867 q += row_i * val;
1868 }
1869 s += normedCoef(gammaSign*fabs(row_i),ii);
1870 }
1871 }
1872
1873
1874 if(modularize){
1875 p -= row_k_.rhs * (1 - row_k_.rhs);
1876 q += row_i_.rhs * row_k_.rhs;
1877 if (gammaSign < 0){
1878 q -= row_i_.rhs;
1879 }
1880 }
1881
1882 int n = gammas_.getNumElements();
1883 if (n==0)
1884 {
1885 resetOriginalTableauRow(basics_[row_i_.num], row_i_, direction);
1886 return -2;
1887 }
1888 gammas_.sortIncrElement();
1889 const int* inds = gammas_.getIndices();
1890 const double * elements = gammas_.getElements();
1891 int bestColumn = -1;
1892 double newSigma = 1e100;
1893 DblEqAssert(sigma_, rhs_weight_*p/r);
1894 bestSigma = sigma_ = rhs_weight_*p/r;
1895 int lastValid = -1;
1896 #ifndef NDEBUG
1897 bool rc_positive=false;
1898 if (M3_.size())
1899 DblEqAssert( gammaSign*(q * r - p * s)/r, chosenReducedCostVal_);
1900 #endif
1901 if ( gammaSign*(q * r - p * s) >= 0)
1902 {
1903 // after recomputing reduced cost (using exact row) it is found to be >=0
1904 resetOriginalTableauRow(basics_[row_i_.num], row_i_, direction);
1905 return -2;
1906 #ifndef NDEBUG
1907 rc_positive = true;
1908 #endif
1909 }
1910 for (int i = 0 ; i < n ; i++)
1911 {
1912 double newRhs = row_k_.rhs + gammaSign * elements[i] * row_i_.rhs;
1913 if (newRhs < rhsTol || newRhs > 1 - rhsTol)
1914 {
1915 // if(i == 0)
1916 break;
1917 }
1918 newSigma = (p + gammaSign * elements[i] * q)*rhs_weight_/(r + gammaSign*elements[i] * s);
1919 #ifdef DEBUG_LAP
1920 double alt = computeCglpObjective(gammaSign*elements[i], false);
1921 DblEqAssert(newSigma, alt);
1922 #endif
1923 if (newSigma > bestSigma - 1e-08*bestSigma)
1924 {
1925 #ifndef NDEBUG
1926
1927 if (!rc_positive)
1928 {
1929 }
1930
1931
1932 if (0 && elements[i] <= 1e-05)
1933 {
1934 }
1935 #endif
1936 break;
1937 }
1938 else if (newSigma <= bestSigma) // && colCandidateToLeave_[inds[i]])
1939 {
1940 bestColumn = inds[i];
1941 bestSigma = newSigma;
1942 lastValid = i;
1943 }
1944
1945
1946 #ifndef NDEBUG
1947 if (rc_positive)
1948 {
1949 break;
1950 }
1951 #endif
1952 int col = nonBasics_[inds[i]];
1953 if (row_i_[col] *gammaSign > 0)
1954 {
1955 p += row_k_[col] * getColsolToCut(col);
1956 q += row_i_[col] * getColsolToCut(col);
1957 r += normedCoef(row_k_[col]*2,col);
1958 s += normedCoef(row_i_[col]*2,col);
1959 }
1960 else
1961 {
1962 p -= row_k_[col] * getColsolToCut(col);
1963 q -= row_i_[col] * getColsolToCut(col);
1964 r -= normedCoef(row_k_[col]*2,col);
1965 s -= normedCoef(row_i_[col]*2,col);
1966 }
1967 if (gammaSign*(q * r - p * s) >= 0) /* function is starting to increase stop here*/
1968 {
1969 #ifndef NDEBUG
1970 if (0 && elements[i] <= 1e-5)
1971 {
1972 }
1973 rc_positive = true;
1974 #endif
1975 break;
1976 }
1977
1978 }
1979 //Get the results did we find a valid pivot? is it degenerate?
1980 if (bestColumn == -1) // Apparently no pivot is within the tolerances
1981 {
1982 resetOriginalTableauRow(basics_[row_i_.num], row_i_, direction);
1983 handler_->message(WarnFailedPivotTol, messages_)<<CoinMessageEol<<CoinMessageEol;
1984 return -1;
1985 }
1986
1987 if (fabs(row_i_[nonBasics_[bestColumn]]) < 1e-05) // Apparently no pivot is within the tolerances
1988 {
1989 resetOriginalTableauRow(basics_[row_i_.num], row_i_, direction);
1990 handler_->message(WarnFailedPivotTol, messages_)<<CoinMessageEol<<CoinMessageEol;
1991 return -2;
1992 }
1993
1994 //std::cout<<"Minimum of f attained at breakpoint "<<lastValid<<std::endl;
1995 assert(bestSigma <= sigma_);
1996 #ifdef DEBUG_LAP
1997 bestSigma_ = bestSigma;
1998 double otherSigma= computeCglpObjective(gammaSign*elements[lastValid], false);
1999 DblEqAssert(bestSigma, otherSigma);
2000 DblEqAssert(bestSigma, computeCglpObjective(gammaSign*elements[lastValid], false, new_row_));
2001 DblEqAssert(bestSigma, computeCglpObjective(new_row_));
2002 assert(row_k_.modularized_ || row_i_.num != row_k_.num);
2003 #endif
2004
2005 #ifdef OLD_COMPUTATION
2006 if (!modularize)
2007 resetOriginalTableauRow(basics_[row_i_.num], row_i_, direction);
2008 #endif
2009 if (bestSigma < sigma_ - 1e-07) //everything has gone ok
2010 {
2011 handler_->message(FoundBestImprovingCol, messages_)<<nonBasics_[bestColumn]<<gammaSign * elements[lastValid]<<bestSigma<<CoinMessageEol<<CoinMessageEol;
2012 inDegenerateSequence_ = false;
2013 assert (bestColumn<0||direction!=0);
2014 return bestColumn;
2015 }
2016 else if (allowDegenerate) //Pivot is degenerate and we allow
2017 {
2018 inDegenerateSequence_ = true;
2019 assert (bestColumn<0||direction!=0);
2020 return bestColumn;
2021 }
2022 else //we don't accept a degenerate pivot
2023 {
2024 handler_->message(WarnFailedBestImprovingCol, messages_)<<chosenReducedCostVal_<<sigma_<<bestSigma<<CoinMessageEol<<CoinMessageEol;
2025 return -1;
2026 }
2027 }
2028
2029
2030 int
findBestPivotColumn(int direction,double pivotTol,bool reducedSpace,bool allowDegenerate,bool modularize)2031 CglLandPSimplex::findBestPivotColumn(int direction,
2032 double pivotTol, bool reducedSpace, bool allowDegenerate, bool modularize)
2033 {
2034 TabRow newRow(this);
2035 newRow.reserve(ncols_ + nrows_);
2036 int varOut=-1;
2037
2038 adjustTableauRow(basics_[row_i_.num], row_i_, direction);
2039
2040 double m = si_->getInfinity();
2041
2042 int j = 0;
2043 double gamma = 0.;
2044
2045 for (; j< ncols_orig_ ; j++)
2046 {
2047 if (reducedSpace &&
2048 !colCandidateToLeave_[j]
2049 )
2050 continue;
2051 if (fabs(row_i_[nonBasics_[j]])< pivotTol)
2052 {
2053 continue;
2054 }
2055 gamma = - row_k_[nonBasics_[j]]/row_i_[nonBasics_[j]];
2056
2057
2058 newRow[basics_[row_k_.num]] = 1.;
2059 newRow.rhs = row_k_.rhs + gamma * row_i_.rhs;
2060 if (newRow.rhs > 1e-5 && newRow.rhs < 1 - 1e-5 )
2061 {
2062 double m_j = computeCglpObjective(gamma, modularize, newRow);
2063 if (m_j < m)
2064 {
2065 varOut = j;
2066 m = m_j;
2067 }
2068 }
2069 }
2070 resetOriginalTableauRow(basics_[row_i_.num], row_i_, direction);
2071
2072 if (m < sigma_ )
2073 {
2074 handler_->message(FoundBestImprovingCol, messages_)<<nonBasics_[varOut]<<gamma<<m<<CoinMessageEol<<CoinMessageEol;
2075 inDegenerateSequence_ = false;
2076 assert (varOut<0||direction!=0);
2077 return varOut;
2078 }
2079 else if (allowDegenerate && m<=sigma_)
2080 {
2081 inDegenerateSequence_ = true;
2082 }
2083 else
2084 {
2085 return -1;
2086 }
2087 return -1;
2088 }
2089
2090 struct reducedCost
2091 {
2092 /** To avoid computing two times the same row direction will have strange
2093 values, direction is -1 or 1 if for only one of the two direction
2094 rc is <0 and is -2 or 2 if the two direction have one <0 rc with the sign
2095 indicating which one of the two directions is the best.<br>
2096 Note that by theory only one reduced cost (for u_i, or v_i)
2097 maybe negative for each direction.
2098 */
2099 int direction;
2100 /** gammSign is the sign of gamma (corresponding to taking rc for u_i or v_i)
2101 for the best of the two rc for this row.*/
2102 int gammaSign;
2103 /** gammaSign2 is the sign of gamma for the worst of the two rc for this row.*/
2104 int gammaSign2;
2105 /** if both reduced costs are <0 value is the smallest of the two.*/
2106 double value;
2107 /** greatest of the two reduced costs */
2108 double value2;
2109 /** index of the row.*/
2110 int row;
operator <LAP::reducedCost2111 bool operator<(const reducedCost & other) const
2112 {
2113 return (value>other.value);
2114 }
2115 };
2116
2117 /** Find incoming and leaving variables which lead to the most violated
2118 adjacent normalized lift-and-project cut.
2119 \remark At this point reduced costs should be already computed.
2120 \return incoming variable variable,
2121 \param leaving variable
2122 \param direction leaving direction
2123 \param numTryRows number rows tried
2124 \param pivotTol pivot tolerance
2125 \param reducedSpace separaration space (reduced or full)
2126 \param allowNonStrictlyImproving wether or not to allow non stricly improving pivots.
2127 */
findBestPivot(int & leaving,int & direction,const CglLandP::Parameters & params)2128 int CglLandPSimplex::findBestPivot(int &leaving, int & direction,
2129 const CglLandP::Parameters & params)
2130 {
2131 // 1. Sort <0 reduced costs in increasing order
2132 // 2. for numTryRows reduced costs call findBestPivotColumn
2133 // if better record
2134
2135 // The reduced cost are already here in rWk1_ is u^l_j
2136 // rwk2_ is u^u_j rWk3_ is v^u_j rWk4_ is v^u_j
2137 // Let's rename not to get too much confused
2138 double * ul_i = &rWk1_[0];
2139 double * uu_i = &rWk2_[0];
2140 double * vl_i = &rWk3_[0];
2141 double * vu_i = &rWk4_[0];
2142
2143 reducedCost * rc = new reducedCost[nNegativeRcRows_];
2144 int k = 0;
2145 rc[k].direction = 0;//initialize first rc
2146 int k2 = 0;
2147 for (int i = 0 ; i < nrows_ ; i++)
2148 {
2149 if (ul_i[i] < -params.pivotTol)
2150 // && rowFlags_[i]) //row has not been flaged
2151 {
2152 rc[k].direction = -1;
2153 rc[k].gammaSign = -1;
2154 rc[k].value = ul_i[i];
2155 rc[k].row = i;
2156 k2++;
2157 }
2158
2159 if (vl_i[i] < -params.pivotTol)
2160 //&& rowFlags_[i]) //row has not been flaged
2161 {
2162 rc[k].direction = -1;
2163 rc[k].gammaSign = 1;
2164 rc[k].value = vl_i[i];
2165 rc[k].row = i;
2166 k2++;
2167 }
2168
2169
2170 if (uu_i[i] < -params.pivotTol)
2171 // && rowFlags_[i]) //row has not been flaged
2172 {
2173 if (rc[k].direction == 0)
2174 {
2175 rc[k].direction = 1;
2176 rc[k].gammaSign = -1;
2177 rc[k].value = uu_i[i];
2178 rc[k].row = i;
2179 }
2180 else
2181 {
2182 if (uu_i[i] <rc[k].value) //this one is better
2183 {
2184 rc[k].direction = 2;
2185 rc[k].gammaSign2 = rc[k].gammaSign;
2186 rc[k].gammaSign = -1;
2187 rc[k].value2 = rc[k].value;
2188 rc[k].value = uu_i[i];
2189 }
2190 else
2191 {
2192 rc[k].direction = -2;
2193 rc[k].gammaSign2 = -1;
2194 rc[k].value2 = uu_i[i];
2195 }
2196 }
2197 k2++;
2198 }
2199 if (vu_i[i] < -params.pivotTol)
2200 //&& rowFlags_[i]) //row has not been flaged
2201 {
2202 if (rc[k].direction==0)
2203 {
2204 rc[k].direction = 1;
2205 rc[k].gammaSign = 1;
2206 rc[k].value = vu_i[i];
2207 rc[k].row = i;
2208 }
2209 else
2210 {
2211 if (vu_i[i] < rc[k].value) //this one is better
2212 {
2213 rc[k].direction = 2;
2214 rc[k].gammaSign2 = rc[k].gammaSign;
2215 rc[k].gammaSign = 1;
2216 rc[k].value2 = rc[k].value;
2217 rc[k].value = vu_i[i];
2218 }
2219 else //the other one is better
2220 {
2221 rc[k].direction = -2;
2222 rc[k].gammaSign2 = 1;
2223 rc[k].value2 = vu_i[i];
2224 }
2225 }
2226 k2++;
2227 }
2228 if (rc[k].direction!=0) //We have added a row with < 0 rc during
2229 //last iteration
2230 {
2231 k++;
2232 if (k<nNegativeRcRows_)
2233 rc[k].direction = 0;
2234 else
2235 break;
2236 }
2237
2238 }
2239
2240 assert(k==nNegativeRcRows_);
2241
2242 //now make a heap
2243 std::make_heap(rc, rc + k);
2244 // assert(rc[0].value==chosenReducedCostVal_);
2245 int bestLeaving = -1;
2246 int bestIncoming = -1;
2247 int bestDirection = 0;
2248
2249 double bestSigma = COIN_DBL_MAX;
2250 double bestRc = COIN_DBL_MAX;
2251 //now scan the heap
2252 #ifndef NDEBUG
2253 int best_l = 0;
2254 #endif
2255 int notImproved = 0;
2256 for (int l = 0; l < k && l < 10 ; l++, notImproved++)
2257 {
2258 if (!rowFlags_[rc[l].row]) continue;//this row has been marked to be skipped
2259 // if(bestLeaving != -1 && rc[l].value > -1e-02) break;
2260 if (rc[l].value > -1e-02) break;
2261 row_i_.num=rc[l].row;
2262 pullTableauRow(row_i_);//Get the tableau row
2263
2264 //compute f+ or f- for the best negative rc corresponding to this row
2265 chosenReducedCostVal_ = rc[l].value;
2266 double sigma;
2267 int incoming =
2268 fastFindBestPivotColumn
2269 (rc[l].direction, rc[l].gammaSign,
2270 params.pivotTol, params.away,
2271 (params.sepSpace==CglLandP::Fractional),
2272 0,
2273 sigma, params.modularize);
2274 if (incoming!=-1 && bestSigma > sigma)
2275 {
2276 // std::cout<<"I found a better pivot "<<sigma - sigma_<< " for indice number "<<l<<std::endl;
2277 #ifndef NDEBUG
2278 best_l = l;
2279 #endif
2280 bestSigma = sigma;
2281 bestIncoming = incoming;
2282 bestLeaving = rc[l].row;
2283 bestDirection = rc[l].direction > 0 ? 1 : -1;
2284 bestRc = rc[l].value;
2285 notImproved = 0;
2286 }
2287
2288 //Now evenutally compute f+ or f- for the other negative rc (if if exists)
2289 if (rc[l].direction == 2 || rc[l].direction == -2)
2290 {
2291 rc[l].direction/= -2;//Reverse the direction
2292 chosenReducedCostVal_ = rc[l].value2;//need to set this for debug double
2293 //checks
2294 incoming = fastFindBestPivotColumn
2295 (rc[l].direction, rc[l].gammaSign2,
2296 params.pivotTol, params.away,
2297 (params.sepSpace==CglLandP::Fractional),
2298 0,
2299 sigma, params.modularize);
2300 if (incoming!=-1 && bestSigma > sigma)
2301 {
2302 // std::cout<<"I found a better pivot "<<sigma - sigma_<<std::endl;
2303 #ifndef NDEBUG
2304 best_l = l;
2305 #endif
2306 bestSigma = sigma;
2307 bestIncoming = incoming;
2308 bestLeaving = rc[l].row;
2309 bestDirection = rc[l].direction;
2310 bestRc = rc[l].value2;
2311 notImproved = 0;
2312 }
2313 }
2314 }
2315
2316 row_i_.num = leaving = bestLeaving;
2317 chosenReducedCostVal_ = bestRc;
2318 assert(best_l <= nNegativeRcRows_);
2319 if (bestLeaving!=-1)
2320 {
2321 // std::cout<<"Best pivot pivot "<<best_l<<std::endl;
2322 pullTableauRow(row_i_);
2323 }
2324 direction = bestDirection;
2325 delete [] rc;
2326 assert (bestIncoming<0||direction!=0);
2327 return bestIncoming;
2328
2329 }
2330
2331 double
computeCglpObjective(const TabRow & row,bool modularize) const2332 CglLandPSimplex::computeCglpObjective(const TabRow &row, bool modularize) const
2333 {
2334 double numerator = -row.rhs * (1 - row.rhs);
2335 double denominator = 1;
2336
2337 const int & n = row.getNumElements();
2338 const int * ind = row.getIndices();
2339 const double * val = row.denseVector();
2340 for (int j = 0 ; j < n ; j++)
2341 {
2342 const int& jj = ind[j];
2343 if (col_in_subspace[jj]==false) continue;
2344 double coeff = val[jj];
2345 if (modularize && isInteger(jj))
2346 coeff = modularizedCoef(coeff, row.rhs);
2347 denominator += normedCoef(fabs(coeff), jj);
2348 numerator += (coeff > 0 ?
2349 coeff *(1- row.rhs):
2350 - coeff * row.rhs)*getColsolToCut(jj);
2351 }
2352 return numerator*rhs_weight_/denominator;
2353 }
2354
2355 double
computeCglpObjective(double gamma,bool strengthen)2356 CglLandPSimplex::computeCglpObjective(double gamma, bool strengthen)
2357 {
2358 double rhs = row_k_.rhs + gamma * row_i_.rhs;
2359 double numerator = - rhs * (1 - rhs);
2360 double denominator = 1;
2361
2362 double coeff = gamma;//newRowCoefficient(basics_[row_i_.num], gamma);
2363 if (strengthen && isInteger(basics_[row_i_.num]))
2364 coeff = modularizedCoef(coeff, rhs);
2365 denominator += normedCoef(fabs(coeff), basics_[row_i_.num]);
2366 numerator += (coeff > 0 ?
2367 coeff *(1- rhs):
2368 - coeff * rhs)*
2369 getColsolToCut(basics_[row_i_.num]);
2370 for (int j = 0 ; j < ncols_ ; j++)
2371 {
2372 if (col_in_subspace[nonBasics_[j]]==false) {
2373 continue;
2374 }
2375 coeff = newRowCoefficient(nonBasics_[j], gamma);
2376 if (strengthen && nonBasics_[j] < ncols_orig_ && isInteger(j))
2377 coeff = modularizedCoef(coeff, rhs);
2378 denominator += normedCoef(fabs(coeff), nonBasics_[j]);
2379 numerator += (coeff > 0 ?
2380 coeff *(1- rhs):
2381 - coeff * rhs)*
2382 getColsolToCut(nonBasics_[j]);
2383 }
2384 return numerator*rhs_weight_/denominator;
2385 }
2386
2387
2388 double
computeCglpObjective(double gamma,bool strengthen,TabRow & newRow)2389 CglLandPSimplex::computeCglpObjective(double gamma, bool strengthen, TabRow & newRow)
2390 {
2391 newRow.clear();
2392 newRow.rhs = row_k_.rhs + gamma * row_i_.rhs;
2393 double numerator = -newRow.rhs * (1 - newRow.rhs);
2394 double denominator = 1;
2395
2396 int * indices = newRow.getIndices();
2397 int k = 0;
2398 {
2399 if (col_in_subspace[basics_[row_i_.num]]==false)
2400 {
2401 DblEqAssert(0.,1.);
2402 }
2403 double & val = newRow[ basics_[row_i_.num]] = gamma;//newRowCoefficient(basics_[row_i_.num], gamma);
2404 indices[k++] = basics_[row_i_.num];
2405 if (strengthen && row_i_.num < ncols_orig_ && isInteger(row_i_.num))
2406 newRow[ basics_[row_i_.num]] = modularizedCoef(newRow[ basics_[row_i_.num]], newRow.rhs);
2407 denominator += normedCoef(fabs(val), basics_[row_i_.num]);
2408 numerator += (val > 0 ?
2409 val *(1- newRow.rhs):
2410 - val * newRow.rhs)*
2411 getColsolToCut(basics_[row_i_.num]);
2412 }
2413 for (int j = 0 ; j < ncols_ ; j++)
2414 {
2415 double & val = newRow[nonBasics_[j]] = newRowCoefficient(nonBasics_[j], gamma);
2416 indices[k++] = nonBasics_[j];
2417 if (strengthen && nonBasics_[j] < ncols_orig_ && isInteger(j))
2418 newRow[ nonBasics_[j]] = modularizedCoef(val, newRow.rhs);
2419 if (col_in_subspace[nonBasics_[j]]==false) continue;
2420 denominator += normedCoef(fabs(val), nonBasics_[j]);
2421 numerator += (val > 0 ?
2422 val *(1- newRow.rhs):
2423 - val * newRow.rhs)*
2424 getColsolToCut(nonBasics_[j]);
2425 }
2426 newRow.setNumElements(k);
2427 //assert (fabs(numerator/denominator - computeCglpObjective(newRow))<1e-04);
2428 return numerator*rhs_weight_/denominator;
2429 }
2430
2431
2432
2433 /** Compute the reduced cost of Cglp */
2434 double
computeCglpRedCost(int direction,int gammaSign,double tau)2435 CglLandPSimplex::computeCglpRedCost(int direction, int gammaSign, double tau)
2436 {
2437 double toBound;
2438 toBound = direction == -1 ? getLoBound(basics_[row_i_.num]) : getUpBound(basics_[row_i_.num]);
2439
2440 double value =0;
2441 int sign = gammaSign * direction;
2442 double tau1 = 0;
2443 double tau2 = 0;
2444 for (unsigned int i = 0 ; i < M3_.size() ; i++)
2445 {
2446 tau1 += fabs(row_i_ [M3_[i]]);
2447 if (sign == 1 && row_i_[M3_[i]] < 0)
2448 {
2449 tau2 += row_i_[M3_[i]] * getColsolToCut(M3_[i]);
2450 }
2451 else if (sign == -1 && row_i_[M3_[i]] > 0)
2452 {
2453 tau2 += row_i_[M3_[i]] * getColsolToCut(M3_[i]);
2454 }
2455 }
2456 double Tau = - sign * (tau + tau2) - tau1 * sigma_;
2457 value = - sigma_ + Tau
2458 + (1 - getColsolToCut(basics_[row_k_.num])) * sign * (row_i_.rhs
2459 - toBound)
2460 + (gammaSign == 1)*direction*(toBound - getColsolToCut(basics_[row_i_.num]));
2461
2462 return value;
2463 }
2464
2465
2466 /** Compute the value of sigma and thau (which are constants for a row i as defined in Mike Perregaard thesis */
2467 double
computeRedCostConstantsInRow()2468 CglLandPSimplex::computeRedCostConstantsInRow()
2469 {
2470 double tau1 = 0; //the part which will be multiplied by sigma
2471 double tau2 = 0;//the rest
2472
2473 for (unsigned int i = 0 ; i < M1_.size() ; i++)
2474 {
2475 tau1 += row_i_[M1_[i]];
2476 }
2477 for (unsigned int i = 0 ; i < M2_.size() ; i++)
2478 {
2479 tau1 -= row_i_[M2_[i]];
2480 tau2 += row_i_[M2_[i]] * getColsolToCut(M2_[i]);
2481 }
2482 return sigma_ * tau1 + tau2;
2483
2484 }
2485
2486 void
updateM1_M2_M3(TabRow & row,double tolerance,bool perturb)2487 CglLandPSimplex::updateM1_M2_M3(TabRow & row, double tolerance, bool perturb)
2488 {
2489 M1_.clear();
2490 M2_.clear();
2491 M3_.clear();
2492 tolerance = 0;
2493 for (int i = 0; i<ncols_ ; i++)
2494 {
2495 const int &ii = nonBasics_[i];
2496 const double &f = row[ii];
2497 if (f< -tolerance)
2498 {
2499 if (col_in_subspace[ii])
2500 {
2501 M1_.push_back(ii);
2502 colCandidateToLeave_[i]=true;
2503 }
2504 else
2505 {
2506 colCandidateToLeave_[i]=false;
2507 }
2508 }
2509 else if (f>tolerance)
2510 {
2511 if (col_in_subspace[ii])
2512 {
2513 M2_.push_back(ii);
2514 colCandidateToLeave_[i]=true;
2515 }
2516 else
2517 {
2518 colCandidateToLeave_[i]=false;
2519 }
2520 }
2521 else
2522 {
2523 if (col_in_subspace[ii])
2524 {
2525 if (perturb) //assign to M1 or M2 at random
2526 {
2527 int sign = CoinDrand48() > 0.5 ? 1 : -1;
2528 if (sign == -1) //put into M1
2529 {
2530 M1_.push_back(ii);
2531 colCandidateToLeave_[i]=true;
2532 }
2533 else //put into M2
2534 {
2535 M2_.push_back(ii);
2536 colCandidateToLeave_[i]=true;
2537 }
2538 }
2539 else
2540 {
2541 M3_.push_back(ii);
2542 colCandidateToLeave_[i] = true;
2543 }
2544 }
2545 else
2546 {
2547 colCandidateToLeave_[i] = false;
2548 }
2549 }
2550 }
2551 //std::cout<<"M3 has "<<M3_.size()<<" variables."<<std::endl;
2552 }
2553
2554 /** Create the intersection cut of row k*/
2555 void
createIntersectionCut(TabRow & row,OsiRowCut & cut) const2556 CglLandPSimplex::createIntersectionCut(TabRow & row, OsiRowCut &cut) const
2557 {
2558 const double * colLower = si_->getColLower();
2559 const double * rowLower = si_->getRowLower();
2560 const double * colUpper = si_->getColUpper();
2561 const double * rowUpper = si_->getRowUpper();
2562 // double f_0 = row.rhs;
2563 //put the row back into original form
2564 for (int j = 0; j < ncols_ ; j++)
2565 {
2566 if ((nonBasics_[j] < ncols_))
2567 {
2568 CoinWarmStartBasis::Status status = getStatus(nonBasics_[j]);
2569
2570 if (status==CoinWarmStartBasis::atLowerBound)
2571 {
2572 // row.rhs += getLoBound(nonBasics_[j]) * row[nonBasics_[j]];
2573 }
2574 else if (status==CoinWarmStartBasis::atUpperBound)
2575 {
2576 row[nonBasics_[j]] = - row[nonBasics_[j]];
2577 // row.rhs += getUpBound(nonBasics_[j]) * row[nonBasics_[j]];
2578 }
2579 else
2580 {
2581 throw;
2582 }
2583 }
2584 }
2585
2586
2587
2588 // return ;
2589
2590 cut.setUb(COIN_DBL_MAX);
2591 double * vec = new double[ncols_orig_+ nrows_orig_ ];
2592 CoinFillN(vec, ncols_orig_ + nrows_orig_, 0.);
2593 double infty = si_->getInfinity();
2594 double cutRhs = row.rhs;
2595 cutRhs = cutRhs * (1 - cutRhs);
2596 for (int j = 0; j < ncols_ ; j++)
2597 {
2598 if (fabs(row[nonBasics_[j]])>1e-10)
2599 {
2600 double value = intersectionCutCoef(row[nonBasics_[j]], row.rhs);
2601
2602 if (nonBasics_[j]<ncols_)
2603 {
2604 CoinWarmStartBasis::Status status = //CoinWarmStartBasis::basic;
2605 basis_->getStructStatus(nonBasics_[j]);
2606 if (status==CoinWarmStartBasis::atUpperBound)
2607 {
2608 value = - intersectionCutCoef(- row[nonBasics_[j]], row.rhs) ;
2609 cutRhs += value * colUpper[nonBasics_[j]];
2610 }
2611 else
2612 cutRhs += value * colLower[nonBasics_[j]];
2613 vec[original_index_[nonBasics_[j]]] += value;
2614 }
2615 else if (nonBasics_[j]>=ncols_)
2616 {
2617 int iRow = nonBasics_[j] - ncols_;
2618
2619 if (rowLower[iRow] > -infty)
2620 {
2621 value = -value;
2622 cutRhs -= value*rowLower[iRow];
2623 assert(basis_->getArtifStatus(iRow)==CoinWarmStartBasis::atUpperBound ||
2624 (fabs(rowLower[iRow] - rowUpper[iRow]) < 1e-08));
2625 }
2626 else
2627 {
2628 cutRhs -= value*rowUpper[iRow];
2629 assert(basis_->getArtifStatus(iRow)==CoinWarmStartBasis::atLowerBound);
2630 }
2631 vec[nonBasics_[j]] = value;
2632 assert(fabs(cutRhs)<1e100);
2633 }
2634 }
2635 }
2636
2637 const CoinPackedMatrix * mat = si_->getMatrixByCol();
2638 const CoinBigIndex * starts = mat->getVectorStarts();
2639 const int * lengths = mat->getVectorLengths();
2640 const double * values = mat->getElements();
2641 const int * indices = mat->getIndices();
2642 for (int j = 0 ; j < ncols_ ; j++)
2643 {
2644 const CoinBigIndex& start = starts[j];
2645 CoinBigIndex end = start + lengths[j];
2646 for (CoinBigIndex k = start ; k < end ; k++)
2647 {
2648 vec[original_index_[j]] -= vec[original_index_[ncols_ + indices[k]]] * values[k];
2649 }
2650
2651 }
2652
2653 //Pack vec into the cut
2654 int * inds = new int [ncols_orig_];
2655 int nelem = 0;
2656 for (int i = 0 ; i < ncols_orig_ ; i++)
2657 {
2658 if (fabs(vec[i]) > COIN_INDEXED_TINY_ELEMENT)
2659 {
2660 vec[nelem] = vec[i];
2661 inds[nelem++] = i;
2662 }
2663 }
2664
2665 cut.setLb(cutRhs);
2666 cut.setRow(nelem, inds, vec, false);
2667 delete [] vec;
2668
2669 }
2670
2671 /** Compute the normalization factor of the cut.*/
2672 double
normalizationFactor(const TabRow & row) const2673 CglLandPSimplex::normalizationFactor(const TabRow & row) const
2674 {
2675 double numerator = rhs_weight_;
2676 double denominator = 1.;
2677 for (int j = 0 ; j < ncols_ ; j++)
2678 {
2679 denominator += fabs(normedCoef(row[nonBasics_[j]], nonBasics_[j]));
2680 }
2681 return numerator/denominator;
2682 }
2683
2684
2685 /** Create MIG cut from row k*/
2686 void
createMIG(TabRow & row,OsiRowCut & cut) const2687 CglLandPSimplex::createMIG( TabRow &row, OsiRowCut &cut) const
2688 {
2689 const double * colLower = si_->getColLower();
2690 const double * rowLower = si_->getRowLower();
2691 const double * colUpper = si_->getColUpper();
2692 const double * rowUpper = si_->getRowUpper();
2693
2694 if (1)
2695 {
2696 double f_0 = row.rhs - floor(row.rhs);
2697 //put the row back into original form
2698 for (int j = 0; j < ncols_ ; j++)
2699 {
2700 if (nonBasics_[j] < ncols_)
2701 {
2702 const CoinWarmStartBasis::Status status = basis_->getStructStatus(nonBasics_[j]);
2703
2704 if (status==CoinWarmStartBasis::atLowerBound)
2705 {
2706 // row.rhs += getLoBound(nonBasics_[j]) * row[nonBasics_[j]];
2707 }
2708 else if (status==CoinWarmStartBasis::atUpperBound)
2709 {
2710 row[nonBasics_[j]] = - row[nonBasics_[j]];
2711 // row.rhs += getUpBound(nonBasics_[j]) * row[nonBasics_[j]];
2712 }
2713 else
2714 {
2715 throw;
2716 }
2717 }
2718 }
2719 //double scaleFactor = normalizationFactor(row);
2720 row.rhs = f_0;
2721
2722 cut.setUb(COIN_DBL_MAX);
2723 double * vec = new double[ncols_orig_ + nrows_orig_];
2724 CoinFillN(vec, ncols_orig_ + nrows_orig_, 0.);
2725 //f_0 = row.rhs - floor(row.rhs);
2726 double infty = si_->getInfinity();
2727 double cutRhs = row.rhs - floor(row.rhs);
2728 cutRhs = cutRhs * (1 - cutRhs);
2729 assert(fabs(cutRhs)<1e100);
2730 for (int j = 0; j < ncols_ ; j++)
2731 {
2732 if (fabs(row[nonBasics_[j]])>0.)
2733 {
2734 {
2735 if (nonBasics_[j]<ncols_orig_)
2736 {
2737 const CoinWarmStartBasis::Status status = basis_->getStructStatus(nonBasics_[j]);
2738 double value;
2739 if (status==CoinWarmStartBasis::atUpperBound)
2740 {
2741 value = - strengthenedIntersectionCutCoef(nonBasics_[j], - row[nonBasics_[j]], row.rhs) ;
2742 cutRhs += value * colUpper[nonBasics_[j]];
2743 }
2744 else if (status==CoinWarmStartBasis::atLowerBound)
2745 {
2746 value = strengthenedIntersectionCutCoef(nonBasics_[j], row[nonBasics_[j]], row.rhs);
2747 cutRhs += value * colLower[nonBasics_[j]];
2748 }
2749 else
2750 {
2751 std::cerr<<"Invalid basis"<<std::endl;
2752 throw -1;
2753 }
2754
2755 assert(fabs(cutRhs)<1e100);
2756 vec[original_index_[nonBasics_[j]]] = value;
2757 }
2758 else
2759 {
2760 int iRow = nonBasics_[j] - ncols_;
2761 double value = strengthenedIntersectionCutCoef(nonBasics_[j], row[nonBasics_[j]], row.rhs);
2762 if (rowUpper[iRow] < infty)
2763 {
2764 cutRhs -= value*rowUpper[iRow];
2765 }
2766 else
2767 {
2768 value = -value;
2769 cutRhs -= value*rowLower[iRow];
2770 assert(basis_->getArtifStatus(iRow)==CoinWarmStartBasis::atUpperBound ||
2771 (rowUpper[iRow] < infty));
2772 }
2773 vec[original_index_[nonBasics_[j]]] = value;
2774 assert(fabs(cutRhs)<1e100);
2775 }
2776 }
2777 }
2778 }
2779
2780 //Eliminate slacks
2781 eliminate_slacks(vec);
2782
2783 //Pack vec into the cut
2784 int * inds = new int [ncols_orig_];
2785 int nelem = 0;
2786 for (int i = 0 ; i < ncols_orig_ ; i++)
2787 {
2788 if (fabs(vec[i]) > COIN_INDEXED_TINY_ELEMENT)
2789 {
2790 vec[nelem] = vec[i];
2791 inds[nelem++] = i;
2792 }
2793 }
2794
2795 cut.setLb(cutRhs);
2796 cut.setRow(nelem, inds, vec, false);
2797 //std::cout<<"Scale factor: "<<scaleFactor<<" rhs weight "<<rhs_weight_<<std::endl;
2798 //scaleCut(cut, scaleFactor);
2799 delete [] vec;
2800 delete [] inds;
2801
2802 }
2803 }
2804
2805 void
eliminate_slacks(double * vec) const2806 CglLandPSimplex::eliminate_slacks(double * vec) const
2807 {
2808 const CoinPackedMatrix * mat = si_->getMatrixByCol();
2809 const CoinBigIndex * starts = mat->getVectorStarts();
2810 const int * lengths = mat->getVectorLengths();
2811 const double * values = mat->getElements();
2812 const int * indices = mat->getIndices();
2813 const double * vecSlacks = vec + ncols_orig_;
2814 for (int j = 0 ; j < ncols_ ; j++)
2815 {
2816 const CoinBigIndex& start = starts[j];
2817 CoinBigIndex end = start + lengths[j];
2818 double & val = vec[original_index_[j]];
2819 for (CoinBigIndex k = start ; k < end ; k++)
2820 {
2821 val -= vecSlacks[indices[k]] * values[k];
2822 }
2823 }
2824 }
2825 void
scaleCut(OsiRowCut & cut,double factor) const2826 CglLandPSimplex::scaleCut(OsiRowCut & cut, double factor) const
2827 {
2828 DblGtAssert(factor, 0.);
2829 cut *= factor;
2830 cut.setLb(cut.lb()*factor);
2831 }
2832
2833 #ifdef APPEND_ROW
2834 void
append_row(int row_num,bool modularize)2835 CglLandPSimplex::append_row(int row_num, bool modularize)
2836 {
2837 int old_idx = basics_[row_num];
2838 int r_idx = nrows_ - 1;
2839
2840 if (basis_->getArtifStatus(nrows_ - 1) == CoinWarmStartBasis::basic)
2841 {
2842 int i = 0;
2843 for (; i < ncols_ ; i++)
2844 {
2845 if (nonBasics_[i] == ncols_ + nrows_ -1)
2846 break;
2847 }
2848 if (i < ncols_)
2849 {
2850 nonBasics_[i] = ncols_ + nrows_ - 1;
2851 assert(basics_[nrows_ - 1] == ncols_ + nrows_ - 1);
2852 basics_[nrows_ - 1] = ncols_ - 1;
2853 }
2854 }
2855 #if 0
2856 if (basis_->getStructStatus(ncols_ - 1) != CoinWarmStartBasis::basic)
2857 {
2858 basis_->setStructStatus(ncols_ - 1, CoinWarmStartBasis::basic);
2859 assert(basis_->getArtifStatus(nrows_ - 1) == CoinWarmStartBasis::basic);
2860 basis_->setArtifStatus(nrows_ - 1, CoinWarmStartBasis::atLowerBound);
2861 si_->pivot(ncols_ - 1, ncols_ + nrows_ - 1, -1);
2862 }
2863 #endif
2864 TabRow row(this);
2865 int rowsize = ncols_orig_ + nrows_orig_ + 1;
2866 row.reserve(rowsize);
2867 row.num = row_num;
2868 pullTableauRow(row);
2869 row.rhs -= floor(row.rhs);
2870 if (basis_->getStructStatus(ncols_ - 1) != CoinWarmStartBasis::basic)
2871 {
2872 basis_->setStructStatus(ncols_ - 1, CoinWarmStartBasis::basic);
2873 assert(basis_->getArtifStatus(nrows_ - 1) == CoinWarmStartBasis::basic);
2874 basis_->setArtifStatus(nrows_ - 1, CoinWarmStartBasis::atLowerBound);
2875 si_->pivot(ncols_ - 1, ncols_ + nrows_ - 1, -1);
2876 }
2877
2878 int n = row.getNumElements();
2879 const int* ind = row.getIndices();
2880 #ifdef COIN_HAS_OSICLP
2881 CoinPackedMatrix * m = clp_->getMutableMatrixByCol();
2882 #endif
2883
2884 #if 1
2885 const double * rowLower = si_->getRowLower();
2886 const double * rowUpper = si_->getRowUpper();
2887 double infty = si_->getInfinity();
2888 double rhs = floor (colsol_[old_idx]);//floor(row.rhs);
2889 std::vector<double> vec(rowsize,0.);
2890 for (int i = 0 ; i < n ; i++)
2891 {
2892 const int &ni = ind[i];
2893 if (ni == old_idx) continue;
2894 if (integers_[ni])
2895 {
2896 double f = modularizedCoef(row[ni],row.rhs);
2897 f = floor(row[ni] - f + 0.5);
2898 row[ni] -= f;
2899 if (ni < ncols_)
2900 {
2901 if (basis_->getStructStatus(ni) == CoinWarmStartBasis::atUpperBound)
2902 {
2903 f *= -1;
2904 rhs += f * getUpBound(ni);
2905 }
2906 if (basis_->getStructStatus(ni) == CoinWarmStartBasis::atLowerBound)
2907 {
2908 //f *= -1;
2909 rhs -= f * getLoBound(ni);
2910 }
2911 vec[ni]=f;
2912 }
2913 else
2914 {
2915 int iRow = ind[i] - ncols_;
2916 if (basis_->getArtifStatus(iRow)==CoinWarmStartBasis::atLowerBound)
2917 {
2918 rhs -= f*rowUpper[iRow];
2919 }
2920 else if (basis_->getArtifStatus(iRow)==CoinWarmStartBasis::atUpperBound)
2921 {
2922 f *= -1;
2923 rhs -= f*rowLower[iRow];
2924 assert(rowLower[iRow] < infty);
2925 }
2926 vec[ni] = f;
2927 }
2928 }
2929 }
2930 vec[old_idx] = 1;
2931
2932 eliminate_slacks(&vec[0]);
2933
2934 for (int i = 0 ; i < rowsize ; i++)
2935 {
2936 if (vec[i])
2937 m->modifyCoefficient(r_idx, i, vec[i], true);
2938 }
2939 si_->setRowBounds(r_idx, rhs, rhs);
2940 #endif
2941
2942 #if 1
2943 si_->messageHandler()->setLogLevel(0);
2944 si_->resolve();
2945 assert(si_->getIterationCount() == 0);
2946 #endif
2947 /* Update the cached info.*/
2948 #ifndef NDEBUG
2949 int * basics = new int[nrows_];
2950 #else
2951 int * basics = basics_;
2952 #endif
2953 si_->getBasics(basics);
2954 #ifndef NDEBUG
2955 for (int i = 0 ; i < r_idx ; i++)
2956 {
2957 assert(basics[i] == basics_[i]);
2958 }
2959 assert(basics[r_idx] == ncols_ - 1);
2960 delete [] basics_;
2961 basics_ = basics;
2962 #endif
2963 assert(basis_->getStructStatus(old_idx) == CoinWarmStartBasis::basic);
2964 assert(basis_->getStructStatus(ncols_ - 1) == CoinWarmStartBasis::basic);
2965
2966 n = 0;
2967 for (int i = 0 ; i < ncols_ ; i++)
2968 {
2969 if (basis_->getStructStatus(i) != CoinWarmStartBasis::basic)
2970 {
2971 nonBasics_[n++] = i;
2972 }
2973 }
2974 for (int i = 0 ; i < nrows_ ; i++)
2975 {
2976 if (basis_->getArtifStatus(i) != CoinWarmStartBasis::basic)
2977 {
2978 nonBasics_[n++] = i + ncols_;
2979 }
2980 }
2981 assert (n == ncols_);
2982 colsol_[ncols_ - 1] = si_->getColSolution()[ncols_-1];
2983 colsolToCut_[ncols_ - 1] = si_->getColSolution()[ncols_-1];
2984 }
2985
2986 void
check_mod_row(TabRow & row)2987 CglLandPSimplex::check_mod_row(TabRow & row)
2988 {
2989 int rowsize = ncols_orig_ + nrows_orig_;
2990 for (int i = 0 ; i < rowsize ; i++)
2991 {
2992 assert(! integers_[i] || ((row[i] <= row.rhs+1e-08) && (row[i] >= row.rhs - 1- 1e-08)));
2993 }
2994 }
2995
2996 void
update_row(TabRow & row)2997 CglLandPSimplex::update_row(TabRow &row)
2998 {
2999 int r_idx = nrows_ - 1;
3000 int c_idx = ncols_ - 1;
3001 assert(basics_[row.num] == ncols_ - 1);
3002 assert(r_idx == row.num);
3003 int rowsize = nrows_ + ncols_;
3004 int n = row.getNumElements();
3005 const int* ind = row.getIndices();
3006 #ifdef COIN_HAS_OSICLP
3007 CoinPackedMatrix * m = clp_->getMutableMatrixByCol();
3008 #endif
3009 //double * m_el = m->getMutableElements();
3010 //const int * m_id = m->getIndices();
3011 //const int * m_le = m->getVectorLengths();
3012 //const CoinBigIndex * m_st = m->getVectorStarts();
3013 const double * rowLower = si_->getRowLower();
3014 const double * rowUpper = si_->getRowUpper();
3015 double infty = si_->getInfinity();
3016 double rhs = floor(row.rhs);
3017 std::vector<double> vec(rowsize,0.);
3018 for (int i = 0 ; i < n ; i++)
3019 {
3020 const int &ni = ind[i];
3021 if (ni == c_idx) continue;
3022 if (integers_[ni])
3023 {
3024 double f = modularizedCoef(row[ni],row.rhs);
3025 f = floor(row[ni] - f + 0.5);
3026 //row[ni] -= f;
3027 if (ni < ncols_)
3028 {
3029 assert(basis_->getStructStatus(ni) != CoinWarmStartBasis::basic);
3030 if (basis_->getStructStatus(ni) == CoinWarmStartBasis::atUpperBound)
3031 {
3032 f *= -1;
3033 rhs += f * getUpBound(ni);
3034 }
3035 if (basis_->getStructStatus(ni) == CoinWarmStartBasis::atLowerBound)
3036 {
3037 rhs -= f * getLoBound(ni);
3038 }
3039 vec[ni]=f;
3040 }
3041 else
3042 {
3043 //continue;
3044 int iRow = ind[i] - ncols_;
3045 if (iRow == r_idx) continue;
3046 assert(basis_->getArtifStatus(iRow) != CoinWarmStartBasis::basic);
3047 if (basis_->getArtifStatus(iRow)==CoinWarmStartBasis::atLowerBound)
3048 {
3049 rhs -= f*rowUpper[iRow];
3050 }
3051 else if (basis_->getArtifStatus(iRow)==CoinWarmStartBasis::atUpperBound)
3052 {
3053 f *= -1;
3054 rhs -= f*rowLower[iRow];
3055 assert(rowUpper[iRow] < infty);
3056 }
3057 vec[ni] = f;
3058 }
3059 }
3060 }
3061 std::copy(vec.begin(),vec.end(), std::ostream_iterator<double>(std::cout, "\t"));
3062 std::cout<<std::endl;
3063
3064 //vec[c_idx] = 1;
3065 eliminate_slacks(&vec[0]);
3066 std::copy(vec.begin(),vec.end(), std::ostream_iterator<double>(std::cout, "\t"));
3067 std::cout<<std::endl;
3068 // assert(vec[c_idx] == 0);
3069 for (int i = 0 ; i < ncols_ ; i++)
3070 {
3071 if (vec[i])
3072 {
3073 double f = m->getCoefficient(r_idx,i);
3074 m->modifyCoefficient(r_idx, i, f + vec[i], true);
3075 }
3076 }
3077 si_->setRowBounds(r_idx, rowLower[r_idx] + rhs, rowUpper[r_idx] + rhs);
3078
3079 #ifdef COIN_HAS_OSICLP
3080 clp_->getModelPtr()->factorize();
3081 #endif
3082
3083 pullTableauRow(row);
3084 assert(row.rhs >= 0. && row.rhs <= 1.);
3085 bool stop = false;
3086 for (int i = 0 ; i < rowsize ; i++)
3087 {
3088 if (integers_[i] && (row[i] > row.rhs || row[i] < row.rhs - 1))
3089 {
3090 stop = true;
3091 }
3092 //assert(! integers_[i] || (row[i] <= row.rhs && row[i] > row.rhs - 1));
3093 }
3094 if (stop) exit(10);
3095
3096 }
3097 #endif
3098
3099 /** Get the row i of the tableau */
3100 void
pullTableauRow(TabRow & row) const3101 CglLandPSimplex::pullTableauRow(TabRow &row) const
3102 {
3103 const double * rowLower = si_->getRowLower();
3104 const double * rowUpper = si_->getRowUpper();
3105
3106 row.clear();
3107 row.modularized_ = false;
3108 double infty = si_->getInfinity();
3109 /* Get the row */
3110 #ifdef COIN_HAS_OSICLP
3111 if (clp_)
3112 {
3113 CoinIndexedVector array2;
3114 array2.borrowVector(nrows_, 0, row.getIndices() + ncols_, row.denseVector() + ncols_);
3115 clp_->getBInvARow(row.num, &row, &array2);
3116 {
3117 int n = array2.getNumElements();
3118 int * indices1 = row.getIndices() + row.getNumElements();
3119 int * indices2 = array2.getIndices();
3120 for ( int i = 0 ; i < n ; i++)
3121 {
3122 *indices1 = indices2[i] + ncols_;
3123 indices1++;
3124 }
3125 row.setNumElements(n + row.getNumElements());
3126 array2.returnVector();
3127 }
3128 }
3129 else
3130 #endif
3131 {
3132 si_->getBInvARow(row.num,row.denseVector(),row.denseVector() + ncols_);
3133 }
3134 //Clear basic element (it is a trouble for most of the computations)
3135 row[basics_[row.num]]=0.;
3136 // row.row[basics_[row.num]]=1;
3137 /* get the rhs */
3138 {
3139 int iCol = basics_[row.num];
3140 if (iCol<ncols_)
3141 row.rhs = si_->getColSolution()[iCol];
3142 else // Osi does not give direct acces to the value of slacks
3143 {
3144 iCol -= ncols_;
3145 row.rhs = - si_->getRowActivity()[iCol];
3146 if (rowLower[iCol]> -infty)
3147 {
3148 row.rhs += rowLower[iCol];
3149 }
3150 else
3151 {
3152 row.rhs+= rowUpper[iCol];
3153 }
3154 }
3155 }
3156 //Now adjust the row of the tableau to reflect non-basic variables activity
3157 for (int j = 0; j < ncols_ ; j++)
3158 {
3159 if (nonBasics_[j]<ncols_)
3160 {
3161 if (basis_->getStructStatus(nonBasics_[j])==CoinWarmStartBasis::atLowerBound)
3162 {
3163 }
3164 else if (basis_->getStructStatus(nonBasics_[j])==CoinWarmStartBasis::atUpperBound)
3165 {
3166 row[nonBasics_[j]] = -row[nonBasics_[j]];
3167 }
3168 else
3169 {
3170 std::cout<<(basis_->getStructStatus(nonBasics_[j])==CoinWarmStartBasis::isFree)<<std::endl;
3171 throw CoinError("Invalid basis","CglLandPSimplex","pullTableauRow");
3172 }
3173 }
3174 else
3175 {
3176 int iRow = nonBasics_[j] - ncols_;
3177
3178 if (basis_->getArtifStatus(iRow)==CoinWarmStartBasis::atUpperBound)
3179 {
3180 row[nonBasics_[j]] = -row[nonBasics_[j]];
3181 }
3182 }
3183 }
3184 // row.clean(1e-30);
3185 }
3186
3187 /** Adjust the row of the tableau to reflect leaving variable direction */
3188 void
adjustTableauRow(int var,TabRow & row,int direction)3189 CglLandPSimplex::adjustTableauRow(int var, TabRow & row, int direction)
3190 {
3191 double bound = 0;
3192 assert(direction != 0);
3193 if (direction > 0)
3194 {
3195 for (int j = 0 ; j < ncols_orig_ ; j++)
3196 row[nonBasics_[j]] = - row[nonBasics_[j]];
3197
3198 row.rhs = -row.rhs;
3199 bound = getUpBound(var);
3200 setColsolToCut(var, bound - getColsolToCut(var));
3201 row.rhs += bound;
3202 }
3203 else if (direction < 0)
3204 {
3205 bound = getLoBound(var);
3206 setColsolToCut(var, getColsolToCut(var) - bound);
3207 row.rhs -= bound;
3208 }
3209 // assert(fabs(row.rhs)<1e100);
3210 }
3211
3212
3213 /** reset the tableau row after a call to adjustTableauRow */
3214 void
resetOriginalTableauRow(int var,TabRow & row,int direction)3215 CglLandPSimplex::resetOriginalTableauRow(int var, TabRow & row, int direction)
3216 {
3217 if (direction > 0)
3218 {
3219 adjustTableauRow(var, row, direction);
3220 }
3221 else
3222 {
3223 double bound = getLoBound(var);
3224 row.rhs += bound;
3225 setColsolToCut(var, getColsolToCut(var) + bound);
3226 }
3227 }
3228
3229 template <class X>
3230 struct SortingOfArray
3231 {
3232 X * array;
SortingOfArrayLAP::SortingOfArray3233 SortingOfArray(X* a): array(a) {}
3234
operator ()LAP::SortingOfArray3235 bool operator()(const int i, const int j)
3236 {
3237 return array[i] < array[j];
3238 }
3239 };
3240
3241 void
removeRows(int nDelete,const int * rowsIdx)3242 CglLandPSimplex::removeRows(int nDelete, const int * rowsIdx)
3243 {
3244 std::vector<int> sortedIdx;
3245 for (int i = 0 ; i < nDelete ; i++)
3246 sortedIdx.push_back(rowsIdx[i]);
3247 std::sort(sortedIdx.end(), sortedIdx.end());
3248 si_->deleteRows(nDelete, rowsIdx);
3249 int k = 1;
3250 int l = sortedIdx[0];
3251 for (int i = sortedIdx[0] + 1 ; k < nDelete ; i++)
3252 {
3253 if (sortedIdx[k] == i)
3254 {
3255 k++;
3256 }
3257 else
3258 {
3259 original_index_[l] = original_index_[i];
3260 l++;
3261 }
3262 }
3263 delete basis_;
3264 basis_ = dynamic_cast<CoinWarmStartBasis *> (si_->getWarmStart());
3265 assert(basis_);
3266
3267 /* Update rowFlags_ */
3268 std::vector<int> order(nrows_);
3269 for (unsigned int i = 0 ; i < order.size() ; i++)
3270 {
3271 order[i] = i;
3272 }
3273 std::sort(order.begin(), order.end(), SortingOfArray<int>(basics_));
3274 k = 0;
3275 l = 0;
3276 for (int i = 0 ; k < nDelete ; i++)
3277 {
3278 if (basics_[order[i]] == sortedIdx[k])
3279 {
3280 basics_[order[i]] = -1;
3281 k++;
3282 }
3283 else
3284 {
3285 order[l] = order[i];
3286 l++;
3287 }
3288 }
3289 k = 0;
3290 for (int i = 0 ; i < nrows_ ; i++)
3291 {
3292 if (basics_[i] == -1)
3293 {
3294 k++;
3295 }
3296 else
3297 {
3298 basics_[l] = basics_[i];
3299 rowFlags_[l] = rowFlags_[i];
3300 rWk1_[l] = rWk1_[i];
3301 rWk2_[l] = rWk2_[i];
3302 rWk4_[l] = rWk3_[i];
3303 rWk4_[l] = rWk4_[i];
3304 if (row_k_.num == i) row_k_.num = l;
3305
3306 l++;
3307 }
3308 }
3309
3310 nrows_ = nrows_ - nDelete;
3311 original_index_.resize(nrows_);
3312
3313 int numStructural = basis_->getNumStructural();
3314 assert(ncols_ = numStructural);
3315 int nNonBasics = 0;
3316 for (int i = 0 ; i < numStructural ; i++)
3317 {
3318 if (basis_->getStructStatus(i) != CoinWarmStartBasis::basic)
3319 {
3320 nonBasics_[nNonBasics++] = i;
3321 }
3322 }
3323
3324 int numArtificial = basis_->getNumArtificial();
3325 assert(nrows_ = numArtificial);
3326 for (int i = 0 ; i < numArtificial ; i++)
3327 {
3328 if (basis_->getArtifStatus(i)!= CoinWarmStartBasis::basic)
3329 {
3330 nonBasics_[nNonBasics++] = i + numStructural;
3331 }
3332 }
3333 assert (nNonBasics == ncols_);
3334 }
3335
3336 void
printEverything()3337 CglLandPSimplex::printEverything(){
3338 row_k_.print(std::cout, 2, nonBasics_, ncols_);
3339 printf("nonBasics_: ");
3340 for(int i = 0 ; i < ncols_ ; i++){
3341 printf("%5i ",nonBasics_[i]);
3342 }
3343 printf("\n");
3344
3345 printf("basics_: ");
3346 for(int i = 0 ; i < nrows_ ; i++){
3347 printf("%5i ",basics_[i]);
3348 }
3349 printf("\n");
3350
3351 printf("source row:");
3352 for(int i = 0 ; i < ncols_ + nrows_ ; i++){
3353 printf("%10.9g ", row_k_[i]);
3354 }
3355 printf("%10.9g", row_k_.rhs);
3356 printf("\n");
3357
3358 printf(" source indices: ");
3359 for(int i = 0 ; i < row_k_.getNumElements() ; i++){
3360 printf("%5i %20.20g ", row_k_.getIndices()[i], row_k_[row_k_.getIndices()[i]]);
3361 }
3362 printf("\n");
3363
3364 printf("colsolToCut: ");
3365 for(int i = 0 ; i < ncols_ + nrows_ ; i++){
3366 printf("%10.6g ", colsolToCut_[i]);
3367 }
3368 printf("\n");
3369
3370 printf("colsol: ");
3371 for(int i = 0 ; i < ncols_ + nrows_ ; i++){
3372 printf("%10.6g ", colsol_[i]);
3373 }
3374 printf("\n");
3375 }
3376 }/* Ends LAP namespace.*/
3377
3378