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

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

  /* 
     optimize wrt to lower-level objective 
     over current feasible lp feasible region 
  */

  MibSModel * model = MibSModel_;

  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 * lColIndices = model->getLowerColInd();
  int * uColIndices = model->getUpperColInd();
  double * lObjCoeffs = model->getLowerObjCoeffs();
  
  //int tCols(lCols + uCols);
  int tCols(oSolver->getNumCols());
  //assert(tCols == oSolver->getNumCols());

  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 * nObjCoeffs = new double[tCols];
  int i(0), index(0);
  
  CoinZeroN(nObjCoeffs, tCols);

  for(i = 0; i < lCols; i++){
    index = lColIndices[i];
    nObjCoeffs[index] = lObjCoeffs[i];
  }

  //MibS objective sense is the opposite of OSI's!
  hSolver->setObjSense(objSense);

  hSolver->setObjective(nObjCoeffs);
 
  //double cutoff(model->getCutoff());
  double cutoff(model->getKnowledgeBroker()->getIncumbentValue());

  if(model->getNumSolutions()){
  
    CoinPackedVector objCon;
    //double rhs(cutoff * objSense);
    //double smlTol(1.0);
    double rhs(cutoff);
       
    for(i = 0; i < tCols; i++){
      objCon.insert(i, oSolver->getObjCoefficients()[i] 
		    * oSolver->getObjSense());
    }
    
    hSolver->addRow(objCon, - hSolver->getInfinity(), rhs);
  }
  
  if(0)
     hSolver->writeLp("lobjheurstic");

  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(hSolver->isProvenOptimal()){

    double upperObjVal(0.0);

    /*****************NEW ******************/

    MibSSolution *mibSol = NULL;

    OsiSolverInterface * lSolver = model->bS_->setUpModel(hSolver, true);
//.........这里部分代码省略.........

//.........这里部分代码省略.........

  int randchoice(0); 
  if(0)
    std::cout << "randchoice " << randchoice << std::endl;


  double cost(0.0);

  //starting from the largest, fix corr upper-level variables
  //then, with these fixed, solve the lower-level problem
  //this yields a feasible solution

  iter = lObjCoeffsOrd.begin(); 
  
  while((usedBudget < intBudget) && (iter != lObjCoeffsOrd.end())){
    
    ind_min_wt = iter->second;
    cost = intCost[ind_min_wt];
    testsol[uColIndices[ind_min_wt]] = 1.0;
    double min_wt = iter->first;
    
    if(0){
      std::cout << "upper: " << ind_min_wt << " " 
		<< uColIndices[ind_min_wt] << " "  
		<< oSolver->getColUpper()[uColIndices[ind_min_wt]] << " " 
		<< oSolver->getColLower()[uColIndices[ind_min_wt]] << std::endl;
      
      std::cout << "lower: " << ind_min_wt << " " 
		<< lColIndices[ind_min_wt] << " "  
		<< oSolver->getColUpper()[lColIndices[ind_min_wt]] << std::endl;
    }

    //if((oSolver->getColUpper()[uColIndices[ind_min_wt]] == 1.0) 
       //&& (oSolver->getColUpper()[lColIndices[ind_min_wt]] > 0)){
    if(oSolver->getColUpper()[uColIndices[ind_min_wt]] > etol_){ 
      
      //if(((usedBudget + cost) <= intBudget) 
      // && checkLowerFeasibility(oSolver, testsol)){
      if((usedBudget + cost) <= intBudget){
	
	//FIXME: SHOULD BE CHECKING FOR CURRENT BOUNDS HERE  
	//fix the corresponding upper-level variable to 1
	randchoice = random() % 2;
	if(0)
	  std::cout << "randchoice " << random << std::endl;
	if(randchoice){
	  fixedVars[ind_min_wt] = 1.0;
	  usedBudget += intCost[ind_min_wt];
	}
      }
    }
    else{
      
      testsol[uColIndices[ind_min_wt]] = 0;
      //break;
      
    }
    iter++;
  }
  
  /*
    now we find a feasible solution by fixing upper-level vars
    and solving the lower-level problem
  */
  
  double * incumbentSol = new double[tCols];
  double * colsol = new double[tCols];

  CoinZeroN(colsol, tCols);

  for(i = 0; i < uCols; i++){
    colsol[uColIndices[i]] = fixedVars[i];
    if(fixedVars[i] == 1.0)
      if(0)
	std::cout << "fixed " << i << std::endl;
  }

  bfSol * sol = getBilevelSolution(colsol, lObjSense * oSolver->getInfinity());

  if(sol){
    double incumbentObjVal = sol->getObjVal();
    CoinCopyN(sol->getColumnSol(), tCols, incumbentSol);
    
    MibSSolution * mibSol = new MibSSolution(tCols,
					     incumbentSol,
					     incumbentObjVal,
					     model);
    
    model->storeSolution(BlisSolutionTypeHeuristic, mibSol);
  }

  //bestObjVal_ = incumbentObjVal;
  //CoinCopyN(incumbentSol, tCols, bestSol_);

  delete [] incumbentSol;
  delete [] testsol;
  //delete [] colsol;
  //delete [] fixedVars;
  //delete sol;
}

