C++ CoinMax函数代码示例

// Apply bounds
void OsiSolverBranch::applyBounds(OsiSolverInterface &solver, int way) const
{
  int base = way + 1;
  assert(way == -1 || way == 1);
  int numberColumns = solver.getNumCols();
  const double *columnLower = solver.getColLower();
  int i;
  for (i = start_[base]; i < start_[base + 1]; i++) {
    int iColumn = indices_[i];
    if (iColumn < numberColumns) {
      double value = CoinMax(bound_[i], columnLower[iColumn]);
      solver.setColLower(iColumn, value);
    } else {
      int iRow = iColumn - numberColumns;
      const double *rowLower = solver.getRowLower();
      double value = CoinMax(bound_[i], rowLower[iRow]);
      solver.setRowLower(iRow, value);
    }
  }
  const double *columnUpper = solver.getColUpper();
  for (i = start_[base + 1]; i < start_[base + 2]; i++) {
    int iColumn = indices_[i];
    if (iColumn < numberColumns) {
      double value = CoinMin(bound_[i], columnUpper[iColumn]);
      solver.setColUpper(iColumn, value);
    } else {
      int iRow = iColumn - numberColumns;
      const double *rowUpper = solver.getRowUpper();
      double value = CoinMin(bound_[i], rowUpper[iRow]);
      solver.setRowUpper(iRow, value);
    }
  }
}

// Return "down" estimate
double
CbcSimpleIntegerPseudoCost::downEstimate() const
{
    OsiSolverInterface * solver = model_->solver();
    const double * solution = model_->testSolution();
    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();
    double value = solution[columnNumber_];
    value = CoinMax(value, lower[columnNumber_]);
    value = CoinMin(value, upper[columnNumber_]);
    if (upper[columnNumber_] == lower[columnNumber_]) {
        // fixed
        return 0.0;
    }
    double integerTolerance =
        model_->getDblParam(CbcModel::CbcIntegerTolerance);
    double below = floor(value + integerTolerance);
    double above = below + 1.0;
    if (above > upper[columnNumber_]) {
        above = below;
        below = above - 1;
    }
    double downCost = CoinMax((value - below) * downPseudoCost_, 0.0);
    return downCost;
}

double
CbcNWay::infeasibility(const OsiBranchingInformation * /*info*/,
                       int &preferredWay) const
{
    int numberUnsatis = 0;
    int j;
    OsiSolverInterface * solver = model_->solver();
    const double * solution = model_->testSolution();
    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();
    double largestValue = 0.0;

    double integerTolerance =
        model_->getDblParam(CbcModel::CbcIntegerTolerance);

    for (j = 0; j < numberMembers_; j++) {
        int iColumn = members_[j];
        double value = solution[iColumn];
        value = CoinMax(value, lower[iColumn]);
        value = CoinMin(value, upper[iColumn]);
        double distance = CoinMin(value - lower[iColumn], upper[iColumn] - value);
        if (distance > integerTolerance) {
            numberUnsatis++;
            largestValue = CoinMax(distance, largestValue);
        }
    }
    preferredWay = 1;
    if (numberUnsatis) {
        return largestValue;
    } else {
        return 0.0; // satisfied
    }
}

// Updates second array for steepest and does devex weights (need not be coded)
void
ClpMatrixBase::subsetTimes2(const ClpSimplex * model,
                            CoinIndexedVector * dj1,
                            const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
                            double referenceIn, double devex,
                            // Array for exact devex to say what is in reference framework
                            unsigned int * reference,
                            double * weights, double scaleFactor)
{
     // get subset which have nonzero tableau elements
     subsetTransposeTimes(model, pi2, dj1, dj2);
     bool killDjs = (scaleFactor == 0.0);
     if (!scaleFactor)
          scaleFactor = 1.0;
     // columns

     int number = dj1->getNumElements();
     const int * index = dj1->getIndices();
     double * updateBy = dj1->denseVector();
     double * updateBy2 = dj2->denseVector();

     for (int j = 0; j < number; j++) {
          double thisWeight;
          double pivot;
          double pivotSquared;
          int iSequence = index[j];
          double value2 = updateBy[j];
          if (killDjs)
               updateBy[j] = 0.0;
          double modification = updateBy2[j];
          updateBy2[j] = 0.0;
          ClpSimplex::Status status = model->getStatus(iSequence);

          if (status != ClpSimplex::basic && status != ClpSimplex::isFixed) {
               thisWeight = weights[iSequence];
               pivot = value2 * scaleFactor;
               pivotSquared = pivot * pivot;

               thisWeight += pivotSquared * devex + pivot * modification;
               if (thisWeight < DEVEX_TRY_NORM) {
                    if (referenceIn < 0.0) {
                         // steepest
                         thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
                    } else {
                         // exact
                         thisWeight = referenceIn * pivotSquared;
                         if (reference(iSequence))
                              thisWeight += 1.0;
                         thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
                    }
               }
               weights[iSequence] = thisWeight;
          }
     }
     dj2->setNumElements(0);
}

