#include <string.h>
#include <float.h>
#include "commonlib.h"
#include "lp_lib.h"
#include "lp_scale.h"
#include "lp_report.h"
#include "lp_simplex.h"
#include "lp_mipbb.h"
#ifdef FORTIFY
# include "lp_fortify.h"
#endif
STATIC BBrec *create_BB(lprec *lp, BBrec *parentBB, MYBOOL dofullcopy)
{
BBrec *newBB;
newBB = (BBrec *) calloc(1, sizeof(*newBB));
if(newBB != NULL) {
if(parentBB == NULL) {
allocREAL(lp, &newBB->upbo, lp->sum + 1, FALSE);
allocREAL(lp, &newBB->lowbo, lp->sum + 1, FALSE);
MEMCOPY(newBB->upbo, lp->orig_upbo, lp->sum + 1);
MEMCOPY(newBB->lowbo, lp->orig_lowbo, lp->sum + 1);
}
else if(dofullcopy) {
allocREAL(lp, &newBB->upbo, lp->sum + 1, FALSE);
allocREAL(lp, &newBB->lowbo, lp->sum + 1, FALSE);
MEMCOPY(newBB->upbo, parentBB->upbo, lp->sum + 1);
MEMCOPY(newBB->lowbo, parentBB->lowbo, lp->sum + 1);
}
else {
newBB->upbo = parentBB->upbo;
newBB->lowbo = parentBB->lowbo;
}
newBB->contentmode = dofullcopy;
newBB->lp = lp;
newBB->parent = parentBB;
}
return( newBB );
}
STATIC BBrec *push_BB(lprec *lp, BBrec *parentBB, int varno, int vartype, int varcus)
{
BBrec *newBB;
if(parentBB == NULL)
parentBB = lp->bb_bounds;
newBB = create_BB(lp, parentBB, FALSE);
if(newBB != NULL) {
newBB->varno = varno;
newBB->vartype = vartype;
newBB->lastvarcus = varcus;
incrementUndoLadder(lp->bb_lowerchange);
newBB->LBtrack++;
incrementUndoLadder(lp->bb_upperchange);
newBB->UBtrack++;
if((parentBB != NULL) && (parentBB->lastrcf > 0)) {
MYBOOL isINT;
int k, ii, nfixed = 0, ntighten = 0;
REAL deltaUL;
for(k = 1; k <= lp->nzdrow[0]; k++) {
ii = lp->nzdrow[k];
#ifdef UseMilpSlacksRCF
isINT = FALSE;
#else
if(ii <= lp->rows)
continue;
isINT = is_int(lp, ii-lp->rows);
#endif
#ifndef UseMilpExpandedRCF
if(!isINT)
continue;
#endif
switch(abs(rcfbound_BB(newBB, ii, isINT, &deltaUL, NULL))) {
case LE: SETMIN(deltaUL, newBB->upbo[ii]);
SETMAX(deltaUL, newBB->lowbo[ii]);
modifyUndoLadder(lp->bb_upperchange, ii, newBB->upbo, deltaUL);
break;
case GE: SETMAX(deltaUL, newBB->lowbo[ii]);
SETMIN(deltaUL, newBB->upbo[ii]);
modifyUndoLadder(lp->bb_lowerchange, ii, newBB->lowbo, deltaUL);
break;
default: continue;
}
if(newBB->upbo[ii] == newBB->lowbo[ii])
nfixed++;
else
ntighten++;
}
if(lp->bb_trace) {
report(lp, DETAILED,
"push_BB: Used reduced cost to fix %d variables and tighten %d bounds\n",
nfixed, ntighten);
}
}
if(parentBB == lp->bb_bounds)
lp->bb_bounds = newBB;
else
newBB->child = parentBB->child;
if(parentBB != NULL)
parentBB->child = newBB;
lp->bb_level++;
if(lp->bb_level > lp->bb_maxlevel)
lp->bb_maxlevel = lp->bb_level;
if(!initbranches_BB(newBB))
newBB = pop_BB(newBB);
else if(MIP_count(lp) > 0) {
if( (lp->bb_level <= 1) && (lp->bb_varactive == NULL) &&
(!allocINT(lp, &lp->bb_varactive, lp->columns+1, TRUE) ||
!initcuts_BB(lp)) )
newBB = pop_BB(newBB);
if(varno > 0) {
lp->bb_varactive[varno-lp->rows]++;
}
}
}
return( newBB );
}
STATIC MYBOOL free_BB(BBrec **BB)
{
MYBOOL parentreturned = FALSE;
if((BB != NULL) && (*BB != NULL)) {
BBrec *parent = (*BB)->parent;
if((parent == NULL) || (*BB)->contentmode) {
FREE((*BB)->upbo);
FREE((*BB)->lowbo);
}
FREE((*BB)->varmanaged);
FREE(*BB);
parentreturned = (MYBOOL) (parent != NULL);
if(parentreturned)
*BB = parent;
}
return( parentreturned );
}
STATIC BBrec *pop_BB(BBrec *BB)
{
int k;
BBrec *parentBB;
lprec *lp = BB->lp;
if(BB == NULL)
return( BB );
parentBB = BB->parent;
if(BB == lp->bb_bounds) {
lp->bb_bounds = parentBB;
if(parentBB != NULL)
parentBB->child = NULL;
}
else {
if(parentBB != NULL)
parentBB->child = BB->child;
if(BB->child != NULL)
BB->child->parent = parentBB;
}
if(lp->bb_upperchange != NULL) {
restoreUndoLadder(lp->bb_upperchange, BB->upbo);
for(; BB->UBtrack > 0; BB->UBtrack--) {
decrementUndoLadder(lp->bb_upperchange);
restoreUndoLadder(lp->bb_upperchange, BB->upbo);
}
}
if(lp->bb_lowerchange != NULL) {
restoreUndoLadder(lp->bb_lowerchange, BB->lowbo);
for(; BB->LBtrack > 0; BB->LBtrack--) {
decrementUndoLadder(lp->bb_lowerchange);
restoreUndoLadder(lp->bb_lowerchange, BB->lowbo);
}
}
lp->bb_level--;
k = BB->varno - lp->rows;
if(lp->bb_level == 0) {
if(lp->bb_varactive != NULL) {
FREE(lp->bb_varactive);
freecuts_BB(lp);
}
if(lp->int_vars+lp->sc_vars > 0)
free_pseudocost(lp);
pop_basis(lp, FALSE);
lp->rootbounds = NULL;
}
else
lp->bb_varactive[k]--;
if(BB->isSOS && (BB->vartype != BB_INT))
SOS_unmark(lp->SOS, 0, k);
else if(BB->isGUB)
SOS_unmark(lp->GUB, 0, k);
if(BB->sc_canset)
lp->sc_lobound[k] *= -1;
#if 1
pop_basis(lp, FALSE);
#else#endif
free_BB(&BB);
return( parentBB );
}
STATIC REAL probe_BB(BBrec *BB)
{
int i, ii;
REAL coefOF, sum = 0;
lprec *lp = BB->lp;
if(lp->solutioncount == 0)
return( lp->infinite );
for(i = 1; i <= lp->columns; i++) {
if(!is_int(lp, i))
continue;
ii = lp->rows + i;
coefOF = lp->obj[i];
if(coefOF < 0) {
if(is_infinite(lp, BB->lowbo[ii]))
return( lp->infinite );
sum += coefOF * (lp->solution[ii]-BB->lowbo[ii]);
}
else {
if(is_infinite(lp, BB->upbo[ii]))
return( lp->infinite );
sum += coefOF * (BB->upbo[ii] - lp->solution[ii]);
}
}
return( sum );
}
STATIC REAL presolve_BB(BBrec *BB)
{
return( 0 );
}
STATIC MYBOOL initbranches_BB(BBrec *BB)
{
REAL new_bound, temp;
int k;
lprec *lp = BB->lp;
BB->nodestatus = NOTRUN;
BB->noderesult = lp->infinite;
push_basis(lp, NULL, NULL, NULL);
if(BB->vartype == BB_REAL)
BB->nodesleft = 1;
else {
BB->nodesleft = 2;
k = BB->varno - lp->rows;
BB->lastsolution = lp->solution[BB->varno];
BB->isSOS = (MYBOOL) ((BB->vartype == BB_SOS) || SOS_is_member(lp->SOS, 0, k));
#ifdef Paranoia
if((BB->vartype == BB_SOS) && !SOS_is_member(lp->SOS, 0, k))
report(lp, SEVERE, "initbranches_BB: Inconsistent identification of SOS variable %s (%d)\n",
get_col_name(lp, k), k);
#endif
BB->isGUB = (MYBOOL) ((BB->vartype == BB_INT) && SOS_can_activate(lp->GUB, 0, k));
if(BB->isGUB) {
BB->varmanaged = SOS_get_candidates(lp->GUB, -1, k, TRUE, BB->upbo, BB->lowbo);
BB->nodesleft++;
}
if(BB->vartype == BB_SOS) {
if(!SOS_can_activate(lp->SOS, 0, k)) {
BB->nodesleft--;
BB->isfloor = TRUE;
}
else
BB->isfloor = (MYBOOL) (BB->lastsolution == 0);
}
else if(lp->bb_usebranch != NULL)
BB->isfloor = (MYBOOL) lp->bb_usebranch(lp, lp->bb_branchhandle, k);
else if(get_var_branch(lp, k) == BRANCH_AUTOMATIC) {
new_bound = modf(BB->lastsolution/get_pseudorange(lp->bb_PseudoCost, k, BB->vartype), &temp);
if(isnan(new_bound))
new_bound = 0;
else if(new_bound < 0)
new_bound += 1.0;
BB->isfloor = (MYBOOL) (new_bound <= 0.5);
if(is_bb_mode(lp, NODE_GREEDYMODE)) {
if(is_bb_mode(lp, NODE_PSEUDOCOSTMODE))
BB->sc_bound = get_pseudonodecost(lp->bb_PseudoCost, k, BB->vartype, BB->lastsolution);
else
BB->sc_bound = mat_getitem(lp->matA, 0, k);
new_bound -= 0.5;
BB->sc_bound *= new_bound;
BB->isfloor = (MYBOOL) (BB->sc_bound > 0);
}
else if(is_bb_mode(lp, NODE_PSEUDOCOSTMODE)) {
BB->isfloor = (MYBOOL) (get_pseudobranchcost(lp->bb_PseudoCost, k, TRUE) >
get_pseudobranchcost(lp->bb_PseudoCost, k, FALSE));
if(is_maxim(lp))
BB->isfloor = !BB->isfloor;
}
if(is_bb_mode(lp, NODE_BRANCHREVERSEMODE))
BB->isfloor = !BB->isfloor;
}
else
BB->isfloor = (MYBOOL) (get_var_branch(lp, k) == BRANCH_FLOOR);
new_bound = fabs(lp->sc_lobound[k]);
BB->sc_bound = new_bound;
BB->sc_canset = (MYBOOL) (new_bound != 0);
new_bound = unscaled_value(lp, new_bound, BB->varno);
if(is_int(lp, k) && ((new_bound > 0) &&
(BB->lastsolution > floor(new_bound)))) {
if(BB->lastsolution < ceil(new_bound))
BB->lastsolution += 1;
modifyUndoLadder(lp->bb_lowerchange, BB->varno, BB->lowbo,
scaled_floor(lp, BB->varno, BB->lastsolution, 1));
}
}
return( fillbranches_BB(BB) );
}
STATIC MYBOOL fillbranches_BB(BBrec *BB)
{
int K, k;
REAL ult_upbo, ult_lowbo;
REAL new_bound, SC_bound, intmargin = BB->lp->epsprimal;
lprec *lp = BB->lp;
MYBOOL OKstatus = FALSE;
if(lp->bb_break || userabort(lp, MSG_MILPSTRATEGY))
return( OKstatus );
K = BB->varno;
if(K > 0) {
k = BB->varno - lp->rows;
ult_upbo = lp->orig_upbo[K];
ult_lowbo = lp->orig_lowbo[K];
SC_bound = unscaled_value(lp, BB->sc_bound, K);
BB->UPbound = lp->infinite;
if((SC_bound > 0) && (fabs(BB->lastsolution) < SC_bound-intmargin)) {
new_bound = 0;
}
else if(BB->vartype == BB_INT) {
#if 1
if(((ult_lowbo >= 0) &&
((floor(BB->lastsolution) <
unscaled_value(lp, MAX(ult_lowbo, fabs(lp->sc_lobound[k])), K)-intmargin))) ||
((ult_upbo <= 0) &&
((floor(BB->lastsolution) >
unscaled_value(lp, MIN(ult_upbo, -fabs(lp->sc_lobound[k])), K)-intmargin)))) {
#else#endif
BB->nodesleft--;
goto SetLB;
}
new_bound = scaled_floor(lp, K, BB->lastsolution, 1);
}
else if(BB->isSOS) {
new_bound = ult_lowbo;
if(is_int(lp, k))
new_bound = scaled_ceil(lp, K, unscaled_value(lp, new_bound, K), -1);
}
else
new_bound = BB->sc_bound;
if(new_bound < BB->lowbo[K])
new_bound = BB->lowbo[K] - my_avoidtiny(new_bound-BB->lowbo[K], intmargin);
if(new_bound < BB->lowbo[K]) {
#ifdef Paranoia
debug_print(lp,
"fillbranches_BB: New upper bound value %g conflicts with old lower bound %g\n",
new_bound, BB->lowbo[K]);
#endif
BB->nodesleft--;
goto SetLB;
}
#ifdef Paranoia
else if(!check_if_less(lp, new_bound, BB->upbo[K], K)) {
BB->nodesleft--;
goto SetLB;
}
#endif
else {
if(fabs(new_bound - BB->lowbo[K]) < intmargin*SCALEDINTFIXRANGE)
new_bound = BB->lowbo[K];
}
BB->UPbound = new_bound;
SetLB:
BB->LObound = -lp->infinite;
if((SC_bound > 0) && (fabs(BB->lastsolution) < SC_bound)) {
if(is_int(lp, k))
new_bound = scaled_ceil(lp, K, SC_bound, 1);
else
new_bound = BB->sc_bound;
}
else if((BB->vartype == BB_INT)) {
if(((ceil(BB->lastsolution) == BB->lastsolution)) ||
(ceil(BB->lastsolution) >
unscaled_value(lp, ult_upbo, K)+intmargin) ||
(BB->isSOS && (BB->lastsolution == 0))) {
BB->nodesleft--;
goto Finish;
}
new_bound = scaled_ceil(lp, K, BB->lastsolution, 1);
}
else if(BB->isSOS) {
if(SOS_is_member_of_type(lp->SOS, k, SOS3))
new_bound = scaled_floor(lp, K, 1, 1);
else {
new_bound = ult_lowbo;
if(is_int(lp, k))
new_bound = scaled_floor(lp, K, unscaled_value(lp, new_bound, K), 1);
if((lp->SOS->maxorder > 2) && (BB->lastsolution == 0) &&
SOS_is_member_of_type(lp->SOS, k, SOSn)) {
BB->isSOS = AUTOMATIC;
}
}
}
else
new_bound = BB->sc_bound;
if(new_bound > BB->upbo[K])
new_bound = BB->upbo[K] + my_avoidtiny(new_bound-BB->upbo[K], intmargin);
if(new_bound > BB->upbo[K]) {
#ifdef Paranoia
debug_print(lp,
"fillbranches_BB: New lower bound value %g conflicts with old upper bound %g\n",
new_bound, BB->upbo[K]);
#endif
BB->nodesleft--;
goto Finish;
}
#ifdef Paranoia
else if(!check_if_less(lp, BB->lowbo[K], new_bound, K)) {
BB->nodesleft--;
goto Finish;
}
#endif
else {
if(fabs(BB->upbo[K]-new_bound) < intmargin*SCALEDINTFIXRANGE)
new_bound = BB->upbo[K];
}
BB->LObound = new_bound;
Finish:
if(BB->nodesleft > 0) {
if(countsUndoLadder(lp->bb_upperchange) > 0) {
incrementUndoLadder(lp->bb_upperchange);
BB->UBtrack++;
}
if(countsUndoLadder(lp->bb_lowerchange) > 0) {
incrementUndoLadder(lp->bb_lowerchange);
BB->LBtrack++;
}
if((BB->vartype != BB_SOS) && (fabs(BB->LObound-BB->UPbound) < intmargin)) {
BB->nodesleft--;
if(fabs(BB->lowbo[K]-BB->LObound) < intmargin)
BB->isfloor = FALSE;
else if(fabs(BB->upbo[K]-BB->UPbound) < intmargin)
BB->isfloor = TRUE;
else
report(BB->lp, IMPORTANT, "fillbranches_BB: Inconsistent equal-valued bounds for %s\n",
get_col_name(BB->lp, k));
}
if((BB->nodesleft == 1) &&
((BB->isfloor && (BB->UPbound >= lp->infinite)) ||
(!BB->isfloor && (BB->LObound <= -lp->infinite))))
BB->isfloor = !BB->isfloor;
BB->isfloor = !BB->isfloor;
while(!OKstatus && lp->spx_status != TIMEOUT && !lp->bb_break && (BB->nodesleft > 0))
OKstatus = nextbranch_BB( BB );
}
if(BB->sc_canset)
lp->sc_lobound[k] *= -1;
}
else {
BB->nodesleft--;
OKstatus = TRUE;
}
return( OKstatus );
}
STATIC MYBOOL nextbranch_BB(BBrec *BB)
{
int k;
lprec *lp = BB->lp;
MYBOOL OKstatus = FALSE;
if(BB->nodessolved > 0) {
restoreUndoLadder(lp->bb_upperchange, BB->upbo);
restoreUndoLadder(lp->bb_lowerchange, BB->lowbo);
}
if(lp->bb_break || userabort(lp, MSG_MILPSTRATEGY)) {
if((lp->bb_level == 1) && (lp->bb_break == AUTOMATIC)) {
lp->bb_break = FALSE;
OKstatus = TRUE;
}
return( OKstatus );
}
if(BB->nodesleft > 0) {
k = BB->varno - lp->rows;
BB->isfloor = !BB->isfloor;
BB->nodesleft--;
if(BB->isSOS && (BB->vartype != BB_INT)) {
if((BB->nodessolved > 0) || ((BB->nodessolved == 0) && (BB->nodesleft == 0))) {
if(BB->isfloor) {
if((BB->nodesleft == 0) && (lp->orig_lowbo[BB->varno] != 0))
return( OKstatus );
}
SOS_unmark(lp->SOS, 0, k);
}
if(BB->isfloor) {
SOS_set_marked(lp->SOS, 0, k, (MYBOOL) (BB->UPbound != 0));
if(BB->isSOS == AUTOMATIC) {
}
}
else {
SOS_set_marked(lp->SOS, 0, k, TRUE);
if(SOS_fix_unmarked(lp->SOS, 0, k, BB->upbo, 0, TRUE,
NULL, lp->bb_upperchange) < 0)
return( OKstatus );
}
}
else if(BB->isGUB) {
if(BB->nodessolved > 0)
SOS_unmark(lp->GUB, 0, k);
if((BB->nodesleft == 0) && !BB->isfloor)
BB->isfloor = !BB->isfloor;
SOS_set_marked(lp->GUB, 0, k, (MYBOOL) !BB->isfloor);
if(BB->isfloor) {
if(SOS_fix_list(lp->GUB, 0, k, BB->upbo,
BB->varmanaged, (MYBOOL) (BB->nodesleft > 0), lp->bb_upperchange) < 0)
return( OKstatus );
}
else {
if(SOS_fix_unmarked(lp->GUB, 0, k, BB->upbo, 0, TRUE,
NULL, lp->bb_upperchange) < 0)
return( OKstatus );
}
}
OKstatus = TRUE;
}
if(OKstatus) {
lp->bb_totalnodes++;
BB->nodestatus = NOTRUN;
BB->noderesult = lp->infinite;
}
return( OKstatus );
}
STATIC MYBOOL initcuts_BB(lprec *lp)
{
return( TRUE );
}
STATIC int updatecuts_BB(lprec *lp)
{
return( 0 );
}
STATIC MYBOOL freecuts_BB(lprec *lp)
{
if(lp->bb_cuttype != NULL)
FREE(lp->bb_cuttype);
return( TRUE );
}
STATIC int solve_LP(lprec *lp, BBrec *BB)
{
int tilted, restored, status;
REAL testOF, *upbo = BB->upbo, *lowbo = BB->lowbo;
BBrec *perturbed = NULL;
if(lp->bb_break)
return(PROCBREAK);
#ifdef Paranoia
debug_print(lp, "solve_LP: Starting solve for iter %.0f, B&B node level %d.\n",
(double) lp->total_iter, lp->bb_level);
if(lp->bb_trace &&
!validate_bounds(lp, upbo, lowbo))
report(lp, SEVERE, "solve_LP: Inconsistent bounds at iter %.0f, B&B node level %d.\n",
(double) lp->total_iter, lp->bb_level);
#endif
impose_bounds(lp, upbo, lowbo);
if(BB->nodessolved > 1)
restore_basis(lp);
status = RUNNING;
tilted = 0;
restored = 0;
while(status == RUNNING) {
status = spx_run(lp, (MYBOOL) (tilted+restored > 0));
lp->bb_status = status;
lp->spx_perturbed = FALSE;
if(tilted < 0)
break;
else if((status == OPTIMAL) && (tilted > 0)) {
if(lp->spx_trace)
report(lp, DETAILED, "solve_LP: Restoring relaxed bounds at level %d.\n",
tilted);
free_BB(&perturbed);
if((perturbed == NULL) || (perturbed == BB)) {
perturbed = NULL;
impose_bounds(lp, upbo, lowbo);
}
else
impose_bounds(lp, perturbed->upbo, perturbed->lowbo);
set_action(&lp->spx_action, ACTION_REBASE | ACTION_RECOMPUTE);
BB->UBzerobased = FALSE;
if(lp->bb_totalnodes == 0)
lp->real_solution = lp->infinite;
status = RUNNING;
tilted--;
restored++;
lp->spx_perturbed = TRUE;
}
else if(((lp->bb_level <= 1) || is_anti_degen(lp, ANTIDEGEN_DURINGBB)) &&
(((status == LOSTFEAS) && is_anti_degen(lp, ANTIDEGEN_LOSTFEAS)) ||
((status == INFEASIBLE) && is_anti_degen(lp, ANTIDEGEN_INFEASIBLE)) ||
((status == NUMFAILURE) && is_anti_degen(lp, ANTIDEGEN_NUMFAILURE)) ||
((status == DEGENERATE) && is_anti_degen(lp, ANTIDEGEN_STALLING)))) {
if((tilted <= DEF_MAXRELAX) &&
!((tilted == 0) && (restored > DEF_MAXRELAX))) {
if(tilted == 0)
perturbed = BB;
perturbed = create_BB(lp, perturbed, TRUE);
#if 1
perturb_bounds(lp, perturbed, TRUE, TRUE, TRUE);
#else#endif
impose_bounds(lp, perturbed->upbo, perturbed->lowbo);
set_action(&lp->spx_action, ACTION_REBASE | ACTION_RECOMPUTE);
BB->UBzerobased = FALSE;
status = RUNNING;
tilted++;
lp->perturb_count++;
lp->spx_perturbed = TRUE;
if(lp->spx_trace)
report(lp, DETAILED, "solve_LP: Starting bound relaxation #%d ('%s')\n",
tilted, get_statustext(lp, status));
}
else {
if(lp->spx_trace)
report(lp, DETAILED, "solve_LP: Relaxation limit exceeded in resolving infeasibility\n");
while((perturbed != NULL) && (perturbed != BB))
free_BB(&perturbed);
perturbed = NULL;
}
}
}
if(status != OPTIMAL) {
if(lp->bb_level <= 1)
lp->bb_parentOF = lp->infinite;
if((status == USERABORT) || (status == TIMEOUT)) {
if((lp->solutioncount == 0) &&
(MIP_count(lp) == 0) &&
((lp->simplex_mode & (SIMPLEX_Phase2_PRIMAL | SIMPLEX_Phase2_DUAL)) > 0)) {
lp->solutioncount++;
construct_solution(lp, NULL);
transfer_solution(lp, TRUE);
}
report(lp, NORMAL, "\nlp_solve optimization was stopped %s.\n",
((status == USERABORT) ? "by the user" : "due to time-out"));
}
else if(BB->varno == 0)
report(lp, NORMAL, "The model %s\n",
(status == UNBOUNDED) ? "is UNBOUNDED" :
((status == INFEASIBLE) ? "is INFEASIBLE" : "FAILED"));
else {
#ifdef Paranoia
if((status != FATHOMED) && (status != INFEASIBLE))
report(lp, SEVERE, "spx_solve: Invalid return code %d during B&B\n", status);
#endif
if(status == FATHOMED)
lp->spx_status = INFEASIBLE;
}
}
else {
construct_solution(lp, NULL);
if((lp->bb_level <= 1) && (restored > 0))
report(lp, DETAILED, "%s numerics encountered; validate accuracy\n",
(restored == 1) ? "Difficult" : "Severe");
if(lp->spx_status != OPTIMAL)
status = lp->spx_status;
else if((lp->bb_totalnodes == 0) && (MIP_count(lp) > 0)) {
if(lp->lag_status != RUNNING) {
report(lp, NORMAL, "\nRelaxed solution " RESULTVALUEMASK " after %10.0f iter is B&B base.\n",
lp->solution[0], (double) lp->total_iter);
report(lp, NORMAL, " \n");
}
if((lp->usermessage != NULL) && (lp->msgmask & MSG_LPOPTIMAL)) {
REAL *best_solution = lp->best_solution;
lp->best_solution = lp->solution;
lp->usermessage(lp, lp->msghandle, MSG_LPOPTIMAL);
lp->best_solution = best_solution;
}
set_var_priority(lp);
}
testOF = my_chsign(is_maxim(lp), my_reldiff(lp->solution[0], lp->real_solution));
if(testOF < -lp->epsprimal) {
report(lp, DETAILED, "solve_LP: A MIP subproblem returned a value better than the base.\n");
status = INFEASIBLE;
lp->spx_status = status;
set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT | ACTION_RECOMPUTE);
}
else if(testOF < 0)
lp->solution[0] = lp->real_solution;
}
return( status );
}
STATIC BBrec *findself_BB(BBrec *BB)
{
int varno = BB->varno, vartype = BB->vartype;
BB = BB->parent;
while((BB != NULL) && (BB->vartype != vartype) && (BB->varno != varno))
BB = BB->parent;
return( BB );
}
STATIC int rcfbound_BB(BBrec *BB, int varno, MYBOOL isINT, REAL *newbound, MYBOOL *isfeasible)
{
int i = FR;
lprec *lp = BB->lp;
REAL deltaRC, rangeLU, deltaOF, lowbo, upbo;
if(lp->is_basic[varno])
return( i );
lowbo = BB->lowbo[varno];
upbo = BB->upbo[varno];
rangeLU = upbo - lowbo;
if(rangeLU > lp->epsprimal) {
#if 1
deltaOF = lp->rhs[0] - lp->bb_workOF;
#elif 0
deltaOF = my_chsign(is_maxim(lp), lp->real_solution) - lp->bb_workOF;
#else#endif
deltaRC = my_chsign(!lp->is_lower[varno], lp->drow[varno]);
if(deltaRC < lp->epspivot)
return( i );
deltaRC = deltaOF / deltaRC;
#ifdef Paranoia
if(deltaRC <= 0)
report(lp, SEVERE, "rcfbound_BB: A negative bound fixing level was encountered after node %.0f\n",
(double) lp->bb_totalnodes);
#endif
if(deltaRC < rangeLU + lp->epsint) {
if(lp->is_lower[varno]) {
if(isINT)
deltaRC = scaled_floor(lp, varno, unscaled_value(lp, deltaRC, varno)+lp->epsprimal, 1);
upbo = lowbo + deltaRC;
deltaRC = upbo;
i = LE;
}
else {
if(isINT)
deltaRC = scaled_ceil(lp, varno, unscaled_value(lp, deltaRC, varno)+lp->epsprimal, 1);
lowbo = upbo - deltaRC;
deltaRC = lowbo;
i = GE;
}
if((isfeasible != NULL) && (upbo - lowbo < -lp->epsprimal))
*isfeasible = FALSE;
else if(fabs(upbo - lowbo) < lp->epsprimal)
i = -i;
if(newbound != NULL) {
my_roundzero(deltaRC, lp->epsprimal);
*newbound = deltaRC;
}
}
}
return( i );
}
STATIC MYBOOL findnode_BB(BBrec *BB, int *varno, int *vartype, int *varcus)
{
int countsossc, countnint, k, reasonmsg = MSG_NONE;
REAL varsol;
MYBOOL is_better = FALSE, is_equal = FALSE, is_feasible = TRUE;
lprec *lp = BB->lp;
*varno = 0;
*vartype = BB_REAL;
*varcus = 0;
countnint = 0;
BB->nodestatus = lp->spx_status;
BB->noderesult = lp->solution[0];
if((lp->bb_limitlevel != 1) && (MIP_count(lp) > 0)) {
countsossc = lp->sos_vars + lp->sc_vars;
if((lp->bb_limitlevel > 0) && (lp->bb_level > lp->bb_limitlevel+countsossc))
return( FALSE );
else if((lp->bb_limitlevel < 0) &&
(lp->bb_level > 2*(lp->int_vars+countsossc)*abs(lp->bb_limitlevel))) {
if(lp->bb_limitlevel == DEF_BB_LIMITLEVEL)
report(lp, IMPORTANT, "findnode_BB: Default B&B limit reached at %d; optionally change strategy or limit.\n\n",
lp->bb_level);
return( FALSE );
}
if(BB->varno == 0) {
varsol = lp->infinite;
if((lp->int_vars+lp->sc_vars > 0) && (lp->bb_PseudoCost == NULL))
lp->bb_PseudoCost = init_pseudocost(lp, get_bb_rule(lp));
}
else {
varsol = lp->solution[BB->varno];
if( ((lp->int_vars > 0) && (BB->vartype == BB_INT)) ||
((lp->sc_vars > 0) && (BB->vartype == BB_SC) && !is_int(lp, BB->varno-lp->rows)) )
update_pseudocost(lp->bb_PseudoCost, BB->varno-lp->rows, BB->vartype, BB->isfloor, varsol);
}
if((lp->bb_totalnodes > 0) && !bb_better(lp, OF_RELAXED, OF_TEST_WE)) {
if(lp->bb_trace)
report(lp, IMPORTANT, "findnode_BB: Simplex failure due to loss of numeric accuracy\n");
lp->spx_status = NUMFAILURE;
return( FALSE );
}
if(((lp->solutioncount == 0) && !bb_better(lp, OF_HEURISTIC, OF_TEST_BE)) ||
((lp->solutioncount > 0) &&
(!bb_better(lp, OF_INCUMBENT | OF_DELTA, OF_TEST_BE | OF_TEST_RELGAP) ||
!bb_better(lp, OF_INCUMBENT | OF_DELTA, OF_TEST_BE)))) {
return( FALSE );
}
if(lp->sc_vars > 0) {
*varno = find_sc_bbvar(lp, &countnint);
if(*varno > 0)
*vartype = BB_SC;
}
if((SOS_count(lp) > 0) && (*varno == 0)) {
*varno = find_sos_bbvar(lp, &countnint, FALSE);
if(*varno < 0)
*varno = 0;
else if(*varno > 0)
*vartype = BB_SOS;
}
if((lp->int_vars > 0) && (*varno == 0)) {
*varno = find_int_bbvar(lp, &countnint, BB, &is_feasible);
if(*varno > 0) {
*vartype = BB_INT;
if((countnint == 1) && !is_feasible) {
BB->lastrcf = 0;
return( FALSE );
}
}
}
#if 1
k = *varno-lp->rows;
if((*varno > 0) && (lp->bb_limitlevel != 0) && (lp->bb_varactive[k] >= abs(lp->bb_limitlevel) )) {
return( FALSE );
}
#endif
if(*varno == 0) {
is_better = (MYBOOL) (lp->solutioncount == 0) || bb_better(lp, OF_INCUMBENT | OF_DELTA, OF_TEST_BT);
#if 1
is_better &= bb_better(lp, OF_INCUMBENT | OF_DELTA, OF_TEST_BT | OF_TEST_RELGAP);
#else#endif
is_equal = !is_better;
if(is_equal) {
if((lp->solutionlimit <= 0) || (lp->solutioncount < lp->solutionlimit)) {
lp->solutioncount++;
SETMIN(lp->bb_solutionlevel, lp->bb_level);
reasonmsg = MSG_MILPEQUAL;
}
}
else if(is_better) {
if(lp->bb_varactive != NULL) {
lp->bb_varactive[0]++;
if((lp->bb_varactive[0] == 1) &&
is_bb_mode(lp, NODE_DEPTHFIRSTMODE) && is_bb_mode(lp, NODE_DYNAMICMODE))
lp->bb_rule &= !NODE_DEPTHFIRSTMODE;
}
if(lp->bb_trace ||
((lp->verbose >= NORMAL) && (lp->print_sol == FALSE) && (lp->lag_status != RUNNING))) {
report(lp, IMPORTANT,
"%s solution " RESULTVALUEMASK " after %10.0f iter, %9.0f nodes (gap %.1f%%)\n",
(lp->bb_improvements == 0) ? "Feasible" : "Improved",
lp->solution[0], (double) lp->total_iter, (double) lp->bb_totalnodes,
100.0*fabs(my_reldiff(lp->solution[0], lp->bb_limitOF)));
}
if(MIP_count(lp) > 0) {
if(lp->bb_improvements == 0)
reasonmsg = MSG_MILPFEASIBLE;
else
reasonmsg = MSG_MILPBETTER;
}
lp->bb_status = FEASFOUND;
lp->bb_solutionlevel = lp->bb_level;
lp->solutioncount = 1;
lp->bb_improvements++;
lp->bb_workOF = lp->rhs[0];
if(lp->bb_breakfirst ||
(!is_infinite(lp, lp->bb_breakOF) && bb_better(lp, OF_USERBREAK, OF_TEST_BE)))
lp->bb_break = TRUE;
}
}
}
else {
is_better = TRUE;
lp->solutioncount = 1;
}
if(is_better || is_equal) {
#ifdef ParanoiaMIP
if((lp->bb_level > 0) &&
(check_solution(lp, lp->columns, lp->solution,
lp->orig_upbo, lp->orig_lowbo, lp->epssolution) != OPTIMAL)) {
lp->solutioncount = 0;
lp->spx_status = NUMFAILURE;
lp->bb_status = lp->spx_status;
lp->bb_break = TRUE;
return( FALSE );
}
#endif
transfer_solution(lp, (MYBOOL) ((lp->do_presolve & PRESOLVE_LASTMASKMODE) != PRESOLVE_NONE));
if((MIP_count(lp) > 0) && (lp->bb_totalnodes > 0)) {
if ((!construct_duals(lp)) ||
(is_presolve(lp, PRESOLVE_SENSDUALS) &&
(!construct_sensitivity_duals(lp) || !construct_sensitivity_obj(lp))
)
) {
}
}
if((reasonmsg != MSG_NONE) && (lp->msgmask & reasonmsg) && (lp->usermessage != NULL))
lp->usermessage(lp, lp->msghandle, reasonmsg);
if(lp->print_sol != FALSE) {
print_objective(lp);
print_solution(lp, 1);
}
}
*varcus = countnint;
if(MIP_count(lp) > 0) {
if((countnint == 0) && (lp->solutioncount == 1) && (lp->solutionlimit == 1) &&
(bb_better(lp, OF_DUALLIMIT, OF_TEST_BE) || bb_better(lp, OF_USERBREAK, OF_TEST_BE | OF_TEST_RELGAP))) {
lp->bb_break = (MYBOOL) (countnint == 0);
return( FALSE );
}
else if(lp->bb_level > 0) {
#ifdef MIPboundWithOF
if((lp->constraintOF > 0) && (countnint == 0))
set_rh(lp, lp->constraintOF, lp->solution[0] + my_chsign(!is_maxim(lp), lp->bb_deltaOF));
#endif
if(lp->spx_trace)
report(lp, DETAILED, "B&B level %5d OPT %16s value " RESULTVALUEMASK "\n",
lp->bb_level, (*varno) ? " " : "INT", lp->solution[0]);
}
return( (MYBOOL) (*varno > 0));
}
else
return( FALSE );
}
STATIC int solve_BB(BBrec *BB)
{
int K, status;
lprec *lp = BB->lp;
status = PROCFAIL;
K = BB->varno;
if(K > 0) {
updatecuts_BB(lp);
if(BB->isfloor)
modifyUndoLadder(lp->bb_upperchange, K, BB->upbo, BB->UPbound);
else
modifyUndoLadder(lp->bb_lowerchange, K, BB->lowbo, BB->LObound);
BB->nodessolved++;
}
status = solve_LP(lp, BB);
#if 1
if((status == OPTIMAL) && (BB->vartype == BB_SOS) && !SOS_is_feasible(lp->SOS, 0, lp->solution))
status = INFEASIBLE;
#endif
return( status );
}
STATIC MYBOOL strongbranch_BB(lprec *lp, BBrec *BB, int varno, int vartype, int varcus)
{
MYBOOL success = FALSE;
int i;
BBrec *strongBB;
lp->is_strongbranch = TRUE;
push_basis(lp, lp->var_basic, lp->is_basic, lp->is_lower);
strongBB = push_BB(lp, BB, lp->rows+varno, vartype, varcus);
if(strongBB == BB)
return( success );
do {
lp->bb_strongbranches++;
if(solve_BB(strongBB) == OPTIMAL) {
success |= 1 << strongBB->isfloor;
strongBB->lastvarcus = 0;
for(i = 1; i <= lp->columns; i++) {
if(is_int(lp, i) && !solution_is_int(lp, lp->rows+i, FALSE))
strongBB->lastvarcus++;
}
update_pseudocost(lp->bb_PseudoCost, varno, strongBB->vartype, strongBB->isfloor,
lp->solution[strongBB->varno]);
}
}
while(nextbranch_BB(strongBB));
strongBB = pop_BB(strongBB);
if(strongBB != BB)
report(lp, SEVERE, "strongbranch_BB: Invalid bound settings restored for variable %d\n",
varno);
pop_basis(lp, TRUE);
set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT | ACTION_RECOMPUTE);
lp->is_strongbranch = FALSE;
return( success );
}
STATIC MYBOOL pre_BB(lprec *lp)
{
return( TRUE );
}
STATIC MYBOOL post_BB(lprec *lp)
{
return( TRUE );
}
STATIC int run_BB(lprec *lp)
{
BBrec *currentBB;
int varno, vartype, varcus, prevsolutions;
int status = NOTRUN;
pre_BB(lp);
prevsolutions = lp->solutioncount;
#ifdef UseMilpSlacksRCF
varno = lp->sum;
#else
varno = lp->columns;
#endif
lp->bb_upperchange = createUndoLadder(lp, varno, 2*MIP_count(lp));
lp->bb_lowerchange = createUndoLadder(lp, varno, 2*MIP_count(lp));
lp->rootbounds = currentBB = push_BB(lp, NULL, 0, BB_REAL, 0);
while(lp->bb_level > 0) {
status = solve_BB(currentBB);
#if 0#endif
if((status == OPTIMAL) && findnode_BB(currentBB, &varno, &vartype, &varcus))
currentBB = push_BB(lp, currentBB, varno, vartype, varcus);
else while((lp->bb_level > 0) && !nextbranch_BB(currentBB))
currentBB = pop_BB(currentBB);
}
freeUndoLadder(&(lp->bb_upperchange));
freeUndoLadder(&(lp->bb_lowerchange));
if(lp->solutioncount > prevsolutions) {
if((status == PROCBREAK) || (status == USERABORT) || (status == TIMEOUT) || userabort(lp, -1))
status = SUBOPTIMAL;
else
status = OPTIMAL;
if(lp->bb_totalnodes > 0)
lp->spx_status = OPTIMAL;
}
post_BB(lp);
return( status );
}