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

void ClpSimplexUnitTest ( const std::string &  mpsDir  )  [friend]

A function that tests the methods in the ClpSimplex class. The only reason for it not to be a member method is that this way it doesn't have to be compiled into the library. And that's a gain, because the library should be compiled with optimization on, but this method should be compiled with debugging.

It also does some testing of ClpFactorization class

Definition at line 652 of file unitTest.cpp.

{
  
  CoinRelFltEq eq(0.000001);

  {
    ClpSimplex solution;
  
    // matrix data
    //deliberate hiccup of 2 between 0 and 1
    CoinBigIndex start[5]={0,4,7,8,9};
    int length[5]={2,3,1,1,1};
    int rows[11]={0,2,-1,-1,0,1,2,0,1,2};
    double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1};
    CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length);
    
    // rim data
    double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0};
    double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10};
    double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10};
    double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0};
    double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0};
    
    // basis 1
    int rowBasis1[3]={-1,-1,-1};
    int colBasis1[5]={1,1,-1,-1,1};
    solution.loadProblem(matrix,colLower,colUpper,objective,
                   rowLower,rowUpper);
    int i;
    solution.createStatus();
    for (i=0;i<3;i++) {
      if (rowBasis1[i]<0) {
      solution.setRowStatus(i,ClpSimplex::atLowerBound);
      } else {
      solution.setRowStatus(i,ClpSimplex::basic);
      }
    }
    for (i=0;i<5;i++) {
      if (colBasis1[i]<0) {
      solution.setColumnStatus(i,ClpSimplex::atLowerBound);
      } else {
      solution.setColumnStatus(i,ClpSimplex::basic);
      }
    }
    solution.setLogLevel(3+4+8+16+32);
    solution.primal();
    for (i=0;i<3;i++) {
      if (rowBasis1[i]<0) {
      solution.setRowStatus(i,ClpSimplex::atLowerBound);
      } else {
      solution.setRowStatus(i,ClpSimplex::basic);
      }
    }
    for (i=0;i<5;i++) {
      if (colBasis1[i]<0) {
      solution.setColumnStatus(i,ClpSimplex::atLowerBound);
      } else {
      solution.setColumnStatus(i,ClpSimplex::basic);
      }
    }
    // intricate stuff does not work with scaling
    solution.scaling(0);
#ifndef NDEBUG
    int returnCode = solution.factorize ( );
    assert(!returnCode);
#else
    solution.factorize ( );
#endif
    const double * colsol = solution.primalColumnSolution();
    const double * rowsol = solution.primalRowSolution();
    solution.getSolution(rowsol,colsol);
#ifndef NDEBUG
    double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0};
    for (i=0;i<5;i++) {
      assert(eq(colsol[i],colsol1[i]));
    }
    // now feed in again without actually doing factorization
#endif
    ClpFactorization factorization2 = *solution.factorization();
    ClpSimplex solution2 = solution;
    solution2.setFactorization(factorization2);
    solution2.createStatus();
    for (i=0;i<3;i++) {
      if (rowBasis1[i]<0) {
      solution2.setRowStatus(i,ClpSimplex::atLowerBound);
      } else {
      solution2.setRowStatus(i,ClpSimplex::basic);
      }
    }
    for (i=0;i<5;i++) {
      if (colBasis1[i]<0) {
      solution2.setColumnStatus(i,ClpSimplex::atLowerBound);
      } else {
      solution2.setColumnStatus(i,ClpSimplex::basic);
      }
    }
    // intricate stuff does not work with scaling
    solution2.scaling(0);
    solution2.getSolution(rowsol,colsol);
    colsol = solution2.primalColumnSolution();
    rowsol = solution2.primalRowSolution();
    for (i=0;i<5;i++) {
      assert(eq(colsol[i],colsol1[i]));
    }
    solution2.setDualBound(0.1);
    solution2.dual();
    objective[2]=-1.0;
    objective[3]=-0.5;
    objective[4]=10.0;
    solution.dual();
    for (i=0;i<3;i++) {
      rowLower[i]=-1.0e20;
      colUpper[i+2]=0.0;
    }
    solution.setLogLevel(3);
    solution.dual();
    double rowObjective[]={1.0,0.5,-10.0};
    solution.loadProblem(matrix,colLower,colUpper,objective,
                   rowLower,rowUpper,rowObjective);
    solution.dual();
    solution.loadProblem(matrix,colLower,colUpper,objective,
                   rowLower,rowUpper,rowObjective);
    solution.primal();
  }
#ifndef COIN_NO_CLP_MESSAGE
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"exmip1";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex solution;
    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
                   m.getObjCoefficients(),
                   m.getRowLower(),m.getRowUpper());
    solution.dual();
    // Test event handling
    MyEventHandler handler;
    solution.passInEventHandler(&handler);
    int numberRows=solution.numberRows();
    // make sure values pass has something to do
    for (int i=0;i<numberRows;i++)
      solution.setRowStatus(i,ClpSimplex::basic);
    solution.primal(1);
    assert (solution.secondaryStatus()==102); // Came out at end of pass
  }
  // Test Message handler
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"exmip1";
    //fn = "Test/subGams4";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex model;
    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
                   m.getObjCoefficients(),
                   m.getRowLower(),m.getRowUpper());
    // Message handler
    MyMessageHandler messageHandler(&model);
    std::cout<<"Testing derived message handler"<<std::endl;
    model.passInMessageHandler(&messageHandler);
    model.messagesPointer()->setDetailMessage(1,102);
    model.setFactorizationFrequency(10);
    model.primal();
    model.primal(0,3);
    model.setObjCoeff(3,-2.9473684210526314);
    model.primal(0,3);
    // Write saved solutions
    int nc = model.getNumCols();
    int s; 
    std::deque<StdVectorDouble> fep = messageHandler.getFeasibleExtremePoints();
    int numSavedSolutions = fep.size();
    for ( s=0; s<numSavedSolutions; ++s ) {
      const StdVectorDouble & solnVec = fep[s];
      for ( int c=0; c<nc; ++c ) {
        if (fabs(solnVec[c])>1.0e-8)
          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
      }
    }
    // Solve again without scaling
    // and maximize then minimize
    messageHandler.clearFeasibleExtremePoints();
    model.scaling(0);
    model.setOptimizationDirection(-1);
    model.primal();
    model.setOptimizationDirection(1);
    model.primal();
    fep = messageHandler.getFeasibleExtremePoints();
    numSavedSolutions = fep.size();
    for ( s=0; s<numSavedSolutions; ++s ) {
      const StdVectorDouble & solnVec = fep[s];
      for ( int c=0; c<nc; ++c ) {
        if (fabs(solnVec[c])>1.0e-8)
          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
      }
    }
  }
