C++ CoinWarmStartBasis类代码示例

//===========================================================================//
CoinWarmStartBasis* UtilAlpsDecodeWarmStart(AlpsEncoded&       encoded,
      AlpsReturnStatus* rc)
{
   //rc not used? not checked?
   int numCols;
   int numRows;
   encoded.readRep(numCols);
   encoded.readRep(numRows);
   int tempInt;
   // Structural
   int nint = (numCols + 15) >> 4;
   char* structuralStatus = new char[4 * nint];
   encoded.readRep(structuralStatus, tempInt);
   assert(tempInt == nint * 4);
   // Artificial
   nint = (numRows + 15) >> 4;
   char* artificialStatus = new char[4 * nint];
   encoded.readRep(artificialStatus, tempInt);
   assert(tempInt == nint * 4);
   CoinWarmStartBasis* ws = new CoinWarmStartBasis();

   if (!ws) {
      throw CoinError("Out of memory", "UtilAlpsDecodeWarmStart", "HELP");
   }

   ws->assignBasisStatus(numCols, numRows,
                         structuralStatus, artificialStatus);
   return ws;
}

CoinWarmStartBasis *CoinPrePostsolveMatrix::getStatus ()

{ int n = ncols_ ;
  int m = nrows_ ;
  CoinWarmStartBasis *wsb = new CoinWarmStartBasis() ;
  wsb->setSize(n,m) ;
  for (int j = 0 ; j < n ; j++)
  { CoinWarmStartBasis::Status statj = 
	CoinWarmStartBasis::Status(getColumnStatus(j)) ;
    wsb->setStructStatus(j,statj) ; }
  for (int i = 0 ; i < m ; i++)
  { CoinWarmStartBasis::Status stati =
	CoinWarmStartBasis::Status(getRowStatus(i)) ;
    wsb->setArtifStatus(i,stati) ; }
  
  return (wsb) ; }

/* returns 0 if no solution, 1 if valid solution
   with better objective value than one passed in
   also returns list of nodes
   This does Fractional Diving
*/
int 
CbcHeuristicDive::fathom(CbcModel * model, int & numberNodes, 
			 CbcSubProblem ** & nodes)
{
  double solutionValue = model->getCutoff();
  numberNodes=0;
  // Get solution array for heuristic solution
  int numberColumns = model_->solver()->getNumCols();
  double * newSolution = new double [4*numberColumns];
  double * lastDjs = newSolution+numberColumns;
  double * originalLower = lastDjs+numberColumns;
  double * originalUpper = originalLower+numberColumns;
  memcpy(originalLower,model_->solver()->getColLower(),
	 numberColumns*sizeof(double));
  memcpy(originalUpper,model_->solver()->getColUpper(),
	 numberColumns*sizeof(double));
  int numberCuts=0;
  OsiRowCut ** cuts = NULL; //new OsiRowCut * [maxIterations_];
  nodes=new CbcSubProblem * [maxIterations_+2];
  int returnCode=solution(solutionValue,numberNodes,numberCuts,
			  cuts,nodes,
			  newSolution);

  if (returnCode==1) {
    // copy to best solution ? or put in solver
    printf("Solution from heuristic fathom\n");
  }
  int numberFeasibleNodes=numberNodes;
  if (returnCode!=1)
    numberFeasibleNodes--;
  if (numberFeasibleNodes>0) {
    CoinWarmStartBasis * basis = nodes[numberFeasibleNodes-1]->status_;
    //double * sort = new double [numberFeasibleNodes];
    //int * whichNode = new int [numberFeasibleNodes];
    //int numberNodesNew=0;
    // use djs on previous unless feasible
    for (int iNode=0;iNode<numberFeasibleNodes;iNode++) {
      CbcSubProblem * sub = nodes[iNode];
      double branchValue = sub->branchValue_;
      int iStatus=sub->problemStatus_;
      int iColumn = sub->branchVariable_;
      bool secondBranch = (iStatus&2)!=0;
      bool branchUp;
      if (!secondBranch)
	branchUp = (iStatus&1)!=0;
      else
	branchUp = (iStatus&1)==0;
      double djValue=lastDjs[iColumn];
      sub->djValue_=fabs(djValue);
      if (!branchUp&&floor(branchValue)==originalLower[iColumn]
	  &&basis->getStructStatus(iColumn) == CoinWarmStartBasis::atLowerBound) {
	if (djValue>0.0) {
	  // naturally goes to LB
	  printf("ignoring branch down on %d (node %d) from value of %g - branch was %s - dj %g\n",
		 iColumn,iNode,branchValue,secondBranch ? "second" : "first",
		 djValue);
	  sub->problemStatus_ |= 4;
	  //} else {
	  // put on list
	  //sort[numberNodesNew]=djValue;
	  //whichNode[numberNodesNew++]=iNode;
	}
      } else if (branchUp&&ceil(branchValue)==originalUpper[iColumn]
	  &&basis->getStructStatus(iColumn) == CoinWarmStartBasis::atUpperBound) {
	if (djValue<0.0) {
	  // naturally goes to UB
	  printf("ignoring branch up on %d (node %d) from value of %g - branch was %s - dj %g\n",
		 iColumn,iNode,branchValue,secondBranch ? "second" : "first",
		 djValue);
	  sub->problemStatus_ |= 4;
	  //} else {
	  // put on list
	  //sort[numberNodesNew]=-djValue;
	  //whichNode[numberNodesNew++]=iNode;
	}
      }
    }
    // use conflict to order nodes
    for (int iCut=0;iCut<numberCuts;iCut++) {
    }
    //CoinSort_2(sort,sort+numberNodesNew,whichNode);
    // create nodes
    // last node will have one way already done
  }
  for (int iCut=0;iCut<numberCuts;iCut++) {
    delete cuts[iCut];
  }
  delete [] cuts;
  delete [] newSolution;
  return returnCode;
}

