prop_obbt.c

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the program and library             */
/*         SCIP --- Solving Constraint Integer Programs                      */
/*                                                                           */
/*  Copyright (c) 2002-2025 Zuse Institute Berlin (ZIB)                      */
/*                                                                           */
/*  Licensed under the Apache License, Version 2.0 (the "License");          */
/*  you may not use this file except in compliance with the License.         */
/*  You may obtain a copy of the License at                                  */
/*                                                                           */
/*      http://www.apache.org/licenses/LICENSE-2.0                           */
/*                                                                           */
/*  Unless required by applicable law or agreed to in writing, software      */
/*  distributed under the License is distributed on an "AS IS" BASIS,        */
/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
/*  See the License for the specific language governing permissions and      */
/*  limitations under the License.                                           */
/*                                                                           */
/*  You should have received a copy of the Apache-2.0 license                */
/*  along with SCIP; see the file LICENSE. If not visit scipopt.org.         */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
/**@file    prop_obbt.c
 * @ingroup DEFPLUGINS_PROP
 * @brief   optimization-based bound tightening propagator
 * @author  Stefan Weltge
 * @author  Benjamin Mueller
 */
 
/**@todo if bound tightenings of other propagators are the reason for lpsolstat != SCIP_LPSOLSTAT_OPTIMAL, resolve LP */
/**@todo only run more than once in root node if primal bound improved or many cuts were added to the LP */
/**@todo filter bounds of a variable already if SCIPisLbBetter()/SCIPisUbBetter() would return FALSE */
/**@todo improve warmstarting of LP solving */
/**@todo include bound value (finite/infinite) into getScore() function */
/**@todo use unbounded ray in filtering */
/**@todo do we want to run if the LP is unbounded, maybe for infinite variable bounds? */
/**@todo add first filter round in direction of objective function */
/**@todo implement conflict resolving callback by calling public method of genvbounds propagator, since the reason are
 *       exactly the variable bounds with nonnegative reduced costs stored in the right-hand side of the generated
 *       generalized variable bound (however, this only makes sense if we run locally)
 */
 
/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
 
#include <assert.h>
#include <string.h>
 
#include "scip/cons_indicator.h"
#include "scip/cons_linear.h"
#include "scip/cons_nonlinear.h"
#include "scip/nlhdlr_bilinear.h"
#include "scip/prop_genvbounds.h"
#include "scip/prop_obbt.h"
#include "scip/pub_cons.h"
#include "scip/pub_lp.h"
#include "scip/pub_message.h"
#include "scip/pub_misc.h"
#include "scip/pub_misc_sort.h"
#include "scip/pub_nlp.h"
#include "scip/pub_prop.h"
#include "scip/pub_tree.h"
#include "scip/pub_var.h"
#include "scip/scip_cons.h"
#include "scip/scip_copy.h"
#include "scip/scip_cut.h"
#include "scip/scip_general.h"
#include "scip/scip_lp.h"
#include "scip/scip_mem.h"
#include "scip/scip_message.h"
#include "scip/scip_nlp.h"
#include "scip/scip_numerics.h"
#include "scip/scip_param.h"
#include "scip/scip_prob.h"
#include "scip/scip_probing.h"
#include "scip/scip_prop.h"
#include "scip/scip_randnumgen.h"
#include "scip/scip_solvingstats.h"
#include "scip/scip_tree.h"
#include "scip/scip_var.h"
 
#define PROP_NAME                       "obbt"
#define PROP_DESC                       "optimization-based bound tightening propagator"
#define PROP_TIMING                     SCIP_PROPTIMING_AFTERLPLOOP
#define PROP_PRIORITY               -1000000 /**< propagator priority */
#define PROP_FREQ                          0 /**< propagator frequency */
#define PROP_DELAY                      TRUE /**< should propagation method be delayed, if other propagators
                                              *   found reductions? */
 
#define DEFAULT_CREATE_GENVBOUNDS       TRUE /**< should obbt try to provide genvbounds if possible? */
#define DEFAULT_FILTERING_NORM          TRUE /**< should coefficients in filtering be normalized w.r.t. the
                                              *   domains sizes? */
#define DEFAULT_APPLY_FILTERROUNDS     FALSE /**< try to filter bounds in so-called filter rounds by solving
                                              *   auxiliary LPs? */
#define DEFAULT_APPLY_TRIVIALFITLERING  TRUE /**< should obbt try to use the LP solution to filter some bounds? */
#define DEFAULT_GENVBDSDURINGFILTER     TRUE /**< try to genrate genvbounds during trivial and aggressive filtering? */
#define DEFAULT_DUALFEASTOL             1e-9 /**< feasibility tolerance for reduced costs used in obbt; this value
                                              *   is used if SCIP's dual feastol is greater */
#define DEFAULT_CONDITIONLIMIT          -1.0 /**< maximum condition limit used in LP solver (-1.0: no limit) */
#define DEFAULT_BOUNDSTREPS            0.001 /**< minimal relative improve for strengthening bounds */
#define DEFAULT_FILTERING_MIN              2 /**< minimal number of filtered bounds to apply another filter
                                              *   round */
#define DEFAULT_ITLIMITFACTOR           10.0 /**< multiple of root node LP iterations used as total LP iteration
                                              *   limit for obbt (<= 0: no limit ) */
#define DEFAULT_MINITLIMIT             5000L /**< minimum LP iteration limit */
#define DEFAULT_ONLYNONCONVEXVARS       TRUE /**< only apply obbt on non-convex variables */
#define DEFAULT_INDICATORS             FALSE /**< apply obbt on variables of indicator constraints? (independent of convexity) */
#define DEFAULT_INDICATORTHRESHOLD       1e6 /**< variables of indicator constraints with smaller upper bound are not considered
                                              *   and upper bound is tightened only if new bound is smaller */
#define DEFAULT_TIGHTINTBOUNDSPROBING   TRUE /**< should bounds of integral variables be tightened during
                                              *   the probing mode? */
#define DEFAULT_TIGHTCONTBOUNDSPROBING FALSE /**< should bounds of continuous variables be tightened during
                                              *   the probing mode? */
#define DEFAULT_ORDERINGALGO               1 /**< which type of ordering algorithm should we use?
                                              *   (0: no, 1: greedy, 2: greedy reverse) */
#define OBBT_SCOREBASE                     5 /**< base that is used to calculate a bounds score value */
#define GENVBOUND_PROP_NAME    "genvbounds"
 
#define DEFAULT_SEPARATESOL            FALSE /**< should the obbt LP solution be separated? note that that by
                                              *   separating solution OBBT will apply all bound tightenings
                                              *   immediatly */
#define DEFAULT_SEPAMINITER                0 /**< minimum number of iteration spend to separate an obbt LP solution */
#define DEFAULT_SEPAMAXITER               10 /**< maximum number of iteration spend to separate an obbt LP solution */
#define DEFAULT_GENVBDSDURINGSEPA       TRUE /**< try to create genvbounds during separation process? */
#define DEFAULT_PROPAGATEFREQ              0 /**< trigger a propagation round after that many bound tightenings
                                              *   (0: no propagation) */
#define DEFAULT_CREATE_BILININEQS       TRUE /**< solve auxiliary LPs in order to find valid inequalities for bilinear terms? */
#define DEFAULT_CREATE_LINCONS         FALSE /**< create linear constraints from inequalities for bilinear terms? */
#define DEFAULT_ITLIMITFAC_BILININEQS    3.0 /**< multiple of OBBT LP limit used as total LP iteration limit for solving bilinear inequality LPs (< 0 for no limit) */
#define DEFAULT_MINNONCONVEXITY         1e-1 /**< minimum nonconvexity for choosing a bilinear term */
#define DEFAULT_RANDSEED                 149 /**< initial random seed */
 
/*
 * Data structures
 */
 
/** bound data */
struct Bound
{
   SCIP_VAR*             var;                /**< variable */
   SCIP_Real             newval;             /**< stores a probably tighter value for this bound */
   SCIP_BOUNDTYPE        boundtype;          /**< type of bound */
   unsigned int          score;              /**< score value that is used to group bounds */
   unsigned int          filtered:1;         /**< thrown out during pre-filtering step */
   unsigned int          found:1;            /**< stores whether a probably tighter value for this bound was found */
   unsigned int          done:1;             /**< has this bound been processed already? */
   unsigned int          nonconvex:1;        /**< is this bound affecting a nonconvex term? */
   unsigned int          indicator:1;        /**< is this bound affecting an indicator constraint? */
   int                   index;              /**< unique index */
};
typedef struct Bound BOUND;
 
/* all possible corners of a rectangular domain */
enum Corner
{
   LEFTBOTTOM  = 1,
   RIGHTBOTTOM = 2,
   RIGHTTOP    = 4,
   LEFTTOP     = 8,
   FILTERED    = 15
};
typedef enum Corner CORNER;
 
/** bilinear bound data */
struct BilinBound
{
   SCIP_EXPR*            expr;               /**< product expression */
   int                   filtered;           /**< corners that could be thrown out during pre-filtering step */
   unsigned int          done:1;             /**< has this bilinear term been processed already? */
   SCIP_Real             score;              /**< score value that is used to group bilinear term bounds */
};
typedef struct BilinBound BILINBOUND;
 
/** propagator data */
struct SCIP_PropData
{
   BOUND**               bounds;             /**< array of interesting bounds */
   BILINBOUND**          bilinbounds;        /**< array of interesting bilinear bounds */
   SCIP_ROW*             cutoffrow;          /**< pointer to current objective cutoff row */
   SCIP_PROP*            genvboundprop;      /**< pointer to genvbound propagator */
   SCIP_RANDNUMGEN*      randnumgen;         /**< random number generator */
   SCIP_Longint          lastnode;           /**< number of last node where obbt was performed */
   SCIP_Longint          npropagatedomreds;  /**< number of domain reductions found during propagation */
   SCIP_Longint          nprobingiterations; /**< number of LP iterations during the probing mode */
   SCIP_Longint          nfilterlpiters;     /**< number of LP iterations spend for filtering */
   SCIP_Longint          minitlimit;         /**< minimum LP iteration limit */
   SCIP_Longint          itlimitbilin;       /**< total LP iterations limit for solving bilinear inequality LPs */
   SCIP_Longint          itusedbilin;        /**< total LP iterations used for solving bilinear inequality LPs */
   SCIP_Real             dualfeastol;        /**< feasibility tolerance for reduced costs used in obbt; this value is
                                              *   used if SCIP's dual feastol is greater */
   SCIP_Real             conditionlimit;     /**< maximum condition limit used in LP solver (-1.0: no limit) */
   SCIP_Real             boundstreps;        /**< minimal relative improve for strengthening bounds */
   SCIP_Real             itlimitfactor;      /**< LP iteration limit for obbt will be this factor times total LP
                                              *   iterations in root node */
   SCIP_Real             itlimitfactorbilin; /**< multiple of OBBT LP limit used as total LP iteration limit for solving bilinear inequality LPs (< 0 for no limit) */
   SCIP_Real             minnonconvexity;    /**< lower bound on minimum absolute value of nonconvex eigenvalues for a bilinear term */
   SCIP_Real             indicatorthreshold; /**< threshold whether upper bounds of vars of indicator conss are considered or tightened */
   SCIP_Bool             applyfilterrounds;  /**< apply filter rounds? */
   SCIP_Bool             applytrivialfilter; /**< should obbt try to use the LP solution to filter some bounds? */
   SCIP_Bool             genvbdsduringfilter;/**< should we try to generate genvbounds during trivial and aggressive
                                              *   filtering? */
   SCIP_Bool             genvbdsduringsepa;  /**< try to create genvbounds during separation process? */
   SCIP_Bool             creategenvbounds;   /**< should obbt try to provide genvbounds if possible? */
   SCIP_Bool             normalize;          /**< should coefficients in filtering be normalized w.r.t. the domains
                                              *   sizes? */
   SCIP_Bool             onlynonconvexvars;  /**< only apply obbt on non-convex variables */
   SCIP_Bool             indicators;         /**< apply obbt on variables of indicator constraints? (independent of convexity) */
   SCIP_Bool             tightintboundsprobing; /**< should bounds of integral variables be tightened during
                                              *   the probing mode? */
   SCIP_Bool             tightcontboundsprobing;/**< should bounds of continuous variables be tightened during
                                              *   the probing mode? */
   SCIP_Bool             separatesol;        /**< should the obbt LP solution be separated? note that that by
                                              *   separating solution OBBT will apply all bound tightenings
                                              *   immediatly */
   SCIP_Bool             createbilinineqs;   /**< solve auxiliary LPs in order to find valid inequalities for bilinear terms? */
   SCIP_Bool             createlincons;      /**< create linear constraints from inequalities for bilinear terms? */
   int                   orderingalgo;       /**< which type of ordering algorithm should we use?
                                              *   (0: no, 1: greedy, 2: greedy reverse) */
   int                   nbounds;            /**< length of interesting bounds array */
   int                   nbilinbounds;       /**< length of interesting bilinear bounds array */
   int                   bilinboundssize;    /**< size of bilinear bounds array */
   int                   boundssize;         /**< size of bounds array */
   int                   nminfilter;         /**< minimal number of filtered bounds to apply another filter round */
   int                   nfiltered;          /**< number of filtered bounds by solving auxiliary variables */
   int                   ntrivialfiltered;   /**< number of filtered bounds because the LP value was equal to the bound */
   int                   nsolvedbounds;      /**< number of solved bounds during the loop in applyObbt() */
   int                   ngenvboundsprobing; /**< number of non-trivial genvbounds generated and added during obbt */
   int                   ngenvboundsaggrfil; /**< number of non-trivial genvbounds found during aggressive filtering */
   int                   ngenvboundstrivfil; /**< number of non-trivial genvbounds found during trivial filtering */
   int                   lastidx;            /**< index to store the last undone and unfiltered bound */
   int                   lastbilinidx;       /**< index to store the last undone and unfiltered bilinear bound */
   int                   sepaminiter;        /**< minimum number of iteration spend to separate an obbt LP solution */
   int                   sepamaxiter;        /**< maximum number of iteration spend to separate an obbt LP solution */
   int                   propagatefreq;      /**< trigger a propagation round after that many bound tightenings
                                              *   (0: no propagation) */
   int                   propagatecounter;   /**< number of bound tightenings since the last propagation round */
};
 
 
/*
 * Local methods
 */
 
/** solves the LP and handles errors */
static
SCIP_RETCODE solveLP(
   SCIP*                 scip,               /**< SCIP data structure */
   int                   itlimit,            /**< maximal number of LP iterations to perform, or -1 for no limit */
   SCIP_Bool*            error,              /**< pointer to store whether an unresolved LP error occurred */
   SCIP_Bool*            optimal             /**< was the LP solved to optimalilty? */
   )
{
   SCIP_LPSOLSTAT lpsolstat;
   SCIP_RETCODE retcode;
 
   assert(scip != NULL);
   assert(itlimit == -1 || itlimit >= 0);
   assert(error != NULL);
   assert(optimal != NULL);
 
   *optimal = FALSE;
   *error = FALSE;
 
   retcode = SCIPsolveProbingLP(scip, itlimit, error, NULL);
 
   lpsolstat = SCIPgetLPSolstat(scip);
 
   /* an error should not kill the overall solving process */
   if( retcode != SCIP_OKAY )
   {
      SCIPwarningMessage(scip, "   error while solving LP in obbt propagator; LP solve terminated with code <%d>\n", retcode);
      SCIPwarningMessage(scip, "   this does not affect the remaining solution procedure --> continue\n");
 
      *error = TRUE;
 
      return SCIP_OKAY;
   }
 
   if( lpsolstat == SCIP_LPSOLSTAT_OPTIMAL )
   {
      assert(!*error);
      *optimal = TRUE;
   }
#ifdef SCIP_DEBUG
   else
   {
      switch( lpsolstat )
      {
      case SCIP_LPSOLSTAT_ITERLIMIT:
         SCIPdebugMsg(scip, "   reached lp iteration limit\n");
         break;
      case SCIP_LPSOLSTAT_TIMELIMIT:
         SCIPdebugMsg(scip, "   reached time limit while solving lp\n");
         break;
      case SCIP_LPSOLSTAT_UNBOUNDEDRAY:
         SCIPdebugMsg(scip, "   lp was unbounded\n");
         break;
      case SCIP_LPSOLSTAT_NOTSOLVED:
         SCIPdebugMsg(scip, "   lp was not solved\n");
         break;
      case SCIP_LPSOLSTAT_ERROR:
         SCIPdebugMsg(scip, "   an error occured during solving lp\n");
         break;
      case SCIP_LPSOLSTAT_INFEASIBLE:
      case SCIP_LPSOLSTAT_OBJLIMIT:
      case SCIP_LPSOLSTAT_OPTIMAL: /* should not appear because it is handled earlier */
      default:
         SCIPdebugMsg(scip, "   received an unexpected solstat during solving lp: %d\n", lpsolstat);
      }
   }
#endif
 
   return SCIP_OKAY;
}
 
/** adds the objective cutoff to the LP; must be in probing mode */
static
SCIP_RETCODE addObjCutoff(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROPDATA*        propdata            /**< data of the obbt propagator */
   )
{
   SCIP_ROW* row;
   SCIP_VAR** vars;
   char rowname[SCIP_MAXSTRLEN];
 
   int nvars;
   int i;
 
   assert(scip != NULL);
   assert(SCIPinProbing(scip));
   assert(propdata != NULL);
   assert(propdata->cutoffrow == NULL);
 
   if( SCIPisInfinity(scip, SCIPgetCutoffbound(scip)) )
   {
      SCIPdebugMsg(scip, "no objective cutoff since there is no cutoff bound\n");
      return SCIP_OKAY;
   }
 
   SCIPdebugMsg(scip, "create objective cutoff and add it to the LP\n");
 
   /* get variables data */
   SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
 
   /* create objective cutoff row; set local flag to FALSE since primal cutoff is globally valid */
   (void) SCIPsnprintf(rowname, SCIP_MAXSTRLEN, "obbt_objcutoff");
   SCIP_CALL( SCIPcreateEmptyRowUnspec(scip, &row, rowname, -SCIPinfinity(scip), SCIPgetCutoffbound(scip), FALSE, FALSE, FALSE) );
   SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
 
   for( i = 0; i < nvars; i++ )
   {
      SCIP_CALL( SCIPaddVarToRow(scip, row, vars[i], SCIPvarGetObj(vars[i])) );
   }
   SCIP_CALL( SCIPflushRowExtensions(scip, row) );
 
   /* add row to the LP */
   SCIP_CALL( SCIPaddRowProbing(scip, row) );
 
   propdata->cutoffrow = row;
   assert(SCIProwIsInLP(propdata->cutoffrow));
 
   return SCIP_OKAY;
}
 
