BCP_lp_user.cpp

Go to the documentation of this file.

// Copyright (C) 2000, International Business Machines
// Corporation and others.  All Rights Reserved.
#include <cmath>
#include <cstdarg>
#include <cstdio>

#include "CoinHelperFunctions.hpp"
#include "CoinTime.hpp"

#include "BCP_math.hpp"
#include "BCP_error.hpp"
#include "BCP_USER.hpp"
#include "BCP_var.hpp"
#include "BCP_cut.hpp"
#include "BCP_problem_core.hpp"
#include "BCP_lp_user.hpp"
#include "BCP_lp.hpp"
#include "BCP_lp_pool.hpp"
#include "BCP_lp_result.hpp"
#include "BCP_lp_node.hpp"
#include "BCP_lp_branch.hpp"
#include "BCP_problem_core.hpp"
#include "BCP_solution.hpp"
#include "BCP_functions.hpp"
#include "BCP_lp_functions.hpp"

//#############################################################################
// Informational methods for the user
double BCP_lp_user::upper_bound() const    { return p->ub(); }
bool BCP_lp_user::over_ub(double lb) const { return p->over_ub(lb); }
int BCP_lp_user::current_phase() const     { return p->phase; }
int BCP_lp_user::current_level() const     { return p->node->level; }
int BCP_lp_user::current_index() const     { return p->node->index; }
int BCP_lp_user::current_iteration() const { return p->node->iteration_count; }
double BCP_lp_user::start_time() const     { return p->start_time; }
BCP_user_data* BCP_lp_user::get_user_data() { return p->node->user_data; }

void BCP_lp_user::print(const bool ifprint, const char * format, ...) const {
  if (ifprint) {
    va_list valist;
    va_start(valist, format);
    vprintf(format, valist);
    va_end(valist);
  }
}

//#############################################################################
// Informational methods for the user
/* Methods to get/set BCP parameters on the fly */
char
BCP_lp_user::get_param(const BCP_lp_par::chr_params key) const
{ return p->par.entry(key); }
int
BCP_lp_user::get_param(const BCP_lp_par::int_params key) const
{ return p->par.entry(key); }
double
BCP_lp_user::get_param(const BCP_lp_par::dbl_params key) const
{ return p->par.entry(key); }
const BCP_string&
BCP_lp_user::get_param(const BCP_lp_par::str_params key) const
{ return p->par.entry(key); }

void BCP_lp_user::set_param(const BCP_lp_par::chr_params key, const bool val)
{ p->par.set_entry(key, val); }
void BCP_lp_user::set_param(const BCP_lp_par::chr_params key, const char val)
{ p->par.set_entry(key, val); }
void BCP_lp_user::set_param(const BCP_lp_par::int_params key, const int val)
{ p->par.set_entry(key, val); }
void BCP_lp_user::set_param(const BCP_lp_par::dbl_params key, const double val)
{ p->par.set_entry(key, val); }
void BCP_lp_user::set_param(const BCP_lp_par::str_params key, const char * val)
{ p->par.set_entry(key, val); }

//#############################################################################

void
BCP_lp_user::send_feasible_solution(const BCP_solution* sol)
{
    // send the feas sol to the tree manager
    p->msg_buf.clear();
    pack_feasible_solution(p->msg_buf, sol);
    p->msg_env->send(p->get_parent() /*tree_manager*/,
                     BCP_Msg_FeasibleSolution, p->msg_buf);

    // update the UB if necessary
    const double obj = sol->objective_value();
    if (p->ub(obj) && p->node->colgen != BCP_GenerateColumns) {
        // FIXME: If we had a flag in the node that indicates not to
        // generate cols in it and in its descendants then the dual obj
        // limit could still be set...
        p->lp_solver->setDblParam(OsiDualObjectiveLimit, obj-p->granularity());
    }
}

//#############################################################################

// These functions are member functions of the VIRTUAL class BCP_lp_user.
// Therefore any concrete class derived from BCP_lp_user can override these
// definitions. These are here as default functions.

//#############################################################################
// a few helper functions for selecting a subset of a double vector
void
BCP_lp_user::select_nonzeros(const double * first, const double * last,
                             const double etol,
                             BCP_vec<int>& nonzeros) const {
    nonzeros.reserve(last - first);
    BCP_vec<double>::const_iterator current = first;
    for ( ; current != last; ++current)
        if (CoinAbs(*current) > etol)
            nonzeros.unchecked_push_back(current - first);
}
//-----------------------------------------------------------------------------
void
BCP_lp_user::select_zeros(const double * first, const double * last,
                          const double etol,
                          BCP_vec<int>& zeros) const {
    zeros.reserve(last - first);
    BCP_vec<double>::const_iterator current = first;
    for ( ; current != last; ++current)
        if (CoinAbs(*current) <= etol)
            zeros.unchecked_push_back(current - first);
}
//-----------------------------------------------------------------------------
void
BCP_lp_user::select_positives(const double * first, const double * last,
                              const double etol,
                              BCP_vec<int>& positives) const {
    positives.reserve(last - first);
    BCP_vec<double>::const_iterator current = first;
    for ( ; current != last; ++current)
        if (*current > etol)
            positives.unchecked_push_back(current - first);
}
//-----------------------------------------------------------------------------
void
BCP_lp_user::select_fractions(const double * first, const double * last,
                              const double etol,
                              BCP_vec<int>& fractions) const {
    fractions.reserve(last - first);
    BCP_vec<double>::const_iterator current = first;
    for ( ; current != last; ++current)
        if (*current-floor(*current) > etol && ceil(*current)-*current > etol)
            fractions.unchecked_push_back(current - first);
}

//#############################################################################
//#############################################################################
// unpack the initial info for the appropriate process
void
BCP_lp_user::unpack_module_data(BCP_buffer & buf)
{
  print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
        "LP: Default unpack_module_data() executed.\n");
}

//#############################################################################

int
BCP_lp_user::process_id() const
{
    return p->get_process_id();
}

int 
BCP_lp_user::parent() const
{
    return p->get_parent();
}

void
BCP_lp_user::send_message(const int target, const BCP_buffer& buf,
                          BCP_message_tag tag)
{
    p->msg_env->send(target, tag, buf);
}    

void
BCP_lp_user::receive_message(const int sender, BCP_buffer& buf,
                             BCP_message_tag tag)
{
    p->msg_env->receive(sender, tag, buf, -1);
}

void
BCP_lp_user::broadcast_message(const BCP_process_t proc_type,
                               const BCP_buffer& buf)
{
    throw BCP_fatal_error("\
BCP_lp_user::broadcast_message: can't broadcast from an LP process.\n");
}