/* Just for debug if odd type matrix.
   Returns number and sum of primal infeasibilities.
*/
int
ClpMatrixBase::checkFeasible(ClpSimplex * model, double & sum) const
{
     int numberRows = model->numberRows();
     double * rhs = new double[numberRows];
     int numberColumns = model->numberColumns();
     int iRow;
     CoinZeroN(rhs, numberRows);
     times(1.0, model->solutionRegion(), rhs, model->rowScale(), model->columnScale());
     int iColumn;
     int logLevel = model->messageHandler()->logLevel();
     int numberInfeasible = 0;
     const double * rowLower = model->lowerRegion(0);
     const double * rowUpper = model->upperRegion(0);
     const double * solution;
     solution = model->solutionRegion(0);
     double tolerance = model->primalTolerance() * 1.01;
     sum = 0.0;
     for (iRow = 0; iRow < numberRows; iRow++) {
          double value = rhs[iRow];
          double value2 = solution[iRow];
          if (logLevel > 3) {
               if (fabs(value - value2) > 1.0e-8)
                    printf("Row %d stored %g, computed %g\n", iRow, value2, value);
          }
          if (value < rowLower[iRow] - tolerance ||
                    value > rowUpper[iRow] + tolerance) {
               numberInfeasible++;
               sum += CoinMax(rowLower[iRow] - value, value - rowUpper[iRow]);
          }
          if (value2 > rowLower[iRow] + tolerance &&
                    value2 < rowUpper[iRow] - tolerance &&
                    model->getRowStatus(iRow) != ClpSimplex::basic) {
               assert (model->getRowStatus(iRow) == ClpSimplex::superBasic);
          }
     }
     const double * columnLower = model->lowerRegion(1);
     const double * columnUpper = model->upperRegion(1);
     solution = model->solutionRegion(1);
     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
          double value = solution[iColumn];
          if (value < columnLower[iColumn] - tolerance ||
                    value > columnUpper[iColumn] + tolerance) {
               numberInfeasible++;
               sum += CoinMax(columnLower[iColumn] - value, value - columnUpper[iColumn]);
          }
          if (value > columnLower[iColumn] + tolerance &&
                    value < columnUpper[iColumn] - tolerance &&
                    model->getColumnStatus(iColumn) != ClpSimplex::basic) {
               assert (model->getColumnStatus(iColumn) == ClpSimplex::superBasic);
          }
     }
     delete [] rhs;
     return numberInfeasible;
}

/** Useful constructor

  Loads actual upper & lower bounds for the specified variable.
*/
CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model,
        int iColumn, double downPseudoCost,
        double upPseudoCost)
        : CbcSimpleInteger(model, iColumn)
{
    downPseudoCost_ = CoinMax(1.0e-10, downPseudoCost);
    upPseudoCost_ = CoinMax(1.0e-10, upPseudoCost);
    breakEven_ = upPseudoCost_ / (upPseudoCost_ + downPseudoCost_);
    upDownSeparator_ = -1.0;
    method_ = 0;
}

inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix,
                        int *low, int *upp, int *fix,
                        CoinDenseVector <double> &b, double *bl, double *bu, double d1, double d2,
                        CoinDenseVector <double> &grad, CoinDenseVector <double> &rL,
                        CoinDenseVector <double> &rU, CoinDenseVector <double> &x,
                        CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
                        CoinDenseVector <double> &y,  CoinDenseVector <double> &z1,
                        CoinDenseVector <double> &z2, CoinDenseVector <double> &r1,
                        CoinDenseVector <double> &r2, double *Pinf, double *Dinf)
{

// Form residuals for the primal and dual equations.
// rL, rU are output, but we input them as full vectors
// initialized (permanently) with any relevant zeros.

// Get some element pointers for efficiency
     double *x_elts  = x.getElements();
     double *r2_elts = r2.getElements();

     for (int k = 0; k < nfix; k++)
          x_elts[fix[k]]  = 0;

     r1.clear();
     r2.clear();
     model->matVecMult( 1, r1, x );
     model->matVecMult( 2, r2, y );
     for (int k = 0; k < nfix; k++)
          r2_elts[fix[k]]  = 0;


     r1      = b    - r1 - d2 * d2 * y;
     r2      = grad - r2 - z1;              // grad includes d1*d1*x
     if (nupp > 0)
          r2    = r2 + z2;

     for (int k = 0; k < nlow; k++)
          rL[low[k]] = bl[low[k]] - x[low[k]] + x1[low[k]];
     for (int k = 0; k < nupp; k++)
          rU[upp[k]] = - bu[upp[k]] + x[upp[k]] + x2[upp[k]];

     double normL = 0.0;
     double normU = 0.0;
     for (int k = 0; k < nlow; k++)
          if (rL[low[k]] > normL) normL = rL[low[k]];
     for (int k = 0; k < nupp; k++)
          if (rU[upp[k]] > normU) normU = rU[upp[k]];

     *Pinf    = CoinMax(normL, normU);
     *Pinf    = CoinMax( r1.infNorm() , *Pinf );
     *Dinf    = r2.infNorm();
     *Pinf    = CoinMax( *Pinf, 1e-99 );
     *Dinf    = CoinMax( *Dinf, 1e-99 );
}

inline void pdxxxresid2(double mu, int nlow, int nupp, int *low, int *upp,
                        CoinDenseVector <double> &cL, CoinDenseVector <double> &cU,
                        CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
                        CoinDenseVector <double> &z1, CoinDenseVector <double> &z2,
                        double *center, double *Cinf, double *Cinf0)
{

// Form residuals for the complementarity equations.
// cL, cU are output, but we input them as full vectors
// initialized (permanently) with any relevant zeros.
// Cinf  is the complementarity residual for X1 z1 = mu e, etc.
// Cinf0 is the same for mu=0 (i.e., for the original problem).

     double maxXz = -1e20;
     double minXz = 1e20;

     double *x1_elts = x1.getElements();
     double *z1_elts = z1.getElements();
     double *cL_elts = cL.getElements();
     for (int k = 0; k < nlow; k++) {
          double x1z1    = x1_elts[low[k]] * z1_elts[low[k]];
          cL_elts[low[k]] = mu - x1z1;
          if (x1z1 > maxXz) maxXz = x1z1;
          if (x1z1 < minXz) minXz = x1z1;
     }

     double *x2_elts = x2.getElements();
     double *z2_elts = z2.getElements();
     double *cU_elts = cU.getElements();
     for (int k = 0; k < nupp; k++) {
          double x2z2    = x2_elts[upp[k]] * z2_elts[upp[k]];
          cU_elts[upp[k]] = mu - x2z2;
          if (x2z2 > maxXz) maxXz = x2z2;
          if (x2z2 < minXz) minXz = x2z2;
     }

     maxXz   = CoinMax( maxXz, 1e-99 );
     minXz   = CoinMax( minXz, 1e-99 );
     *center  = maxXz / minXz;

     double normL = 0.0;
     double normU = 0.0;
     for (int k = 0; k < nlow; k++)
          if (cL_elts[low[k]] > normL) normL = cL_elts[low[k]];
     for (int k = 0; k < nupp; k++)
          if (cU_elts[upp[k]] > normU) normU = cU_elts[upp[k]];
     *Cinf    = CoinMax( normL, normU);
     *Cinf0   = maxXz;
}