/** determines, whether a variable is already locally fixed */
static
SCIP_Bool varIsFixedLocal(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_VAR*             var                 /**< variable to check */
   )
{
   return SCIPisFeasEQ(scip, SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var));
}
 
/** sets objective to minimize or maximize a single variable */
static
SCIP_RETCODE setObjProbing(
   SCIP*                 scip,
   SCIP_PROPDATA*        propdata,
   BOUND*                bound,
   SCIP_Real             coef
   )
{
#ifdef SCIP_DEBUG
   SCIP_VAR** vars;
   int nvars;
   int counter;
   int i;
#endif
 
   assert( scip != NULL );
   assert( propdata != NULL );
   assert( bound != NULL );
 
   /* set the objective for bound->var */
   if( bound->boundtype == SCIP_BOUNDTYPE_LOWER )
   {
      SCIP_CALL( SCIPchgVarObjProbing(scip, bound->var, coef) );
   }
   else
   {
      SCIP_CALL( SCIPchgVarObjProbing(scip, bound->var, -coef) );
   }
 
#ifdef SCIP_DEBUG
   vars = SCIPgetVars(scip);
   nvars = SCIPgetNVars(scip);
   counter = 0;
 
   for( i = 0; i < nvars; ++i )
   {
      if( SCIPgetVarObjProbing(scip, vars[i]) != 0.0 )
         ++counter;
   }
 
   assert((counter == 0 && coef == 0.0) || (counter == 1 && coef != 0.0));
#endif
 
   return SCIP_OKAY;
}
 
/** determines whether variable should be included in the right-hand side of the generalized variable bound */
static
SCIP_Bool includeVarGenVBound(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_VAR*             var                 /**< variable to check */
   )
{
   SCIP_Real redcost;
 
   assert(scip != NULL);
   assert(var != NULL);
 
   if( SCIPvarGetStatus(var) != SCIP_VARSTATUS_COLUMN )
      return FALSE;
 
   redcost = SCIPgetVarRedcost(scip, var);
   assert(redcost != SCIP_INVALID); /*lint !e777 */
 
   if( redcost == SCIP_INVALID ) /*lint !e777 */
      return FALSE;
 
   if( redcost < SCIPdualfeastol(scip) && redcost > -SCIPdualfeastol(scip) )
      return FALSE;
 
   return TRUE;
}
 
/** returns number of LP iterations left (-1: no limit ) */
static
int getIterationsLeft(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_Longint          nolditerations,     /**< iterations count at the beginning of the corresponding function */
   SCIP_Longint          itlimit             /**< LP iteration limit (-1: no limit) */
   )
{
   SCIP_Longint itsleft;
 
   assert(scip != NULL);
   assert(nolditerations >= 0);
   assert(itlimit == -1 || itlimit >= 0);
 
   if( itlimit == -1 )
   {
      SCIPdebugMsg(scip, "iterations left: unlimited\n");
      return -1;
   }
   else
   {
      itsleft = itlimit - ( SCIPgetNLPIterations(scip) - nolditerations );
      itsleft = MAX(itsleft, 0);
      itsleft = MIN(itsleft, INT_MAX);
 
      SCIPdebugMsg(scip, "iterations left: %d\n", (int) itsleft);
      return (int) itsleft;
   }
}
 
/** returns the objective coefficient for a variable's bound that will be chosen during filtering */
static
SCIP_Real getFilterCoef(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROPDATA*        propdata,           /**< data of the obbt propagator */
   SCIP_VAR*             var,                /**< variable */
   SCIP_BOUNDTYPE        boundtype           /**< boundtype to be filtered? */
   )
{
   SCIP_Real lb;
   SCIP_Real ub;
 
   assert(scip != NULL);
   assert(propdata != NULL);
   assert(var != NULL);
 
   lb = SCIPvarGetLbLocal(var);
   ub = SCIPvarGetUbLocal(var);
 
   /* this function should not be called for fixed variables */
   assert(!varIsFixedLocal(scip, var));
 
   /* infinite bounds will not be reached */
   if( boundtype == SCIP_BOUNDTYPE_LOWER && SCIPisInfinity(scip, -lb) )
      return 0.0;
   if( boundtype == SCIP_BOUNDTYPE_UPPER && SCIPisInfinity(scip, ub) )
      return 0.0;
 
   if( propdata->normalize )
   {
      /* if the length of the domain is too large then the coefficient should be set to +/- 1.0 */
      if( boundtype == SCIP_BOUNDTYPE_LOWER && SCIPisInfinity(scip, ub) )
         return 1.0;
      if( boundtype == SCIP_BOUNDTYPE_UPPER && SCIPisInfinity(scip, -lb) )
         return -1.0;
 
      /* otherwise the coefficient is +/- 1.0 / ( ub - lb ) */
      return boundtype == SCIP_BOUNDTYPE_LOWER ? 1.0 / (ub - lb) : -1.0 / (ub - lb);
   }
   else
   {
      return boundtype == SCIP_BOUNDTYPE_LOWER ? 1.0 : -1.0;
   }
}
 
/** creates a genvbound if the dual LP solution provides such information
 *
 *  Consider the problem
 *
 *     min { +/- x_i : obj * x <= z, lb <= Ax <= ub, l <= x <= u },
 *
 *  where z is the current cutoff bound. Let (mu, nu, gamma, alpha, beta) >= 0 be the optimal solution of the dual of
 *  problem (P), where the variables correspond to the primal inequalities in the following way:
 *
 *           Ax >=  lb    <->   mu
 *          -Ax >= -ub    <->   nu
 *     -obj * x >=  -z    <->   gamma
 *            x >=   l    <->   alpha
 *           -x >=  -u    <->   beta
 *
 *  Fixing these multipliers, by weak duality, we obtain the inequality
 *
 *     +/- x_i >= lb*mu - ub*nu - z*gamma + l*alpha - u*beta
 *
 *  that holds for all primal feasible points x with objective value at least z. Setting
 *
 *     c = lb*mu - ub*nu, redcost_k = alpha_k - beta_k
 *
 *  we obtain the inequality
 *
 *     +/- x_i >= sum ( redcost_k * x_k ) + (-gamma) * cutoff_bound + c,
 *
 *  that holds for all primal feasible points with objective value at least cutoff_bound. Therefore, the latter
 *  inequality can be added as a generalized variable bound.
 */
static
SCIP_RETCODE createGenVBound(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROPDATA*        propdata,           /**< data of the obbt propagator */
   BOUND*                bound,              /**< bound of x_i */
   SCIP_Bool*            found               /**< pointer to store if we have found a non-trivial genvbound */
   )
{
   assert(scip != NULL);
   assert(bound != NULL);
   assert(propdata != NULL);
   assert(propdata->genvboundprop != NULL);
   assert(found != NULL);
 
   *found = FALSE;
 
   /* make sure we are in probing mode having an optimal LP solution */
   assert(SCIPinProbing(scip));
 
   assert(SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL);
 
   /* only genvbounds created in the root node are globally valid
    *
    * note: depth changes to one if we use the probing mode to solve the obbt LPs
    */
   assert(SCIPgetDepth(scip) == 0 || (SCIPinProbing(scip) && SCIPgetDepth(scip) == 1));
 
   SCIPdebugMsg(scip, "      try to create a genvbound for <%s>...\n", SCIPvarGetName(bound->var));
 
   /* a genvbound with a multiplier for x_i would not help us */
   if( SCIPisZero(scip, SCIPgetVarRedcost(scip, bound->var)) )
   {
      SCIP_VAR** vars;                          /* global variables array */
      SCIP_VAR** genvboundvars;                 /* genvbound variables array */
 
      SCIP_VAR* xi;                             /* variable x_i */
 
      SCIP_Real* genvboundcoefs;                /* genvbound coefficients array */
 
      SCIP_Real gamma_dual;                     /* dual multiplier of objective cutoff */
 
      int k;                                    /* variable for indexing global variables array */
      int ncoefs;                               /* number of nonzero coefficients in genvbound */
      int nvars;                                /* number of global variables */
 
      /* set x_i */
      xi = bound->var;
 
      /* get variable data */
      SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
 
      /* count nonzero coefficients in genvbound */
      ncoefs = 0;
      for( k = 0; k < nvars; k++ )
      {
         if( includeVarGenVBound(scip, vars[k]) )
         {
            assert(vars[k] != xi);
            ncoefs++;
         }
      }
 
      /* get dual multiplier for the objective cutoff (set to zero if there is no) */
      if( propdata->cutoffrow == NULL )
      {
         gamma_dual = 0.0;
      }
      else
      {
         assert(!SCIPisInfinity(scip, SCIPgetCutoffbound(scip)));
 
         /* note that the objective cutoff is of the form
          *    -inf <= obj * x <= cutoff_bound
          * but we want the positive dual multiplier!
          */
         gamma_dual = -SCIProwGetDualsol(propdata->cutoffrow);
 
         /* we need to treat gamma to be exactly 0 if it is below the dual feasibility tolerance, see #2914 */
         if( EPSZ(gamma_dual, SCIPdualfeastol(scip)) )
            gamma_dual = 0.0;
      }
 
      /* we need at least one nonzero coefficient or a nonzero dual multiplier for the objective cutoff */
      if( ncoefs > 0 || gamma_dual != 0.0 )
      {
         SCIP_Bool addgenvbound;                /* if everything is fine with the redcosts and the bounds, add the genvbound */
         SCIP_Real c;                           /* helper variable to calculate constant term in genvbound */
         int idx;                               /* variable for indexing genvbound's coefficients array */
 
         /* add the bound if the bool is still TRUE after the loop */
         addgenvbound = TRUE;
 
         /* there should be no coefficient for x_i */
         assert(SCIPisZero(scip, SCIPgetVarRedcost(scip, xi)));
 
         /* allocate memory for storing the genvbounds right-hand side variables and coefficients */
         SCIP_CALL( SCIPallocBufferArray(scip, &(genvboundvars), ncoefs) );
         SCIP_CALL( SCIPallocBufferArray(scip, &(genvboundcoefs), ncoefs) );
 
         /* set c = lb*mu - ub*nu - z*gamma + l*alpha - u*beta */
         c = SCIPgetLPObjval(scip);
 
         /* subtract ( - z * gamma ) from c */
         c += SCIPgetCutoffbound(scip) * gamma_dual;
 
         /* subtract ( l*alpha - u*beta ) from c and set the coefficients of the variables */
         idx = 0;
         for( k = 0; k < nvars; k++ )
         {
            SCIP_VAR* xk;
 
            xk = vars[k];
 
            if( includeVarGenVBound(scip, xk) )
            {
               SCIP_Real redcost;
 
               redcost = SCIPgetVarRedcost(scip, xk);
 
               assert(redcost != SCIP_INVALID); /*lint !e777 */
               assert(xk != xi);
 
               /* in this case dont add a genvbound */
               if( ( (redcost > SCIPdualfeastol(scip)) && SCIPisInfinity(scip, -SCIPvarGetLbLocal(xk)) ) ||
                  ( (redcost < -SCIPdualfeastol(scip)) && SCIPisInfinity(scip, SCIPvarGetUbLocal(xk)) ) )
               {
                  addgenvbound = FALSE;
                  break;
               }
 
               /* store coefficients */
               assert(idx < ncoefs);
               genvboundvars[idx] = xk;
               genvboundcoefs[idx] = redcost;
               idx++;
 
               /* if redcost > 0, then redcost = alpha_k, otherwise redcost = - beta_k */
               assert(redcost <= 0 || !SCIPisInfinity(scip, -SCIPvarGetLbLocal(xk)));
               assert(redcost >= 0 || !SCIPisInfinity(scip, SCIPvarGetUbLocal(xk)));
               c -= redcost > 0 ? redcost * SCIPvarGetLbLocal(xk) : redcost * SCIPvarGetUbLocal(xk);
            }
         }
 
         assert(!addgenvbound || idx == ncoefs);
 
         /* add genvbound */
         if( addgenvbound && !SCIPisInfinity(scip, -c) )
         {
#ifndef NDEBUG
            /* check whether the activity of the LVB in the optimal solution of the LP is equal to the LP objective value */
            SCIP_Real activity = c - gamma_dual * SCIPgetCutoffbound(scip);
 
            for( k = 0; k < ncoefs; ++k )
               activity += genvboundcoefs[k] * SCIPvarGetLPSol(genvboundvars[k]);
 
            SCIPdebugMsg(scip, "LVB activity = %g lpobj = %g\n", activity, SCIPgetLPObjval(scip));
            assert(EPSZ(SCIPrelDiff(activity, SCIPgetLPObjval(scip)), 18.0 * SCIPdualfeastol(scip)));
#endif
 
            SCIPdebugMsg(scip, "         adding genvbound\n");
            SCIP_CALL( SCIPgenVBoundAdd(scip, propdata->genvboundprop, genvboundvars, xi, genvboundcoefs, ncoefs,
                  gamma_dual < SCIPdualfeastol(scip) ? 0.0 : -gamma_dual, c, bound->boundtype) );
            *found = TRUE;
         }
 
         /* free arrays */
         SCIPfreeBufferArray(scip, &genvboundcoefs);
         SCIPfreeBufferArray(scip, &genvboundvars);
      }
      else
      {
         SCIPdebugMsg(scip, "         trivial genvbound, skipping\n");
      }
   }
   else
   {
      SCIPdebugMsg(scip, "         found multiplier for <%s>: %g, skipping\n",
         SCIPvarGetName(bound->var), SCIPgetVarRedcost(scip, bound->var));
   }
 
   return SCIP_OKAY;
}
 
/** exchange a bound which has been processed and updates the last undone and unfiltered bound index
 *  NOTE: this method has to be called after filtering or processing a bound
 */
static
void exchangeBounds(
   SCIP_PROPDATA*        propdata,           /**< propagator data */
   int                   i                   /**< bound that was filtered or processed */
   )
{
   assert(i >= 0 && i < propdata->nbounds);
   assert(propdata->lastidx >= 0 && propdata->lastidx < propdata->nbounds);
 
   /* exchange the bounds */
   if( propdata->lastidx != i )
   {
      BOUND* tmp;
 
      tmp = propdata->bounds[i];
      propdata->bounds[i] = propdata->bounds[propdata->lastidx];
      propdata->bounds[propdata->lastidx] = tmp;
   }
 
   propdata->lastidx -= 1;
}
 
/** helper function to return a corner of the domain of two variables */
static
void getCorner(
   SCIP_VAR*             x,                  /**< first variable */
   SCIP_VAR*             y,                  /**< second variable */
   CORNER                corner,             /**< corner */
   SCIP_Real*            px,                 /**< buffer to store point for x */
   SCIP_Real*            py                  /**< buffer to store point for y */
   )
{
   assert(x != NULL);
   assert(y != NULL);
   assert(px != NULL);
   assert(py != NULL);
 
   switch( corner )
   {
      case LEFTBOTTOM:
         *px = SCIPvarGetLbGlobal(x);
         *py = SCIPvarGetLbGlobal(y);
         break;
      case RIGHTBOTTOM:
         *px = SCIPvarGetUbGlobal(x);
         *py = SCIPvarGetLbGlobal(y);
         break;
      case LEFTTOP:
         *px = SCIPvarGetLbGlobal(x);
         *py = SCIPvarGetUbGlobal(y);
         break;
      case RIGHTTOP:
         *px = SCIPvarGetUbGlobal(x);
         *py = SCIPvarGetUbGlobal(y);
         break;
      case FILTERED:
         SCIPABORT();
   }
}
 
/** helper function to return the two end points of a diagonal */
static
void getCorners(
   SCIP_VAR*             x,                  /**< first variable */
   SCIP_VAR*             y,                  /**< second variable */
   CORNER                corner,             /**< corner */
   SCIP_Real*            xs,                 /**< buffer to store start point for x */
   SCIP_Real*            ys,                 /**< buffer to store start point for y */
   SCIP_Real*            xt,                 /**< buffer to store end point for x */
   SCIP_Real*            yt                  /**< buffer to store end point for y */
   )
{
   assert(x != NULL);
   assert(y != NULL);
   assert(xs != NULL);
   assert(ys != NULL);
   assert(xt != NULL);
   assert(yt != NULL);
 
   /* get end point */
   getCorner(x,y, corner, xt, yt);
 
   /* get start point */
   switch( corner )
   {
      case LEFTBOTTOM:
         getCorner(x,y, RIGHTTOP, xs, ys);
         break;
      case RIGHTBOTTOM:
         getCorner(x,y, LEFTTOP, xs, ys);
         break;
      case LEFTTOP:
         getCorner(x,y, RIGHTBOTTOM, xs, ys);
         break;
      case RIGHTTOP:
         getCorner(x,y, LEFTBOTTOM, xs, ys);
         break;
      case FILTERED:
         SCIPABORT();
   }
}
 
/** returns the first variable of a bilinear bound */
static
SCIP_VAR* bilinboundGetX(
   BILINBOUND*           bilinbound          /**< bilinear bound */
   )
{
   assert(bilinbound->expr != NULL);
   assert(SCIPexprGetNChildren(bilinbound->expr) == 2);
 
   return SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(bilinbound->expr)[0]);
}
 