void
BCP_lp_user::process_message(BCP_buffer& buf)
{
    throw BCP_fatal_error("\
BCP_lp_user::process_message() invoked but not overridden!\n");
}

//#############################################################################

OsiSolverInterface * 
BCP_lp_user::initialize_solver_interface()
{
    throw BCP_fatal_error("\
BCP_lp_user::initialize_solver_interface() invoked but not overridden!\n");
    return 0;
}

//#############################################################################

void 
BCP_lp_user::initialize_int_and_sos_list(std::vector<OsiObject *>& intAndSosObjects)
{
    return;
}

//#############################################################################
void
BCP_lp_user::initialize_new_search_tree_node(const BCP_vec<BCP_var*>& vars,
                                             const BCP_vec<BCP_cut*>& cuts,
                                             const BCP_vec<BCP_obj_status>& vs,
                                             const BCP_vec<BCP_obj_status>& cs,
                                             BCP_vec<int>& var_changed_pos,
                                             BCP_vec<double>& var_new_bd,
                                             BCP_vec<int>& cut_changed_pos,
                                             BCP_vec<double>& cut_new_bd)
{
  print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
        "LP: Default initialize_new_search_tree_node() executed.\n");
}

//#############################################################################

void
BCP_lp_user::load_problem(OsiSolverInterface& osi, BCP_problem_core* core,
                          BCP_var_set& vars, BCP_cut_set& cuts)
{
  const int varnum = vars.size();
  const int cutnum = cuts.size();

  if (varnum == 0) {
    // *FIXME* : kill all the processes
    throw BCP_fatal_error("There are no vars in the description for node %i!\n",
                          current_index());
  }

  if (cutnum == 0) {
    // *FIXME* : kill all the processes
    throw BCP_fatal_error("There are no cuts in the description for node %i!\n",
                          current_index());
  }

  BCP_vec<BCP_col*> cols;
  BCP_vec<BCP_row*> rows;

  int bvarnum = core->varnum();
  int bcutnum = core->cutnum();

  if (varnum == 0 || cutnum == 0){
    throw BCP_fatal_error("No rows or no cols to create a matrix from!\n");
  }

  BCP_lp_relax* m = 0;
  if (bvarnum == 0) {
    // no core vars. doesn't matter if there're any core cuts, the starting
    // matrix is computed the same way
    cols.reserve(varnum);
    vars_to_cols(cuts, vars, cols,
                 *p->lp_result, BCP_Object_FromTreeManager, false);
    BCP_vec<double> RLB;
    BCP_vec<double> RUB;
    RLB.reserve(cutnum);
    RUB.reserve(cutnum);
    BCP_cut_set::const_iterator ci = cuts.begin();
    BCP_cut_set::const_iterator lastci = cuts.end();
    for ( ; ci != lastci; ++ci) {
      RLB.unchecked_push_back((*ci)->lb());
      RUB.unchecked_push_back((*ci)->ub());
    }
    m = new BCP_lp_relax(cols, RLB, RUB);
    purge_ptr_vector(cols);
  } else {
    if (bcutnum == 0) {
      rows.reserve(cutnum);
      cuts_to_rows(vars, cuts, rows,
                   *p->lp_result, BCP_Object_FromTreeManager, false);
      BCP_vec<double> CLB;
      BCP_vec<double> CUB;
      BCP_vec<double> OBJ;
      CLB.reserve(varnum);
      CUB.reserve(varnum);
      OBJ.reserve(varnum);
      BCP_var_set::const_iterator vi = vars.begin();
      BCP_var_set::const_iterator lastvi = vars.end();
      for ( ; vi != lastvi; ++vi) {
        CLB.unchecked_push_back((*vi)->lb());
        CUB.unchecked_push_back((*vi)->ub());
        OBJ.unchecked_push_back((*vi)->obj());
      }
      m = new BCP_lp_relax(rows, CLB, CUB, OBJ);
      purge_ptr_vector(rows);
    } else {
      // has core vars and cuts, the starting matrix is the core matrix
      m = core->matrix;
    }
  }
   
  // We have the description in p.node. Load it into the lp solver.
  // First load the core matrix
  osi.loadProblem(*m,
                  m->ColLowerBound().begin(), m->ColUpperBound().begin(),
                  m->Objective().begin(),
                  m->RowLowerBound().begin(), m->RowUpperBound().begin());

  if (bvarnum > 0 && bcutnum > 0) {
    //-----------------------------------------------------------------------
    // We have to add the 'added' stuff only if we had a core matrix
    //-----------------------------------------------------------------------
    // Add the Named and Algo cols if there are any (and if we have any cols)
    if (varnum > bvarnum && bcutnum > 0) {
      BCP_vec<BCP_var*> added_vars(vars.entry(bvarnum), vars.end());
      BCP_vec<BCP_cut*> core_cuts(cuts.begin(), cuts.entry(bcutnum));
      cols.reserve(vars.size());
      vars_to_cols(core_cuts, added_vars, cols,
                   *p->lp_result, BCP_Object_FromTreeManager, false);
      BCP_lp_add_cols_to_lp(cols, &osi);
      purge_ptr_vector(cols);
    }
    //-----------------------------------------------------------------------
    // Add the Named and Algo rows if there are any (and if we have any such
    // rows AND if there are core vars)
    if (cutnum > bcutnum) {
      BCP_vec<BCP_cut*> added_cuts(cuts.entry(bcutnum), cuts.end());
      rows.reserve(added_cuts.size());
      BCP_fatal_error::abort_on_error = false;
      try {
         cuts_to_rows(vars, added_cuts, rows,
                      *p->lp_result, BCP_Object_FromTreeManager, false);
      }
      catch (...) {
      }
      BCP_fatal_error::abort_on_error = true;
      BCP_lp_add_rows_to_lp(rows, &osi);
      purge_ptr_vector(rows);
    }
  } else {
    // Otherwise (i.e., if we had no core matrix) we just have to get rid of
    // 'm'.
    delete m;
  }
  const bool dumpcuts = get_param(BCP_lp_par::Lp_DumpNodeDescCuts);
  const bool dumpvars = get_param(BCP_lp_par::Lp_DumpNodeDescVars);
  if (dumpcuts || dumpvars) {
    printf("LP: NEW ACTIVE NODE %i DUMP START ===========================\n",
           current_index());
    if (dumpvars) {
      printf("    Core Vars (bvarnum %i)\n", bvarnum);
      for (int i = 0; i < bvarnum; ++i) {
        printf("        var %4i: bcpind %i,   status %i,   lb %f,   ub %f\n",
               i, vars[i]->bcpind(), vars[i]->status(),
               vars[i]->lb(), vars[i]->ub());
      }
    }
    if (dumpcuts) {
      printf("    Core Cuts (bcutnum %i)\n", bcutnum);
      for (int i = 0; i < bcutnum; ++i) {
        printf("        cut %4i: bcpind %i,   status %i,   lb %f,   ub %f\n",
               i, cuts[i]->bcpind(), cuts[i]->status(),
               cuts[i]->lb(), cuts[i]->ub());
      }
    }
    if (dumpvars) {
      printf("    Extra Vars (extra varnum %i)\n", varnum - bvarnum);
      for (int i = bvarnum; i < varnum; ++i) {
        printf("        var %4i: bcpind %i,   status %i,   lb %f,   ub %f\n",
               i, vars[i]->bcpind(), vars[i]->status(),
               vars[i]->lb(), vars[i]->ub());
      }
    }
    if (dumpcuts) {
      printf("    Extra Cuts (extra cutnum %i)\n", cutnum - bcutnum);
      for (int i = bcutnum; i < cutnum; ++i) {
        printf("        cut %4i: bcpind %i,   status %i,   lb %f,   ub %f\n",
               i, cuts[i]->bcpind(), cuts[i]->status(),
               cuts[i]->lb(), cuts[i]->ub());
      }
    }
    printf("LP: NEW ACTIVE NODE %i DUMP END ===========================\n",
           current_index());
  }
}

