1 #include <cstdlib>
2 #include <cmath>
3 #include <ctime>
4 #include <cassert>
5 #include <cstdio>
6 #include <cstring>
7 #include <algorithm>
8 #include <vector>
9 #include <string>
10 #include <map>
11 #include <OsiSolverInterface.hpp>
12 #include "CbcMessage.hpp"
13 #include "CbcHeuristic.hpp"
14 #include <CbcModel.hpp>
15 #include "CbcMipStartIO.hpp"
16 #include "CbcSOS.hpp"
17 #include "CoinTime.hpp"
18
19 using namespace std;
20
isNumericStr(const char * str)21 bool isNumericStr(const char *str)
22 {
23 const size_t l = strlen(str);
24
25 for (size_t i = 0; i < l; ++i)
26 if (!(isdigit(str[i]) || (str[i] == '.') || (str[i] == '-') || (str[i] == '+') || (str[i] == 'e')))
27 return false;
28
29 return true;
30 }
31
readMIPStart(CbcModel * model,const char * fileName,vector<pair<string,double>> & colValues,double &)32 int readMIPStart(CbcModel *model, const char *fileName,
33 vector< pair< string, double > > &colValues,
34 double & /*solObj*/)
35 {
36 #define STR_SIZE 256
37 FILE *f = fopen(fileName, "r");
38 if (!f)
39 return 1;
40 char line[STR_SIZE];
41
42 int nLine = 0;
43 char printLine[STR_SIZE];
44 while (fgets(line, STR_SIZE, f)) {
45 ++nLine;
46 char col[4][STR_SIZE];
47 int nread = sscanf(line, "%s %s %s %s", col[0], col[1], col[2], col[3]);
48 if (!nread)
49 continue;
50 /* line with variable value */
51 if (strlen(col[0]) && isdigit(col[0][0]) && (nread >= 3)) {
52 if (!isNumericStr(col[0])) {
53 sprintf(printLine, "Reading: %s, line %d - first column in mipstart file should be numeric, ignoring.", fileName, nLine);
54 model->messageHandler()->message(CBC_GENERAL, model->messages()) << printLine << CoinMessageEol;
55 continue;
56 }
57 if (!isNumericStr(col[2])) {
58 sprintf(printLine, "Reading: %s, line %d - Third column in mipstart file should be numeric, ignoring.", fileName, nLine);
59 model->messageHandler()->message(CBC_GENERAL, model->messages()) << printLine << CoinMessageEol;
60 continue;
61 }
62
63 char *name = col[1];
64 double value = atof(col[2]);
65
66 colValues.push_back(pair< string, double >(string(name), value));
67 }
68 }
69
70 if (colValues.size()) {
71 sprintf(printLine, "MIPStart values read for %d variables.", static_cast< int >(colValues.size()));
72 model->messageHandler()->message(CBC_GENERAL, model->messages()) << printLine << CoinMessageEol;
73 if (colValues.size() < model->getNumCols()) {
74 int numberColumns = model->getNumCols();
75 OsiSolverInterface *solver = model->solver();
76 vector< pair< string, double > > fullValues;
77 /* for fast search of column names */
78 map< string, int > colIdx;
79 for (int i = 0; i < numberColumns; i++) {
80 fullValues.push_back(pair< string, double >(solver->getColName(i), 0.0));
81 colIdx[solver->getColName(i)] = i;
82 }
83 for (int i = 0; (i < static_cast< int >(colValues.size())); ++i) {
84 map< string, int >::const_iterator mIt = colIdx.find(colValues[i].first);
85 if (mIt != colIdx.end()) {
86 const int idx = mIt->second;
87 double v = colValues[i].second;
88 fullValues[idx].second = v;
89 }
90 }
91 colValues = fullValues;
92 }
93 } else {
94 sprintf(printLine, "No mipstart solution read from %s", fileName);
95 model->messageHandler()->message(CBC_GENERAL, model->messages()) << printLine << CoinMessageEol;
96 fclose(f);
97 return 1;
98 }
99
100 fclose(f);
101 return 0;
102 }
103
computeCompleteSolution(CbcModel * model,const vector<string> colNames,const std::vector<std::pair<std::string,double>> & colValues,double * sol,double & obj)104 int computeCompleteSolution(CbcModel *model,
105 const vector< string > colNames,
106 const std::vector< std::pair< std::string, double > > &colValues,
107 double *sol, double &obj)
108 {
109 if (!model->getNumCols())
110 return 0;
111
112 int status = 0;
113 double compObj = COIN_DBL_MAX;
114 bool foundIntegerSol = false;
115 OsiSolverInterface *lp = model->solver()->clone();
116 map< string, int > colIdx;
117 assert((static_cast< int >(colNames.size())) == lp->getNumCols());
118 /* for fast search of column names */
119 for (int i = 0; (i < static_cast< int >(colNames.size())); ++i)
120 colIdx[colNames[i]] = i;
121
122 char printLine[STR_SIZE];
123 int fixed = 0;
124 int notFound = 0;
125 char colNotFound[256] = "";
126 int nContinuousFixed = 0;
127 double *realObj = new double[lp->getNumCols()];
128 memcpy(realObj, lp->getObjCoefficients(), sizeof(double)*lp->getNumCols());
129
130 #ifndef JUST_FIX_INTEGER
131 #define JUST_FIX_INTEGER 0
132 #endif
133
134 #if JUST_FIX_INTEGER > 1
135 // all not mentioned are at zero
136 for (int i = 0; (i < lp->getNumCols()); ++i) {
137 if (lp->isInteger(i))
138 lp->setColBounds(i, 0.0, 0.0);
139 }
140 #endif
141 for (int i = 0; (i < static_cast< int >(colValues.size())); ++i) {
142 map< string, int >::const_iterator mIt = colIdx.find(colValues[i].first);
143 if (mIt == colIdx.end()) {
144 if (!notFound)
145 strcpy(colNotFound, colValues[i].first.c_str());
146 notFound++;
147 } else {
148 const int idx = mIt->second;
149 double v = colValues[i].second;
150 #if JUST_FIX_INTEGER
151 if (!lp->isInteger(idx))
152 continue;
153 #endif
154 if (fabs(v) < 1e-8)
155 v = 0.0;
156 if (lp->isInteger(idx)) // just to avoid small
157 v = floor(v + 0.5); // fractional garbage
158 else
159 nContinuousFixed++;
160
161 lp->setColBounds(idx, v, v);
162 ++fixed;
163 }
164 }
165
166 if (!fixed) {
167 model->messageHandler()->message(CBC_GENERAL, model->messages())
168 << "Warning: MIPstart solution is not valid, column names do not match, ignoring it."
169 << CoinMessageEol;
170 goto TERMINATE;
171 }
172
173 if (notFound >= ((static_cast< double >(colNames.size())) * 0.5)) {
174 sprintf(printLine, "Warning: %d column names were not found (e.g. %s) while filling solution.", notFound, colNotFound);
175 model->messageHandler()->message(CBC_GENERAL, model->messages())
176 << printLine << CoinMessageEol;
177 }
178 #if JUST_FIX_INTEGER
179 lp->setHintParam(OsiDoPresolveInInitial, true, OsiHintDo);
180 #endif
181 lp->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
182 lp->initialSolve();
183
184 if ((lp->isProvenPrimalInfeasible()) || (lp->isProvenDualInfeasible())) {
185 if (nContinuousFixed) {
186 model->messageHandler()->message(CBC_GENERAL, model->messages())
187 << "Trying just fixing integer variables (and fixingish SOS)." << CoinMessageEol;
188 int numberColumns = lp->getNumCols();
189 const double *oldLower = model->solver()->getColLower();
190 const double *oldUpper = model->solver()->getColUpper();
191 double *savedSol = CoinCopyOfArray(lp->getColLower(), numberColumns);
192 for (int i = 0; i < numberColumns; ++i) {
193 if (!lp->isInteger(i)) {
194 lp->setColLower(i, oldLower[i]);
195 lp->setColUpper(i, oldUpper[i]);
196 }
197 }
198 // but look at SOS
199 int numberObjects = model->numberObjects();
200 for (int i = 0; i < numberObjects; i++) {
201 const CbcSOS *object = dynamic_cast< const CbcSOS * >(model->object(i));
202 if (object) {
203 int n = object->numberMembers();
204 const int *members = object->members();
205 int sosType = object->sosType();
206 if (sosType == 1) {
207 // non zero can take any value - others zero
208 int iColumn = -1;
209 for (int j = 0; j < n; j++) {
210 int jColumn = members[j];
211 if (savedSol[jColumn])
212 iColumn = jColumn;
213 }
214 for (int j = 0; j < n; j++) {
215 int jColumn = members[j];
216 if (jColumn != iColumn) {
217 lp->setColLower(jColumn, 0.0);
218 lp->setColUpper(jColumn, 0.0);
219 }
220 }
221 } else if (sosType == 2) {
222 // SOS 2 - make a guess if just one nonzero
223 int jA = -1;
224 int jB = -1;
225 for (int j = 0; j < n; j++) {
226 int jColumn = members[j];
227 if (savedSol[jColumn]) {
228 if (jA == -1)
229 jA = j;
230 jB = j;
231 }
232 }
233 if (jB > jA + 1) {
234 jB = jA + 1;
235 } else if (jA == jB) {
236 if (jA == n - 1)
237 jA--;
238 else
239 jB++;
240 }
241 for (int j = 0; j < n; j++) {
242 if (j != jA && j != jB) {
243 int jColumn = members[j];
244 lp->setColLower(jColumn, 0.0);
245 lp->setColUpper(jColumn, 0.0);
246 }
247 }
248 }
249 }
250 }
251 delete[] savedSol;
252 lp->initialSolve();
253 } else {
254 model->messageHandler()->message(CBC_GENERAL, model->messages())
255 << "Fixing only non-zero variables." << CoinMessageEol;
256 /* unfix all variables which are zero */
257 int notZeroAnymore = 0;
258 for (int i = 0; (i < lp->getNumCols()); ++i)
259 if (((fabs(lp->getColLower()[i])) <= 1e-8) && (fabs(lp->getColLower()[i] - lp->getColUpper()[i]) <= 1e-8)) {
260 const double *oldLower = model->solver()->getColLower();
261 const double *oldUpper = model->solver()->getColUpper();
262 lp->setColLower(i, oldLower[i]);
263 lp->setColUpper(i, oldUpper[i]);
264 notZeroAnymore++;
265 }
266 if (notZeroAnymore)
267 lp->initialSolve();
268 }
269 }
270
271 if (!lp->isProvenOptimal()) {
272 model->messageHandler()->message(CBC_GENERAL, model->messages())
273 << "Warning: mipstart values could not be used to build a solution." << CoinMessageEol;
274 status = 1;
275 goto TERMINATE;
276 }
277
278 /* some additional effort is needed to provide an integer solution */
279 if (lp->getFractionalIndices().size() > 0) {
280 sprintf(printLine, "MIPStart solution provided values for %d of %d integer variables, %d variables are still fractional.", fixed, lp->getNumIntegers(), static_cast< int >(lp->getFractionalIndices().size()));
281 model->messageHandler()->message(CBC_GENERAL, model->messages())
282 << printLine << CoinMessageEol;
283 double start = CoinCpuTime();
284 #if 1
285 CbcSerendipity heuristic(*model);
286 heuristic.setFractionSmall(2.0);
287 heuristic.setFeasibilityPumpOptions(1008013);
288 int returnCode = heuristic.smallBranchAndBound(lp,
289 1000, sol,
290 compObj,
291 model->getCutoff(),
292 "ReduceInMIPStart");
293 if ((returnCode & 1) != 0) {
294 sprintf(printLine, "Mini branch and bound defined values for remaining variables in %.2f seconds.",
295 CoinCpuTime() - start);
296 model->messageHandler()->message(CBC_GENERAL, model->messages())
297 << printLine << CoinMessageEol;
298 foundIntegerSol = true;
299 obj = compObj;
300 }
301 #else
302 CbcModel babModel(*lp);
303 lp->writeLp("lessFix");
304 babModel.setLogLevel(2);
305 babModel.setMaximumNodes(1000);
306 babModel.setMaximumSeconds(60);
307 babModel.branchAndBound();
308 if (babModel.bestSolution()) {
309 sprintf(printLine, "Mini branch and bound defined values for remaining variables in %.2f seconds.",
310 CoinCpuTime() - start);
311 model->messageHandler()->message(CBC_GENERAL, model->messages())
312 << printLine << CoinMessageEol;
313 copy(babModel.bestSolution(), babModel.bestSolution() + babModel.getNumCols(), sol);
314 foundIntegerSol = true;
315 obj = compObj = babModel.getObjValue();
316 }
317 #endif
318 else {
319 model->messageHandler()->message(CBC_GENERAL, model->messages())
320 << "Warning: mipstart values could not be used to build a solution." << CoinMessageEol;
321 status = 1;
322 goto TERMINATE;
323 }
324 } else {
325 foundIntegerSol = true;
326
327 obj = 0.0;
328 for ( int i=0 ; (i<lp->getNumCols()) ; ++i )
329 obj += realObj[i]*lp->getColSolution()[i];
330 compObj = obj;
331 copy(lp->getColSolution(), lp->getColSolution() + lp->getNumCols(), sol);
332 }
333
334 if (foundIntegerSol) {
335 sprintf(printLine, "MIPStart provided solution with cost %g", compObj);
336 model->messageHandler()->message(CBC_GENERAL, model->messages())
337 << printLine << CoinMessageEol;
338 #if 0
339 {
340 int numberColumns=lp->getNumCols();
341 double largestInfeasibility = 0.0;
342 double primalTolerance ;
343 double offset;
344 lp->getDblParam(OsiObjOffset, offset);
345 lp->getDblParam(OsiPrimalTolerance, primalTolerance) ;
346 const double *objective = lp->getObjCoefficients() ;
347 const double * rowLower = lp->getRowLower() ;
348 const double * rowUpper = lp->getRowUpper() ;
349 const double * columnLower = lp->getColLower() ;
350 const double * columnUpper = lp->getColUpper() ;
351 int numberRows = lp->getNumRows() ;
352 double *rowActivity = new double[numberRows] ;
353 memset(rowActivity, 0, numberRows*sizeof(double)) ;
354 double *rowSum = new double[numberRows] ;
355 memset(rowSum, 0, numberRows*sizeof(double)) ;
356 const double * element = lp->getMatrixByCol()->getElements();
357 const int * row = lp->getMatrixByCol()->getIndices();
358 const CoinBigIndex * columnStart = lp->getMatrixByCol()->getVectorStarts();
359 const int * columnLength = lp->getMatrixByCol()->getVectorLengths();
360 const CoinPackedMatrix * rowCopy = lp->getMatrixByRow();
361 const int * column = rowCopy->getIndices();
362 const int * rowLength = rowCopy->getVectorLengths();
363 const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
364 const double * elementByRow = rowCopy->getElements();
365 double objValue=-offset;
366 for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
367 double value = sol[iColumn];
368 if (lp->isInteger(iColumn))
369 assert (fabs(value-floor(value+0.5))<1.0e-6);
370 objValue += value*objective[iColumn];
371 if (value>columnUpper[iColumn]) {
372 if (value-columnUpper[iColumn]>1.0e-8)
373 printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]);
374 value=columnUpper[iColumn];
375 } else if (value<columnLower[iColumn]) {
376 if (value-columnLower[iColumn]<-1.0e-8)
377 printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]);
378 value=columnLower[iColumn];
379 }
380 if (value) {
381 CoinBigIndex start = columnStart[iColumn];
382 CoinBigIndex end = start + columnLength[iColumn];
383 for (CoinBigIndex j = start; j < end; j++) {
384 int iRow = row[j];
385 if (fabs(value)<1.0e-6&&fabs(value*element[j])>1.0e-5)
386 printf("Column %d row %d value %.8g element %g %s\n",
387 iColumn,iRow,value,element[j],lp->isInteger(iColumn) ? "integer" : "");
388 rowActivity[iRow] += value * element[j];
389 rowSum[iRow] += fabs(value * element[j]);
390 }
391 }
392 }
393 for (int i = 0 ; i < numberRows ; i++) {
394 #if 0 //def CLP_INVESTIGATE
395 double inf;
396 inf = rowLower[i] - rowActivity[i];
397 if (inf > primalTolerance)
398 printf("Row %d inf %g sum %g %g <= %g <= %g\n",
399 i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
400 inf = rowActivity[i] - rowUpper[i];
401 if (inf > primalTolerance)
402 printf("Row %d inf %g sum %g %g <= %g <= %g\n",
403 i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
404 #endif
405 double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
406 rowLower[i]-rowActivity[i]);
407 // but allow for errors
408 double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
409 if (infeasibility>largestInfeasibility*factor) {
410 largestInfeasibility = infeasibility/factor;
411 printf("Cinf of %g on row %d sum %g scaled %g\n",
412 infeasibility,i,rowSum[i],largestInfeasibility);
413 if (infeasibility>1.0e10) {
414 for (CoinBigIndex j=rowStart[i];
415 j<rowStart[i]+rowLength[i];j++) {
416 printf("col %d element %g\n",
417 column[j],elementByRow[j]);
418 }
419 }
420 }
421 }
422 delete [] rowActivity ;
423 delete [] rowSum;
424 if (largestInfeasibility > 10.0*primalTolerance)
425 printf("Clargest infeasibility is %g - obj %g\n", largestInfeasibility,objValue);
426 else
427 printf("Cfeasible (%g) - obj %g\n", largestInfeasibility,objValue);
428 }
429 #endif
430 for (int i = 0; (i < lp->getNumCols()); ++i) {
431 #if 0
432 if (sol[i]<1e-8)
433 sol[i] = 0.0;
434 else
435 if (lp->isInteger(i))
436 sol[i] = floor( sol[i]+0.5 );
437 #else
438 if (lp->isInteger(i)) {
439 //if (fabs(sol[i] - floor( sol[i]+0.5 ))>1.0e-8)
440 //printf("bad sol for %d - %.12g\n",i,sol[i]);
441 sol[i] = floor(sol[i] + 0.5);
442 }
443 #endif
444 }
445 #if 0
446 {
447 int numberColumns=lp->getNumCols();
448 double largestInfeasibility = 0.0;
449 double primalTolerance ;
450 double offset;
451 lp->getDblParam(OsiObjOffset, offset);
452 lp->getDblParam(OsiPrimalTolerance, primalTolerance) ;
453 const double *objective = lp->getObjCoefficients() ;
454 const double * rowLower = lp->getRowLower() ;
455 const double * rowUpper = lp->getRowUpper() ;
456 const double * columnLower = lp->getColLower() ;
457 const double * columnUpper = lp->getColUpper() ;
458 int numberRows = lp->getNumRows() ;
459 double *rowActivity = new double[numberRows] ;
460 memset(rowActivity, 0, numberRows*sizeof(double)) ;
461 double *rowSum = new double[numberRows] ;
462 memset(rowSum, 0, numberRows*sizeof(double)) ;
463 const double * element = lp->getMatrixByCol()->getElements();
464 const int * row = lp->getMatrixByCol()->getIndices();
465 const CoinBigIndex * columnStart = lp->getMatrixByCol()->getVectorStarts();
466 const int * columnLength = lp->getMatrixByCol()->getVectorLengths();
467 const CoinPackedMatrix * rowCopy = lp->getMatrixByRow();
468 const int * column = rowCopy->getIndices();
469 const int * rowLength = rowCopy->getVectorLengths();
470 const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
471 const double * elementByRow = rowCopy->getElements();
472 double objValue=-offset;
473 for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
474 double value = sol[iColumn];
475 if (lp->isInteger(iColumn))
476 assert (fabs(value-floor(value+0.5))<1.0e-6);
477 objValue += value*objective[iColumn];
478 if (value>columnUpper[iColumn]) {
479 if (value-columnUpper[iColumn]>1.0e-8)
480 printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]);
481 value=columnUpper[iColumn];
482 } else if (value<columnLower[iColumn]) {
483 if (value-columnLower[iColumn]<-1.0e-8)
484 printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]);
485 value=columnLower[iColumn];
486 }
487 if (value) {
488 CoinBigIndex start = columnStart[iColumn];
489 CoinBigIndex end = start + columnLength[iColumn];
490 for (CoinBigIndex j = start; j < end; j++) {
491 int iRow = row[j];
492 rowActivity[iRow] += value * element[j];
493 rowSum[iRow] += fabs(value * element[j]);
494 }
495 }
496 }
497 for (int i = 0 ; i < numberRows ; i++) {
498 #if 0 //def CLP_INVESTIGATE
499 double inf;
500 inf = rowLower[i] - rowActivity[i];
501 if (inf > primalTolerance)
502 printf("Row %d inf %g sum %g %g <= %g <= %g\n",
503 i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
504 inf = rowActivity[i] - rowUpper[i];
505 if (inf > primalTolerance)
506 printf("Row %d inf %g sum %g %g <= %g <= %g\n",
507 i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
508 #endif
509 double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
510 rowLower[i]-rowActivity[i]);
511 // but allow for errors
512 double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
513 if (infeasibility>largestInfeasibility*factor) {
514 largestInfeasibility = infeasibility/factor;
515 printf("Dinf of %g on row %d sum %g scaled %g\n",
516 infeasibility,i,rowSum[i],largestInfeasibility);
517 if (infeasibility>1.0e10) {
518 for (CoinBigIndex j=rowStart[i];
519 j<rowStart[i]+rowLength[i];j++) {
520 printf("col %d element %g\n",
521 column[j],elementByRow[j]);
522 }
523 }
524 }
525 }
526 delete [] rowActivity ;
527 delete [] rowSum;
528 if (largestInfeasibility > 10.0*primalTolerance)
529 printf("Dlargest infeasibility is %g - obj %g\n", largestInfeasibility,objValue);
530 else
531 printf("Dfeasible (%g) - obj %g\n", largestInfeasibility,objValue);
532 }
533 #endif
534 #if JUST_FIX_INTEGER
535 const double *oldLower = model->solver()->getColLower();
536 const double *oldUpper = model->solver()->getColUpper();
537 const double *dj = lp->getReducedCost();
538 int nNaturalLB = 0;
539 int nMaybeLB = 0;
540 int nForcedLB = 0;
541 int nNaturalUB = 0;
542 int nMaybeUB = 0;
543 int nForcedUB = 0;
544 int nOther = 0;
545 for (int i = 0; i < lp->getNumCols(); ++i) {
546 if (lp->isInteger(i)) {
547 if (sol[i] == oldLower[i]) {
548 if (dj[i] > 1.0e-5)
549 nNaturalLB++;
550 else if (dj[i] < -1.0e-5)
551 nForcedLB++;
552 else
553 nMaybeLB++;
554 } else if (sol[i] == oldUpper[i]) {
555 if (dj[i] < -1.0e-5)
556 nNaturalUB++;
557 else if (dj[i] > 1.0e-5)
558 nForcedUB++;
559 else
560 nMaybeUB++;
561 } else {
562 nOther++;
563 }
564 }
565 }
566 printf("%d other, LB %d natural, %d neutral, %d forced, UB %d natural, %d neutral, %d forced\n",
567 nOther, nNaturalLB, nMaybeLB, nForcedLB,
568 nNaturalUB, nMaybeUB, nForcedUB = 0);
569 #endif
570 }
571
572 TERMINATE:
573 delete[] realObj;
574 delete lp;
575 return status;
576 }
577 #undef STR_SIZE
578
579 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
580 */
581