/** returns the second variable of a bilinear bound */
static
SCIP_VAR* bilinboundGetY(
   BILINBOUND*           bilinbound          /**< bilinear bound */
   )
{
   assert(bilinbound->expr != NULL);
   assert(SCIPexprGetNChildren(bilinbound->expr) == 2);
 
   return SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(bilinbound->expr)[1]);
}
 
/** returns the negative locks of the expression in a bilinear bound */
static
int bilinboundGetLocksNeg(
   BILINBOUND*           bilinbound          /**< bilinear bound */
   )
{
   assert(bilinbound->expr != NULL);
 
   return SCIPgetExprNLocksNegNonlinear(bilinbound->expr);
}
 
/** returns the positive locks of the expression in a bilinear bound */
static
int bilinboundGetLocksPos(
   BILINBOUND*           bilinbound          /**< bilinear bound */
   )
{
   assert(bilinbound->expr != NULL);
 
   return SCIPgetExprNLocksPosNonlinear(bilinbound->expr);
}
 
/** computes the score of a bilinear term bound */
static
SCIP_Real bilinboundGetScore(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_RANDNUMGEN*      randnumgen,         /**< random number generator */
   BILINBOUND*           bilinbound          /**< bilinear bound */
   )
{
   SCIP_VAR* x = bilinboundGetX(bilinbound);
   SCIP_VAR* y = bilinboundGetY(bilinbound);
   SCIP_Real lbx = SCIPvarGetLbLocal(x);
   SCIP_Real ubx = SCIPvarGetUbLocal(x);
   SCIP_Real lby = SCIPvarGetLbLocal(y);
   SCIP_Real uby = SCIPvarGetUbLocal(y);
   SCIP_Real score;
 
   assert(scip != NULL);
   assert(randnumgen != NULL);
   assert(bilinbound != NULL);
 
   /* consider how often a bilinear term is present in the problem */
   score = bilinboundGetLocksNeg(bilinbound) + bilinboundGetLocksPos(bilinbound);
 
   /* penalize small variable domains; TODO tune the factor in the logarithm, maybe add a parameter for it */
   if( ubx - lbx < 0.5 )
      score += log(2.0*(ubx-lbx) + SCIPepsilon(scip));
   if( uby - lby < 0.5 )
      score += log(2.0*(uby-lby) + SCIPepsilon(scip));
 
   /* consider interiority of variables in the LP solution */
   if( SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL )
   {
      SCIP_Real solx = SCIPvarGetLPSol(x);
      SCIP_Real soly = SCIPvarGetLPSol(y);
      SCIP_Real interiorityx = MIN(solx-lbx, ubx-solx) / MAX(ubx-lbx, SCIPepsilon(scip)); /*lint !e666*/
      SCIP_Real interiorityy = MIN(soly-lby, uby-soly) / MAX(uby-lby, SCIPepsilon(scip)); /*lint !e666*/
 
      score += interiorityx + interiorityy;
   }
 
   /* randomize score */
   score *= 1.0 + SCIPrandomGetReal(randnumgen, -SCIPepsilon(scip), SCIPepsilon(scip));
 
   return score;
}
 
/** determines whether a variable of an indicator constraint is (still) interesting
 *
 * A variable is interesting if it is not only part of indicator constraints or if the upper bound is greater than the given threshold.
 */
static
SCIP_Bool indicatorVarIsInteresting(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_VAR*             var,                /**< variable to check */
   int                   nlcount,            /**< number of nonlinear constraints containing the variable
                                              *   or number of non-convex terms containing the variable
                                              *  (depends on propdata->onlynonconvexvars)  */
   int                   nindcount,          /**< number of indicator constraints containing the variable
                                              *   or 0 (depends on propdata->indicators) */
   SCIP_Real             threshold           /**< variables with smaller upper bound are not interesting */
   )
{
   /* if variable is only part of indicator constraints, consider current upper bound */
   if( nlcount == 0 && nindcount > 0 )
   {
      if( SCIPisLE(scip, SCIPvarGetUbLocal(var), threshold) )
         return FALSE;
   }
 
   return TRUE;
}
 
/** trying to filter some bounds using the existing LP solution */
static
SCIP_RETCODE filterExistingLP(
   SCIP*                 scip,               /**< original SCIP data structure */
   SCIP_PROPDATA*        propdata,           /**< data of the obbt propagator */
   int*                  nfiltered,          /**< how many bounds were filtered this round? */
   BOUND*                currbound           /**< bound for which OBBT LP was solved (Note: might be NULL) */
   )
{
   int i;
 
   assert(scip != NULL);
   assert(propdata != NULL);
   assert(nfiltered != NULL);
 
   *nfiltered = 0;
 
   /* only apply filtering if an LP solution is at hand */
   if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
   {
      SCIPdebugMsg(scip, "can't filter using existing lp solution since it was not solved to optimality\n");
      return SCIP_OKAY;
   }
 
   /* check if a bound is tight */
   for( i = propdata->nbounds - 1; i >= 0; --i )
   {
      BOUND* bound;                          /* shortcut for current bound */
 
      SCIP_Real solval;                      /* the variables value in the current solution */
      SCIP_Real boundval;                    /* current local bound for the variable */
 
      bound = propdata->bounds[i];
      if( bound->filtered || bound->done )
         continue;
 
      boundval = bound->boundtype == SCIP_BOUNDTYPE_UPPER ?
         SCIPvarGetUbLocal(bound->var) : SCIPvarGetLbLocal(bound->var);
      solval = SCIPvarGetLPSol(bound->var);
 
      /* bound is tight; since this holds for all fixed variables, those are filtered here automatically; if the lp solution
       * is infinity, then also the bound is tight */
      if( (bound->boundtype == SCIP_BOUNDTYPE_UPPER &&
               (SCIPisInfinity(scip, solval) || SCIPisFeasGE(scip, solval, boundval)))
            || (bound->boundtype == SCIP_BOUNDTYPE_LOWER &&
               (SCIPisInfinity(scip, -solval) || SCIPisFeasLE(scip, solval, boundval))) )
      {
         SCIP_BASESTAT basestat;
 
         /* mark bound as filtered */
         bound->filtered = TRUE;
         SCIPdebugMsg(scip, "trivial filtered var: %s boundval=%e solval=%e\n", SCIPvarGetName(bound->var), boundval, solval);
 
         /* get the basis status of the variable */
         basestat = SCIPcolGetBasisStatus(SCIPvarGetCol(bound->var));
 
         /* solve corresponding OBBT LP and try to generate a nontrivial genvbound */
         if( propdata->genvbdsduringfilter && currbound != NULL && basestat == SCIP_BASESTAT_BASIC )
         {
#ifndef NDEBUG
            int j;
#endif
            SCIP_Bool optimal;
            SCIP_Bool error;
 
            /* set objective coefficient of the bound */
            SCIP_CALL( SCIPchgVarObjProbing(scip, currbound->var, 0.0) );
            SCIP_CALL( setObjProbing(scip, propdata, bound, 1.0) );
 
#ifndef NDEBUG
            for( j = 0; j < SCIPgetNVars(scip); ++j )
            {
               SCIP_VAR* var;
 
               var = SCIPgetVars(scip)[j];
               assert(var != NULL);
               assert(SCIPisZero(scip, SCIPgetVarObjProbing(scip, var)) || var == bound->var);
            }
#endif
 
            /* solve the OBBT LP */
            propdata->nprobingiterations -= SCIPgetNLPIterations(scip);
            SCIP_CALL( solveLP(scip, -1, &error, &optimal) );
            propdata->nprobingiterations += SCIPgetNLPIterations(scip);
            assert(propdata->nprobingiterations >= 0);
 
            /* try to generate a genvbound if we have solved the OBBT LP */
            if( optimal && propdata->genvboundprop != NULL
                  && (SCIPgetDepth(scip) == 0 || (SCIPinProbing(scip) && SCIPgetDepth(scip) == 1)) )
            {
               SCIP_Bool found;
 
               assert(!error);
               SCIP_CALL( createGenVBound(scip, propdata, bound, &found) );
 
               if( found )
               {
                  propdata->ngenvboundstrivfil += 1;
                  SCIPdebugMsg(scip, "found genvbound during trivial filtering\n");
               }
            }
 
            /* restore objective function */
            SCIP_CALL( setObjProbing(scip, propdata, bound, 0.0) );
            SCIP_CALL( setObjProbing(scip, propdata, currbound, 1.0) );
         }
 
         /* exchange bound i with propdata->bounds[propdata->lastidx] */
         if( propdata->lastidx >= 0 )
            exchangeBounds(propdata, i);
 
         /* increase number of filtered variables */
         (*nfiltered)++;
      }
   }
 
   /* try to filter bilinear bounds */
   for( i = propdata->lastbilinidx; i < propdata->nbilinbounds; ++i )
   {
      CORNER corners[4] = {LEFTTOP, LEFTBOTTOM, RIGHTTOP, RIGHTBOTTOM};
      BILINBOUND* bilinbound = propdata->bilinbounds[i];
      SCIP_Real solx;
      SCIP_Real soly;
      SCIPdebug(int oldfiltered;)
      int j;
 
      /* skip processed and filtered bounds */
      if( bilinbound->done || bilinbound->filtered == FILTERED ) /*lint !e641*/
         continue;
 
      SCIPdebug(oldfiltered = bilinbound->filtered;)
      solx = SCIPvarGetLPSol(bilinboundGetX(bilinbound));
      soly = SCIPvarGetLPSol(bilinboundGetY(bilinbound));
 
      /* check cases of unbounded solution values */
      if( SCIPisInfinity(scip, solx) )
         bilinbound->filtered = bilinbound->filtered | RIGHTTOP | RIGHTBOTTOM; /*lint !e641*/
      else if( SCIPisInfinity(scip, -solx) )
         bilinbound->filtered = bilinbound->filtered | LEFTTOP | LEFTBOTTOM; /*lint !e641*/
 
      if( SCIPisInfinity(scip, soly) )
         bilinbound->filtered = bilinbound->filtered | RIGHTTOP | LEFTTOP; /*lint !e641*/
      else if( SCIPisInfinity(scip, -soly) )
         bilinbound->filtered = bilinbound->filtered | RIGHTBOTTOM | LEFTBOTTOM; /*lint !e641*/
 
      /* check all corners */
      for( j = 0; j < 4; ++j )
      {
         SCIP_Real xt = SCIP_INVALID;
         SCIP_Real yt = SCIP_INVALID;
 
         getCorner(bilinboundGetX(bilinbound), bilinboundGetY(bilinbound), corners[j], &xt, &yt);
 
         if( (SCIPisInfinity(scip, REALABS(solx)) || SCIPisFeasEQ(scip, xt, solx))
            && (SCIPisInfinity(scip, REALABS(soly)) || SCIPisFeasEQ(scip, yt, soly)) )
            bilinbound->filtered = bilinbound->filtered | corners[j]; /*lint !e641*/
      }
 
#ifdef SCIP_DEBUG
      if( oldfiltered != bilinbound->filtered )
      {
         SCIP_VAR* x = bilinboundGetX(bilinbound);
         SCIP_VAR* y = bilinboundGetY(bilinbound);
         SCIPdebugMessage("filtered corners %d for (%s,%s) = (%g,%g) in [%g,%g]x[%g,%g]\n",
            bilinbound->filtered - oldfiltered, SCIPvarGetName(x), SCIPvarGetName(y), solx, soly,
            SCIPvarGetLbGlobal(x), SCIPvarGetUbGlobal(x), SCIPvarGetLbGlobal(y), SCIPvarGetUbGlobal(y));
      }
#endif
   }
 
   return SCIP_OKAY;
}
 
/** enforces one round of filtering */
static
SCIP_RETCODE filterRound(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROPDATA*        propdata,           /**< data of the obbt propagator */
   int                   itlimit,            /**< LP iteration limit (-1: no limit) */
   int*                  nfiltered,          /**< how many bounds were filtered this round */
   SCIP_Real*            objcoefs,           /**< array to store the nontrivial objective coefficients */
   int*                  objcoefsinds,       /**< array to store bound indices for which their corresponding variables
                                              *   has a nontrivial objective coefficient */
   int                   nobjcoefs           /**< number of nontrivial objective coefficients */
   )
{
   SCIP_VAR** vars;                          /* array of the problems variables */
   SCIP_Bool error;
   SCIP_Bool optimal;
 
   int nvars;                                /* number of the problems variables */
   int i;
 
   assert(scip != NULL);
   assert(SCIPinProbing(scip));
   assert(propdata != NULL);
   assert(itlimit == -1 || itlimit >= 0);
   assert(nfiltered != NULL);
   assert(objcoefs != NULL);
   assert(objcoefsinds != NULL);
   assert(nobjcoefs >= 0);
 
   *nfiltered = 0;
 
   /* get variable data */
   SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
 
   /* solve LP */
   propdata->nfilterlpiters -= (int) SCIPgetNLPIterations(scip);
   SCIP_CALL( solveLP(scip, itlimit, &error, &optimal) );
   propdata->nfilterlpiters += (int) SCIPgetNLPIterations(scip);
   assert(propdata->nfilterlpiters >= 0);
 
   if( !optimal )
   {
      SCIPdebugMsg(scip, "skipping filter round since the LP was not solved to optimality\n");
      return SCIP_OKAY;
   }
 
   assert(!error);
 
   /* check if a bound is tight */
   for( i = 0; i < propdata->nbounds; i++ )
   {
      BOUND* bound;                          /* shortcut for current bound */
 
      SCIP_Real solval;                      /* the variables value in the current solution */
      SCIP_Real boundval;                    /* current local bound for the variable */
 
      bound = propdata->bounds[i];
 
      /* if bound is filtered it was handled already before */
      if( bound->filtered )
         continue;
 
      boundval = bound->boundtype == SCIP_BOUNDTYPE_UPPER ?
         SCIPvarGetUbLocal(bound->var) : SCIPvarGetLbLocal(bound->var);
      solval = SCIPvarGetLPSol(bound->var);
 
      /* bound is tight */
      if( (bound->boundtype == SCIP_BOUNDTYPE_UPPER && SCIPisFeasGE(scip, solval, boundval))
         || (bound->boundtype == SCIP_BOUNDTYPE_LOWER && SCIPisFeasLE(scip, solval, boundval)) )
      {
         SCIP_Real objcoef;
         SCIP_BASESTAT basestat;
 
         /* mark bound as filtered */
         bound->filtered = TRUE;
 
         /* get the basis status of the variable */
         basestat = SCIPcolGetBasisStatus(SCIPvarGetCol(bound->var));
 
         /* increase number of filtered variables */
         (*nfiltered)++;
 
         /* solve corresponding OBBT LP and try to generate a nontrivial genvbound */
         if( propdata->genvbdsduringfilter && basestat == SCIP_BASESTAT_BASIC )
         {
            int j;
 
            /* set all objective coefficients to zero */
            for( j = 0; j < nobjcoefs; ++j )
            {
               BOUND* filterbound;
 
               filterbound = propdata->bounds[ objcoefsinds[j] ];
               assert(filterbound != NULL);
 
               SCIP_CALL( SCIPchgVarObjProbing(scip, filterbound->var, 0.0) );
            }
 
#ifndef NDEBUG
            for( j = 0; j < nvars; ++j )
               assert(SCIPisZero(scip, SCIPgetVarObjProbing(scip, vars[j])));
#endif
 
            /* set objective coefficient of the bound */
            SCIP_CALL( setObjProbing(scip, propdata, bound, 1.0) );
 
            /* solve the OBBT LP */
            propdata->nfilterlpiters -= (int) SCIPgetNLPIterations(scip);
            SCIP_CALL( solveLP(scip, -1, &error, &optimal) );
            propdata->nfilterlpiters += (int) SCIPgetNLPIterations(scip);
            assert(propdata->nfilterlpiters >= 0);
 
            /* try to generate a genvbound if we have solved the OBBT LP */
            if( optimal && propdata->genvboundprop != NULL
                  && (SCIPgetDepth(scip) == 0 || (SCIPinProbing(scip) && SCIPgetDepth(scip) == 1)) )
            {
               SCIP_Bool found;
 
               assert(!error);
               SCIP_CALL( createGenVBound(scip, propdata, bound, &found) );
 
               if( found )
               {
                  propdata->ngenvboundsaggrfil += 1;
                  SCIPdebugMsg(scip, "found genvbound during aggressive filtering\n");
               }
            }
 
            /* restore objective function */
            for( j = 0; j < nobjcoefs; ++j )
            {
               BOUND* filterbound;
 
               filterbound = propdata->bounds[ objcoefsinds[j] ];
               assert(filterbound != NULL);
 
               /* NOTE: only restore coefficients of nonfiltered bounds */
               if( !filterbound->filtered )
               {
                  assert(!SCIPisZero(scip, objcoefs[j]));
                  SCIP_CALL( SCIPchgVarObjProbing(scip, propdata->bounds[ objcoefsinds[j] ]->var, objcoefs[j]) );
               }
            }
         }
 
         /* get the corresponding variable's objective coefficient */
         objcoef = SCIPgetVarObjProbing(scip, bound->var);
 
         /* change objective coefficient if it was set up for this bound */
          if( (bound->boundtype == SCIP_BOUNDTYPE_UPPER && SCIPisNegative(scip, objcoef))
             || (bound->boundtype == SCIP_BOUNDTYPE_LOWER && SCIPisPositive(scip, objcoef)) )
          {
             SCIP_CALL( SCIPchgVarObjProbing(scip, bound->var, 0.0) );
          }
      }
   }
 
   return SCIP_OKAY;
}
 