//#############################################################################
// Opportunity to reset things before optimization
void
BCP_lp_user::modify_lp_parameters(OsiSolverInterface* lp, const int changeType,
                                  bool in_strong_branching)
{
  print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
        "LP: Default prepare_for_optimization() executed.\n");

  if (changeType == 1) { // primal feas was affected
    lp->setHintParam(OsiDoDualInResolve, true, OsiHintTry);
    return;
  }
  if (changeType == 2) { // dual feas was affected
    lp->setHintParam(OsiDoDualInResolve, false, OsiHintTry);
    return;
  }
  // Neither or both, so let the solver decide
  lp->setHintParam(OsiDoDualInResolve, true, OsiHintIgnore);
}

//#############################################################################
void
BCP_lp_user::process_lp_result(const BCP_lp_result& lpres,
                               const BCP_vec<BCP_var*>& vars,
                               const BCP_vec<BCP_cut*>& cuts,
                               const double old_lower_bound,
                               double& true_lower_bound,
                               BCP_solution*& sol,
                               BCP_vec<BCP_cut*>& new_cuts,
                               BCP_vec<BCP_row*>& new_rows,
                               BCP_vec<BCP_var*>& new_vars,
                               BCP_vec<BCP_col*>& new_cols)
{
    p->user_has_lp_result_processing = false;
}

//#############################################################################
// Generating a true lower bound
double
BCP_lp_user::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)
{
    // If columns are to be generated then we can't say anything, just return
    // the current lower bound
    if (p->node->colgen != BCP_DoNotGenerateColumns_Fathom)
        return old_lower_bound;

    // Otherwise we got the examine the termination code and the objective
    // value of the LP solution
    const int tc = lpres.termcode();
    if (tc & BCP_ProvenOptimal)
        return lpres.objval();

    // The limit (the upper bound) on the dual objective is proven to be
    // reached, but the objval might not reflect this! (the LP solver may not
    // make the last iteration that pushes objval over the limit). So we return
    // a high value ourselves.
    if (tc & BCP_DualObjLimReached)
        return p->ub() + 1e-5;

    // We can't say anything in any other case
    // (BCP_ProvenPrimalInf | BCP_ProvenDualInf | BCP_PrimalObjLimReached |
    //  BCP_TimeLimit | BCP_Abandoned), not to mention that some of these are
    //  impossible. Just return the current bound.
    return old_lower_bound;
}

//#############################################################################
// Feasibility testing
BCP_solution*
BCP_lp_user::test_feasibility(const BCP_lp_result& lpres,
                              const BCP_vec<BCP_var*>& vars,
                              const BCP_vec<BCP_cut*>& cuts)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default test_feasibility() executed.\n");

    const double etol = p->param(BCP_lp_par::IntegerTolerance);
    BCP_feasibility_test test =
        static_cast<BCP_feasibility_test>(p->param(BCP_lp_par::FeasibilityTest));

    BCP_solution* sol = NULL;
    switch (test){
    case BCP_Binary_Feasible:   sol = test_binary(lpres, vars, etol); break;
    case BCP_Integral_Feasible: sol = test_integral(lpres, vars, etol); break;
    case BCP_FullTest_Feasible: sol = test_full(lpres, vars, etol); break;
    default:
        throw BCP_fatal_error("LP: unknown test_feasibility() rule.\n");
    }

    return sol;
}

//-----------------------------------------------------------------------------