double
CbcBranchToFixLots::infeasibility(const OsiBranchingInformation * /*info*/,
                                  int &preferredWay) const
{
    preferredWay = -1;
    CbcNode * node = model_->currentNode();
    int depth;
    if (node)
        depth = CoinMax(node->depth(), 0);
    else
        return 0.0;
    if (depth_ < 0) {
        return 0.0;
    } else if (depth_ > 0) {
        if ((depth % depth_) != 0)
            return 0.0;
    }
    if (djTolerance_ != -1.234567) {
        if (!shallWe())
            return 0.0;
        else
            return 1.0e20;
    } else {
        // See if 3 in same row and sum <FIX_IF_LESS?
        int numberRows = matrixByRow_.getNumRows();
        const double * solution = model_->testSolution();
        const int * column = matrixByRow_.getIndices();
        const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
        const int * rowLength = matrixByRow_.getVectorLengths();
        double bestSum = 1.0;
        int nBest = -1;
        OsiSolverInterface * solver = model_->solver();
        for (int i = 0; i < numberRows; i++) {
            int numberUnsatisfied = 0;
            double sum = 0.0;
            for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
                int iColumn = column[j];
                if (solver->isInteger(iColumn)) {
                    double solValue = solution[iColumn];
                    if (solValue > 1.0e-5 && solValue < FIX_IF_LESS) {
                        numberUnsatisfied++;
                        sum += solValue;
                    }
                }
            }
            if (numberUnsatisfied >= 3 && sum < FIX_IF_LESS) {
                // possible
                if (numberUnsatisfied > nBest ||
                        (numberUnsatisfied == nBest && sum < bestSum)) {
                    nBest = numberUnsatisfied;
                    bestSum = sum;
                }
            }
        }
        if (nBest > 0)
            return 1.0e20;
        else
            return 0.0;
    }
}

void
OsiTestSolverInterface::rowRimResize_(const int newSize)
{
   if (newSize > maxNumrows_) {
      double* rub   = rowupper_;
      double* rlb   = rowlower_;
      char*   sense = rowsense_;
      double* right = rhs_;
      double* range = rowrange_;
      double* dual  = rowprice_;
      double* left  = lhs_;
      maxNumrows_ = CoinMax(1000, (newSize * 5) / 4);
      rowRimAllocator_();
      const int rownum = getNumRows();
      CoinDisjointCopyN(rub  , rownum, rowupper_);
      CoinDisjointCopyN(rlb  , rownum, rowlower_);
      CoinDisjointCopyN(sense, rownum, rowsense_);
      CoinDisjointCopyN(right, rownum, rhs_);
      CoinDisjointCopyN(range, rownum, rowrange_);
      CoinDisjointCopyN(dual , rownum, rowprice_);
      CoinDisjointCopyN(left , rownum, lhs_);
      delete[] rub;
      delete[] rlb;
      delete[] sense;
      delete[] right;
      delete[] range;
      delete[] dual;
      delete[] left;
   }
}

// Creates a branching object
CbcBranchingObject * 
CbcSimpleIntegerFixed::createBranch(OsiSolverInterface * solver,
					    const OsiBranchingInformation * info, int way)  
{
  const double * solution = model_->testSolution();
  const double * lower = solver->getColLower();
  const double * upper = solver->getColUpper();
  double value = solution[columnNumber_];
  value = CoinMax(value, lower[columnNumber_]);
  value = CoinMin(value, upper[columnNumber_]);
  assert (upper[columnNumber_]>lower[columnNumber_]);
  if (!model_->hotstartSolution()) {
    double nearest = floor(value+0.5);
    double integerTolerance = 
    model_->getDblParam(CbcModel::CbcIntegerTolerance);
    if (fabs(value-nearest)<integerTolerance) {
      // adjust value
      if (nearest!=upper[columnNumber_])
	value = nearest+2.0*integerTolerance;
      else
	value = nearest-2.0*integerTolerance;
    }
  } else {
    const double * hotstartSolution = model_->hotstartSolution();
    double targetValue = hotstartSolution[columnNumber_];
    if (way>0)
      value = targetValue-0.1;
    else
      value = targetValue+0.1;
  }
  CbcBranchingObject * branch = new CbcIntegerBranchingObject(model_,columnNumber_,way,
					     value);
  branch->setOriginalObject(this);
  return branch;
}

// Returns true if y better than x
bool 
CbcCompareUser::test (CbcNode * x, CbcNode * y)
{
  if (weight_==-1.0&&(y->depth()>7||x->depth()>7)) {
    // before solution
    /* printf("x %d %d %g, y %d %d %g\n",
       x->numberUnsatisfied(),x->depth(),x->objectiveValue(),
       y->numberUnsatisfied(),y->depth(),y->objectiveValue()); */
    if (x->numberUnsatisfied() > y->numberUnsatisfied()) {
      return true;
    } else if (x->numberUnsatisfied() < y->numberUnsatisfied()) {
      return false;
    } else {
      int testX = x->depth();
      int testY = y->depth();
      if (testX!=testY)
	return testX < testY;
      else
	return equalityTest(x,y); // so ties will be broken in consistent manner
    }
  } else {
    // after solution
    double weight = CoinMax(weight_,0.0);
    double testX =  x->objectiveValue()+ weight*x->numberUnsatisfied();
    double testY = y->objectiveValue() + weight*y->numberUnsatisfied();
    if (testX!=testY)
      return testX > testY;
    else
      return equalityTest(x,y); // so ties will be broken in consistent manner
  }
}

