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

Last change on this file since 1943 was 1943, checked in by forrest, 8 years ago

more options, copy statistics structure analysis
start coding of "switch" variables i.e. badly scaled ints or hi/lo
changes to allow more influence on small branch and bound
changes to get correct printout with threads

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