// Apply subproblem
void
CbcSubProblem::apply(OsiSolverInterface * solver, int what) const
{
    int i;
    if ((what&1) != 0) {
#ifndef NDEBUG
        int nSame = 0;
#endif
        for (i = 0; i < numberChangedBounds_; i++) {
            int variable = variables_[i];
            int k = variable & 0x3fffffff;
            if ((variable&0x80000000) == 0) {
                // lower bound changing
                //#define CBC_PRINT2
#ifdef CBC_PRINT2
                if (solver->getColLower()[k] != newBounds_[i])
                    printf("lower change for column %d - from %g to %g\n", k, solver->getColLower()[k], newBounds_[i]);
#endif
#ifndef NDEBUG
                if ((variable&0x40000000) == 0 && true) {
                    double oldValue = solver->getColLower()[k];
                    assert (newBounds_[i] > oldValue - 1.0e-8);
                    if (newBounds_[i] < oldValue + 1.0e-8) {
#ifdef CBC_PRINT2
                        printf("bad null lower change for column %d - bound %g\n", k, oldValue);
#endif
                        if (newBounds_[i] == oldValue)
                            nSame++;
                    }
                }
#endif
                solver->setColLower(k, newBounds_[i]);
            } else {
                // upper bound changing
#ifdef CBC_PRINT2
                if (solver->getColUpper()[k] != newBounds_[i])
                    printf("upper change for column %d - from %g to %g\n", k, solver->getColUpper()[k], newBounds_[i]);
#endif
#ifndef NDEBUG
                if ((variable&0x40000000) == 0 && true) {
                    double oldValue = solver->getColUpper()[k];
                    assert (newBounds_[i] < oldValue + 1.0e-8);
                    if (newBounds_[i] > oldValue - 1.0e-8) {
#ifdef CBC_PRINT2
                        printf("bad null upper change for column %d - bound %g\n", k, oldValue);
#endif
                        if (newBounds_[i] == oldValue)
                            nSame++;
                    }
                }
#endif
                solver->setColUpper(k, newBounds_[i]);
            }
        }
#ifndef NDEBUG
#ifdef CBC_PRINT2
        if (nSame && (nSame < numberChangedBounds_ || (what&3) != 3))
            printf("%d changes out of %d redundant %d\n",
                   nSame, numberChangedBounds_, what);
        else if (numberChangedBounds_ && what == 7 && !nSame)
            printf("%d good changes %d\n",
                   numberChangedBounds_, what);
#endif
#endif
    }
#ifdef JJF_ZERO
    if ((what&2) != 0) {
        OsiClpSolverInterface * clpSolver
        = dynamic_cast<OsiClpSolverInterface *> (solver);
        assert (clpSolver);
        //assert (clpSolver->getNumRows()==numberRows_);
        //clpSolver->setBasis(*status_);
        // Current basis
        CoinWarmStartBasis * basis = clpSolver->getPointerToWarmStart();
        printf("BBBB\n");
        basis->print();
        assert (basis->fullBasis());
        basis->applyDiff(status_);
        printf("diff applied %x\n", status_);
        printf("CCCC\n");
        basis->print();
        assert (basis->fullBasis());
#ifndef NDEBUG
        if (!basis->fullBasis())
            printf("Debug this basis!!\n");
#endif
        clpSolver->setBasis(*basis);
    }
#endif
    if ((what&8) != 0) {
        OsiClpSolverInterface * clpSolver
        = dynamic_cast<OsiClpSolverInterface *> (solver);
        assert (clpSolver);
        clpSolver->setBasis(*status_);
        delete status_;
        status_ = NULL;
    }
}

