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

Last change on this file since 2013 was 2013, checked in by forrest, 5 years ago

mods for mipstart

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