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

//#############################################################################
void 
MibSHeuristic::objCutHeuristic()
{

  /* Solve the LP relaxation with the new constraint d^2 y <= d^y* */

  MibSModel * model = MibSModel_;

  //OsiSolverInterface * oSolver = model->origLpSolver_;
  OsiSolverInterface * oSolver = model->getSolver();
  //OsiSolverInterface * hSolver = new OsiCbcSolverInterface();
  OsiSolverInterface * hSolver = new OsiSymSolverInterface();

  double objSense(model->getLowerObjSense());  
  int lCols(model->getLowerDim());
  int uCols(model->getUpperDim());
  int tCols(lCols + uCols);
  int * lColIndices = model->getLowerColInd();
  int * uColIndices = model->getUpperColInd();
  double * lObjCoeffs = model->getLowerObjCoeffs();

  hSolver->loadProblem(*oSolver->getMatrixByCol(),
		       oSolver->getColLower(), oSolver->getColUpper(),
		       oSolver->getObjCoefficients(),
		       oSolver->getRowLower(), oSolver->getRowUpper());

  int j(0);
  for(j = 0; j < tCols; j++){
    if(oSolver->isInteger(j))
      hSolver->setInteger(j);
  }

  double * optLowerSolutionOrd = model->bS_->optLowerSolutionOrd_;

  CoinPackedVector objCon;
  int i(0), index(0);
  double rhs(0.0);

  for(i = 0; i < lCols; i++){
    index = lColIndices[i];
    objCon.insert(index, lObjCoeffs[i] * objSense);
    //should this be ordered? and should lObjCoeffs by at index?
    //rhs += optLowerSolutionOrd_[i] * lObjCoeffs[i] * objSense;
    rhs += optLowerSolutionOrd[i] * lObjCoeffs[i] * objSense;
  }

  //Hmm, I think this was wrong before...?
  //  hSolver->addRow(objCon, - hSolver->getInfinity(), rhs);
  hSolver->addRow(objCon, rhs, hSolver->getInfinity());

  /* optimize w.r.t. to the UL objective with the new row */
  if(0){
    dynamic_cast<OsiCbcSolverInterface *> 
      (hSolver)->getModelPtr()->messageHandler()->setLogLevel(0);
  }
  else{
    dynamic_cast<OsiSymSolverInterface *> 
      (hSolver)->setSymParam("prep_level", -1);
    
    dynamic_cast<OsiSymSolverInterface *> 
      (hSolver)->setSymParam("verbosity", -2);

    dynamic_cast<OsiSymSolverInterface *> 
      (hSolver)->setSymParam("max_active_nodes", 1);
  }

  hSolver->branchAndBound();

  if(0)
    hSolver->writeLp("objcutheuristic");

  if(hSolver->isProvenOptimal()){

    MibSSolution *mibSol = NULL;

    OsiSolverInterface * lSolver = model->bS_->setUpModel(hSolver, true);

    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();

    const double * sol = hSolver->getColSolution();
    double objVal(lSolver->getObjValue() * objSense);
    double etol(etol_);
    double lowerObj = getLowerObj(sol, objSense);  
    
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      
      /*
	fix upper-level objective to current value and 
	reoptimize wrt to lower-level objective
      */
      
      //OsiSolverInterface * nSolver = new OsiCbcSolverInterface();
      OsiSolverInterface * nSolver = new OsiSymSolverInterface();
      nSolver->loadProblem(*oSolver->getMatrixByCol(),
			   oSolver->getColLower(), oSolver->getColUpper(),
			   oSolver->getObjCoefficients(),
			   oSolver->getRowLower(), oSolver->getRowUpper());
      for(j = 0; j < tCols; j++){
	if(oSolver->isInteger(j))
	  nSolver->setInteger(j);
      }
      

      CoinZeroN(nObjCoeffs, tCols);
      
      for(i = 0; i < lCols; i++){
	index = lColIndices[i];
	nObjCoeffs[index] = lObjCoeffs[i] * lObjSense;
      }
      
      nSolver->setObjective(nObjCoeffs);
      
      CoinPackedVector objCon;
      
      for(i = 0; i < tCols; i++){
	objCon.insert(i, uObjCoeffs[i] * uObjSense);
      }
      
      nSolver->addRow(objCon, upperObjVal, upperObjVal);
      nSolver->writeLp("beta1");
      if(0){
	dynamic_cast<OsiCbcSolverInterface *> 
	  (nSolver)->getModelPtr()->messageHandler()->setLogLevel(0);
      }
      else{
	 dynamic_cast<OsiSymSolverInterface *> 
	    (nSolver)->setSymParam("prep_level", -1);

	 dynamic_cast<OsiSymSolverInterface *> 
	    (nSolver)->setSymParam("verbosity", -2);

	 dynamic_cast<OsiSymSolverInterface *> 
	    (nSolver)->setSymParam("max_active_nodes", 1);
      }

      nSolver->branchAndBound();
     

      double * colsol = new double[tCols];

      if(nSolver->isProvenOptimal()){
	lowerObjVal = nSolver->getObjValue();
	CoinCopyN(nSolver->getColSolution(), tCols, colsol);
      }
      else{
	//just take the current solution
	lowerObjVal = sSolver->getObjValue();
	CoinCopyN(sSolver->getColSolution(), tCols, colsol);
      }

      delete[] nObjCoeffs;

//.........这里部分代码省略.........
    model.setDblParam(CbcModel::CbcMaximumSeconds,60.0*minutes);
  }
  // Switch off most output
  if (model.getNumCols()<3000) {
    model.messageHandler()->setLogLevel(1);
    //model.solver()->messageHandler()->setLogLevel(0);
  } else {
    model.messageHandler()->setLogLevel(2);
    model.solver()->messageHandler()->setLogLevel(1);
  }
  // Default strategy will leave cut generators as they exist already
  // so cutsOnlyAtRoot (1) ignored
  // numberStrong (2) is 5 (default)
  // numberBeforeTrust (3) is 5 (default is 0)
  // printLevel (4) defaults (0)
  CbcStrategyDefault strategy(true,5,5);
  // Set up pre-processing to find sos if wanted
  if (preProcess)
    strategy.setupPreProcessing(2);
  model.setStrategy(strategy);

  // Go round adding cuts to cutoff last solution
  // Stop after finding 20 best solutions
  for (int iPass=0;iPass<20;iPass++) {
    time1 = CoinCpuTime();
    // Do complete search
    model.branchAndBound();
    
    std::cout<<mpsFileName<<" took "<<CoinCpuTime()-time1<<" seconds, "
             <<model.getNodeCount()<<" nodes with objective "
             <<model.getObjValue()
             <<(!model.status() ? " Finished" : " Not finished")
             <<std::endl;
    // Stop if infeasible
    if (model.isProvenInfeasible())
      break;
    // Print solution if finished - we can't get names from Osi! - so get from OsiClp
    
    assert (model.getMinimizationObjValue()<1.0e50);
    OsiSolverInterface * solver = model.solver();
    int numberColumns = solver->getNumCols();
    
    const double * solution = model.bestSolution();
    //const double * lower = solver->getColLower();
    //const double * upper = solver->getColUpper();

    // Get names from solver1 (as OsiSolverInterface may lose)
    std::vector<std::string> columnNames = *solver1.getModelPtr()->columnNames();
    
    int iColumn;
    std::cout<<std::setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14);
    
    std::cout<<"--------------------------------------"<<std::endl;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      double value=solution[iColumn];
      if (fabs(value)>1.0e-7&&solver->isInteger(iColumn)) 
	std::cout<<std::setw(6)<<iColumn<<" "
                 <<columnNames[iColumn]<<" "
                 <<value
          //<<" "<<lower[iColumn]<<" "<<upper[iColumn]
                 <<std::endl;
    }
    std::cout<<"--------------------------------------"<<std::endl;
  
    std::cout<<std::resetiosflags(std::ios::fixed|std::ios::showpoint|std::ios::scientific);
    /* Now add cut to reference copy.
       resetting to reference copy also gets rid of best solution so we
       should either save best solution, reset, add cut OR
       add cut to reference copy then reset - this is doing latter
    */
    OsiSolverInterface * refSolver = model.referenceSolver();
    const double * bestSolution = model.bestSolution();
    const double * originalLower = refSolver->getColLower();
    const double * originalUpper = refSolver->getColUpper();
    CoinPackedVector cut;
    double rhs = 1.0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      double value=bestSolution[iColumn];
      if (solver->isInteger(iColumn)) {
        // only works for 0-1 variables
        assert (originalLower[iColumn]==0.0&&
                originalUpper[iColumn]==1.0);
        // double check integer
        assert (fabs(floor(value+0.5)-value)<1.0e-5);
        if (value>0.5) {
          // at 1.0
          cut.insert(iColumn,-1.0);
          rhs -= 1.0;
        } else {
          // at 0.0
          cut.insert(iColumn,1.0);
        }
      }
    }
    // now add cut
    refSolver->addRow(cut,rhs,COIN_DBL_MAX);
    model.resetToReferenceSolver();
  }
  return 0;
}    