/** filter some bounds that are not improvable by solving auxiliary LPs */
static
SCIP_RETCODE filterBounds(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROPDATA*        propdata,           /**< data of the obbt propagator */
   SCIP_Longint          itlimit             /**< LP iteration limit (-1: no limit) */
   )
{
   SCIP_VAR** vars;
   SCIP_Longint nolditerations;
   SCIP_Real* objcoefs;               /* array to store the nontrivial objective coefficients */
   int* objcoefsinds;                 /* array to store bound indices for which the corresponding variable
                                       * has a nontrivial objective coefficient */
   int nobjcoefs;                     /* number of nontrivial objective coefficients */
   int nleftiterations;
   int i;
   int nfiltered;
   int ntotalfiltered;
   int nvars;
 
   assert(scip != NULL);
   assert(SCIPinProbing(scip));
   assert(propdata != NULL);
   assert(itlimit == -1 || itlimit >= 0);
 
   ntotalfiltered = 0;
   nolditerations = SCIPgetNLPIterations(scip);
   nleftiterations = getIterationsLeft(scip, nolditerations, itlimit);
 
   /* get variable data */
   SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
 
   SCIPdebugMsg(scip, "start filter rounds\n");
 
   SCIP_CALL( SCIPallocBufferArray(scip, &objcoefs, propdata->nbounds) );
   SCIP_CALL( SCIPallocBufferArray(scip, &objcoefsinds, propdata->nbounds) );
   nobjcoefs = 0;
 
   /*
    * 1.) Filter bounds of variables that are only part of indicator constraints if they are not interesting any more
    */
   for( i = 0; i < propdata->nbounds; i++ )
   {
      if( !propdata->bounds[i]->filtered && !propdata->bounds[i]->done && propdata->bounds[i]->indicator && !propdata->bounds[i]->nonconvex )
      {
         if( !indicatorVarIsInteresting(scip, vars[i], (int)propdata->bounds[i]->nonconvex, (int)propdata->bounds[i]->indicator, propdata->indicatorthreshold) )
         {
            /* mark bound as filtered */
            propdata->bounds[i]->filtered = TRUE;
 
            /* increase number of filtered variables */
            ntotalfiltered++;
         }
      }
   }
 
   /*
    * 2.) Try first to filter lower bounds of interesting variables, whose bounds are not already filtered
    */
 
   for( i = 0; i < nvars; i++ )
   {
      SCIP_CALL( SCIPchgVarObjProbing(scip, vars[i], 0.0) );
   }
 
   for( i = 0; i < propdata->nbounds; i++ )
   {
      if( propdata->bounds[i]->boundtype == SCIP_BOUNDTYPE_LOWER && !propdata->bounds[i]->filtered
            && !propdata->bounds[i]->done )
      {
         SCIP_Real objcoef;
 
         objcoef = getFilterCoef(scip, propdata, propdata->bounds[i]->var, SCIP_BOUNDTYPE_LOWER);
 
         if( !SCIPisZero(scip, objcoef) )
         {
            SCIP_CALL( SCIPchgVarObjProbing(scip, propdata->bounds[i]->var, objcoef) );
 
            /* store nontrivial objective coefficients */
            objcoefs[nobjcoefs] = objcoef;
            objcoefsinds[nobjcoefs] = i;
            ++nobjcoefs;
         }
      }
   }
 
   do
   {
      SCIPdebugMsg(scip, "doing a lower bounds round\n");
      SCIP_CALL( filterRound(scip, propdata, nleftiterations, &nfiltered, objcoefs, objcoefsinds, nobjcoefs) );
      ntotalfiltered += nfiltered;
      SCIPdebugMsg(scip, "filtered %d more bounds in lower bounds round\n", nfiltered);
 
      /* update iterations left */
      nleftiterations = getIterationsLeft(scip, nolditerations, itlimit);
   }
   while( nfiltered >= propdata->nminfilter && ( nleftiterations == -1 ||  nleftiterations > 0 ) );
 
   /*
    * 3.) Now try to filter the remaining upper bounds of interesting variables, whose bounds are not already filtered
    */
 
   /* set all objective coefficients to zero */
   for( i = 0; i < nobjcoefs; i++ )
   {
      BOUND* bound;
 
      assert(objcoefsinds[i] >= 0 && objcoefsinds[i] < propdata->nbounds);
      bound = propdata->bounds[ objcoefsinds[i] ];
      assert(bound != NULL);
      SCIP_CALL( SCIPchgVarObjProbing(scip, bound->var, 0.0) );
   }
 
   /* reset number of nontrivial objective coefficients */
   nobjcoefs = 0;
 
#ifndef NDEBUG
   for( i = 0; i < nvars; ++i )
      assert(SCIPisZero(scip, SCIPgetVarObjProbing(scip, vars[i])));
#endif
 
   for( i = 0; i < propdata->nbounds; i++ )
   {
      if( propdata->bounds[i]->boundtype == SCIP_BOUNDTYPE_UPPER && !propdata->bounds[i]->filtered )
      {
         SCIP_Real objcoef;
 
         objcoef = getFilterCoef(scip, propdata, propdata->bounds[i]->var, SCIP_BOUNDTYPE_UPPER);
 
         if( !SCIPisZero(scip, objcoef) )
         {
            SCIP_CALL( SCIPchgVarObjProbing(scip, propdata->bounds[i]->var, objcoef) );
 
            /* store nontrivial objective coefficients */
            objcoefs[nobjcoefs] = objcoef;
            objcoefsinds[nobjcoefs] = i;
            ++nobjcoefs;
         }
      }
   }
 
   do
   {
      SCIPdebugMsg(scip, "doing an upper bounds round\n");
      SCIP_CALL( filterRound(scip, propdata, nleftiterations, &nfiltered, objcoefs, objcoefsinds, nobjcoefs) );
      SCIPdebugMsg(scip, "filtered %d more bounds in upper bounds round\n", nfiltered);
      ntotalfiltered += nfiltered;
      /* update iterations left */
      nleftiterations = getIterationsLeft(scip, nolditerations, itlimit);
   }
   while( nfiltered >= propdata->nminfilter && ( nleftiterations == -1 ||  nleftiterations > 0 ) );
 
   SCIPdebugMsg(scip, "filtered %d this round\n", ntotalfiltered);
   propdata->nfiltered += ntotalfiltered;
 
   /* free array */
   SCIPfreeBufferArray(scip, &objcoefsinds);
   SCIPfreeBufferArray(scip, &objcoefs);
 
   return SCIP_OKAY;
}
 
/** applies possible bound changes that were found */
static
SCIP_RETCODE applyBoundChgs(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROPDATA*        propdata,           /**< data of the obbt propagator */
   SCIP_RESULT*          result              /**< result pointer */
   )
{
#ifdef SCIP_DEBUG
   int ntightened;                           /* stores the number of successful bound changes */
#endif
   int i;
 
   assert(scip != NULL);
   assert(!SCIPinProbing(scip));
   assert(propdata != NULL);
   assert(result != NULL);
   assert(*result == SCIP_DIDNOTFIND);
 
   SCIPdebug( ntightened = 0 );
 
   for( i = 0; i < propdata->nbounds; i++ )
   {
      BOUND* bound;                          /* shortcut to the current bound */
      SCIP_Bool infeas;                      /* stores wether a tightening approach forced an infeasibilty */
      SCIP_Bool tightened;                   /* stores wether a tightening approach was successful */
 
      bound = propdata->bounds[i];
      infeas = FALSE;
 
      if( bound->found )
      {
         SCIPdebug( double oldbound = (bound->boundtype == SCIP_BOUNDTYPE_LOWER)
            ? SCIPvarGetLbLocal(bound->var)
            : SCIPvarGetUbLocal(bound->var) );
 
         if( bound->boundtype == SCIP_BOUNDTYPE_LOWER )
         {
            SCIP_CALL( SCIPtightenVarLb(scip, bound->var, bound->newval, FALSE, &infeas, &tightened) );
         }
         else
         {
            /* tighten only if new bound is small enough due to numerical reasons */
            if( SCIPisLE(scip, bound->newval, propdata->indicatorthreshold) )
            {
               SCIP_CALL( SCIPtightenVarUb(scip, bound->var, bound->newval, FALSE, &infeas, &tightened) );
            }
            else
               tightened = FALSE;
         }
 
         /* handle information about the success */
         if( infeas )
         {
            *result = SCIP_CUTOFF;
            SCIPdebugMsg(scip, "cut off\n");
            break;
         }
 
         if( tightened )
         {
            SCIPdebug( SCIPdebugMsg(scip, "tightended: %s old: %e new: %e\n" , SCIPvarGetName(bound->var), oldbound,
                  bound->newval) );
 
            *result = SCIP_REDUCEDDOM;
            SCIPdebug( ntightened++ );
         }
      }
   }
 
   SCIPdebug( SCIPdebugMsg(scip, "tightened bounds: %d\n", ntightened) );
 
   return SCIP_OKAY;
}
 
/** tries to tighten a bound in probing mode  */
static
SCIP_RETCODE tightenBoundProbing(
   SCIP*                 scip,               /**< SCIP data structure */
   BOUND*                bound,              /**< bound that could be tightened */
   SCIP_Real             newval,             /**< new bound value */
   SCIP_Bool*            tightened           /**< was tightening successful? */
   )
{
   SCIP_Real lb;
   SCIP_Real ub;
 
   assert(scip != NULL);
   assert(SCIPinProbing(scip));
   assert(bound != NULL);
   assert(tightened != NULL);
 
   *tightened = FALSE;
 
   /* get old bounds */
   lb = SCIPvarGetLbLocal(bound->var);
   ub = SCIPvarGetUbLocal(bound->var);
 
   if( bound->boundtype == SCIP_BOUNDTYPE_LOWER )
   {
      /* round bounds new value if variable is integral */
      if( SCIPvarIsIntegral(bound->var) )
         newval = SCIPceil(scip, newval);
 
      /* ensure that we give consistent bounds to the LP solver */
      if( newval > ub )
         newval = ub;
 
      /* tighten if really better */
      if( SCIPisLbBetter(scip, newval, lb, ub) )
      {
         SCIP_CALL( SCIPchgVarLbProbing(scip, bound->var, newval) );
         *tightened = TRUE;
      }
   }
   else
   {
      /* round bounds new value if variable is integral */
      if( SCIPvarIsIntegral(bound->var) )
         newval = SCIPfloor(scip, newval);
 
      /* ensure that we give consistent bounds to the LP solver */
      if( newval < lb )
         newval = lb;
 
      /* tighten if really better */
      if( SCIPisUbBetter(scip, newval, lb, ub) )
      {
         SCIP_CALL( SCIPchgVarUbProbing(scip, bound->var, newval) );
         *tightened = TRUE;
      }
   }
 
   return SCIP_OKAY;
}
 
/** comparison method for two bounds w.r.t. their scores */
static
SCIP_DECL_SORTPTRCOMP(compBoundsScore)
{
   BOUND* bound1 = (BOUND*) elem1;
   BOUND* bound2 = (BOUND*) elem2;
 
   return bound1->score == bound2->score ? 0 : ( bound1->score > bound2->score ? 1 : -1 );
}
 
/** comparison method for two bilinear term bounds w.r.t. their scores */
static
SCIP_DECL_SORTPTRCOMP(compBilinboundsScore)
{
   BILINBOUND* bound1 = (BILINBOUND*) elem1;
   BILINBOUND* bound2 = (BILINBOUND*) elem2;
 
   return bound1->score == bound2->score ? 0 : ( bound1->score > bound2->score ? 1 : -1 ); /*lint !e777*/
}
 
/** comparison method for two bounds w.r.t. their boundtype */
static
SCIP_DECL_SORTPTRCOMP(compBoundsBoundtype)
{
   int diff;
   BOUND* bound1 = (BOUND*) elem1;
   BOUND* bound2 = (BOUND*) elem2;
 
   /* prioritize undone bounds */
   diff = (!bound1->done ? 1 : 0) - (!bound2->done ? 1 : 0);
   if( diff != 0 )
      return diff;
 
   /* prioritize unfiltered bounds */
   diff = (!bound1->filtered ? 1 : 0) - (!bound2->filtered ? 1 : 0);
   if( diff != 0 )
      return diff;
 
   diff = (bound1->boundtype == SCIP_BOUNDTYPE_LOWER ? 1 : 0) - (bound2->boundtype == SCIP_BOUNDTYPE_LOWER ? 1 : 0);
   if( diff != 0 )  /* cppcheck-suppress knownConditionTrueFalse */
      return diff;
 
   return (bound1->score == bound2->score) ? 0 : (bound1->score > bound2->score ? 1 : -1);
}
 
/** sort the propdata->bounds array with their distance or their boundtype key */
static
SCIP_RETCODE sortBounds(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROPDATA*        propdata            /**< propagator data */
   )
{
   assert(scip != NULL);
   assert(propdata != NULL);
 
   SCIPdebugMsg(scip, "sort bounds\n");
   SCIPsortDownPtr((void**) propdata->bounds, compBoundsBoundtype, propdata->nbounds);
 
   return SCIP_OKAY;
}
 
/** evaluates a bound for the current LP solution */
static
SCIP_Real evalBound(
   SCIP*                 scip,
   BOUND*                bound
   )
{
   assert(scip != NULL);
   assert(bound != NULL);
 
   if( bound->boundtype == SCIP_BOUNDTYPE_LOWER )
      return REALABS( SCIPvarGetLPSol(bound->var) - SCIPvarGetLbLocal(bound->var) );
   else
      return REALABS( SCIPvarGetUbLocal(bound->var) - SCIPvarGetLPSol(bound->var) );
}
 
/** returns the index of the next undone and unfiltered bound with the smallest distance */
static
int nextBound(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROPDATA*        propdata,           /**< data of the obbt propagator */
   SCIP_Bool             convexphase         /**< consider only convex variables? */
   )
{
   SCIP_Real bestval;
   int bestidx;
   int k;
 
   assert(scip != NULL);
   assert(propdata != NULL);
 
   bestidx = -1;
   bestval = SCIPinfinity(scip);
 
   for( k = 0; k <= propdata->lastidx; ++k )
   {
      BOUND* tmpbound;
      tmpbound = propdata->bounds[k];
 
      assert(tmpbound != NULL);
 
      /* variables of indicator constraints are considered as nonconvex */
      if( !tmpbound->filtered && !tmpbound->done && (tmpbound->nonconvex == !convexphase || tmpbound->indicator == !convexphase) )
      {
         SCIP_Real boundval;
 
         /* return the next bound which is not done or unfiltered yet */
         if( propdata->orderingalgo == 0 )
            return k;
 
         boundval = evalBound(scip, tmpbound);
 
         /* negate boundval if we use the reverse greedy algorithm */
         boundval = (propdata->orderingalgo == 2) ? -1.0 * boundval : boundval;
 
         if( bestidx == -1 || boundval < bestval )
         {
            bestidx = k;
            bestval = boundval;
         }
      }
   }
 
   return bestidx;  /*lint !e438*/
}
 
/** try to separate the solution of the last OBBT LP in order to learn better variable bounds; we apply additional
 *  separation rounds as long as the routine finds better bounds; because of dual degeneracy we apply a minimum number of
 *  separation rounds
 */
static
SCIP_RETCODE applySeparation(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROPDATA*        propdata,           /**< data of the obbt propagator */
   BOUND*                currbound,          /**< current bound */
   SCIP_Longint*         nleftiterations,    /**< number of left iterations (-1 for no limit) */
   SCIP_Bool*            success             /**< pointer to store if we have found a better bound */
   )
{
   SCIP_Bool inroot;
   int i;
 
   assert(nleftiterations != NULL);
   assert(success != NULL);
   assert(SCIPinProbing(scip));
 
   *success = FALSE;
 
   /* check if we are originally in the root node */
   inroot = SCIPgetDepth(scip) == 1;
 
   for( i = 0; i <= propdata->sepamaxiter; ++i )
   {
      SCIPdebug( SCIP_Longint nlpiter; )
      SCIP_Real oldval;
      SCIP_Bool cutoff;
      SCIP_Bool delayed;
      SCIP_Bool error;
      SCIP_Bool optimal;
      SCIP_Bool tightened;
 
      oldval = SCIPvarGetLPSol(currbound->var);
 
      /* find and store cuts to separate the current LP solution */
      SCIP_CALL( SCIPseparateSol(scip, NULL, inroot, TRUE, FALSE, &delayed, &cutoff) );
      SCIPdebugMsg(scip, "applySeparation() - ncuts = %d\n", SCIPgetNCuts(scip));
 
      /* leave if we did not found any cut */
      if( SCIPgetNCuts(scip) == 0 )
         break;
 
      /* apply cuts and resolve LP */
      SCIP_CALL( SCIPapplyCutsProbing(scip, &cutoff) );
      assert(SCIPgetNCuts(scip) == 0);
      SCIPdebug( nlpiter = SCIPgetNLPIterations(scip); )
      SCIP_CALL( solveLP(scip, (int) *nleftiterations, &error, &optimal) );
      SCIPdebug( nlpiter = SCIPgetNLPIterations(scip) - nlpiter; )
      SCIPdebug( SCIPdebugMsg(scip, "applySeparation() - optimal=%u error=%u lpiter=%" SCIP_LONGINT_FORMAT "\n", optimal, error, nlpiter); )
      SCIPdebugMsg(scip, "oldval = %e newval = %e\n", oldval, SCIPvarGetLPSol(currbound->var));
 
      /* leave if we did not solve the LP to optimality or an error occured */
      if( error || !optimal )
         break;
 
      /* try to generate a genvbound */
      if( inroot && propdata->genvboundprop != NULL && propdata->genvbdsduringsepa )
      {
         SCIP_Bool found;
         SCIP_CALL( createGenVBound(scip, propdata, currbound, &found) );
         propdata->ngenvboundsprobing += found ? 1 : 0;
      }
 
      /* try to tight the variable bound */
      tightened = FALSE;
      if( !SCIPisEQ(scip, oldval, SCIPvarGetLPSol(currbound->var)) )
      {
         SCIP_CALL( tightenBoundProbing(scip, currbound, SCIPvarGetLPSol(currbound->var), &tightened) );
         SCIPdebugMsg(scip, "apply separation - tightened=%u oldval=%e newval=%e\n", tightened, oldval,
            SCIPvarGetLPSol(currbound->var));
 
         *success |= tightened;
      }
 
      /* leave the separation if we did not tighten the bound and proceed at least propdata->sepaminiter iterations */
      if( !tightened && i >= propdata->sepaminiter )
         break;
   }
 
   return SCIP_OKAY;
}
 
