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

  TNLPSolver::ReturnStatus LpBranchingSolver::
  solveFromHotStart(OsiTMINLPInterface* tminlp_interface)
  {
    TNLPSolver::ReturnStatus retstatus = TNLPSolver::solvedOptimal;

    // updated the bounds of the linear solver
    std::vector<int> diff_low_bnd_index;
    std::vector<double> diff_low_bnd_value;
    std::vector<int> diff_up_bnd_index;
    std::vector<double> diff_up_bnd_value;

    // Get the bounds.  We assume that the bounds in the linear solver
    // are always the original ones
    const int numCols = tminlp_interface->getNumCols();
    const double* colLow_orig = lin_->getColLower();
    const double* colUp_orig = lin_->getColUpper();
    const double* colLow = tminlp_interface->getColLower();
    const double* colUp = tminlp_interface->getColUpper();

    OsiSolverInterface * lin = lin_;
    // eventualy clone lin_
    if(warm_start_mode_ == Clone){
      lin = lin_->clone();
//      std::cout<<"Cloning it"<<std::endl;
    }
    // Set the bounds on the LP solver according to the changes in
    // tminlp_interface
    for (int i=0; i<numCols; i++) {
      const double& lo = colLow[i];
      if (colLow_orig[i] < lo) {
        if(warm_start_mode_ == Basis){
          diff_low_bnd_value.push_back(colLow_orig[i]);
          diff_low_bnd_index.push_back(i);
        }
        lin->setColLower(i,lo);
      }
      const double& up = colUp[i];
      if (colUp_orig[i] > up) {
        if(warm_start_mode_ == Basis){
          diff_up_bnd_index.push_back(i);
          diff_up_bnd_value.push_back(colUp_orig[i]);
        }
        lin->setColUpper(i,lo);
      }
    }

    if(warm_start_mode_ == Basis){
      lin->setWarmStart(warm_);
    }

    lin->resolve();

    double obj = lin->getObjValue();
    bool go_on = true;
    if (lin->isProvenPrimalInfeasible() || 
        lin->isDualObjectiveLimitReached()) {
      retstatus = TNLPSolver::provenInfeasible;
      go_on = false;
    }
    else if (lin->isIterationLimitReached()) {
      retstatus = TNLPSolver::iterationLimit;
      go_on = false;
    }
    else {
      if (maxCuttingPlaneIterations_ > 0 && go_on) {
        double violation;
        obj = ecp_->doEcpRounds(*lin, true, &violation);
        if (obj == COIN_DBL_MAX) {
          retstatus = TNLPSolver::provenInfeasible;
        }
        else if (violation <= 1e-8) {
          retstatus = TNLPSolver::solvedOptimal;
        }
      }
    }
    tminlp_interface->problem()->set_obj_value(obj);
    tminlp_interface->problem()->Set_x_sol(numCols, lin_->getColSolution());

    //restore the original bounds
    if(warm_start_mode_ == Basis){
      for (unsigned int i = 0; i < diff_low_bnd_index.size(); i++) {
        lin_->setColLower(diff_low_bnd_index[i],diff_low_bnd_value[i]);
      }
      for (unsigned int i = 0; i < diff_up_bnd_index.size(); i++) {
        lin_->setColUpper(diff_up_bnd_index[i],diff_up_bnd_value[i]);
      }
    }
    else {
      delete lin;
    }
    return retstatus;
  }

