Logo Search packages:      
Sourcecode: rsymphony version File versions  Download package

OsiSolverInterface * CglPreProcess::preProcessNonDefault ( OsiSolverInterface model,
int  makeEquality = 0,
int  numberPasses = 5,
int  tuning = 0 
)

preProcess problem - returning new problem. If makeEquality true then <= cliques converted to ==. Presolve will be done numberPasses times.

Returns NULL if infeasible

This version assumes user has added cut generators to CglPreProcess object before calling it. As an example use coding in preProcess If makeEquality is 1 add slacks to get cliques, if 2 add slacks to get sos (but only if looks plausible) and keep sos info

Definition at line 1046 of file CglPreProcess.cpp.

References flopc::abs(), OsiSolverInterface::addCol(), OsiSolverInterface::addCols(), CglStored::addCut(), addCutGenerator(), OsiSolverInterface::addRows(), CoinWarmStartBasis::atLowerBound, CoinWarmStartBasis::basic, flopc::ceil(), CoinPackedMatrix::cleanMatrix(), OsiSolverInterface::clone(), cuts_, OsiSolverInterface::deleteCols(), OsiSolverInterface::deleteRows(), flopc::floor(), generator_, OsiSolverInterface::getColLower(), OsiSolverInterface::getColSolution(), OsiSolverInterface::getColUpper(), OsiSolverInterface::getDblParam(), CoinPackedMatrix::getElements(), OsiSolverInterface::getEmptyWarmStart(), OsiSolverInterface::getHintParam(), CoinPackedMatrix::getIndices(), OsiSolverInterface::getIterationCount(), OsiSolverInterface::getMatrixByCol(), OsiSolverInterface::getMatrixByRow(), CoinPackedMatrix::getMutableElements(), OsiSolverInterface::getNumCols(), OsiSolverInterface::getNumElements(), CoinPackedMatrix::getNumRows(), OsiSolverInterface::getNumRows(), OsiSolverInterface::getObjCoefficients(), OsiSolverInterface::getObjSense(), OsiSolverInterface::getRowLower(), OsiSolverInterface::getRowUpper(), OsiSolverInterface::getStrParam(), CoinPackedMatrix::getVectorLengths(), CoinPackedMatrix::getVectorStarts(), OsiSolverInterface::getWarmStart(), handler_, OsiSolverInterface::initialSolve(), OsiSolverInterface::isBinary(), OsiSolverInterface::isInteger(), OsiSolverInterface::isProvenDualInfeasible(), OsiSolverInterface::isProvenOptimal(), CoinMessageHandler::logLevel(), makeInteger(), CoinMessageHandler::message(), OsiSolverInterface::messageHandler(), messages_, model_, modified(), modifiedModel_, numberCutGenerators_, numberIterationsPre_, numberProhibited_, numberRowType_, numberSolvers_, numberSOS(), numberSOS_, options_, originalColumn_, originalModel_, originalRow_, presolve_, OsiPresolve::presolvedModel(), prohibited_, CoinThreadRandom::randomDouble(), reducedCostFix(), OsiSolverInterface::resolve(), rowType_, CoinWarmStartBasis::setArtifStatus(), OsiSolverInterface::setColLower(), OsiSolverInterface::setColSolution(), OsiSolverInterface::setColUpper(), OsiSolverInterface::setDblParam(), OsiSolverInterface::setHintParam(), OsiSolverInterface::setInteger(), CoinMessageHandler::setLogLevel(), OsiSolverInterface::setObjective(), OsiPresolve::setPresolveActions(), OsiSolverInterface::setRowLower(), OsiSolverInterface::setRowUpper(), CoinWarmStartBasis::setStructStatus(), OsiSolverInterface::setWarmStart(), startModel_, startSOS_, tightenPrimalBounds(), typeSOS_, update(), weightSOS_, whichSOS_, and OsiSolverInterface::writeMps().

Referenced by preProcess().