#endif
  // Test dual ranging
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"exmip1";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex model;
    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
                   m.getObjCoefficients(),
                   m.getRowLower(),m.getRowUpper());
    model.primal();
    int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
    double costIncrease[13];
    int sequenceIncrease[13];
    double costDecrease[13];
    int sequenceDecrease[13];
    // ranging
    model.dualRanging(13,which,costIncrease,sequenceIncrease,
                  costDecrease,sequenceDecrease);
    int i;
    for ( i=0;i<13;i++)
      printf("%d increase %g %d, decrease %g %d\n",
           i,costIncrease[i],sequenceIncrease[i],
           costDecrease[i],sequenceDecrease[i]);
    assert (fabs(costDecrease[3])<1.0e-4);
    assert (fabs(costIncrease[7]-1.0)<1.0e-4);
    model.setOptimizationDirection(-1);
    {
      int j;
      double * obj = model.objective();
      int n=model.numberColumns();
      for (j=0;j<n;j++) 
      obj[j] *= -1.0;
    }
    double costIncrease2[13];
    int sequenceIncrease2[13];
    double costDecrease2[13];
    int sequenceDecrease2[13];
    // ranging
    model.dualRanging(13,which,costIncrease2,sequenceIncrease2,
                  costDecrease2,sequenceDecrease2);
    for (i=0;i<13;i++) {
      assert (fabs(costIncrease[i]-costDecrease2[i])<1.0e-6);
      assert (fabs(costDecrease[i]-costIncrease2[i])<1.0e-6);
      assert (sequenceIncrease[i]==sequenceDecrease2[i]);
      assert (sequenceDecrease[i]==sequenceIncrease2[i]);
    }
    // Now delete all rows and see what happens
    model.deleteRows(model.numberRows(),which);
    model.primal();
    // ranging
    if (!model.dualRanging(8,which,costIncrease,sequenceIncrease,
                           costDecrease,sequenceDecrease)) {
      for (i=0;i<8;i++) {
        printf("%d increase %g %d, decrease %g %d\n",
               i,costIncrease[i],sequenceIncrease[i],
               costDecrease[i],sequenceDecrease[i]);
      }
    }
  }
  // Test primal ranging
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"exmip1";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex model;
    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
                   m.getObjCoefficients(),
                   m.getRowLower(),m.getRowUpper());
    model.primal();
    int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
    double valueIncrease[13];
    int sequenceIncrease[13];
    double valueDecrease[13];
    int sequenceDecrease[13];
    // ranging
    model.primalRanging(13,which,valueIncrease,sequenceIncrease,
                  valueDecrease,sequenceDecrease);
    int i;
    for ( i=0;i<13;i++)
      printf("%d increase %g %d, decrease %g %d\n",
           i,valueIncrease[i],sequenceIncrease[i],
           valueDecrease[i],sequenceDecrease[i]);
    assert (fabs(valueIncrease[3]-0.642857)<1.0e-4);
    assert (fabs(valueIncrease[8]-2.95113)<1.0e-4);
#if 0
    // out until I find optimization bug
    // Test parametrics
    ClpSimplexOther * model2 = (ClpSimplexOther *) (&model);
    double rhs[]={ 1.0,2.0,3.0,4.0,5.0};
    double endingTheta=1.0;
    model2->scaling(0);
    model2->setLogLevel(63);
    model2->parametrics(0.0,endingTheta,0.1,
                        NULL,NULL,rhs,rhs,NULL);
