C++ OsiSolverInterface::getRowActivity方法代码示例

// This looks at solution and sets bounds to contain solution
void 
CbcLink::feasibleRegion()
{
  int j;
  int firstNonZero=-1;
  int lastNonZero = -1;
  OsiSolverInterface * solver = model_->solver();
  const double * solution = model_->testSolution();
  const double * upper = solver->getColUpper();
  double integerTolerance = 
    model_->getDblParam(CbcModel::CbcIntegerTolerance);
  double weight = 0.0;
  double sum =0.0;

  int base=0;
  for (j=0;j<numberMembers_;j++) {
    for (int k=0;k<numberLinks_;k++) {
      int iColumn = which_[base+k];
      double value = CoinMax(0.0,solution[iColumn]);
      sum += value;
      if (value>integerTolerance&&upper[iColumn]) {
        weight += weights_[j]*value;
        if (firstNonZero<0)
          firstNonZero=j;
        lastNonZero=j;
      }
    }
    base += numberLinks_;
  }
#ifdef DISTANCE
  if (lastNonZero-firstNonZero>sosType_-1) {
    /* may still be satisfied.
       For LOS type 2 we might wish to move coding around
       and keep initial info in model_ for speed
    */
    int iWhere;
    bool possible=false;
    for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) {
      if (fabs(weight-weights_[iWhere])<1.0e-8) {
	possible=true;
	break;
      }
    }
    if (possible) {
      // One could move some of this (+ arrays) into model_
      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
      const double * element = matrix->getMutableElements();
      const int * row = matrix->getIndices();
      const CoinBigIndex * columnStart = matrix->getVectorStarts();
      const int * columnLength = matrix->getVectorLengths();
      const double * rowSolution = solver->getRowActivity();
      const double * rowLower = solver->getRowLower();
      const double * rowUpper = solver->getRowUpper();
      int numberRows = matrix->getNumRows();
      double * array = new double [numberRows];
      CoinZeroN(array,numberRows);
      int * which = new int [numberRows];
      int n=0;
      int base=numberLinks_*firstNonZero;
      for (j=firstNonZero;j<=lastNonZero;j++) {
	for (int k=0;k<numberLinks_;k++) {
	  int iColumn = which_[base+k];
	  double value = CoinMax(0.0,solution[iColumn]);
	  if (value>integerTolerance&&upper[iColumn]) {
	    value = CoinMin(value,upper[iColumn]);
	    for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
	      int iRow = row[j];
	      double a = array[iRow];
	      if (a) {
		a += value*element[j];
		if (!a)
		  a = 1.0e-100;
	      } else {
		which[n++]=iRow;
		a=value*element[j];
		assert (a);
	      }
	      array[iRow]=a;
	    }
	  }
	}
	base += numberLinks_;
      }
      base=numberLinks_*iWhere;
      for (int k=0;k<numberLinks_;k++) {
	int iColumn = which_[base+k];
	const double value = 1.0;
	for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
	  int iRow = row[j];
	  double a = array[iRow];
	  if (a) {
	    a -= value*element[j];
	    if (!a)
	      a = 1.0e-100;
	  } else {
	    which[n++]=iRow;
	    a=-value*element[j];
	    assert (a);
	  }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    //   int numRows = m->getNumRows();
    CoinFillN(integers_ ,n, true);
    for (int i = 0 ;  i < numCols ; i++)
    {
        if (si.isContinuous(i))
            integers_[i] = false;
    }
    bool * integerSlacks = integers_ + numCols;
    for (int i = 0 ; i < nBasics ; i++)
    {
        if (rowLower[i] > -1e50 && INT_INFEAS(rowLower[i]) > 1e-15)
            integerSlacks[i] = false;
        if (rowUpper[i] < 1e50 && INT_INFEAS(rowUpper[i]) > 1e-15)
            integerSlacks[i] = false;
    }
    for (int i = 0 ;  i < numCols ; i++)
    {
        CoinBigIndex end = starts[i] + lengths[i];
        if (integers_[i])
        {
            for (CoinBigIndex k=starts[i] ; k < end; k++)
            {
                if (integerSlacks[inds[k]] && INT_INFEAS(elems[k])>1e-15 )
                    integerSlacks[inds[k]] = false;
            }
        }
        else
        {
            for (CoinBigIndex k=starts[i] ; k < end; k++)
            {
                if (integerSlacks[inds[k]])
                    integerSlacks[inds[k]] = false;
            }
        }
    }

    CoinCopyN(si.getColSolution(), si.getNumCols(), colsol_);
    CoinCopyN(si.getRowActivity(), si.getNumRows(), slacks_);
    for (int i = 0 ; i < si.getNumRows() ; i++)
    {
        slacks_[i]*=-1;
        if (rowLower[i]>-1e50)
        {
            slacks_[i] += rowLower[i];
        }
        else
        {
            slacks_[i] += rowUpper[i];
        }
    }
    //Now get the fill the arrays;
    nNonBasics = 0;
    nBasics = 0;



    //For having the index variables correctly ordered we need to access to OsiSimplexInterface
    {
        OsiSolverInterface * ncSi = (const_cast<OsiSolverInterface *>(&si));
        ncSi->enableSimplexInterface(0);
        ncSi->getBasics(basics_);
	// Save enabled solver
	solver_ = si.clone();
#ifdef COIN_HAS_OSICLP
	OsiClpSolverInterface * clpSi = dynamic_cast<OsiClpSolverInterface *>(solver_);
	const OsiClpSolverInterface * clpSiRhs = dynamic_cast<const OsiClpSolverInterface *>(&si);
	if (clpSi)
	  clpSi->getModelPtr()->copyEnabledStuff(clpSiRhs->getModelPtr());;
#endif
        ncSi->disableSimplexInterface();
    }

    int numStructural = basis_->getNumStructural();
    for (int i = 0 ; i < numStructural ; i++)
    {
        if (basis_->getStructStatus(i)== CoinWarmStartBasis::basic)
        {
            nBasics++;
            //Basically do nothing
        }
        else
        {
            nonBasics_[nNonBasics++] = i;
        }
    }

    int numArtificial = basis_->getNumArtificial();
    for (int i = 0 ; i < numArtificial ; i++)
    {
        if (basis_->getArtifStatus(i)== CoinWarmStartBasis::basic)
        {
            //Just check number of basics
            nBasics++;
        }
        else
        {
            nonBasics_[nNonBasics++] = i + basis_->getNumStructural();
        }
    }
}

// Infeasibility - large is 0.5
double 
CbcLink::infeasibility(int & preferredWay) const
{
  int j;
  int firstNonZero=-1;
  int lastNonZero = -1;
  OsiSolverInterface * solver = model_->solver();
  const double * solution = model_->testSolution();
  //const double * lower = solver->getColLower();
  const double * upper = solver->getColUpper();
  double integerTolerance = 
    model_->getDblParam(CbcModel::CbcIntegerTolerance);
  double weight = 0.0;
  double sum =0.0;

  // check bounds etc
  double lastWeight=-1.0e100;
  int base=0;
  for (j=0;j<numberMembers_;j++) {
    for (int k=0;k<numberLinks_;k++) {
      int iColumn = which_[base+k];
      //if (lower[iColumn])
      //throw CoinError("Non zero lower bound in CBCLink","infeasibility","CbcLink");
      if (lastWeight>=weights_[j]-1.0e-7)
        throw CoinError("Weights too close together in CBCLink","infeasibility","CbcLink");
      double value = CoinMax(0.0,solution[iColumn]);
      sum += value;
      if (value>integerTolerance&&upper[iColumn]) {
        // Possibly due to scaling a fixed variable might slip through
        if (value>upper[iColumn]+1.0e-8) {
          // Could change to #ifdef CBC_DEBUG
#ifndef NDEBUG
          if (model_->messageHandler()->logLevel()>1)
            printf("** Variable %d (%d) has value %g and upper bound of %g\n",
                   iColumn,j,value,upper[iColumn]);
#endif
        } 
	value = CoinMin(value,upper[iColumn]);
        weight += weights_[j]*value;
        if (firstNonZero<0)
          firstNonZero=j;
        lastNonZero=j;
      }
    }
    base += numberLinks_;
  }
  double valueInfeasibility;
  preferredWay=1;
  if (lastNonZero-firstNonZero>=sosType_) {
    // find where to branch
    assert (sum>0.0);
    weight /= sum;
    valueInfeasibility = lastNonZero-firstNonZero+1;
    valueInfeasibility *= 0.5/((double) numberMembers_);
    //#define DISTANCE
#ifdef DISTANCE
    assert (sosType_==1); // code up
    /* may still be satisfied.
       For LOS type 2 we might wish to move coding around
       and keep initial info in model_ for speed
    */
    int iWhere;
    bool possible=false;
    for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) {
      if (fabs(weight-weights_[iWhere])<1.0e-8) {
	possible=true;
	break;
      }
    }
    if (possible) {
      // One could move some of this (+ arrays) into model_
      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
      const double * element = matrix->getMutableElements();
      const int * row = matrix->getIndices();
      const CoinBigIndex * columnStart = matrix->getVectorStarts();
      const int * columnLength = matrix->getVectorLengths();
      const double * rowSolution = solver->getRowActivity();
      const double * rowLower = solver->getRowLower();
      const double * rowUpper = solver->getRowUpper();
      int numberRows = matrix->getNumRows();
      double * array = new double [numberRows];
      CoinZeroN(array,numberRows);
      int * which = new int [numberRows];
      int n=0;
      int base=numberLinks_*firstNonZero;
      for (j=firstNonZero;j<=lastNonZero;j++) {
	for (int k=0;k<numberLinks_;k++) {
	  int iColumn = which_[base+k];
	  double value = CoinMax(0.0,solution[iColumn]);
	  if (value>integerTolerance&&upper[iColumn]) {
	    value = CoinMin(value,upper[iColumn]);
	    for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
	      int iRow = row[j];
	      double a = array[iRow];
	      if (a) {
		a += value*element[j];
		if (!a)
		  a = 1.0e-100;
	      } else {
//.........这里部分代码省略.........