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

Last change on this file since 2349 was 2349, checked in by forrest, 3 years ago

put back tolerance a bit earlier and fix minor bug in mipstart

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