//.........这里部分代码省略.........

  int i, ij, k=0;
  int nPlus1=n+1;
  int offset = AtildeStarts[m]+m;
  for (i=0; i<m; i++){
    for (ij=AtildeStarts[i];ij<AtildeStarts[i]+AtildeLengths[i];ij++){
      BElements[k]=AtildeElements[ij];
      BElements[k+offset]=-AtildeElements[ij];
      BIndices[k]= AtildeIndices[ij];
      BIndices[k+offset]= AtildeIndices[ij];

      k++;
    }
    BElements[k]=btilde[i];
    BElements[k+offset]=btilde[i];
    BIndices[k]=n;
    BIndices[k+offset]=nPlus1;
    BStarts[i]= AtildeStarts[i]+i;
    BStarts[i+m]=offset+BStarts[i];// = AtildeStarts[m]+m+AtildeStarts[i]+i
    BLengths[i]= AtildeLengths[i]+1;
    BLengths[i+m]= AtildeLengths[i]+1;
    k++;
  }

  BStarts[twoM]=BStarts[twoM-1]+BLengths[twoM-1];

  // Cols that will be deleted each iteration
  int BNumColsLessOne=BNumCols-1;
  int BNumColsLessTwo=BNumCols-2;
  const int delCols[2] = {BNumColsLessOne, BNumColsLessTwo};

  // Set lower bound on u and v
  // u_0, v_0 will be reset as free
  const double solverINFINITY = si.getInfinity();
  double * BColLowers = new double[BNumCols];
  double * BColUppers = new double[BNumCols];
  CoinFillN(BColLowers,BNumCols,0.0);  
  CoinFillN(BColUppers,BNumCols,solverINFINITY); 

  // Set row lowers and uppers.
  // The rhs is zero, for but the last two rows.
  // For these the rhs is beta_
  double * BRowLowers = new double[BNumRows];
  double * BRowUppers = new double[BNumRows];
  CoinFillN(BRowLowers,BNumRows,0.0);  
  CoinFillN(BRowUppers,BNumRows,0.0);
  BRowLowers[BNumRows-2]=beta_;
  BRowUppers[BNumRows-2]=beta_;
  BRowLowers[BNumRows-1]=beta_;
  BRowUppers[BNumRows-1]=beta_;


  // Calculate base objective <<x^T,Atilde^T>,u>
  // Note: at each iteration coefficient u_0
  //       changes to <x^T,e_j>
  //       w=(u,v,beta,v_0,u_0) size 2m+3
  //       So, BOjective[2m+2]=x[j]
  double * BObjective= new double[BNumCols];
  double * Atildex = new double[m];
  CoinFillN(BObjective,BNumCols,0.0);
  Atilde->times(x,Atildex); // Atildex is size m, x is size n
  CoinDisjointCopyN(Atildex,m,BObjective); 

  // Number of cols and size of Elements vector
  // in B without the v_0 and u_0 cols
  int BFullSizeLessThree = BFullSize-3;

