source: trunk/Cbc/src/CbcMipStartIO.cpp @ 1957

Last change on this file since 1957 was 1957, checked in by forrest, 6 years ago

minor fixes and take out printout

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.6 KB
Line 
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 "CoinTime.hpp"
17
18using namespace std;
19
20
21bool 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]=='e')))
27         return false;
28
29   return true;
30}
31
32int 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   {
46      ++nLine;
47      char col[4][STR_SIZE];
48      int nread = sscanf( line, "%s %s %s %s", col[0], col[1], col[2], col[3] );
49      if (!nread)
50         continue;
51      /* line with variable value */
52      if (strlen(col[0])&&isdigit(col[0][0])&&(nread>=3))
53      {
54         if (!isNumericStr(col[0]))
55         {
56            sprintf( printLine, "Reading: %s, line %d - first column in mipstart file should be numeric, ignoring.", fileName, nLine );
57            model->messageHandler()->message(CBC_GENERAL, model->messages())
58              << printLine << CoinMessageEol;
59            continue;
60         }
61         if (!isNumericStr(col[2]))
62         {
63            sprintf( printLine, "Reading: %s, line %d - Third column in mipstart file should be numeric, ignoring.", fileName, nLine  );
64            model->messageHandler()->message(CBC_GENERAL, model->messages())
65              << printLine << CoinMessageEol;
66            continue;
67         }
68
69         //int idx = atoi( col[0] );
70         char *name = col[1];
71         double value = atof( col[2] );
72         //double obj = 0.0;
73//         if (nread >= 4)
74//            obj = atof( col[3] );
75
76         colValues.push_back( pair<string, double>(string(name),value) );
77      }
78   }
79
80   if (colValues.size()) {
81      sprintf( printLine,"mipstart values read for %d variables.", (int)colValues.size());
82      model->messageHandler()->message(CBC_GENERAL, model->messages())
83        << printLine << CoinMessageEol;
84   } else
85   {
86      sprintf( printLine, "No mipstart solution read from %s", fileName );
87      model->messageHandler()->message(CBC_GENERAL, model->messages())
88        << printLine << CoinMessageEol;
89      return 1;
90   }
91
92   fclose(f);
93   return 0;
94}
95
96int computeCompleteSolution( CbcModel * model,
97                             const vector< string > colNames,
98                             const std::vector< std::pair< std::string, double > > &colValues,
99                             double *sol, double &obj )
100{
101   int status = 0;
102   double compObj = COIN_DBL_MAX;
103   bool foundIntegerSol = false;
104   OsiSolverInterface *lp = model->solver()->clone();
105   map< string, int > colIdx;
106   assert( ((int)colNames.size()) == lp->getNumCols() );
107
108   /* for fast search of column names */
109   for ( int i=0 ; (i<(int)colNames.size()) ; ++i )
110      colIdx[colNames[i]] = i;
111
112   char printLine[STR_SIZE];
113   int fixed = 0;
114   int notFound = 0;
115   char colNotFound[256] = "";
116   int nContinuousFixed = 0;
117#ifndef JUST_FIX_INTEGER
118#define JUST_FIX_INTEGER 0
119#endif
120#if JUST_FIX_INTEGER > 1
121   // all not mentioned are at zero
122   for ( int i=0 ; (i<lp->getNumCols()) ; ++i )
123     {
124       if (lp->isInteger(i))
125         lp->setColBounds( i, 0.0, 0.0 );
126     }
127#endif
128   for ( int i=0 ; (i<(int)colValues.size()) ; ++i )
129   {
130      map< string, int >::const_iterator mIt = colIdx.find( colValues[i].first );
131      if ( mIt == colIdx.end() )
132      {
133         if (!notFound)
134            strcpy( colNotFound, colValues[i].first.c_str() );
135         notFound++;
136      }
137      else
138      {
139         const int idx = mIt->second;
140         double v = colValues[i].second;
141#if JUST_FIX_INTEGER
142         if (!lp->isInteger(idx))
143           continue;
144#endif
145         if (v<1e-8)
146            v = 0.0;
147         if (lp->isInteger(idx))  // just to avoid small
148            v = floor( v+0.5 );   // fractional garbage
149         else
150           nContinuousFixed++;
151         lp->setColBounds( idx, v, v );
152         ++fixed;
153      }
154   }
155
156   if (!fixed)
157   {
158      model->messageHandler()->message(CBC_GENERAL, model->messages())
159        << "Warning: MIPstart solution is not valid, ignoring it."
160        << CoinMessageEol;
161      goto TERMINATE;
162   }
163
164   if ( notFound >= ( ((double)colNames.size()) * 0.5 ) ) {
165      sprintf( printLine, "Warning: %d column names were not found (e.g. %s) while filling solution.", notFound, colNotFound );
166      model->messageHandler()->message(CBC_GENERAL, model->messages())
167        << printLine << CoinMessageEol;
168   }
169#if JUST_FIX_INTEGER
170   lp->setHintParam(OsiDoPresolveInInitial, true, OsiHintDo) ;
171#endif
172   lp->initialSolve();
173   //lp->writeMps("fixed","mps");
174   if (!lp->isProvenOptimal())
175   {
176      model->messageHandler()->message(CBC_GENERAL, model->messages())
177        << "Warning: mipstart values could not be used to build a solution." << CoinMessageEol;
178      if (nContinuousFixed) {
179        model->messageHandler()->message(CBC_GENERAL, model->messages())
180          << "Trying just fixing integer variables." << CoinMessageEol;
181        int numberColumns = lp->getNumCols();
182        const double * oldLower = model->solver()->getColLower();
183        const double * oldUpper = model->solver()->getColUpper();
184        for ( int i=0 ; i<numberColumns ; ++i ) {
185          if (!lp->isInteger(i)) {
186            lp->setColLower(i,oldLower[i]);
187            lp->setColUpper(i,oldUpper[i]);
188          }
189        }
190        lp->initialSolve();
191        if (!lp->isProvenOptimal())
192          model->messageHandler()->message(CBC_GENERAL, model->messages())
193            << "Still no good." << CoinMessageEol;
194      }
195      if (!lp->isProvenOptimal()) {
196        status = 1;
197        goto TERMINATE;
198      }
199   }
200
201   /* some additional effort is needed to provide an integer solution */
202   if ( lp->getFractionalIndices().size() > 0 )
203   {
204      sprintf( printLine,"MIPStart solution provided values for %d of %d integer variables, %d variables are still fractional.", fixed, lp->getNumIntegers(), (int)lp->getFractionalIndices().size() );
205      model->messageHandler()->message(CBC_GENERAL, model->messages())
206        << printLine << CoinMessageEol;
207      double start = CoinCpuTime();
208#if 1
209      CbcSerendipity heuristic(*model);
210      heuristic.setFractionSmall(2.0);
211      heuristic.setFeasibilityPumpOptions(1008013);
212      int returnCode = heuristic.smallBranchAndBound(lp,
213                                                     1000, sol,
214                                                     compObj,
215                                                     model->getCutoff(), 
216                                                     "ReduceInMIPStart");
217      if ((returnCode&1) != 0) {
218         sprintf( printLine,"Mini branch and bound defined values for remaining variables in %.2f seconds.", 
219                  CoinCpuTime()-start);
220         model->messageHandler()->message(CBC_GENERAL, model->messages())
221           << printLine << CoinMessageEol;
222         foundIntegerSol = true;
223         obj = compObj;
224      }
225#else
226      CbcModel babModel( *lp );
227      babModel.setLogLevel( 0 );
228      babModel.setMaximumNodes( 500 );
229      babModel.setMaximumSeconds( 60 );
230      babModel.branchAndBound();
231      if (babModel.bestSolution())
232      {
233         sprintf( printLine,"Mini branch and bound defined values for remaining variables in %.2f seconds.", 
234                  CoinCpuTime()-start);
235         model->messageHandler()->message(CBC_GENERAL, model->messages())
236           << printLine << CoinMessageEol;
237         copy( babModel.bestSolution(), babModel.bestSolution()+babModel.getNumCols(), sol );
238         foundIntegerSol = true;
239         obj = compObj = babModel.getObjValue();
240      }
241#endif
242      else
243      {
244         model->messageHandler()->message(CBC_GENERAL, model->messages())
245           << "Warning: mipstart values could not be used to build a solution." << CoinMessageEol;
246         status = 1;
247         goto TERMINATE;
248      }
249   }
250   else
251   {
252      foundIntegerSol = true;
253      obj = compObj = lp->getObjValue();
254      copy( lp->getColSolution(), lp->getColSolution()+lp->getNumCols(), sol );
255   }
256
257   if ( foundIntegerSol )
258   {
259      sprintf( printLine,"mipstart provided solution with cost %g", compObj);
260      model->messageHandler()->message(CBC_GENERAL, model->messages())
261        << printLine << CoinMessageEol;
262#if 0
263      {
264        int numberColumns=lp->getNumCols();
265        double largestInfeasibility = 0.0;
266        double primalTolerance ;
267        double offset;
268        lp->getDblParam(OsiObjOffset, offset);
269        lp->getDblParam(OsiPrimalTolerance, primalTolerance) ;
270        const double *objective = lp->getObjCoefficients() ;
271        const double * rowLower = lp->getRowLower() ;
272        const double * rowUpper = lp->getRowUpper() ;
273        const double * columnLower = lp->getColLower() ;
274        const double * columnUpper = lp->getColUpper() ;
275        int numberRows = lp->getNumRows() ;
276        double *rowActivity = new double[numberRows] ;
277        memset(rowActivity, 0, numberRows*sizeof(double)) ;
278        double *rowSum = new double[numberRows] ;
279        memset(rowSum, 0, numberRows*sizeof(double)) ;
280        const double * element = lp->getMatrixByCol()->getElements();
281        const int * row = lp->getMatrixByCol()->getIndices();
282        const CoinBigIndex * columnStart = lp->getMatrixByCol()->getVectorStarts();
283        const int * columnLength = lp->getMatrixByCol()->getVectorLengths();
284        const CoinPackedMatrix * rowCopy = lp->getMatrixByRow();
285        const int * column = rowCopy->getIndices();
286        const int * rowLength = rowCopy->getVectorLengths();
287        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
288        const double * elementByRow = rowCopy->getElements();
289        double objValue=-offset;
290        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
291          double value = sol[iColumn];
292          if (lp->isInteger(iColumn))
293            assert (fabs(value-floor(value+0.5))<1.0e-6);
294          objValue += value*objective[iColumn];
295          if (value>columnUpper[iColumn]) {
296            if (value-columnUpper[iColumn]>1.0e-8)
297              printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]);
298            value=columnUpper[iColumn];
299          } else if (value<columnLower[iColumn]) {
300            if (value-columnLower[iColumn]<-1.0e-8)
301              printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]);
302            value=columnLower[iColumn];
303          }
304          if (value) {
305            CoinBigIndex start = columnStart[iColumn];
306            CoinBigIndex end = start + columnLength[iColumn];
307            for (CoinBigIndex j = start; j < end; j++) {
308              int iRow = row[j];
309              if (fabs(value)<1.0e-6&&fabs(value*element[j])>1.0e-5)
310                printf("Column %d row %d value %.8g element %g %s\n",
311                       iColumn,iRow,value,element[j],lp->isInteger(iColumn) ? "integer" : "");
312              rowActivity[iRow] += value * element[j];
313              rowSum[iRow] += fabs(value * element[j]);
314            }
315          }
316        }
317        for (int i = 0 ; i < numberRows ; i++) {
318#if 0 //def CLP_INVESTIGATE
319          double inf;
320          inf = rowLower[i] - rowActivity[i];
321          if (inf > primalTolerance)
322            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
323                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
324          inf = rowActivity[i] - rowUpper[i];
325          if (inf > primalTolerance)
326            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
327                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
328#endif
329          double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
330                                         rowLower[i]-rowActivity[i]);
331          // but allow for errors
332          double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
333          if (infeasibility>largestInfeasibility*factor) {
334            largestInfeasibility = infeasibility/factor;
335            printf("Cinf of %g on row %d sum %g scaled %g\n",
336                   infeasibility,i,rowSum[i],largestInfeasibility);
337            if (infeasibility>1.0e10) {
338              for (CoinBigIndex j=rowStart[i];
339                   j<rowStart[i]+rowLength[i];j++) {
340                printf("col %d element %g\n",
341                       column[j],elementByRow[j]);
342              }
343            }
344          }
345        }
346        delete [] rowActivity ;
347        delete [] rowSum;
348        if (largestInfeasibility > 10.0*primalTolerance)
349          printf("Clargest infeasibility is %g - obj %g\n", largestInfeasibility,objValue);
350        else
351          printf("Cfeasible (%g) - obj %g\n", largestInfeasibility,objValue);
352      }
353#endif
354      for ( int i=0 ; (i<lp->getNumCols()) ; ++i )
355      {
356#if 0
357         if (sol[i]<1e-8)
358            sol[i] = 0.0;
359         else
360            if (lp->isInteger(i))
361               sol[i] = floor( sol[i]+0.5 );
362#else
363         if (lp->isInteger(i)) {
364           //if (fabs(sol[i] - floor( sol[i]+0.5 ))>1.0e-8)
365           //printf("bad sol for %d - %.12g\n",i,sol[i]);
366           sol[i] = floor( sol[i]+0.5 );
367         }
368#endif
369      }
370#if 0
371      {
372        int numberColumns=lp->getNumCols();
373        double largestInfeasibility = 0.0;
374        double primalTolerance ;
375        double offset;
376        lp->getDblParam(OsiObjOffset, offset);
377        lp->getDblParam(OsiPrimalTolerance, primalTolerance) ;
378        const double *objective = lp->getObjCoefficients() ;
379        const double * rowLower = lp->getRowLower() ;
380        const double * rowUpper = lp->getRowUpper() ;
381        const double * columnLower = lp->getColLower() ;
382        const double * columnUpper = lp->getColUpper() ;
383        int numberRows = lp->getNumRows() ;
384        double *rowActivity = new double[numberRows] ;
385        memset(rowActivity, 0, numberRows*sizeof(double)) ;
386        double *rowSum = new double[numberRows] ;
387        memset(rowSum, 0, numberRows*sizeof(double)) ;
388        const double * element = lp->getMatrixByCol()->getElements();
389        const int * row = lp->getMatrixByCol()->getIndices();
390        const CoinBigIndex * columnStart = lp->getMatrixByCol()->getVectorStarts();
391        const int * columnLength = lp->getMatrixByCol()->getVectorLengths();
392        const CoinPackedMatrix * rowCopy = lp->getMatrixByRow();
393        const int * column = rowCopy->getIndices();
394        const int * rowLength = rowCopy->getVectorLengths();
395        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
396        const double * elementByRow = rowCopy->getElements();
397        double objValue=-offset;
398        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
399          double value = sol[iColumn];
400          if (lp->isInteger(iColumn))
401            assert (fabs(value-floor(value+0.5))<1.0e-6);
402          objValue += value*objective[iColumn];
403          if (value>columnUpper[iColumn]) {
404            if (value-columnUpper[iColumn]>1.0e-8)
405              printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]);
406            value=columnUpper[iColumn];
407          } else if (value<columnLower[iColumn]) {
408            if (value-columnLower[iColumn]<-1.0e-8)
409              printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]);
410            value=columnLower[iColumn];
411          }
412          if (value) {
413            CoinBigIndex start = columnStart[iColumn];
414            CoinBigIndex end = start + columnLength[iColumn];
415            for (CoinBigIndex j = start; j < end; j++) {
416              int iRow = row[j];
417              rowActivity[iRow] += value * element[j];
418              rowSum[iRow] += fabs(value * element[j]);
419            }
420          }
421        }
422        for (int i = 0 ; i < numberRows ; i++) {
423#if 0 //def CLP_INVESTIGATE
424          double inf;
425          inf = rowLower[i] - rowActivity[i];
426          if (inf > primalTolerance)
427            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
428                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
429          inf = rowActivity[i] - rowUpper[i];
430          if (inf > primalTolerance)
431            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
432                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
433#endif
434          double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
435                                         rowLower[i]-rowActivity[i]);
436          // but allow for errors
437          double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
438          if (infeasibility>largestInfeasibility*factor) {
439            largestInfeasibility = infeasibility/factor;
440            printf("Dinf of %g on row %d sum %g scaled %g\n",
441                   infeasibility,i,rowSum[i],largestInfeasibility);
442            if (infeasibility>1.0e10) {
443              for (CoinBigIndex j=rowStart[i];
444                   j<rowStart[i]+rowLength[i];j++) {
445                printf("col %d element %g\n",
446                       column[j],elementByRow[j]);
447              }
448            }
449          }
450        }
451        delete [] rowActivity ;
452        delete [] rowSum;
453        if (largestInfeasibility > 10.0*primalTolerance)
454          printf("Dlargest infeasibility is %g - obj %g\n", largestInfeasibility,objValue);
455        else
456          printf("Dfeasible (%g) - obj %g\n", largestInfeasibility,objValue);
457      }
458#endif
459#if JUST_FIX_INTEGER
460      const double * oldLower = model->solver()->getColLower();
461      const double * oldUpper = model->solver()->getColUpper();
462      const double * dj = lp->getReducedCost();
463      int nNaturalLB=0;
464      int nMaybeLB=0;
465      int nForcedLB=0;
466      int nNaturalUB=0;
467      int nMaybeUB=0;
468      int nForcedUB=0;
469      int nOther=0;
470      for ( int i=0 ; i<lp->getNumCols() ; ++i ) {
471        if (lp->isInteger(i)) {
472          if (sol[i]==oldLower[i]) {
473            if (dj[i]>1.0e-5)
474              nNaturalLB++;
475            else if (dj[i]<-1.0e-5)
476              nForcedLB++;
477            else
478              nMaybeLB++;
479          } else if (sol[i]==oldUpper[i]) {
480            if (dj[i]<-1.0e-5)
481              nNaturalUB++;
482            else if (dj[i]>1.0e-5)
483              nForcedUB++;
484            else
485              nMaybeUB++;
486          } else {
487            nOther++;
488          }
489        }
490      }
491      printf("%d other, LB %d natural, %d neutral, %d forced, UB %d natural, %d neutral, %d forced\n",
492             nOther,nNaturalLB,nMaybeLB,nForcedLB,
493             nNaturalUB,nMaybeUB,nForcedUB=0);
494#endif
495   }
496
497TERMINATE:
498   delete lp;
499   return status;
500}
501#undef STR_SIZE
Note: See TracBrowser for help on using the repository browser.