void
CoinPackedVector::append(const CoinPackedVectorBase & caboose)
{
   const int cs = caboose.getNumElements();
   if (cs == 0) {
       return;
   }
   if (testForDuplicateIndex()) {
       // Just to initialize the index heap
       indexSet("append (1st call)", "CoinPackedVector");
   }
   const int s = nElements_;
   // Make sure there is enough room for the caboose
   if ( capacity_ < s + cs)
      reserve(CoinMax(s + cs, 2 * capacity_));

   const int * cind = caboose.getIndices();
   const double * celem = caboose.getElements();
   CoinDisjointCopyN(cind, cs, indices_ + s);
   CoinDisjointCopyN(celem, cs, elements_ + s);
   CoinIotaN(origIndices_ + s, cs, s);
   nElements_ += cs;
   if (testForDuplicateIndex()) {
      std::set<int>& is = *indexSet("append (2nd call)", "CoinPackedVector");
      for (int i = 0; i < cs; ++i) {
	 if (!is.insert(cind[i]).second)
	    throw CoinError("duplicate index", "append", "CoinPackedVector");
      }
   }
}

//-------------------------------------------------------------------
// Generate Stored cuts
//------------------------------------------------------------------- 
void 
CglStoredUser::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
			     const CglTreeInfo info) const
{
  // Get basic problem information
  const double * solution = si.getColSolution();
  if (info.inTree&&info.pass>numberPasses_) {
    // only continue if integer feasible
    int numberColumns=si.getNumCols(); 
    int i;
    const double * colUpper = si.getColUpper();
    const double * colLower = si.getColLower();
    int numberAway=0;
    for (i=0;i<numberColumns;i++) {
      double value = solution[i];
      // In case slightly away from bounds
      value = CoinMax(colLower[i],value);
      value = CoinMin(colUpper[i],value);
      if (si.isInteger(i)&&fabs(value-fabs(value+0.5))>1.0e-5) 
	numberAway++;
    }
    if (numberAway)
      return; // let code branch
  }
  int numberRowCuts = cuts_.sizeRowCuts();
  for (int i=0;i<numberRowCuts;i++) {
    const OsiRowCut * rowCutPointer = cuts_.rowCutPtr(i);
    double violation = rowCutPointer->violated(solution);
    if (violation>=requiredViolation_)
      cs.insert(*rowCutPointer);
  }
}

// Creates a branching object
CbcBranchingObject * 
CbcLink::createBranch(int way) 
{
  int j;
  const double * solution = model_->testSolution();
  double integerTolerance = 
      model_->getDblParam(CbcModel::CbcIntegerTolerance);
  OsiSolverInterface * solver = model_->solver();
  const double * upper = solver->getColUpper();
  int firstNonFixed=-1;
  int lastNonFixed=-1;
  int firstNonZero=-1;
  int lastNonZero = -1;
  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];
      if (upper[iColumn]) {
        double value = CoinMax(0.0,solution[iColumn]);
        sum += value;
        if (firstNonFixed<0)
          firstNonFixed=j;
        lastNonFixed=j;
        if (value>integerTolerance) {
          weight += weights_[j]*value;
          if (firstNonZero<0)
            firstNonZero=j;
          lastNonZero=j;
        }
      }
    }
    base += numberLinks_;
  }
  assert (lastNonZero-firstNonZero>=sosType_) ;
  // find where to branch
  assert (sum>0.0);
  weight /= sum;
  int iWhere;
  double separator=0.0;
  for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) 
    if (weight<weights_[iWhere+1])
      break;
  if (sosType_==1) {
    // SOS 1
    separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]);
  } else {
    // SOS 2
    if (iWhere==firstNonFixed)
      iWhere++;;
    if (iWhere==lastNonFixed-1)
      iWhere = lastNonFixed-2;
    separator = weights_[iWhere+1];
  }
  // create object
  CbcBranchingObject * branch;
  branch = new CbcLinkBranchingObject(model_,this,way,separator);
  return branch;
}