{
  originalModel_ = & model;
  numberSolvers_ = numberPasses;
  model_ = new OsiSolverInterface * [numberSolvers_];
  modifiedModel_ = new OsiSolverInterface * [numberSolvers_];
  presolve_ = new OsiPresolve * [numberSolvers_];
  for (int i=0;i<numberSolvers_;i++) {
    model_[i]=NULL;
    modifiedModel_[i]=NULL;
    presolve_[i]=NULL;
  }
  // Put presolve option on
  tuning |= 8;
  // clear original
  delete [] originalColumn_;
  delete [] originalRow_;
  originalColumn_=NULL;
  originalRow_=NULL;
  //startModel_=&model;
  // make clone
  delete startModel_;
  startModel_ = originalModel_->clone();
  CoinPackedMatrix matrixByRow(*originalModel_->getMatrixByRow());
  int numberRows = originalModel_->getNumRows();
  if (rowType_)
    assert (numberRowType_==numberRows);
  int numberColumns = originalModel_->getNumCols();
  //int originalNumberColumns=numberColumns;
  int minimumLength = 5;
  int numberModifiedPasses=10;
  if (numberPasses<=1)
    numberModifiedPasses=1; // lightweight preprocessing
  else if (numberPasses<=2)
    numberModifiedPasses=2; // fairly lightweight preprocessing
  if (tuning>=10000) {
    numberModifiedPasses=(tuning-10000)/10000;
    tuning %= 10000;
    //minimumLength = tuning;
  }
  //bool heavyProbing = (tuning&1)!=0;
  int makeIntegers = (tuning&6)>>1;
  // See if we want to do initial presolve
  int doInitialPresolve = (tuning&8)>>3;
  if (numberSolvers_<2)
    doInitialPresolve=0;
  // We want to add columns
  int numberSlacks=0;
  int * rows = new int[numberRows];
  double * element =new double[numberRows];
  
  int iRow;
  
  int numberCliques=0;
  int * which = new int[numberColumns];

  // Statistics
  int totalP1=0,totalM1=0;
  int numberFixed=0;
  // May just find it is infeasible
  bool feasible=true;
  
  // Row copy
  const double * elementByRow = matrixByRow.getElements();
  const int * column = matrixByRow.getIndices();
  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
  const int * rowLength = matrixByRow.getVectorLengths();
  
  const double * lower = originalModel_->getColLower();
  const double * upper = originalModel_->getColUpper();
  const double * rowLower = originalModel_->getRowLower();
  const double * rowUpper = originalModel_->getRowUpper();
  // Clean bounds
  int iColumn;
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    if (originalModel_->isInteger(iColumn)) {
      double lo = CoinMax(lower[iColumn],ceil(lower[iColumn]-1.0e-6));
      if (lo>lower[iColumn])
      originalModel_->setColLower(iColumn,lo);
      double up = CoinMin(upper[iColumn],floor(upper[iColumn]+1.0e-6));
      if (up<upper[iColumn])
      originalModel_->setColUpper(iColumn,up);
      if (lo>up)
      feasible=false;
    }
  }
  bool allToGub = makeEquality==5;
  if (allToGub)
    makeEquality=3;
  // Initialize random seed
  CoinThreadRandom randomGenerator(987654321);
  bool justOnesWithObj=false;
  if (makeEquality==2||makeEquality==3||makeEquality==4) {
    int iRow, iColumn;
    int numberIntegers = 0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      if (originalModel_->isInteger(iColumn))
        numberIntegers++;
    }
    // Look for possible SOS
    int numberSOS=0;
    int * mark = new int[numberColumns];
    CoinFillN(mark,numberColumns,-1);
    int numberOverlap=0;
    int numberInSOS=0;
    // See if worthwhile creating accumulation variables
    int firstOther=numberRows;
    int * whichRow = new int[numberRows];
    for (iRow=0;iRow<numberRows;iRow++) {
      if (rowUpper[iRow]==1.0) {
        if (rowLength[iRow]<5)
          continue;
        bool goodRow=true;
      bool overlap=false;
        for (int j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
          int iColumn = column[j];
          if (elementByRow[j]!=1.0||!originalModel_->isInteger(iColumn)||lower[iColumn]) {
            goodRow=false;
            break;
          }
          if (mark[iColumn]>=0) {
            overlap=true;
            numberOverlap++;
          }
        }
        if (goodRow) {
        if (!overlap) {
          // mark all
          for (int j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
            int iColumn = column[j];
            mark[iColumn]=numberSOS;
          }
          numberSOS++;
          numberInSOS += rowLength[iRow];
        }
        // May still be interesting even if overlap
        if (rowLength[iRow]>=5) {
          firstOther--;
          whichRow[firstOther]=iRow;
        }
        }
      }
    }
    if (makeEquality==2&&false) {
      if(numberOverlap||numberIntegers>numberInSOS+1) {
      // try just ones with costs
      CoinFillN(mark,numberColumns,-1);
      numberOverlap=0;
      numberInSOS=0;
      bool allCostsInSOS=true;
      const double *objective = originalModel_->getObjCoefficients() ;
      for (iRow=0;iRow<numberRows;iRow++) {
        if (rowUpper[iRow]==1.0&&rowLength[iRow]>=5) {
          bool goodRow=true;
          bool overlap=false;
          int nObj=0;
          for (int j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
            int iColumn = column[j];
            if (elementByRow[j]!=1.0||!originalModel_->isInteger(iColumn)||lower[iColumn]) {
            goodRow=false;
            }
            if (objective[iColumn])
            nObj++;
            if (mark[iColumn]>=0) {
            overlap=true;
            numberOverlap++;
            }
          }
          if (nObj&&nObj>=rowLength[iRow]-1) {
            if (goodRow) {
            if (!overlap) {
              // mark all
              for (int j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
                int iColumn = column[j];
                mark[iColumn]=numberSOS;
              }
              numberSOS++;
              numberInSOS += rowLength[iRow];
            }
            } else {
            // no good
            allCostsInSOS=false;
            }
          }
        }
      }
      if (numberInSOS&&allCostsInSOS) {
        int nGoodObj=0;
        int nBadObj=0;
        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
          if (objective[iColumn]) {
            if (mark[iColumn]>=0)
            nGoodObj++;
            else
            nBadObj++;
          }
        }
        if (nBadObj*10<nGoodObj) {
          justOnesWithObj=true;
          makeEquality=3;
#ifdef CLP_INVESTIGATE
          printf("trying SOS as all costs there\n");
#endif
        }
      }
      }
    }
    if (firstOther<numberRows&&makeEquality==4) {
      CoinPackedMatrix * matrixByColumn = const_cast<CoinPackedMatrix *>(startModel_->getMatrixByCol());
      // Column copy
      const int * row = matrixByColumn->getIndices();
      const CoinBigIndex * columnStart = matrixByColumn->getVectorStarts();
      const int * columnLength = matrixByColumn->getVectorLengths(); 
      double * columnElements = matrixByColumn->getMutableElements();
      int * rowCount = new int[numberRows];
      memset(rowCount,0,numberRows*sizeof(int));
      double * rowValue = new double [numberRows];
      int numberY=0;
      int numberElements=0;
      int numberSOS=0;
      for (int kRow=firstOther;kRow<numberRows;kRow++) {
      int iRow=whichRow[kRow];
      int n=0;
      int j;
      for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
        int iColumn = column[j];
        for (int k=columnStart[iColumn];k<columnStart[iColumn]+columnLength[iColumn];
             k++) {
          int jRow = row[k];
          double value = columnElements[k];
          if (jRow!=iRow) {
            if (rowCount[jRow]>0) {
            if (value!=rowValue[jRow])
              rowCount[jRow]=-1; // no good
            else
              rowCount[jRow]++;
            } else if (!rowCount[jRow]) {
            whichRow[n++]=jRow;
            rowCount[jRow]=1;
            rowValue[jRow]=value;
            }
          }
        }
      }
      int bestRow=-1;
      int bestCount=4;
      for (j=0;j<n;j++) {
        int jRow = whichRow[j];
        int count=rowCount[jRow];
        rowCount[jRow]=0;
        if (count>=5) {
          numberY++;
          numberElements+=count;
        }
        if (count>bestCount) {
          // possible
          bestRow=jRow;
          bestCount=count;
        }
      }
      if (bestRow>=0) {
        numberSOS++;
        numberY++;
        numberElements+=bestCount;
      }
      }
      if (numberY) {
      // Some may be duplicates
      // make sure ordered
      matrixByRow.orderMatrix();
      elementByRow = matrixByRow.getElements();
      column = matrixByRow.getIndices();
      rowStart = matrixByRow.getVectorStarts();
      rowLength = matrixByRow.getVectorLengths();
      CoinBigIndex * newStart = new CoinBigIndex[numberY+1];
      int * newColumn = new int [numberElements];
      double * newValue = new double [numberElements];
      double * hash = new double [numberY];
      double * hashColumn = new double [numberColumns];
      int i;
      for (i=0;i<numberColumns;i++)
        hashColumn[i]=randomGenerator.randomDouble();
      double * valueY = new double [3*numberY];
      int * rowY = new int [3*numberY];
      int * columnY = new int[3*numberY];
      // For new solution
      double * newSolution = new double [numberColumns+numberY];
      memcpy(newSolution,startModel_->getColSolution(),numberColumns*sizeof(double));
      memset(rowCount,0,numberRows*sizeof(int));
      // List of SOS entries to zero out
      CoinBigIndex * where = new CoinBigIndex[numberColumns];
      numberY=0;
      numberElements=0;
      int numberElementsY=0;
      newStart[0]=0;
      for (int kRow=firstOther;kRow<numberRows;kRow++) {
        int iRow=whichRow[kRow];
        int n=0;
        int j;
        int saveNumberY=numberY;
        for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
          int iColumn = column[j];
          for (int k=columnStart[iColumn];k<columnStart[iColumn]+columnLength[iColumn];
             k++) {
            int jRow = row[k];
            double value = columnElements[k];
            if (jRow!=iRow) {
            if (rowCount[jRow]>0) {
              if (value!=rowValue[jRow])
                rowCount[jRow]=-1; // no good
              else
                rowCount[jRow]++;
            } else if (!rowCount[jRow]) {
              whichRow[n++]=jRow;
              rowCount[jRow]=1;
              rowValue[jRow]=value;
              assert (value);
            }
            }
          }
        }
        for (i=0;i<n;i++) {
          // Sort so fewest first
          std::sort(whichRow,whichRow+n);
          int jRow = whichRow[i];
          int count=rowCount[jRow];
          rowCount[jRow]=0;
          if (count>=5) {
            //assert (count<rowLength[jRow]); // not error - just need to think
            // mark so not looked at again
            rowCount[jRow]=-count;
            // form new row
            double value=0.0;
            double hashValue=0.0;
            int nInSOS=0;
            double valueOfY=0.0;
            for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
            int iColumn = column[j];
            for (int k=columnStart[iColumn];k<columnStart[iColumn]+columnLength[iColumn];
                 k++) {
              if (row[k]==jRow) {
                value = columnElements[k];
                newColumn[numberElements]=iColumn;
                newValue[numberElements++]=1.0;
                hashValue += hashColumn[iColumn];
                columnElements[k]=0.0;
                valueOfY += newSolution[iColumn];
              } else if (row[k]==iRow) {
                if (columnElements[k])
                  where[nInSOS++]=k;
              }
            }
            }
            // See if already exists
            int n=numberElements-newStart[numberY];
            for (j=0;j<numberY;j++) {
            if (hashValue==hash[j]) {
              // Double check
              int offset=newStart[numberY]-newStart[j];
              if (n==newStart[j+1]-newStart[j]) {
                int k;
                for (k=newStart[j];k<newStart[j]+n;k++) {
                  if (newColumn[k]!=newColumn[k+offset])
                  break;
                }
                if (k==newStart[j+1])
                  break;
              }
            }
            }
            if (j==numberY) {
            // not duplicate
            newSolution[numberY+numberColumns]=valueOfY;
            numberY++;
            newStart[numberY]=numberElements;
            hash[j]=hashValue;
            // Now do -1
            rowY[numberElementsY]=j+numberRows;
            columnY[numberElementsY]=j;
            valueY[numberElementsY++]=-1;
            if (n==nInSOS) {
              // SOS entry
              rowY[numberElementsY]=iRow;
              columnY[numberElementsY]=j;
              valueY[numberElementsY++]=1;
              for (int i=0;i<n;i++) {
                int iEl = where[i];
                columnElements[iEl]=0.0;
              }
            }
            } else {
            // duplicate
            numberElements=newStart[numberY];
            }
            // Now do 
            rowY[numberElementsY]=jRow;
            columnY[numberElementsY]=j;
            valueY[numberElementsY++]=value;
          }
        }
        if (numberY>saveNumberY)
          rowCount[iRow]=-1000;
      }
      delete [] hash;
      delete [] hashColumn;
      matrixByColumn->cleanMatrix();
      // Now add rows
      double * rhs = new double[numberY];
      memset(rhs,0,numberY*sizeof(double));
      startModel_->addRows(numberY,newStart,newColumn,newValue,rhs,rhs);
      delete [] rhs;
      delete [] newStart;
      delete [] newColumn;
      delete [] newValue;
      delete [] where;
      // Redo matrix
      CoinPackedMatrix add(true,rowY,columnY,valueY,numberElementsY);
      delete [] valueY;
      delete [] rowY;
      delete [] columnY;
      const int * row = add.getIndices();
      const CoinBigIndex * columnStart = add.getVectorStarts();
      //const int * columnLength = add.getVectorLengths(); 
      double * columnElements = add.getMutableElements();
      double * lo = new double [numberY];
      double * up = new double [numberY];
      for (i=0;i<numberY;i++) {
        lo[i]=0.0;
        up[i]=1.0;
      }
      startModel_->addCols(numberY,columnStart,row,columnElements,lo,up,NULL);
      delete [] lo;
      delete [] up;
      for (i=0;i<numberY;i++) 
        startModel_->setInteger(i+numberColumns);
      CoinWarmStartBasis* basis =
        dynamic_cast <CoinWarmStartBasis*>(startModel_->getWarmStart()) ;
      if (basis) {
        for (i=0;i<numberY;i++) {
          basis->setArtifStatus(i+numberRows,CoinWarmStartBasis::atLowerBound);
          basis->setStructStatus(i+numberColumns,CoinWarmStartBasis::basic);
        }
        startModel_->setWarmStart(basis);
        delete basis;
      }
      startModel_->setColSolution(newSolution);
      delete [] newSolution;
      if (numberElements<10*CoinMin(numberColumns,100*numberY)) {
        handler_->message(CGL_ADDED_INTEGERS,messages_)
          <<numberY<<numberSOS<<numberElements
          <<CoinMessageEol;
        //startModel_->writeMps("new");
        numberColumns += numberY;
        bool saveTakeHint;
        OsiHintStrength saveStrength;
        startModel_->getHintParam(OsiDoDualInResolve,
                        saveTakeHint,saveStrength);
        startModel_->setHintParam(OsiDoDualInResolve,false,OsiHintTry);
        startModel_->resolve();
        numberIterationsPre_ += startModel_->getIterationCount();
        startModel_->setHintParam(OsiDoDualInResolve,saveTakeHint,saveStrength);
      } else {
        // not such a good idea?
        delete startModel_;
        startModel_=NULL;
      }
      }
      delete [] rowValue;
      delete [] rowCount;
    }
    if (makeEquality==4) {
      makeEquality=0;
#if 1
      // Try and make continuous variables integer
      // make clone
      if (!startModel_)
      startModel_ = originalModel_->clone();
      makeInteger();
#endif
    }
    delete [] whichRow;
    delete [] mark;
    if (numberSOS) {
      if (makeEquality==2) {
      if(numberOverlap||numberIntegers>numberInSOS+1) {
        handler_->message(CGL_PROCESS_SOS2,messages_)
          <<numberSOS<<numberInSOS<<numberIntegers<<numberOverlap
          <<CoinMessageEol;
        makeEquality=0;
      }
      }
    } else {
      // no sos
      makeEquality=0;
    }
  }
  if (startModel_) {
    lower = originalModel_->getColLower();
    upper = originalModel_->getColUpper();
    rowLower = originalModel_->getRowLower();
    rowUpper = originalModel_->getRowUpper();
  }
  // See if all + 1
  bool allPlusOnes=true;
  int nPossible=0;
  int numberMadeEquality=0;
  for (iRow=0;iRow<numberRows;iRow++) {
    int numberP1=0, numberM1=0;
    int numberTotal=0;
    int j;
    double upperValue=rowUpper[iRow];
    double lowerValue=rowLower[iRow];
    bool good=true;
    bool possibleSlack=true;
    bool allPlus=true;
    for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
      int iColumn = column[j];
      double value = elementByRow[j];
      if (upper[iColumn]-lower[iColumn]<1.0e-8) {
        // fixed
        upperValue -= lower[iColumn]*value;
        lowerValue -= lower[iColumn]*value;
        continue;
      } else if (!originalModel_->isBinary(iColumn)) {
        good = false;
      possibleSlack=false;
        //break;
      } else {
      numberTotal++;
      }
      
      if (fabs(value-floor(value+0.5))>1.0e-12)
      possibleSlack=false;;
      if (fabs(value)!=1.0) {
        good=false;
        allPlus=false;
      } else if (value>0.0) {
        which[numberP1++]=iColumn;
      } else {
        numberM1++;
        which[numberColumns-numberM1]=iColumn;
        allPlus=false;
      }
    }
    if (possibleSlack) {
      if(upperValue>1.0e20&&lowerValue>-1.0e12) {
      possibleSlack =  (fabs(lowerValue-floor(lowerValue+0.5))<1.0e-12);
      } else if(lowerValue<-1.0e20&&upperValue<1.0e12) {
      possibleSlack =  (fabs(upperValue-floor(upperValue+0.5))<1.0e-12);
      } else {
      possibleSlack=false;
      }
    }
    if (allPlus)
      nPossible++;
    int iUpper = static_cast<int> (floor(upperValue+1.0e-5));
    int iLower = static_cast<int> (ceil(lowerValue-1.0e-5));
    int state=0;
    if (upperValue<1.0e6) {
      if (iUpper==1-numberM1)
        state=1;
      else if (iUpper==-numberM1)
        state=2;
      else if (iUpper<-numberM1)
        state=3;
      if (fabs((static_cast<double> (iUpper))-upperValue)>1.0e-9)
        state =-1;
    }
    if (!state&&lowerValue>-1.0e6) {
      if (-iLower==1-numberP1)
        state=-1;
      else if (-iLower==-numberP1)
        state=-2;
      else if (-iLower<-numberP1)
        state=-3;
      if (fabs((static_cast<double> (iLower))-lowerValue)>1.0e-9)
        state =-1;
    }
    if (good&&state>0) {
      if (abs(state)==3) {
        // infeasible
        feasible=false;
        break;
      } else if (abs(state)==2) {
        // we can fix all
        numberFixed += numberP1+numberM1;
        int i;
        if (state>0) {
          // fix all +1 at 0, -1 at 1
          for (i=0;i<numberP1;i++)
            originalModel_->setColUpper(which[i],0.0);
          for (i=0;i<numberM1;i++)
            originalModel_->setColLower(which[numberColumns-i-1],1.0);
        } else {
          // fix all +1 at 1, -1 at 0
          for (i=0;i<numberP1;i++)
            originalModel_->setColLower(which[i],1.0);
          for (i=0;i<numberM1;i++)
            originalModel_->setColUpper(which[numberColumns-i-1],0.0);
        }
      } else {
        if (!makeEquality||(makeEquality==-1&&numberM1+numberP1<minimumLength))
          continue;
        if (makeEquality==2||makeEquality==3) {
          if (numberM1||numberP1<minimumLength) 
            continue;
        }
        numberCliques++;
        if (iLower!=iUpper) {
          element[numberSlacks]=state;
          rows[numberSlacks++]=iRow;
        }
        if (state>0) {
          totalP1 += numberP1;
          totalM1 += numberM1;
        } else {
          totalP1 += numberM1;
          totalM1 += numberP1;
        }
      }
    }
    if (possibleSlack&&makeEquality==-2&&(!good||state<=0)) {
      if (numberTotal<minimumLength)
      continue;
      numberMadeEquality++;
      element[numberSlacks] = (upperValue<1.0e10) ? 1.0 : -1.0; 
      rows[numberSlacks++]=iRow+numberRows;
    }
  }
  // allow if some +1's
  allPlusOnes = 10*nPossible>numberRows;
  delete [] which;
  if (!feasible) {
    handler_->message(CGL_INFEASIBLE,messages_)
      <<CoinMessageEol;
    delete [] rows;
    delete [] element;
    return NULL;
  } else {
    if (numberCliques) {
      handler_->message(CGL_CLIQUES,messages_)
        <<numberCliques
        << (static_cast<double>(totalP1+totalM1))/
      (static_cast<double> (numberCliques))
        <<CoinMessageEol;
      //printf("%d of these could be converted to equality constraints\n",
      //     numberSlacks);
    }
    if (numberFixed)
      handler_->message(CGL_FIXED,messages_)
        <<numberFixed
        <<CoinMessageEol;
  }
  if (numberSlacks&&makeEquality&&!justOnesWithObj) {
    handler_->message(CGL_SLACKS,messages_)
      <<numberSlacks
      <<CoinMessageEol;
    // add variables to make equality rows
    // Get new model
    if (!startModel_) {
      assert (!startModel_);
      startModel_ = originalModel_->clone();
    }
    for (int i=0;i<numberSlacks;i++) {
      int iRow = rows[i];
      double value = element[i];
      double lowerValue = 0.0;
      double upperValue = 1.0;
      double objValue  = 0.0;
      if (iRow>=numberRows) {
      // just a slack not a clique
      upperValue=COIN_DBL_MAX;
      iRow -= numberRows;
      }
      CoinPackedVector column(1,&iRow,&value);
      startModel_->addCol(column,lowerValue,upperValue,objValue);
      // set integer
      startModel_->setInteger(numberColumns+i);
      if (value >0)
      startModel_->setRowLower(iRow,rowUpper[iRow]);
      else
      startModel_->setRowUpper(iRow,rowLower[iRow]);
    }
  } else if (!startModel_) {
    // make clone anyway so can tighten bounds
    startModel_ = originalModel_->clone();
  }
  // move objective to integers or to aggregated
  lower = startModel_->getColLower();
  upper = startModel_->getColUpper();
  rowLower = startModel_->getRowLower();
  rowUpper = startModel_->getRowUpper();
  matrixByRow = CoinPackedMatrix(*startModel_->getMatrixByRow());
  elementByRow = matrixByRow.getElements();
  column = matrixByRow.getIndices();
  rowStart = matrixByRow.getVectorStarts();
  rowLength = matrixByRow.getVectorLengths();
  char * marked = new char [numberColumns];
  memset(marked,0,numberColumns);
  numberRows=startModel_->getNumRows();
  numberColumns=startModel_->getNumCols();
  //CoinPackedMatrix * matrixByColumn = const_cast<CoinPackedMatrix *>(startModel_->getMatrixByCol());
  // Column copy
  //const int * row = matrixByColumn->getIndices();
  //const CoinBigIndex * columnStart = matrixByColumn->getVectorStarts();
  //const int * columnLength = startModel_->getMatrixByCol()->getVectorLengths(); 
  //const double * columnElements = matrixByColumn->getElements();
  double * obj = CoinCopyOfArray(startModel_->getObjCoefficients(),numberColumns);
  double offset;
  int numberMoved=0;
  startModel_->getDblParam(OsiObjOffset,offset);
  for (iRow=0;iRow<numberRows;iRow++) {
    //int slack = -1;
    int nPlus=0;
    int nMinus=0;
    int iPlus=-1;
    int iMinus=-1;
    double valuePlus=0;
    double valueMinus=0;
    //bool allInteger=true;
    double rhs = rowLower[iRow];
    if (rhs!=rowUpper[iRow])
      continue;
    //int multiple=0;
    //int iSlack=-1;
    int numberContinuous=0;
    for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
      int iColumn = column[j];
      double value = elementByRow[j];
      if (upper[iColumn]>lower[iColumn]) {
      if (startModel_->isInteger(iColumn)) {
#if 0
        if (columnLength[iColumn]==1) {
          if (value==1.0) {
          }
        }
        if (value!=floor(value+0.5))
          allInteger=false;
        if (allInteger&&fabs(value)<1.0e8) {
          if (!multiple)
            multiple = static_cast<int> (fabs(value));
          else if (multiple>0)
            multiple = gcd(multiple,static_cast<int> (fabs(value)));
        } else {
          allInteger=false;
        }
#endif
      } else {
        numberContinuous++;
      }
      if (value>0.0) {
        if (nPlus>0&&value!=valuePlus) {
          nPlus = - numberColumns;
        } else if (!nPlus) {
          nPlus=1;
          iPlus=iColumn;
          valuePlus=value;
        } else {
          nPlus++;
        }
      } else {
        if (nMinus>0&&value!=valueMinus) {
          nMinus = - numberColumns;
        } else if (!nMinus) {
          nMinus=1;
          iMinus=iColumn;
          valueMinus=value;
        } else {
          nMinus++;
        }
      }
      } else {
      rhs -= lower[iColumn]*value;
      }
    }
    if (((nPlus==1&&startModel_->isInteger(iPlus)&&nMinus>0)||
       (nMinus==1&&startModel_->isInteger(iMinus)&&nPlus>0))&&numberContinuous&&true) {
      int jColumn;
      double multiplier;
      if (nPlus==1) {
      jColumn = iPlus;
      multiplier = fabs(valuePlus/valueMinus);
      rhs /= -valueMinus;
      } else {
      jColumn = iMinus;
      multiplier = fabs(valueMinus/valuePlus);
      rhs /= valuePlus;
      }
      double smallestPos=COIN_DBL_MAX;
      double smallestNeg=-COIN_DBL_MAX;
      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
      int iColumn = column[j];
      if (iColumn!=jColumn) {
        double objValue = obj[iColumn];
        if (upper[iColumn]>lower[iColumn]) {
          if (objValue>=0.0)
            smallestPos=CoinMin(smallestPos,objValue);
          else
            smallestNeg=CoinMax(smallestNeg,objValue);
        }
      }
      }
      if (smallestPos>0.0) {
      double move=0.0;
      if(smallestNeg==-COIN_DBL_MAX)
        move=smallestPos;
      else if (smallestPos==COIN_DBL_MAX)
        move=smallestNeg;
      if (move) {
        // can move objective
        numberMoved++;
#ifdef COIN_DEVELOP
        if (rhs)
          printf("ZZZ on col %d move %g offset %g\n",
               jColumn,move,move*rhs);
#endif
        offset += move*rhs;
        for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
          int iColumn = column[j];
          if (iColumn!=jColumn) {
            if (upper[iColumn]>lower[iColumn]) {
            obj[iColumn] -= move;
            }
          } else {
            obj[jColumn] += move*multiplier;
          }
        }
      }
      }
    }
  }
