source: stable/2.9/Cbc/src/CbcMipStartIO.cpp @ 2187

Last change on this file since 2187 was 2187, checked in by stefan, 3 years ago

more sync with trunk

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