void 
HeuristicInnerApproximation::extractInnerApproximation(Bonmin::OsiTMINLPInterface & nlp, OsiSolverInterface &si,
  const double * x, bool getObj) {
   printf("************  Start extracting inner approx");
   int n;
   int m;
   int nnz_jac_g;
   int nnz_h_lag;
   Ipopt::TNLP::IndexStyleEnum index_style;
   Bonmin::TMINLP2TNLP * problem = nlp.problem(); 
   //Get problem information
   problem->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
   
   Bonmin::vector<int> jRow(nnz_jac_g);
   Bonmin::vector<int> jCol(nnz_jac_g);
   Bonmin::vector<double> jValues(nnz_jac_g);
   problem->eval_jac_g(n, NULL, 0, m, nnz_jac_g, jRow(), jCol(), NULL);
   if(index_style == Ipopt::TNLP::FORTRAN_STYLE)//put C-style
   {
     for(int i = 0 ; i < nnz_jac_g ; i++){
       jRow[i]--;
       jCol[i]--;
     }
   }
   
   //get Jacobian
   problem->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL,
       jValues());
   
   Bonmin::vector<double> g(m);
   problem->eval_g(n, x, 1, m, g());
   
   Bonmin::vector<int> nonLinear(m);
   //store non linear constraints (which are to be removed from IA)
   int numNonLinear = 0;
   const double * rowLower = nlp.getRowLower();
   const double * rowUpper = nlp.getRowUpper();
   const double * colLower = nlp.getColLower();
   const double * colUpper = nlp.getColUpper();
   assert(m == nlp.getNumRows());
   double infty = si.getInfinity();
   double nlp_infty = nlp.getInfinity();
   Bonmin::vector<Ipopt::TNLP::LinearityType>  constTypes(m);
   Bonmin::vector<Ipopt::TNLP::LinearityType>  varTypes(n);
   problem->get_constraints_linearity(m, constTypes());
   problem->get_variables_linearity(n, varTypes());
   for (int i = 0; i < m; i++) {
     if (constTypes[i] == Ipopt::TNLP::NON_LINEAR) {
       nonLinear[numNonLinear++] = i;
     }
   }
   Bonmin::vector<double> rowLow(m - numNonLinear);
   Bonmin::vector<double> rowUp(m - numNonLinear);
   int ind = 0;
   for (int i = 0; i < m; i++) {
     if (constTypes[i] != Ipopt::TNLP::NON_LINEAR) {
       if (rowLower[i] > -nlp_infty) {
         //   printf("Lower %g ", rowLower[i]);
         rowLow[ind] = (rowLower[i]);
       } else
         rowLow[ind] = -infty;
       if (rowUpper[i] < nlp_infty) {
         //   printf("Upper %g ", rowUpper[i]);
         rowUp[ind] = (rowUpper[i]);
       } else
         rowUp[ind] = infty;
       ind++;
     }
   
   }
   
   CoinPackedMatrix mat(true, jRow(), jCol(), jValues(), nnz_jac_g);
   mat.setDimensions(m, n); // In case matrix was empty, this should be enough
   
   //remove non-linear constraints
   mat.deleteRows(numNonLinear, nonLinear());
   
   int numcols = nlp.getNumCols();
   Bonmin::vector<double> obj(numcols);
   for (int i = 0; i < numcols; i++)
     obj[i] = 0.;
   
   si.loadProblem(mat, nlp.getColLower(), nlp.getColUpper(), 
                  obj(), rowLow(), rowUp());
   const Bonmin::TMINLP::VariableType* variableType = problem->var_types();
   for (int i = 0; i < n; i++) {
     if ((variableType[i] == Bonmin::TMINLP::BINARY) || (variableType[i] == Bonmin::TMINLP::INTEGER))
       si.setInteger(i);
   }
   if (getObj) {
     bool addObjVar = false;
     if (problem->hasLinearObjective()) {
       double zero;
       Bonmin::vector<double> x0(n, 0.);
       problem->eval_f(n, x0(), 1, zero);
       si.setDblParam(OsiObjOffset, -zero);
       //Copy the linear objective and don't create a dummy variable.
       problem->eval_grad_f(n, x, 1, obj());
       si.setObjective(obj());
     } else {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	if ( fabs( yCoef[i] ) > EPSILON_ ) {

	  if ( sign[i] == CGLFLOW_COL_CONTPOS )
	    cutCoef[cutLen] = coef[i] * yCoef[i];
	  else
	    cutCoef[cutLen] = -coef[i] * yCoef[i];

	  cutInd[cutLen++] = ind[i];
	}

	if ( fabs( xCoef[i] ) > EPSILON_ ) {
	  if ( VUB.getVar() != UNDEFINED_ ) {
	    cutCoef[cutLen] = xCoef[i];
	    cutInd[cutLen++] = VUB.getVar();
	  }
	  else
	    cutRHS -= xCoef[i];
	}
      }

      if ( ( sign[i] == CGLFLOW_COL_BINPOS ) ||
	   ( sign[i] == CGLFLOW_COL_BINNEG ) ) {
	if (fabs(yCoef[i]) > EPSILON_ || fabs(xCoef[i]) > EPSILON_) {
	  if (sign[i] == CGLFLOW_COL_BINPOS)
	    cutCoef[cutLen] = coef[i] * yCoef[i] + xCoef[i];
	  else
	    cutCoef[cutLen] = -coef[i] * yCoef[i] + xCoef[i];
	  cutInd[cutLen++] = ind[i];
	}
      }
    }

    for ( i = 0; i < cutLen; ++i ) {
      for ( j = 0; j < i; j++ ) {
	if ( cutInd[j] == cutInd[i] ) { /* Duplicate*/
	  cutCoef[j] += cutCoef[i];
	  cutInd[i] = -1;
	}
      }
    }

    for ( j = 0, i = 0; i < cutLen; ++i ) {
      if ( ( cutInd[i] == -1 ) || ( fabs( cutCoef[i]) < EPSILON_ ) ){
	/* Small coeff*/
      }
      else {
	cutCoef[j] = cutCoef[i];
	cutInd[j] = cutInd[i];
	j++;
      }
    }

    cutLen = j;
    // Skip if no elements ? - bug somewhere
    assert (cutLen);

    // Recheck the violation.
    violation = 0.0;
    for (i = 0; i < cutLen; ++i)
      violation += cutCoef[i] * xlp[cutInd[i]];

    violation -= cutRHS;

    if ( violation > TOLERANCE_ ) {
      flowCut.setRow(cutLen, cutInd, cutCoef);
      flowCut.setLb(-1.0 * si.getInfinity());
      flowCut.setUb(cutRHS);
      flowCut.setEffectiveness(violation);
      generated = true;

      if(CGLFLOW_DEBUG) {
	std::cout << "generateOneFlowCover(): Found a cut" << std::endl;
      }
    }
    else {
      if(CGLFLOW_DEBUG) {
	std::cout << "generateOneFlowCover(): Lost a cut" << std::endl;
      }
    }
  }

  //-------------------------------------------------------------------------
  delete [] sign;
  delete [] up;
  delete [] x;
  delete [] y;
  delete [] candidate;
  delete [] label;
  delete [] ratio;
  delete [] rho;
  delete [] xCoef;
  delete [] yCoef;
  delete [] mt;
  delete [] M;
  delete [] order;
  delete [] cutInd;
  delete [] cutCoef;

  return generated;
}

