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

Last change on this file since 2186 was 2186, checked in by stefan, 4 years ago

sync with trunk

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