AAP_lp.cpp

// $Id: AAP_lp.cpp,v 1.3 2003/12/03 01:52:30 magh Exp $

/* ------------------------------------------------------------------------
 Author: Matthew Galati (magh@lehigh.edu)

 (c) Copyright 2003 Lehigh University. All Rights Reserved.

 This software is licensed under the Common Public License. Please see 
 accompanying file for terms.    
---------------------------------------------------------------------------*/

#ifdef COIN_USE_CPX
#include "OsiCpxSolverInterface.hpp"
#include "cplex.h"
typedef OsiCpxSolverInterface OsiXxxSolverInterface;
#endif
#ifdef COIN_USE_CLP
#include "OsiClpSolverInterface.hpp"
typedef OsiClpSolverInterface OsiXxxSolverInterface;
#endif
#ifdef COIN_USE_OSL
#include "OsiOslSolverInterface.hpp"
typedef OsiOslSolverInterface OsiXxxSolverInterface;
#endif

#include "BCP_lp.hpp"
#include "BCP_lp_node.hpp"
#include "AAP_lp.hpp"
#include "AAP_var.hpp"
#include "AAP_user_data.hpp"

/*--------------------------------------------------------------------------*/
void AAP_lp::unpack_module_data(BCP_buffer& buf) {
  //Unpack the information passed from the TM process
  lp_par.unpack(buf); //parameter list for LP
  aap = new AAP(buf); //the AAP instance
}

/*--------------------------------------------------------------------------*/
void AAP_lp::pack_var_algo(const BCP_var_algo* var, BCP_buffer& buf){
  //Pack an algorithmic variable (AAP_var is derived from BCP_var_algo)
  {
    const AAP_var * mvar = dynamic_cast<const AAP_var*>(var);
    if(mvar){
      mvar->pack(buf);
      return;
    }
  }
  throw BCP_fatal_error("AAP_lp::pack_var_algo() : unknown cut type!\n");
}

/*--------------------------------------------------------------------------*/
BCP_var_algo* AAP_lp::unpack_var_algo(BCP_buffer& buf){
  //Unpack an algorithmic variable (AAP_var is derived from BCP_var_algo)
  return new AAP_var(buf);
}

/*--------------------------------------------------------------------------*/
void AAP_lp::pack_user_data(const BCP_user_data* ud, BCP_buffer& buf){
  //Pack user data
  {
    const AAP_user_data * data = dynamic_cast<const AAP_user_data*>(ud);
    if(data){
      data->pack(buf);
      return;
    }
  }
  throw BCP_fatal_error("AAP_tm::pack_user_data() : cast failed!\n");
}

/*--------------------------------------------------------------------------*/
BCP_user_data* AAP_lp::unpack_user_data(BCP_buffer& buf){
  return new AAP_user_data(buf);
}

/*--------------------------------------------------------------------------*/
OsiSolverInterface* AAP_lp::initialize_solver_interface(){
  //Initialize the OSI Xxx solver interface
  OsiXxxSolverInterface * si = new OsiXxxSolverInterface();

#ifdef COIN_USE_CPX
  //Using cplex, we need to shut off the presolver to make sure get a dual 
  //ray in the case of primal infeasibility.
  CPXsetintparam(si->getEnvironmentPtr(),CPX_PARAM_PREIND, 0);
#endif

  return si;
}

/*--------------------------------------------------------------------------*/
void AAP_lp::display_lp_solution(const BCP_lp_result& lpres,
                                 const BCP_vec<BCP_var*>& vars,
                                 const BCP_vec<BCP_cut*>& cuts,
                                 const bool final_lp_solution){

#ifdef DBG_PRINT 
  //Display the current LP solution
  cout << "UB: " << upper_bound() << " LB: " 
       << getLpProblemPointer()->node->true_lower_bound << endl;

  const int varnum = vars.size();
  for (int c = 0; c < varnum; ++c) {   
    if(lpres.x()[c] > tolerance){
      const AAP_var* mvar = dynamic_cast<const AAP_var*>(vars[c]);
      if(mvar){
        cout << "Column " << c << " : " << lpres.x()[c] << endl; 
        mvar->print(aap->getDimension());
      }
    }
  }
#endif

}