void 
HeuristicInnerApproximation::extractInnerApproximation(OsiTMINLPInterface & nlp, OsiSolverInterface &si,
                                                       const double * x, bool getObj) {
   int n;
   int m;
   int nnz_jac_g;
   int nnz_h_lag;
   Ipopt::TNLP::IndexStyleEnum index_style;
   TMINLP2TNLP * problem = nlp.problem(); 
   //Get problem information
   problem->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
   
   vector<int> jRow(nnz_jac_g);
   vector<int> jCol(nnz_jac_g);
   vector<double> jValues(nnz_jac_g);
   problem->eval_jac_g(n, NULL, 0, m, nnz_jac_g, jRow(), jCol(), NULL);
   if(index_style == Ipopt::TNLP::FORTRAN_STYLE)//put C-style
   {
     for(int i = 0 ; i < nnz_jac_g ; i++){
       jRow[i]--;
       jCol[i]--;
     }
   }
   
   //get Jacobian
   problem->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL,
       jValues());
   
   vector<double> g(m);
   problem->eval_g(n, x, 1, m, g());
   
   vector<int> nonLinear(m);
   //store non linear constraints (which are to be removed from IA)
   int numNonLinear = 0;
   const double * rowLower = nlp.getRowLower();
   const double * rowUpper = nlp.getRowUpper();
   const double * colLower = nlp.getColLower();
   const double * colUpper = nlp.getColUpper();
   assert(m == nlp.getNumRows());
   double infty = si.getInfinity();
   double nlp_infty = nlp.getInfinity();
   vector<Ipopt::TNLP::LinearityType>  constTypes(m);
   problem->get_constraints_linearity(m, constTypes());
   for (int i = 0; i < m; i++) {
     if (constTypes[i] == Ipopt::TNLP::NON_LINEAR) {
       nonLinear[numNonLinear++] = i;
     }
   }
   vector<double> rowLow(m - numNonLinear);
   vector<double> rowUp(m - numNonLinear);
   int ind = 0;
   for (int i = 0; i < m; i++) {
     if (constTypes[i] != Ipopt::TNLP::NON_LINEAR) {
       if (rowLower[i] > -nlp_infty) {
         //   printf("Lower %g ", rowLower[i]);
         rowLow[ind] = (rowLower[i]);
       } else
         rowLow[ind] = -infty;
       if (rowUpper[i] < nlp_infty) {
         //   printf("Upper %g ", rowUpper[i]);
         rowUp[ind] = (rowUpper[i]);
       } else
         rowUp[ind] = infty;
       ind++;
     }
   
   }
   
   CoinPackedMatrix mat(true, jRow(), jCol(), jValues(), nnz_jac_g);
   mat.setDimensions(m, n); // In case matrix was empty, this should be enough
   
   //remove non-linear constraints
   mat.deleteRows(numNonLinear, nonLinear());
   
   int numcols = nlp.getNumCols();
   vector<double> obj(numcols);
   for (int i = 0; i < numcols; i++)
     obj[i] = 0.;
   
   si.loadProblem(mat, nlp.getColLower(), nlp.getColUpper(), 
                  obj(), rowLow(), rowUp());
   const Bonmin::TMINLP::VariableType* variableType = problem->var_types();
   for (int i = 0; i < n; i++) {
     if ((variableType[i] == TMINLP::BINARY) || (variableType[i]
         == TMINLP::INTEGER))
       si.setInteger(i);
   }
   if (getObj) {
     bool addObjVar = false;
     if (problem->hasLinearObjective()) {
       double zero;
       vector<double> x0(n, 0.);
       problem->eval_f(n, x0(), 1, zero);
       si.setDblParam(OsiObjOffset, -zero);
       //Copy the linear objective and don't create a dummy variable.
       problem->eval_grad_f(n, x, 1, obj());
       si.setObjective(obj());
     } else {
       addObjVar = true;
     }
//.........这里部分代码省略.........

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

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

   // Build our own instance from scratch

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

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

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

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

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

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

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

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

  

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

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

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

   return 0;
}

