C++ CoinPackedMatrix::appendRow方法代码示例

void OsiIF::loadDummyRow(OsiSolverInterface* s2, const double* lbounds, const double* ubounds, const double* objectives)
{
	CoinPackedVector *coinrow = new CoinPackedVector();
	CoinPackedMatrix *matrix =  new CoinPackedMatrix(false,0,0);
	matrix->setDimensions(0, numCols_);
	ArrayBuffer<int> dummy(1,false);
	dummy.push(0);
	char *senses = new char[1];
	double *rhs = new double[1];
	double *ranges = new double[1];
	coinrow->insert(0, 1.);
	matrix->appendRow(*coinrow);
	senses[0] = 'E';
	rhs[0] = 1.;
	ranges[0] = 0.;

	lpSolverTime_.start();
	s2->loadProblem(*matrix, lbounds, ubounds, objectives, senses, rhs, ranges);
	lpSolverTime_.stop();
	_remRows(dummy);

	delete coinrow;
	delete matrix;
	freeChar(senses);
	freeDouble(rhs);
	freeDouble(ranges);

	return;
}

ProblemStatus ColaModel::solve_numeric() {
  solve(false);
  // 1. determine original binding constraints
  std::vector<int> binding;
  int num_rows = getNumRows();
  int orig_num_rows = num_rows - total_num_cuts_;
  int num_cols = getNumCols();
  int * cstat = new int[num_cols];
  int * rstat = new int[num_rows];
  getBasisStatus(cstat, rstat);
  for (int i=0; i<orig_num_rows; ++i) {
    if (rstat[i]==2 or rstat[i]==3) {
      binding.push_back(i);
    }
  }
  int num_binding = binding.size();
  delete[] cstat;
  delete[] rstat;
  // 2. compute AA^T
  // A is sparse and instance of CoinPackedMatrix
  // 2.1 get A
  CoinPackedMatrix const * mat = getMatrixByRow();
  CoinPackedMatrix * A = new CoinPackedMatrix(false, 500, 500);
  std::vector<int>::const_iterator it;
  for (it=binding.begin(); it!=binding.end(); ++it) {
    int first = mat->getVectorFirst(*it);
    int last = mat->getVectorLast(*it);
    int num_elem  = last-first;
    int const * mat_cols = mat->getIndices() + first;
    double const * mat_value = mat->getElements() + first;
    A->appendRow(CoinPackedVector(num_elem, mat_cols, mat_value));
  }
  // 2.2 Compute AAt
  CoinPackedMatrix * AAt = new CoinPackedMatrix();
  double * Aa_i = new double[num_binding];
  int * Aa_i_cols = new int[num_binding];
  double * Aa_i_vals = new double[num_binding];
  for (it=binding.begin(); it!=binding.end(); ++it) {
    // A times row i of A
    int first = mat->getVectorFirst(*it);
    int last = mat->getVectorLast(*it);
    int num_elem  = last-first;
    int const * mat_cols = mat->getIndices() + first;
    double const * mat_value = mat->getElements() + first;
    A->times(CoinPackedVector(num_elem, mat_cols, mat_value), Aa_i);
    // sparsify and insert Aa_i
    int Aa_i_size = 0;
    for (int i=0; i<num_binding; ++i) {
      if (Aa_i[i]!=0.0) {
        Aa_i_cols[Aa_i_size] = i;
        Aa_i_vals[Aa_i_size] = Aa_i[i];
        Aa_i_size++;
      }
    }
    AAt->appendCol(CoinPackedVector(Aa_i_size, Aa_i_cols, Aa_i_vals));
  }
  delete[] Aa_i;
  delete[] Aa_i_cols;
  delete[] Aa_i_vals;
  //AAt->dumpMatrix();
  // 3. compute Ac
  double * Ac = new double[num_binding];
  double const * obj = getObjCoefficients();
  int obj_size;
  int * obj_cols = new int[num_cols];
  double * obj_vals = new double[num_cols];
  for (int i=0; i<num_cols; ++i) {
    if (obj[i]!=0) {
      obj_cols[obj_size] = i;
      obj_vals[obj_size] = obj[i];
      obj_size++;
    }
  }
  A->times(CoinPackedVector(obj_size, obj_cols, obj_vals), Ac);
  delete[] obj_cols;
  delete[] obj_vals;
  // 4. compute (b-Ac)
  double * b = new double[num_binding];
  int k=0;
  for (it=binding.begin(); it!=binding.end(); ++it) {
    b[k] = getRightHandSide()[*it];
    k++;
  }
  double * b_Ac = new double[num_binding];
  for (int i=0; i<num_binding; ++i) {
    b_Ac[i] = b[i] - Ac[i];
  }
  // 5. solve AA^Ty=(b-Ac)
  // 5.1 get AAt in lower triangular format
  double ** AAt_dense = new double*[num_binding];
  for (int i=0; i<num_binding; ++i) {
    AAt_dense[i] = new double[num_binding]();
  }
  int const * AAt_cols = AAt->getIndices();
  double const * AAt_value = AAt->getElements();
  for (int i=0; i<num_binding; ++i) {
    // get row i
    int first = AAt->getVectorFirst(i);
    int last = AAt->getVectorLast(i);
    //int num_elem  = last-first;
//.........这里部分代码省略.........

bool OSI_X_SolverWrapper<SOLVERINTERFACE>::setup(const LP_Constraints & cstraints) //cstraints <-> constraints
{
  bool bOk = true;
  if ( si == NULL )
  {
    return false;
  }
  assert(_nbParams == cstraints._nbParams);

  
  int NUMVAR = cstraints._constraintMat.cols();
  int NUMCON = cstraints._constraintMat.rows();
  int NUMANZ = cstraints._constraintMat.cols() * cstraints._constraintMat.rows(); //DENSE MATRIX
  
  std::vector<double> col_lb(NUMVAR);//the column lower bounds
  std::vector<double> col_ub(NUMVAR);//the column upper bounds
  
  this->_nbParams = NUMVAR;

  if (cstraints._bminimize) 
  {
    si->setObjSense( 1 );
  }
  else  
  {
    si->setObjSense( -1 );
  }

  const Mat & A = cstraints._constraintMat;

  //Equality constraint will be handked by two constraintsdue to the API limitation. 
  size_t nbLine = A.rows() + std::count(cstraints._vec_sign.begin(), cstraints._vec_sign.end(), EQ);
  
  std::vector<double> row_lb(nbLine);//the row lower bounds
  std::vector<double> row_ub(nbLine);//the row upper bounds
  
  CoinPackedMatrix * matrix = new CoinPackedMatrix(false,0,0);
  matrix->setDimensions(0, NUMVAR);

  
  //-- Add row-wise constraint
  size_t indexRow = 0;
  for (int i=0; i < A.rows(); ++i) 
  {
    Vec temp = A.row(i);
    
    CoinPackedVector row;
    if ( cstraints._vec_sign[i] == EQ || cstraints._vec_sign[i] == LE )
    {
      int coef = 1;
      for ( int j = 0; j < A.cols() ; j++ )
      {
	row.insert(j, coef * temp.data()[j]);
      }
      row_lb[indexRow] = -1.0 * si->getInfinity();
      row_ub[indexRow] = coef * cstraints._Cst_objective(i);
      matrix->appendRow(row);
      indexRow++;
    }
    if ( cstraints._vec_sign[i] == EQ || cstraints._vec_sign[i] == GE )
    {
      int coef = -1;
      for ( int j = 0; j < A.cols() ; j++ )
      {
	      row.insert(j, coef * temp.data()[j]);
      }
      row_lb[indexRow] = -1.0 * si->getInfinity();
      row_ub[indexRow] = coef * cstraints._Cst_objective(i);
      matrix->appendRow(row);      
      indexRow++;
    }
  }
  
  //-- Setup bounds
  if (cstraints._vec_bounds.size() == 1) 
  {
    // Setup the same bound for all the parameter
    for (int i=0; i < this->_nbParams; ++i)
    {
      col_lb[i] = cstraints._vec_bounds[0].first;
      col_ub[i] = cstraints._vec_bounds[0].second;
    }
  }
  else  
  {

    for (int i=0; i < this->_nbParams; ++i)
    {      
      col_lb[i] = cstraints._vec_bounds[i].first;
      col_ub[i] = cstraints._vec_bounds[i].second;
    }
  }  
 
  si->loadProblem(*matrix, &col_lb[0], &col_ub[0], cstraints._vec_cost.empty() ? NULL : &cstraints._vec_cost[0], &row_lb[0], &row_ub[0] );
  
  delete matrix;
  
  return bOk;
}

int
main(void)
{
   // Create a problem pointer.  We use the base class here.
   OsiSolverInterface *si;

   // When we instantiate the object, we need a specific derived class.
   si = new OSIXXX;

   // Build our own instance from scratch

   /*
    * This section adapted from Matt Galati's example 
    * on the COIN-OR Tutorial website.
    *
    * Problem from Bertsimas, Tsitsiklis page 21
    *  
    *  optimal solution: x* = (1,1)
    *  
    *  minimize -1 x0 - 1 x1
    *  s.t       1 x0 + 2 x1 <= 3
    *            2 x0 + 1 x1 <= 3
    *              x0        >= 0
    *              x1        >= 0
    */

   int n_cols = 2;
   double *objective    = new double[n_cols];//the objective coefficients
   double *col_lb       = new double[n_cols];//the column lower bounds
   double *col_ub       = new double[n_cols];//the column upper bounds

   //Define the objective coefficients.
   //minimize -1 x0 - 1 x1
   objective[0] = -1.0;
   objective[1] = -1.0;

   //Define the variable lower/upper bounds.
   // x0 >= 0   =>  0 <= x0 <= infinity
   // x1 >= 0   =>  0 <= x1 <= infinity
   col_lb[0] = 0.0;
   col_lb[1] = 0.0;
   col_ub[0] = si->getInfinity();
   col_ub[1] = si->getInfinity();
     
   int n_rows = 2;
   double *row_lb = new double[n_rows]; //the row lower bounds
   double *row_ub = new double[n_rows]; //the row upper bounds
     
   //Define the constraint matrix.
   CoinPackedMatrix *matrix =  new CoinPackedMatrix(false,0,0);
   matrix->setDimensions(0, n_cols);

   //1 x0 + 2 x1 <= 3  =>  -infinity <= 1 x0 + 2 x2 <= 3
   CoinPackedVector row1;
   row1.insert(0, 1.0);
   row1.insert(1, 2.0);
   row_lb[0] = -1.0 * si->getInfinity();
   row_ub[0] = 3.0;
   matrix->appendRow(row1);

   //2 x0 + 1 x1 <= 3  =>  -infinity <= 2 x0 + 1 x1 <= 3
   CoinPackedVector row2;
   row2.insert(0, 2.0);
   row2.insert(1, 1.0);
   row_lb[1] = -1.0 * si->getInfinity();
   row_ub[1] = 3.0;
   matrix->appendRow(row2);

   //load the problem to OSI
   si->loadProblem(*matrix, col_lb, col_ub, objective, row_lb, row_ub);

   //write the MPS file to a file called example.mps
   si->writeMps("example");

  

   // Solve the (relaxation of the) problem
   si->initialSolve();

   // Check the solution
   if ( si->isProvenOptimal() ) { 
      std::cout << "Found optimal solution!" << std::endl; 
      std::cout << "Objective value is " << si->getObjValue() << std::endl;

      int n = si->getNumCols();
      const double *solution;
      solution = si->getColSolution();
      // We could then print the solution or examine it.
   } else {
      std::cout << "Didn't find optimal solution." << std::endl;
      // Could then check other status functions.
   }

   return 0;
}

//.........这里部分代码省略.........
     //CoinZeroN(objCoeffs, lCols);
     
     int intCnt(0);
  
     /** Fill in array of lower-level integer variables **/
     
     for(i = 0; i < lCols; i++){
	index1 = lColIndices[i];
	if(oSolver->isInteger(index1)){
	   integerVars[intCnt] = i;
	   intCnt++;
	}
     }
     
     CoinDisjointCopyN(lObjCoeffs, lCols, objCoeffs);
     
     CoinPackedMatrix * newMat = new CoinPackedMatrix(false, 0, 0);
     newMat->setDimensions(0, lCols);
     double tmp(0.0);
     
     /*
       for(i = 0; i < lRows; i++){
       CoinPackedVector row;
       index1 = lRowIndices[i];
       start = matStarts[index1];
       end = start + matrix->getVectorSize(index1);
       for(j = start; j < end; j++){
       index2 = matIndices[j];
       //tmp = findIndex(index, lCols, lColIndices);
       tmp = binarySearch(0, lCols - 1, index2, lColIndices);
       if(tmp > -1)
       row.insert(tmp, matElements[j]);
       }
       newMat->appendRow(row);
       }
     */
     for(i = 0; i < lRows; i++){
	CoinPackedVector row;
	index1 = lRowIndices[i];
	for(j = 0; j < lCols; j++){
	   index2 = lColIndices[j];
	   tmp = matrix->getCoefficient(index1, index2);
	   row.insert(j, tmp);
	}
	newMat->appendRow(row);
     }

     /*
       nSolver->assignProblem(newMat, colLb, colUb,
       objCoeffs, rowLb, rowUb);
     */
     
     nSolver->loadProblem(*newMat, colLb, colUb,
			  objCoeffs, rowLb, rowUb);
     
     for(i = 0; i < intCnt; i++){
	nSolver->setInteger(integerVars[i]);
     }
     //nSolver->setInteger(integerVars, intCnt);
     
     nSolver->setObjSense(objSense); //1 min; -1 max
     
     nSolver->setHintParam(OsiDoReducePrint, true, OsiHintDo);

#if 0
     if(0){

void OsiIF::_initialize(
	OptSense sense,
	int nRow,
	int maxRow,
	int nCol,
	int maxCol,
	Array<double> &obj,
	Array<double> &lBound,
	Array<double> &uBound,
	Array<Row*> &rows)
{
	osiLP_ = getDefaultInterface();
	currentSolverType_ = Exact;

	// switch off output from the solver
	// can be reset in setSolverParameters
	osiLP_->setHintParam(OsiDoReducePrint, true, OsiHintDo);
	osiLP_->messageHandler()->setLogLevel(0);
	master_->setSolverParameters(osiLP_, currentSolverType() == Approx);

	numRows_ = nRow;
	numCols_ = nCol;
	double *lbounds = new double[numCols_];
	double *ubounds = new double[numCols_];
	double *objectives = new double[numCols_];

	CoinPackedVector *coinrow = new CoinPackedVector();
	CoinPackedMatrix *matrix =  new CoinPackedMatrix(false,0,0);
	matrix->setDimensions(0, numCols_);

	for (int i = 0; i < numCols_; i++){
		lbounds[i] = lBound[i];
		ubounds[i] = uBound[i];
		objectives[i] = obj[i];
	}

	if (currentSolverType() == Exact && numRows_ == 0 && master_->defaultLpSolver() == Master::CPLEX) {
		loadDummyRow(osiLP_, lbounds, ubounds, objectives);
	}
	else {
		char *senses = new char[numRows_];
		double *rhs = new double[numRows_];
		double *ranges = new double[numRows_];

		for (int i = 0; i < numRows_; i++){
			coinrow->clear();
			for (int j = 0; j < rows[i]->nnz(); j++){
				coinrow->insert(rows[i]->support(j), rows[i]->coeff(j));
			}
			matrix->appendRow(*coinrow);
			senses[i] = csense2osi(rows[i]->sense());
			rhs[i] = rows[i]->rhs();
			ranges[i] = 0.0;
		}
		lpSolverTime_.start();
		osiLP_->loadProblem(*matrix, lbounds, ubounds, objectives, senses, rhs, ranges);
		lpSolverTime_.stop();

		freeChar(senses);
		freeDouble(rhs);
		freeDouble(ranges);
	}

	// set the sense of the optimization
	_sense(sense);

	// get the pointers to the solution, reduced costs etc.
	lpSolverTime_.start();
	numRows_ = osiLP_->getNumRows();
	numCols_ = osiLP_->getNumCols();
	rhs_ = osiLP_->getRightHandSide();
	rowsense_ = osiLP_->getRowSense();
	colupper_ = osiLP_->getColUpper();
	collower_ = osiLP_->getColLower();
	objcoeff_ = osiLP_->getObjCoefficients();
	if( ws_ != nullptr )
		delete ws_;
	//ws_ = dynamic_cast<CoinWarmStartBasis *>(osiLP_->getWarmStart());
	ws_=nullptr;

	xValStatus_ = recoStatus_ = yValStatus_ = slackStatus_ = basisStatus_ = Missing;
	lpSolverTime_.stop();

	delete coinrow;
	delete matrix;
	freeDouble(lbounds);
	freeDouble(ubounds);
	freeDouble(objectives);
}