BCP_solution_generic*
BCP_lp_user::test_binary(const BCP_lp_result& lpres,
                         const BCP_vec<BCP_var*>& vars,
                         const double etol) const
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default test_binary() executed.\n");

    // Do anything only if the termination code is sensible
    const int tc = lpres.termcode();
    if (! (tc & BCP_ProvenOptimal))
        return 0;

    const double * x = lpres.x();
    const int varnum = vars.size();
    const double etol1 = 1 - etol;
    int i;

    for (i = 0 ; i < varnum; ++i) {
        const double xi = x[i];
        if (xi > etol && xi < etol1)
            return(NULL);
    }

    // This solution does not own the pointers to the variables
    BCP_solution_generic* sol = new BCP_solution_generic(false); 
    double obj = 0;
    for (i = 0 ; i < varnum; ++i) {
        if (x[i] > etol1) {
            sol->add_entry(vars[i], 1.0);
            obj += vars[i]->obj();
        }
    }
    sol->_objective = obj;
    return sol;
}
//-----------------------------------------------------------------------------
BCP_solution_generic*
BCP_lp_user::test_integral(const BCP_lp_result& lpres,
                           const BCP_vec<BCP_var*>& vars,
                           const double etol) const
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default test_integral() executed.\n");

    // Do anything only if the termination code is sensible
    const int tc = lpres.termcode();
    if (! (tc & BCP_ProvenOptimal))
        return 0;

    const double * x = lpres.x();
    const int varnum = vars.size();
    const double etol1 = 1 - etol;
    int i;

    for (i = 0 ; i < varnum; ++i) {
        const double val = x[i] - floor(x[i]);
        if (val > etol && val < etol1)
            return(NULL);
    }
  
    // This solution does not own the pointers to the variables
    BCP_solution_generic* sol = new BCP_solution_generic(false);
    double obj = 0;
    for (i = 0 ; i < varnum; ++i) {
        const double val = floor(x[i] + 0.5);
        if (val != 0.0) { // it's OK to test the double against 0,
                          // since we took the floor()
            sol->add_entry(vars[i], val);
            obj += val * vars[i]->obj();
        }
    }
    sol->_objective = obj;
    return sol;
}
//-----------------------------------------------------------------------------
BCP_solution_generic*
BCP_lp_user::test_full(const BCP_lp_result& lpres,
                       const BCP_vec<BCP_var*>& vars,
                       const double etol) const
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default test_full() executed.\n");

    // Do anything only if the termination code is sensible
    const int tc = lpres.termcode();
    if (! (tc & BCP_ProvenOptimal))
        return 0;

    const double * x = lpres.x();
    const int varnum = vars.size();
    const double etol1 = 1 - etol;
    int i;

    for (i = 0 ; i < varnum; ++i) {
        switch (vars[i]->var_type()){
        case BCP_BinaryVar:
            {
                const double val = x[i];
                if (val > etol && val < etol1)
                    return(NULL);
            }
            break;
        case BCP_IntegerVar:
            {
                const double val = x[i] - floor(x[i]);
                if (val > etol && val < etol1)
                    return(NULL);
            }
            break;
        case BCP_ContinuousVar:
            break;
        }
    }

    // This solution does not own the pointers to the variables
    BCP_solution_generic* sol = new BCP_solution_generic(false);
    double obj = 0;
    for (i = 0 ; i < varnum; ++i) {
        // round the integer vars
        const double val = (vars[i]->var_type() == BCP_ContinuousVar) ?
            x[i] : floor(x[i] + 0.5);
        if (val < -etol || val > etol) {
            sol->add_entry(vars[i], val);
            obj += val * vars[i]->obj();
        }
    }
    sol->_objective = obj;
    return sol;
}

//#############################################################################
// Heuristic solution generation
BCP_solution*
BCP_lp_user::generate_heuristic_solution(const BCP_lp_result& lpres,
                                         const BCP_vec<BCP_var*>& vars,
                                         const BCP_vec<BCP_cut*>& cuts)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default generate_heuristic_solution() executed.\n");
    return NULL;
}
//#############################################################################
// Feasible solution packing
void
BCP_lp_user::pack_feasible_solution(BCP_buffer& buf, const BCP_solution* sol)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default pack_feasible_solution() executed.\n");

    const BCP_solution_generic* gensol =
        dynamic_cast<const BCP_solution_generic*>(sol);
    const BCP_vec<BCP_var*> solvars = gensol->_vars;
    const BCP_vec<double> values = gensol->_values;
    const int size = solvars.size();
    buf.pack(size);
    for (int i = 0; i < size; ++i) {
        buf.pack(values[i]);
        p->pack_var(*solvars[i]);
    }
    buf.pack(gensol->objective_value());
}

//#############################################################################
// pack message for CG/CP

void
BCP_lp_user::pack_primal_solution(BCP_buffer& buf,
                                  const BCP_lp_result& lpres,
                                  const BCP_vec<BCP_var*>& vars,
                                  const BCP_vec<BCP_cut*>& cuts)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default pack_for_cg() executed.\n");

    BCP_vec<int> coll;

    const double petol = lpres.primalTolerance();
    const double * x = lpres.x();
    const int varnum = vars.size();

    switch (p->param(BCP_lp_par::InfoForCG)){
    case BCP_PrimalSolution_Nonzeros:
        select_nonzeros(x, x + varnum, petol, coll);
        buf.set_msgtag(BCP_Msg_ForCG_PrimalNonzeros);
        break;
    case BCP_PrimalSolution_Fractions:
        select_fractions(x, x+varnum, petol, coll);
        buf.set_msgtag(BCP_Msg_ForCG_PrimalFractions);
        break;
    case BCP_PrimalSolution_Full:
        coll.reserve(varnum);
        for (int i = 0; i < varnum; ++i) {
            coll.unchecked_push_back(i);
        }
        buf.set_msgtag(BCP_Msg_ForCG_PrimalFull);
        break;
    default:
        throw BCP_fatal_error("Incorrect msgtag in pack_for_cg() !\n");
    }

    const int size = coll.size();
    buf.pack(size);
    if (size > 0){
        BCP_vec<int>::const_iterator pos = coll.begin() - 1;
        const BCP_vec<int>::const_iterator last_pos = coll.end();
        while (++pos != last_pos) {
            buf.pack(x[*pos]);
            p->pack_var(*vars[*pos]);
        }
    }
}

//=============================================================================
// pack message for VG/VP

void
BCP_lp_user::pack_dual_solution(BCP_buffer& buf,
                                const BCP_lp_result& lpres,
                                const BCP_vec<BCP_var*>& vars,
                                const BCP_vec<BCP_cut*>& cuts)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default pack_for_vg() executed.\n");

    BCP_vec<int> coll;

    const double detol = lpres.dualTolerance();
    const double * pi = lpres.pi();
    const int cutnum = cuts.size();

    switch (p->param(BCP_lp_par::InfoForVG)){
    case BCP_DualSolution_Nonzeros:
        select_nonzeros(pi, pi+cutnum, detol, coll);
        buf.set_msgtag(BCP_Msg_ForVG_DualNonzeros);
        break;
    case BCP_DualSolution_Full:
        coll.reserve(cutnum);
        for (int i = 0; i < cutnum; ++i) {
            coll.unchecked_push_back(i);
        }
        buf.set_msgtag(BCP_Msg_ForVG_DualFull);
        break;
    default:
        throw BCP_fatal_error("Incorrect msgtag in pack_lp_solution() !\n");
    }

    const int size = coll.size();
    buf.pack(size);
    if (size > 0){
        BCP_vec<int>::const_iterator pos = coll.begin() - 1;
        const BCP_vec<int>::const_iterator last_pos = coll.end();
        while (++pos != last_pos) {
            buf.pack(pi[*pos]);
            p->pack_cut(*cuts[*pos]);
        }
    }
}