/** finds new variable bounds until no iterations left or all bounds have been checked */
static
SCIP_RETCODE findNewBounds(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROPDATA*        propdata,           /**< data of the obbt propagator */
   SCIP_Longint*         nleftiterations,    /**< pointer to store the number of left iterations */
   SCIP_Bool             convexphase         /**< consider only convex variables? */
   )
{
   SCIP_Longint nolditerations;
   SCIP_Bool iterationsleft;
   BOUND* currbound;
   SCIP_Longint itlimit;
   int nextboundidx;
 
   assert(scip != NULL);
   assert(propdata != NULL);
   assert(nleftiterations != NULL);
 
   /* update the number of left iterations */
   nolditerations = SCIPgetNLPIterations(scip);
   itlimit = *nleftiterations;
   assert(*nleftiterations == getIterationsLeft(scip, nolditerations, itlimit));
   iterationsleft = (*nleftiterations == -1) || (*nleftiterations > 0);
 
   /* To improve the performance we sort the bound in such a way that the undone and
    * unfiltered bounds are at the end of propdata->bounds. We calculate and update
    * the position of the last unfiltered and undone bound in propdata->lastidx
    */
   if( !convexphase )
   {
      /* sort bounds */
      SCIP_CALL( sortBounds(scip, propdata) );
 
      /* if the first bound is filtered or done then there is no bound left */
      if( propdata->bounds[0]->done || propdata->bounds[0]->filtered )
      {
         SCIPdebugMsg(scip, "no unprocessed/unfiltered bound left\n");
         return SCIP_OKAY;
      }
 
      /* compute the last undone and unfiltered node */
      propdata->lastidx = 0;
      while( propdata->lastidx < propdata->nbounds - 1 && !propdata->bounds[propdata->lastidx]->done &&
            !propdata->bounds[propdata->lastidx]->filtered )
         ++propdata->lastidx;
 
      SCIPdebugMsg(scip, "lastidx = %d\n", propdata->lastidx);
   }
 
   /* find the first unprocessed bound */
   nextboundidx = nextBound(scip, propdata, convexphase);
 
   /* skip if there is no bound left */
   if( nextboundidx == -1 )
   {
      SCIPdebugMsg(scip, "no unprocessed/unfiltered bound left\n");
      return SCIP_OKAY;
   }
 
   currbound = propdata->bounds[nextboundidx];
   assert(!currbound->done && !currbound->filtered);
 
   /* main loop */
   while( iterationsleft &&  !SCIPisStopped(scip) )
   {
      SCIP_Bool optimal;
      SCIP_Bool error;
      int nfiltered;
 
      assert(currbound != NULL);
      assert(currbound->done == FALSE);
      assert(currbound->filtered == FALSE);
 
      /* do not visit currbound more than once */
      currbound->done = TRUE;
      exchangeBounds(propdata, nextboundidx);
 
      /* set objective for curr */
      SCIP_CALL( setObjProbing(scip, propdata, currbound, 1.0) );
 
      SCIPdebugMsg(scip, "before solving      Boundtype: %d , LB: %e , UB: %e\n",
         currbound->boundtype == SCIP_BOUNDTYPE_LOWER, SCIPvarGetLbLocal(currbound->var),
         SCIPvarGetUbLocal(currbound->var) );
      SCIPdebugMsg(scip, "before solving      var <%s>, LP value: %f\n",
         SCIPvarGetName(currbound->var), SCIPvarGetLPSol(currbound->var));
 
      SCIPdebugMsg(scip, "probing iterations before solve: %lld \n", SCIPgetNLPIterations(scip));
 
      propdata->nprobingiterations -= SCIPgetNLPIterations(scip);
 
      /* now solve the LP */
      SCIP_CALL( solveLP(scip, (int) *nleftiterations, &error, &optimal) );
 
      propdata->nprobingiterations += SCIPgetNLPIterations(scip);
      propdata->nsolvedbounds++;
 
      SCIPdebugMsg(scip, "probing iterations after solve: %lld \n", SCIPgetNLPIterations(scip));
      SCIPdebugMsg(scip, "OPT: %u ERROR: %u\n" , optimal, error);
      SCIPdebugMsg(scip, "after solving      Boundtype: %d , LB: %e , UB: %e\n",
         currbound->boundtype == SCIP_BOUNDTYPE_LOWER, SCIPvarGetLbLocal(currbound->var),
         SCIPvarGetUbLocal(currbound->var) );
      SCIPdebugMsg(scip, "after solving      var <%s>, LP value: %f\n",
         SCIPvarGetName(currbound->var), SCIPvarGetLPSol(currbound->var));
 
      /* update nleftiterations */
      *nleftiterations = getIterationsLeft(scip, nolditerations, itlimit);
      iterationsleft = (*nleftiterations == -1) || (*nleftiterations > 0);
 
      if( error )
      {
         SCIPdebugMsg(scip, "ERROR during LP solving\n");
 
         /* set the objective of currbound to zero to null the whole objective; otherwise the objective is wrong when
          * we call findNewBounds() for the convex phase
          */
         SCIP_CALL( SCIPchgVarObjProbing(scip, currbound->var, 0.0) );
 
         return SCIP_OKAY;
      }
 
      if( optimal )
      {
         SCIP_Bool success;
 
         currbound->newval = SCIPvarGetLPSol(currbound->var);
         currbound->found = TRUE;
 
         /* in root node we may want to create a genvbound (independent of tightening success) */
         if( (SCIPgetDepth(scip) == 0 || (SCIPinProbing(scip) && SCIPgetDepth(scip) == 1))
               && propdata->genvboundprop != NULL )
         {
            SCIP_Bool found;
 
            SCIP_CALL( createGenVBound(scip, propdata, currbound, &found) );
 
            if( found )
               propdata->ngenvboundsprobing += 1;
         }
 
         /* try to tighten bound in probing mode */
         success = FALSE;
         if( propdata->tightintboundsprobing && SCIPvarIsIntegral(currbound->var) )
         {
            SCIPdebugMsg(scip, "tightening bound %s = %e bounds: [%e, %e]\n", SCIPvarGetName(currbound->var),
                currbound->newval, SCIPvarGetLbLocal(currbound->var), SCIPvarGetUbLocal(currbound->var) );
            SCIP_CALL( tightenBoundProbing(scip, currbound, currbound->newval, &success) );
            SCIPdebugMsg(scip, "tightening bound %s\n", success ? "successful" : "not successful");
         }
         else if( propdata->tightcontboundsprobing && !SCIPvarIsIntegral(currbound->var) )
         {
            SCIPdebugMsg(scip, "tightening bound %s = %e bounds: [%e, %e]\n", SCIPvarGetName(currbound->var),
               currbound->newval, SCIPvarGetLbLocal(currbound->var), SCIPvarGetUbLocal(currbound->var) );
            SCIP_CALL( tightenBoundProbing(scip, currbound, currbound->newval, &success) );
            SCIPdebugMsg(scip, "tightening bound %s\n", success ? "successful" : "not successful");
         }
 
         /* separate current OBBT LP solution */
         if( iterationsleft && propdata->separatesol )
         {
            propdata->nprobingiterations -= SCIPgetNLPIterations(scip);
            SCIP_CALL( applySeparation(scip, propdata, currbound, nleftiterations, &success) );
            propdata->nprobingiterations += SCIPgetNLPIterations(scip);
 
            /* remember best solution value after solving additional separations LPs */
            if( success )
            {
#ifndef NDEBUG
               SCIP_Real newval = SCIPvarGetLPSol(currbound->var);
 
               /* round new bound if the variable is integral */
               if( SCIPvarIsIntegral(currbound->var) )
                  newval = currbound->boundtype == SCIP_BOUNDTYPE_LOWER ?
                     SCIPceil(scip, newval) : SCIPfloor(scip, newval);
 
               assert((currbound->boundtype == SCIP_BOUNDTYPE_LOWER &&
                     SCIPisGT(scip, newval, currbound->newval))
                  || (currbound->boundtype == SCIP_BOUNDTYPE_UPPER &&
                     SCIPisLT(scip, newval, currbound->newval)));
#endif
 
               currbound->newval = SCIPvarGetLPSol(currbound->var);
            }
         }
 
         /* filter bound candidates by using the current LP solution */
         if( propdata->applytrivialfilter )
         {
            SCIP_CALL( filterExistingLP(scip, propdata, &nfiltered, currbound) );
            SCIPdebugMsg(scip, "filtered %d bounds via inspecting present LP solution\n", nfiltered);
            propdata->ntrivialfiltered += nfiltered;
         }
 
         propdata->propagatecounter += success ? 1 : 0;
 
         /* propagate if we have found enough bound tightenings */
         if( propdata->propagatefreq != 0 && propdata->propagatecounter >= propdata->propagatefreq )
         {
            SCIP_Longint ndomredsfound;
            SCIP_Bool cutoff;
 
            SCIP_CALL( SCIPpropagateProbing(scip, 0, &cutoff, &ndomredsfound) );
            SCIPdebugMsg(scip, "propagation - cutoff %u  ndomreds %" SCIP_LONGINT_FORMAT "\n", cutoff, ndomredsfound);
 
            propdata->npropagatedomreds += ndomredsfound;
            propdata->propagatecounter = 0;
         }
      }
 
      /* set objective to zero */
      SCIP_CALL( setObjProbing(scip, propdata, currbound, 0.0) );
 
      /* find the first unprocessed bound */
      nextboundidx = nextBound(scip, propdata, convexphase);
 
      /* check if there is no unprocessed and unfiltered node left */
      if( nextboundidx == -1 )
      {
         SCIPdebugMsg(scip, "NO unvisited/unfiltered bound left!\n");
         break;
      }
 
      currbound = propdata->bounds[nextboundidx];
      assert(!currbound->done && !currbound->filtered);
   }
 
   if( iterationsleft )
   {
      SCIPdebugMsg(scip, "still iterations left: %" SCIP_LONGINT_FORMAT "\n", *nleftiterations);
   }
   else
   {
      SCIPdebugMsg(scip, "no iterations left\n");
   }
 
   return SCIP_OKAY;
}
 
 
/** main function of obbt */
static
SCIP_RETCODE applyObbt(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROPDATA*        propdata,           /**< data of the obbt propagator */
   SCIP_Longint          itlimit,            /**< LP iteration limit (-1: no limit) */
   SCIP_RESULT*          result              /**< result pointer */
   )
{
   SCIP_VAR** vars;
   SCIP_Real* oldlbs;
   SCIP_Real* oldubs;
   SCIP_Longint lastnpropagatedomreds;
   SCIP_Longint nleftiterations;
   SCIP_Real oldconditionlimit;
   SCIP_Real oldboundstreps;
   SCIP_Real olddualfeastol;
   SCIP_Bool hasconditionlimit;
   SCIP_Bool continuenode;
   SCIP_Bool boundleft;
   int oldpolishing;
   int nfiltered;
   int nvars;
   int i;
 
   assert(scip != NULL);
   assert(propdata != NULL);
   assert(itlimit == -1 || itlimit >= 0);
 
   SCIPdebugMsg(scip, "apply obbt\n");
 
   oldlbs = NULL;
   oldubs = NULL;
   lastnpropagatedomreds = propdata->npropagatedomreds;
   nleftiterations = itlimit;
   continuenode = SCIPnodeGetNumber(SCIPgetCurrentNode(scip)) == propdata->lastnode;
   propdata->lastidx = -1;
   boundleft = FALSE;
   *result = SCIP_DIDNOTFIND;
 
   /* store old variable bounds if we use propagation during obbt */
   if( propdata->propagatefreq > 0 )
   {
      SCIP_CALL( SCIPallocBufferArray(scip, &oldlbs, propdata->nbounds) );
      SCIP_CALL( SCIPallocBufferArray(scip, &oldubs, propdata->nbounds) );
   }
 
   /* reset bound data structure flags; fixed variables are marked as filtered */
   for( i = 0; i < propdata->nbounds; i++ )
   {
      BOUND* bound = propdata->bounds[i];
      bound->found = FALSE;
 
      /* store old variable bounds */
      if( oldlbs != NULL && oldubs != NULL )
      {
         oldlbs[bound->index] = SCIPvarGetLbLocal(bound->var);
         oldubs[bound->index] = SCIPvarGetUbLocal(bound->var);
      }
 
      /* reset 'done' and 'filtered' flag in a new B&B node */
      if( !continuenode )
      {
         bound->done = FALSE;
         bound->filtered = FALSE;
      }
 
      /* mark fixed variables as filtered */
      bound->filtered |= varIsFixedLocal(scip, bound->var);
 
      /* check for an unprocessed bound */
      if( !bound->filtered && !bound->done )
         boundleft = TRUE;
   }
 
   /* no bound left to check */
   if( !boundleft )
      goto TERMINATE;
 
   /* filter variables via inspecting present LP solution */
   if( propdata->applytrivialfilter && !continuenode )
   {
      SCIP_CALL( filterExistingLP(scip, propdata, &nfiltered, NULL) );
      SCIPdebugMsg(scip, "filtered %d bounds via inspecting present LP solution\n", nfiltered);
      propdata->ntrivialfiltered += nfiltered;
   }
 
   /* store old dualfeasibletol */
   olddualfeastol = SCIPdualfeastol(scip);
 
   /* start probing */
   SCIP_CALL( SCIPstartProbing(scip) );
   SCIPdebugMsg(scip, "start probing\n");
 
   /* tighten dual feastol */
   if( propdata->dualfeastol < olddualfeastol )
   {
      SCIP_CALL( SCIPchgDualfeastol(scip, propdata->dualfeastol) );
   }
 
   /* tighten condition limit */
   hasconditionlimit = (SCIPgetRealParam(scip, "lp/conditionlimit", &oldconditionlimit) == SCIP_OKAY);
   if( !hasconditionlimit )
   {
      SCIPwarningMessage(scip, "obbt propagator could not set condition limit in LP solver - running without\n");
   }
   else if( propdata->conditionlimit > 0.0 && (oldconditionlimit < 0.0 || propdata->conditionlimit < oldconditionlimit) )
   {
      SCIP_CALL( SCIPsetRealParam(scip, "lp/conditionlimit", propdata->conditionlimit) );
   }
 
   /* tighten relative bound improvement limit */
   SCIP_CALL( SCIPgetRealParam(scip, "numerics/boundstreps", &oldboundstreps) );
   if( !SCIPisEQ(scip, oldboundstreps, propdata->boundstreps) )
   {
     SCIP_CALL( SCIPsetRealParam(scip, "numerics/boundstreps", propdata->boundstreps) );
   }
 
   /* add objective cutoff */
   SCIP_CALL( addObjCutoff(scip, propdata) );
 
   /* deactivate LP polishing */
   SCIP_CALL( SCIPgetIntParam(scip, "lp/solutionpolishing", &oldpolishing) );
   SCIP_CALL( SCIPsetIntParam(scip, "lp/solutionpolishing", 0) );
 
   /* apply filtering */
   if( propdata->applyfilterrounds )
   {
      SCIP_CALL( filterBounds(scip, propdata, nleftiterations) );
   }
 
   /* set objective coefficients to zero */
   vars = SCIPgetVars(scip);
   nvars = SCIPgetNVars(scip);
   for( i = 0; i < nvars; ++i )
   {
      /* note that it is not possible to change the objective of non-column variables during probing; we have to take
       * care of the objective contribution of loose variables in createGenVBound()
       */
      if( SCIPvarGetObj(vars[i]) != 0.0 && SCIPvarGetStatus(vars[i]) == SCIP_VARSTATUS_COLUMN )
      {
         SCIP_CALL( SCIPchgVarObjProbing(scip, vars[i], 0.0) );
      }
   }
 
   /* find new bounds for the variables */
   SCIP_CALL( findNewBounds(scip, propdata, &nleftiterations, FALSE) );
 
   if( nleftiterations > 0 || itlimit < 0 )
   {
      SCIP_CALL( findNewBounds(scip, propdata, &nleftiterations, TRUE) );
   }
 
   /* reset dual feastol and condition limit */
   SCIP_CALL( SCIPchgDualfeastol(scip, olddualfeastol) );
   if( hasconditionlimit )
   {
      SCIP_CALL( SCIPsetRealParam(scip, "lp/conditionlimit", oldconditionlimit) );
   }
 
   /* update bound->newval if we have learned additional bound tightenings during SCIPpropagateProbing() */
   if( oldlbs != NULL && oldubs != NULL && propdata->npropagatedomreds - lastnpropagatedomreds > 0 )
   {
      assert(propdata->propagatefreq > 0);
      for( i = 0; i < propdata->nbounds; ++i )
      {
         BOUND* bound = propdata->bounds[i];
 
         /* it might be the case that a bound found by the additional propagation is better than the bound found after solving an OBBT
          * LP
          */
         if( bound->found )
         {
            if( bound->boundtype == SCIP_BOUNDTYPE_LOWER )
               bound->newval = MAX(bound->newval, SCIPvarGetLbLocal(bound->var)); /*lint !e666*/
            else
               bound->newval = MIN(bound->newval, SCIPvarGetUbLocal(bound->var)); /*lint !e666*/
         }
         else
         {
            SCIP_Real oldlb;
            SCIP_Real oldub;
 
            oldlb = oldlbs[bound->index];
            oldub = oldubs[bound->index];
 
            if( bound->boundtype == SCIP_BOUNDTYPE_LOWER && SCIPisLbBetter(scip, SCIPvarGetLbLocal(bound->var), oldlb, oldub) )
            {
               SCIPdebugMsg(scip, "tighter lower bound due to propagation: %d - %e -> %e\n", i, oldlb, SCIPvarGetLbLocal(bound->var));
               bound->newval = SCIPvarGetLbLocal(bound->var);
               bound->found = TRUE;
            }
 
            if( bound->boundtype == SCIP_BOUNDTYPE_UPPER && SCIPisUbBetter(scip, SCIPvarGetUbLocal(bound->var), oldlb, oldub) )
            {
               SCIPdebugMsg(scip, "tighter upper bound due to propagation: %d - %e -> %e\n", i, oldub, SCIPvarGetUbLocal(bound->var));
               bound->newval = SCIPvarGetUbLocal(bound->var);
               bound->found = TRUE;
            }
         }
      }
   }
 
   /* reset relative bound improvement limit */
   SCIP_CALL( SCIPsetRealParam(scip, "numerics/boundstreps", oldboundstreps) );
 
   /* reset original LP polishing setting */
   SCIP_CALL( SCIPsetIntParam(scip, "lp/solutionpolishing", oldpolishing) );
 
   /* end probing */
   SCIP_CALL( SCIPendProbing(scip) );
   SCIPdebugMsg(scip, "end probing!\n");
 
   /* release cutoff row if there is one */
   if( propdata->cutoffrow != NULL )
   {
      assert(!SCIProwIsInLP(propdata->cutoffrow));
      SCIP_CALL( SCIPreleaseRow(scip, &(propdata->cutoffrow)) );
   }
 
   /* apply buffered bound changes */
   SCIP_CALL( applyBoundChgs(scip, propdata, result) );
 
TERMINATE:
   SCIPfreeBufferArrayNull(scip, &oldubs);
   SCIPfreeBufferArrayNull(scip, &oldlbs);
 
   return SCIP_OKAY;
}
 