#endif
  }
  // Test binv etc
  {    
    /* 
       Wolsey : Page 130
       max 4x1 -  x2
       7x1 - 2x2    <= 14
       x2    <= 3
       2x1 - 2x2    <= 3
       x1 in Z+, x2 >= 0

       note slacks are -1 in Clp so signs may be different
    */
    
    int n_cols = 2;
    int n_rows = 3;
    
    double obj[2] = {-4.0, 1.0};
    double collb[2] = {0.0, 0.0};
    double colub[2] = {COIN_DBL_MAX, COIN_DBL_MAX};
    double rowlb[3] = {-COIN_DBL_MAX, -COIN_DBL_MAX, -COIN_DBL_MAX};
    double rowub[3] = {14.0, 3.0, 3.0};
    
    int rowIndices[5] =  {0,     2,    0,    1,    2};
    int colIndices[5] =  {0,     0,    1,    1,    1};
    double elements[5] = {7.0, 2.0, -2.0,  1.0, -2.0};
    CoinPackedMatrix M(true, rowIndices, colIndices, elements, 5);

    ClpSimplex model;
    model.loadProblem(M, collb, colub, obj, rowlb, rowub);
    model.dual(0,1); // keep factorization 
    
    //check that the tableau matches wolsey (B-1 A)
    // slacks in second part of binvA
    double * binvA = reinterpret_cast<double*> (malloc((n_cols+n_rows) * sizeof(double)));
    
    printf("B-1 A by row\n");
    int i;
    for( i = 0; i < n_rows; i++){
      model.getBInvARow(i, binvA,binvA+n_cols);
      printf("row: %d -> ",i);
      for(int j=0; j < n_cols+n_rows; j++){
      printf("%g, ", binvA[j]);
      }
      printf("\n");
    }
    // See if can re-use factorization AND arrays
    model.primal(0,3+4); // keep factorization
    // And do by column
    printf("B-1 A by column\n");
    for( i = 0; i < n_rows+n_cols; i++){
      model.getBInvACol(i, binvA);
      printf("column: %d -> ",i);
      for(int j=0; j < n_rows; j++){
      printf("%g, ", binvA[j]);
      }
      printf("\n");
    }
    /* Do twice -
       without and with scaling
    */
    // set scaling off
    model.scaling(0);
    for (int iPass=0;iPass<2;iPass++) {
      model.primal(0,3+4); // keep factorization
      const double * rowScale = model.rowScale();
      const double * columnScale = model.columnScale();
      if (!iPass)
        assert (!rowScale);
      else
        assert (rowScale); // only true for this example
      /* has to be exactly correct as in OsiClpsolverInterface.cpp
         (also redo each pass as may change
      */
      printf("B-1 A");
      for( i = 0; i < n_rows; i++){
        model.getBInvARow(i, binvA,binvA+n_cols);
        printf("\nrow: %d -> ",i);
        int j;
        // First columns
        for(j=0; j < n_cols; j++){
        if (binvA[j]) {
          printf("(%d %g), ", j, binvA[j]);
        }
        }
        // now rows
        for(j=0; j < n_rows; j++){
        if (binvA[j+n_cols]) {
          printf("(%d %g), ", j+n_cols, binvA[j+n_cols]);
        }
        }
      }
      printf("\n");
      printf("And by column (trickier)");
      const int * pivotVariable = model.pivotVariable();
      for( i = 0; i < n_cols+n_rows; i++){
        model.getBInvACol(i, binvA);
        printf("\ncolumn: %d -> ",i);
        for(int j=0; j < n_rows; j++){
        if (binvA[j]) {
          // need to know pivot variable for +1/-1 (slack) and row/column scaling
          int pivot = pivotVariable[j];
          if (pivot<n_cols) {
            // scaled coding is in just in case 
            if (!columnScale) {
            printf("(%d %g), ", j, binvA[j]);
            } else {
            printf("(%d %g), ", j, binvA[j]*columnScale[pivot]);
            }
          } else {
            if (!rowScale) {
            printf("(%d %g), ", j, binvA[j]);
            } else {
            printf("(%d %g), ", j, binvA[j]/rowScale[pivot-n_cols]);
            }
          }
        }
        }
      }
      printf("\n");
      printf("binvrow");
      for( i = 0; i < n_rows; i++){
      model.getBInvRow(i, binvA);
      printf("\nrow: %d -> ",i);
      int j;
      for (j=0;j<n_rows;j++) {
        if (binvA[j])
          printf("(%d %g), ", j, binvA[j]);
      }
      }
      printf("\n");
      printf("And by column ");
      for( i = 0; i < n_rows; i++){
      model.getBInvCol(i, binvA);
      printf("\ncol: %d -> ",i);
      int j;
      for (j=0;j<n_rows;j++) {
        if (binvA[j])
          printf("(%d %g), ", j, binvA[j]);
      }
      }
      printf("\n");
      // now deal with next pass
      if (!iPass) {
      // get scaling for testing
      model.scaling(1);
      }
    }
    free(binvA);
    model.setColUpper(1,2.0);
    model.dual(0,2+4); // use factorization and arrays
    model.dual(0,2); // hopefully will not use factorization
    model.primal(0,3+4); // keep factorization
    // but say basis has changed
    model.setWhatsChanged(model.whatsChanged()&(~512));
    model.dual(0,2); // hopefully will not use factorization
  }
  // test steepest edge
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"finnis";
    int returnCode = m.readMps(fn.c_str(),"mps");
    if (returnCode) {
      // probable cause is that gz not there
      fprintf(stderr,"Unable to open finnis.mps in %s\n", dirSample.c_str());
      fprintf(stderr,"Most probable cause is finnis.mps is gzipped i.e. finnis.mps.gz and libz has not been activated\n");
      fprintf(stderr,"Either gunzip files or edit Makefiles/Makefile.location to get libz\n");
      exit(999);
    }
    ClpModel model;
    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),
                m.getColUpper(),
                m.getObjCoefficients(),
                m.getRowLower(),m.getRowUpper());
    ClpSimplex solution(model);

    solution.scaling(1); 
    solution.setDualBound(1.0e8);
    //solution.factorization()->maximumPivots(1);
    //solution.setLogLevel(3);
    solution.setDualTolerance(1.0e-7);
    // set objective sense,
    ClpDualRowSteepest steep;
    solution.setDualRowPivotAlgorithm(steep);
    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
    solution.dual();
  }
  // test normal solution
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"afiro";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex solution;
    ClpModel model;
    // do twice - without and with scaling
    int iPass;
    for (iPass=0;iPass<2;iPass++) {
      // explicit row objective for testing
      int nr = m.getNumRows();
      double * rowObj = new double[nr];
      CoinFillN(rowObj,nr,0.0);
      model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
                  m.getObjCoefficients(),
                  m.getRowLower(),m.getRowUpper(),rowObj);
      delete [] rowObj;
      solution = ClpSimplex(model);
      if (iPass) {
      solution.scaling();
      }
      solution.dual();
      solution.dual();
      // test optimal
      assert (solution.status()==0);
      int numberColumns = solution.numberColumns();
      int numberRows = solution.numberRows();
      CoinPackedVector colsol(numberColumns,solution.primalColumnSolution());
      double * objective = solution.objective();