/*--------------------------------------------------------------------------*/
double AAP_lp::compute_lower_bound(const double old_lower_bound,
                                   const BCP_lp_result& lpres,
                                   const BCP_vec<BCP_var*>& vars,
                                   const BCP_vec<BCP_cut*>& cuts){

  // In order to get the bound we are solving the subproblem over the 
  // set of reduced costs. So, if we find any columns with negative 
  // reduced cost, we should save them for generate_vars_in_lp.  

  // Only try to generate_vars here if we have good duals 
  if(lpres.termcode() & BCP_ProvenPrimalInf)
    return -DBL_MAX;

  bool new_var = generate_vars(lpres.pi(), false);

  if(new_var){
    return old_lower_bound;
    //return max(old_lower_bound, lpres.objval() + most_reduced_cost);
  }
  else{
    const int tc = lpres.termcode();
    if (tc & BCP_ProvenOptimal)
      return lpres.objval();
    else return old_lower_bound;
  }
}

/*--------------------------------------------------------------------------*/
void AAP_lp::generate_vars_in_lp(const BCP_lp_result& lpres,
                                 const BCP_vec<BCP_var*>& vars,
                                 const BCP_vec<BCP_cut*>& cuts,
                                 const bool before_fathom,
                                 BCP_vec<BCP_var*>& new_vars,
                                 BCP_vec<BCP_col*>& new_cols){

  /*
    Generate variables with negative reduced cost to be added 
    to the LP by solving an assignment problem. This work has
    already been done in compute_lower_bound and was saved in 
    generated_vars.

    The reduced cost of an assignment point is:
    rc_s = sum{ijk in H} s_ijk c_ijk                                    [c_s]
    - sum{i in I} sum{jk in J x K} s_ijk uhat_i       [dual of I constraints]
    - 1 * alpha                                [dual of convexity constraint]
  */
    
  new_vars.reserve(new_vars.size() + generated_vars.size());
  for (unsigned int i = 0; i < generated_vars.size(); ++i)
    new_vars.push_back(generated_vars[i]);
  generated_vars.clear();
}

/*--------------------------------------------------------------------------*/
void AAP_lp::restore_feasibility(const BCP_lp_result& lpres,
                                 const std::vector<double*> dual_rays,
                                 const BCP_vec<BCP_var*>& vars,
                                 const BCP_vec<BCP_cut*>& cuts,
                                 BCP_vec<BCP_var*>& vars_to_add,
                                 BCP_vec<BCP_col*>& cols_to_add){

  // Generate a variable that breaks the certificate of infeasibility 
  // using the dual ray as the cost vector and solving an assignment problem. 

  for(unsigned int r = 0; r < dual_rays.size(); r++){
    generate_vars(dual_rays[r], true);
    vars_to_add.reserve(vars_to_add.size() + generated_vars.size());
    for (unsigned int i = 0; i < generated_vars.size(); ++i)
      vars_to_add.push_back(generated_vars[i]);
    
    generated_vars.clear();
  }
}