//.........这里部分代码省略.........

/*
  This admittedly looks bogus, but it's the equivalent operation on the matrix
  for inserting a cut of the form -Inf <= +Inf (i.e., a cut with no
  coefficients).
*/
      CoinPackedMatrix exmip1Mtx ;
      exmip1Mtx.reverseOrderedCopyOf(BuildExmip1Mtx()) ;
      CoinPackedVector freeRow ;
      exmip1Mtx.appendRow(freeRow) ;
      OSIUNITTEST_ASSERT_ERROR(exmip1Mtx.isEquivalent(*lhsmbr), {}, "cbc", "matrix by row after assignment");
#endif	// OSICBC_TEST_MTX_STRUCTURE
    }
  }

  // Test add/delete columns
  {    
    OsiCbcSolverInterface m;
    std::string fn = mpsDir+"p0033";
    m.readMps(fn.c_str(),"mps");
    double inf = m.getInfinity();

    CoinPackedVector c0;
    c0.insert(0, 4);
    c0.insert(1, 1);
    m.addCol(c0, 0, inf, 3);
    m.initialSolve();
    double objValue = m.getObjValue();
    CoinRelFltEq eq(1.0e-2);
    OSIUNITTEST_ASSERT_ERROR(eq(objValue,2520.57), {}, "cbc", "objvalue after adding col");

    // Try deleting first column that's nonbasic at lower bound (0).
    int * d = new int[1];
    CoinWarmStartBasis *cwsb = dynamic_cast<CoinWarmStartBasis *>(m.getWarmStart()) ;
    OSIUNITTEST_ASSERT_ERROR(cwsb != NULL, {}, "cbc", "get warmstart basis");
    CoinWarmStartBasis::Status stati ;
    int iCol ;
    for (iCol = 0 ;  iCol < cwsb->getNumStructural() ; iCol++)
    { stati = cwsb->getStructStatus(iCol) ;
      if (stati == CoinWarmStartBasis::atLowerBound) break ; }
    d[0]=iCol;
    m.deleteCols(1,d);
    delete [] d;
    delete cwsb;
    d=NULL;
    m.resolve();
    objValue = m.getObjValue();
    OSIUNITTEST_ASSERT_ERROR(eq(objValue,2520.57), {}, "clp", "objvalue after deleting first col");

    // Try deleting column we added. If basic, go to initialSolve as deleting
    // basic variable trashes basis required for warm start.
    iCol = m.getNumCols()-1;
    cwsb = dynamic_cast<CoinWarmStartBasis *>(m.getWarmStart()) ;
    stati =  cwsb->getStructStatus(iCol) ;
    delete cwsb;
    m.deleteCols(1,&iCol);
    if (stati == CoinWarmStartBasis::basic)
    { m.initialSolve() ; }
    else
    { m.resolve(); }
    objValue = m.getObjValue();
    OSIUNITTEST_ASSERT_ERROR(eq(objValue,2520.57), {}, "clp", "objvalue after deleting added col");
  }

  // Build a model
  {    

//.........这里部分代码省略.........
                        std::cout << "x[40] = " << saveSolution[40] 
                                  << std::endl;
                    }