/** computes a valid inequality from the current LP relaxation for a bilinear term xy only involving x and y; the
 *  inequality is found by optimizing along the line connecting the points (xs,ys) and (xt,yt) over the currently given
 *  linear relaxation of the problem; this optimization problem is an LP
 *
 *  max lambda
 *  s.t. Ax <= b
 *       (x,y) = (xs,ys) + lambda ((xt,yt) - (xs,ys))
 *       lambda in [0,1]
 *
 *  which is equivalent to
 *
 *  max x
 *  s.t. (1) Ax <= b
 *       (2) (x - xs) / (xt - xs) = (y - ys) / (yt - ys)
 *
 *  Let x* be the optimal primal and (mu,theta) be the optimal dual solution of this LP. The KKT conditions imply that
 *  the aggregation of the linear constraints mu*Ax <= mu*b can be written as
 *
 *  x * (1 - theta) / (xt - xs) + y * theta / (yt - ys) = mu * Ax <= mu * b
 *
 *  <=> alpha * x + beta * y <= mu * b = alpha * (x*) + beta * (y*)
 *
 *  which is a valid inequality in the (x,y)-space; in order to avoid numerical difficulties when (xs,ys) is too close
 *  to (xt,yt), we scale constraint (1) by max{1,|xt-xs|,|yt-ys|} beforehand
 */
static
SCIP_RETCODE solveBilinearLP(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_VAR*             x,                  /**< first variable */
   SCIP_VAR*             y,                  /**< second variable */
   SCIP_Real             xs,                 /**< x-coordinate of the first point */
   SCIP_Real             ys,                 /**< y-coordinate of the first point */
   SCIP_Real             xt,                 /**< x-coordinate of the second point */
   SCIP_Real             yt,                 /**< y-coordinate of the second point */
   SCIP_Real*            xcoef,              /**< pointer to store the coefficient of x */
   SCIP_Real*            ycoef,              /**< pointer to store the coefficient of y */
   SCIP_Real*            constant,           /**< pointer to store the constant */
   SCIP_Longint          iterlim,            /**< iteration limit (-1: for no limit) */
   int*                  nnonzduals          /**< buffer to store the number of non-zero dual multipliers except for
                                              *   the auxiliary row (NULL if not needed) */
   )
{
   SCIP_ROW* row;
   SCIP_Real signx;
   SCIP_Real scale;
   SCIP_Real side;
   SCIP_Bool lperror;
 
   assert(xcoef != NULL);
   assert(ycoef != NULL);
   assert(constant != NULL);
   assert(SCIPinProbing(scip));
 
   *xcoef = SCIP_INVALID;
   *ycoef = SCIP_INVALID;
   *constant= SCIP_INVALID;
   if( nnonzduals != NULL )
      *nnonzduals = 0;
 
   SCIPdebugMsg(scip, "   solve bilinear LP for (%s,%s) from (%g,%g) to (%g,%g)\n", SCIPvarGetName(x), SCIPvarGetName(y), xs,
      ys, xt, yt);
 
   /* skip computations if (xs,ys) and (xt,yt) are too close to each other or contain too large values */
   if( SCIPisFeasEQ(scip, xs, xt) || SCIPisFeasEQ(scip, ys, yt)
      || SCIPisHugeValue(scip, REALABS(xs)) || SCIPisHugeValue(scip, REALABS(xt))
      || SCIPisHugeValue(scip, REALABS(ys)) || SCIPisHugeValue(scip, REALABS(yt)) )
   {
      SCIPdebugMsg(scip, "   -> skip: bounds are too close/large\n");
      return SCIP_OKAY;
   }
 
   /* compute scaler for the additional linear constraint */
   scale = MIN(MAX3(1.0, REALABS(xt-xs), REALABS(yt-ys)), 100.0); /*lint !e666*/
 
   /* set objective function */
   signx = (xs > xt) ? 1.0 : -1.0;
   SCIP_CALL( SCIPchgVarObjProbing(scip, x, signx) );
 
   /* create new probing node to remove the added LP row afterwards */
   SCIP_CALL( SCIPnewProbingNode(scip) );
 
   /* create row */
   side = scale * (xs/(xt-xs) - ys/(yt-ys));
   SCIP_CALL( SCIPcreateEmptyRowUnspec(scip, &row, "bilinrow", side, side, FALSE, FALSE, TRUE) );
   SCIP_CALL( SCIPaddVarToRow(scip, row, x, scale/(xt-xs)) );
   SCIP_CALL( SCIPaddVarToRow(scip, row, y, -scale/(yt-ys)) );
   SCIP_CALL( SCIPaddRowProbing(scip, row) );
 
   /* solve probing LP */
#ifdef NDEBUG
   {
      SCIP_RETCODE retstat;
      retstat = SCIPsolveProbingLP(scip, iterlim, &lperror, NULL);
      if( retstat != SCIP_OKAY )
      {
         SCIPwarningMessage(scip, "Error while solving LP in quadratic constraint handler; LP solve terminated with" \
            "code <%d>\n", retstat);
      }
   }
#else
   SCIP_CALL( SCIPsolveProbingLP(scip, (int)iterlim, &lperror, NULL) ); /*lint !e712*/
#endif
 
   SCIPdebugMsg(scip, "   solved probing LP -> lperror=%u lpstat=%d\n", lperror, SCIPgetLPSolstat(scip));
 
   /* collect dual and primal solution entries */
   if( !lperror  && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL )
   {
      SCIP_Real xval = SCIPvarGetLPSol(x);
      SCIP_Real yval = SCIPvarGetLPSol(y);
      SCIP_Real mu = -SCIProwGetDualsol(row);
 
      SCIPdebugMsg(scip, "   primal=(%g,%g) dual=%g\n", xval, yval, mu);
 
      /* xcoef x + ycoef y <= constant */
      *xcoef  = -signx - (mu * scale) / (xt - xs);
      *ycoef = (mu * scale) / (yt - ys);
      *constant = (*xcoef) * xval + (*ycoef) * yval;
 
      /* xcoef x <= -ycoef y + constant */
      *ycoef = -(*ycoef);
 
      /* inequality is only useful when both coefficients are different from zero; normalize inequality if possible */
      if( !SCIPisFeasZero(scip, *xcoef) && !SCIPisFeasZero(scip, *ycoef) )
      {
         SCIP_Real val = REALABS(*xcoef);
         int r;
 
         *xcoef /= val;
         *ycoef /= val;
         *constant /= val;
 
         if( SCIPisZero(scip, *constant) )
            *constant = 0.0;
 
         if( nnonzduals != NULL )
         {
            /* count the number of non-zero dual multipliers except for the added row */
            for( r = 0; r < SCIPgetNLPRows(scip); ++r )
            {
               if( SCIPgetLPRows(scip)[r] != row && !SCIPisFeasZero(scip, SCIProwGetDualsol(SCIPgetLPRows(scip)[r])) )
                  ++(*nnonzduals);
            }
         }
      }
      else
      {
         *xcoef = SCIP_INVALID;
         *ycoef = SCIP_INVALID;
         *constant = SCIP_INVALID;
      }
   }
 
   /* release row and backtrack probing node */
   SCIP_CALL( SCIPreleaseRow(scip, &row) );
   SCIP_CALL( SCIPbacktrackProbing(scip, 0) );
 
   /* reset objective function */
   SCIP_CALL( SCIPchgVarObjProbing(scip, x, 0.0) );
 
   return SCIP_OKAY;
}
 
/* applies obbt for finding valid inequalities for bilinear terms; function works as follows:
 *
 *  1. start probing mode
 *  2. add objective cutoff (if possible)
 *  3. set objective function to zero
 *  4. set feasibility, optimality, and relative bound improvement tolerances of SCIP
 *  5. main loop
 *  6. restore old tolerances
 *
 */
static
SCIP_RETCODE applyObbtBilinear(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROPDATA*        propdata,           /**< data of the obbt propagator */
   SCIP_Longint          itlimit,            /**< LP iteration limit (-1: no limit) */
   SCIP_RESULT*          result              /**< result pointer */
   )
{
   SCIP_VAR** vars;
   SCIP_Real oldfeastol;
   SCIP_Bool lperror;
   SCIP_Longint nolditerations;
   SCIP_Longint nleftiterations;
   SCIP_CONSHDLR* conshdlr;
   SCIP_NLHDLR* bilinearnlhdlr;
   int nvars;
   int i;
 
   assert(scip != NULL);
   assert(propdata != NULL);
   assert(itlimit == -1 || itlimit >= 0);
   assert(result != NULL);
 
   if( propdata->nbilinbounds <= 0 || SCIPgetDepth(scip) != 0 || propdata->lastbilinidx >= propdata->nbilinbounds )
      return SCIP_OKAY;
 
   SCIPdebugMsg(scip, "call applyObbtBilinear starting from %d\n", propdata->lastbilinidx);
 
   /* find nonlinear handler for bilinear terms */
   conshdlr = SCIPfindConshdlr(scip, "nonlinear");
   bilinearnlhdlr = conshdlr != NULL ? SCIPfindNlhdlrNonlinear(conshdlr, "bilinear") : NULL;
 
   /* no nonlinear handler available -> skip */
   if( bilinearnlhdlr == NULL )
      return SCIP_OKAY;
 
   vars = SCIPgetVars(scip);
   nvars = SCIPgetNVars(scip);
 
   nolditerations = SCIPgetNLPIterations(scip);
   nleftiterations = getIterationsLeft(scip, nolditerations, itlimit);
   SCIPdebugMsg(scip, "iteration limit: %lld\n", nleftiterations);
 
   /* 1. start probing */
   SCIP_CALL( SCIPstartProbing(scip) );
 
   /* 2. add objective cutoff */
   SCIP_CALL( addObjCutoff(scip, propdata) );
 
   /* 3. set objective function to zero */
   for( i = 0; i < nvars; ++i )
   {
      SCIP_CALL( SCIPchgVarObjProbing(scip, vars[i], 0.0) );
   }
 
   /* 4. tighten LP feasibility tolerance to be at most feastol/10.0 */
   oldfeastol = SCIPchgRelaxfeastol(scip, SCIPfeastol(scip) / 10.0);
 
   /* we need to solve the probing LP before creating new probing nodes in solveBilinearLP() */
   SCIP_CALL( SCIPsolveProbingLP(scip, (int)nleftiterations, &lperror, NULL) );
 
   if( lperror )
      goto TERMINATE;
 
   /* 5. main loop */
   for( i = propdata->lastbilinidx; i < propdata->nbilinbounds
      && (nleftiterations > 0 || nleftiterations == -1)
      && (propdata->itlimitbilin < 0 || propdata->itlimitbilin > propdata->itusedbilin )
      && !SCIPisStopped(scip); ++i )
   {
      CORNER corners[4] = {LEFTBOTTOM, LEFTTOP, RIGHTTOP, RIGHTBOTTOM};
      BILINBOUND* bilinbound;
      int k;
 
      bilinbound = propdata->bilinbounds[i];
      assert(bilinbound != NULL);
 
      SCIPdebugMsg(scip, "process %d: %s %s done=%u filtered=%d nunderest=%d noverest=%d\n", i,
         SCIPvarGetName(bilinboundGetX(bilinbound)), SCIPvarGetName(bilinboundGetY(bilinbound)), bilinbound->done,
         bilinbound->filtered, bilinboundGetLocksNeg(bilinbound), bilinboundGetLocksPos(bilinbound));
 
      /* we already solved LPs for this bilinear term */
      if( bilinbound->done || bilinbound->filtered == (int)FILTERED )
         continue;
 
      /* iterate through all corners
       *
       * 0: (xs,ys)=(ubx,lby) (xt,yt)=(lbx,uby) -> underestimate
       * 1: (xs,ys)=(ubx,uby) (xt,yt)=(lbx,lby) -> overestimate
       * 2: (xs,ys)=(lbx,uby) (xt,yt)=(ubx,lby) -> underestimate
       * 3: (xs,ys)=(lbx,lby) (xt,yt)=(ubx,uby) -> overestimate
       */
      for( k = 0; k < 4; ++k )
      {
         CORNER corner = corners[k];
         SCIP_VAR* x = bilinboundGetX(bilinbound);
         SCIP_VAR* y = bilinboundGetY(bilinbound);
         SCIP_Real xcoef;
         SCIP_Real ycoef;
         SCIP_Real constant;
         SCIP_Real xs = SCIP_INVALID;
         SCIP_Real ys = SCIP_INVALID;
         SCIP_Real xt = SCIP_INVALID;
         SCIP_Real yt = SCIP_INVALID;
         int nnonzduals = 0;
 
         /* skip corners that lead to an under- or overestimate that is not needed */
         if( ((corner == LEFTTOP || corner == RIGHTBOTTOM) && bilinboundGetLocksPos(bilinbound) == 0)
            || ((corner == LEFTBOTTOM || corner == RIGHTTOP) && bilinboundGetLocksNeg(bilinbound) == 0) )
            continue;
 
         /* check whether corner has been filtered already */
         if( (bilinbound->filtered & corner) != 0 ) /*lint !e641*/
            continue;
 
         /* get corners (xs,ys) and (xt,yt) */
         getCorners(x, y, corner, &xs, &ys, &xt, &yt);
 
         /* skip target corner points with too large values */
         if( SCIPisHugeValue(scip, REALABS(xt)) || SCIPisHugeValue(scip, REALABS(yt)) )
            continue;
 
         /* compute inequality */
         propdata->itusedbilin -= SCIPgetNLPIterations(scip);
         SCIP_CALL( solveBilinearLP(scip, x, y, xs, ys, xt, yt, &xcoef, &ycoef, &constant, -1L,
            propdata->createlincons ? &nnonzduals : NULL) ); /*lint !e826*/
         propdata->itusedbilin += SCIPgetNLPIterations(scip);
 
         /* update number of LP iterations */
         nleftiterations = getIterationsLeft(scip, nolditerations, itlimit);
         SCIPdebugMsg(scip, "LP iterations left: %lld\n", nleftiterations);
 
         /* add inequality to quadratic constraint handler if it separates (xt,yt) */
         if( !SCIPisHugeValue(scip, xcoef)  && !SCIPisFeasZero(scip, xcoef)
            && REALABS(ycoef) < 1e+3 && REALABS(ycoef) > 1e-3 /* TODO add a parameter for this */
            && SCIPisFeasGT(scip, (xcoef*xt - ycoef*yt - constant) / sqrt(SQR(xcoef) + SQR(ycoef) + SQR(constant)), 1e-2) )
         {
            SCIP_Bool success;
 
            /* add inequality to the associated product expression */
            SCIP_CALL( SCIPaddIneqBilinear(scip, bilinearnlhdlr, bilinbound->expr, xcoef, ycoef,
               constant, &success) );
 
            /* check whether the inequality has been accepted */
            if( success )
            {
               *result = SCIP_REDUCEDDOM;
               SCIPdebugMsg(scip, "   found %g x <= %g y + %g with violation %g\n", xcoef, ycoef, constant,
                  (xcoef*xt - ycoef*yt - constant) / sqrt(SQR(xcoef) + SQR(ycoef) + SQR(constant)));
 
               /* create a linear constraint that is only used for propagation */
               if( propdata->createlincons && nnonzduals > 1 )
               {
                  SCIP_CONS* cons;
                  char name[SCIP_MAXSTRLEN];
                  SCIP_VAR* linvars[2] = {x, y};
                  SCIP_Real linvals[2] = {xcoef, -ycoef};
                  SCIP_Real rhs = constant;
 
                  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "bilincons_%s_%s", SCIPvarGetName(x), SCIPvarGetName(y));
                  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, name, 2, linvars, linvals, -SCIPinfinity(scip), rhs,
                     FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE) );
 
                  SCIP_CALL( SCIPaddCons(scip, cons) );
                  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
               }
            }
         }
      }
 
      /* mark the bound as processed */
      bilinbound->done = TRUE;
   }
 
   /* remember last unprocessed bilinear term */
   propdata->lastbilinidx = i;
 
  TERMINATE:
   /* end probing */
   SCIP_CALL( SCIPendProbing(scip) );
 
   /* release cutoff row if there is one */
   if( propdata->cutoffrow != NULL )
   {
      assert(!SCIProwIsInLP(propdata->cutoffrow));
      SCIP_CALL( SCIPreleaseRow(scip, &(propdata->cutoffrow)) );
   }
 
   /* 6. restore old tolerance */
   (void) SCIPchgRelaxfeastol(scip, oldfeastol);
 
   return SCIP_OKAY;
}
 