#ifndef NDEBUG
      double objValue = colsol.dotProduct(objective);
#endif
      CoinRelFltEq eq(1.0e-8);
      assert(eq(objValue,-4.6475314286e+02));
      solution.dual();
      assert(eq(solution.objectiveValue(),-4.6475314286e+02));
      double * lower = solution.columnLower();
      double * upper = solution.columnUpper();
      double * sol = solution.primalColumnSolution();
      double * result = new double[numberColumns];
      CoinFillN ( result, numberColumns,0.0);
      solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
      int iRow , iColumn;
      // see if feasible and dual feasible
      for (iColumn=0;iColumn<numberColumns;iColumn++) {
      double value = sol[iColumn];
      assert(value<upper[iColumn]+1.0e-8);
      assert(value>lower[iColumn]-1.0e-8);
      value = objective[iColumn]-result[iColumn];
      assert (value>-1.0e-5);
      if (sol[iColumn]>1.0e-5)
        assert (value<1.0e-5);
      }
      delete [] result;
      result = new double[numberRows];
      CoinFillN ( result, numberRows,0.0);
      solution.matrix()->times(colsol, result);
      lower = solution.rowLower();
      upper = solution.rowUpper();
      sol = solution.primalRowSolution();
#ifndef NDEBUG
      for (iRow=0;iRow<numberRows;iRow++) {
      double value = result[iRow];
      assert(eq(value,sol[iRow]));
      assert(value<upper[iRow]+1.0e-8);
      assert(value>lower[iRow]-1.0e-8);
      }
#endif
      delete [] result;
      // test row objective
      double * rowObjective = solution.rowObjective();
      CoinDisjointCopyN(solution.dualRowSolution(),numberRows,rowObjective);
      CoinDisjointCopyN(solution.dualColumnSolution(),numberColumns,objective);
      // this sets up all slack basis
      solution.createStatus();
      solution.dual();
      CoinFillN(rowObjective,numberRows,0.0);
      CoinDisjointCopyN(m.getObjCoefficients(),numberColumns,objective);
      solution.dual();
    }
  }
  // test unbounded
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"brandy";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex solution;
    // do twice - without and with scaling
    int iPass;
    for (iPass=0;iPass<2;iPass++) {
      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
                  m.getObjCoefficients(),
                  m.getRowLower(),m.getRowUpper());
      if (iPass)
      solution.scaling();
      solution.setOptimizationDirection(-1);
      // test unbounded and ray
#ifdef DUAL
      solution.setDualBound(100.0);
      solution.dual();
#else
      solution.primal();
#endif
      assert (solution.status()==2);
      int numberColumns = solution.numberColumns();
      int numberRows = solution.numberRows();
      double * lower = solution.columnLower();
      double * upper = solution.columnUpper();
      double * sol = solution.primalColumnSolution();
      double * ray = solution.unboundedRay();
      double * objective = solution.objective();
      double objChange=0.0;
      int iRow , iColumn;
      // make sure feasible and columns form ray
      for (iColumn=0;iColumn<numberColumns;iColumn++) {
      double value = sol[iColumn];
      assert(value<upper[iColumn]+1.0e-8);
      assert(value>lower[iColumn]-1.0e-8);
      value = ray[iColumn];
      if (value>0.0)
        assert(upper[iColumn]>1.0e30);
      else if (value<0.0)
        assert(lower[iColumn]<-1.0e30);
      objChange += value*objective[iColumn];
      }
      // make sure increasing objective
      assert(objChange>0.0);
      double * result = new double[numberRows];
      CoinFillN ( result, numberRows,0.0);
      solution.matrix()->times(sol, result);
      lower = solution.rowLower();
      upper = solution.rowUpper();
      sol = solution.primalRowSolution();
#ifndef NDEBUG
      for (iRow=0;iRow<numberRows;iRow++) {
      double value = result[iRow];
      assert(eq(value,sol[iRow]));
      assert(value<upper[iRow]+2.0e-8);
      assert(value>lower[iRow]-2.0e-8);
      }
#endif
      CoinFillN ( result, numberRows,0.0);
      solution.matrix()->times(ray, result);
      // there may be small differences (especially if scaled)
      for (iRow=0;iRow<numberRows;iRow++) {
      double value = result[iRow];
      if (value>1.0e-8)
        assert(upper[iRow]>1.0e30);
      else if (value<-1.0e-8)
        assert(lower[iRow]<-1.0e30);
      }
      delete [] result;
      delete [] ray;
    }
  }
  // test infeasible
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"brandy";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex solution;
    // do twice - without and with scaling
    int iPass;
    for (iPass=0;iPass<2;iPass++) {
      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
                  m.getObjCoefficients(),
                  m.getRowLower(),m.getRowUpper());
      if (iPass)
      solution.scaling();
      // test infeasible and ray
      solution.columnUpper()[0]=0.0;
#ifdef DUAL
      solution.setDualBound(100.0);
      solution.dual();
#else
      solution.primal();
