Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/src/CbcMipStartIO.cpp

    r2012 r2013  
    2424
    2525   for ( size_t i=0 ; i<l ; ++i )
    26      if (!(isdigit(str[i])||(str[i]=='.')||(str[i]=='-')))
     26     if (!(isdigit(str[i])||(str[i]=='.')||(str[i]=='-')||(str[i]=='e')))
    2727         return false;
    2828
     
    126126   map< string, int > colIdx;
    127127   assert( ((int)colNames.size()) == lp->getNumCols() );
    128 
    129128   /* for fast search of column names */
    130129   for ( int i=0 ; (i<(int)colNames.size()) ; ++i )
     
    136135   char colNotFound[256] = "";
    137136   int nContinuousFixed = 0;
     137#ifndef JUST_FIX_INTEGER
     138#define JUST_FIX_INTEGER 0
     139#endif
     140#if JUST_FIX_INTEGER > 1
     141   // all not mentioned are at zero
     142   for ( int i=0 ; (i<lp->getNumCols()) ; ++i )
     143     {
     144       if (lp->isInteger(i))
     145         lp->setColBounds( i, 0.0, 0.0 );
     146     }
     147#endif
    138148   for ( int i=0 ; (i<(int)colValues.size()) ; ++i )
    139149   {
     
    149159         const int idx = mIt->second;
    150160         double v = colValues[i].second;
     161#if JUST_FIX_INTEGER
     162         if (!lp->isInteger(idx))
     163           continue;
     164#endif
    151165         if (v<1e-8)
    152166            v = 0.0;
     
    173187        << printLine << CoinMessageEol;
    174188   }
    175 
     189#if JUST_FIX_INTEGER
     190   lp->setHintParam(OsiDoPresolveInInitial, true, OsiHintDo) ;
     191#endif
     192   lp->setDblParam(OsiDualObjectiveLimit,COIN_DBL_MAX);
    176193   lp->initialSolve();
     194   //lp->writeMps("fixed","mps");
    177195   if (!lp->isProvenOptimal())
    178196   {
     
    263281      model->messageHandler()->message(CBC_GENERAL, model->messages())
    264282        << printLine << CoinMessageEol;
     283#if 0
     284      {
     285        int numberColumns=lp->getNumCols();
     286        double largestInfeasibility = 0.0;
     287        double primalTolerance ;
     288        double offset;
     289        lp->getDblParam(OsiObjOffset, offset);
     290        lp->getDblParam(OsiPrimalTolerance, primalTolerance) ;
     291        const double *objective = lp->getObjCoefficients() ;
     292        const double * rowLower = lp->getRowLower() ;
     293        const double * rowUpper = lp->getRowUpper() ;
     294        const double * columnLower = lp->getColLower() ;
     295        const double * columnUpper = lp->getColUpper() ;
     296        int numberRows = lp->getNumRows() ;
     297        double *rowActivity = new double[numberRows] ;
     298        memset(rowActivity, 0, numberRows*sizeof(double)) ;
     299        double *rowSum = new double[numberRows] ;
     300        memset(rowSum, 0, numberRows*sizeof(double)) ;
     301        const double * element = lp->getMatrixByCol()->getElements();
     302        const int * row = lp->getMatrixByCol()->getIndices();
     303        const CoinBigIndex * columnStart = lp->getMatrixByCol()->getVectorStarts();
     304        const int * columnLength = lp->getMatrixByCol()->getVectorLengths();
     305        const CoinPackedMatrix * rowCopy = lp->getMatrixByRow();
     306        const int * column = rowCopy->getIndices();
     307        const int * rowLength = rowCopy->getVectorLengths();
     308        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     309        const double * elementByRow = rowCopy->getElements();
     310        double objValue=-offset;
     311        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     312          double value = sol[iColumn];
     313          if (lp->isInteger(iColumn))
     314            assert (fabs(value-floor(value+0.5))<1.0e-6);
     315          objValue += value*objective[iColumn];
     316          if (value>columnUpper[iColumn]) {
     317            if (value-columnUpper[iColumn]>1.0e-8)
     318              printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]);
     319            value=columnUpper[iColumn];
     320          } else if (value<columnLower[iColumn]) {
     321            if (value-columnLower[iColumn]<-1.0e-8)
     322              printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]);
     323            value=columnLower[iColumn];
     324          }
     325          if (value) {
     326            CoinBigIndex start = columnStart[iColumn];
     327            CoinBigIndex end = start + columnLength[iColumn];
     328            for (CoinBigIndex j = start; j < end; j++) {
     329              int iRow = row[j];
     330              if (fabs(value)<1.0e-6&&fabs(value*element[j])>1.0e-5)
     331                printf("Column %d row %d value %.8g element %g %s\n",
     332                       iColumn,iRow,value,element[j],lp->isInteger(iColumn) ? "integer" : "");
     333              rowActivity[iRow] += value * element[j];
     334              rowSum[iRow] += fabs(value * element[j]);
     335            }
     336          }
     337        }
     338        for (int i = 0 ; i < numberRows ; i++) {
     339#if 0 //def CLP_INVESTIGATE
     340          double inf;
     341          inf = rowLower[i] - rowActivity[i];
     342          if (inf > primalTolerance)
     343            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     344                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     345          inf = rowActivity[i] - rowUpper[i];
     346          if (inf > primalTolerance)
     347            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     348                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     349#endif
     350          double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
     351                                         rowLower[i]-rowActivity[i]);
     352          // but allow for errors
     353          double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
     354          if (infeasibility>largestInfeasibility*factor) {
     355            largestInfeasibility = infeasibility/factor;
     356            printf("Cinf of %g on row %d sum %g scaled %g\n",
     357                   infeasibility,i,rowSum[i],largestInfeasibility);
     358            if (infeasibility>1.0e10) {
     359              for (CoinBigIndex j=rowStart[i];
     360                   j<rowStart[i]+rowLength[i];j++) {
     361                printf("col %d element %g\n",
     362                       column[j],elementByRow[j]);
     363              }
     364            }
     365          }
     366        }
     367        delete [] rowActivity ;
     368        delete [] rowSum;
     369        if (largestInfeasibility > 10.0*primalTolerance)
     370          printf("Clargest infeasibility is %g - obj %g\n", largestInfeasibility,objValue);
     371        else
     372          printf("Cfeasible (%g) - obj %g\n", largestInfeasibility,objValue);
     373      }
     374#endif
    265375      for ( int i=0 ; (i<lp->getNumCols()) ; ++i )
    266376      {
     377#if 0
    267378         if (sol[i]<1e-8)
    268379            sol[i] = 0.0;
     
    270381            if (lp->isInteger(i))
    271382               sol[i] = floor( sol[i]+0.5 );
    272       }
     383#else
     384         if (lp->isInteger(i)) {
     385           //if (fabs(sol[i] - floor( sol[i]+0.5 ))>1.0e-8)
     386           //printf("bad sol for %d - %.12g\n",i,sol[i]);
     387           sol[i] = floor( sol[i]+0.5 );
     388         }
     389#endif
     390      }
     391#if 0
     392      {
     393        int numberColumns=lp->getNumCols();
     394        double largestInfeasibility = 0.0;
     395        double primalTolerance ;
     396        double offset;
     397        lp->getDblParam(OsiObjOffset, offset);
     398        lp->getDblParam(OsiPrimalTolerance, primalTolerance) ;
     399        const double *objective = lp->getObjCoefficients() ;
     400        const double * rowLower = lp->getRowLower() ;
     401        const double * rowUpper = lp->getRowUpper() ;
     402        const double * columnLower = lp->getColLower() ;
     403        const double * columnUpper = lp->getColUpper() ;
     404        int numberRows = lp->getNumRows() ;
     405        double *rowActivity = new double[numberRows] ;
     406        memset(rowActivity, 0, numberRows*sizeof(double)) ;
     407        double *rowSum = new double[numberRows] ;
     408        memset(rowSum, 0, numberRows*sizeof(double)) ;
     409        const double * element = lp->getMatrixByCol()->getElements();
     410        const int * row = lp->getMatrixByCol()->getIndices();
     411        const CoinBigIndex * columnStart = lp->getMatrixByCol()->getVectorStarts();
     412        const int * columnLength = lp->getMatrixByCol()->getVectorLengths();
     413        const CoinPackedMatrix * rowCopy = lp->getMatrixByRow();
     414        const int * column = rowCopy->getIndices();
     415        const int * rowLength = rowCopy->getVectorLengths();
     416        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     417        const double * elementByRow = rowCopy->getElements();
     418        double objValue=-offset;
     419        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     420          double value = sol[iColumn];
     421          if (lp->isInteger(iColumn))
     422            assert (fabs(value-floor(value+0.5))<1.0e-6);
     423          objValue += value*objective[iColumn];
     424          if (value>columnUpper[iColumn]) {
     425            if (value-columnUpper[iColumn]>1.0e-8)
     426              printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]);
     427            value=columnUpper[iColumn];
     428          } else if (value<columnLower[iColumn]) {
     429            if (value-columnLower[iColumn]<-1.0e-8)
     430              printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]);
     431            value=columnLower[iColumn];
     432          }
     433          if (value) {
     434            CoinBigIndex start = columnStart[iColumn];
     435            CoinBigIndex end = start + columnLength[iColumn];
     436            for (CoinBigIndex j = start; j < end; j++) {
     437              int iRow = row[j];
     438              rowActivity[iRow] += value * element[j];
     439              rowSum[iRow] += fabs(value * element[j]);
     440            }
     441          }
     442        }
     443        for (int i = 0 ; i < numberRows ; i++) {
     444#if 0 //def CLP_INVESTIGATE
     445          double inf;
     446          inf = rowLower[i] - rowActivity[i];
     447          if (inf > primalTolerance)
     448            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     449                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     450          inf = rowActivity[i] - rowUpper[i];
     451          if (inf > primalTolerance)
     452            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     453                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     454#endif
     455          double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
     456                                         rowLower[i]-rowActivity[i]);
     457          // but allow for errors
     458          double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
     459          if (infeasibility>largestInfeasibility*factor) {
     460            largestInfeasibility = infeasibility/factor;
     461            printf("Dinf of %g on row %d sum %g scaled %g\n",
     462                   infeasibility,i,rowSum[i],largestInfeasibility);
     463            if (infeasibility>1.0e10) {
     464              for (CoinBigIndex j=rowStart[i];
     465                   j<rowStart[i]+rowLength[i];j++) {
     466                printf("col %d element %g\n",
     467                       column[j],elementByRow[j]);
     468              }
     469            }
     470          }
     471        }
     472        delete [] rowActivity ;
     473        delete [] rowSum;
     474        if (largestInfeasibility > 10.0*primalTolerance)
     475          printf("Dlargest infeasibility is %g - obj %g\n", largestInfeasibility,objValue);
     476        else
     477          printf("Dfeasible (%g) - obj %g\n", largestInfeasibility,objValue);
     478      }
     479#endif
     480#if JUST_FIX_INTEGER
     481      const double * oldLower = model->solver()->getColLower();
     482      const double * oldUpper = model->solver()->getColUpper();
     483      const double * dj = lp->getReducedCost();
     484      int nNaturalLB=0;
     485      int nMaybeLB=0;
     486      int nForcedLB=0;
     487      int nNaturalUB=0;
     488      int nMaybeUB=0;
     489      int nForcedUB=0;
     490      int nOther=0;
     491      for ( int i=0 ; i<lp->getNumCols() ; ++i ) {
     492        if (lp->isInteger(i)) {
     493          if (sol[i]==oldLower[i]) {
     494            if (dj[i]>1.0e-5)
     495              nNaturalLB++;
     496            else if (dj[i]<-1.0e-5)
     497              nForcedLB++;
     498            else
     499              nMaybeLB++;
     500          } else if (sol[i]==oldUpper[i]) {
     501            if (dj[i]<-1.0e-5)
     502              nNaturalUB++;
     503            else if (dj[i]>1.0e-5)
     504              nForcedUB++;
     505            else
     506              nMaybeUB++;
     507          } else {
     508            nOther++;
     509          }
     510        }
     511      }
     512      printf("%d other, LB %d natural, %d neutral, %d forced, UB %d natural, %d neutral, %d forced\n",
     513             nOther,nNaturalLB,nMaybeLB,nForcedLB,
     514             nNaturalUB,nMaybeUB,nForcedUB=0);
     515#endif
    273516   }
    274517
Note: See TracChangeset for help on using the changeset viewer.