//.........这里部分代码省略.........
	   (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);
     }
#endif
     delete [] integerVars;

  }else{
     nSolver = solver_;
  }

#define SYM_VERSION_IS_WS strcmp(SYMPHONY_VERSION, "WS")  

#if SYMPHONY_VERSION_IS_WS
  if (feasCheckSolver == "SYMPHONY" && probType == 1 && warmStartLL &&
      !newOsi && doDualFixing){ //Interdiction

     /** Get upper bound from best known (feasible) lower level solution and try 
	 to fix additional variables by sensitivity analysis **/

     std::vector<std::pair<AlpsKnowledge*, double> > solutionPool;
     model_->getKnowledgeBroker()->
	getAllKnowledges(AlpsKnowledgeTypeSolution, solutionPool);

     const double * sol; 
     double objval, Ub(objSense*nSolver->getInfinity());
     BlisSolution* blisSol;
     std::vector<std::pair<AlpsKnowledge*, double> >::const_iterator si;
     for (si = solutionPool.begin(); si != solutionPool.end(); ++si){
	blisSol = dynamic_cast<BlisSolution*>(si->first);
	sol = blisSol->getValues();
	for (i = 0; i < uCols; i++){
	   if (lpSol[uColIndices[i]] > 1 - etol &&
	       sol[lColIndices[i]] > 1-etol){
	      break;
	   }
	}
	if (i == uCols && -objSense*blisSol->getQuality() < Ub){
	   Ub = -objSense*blisSol->getQuality();
	}
     }

     /** Figure out which variables get fixed by the upper level solution **/

     int *newUbInd = new int[uCols];
     int *newLbInd = new int[uCols];
     double *newUbVal = new double[uCols];
     double *newLbVal = new double[uCols];
     double newLb;
     for (i = 0; i < uCols; i++){
	newUbInd[i] = uColIndices[i];
	newLbInd[i] = uColIndices[i];
	newLbVal[i] = 0;
	if (lpSol[uColIndices[i]] > 1 - etol){
	   newUbVal[i] = 0;
	}else{
	   newUbVal[i] = 1;
	}

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;
}