#endif
      assert (solution.status()==1);
      int numberColumns = solution.numberColumns();
      int numberRows = solution.numberRows();
      double * lower = solution.rowLower();
      double * upper = solution.rowUpper();
      double * ray = solution.infeasibilityRay();
      assert(ray);
      // construct proof of infeasibility
      int iRow , iColumn;
      double lo=0.0,up=0.0;
      int nl=0,nu=0;
      for (iRow=0;iRow<numberRows;iRow++) {
      if (lower[iRow]>-1.0e20) {
        lo += ray[iRow]*lower[iRow];
      } else {
        if (ray[iRow]>1.0e-8) 
          nl++;
      }
      if (upper[iRow]<1.0e20) {
        up += ray[iRow]*upper[iRow];
      } else {
        if (ray[iRow]>1.0e-8) 
          nu++;
      }
      }
      if (nl)
      lo=-1.0e100;
      if (nu)
      up=1.0e100;
      double * result = new double[numberColumns];
      double lo2=0.0,up2=0.0;
      CoinFillN ( result, numberColumns,0.0);
      solution.matrix()->transposeTimes(ray, result);
      lower = solution.columnLower();
      upper = solution.columnUpper();
      nl=nu=0;
      for (iColumn=0;iColumn<numberColumns;iColumn++) {
      if (result[iColumn]>1.0e-8) {
        if (lower[iColumn]>-1.0e20)
          lo2 += result[iColumn]*lower[iColumn];
        else
          nl++;
        if (upper[iColumn]<1.0e20)
          up2 += result[iColumn]*upper[iColumn];
        else
          nu++;
      } else if (result[iColumn]<-1.0e-8) {
        if (lower[iColumn]>-1.0e20)
          up2 += result[iColumn]*lower[iColumn];
        else
          nu++;
        if (upper[iColumn]<1.0e20)
          lo2 += result[iColumn]*upper[iColumn];
        else
          nl++;
      }
      }
      if (nl)
      lo2=-1.0e100;
      if (nu)
      up2=1.0e100;
      // make sure inconsistency
      assert(lo2>up||up2<lo);
      delete [] result;
      delete [] ray;
    }
  }
  // test delete and add
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"brandy";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex solution;
    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
                   m.getObjCoefficients(),
                   m.getRowLower(),m.getRowUpper());
    solution.dual();
    CoinRelFltEq eq(1.0e-8);
    assert(eq(solution.objectiveValue(),1.5185098965e+03));

    int numberColumns = solution.numberColumns();
    int numberRows = solution.numberRows();
    double * saveObj = new double [numberColumns];
    double * saveLower = new double[numberRows+numberColumns];
    double * saveUpper = new double[numberRows+numberColumns];
    int * which = new int [numberRows+numberColumns];

    int numberElements = m.getMatrixByCol()->getNumElements();
    int * starts = new int[numberRows+numberColumns];
    int * index = new int[numberElements];
    double * element = new double[numberElements];

    const CoinBigIndex * startM;
    const int * lengthM;
    const int * indexM;
    const double * elementM;

    int n,nel;

    // delete non basic columns
    n=0;
    nel=0;
    int iRow , iColumn;
    const double * lower = m.getColLower();
    const double * upper = m.getColUpper();
    const double * objective = m.getObjCoefficients();
    startM = m.getMatrixByCol()->getVectorStarts();
    lengthM = m.getMatrixByCol()->getVectorLengths();
    indexM = m.getMatrixByCol()->getIndices();
    elementM = m.getMatrixByCol()->getElements();
    starts[0]=0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      if (solution.getColumnStatus(iColumn)!=ClpSimplex::basic) {
      saveObj[n]=objective[iColumn];
      saveLower[n]=lower[iColumn];
      saveUpper[n]=upper[iColumn];
      int j;
      for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
        index[nel]=indexM[j];
        element[nel++]=elementM[j];
      }
      which[n++]=iColumn;
      starts[n]=nel;
      }
    }
    solution.deleteColumns(n,which);
    solution.dual();
    // Put back
    solution.addColumns(n,saveLower,saveUpper,saveObj,
                  starts,index,element);
    solution.dual();
    assert(eq(solution.objectiveValue(),1.5185098965e+03));
    // Delete all columns and add back
    n=0;
    nel=0;
    starts[0]=0;
    lower = m.getColLower();
    upper = m.getColUpper();
    objective = m.getObjCoefficients();
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      saveObj[n]=objective[iColumn];
      saveLower[n]=lower[iColumn];
      saveUpper[n]=upper[iColumn];
      int j;
      for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
      index[nel]=indexM[j];
      element[nel++]=elementM[j];
      }
      which[n++]=iColumn;
      starts[n]=nel;
    }
    solution.deleteColumns(n,which);
    solution.dual();
    // Put back
    solution.addColumns(n,saveLower,saveUpper,saveObj,
                  starts,index,element);
    solution.dual();
    assert(eq(solution.objectiveValue(),1.5185098965e+03));

    // reload with original
    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
                   m.getObjCoefficients(),
                   m.getRowLower(),m.getRowUpper());
    // delete half rows
    n=0;
    nel=0;
    lower = m.getRowLower();
    upper = m.getRowUpper();
    startM = m.getMatrixByRow()->getVectorStarts();
    lengthM = m.getMatrixByRow()->getVectorLengths();
    indexM = m.getMatrixByRow()->getIndices();
    elementM = m.getMatrixByRow()->getElements();
    starts[0]=0;
    for (iRow=0;iRow<numberRows;iRow++) {
      if ((iRow&1)==0) {
      saveLower[n]=lower[iRow];
      saveUpper[n]=upper[iRow];
      int j;
      for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
        index[nel]=indexM[j];
        element[nel++]=elementM[j];
      }
      which[n++]=iRow;
      starts[n]=nel;
      }
    }
    solution.deleteRows(n,which);
    solution.dual();
    // Put back
    solution.addRows(n,saveLower,saveUpper,
                  starts,index,element);
    solution.dual();
    assert(eq(solution.objectiveValue(),1.5185098965e+03));
    solution.writeMps("yy.mps");
    // Delete all rows
    n=0;
    nel=0;
    lower = m.getRowLower();
    upper = m.getRowUpper();
    starts[0]=0;
    for (iRow=0;iRow<numberRows;iRow++) {
      saveLower[n]=lower[iRow];
      saveUpper[n]=upper[iRow];
      int j;
      for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
      index[nel]=indexM[j];
      element[nel++]=elementM[j];
      }
      which[n++]=iRow;
      starts[n]=nel;
    }
    solution.deleteRows(n,which);
    solution.dual();
    // Put back
    solution.addRows(n,saveLower,saveUpper,
                  starts,index,element);
    solution.dual();
    solution.writeMps("xx.mps");
    assert(eq(solution.objectiveValue(),1.5185098965e+03));
    // Zero out status array to give some interest
    memset(solution.statusArray()+numberColumns,0,numberRows);
    solution.primal(1);
    assert(eq(solution.objectiveValue(),1.5185098965e+03));
    // Delete all columns and rows
    n=0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      which[n++]=iColumn;
      starts[n]=nel;
    }
    solution.deleteColumns(n,which);
    n=0;
    for (iRow=0;iRow<numberRows;iRow++) {
      which[n++]=iRow;
      starts[n]=nel;
    }
    solution.deleteRows(n,which);

    delete [] saveObj;
    delete [] saveLower;
    delete [] saveUpper;
    delete [] which;
    delete [] starts;
    delete [] index;
    delete [] element;
  }