#ifdef COIN_DEVELOP
  if (numberMoved)
    printf("ZZZ %d costs moved\n",numberMoved);
#endif
  startModel_->setDblParam(OsiObjOffset,offset);
  startModel_->setObjective(obj);
  delete [] obj;
  delete [] marked;
  delete [] rows;
  delete [] element;
  if (makeIntegers) {
    makeIntegers2(startModel_,makeIntegers);
  }
  int infeas=0;
  // Do we want initial presolve
  if (doInitialPresolve) {
    assert (doInitialPresolve==1);
    OsiSolverInterface * presolvedModel;
    OsiSolverInterface * oldModel = startModel_;
    OsiPresolve * pinfo = new OsiPresolve();
    int presolveActions=0;
    // Allow dual stuff on integers
    // Allow stuff which may not unroll cleanly
    presolveActions=1+16;
    if ((tuning&32)!=0)
      presolveActions |= 32;
    // Do not allow all +1 to be tampered with
    //if (allPlusOnes)
    //presolveActions |= 2;
    // allow transfer of costs
    // presolveActions |= 4;
    // If trying for SOS don't allow some transfers
    if (makeEquality==2||makeEquality==3)
      presolveActions |= 8;
    pinfo->setPresolveActions(presolveActions);
    if (prohibited_)
      assert (numberProhibited_==oldModel->getNumCols());
    int saveLogLevel = oldModel->messageHandler()->logLevel();
    if (saveLogLevel==1)
      oldModel->messageHandler()->setLogLevel(0);
    std::string solverName;
    oldModel->getStrParam(OsiSolverName,solverName);
    // Extend if you want other solvers to keep solution
    bool keepSolution=solverName=="clp";
    presolvedModel = pinfo->presolvedModel(*oldModel,1.0e-7,true,5,prohibited_,keepSolution,rowType_);
    oldModel->messageHandler()->setLogLevel(saveLogLevel);
    if (presolvedModel) {
      presolvedModel->messageHandler()->setLogLevel(saveLogLevel);
      // update prohibited and rowType
      update(pinfo,presolvedModel);
      if (!presolvedModel->getNumRows()) {
      doInitialPresolve=0;
      delete presolvedModel;
      delete pinfo;
      } else {
      model_[0]=presolvedModel;
      presolve_[0]=pinfo;
      modifiedModel_[0]=presolvedModel->clone();
      }
    } else {
      infeas=1;
      doInitialPresolve=0;
      delete presolvedModel;
      delete pinfo;
    }
  }
  // tighten bounds