//.........这里部分代码省略.........
                    for (int j = 0; j < numCols; j++)
                        objectiveValue += solutions.solution[k][j] * originalObjCoeff[j];
                    cutoff = objectiveValue;
                    clpSolver->addRow(numCols, addRowIndex, originalObjCoeff, -COIN_DBL_MAX, cutoff);

                    // Todo: pick up the best solution in the block (not
                    // the last).
                    solutionValue = objectiveValue;
                    for (int m = 0; m < numCols; m++)
                        betterSolution[m] = solutions.solution[k][m];
                    numFeasibles++;
                }
            }

        // Go through the block of solution and decide if to call smallBB
        for (int k = startIndex; k < solutions.numberSolutions; k++) {
            if (solutions.numberUnsatisfied[k] <= fixThreshold) {
                // get new copy
                OsiSolverInterface * newSolver;
                newSolver = new OsiClpSolverInterface(*clpSolver);
                newSolver->setObjSense(1);
                newSolver->setObjective(originalObjCoeff);
                int numberColumns = newSolver->getNumCols();
                int numFixed = 0;

                // Fix the first fixThreshold number of integer vars
                // that are satisfied
                for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
                    if (newSolver->isInteger(iColumn)) {
                        double value = solutions.solution[k][iColumn];
                        double intValue = floor(value + 0.5);
                        if (fabs(value - intValue) < 1.0e-5) {
                            newSolver->setColLower(iColumn, intValue);
                            newSolver->setColUpper(iColumn, intValue);
                            numFixed++;
                            if (numFixed > numInt - fixThreshold)
                                break;
                        }
                    }
                }
                COIN_DETAIL_PRINT(printf("numFixed: %d\n", numFixed));
                COIN_DETAIL_PRINT(printf("fixThreshold: %f\n", fixThreshold));
		COIN_DETAIL_PRINT(printf("numInt: %d\n", numInt));
                double *newSolution = new double[numCols];
                double newSolutionValue;

                // Call smallBB on the modified problem
                int returnCode = smallBranchAndBound(newSolver, maxNode, newSolution,
                                                     newSolutionValue, cutoff, "mini");

                // If smallBB found a solution, update the better
                // solution and solutionValue (we gave smallBB our
                // cutoff, so it only finds improving solutions)
                if (returnCode == 1 || returnCode == 3) {
                    numFeasibles ++;
                    solutionValue = newSolutionValue;
                    for (int m = 0; m < numCols; m++)
                        betterSolution[m] = newSolution[m];
                    COIN_DETAIL_PRINT(printf("cutoff: %f\n", newSolutionValue));
                    COIN_DETAIL_PRINT(printf("time: %.2lf\n", CoinCpuTime() - start));
                }
                didMiniBB = 1;
                COIN_DETAIL_PRINT(printf("returnCode: %d\n", returnCode));

                //Update sumReturnCode array
                for (int iRC = 0; iRC < 6; iRC++) {

//.........这里部分代码省略.........
                // Both branch can be fathomed
                bStatus = -2;
            }
            else if (!downKeep) {
                // Down branch can be fathomed.
                lbInd[numLowerTightens] = colInd;
                newLB[numLowerTightens++] = ceil(lpX);
                //break;
            }
            else if (!upKeep) {
                // Up branch can be fathomed.
                ubInd[numUpperTightens] = colInd;
                newUB[numUpperTightens++] = floor(lpX);
                // break;
            }

            // Update pseudocost.
            if(downGood) {
                firstObjects[i]->pseudocost().update(-1, downDeg, lpX);
            }
            if(downGood) {
                firstObjects[i]->pseudocost().update(1, upDeg, lpX);
            }
        }

        //--------------------------------------------------
        // Set new bounds in lp solver for resolving
        //--------------------------------------------------

        if (bStatus != -2) {
            if (numUpperTightens > 0) {
                bStatus = -1;
                for (i = 0; i < numUpperTightens; ++i) {
                    solver->setColUpper(ubInd[i], newUB[i]);
                }
            }
            if (numLowerTightens > 0) {
                bStatus = -1;
                for (i = 0; i < numLowerTightens; ++i) {
                    solver->setColLower(lbInd[i], newLB[i]);
                }
            }
        }

        //--------------------------------------------------
        // Unmark hotstart and recover LP solver.
        //--------------------------------------------------

        solver->unmarkHotStart();
        solver->setColSolution(saveSolution);
        solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit);
        solver->setWarmStart(ws);
        delete ws;
    }

    //std::cout << "REL: bStatus = " << bStatus << std::endl;

    if (bStatus < 0) {
        // Infeasible or monotone.
        goto TERM_CREATE;
    }
    else {
        // All object's pseudocost have been initialized.
        // Sort them, and do strong branch for the unreliable one
        // NOTE: it set model->savedLpSolution.
        // model->feasibleSolution(numIntegerInfs, numObjectInfs);

//.........这里部分代码省略.........
            colInd = firstObjects[i]->columnIndex();
            
            lpX = saveSolution[colInd];
            
            BlisStrongBranch(model, objValue, colInd, lpX,
                             saveLower, saveUpper,
                             downKeep, downGood, downDeg,
                             upKeep, upGood, upDeg);
            
            if(!downKeep && !upKeep) {
                // Both branch can be fathomed
                bStatus = -2;
            }
            else if (!downKeep) {
                // Down branch can be fathomed.
                lbInd[numLowerTightens] = colInd;
                newLB[numLowerTightens++] = ceil(lpX);
            }
            else if (!upKeep) {
                // Up branch can be fathomed.
                ubInd[numUpperTightens] = colInd;
                newUB[numUpperTightens++] = floor(lpX);
            }
        }

        //--------------------------------------------------
        // Set new bounds in lp solver for resolving
        //--------------------------------------------------
        
        if (bStatus != -2) {
            if (numUpperTightens > 0) {
                bStatus = -1;
                for (i = 0; i < numUpperTightens; ++i) {
                    solver->setColUpper(ubInd[i], newUB[i]);
                }
            }
            if (numLowerTightens > 0) {
                bStatus = -1;
                for (i = 0; i < numLowerTightens; ++i) {
                    solver->setColLower(lbInd[i], newLB[i]);
                }
            }
        }
	
        //--------------------------------------------------
        // Unmark hotstart and recover LP solver.
        //--------------------------------------------------
        
        solver->unmarkHotStart();
        solver->setColSolution(saveSolution);
        solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit);
        solver->setWarmStart(ws);
        delete ws;
    }
    
    if (bStatus < 0) {
	goto TERM_CREATE;
    }
    else {
        // Create a set of candidate branching objects. 
        numBranchObjects_ = numInfs;
        branchObjects_ = new BcpsBranchObject* [numInfs];        
        
        // NOTE: it set model->savedLpSolution.
        
        sumDeg = 0.0;

//.........这里部分代码省略.........
                    } else {
                        if (fabs(value - floor(value + 0.5)) > 1.0e-12)
                            allIntegerCoeff = false;
                    }
                } else {
                    sumFixed += lower[jColumn] * value;
                }
            }
            double low = rowLower[iRow];
            if (low > -1.0e20) {
                low -= sumFixed;
                if (fabs(low - floor(low + 0.5)) > 1.0e-12)
                    allIntegerCoeff = false;
            }
            double up = rowUpper[iRow];
            if (up < 1.0e20) {
                up -= sumFixed;
                if (fabs(up - floor(up + 0.5)) > 1.0e-12)
                    allIntegerCoeff = false;
            }
            if (!allIntegerCoeff)
                continue; // can't do
            if (numberContinuous == 1) {
                // see if really integer
                // This does not allow for complicated cases
                if (low == up) {
                    if (fabs(value1) > 1.0e-3) {
                        value1 = 1.0 / value1;
                        if (fabs(value1 - floor(value1 + 0.5)) < 1.0e-12) {
                            // integer
                            changed[numberChanged++] = jColumn1;
                            solver->setInteger(jColumn1);
                            if (upper[jColumn1] > 1.0e20)
                                solver->setColUpper(jColumn1, 1.0e20);
                            if (lower[jColumn1] < -1.0e20)
                                solver->setColLower(jColumn1, -1.0e20);
                        }
                    }
                } else {
                    if (fabs(value1) > 1.0e-3) {
                        value1 = 1.0 / value1;
                        if (fabs(value1 - floor(value1 + 0.5)) < 1.0e-12) {
                            // This constraint will not stop it being integer
                            ignore[iRow] = 1;
                        }
                    }
                }
            } else if (numberContinuous == 2) {
                if (low == up) {
                    /* need general theory - for now just look at 2 cases -
                       1 - +- 1 one in column and just costs i.e. matching objective
                       2 - +- 1 two in column but feeds into G/L row which will try and minimize
                    */
                    if (fabs(value1) == 1.0 && value1*value2 == -1.0 && !lower[jColumn1]
                            && !lower[jColumn2]) {
                        int n = 0;
                        int i;
                        double objChange = direction * (objective[jColumn1] + objective[jColumn2]);
                        double bound = CoinMin(upper[jColumn1], upper[jColumn2]);
                        bound = CoinMin(bound, 1.0e20);
                        for ( i = columnStart[jColumn1]; i < columnStart[jColumn1] + columnLength[jColumn1]; i++) {
                            int jRow = row[i];
                            double value = element[i];
                            if (jRow != iRow) {
                                which[n++] = jRow;
                                changeRhs[jRow] = value;

bool
BlisConGenerator::generateCons(OsiCuts & coinCuts , bool fullScan)
{
    bool status = false;

    if (strategy_ == -2) {
        // This con generator has been disabled.
        return false;
    }

    OsiSolverInterface * solver = model_->solver();

#if defined(BLIS_DEBUG_MORE)
    std::cout << "model_->getNodeCount() = " << model_->getNodeCount()
              << std::endl;
#endif

    if ( fullScan ||
         ((strategy_ > 0) && (model_->getNumNodes() % strategy_) == 0) ) {

        //--------------------------------------------------
        // Start to generate cons ...
        //--------------------------------------------------

	int j;
        double start = CoinCpuTime();
        int numConsBefore = coinCuts.sizeCuts();
        int numRowsBefore = coinCuts.sizeRowCuts();

        assert(generator_ != NULL);

        CglProbing* generator = dynamic_cast<CglProbing *>(generator_);

        if (!generator) {
            generator_->generateCuts(*solver, coinCuts);
        }
        else {
            // It is probing - return tight column bound
 	    CglTreeInfo info;
            generator->generateCutsAndModify(*solver, coinCuts, &info);
            const double * tightLower = generator->tightLower();
            const double * lower = solver->getColLower();
            const double * tightUpper = generator->tightUpper();
            const double * upper = solver->getColUpper();
            const double * solution = solver->getColSolution();

            int numberColumns = solver->getNumCols();
            double primalTolerance = 1.0e-8;
            for (j = 0; j < numberColumns; ++j) {
                if ( (tightUpper[j] == tightLower[j]) &&
                     (upper[j] > lower[j]) ) {
                    // fix column j
                    solver->setColLower(j, tightLower[j]);
                    solver->setColUpper(j, tightUpper[j]);
                    if ( (tightLower[j] > solution[j] + primalTolerance) ||
                         (tightUpper[j] < solution[j] - primalTolerance) ) {
                        status = true;
                    }
                }
            }
        } // EOF probing.

        //--------------------------------------------------
        // Remove zero length row cuts.
        //--------------------------------------------------

	int numRowCons = coinCuts.sizeRowCuts();
	for (j = numRowsBefore; j < numRowCons; ++j) {
	    OsiRowCut & rCut = coinCuts.rowCut(j);
	    int len = rCut.row().getNumElements();
#ifdef BLIS_DEBUG_MORE
	    std::cout << "Cut " << j<<": length = " << len << std::endl;
#endif
	    if (len == 0) {
		// Empty cuts
		coinCuts.eraseRowCut(j);
		--j;
		--numRowCons;
#ifdef BLIS_DEBUG
		std::cout << "WARNING: Empty cut from " << name_ << std::endl;
#endif
	    }
	    else if (len < 0) {
#ifdef BLIS_DEBUG
		std::cout << "ERROR: Cut length = " << len << std::endl;
#endif
		// Error
		assert(0);
	    }
	}

        //--------------------------------------------------
        // Update statistics.
        //--------------------------------------------------

        ++calls_;
        numConsGenerated_ += (coinCuts.sizeCuts() - numConsBefore);
        time_ += (CoinCpuTime() - start);
        if (numConsGenerated_ == 0) {
            ++noConsCalls_;
//.........这里部分代码省略.........

// Generate constraints for the model data contained in si.
// The generated constraints are inserted into and returned in the 
// constraint pool.
// Default implementation use Cgl cut generators.
bool
BlisConGenerator::generateConstraints(BcpsConstraintPool &conPool)
{
    bool status = false;
    
    OsiSolverInterface * solver = model_->solver();
    
#if defined(BLIS_DEBUG_MORE)
    std::cout << "model_->getNodeCount() = " << model_->getNodeCount()
              << std::endl;
#endif
    
    //--------------------------------------------------
    // Start to generate constraints...
    //--------------------------------------------------

    assert(generator_ != NULL);
    
    int j;
    OsiCuts newOsiCuts;
    CglProbing* generator = dynamic_cast<CglProbing *>(generator_);
    
    if (generator) {
	// It is CglProbing - return tight column bounds
        CglTreeInfo info;
	generator->generateCutsAndModify(*solver, newOsiCuts, &info);
	const double * tightLower = generator->tightLower();
	const double * lower = solver->getColLower();
	const double * tightUpper = generator->tightUpper();
	const double * upper = solver->getColUpper();
	const double * solution = solver->getColSolution();
	
	int numberColumns = solver->getNumCols();
	double primalTolerance = 1.0e-8;
	for (j = 0; j < numberColumns; ++j) {
	    if ( (tightUpper[j] == tightLower[j]) &&
		 (upper[j] > lower[j]) ) {
		// Fix column j
		solver->setColLower(j, tightLower[j]);
		solver->setColUpper(j, tightUpper[j]);
		if ( (tightLower[j] > solution[j] + primalTolerance) ||
		     (tightUpper[j] < solution[j] - primalTolerance) ) {
		    status = true;
		}
	    }
	}
    }
    else {
	// Other Cgl cut generators
	generator_->generateCuts(*solver, newOsiCuts);
    }
    
    //--------------------------------------------------
    // Create blis constraints and remove zero length row cuts.
    //--------------------------------------------------
    
    int numNewConstraints = newOsiCuts.sizeRowCuts();
    for (j = 0; j < numNewConstraints; ++j) {
	OsiRowCut & rCut = newOsiCuts.rowCut(j);
	int len = rCut.row().getNumElements();

#ifdef BLIS_DEBUG_MORE
	std::cout << "Cut " << j<<": length = " << len << std::endl;
#endif
	if (len > 0) {
	    // Create BlisConstraints from OsiCuts.
	    BlisConstraint *blisCon = BlisOsiCutToConstraint(&rCut);
	    conPool.addConstraint(blisCon);
	}
	else if (len == 0) {
	    // Empty cuts
#ifdef BLIS_DEBUG
	    std::cout << "WARNING: Empty cut from " << name_ << std::endl;
#endif
	}
	else {
#ifdef BLIS_DEBUG
	    std::cout << "ERROR: Cut length = " << len << std::endl;
#endif
	    // Error
	    assert(0);
	}
    }

    // Adjust cut strategy.
    if ( (strategy_ == BlisCutStrategyAuto) &&
	 (noConsCalls_ > BLIS_CUT_DISABLE) ) {
	strategy_ = BlisCutStrategyNone;
    }

    return status;
}

//.........这里部分代码省略.........
                nelLargest = CoinMax(nelLargest, nelCreate);
            }
            if (nTotal > maxTotal)
                badModel = true;
            if (!badModel) {
                // Now arrays for building
                nelLargest = CoinMax(nelLargest, nLargest) + 1;
                double * buildObj = new double [nLargest];
                double * buildElement = new double [nelLargest];
                int * buildStart = new int[nLargest+1];
                int * buildRow = new int[nelLargest];
                // alow for integers in knapsacks
                OsiObject ** object = new OsiObject * [numberKnapsack+nTotal];
                int nSOS = 0;
                int nObj = numberKnapsack;
                for (iKnapsack = 0; iKnapsack < numberKnapsack; iKnapsack++) {
                    knapsackStart[iKnapsack] = finalModel->getNumCols();
                    iRow = knapsackRow[iKnapsack];
                    int nCreate = 10000;
                    coinModel.expandKnapsack(iRow, nCreate, buildObj, buildStart, buildRow, buildElement);
                    // Redo row numbers
                    for (iColumn = 0; iColumn < nCreate; iColumn++) {
                        for (int j = buildStart[iColumn]; j < buildStart[iColumn+1]; j++) {
                            int jRow = buildRow[j];
                            jRow = lookupRow[jRow];
                            assert (jRow >= 0 && jRow < nRow);
                            buildRow[j] = jRow;
                        }
                    }
                    finalModel->addCols(nCreate, buildStart, buildRow, buildElement, NULL, NULL, buildObj);
                    int numberFinal = finalModel->getNumCols();
                    for (iColumn = numberOther; iColumn < numberFinal; iColumn++) {
                        if (markKnapsack[iKnapsack] < 0) {
                            finalModel->setColUpper(iColumn, maxCoefficient);
                            finalModel->setInteger(iColumn);
                        } else {
                            finalModel->setColUpper(iColumn, maxCoefficient + 1.0);
                            finalModel->setInteger(iColumn);
                        }
                        OsiSimpleInteger * sosObject = new OsiSimpleInteger(finalModel, iColumn);
                        sosObject->setPriority(1000000);
                        object[nObj++] = sosObject;
                        buildRow[iColumn-numberOther] = iColumn;
                        buildElement[iColumn-numberOther] = 1.0;
                    }
                    if (markKnapsack[iKnapsack] < 0) {
                        // convexity row
                        finalModel->addRow(numberFinal - numberOther, buildRow, buildElement, 1.0, 1.0);
                    } else {
                        int iColumn = markKnapsack[iKnapsack];
                        int n = numberFinal - numberOther;
                        buildRow[n] = iColumn;
                        buildElement[n++] = -fabs(coefficient[iKnapsack]);
                        // convexity row (sort of)
                        finalModel->addRow(n, buildRow, buildElement, 0.0, 0.0);
                        OsiSOS * sosObject = new OsiSOS(finalModel, n - 1, buildRow, NULL, 1);
                        sosObject->setPriority(iKnapsack + SOSPriority);
                        // Say not integral even if is (switch off heuristics)
                        sosObject->setIntegerValued(false);
                        object[nSOS++] = sosObject;
                    }
                    numberOther = numberFinal;
                }
                finalModel->addObjects(nObj, object);
                for (iKnapsack = 0; iKnapsack < nObj; iKnapsack++)
                    delete object[iKnapsack];

//.........这里部分代码省略.........
      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);
	  }
	  array[iRow]=a;
	}
      }
      for (j=0;j<n;j++) {
	int iRow = which[j];
	// moving to point will increase row solution by this
	double distance = array[iRow];
	if (distance>1.0e-8) {
	  if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e-8) {
	    possible=false;
	    break;
	  }
	} else if (distance<-1.0e-8) {
	  if (distance+rowSolution[iRow]<rowLower[iRow]-1.0e-8) {
	    possible=false;
	    break;
	  } 
	}
      }
      for (j=0;j<n;j++)
	array[which[j]]=0.0;
      delete [] array;
      delete [] which;
      if (possible) {
	printf("possible feas region %d %d %d\n",firstNonZero,lastNonZero,iWhere);
	firstNonZero=iWhere;
	lastNonZero=iWhere;
      }
    }
  }
#else
  assert (lastNonZero-firstNonZero<sosType_) ;
#endif
  base=0;
  for (j=0;j<firstNonZero;j++) {
    for (int k=0;k<numberLinks_;k++) {
      int iColumn = which_[base+k];
      solver->setColUpper(iColumn,0.0);
    }
    base += numberLinks_;
  }
  // skip
  base += numberLinks_;
  for (j=lastNonZero+1;j<numberMembers_;j++) {
    for (int k=0;k<numberLinks_;k++) {
      int iColumn = which_[base+k];
      solver->setColUpper(iColumn,0.0);
    }
    base += numberLinks_;
  }
}