int main(int argc, char **argv)
{
    VRPH_version();

    int i, j, k, n, status, num_attempts, *sol_buff, *IP_sol_buff;
    char in_file[200];
    double lambda, best_heur_sol=VRP_INFINITY;
    bool first_sol=false, bootstrap=false;;
    VRPSolution *fresh_solution;
    OsiSolverInterface *si;
    const double *x;
    int last_num_cols=0, route_id=0;
    time_t start, stop;
    int *orderings[MAX_ROUTES];
    for(i=0;i<MAX_ROUTES;i++)
        orderings[i]=NULL;
    

    // Set timing counters to 0
    heur_time=mip_time=0;

    // Check arguments
    if(argc<5)
    {
        fprintf(stderr,"Usage: %s -f input_file -n num_runs [-v,-b,-c max_columns -d cols_to_delete]\n",
            argv[0]);
        fprintf(stderr,"\t Will solve the problem num_solutions times and add the routes\n");
        fprintf(stderr,"\t to a set partitioning problem.\n");
        fprintf(stderr,"\t Other options:\n");
        fprintf(stderr,"\t -v runs in verbose mode\n");
        fprintf(stderr,"\t -b will use bootstrapping where we send the set partitioning\n"
                       "\t    solution back to the metaheuristic solver\n");
        fprintf(stderr,"\t -c max_columns will allow this many active columns/variables in the IP.\n");
        fprintf(stderr,"\t    Default value is max_columns=500\n");
        fprintf(stderr,"\t -d num_cols_to_delete will delete this many columns once we have too many\n");
        fprintf(stderr,"\t    in the IP. Default value is num_cols_to_delete=100\n");
        exit(-1);
    }

    // Set defaults
    verbose=false;
    max_columns=500;
    num_cols_to_delete=100;

    // Parse command line
    for(i=0;i<argc;i++)
    {
        if(strcmp(argv[i],"-f")==0)
            strcpy(in_file,argv[i+1]);
        if(strcmp(argv[i],"-n")==0)
            num_attempts=atoi(argv[i+1]);
        if(strcmp(argv[i],"-v")==0)
            verbose=true;
         if(strcmp(argv[i],"-b")==0)
            bootstrap=true;
         if(strcmp(argv[i],"-c")==0)
             max_columns=atoi(argv[i+1]);
         if(strcmp(argv[i],"-d")==0)
             num_cols_to_delete=atoi(argv[i+1]);
    }

    // This is the # of non-VRPH_DEPOT nodes
    n=VRPGetDimension(in_file);
    // This will be used to import/export solutions
    fresh_solution = new VRPSolution(n);
    // Create buffers for importing solutions
    sol_buff= new int[n+2];
    IP_sol_buff = new int[n+2];

    // Declare an OSI interface
    si=new OsiGlpkSolverInterface;
    si->setIntParam(OsiNameDiscipline,2);

    for(i=0;i<n;i++)
    {
        si->addRow(0,NULL,NULL,1,1);
    }

    // Declare a VRP of the right size and import the file
    VRP V(n);
    ClarkeWright CW(n);
    VRPRoute route(n);
    V.read_TSPLIB_file(in_file);
    // Set up a "route warehouse" to store the routes to be added to the IP
    V.route_wh=new VRPRouteWarehouse(HASH_TABLE_SIZE);
  
    // Set up a minimization problem
    si->setObjSense(1);
    // Set to error only output
    si->setHintParam(OsiDoReducePrint,true, OsiHintDo);
    // Unfortunately GLPK still prints out something regarding the conflict graph

    for(i=0;i<num_attempts;i++)
    {
        if(i==0 || !bootstrap) 
        {
            lambda=.5+1.5*lcgrand(0);
            // Start with a clean VRP object
            V.reset();
            CW.Construct(&V, lambda, false);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                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];
                delete [] object;
                // Can we move any rows to cuts
                const int * cutMarker = coinModel.cutMarker();
                if (cutMarker && 0) {
                    printf("AMPL CUTS OFF until global cuts fixed\n");
                    cutMarker = NULL;
                }
                if (cutMarker) {
                    // Row copy
                    const CoinPackedMatrix * matrixByRow = finalModel->getMatrixByRow();
                    const double * elementByRow = matrixByRow->getElements();
                    const int * column = matrixByRow->getIndices();
                    const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
                    const int * rowLength = matrixByRow->getVectorLengths();

OsiSolverInterface* cpropagation_preprocess(CPropagation *cp, int nindexes[])
{
    if(cp->varsToFix == 0)
    {
        /* printf("There are no variables to remove from the problem!\n"); */
        return NULL; /* returns a pointer to original solver */
    }

    const double *colLb = problem_vars_lower_bound(cp->problem), *colUb = problem_vars_upper_bound(cp->problem);
    const double *objCoef = problem_vars_obj_coefs(cp->problem);
    const char *ctype = problem_vars_type(cp->problem);

    double sumFixedObj = 0.0; /* stores the sum of objective coefficients of all variables fixed to 1 */

    OsiSolverInterface *preProcSolver = new OsiClpSolverInterface();
    preProcSolver->setIntParam(OsiNameDiscipline, 2);
    preProcSolver->messageHandler()->setLogLevel(0);
    preProcSolver->setHintParam(OsiDoReducePrint,true,OsiHintTry);
    //preProcSolver->setObjName(cp->solver->getObjName());

    for(int i = 0, j = 0; i < problem_num_cols(cp->problem); i++)
    {
        nindexes[i] = -1;
        if(cp->isToFix[i] == UNFIXED)
        {
            preProcSolver->addCol(0, NULL, NULL, colLb[i], colUb[i], objCoef[i]);
            preProcSolver->setColName(j, problem_var_name(cp->problem, i));
            if(problem_var_type(cp->problem, i) == CONTINUOUS)
                preProcSolver->setContinuous(j);
            else 
                preProcSolver->setInteger(j);
            nindexes[i] = j++;
        }
        else if(cp->isToFix[i] == ACTIVATE)
            sumFixedObj += objCoef[i];
    }

    if(fabs(sumFixedObj) > EPS)
    {
        /* adding a variable with cost equals to the sum of all coefficients of variables fixed to 1 */
        preProcSolver->addCol(0, NULL, NULL, 1.0, 1.0, sumFixedObj);
        preProcSolver->setColName(preProcSolver->getNumCols()-1, "sumFixedObj");
        preProcSolver->setInteger(preProcSolver->getNumCols()-1);
    }

    for(int idxRow = 0; idxRow < problem_num_rows(cp->problem); idxRow++)
    {
        const int nElements = problem_row_size(cp->problem, idxRow);
        const int *idxs = problem_row_idxs(cp->problem, idxRow);
        const double *coefs = problem_row_coefs(cp->problem, idxRow);
        vector< int > vidx; vidx.reserve(problem_num_cols(cp->problem));
        vector< double > vcoef; vcoef.reserve(problem_num_cols(cp->problem));
        double activeCoefs = 0.0;

        for(int i = 0; i < nElements; i++)
        {
            if(cp->isToFix[idxs[i]] == UNFIXED)
            {
                assert(nindexes[idxs[i]] >= 0 && nindexes[idxs[i]] < problem_num_cols(cp->problem));
                vidx.push_back(nindexes[idxs[i]]);
                vcoef.push_back(coefs[i]);
            }
            else if(cp->isToFix[idxs[i]] == ACTIVATE)
            	activeCoefs += coefs[i];
        }

        if(!vidx.empty())
        {
        	double rlb, rub;
        	const char sense = problem_row_sense(cp->problem, idxRow);
        	
        	if(sense == 'E')
            {
                rlb = problem_row_rhs(cp->problem, idxRow) - activeCoefs;
                rub = problem_row_rhs(cp->problem, idxRow) - activeCoefs;
            }
        	else if(sense == 'L')
            {
                rlb = preProcSolver->getInfinity();
                rub = problem_row_rhs(cp->problem, idxRow) - activeCoefs;
            }
        	else if(sense == 'G')
            {
                rlb = problem_row_rhs(cp->problem, idxRow) - activeCoefs;
                rub = preProcSolver->getInfinity();
            }
        	else
        	{
        		fprintf(stderr, "Error: invalid type of constraint!\n");
        		exit(EXIT_FAILURE);
        	}

        	preProcSolver->addRow((int)vcoef.size(), &vidx[0], &vcoef[0], rlb, rub);
            preProcSolver->setRowName(idxRow, problem_row_name(cp->problem, idxRow));
        }
	}

    return preProcSolver;
}