/*

  Virtuous solvers may require a refresh via initialSolve if this
  call is ever changed to give a nonzero value to the (default) second
  parameter. Previous actions may have made significant changes to the
  constraint system. Safe as long as tightenPrimalBounds doesn't ask for
  the current solution.
*/
  if (!infeas) {
    //startModel_->writeMps("before-tighten");
    infeas = tightenPrimalBounds(*startModel_);
    //startModel_->writeMps("after-tighten");
  }
  if (infeas) {
    handler_->message(CGL_INFEASIBLE,messages_)
      <<CoinMessageEol;
    return NULL;
  }
  OsiSolverInterface * returnModel=NULL;
  int numberChanges;
  {
    // Give a hint to do dual
    bool saveTakeHint;
    OsiHintStrength saveStrength;
    startModel_->getHintParam(OsiDoDualInInitial,
                        saveTakeHint,saveStrength);
    startModel_->setHintParam(OsiDoDualInInitial,true,OsiHintTry);
    startModel_->initialSolve();
    numberIterationsPre_ += startModel_->getIterationCount();
    // double check
    if (!startModel_->isProvenOptimal()) {
      if (!startModel_->isProvenDualInfeasible()) {
      // Do presolves
      bool saveHint;
      OsiHintStrength saveStrength;
      startModel_->getHintParam(OsiDoPresolveInInitial,saveHint,saveStrength);
      startModel_->setHintParam(OsiDoPresolveInInitial,true,OsiHintTry);
      startModel_->setHintParam(OsiDoDualInInitial,false,OsiHintTry);
      startModel_->initialSolve();
      numberIterationsPre_ += startModel_->getIterationCount();
      if (!startModel_->isProvenDualInfeasible()) {
        CoinWarmStart * empty = startModel_->getEmptyWarmStart();
        startModel_->setWarmStart(empty);
        delete empty;
        startModel_->setHintParam(OsiDoDualInInitial,true,OsiHintTry);
        startModel_->initialSolve();
        numberIterationsPre_ += startModel_->getIterationCount();
      }
      startModel_->setHintParam(OsiDoPresolveInInitial,saveHint,saveStrength);
      }
    }
    startModel_->setHintParam(OsiDoDualInInitial,saveTakeHint,saveStrength);
  }
  if (!startModel_->isProvenOptimal()) {
    if (!startModel_->isProvenDualInfeasible()) {
      handler_->message(CGL_INFEASIBLE,messages_)<< CoinMessageEol ;
#ifdef COIN_DEVELOP
      startModel_->writeMps("infeas");
#endif
    } else {
      handler_->message(CGL_UNBOUNDED,messages_)<< CoinMessageEol ;
    }
    return NULL;
  }
  reducedCostFix(*startModel_);
  if (!numberSolvers_) {
    // just fix
    OsiSolverInterface * newModel = modified(startModel_,false,numberChanges,0,numberModifiedPasses);
    if (startModel_!=originalModel_)
      delete startModel_;
    startModel_=newModel;
    returnModel=startModel_;
  } else {
    OsiSolverInterface * presolvedModel;
    OsiSolverInterface * oldModel = startModel_;
    if (doInitialPresolve)
      oldModel = modifiedModel_[0];
    //CglDuplicateRow dupCuts(oldModel);
    //dupCuts.setLogLevel(1);
    // If +1 try duplicate rows
    if (allPlusOnes) {
#if 0 
      // put at beginning
      CglCutGenerator ** temp = generator_;
      generator_ = new CglCutGenerator * [numberCutGenerators_+1];
      memcpy(generator_+1,temp,numberCutGenerators_*sizeof(CglCutGenerator *));
      delete[] temp ;
      numberCutGenerators_++;
      generator_[0]=new CglDuplicateRow(oldModel);
#else
      CglDuplicateRow dupCuts(oldModel);
      addCutGenerator(&dupCuts);
#endif
    }
    for (int iPass=doInitialPresolve;iPass<numberSolvers_;iPass++) {
      // Look at Vubs
      {
        const double * columnLower = oldModel->getColLower();
        const double * columnUpper = oldModel->getColUpper();
        const CoinPackedMatrix * rowCopy = oldModel->getMatrixByRow();
        const int * column = rowCopy->getIndices();
        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
        const int * rowLength = rowCopy->getVectorLengths(); 
        const double * rowElements = rowCopy->getElements();
        const CoinPackedMatrix * columnCopy = oldModel->getMatrixByCol();
        //const int * row = columnCopy->getIndices();
        //const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
        const int * columnLength = columnCopy->getVectorLengths(); 
        //const double * columnElements = columnCopy->getElements();
        const double * rowLower = oldModel->getRowLower();
        const double * rowUpper = oldModel->getRowUpper();
        const double * objective = oldModel->getObjCoefficients();
        double direction = oldModel->getObjSense();
        int numberRows = oldModel->getNumRows();
        for (int iRow=0;iRow<numberRows;iRow++) {
          if (rowLength[iRow]==2&&(rowLower[iRow]<-1.0e20||rowUpper[iRow]>1.0e20)) {
            CoinBigIndex start = rowStart[iRow];
            int iColumn1 = column[start];
            int iColumn2 = column[start+1];
            double value1 = rowElements[start];
            double value2 = rowElements[start+1];
            double upper;
            if (rowLower[iRow]<-1.0e20) {
              if (rowUpper[iRow]<1.0e20)
                upper = rowUpper[iRow];
              else
                continue; // free row
            } else {
              upper = - rowLower[iRow];
              value1=-value1;
              value2=-value2;
            }
            //for now just singletons
            bool integer1 = oldModel->isInteger(iColumn1);
            bool integer2 = oldModel->isInteger(iColumn2);
            int debug=0;
            if (columnLength[iColumn1]==1) {
              if (integer1) {
                debug=0;// no good
              } else if (integer2) {
                // possible
                debug=1;
              }
            } else if (columnLength[iColumn2]==1) {
              if (integer2) {
                debug=-1; // print and skip
              } else if (integer1) {
                // possible
                debug=1;
                double valueD = value1;
                value1 = value2;
                value2 = valueD;
                int valueI = iColumn1;
                iColumn1 = iColumn2;
                iColumn2 = valueI;
                bool valueB = integer1;
                integer1 = integer2;
                integer2 = valueB;
              }
            }
            if (debug&&0) {
              printf("%d %d elements%selement %g and %d %d elements%selement %g <= %g\n",
                     iColumn1,columnLength[iColumn1],integer1 ? " (integer) " : " ",value1,
                     iColumn2,columnLength[iColumn2],integer2 ? " (integer) " : " ",value2,
                     upper);
            }
            if (debug>0) {
              if (value1>0.0&&objective[iColumn1]*direction<0.0) {
                // will push as high as possible so make ==
                // highest effective rhs
                if (value2>0) 
                  upper -= value2 * columnLower[iColumn2];
                else
                  upper -= value2 * columnUpper[iColumn2];
                if (columnUpper[iColumn1]>1.0e20||
                    columnUpper[iColumn1]*value1>=upper) {
                  //printf("looks possible\n");
                  // make equality
                  if (rowLower[iRow]<-1.0e20) 
                    oldModel->setRowLower(iRow,rowUpper[iRow]);
                  else
                    oldModel->setRowUpper(iRow,rowLower[iRow]);
                } else {
                  // may be able to make integer
                  // may just be better to use to see objective integral
                  if (upper==floor(upper)&&value2==floor(value2)&&
                      value1==floor(value1)&&objective[iColumn1]==floor(objective[iColumn1]))
                    oldModel->setInteger(iColumn1);
                  //printf("odd3\n");
                }
              } else if (value1<0.0&&objective[iColumn1]*direction>0.0) {
                //printf("odd4\n");
              } else {
                //printf("odd2\n");
              }
            } else if (debug<0) {
              //printf("odd1\n");
            }
          }
        }
      }
      OsiPresolve * pinfo = new OsiPresolve();
      int presolveActions=0;
      // Allow dual stuff on integers
      // Allow stuff which may not unroll cleanly
      presolveActions=1+16;
      // Do not allow all +1 to be tampered with
      //if (allPlusOnes)
      //presolveActions |= 2;
      // allow transfer of costs
      // presolveActions |= 4;
      // If trying for SOS don't allow some transfers
      if (makeEquality==2||makeEquality==3)
        presolveActions |= 8;
      pinfo->setPresolveActions(presolveActions);
      if (prohibited_)
        assert (numberProhibited_==oldModel->getNumCols());
/*
  VIRTUOUS but possible bad for performance 
  
  At this point, the solution is most likely stale: we may have added cuts as
  we left the previous call to modified(), or we may have changed row bounds
  in VUB analysis just above. Continuous presolve doesn't need a solution
  unless we want it to transform the current solution to match the presolved
  model.
*/
      int saveLogLevel = oldModel->messageHandler()->logLevel();
      if (saveLogLevel==1)
      oldModel->messageHandler()->setLogLevel(0);
      std::string solverName;
      oldModel->getStrParam(OsiSolverName,solverName);
      // Extend if you want other solvers to keep solution
      bool keepSolution=solverName=="clp";
      presolvedModel = pinfo->presolvedModel(*oldModel,1.0e-7,true,5,
                                   prohibited_,keepSolution,rowType_);
      oldModel->messageHandler()->setLogLevel(saveLogLevel);
      if (!presolvedModel) {
        returnModel=NULL;
      delete pinfo;
        break;
      }
      presolvedModel->messageHandler()->setLogLevel(saveLogLevel);
      // update prohibited and rowType
      update(pinfo,presolvedModel);
      //char name[20];
      //sprintf(name,"prex%2.2d.mps",iPass);
      //presolvedModel->writeMpsNative(name, NULL, NULL,0,1,0);
      model_[iPass]=presolvedModel;
      presolve_[iPass]=pinfo;
      if (!presolvedModel->getNumRows()) {
        // was returnModel=oldModel;
        returnModel=presolvedModel;
        numberSolvers_=iPass+1;
        break; // model totally solved
      }
      bool constraints = iPass<numberPasses-1;
      // Give a hint to do primal
      bool saveTakeHint;
      OsiHintStrength saveStrength;
      presolvedModel->getHintParam(OsiDoDualInInitial,
                                   saveTakeHint,saveStrength);
      //if (iPass)
      presolvedModel->setHintParam(OsiDoDualInInitial,false,OsiHintTry);
      presolvedModel->initialSolve();
      numberIterationsPre_ += presolvedModel->getIterationCount();
      presolvedModel->setHintParam(OsiDoDualInInitial,saveTakeHint,saveStrength);
      if (!presolvedModel->isProvenOptimal()) {
        returnModel=NULL;
        break;
      }
      // maybe we can fix some
      int numberFixed = 
      reducedCostFix(*presolvedModel);
#ifdef COIN_DEVELOP
      if (numberFixed)
      printf("%d variables fixed on reduced cost\n",numberFixed);
#endif
      OsiSolverInterface * newModel = modified(presolvedModel,constraints,numberChanges,iPass-doInitialPresolve,numberModifiedPasses);
      returnModel=newModel;
      if (!newModel) {
        break;
      }
      modifiedModel_[iPass]=newModel;
      oldModel=newModel;
      //sprintf(name,"pre%2.2d.mps",iPass);
      //newModel->writeMpsNative(name, NULL, NULL,0,1,0);
      if (!numberChanges&&!numberFixed) {
#ifdef COIN_DEVELOP
      printf("exiting after pass %d of %d\n",iPass,numberSolvers_);
#endif
        numberSolvers_=iPass+1;
        break;
      }
    }
  }
  if (returnModel) {
    if (returnModel->getNumRows()) {
      // tighten bounds
      int infeas = tightenPrimalBounds(*returnModel);
      if (infeas) {
        delete returnModel;
      for (int iPass=0;iPass<numberSolvers_;iPass++) {
        if (returnModel==modifiedModel_[iPass])
          modifiedModel_[iPass]=NULL;
      }
        if (returnModel==startModel_&&startModel_!=originalModel_)
          startModel_=NULL;
        returnModel=NULL;
      }
    }
  } else {
    handler_->message(CGL_INFEASIBLE,messages_)
      <<CoinMessageEol;
  }
  int numberIntegers=0;
  if (returnModel) {
    int iColumn;
    int numberColumns = returnModel->getNumCols();
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      if (returnModel->isInteger(iColumn))
        numberIntegers++;
    }
  }
  if ((makeEquality==2||makeEquality==3)&&numberCliques&&returnModel) {
    int iRow, iColumn;
    int numberColumns = returnModel->getNumCols();
    int numberRows = returnModel->getNumRows();
    const double * objective = returnModel->getObjCoefficients();
    // get row copy
    const CoinPackedMatrix * matrix = returnModel->getMatrixByRow();
    const double * element = matrix->getElements();
    const int * column = matrix->getIndices();
    const CoinBigIndex * rowStart = matrix->getVectorStarts();
    const int * rowLength = matrix->getVectorLengths();
    const double * rowLower = returnModel->getRowLower();
    const double * rowUpper = returnModel->getRowUpper();
    const double * columnLower = returnModel->getColLower();
    
    // Look for possible SOS
    int numberSOS=0;
    int * mark = new int[numberColumns];
    int * sosRow = new int [numberRows];
    CoinZeroN(sosRow,numberRows);
    CoinFillN(mark,numberColumns,-1);
    int numberOverlap=0;
    int numberInSOS=0;
    for (iRow=0;iRow<numberRows;iRow++) {
      if (rowLower[iRow]==1.0&&rowUpper[iRow]==1.0) {
        if ((rowLength[iRow]<5&&!justOnesWithObj)||(rowLength[iRow]<20&&allToGub))
          continue;
        bool goodRow=true;
      int nObj=0;
        for (int j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
          iColumn = column[j];
          if (element[j]!=1.0||!returnModel->isInteger(iColumn)||columnLower[iColumn]) {
            goodRow=false;
            break;
          }
          if (mark[iColumn]>=0&&!allToGub) {
            goodRow=false;
            numberOverlap++;
          }
        if (objective[iColumn])
          nObj++;
        }
      if (goodRow&&justOnesWithObj) {
        if (!nObj||nObj<rowLength[iRow]-1) 
          goodRow=false;
      }
        if (goodRow) {
          // mark all
          for (int j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
            int iColumn = column[j];
            mark[iColumn]=numberSOS;
          }
          sosRow[numberSOS++]=iRow;
          numberInSOS += rowLength[iRow];
        }
      }
    }
    if (numberSOS) {
      if (makeEquality==2&&(numberOverlap||numberIntegers>numberInSOS+1)) {
        handler_->message(CGL_PROCESS_SOS2,messages_)
          <<numberSOS<<numberInSOS<<numberIntegers<<numberOverlap
          <<CoinMessageEol;
      } else {
        handler_->message(CGL_PROCESS_SOS1,messages_)
          <<numberSOS<<numberInSOS
          <<CoinMessageEol;
        numberSOS_=numberSOS;
        typeSOS_ = new int[numberSOS_];
        startSOS_ = new int[numberSOS_+1];
        whichSOS_ = new int[numberInSOS];
        weightSOS_ = new double[numberInSOS];
        numberInSOS=0;
        startSOS_[0]=0;
        const CoinPackedMatrix * columnCopy = returnModel->getMatrixByCol();
        const int * row = columnCopy->getIndices();
        const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
        const int * columnLength = columnCopy->getVectorLengths(); 
        const double * columnElements = columnCopy->getElements();
        const double * objective = returnModel->getObjCoefficients();
        int * numberInRow = new int [numberRows];
        double * sort = new double[numberColumns];
        int * which = new int[numberColumns];
        for (int iSOS =0;iSOS<numberSOS_;iSOS++) {
          int n=0;
          int numberObj=0;
          CoinZeroN(numberInRow,numberRows);
        int kRow = sosRow[iSOS];
          for (int j=rowStart[kRow];j<rowStart[kRow]+rowLength[kRow];j++) {
            int iColumn = column[j];
          whichSOS_[numberInSOS]=iColumn;
          weightSOS_[numberInSOS]=n;
          numberInSOS++;
          n++;
          if (objective[iColumn])
            numberObj++;
          for (CoinBigIndex j=columnStart[iColumn];
             j<columnStart[iColumn]+columnLength[iColumn];j++) {
            int iRow = row[j];
            if (!sosRow[iRow])
            numberInRow[iRow]++;
          }
        }
          // See if any rows look good
          int bestRow=-1;
          int numberDifferent=1;
          int start = startSOS_[iSOS];
          for (int iRow=0;iRow<numberRows;iRow++) {
            if (numberInRow[iRow]>=n-1) {
              // See how many different
              int i;
              for ( i=0;i<n;i++) {
                int iColumn = whichSOS_[i+start];
                sort[i]=0.0;
                which[i]=iColumn;
                for (CoinBigIndex j=columnStart[iColumn];
                     j<columnStart[iColumn]+columnLength[iColumn];j++) {
                  int jRow = row[j];
                  if (jRow==iRow) {
                    sort[i]=columnElements[j];
                    break;
                  }
                }
              }
              // sort
              CoinSort_2(sort,sort+n,which);
              double last = sort[0];
              int nDiff=1;
              for ( i=1;i<n;i++) {
                if (sort[i]>last+CoinMax(fabs(last)*1.0e-8,1.0e-5)) {
                  nDiff++;
                }
                last = sort[i];
              }
              if (nDiff>numberDifferent) {
                numberDifferent = nDiff;
                bestRow=iRow;
              }
            }
          }
          if (numberObj>=n-1||bestRow<0) {
            int i;
            for ( i=0;i<n;i++) {
              int iColumn = whichSOS_[i+start];
              sort[i]=objective[iColumn];
              which[i]=iColumn;
            }
            // sort
            CoinSort_2(sort,sort+n,which);
            double last = sort[0];
            int nDiff=1;
            for ( i=1;i<n;i++) {
              if (sort[i]>last+CoinMax(fabs(last)*1.0e-8,1.0e-5)) {
                nDiff++;
              }
              last = sort[i];
            }
            if (nDiff>numberDifferent) {
              numberDifferent = nDiff;
              bestRow=numberRows;
            }
          }
          if (bestRow>=0) {
            // if not objective - recreate
            if (bestRow<numberRows) {
              int i;
              for ( i=0;i<n;i++) {
                int iColumn = whichSOS_[i+start];
                sort[i]=0.0;
                which[i]=iColumn;
                for (CoinBigIndex j=columnStart[iColumn];
                     j<columnStart[iColumn]+columnLength[iColumn];j++) {
                  int jRow = row[j];
                  if (jRow==bestRow) {
                    sort[i]=columnElements[j];
                    break;
                  }
                }
              }
              // sort
              CoinSort_2(sort,sort+n,which);
            }
            // make sure gaps OK
            double last = sort[0];
            for (int i=1;i<n;i++) {
              double next = last+CoinMax(fabs(last)*1.0e-8,1.0e-5);
              sort[i]=CoinMax(sort[i],next);
              last = sort[i];
            }
            //CoinCopyN(sort,n,weightSOS_+start);
            //CoinCopyN(which,n,whichSOS_+start);
          }
          typeSOS_[iSOS]=1;
          startSOS_[iSOS+1]=numberInSOS;
        }
        delete [] numberInRow;
        delete [] sort;
        delete [] which;
      }
    }
    delete [] mark;
    delete [] sosRow;
  }
  if (returnModel) {
    if (makeIntegers) 
      makeIntegers2(returnModel,makeIntegers);
    handler_->message(CGL_PROCESS_STATS2,messages_)
      <<returnModel->getNumRows()<<returnModel->getNumCols()
      <<numberIntegers<<returnModel->getNumElements()
      <<CoinMessageEol;
    // If can make some cuts then do so
    if (rowType_) {
      int numberRows = returnModel->getNumRows();
      int numberCuts=0;
      for (int i=0;i<numberRows;i++) {
      if (rowType_[i]>0)
        numberCuts++;
      }
      if (numberCuts) {
      CglStored stored;
      
      int * whichRow = new int[numberRows];
      // get row copy
      const CoinPackedMatrix * rowCopy = returnModel->getMatrixByRow();
      const int * column = rowCopy->getIndices();
      const int * rowLength = rowCopy->getVectorLengths();
      const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
      const double * rowLower = returnModel->getRowLower();
      const double * rowUpper = returnModel->getRowUpper();
      const double * element = rowCopy->getElements();
      int iRow,nDelete=0;
      for (iRow=0;iRow<numberRows;iRow++) {
        if (rowType_[iRow]==1) {
          // take out
          whichRow[nDelete++]=iRow;
        }
      }
      for (int jRow=0;jRow<nDelete;jRow++) {
        iRow=whichRow[jRow];
        int start = rowStart[iRow];
        stored.addCut(rowLower[iRow],rowUpper[iRow],rowLength[iRow],
                  column+start,element+start);
      }
      returnModel->deleteRows(nDelete,whichRow);
      delete [] whichRow;
      cuts_ = stored;
      }
    }
  }