//#############################################################################
// lp solution displaying
void
BCP_lp_user::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)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default display_lp_solution() executed.\n");

    if (final_lp_solution) {
        if (! p->param(BCP_lp_par::LpVerb_FinalRelaxedSolution))
            return;
        print(true,"  LP : Displaying LP solution (FinalRelaxedSolution) :\n");
    } else {
        if (! p->param(BCP_lp_par::LpVerb_RelaxedSolution))
            return;
        print(true,"  LP : Displaying LP solution (RelaxedSolution) :\n");
    }

    const double ietol = p->param(BCP_lp_par::IntegerTolerance);

    print(true, "  LP : Displaying solution :\n");
    BCP_vec<int> coll;
    const double * x = lpres.x();
    select_nonzeros(x, x+vars.size(), ietol, coll);
    const int size = coll.size();
    for (int i = 0; i < size; ++i) {
        const int ind = coll[i];
        vars[ind]->display(x[ind]);
    }
}

//#############################################################################
// restoring feasibility
void
BCP_lp_user::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)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default restore_feasibility() executed.\n");
}

//#############################################################################
//#############################################################################

void
BCP_lp_user::cuts_to_rows(const BCP_vec<BCP_var*>& vars, // on what to expand
                          BCP_vec<BCP_cut*>& cuts,       // what to expand
                          BCP_vec<BCP_row*>& rows,       // the expanded rows
                          // things that the user can use for lifting cuts
                          // if allowed
                          const BCP_lp_result& lpres,
                          BCP_object_origin origin, bool allow_multiple)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default cuts_to_rows() executed.\n");
    throw BCP_fatal_error("cuts_to_rows() missing.\n");
}
//-----------------------------------------------------------------------------
void
BCP_lp_user::vars_to_cols(const BCP_vec<BCP_cut*>& cuts, // on what to expand
                          BCP_vec<BCP_var*>& vars,       // what to expand
                          BCP_vec<BCP_col*>& cols,       // the expanded cols
                          // things that the user can use for lifting vars
                          // if allowed
                          const BCP_lp_result& lpres,
                          BCP_object_origin origin, bool allow_multiple)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
        "LP: Default vars_to_cols() executed.\n");
    throw BCP_fatal_error("vars_to_cols() missing.\n");
}

//#############################################################################

void
BCP_lp_user::generate_cuts_in_lp(const BCP_lp_result& lpres,
                                 const BCP_vec<BCP_var*>& vars,
                                 const BCP_vec<BCP_cut*>& cuts,
                                 BCP_vec<BCP_cut*>& new_cuts,
                                 BCP_vec<BCP_row*>& new_rows)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default generate_cuts_in_lp() executed.\n");
}
//-----------------------------------------------------------------------------
void
BCP_lp_user::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)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default generate_vars_in_lp() executed.\n");
}
//-----------------------------------------------------------------------------
BCP_object_compare_result
BCP_lp_user::compare_cuts(const BCP_cut* c0, const BCP_cut* c1)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default compare_cuts() executed.\n");
    return BCP_DifferentObjs;
}
//-----------------------------------------------------------------------------
BCP_object_compare_result
BCP_lp_user::compare_vars(const BCP_var* v0, const BCP_var* v1)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default compare_vars() executed.\n");
    return BCP_DifferentObjs;
}

//#############################################################################

void
BCP_lp_user::select_vars_to_delete(const BCP_lp_result& lpres,
                                   const BCP_vec<BCP_var*>& vars,
                                   const BCP_vec<BCP_cut*>& cuts,
                                   const bool before_fathom,
                                   BCP_vec<int>& deletable)
{
    if (before_fathom && p->param(BCP_lp_par::NoCompressionAtFathom))
        return;
    const int varnum = vars.size();
    deletable.reserve(varnum);
    for (int i = p->core->varnum(); i < varnum; ++i) {
        BCP_var *var = vars[i];
        if (var->is_to_be_removed() ||
            (! var->is_non_removable() && var->lb() == 0 && var->ub() == 0)) {
            deletable.unchecked_push_back(i);
        }
    }
}

//=============================================================================

void
BCP_lp_user::select_cuts_to_delete(const BCP_lp_result& lpres,
                                   const BCP_vec<BCP_var*>& vars,
                                   const BCP_vec<BCP_cut*>& cuts,
                                   const bool before_fathom,
                                   BCP_vec<int>& deletable)
{
    if (before_fathom && p->param(BCP_lp_par::NoCompressionAtFathom))
        return;
    const int cutnum = cuts.size();
    const int ineff_to_delete = p->param(BCP_lp_par::IneffectiveBeforeDelete);
    const double lb = lpres.objval();
    const BCP_vec<double>& lb_at_cutgen = p->node->lb_at_cutgen;
    deletable.reserve(cutnum);
    for (int i = p->core->cutnum(); i < cutnum; ++i) {
        BCP_cut *cut = cuts[i];
        if (cut->is_to_be_removed() ||
            (! cut->is_non_removable() &&
             cut->effective_count() <= -ineff_to_delete &&
             lb_at_cutgen[i] < lb - 0.0001)) {
            deletable.unchecked_push_back(i);
        }
    }
}

//#############################################################################

void
BCP_lp_user::logical_fixing(const BCP_lp_result& lpres,
                            const BCP_vec<BCP_var*>& vars,
                            const BCP_vec<BCP_cut*>& cuts,
                            const BCP_vec<BCP_obj_status>& var_status,
                            const BCP_vec<BCP_obj_status>& cut_status,
                            const int var_bound_changes_since_logical_fixing,
                            BCP_vec<int>& changed_pos,
                            BCP_vec<double>& new_bd)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default logical_fixing() executed.\n");
}