#if 1
  // Test barrier
  {
    CoinMpsIO m;
    std::string fn = dirSample+"exmip1";
    m.readMps(fn.c_str(),"mps");
    ClpInterior solution;
    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
                   m.getObjCoefficients(),
                   m.getRowLower(),m.getRowUpper());
    solution.primalDual();
  }
#endif
  // test network 
#define QUADRATIC
  if (1) {    
    std::string fn = dirSample+"input.130";
    int numberColumns;
    int numberRows;
    
    FILE * fp = fopen(fn.c_str(),"r");
    if (!fp) {
      // Try in Data/Sample
      fn = "Data/Sample/input.130";
      fp = fopen(fn.c_str(),"r");
    }
    if (!fp) {
      fprintf(stderr,"Unable to open file input.130 in dirSample or Data/Sample directory\n");
    } else {
      int problem;
      char temp[100];
      // read and skip 
      int x=fscanf(fp,"%s",temp);
      if (x<0)
      throw("bad fscanf");
      assert (!strcmp(temp,"BEGIN"));
      x=fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows, 
           &numberColumns);
      if (x<0)
      throw("bad fscanf");
      // scan down to SUPPLY
      while (fgets(temp,100,fp)) {
      if (!strncmp(temp,"SUPPLY",6))
        break;
      }
      if (strncmp(temp,"SUPPLY",6)) {
      fprintf(stderr,"Unable to find SUPPLY\n");
      exit(2);
      }
      // get space for rhs
      double * lower = new double[numberRows];
      double * upper = new double[numberRows];
      int i;
      for (i=0;i<numberRows;i++) {
      lower[i]=0.0;
      upper[i]=0.0;
      }
      // ***** Remember to convert to C notation
      while (fgets(temp,100,fp)) {
      int row;
      int value;
      if (!strncmp(temp,"ARCS",4))
        break;
      sscanf(temp,"%d %d",&row,&value);
      upper[row-1]=-value;
      lower[row-1]=-value;
      }
      if (strncmp(temp,"ARCS",4)) {
      fprintf(stderr,"Unable to find ARCS\n");
      exit(2);
      }
      // number of columns may be underestimate
      int * head = new int[2*numberColumns];
      int * tail = new int[2*numberColumns];
      int * ub = new int[2*numberColumns];
      int * cost = new int[2*numberColumns];
      // ***** Remember to convert to C notation
      numberColumns=0;
      while (fgets(temp,100,fp)) {
      int iHead;
      int iTail;
      int iUb;
      int iCost;
      if (!strncmp(temp,"DEMAND",6))
        break;
      sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb);
      iHead--;
      iTail--;
      head[numberColumns]=iHead;
      tail[numberColumns]=iTail;
      ub[numberColumns]=iUb;
      cost[numberColumns]=iCost;
      numberColumns++;
      }
      if (strncmp(temp,"DEMAND",6)) {
      fprintf(stderr,"Unable to find DEMAND\n");
      exit(2);
      }
      // ***** Remember to convert to C notation
      while (fgets(temp,100,fp)) {
      int row;
      int value;
      if (!strncmp(temp,"END",3))
        break;
      sscanf(temp,"%d %d",&row,&value);
      upper[row-1]=value;
      lower[row-1]=value;
      }
      if (strncmp(temp,"END",3)) {
      fprintf(stderr,"Unable to find END\n");
      exit(2);
      }
      printf("Problem %d has %d rows and %d columns\n",problem,
           numberRows,numberColumns);
      fclose(fp);
      ClpSimplex  model;
      // now build model
      
      double * objective =new double[numberColumns];
      double * lowerColumn = new double[numberColumns];
      double * upperColumn = new double[numberColumns];
      
      double * element = new double [2*numberColumns];
      CoinBigIndex * start = new CoinBigIndex [numberColumns+1];
      int * row = new int[2*numberColumns];
      start[numberColumns]=2*numberColumns;
      for (i=0;i<numberColumns;i++) {
      start[i]=2*i;
      element[2*i]=-1.0;
      element[2*i+1]=1.0;
      row[2*i]=head[i];
      row[2*i+1]=tail[i];
      lowerColumn[i]=0.0;
      upperColumn[i]=ub[i];
      objective[i]=cost[i];
      }
      // Create Packed Matrix
      CoinPackedMatrix matrix;
      int * lengths = NULL;
      matrix.assignMatrix(true,numberRows,numberColumns,
                    2*numberColumns,element,row,start,lengths);
      // load model
      model.loadProblem(matrix,
                  lowerColumn,upperColumn,objective,
                  lower,upper);
      model.factorization()->maximumPivots(200+model.numberRows()/100);
      model.createStatus();
      double time1 = CoinCpuTime();
      model.dual();
      std::cout<<"Network problem, ClpPackedMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
      ClpPlusMinusOneMatrix * plusMinus = new ClpPlusMinusOneMatrix(matrix);
      assert (plusMinus->getIndices()); // would be zero if not +- one
      //ClpPlusMinusOneMatrix *plusminus_matrix;

      //plusminus_matrix = new ClpPlusMinusOneMatrix;

      //plusminus_matrix->passInCopy(numberRows, numberColumns, true, plusMinus->getMutableIndices(),
      //                         plusMinus->startPositive(),plusMinus->startNegative());
      model.loadProblem(*plusMinus,
                  lowerColumn,upperColumn,objective,
                  lower,upper);
      //model.replaceMatrix( plusminus_matrix , true);
      delete plusMinus;
      //model.createStatus();
      //model.initialSolve();
      //model.writeMps("xx.mps");
      
      model.factorization()->maximumPivots(200+model.numberRows()/100);
      model.createStatus();
      time1 = CoinCpuTime();
      model.dual();
      std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
      ClpNetworkMatrix network(numberColumns,head,tail);
      model.loadProblem(network,
                  lowerColumn,upperColumn,objective,
                  lower,upper);
      
      model.factorization()->maximumPivots(200+model.numberRows()/100);
      model.createStatus();
      time1 = CoinCpuTime();
      model.dual();
      std::cout<<"Network problem, ClpNetworkMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
      delete [] lower;
      delete [] upper;
      delete [] head;
      delete [] tail;
      delete [] ub;
      delete [] cost;
      delete [] objective;
      delete [] lowerColumn;
      delete [] upperColumn;
    }
  }