#if 0
  if (returnModel) {
    int numberColumns = returnModel->getNumCols();
    int numberRows = returnModel->getNumRows();
    int * del = new int [CoinMax(numberColumns,numberRows)];
    int * original = new int [numberColumns];
    int nDel=0;
    for (int i=0;i<numberColumns;i++) {
      original[i]=i;
      if (returnModel->isInteger(i))
      del[nDel++]=i;
    }
    int nExtra=0;
    if (nDel&&nDel!=numberColumns&&(options_&1)!=0&&false) {
      OsiSolverInterface * yyyy = returnModel->clone();
      int nPass=0;
      while (nDel&&nPass<10) {
      nPass++;
      OsiSolverInterface * xxxx = yyyy->clone();
      int nLeft=0;
      for (int i=0;i<nDel;i++) 
        original[del[i]]=-1;
      for (int i=0;i<numberColumns;i++) {
        int kOrig=original[i];
        if (kOrig>=0)
          original[nLeft++]=kOrig;
      }
      assert (nLeft==numberColumns-nDel);
      xxxx->deleteCols(nDel,del);
      numberColumns = xxxx->getNumCols();
      const CoinPackedMatrix * rowCopy = xxxx->getMatrixByRow();
      numberRows = rowCopy->getNumRows();
      const int * column = rowCopy->getIndices();
      const int * rowLength = rowCopy->getVectorLengths();
      const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
      const double * rowLower = xxxx->getRowLower();
      const double * rowUpper = xxxx->getRowUpper();
      const double * element = rowCopy->getElements();
        const CoinPackedMatrix * columnCopy = xxxx->getMatrixByCol();
        const int * columnLength = columnCopy->getVectorLengths(); 
      nDel=0;
      // Could do gcd stuff on ones with costs
      for (int i=0;i<numberRows;i++) {
        if (!rowLength[i]) {
          del[nDel++]=i;
        } else if (rowLength[i]==1) {
          int k=rowStart[i];
          int iColumn = column[k];
          if (!xxxx->isInteger(iColumn)) {
            double mult =1.0/fabs(element[k]);
            if (rowLower[i]<-1.0e20) {
            double value = rowUpper[i]*mult;
            if (fabs(value-floor(value+0.5))<1.0e-8) {
              del[nDel++]=i;
              if (columnLength[iColumn]==1) {
                xxxx->setInteger(iColumn);
                int kOrig=original[iColumn];
                returnModel->setInteger(kOrig);
              }
            }
            } else if (rowUpper[i]>1.0e20) {
            double value = rowLower[i]*mult;
            if (fabs(value-floor(value+0.5))<1.0e-8) {
              del[nDel++]=i;
              if (columnLength[iColumn]==1) {
                xxxx->setInteger(iColumn);
                int kOrig=original[iColumn];
                returnModel->setInteger(kOrig);
              }
            }
            } else {
            double value = rowUpper[i]*mult;
            if (rowLower[i]==rowUpper[i]&&
                fabs(value-floor(value+0.5))<1.0e-8) {
              del[nDel++]=i;
              xxxx->setInteger(iColumn);
              int kOrig=original[iColumn];
              returnModel->setInteger(kOrig);
            }
            }
          }
        } else {
          // only if all singletons
          bool possible=false;
          if (rowLower[i]<-1.0e20) {
            double value = rowUpper[i];
            if (fabs(value-floor(value+0.5))<1.0e-8) 
            possible=true;
          } else if (rowUpper[i]>1.0e20) {
            double value = rowLower[i];
            if (fabs(value-floor(value+0.5))<1.0e-8) 
            possible=true;
          } else {
            double value = rowUpper[i];
            if (rowLower[i]==rowUpper[i]&&
              fabs(value-floor(value+0.5))<1.0e-8)
            possible=true;
          }
          if (possible) {
            for (CoinBigIndex j=rowStart[i];
               j<rowStart[i]+rowLength[i];j++) {
            int iColumn = column[j];
            if (columnLength[iColumn]!=1||fabs(element[j])!=1.0) {
              possible=false;
              break;
            }
            }
            if (possible) {
            for (CoinBigIndex j=rowStart[i];
                 j<rowStart[i]+rowLength[i];j++) {
              int iColumn = column[j];
              if (!xxxx->isInteger(iColumn)) {
                xxxx->setInteger(iColumn);
                int kOrig=original[iColumn];
                returnModel->setInteger(kOrig);
              }
            }
            del[nDel++]=i;
            }
          }
        }
      }
      if (nDel) {
        xxxx->deleteRows(nDel,del);
      }
      if (nDel!=numberRows) {
        nDel=0;
        for (int i=0;i<numberColumns;i++) {
          if (xxxx->isInteger(i)) {
            del[nDel++]=i;
            nExtra++;
          }
        }
      } 
      delete yyyy;
      yyyy=xxxx->clone();
      }
      numberColumns = yyyy->getNumCols();
      numberRows = yyyy->getNumRows();
      if (!numberColumns||!numberRows) {
      printf("All gone\n");
      int numberColumns = returnModel->getNumCols();
      for (int i=0;i<numberColumns;i++)
        assert(returnModel->isInteger(i));
      }
      // Would need to check if original bounds integer
      //yyyy->writeMps("noints");
      delete yyyy;
      printf("Creating simplified model with %d rows and %d columns - %d extra integers\n",
           numberRows,numberColumns,nExtra);
    }
    delete [] del;
    delete [] original;
    //exit(2);
  }
#endif
  return returnModel;
}


Generated by  Doxygen 1.6.0   Back to index