//#############################################################################
void
BCP_lp_user::reduced_cost_fixing(const double* dj, const double* x,
                                 const double gap,
                                 BCP_vec<BCP_var*>& vars, int& newly_changed)
{
    OsiBabSolver* babSolver = getOsiBabSolver();
    if (babSolver && !babSolver->reducedCostsAccurate())
        return;

    newly_changed = 0;
    const bool atZero = get_param(BCP_lp_par::DoReducedCostFixingAtZero);
    const bool atAny = get_param(BCP_lp_par::DoReducedCostFixingAtAnything);

    if (! atZero && ! atAny)
        return;

    double petol = 0.0;
    p->lp_solver->getDblParam(OsiPrimalTolerance, petol);
    double detol = 0.0;
    p->lp_solver->getDblParam(OsiDualTolerance, detol);

    // If the gap is negative that means that we are above the limit, so
    // don't do anything.
    if (gap < 0)
        return;

    const int varnum = vars.size();
    BCP_vec<int> changed_indices;
    changed_indices.reserve(varnum);
    BCP_vec<double> changed_bounds;
    changed_bounds.reserve(2 * varnum);

    // Note that when this function is called, we must have a dual
    // feasible dual solution. Therefore we can use the lagrangean
    // relaxation to fix variables.

    // *FIXME* : If we knew that there are integral vars only, then
    // we could leave out the test for BCP_ContinuousVar...
    for (int i = 0; i < varnum; ++i) {
        BCP_var* var = vars[i];
        if (! var->is_fixed() && var->var_type() != BCP_ContinuousVar){
            if (dj[i] > detol) {
                const double lb = var->lb();
                const double new_ub = lb + floor(gap / dj[i]);
                if (new_ub < var->ub() && (atAny || CoinAbs(x[i])<petol) ) {
                    vars[i]->set_ub(new_ub);
                    changed_indices.unchecked_push_back(i);
                    changed_bounds.unchecked_push_back(lb);
                    changed_bounds.unchecked_push_back(new_ub);
                }
            } else if (dj[i] < -detol) {
                const double ub = var->ub();
                const double new_lb = ub - floor(gap / (-dj[i]));
                if (new_lb > var->lb() && (atAny || CoinAbs(x[i])<petol) ) {
                    vars[i]->set_lb(new_lb);
                    changed_indices.unchecked_push_back(i);
                    changed_bounds.unchecked_push_back(new_lb);
                    changed_bounds.unchecked_push_back(ub);
                }
            }
        }
    }
    newly_changed = changed_indices.size();
    if (newly_changed > 0) {
        p->lp_solver->setColSetBounds(changed_indices.begin(),
                                      changed_indices.end(),
                                      changed_bounds.begin());
    }
}

//#############################################################################
//#############################################################################

int
BCP_lp_user::try_to_branch(OsiBranchingInformation& branchInfo,
                           OsiSolverInterface* solver,
                           OsiChooseVariable* choose,
                           OsiBranchingObject*& branchObject,
                           bool allowVarFix)
{
    int returnStatus = 0;
    int numUnsatisfied = choose->setupList(&branchInfo, true);
    choose->setBestObjectIndex(-1);
    if (numUnsatisfied > 0) {
        if (choose->numberOnList() == 0) {
            // Nothing left for strong branching to evaluate
            if (choose->numberOnList() > 0 || choose->numberStrong() == 0) {
                // There is something on the list
                choose->setBestObjectIndex(choose->candidates()[0]);
            } else {
                // There is nothing on the list
                numUnsatisfied = choose->setupList(&branchInfo, false);
                if (numUnsatisfied > 0) {
                    choose->setBestObjectIndex(choose->candidates()[0]);
                }
            }
        } else {
            // Do the strong branching
            int ret = choose->chooseVariable(solver, &branchInfo, allowVarFix);
            /* Check if SB has fixed anything, and if so, apply the same
               changes to the vars */
            const double * clb = solver->getColLower();
            const double * cub = solver->getColUpper();
            BCP_vec<BCP_var*>& vars = p->node->vars;
            for (int i = numUnsatisfied - 1; i >= 0; --i) {
                const int ind =
                  solver->object(choose->candidates()[i])->columnNumber();
                if (ind >= 0) {
                  assert(vars[ind]->lb() <= clb[ind]);
                  assert(vars[ind]->ub() >= cub[ind]);
                  vars[ind]->set_lb_ub(clb[ind], cub[ind]);
                }
            }
                
            /* update number of strong iterations etc
            model->incrementStrongInfo(choose->numberStrongDone(),
                                       choose->numberStrongIterations(),
                                       ret==-1 ? 0:choose->numberStrongFixed(),
                                       ret==-1);
            */
            if (ret > 1) {
                // has fixed some
                returnStatus = -1;
            } else if (ret == -1) {
                // infeasible
                returnStatus = -2;
            } else if (ret == 0) {
                // normal
                returnStatus = 0;
                numUnsatisfied = 1;
            } else {
                // ones on list satisfied - double check
                numUnsatisfied = choose->setupList(&branchInfo, false);
                if (numUnsatisfied > 0) {
                    choose->setBestObjectIndex(choose->candidates()[0]);
                }
            }
        }
    }
    if (! returnStatus) {
        if (numUnsatisfied > 0) {
            // create branching object
            /* FIXME: how the objects are created? */
            const OsiObject * obj = solver->object(choose->bestObjectIndex());
            branchObject = obj->createBranch(solver,
                                             &branchInfo,
                                             obj->whichWay());
        }
    }

    return returnStatus;
}

//#############################################################################

BCP_branching_decision BCP_lp_user::
select_branching_candidates(const BCP_lp_result& lpres,
                            const BCP_vec<BCP_var*>& vars,
                            const BCP_vec<BCP_cut*>& cuts,
                            const BCP_lp_var_pool& local_var_pool,
                            const BCP_lp_cut_pool& local_cut_pool,
                            BCP_vec<BCP_lp_branching_object*>& cands,
                            bool force_branch)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default select_branching_candidates() executed.\n");

    if (lpres.termcode() & BCP_Abandoned) {
      // *THINK*: maybe we should branch through...
      print(true, "\
LP: ############ LP solver abandoned. Branching through...\n");
    }

    // *THINK* : this branching object selection is very primitive
    // *THINK* : should check for tail-off, could check for branching cuts, etc
    if (! force_branch &&
        (local_var_pool.size() > 0 || local_cut_pool.size() > 0))
        return BCP_DoNotBranch;

    OsiSolverInterface* lp = p->lp_solver;

    /* The last true: make a copy in brInfo of the sol of solver
       instead of just pointing into the solver's copy */
    OsiBranchingInformation brInfo(lp, true, true);
    lp->getDblParam(OsiDualObjectiveLimit, brInfo.cutoff_);
    brInfo.integerTolerance_ = p->param(BCP_lp_par::IntegerTolerance);
    brInfo.timeRemaining_ = get_param(BCP_lp_par::MaxRunTime) - CoinCpuTime();
    brInfo.numberSolutions_ = 0; /*FIXME*/
    brInfo.numberBranchingSolutions_ = 0; /*FIXME numBranchingSolutions_;*/
    brInfo.depth_ = current_level();

    OsiChooseStrong* strong = new OsiChooseStrong(lp);
    strong->setNumberBeforeTrusted(5); // the default in Cbc
    strong->setNumberStrong(p->param(BCP_lp_par::StrongBranchNum));
    strong->setTrustStrongForSolution(false);
    strong->setShadowPriceMode(0);

    OsiChooseVariable * choose = strong;
    OsiBranchingObject* brObj = NULL;

    bool allowVarFix = true;
    /*
      <li>  0: A branching object has been installed
      <li> -1: A monotone object was discovered
      <li> -2: An infeasible object was discovered
    */
    int brResult =
        try_to_branch(brInfo, lp, choose, brObj, allowVarFix);
    const int bestWhichWay = choose->bestWhichWay();