#ifdef QUADRATIC
  // Test quadratic to solve linear
  if (1) {    
    CoinMpsIO m;
    std::string fn = dirSample+"exmip1";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex solution;
    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
                   m.getObjCoefficients(),
                   m.getRowLower(),m.getRowUpper());
    //solution.dual();
    // get quadratic part
    int numberColumns=solution.numberColumns();
    int * start=new int [numberColumns+1];
    int * column = new int[numberColumns];
    double * element = new double[numberColumns];
    int i;
    start[0]=0;
    int n=0;
    int kk=numberColumns-1;
    int kk2=numberColumns-1;
    for (i=0;i<numberColumns;i++) {
      if (i>=kk) {
      column[n]=i;
      if (i>=kk2)
        element[n]=1.0e-1;
      else
        element[n]=0.0;
      n++;
      }
      start[i+1]=n;
    }
    // Load up objective
    solution.loadQuadraticObjective(numberColumns,start,column,element);
    delete [] start;
    delete [] column;
    delete [] element;
    //solution.quadraticSLP(50,1.0e-4);
    CoinRelFltEq eq(1.0e-4);
    //assert(eq(objValue,820.0));
    //solution.setLogLevel(63);
    solution.primal();
    printSol(solution);
    //assert(eq(objValue,3.2368421));
    //exit(77);
  }
  // Test quadratic
  if (1) {    
    CoinMpsIO m;
    std::string fn = dirSample+"share2qp";
    //fn = "share2qpb";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex model;
    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
                   m.getObjCoefficients(),
                   m.getRowLower(),m.getRowUpper());
    model.dual();
    // get quadratic part
    int * start=NULL;
    int * column = NULL;
    double * element = NULL;
    m.readQuadraticMps(NULL,start,column,element,2);
    int column2[200];
    double element2[200];
    int start2[80];
    int j;
    start2[0]=0;
    int nel=0;
    bool good=false;
    for (j=0;j<79;j++) {
      if (start[j]==start[j+1]) {
      column2[nel]=j;
      element2[nel]=0.0;
      nel++;
      } else {
      int i;
      for (i=start[j];i<start[j+1];i++) {
        column2[nel]=column[i];
        element2[nel++]=element[i];
      }
      }
      start2[j+1]=nel;
    }
    // Load up objective
    if (good)
      model.loadQuadraticObjective(model.numberColumns(),start2,column2,element2);
    else
      model.loadQuadraticObjective(model.numberColumns(),start,column,element);
    delete [] start;
    delete [] column;
    delete [] element;
    int numberColumns=model.numberColumns();
    model.scaling(0);
#if 0
    model.nonlinearSLP(50,1.0e-4);
#else
    // Get feasible
    ClpObjective * saveObjective = model.objectiveAsObject()->clone();
    ClpLinearObjective zeroObjective(NULL,numberColumns);
    model.setObjective(&zeroObjective);
    model.dual();
    model.setObjective(saveObjective);
    delete saveObjective;
#endif
    //model.setLogLevel(63);
    //exit(77);
    model.setFactorizationFrequency(10);
    model.primal();
    printSol(model);
