1 /*===========================================================================*
2 * This file is part of the Bcps Linear Solver (BLIS). *
3 * *
4 * BLIS is distributed under the Eclipse Public License as part of the *
5 * COIN-OR repository (http://www.coin-or.org). *
6 * *
7 * Authors: *
8 * *
9 * Yan Xu, Lehigh University *
10 * Ted Ralphs, Lehigh University *
11 * *
12 * Conceptual Design: *
13 * *
14 * Yan Xu, Lehigh University *
15 * Ted Ralphs, Lehigh University *
16 * Laszlo Ladanyi, IBM T.J. Watson Research Center *
17 * Matthew Saltzman, Clemson University *
18 * *
19 * *
20 * Copyright (C) 2001-2017, Lehigh University, Yan Xu, and Ted Ralphs. *
21 * All Rights Reserved. *
22 *===========================================================================*/
23
24 #include "CoinHelperFunctions.hpp"
25 #include "CoinPackedVector.hpp"
26 #include "CoinWarmStartBasis.hpp"
27
28 #include "OsiRowCut.hpp"
29
30 #include "AlpsKnowledgeBroker.h"
31
32 #include "BlisHelp.h"
33 #include "BlisConstraint.h"
34 #include "BlisModel.h"
35 #include "BlisSolution.h"
36
37 //#############################################################################
38
39 /** Convert a OsiRowCut to a Blis Contraint. */
BlisOsiCutToConstraint(const OsiRowCut * rowCut)40 BlisConstraint * BlisOsiCutToConstraint(const OsiRowCut *rowCut)
41 {
42 int size = rowCut->row().getNumElements();
43 assert(size > 0);
44
45 const int *ind = rowCut->row().getIndices();
46 const double *val = rowCut->row().getElements();
47
48 double lower = rowCut->lb();
49 double upper = rowCut->ub();
50
51 BlisConstraint *con = new BlisConstraint(lower, upper,
52 lower, upper,
53 size, ind, val);
54
55 if (!con) {
56 // No memory
57 throw CoinError("Out of Memory", "Blis_OsiCutToConstraint", "NONE");
58 }
59
60 return con;
61 }
62
63 //#############################################################################
64
65 /** Convert a Blis constraint to a OsiRowCut. */
BlisConstraintToOsiCut(const BlisConstraint * con)66 OsiRowCut * BlisConstraintToOsiCut(const BlisConstraint * con)
67 {
68 double lower = CoinMax(con->getLbHard(), con->getLbSoft());
69 double upper = CoinMin(con->getUbHard(), con->getUbSoft());
70
71 OsiRowCut * cut = new OsiRowCut;
72 if (!cut) {
73 /* Out of memory. */
74 throw CoinError("Out of Memory", "Blis_constraintToOsiCut", "NONE");
75 }
76
77 assert(con->getSize() > 0);
78
79 cut->setLb(lower);
80 cut->setUb(upper);
81 cut->setRow(con->getSize(),
82 con->getIndices(),
83 con->getValues());
84
85 return cut;
86 }
87
88 //#############################################################################
89
BlisStrongBranch(BlisModel * model,double objValue,int colInd,double x,const double * saveLower,const double * saveUpper,bool & downKeep,bool & downFinished,double & downDeg,bool & upKeep,bool & upFinished,double & upDeg)90 int BlisStrongBranch(BlisModel *model, double objValue, int colInd, double x,
91 const double *saveLower, const double *saveUpper,
92 bool &downKeep, bool &downFinished, double &downDeg,
93 bool &upKeep, bool &upFinished, double &upDeg)
94 {
95 int status = BLIS_OK;
96 int lpStatus = 0;
97
98 int j, numIntInfDown, numObjInfDown;
99
100 double newObjValue;
101
102 OsiSolverInterface * solver = model->solver();
103
104 int numCols = solver->getNumCols();
105 const double * lower = solver->getColLower();
106 const double * upper = solver->getColUpper();
107
108
109 // restore bounds
110 int numDiff = 0;
111
112 #ifdef BLIS_DEBUG_MORE
113 for (j = 0; j < numCols; ++j) {
114 if (saveLower[j] != lower[j]) {
115 //solver->setColLower(j, saveLower[j]);
116 ++numDiff;
117 }
118 if (saveUpper[j] != upper[j]) {
119 //solver->setColUpper(j, saveUpper[j]);
120 ++numDiff;
121 }
122 }
123 std::cout << "BEFORE: numDiff = " << numDiff << std::endl;
124 #endif
125
126 //------------------------------------------------------
127 // Branching down.
128 //------------------------------------------------------
129
130 solver->setColUpper(colInd, floor(x));
131 solver->solveFromHotStart();
132
133 newObjValue = solver->getObjSense() * solver->getObjValue();
134 downDeg = newObjValue - objValue;
135
136 if (solver->isProvenOptimal()) {
137 lpStatus = 0; // optimal
138 #ifdef BLIS_DEBUG_MORE
139 printf("STRONG: COL[%d]: downDeg=%g, x=%g\n", colInd, downDeg, x);
140 #endif
141
142 if (model->feasibleSolution(numIntInfDown, numObjInfDown)) {
143 #ifdef BLIS_DEBUG_MORE
144 printf("STRONG:down:found a feasible solution\n");
145 #endif
146
147 model->setBestSolution(BLIS_SOL_STRONG,
148 newObjValue,
149 solver->getColSolution());
150
151 BlisSolution* ksol = new BlisSolution(solver->getNumCols(),
152 solver->getColSolution(),
153 newObjValue);
154
155 model->broker()->addKnowledge(AlpsKnowledgeTypeSolution,
156 ksol,
157 newObjValue);
158
159 downKeep = false;
160 }
161 else {
162 downKeep = true;
163 }
164 downFinished = true;
165 }
166 else if (solver->isIterationLimitReached() &&
167 !solver->isDualObjectiveLimitReached()) {
168 lpStatus = 2; // unknown
169 downKeep = true;
170 downFinished = false;
171 }
172 else {
173 lpStatus = 1; // infeasible
174 downKeep = false;
175 downFinished = false;
176 }
177
178 #ifdef BLIS_DEBUG_MORE
179 std::cout << "Down: lpStatus = " << lpStatus << std::endl;
180 #endif
181
182 // restore bounds
183 numDiff = 0;
184 for (j = 0; j < numCols; ++j) {
185 if (saveLower[j] != lower[j]) {
186 solver->setColLower(j, saveLower[j]);
187 ++numDiff;
188 }
189 if (saveUpper[j] != upper[j]) {
190 solver->setColUpper(j, saveUpper[j]);
191 ++numDiff;
192 }
193 }
194 #ifdef BLIS_DEBUG
195 assert(numDiff > 0);
196 //std::cout << "numDiff = " << numDiff << std::endl;
197 #endif
198
199 //----------------------------------------------
200 // Branching up.
201 //----------------------------------------------
202
203 solver->setColLower(colInd, ceil(x));
204 solver->solveFromHotStart();
205
206 newObjValue = solver->getObjSense() * solver->getObjValue();
207 upDeg = newObjValue - objValue;
208
209 if (solver->isProvenOptimal()) {
210 lpStatus = 0; // optimal
211
212 #ifdef BLIS_DEBUG_MORE
213 printf("STRONG: COL[%d]: upDeg=%g, x=%g\n", colInd, upDeg, x);
214 #endif
215
216 if (model->feasibleSolution(numIntInfDown, numObjInfDown)) {
217 #ifdef BLIS_DEBUG_MORE
218 printf("STRONG: up:found a feasible solution\n");
219 #endif
220
221 model->setBestSolution(BLIS_SOL_STRONG,
222 newObjValue,
223 solver->getColSolution());
224
225 BlisSolution* ksol = new BlisSolution(solver->getNumCols(),
226 solver->getColSolution(),
227 newObjValue);
228
229 model->broker()->addKnowledge(AlpsKnowledgeTypeSolution,
230 ksol,
231 newObjValue);
232 // FIXME: should not keep this branch.
233 upKeep = false;
234 }
235 else {
236 upKeep = true;
237 }
238 upFinished = true;
239 }
240 else if (solver->isIterationLimitReached()
241 &&!solver->isDualObjectiveLimitReached()) {
242 lpStatus = 2; // unknown
243 upKeep = true;
244 upFinished = false;
245 }
246 else {
247 lpStatus = 1; // infeasible
248 upKeep = false;
249 upFinished = false;
250 }
251
252 #ifdef BLIS_DEBUG_MORE
253 std::cout << "STRONG: Up: lpStatus = " << lpStatus << std::endl;
254 #endif
255
256 // restore bounds
257 for (j = 0; j < numCols; ++j) {
258 if (saveLower[j] != lower[j]) {
259 solver->setColLower(j,saveLower[j]);
260 }
261 if (saveUpper[j] != upper[j]) {
262 solver->setColUpper(j,saveUpper[j]);
263 }
264 }
265
266 return status;
267 }
268
269 //#############################################################################
270
BlisEncodeWarmStart(AlpsEncoded * encoded,const CoinWarmStartBasis * ws)271 int BlisEncodeWarmStart(AlpsEncoded *encoded, const CoinWarmStartBasis *ws)
272 {
273
274 int status = BLIS_OK;
275 int numCols = ws->getNumStructural();
276 int numRows = ws->getNumArtificial();
277
278 encoded->writeRep(numCols);
279 encoded->writeRep(numRows);
280
281 // Pack structural.
282 int nint = (ws->getNumStructural() + 15) >> 4;
283 encoded->writeRep(ws->getStructuralStatus(), nint * 4);
284
285 // Pack artificial.
286 nint = (ws->getNumArtificial() + 15) >> 4;
287 encoded->writeRep(ws->getArtificialStatus(), nint * 4);
288
289 return status;
290 }
291
292 //#############################################################################
293
BlisDecodeWarmStart(AlpsEncoded & encoded,AlpsReturnStatus * rc)294 CoinWarmStartBasis *BlisDecodeWarmStart(AlpsEncoded &encoded,
295 AlpsReturnStatus *rc)
296 {
297 int numCols;
298 int numRows;
299
300 encoded.readRep(numCols);
301 encoded.readRep(numRows);
302
303 int tempInt;
304
305 // Structural
306 int nint = (numCols + 15) >> 4;
307 char *structuralStatus = new char[4 * nint];
308 encoded.readRep(structuralStatus, tempInt);
309 assert(tempInt == nint*4);
310
311 // Artificial
312 nint = (numRows + 15) >> 4;
313 char *artificialStatus = new char[4 * nint];
314 encoded.readRep(artificialStatus, tempInt);
315 assert(tempInt == nint*4);
316
317 CoinWarmStartBasis *ws = new CoinWarmStartBasis();
318 if (!ws) {
319 throw CoinError("Out of memory", "BlisDecodeWarmStart", "HELP");
320 }
321
322 ws->assignBasisStatus(numCols, numRows,
323 structuralStatus, artificialStatus);
324
325 return ws;
326
327 }
328
329 //#############################################################################
330
331 /** Compute and return a hash value of an Osi row cut. */
BlisHashingOsiRowCut(const OsiRowCut * rowCut,const BlisModel * model)332 double BlisHashingOsiRowCut(const OsiRowCut *rowCut,
333 const BlisModel *model)
334 {
335 int size = rowCut->row().getNumElements();
336 assert(size > 0);
337
338 int ind, k;
339
340 const int *indices = rowCut->row().getIndices();
341 const double * randoms = model->getConRandoms();
342
343 double hashValue_ = 0.0;
344 for (k = 0; k < size; ++k) {
345 ind = indices[k];
346 hashValue_ += randoms[ind] * ind;
347 }
348 #ifdef BLIS_DEBUG_MORE
349 std::cout << "hashValue_=" << hashValue_ << std::endl;
350 #endif
351 return hashValue_;
352 }
353
354 //#############################################################################
355
356 /** Check if a row cut parallel with another row cut. */
BlisParallelCutCut(OsiRowCut * rowCut1,OsiRowCut * rowCut2,double threshold)357 bool BlisParallelCutCut(OsiRowCut * rowCut1,
358 OsiRowCut * rowCut2,
359 double threshold)
360 {
361 int size1 = rowCut1->row().getNumElements();
362 int size2 = rowCut2->row().getNumElements();
363 assert(size1 > 0 && size2 > 0);
364
365 int i, j, k;
366
367 //------------------------------------------------------
368 // Assume no duplicate indices.
369 // rowCut->setRow() tests duplicate indices.
370 //------------------------------------------------------
371
372 rowCut1->sortIncrIndex();
373 rowCut2->sortIncrIndex();
374
375 const int *indices1 = rowCut1->row().getIndices();
376 const double *elems1 = rowCut1->row().getElements();
377
378 const int *indices2 = rowCut2->row().getIndices();
379 const double *elems2 = rowCut2->row().getElements();
380
381 //------------------------------------------------------
382 // Compute norms.
383 //------------------------------------------------------
384
385 double norm1 = 0.0;
386 double norm2 = 0.0;
387
388 for (k = 0; k < size1; ++k) {
389 norm1 += elems1[k]*elems1[k];
390 }
391 for (k = 0; k < size2; ++k) {
392 norm2 += elems2[k]*elems2[k];
393 }
394 norm1 = sqrt(norm1);
395 norm2 = sqrt(norm2);
396
397 //------------------------------------------------------
398 // Compute angel cut1 * cut2
399 //------------------------------------------------------
400
401 double denorm = 0.0;
402 double angle;
403 int index1, index2;
404 i = 0; j = 0;
405 while(true) {
406 index1 = indices1[i];
407 index2 = indices2[j];
408 if (index1 == index2) {
409 denorm += (elems1[i] * elems2[j]);
410 ++i;
411 ++j;
412 if ((i >= size1) || (j >= size2)) {
413 break;
414 }
415 }
416 else if (index1 > index2) {
417 ++j;
418 if (j >= size2) break;
419 }
420 else {
421 ++i;
422 if (i >= size1) break;
423 }
424 }
425 denorm = fabs(denorm);
426 angle = denorm/(norm1 * norm2);
427 assert(angle >= 0.0 && angle <= 1.000001);
428
429 if (angle >= threshold) {
430 return true;
431 }
432 else {
433 return false;
434 }
435 }
436
437 //#############################################################################
438
439 /** Check if a row cut parallel with a constraint. */
BlisParallelCutCon(OsiRowCut * rowCut,BlisConstraint * con,double threshold)440 bool BlisParallelCutCon(OsiRowCut * rowCut,
441 BlisConstraint * con,
442 double threshold)
443 {
444 bool parallel;
445 // Convert con to row cut
446 OsiRowCut * rowCut2 = BlisConstraintToOsiCut(con);
447
448 parallel = BlisParallelCutCut(rowCut,
449 rowCut2,
450 threshold);
451
452 return parallel;
453 }
454
455 //#############################################################################
456