#if 0
    /* FIXME:before doing anything check if we have found a new solution */
    if (choose->goodSolution()
        &&model->problemFeasibility()->feasible(model,-1)>=0) {
        // yes
        double objValue = choose->goodObjectiveValue();
        model->setBestSolution(CBC_STRONGSOL,
                               objValue,
                               choose->goodSolution()) ;
        model->setLastHeuristic(NULL);
        model->incrementUsed(choose->goodSolution());
        choose->clearGoodSolution();
    }
#endif

    delete choose;

    switch (brResult) {
    case -2:
        // when doing strong branching a candidate has proved that the
        // problem is infeasible
        delete brObj;
        return BCP_DoNotBranch_Fathomed;
    case -1:
        // OsiChooseVariable::chooseVariable() returned 2, 3, or 4
        if (!brObj) {
            // just go back and resolve
            return BCP_DoNotBranch;
        }
        // otherwise might as well branch. The forced variable is
        // unlikely to jump up one more (though who knows...)
        break;
    case 0:
        if (!brObj) {
            // nothing got fixed, yet couldn't find something to branch on
            throw BCP_fatal_error("BM: Couldn't branch!\n");
        }
        // we've got a branching object
        break;
    default:
        throw BCP_fatal_error("\
BCP: BCP_lp_user::try_to_branch returned with unknown return code.\n");
    }

    // If there were some fixings (brResult < 0) then propagate them where
    // needed
    if (allowVarFix && brResult < 0) {
        const double* clb = lp->getColLower();
        const double* cub = lp->getColUpper();
        /* There may or may not have been changes, but faster to set then to
           test... */
        BCP_vec<BCP_var*>& vars = p->node->vars;
        int numvars = vars.size();
        for (int i = 0; i < numvars; ++i) {
            vars[i]->change_bounds(clb[i], cub[i]);
        }
    }
    
    // all possibilities are 2-way branches
    int order[2] = {0, 1};
    if (bestWhichWay == 1) {
        order[0] = 1;
        order[1] = 0;
    }

    // Now interpret the result (at this point we must have a brObj
    OsiIntegerBranchingObject* intBrObj =
        dynamic_cast<OsiIntegerBranchingObject*>(brObj);
    if (intBrObj) {
        BCP_lp_integer_branching_object o(intBrObj);
        cands.push_back(new BCP_lp_branching_object(o, order));
    }
    OsiSOSBranchingObject* sosBrObj =
        dynamic_cast<OsiSOSBranchingObject*>(brObj);
    if (sosBrObj) {
        BCP_lp_sos_branching_object o(sosBrObj);
        cands.push_back(new BCP_lp_branching_object(lp, o, order));
    }

    delete brObj;
    
    if (cands.size() == 0) {
        throw BCP_fatal_error("\
 LP : No var/cut in pool but couldn't select branching object.\n");
    }
    return BCP_DoBranch;
}

//-----------------------------------------------------------------------------
void
BCP_lp_user::append_branching_vars(const double* x,
                                   const BCP_vec<BCP_var*>& vars,
                                   const BCP_vec<int>& select_pos,
                                   BCP_vec<BCP_lp_branching_object*>& cans)
{
    BCP_vec<int>::const_iterator spi;
    BCP_vec<double> vbd(4, 0.0);
    BCP_vec<int> vpos(1, 0);

    for (spi = select_pos.begin(); spi != select_pos.end(); ++spi) {
        const int pos = *spi;
        vpos[0] = pos;
        vbd[0] = vars[pos]->lb();
        vbd[1] = floor(x[pos]);
        vbd[2] = ceil(x[pos]);
        vbd[3] = vars[pos]->ub();
        cans.push_back
            (new  BCP_lp_branching_object(2,
                                          0, 0, /* vars/cuts_to_add */
                                          &vpos, 0, &vbd, 0, /* forced parts */
                                          0, 0, 0, 0 /* implied parts */));
    }
}

