decompose.cpp

// Copyright (C) 2003, International Business Machines
// Corporation and others.  All Rights Reserved.

#include "ClpSimplex.hpp"
#include "CoinMpsIO.hpp"
#include <iomanip>

int main (int argc, const char *argv[])
{
  ClpSimplex  model;
  int status;
  // Keep names
  if (argc<2) {
    //status=model.readMps("/home/forrest/data/ken-18.mps.gz",true);
    status=model.readMps("small.mps",true);
  } else {
    status=model.readMps(argv[1],true);
  }
  if (status)
    exit(10);
  /*
    This driver does a simple Dantzig Wolfe decomposition
  */
  // Get master rows in some mysterious way
  int numberRows = model.numberRows();
  int * rowBlock = new int[numberRows];
  int iRow;
  for (iRow=0;iRow<numberRows;iRow++)
    rowBlock[iRow]=-2;
  // these are master rows
  if (numberRows==105127) {
    // ken-18
    for (iRow=104976;iRow<numberRows;iRow++)
      rowBlock[iRow]=-1;
  } else if (numberRows==810) {
    for (iRow=81;iRow<84;iRow++)
      rowBlock[iRow]=-1;
  } else if (numberRows==5418) {
    for (iRow=564;iRow<603;iRow++)
      rowBlock[iRow]=-1;
  } else if (numberRows==1503) {
    // degen3
    for (iRow=0;iRow<561;iRow++)
      rowBlock[iRow]=-1;
  }
  CoinPackedMatrix * matrix = model.matrix();
  // get row copy
  CoinPackedMatrix rowCopy = *matrix;
  rowCopy.reverseOrdering();
  const int * row = matrix->getIndices();
  const int * columnLength = matrix->getVectorLengths();
  const CoinBigIndex * columnStart = matrix->getVectorStarts();
  //const double * elementByColumn = matrix->getElements();
  const int * column = rowCopy.getIndices();
  const int * rowLength = rowCopy.getVectorLengths();
  const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
  //const double * elementByRow = rowCopy.getElements();
  int numberBlocks = 0;
  int * stack = new int [numberRows];
  // to say if column looked at
  int numberColumns = model.numberColumns();
  int * columnBlock = new int[numberColumns];
  int iColumn;
  for (iColumn=0;iColumn<numberColumns;iColumn++)
    columnBlock[iColumn]=-2;
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    int kstart = columnStart[iColumn];
    int kend = columnStart[iColumn]+columnLength[iColumn];
    if (columnBlock[iColumn]==-2) {
      // column not allocated
      int j;
      int nstack=0;
      for (j=kstart;j<kend;j++) {
        int iRow= row[j];
        if (rowBlock[iRow]!=-1) {
          assert(rowBlock[iRow]==-2);
          rowBlock[iRow]=numberBlocks; // mark
          stack[nstack++] = iRow;
        }
      }
      if (nstack) {
        // new block - put all connected in
        numberBlocks++;
        columnBlock[iColumn]=numberBlocks-1;
        while (nstack) {
          int iRow = stack[--nstack];
          int k;
          for (k=rowStart[iRow];k<rowStart[iRow]+rowLength[iRow];k++) {
            int iColumn = column[k];
            int kkstart = columnStart[iColumn];
            int kkend = kkstart + columnLength[iColumn];
            if (columnBlock[iColumn]==-2) {
              columnBlock[iColumn]=numberBlocks-1; // mark
              // column not allocated
              int jj;
              for (jj=kkstart;jj<kkend;jj++) {
                int jRow= row[jj];
                if (rowBlock[jRow]==-2) {
                  rowBlock[jRow]=numberBlocks-1;
                  stack[nstack++]=jRow;
                }
              }
            } else {
              assert (columnBlock[iColumn]==numberBlocks-1);
            }
          }
        }
      } else {
        // Only in master
        columnBlock[iColumn]=-1;
      }
    }
  }
  printf("%d blocks found\n",numberBlocks);
  delete [] stack;
  // make up problems
  CoinPackedMatrix * top = new CoinPackedMatrix [numberBlocks];
  ClpSimplex * sub = new ClpSimplex [numberBlocks];
  ClpSimplex master;
  // Create all sub problems
  // Could do much faster - but do that later
  int * whichRow = new int [numberRows];
  int * whichColumn = new int [numberColumns];
  // get top matrix
  CoinPackedMatrix topMatrix = *model.matrix();
  int numberRow2,numberColumn2;
  numberRow2=0;
  for (iRow=0;iRow<numberRows;iRow++)
    if (rowBlock[iRow]>=0)
      whichRow[numberRow2++]=iRow;
  topMatrix.deleteRows(numberRow2,whichRow);
  int iBlock;
  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
    numberRow2=0;
    numberColumn2=0;
    for (iRow=0;iRow<numberRows;iRow++)
      if (iBlock==rowBlock[iRow])
        whichRow[numberRow2++]=iRow;
    for (iColumn=0;iColumn<numberColumns;iColumn++)
      if (iBlock==columnBlock[iColumn])
        whichColumn[numberColumn2++]=iColumn;
    sub[iBlock]=ClpSimplex(&model,numberRow2,whichRow,
                           numberColumn2,whichColumn);