/*--------------------------------------------------------------------------*/
void AAP_lp::construct_cost(const double * pi, const bool from_restore,
                            double * ap_costv, map<int,int> & jk_mini){
  
  /*
    Assignment Problem objective:
    min {jk in J x K} c'_jk x_jk    
    
    c_'jk = min{i in I} c_ijk - uhat_i OR
    c_'jk = min{i in I} ray_i
  */
  
#ifdef DBG_PRINT
  const AAP_user_data * datatmp 
    = dynamic_cast<const AAP_user_data*>(get_user_data());
  cout << "ONES: ";
  for(unsigned int i = 0; i < datatmp->ones.size(); i++)
    cout << datatmp->ones[i] << " ";
  cout << endl;
  cout << "ZEROS: ";
  for(unsigned int i = 0; i < datatmp->zeros.size(); i++)
    cout << datatmp->zeros[i] << " ";
  cout << endl;
#endif

  const int dimension = aap->getDimension();
  
  int index, index0, i, j, k, mult;

  //Construct the full cost vector across all ijk
  double * rc = new double[aap->n_cols()];
  if(from_restore){
    CoinFillN(rc, aap->n_cols(), 0.0);
    mult = 1;
  }
  else{
    CoinDisjointCopyN(aap->getAssignCost(), aap->n_cols(), rc);
    mult = -1;
  }

  int c = 0;
  for(i = 0; i < dimension; i++){
    for(j = 0; j < dimension; j++){
      for(k = 0; k < dimension; k++){
        rc[c] += (mult * pi[i]);
        c++;
      }
    }
  }

  //Enforce branching rules using a bigM cost
  const AAP_user_data * data 
    = dynamic_cast<const AAP_user_data*>(get_user_data());

  //For x_ijk = 1, for all s in F' such that s_ijk != 1, set lambda_s = 0
  for(index = 0; index < (int)data->ones.size(); index++)
    rc[data->ones[index]] += -bigM;

  //For x_ijk = 0, for all s in F' such that s_ijk != 1, set lambda_s = 0
  for(index = 0; index < (int)data->zeros.size(); index++)
    rc[data->zeros[index]] += bigM;


  //Now construct the cost vector for the AP problem
  //(j,k) -> (i*,j,k) for i* arg min{i in I} rc_ijk
  c = 0;
  for(j = 0; j < dimension; j++){
    for(k = 0; k < dimension; k++){
      index0 = index3(0,j,k,dimension);
      pair<int,double> mini = make_pair(index0,rc[index0]);
      for(i = 1; i < dimension; i++){
        index = index3(i,j,k,dimension);
        if(rc[index] < mini.second)
          mini = make_pair(index, rc[index]);
      }
      jk_mini.insert(make_pair(c, mini.first));
      ap_costv[c] = mini.second;
      c++;
    }
  } 

  delete [] rc;
}

/*--------------------------------------------------------------------------*/
bool AAP_lp::generate_vars(const double * pi, const bool from_restore){

  const int dimension = aap->getDimension();

  double objective_value;
  double * solution  = new double[dimension * dimension];
  double * ap_costv  = new double[dimension * dimension];//BCP will delete
  map<int,int> jk_mini;

  construct_cost(pi, from_restore, ap_costv, jk_mini);
  
  int return_code = solve_assignment_problem(ap_costv, solution, 
                                             objective_value); 
  
  if(return_code){
    delete [] solution;
    return false;
  }

  double most_reduced_cost;
  const AAP_user_data * data 
    = dynamic_cast<const AAP_user_data*>(get_user_data());
  const int n_ones = data->ones.size();

  if(from_restore)
    most_reduced_cost = objective_value + pi[dimension] + (n_ones * bigM);
  else
    most_reduced_cost = objective_value - pi[dimension] + (n_ones * bigM);
  
  if(most_reduced_cost >= -tolerance){
    delete [] solution;
    return false;
  }

  //If the column has negative reduced cost, then transform the AP point 
  //back into AAP indices and create an AAP_var.
  double ap_cost = 0.0;
  vector<int> ap; ap.reserve(dimension);
  int index_aap, j, k, c = 0;
  for(j = 0; j < dimension; j++){
    for(k = 0; k < dimension; k++){
      if(solution[c] > tolerance){
        index_aap = jk_mini[c];
        ap.push_back(index_aap);
        ap_cost += aap->getAssignCost()[index_aap];
      }
      c++;
    }
  }
  AAP_var * new_var = new AAP_var(ap, ap_cost);
  generated_vars.push_back(new_var);
  delete [] solution;
  return true;
}