#ifndef NDEBUG
    double objValue = model.getObjValue();
#endif
    CoinRelFltEq eq(1.0e-4);
    assert(eq(objValue,-400.92));
    // and again for barrier
    model.barrier(false);
    //printSol(model);
    model.allSlackBasis();
    model.primal();
    //printSol(model);
  }
  if (0) {    
    CoinMpsIO m;
    std::string fn = "./beale";
    //fn = "./jensen";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex solution;
    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
                   m.getObjCoefficients(),
                   m.getRowLower(),m.getRowUpper());
    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
    solution.dual();
    // get quadratic part
    int * start=NULL;
    int * column = NULL;
    double * element = NULL;
    m.readQuadraticMps(NULL,start,column,element,2);
    // Load up objective
    solution.loadQuadraticObjective(solution.numberColumns(),start,column,element);
    delete [] start;
    delete [] column;
    delete [] element;
    solution.primal(1);
    solution.nonlinearSLP(50,1.0e-4);
    double objValue = solution.getObjValue();
    CoinRelFltEq eq(1.0e-4);
    assert(eq(objValue,0.5));
    solution.primal();
    objValue = solution.getObjValue();
    assert(eq(objValue,0.5));
  }
#endif
  // Test CoinStructuredModel
  {

    // Sub block
    CoinModel sub;
    {
      // matrix data
      //deliberate hiccup of 2 between 0 and 1
      CoinBigIndex start[5]={0,4,7,8,9};
      int length[5]={2,3,1,1,1};
      int rows[11]={0,2,-1,-1,0,1,2,0,1,2};
      double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1};
      CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length);
      // by row
      matrix.reverseOrdering();
      const double * element = matrix.getElements();
      const int * column = matrix.getIndices();
      const CoinBigIndex * rowStart = matrix.getVectorStarts();
      const int * rowLength = matrix.getVectorLengths();
      
      // rim data
      //double objective[5]={-4.0,1.0,0.0,0.0,0.0};
      double rowLower[3]={14.0,3.0,3.0};
      double rowUpper[3]={14.0,3.0,3.0};
      //double columnLower[5]={0.0,0.0,0.0,0.0,0.0};
      //double columnUpper[5]={100.0,100.0,100.0,100.0,100.0};
      
      for (int i=0;i<3;i++) {
      sub.addRow(rowLength[i],column+rowStart[i],
               element+rowStart[i],rowLower[i],rowUpper[i]);
      }
      //for (int i=0;i<5;i++) {
      //sub.setColumnBounds(i,columnLower[i],columnUpper[i]);
      //sub.setColumnObjective(i,objective[i]);
      //}
      sub.convertMatrix();
    }
    // Top
    CoinModel top;
    {
      // matrix data
      CoinBigIndex start[5]={0,2,4,6,8};
      int length[5]={2,2,2,2,2};
      int rows[10]={0,1,0,1,0,1,0,1,0,1};
      double elements[10]={1,1,1,1,1,1,1,1,1,1};
      CoinPackedMatrix matrix(true,2,5,8,elements,rows,start,length);
      // by row
      matrix.reverseOrdering();
      const double * element = matrix.getElements();
      const int * column = matrix.getIndices();
      const CoinBigIndex * rowStart = matrix.getVectorStarts();
      const int * rowLength = matrix.getVectorLengths();
      
      // rim data
      double objective[5]={-4.0,1.0,0.0,0.0,0.0};
      //double rowLower[3]={14.0,3.0,3.0};
      //double rowUpper[3]={14.0,3.0,3.0};
      double columnLower[5]={0.0,0.0,0.0,0.0,0.0};
      double columnUpper[5]={100.0,100.0,100.0,100.0,100.0};
      
      for (int i=0;i<2;i++) {
      top.addRow(rowLength[i],column+rowStart[i],
               element+rowStart[i],
               -COIN_DBL_MAX,COIN_DBL_MAX);
      }
      for (int i=0;i<5;i++) {
      top.setColumnBounds(i,columnLower[i],columnUpper[i]);
      top.setColumnObjective(i,objective[i]);
      }
      top.convertMatrix();
    }
    // Create a structured model
    CoinStructuredModel structured;
    int numberBlocks=5;
    for (int i=0;i<numberBlocks;i++) {
      std::string topName="row_master";
      std::string blockName="block_";
      char bName = static_cast<char>('a'+static_cast<char>(i));
      blockName.append(1,bName);
      structured.addBlock(topName,blockName,top);
      structured.addBlock(blockName,blockName,sub);
    }
    // Set rhs on first block
    CoinModel * first = structured.coinBlock(0);
    for (int i=0;i<2;i++) {
      first->setRowLower(i,0.0);
      first->setRowUpper(i,100.0);
    }
    // Refresh whats set
    structured.refresh(0);
    // Could perturb stuff, but for first go don't bother
    ClpSimplex fullModel;
    // There is no original stuff set - think
    fullModel.loadProblem(structured,false);
    fullModel.dual();
    fullModel.dropNames();
    fullModel.writeMps("test.mps");
    // Make up very simple nested model - not realistic
    // Create a structured model
    CoinStructuredModel structured2;
    numberBlocks=3;
    for (int i=0;i<numberBlocks;i++) {
      std::string blockName="block_";
      char bName = static_cast<char>('a'+static_cast<char>(i));
      blockName.append(1,bName);
      structured2.addBlock(blockName,blockName,structured);
    }
    fullModel.loadProblem(structured2,false);
    fullModel.dual();
    fullModel.dropNames();
    fullModel.writeMps("test2.mps");
  }
}


Generated by  Doxygen 1.6.0   Back to index