//-----------------------------------------------------------------------------
BCP_branching_object_relation BCP_lp_user::
compare_branching_candidates(BCP_presolved_lp_brobj* new_presolved,
                             BCP_presolved_lp_brobj* old_presolved)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default compare_presolved_branching_objects() executed.\n");

    // change the objvals according to the termcodes (in order to be able to
    // make decisions based on objval, no matter what the termcode is).
    new_presolved->fake_objective_values(p->lp_result->objval());

    // If using this branching object we can fathom all children, then choose
    // this object
    if (new_presolved->fathomable(p->ub() - p->granularity()))
        return(BCP_NewPresolvedIsBetter_BranchOnIt);

    // If this is the first, keep it
    if (! old_presolved)
        return(BCP_NewPresolvedIsBetter);

    // If something had gone wrong with at least one descendant in the new then
    // we prefer to keep the old.
    if (new_presolved->had_numerical_problems())
        return(BCP_OldPresolvedIsBetter);

    // OK, so all descendants finished just fine. Do whatever the parameter
    // says

    if (p->param(BCP_lp_par::BranchingObjectComparison)&BCP_Comparison_Objval){
        const double objetol = 1e-7;
      
        BCP_vec<double> new_obj;
        new_presolved->get_lower_bounds(new_obj);
        std::sort(new_obj.begin(), new_obj.end());
        const int new_not_fathomed =
            new_obj.index(std::lower_bound(new_obj.begin(), new_obj.end(),
                                           BCP_DBL_MAX / 10));

        BCP_vec<double> old_obj;
        old_presolved->get_lower_bounds(old_obj);
        std::sort(old_obj.begin(), old_obj.end());
        const int old_not_fathomed =
            old_obj.index(std::lower_bound(old_obj.begin(), old_obj.end(),
                                           BCP_DBL_MAX / 10));

        if (new_not_fathomed < old_not_fathomed)
            return BCP_NewPresolvedIsBetter;

        if (new_not_fathomed > old_not_fathomed)
            return BCP_OldPresolvedIsBetter;

        BCP_vec<double>::const_iterator new_first = new_obj.begin();
        BCP_vec<double>::const_iterator new_last = new_obj.end();
        BCP_vec<double>::const_iterator old_first = old_obj.begin();
        BCP_vec<double>::const_iterator old_last = old_obj.end();
   
        switch(p->param(BCP_lp_par::BranchingObjectComparison)){
        case BCP_LowestAverageObjval:
        case BCP_HighestAverageObjval:
            {
                double newavg = 0;
                for ( ; new_first != new_last; ++new_first) {
                    if (*new_first < BCP_DBL_MAX / 10)
                        newavg += *new_first;
                }
                newavg /= new_not_fathomed;
                double oldavg = 0;
                for ( ; old_first != old_last; ++old_first) {
                    if (*old_first < BCP_DBL_MAX / 10)
                        oldavg += *old_first;
                }
                oldavg /= old_not_fathomed;
                const bool high =
                    (p->param(BCP_lp_par::BranchingObjectComparison)
                     == BCP_HighestAverageObjval);
                return (high && newavg > oldavg) || (!high && newavg < oldavg)?
                    BCP_NewPresolvedIsBetter : BCP_OldPresolvedIsBetter;
            }

        case BCP_LowestLowObjval:
            for ( ; new_first != new_last && old_first != old_last;
                  ++new_first, ++old_first) {
                if (*new_first < *old_first - objetol)
                    return BCP_NewPresolvedIsBetter;
                if (*new_first > *old_first + objetol)
                    return BCP_OldPresolvedIsBetter;
            }
            return old_first == old_last ?
                BCP_OldPresolvedIsBetter : BCP_NewPresolvedIsBetter;

        case BCP_HighestLowObjval:
            for ( ; new_first != new_last && old_first != old_last;
                  ++new_first, ++old_first) {
                if (*new_first > *old_first + objetol)
                    return BCP_NewPresolvedIsBetter;
                if (*new_first < *old_first - objetol)
                    return BCP_OldPresolvedIsBetter;
            }
            return old_first == old_last ?
                BCP_OldPresolvedIsBetter : BCP_NewPresolvedIsBetter;

        case BCP_LowestHighObjval:
            while (new_first != new_last && old_first != old_last) {
                --new_last;
                --old_last;
                if (*new_last < *old_last - objetol)
                    return BCP_NewPresolvedIsBetter;
                if (*new_last > *old_last + objetol)
                    return BCP_OldPresolvedIsBetter;
            }
            return old_first == old_last ?
                BCP_OldPresolvedIsBetter : BCP_NewPresolvedIsBetter;

        case BCP_HighestHighObjval:
            while (new_first != new_last && old_first != old_last) {
                --new_last;
                --old_last;
                if (*new_last > *old_last + objetol)
                    return BCP_NewPresolvedIsBetter;
                if (*new_last < *old_last - objetol)
                    return BCP_OldPresolvedIsBetter;
            }
            return old_first == old_last ?
                BCP_OldPresolvedIsBetter : BCP_NewPresolvedIsBetter;
        default:
            throw BCP_fatal_error("\
Unknown branching object comparison rule.\n");
        }
    }

    // shouldn't ever get here
    throw BCP_fatal_error("No branching object comparison rule is chosen.\n");

    // fake return
    return BCP_NewPresolvedIsBetter;
}

//-----------------------------------------------------------------------------
void
BCP_lp_user::set_actions_for_children(BCP_presolved_lp_brobj* best)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default set_actions_for_children() executed.\n");

    // by default every action is set to BCP_ReturnChild
    if (p->node->dive == BCP_DoNotDive)
        return;
   
    BCP_vec<BCP_child_action>& action = best->action();
    int i, ind;

    switch (p->param(BCP_lp_par::ChildPreference)){
    case BCP_PreferDiveDown:
        ind = 0;
        break;

    case BCP_PreferDiveUp:
        ind = best->candidate()->child_num - 1;
        break;
                
    case BCP_PreferChild_LowBound:
        // NOTE: if the lowest objval child is fathomed then everything is
        for (ind = 0, i = best->candidate()->child_num - 1; i > 0; --i){
            if (best->lpres(i).objval() < best->lpres(ind).objval())
                ind = i;
        }
        break;

    case BCP_PreferChild_HighBound:
        // NOTE: this selects the highest objval child NOT FATHOMED, thus if
        // the highest objval child is fathomed then everything is
        for (ind = 0, i = best->candidate()->child_num - 1; i > 0; --i) {
            if (! p->over_ub(best->lpres(i).objval()) &&
                best->lpres(i).objval() < best->lpres(ind).objval())
                ind = i;
        }
        break;
    default:
        // *THINK* : fractional branching
        throw BCP_fatal_error("LP : Unrecognized child preference.\n");
    }
    // mark the ind-th to be kept (if not over ub)
    if (! p->over_ub(best->lpres(ind).objval()))
        action[ind] = BCP_KeepChild;
}

//#############################################################################

void
BCP_lp_user::set_user_data_for_children(BCP_presolved_lp_brobj* best, 
                                        const int selected)
{
    using_deprecated_set_user_data_for_children = true;
    set_user_data_for_children(best);
    print(using_deprecated_set_user_data_for_children,
          "\
*** WARNING *** WARNING *** WARNING *** WARNING *** WARNING *** WARNING ***\n\
You have overridden\n\
  BCP_lp_user::set_user_data_for_children(BCP_presolved_lp_brobj* best)\n\
which is a deprecated virtual function. Please override\n\
  BCP_lp_user::set_user_data_for_children(BCP_presolved_lp_brobj* best,\n\
                                          const int selected)\n\
instead. The old version will go away, your code will still compile, but it\n\
will not do what you intend it to be doing.\n\
*** WARNING *** WARNING *** WARNING *** WARNING *** WARNING *** WARNING ***\n"\
               );
}

//#############################################################################

void
BCP_lp_user::set_user_data_for_children(BCP_presolved_lp_brobj* best)
{
    using_deprecated_set_user_data_for_children = false;
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default set_user_data_for_children() executed.\n");
}

//#############################################################################
// purging the slack cut pool (candidates for branching on cut)
void
BCP_lp_user::purge_slack_pool(const BCP_vec<BCP_cut*>& slack_pool,
                              BCP_vec<int>& to_be_purged)
{
    print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
          "LP: Default purge_slack_pool() executed.\n");

    switch (p->param(BCP_lp_par::SlackCutDiscardingStrategy)) {
    case BCP_DiscardSlackCutsAtNewNode:
        if (current_iteration() != 1)
            break;
    case BCP_DiscardSlackCutsAtNewIteration:
        {
            const int size = slack_pool.size();
            if (size > 0) {
                to_be_purged.reserve(size);
                for (int i = 0; i < size; ++i)
                    to_be_purged.unchecked_push_back(i);
            }
        }
        break;
    }
}

//#############################################################################