BlisReturnStatus
BlisStrongBranch(BlisModel *model, double objValue, int colInd, double x,
		 const double *saveLower, const double *saveUpper,
		 bool &downKeep, bool &downFinished, double &downDeg,
		 bool &upKeep, bool &upFinished, double &upDeg)
{
    BlisReturnStatus status = BlisReturnStatusOk;
    int lpStatus = 0;

    int j, numIntInfDown, numObjInfDown;

    double newObjValue;
    
    OsiSolverInterface * solver = model->solver();
    
    int numCols = solver->getNumCols();
    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();

    // Restore bounds
    int numDiff = 0;

    BlisSolution* ksol = NULL;

    int ind = model->getIntObjIndices()[colInd];
    BlisObjectInt *intObj = dynamic_cast<BlisObjectInt *>(model->objects(ind));
    
#ifdef BLIS_DEBUG_MORE
    for (j = 0; j < numCols; ++j) {
	if (saveLower[j] != lower[j]) {
	    //solver->setColLower(j, saveLower[j]);
            ++numDiff;
	}
	if (saveUpper[j] != upper[j]) {
	    //solver->setColUpper(j, saveUpper[j]);
            ++numDiff;
	}
    }
    std::cout << "BEFORE: numDiff = " << numDiff << std::endl;
#endif	 
   
    //------------------------------------------------------
    // Branching down.
    //------------------------------------------------------

    solver->setColUpper(colInd, floor(x));
    solver->solveFromHotStart();
    
    newObjValue = solver->getObjSense() * solver->getObjValue();
    downDeg = newObjValue - objValue;
    
    if (solver->isProvenOptimal()) {
	lpStatus = 0; // optimal
#ifdef BLIS_DEBUG_MORE
        printf("STRONG: COL[%d]: downDeg=%g, x=%g\n", colInd, downDeg, x);
#endif
        // Update pseudocost
        intObj->pseudocost().update(-1, downDeg, x);
        model->setSharedObjectMark(ind);        

        // Check if ip feasible
        ksol = model->feasibleSolution(numIntInfDown, numObjInfDown);
        if (ksol) {
#ifdef BLIS_DEBUG_MORE
            printf("STRONG:Down:found a feasible solution\n");
#endif
            
            model->storeSolution(BlisSolutionTypeStrong, ksol);
	    downKeep = false;
        }
	else {
	    downKeep = true;
	}
	downFinished = true;
    }
    else if (solver->isIterationLimitReached() && 
	     !solver->isDualObjectiveLimitReached()) {
	lpStatus = 2;      // unknown 
	downKeep = true;
	downFinished = false;
    }
    else {
        downDeg = 1.0e20;
	lpStatus = 1; // infeasible
	downKeep = false;
	downFinished = false;
    }       
            
#ifdef BLIS_DEBUG_MORE
    std::cout << "Down: lpStatus = " << lpStatus << std::endl;
#endif
    
    // restore bounds
    numDiff = 0;
    for (j = 0; j < numCols; ++j) {
	if (saveLower[j] != lower[j]) {
	    solver->setColLower(j, saveLower[j]);
            ++numDiff;
	}
	if (saveUpper[j] != upper[j]) {
//.........这里部分代码省略.........

int
CbcHeuristicNaive::solution(double & solutionValue,
                            double * betterSolution)
{
    numCouldRun_++;
    // See if to do
    bool atRoot = model_->getNodeCount() == 0;
    int passNumber = model_->getCurrentPassNumber();
    if (!when() || (when() == 1 && model_->phase() != 1) || !atRoot || passNumber != 1)
        return 0; // switched off
    // Don't do if it was this heuristic which found solution!
    if (this == model_->lastHeuristic())
        return 0;
    numRuns_++;
    double cutoff;
    model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
    double direction = model_->solver()->getObjSense();
    cutoff *= direction;
    cutoff = CoinMin(cutoff, solutionValue);
    OsiSolverInterface * solver = model_->continuousSolver();
    if (!solver)
        solver = model_->solver();
    const double * colLower = solver->getColLower();
    const double * colUpper = solver->getColUpper();
    const double * objective = solver->getObjCoefficients();

    int numberColumns = model_->getNumCols();
    int numberIntegers = model_->numberIntegers();
    const int * integerVariable = model_->integerVariable();

    int i;
    bool solutionFound = false;
    CoinWarmStartBasis saveBasis;
    CoinWarmStartBasis * basis =
        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
    if (basis) {
        saveBasis = * basis;
        delete basis;
    }
    // First just fix all integers as close to zero as possible
    OsiSolverInterface * newSolver = cloneBut(7); // wassolver->clone();
    for (i = 0; i < numberIntegers; i++) {
        int iColumn = integerVariable[i];
        double lower = colLower[iColumn];
        double upper = colUpper[iColumn];
        double value;
        if (lower > 0.0)
            value = lower;
        else if (upper < 0.0)
            value = upper;
        else
            value = 0.0;
        newSolver->setColLower(iColumn, value);
        newSolver->setColUpper(iColumn, value);
    }
    newSolver->initialSolve();
    if (newSolver->isProvenOptimal()) {
        double solValue = newSolver->getObjValue() * direction ;
        if (solValue < cutoff) {
            // we have a solution
            solutionFound = true;
            solutionValue = solValue;
            memcpy(betterSolution, newSolver->getColSolution(),
                   numberColumns*sizeof(double));
            COIN_DETAIL_PRINT(printf("Naive fixing close to zero gave solution of %g\n", solutionValue));
            cutoff = solValue - model_->getCutoffIncrement();
        }
    }
    // Now fix all integers as close to zero if zero or large cost
    int nFix = 0;
    for (i = 0; i < numberIntegers; i++) {
        int iColumn = integerVariable[i];
        double lower = colLower[iColumn];
        double upper = colUpper[iColumn];
        double value;
        if (fabs(objective[i]) > 0.0 && fabs(objective[i]) < large_) {
            nFix++;
            if (lower > 0.0)
                value = lower;
            else if (upper < 0.0)
                value = upper;
            else
                value = 0.0;
            newSolver->setColLower(iColumn, value);
            newSolver->setColUpper(iColumn, value);
        } else {
            // set back to original
            newSolver->setColLower(iColumn, lower);
            newSolver->setColUpper(iColumn, upper);
        }
    }
    const double * solution = solver->getColSolution();
    if (nFix) {
        newSolver->setWarmStart(&saveBasis);
        newSolver->setColSolution(solution);
        newSolver->initialSolve();
        if (newSolver->isProvenOptimal()) {
            double solValue = newSolver->getObjValue() * direction ;
            if (solValue < cutoff) {
                // try branch and bound
//.........这里部分代码省略.........

// This version fixes stuff and does IP
int
CbcHeuristicLocal::solutionFix(double & objectiveValue,
                               double * newSolution,
                               const int * /*keep*/)
{
/*
  If when is set to off (0), or set to root (1) and we're not at the root,
  return. If this heuristic discovered the current solution, don't continue.
*/

    numCouldRun_++;
    // See if to do
    if (!when() || (when() == 1 && model_->phase() != 1))
        return 0; // switched off
    // Don't do if it was this heuristic which found solution!
    if (this == model_->lastHeuristic())
        return 0;
/*
  Load up a new solver with the solution.

  Why continuousSolver(), as opposed to solver()?
*/
    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
    const double * colLower = newSolver->getColLower();
    //const double * colUpper = newSolver->getColUpper();

    int numberIntegers = model_->numberIntegers();
    const int * integerVariable = model_->integerVariable();
/*
  The net effect here is that anything that hasn't moved from its lower bound
  will be fixed at lower bound.

  See comments in solution() w.r.t. asymmetric treatment of upper and lower
  bounds.
*/

    int i;
    int nFix = 0;
    for (i = 0; i < numberIntegers; i++) {
        int iColumn = integerVariable[i];
        const OsiObject * object = model_->object(i);
        // get original bounds
        double originalLower;
        double originalUpper;
        getIntegerInformation( object, originalLower, originalUpper);
        newSolver->setColLower(iColumn, CoinMax(colLower[iColumn], originalLower));
        if (!used_[iColumn]) {
            newSolver->setColUpper(iColumn, colLower[iColumn]);
            nFix++;
        }
    }
/*
  Try a `small' branch-and-bound search. The notion here is that we've fixed a
  lot of variables and reduced the amount of `free' problem to a point where a
  small BaB search will suffice to fully explore the remaining problem. This
  routine will execute integer presolve, then call branchAndBound to do the
  actual search.
*/
    int returnCode = 0;
#ifdef CLP_INVESTIGATE2
    printf("Fixing %d out of %d (%d continuous)\n",
           nFix, numberIntegers, newSolver->getNumCols() - numberIntegers);
#endif
    if (nFix*10 <= numberIntegers) {
        // see if we can fix more
        int * which = new int [2*(numberIntegers-nFix)];
        int * sort = which + (numberIntegers - nFix);
        int n = 0;
        for (i = 0; i < numberIntegers; i++) {
            int iColumn = integerVariable[i];
            if (used_[iColumn]) {
                which[n] = iColumn;
                sort[n++] = used_[iColumn];
            }
        }
        CoinSort_2(sort, sort + n, which);
        // only half fixed in total
        n = CoinMin(n, numberIntegers / 2 - nFix);
        int allow = CoinMax(numberSolutions_ - 2, sort[0]);
        int nFix2 = 0;
        for (i = 0; i < n; i++) {
            int iColumn = integerVariable[i];
            if (used_[iColumn] <= allow) {
                newSolver->setColUpper(iColumn, colLower[iColumn]);
                nFix2++;
            } else {
                break;
            }
        }
        delete [] which;
        nFix += nFix2;
#ifdef CLP_INVESTIGATE2
        printf("Number fixed increased from %d to %d\n",
               nFix - nFix2, nFix);
#endif
    }
    if (nFix*10 > numberIntegers) {
        returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, objectiveValue,
                                         objectiveValue, "CbcHeuristicLocal");
//.........这里部分代码省略.........

//#############################################################################
bfSol*
MibSHeuristic::getBilevelSolution1(const double * sol)
{

  /* 
     Find a bilevel feasible solution by solving the LL problem
     for a fixed UL solution, given by the UL portion of sol
  */


  MibSModel * model = MibSModel_;
  OsiSolverInterface * oSolver = model->getSolver();
  OsiSolverInterface * lSolver = new OsiCbcSolverInterface(oSolver);
  
  //double uObjSense(model->getSolver()->getObjSense());
  double lObjSense(model->getLowerObjSense());  
  int lCols(model->getLowerDim());
  int uCols(model->getUpperDim());
  int * lColIndices = model->getLowerColInd();
  int uRowNum = model->getUpperRowNum();
  int lRowNum = model->getLowerRowNum();
  int * uRowIndices = model->getUpperRowInd();
  int * lRowIndices = model->getLowerRowInd();
  int * uColIndices = model->getUpperColInd();
  double * lObjCoeffs = model->getLowerObjCoeffs();

  int tCols(uCols + lCols); 
  int i(0), index(0);

  /* delete the UL rows */
  lSolver->deleteRows(uRowNum, uRowIndices);

  /* Fix the UL variables to their current value in sol */

  for(i = 0; i < uCols; i++){
    index = uColIndices[i];
    lSolver->setColLower(index, sol[index]);
    lSolver->setColUpper(index, sol[index]);
  }

  /* Set the objective to the LL objective coefficients */

  double * nObjCoeffs = new double[tCols];
  CoinZeroN(nObjCoeffs, tCols);
      
  for(i = 0; i < lCols; i++){
    index = lColIndices[i];
    nObjCoeffs[index] = lObjCoeffs[i] * lObjSense;
  }
  
  lSolver->setObjective(nObjCoeffs);

  if(0){
    dynamic_cast<OsiCbcSolverInterface *> 
      (lSolver)->getModelPtr()->messageHandler()->setLogLevel(0);
  }
  else{
     dynamic_cast<OsiSymSolverInterface *> 
	(lSolver)->setSymParam("prep_level", -1);
     
     dynamic_cast<OsiSymSolverInterface *> 
	(lSolver)->setSymParam("verbosity", -2);
     
     dynamic_cast<OsiSymSolverInterface *> 
	(lSolver)->setSymParam("max_active_nodes", 1);
  }

  lSolver->branchAndBound();

  if(lSolver->isProvenOptimal()){

    double objVal(0.0);

    for(i = 0; i < tCols; i++)
      objVal += lSolver->getColSolution()[i] * oSolver->getObjCoefficients()[i];
    
    double * colsol = new double[tCols];
    CoinCopyN(lSolver->getColSolution(), tCols, colsol);
 
    bfSol * bfsol = 
      new bfSol(objVal, colsol);
    delete lSolver;
    return bfsol;
  }
  else{
    delete lSolver;
    return NULL;
  }

}

//.........这里部分代码省略.........
	
     /** If we found an upper bound, then do the dual fixing **/

     if (Ub < objSense*nSolver->getInfinity()){
	sym_environment *env = 
	   dynamic_cast<OsiSymSolverInterface *>(nSolver)->getSymphonyEnvironment();

	for (i = 0; i < uCols; i++){
	   if (newUbVal[i] == 1){
	      // Try fixing it to zero
	      newUbVal[i] = 0; 
	      sym_get_lb_for_new_rhs(env, 0, NULL, NULL, uCols, newLbInd,
				     newLbVal, uCols, newUbInd, newUbVal,
				     &newLb);
	      if (objSense*newLb > Ub + etol){
		 //Victory! This variable can be fixed to 1 permanently
		 newLbVal[i] = 1;
	      }
	      //Set upper bound back to 1
	      newUbVal[i] = 1;
	      
	      if (newLbVal[i] == 0){
		 // Try fixing it to one
		 newLbVal[i] = 1;
		 sym_get_lb_for_new_rhs(env, 0, NULL, NULL, uCols, newLbInd,
					newLbVal, uCols, newUbInd, newUbVal,
					&newLb);
		 if (objSense*newLb > Ub + etol){
		    //Victory! This variable can be fixed to 0 permanently
		    newUbVal[i] = 0;
		 }
		 newLbVal[i] = 0;
	      }
	   }
	   
	}
     }
     
     /** Now set the row bounds to account for fixings **/
     
     /** This is probably very frgaile. Assuming interdiction 
	 rows come last. It would be better to set variable 
	 bounds directly, but this doesn't seem to work right now. **/
     int iRowStart = lRows-uCols; 

#if 1
     for(i = iRowStart; i < lRows; i++){
	nSolver->setRowLower(i, newLbVal[i-iRowStart]);
	nSolver->setRowUpper(i, newUbVal[i-iRowStart]);
     }
#else
     for(i = 0; i < uCols; i++){
	nSolver->setColLower(i, newLbVal[i]);
	nSolver->setColUpper(i, newUbVal[i]);
     }
#endif     
     delete[] newUbInd;
     delete[] newLbInd;
     delete[] newUbVal;
     delete[] newLbVal;
  }else{
#endif
     //FIXME: NEED TO GET ROW SENSE HERE
     
     /** Get contribution of upper-level columns **/
     
     double * upComp = new double[lRows];
     CoinFillN(upComp, lRows, 0.0);
     
     for(i = 0; i < lRows; i++){
	index1 = lRowIndices[i];
	for(j = 0; j < uCols; j++){
	   index2 = uColIndices[j];
	   coeff = matrix->getCoefficient(index1, index2);
	   if (coeff != 0){
	      upComp[i] += coeff * lpSol[index2];
	   }
	}
     }
     
     /** Correct the row bounds to account for fixed upper-level vars **/
     
     for(i = 0; i < lRows; i++){
	nSolver->setRowLower(i, rowLb[i] - upComp[i]);
	nSolver->setRowUpper(i, rowUb[i] - upComp[i]);
     }
     
     delete [] upComp;

#if SYMPHONY_VERSION_IS_WS
  }
#endif

  //I don't think this is needed
  //if(!getWarmStart())
  //  setWarmStart(nSolver->getWarmStart());

  return nSolver;

}

//.........这里部分代码省略.........
        if (numberNodes < lastNode_ + 12)
            return 0;
        // Do at 50 and 100
        if ((numberNodes > 40 && numberNodes <= 50) || (numberNodes > 90 && numberNodes < 100))
            numberNodes = howOften_;
    }
    if ((numberNodes % howOften_) == 0 && (model_->getCurrentPassNumber() <= 1 ||
                                           model_->getCurrentPassNumber() == 999999)) {
        lastNode_ = model_->getNodeCount();
        OsiSolverInterface * solver = model_->solver();

        int numberIntegers = model_->numberIntegers();
        const int * integerVariable = model_->integerVariable();

        const double * currentSolution = solver->getColSolution();
        OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
        //const double * colLower = newSolver->getColLower();
        //const double * colUpper = newSolver->getColUpper();

        double primalTolerance;
        solver->getDblParam(OsiPrimalTolerance, primalTolerance);

        // Sort on distance
        double * distance = new double [numberIntegers];
        int * which = new int [numberIntegers];

        int i;
        int nFix = 0;
        double tolerance = 10.0 * primalTolerance;
        for (i = 0; i < numberIntegers; i++) {
            int iColumn = integerVariable[i];
            const OsiObject * object = model_->object(i);
            // get original bounds
            double originalLower;
            double originalUpper;
            getIntegerInformation( object, originalLower, originalUpper);
            double valueInt = bestSolution[iColumn];
            if (valueInt < originalLower) {
                valueInt = originalLower;
            } else if (valueInt > originalUpper) {
                valueInt = originalUpper;
            }
            baseSolution_[iColumn] = currentSolution[iColumn];
            distance[i] = fabs(currentSolution[iColumn] - valueInt);
            which[i] = i;
            if (fabs(currentSolution[iColumn] - valueInt) < tolerance)
                nFix++;
        }
        CoinSort_2(distance, distance + numberIntegers, which);
        nDifferent_ = numberIntegers - nFix;
        stepSize_ = nDifferent_ / 10;
        k_ = stepSize_;
        //nFix = numberIntegers-stepSize_;
        for (i = 0; i < nFix; i++) {
            int j = which[i];
            int iColumn = integerVariable[j];
            const OsiObject * object = model_->object(i);
            // get original bounds
            double originalLower;
            double originalUpper;
            getIntegerInformation( object, originalLower, originalUpper);
            double valueInt = bestSolution[iColumn];
            if (valueInt < originalLower) {
                valueInt = originalLower;
            } else if (valueInt > originalUpper) {
                valueInt = originalUpper;
            }
            double nearest = floor(valueInt + 0.5);
            newSolver->setColLower(iColumn, nearest);
            newSolver->setColUpper(iColumn, nearest);
        }
        delete [] distance;
        delete [] which;
        if (nFix > numberIntegers / 5) {
            //printf("%d integers have samish value\n",nFix);
            returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
                                             model_->getCutoff(), "CbcHeuristicVND");
            if (returnCode < 0)
                returnCode = 0; // returned on size
            else
                numRuns_++;
            if ((returnCode&1) != 0)
                numberSuccesses_++;
            //printf("return code %d",returnCode);
            if ((returnCode&2) != 0) {
                // could add cut
                returnCode &= ~2;
                //printf("could add cut with %d elements (if all 0-1)\n",nFix);
            } else {
                //printf("\n");
            }
            numberTries_++;
            if ((numberTries_ % 10) == 0 && numberSuccesses_*3 < numberTries_)
                howOften_ += static_cast<int> (howOften_ * decayFactor_);
        }

        delete newSolver;
    }
    return returnCode;
}