#if 1
    // temp
    double * upper = sub[iBlock].columnUpper();
    for (iColumn=0;iColumn<numberColumn2;iColumn++)
      upper[iColumn]=100.0;
#endif
    // and top matrix
    CoinPackedMatrix matrix = topMatrix;
    // and delete bits
    numberColumn2=0;
    for (iColumn=0;iColumn<numberColumns;iColumn++)
      if (iBlock!=columnBlock[iColumn])
        whichColumn[numberColumn2++]=iColumn;
    matrix.deleteCols(numberColumn2,whichColumn);
    top[iBlock]=matrix;
  }
  // and master
  numberRow2=0;
  numberColumn2=0;
  for (iRow=0;iRow<numberRows;iRow++)
    if (rowBlock[iRow]==-1)
      whichRow[numberRow2++]=iRow;
  for (iColumn=0;iColumn<numberColumns;iColumn++)
    if (columnBlock[iColumn]==-1)
      whichColumn[numberColumn2++]=iColumn;
  ClpModel masterModel(&model,numberRow2,whichRow,
                    numberColumn2,whichColumn);
  master=ClpSimplex(masterModel);
  delete [] whichRow;
  delete [] whichColumn;

  // Overkill in terms of space
  int numberMasterRows = master.numberRows();
  int * columnAdd = new int[numberBlocks+1];
  int * rowAdd = new int[numberBlocks*(numberMasterRows+1)];
  double * elementAdd = new double[numberBlocks*(numberMasterRows+1)];
  double * objective = new double[numberBlocks];
  int maxPass=500;
  int iPass;
  double lastObjective=1.0e31;
  // Create convexity rows for proposals 
  int numberMasterColumns = master.numberColumns();
  master.resize(numberMasterRows+numberBlocks,numberMasterColumns);
  // Arrays to say which block and when created
  int maximumColumns = 2*numberMasterRows+10*numberBlocks;
  int * whichBlock = new int[maximumColumns];
  int * when = new int[maximumColumns];
  int numberColumnsGenerated=numberBlocks;
  // fill in rhs and add in artificials
  {
    double * rowLower = master.rowLower();
    double * rowUpper = master.rowUpper();
    int iBlock;
    columnAdd[0] = 0;
    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
      int iRow = iBlock + numberMasterRows;;
      rowLower[iRow]=1.0;
      rowUpper[iRow]=1.0;
      rowAdd[iBlock] = iRow;
      elementAdd[iBlock] = 1.0;
      objective[iBlock] = 1.0e9;
      columnAdd[iBlock+1] = iBlock+1;
      when[iBlock]=-1;
      whichBlock[iBlock] = iBlock;
    }
    master.addColumns(numberBlocks,NULL,NULL,objective,
                      columnAdd, rowAdd, elementAdd);
  }
  // and resize matrix to double check clp will be happy
  //master.matrix()->setDimensions(numberMasterRows+numberBlocks,
  //                     numberMasterColumns+numberBlocks);

  for (iPass=0;iPass<maxPass;iPass++) {
    printf("Start of pass %d\n",iPass);
    // Solve master - may be infeasible
    master.scaling(false);
    if (0) {
      CoinMpsIO writer;
      writer.setMpsData(*master.matrix(), COIN_DBL_MAX,
                        master.getColLower(), master.getColUpper(),
                        master.getObjCoefficients(),
                        (const char*) 0 /*integrality*/,
                        master.getRowLower(), master.getRowUpper(),
                        NULL,NULL);
      writer.writeMps("yy.mps");
    }

    master.primal();
    if (master.numberIterations()==0&&iPass)
      break; // finished
    if (master.objectiveValue()>lastObjective-1.0e-7&&iPass>555)
      break; // finished
    lastObjective = master.objectiveValue();
    // mark basic ones and delete if necessary
    int iColumn;
    numberColumnsGenerated=master.numberColumns()-numberMasterColumns;
    for (iColumn=0;iColumn<numberColumnsGenerated;iColumn++) {
      if (master.getStatus(iColumn+numberMasterColumns)==ClpSimplex::basic)
        when[iColumn]=iPass;
    }
    if (numberColumnsGenerated+numberBlocks>maximumColumns) {
      // delete
      int numberKeep=0;
      int numberDelete=0;
      int * whichDelete = new int[numberColumns];
      for (iColumn=0;iColumn<numberColumnsGenerated;iColumn++) {
        if (when[iColumn]>iPass-7) {
          // keep
          when[numberKeep]=when[iColumn];
          whichBlock[numberKeep++]=whichBlock[iColumn];
        } else {
          // delete
          whichDelete[numberDelete++]=iColumn+numberMasterColumns;
        }
      }
      numberColumnsGenerated -= numberDelete;
      master.deleteColumns(numberDelete,whichDelete);
      delete [] whichDelete;
    }
    const double * dual=NULL;
    bool deleteDual=false;
    if (master.isProvenOptimal()) {
      dual = master.dualRowSolution();
    } else if (master.isProvenPrimalInfeasible()) {
      // could do composite objective
      dual = master.infeasibilityRay();
      deleteDual = true;
      printf("The sum of infeasibilities is %g\n",
             master.sumPrimalInfeasibilities());
    } else if (!master.numberColumns()) {
      assert(!iPass);
      dual = master.dualRowSolution();
      memset(master.dualRowSolution(),
             0,(numberMasterRows+numberBlocks)*sizeof(double));
    } else {
      abort();
    }
    // Create objective for sub problems and solve
    columnAdd[0]=0;
    int numberProposals=0;
    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
      int numberColumns2 = sub[iBlock].numberColumns();
      double * saveObj = new double [numberColumns2];
      double * objective2 = sub[iBlock].objective();
      memcpy(saveObj,objective2,numberColumns2*sizeof(double));
      // new objective
      top[iBlock].transposeTimes(dual,objective2);
      int i;
      if (master.isProvenOptimal()) {
        for (i=0;i<numberColumns2;i++)
          objective2[i] = saveObj[i]-objective2[i];
      } else {
        for (i=0;i<numberColumns2;i++)
          objective2[i] = -objective2[i];
      }
      sub[iBlock].primal();
      memcpy(objective2,saveObj,numberColumns2*sizeof(double));
      // get proposal
      if (sub[iBlock].numberIterations()||!iPass) {
        double objValue =0.0;
        int start = columnAdd[numberProposals];
        // proposal
        if (sub[iBlock].isProvenOptimal()) {
          const double * solution = sub[iBlock].primalColumnSolution();
          top[iBlock].times(solution,elementAdd+start);
          for (i=0;i<numberColumns2;i++)
            objValue += solution[i]*saveObj[i];
          // See if good dj and pack down
          int number=start;
          double dj = objValue;
          if (!master.isProvenOptimal()) 
            dj=0.0;
          double smallest=1.0e100;
          double largest=0.0;
          for (i=0;i<numberMasterRows;i++) {
            double value = elementAdd[start+i];
            if (fabs(value)>1.0e-15) {
              dj -= dual[i]*value;
              smallest = min(smallest,fabs(value));
              largest = max(largest,fabs(value));
              rowAdd[number]=i;
              elementAdd[number++]=value;
            }
          }
          // and convexity
          dj -= dual[numberMasterRows+iBlock];
          rowAdd[number]=numberMasterRows+iBlock;
          elementAdd[number++]=1.0;
          // if elements large then scale?
          //if (largest>1.0e8||smallest<1.0e-8)
          printf("For subproblem %d smallest - %g, largest %g - dj %g\n",
                 iBlock,smallest,largest,dj);
          if (dj<-1.0e-6||!iPass) {
            // take
            objective[numberProposals]=objValue;
            columnAdd[++numberProposals]=number;
            when[numberColumnsGenerated]=iPass;
            whichBlock[numberColumnsGenerated++]=iBlock;
          }
        } else if (sub[iBlock].isProvenDualInfeasible()) {
          // use ray
          const double * solution = sub[iBlock].unboundedRay();
          top[iBlock].times(solution,elementAdd+start);
          for (i=0;i<numberColumns2;i++)
            objValue += solution[i]*saveObj[i];
          // See if good dj and pack down
          int number=start;
          double dj = objValue;
          double smallest=1.0e100;
          double largest=0.0;
          for (i=0;i<numberMasterRows;i++) {
            double value = elementAdd[start+i];
            if (fabs(value)>1.0e-15) {
              dj -= dual[i]*value;
              smallest = min(smallest,fabs(value));
              largest = max(largest,fabs(value));
              rowAdd[number]=i;
              elementAdd[number++]=value;
            }
          }
          // if elements large or small then scale?
          //if (largest>1.0e8||smallest<1.0e-8)
            printf("For subproblem ray %d smallest - %g, largest %g - dj %g\n",
                   iBlock,smallest,largest,dj);
          if (dj<-1.0e-6) {
            // take
            objective[numberProposals]=objValue;
            columnAdd[++numberProposals]=number;
            when[numberColumnsGenerated]=iPass;
            whichBlock[numberColumnsGenerated++]=iBlock;
          }
        } else {
          abort();
        }
      }
    }
    if (deleteDual)
      delete [] dual;
    if (numberProposals) 
      master.addColumns(numberProposals,NULL,NULL,objective,
                        columnAdd,rowAdd,elementAdd);
  }
  // now put back a good solution
  double * lower = new double[numberMasterRows+numberBlocks];
  double * upper = new double[numberMasterRows+numberBlocks];
  numberColumnsGenerated  += numberMasterColumns;
  double * sol = new double[numberColumnsGenerated];
  const double * solution = master.primalColumnSolution();
  const double * masterLower = master.rowLower();
  const double * masterUpper = master.rowUpper();
  double * fullSolution = model.primalColumnSolution();
  const double * fullLower = model.columnLower();
  const double * fullUpper = model.columnUpper();
  int kColumn=0;
  for (iRow=0;iRow<numberRows;iRow++)
    model.setRowStatus(iRow,ClpSimplex::superBasic);
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    if (columnBlock[iColumn]==-1) {
      model.setStatus(iColumn,master.getStatus(kColumn));
      fullSolution[iColumn]=solution[kColumn++];
    }
  }
  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
    // convert top bit to by rows
    top[iBlock].reverseOrdering();
    // zero solution
    memset(sol,0,numberColumnsGenerated*sizeof(double));
    int i;
    for (i=numberMasterColumns;i<numberColumnsGenerated;i++)
      if (whichBlock[i-numberMasterColumns]==iBlock)
        sol[i] = solution[i];
    memset(lower,0,numberMasterRows*sizeof(double));
    master.times(1.0,sol,lower);
    for (iRow=0;iRow<numberMasterRows;iRow++) {
      double value = lower[iRow];
      if (masterUpper[iRow]<1.0e20) 
        upper[iRow] = value;
      else
        upper[iRow]=COIN_DBL_MAX;
      if (masterLower[iRow]>-1.0e20) 
        lower[iRow] = value;
      else
        lower[iRow]=-COIN_DBL_MAX;
    }
    sub[iBlock].addRows(numberMasterRows,lower,upper,
                        top[iBlock].getVectorStarts(),
                        top[iBlock].getVectorLengths(),
                        top[iBlock].getIndices(),
                        top[iBlock].getElements());
    sub[iBlock].primal();
    const double * subSolution = sub[iBlock].primalColumnSolution();
    // move solution
    kColumn=0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      if (columnBlock[iColumn]==iBlock) {
        model.setStatus(iColumn,sub[iBlock].getStatus(kColumn));
        fullSolution[iColumn]=subSolution[kColumn++];
      }
    }
    assert(kColumn==sub[iBlock].numberColumns());
  }
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    if (fullSolution[iColumn]<fullUpper[iColumn]-1.0e-8&&
        fullSolution[iColumn]>fullLower[iColumn]+1.0e-8) {
      model.setStatus(iColumn,ClpSimplex::basic);
    }
  }
  model.primal(1);
  delete [] sol;
  delete [] lower;
  delete [] upper;
  delete [] whichBlock;
  delete [] when;
  delete [] columnAdd;
  delete [] rowAdd;
  delete [] elementAdd;
  delete [] objective;
  delete [] top;
  delete [] sub;
  delete [] rowBlock;
  delete [] columnBlock;
  return 0;
}    

---