#endif

                    intObject = NULL;
                }
                else {
                    // TODO: currently all are integer objects.
#ifdef BLIS_DEBU
                    assert(0);
#endif
                }
                
            }
        }
            
        if (numInfs) {
#if 0
            std::cout << "PSEUDO: numInfs = " << numInfs
                      << std::endl;
#endif
            break;
        }
        else if (pass == 0) {
            // The first pass and is IP feasible.
            
#if 1
            std::cout << "ERROR: PSEUDO: given a integer feasible sol, no fraction variable" << std::endl;
            assert(0);
#endif      
            
            roundAgain = false;
            CoinWarmStartBasis * ws = 
                dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
            if (!ws) break;
            
            // Force solution values within bounds
            for (i = 0; i < numCols; ++i) {
                lpX = saveSolution[i];
                if (lpX < lower[i]) {
                    saveSolution[i] = lower[i];
                    roundAgain = true;
                    ws->setStructStatus(i, CoinWarmStartBasis::atLowerBound);
                } 
                else if (lpX > upper[i]) {
                    saveSolution[i] = upper[i];
                    roundAgain = true;
                    ws->setStructStatus(i, CoinWarmStartBasis::atUpperBound);
                } 
            }
            
            if (roundAgain) {
                // Need resolve and do the second round selection.
                solver->setWarmStart(ws);
                delete ws;
                
                // Resolve.
                solver->resolve();
		
                if (!solver->isProvenOptimal()) {
                    // Become infeasible, can do nothing. 
                    bStatus = -2;
                    goto TERM_CREATE;
                }
                else {

//.........这里部分代码省略.........
# endif
	if (basicCnt>prob.nrows_) {
	  // Take out slacks
	  double * acts = prob.acts_;
	  double * rlo = prob.rlo_;
	  double * rup = prob.rup_;
	  double infinity = si.getInfinity();
	  for (i = 0 ; i < prob.nrows_ ; i++) {
	    status = prob.getRowStatus(i);
	    if (status == CoinPrePostsolveMatrix::basic) {
	      basicCnt-- ;
	      double down = acts[i]-rlo[i];
	      double up = rup[i]-acts[i];
	      if (CoinMin(up,down)<infinity) {
		if (down<=up)
		  prob.setRowStatus(i,CoinPrePostsolveMatrix::atLowerBound);
		else
		  prob.setRowStatus(i,CoinPrePostsolveMatrix::atUpperBound);
	      } else {
		prob.setRowStatus(i,CoinPrePostsolveMatrix::isFree);
	      }
	    }
	    if (basicCnt==prob.nrows_)
	      break;
	  }
	}
      }
#endif

/*
  Install the status and primal solution, if we've been carrying them along.

  The code that copies status is efficient but brittle. The current definitions
  for CoinWarmStartBasis::Status and CoinPrePostsolveMatrix::Status are in
  one-to-one correspondence. This code will fail if that ever changes.
*/
      if (doStatus) {
	presolvedModel_->setColSolution(prob.sol_);
	CoinWarmStartBasis *basis =
	  dynamic_cast<CoinWarmStartBasis *>(presolvedModel_->getEmptyWarmStart());
	basis->resize(prob.nrows_,prob.ncols_);
	int i;
	for (i=0;i<prob.ncols_;i++) {
	  CoinWarmStartBasis::Status status =
	    static_cast<CoinWarmStartBasis::Status> (prob.getColumnStatus(i));
	  basis->setStructStatus(i,status);
	}
	for (i=0;i<prob.nrows_;i++) {
	  CoinWarmStartBasis::Status status =
	    static_cast<CoinWarmStartBasis::Status> (prob.getRowStatus(i));
	  basis->setArtifStatus(i,status);
	}
	presolvedModel_->setWarmStart(basis);
	delete basis ;
	delete [] prob.sol_;
	delete [] prob.acts_;
	delete [] prob.colstat_;
	prob.sol_=NULL;
	prob.acts_=NULL;
	prob.colstat_=NULL;
      }
/*
  Copy original column and row information from the CoinPresolveMatrix object
  so it'll be available for postsolve.
*/
      int ncolsNow = presolvedModel_->getNumCols();

//.........这里部分代码省略.........
                    if (intObject->columnIndex() == 15) {
                        std::cout << "x[15] = " << saveSolution[15]
                                  << std::endl;
                    }
#endif

                    intObject = NULL;
                }
                else {
                    // TODO: currently all are integer objects.
#ifdef BLIS_DEBU
                    assert(0);
#endif
                }

            }
        }

        if (numInfs) {
#ifdef BLIS_DEBUG_MORE
            std::cout << "REL: numInfs = " << numInfs
                      << std::endl;
#endif
            break;
        }
        else if (pass == 0) {
            // The first pass and is IP feasible.

#ifdef BLIS_DEBUG
            std::cout << "REL: given a feasible sol" << std::endl;
#endif

            roundAgain = false;
            CoinWarmStartBasis * ws =
                dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
            if (!ws) break;

            // Force solution values within bounds
            for (i = 0; i < numCols; ++i) {
                lpX = saveSolution[i];
                if (lpX < lower[i]) {
                    saveSolution[i] = lower[i];
                    roundAgain = true;
                    ws->setStructStatus(i, CoinWarmStartBasis::atLowerBound);
                }
                else if (lpX > upper[i]) {
                    saveSolution[i] = upper[i];
                    roundAgain = true;
                    ws->setStructStatus(i, CoinWarmStartBasis::atUpperBound);
                }
            }

            if (roundAgain) {
                // Need resolve and do the second round selection.
                solver->setWarmStart(ws);
                delete ws;

                // Resolve.
                solver->resolve();

                if (!solver->isProvenOptimal()) {
                    // Become infeasible, can do nothing.
                    bStatus = -2;
                    goto TERM_CREATE;
                }
                else {

void
CbcTreeArray::cleanTree(CbcModel * model, double cutoff, double & bestPossibleObjective)
{
    int j;
    int nNodes = size();
    int lastNode = nNodes + 1;
    CbcNode ** nodeArray = new CbcNode * [lastNode];
    int * depth = new int [lastNode];
    int k = 0;
    int kDelete = lastNode;
    bestPossibleObjective = 1.0e100 ;
    /*
        Destructively scan the heap. Nodes to be retained go into the front of
        nodeArray, nodes to be deleted into the back. Store the depth in a
        correlated array for nodes to be deleted.
    */
    for (j = 0; j < nNodes; j++) {
        CbcNode * node = nodes_.front();
        nodes_.front()->setOnTree(false);
        std::pop_heap(nodes_.begin(), nodes_.end(), comparison_);
        nodes_.pop_back();
        double value = node ? node->objectiveValue() : COIN_DBL_MAX;
        if (node && value >= cutoff) {
            // double check in case node can change its mind!
            value = node->checkIsCutoff(cutoff);
        }
        if (value >= cutoff || !node->active()) {
            if (node) {
                nodeArray[--kDelete] = node;
                depth[kDelete] = node->depth();
            }
        } else {
            bestPossibleObjective = CoinMin(bestPossibleObjective, value);
            nodeArray[k++] = node;
        }
    }
    if (lastNode_) {
        double value = lastNode_->objectiveValue();
        bestPossibleObjective = CoinMin(bestPossibleObjective, value);
        if (value >= cutoff || !lastNode_->active()) {
            nodeArray[--kDelete] = lastNode_;
            depth[kDelete] = lastNode_->depth();
            lastNode_ = NULL;
        }
    }
    CbcCompareDefault * compareDefault
    = dynamic_cast<CbcCompareDefault *> (comparison_.test_);
    assert (compareDefault);
    compareDefault->setBestPossible(bestPossibleObjective);
    compareDefault->setCutoff(cutoff);
    /*
      Rebuild the heap using the retained nodes.
    */
    for (j = 0; j < k; j++) {
        CbcNode * node = nodeArray[j];
        node->setOnTree(true);
        nodes_.push_back(node);
        std::push_heap(nodes_.begin(), nodes_.end(), comparison_);
    }
    /*
      Sort the list of nodes to be deleted, nondecreasing.
    */
    CoinSort_2(depth + kDelete, depth + lastNode, nodeArray + kDelete);
    /*
      Work back from deepest to shallowest. In spite of the name, addCuts1 is
      just a preparatory step. When it returns, the following will be true:
        * all cuts are removed from the solver's copy of the constraint system;
        * lastws will be a basis appropriate for the specified node;
        * variable bounds will be adjusted to be appropriate for the specified
          node;
        * addedCuts_ (returned via addedCuts()) will contain a list of cuts that
          should be added to the constraint system at this node (but they have
          not actually been added).
      Then we scan the cut list for the node. Decrement the reference count
      for the cut, and if it's gone to 0, really delete it.

      I don't yet see why the checks for status != basic and addedCuts_[i] != 0
      are necessary. When reconstructing a node, these checks are used to skip
      over loose cuts, excluding them from the reconstituted basis. But here
      we're just interested in correcting the reference count. Tight/loose
      should make no difference.

      Arguably a separate routine should be used in place of addCuts1. It's
      doing more work than needed, modifying the model to match a subproblem
      at a node that will be discarded.  Then again, we seem to need the basis.
    */
    for (j = lastNode - 1; j >= kDelete; j--) {
        CbcNode * node = nodeArray[j];
        CoinWarmStartBasis *lastws = model->getEmptyBasis() ;

        model->addCuts1(node, lastws);
        // Decrement cut counts
        assert (node);
        //assert (node->nodeInfo());
        int numberLeft = (node->nodeInfo()) ? node->nodeInfo()->numberBranchesLeft() : 0;
        int i;
        for (i = 0; i < model->currentNumberCuts(); i++) {
            // take off node
            CoinWarmStartBasis::Status status =
                lastws->getArtifStatus(i + model->numberRowsAtContinuous());
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      OSIUNITTEST_ASSERT_ERROR(ei[ 2] == 3, {}, "glpk", "matrix by row after assignment: indices");
      OSIUNITTEST_ASSERT_ERROR(ei[ 1] == 4, {}, "glpk", "matrix by row after assignment: indices");
      OSIUNITTEST_ASSERT_ERROR(ei[ 0] == 7, {}, "glpk", "matrix by row after assignment: indices");
      OSIUNITTEST_ASSERT_ERROR(ei[ 6] == 1, {}, "glpk", "matrix by row after assignment: indices");
      OSIUNITTEST_ASSERT_ERROR(ei[ 5] == 2, {}, "glpk", "matrix by row after assignment: indices");
      OSIUNITTEST_ASSERT_ERROR(ei[ 8] == 2, {}, "glpk", "matrix by row after assignment: indices");
      OSIUNITTEST_ASSERT_ERROR(ei[ 7] == 5, {}, "glpk", "matrix by row after assignment: indices");
      OSIUNITTEST_ASSERT_ERROR(ei[10] == 3, {}, "glpk", "matrix by row after assignment: indices");
      OSIUNITTEST_ASSERT_ERROR(ei[ 9] == 6, {}, "glpk", "matrix by row after assignment: indices");
      OSIUNITTEST_ASSERT_ERROR(ei[13] == 0, {}, "glpk", "matrix by row after assignment: indices");
      OSIUNITTEST_ASSERT_ERROR(ei[12] == 4, {}, "glpk", "matrix by row after assignment: indices");
      OSIUNITTEST_ASSERT_ERROR(ei[11] == 7, {}, "glpk", "matrix by row after assignment: indices");
    }
    
  }

  // Test add/delete columns
  {    
    OsiGlpkSolverInterface m;
    std::string fn = mpsDir+"p0033";
    m.readMps(fn.c_str(),"mps");
    double inf = m.getInfinity();

    CoinPackedVector c0;
    c0.insert(0, 4);
    c0.insert(1, 1);
    m.addCol(c0, 0, inf, 3);
    m.initialSolve();
    double objValue = m.getObjValue();
    CoinRelFltEq eq(1.0e-2);
    OSIUNITTEST_ASSERT_ERROR(eq(objValue,2520.57), {}, "glpk", "add/delete columns: first optimal value");
    // Try deleting first column that's nonbasic at lower bound (0).
    int * d = new int[1];
    CoinWarmStartBasis *cwsb = dynamic_cast<CoinWarmStartBasis *>(m.getWarmStart()) ;
    OSIUNITTEST_ASSERT_ERROR(cwsb != NULL, {}, "glpk", "add/delete columns: have warm start basis");
    CoinWarmStartBasis::Status stati ;
    int iCol ;
    for (iCol = 0 ;  iCol < cwsb->getNumStructural() ; iCol++)
    { stati = cwsb->getStructStatus(iCol) ;
      if (stati == CoinWarmStartBasis::atLowerBound) break ; }
    d[0]=iCol;
    m.deleteCols(1,d);
    delete [] d;
    delete cwsb;
    d=NULL;
    m.resolve();
    objValue = m.getObjValue();
    OSIUNITTEST_ASSERT_ERROR(eq(objValue,2520.57), {}, "glpk", "add/delete columns: optimal value after deleting nonbasic column");
    // Try deleting column we added. If basic, go to initialSolve as deleting
    // basic variable trashes basis required for warm start.
    iCol = m.getNumCols()-1;
    cwsb = dynamic_cast<CoinWarmStartBasis *>(m.getWarmStart()) ;
    stati =  cwsb->getStructStatus(iCol) ;
    delete cwsb;
    m.deleteCols(1,&iCol);
    if (stati == CoinWarmStartBasis::basic)
    { m.initialSolve() ; }
    else
    { m.resolve(); }
    objValue = m.getObjValue();
    OSIUNITTEST_ASSERT_ERROR(eq(objValue,2520.57), {}, "glpk", "add/delete columns: optimal value after deleting added column");
  }

#if 0
  // ??? Simplex routines not adapted to OsiGlpk yet

//--------------------------------------------------------------------------
// ** At present this does not use any solver
void
CglGomoryUnitTest(
  const OsiSolverInterface * baseSiP,
  const std::string mpsDir )
{
  CoinRelFltEq eq(0.000001);

  // Test default constructor
  {
    CglGomory aGenerator;
    assert (aGenerator.getLimit()==50);
    assert (aGenerator.getAway()==0.05);
  }
  
  // Test copy & assignment etc
  {
    CglGomory rhs;
    {
      CglGomory bGenerator;
      bGenerator.setLimit(99);
      bGenerator.setAway(0.2);
      CglGomory cGenerator(bGenerator);
      rhs=bGenerator;
      assert (rhs.getLimit()==99);
      assert (rhs.getAway()==0.2);
    }
  }

  // Test explicit form - all integer (pg 125 Wolsey)
  if (1) {
    OsiCuts osicuts;
    CglGomory test1;
    int i;
    int nOldCuts=0,nRowCuts;
 
    // 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 (objective not used just yet)
    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};
  
    // integer
    char intVar[7]={2,2,2,2,2,2,2};

    // basis 1
    int rowBasis1[3]={-1,-1,-1};
    int colBasis1[5]={1,1,-1,-1,1};
    CoinWarmStartBasis warm;
    warm.setSize(5,3);
    for (i=0;i<3;i++) {
      if (rowBasis1[i]<0) {
	warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound);
      } else {
	warm.setArtifStatus(i,CoinWarmStartBasis::basic);
      }
    }
    for (i=0;i<5;i++) {
      if (colBasis1[i]<0) {
	warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound);
      } else {
	warm.setStructStatus(i,CoinWarmStartBasis::basic);
      }
    }

    // solution 1
    double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0};
    test1.generateCuts(NULL, osicuts, matrix,
		       /*objective,*/ colsol1,
		 colLower, colUpper,
		 rowLower, rowUpper, intVar, &warm);
    nRowCuts = osicuts.sizeRowCuts();
    std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl;
    assert (nRowCuts==2);
    // cuts always <=
    int testCut=0; // test first cut as stronger
    double rhs=-6.0;
    double testCut1[5]={0.0,0.0,-1.0,-2.0,0.0};
    double * cut = testCut1;
    double * colsol = colsol1;
    for (i=nOldCuts; i<nRowCuts; i++){
      OsiRowCut rcut;
      CoinPackedVector rpv;
      rcut = osicuts.rowCut(i);
      rpv = rcut.row();
      const int n = rpv.getNumElements();
      const int * indices = rpv.getIndices();
      double* elements = rpv.getElements();
      double sum2=0.0;
      int k=0;
      for (k=0; k<n; k++){
//.........这里部分代码省略.........

void OsiIF::_loadBasis(
	Array<LPVARSTAT::STATUS> &lpVarStat,
	Array<SlackStat::STATUS> &slackStat)
{
	int lps = lpVarStat.size();
	int sls = slackStat.size();

	CoinWarmStartBasis *ws = nullptr;
	ws = new CoinWarmStartBasis();
	ws->setSize(numCols_, numRows_);

	if (osiLP_->getNumCols() > lps) {
		Logger::ifout() << "OsiIF::_loadBasis: mismatch in number of columns: OSI " << osiLP_->getNumCols() << ", Abacus: " << lps << "\n";
		OGDF_THROW_PARAM(AlgorithmFailureException, ogdf::afcOsiIf);
	}
	for (int i = 0; i < numCols_; i++)
		ws->setStructStatus(i, lpVarStat2osi(lpVarStat[i]));

	if (osiLP_->getNumRows() > sls) {
		Logger::ifout() << "OsiIF::_loadBasis: mismatch in number of rows: OSI " << osiLP_->getNumCols() << ", Abacus: " << sls << "\n";
		OGDF_THROW_PARAM(AlgorithmFailureException, ogdf::afcOsiIf);
	}
	for (int i = 0; i < numRows_; i++)
		ws->setArtifStatus(i, slackStat2osi(slackStat[i]));

	lpSolverTime_.start();
	slackStatus_ = basisStatus_ = Missing;
	int status = 0;
	// FIXME loadBasis
	// better test whether the number of basic structurals is correct?
	if (ws->numberBasicStructurals() > 0) {
		status = osiLP_->setWarmStart(dynamic_cast<CoinWarmStart *> (ws));
		if (ws_ != nullptr) delete ws_;
		ws_ = dynamic_cast<CoinWarmStartBasis*> (osiLP_->getWarmStart());
		if (ws_ != nullptr) {
			delete[] cStat_;
			int nStructBytes = (int) ceil( ws_->getNumStructural() / 4.0);
			cStat_ = new char[nStructBytes];
			for(int i = 0; i < nStructBytes; i++) {
				cStat_[i] = ws_->getStructuralStatus()[i];
			}

			delete[] rStat_;
			int nArtBytes = (int) ceil( ws_->getNumArtificial() / 4.0 );
			rStat_ = new char[nArtBytes];
			for(int i = 0; i < nArtBytes; i++) {
				rStat_[i] = ws_->getArtificialStatus()[i];
			}

			basisStatus_ = Available;
		} else
			basisStatus_ = Missing;
	} else
		status = 2;
	lpSolverTime_.stop();

	delete ws;

	if (status == 0) {
		Logger::ifout()
			<< "OsiIF::_loadBasis(): loading the new basis has failed. Status "
			<< status << endl;
		// FIXME loadBasis
		return;
	} else
		return;
}