/*--------------------------------------------------------------------------*/
int AAP_lp::solve_assignment_problem(double * ap_costv,
                                     double * solution,
                                     double & objective_value){
  /*
    Construct the Assignment Problem as an OSI instance

    Assignment Problem - relaxation of I:
    min {jk in J x K} c'_jk x_jk
    sum {k in K}            x_jk = 1, forall j in J
    sum {j in J}            x_jk = 1, forall k in K
    x_jk in {0,1}, for all jk in J x K

    c_'jk = min{i in I} c_ijk - uhat_i OR
    c_'jk = min{i in I} ray_i
  */

  OsiXxxSolverInterface * siAP = new OsiXxxSolverInterface();
  
  const int dimension = aap->getDimension();
  const int n_cols    = dimension * dimension;
  const int n_rows    = 2 * dimension;
  
  int i, j, k, index;
  double * col_lb    = new double[n_cols];
  double * col_ub    = new double[n_cols];
  double * row_lb_ub = new double[n_rows];
  CoinFillN(col_lb, n_cols, 0.0);
  CoinFillN(col_ub, n_cols, 1.0);
  CoinFillN(row_lb_ub, n_rows, 1.0);  


  //Adjust for branching:
  const AAP_user_data * data 
    = dynamic_cast<const AAP_user_data*>(get_user_data());

  //This should have been already taken care of using bigM, but using 
  //bounds should be used to help the LP solver.

  //For x_ijk = 1, for all s in F' such that s_ijk != 1, set lambda_s = 0
  for(index = 0; index < (int)data->ones.size(); index++){
    ijk(data->ones[index], dimension, i, j, k);
    col_lb[index2(j,k,dimension)] = 1.0;
  }


  CoinPackedMatrix * M =  new CoinPackedMatrix(false,0,0);
  M->setDimensions(0, n_cols);
  
  //sum {k in K}            x_jk = 1, forall j in J
  for(j = 0; j < dimension; j++){
    CoinPackedVector row;
    for(k = 0; k < dimension; k++){
      index = index2(j,k,dimension);
      row.insert(index, 1.0);
    }
    M->appendRow(row);
  }

  //sum {j in J}            x_jk = 1, forall k in K
  for(k = 0; k < dimension; k++){
    CoinPackedVector row;
    for(j = 0; j < dimension; j++){
      index = index2(j,k,dimension);
      row.insert(index, 1.0);
    }
    M->appendRow(row);
  }

  //Note on OSI usage:
  //assignProblem assumes ownership of the arguments (and deletes memory)
  //loadProblem makes a copy of the arguments
  siAP->assignProblem(M, col_lb, col_ub, ap_costv, row_lb_ub, row_lb_ub);

  //AP is Totally Unimodular, so LP is sufficient for integral solutions
  siAP->initialSolve();

  int return_code = 0;
  if(siAP->isProvenOptimal()){
    objective_value = siAP->getObjValue();
    CoinDisjointCopyN(siAP->getColSolution(), n_cols, solution);
  }
  else
    return_code = 1;

  delete siAP;

  return return_code;
}

/*--------------------------------------------------------------------------*/
void AAP_lp::vars_to_cols(const BCP_vec<BCP_cut*>& cuts,
                          BCP_vec<BCP_var*>& vars,       
                          BCP_vec<BCP_col*>& cols,
                          const BCP_lp_result& lpres,
                          BCP_object_origin origin, bool allow_multiple){

  /*
    Expand the variable (AP point) into a column in the LP:
    min {s in F'}                c_s lambda_s             
    sum {s in F', jk in J x K} s_ijk lambda_s = 1, forall i in I (5)
    sum {s in F'}                    lambda_s = 1                (6)
  */
  const int varnum = vars.size();
  if (varnum == 0)
    return; 

  cols.reserve(varnum);
  for (int c = 0; c < varnum; ++c) {   
    const AAP_var* mvar = dynamic_cast<const AAP_var*>(vars[c]);
    if(mvar){

      CoinPackedVector col;

      //coefficient = number of i's in AP point (can be > 1)
      vector<int> coeff(aap->getDimension(), 0);
      for(unsigned int p = 0; p < mvar->ap.size(); p++)
        coeff[ijk_i(mvar->ap[p], aap->getDimension())]++;
  
      for(int i = 0; i < aap->getDimension(); i++)
        if(coeff[i] > 0) 
          col.insert(i, static_cast<double>(coeff[i])); //const           (5)
      col.insert(aap->getDimension(), 1.0);             //convexity const (6)

      cols.unchecked_push_back(new BCP_col(col, mvar->cost, 
                                           vars[c]->lb(), vars[c]->ub()));
    }
  }
}

---