/** computes the score of a bound */
static
unsigned int getScore(
   SCIP*                 scip,               /**< SCIP data structure */
   BOUND*                bound,              /**< pointer of bound */
   int                   nlcount,            /**< number of nonlinear constraints containing the bounds variable */
   int                   nindcount,          /**< number of indicator constraints containing the bounds variable */
   int                   maxnlcount,         /**< maximal number of nonlinear and indicator constraints a variable appears in */
   SCIP_Real             smallub             /**< variables with upper bound smaller than this value are counted in half iff part of indicator constraints */
   )
{
   SCIP_Real counter;
   unsigned int score;                       /* score to be computed */
 
   assert(scip != NULL);
   assert(bound != NULL);
   assert(nlcount >= 0);
   assert(nindcount >= 0);
   assert(maxnlcount >= nlcount + nindcount);
 
   counter = nlcount;
   if( nindcount > 0 )
   {
      /* variables with small upper bound are counted in half
       * since the probability is high that the corresponding indicator constraint is already reformulated as bigM constraint
       */
      if( SCIPvarGetUbLocal(bound->var) <= smallub )
         counter += 0.5 * nindcount;
      else
         counter += nindcount;
   }
 
   /* score = ( nlcount * ( BASE - 1 ) / maxnlcount ) * BASE^2 + vartype * BASE + boundtype */
   score = (unsigned int) ( counter > 0 ? (OBBT_SCOREBASE * counter * ( OBBT_SCOREBASE - 1 )) / maxnlcount : 0 ); /*lint !e414*/
   switch( SCIPvarGetType(bound->var) )
   {
   case SCIP_VARTYPE_INTEGER:
      score += 1;
      break;
   case SCIP_VARTYPE_IMPLINT:
      score += 2;
      break;
   case SCIP_VARTYPE_CONTINUOUS:
      score += 3;
      break;
   case SCIP_VARTYPE_BINARY:
      score += 4;
      break;
   default:
      break;
   }
 
   score *= OBBT_SCOREBASE;
   if( bound->boundtype == SCIP_BOUNDTYPE_UPPER )
      score += 1;
 
   return score;
}
 
/** count how often each variable is used in a nonconvex term */
static
SCIP_RETCODE getNLPVarsNonConvexity(
   SCIP*                 scip,               /**< SCIP data structure */
   unsigned int*         nccounts            /**< store the number each variable appears in a
                                              *   non-convex term */
   )
{
   SCIP_CONSHDLR* conshdlr;
   SCIP_HASHMAP* var2expr;
   int nvars;
   int i;
 
   assert(scip != NULL);
   assert(nccounts != NULL);
 
   nvars = SCIPgetNVars(scip);
 
   /* initialize nccounts to zero */
   BMSclearMemoryArray(nccounts, nvars);
 
   /* get nonlinear constraint handler */
   conshdlr = SCIPfindConshdlr(scip, "nonlinear");
   if( conshdlr == NULL || SCIPconshdlrGetNConss(conshdlr) == 0 )
      return SCIP_OKAY;
 
   var2expr = SCIPgetVarExprHashmapNonlinear(conshdlr);
   assert(var2expr != NULL);
 
   for( i = 0; i < SCIPgetNVars(scip); ++i )
   {
      SCIP_VAR* var;
 
      var = SCIPgetVars(scip)[i];
      assert(var != NULL);
 
      if( SCIPhashmapExists(var2expr, (void*) var) )
      {
         SCIP_EXPR* expr = (SCIP_EXPR*)SCIPhashmapGetImage(var2expr, (void*) var);
         assert(expr != NULL);
         assert(SCIPisExprVar(scip, expr));
 
         nccounts[SCIPvarGetProbindex(var)] = SCIPgetExprNSepaUsesActivityNonlinear(expr);
      }
   }
 
#ifdef SCIP_DEBUG
   for( i = 0; i < SCIPgetNVars(scip); ++i)
   {
      SCIP_VAR* var = SCIPgetVars(scip)[i];
      assert(var != NULL);
      SCIPdebugMsg(scip, "nccounts[%s] = %u\n", SCIPvarGetName(var), nccounts[SCIPvarGetProbindex(var)]);
   }
#endif
 
   return SCIP_OKAY;
}
 
/** computes for each variable the number of indicator constraints in which the variable appears */
static
SCIP_RETCODE getNVarsIndicators(
   SCIP*                 scip,               /**< SCIP data structure */
   int*                  nindcount           /**< array that stores in how many indicator conss each variable appears */
   )
{
   SCIP_CONSHDLR* conshdlr;
   SCIP_CONS** indconss;
   int nvars;
   int nindconss;
 
   assert(scip != NULL);
   assert(nindcount != NULL);
 
   nvars = SCIPgetNVars(scip);
 
   /* initialize nindcount to zero */
   BMSclearMemoryArray(nindcount, nvars);
 
   /* get indicator constraint handler */
   conshdlr = SCIPfindConshdlr(scip, "indicator");
   if( conshdlr == NULL || SCIPconshdlrGetNConss(conshdlr) == 0 )
      return SCIP_OKAY;
 
   nindconss = SCIPconshdlrGetNConss(conshdlr);
   indconss = SCIPconshdlrGetConss(conshdlr);
 
   for( int i = 0; i < nindconss; ++i )
   {
      SCIP_VAR** consvars;
      SCIP_VAR* slackvar;
      SCIP_CONS* lincons;
      int nconsvars;
 
      lincons = SCIPgetLinearConsIndicator(indconss[i]);
      assert(lincons!=NULL);
 
      nconsvars = SCIPgetNVarsLinear(scip, lincons);
      consvars = SCIPgetVarsLinear(scip, lincons);
      assert(consvars != NULL);
      slackvar = SCIPgetSlackVarIndicator(indconss[i]);
 
      for( int v = 0; v < nconsvars; ++v )
      {
         /* we should skip the slackvariable */
         if( consvars[v] == slackvar )
            continue;
 
         /* consider only active variables */
         if( SCIPvarGetStatus(consvars[v]) == SCIP_VARSTATUS_COLUMN )
         {
            nindcount[SCIPvarGetProbindex(consvars[v])] += 1;
         }
      }
   }
 
   return SCIP_OKAY;
}
 
/** determines whether a variable is interesting */
static
SCIP_Bool varIsInteresting(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_VAR*             var,                /**< variable to check */
   int                   nlcount,            /**< number of nonlinear constraints containing the variable
                                              *   or number of non-convex terms containing the variable
                                              *  (depends on propdata->onlynonconvexvars)  */
   int                   nindcount           /**< number of indicator constraints containing the variable
                                              *   or 0 (depends on propdata->indicators) */
   )
{
   assert(SCIPgetDepth(scip) == 0);
 
   return !SCIPvarIsBinary(var) && SCIPvarGetStatus(var) == SCIP_VARSTATUS_COLUMN && (nlcount > 0 || nindcount > 0)
      && !varIsFixedLocal(scip, var);
}
 
/** initializes interesting bounds */
static
SCIP_RETCODE initBounds(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROPDATA*        propdata            /**< data of the obbt propagator */
   )
{
   SCIP_CONSHDLR* conshdlr;
   SCIP_VAR** vars;                          /* array of the problems variables */
   int* nlcount;                             /* array that stores in how many nonlinearities each variable appears */
   int* nindcount;                           /* array that stores in how many indicator constraints each variable appears */
   unsigned int* nccount;                    /* array that stores in how many nonconvexities each variable appears */
   SCIP_Real maxcouplingvalue;
   SCIP_Real sepacouplingvalue;
   SCIP_Real smallub;
 
   int bdidx;                                /* bound index inside propdata->bounds */
   int maxnlcount;                           /* maximal number of nonlinear and indicator constraints a variable appears in */
   int nvars;                                /* number of the problems variables */
   int i;
 
   assert(scip != NULL);
   assert(propdata != NULL);
   assert(SCIPisNLPConstructed(scip) || propdata->indicators);
 
   SCIPdebugMsg(scip, "initialize bounds\n");
 
   /* get variable data */
   SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
 
   SCIP_CALL( SCIPallocBufferArray(scip, &nlcount, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &nccount, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &nindcount, nvars) );
 
   /* count nonlinearities */
   if( SCIPisNLPConstructed(scip) )
   {
      assert(SCIPgetNNLPVars(scip) == nvars);
      SCIP_CALL( SCIPgetNLPVarsNonlinearity(scip, nlcount) );
      SCIP_CALL( getNLPVarsNonConvexity(scip, nccount) );
   }
   else
   {
      /* initialize to zero */
      BMSclearMemoryArray(nlcount, nvars);
      BMSclearMemoryArray(nccount, nvars);
   }
 
   /* count indicators */
   if( propdata->indicators )
   {
      SCIP_CALL( getNVarsIndicators(scip, nindcount) );
   }
   else
   {
      /* initialize to zero */
      BMSclearMemoryArray(nindcount, nvars);
   }
 
   maxnlcount = 0;
   for( i = 0; i < nvars; i++ )
   {
      if( maxnlcount < (nlcount[i] + nindcount[i]) )
         maxnlcount = nlcount[i] + nindcount[i];
   }
 
   /* allocate interesting bounds array */
   propdata->boundssize = 2 * nvars;
   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(propdata->bounds), 2 * nvars) );
 
   SCIP_CALL( SCIPgetRealParam(scip, "constraints/indicator/maxcouplingvalue", &maxcouplingvalue) );
   SCIP_CALL( SCIPgetRealParam(scip, "constraints/indicator/sepacouplingvalue", &sepacouplingvalue) );
 
   smallub = MIN(maxcouplingvalue, sepacouplingvalue);
 
   /* get all interesting variables and their bounds */
   bdidx = 0;
   for( i = 0; i < nvars; i++ )
   {
      if( varIsInteresting(scip, vars[i], (propdata->onlynonconvexvars ? (int)nccount[i] : nlcount[i]), (propdata->indicators ? nindcount[i] : 0))
         && indicatorVarIsInteresting(scip, vars[i], (propdata->onlynonconvexvars ? (int)nccount[i] : nlcount[i]), (propdata->indicators ? nindcount[i] : 0), propdata->indicatorthreshold) )
      {
         BOUND** bdaddress;
 
         /* create lower bound */
         bdaddress = &(propdata->bounds[bdidx]);
         SCIP_CALL( SCIPallocBlockMemory(scip, bdaddress) );
         propdata->bounds[bdidx]->boundtype = SCIP_BOUNDTYPE_LOWER;
         propdata->bounds[bdidx]->var = vars[i];
         propdata->bounds[bdidx]->found = FALSE;
         propdata->bounds[bdidx]->filtered = FALSE;
         propdata->bounds[bdidx]->newval = 0.0;
         propdata->bounds[bdidx]->score = getScore(scip, propdata->bounds[bdidx], nlcount[i], nindcount[i], maxnlcount, smallub);
         propdata->bounds[bdidx]->done = FALSE;
         propdata->bounds[bdidx]->nonconvex = (nccount[i] > 0);
         propdata->bounds[bdidx]->indicator = (nindcount[i] > 0);
         propdata->bounds[bdidx]->index = bdidx;
         bdidx++;
 
         /* create upper bound */
         bdaddress = &(propdata->bounds[bdidx]);
         SCIP_CALL( SCIPallocBlockMemory(scip, bdaddress) );
         propdata->bounds[bdidx]->boundtype = SCIP_BOUNDTYPE_UPPER;
         propdata->bounds[bdidx]->var = vars[i];
         propdata->bounds[bdidx]->found = FALSE;
         propdata->bounds[bdidx]->filtered = FALSE;
         propdata->bounds[bdidx]->newval = 0.0;
         propdata->bounds[bdidx]->score = getScore(scip, propdata->bounds[bdidx], nlcount[i], nindcount[i], maxnlcount, smallub);
         propdata->bounds[bdidx]->done = FALSE;
         propdata->bounds[bdidx]->nonconvex = (nccount[i] > 0);
         propdata->bounds[bdidx]->indicator = (nindcount[i] > 0);
         propdata->bounds[bdidx]->index = bdidx;
         bdidx++;
      }
   }
 
   /* set number of interesting bounds */
   propdata->nbounds = bdidx;
 
   conshdlr = SCIPfindConshdlr(scip, "nonlinear");
 
   /* get all product expressions from nonlinear constraint handler */
   if( propdata->nbounds > 0 && conshdlr != NULL && propdata->createbilinineqs )
   {
      SCIP_NLHDLR* bilinnlhdlr;
      SCIP_EXPR** exprs;
      int nexprs;
 
      /* find nonlinear handler for bilinear terms */
      bilinnlhdlr = SCIPfindNlhdlrNonlinear(conshdlr, "bilinear");
      assert(bilinnlhdlr != NULL);
 
      /* collect all bilinear product in all nonlinear constraints */
      exprs = SCIPgetExprsBilinear(bilinnlhdlr);
      nexprs = SCIPgetNExprsBilinear(bilinnlhdlr);
 
      if( nexprs > 0 )
      {
         SCIP_CALL( SCIPallocBlockMemoryArray(scip, &propdata->bilinbounds, nexprs) );
         propdata->bilinboundssize = nexprs;
         propdata->nbilinbounds = 0;
 
         /* store candidates as bilinear bounds */
         for( i = 0; i < nexprs; ++i )
         {
            BILINBOUND* bilinbound;
            SCIP_VAR* x;
            SCIP_VAR* y;
 
            assert(exprs[i] != NULL);
            assert(SCIPexprGetNChildren(exprs[i]) == 2);
 
            x = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(exprs[i])[0]);
            y = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(exprs[i])[1]);
            assert(x != NULL);
            assert(y != NULL);
            assert(x != y);
 
            /* skip almost fixed variables */
            if( !varIsInteresting(scip, x, 1, 0) || !varIsInteresting(scip, y, 1, 0) )
               continue;
 
            /* create bilinear bound */
            SCIP_CALL( SCIPallocBlockMemory(scip, &propdata->bilinbounds[propdata->nbilinbounds]) ); /*lint !e866*/
            bilinbound = propdata->bilinbounds[propdata->nbilinbounds];
            BMSclearMemory(bilinbound);
 
            /* store and capture expression */
            bilinbound->expr = exprs[i];
            SCIPcaptureExpr(bilinbound->expr);
 
            /* compute a descent score */
            bilinbound->score = bilinboundGetScore(scip, propdata->randnumgen, bilinbound);
 
            /* increase the number of bilinear bounds */
            ++(propdata->nbilinbounds);
 
            SCIPdebugMsg(scip, "added bilinear bound for %s %s\n", SCIPvarGetName(x), SCIPvarGetName(y));
         }
      }
 
      /* sort bounds according to decreasing score */
      if( propdata->nbilinbounds > 1 )
      {
         SCIPsortDownPtr((void**) propdata->bilinbounds, compBilinboundsScore, propdata->nbilinbounds);
      }
   }
 
   /* free memory for buffering nonlinearities */
   assert(nlcount != NULL);
   assert(nccount != NULL);
   assert(nindcount != NULL);
   SCIPfreeBufferArray(scip, &nindcount);
   SCIPfreeBufferArray(scip, &nccount);
   SCIPfreeBufferArray(scip, &nlcount);
 
   /*  propdata->bounds array if empty */
   if( propdata->nbounds <= 0 )
   {
      assert(propdata->nbounds == 0);
      assert(propdata->boundssize >= 0 );
      SCIPfreeBlockMemoryArray(scip, &(propdata->bounds), propdata->boundssize);
   }
 
   SCIPdebugMsg(scip, "problem has %d/%d interesting bounds\n", propdata->nbounds, 2 * nvars);
 
   if( propdata->nbounds > 0 )
   {
      /* sort bounds according to decreasing score; although this initial order will be overruled by the distance
       * criterion later, gives a more well-defined starting situation for OBBT and might help to reduce solver
       * variability
       */
      SCIPsortDownPtr((void**) propdata->bounds, compBoundsScore, propdata->nbounds);
   }
 
   return SCIP_OKAY;
}
 
/*
 * Callback methods of propagator
 */
 
/** copy method for propagator plugins (called when SCIP copies plugins)
 *
 *  @note The UG framework assumes that all default plug-ins of SCIP implement a copy callback. We check
 *  SCIPgetSubscipDepth() in PROPEXEC to prevent the propagator to run in a sub-SCIP.
 */
static
SCIP_DECL_PROPCOPY(propCopyObbt)
{  /*lint --e{715}*/
   assert(scip != NULL);
   assert(prop != NULL);
   assert(strcmp(SCIPpropGetName(prop), PROP_NAME) == 0);
 
   /* call inclusion method of constraint handler */
   SCIP_CALL( SCIPincludePropObbt(scip) );
 
   return SCIP_OKAY;
}
 
/** solving process initialization method of propagator (called when branch and bound process is about to begin) */
static
SCIP_DECL_PROPINITSOL(propInitsolObbt)
{  /*lint --e{715}*/
   SCIP_PROPDATA* propdata;
 
   assert(scip != NULL);
   assert(prop != NULL);
   assert(strcmp(SCIPpropGetName(prop), PROP_NAME) == 0);
 
   /* get propagator data */
   propdata = SCIPpropGetData(prop);
   assert(propdata != NULL);
 
   propdata->bounds = NULL;
   propdata->nbounds = -1;
   propdata->boundssize = 0;
   propdata->cutoffrow = NULL;
   propdata->lastnode = -1;
 
   /* if genvbounds propagator is not available, we cannot create genvbounds */
   propdata->genvboundprop = propdata->creategenvbounds ? SCIPfindProp(scip, GENVBOUND_PROP_NAME) : NULL;
 
   SCIPdebugMsg(scip, "creating genvbounds: %s\n", propdata->genvboundprop != NULL ? "true" : "false");
 
   /* create random number generator */
   SCIP_CALL( SCIPcreateRandom(scip, &propdata->randnumgen, DEFAULT_RANDSEED, TRUE) );
 
   return SCIP_OKAY;
}
 
