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

Last change on this file since 2385 was 2385, checked in by forrest, 2 years ago

mip start sos

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.9 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 (and fixingish SOS)." << CoinMessageEol;
198         int numberColumns = lp->getNumCols();
199         const double *oldLower = model->solver()->getColLower();
200         const double *oldUpper = model->solver()->getColUpper();
201         double * savedSol = CoinCopyOfArray(lp->getColLower(),numberColumns);
202         for ( int i=0 ; i<numberColumns ; ++i ) {
203            if (!lp->isInteger(i)) { 
204               lp->setColLower(i,oldLower[i]);
205               lp->setColUpper(i,oldUpper[i]);
206            }
207         }
208         // but look at SOS
209         int numberObjects = model->numberObjects();
210         for (int i=0;i<numberObjects;i++) {
211           const OsiSOS * object =
212             dynamic_cast<const OsiSOS *>(model->object(i));
213           if (object) {
214             int n=object->numberMembers();
215             const int * members = object->members();
216             int sosType = object->sosType();
217             if (sosType==1) {
218               // non zero can take any value - others zero
219               int iColumn=-1;
220               for (int j=0;j<n;j++) {
221                 int jColumn=members[j];
222                 if (savedSol[jColumn]) 
223                   iColumn = jColumn;
224               }
225               for (int j=0;j<n;j++) {
226                 int jColumn=members[j];
227                 if (jColumn!=iColumn) {
228                   lp->setColLower(jColumn,0.0);
229                   lp->setColUpper(jColumn,0.0);
230                 }
231               }
232             } else if (sosType==2) {
233               // SOS 2 - make a guess if just one nonzero
234               int jA=-1;
235               int jB=-1;
236               for (int j=0;j<n;j++) {
237                 int jColumn=members[j];
238                 if (savedSol[jColumn]) {
239                   if (jA==-1)
240                     jA = j;
241                   jB = j;
242                 }
243               }
244               if (jB>jA+1) {
245                 jB = jA+1;
246               } else if (jA==jB) {
247                 if (jA==n-1)
248                   jA--;
249                 else
250                   jB++;
251               }
252               for (int j=0;j<n;j++) {
253                 if (j!=jA&&j!=jB) {
254                   int jColumn=members[j];
255                   lp->setColLower(jColumn,0.0);
256                   lp->setColUpper(jColumn,0.0);
257                 }
258               }
259             }
260           }
261         }
262         delete [] savedSol;
263         lp->initialSolve();
264      }
265      else
266      {
267         model->messageHandler()->message(CBC_GENERAL, model->messages())
268             << "Fixing only non-zero variables." << CoinMessageEol;
269         /* unfix all variables which are zero */
270         int notZeroAnymore = 0;
271         for ( int i=0 ; (i<lp->getNumCols()) ; ++i )
272             if ( ((fabs(lp->getColLower()[i])) <= 1e-8) && (fabs(lp->getColLower()[i]-lp->getColUpper()[i]) <= 1e-8) )
273             {
274                const double *oldLower = model->solver()->getColLower();
275                const double *oldUpper = model->solver()->getColUpper();
276                lp->setColLower(i,oldLower[i]);
277                lp->setColUpper(i,oldUpper[i]);
278                notZeroAnymore++;
279             }
280         if (notZeroAnymore)
281             lp->initialSolve();
282      }
283   }
284
285   if (!lp->isProvenOptimal())
286   {
287      model->messageHandler()->message(CBC_GENERAL, model->messages())
288           << "Warning: mipstart values could not be used to build a solution." << CoinMessageEol;
289      status = 1;
290      goto TERMINATE;
291   }
292   
293   /* some additional effort is needed to provide an integer solution */
294   if ( lp->getFractionalIndices().size() > 0 )
295   {
296      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()) );
297      model->messageHandler()->message(CBC_GENERAL, model->messages())
298        << printLine << CoinMessageEol;
299      double start = CoinCpuTime();
300#if 1
301      CbcSerendipity heuristic(*model);
302      heuristic.setFractionSmall(2.0);
303      heuristic.setFeasibilityPumpOptions(1008013);
304      int returnCode = heuristic.smallBranchAndBound(lp,
305                                                     1000, sol,
306                                                     compObj,
307                                                     model->getCutoff(), 
308                                                     "ReduceInMIPStart");
309      if ((returnCode&1) != 0) {
310         sprintf( printLine,"Mini branch and bound defined values for remaining variables in %.2f seconds.", 
311                  CoinCpuTime()-start);
312         model->messageHandler()->message(CBC_GENERAL, model->messages())
313           << printLine << CoinMessageEol;
314         foundIntegerSol = true;
315         obj = compObj;
316      }
317#else
318      CbcModel babModel( *lp );
319      lp->writeLp("lessFix");
320      babModel.setLogLevel( 2 );
321      babModel.setMaximumNodes( 1000 );
322      babModel.setMaximumSeconds( 60 );
323      babModel.branchAndBound();
324      if (babModel.bestSolution())
325      {
326         sprintf( printLine,"Mini branch and bound defined values for remaining variables in %.2f seconds.", 
327                  CoinCpuTime()-start);
328         model->messageHandler()->message(CBC_GENERAL, model->messages())
329           << printLine << CoinMessageEol;
330         copy( babModel.bestSolution(), babModel.bestSolution()+babModel.getNumCols(), sol );
331         foundIntegerSol = true;
332         obj = compObj = babModel.getObjValue();
333      }
334#endif
335      else
336      {
337          model->messageHandler()->message(CBC_GENERAL, model->messages())
338              << "Warning: mipstart values could not be used to build a solution." << CoinMessageEol;
339         status = 1;
340         goto TERMINATE;
341      }
342   }
343   else
344   {
345      foundIntegerSol = true;
346      obj = compObj = lp->getObjValue();
347      copy( lp->getColSolution(), lp->getColSolution()+lp->getNumCols(), sol );
348   }
349
350   if ( foundIntegerSol )
351   {
352      sprintf( printLine,"MIPStart provided solution with cost %g", compObj);
353      model->messageHandler()->message(CBC_GENERAL, model->messages())
354           << printLine << CoinMessageEol;
355#if 0
356      {
357        int numberColumns=lp->getNumCols();
358        double largestInfeasibility = 0.0;
359        double primalTolerance ;
360        double offset;
361        lp->getDblParam(OsiObjOffset, offset);
362        lp->getDblParam(OsiPrimalTolerance, primalTolerance) ;
363        const double *objective = lp->getObjCoefficients() ;
364        const double * rowLower = lp->getRowLower() ;
365        const double * rowUpper = lp->getRowUpper() ;
366        const double * columnLower = lp->getColLower() ;
367        const double * columnUpper = lp->getColUpper() ;
368        int numberRows = lp->getNumRows() ;
369        double *rowActivity = new double[numberRows] ;
370        memset(rowActivity, 0, numberRows*sizeof(double)) ;
371        double *rowSum = new double[numberRows] ;
372        memset(rowSum, 0, numberRows*sizeof(double)) ;
373        const double * element = lp->getMatrixByCol()->getElements();
374        const int * row = lp->getMatrixByCol()->getIndices();
375        const CoinBigIndex * columnStart = lp->getMatrixByCol()->getVectorStarts();
376        const int * columnLength = lp->getMatrixByCol()->getVectorLengths();
377        const CoinPackedMatrix * rowCopy = lp->getMatrixByRow();
378        const int * column = rowCopy->getIndices();
379        const int * rowLength = rowCopy->getVectorLengths();
380        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
381        const double * elementByRow = rowCopy->getElements();
382        double objValue=-offset;
383        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
384          double value = sol[iColumn];
385          if (lp->isInteger(iColumn))
386            assert (fabs(value-floor(value+0.5))<1.0e-6);
387          objValue += value*objective[iColumn];
388          if (value>columnUpper[iColumn]) {
389            if (value-columnUpper[iColumn]>1.0e-8)
390              printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]);
391            value=columnUpper[iColumn];
392          } else if (value<columnLower[iColumn]) {
393            if (value-columnLower[iColumn]<-1.0e-8)
394              printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]);
395            value=columnLower[iColumn];
396          }
397          if (value) {
398            CoinBigIndex start = columnStart[iColumn];
399            CoinBigIndex end = start + columnLength[iColumn];
400            for (CoinBigIndex j = start; j < end; j++) {
401              int iRow = row[j];
402              if (fabs(value)<1.0e-6&&fabs(value*element[j])>1.0e-5)
403                printf("Column %d row %d value %.8g element %g %s\n",
404                       iColumn,iRow,value,element[j],lp->isInteger(iColumn) ? "integer" : "");
405              rowActivity[iRow] += value * element[j];
406              rowSum[iRow] += fabs(value * element[j]);
407            }
408          }
409        }
410        for (int i = 0 ; i < numberRows ; i++) {
411#if 0 //def CLP_INVESTIGATE
412          double inf;
413          inf = rowLower[i] - rowActivity[i];
414          if (inf > primalTolerance)
415            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
416                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
417          inf = rowActivity[i] - rowUpper[i];
418          if (inf > primalTolerance)
419            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
420                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
421#endif
422          double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
423                                         rowLower[i]-rowActivity[i]);
424          // but allow for errors
425          double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
426          if (infeasibility>largestInfeasibility*factor) {
427            largestInfeasibility = infeasibility/factor;
428            printf("Cinf of %g on row %d sum %g scaled %g\n",
429                   infeasibility,i,rowSum[i],largestInfeasibility);
430            if (infeasibility>1.0e10) {
431              for (CoinBigIndex j=rowStart[i];
432                   j<rowStart[i]+rowLength[i];j++) {
433                printf("col %d element %g\n",
434                       column[j],elementByRow[j]);
435              }
436            }
437          }
438        }
439        delete [] rowActivity ;
440        delete [] rowSum;
441        if (largestInfeasibility > 10.0*primalTolerance)
442          printf("Clargest infeasibility is %g - obj %g\n", largestInfeasibility,objValue);
443        else
444          printf("Cfeasible (%g) - obj %g\n", largestInfeasibility,objValue);
445      }
446#endif
447      for ( int i=0 ; (i<lp->getNumCols()) ; ++i )
448      {
449#if 0
450         if (sol[i]<1e-8)
451            sol[i] = 0.0;
452         else
453            if (lp->isInteger(i))
454               sol[i] = floor( sol[i]+0.5 );
455#else
456         if (lp->isInteger(i)) {
457           //if (fabs(sol[i] - floor( sol[i]+0.5 ))>1.0e-8)
458           //printf("bad sol for %d - %.12g\n",i,sol[i]);
459           sol[i] = floor( sol[i]+0.5 );
460         }
461#endif
462      }
463#if 0
464      {
465        int numberColumns=lp->getNumCols();
466        double largestInfeasibility = 0.0;
467        double primalTolerance ;
468        double offset;
469        lp->getDblParam(OsiObjOffset, offset);
470        lp->getDblParam(OsiPrimalTolerance, primalTolerance) ;
471        const double *objective = lp->getObjCoefficients() ;
472        const double * rowLower = lp->getRowLower() ;
473        const double * rowUpper = lp->getRowUpper() ;
474        const double * columnLower = lp->getColLower() ;
475        const double * columnUpper = lp->getColUpper() ;
476        int numberRows = lp->getNumRows() ;
477        double *rowActivity = new double[numberRows] ;
478        memset(rowActivity, 0, numberRows*sizeof(double)) ;
479        double *rowSum = new double[numberRows] ;
480        memset(rowSum, 0, numberRows*sizeof(double)) ;
481        const double * element = lp->getMatrixByCol()->getElements();
482        const int * row = lp->getMatrixByCol()->getIndices();
483        const CoinBigIndex * columnStart = lp->getMatrixByCol()->getVectorStarts();
484        const int * columnLength = lp->getMatrixByCol()->getVectorLengths();
485        const CoinPackedMatrix * rowCopy = lp->getMatrixByRow();
486        const int * column = rowCopy->getIndices();
487        const int * rowLength = rowCopy->getVectorLengths();
488        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
489        const double * elementByRow = rowCopy->getElements();
490        double objValue=-offset;
491        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
492          double value = sol[iColumn];
493          if (lp->isInteger(iColumn))
494            assert (fabs(value-floor(value+0.5))<1.0e-6);
495          objValue += value*objective[iColumn];
496          if (value>columnUpper[iColumn]) {
497            if (value-columnUpper[iColumn]>1.0e-8)
498              printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]);
499            value=columnUpper[iColumn];
500          } else if (value<columnLower[iColumn]) {
501            if (value-columnLower[iColumn]<-1.0e-8)
502              printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]);
503            value=columnLower[iColumn];
504          }
505          if (value) {
506            CoinBigIndex start = columnStart[iColumn];
507            CoinBigIndex end = start + columnLength[iColumn];
508            for (CoinBigIndex j = start; j < end; j++) {
509              int iRow = row[j];
510              rowActivity[iRow] += value * element[j];
511              rowSum[iRow] += fabs(value * element[j]);
512            }
513          }
514        }
515        for (int i = 0 ; i < numberRows ; i++) {
516#if 0 //def CLP_INVESTIGATE
517          double inf;
518          inf = rowLower[i] - rowActivity[i];
519          if (inf > primalTolerance)
520            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
521                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
522          inf = rowActivity[i] - rowUpper[i];
523          if (inf > primalTolerance)
524            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
525                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
526#endif
527          double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
528                                         rowLower[i]-rowActivity[i]);
529          // but allow for errors
530          double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
531          if (infeasibility>largestInfeasibility*factor) {
532            largestInfeasibility = infeasibility/factor;
533            printf("Dinf of %g on row %d sum %g scaled %g\n",
534                   infeasibility,i,rowSum[i],largestInfeasibility);
535            if (infeasibility>1.0e10) {
536              for (CoinBigIndex j=rowStart[i];
537                   j<rowStart[i]+rowLength[i];j++) {
538                printf("col %d element %g\n",
539                       column[j],elementByRow[j]);
540              }
541            }
542          }
543        }
544        delete [] rowActivity ;
545        delete [] rowSum;
546        if (largestInfeasibility > 10.0*primalTolerance)
547          printf("Dlargest infeasibility is %g - obj %g\n", largestInfeasibility,objValue);
548        else
549          printf("Dfeasible (%g) - obj %g\n", largestInfeasibility,objValue);
550      }
551#endif
552#if JUST_FIX_INTEGER
553      const double * oldLower = model->solver()->getColLower();
554      const double * oldUpper = model->solver()->getColUpper();
555      const double * dj = lp->getReducedCost();
556      int nNaturalLB=0;
557      int nMaybeLB=0;
558      int nForcedLB=0;
559      int nNaturalUB=0;
560      int nMaybeUB=0;
561      int nForcedUB=0;
562      int nOther=0;
563      for ( int i=0 ; i<lp->getNumCols() ; ++i ) {
564        if (lp->isInteger(i)) {
565          if (sol[i]==oldLower[i]) {
566            if (dj[i]>1.0e-5)
567              nNaturalLB++;
568            else if (dj[i]<-1.0e-5)
569              nForcedLB++;
570            else
571              nMaybeLB++;
572          } else if (sol[i]==oldUpper[i]) {
573            if (dj[i]<-1.0e-5)
574              nNaturalUB++;
575            else if (dj[i]>1.0e-5)
576              nForcedUB++;
577            else
578              nMaybeUB++;
579          } else {
580            nOther++;
581          }
582        }
583      }
584      printf("%d other, LB %d natural, %d neutral, %d forced, UB %d natural, %d neutral, %d forced\n",
585             nOther,nNaturalLB,nMaybeLB,nForcedLB,
586             nNaturalUB,nMaybeUB,nForcedUB=0);
587#endif
588   }
589
590TERMINATE:
591   delete lp;
592   return status;
593}
594#undef STR_SIZE
Note: See TracBrowser for help on using the repository browser.