/** execution method of propagator */
static
SCIP_DECL_PROPEXEC(propExecObbt)
{  /*lint --e{715}*/
   SCIP_PROPDATA* propdata;
   SCIP_Longint itlimit;
 
   assert(scip != NULL);
   assert(prop != NULL);
   assert(strcmp(SCIPpropGetName(prop), PROP_NAME) == 0);
 
   *result = SCIP_DIDNOTRUN;
 
   /* do not run in: presolving, repropagation, probing mode, if no objective propagation is allowed  */
   if( SCIPgetStage(scip) != SCIP_STAGE_SOLVING || SCIPinRepropagation(scip) || SCIPinProbing(scip) || !SCIPallowWeakDualReds(scip) )
      return SCIP_OKAY;
 
   /* do not run propagator in a sub-SCIP */
   if( SCIPgetSubscipDepth(scip) > 0 )
      return SCIP_OKAY;
 
   /* get propagator data */
   propdata = SCIPpropGetData(prop);
   assert(propdata != NULL);
 
   /* only run for nonlinear problems, i.e., if NLP is constructed
    * or if indicator constraints exists and should be considered
    */
   if( !SCIPisNLPConstructed(scip)
      && (!propdata->indicators || SCIPfindConshdlr(scip, "indicator") == NULL || SCIPconshdlrGetNConss(SCIPfindConshdlr(scip, "indicator")) == 0) )
   {
      SCIPdebugMsg(scip, "NLP not constructed and no indicator constraints available, skipping obbt\n");
      return SCIP_OKAY;
   }
 
   /* only run if LP all columns are in the LP, i.e., the LP is a relaxation; e.g., do not run if pricers are active
    * since pricing is not performed in probing mode
    */
   if( !SCIPallColsInLP(scip) )
   {
      SCIPdebugMsg(scip, "not all columns in LP, skipping obbt\n");
      return SCIP_OKAY;
   }
 
   /* ensure that bounds are initialized */
   if( propdata->nbounds == -1 )
   {
      /* bounds must be initialized at root node */
      if( SCIPgetDepth(scip) == 0 )
      {
         SCIP_CALL( initBounds(scip, propdata) );
      }
      else
      {
         assert(!SCIPinProbing(scip));
         return SCIP_OKAY;
      }
   }
   assert(propdata->nbounds >= 0);
 
   /* do not run if there are no interesting bounds */
   /**@todo disable */
   if( propdata->nbounds <= 0 )
   {
      SCIPdebugMsg(scip, "there are no interesting bounds\n");
      return SCIP_OKAY;
   }
 
   /* only run once in a node != root */
   if( SCIPgetDepth(scip) > 0 && SCIPnodeGetNumber(SCIPgetCurrentNode(scip)) == propdata->lastnode )
   {
      return SCIP_OKAY;
   }
 
   SCIPdebugMsg(scip, "applying obbt for problem <%s> at depth %d\n", SCIPgetProbName(scip), SCIPgetDepth(scip));
 
   /* without an optimal LP solution we don't want to run; this may be because propagators with higher priority have
    * already found reductions or numerical troubles occured during LP solving
    */
   if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL && SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_UNBOUNDEDRAY )
   {
      SCIPdebugMsg(scip, "aborting since no optimal LP solution is at hand\n");
      return SCIP_OKAY;
   }
 
   /* compute iteration limit */
   if( propdata->itlimitfactor > 0.0 )
      itlimit = (SCIP_Longint) MAX(propdata->itlimitfactor * SCIPgetNRootLPIterations(scip),
         propdata->minitlimit); /*lint !e666*/
   else
      itlimit = -1;
 
   /* apply obbt */
   SCIP_CALL( applyObbt(scip, propdata, itlimit, result) );
   assert(*result != SCIP_DIDNOTRUN);
 
   /* compute globally inequalities for bilinear terms */
   if( propdata->createbilinineqs )
   {
      /* set LP iteration limit */
      if( propdata->itlimitbilin == 0L )
      {
         /* no iteration limit if itlimit < 0 or itlimitfactorbilin < 0 */
         propdata->itlimitbilin = (itlimit < 0 || propdata->itlimitfactorbilin < 0)
            ? -1L : (SCIP_Longint)(itlimit * propdata->itlimitfactorbilin);
      }
 
      SCIP_CALL( applyObbtBilinear(scip, propdata, itlimit, result) );
   }
 
   /* set current node as last node */
   propdata->lastnode = SCIPnodeGetNumber(SCIPgetCurrentNode(scip));
 
   return SCIP_OKAY;
}
 
/** propagation conflict resolving method of propagator */
static
SCIP_DECL_PROPRESPROP(propRespropObbt)
{  /*lint --e{715}*/
   *result = SCIP_DIDNOTFIND;
 
   return SCIP_OKAY;
}
 
/** solving process deinitialization method of propagator (called before branch and bound process data is freed) */
static
SCIP_DECL_PROPEXITSOL(propExitsolObbt)
{  /*lint --e{715}*/
   SCIP_PROPDATA* propdata;
   int i;
 
   assert(scip != NULL);
   assert(prop != NULL);
   assert(strcmp(SCIPpropGetName(prop), PROP_NAME) == 0);
 
   /* get propagator data */
   propdata = SCIPpropGetData(prop);
   assert(propdata != NULL);
 
   /* free random number generator */
   SCIPfreeRandom(scip, &propdata->randnumgen);
   propdata->randnumgen = NULL;
 
   /* note that because we reset filtered flags to false at each call to obbt, the same bound may be filtered multiple
    * times
    */
   SCIPstatisticMessage("DIVE-LP: %" SCIP_LONGINT_FORMAT "  NFILTERED: %d NTRIVIALFILTERED: %d NSOLVED: %d "
      "FILTER-LP: %" SCIP_LONGINT_FORMAT " NGENVB(dive): %d NGENVB(aggr.): %d NGENVB(triv.) %d\n",
      propdata->nprobingiterations, propdata->nfiltered, propdata->ntrivialfiltered, propdata->nsolvedbounds,
      propdata->nfilterlpiters, propdata->ngenvboundsprobing, propdata->ngenvboundsaggrfil, propdata->ngenvboundstrivfil);
 
   /* free bilinear bounds */
   if( propdata->bilinboundssize > 0 )
   {
      for( i = propdata->nbilinbounds - 1; i >= 0; --i )
      {
         assert(propdata->bilinbounds[i] != NULL);
         assert(propdata->bilinbounds[i]->expr != NULL);
 
         /* release expression */
         SCIP_CALL( SCIPreleaseExpr(scip, &propdata->bilinbounds[i]->expr) );
 
         SCIPfreeBlockMemory(scip, &propdata->bilinbounds[i]); /*lint !e866*/
      }
      SCIPfreeBlockMemoryArray(scip, &propdata->bilinbounds, propdata->bilinboundssize);
      propdata->bilinboundssize = 0;
      propdata->nbilinbounds = 0;
   }
 
   /* free memory allocated for the bounds */
   if( propdata->nbounds > 0 )
   {
      /* free bounds */
      for( i = propdata->nbounds - 1; i >= 0; i-- )
      {
         SCIPfreeBlockMemory(scip, &(propdata->bounds[i])); /*lint !e866*/
      }
      SCIPfreeBlockMemoryArray(scip, &(propdata->bounds), propdata->boundssize);
   }
 
   propdata->nbounds = -1;
   propdata->itlimitbilin = 0;
   propdata->itusedbilin = 0;
 
   return SCIP_OKAY;
}
 
/** destructor of propagator to free user data (called when SCIP is exiting) */
static
SCIP_DECL_PROPFREE(propFreeObbt)
{  /*lint --e{715}*/
   SCIP_PROPDATA* propdata;
 
   assert(strcmp(SCIPpropGetName(prop), PROP_NAME) == 0);
 
   /* free propagator data */
   propdata = SCIPpropGetData(prop);
   assert(propdata != NULL);
 
   SCIPfreeBlockMemory(scip, &propdata);
 
   SCIPpropSetData(prop, NULL);
 
   return SCIP_OKAY;
}
 
/*
 * propagator specific interface methods
 */
 
/** creates the obbt propagator and includes it in SCIP */
SCIP_RETCODE SCIPincludePropObbt(
   SCIP*                 scip                /**< SCIP data structure */
   )
{
   SCIP_PROPDATA* propdata;
   SCIP_PROP* prop;
 
   /* create obbt propagator data */
   SCIP_CALL( SCIPallocBlockMemory(scip, &propdata) );
   BMSclearMemory(propdata);
 
   /* initialize statistic variables */
   propdata->nprobingiterations = 0;
   propdata->nfiltered = 0;
   propdata->ntrivialfiltered = 0;
   propdata->nsolvedbounds = 0;
   propdata->ngenvboundsprobing = 0;
   propdata->ngenvboundsaggrfil = 0;
   propdata->ngenvboundstrivfil = 0;
   propdata->nfilterlpiters = 0;
   propdata->lastidx = -1;
   propdata->propagatecounter = 0;
   propdata->npropagatedomreds = 0;
 
   /* include propagator */
   SCIP_CALL( SCIPincludePropBasic(scip, &prop, PROP_NAME, PROP_DESC, PROP_PRIORITY, PROP_FREQ, PROP_DELAY, PROP_TIMING,
         propExecObbt, propdata) );
 
   SCIP_CALL( SCIPsetPropCopy(scip, prop, propCopyObbt) );
   SCIP_CALL( SCIPsetPropFree(scip, prop, propFreeObbt) );
   SCIP_CALL( SCIPsetPropExitsol(scip, prop, propExitsolObbt) );
   SCIP_CALL( SCIPsetPropInitsol(scip, prop, propInitsolObbt) );
   SCIP_CALL( SCIPsetPropResprop(scip, prop, propRespropObbt) );
 
   SCIP_CALL( SCIPaddBoolParam(scip, "propagating/" PROP_NAME "/creategenvbounds",
         "should obbt try to provide genvbounds if possible?",
         &propdata->creategenvbounds, TRUE, DEFAULT_CREATE_GENVBOUNDS, NULL, NULL) );
 
   SCIP_CALL( SCIPaddBoolParam(scip, "propagating/" PROP_NAME "/normalize",
         "should coefficients in filtering be normalized w.r.t. the domains sizes?",
         &propdata->normalize, TRUE, DEFAULT_FILTERING_NORM, NULL, NULL) );
 
   SCIP_CALL( SCIPaddBoolParam(scip, "propagating/" PROP_NAME "/applyfilterrounds",
         "try to filter bounds in so-called filter rounds by solving auxiliary LPs?",
         &propdata->applyfilterrounds, TRUE, DEFAULT_APPLY_FILTERROUNDS, NULL, NULL) );
 
   SCIP_CALL( SCIPaddBoolParam(scip, "propagating/" PROP_NAME "/applytrivialfilter",
         "try to filter bounds with the LP solution after each solve?",
         &propdata->applytrivialfilter, TRUE, DEFAULT_APPLY_TRIVIALFITLERING, NULL, NULL) );
 
   SCIP_CALL( SCIPaddBoolParam(scip, "propagating/" PROP_NAME "/genvbdsduringfilter",
         "should we try to generate genvbounds during trivial and aggressive filtering?",
         &propdata->genvbdsduringfilter, TRUE, DEFAULT_GENVBDSDURINGFILTER, NULL, NULL) );
 
  SCIP_CALL( SCIPaddBoolParam(scip, "propagating/" PROP_NAME "/genvbdsduringsepa",
         "try to create genvbounds during separation process?",
         &propdata->genvbdsduringsepa, TRUE, DEFAULT_GENVBDSDURINGSEPA, NULL, NULL) );
 
   SCIP_CALL( SCIPaddIntParam(scip, "propagating/" PROP_NAME "/minfilter",
         "minimal number of filtered bounds to apply another filter round",
         &propdata->nminfilter, TRUE, DEFAULT_FILTERING_MIN, 1, INT_MAX, NULL, NULL) );
 
   SCIP_CALL( SCIPaddRealParam(scip, "propagating/" PROP_NAME "/itlimitfactor",
         "multiple of root node LP iterations used as total LP iteration limit for obbt (<= 0: no limit )",
         &propdata->itlimitfactor, FALSE, DEFAULT_ITLIMITFACTOR, SCIP_REAL_MIN, SCIP_REAL_MAX, NULL, NULL) );
 
   SCIP_CALL( SCIPaddRealParam(scip, "propagating/" PROP_NAME "/itlimitfactorbilin",
         "multiple of OBBT LP limit used as total LP iteration limit for solving bilinear inequality LPs (< 0 for no limit)",
         &propdata->itlimitfactorbilin, FALSE, DEFAULT_ITLIMITFAC_BILININEQS, SCIP_REAL_MIN, SCIP_REAL_MAX, NULL, NULL) );
 
   SCIP_CALL( SCIPaddRealParam(scip, "propagating/" PROP_NAME "/minnonconvexity",
         "minimum absolute value of nonconvex eigenvalues for a bilinear term",
         &propdata->minnonconvexity, FALSE, DEFAULT_MINNONCONVEXITY, 0.0, SCIP_REAL_MAX, NULL, NULL) );
 
   SCIP_CALL( SCIPaddLongintParam(scip, "propagating/" PROP_NAME "/minitlimit",
         "minimum LP iteration limit",
         &propdata->minitlimit, FALSE, DEFAULT_MINITLIMIT, 0L, SCIP_LONGINT_MAX, NULL, NULL) );
 
   SCIP_CALL( SCIPaddRealParam(scip, "propagating/" PROP_NAME "/dualfeastol",
         "feasibility tolerance for reduced costs used in obbt; this value is used if SCIP's dual feastol is greater",
         &propdata->dualfeastol, FALSE, DEFAULT_DUALFEASTOL, 0.0, SCIP_REAL_MAX, NULL, NULL) );
 
   SCIP_CALL( SCIPaddRealParam(scip, "propagating/" PROP_NAME "/conditionlimit",
         "maximum condition limit used in LP solver (-1.0: no limit)",
         &propdata->conditionlimit, FALSE, DEFAULT_CONDITIONLIMIT, -1.0, SCIP_REAL_MAX, NULL, NULL) );
 
   SCIP_CALL( SCIPaddRealParam(scip, "propagating/" PROP_NAME "/boundstreps",
         "minimal relative improve for strengthening bounds",
         &propdata->boundstreps, FALSE, DEFAULT_BOUNDSTREPS, 0.0, 1.0, NULL, NULL) );
 
   SCIP_CALL( SCIPaddRealParam(scip, "propagating/" PROP_NAME "/indicatorthreshold",
         "threshold whether upper bounds of vars of indicator conss are considered or tightened",
         &propdata->indicatorthreshold, TRUE, DEFAULT_INDICATORTHRESHOLD, 0.0, SCIP_REAL_MAX, NULL, NULL) );
 
  SCIP_CALL( SCIPaddBoolParam(scip, "propagating/" PROP_NAME "/onlynonconvexvars",
         "only apply obbt on non-convex variables",
         &propdata->onlynonconvexvars, TRUE, DEFAULT_ONLYNONCONVEXVARS, NULL, NULL) );
 
  SCIP_CALL( SCIPaddBoolParam(scip, "propagating/" PROP_NAME "/indicators",
         "apply obbt on variables of indicator constraints? (independent of convexity)",
         &propdata->indicators, TRUE, DEFAULT_INDICATORS, NULL, NULL) );
 
  SCIP_CALL( SCIPaddBoolParam(scip, "propagating/" PROP_NAME "/tightintboundsprobing",
         "should integral bounds be tightened during the probing mode?",
         &propdata->tightintboundsprobing, TRUE, DEFAULT_TIGHTINTBOUNDSPROBING, NULL, NULL) );
 
  SCIP_CALL( SCIPaddBoolParam(scip, "propagating/" PROP_NAME "/tightcontboundsprobing",
         "should continuous bounds be tightened during the probing mode?",
         &propdata->tightcontboundsprobing, TRUE, DEFAULT_TIGHTCONTBOUNDSPROBING, NULL, NULL) );
 
  SCIP_CALL( SCIPaddBoolParam(scip, "propagating/" PROP_NAME "/createbilinineqs",
         "solve auxiliary LPs in order to find valid inequalities for bilinear terms?",
         &propdata->createbilinineqs, TRUE, DEFAULT_CREATE_BILININEQS, NULL, NULL) );
 
  SCIP_CALL( SCIPaddBoolParam(scip, "propagating/" PROP_NAME "/createlincons",
         "create linear constraints from inequalities for bilinear terms?",
         &propdata->createlincons, TRUE, DEFAULT_CREATE_LINCONS, NULL, NULL) );
 
  SCIP_CALL( SCIPaddIntParam(scip, "propagating/" PROP_NAME "/orderingalgo",
        "select the type of ordering algorithm which should be used (0: no special ordering, 1: greedy, 2: greedy reverse)",
        &propdata->orderingalgo, TRUE, DEFAULT_ORDERINGALGO, 0, 2, NULL, NULL) );
 
  SCIP_CALL( SCIPaddBoolParam(scip, "propagating/" PROP_NAME "/separatesol",
         "should the obbt LP solution be separated?",
         &propdata->separatesol, TRUE, DEFAULT_SEPARATESOL, NULL, NULL) );
 
  SCIP_CALL( SCIPaddIntParam(scip, "propagating/" PROP_NAME "/sepaminiter",
        "minimum number of iteration spend to separate an obbt LP solution",
        &propdata->sepaminiter, TRUE, DEFAULT_SEPAMINITER, 0, INT_MAX, NULL, NULL) );
 
  SCIP_CALL( SCIPaddIntParam(scip, "propagating/" PROP_NAME "/sepamaxiter",
        "maximum number of iteration spend to separate an obbt LP solution",
        &propdata->sepamaxiter, TRUE, DEFAULT_SEPAMAXITER, 0, INT_MAX, NULL, NULL) );
 
  SCIP_CALL( SCIPaddIntParam(scip, "propagating/" PROP_NAME "/propagatefreq",
        "trigger a propagation round after that many bound tightenings (0: no propagation)",
        &propdata->propagatefreq, TRUE, DEFAULT_PROPAGATEFREQ, 0, INT_MAX, NULL, NULL) );
 
   return SCIP_OKAY;
}