cons_nonlinear.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-2022 Konrad-Zuse-Zentrum                            */
/*                            fuer Informationstechnik Berlin                */
/*                                                                           */
/*  SCIP is distributed under the terms of the ZIB Academic License.         */
/*                                                                           */
/*  You should have received a copy of the ZIB Academic License              */
/*  along with SCIP; see the file COPYING. If not visit scipopt.org.         */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/**@file   cons_nonlinear.c
 * @ingroup DEFPLUGINS_CONS
 * @brief  constraint handler for nonlinear constraints specified by algebraic expressions
 * @author Ksenia Bestuzheva
 * @author Benjamin Mueller
 * @author Felipe Serrano
 * @author Stefan Vigerske
 */

/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/

#ifdef SCIP_DEBUG
#define ENFO_LOGGING
#endif

/* enable to get log output for enforcement */
/* #define ENFO_LOGGING */
/* define to get enforcement logging into file */
/* #define ENFOLOGFILE "consexpr_enfo.log" */

/* define to get more debug output from domain propagation */
/* #define DEBUG_PROP */

/*lint -e440*/
/*lint -e441*/
/*lint -e528*/
/*lint -e666*/
/*lint -e777*/
/*lint -e866*/

#include <assert.h>
#include <ctype.h>

#include "scip/cons_nonlinear.h"
#include "scip/nlhdlr.h"
#include "scip/expr_var.h"
#include "scip/expr_sum.h"
#include "scip/expr_value.h"
#include "scip/expr_pow.h"
#include "scip/nlhdlr_convex.h"
#include "scip/cons_linear.h"
#include "scip/cons_varbound.h"
#include "scip/cons_and.h"
#include "scip/cons_bounddisjunction.h"
#include "scip/heur_subnlp.h"
#include "scip/heur_trysol.h"
#include "scip/nlpi_ipopt.h"  /* for SCIPsolveLinearEquationsIpopt */
#include "scip/debug.h"
#include "scip/dialog_default.h"

/* fundamental constraint handler properties */
#define CONSHDLR_NAME          "nonlinear"
#define CONSHDLR_DESC          "handler for nonlinear constraints specified by algebraic expressions"
#define CONSHDLR_ENFOPRIORITY       -60 /**< priority of the constraint handler for constraint enforcing */
#define CONSHDLR_CHECKPRIORITY -4000010 /**< priority of the constraint handler for checking feasibility */
#define CONSHDLR_EAGERFREQ          100 /**< frequency for using all instead of only the useful constraints in separation,
                                         *   propagation and enforcement, -1 for no eager evaluations, 0 for first only */
#define CONSHDLR_NEEDSCONS         TRUE /**< should the constraint handler be skipped, if no constraints are available? */

/* optional constraint handler properties */
#define CONSHDLR_SEPAPRIORITY        10 /**< priority of the constraint handler for separation */
#define CONSHDLR_SEPAFREQ             1 /**< frequency for separating cuts; zero means to separate only in the root node */
#define CONSHDLR_DELAYSEPA        FALSE /**< should separation method be delayed, if other separators found cuts? */

#define CONSHDLR_PROPFREQ             1 /**< frequency for propagating domains; zero means only preprocessing propagation */
#define CONSHDLR_DELAYPROP        FALSE /**< should propagation method be delayed, if other propagators found reductions? */
#define CONSHDLR_PROP_TIMING     SCIP_PROPTIMING_BEFORELP /**< propagation timing mask of the constraint handler*/

#define CONSHDLR_PRESOLTIMING    SCIP_PRESOLTIMING_ALWAYS /**< presolving timing of the constraint handler (fast, medium, or exhaustive) */
#define CONSHDLR_MAXPREROUNDS        -1 /**< maximal number of presolving rounds the constraint handler participates in (-1: no limit) */

/* properties of the nonlinear constraint handler statistics table */
#define TABLE_NAME_NONLINEAR           "cons_nonlinear"
#define TABLE_DESC_NONLINEAR           "nonlinear constraint handler statistics"
#define TABLE_POSITION_NONLINEAR       14600                  /**< the position of the statistics table */
#define TABLE_EARLIEST_STAGE_NONLINEAR SCIP_STAGE_TRANSFORMED /**< output of the statistics table is only printed from this stage onwards */

/* properties of the nonlinear handler statistics table */
#define TABLE_NAME_NLHDLR              "nlhdlr"
#define TABLE_DESC_NLHDLR              "nonlinear handler statistics"
#define TABLE_POSITION_NLHDLR          14601                  /**< the position of the statistics table */
#define TABLE_EARLIEST_STAGE_NLHDLR    SCIP_STAGE_PRESOLVING  /**< output of the statistics table is only printed from this stage onwards */

#define DIALOG_NAME            "nlhdlrs"
#define DIALOG_DESC            "display nonlinear handlers"
#define DIALOG_ISSUBMENU          FALSE

#define VERTEXPOLY_MAXPERTURBATION      1e-3 /**< maximum perturbation */
#define VERTEXPOLY_USEDUALSIMPLEX       TRUE /**< use dual or primal simplex algorithm? */
#define VERTEXPOLY_RANDNUMINITSEED  20181029 /**< seed for random number generator, which is used to move points away from the boundary */
#define VERTEXPOLY_ADJUSTFACETFACTOR     1e1 /**< adjust resulting facets in checkRikun() up to a violation of this value times lpfeastol */

#define BRANCH_RANDNUMINITSEED      20191229 /**< seed for random number generator, which is used to select from several similar good branching candidates */

#define BILIN_MAXNAUXEXPRS                10 /**< maximal number of auxiliary expressions per bilinear term */

/** translate from one value of infinity to another
 *
 *  if val is &ge; infty1, then give infty2, else give val
 */
#define infty2infty(infty1, infty2, val) ((val) >= (infty1) ? (infty2) : (val))

/** translates x to 2^x for non-negative integer x */
#define POWEROFTWO(x) (0x1u << (x))

#ifdef ENFO_LOGGING
#define ENFOLOG(x) if( SCIPgetSubscipDepth(scip) == 0 && SCIPgetVerbLevel(scip) >= SCIP_VERBLEVEL_NORMAL ) { x }
FILE* enfologfile = NULL;
#else
#define ENFOLOG(x)
#endif

/*
 * Data structures
 */

/** enforcement data of an expression */
typedef struct
{
   SCIP_NLHDLR*          nlhdlr;             /**< nonlinear handler */
   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata;     /**< data of nonlinear handler */
   SCIP_NLHDLR_METHOD nlhdlrparticipation; /**< methods where nonlinear handler participates */
   SCIP_Bool             issepainit;         /**< was the initsepa callback of nlhdlr called */
   SCIP_Real             auxvalue;           /**< auxiliary value of expression w.r.t. currently enforced solution */
   SCIP_Bool             sepabelowusesactivity;/**< whether sepabelow uses activity of some expression */
   SCIP_Bool             sepaaboveusesactivity;/**< whether sepaabove uses activity of some expression */
} EXPRENFO;

/** data stored by constraint handler in an expression that belongs to a nonlinear constraint */
struct SCIP_Expr_OwnerData
{
   SCIP_CONSHDLR*        conshdlr;           /** nonlinear constraint handler */

   /* locks and monotonicity */
   int                   nlockspos;          /**< positive locks counter */
   int                   nlocksneg;          /**< negative locks counter */
   SCIP_MONOTONE*        monotonicity;       /**< array containing monotonicity of expression w.r.t. each child */
   int                   monotonicitysize;   /**< length of monotonicity array */

   /* propagation (in addition to activity that is stored in expr) */
   SCIP_INTERVAL         propbounds;         /**< bounds to propagate in reverse propagation */
   unsigned int          propboundstag;      /**< tag to indicate whether propbounds are valid for the current propagation rounds */
   SCIP_Bool             inpropqueue;        /**< whether expression is queued for propagation */

   /* enforcement of expr == auxvar (or expr <= auxvar, or expr >= auxvar) */
   EXPRENFO**            enfos;              /**< enforcements */
   int                   nenfos;             /**< number of enforcements, or -1 if not initialized */
   unsigned int          lastenforced;       /**< last enforcement round where expression was enforced successfully */
   unsigned int          nactivityusesprop;  /**< number of nonlinear handlers whose activity computation (or domain propagation) depends on the activity of the expression */
   unsigned int          nactivityusessepa;  /**< number of nonlinear handlers whose separation (estimate or enfo) depends on the activity of the expression */
   unsigned int          nauxvaruses;        /**< number of nonlinear handlers whose separation uses an auxvar in the expression */
   SCIP_VAR*             auxvar;             /**< auxiliary variable used for outer approximation cuts */

   /* branching */
   SCIP_Real             violscoresum;       /**< sum of violation scores for branching stored for this expression */
   SCIP_Real             violscoremax;       /**< max of violation scores for branching stored for this expression */
   int                   nviolscores;        /**< number of violation scores stored for this expression */
   unsigned int          violscoretag;       /**< tag to decide whether a violation score of an expression needs to be initialized */

   /* additional data for variable expressions (TODO move into sub-struct?) */
   SCIP_CONS**           conss;              /**< constraints in which this variable appears */
   int                   nconss;             /**< current number of constraints in conss */
   int                   consssize;          /**< length of conss array */
   SCIP_Bool             consssorted;        /**< is the array of constraints sorted */

   int                   filterpos;          /**< position of eventdata in SCIP's event filter, -1 if not catching events */
};

/** constraint data for nonlinear constraints */
struct SCIP_ConsData
{
   /* data that defines the constraint: expression and sides */
   SCIP_EXPR*            expr;               /**< expression that represents this constraint */
   SCIP_Real             lhs;                /**< left-hand side */
   SCIP_Real             rhs;                /**< right-hand side */

   /* variables */
   SCIP_EXPR**           varexprs;           /**< array containing all variable expressions */
   int                   nvarexprs;          /**< total number of variable expressions */
   SCIP_Bool             catchedevents;      /**< do we catch events on variables? */

   /* constraint violation */
   SCIP_Real             lhsviol;            /**< violation of left-hand side by current solution */
   SCIP_Real             rhsviol;            /**< violation of right-hand side by current solution */
   SCIP_Real             gradnorm;           /**< norm of gradient of constraint function in current solution (if evaluated) */
   SCIP_Longint          gradnormsoltag;     /**< tag of solution used that gradnorm corresponds to */

   /* status flags */
   unsigned int          ispropagated:1;     /**< did we propagate the current bounds already? */
   unsigned int          issimplified:1;     /**< did we simplify the expression tree already? */

   /* locks */
   int                   nlockspos;          /**< number of positive locks */
   int                   nlocksneg;          /**< number of negative locks */

   /* repair infeasible solutions */
   SCIP_VAR*             linvardecr;         /**< variable that may be decreased without making any other constraint infeasible, or NULL if none */
   SCIP_VAR*             linvarincr;         /**< variable that may be increased without making any other constraint infeasible, or NULL if none */
   SCIP_Real             linvardecrcoef;     /**< linear coefficient of linvardecr */
   SCIP_Real             linvarincrcoef;     /**< linear coefficient of linvarincr */

   /* miscellaneous */
   SCIP_EXPRCURV         curv;               /**< curvature of the root expression w.r.t. the original variables */
   SCIP_NLROW*           nlrow;              /**< a nonlinear row representation of this constraint */
   int                   consindex;          /**< an index of the constraint that is unique among all expr-constraints in this SCIP instance and is constant */
};

/** constraint upgrade method */
typedef struct
{
   SCIP_DECL_NONLINCONSUPGD((*consupgd));    /**< method to call for upgrading nonlinear constraint */
   int                   priority;           /**< priority of upgrading method */
   SCIP_Bool             active;             /**< is upgrading enabled */
} CONSUPGRADE;

/** constraint handler data */
struct SCIP_ConshdlrData
{
   /* nonlinear handler */
   SCIP_NLHDLR**         nlhdlrs;            /**< nonlinear handlers */
   int                   nnlhdlrs;           /**< number of nonlinear handlers */
   int                   nlhdlrssize;        /**< size of nlhdlrs array */
   SCIP_Bool             indetect;           /**< whether we are currently in detectNlhdlr */
   SCIP_Bool             registerusesactivitysepabelow; /**< a flag that is used only during \ref @detectNlhdlr() */
   SCIP_Bool             registerusesactivitysepaabove; /**< a flag that is used only during \ref @detectNlhdlr() */

   /* constraint upgrades */
   CONSUPGRADE**         consupgrades;       /**< constraint upgrade methods for specializing nonlinear constraints */
   int                   consupgradessize;   /**< size of consupgrades array */
   int                   nconsupgrades;      /**< number of constraint upgrade methods */

   /* other plugins */
   SCIP_EVENTHDLR*       eventhdlr;          /**< handler for variable bound change events */
   SCIP_HEUR*            subnlpheur;         /**< a pointer to the subnlp heuristic, if available */
   SCIP_HEUR*            trysolheur;         /**< a pointer to the trysol heuristic, if available */

   /* tags and counters */
   int                   auxvarid;           /**< unique id for the next auxiliary variable */
   SCIP_Longint          curboundstag;       /**< tag indicating current variable bounds */
   SCIP_Longint          lastboundrelax;     /**< tag when bounds where most recently relaxed */
   SCIP_Longint          lastvaractivitymethodchange; /**< tag when method used to evaluate activity of variables changed last */
   unsigned int          enforound;          /**< total number of enforcement calls, including current one */
   int                   lastconsindex;      /**< last used consindex, plus one */

   /* activity intervals and domain propagation */
   SCIP_DECL_EXPR_INTEVALVAR((*intevalvar)); /**< method currently used for activity calculation of variable expressions */
   SCIP_Bool             globalbounds;       /**< whether global variable bounds should be used for activity calculation */
   SCIP_QUEUE*           reversepropqueue;   /**< expression queue to be used in reverse propagation, filled by SCIPtightenExprIntervalNonlinear */
   SCIP_Bool             forceboundtightening; /**< whether bound change passed to SCIPtightenExprIntervalNonlinear should be forced */
   unsigned int          curpropboundstag;   /**< tag indicating current propagation rounds, to match with expr->propboundstag */

   /* parameters */
   int                   maxproprounds;      /**< limit on number of propagation rounds for a set of constraints within one round of SCIP propagation */
   SCIP_Bool             propauxvars;        /**< whether to check bounds of all auxiliary variable to seed reverse propagation */
   char                  varboundrelax;      /**< strategy on how to relax variable bounds during bound tightening */
   SCIP_Real             varboundrelaxamount; /**< by how much to relax variable bounds during bound tightening */
   SCIP_Real             conssiderelaxamount; /**< by how much to relax constraint sides during bound tightening */
   SCIP_Real             vp_maxperturb;      /**< maximal relative perturbation of reference point */
   SCIP_Real             vp_adjfacetthreshold; /**< adjust computed facet up to a violation of this value times lpfeastol */
   SCIP_Bool             vp_dualsimplex;     /**< whether to use dual simplex instead of primal simplex for facet computing LP */
   SCIP_Bool             reformbinprods;     /**< whether to reformulate products of binary variables during presolving */
   SCIP_Bool             reformbinprodsand;  /**< whether to use the AND constraint handler for reformulating binary products */
   int                   reformbinprodsfac;  /**< minimum number of terms to reformulate bilinear binary products by factorizing variables (<= 1: disabled) */
   SCIP_Bool             forbidmultaggrnlvar; /**< whether to forbid multiaggregation of variables that appear in a nonlinear term of a constraint */
   SCIP_Bool             tightenlpfeastol;   /**< whether to tighten LP feasibility tolerance during enforcement, if it seems useful */
   SCIP_Bool             propinenforce;      /**< whether to (re)run propagation in enforcement */
   SCIP_Real             weakcutthreshold;   /**< threshold for when to regard a cut from an estimator as weak */
   SCIP_Real             strongcutmaxcoef;   /**< "strong" cuts will be scaled to have their maximal coef in [1/strongcutmaxcoef,strongcutmaxcoef] */
   SCIP_Bool             strongcutefficacy;  /**< consider efficacy requirement when deciding whether a cut is "strong" */
   SCIP_Bool             forcestrongcut;     /**< whether to force "strong" cuts in enforcement */
   SCIP_Real             enfoauxviolfactor;  /**< an expression will be enforced if the "auxiliary" violation is at least enfoauxviolfactor times the "original" violation */
   SCIP_Real             weakcutminviolfactor; /**< retry with weak cuts for constraints with violation at least this factor of maximal violated constraints */
   char                  rownotremovable;    /**< whether to make rows to be non-removable in the node where they are added (can prevent some cycling): 'o'ff, in 'e'nforcement only, 'a'lways */
   char                  violscale;          /**< method how to scale violations to make them comparable (not used for feasibility check) */
   char                  checkvarlocks;      /**< whether variables contained in a single constraint should be forced to be at their lower or upper bounds ('d'isable, change 't'ype, add 'b'ound disjunction) */
   int                   branchauxmindepth;  /**< from which depth on to allow branching on auxiliary variables */
   SCIP_Bool             branchexternal;     /**< whether to use external branching candidates for branching */
   SCIP_Real             branchhighviolfactor; /**< consider a constraint highly violated if its violation is >= this factor * maximal violation among all constraints */
   SCIP_Real             branchhighscorefactor; /**< consider a variable branching score high if its branching score >= this factor * maximal branching score among all variables */
   SCIP_Real             branchviolweight;   /**< weight by how much to consider the violation assigned to a variable for its branching score */
   SCIP_Real             branchdualweight;   /**< weight by how much to consider the dual values of rows that contain a variable for its branching score */
   SCIP_Real             branchpscostweight; /**< weight by how much to consider the pseudo cost of a variable for its branching score */
   SCIP_Real             branchdomainweight; /**< weight by how much to consider the domain width in branching score */
   SCIP_Real             branchvartypeweight;/**< weight by how much to consider variable type in branching score */
   char                  branchscoreagg;     /**< how to aggregate several branching scores given for the same expression ('a'verage, 'm'aximum, or 's'um) */
   char                  branchviolsplit;    /**< method used to split violation in expression onto variables ('u'niform, 'm'idness of solution, 'd'omain width, 'l'ogarithmic domain width) */
   SCIP_Real             branchpscostreliable; /**< minimum pseudo-cost update count required to consider pseudo-costs reliable */
   char                  linearizeheursol;   /**< whether tight linearizations of nonlinear constraints should be added to cutpool when some heuristics finds a new solution ('o'ff, on new 'i'ncumbents, on 'e'very solution) */
   SCIP_Bool             assumeconvex;       /**< whether to assume that any constraint is convex */

   /* statistics */
   SCIP_Longint          nweaksepa;          /**< number of times we used "weak" cuts for enforcement */
   SCIP_Longint          ntightenlp;         /**< number of times we requested solving the LP with a smaller feasibility tolerance when enforcing */
   SCIP_Longint          ndesperatetightenlp; /**< number of times we requested solving the LP with a smaller feasibility tolerance when enforcing because we didn't know anything better */
   SCIP_Longint          ndesperatebranch;   /**< number of times we branched on some variable because normal enforcement was not successful */
   SCIP_Longint          ndesperatecutoff;   /**< number of times we cut off a node in enforcement because no branching candidate could be found */
   SCIP_Longint          nforcelp;           /**< number of times we forced solving the LP when enforcing a pseudo solution */
   SCIP_CLOCK*           canonicalizetime;   /**< time spend for canonicalization */
   SCIP_Longint          ncanonicalizecalls; /**< number of times we called canonicalization */

   /* facets of envelops of vertex-polyhedral functions */
   SCIP_RANDNUMGEN*      vp_randnumgen;      /**< random number generator used to perturb reference point */
   SCIP_LPI*             vp_lp[SCIP_MAXVERTEXPOLYDIM+1];  /**< LPs used to compute facets for functions of different dimension */

   /* hashing of bilinear terms */
   SCIP_HASHTABLE*       bilinhashtable;     /**< hash table for bilinear terms */
   SCIP_CONSNONLINEAR_BILINTERM* bilinterms; /**< bilinear terms */
   int                   nbilinterms;        /**< total number of bilinear terms */
   int                   bilintermssize;     /**< size of bilinterms array */
   int                   bilinmaxnauxexprs;  /**< maximal number of auxiliary expressions per bilinear term */

   /* branching */
   SCIP_RANDNUMGEN*      branchrandnumgen;   /**< random number generated used in branching variable selection */
   char                  branchpscostupdatestrategy; /**< value of parameter branching/lpgainnormalize */

   /* misc */
   SCIP_Bool             checkedvarlocks;    /**< whether variables contained in a single constraint have been already considered */
   SCIP_HASHMAP*         var2expr;           /**< hashmap to map SCIP variables to variable-expressions */
   int                   newsoleventfilterpos; /**< filter position of new solution event handler, if caught */
};

/** branching candidate with various scores */
typedef struct
{
   SCIP_EXPR*            expr;               /**< expression that holds branching candidate */
   SCIP_Real             auxviol;            /**< aux-violation score of candidate */
   SCIP_Real             domain;             /**< domain score of candidate */
   SCIP_Real             dual;               /**< dual score of candidate */
   SCIP_Real             pscost;             /**< pseudo-cost score of candidate */
   SCIP_Real             vartype;            /**< variable type score of candidate */
   SCIP_Real             weighted;           /**< weighted sum of other scores, see scoreBranchingCandidates() */
} BRANCHCAND;

/*
 * Local methods
 */

/* forward declaration */
static
SCIP_RETCODE forwardPropExpr(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_EXPR*            rootexpr,           /**< expression */
   SCIP_Bool             tightenauxvars,     /**< should the bounds of auxiliary variables be tightened? */
   SCIP_Bool*            infeasible,         /**< buffer to store whether the problem is infeasible (NULL if not needed) */
   int*                  ntightenings        /**< buffer to store the number of auxiliary variable tightenings (NULL if not needed) */
   );

/** frees auxiliary variables of expression, if any */
static
SCIP_RETCODE freeAuxVar(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EXPR*            expr                /**< expression which auxvar to free, if any */
   )
{
   SCIP_EXPR_OWNERDATA* mydata;

   assert(scip != NULL);
   assert(expr != NULL);

   mydata = SCIPexprGetOwnerData(expr);
   assert(mydata != NULL);

   if( mydata->auxvar == NULL )
      return SCIP_OKAY;

   SCIPdebugMsg(scip, "remove auxiliary variable <%s> for expression %p\n", SCIPvarGetName(mydata->auxvar), (void*)expr);

   /* remove variable locks
    * as this is a relaxation-only variable, no other plugin should use it for deducing any type of reductions or cutting planes
    */
   SCIP_CALL( SCIPaddVarLocks(scip, mydata->auxvar, -1, -1) );

   /* release auxiliary variable */
   SCIP_CALL( SCIPreleaseVar(scip, &mydata->auxvar) );
   assert(mydata->auxvar == NULL);

   return SCIP_OKAY;
}

/** frees data used for enforcement of expression, that is, nonlinear handlers
 *
 * can also clear indicators whether expr needs enforcement methods, that is,
 * free an associated auxiliary variable and reset the nactivityuses counts
 */
static
SCIP_RETCODE freeEnfoData(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EXPR*            expr,               /**< expression whose enforcement data will be released */
   SCIP_Bool             freeauxvar          /**< whether aux var should be released and activity usage counts be reset */
   )
{
   SCIP_EXPR_OWNERDATA* mydata;
   int e;

   mydata = SCIPexprGetOwnerData(expr);
   assert(mydata != NULL);

   if( freeauxvar )
   {
      /* free auxiliary variable */
      SCIP_CALL( freeAuxVar(scip, expr) );
      assert(mydata->auxvar == NULL);

      /* reset count on activity and auxvar usage */
      mydata->nactivityusesprop = 0;
      mydata->nactivityusessepa = 0;
      mydata->nauxvaruses = 0;
   }

   /* free data stored by nonlinear handlers */
   for( e = 0; e < mydata->nenfos; ++e )
   {
      SCIP_NLHDLR* nlhdlr;

      assert(mydata->enfos[e] != NULL);

      nlhdlr = mydata->enfos[e]->nlhdlr;
      assert(nlhdlr != NULL);

      if( mydata->enfos[e]->issepainit )
      {
         /* call the separation deinitialization callback of the nonlinear handler */
         SCIP_CALL( SCIPnlhdlrExitsepa(scip, nlhdlr, expr, mydata->enfos[e]->nlhdlrexprdata) );
         mydata->enfos[e]->issepainit = FALSE;
      }

      /* free nlhdlr exprdata, if there is any and there is a method to free this data */
      if( mydata->enfos[e]->nlhdlrexprdata != NULL )
      {
         SCIP_CALL( SCIPnlhdlrFreeexprdata(scip, nlhdlr, expr, &mydata->enfos[e]->nlhdlrexprdata) );
         assert(mydata->enfos[e]->nlhdlrexprdata == NULL);
      }

      /* free enfo data */
      SCIPfreeBlockMemory(scip, &mydata->enfos[e]);
   }

   /* free array with enfo data */
   SCIPfreeBlockMemoryArrayNull(scip, &mydata->enfos, mydata->nenfos);

   /* we need to look at this expression in detect again */
   mydata->nenfos = -1;

   return SCIP_OKAY;
}

/** callback that frees data that this conshdlr stored in an expression */
static
SCIP_DECL_EXPR_OWNERFREE(exprownerFree)
{
   assert(scip != NULL);
   assert(expr != NULL);
   assert(ownerdata != NULL);
   assert(*ownerdata != NULL);

   /* expression should not be locked anymore */
   assert((*ownerdata)->nlockspos == 0);
   assert((*ownerdata)->nlocksneg == 0);

   SCIP_CALL( freeEnfoData(scip, expr, TRUE) );

   /* expression should not be enforced anymore */
   assert((*ownerdata)->nenfos <= 0);
   assert((*ownerdata)->auxvar == NULL);

   if( SCIPisExprVar(scip, expr) )
   {
      SCIP_CONSHDLRDATA* conshdlrdata;
      SCIP_VAR* var;

      /* there should be no constraints left that still use this variable */
      assert((*ownerdata)->nconss == 0);
      /* thus, there should also be no variable event catched (via this exprhdlr) */
      assert((*ownerdata)->filterpos == -1);

      SCIPfreeBlockMemoryArrayNull(scip, &(*ownerdata)->conss, (*ownerdata)->consssize);

      /* update var2expr hashmap in conshdlrdata */
      conshdlrdata = SCIPconshdlrGetData((*ownerdata)->conshdlr);
      assert(conshdlrdata != NULL);

      var = SCIPgetVarExprVar(expr);
      assert(var != NULL);

      /* remove var -> expr map from hashmap if present
       *  (if no variable-expression stored for var hashmap, then the var hasn't been used in any constraint, so do nothing
       *   if variable-expression stored for var is different, then also do nothing)
       */
      if( SCIPhashmapGetImage(conshdlrdata->var2expr, var) == (void*)expr )
      {
         SCIP_CALL( SCIPhashmapRemove(conshdlrdata->var2expr, var) );
      }
   }

   SCIPfreeBlockMemory(scip, ownerdata);

   return SCIP_OKAY;
}

static
SCIP_DECL_EXPR_OWNERPRINT(exprownerPrint)
{  /*lint --e{715}*/
   assert(ownerdata != NULL);

   /* print nl handlers associated to expr */
   if( ownerdata->nenfos > 0 )
   {
      int i;
      SCIPinfoMessage(scip, file, "   {");

      for( i = 0; i < ownerdata->nenfos; ++i )
      {
         SCIPinfoMessage(scip, file, "%s:", SCIPnlhdlrGetName(ownerdata->enfos[i]->nlhdlr));
         if( ownerdata->enfos[i]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_ACTIVITY )
            SCIPinfoMessage(scip, file, "a");
         if( ownerdata->enfos[i]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_SEPABELOW )
            SCIPinfoMessage(scip, file, "u");
         if( ownerdata->enfos[i]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_SEPAABOVE )
            SCIPinfoMessage(scip, file, "o");
         if( i < ownerdata->nenfos-1 )
            SCIPinfoMessage(scip, file, ", ");
      }

      SCIPinfoMessage(scip, file, "}");
   }

   /* print aux var associated to expr */
   if( ownerdata->auxvar != NULL )
   {
      SCIPinfoMessage(scip, file, "  (<%s> in [%g, %g])", SCIPvarGetName(ownerdata->auxvar), SCIPvarGetLbLocal(ownerdata->auxvar), SCIPvarGetUbLocal(ownerdata->auxvar));
   }
   SCIPinfoMessage(scip, file, "\n");

   return SCIP_OKAY;
}

/** possibly reevaluates and then returns the activity of the expression
 *
 * Reevaluate activity if currently stored is not up to date (some bound was changed since last evaluation).
 */
static
SCIP_DECL_EXPR_OWNEREVALACTIVITY(exprownerEvalactivity)
{
   SCIP_CONSHDLRDATA* conshdlrdata;

   assert(scip != NULL);
   assert(expr != NULL);
   assert(ownerdata != NULL);

   conshdlrdata = SCIPconshdlrGetData(ownerdata->conshdlr);
   assert(conshdlrdata != NULL);

   if( SCIPexprGetActivityTag(expr) < conshdlrdata->curboundstag )
   {
      /* update activity of expression */
      SCIP_CALL( forwardPropExpr(scip, ownerdata->conshdlr, expr, FALSE, NULL, NULL) );

      assert(SCIPexprGetActivityTag(expr) == conshdlrdata->curboundstag);
   }

   return SCIP_OKAY;
}

/** callback that creates data that this conshdlr wants to store in an expression */
static
SCIP_DECL_EXPR_OWNERCREATE(exprownerCreate)
{
   assert(scip != NULL);
   assert(expr != NULL);
   assert(ownerdata != NULL);

   SCIP_CALL( SCIPallocClearBlockMemory(scip, ownerdata) );
   (*ownerdata)->nenfos = -1;
   (*ownerdata)->conshdlr = (SCIP_CONSHDLR*)ownercreatedata;

   if( SCIPisExprVar(scip, expr) )
   {
      SCIP_CONSHDLRDATA* conshdlrdata;
      SCIP_VAR* var;

      (*ownerdata)->filterpos = -1;

      /* add to var2expr hashmap if not having expr for var yet */

      conshdlrdata = SCIPconshdlrGetData((*ownerdata)->conshdlr);
      assert(conshdlrdata != NULL);

      var = SCIPgetVarExprVar(expr);

      if( !SCIPhashmapExists(conshdlrdata->var2expr, (void*)var) )
      {
         /* store the variable expression in the hashmap */
         SCIP_CALL( SCIPhashmapInsert(conshdlrdata->var2expr, (void*)var, (void*)expr) );
      }
      else
      {
         /* if expr was just created, then it shouldn't already be stored as image of var */
         assert(SCIPhashmapGetImage(conshdlrdata->var2expr, (void*)var) != (void*)expr);
      }
   }
   else
   {
      /* just so that we can use filterpos to recognize whether an expr is a varexpr if not having a SCIP pointer around */
      (*ownerdata)->filterpos = -2;
   }

   *ownerfree = exprownerFree;
   *ownerprint = exprownerPrint;
   *ownerevalactivity = exprownerEvalactivity;

   return SCIP_OKAY;
}

/** creates a variable expression or retrieves from hashmap in conshdlr data */
static
SCIP_RETCODE createExprVar(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< nonlinear constraint handler */
   SCIP_EXPR**           expr,               /**< pointer where to store expression */
   SCIP_VAR*             var                 /**< variable to be stored */
   )
{
   assert(conshdlr != NULL);
   assert(expr != NULL);
   assert(var != NULL);

   /* get variable expression representing the given variable if there is one already */
   *expr = (SCIP_EXPR*) SCIPhashmapGetImage(SCIPconshdlrGetData(conshdlr)->var2expr, (void*) var);

   if( *expr == NULL )
   {
      /* create a new variable expression; this also captures the expression */
      SCIP_CALL( SCIPcreateExprVar(scip, expr, var, exprownerCreate, (void*)conshdlr) );
      assert(*expr != NULL);
      /* exprownerCreate should have added var->expr to var2expr */
      assert(SCIPhashmapGetImage(SCIPconshdlrGetData(conshdlr)->var2expr, (void*)var) == (void*)*expr);
   }
   else
   {
      /* only capture already existing expr to get a consistent uses-count */
      SCIPcaptureExpr(*expr);
   }

   return SCIP_OKAY;
}

/* map var exprs to var-expr from var2expr hashmap */
static
SCIP_DECL_EXPR_MAPEXPR(mapexprvar)
{  /*lint --e{715}*/
   SCIP_CONSHDLR* conshdlr = (SCIP_CONSHDLR*)mapexprdata;

   assert(sourcescip != NULL);
   assert(targetscip != NULL);
   assert(sourceexpr != NULL);
   assert(targetexpr != NULL);
   assert(*targetexpr == NULL);
   assert(mapexprdata != NULL);

   /* do not provide map if not variable */
   if( !SCIPisExprVar(sourcescip, sourceexpr) )
      return SCIP_OKAY;

   SCIP_CALL( createExprVar(targetscip, conshdlr, targetexpr, SCIPgetVarExprVar(sourceexpr)) );

   return SCIP_OKAY;
}

/* map var exprs to var-expr from var2expr hashmap corresponding to transformed var */
static
SCIP_DECL_EXPR_MAPEXPR(mapexprtransvar)
{  /*lint --e{715}*/
   SCIP_CONSHDLR* conshdlr = (SCIP_CONSHDLR*)mapexprdata;
   SCIP_VAR* var;

   assert(sourcescip != NULL);
   assert(targetscip != NULL);
   assert(sourceexpr != NULL);
   assert(targetexpr != NULL);
   assert(*targetexpr == NULL);
   assert(mapexprdata != NULL);

   /* do not provide map if not variable */
   if( !SCIPisExprVar(sourcescip, sourceexpr) )
      return SCIP_OKAY;

   var = SCIPgetVarExprVar(sourceexpr);
   assert(var != NULL);

   /* transform variable */
   SCIP_CALL( SCIPgetTransformedVar(sourcescip, var, &var) );
   assert(var != NULL);

   SCIP_CALL( createExprVar(targetscip, conshdlr, targetexpr, var) );

   return SCIP_OKAY;
}

/** stores all variable expressions into a given constraint */
static
SCIP_RETCODE storeVarExprs(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONSDATA*        consdata            /**< constraint data */
   )
{
   SCIP_CONSHDLRDATA* conshdlrdata;
   int varexprssize;
   int i;

   assert(consdata != NULL);

   /* skip if we have stored the variable expressions already */
   if( consdata->varexprs != NULL )
      return SCIP_OKAY;

   assert(consdata->varexprs == NULL);
   assert(consdata->nvarexprs == 0);

   /* get an upper bound on number of variable expressions */
   if( consdata->issimplified )
   {
      /* if simplified, then we should have removed inactive variables and replaced common subexpressions,
       * so we cannot have more variable expression than the number of active variables
       */
      varexprssize = SCIPgetNVars(scip);
   }
   else
   {
      SCIP_CALL( SCIPgetExprNVars(scip, consdata->expr, &varexprssize) );
   }

   /* create array to store all variable expressions */
   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->varexprs, varexprssize) );

   SCIP_CALL( SCIPgetExprVarExprs(scip, consdata->expr, consdata->varexprs, &(consdata->nvarexprs)) );
   assert(varexprssize >= consdata->nvarexprs);

   /* shrink array if there are less variables in the expression than in the problem */
   if( varexprssize > consdata->nvarexprs )
   {
      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->varexprs, varexprssize, consdata->nvarexprs) );
   }

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);
   assert(conshdlrdata->var2expr != NULL);

   /* ensure that for every variable an entry exists in the var2expr hashmap
    * when removing duplicate subexpressions it can happen that a var->varexpr map was removed from the hashmap
    */
   for( i = 0; i < consdata->nvarexprs; ++i )
   {
      if( !SCIPhashmapExists(conshdlrdata->var2expr, SCIPgetVarExprVar(consdata->varexprs[i])) )
      {
         SCIP_CALL( SCIPhashmapInsert(conshdlrdata->var2expr, SCIPgetVarExprVar(consdata->varexprs[i]), consdata->varexprs[i]) );
      }
   }

   return SCIP_OKAY;
}

/** frees all variable expression stored in storeVarExprs() */
static
SCIP_RETCODE freeVarExprs(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSDATA*        consdata            /**< constraint data */
   )
{
   int i;

   assert(consdata != NULL);

   /* skip if we have stored the variable expressions already*/
   if( consdata->varexprs == NULL )
      return SCIP_OKAY;

   assert(consdata->varexprs != NULL);
   assert(consdata->nvarexprs >= 0);

   /* release variable expressions */
   for( i = 0; i < consdata->nvarexprs; ++i )
   {
      assert(consdata->varexprs[i] != NULL);
      SCIP_CALL( SCIPreleaseExpr(scip, &consdata->varexprs[i]) );
      assert(consdata->varexprs[i] == NULL);
   }

   /* free variable expressions */
   SCIPfreeBlockMemoryArrayNull(scip, &consdata->varexprs, consdata->nvarexprs);
   consdata->varexprs = NULL;
   consdata->nvarexprs = 0;

   return SCIP_OKAY;
}

/** interval evaluation of variables as used in bound tightening
 *
 * Returns slightly relaxed local variable bounds of a variable as interval.
 * Does not relax beyond integer values, thus does not relax bounds on integer variables at all.
 */
static
SCIP_DECL_EXPR_INTEVALVAR(intEvalVarBoundTightening)
{
   SCIP_INTERVAL interval;
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_Real lb;
   SCIP_Real ub;

   assert(scip != NULL);
   assert(var != NULL);

   conshdlrdata = (SCIP_CONSHDLRDATA*)intevalvardata;
   assert(conshdlrdata != NULL);

   if( conshdlrdata->globalbounds )
   {
      lb = SCIPvarGetLbGlobal(var);
      ub = SCIPvarGetUbGlobal(var);
   }
   else
   {
      lb = SCIPvarGetLbLocal(var);
      ub = SCIPvarGetUbLocal(var);
   }
   assert(lb <= ub);  /* SCIP should ensure that variable bounds are not contradicting */

   /* implicit integer variables may have non-integer bounds, apparently (run space25a) */
   if( SCIPvarGetType(var) == SCIP_VARTYPE_IMPLINT )
   {
      lb = EPSROUND(lb, 0.0); /*lint !e835*/
      ub = EPSROUND(ub, 0.0); /*lint !e835*/
   }

   /* integer variables should always have integral bounds in SCIP */
   assert(EPSFRAC(lb, 0.0) == 0.0 || !SCIPvarIsIntegral(var));  /*lint !e835*/
   assert(EPSFRAC(ub, 0.0) == 0.0 || !SCIPvarIsIntegral(var));  /*lint !e835*/

   switch( conshdlrdata->varboundrelax )
   {
      case 'n' : /* no relaxation */
         break;

      case 'a' : /* relax by absolute value */
      {
         /* do not look at integer variables, they already have integral bounds, so wouldn't be relaxed */
         if( SCIPvarIsIntegral(var) )
            break;

         if( !SCIPisInfinity(scip, -lb) )
         {
            /* reduce lb by epsilon, or to the next integer value, which ever is larger */
            SCIP_Real bnd = floor(lb);
            lb = MAX(bnd, lb - conshdlrdata->varboundrelaxamount);
         }

         if( !SCIPisInfinity(scip, ub) )
         {
            /* increase ub by epsilon, or to the next integer value, which ever is smaller */
            SCIP_Real bnd = ceil(ub);
            ub = MIN(bnd, ub + conshdlrdata->varboundrelaxamount);
         }

         break;
      }

      case 'b' : /* relax always by absolute value */
      {
         /* do not look at integer variables, they already have integral bounds, so wouldn't be relaxed */
         if( SCIPvarIsIntegral(var) )
            break;

         if( !SCIPisInfinity(scip, -lb) )
            lb -= conshdlrdata->varboundrelaxamount;

         if( !SCIPisInfinity(scip, ub) )
            ub += conshdlrdata->varboundrelaxamount;

         break;
      }

      case 'r' : /* relax by relative value */
      {
         /* do not look at integer variables, they already have integral bounds, so wouldn't be relaxed */
         if( SCIPvarIsIntegral(var) )
            break;

         /* relax bounds by epsilon*max(1,|bnd|), instead of just epsilon as in case 'a', thus we trust the first log(epsilon) digits
          * however, when domains get small, relaxing can excessively weaken bound tightening, thus do only fraction of |ub-lb| if that is smaller
          * further, do not relax beyond next integer value
          */
         if( !SCIPisInfinity(scip, -lb) )
         {
            SCIP_Real bnd = floor(lb);
            lb = MAX(bnd, lb - MIN(conshdlrdata->varboundrelaxamount * MAX(1.0, REALABS(lb)), 0.001 * REALABS(ub-lb)));
         }

         if( !SCIPisInfinity(scip, ub) )
         {
            SCIP_Real bnd = ceil(ub);
            ub = MIN(bnd, ub + MIN(conshdlrdata->varboundrelaxamount * MAX(1.0, REALABS(ub)), 0.001 * REALABS(ub-lb)));
         }

         break;
      }

      default :
      {
         SCIPerrorMessage("Unsupported value '%c' for varboundrelax option.\n", conshdlrdata->varboundrelax);
         SCIPABORT();
         break;
      }
   }

   /* convert SCIPinfinity() to SCIP_INTERVAL_INFINITY */
   lb = -infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, -lb);
   ub =  infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, ub);
   assert(lb <= ub);

   SCIPintervalSetBounds(&interval, lb, ub);

   return interval;
}

/** compares two nonlinear constraints by its index
 *
 * Usable as compare operator in array sort functions.
 */
static
SCIP_DECL_SORTPTRCOMP(compIndexConsNonlinear)
{
   SCIP_CONSDATA* consdata1 = SCIPconsGetData((SCIP_CONS*)elem1);
   SCIP_CONSDATA* consdata2 = SCIPconsGetData((SCIP_CONS*)elem2);

   assert(consdata1 != NULL);
   assert(consdata2 != NULL);

   return consdata1->consindex - consdata2->consindex;
}

/** processes variable fixing or bound change event */
static
SCIP_DECL_EVENTEXEC(processVarEvent)
{  /*lint --e{715}*/
   SCIP_EVENTTYPE eventtype;
   SCIP_EXPR* expr;
   SCIP_EXPR_OWNERDATA* ownerdata;

   eventtype = SCIPeventGetType(event);
   assert(eventtype & (SCIP_EVENTTYPE_BOUNDCHANGED | SCIP_EVENTTYPE_VARFIXED));

   assert(eventdata != NULL);
   expr = (SCIP_EXPR*) eventdata;
   assert(SCIPisExprVar(scip, expr));

   SCIPdebugMsg(scip, "  exec event %" SCIP_EVENTTYPE_FORMAT " for variable <%s> (local [%g,%g], global [%g,%g])\n", eventtype,
         SCIPvarGetName(SCIPeventGetVar(event)),
         SCIPvarGetLbLocal(SCIPeventGetVar(event)), SCIPvarGetUbLocal(SCIPeventGetVar(event)),
         SCIPvarGetLbGlobal(SCIPeventGetVar(event)), SCIPvarGetUbGlobal(SCIPeventGetVar(event)));

   ownerdata = SCIPexprGetOwnerData(expr);
   assert(ownerdata != NULL);
   /* we only catch varevents for variables in constraints, so there should be constraints */
   assert(ownerdata->nconss > 0);
   assert(ownerdata->conss != NULL);

   /* notify constraints that use this variable expression (expr) to repropagate and possibly resimplify
    * - propagation can only find something new if a bound was tightened
    * - simplify can only find something new if a var is fixed (or maybe a bound is tightened)
    *   and we look at global changes (that is, we are not looking at boundchanges in probing)
    */
   if( eventtype & (SCIP_EVENTTYPE_BOUNDTIGHTENED | SCIP_EVENTTYPE_VARFIXED) )
   {
      SCIP_CONSDATA* consdata;
      int c;

      for( c = 0; c < ownerdata->nconss; ++c )
      {
         assert(ownerdata->conss[c] != NULL);
         consdata = SCIPconsGetData(ownerdata->conss[c]);

         /* if bound tightening, then mark constraints to be propagated again
          * TODO we could try be more selective here and only trigger a propagation if a relevant bound has changed,
          *   that is, we don't need to repropagate x + ... <= rhs if only the upper bound of x has been tightened
          *   the locks don't help since they are not available separately for each constraint
          */
         if( eventtype & SCIP_EVENTTYPE_BOUNDTIGHTENED )
         {
            consdata->ispropagated = FALSE;
            SCIPdebugMsg(scip, "  marked <%s> for propagate\n", SCIPconsGetName(ownerdata->conss[c]));
         }

         /* if still in presolve (but not probing), then mark constraints to be unsimplified */
         if( SCIPgetStage(scip) == SCIP_STAGE_PRESOLVING && !SCIPinProbing(scip) )
         {
            consdata->issimplified = FALSE;
            SCIPdebugMsg(scip, "  marked <%s> for simplify\n", SCIPconsGetName(ownerdata->conss[c]));
         }
      }
   }

   /* update curboundstag, lastboundrelax, and expr activity */
   if( eventtype & SCIP_EVENTTYPE_BOUNDCHANGED )
   {
      SCIP_CONSHDLRDATA* conshdlrdata;
      SCIP_INTERVAL activity;

      conshdlrdata = SCIPconshdlrGetData(ownerdata->conshdlr);
      assert(conshdlrdata != NULL);

      /* increase tag on bounds */
      ++conshdlrdata->curboundstag;
      assert(conshdlrdata->curboundstag > 0);

      /* remember also if we relaxed bounds now */
      if( eventtype & SCIP_EVENTTYPE_BOUNDRELAXED )
         conshdlrdata->lastboundrelax = conshdlrdata->curboundstag;

      /* update the activity of the var-expr here immediately
       * (we could call expr->activity = intevalvar(var, consdhlr) directly, but then the exprhdlr statistics are not updated)
       */
      SCIP_CALL( SCIPcallExprInteval(scip, expr, &activity, conshdlrdata->intevalvar, conshdlrdata) );
      /* activity = conshdlrdata->intevalvar(scip, SCIPgetVarExprVar(expr), conshdlrdata); */
#ifdef DEBUG_PROP
      SCIPdebugMsg(scip, "  var-exprhdlr::inteval = [%.20g, %.20g]\n", activity.inf, activity.sup);
#endif
      SCIPexprSetActivity(expr, activity, conshdlrdata->curboundstag);
   }

   return SCIP_OKAY;
}

/** registers event handler to catch variable events on variable
 *
 * Additionally, the given constraint is stored in the ownerdata of the variable-expression.
 * When an event occurs, all stored constraints are notified.
 */
static
SCIP_RETCODE catchVarEvent(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EVENTHDLR*       eventhdlr,          /**< event handler */
   SCIP_EXPR*            expr,               /**< variable expression */
   SCIP_CONS*            cons                /**< nonlinear constraint */
   )
{
   SCIP_EXPR_OWNERDATA* ownerdata;

   assert(eventhdlr != NULL);
   assert(expr != NULL);
   assert(SCIPisExprVar(scip, expr));
   assert(cons != NULL);

   ownerdata = SCIPexprGetOwnerData(expr);
   assert(ownerdata != NULL);

#ifndef NDEBUG
   /* assert that constraint does not double-catch variable */
   {
      int i;
      for( i = 0; i < ownerdata->nconss; ++i )
         assert(ownerdata->conss[i] != cons);
   }
#endif

   /* append cons to ownerdata->conss */
   SCIP_CALL( SCIPensureBlockMemoryArray(scip, &ownerdata->conss, &ownerdata->consssize, ownerdata->nconss + 1) );
   ownerdata->conss[ownerdata->nconss++] = cons;
   /* we're not capturing the constraint here to avoid circular references */

   /* updated sorted flag */
   if( ownerdata->nconss <= 1 )
      ownerdata->consssorted = TRUE;
   else if( ownerdata->consssorted )
      ownerdata->consssorted = compIndexConsNonlinear(ownerdata->conss[ownerdata->nconss-2], ownerdata->conss[ownerdata->nconss-1]) > 0;

   /* catch variable events, if not done so yet (first constraint) */
   if( ownerdata->filterpos < 0 )
   {
      SCIP_EVENTTYPE eventtype;

      assert(ownerdata->nconss == 1);

      eventtype = SCIP_EVENTTYPE_BOUNDCHANGED | SCIP_EVENTTYPE_VARFIXED;

      SCIP_CALL( SCIPcatchVarEvent(scip, SCIPgetVarExprVar(expr), eventtype, eventhdlr, (SCIP_EVENTDATA*)expr, &ownerdata->filterpos) );
      assert(ownerdata->filterpos >= 0);
   }

   return SCIP_OKAY;
}

/** catch variable events */
static
SCIP_RETCODE catchVarEvents(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EVENTHDLR*       eventhdlr,          /**< event handler */
   SCIP_CONS*            cons                /**< constraint for which to catch bound change events */
   )
{
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_CONSDATA* consdata;
   SCIP_EXPR* expr;
   int i;

   assert(eventhdlr != NULL);
   assert(cons != NULL);

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);
   assert(consdata->varexprs != NULL);
   assert(consdata->nvarexprs >= 0);

   /* check if we have catched variable events already */
   if( consdata->catchedevents )
      return SCIP_OKAY;

   conshdlrdata = SCIPconshdlrGetData(SCIPconsGetHdlr(cons));
   assert(conshdlrdata != NULL);
   assert(conshdlrdata->intevalvar == intEvalVarBoundTightening);

   SCIPdebugMsg(scip, "catchVarEvents for %s\n", SCIPconsGetName(cons));

   for( i = 0; i < consdata->nvarexprs; ++i )
   {
      expr = consdata->varexprs[i];

      assert(expr != NULL);
      assert(SCIPisExprVar(scip, expr));

      SCIP_CALL( catchVarEvent(scip, eventhdlr, expr, cons) );

      /* from now on, activity of var-expr will usually be updated in processVarEvent if variable bound is changing
       * since we just registered this eventhdlr, we should make sure that the activity is also up to date now
       */
      if( SCIPexprGetActivityTag(expr) < conshdlrdata->curboundstag )
      {
         SCIP_INTERVAL activity;
         SCIP_CALL( SCIPcallExprInteval(scip, expr, &activity, intEvalVarBoundTightening, conshdlrdata) );
         /* activity = intEvalVarBoundTightening(scip, SCIPgetVarExprVar(expr), conshdlrdata); */
         SCIPexprSetActivity(expr, activity, conshdlrdata->curboundstag);
#ifdef DEBUG_PROP
         SCIPdebugMsg(scip, "var-exprhdlr::inteval for var <%s> = [%.20g, %.20g]\n", SCIPvarGetName(SCIPgetVarExprVar(expr)), activity.inf, activity.sup);
#endif
      }
   }

   consdata->catchedevents = TRUE;

   return SCIP_OKAY;
}

/** unregisters event handler to catch variable events on variable
 *
 * The given constraint is removed from the constraints array in the ownerdata of the variable-expression.
 * If this was the last constraint, then the event handler is unregistered for this variable.
 */
static
SCIP_RETCODE dropVarEvent(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EVENTHDLR*       eventhdlr,          /**< event handler */
   SCIP_EXPR*            expr,               /**< variable expression */
   SCIP_CONS*            cons                /**< expr constraint */
   )
{
   SCIP_EXPR_OWNERDATA* ownerdata;
   int pos;

   assert(eventhdlr != NULL);
   assert(expr != NULL);
   assert(SCIPisExprVar(scip, expr));
   assert(cons != NULL);

   ownerdata = SCIPexprGetOwnerData(expr);
   assert(ownerdata != NULL);
   assert(ownerdata->nconss > 0);

   if( ownerdata->conss[ownerdata->nconss-1] == cons )
   {
      pos = ownerdata->nconss-1;
   }
   else
   {
      if( !ownerdata->consssorted )
      {
         SCIPsortPtr((void**)ownerdata->conss, compIndexConsNonlinear, ownerdata->nconss);
         ownerdata->consssorted = TRUE;
      }

      if( !SCIPsortedvecFindPtr((void**)ownerdata->conss, compIndexConsNonlinear, cons, ownerdata->nconss, &pos) )
      {
         SCIPerrorMessage("Constraint <%s> not in constraint array of expression for variable <%s>\n", SCIPconsGetName(cons), SCIPvarGetName(SCIPgetVarExprVar(expr)));
         return SCIP_ERROR;
      }
      assert(pos >= 0 && pos < ownerdata->nconss);
   }
   assert(ownerdata->conss[pos] == cons);

   /* move last constraint into position of removed constraint */
   if( pos < ownerdata->nconss-1 )
   {
      ownerdata->conss[pos] = ownerdata->conss[ownerdata->nconss-1];
      ownerdata->consssorted = FALSE;
   }
   --ownerdata->nconss;

   /* drop variable events if that was the last constraint */
   if( ownerdata->nconss == 0 )
   {
      SCIP_EVENTTYPE eventtype;

      assert(ownerdata->filterpos >= 0);

      eventtype = SCIP_EVENTTYPE_BOUNDCHANGED | SCIP_EVENTTYPE_VARFIXED;

      SCIP_CALL( SCIPdropVarEvent(scip, SCIPgetVarExprVar(expr), eventtype, eventhdlr, (SCIP_EVENTDATA*)expr, ownerdata->filterpos) );
      ownerdata->filterpos = -1;
   }

   return SCIP_OKAY;
}

/** drop variable events */
static
SCIP_RETCODE dropVarEvents(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EVENTHDLR*       eventhdlr,          /**< event handler */
   SCIP_CONS*            cons                /**< constraint for which to drop bound change events */
   )
{
   SCIP_CONSDATA* consdata;
   int i;

   assert(eventhdlr != NULL);
   assert(cons != NULL);

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);

   /* check if we have catched variable events already */
   if( !consdata->catchedevents )
      return SCIP_OKAY;

   assert(consdata->varexprs != NULL);
   assert(consdata->nvarexprs >= 0);

   SCIPdebugMsg(scip, "dropVarEvents for %s\n", SCIPconsGetName(cons));

   for( i = consdata->nvarexprs - 1; i >= 0; --i )
   {
      assert(consdata->varexprs[i] != NULL);

      SCIP_CALL( dropVarEvent(scip, eventhdlr, consdata->varexprs[i], cons) );
   }

   consdata->catchedevents = FALSE;

   return SCIP_OKAY;
}

/** creates and captures a nonlinear constraint
 *
 * @attention Use copyexpr=FALSE only if expr is already "owned" by conshdlr, that is, if expressions were created with exprownerCreate() and ownerdata passed in the last two arguments
 */
static
SCIP_RETCODE createCons(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS**           cons,               /**< pointer to hold the created constraint */
   const char*           name,               /**< name of constraint */
   SCIP_EXPR*            expr,               /**< expression of constraint (must not be NULL) */
   SCIP_Real             lhs,                /**< left hand side of constraint */
   SCIP_Real             rhs,                /**< right hand side of constraint */
   SCIP_Bool             copyexpr,           /**< whether to copy the expression or reuse the given expr (capture it) */
   SCIP_Bool             initial,            /**< should the LP relaxation of constraint be in the initial LP?
                                              *   Usually set to TRUE. Set to FALSE for 'lazy constraints'. */
   SCIP_Bool             separate,           /**< should the constraint be separated during LP processing?
                                              *   Usually set to TRUE. */
   SCIP_Bool             enforce,            /**< should the constraint be enforced during node processing?
                                              *   TRUE for model constraints, FALSE for additional, redundant constraints. */
   SCIP_Bool             check,              /**< should the constraint be checked for feasibility?
                                              *   TRUE for model constraints, FALSE for additional, redundant constraints. */
   SCIP_Bool             propagate,          /**< should the constraint be propagated during node processing?
                                              *   Usually set to TRUE. */
   SCIP_Bool             local,              /**< is constraint only valid locally?
                                              *   Usually set to FALSE. Has to be set to TRUE, e.g., for branching constraints. */
   SCIP_Bool             modifiable,         /**< is constraint modifiable (subject to column generation)?
                                              *   Usually set to FALSE. In column generation applications, set to TRUE if pricing
                                              *   adds coefficients to this constraint. */
   SCIP_Bool             dynamic,            /**< is constraint subject to aging?
                                              *   Usually set to FALSE. Set to TRUE for own cuts which
                                              *   are separated as constraints. */
   SCIP_Bool             removable           /**< should the relaxation be removed from the LP due to aging or cleanup?
                                              *   Usually set to FALSE. Set to TRUE for 'lazy constraints' and 'user cuts'. */
   )
{
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_CONSDATA* consdata;

   assert(conshdlr != NULL);
   assert(expr != NULL);

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);

   if( local && SCIPgetDepth(scip) != 0 )
   {
      SCIPerrorMessage("Locally valid nonlinear constraints are not supported, yet.\n");
      return SCIP_INVALIDCALL;
   }

   /* TODO we should allow for non-initial nonlinear constraints */
   if( !initial )
   {
      SCIPerrorMessage("Non-initial nonlinear constraints are not supported, yet.\n");
      return SCIP_INVALIDCALL;
   }

   /* create constraint data */
   SCIP_CALL( SCIPallocClearBlockMemory(scip, &consdata) );

   if( copyexpr )
   {
      /* copy expression, thereby map variables expressions to already existing variables expressions in var2expr map, or augment var2expr map */
      SCIP_CALL( SCIPduplicateExpr(scip, expr, &consdata->expr, mapexprvar, conshdlr, exprownerCreate, (void*)conshdlr) );
   }
   else
   {
      consdata->expr = expr;
      SCIPcaptureExpr(consdata->expr);
   }
   consdata->lhs = lhs;
   consdata->rhs = rhs;
   consdata->consindex = conshdlrdata->lastconsindex++;
   consdata->curv = SCIP_EXPRCURV_UNKNOWN;

   /* create constraint */
   SCIP_CALL( SCIPcreateCons(scip, cons, name, conshdlr, consdata, initial, separate, enforce, check, propagate,
         local, modifiable, dynamic, removable, FALSE) );

   return SCIP_OKAY;
}

/** returns absolute violation for auxvar relation in an expression w.r.t. original variables
 *
 * Assume the expression is f(x), where x are original (i.e., not auxiliary) variables.
 * Assume that f(x) is associated with auxiliary variable z.
 *
 * If there are negative locks, then return the violation of z &le; f(x) and sets `violover` to TRUE.
 * If there are positive locks, then return the violation of z &ge; f(x) and sets `violunder` to TRUE.
 * Of course, if there both negative and positive locks, then return the violation of z = f(x).
 * If f could not be evaluated, then return SCIPinfinity() and set both `violover` and `violunder` to TRUE.
 *
 * @note This does not reevaluate the violation, but assumes that the expression has been evaluated
 */
static
SCIP_Real getExprAbsOrigViolation(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EXPR*            expr,               /**< expression */
   SCIP_SOL*             sol,                /**< solution that has been evaluated */
   SCIP_Bool*            violunder,          /**< buffer to store whether z >= f(x) is violated, or NULL */
   SCIP_Bool*            violover            /**< buffer to store whether z <= f(x) is violated, or NULL */
   )
{
   SCIP_EXPR_OWNERDATA* ownerdata;
   SCIP_Real auxvarvalue;

   assert(expr != NULL);

   ownerdata = SCIPexprGetOwnerData(expr);
   assert(ownerdata != NULL);
   assert(ownerdata->auxvar != NULL);

   if( SCIPexprGetEvalValue(expr) == SCIP_INVALID )
   {
      if( violunder != NULL )
         *violunder = TRUE;
      if( violover != NULL )
         *violover = TRUE;
      return SCIPinfinity(scip);
   }

   auxvarvalue = SCIPgetSolVal(scip, sol, ownerdata->auxvar);

   if( ownerdata->nlocksneg > 0 && auxvarvalue > SCIPexprGetEvalValue(expr) )
   {
      if( violunder != NULL )
         *violunder = FALSE;
      if( violover != NULL )
         *violover = TRUE;
      return auxvarvalue - SCIPexprGetEvalValue(expr);
   }

   if( ownerdata->nlockspos > 0 && SCIPexprGetEvalValue(expr) > auxvarvalue )
   {
      if( violunder != NULL )
         *violunder = TRUE;
      if( violover != NULL )
         *violover = FALSE;
      return SCIPexprGetEvalValue(expr) - auxvarvalue;
   }

   if( violunder != NULL )
      *violunder = FALSE;
   if( violover != NULL )
      *violover = FALSE;
   return 0.0;
}

/** returns absolute violation for auxvar relation in an expression w.r.t. auxiliary variables
 *
 * Assume the expression is f(w), where w are auxiliary variables that were introduced by some nlhdlr.
 * Assume that f(w) is associated with auxiliary variable z.
 *
 * If there are negative locks, then return the violation of z &le; f(w) and sets `violover` to TRUE.
 * If there are positive locks, then return the violation of z &ge; f(w) and sets `violunder` to TRUE.
 * Of course, if there both negative and positive locks, then return the violation of z = f(w).
 * If f could not be evaluated, then return SCIPinfinity() and set both `violover` and `violunder` to TRUE.
 *
 * @note This does not reevaluate the violation, but assumes that f(w) is passed in with auxvalue.
 */
static
SCIP_Real getExprAbsAuxViolation(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EXPR*            expr,               /**< expression */
   SCIP_Real             auxvalue,           /**< value of f(w) */
   SCIP_SOL*             sol,                /**< solution that has been evaluated */
   SCIP_Bool*            violunder,          /**< buffer to store whether z >= f(w) is violated, or NULL */
   SCIP_Bool*            violover            /**< buffer to store whether z <= f(w) is violated, or NULL */
   )
{
   SCIP_EXPR_OWNERDATA* ownerdata;
   SCIP_Real auxvarvalue;

   assert(expr != NULL);

   ownerdata = SCIPexprGetOwnerData(expr);
   assert(ownerdata != NULL);
   assert(ownerdata->auxvar != NULL);

   if( auxvalue == SCIP_INVALID )
   {
      if( violunder != NULL )
         *violunder = TRUE;
      if( violover != NULL )
         *violover = TRUE;
      return SCIPinfinity(scip);
   }

   auxvarvalue = SCIPgetSolVal(scip, sol, ownerdata->auxvar);

   if( ownerdata->nlocksneg > 0 && auxvarvalue > auxvalue )
   {
      if( violunder != NULL )
         *violunder = FALSE;
      if( violover != NULL )
         *violover = TRUE;
      return auxvarvalue - auxvalue;
   }

   if( ownerdata->nlockspos > 0 && auxvalue > auxvarvalue )
   {
      if( violunder != NULL )
         *violunder = TRUE;
      if( violover != NULL )
         *violover = FALSE;
      return auxvalue - auxvarvalue;
   }

   if( violunder != NULL )
      *violunder = FALSE;
   if( violover != NULL )
      *violover = FALSE;

   return 0.0;
}

/** computes violation of a constraint */
static
SCIP_RETCODE computeViolation(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons,               /**< constraint */
   SCIP_SOL*             sol,                /**< solution or NULL if LP solution should be used */
   SCIP_Longint          soltag              /**< tag that uniquely identifies the solution (with its values), or 0. */
   )
{
   SCIP_CONSDATA* consdata;
   SCIP_Real activity;

   assert(scip != NULL);
   assert(cons != NULL);

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);

   SCIP_CALL( SCIPevalExpr(scip, consdata->expr, sol, soltag) );
   activity = SCIPexprGetEvalValue(consdata->expr);

   /* consider constraint as violated if it is undefined in the current point */
   if( activity == SCIP_INVALID )
   {
      consdata->lhsviol = SCIPinfinity(scip);
      consdata->rhsviol = SCIPinfinity(scip);
      return SCIP_OKAY;
   }

   /* compute violations */
   consdata->lhsviol = SCIPisInfinity(scip, -consdata->lhs) ? -SCIPinfinity(scip) : consdata->lhs  - activity;
   consdata->rhsviol = SCIPisInfinity(scip,  consdata->rhs) ? -SCIPinfinity(scip) : activity - consdata->rhs;

   return SCIP_OKAY;
}

/** returns absolute violation of a constraint
 *
 * @note This does not reevaluate the violation, but assumes that computeViolation() has been called before.
 */
static
SCIP_Real getConsAbsViolation(
   SCIP_CONS*            cons                /**< constraint */
   )
{
   SCIP_CONSDATA* consdata;

   assert(cons != NULL);

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);

   return MAX3(0.0, consdata->lhsviol, consdata->rhsviol);
}

/** computes relative violation of a constraint
 *
 * @note This does not reevaluate the violation, but assumes that computeViolation() has been called before.
 */
static
SCIP_RETCODE getConsRelViolation(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons,               /**< constraint */
   SCIP_Real*            viol,               /**< buffer to store violation */
   SCIP_SOL*             sol,                /**< solution or NULL if LP solution should be used */
   SCIP_Longint          soltag              /**< tag that uniquely identifies the solution (with its values), or 0 */
   )
{
   SCIP_CONSHDLR* conshdlr;
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_CONSDATA* consdata;
   SCIP_Real scale;

   assert(cons != NULL);
   assert(viol != NULL);

   conshdlr = SCIPconsGetHdlr(cons);
   assert(conshdlr != NULL);

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);

   *viol = getConsAbsViolation(cons);

   if( conshdlrdata->violscale == 'n' )
      return SCIP_OKAY;

   if( SCIPisInfinity(scip, *viol) )
      return SCIP_OKAY;

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);

   if( conshdlrdata->violscale == 'a' )
   {
      scale = MAX(1.0, REALABS(SCIPexprGetEvalValue(consdata->expr)));

      /* consider value of side that is violated for scaling, too */
      if( consdata->lhsviol > 0.0 && REALABS(consdata->lhs) > scale )
      {
         assert(!SCIPisInfinity(scip, -consdata->lhs));
         scale = REALABS(consdata->lhs);
      }
      else if( consdata->rhsviol > 0.0 && REALABS(consdata->rhs) > scale )
      {
         assert(!SCIPisInfinity(scip,  consdata->rhs));
         scale = REALABS(consdata->rhs);
      }

      *viol /= scale;
      return SCIP_OKAY;
   }

   /* if not 'n' or 'a', then it has to be 'g' at the moment */
   assert(conshdlrdata->violscale == 'g');
   if( soltag == 0L || consdata->gradnormsoltag != soltag )
   {
      /* we need the varexprs to conveniently access the gradient */
      SCIP_CALL( storeVarExprs(scip, conshdlr, consdata) );

      /* update cached value of norm of gradient */
      consdata->gradnorm = 0.0;

      /* compute gradient */
      SCIP_CALL( SCIPevalExprGradient(scip, consdata->expr, sol, soltag) );

      /* gradient evaluation error -> no scaling */
      if( SCIPexprGetDerivative(consdata->expr) != SCIP_INVALID )
      {
         int i;
         for( i = 0; i < consdata->nvarexprs; ++i )
         {
            SCIP_Real deriv;

            assert(SCIPexprGetDiffTag(consdata->expr) == SCIPexprGetDiffTag(consdata->varexprs[i]));
            deriv = SCIPexprGetDerivative(consdata->varexprs[i]);
            if( deriv == SCIP_INVALID )
            {
               /* SCIPdebugMsg(scip, "gradient evaluation error for component %d\n", i); */
               consdata->gradnorm = 0.0;
               break;
            }

            consdata->gradnorm += deriv*deriv;
         }
      }
      consdata->gradnorm = sqrt(consdata->gradnorm);
      consdata->gradnormsoltag = soltag;
   }

   *viol /= MAX(1.0, consdata->gradnorm);

   return SCIP_OKAY;
}

/** returns whether constraint is currently violated
 *
 * @note This does not reevaluate the violation, but assumes that computeViolation() has been called before.
 */
static
SCIP_Bool isConsViolated(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons                /**< constraint */
   )
{
   return getConsAbsViolation(cons) > SCIPfeastol(scip);
}

/** checks for a linear variable that can be increased or decreased without harming feasibility */
static
void findUnlockedLinearVar(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons                /**< constraint */
   )
{
   SCIP_CONSDATA* consdata;
   int poslock;
   int neglock;
   int i;

   assert(cons != NULL);

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);

   consdata->linvarincr = NULL;
   consdata->linvardecr = NULL;
   consdata->linvarincrcoef = 0.0;
   consdata->linvardecrcoef = 0.0;

   /* root expression is not a sum -> no unlocked linear variable available */
   if( !SCIPisExprSum(scip, consdata->expr) )
      return;

   for( i = 0; i < SCIPexprGetNChildren(consdata->expr); ++i )
   {
      SCIP_EXPR* child;

      child = SCIPexprGetChildren(consdata->expr)[i];
      assert(child != NULL);

      /* check whether the child is a variable expression */
      if( SCIPisExprVar(scip, child) )
      {
         SCIP_VAR* var = SCIPgetVarExprVar(child);
         SCIP_Real coef = SCIPgetCoefsExprSum(consdata->expr)[i];

         if( coef > 0.0 )
         {
            poslock = !SCIPisInfinity(scip,  consdata->rhs) ? 1 : 0;
            neglock = !SCIPisInfinity(scip, -consdata->lhs) ? 1 : 0;
         }
         else
         {
            poslock = !SCIPisInfinity(scip, -consdata->lhs) ? 1 : 0;
            neglock = !SCIPisInfinity(scip,  consdata->rhs) ? 1 : 0;
         }
         SCIPdebugMsg(scip, "child <%s> locks: %d %d\n", SCIPvarGetName(var), SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL), SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL));

         if( SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) - neglock == 0 )
         {
            /* for a*x + f(y) \in [lhs, rhs], we can decrease x without harming other constraints
             * if we have already one candidate, then take the one where the loss in the objective function is less
             */
            if( (consdata->linvardecr == NULL) ||
               (SCIPvarGetObj(consdata->linvardecr) / consdata->linvardecrcoef > SCIPvarGetObj(var) / coef) )
            {
               consdata->linvardecr = var;
               consdata->linvardecrcoef = coef;
            }
         }

         if( SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) - poslock == 0 )
         {
            /* for a*x + f(y) \in [lhs, rhs], we can increase x without harm
             * if we have already one candidate, then take the one where the loss in the objective function is less
             */
            if( (consdata->linvarincr == NULL) ||
               (SCIPvarGetObj(consdata->linvarincr) / consdata->linvarincrcoef > SCIPvarGetObj(var) / coef) )
            {
               consdata->linvarincr = var;
               consdata->linvarincrcoef = coef;
            }
         }
      }
   }

   assert(consdata->linvarincr == NULL || consdata->linvarincrcoef != 0.0);
   assert(consdata->linvardecr == NULL || consdata->linvardecrcoef != 0.0);

   if( consdata->linvarincr != NULL )
   {
      SCIPdebugMsg(scip, "may increase <%s> to become feasible\n", SCIPvarGetName(consdata->linvarincr));
   }
   if( consdata->linvardecr != NULL )
   {
      SCIPdebugMsg(scip, "may decrease <%s> to become feasible\n", SCIPvarGetName(consdata->linvardecr));
   }
}

/** Given a solution where every nonlinear constraint is either feasible or can be made feasible by
 *  moving a linear variable, construct the corresponding feasible solution and pass it to the trysol heuristic.
 *
 *  The method assumes that this is always possible and that not all constraints are feasible already.
 */
static
SCIP_RETCODE proposeFeasibleSolution(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS**           conss,              /**< constraints to process */
   int                   nconss,             /**< number of constraints */
   SCIP_SOL*             sol,                /**< solution to process */
   SCIP_Bool*            success             /**< buffer to store whether we succeeded to construct a solution that satisfies all provided constraints */
   )
{
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_SOL* newsol;
   int c;

   assert(scip  != NULL);
   assert(conshdlr != NULL);
   assert(conss != NULL || nconss == 0);
   assert(success != NULL);

   *success = FALSE;

   /* don't propose new solutions if not in presolve or solving */
   if( SCIPgetStage(scip) < SCIP_STAGE_INITPRESOLVE || SCIPgetStage(scip) >= SCIP_STAGE_SOLVED )
      return SCIP_OKAY;

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);

   if( sol != NULL )
   {
      SCIP_CALL( SCIPcreateSolCopy(scip, &newsol, sol) );
   }
   else
   {
      SCIP_CALL( SCIPcreateLPSol(scip, &newsol, NULL) );
   }
   SCIP_CALL( SCIPunlinkSol(scip, newsol) );
   SCIPdebugMsg(scip, "attempt to make solution from <%s> feasible by shifting linear variable\n",
      sol != NULL ? (SCIPsolGetHeur(sol) != NULL ? SCIPheurGetName(SCIPsolGetHeur(sol)) : "tree") : "LP");

   for( c = 0; c < nconss; ++c )
   {
      SCIP_CONSDATA* consdata = SCIPconsGetData(conss[c]);  /*lint !e613*/
      SCIP_Real viol = 0.0;
      SCIP_Real delta;
      SCIP_Real gap;

      assert(consdata != NULL);

      /* get absolute violation and sign */
      if( consdata->lhsviol > SCIPfeastol(scip) )
         viol = consdata->lhsviol; /* lhs - activity */
      else if( consdata->rhsviol > SCIPfeastol(scip) )
         viol = -consdata->rhsviol; /* rhs - activity */
      else
         continue; /* constraint is satisfied */

      if( consdata->linvarincr != NULL &&
         ((viol > 0.0 && consdata->linvarincrcoef > 0.0) || (viol < 0.0 && consdata->linvarincrcoef < 0.0)) )
      {
         SCIP_VAR* var = consdata->linvarincr;

         /* compute how much we would like to increase var */
         delta = viol / consdata->linvarincrcoef;
         assert(delta > 0.0);

         /* if var has an upper bound, may need to reduce delta */
         if( !SCIPisInfinity(scip, SCIPvarGetUbGlobal(var)) )
         {
            gap = SCIPvarGetUbGlobal(var) - SCIPgetSolVal(scip, newsol, var);
            delta = MIN(MAX(0.0, gap), delta);
         }
         if( SCIPisPositive(scip, delta) )
         {
            /* if variable is integral, round delta up so that it will still have an integer value */
            if( SCIPvarIsIntegral(var) )
               delta = SCIPceil(scip, delta);

            SCIP_CALL( SCIPincSolVal(scip, newsol, var, delta) );
            SCIPdebugMsg(scip, "increase <%s> by %g to %g to remedy lhs-violation %g of cons <%s>\n",
               SCIPvarGetName(var), delta, SCIPgetSolVal(scip, newsol, var), viol, SCIPconsGetName(conss[c]));  /*lint !e613*/

            /* adjust constraint violation, if satisfied go on to next constraint */
            viol -= consdata->linvarincrcoef * delta;
            if( SCIPisZero(scip, viol) )
               continue;
         }
      }

      assert(viol != 0.0);
      if( consdata->linvardecr != NULL &&
         ((viol > 0.0 && consdata->linvardecrcoef < 0.0) || (viol < 0.0 && consdata->linvardecrcoef > 0.0)) )
      {
         SCIP_VAR* var = consdata->linvardecr;

         /* compute how much we would like to decrease var */
         delta = viol / consdata->linvardecrcoef;
         assert(delta < 0.0);

         /* if var has a lower bound, may need to reduce delta */
         if( !SCIPisInfinity(scip, -SCIPvarGetLbGlobal(var)) )
         {
            gap = SCIPgetSolVal(scip, newsol, var) - SCIPvarGetLbGlobal(var);
            delta = MAX(MIN(0.0, gap), delta);
         }
         if( SCIPisNegative(scip, delta) )
         {
            /* if variable is integral, round delta down so that it will still have an integer value */
            if( SCIPvarIsIntegral(var) )
               delta = SCIPfloor(scip, delta);
            SCIP_CALL( SCIPincSolVal(scip, newsol, consdata->linvardecr, delta) );
            /*lint --e{613} */
            SCIPdebugMsg(scip, "increase <%s> by %g to %g to remedy rhs-violation %g of cons <%s>\n",
               SCIPvarGetName(var), delta, SCIPgetSolVal(scip, newsol, var), viol, SCIPconsGetName(conss[c]));

            /* adjust constraint violation, if satisfied go on to next constraint */
            viol -= consdata->linvardecrcoef * delta;
            if( SCIPisZero(scip, viol) )
               continue;
         }
      }

      /* still here... so probably we could not make constraint feasible due to variable bounds, thus give up */
      break;
   }

   /* if we have a solution that should satisfy all quadratic constraints and has a better objective than the current upper bound,
    * then pass it to the trysol heuristic
    */
   if( c == nconss && (SCIPisInfinity(scip, SCIPgetUpperbound(scip)) || SCIPisSumLT(scip, SCIPgetSolTransObj(scip, newsol), SCIPgetUpperbound(scip))) )
   {
      SCIPdebugMsg(scip, "pass solution with objective val %g to trysol heuristic\n", SCIPgetSolTransObj(scip, newsol));

      assert(conshdlrdata->trysolheur != NULL);
      SCIP_CALL( SCIPheurPassSolTrySol(scip, conshdlrdata->trysolheur, newsol) );

      *success = TRUE;
   }

   SCIP_CALL( SCIPfreeSol(scip, &newsol) );

   return SCIP_OKAY;
}

/** adds globally valid tight estimators in a given solution as cut to cutpool
 *
 * Called by addTightEstimatorCuts() for a specific expression, nlhdlr, and estimate-direction (over or under).
 */
static
SCIP_RETCODE addTightEstimatorCut(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS*            cons,               /**< constraint */
   SCIP_EXPR*            expr,               /**< expression */
   EXPRENFO*             exprenfo,           /**< expression enfo data, e.g., nlhdlr to use */
   SCIP_SOL*             sol,                /**< reference point where to estimate */
   SCIP_Bool             overestimate,       /**< whether to overestimate */
   SCIP_PTRARRAY*        rowpreps            /**< array for rowpreps */
   )
{
   SCIP_Bool estimatesuccess = FALSE;
   SCIP_Bool branchscoresuccess = FALSE;
   int minidx;
   int maxidx;
   int r;

   assert(scip != NULL);
   assert(expr != NULL);
   assert(exprenfo != NULL);
   assert(rowpreps != NULL);

   ENFOLOG( SCIPinfoMessage(scip, enfologfile, "   %sestimate using nlhdlr <%s> for expr %p (%s)\n",
      overestimate ? "over" : "under", SCIPnlhdlrGetName(exprenfo->nlhdlr), (void*)expr, SCIPexprhdlrGetName(SCIPexprGetHdlr(expr))); )

   SCIP_CALL( SCIPnlhdlrEstimate(scip, conshdlr, exprenfo->nlhdlr, expr, exprenfo->nlhdlrexprdata, sol,
      exprenfo->auxvalue, overestimate, SCIPinfinity(scip), FALSE, rowpreps, &estimatesuccess, &branchscoresuccess) );

   minidx = SCIPgetPtrarrayMinIdx(scip, rowpreps);
   maxidx = SCIPgetPtrarrayMaxIdx(scip, rowpreps);
   assert(estimatesuccess == (minidx <= maxidx));

   if( !estimatesuccess )
   {
      ENFOLOG( SCIPinfoMessage(scip, enfologfile, "    estimate of nlhdlr %s failed\n", SCIPnlhdlrGetName(exprenfo->nlhdlr)); )
      return SCIP_OKAY;
   }

   for( r = minidx; r <= maxidx; ++r )
   {
      SCIP_ROWPREP* rowprep;
      SCIP_ROW* row;
      SCIP_Real estimateval;
      int i;

      rowprep = (SCIP_ROWPREP*) SCIPgetPtrarrayVal(scip, rowpreps, r);
      assert(rowprep != NULL);
      assert(SCIProwprepGetSidetype(rowprep) == (overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT));

      /* if estimators is only local valid, then skip */
      if( SCIProwprepIsLocal(rowprep) )
      {
         ENFOLOG( SCIPinfoMessage(scip, enfologfile, "    skip local estimator\n"); )
         SCIPfreeRowprep(scip, &rowprep);
         continue;
      }

      /* compute value of estimator */
      estimateval = -SCIProwprepGetSide(rowprep);
      for( i = 0; i < SCIProwprepGetNVars(rowprep); ++i )
         estimateval += SCIProwprepGetCoefs(rowprep)[i] * SCIPgetSolVal(scip, sol, SCIProwprepGetVars(rowprep)[i]);

      /* if estimator value is not tight (or even "more than tight", e.g., when estimating in integer vars), then skip */
      if( (overestimate && !SCIPisFeasLE(scip, estimateval, SCIPexprGetEvalValue(expr))) ||
         (!overestimate && !SCIPisFeasGE(scip, estimateval, SCIPexprGetEvalValue(expr))) )
      {
         ENFOLOG( SCIPinfoMessage(scip, enfologfile, "    skip non-tight estimator with value %g, expr value %g\n", estimateval, SCIPexprGetEvalValue(expr)); )
         SCIPfreeRowprep(scip, &rowprep);
         continue;
      }

      /* complete estimator to cut and clean it up */
      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPgetExprAuxVarNonlinear(expr), -1.0) );
      SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, sol, SCIPinfinity(scip), &estimatesuccess) );

      /* if cleanup failed or rowprep is local now, then skip */
      if( !estimatesuccess || SCIProwprepIsLocal(rowprep) )
      {
         ENFOLOG( SCIPinfoMessage(scip, enfologfile, "    skip after cleanup failed or made estimator locally valid\n"); )
         SCIPfreeRowprep(scip, &rowprep);
         continue;
      }

      /* generate row and add to cutpool */
      SCIP_CALL( SCIPgetRowprepRowCons(scip, &row, rowprep, cons) );

      ENFOLOG( SCIPinfoMessage(scip, enfologfile, "    adding cut ");
      SCIP_CALL( SCIPprintRow(scip, row, enfologfile) ); )

      SCIP_CALL( SCIPaddPoolCut(scip, row) );
      /* SCIPnlhdlrIncrementNSeparated(nlhdlr); */

      SCIP_CALL( SCIPreleaseRow(scip, &row) );
      SCIPfreeRowprep(scip, &rowprep);
   }

   SCIP_CALL( SCIPclearPtrarray(scip, rowpreps) );

   return SCIP_OKAY;
}

/** adds globally valid tight estimators in a given solution as cuts to cutpool
 *
 * Essentially we want to ensure that the LP relaxation is tight in the new solution, if possible.
 * For convex constraints, we would achieve this by linearizing.
 * To avoid checking explicitly for convexity, we compute estimators via any nlhdlr that didn't say it would
 * use bound information and check whether the estimator is tight.
 *
 * Since linearization may happen in auxiliary variables, we ensure that auxiliary variables are set
 * to the eval-value of its expression, i.e., we change sol so it is also feasible in the extended formulation.
 */
static
SCIP_RETCODE addTightEstimatorCuts(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS**           conss,              /**< constraints */
   int                   nconss,             /**< number of constraints */
   SCIP_SOL*             sol                 /**< reference point where to estimate */
   )
{
   SCIP_CONSDATA* consdata;
   SCIP_Longint soltag;
   SCIP_EXPRITER* it;
   SCIP_EXPR* expr;
   SCIP_PTRARRAY* rowpreps;
   int c, e;

   assert(scip != NULL);
   assert(conshdlr != NULL);
   assert(conss != NULL || nconss == 0);

   ENFOLOG( SCIPinfoMessage(scip, enfologfile, "add tight estimators in new solution from <%s> to cutpool\n", SCIPheurGetName(SCIPsolGetHeur(sol))); )

   /* TODO probably we just evaluated all expressions when checking the sol before it was added
    * would be nice to recognize this and skip reevaluating
    */
   soltag = SCIPgetExprNewSoltag(scip);

   SCIP_CALL( SCIPcreatePtrarray(scip, &rowpreps) );

   SCIP_CALL( SCIPcreateExpriter(scip, &it) );
   SCIP_CALL( SCIPexpriterInit(it, NULL, SCIP_EXPRITER_DFS, FALSE) );
   SCIPexpriterSetStagesDFS(it, SCIP_EXPRITER_LEAVEEXPR);

   for( c = 0; c < nconss; ++c )
   {
      /* skip constraints that are not enabled or deleted or have separation disabled */
      if( !SCIPconsIsEnabled(conss[c]) || SCIPconsIsDeleted(conss[c]) || !SCIPconsIsSeparationEnabled(conss[c]) )
         continue;
      assert(SCIPconsIsActive(conss[c]));

      consdata = SCIPconsGetData(conss[c]);
      assert(consdata != NULL);

      /* TODO we could remember for which constraints there is a chance that we would add anything,
       * i.e., there is some convex-like expression, and skip other constraints
       */

      ENFOLOG(
      {
         int i;
         SCIPinfoMessage(scip, enfologfile, " constraint ");
         SCIP_CALL( SCIPprintCons(scip, conss[c], enfologfile) );
         SCIPinfoMessage(scip, enfologfile, "\n and point\n");
         for( i = 0; i < consdata->nvarexprs; ++i )
         {
            SCIP_VAR* var;
            var = SCIPgetVarExprVar(consdata->varexprs[i]);
            SCIPinfoMessage(scip, enfologfile, "  %-10s = %15g bounds: [%15g,%15g]\n", SCIPvarGetName(var),
                  SCIPgetSolVal(scip, sol, var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var));
         }
      })

      SCIP_CALL( SCIPevalExpr(scip, consdata->expr, sol, soltag) );
      assert(SCIPexprGetEvalValue(consdata->expr) != SCIP_INVALID);

      for( expr = SCIPexpriterRestartDFS(it, consdata->expr); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
      {
         SCIP_EXPR_OWNERDATA* ownerdata;

         ownerdata = SCIPexprGetOwnerData(expr);
         assert(ownerdata != NULL);

         /* we can only generate a cut from an estimator if there is an auxvar */
         if( ownerdata->auxvar == NULL )
            continue;

         /* set value for auxvar in sol to value of expr, in case it is used to compute estimators higher up of this expression */
         assert(SCIPexprGetEvalTag(expr) == soltag);
         assert(SCIPexprGetEvalValue(expr) != SCIP_INVALID);
         SCIP_CALL( SCIPsetSolVal(scip, sol, ownerdata->auxvar, SCIPexprGetEvalValue(expr)) );

         /* generate cuts from estimators of each nonlinear handler that provides estimates */
         for( e = 0; e < ownerdata->nenfos; ++e )
         {
            SCIP_NLHDLR* nlhdlr;

            nlhdlr = ownerdata->enfos[e]->nlhdlr;
            assert(nlhdlr != NULL);

            /* skip nlhdlr that does not implement estimate (so it does enfo) */
            if( !SCIPnlhdlrHasEstimate(nlhdlr) )
               continue;

            /* skip nlhdlr that does not participate in separation or looks like it would give only locally-valid estimators
             * (because it uses activities on vars/auxvars)
             */
            if( ((ownerdata->enfos[e]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_SEPAABOVE) == 0 || ownerdata->enfos[e]->sepaaboveusesactivity) &&
                ((ownerdata->enfos[e]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_SEPABELOW) == 0 || ownerdata->enfos[e]->sepabelowusesactivity) )
               continue;

            /* skip nlhdlr_default on sum, as the estimator doesn't depend on the reference point (expr is linear in auxvars) */
            if( SCIPisExprSum(scip, expr) && strcmp(SCIPnlhdlrGetName(nlhdlr), "default") == 0 )
               continue;

            /* evaluate the expression w.r.t. the nlhdlrs auxiliary variables, since some nlhdlr expect this before their estimate is called */
            SCIP_CALL( SCIPnlhdlrEvalaux(scip, nlhdlr, expr, ownerdata->enfos[e]->nlhdlrexprdata, &ownerdata->enfos[e]->auxvalue, sol) );
            ENFOLOG(
               SCIPinfoMessage(scip, enfologfile, "  expr ");
               SCIPprintExpr(scip, expr, enfologfile);
               SCIPinfoMessage(scip, enfologfile, " (%p): evalvalue %.15g auxvarvalue %.15g, nlhdlr <%s> auxvalue: %.15g\n",
                  (void*)expr, SCIPexprGetEvalValue(expr), SCIPgetSolVal(scip, sol, ownerdata->auxvar), SCIPnlhdlrGetName(nlhdlr), ownerdata->enfos[e]->auxvalue);
            )
            /* due to setting values of auxvars to expr values in sol, the auxvalue should equal to expr evalvalue */
            assert(SCIPisEQ(scip, ownerdata->enfos[e]->auxvalue, SCIPexprGetEvalValue(expr)));

            /* if nlhdlr wants to be called for overestimate and does not use local bounds, then call estimate of nlhdlr */
            if( (ownerdata->enfos[e]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_SEPAABOVE) && !ownerdata->enfos[e]->sepaaboveusesactivity )
            {
               SCIP_CALL( addTightEstimatorCut(scip, conshdlr, conss[c], expr, ownerdata->enfos[e], sol, TRUE, rowpreps) );
            }

            /* if nlhdlr wants to be called for underestimate and does not use local bounds, then call estimate of nlhdlr */
            if( (ownerdata->enfos[e]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_SEPABELOW) && !ownerdata->enfos[e]->sepabelowusesactivity )
            {
               SCIP_CALL( addTightEstimatorCut(scip, conshdlr, conss[c], expr, ownerdata->enfos[e], sol, FALSE, rowpreps) );
            }
         }
      }
   }

   SCIPfreeExpriter(&it);
   SCIP_CALL( SCIPfreePtrarray(scip, &rowpreps) );

   return SCIP_OKAY;
}

/** processes the event that a new primal solution has been found */
static
SCIP_DECL_EVENTEXEC(processNewSolutionEvent)
{
   SCIP_CONSHDLR* conshdlr;
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_SOL* sol;

   assert(scip != NULL);
   assert(event != NULL);
   assert(eventdata != NULL);
   assert(eventhdlr != NULL);
   assert(SCIPeventGetType(event) & SCIP_EVENTTYPE_SOLFOUND);

   conshdlr = (SCIP_CONSHDLR*)eventdata;

   if( SCIPconshdlrGetNConss(conshdlr) == 0 )
      return SCIP_OKAY;

   sol = SCIPeventGetSol(event);
   assert(sol != NULL);

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);

   /* we are only interested in solution coming from some heuristic other than trysol, but not from the tree
    * the reason for ignoring trysol solutions is that they may come ~~from an NLP solve in sepalp, where we already added linearizations, or are~~
    * from the tree, but postprocessed via proposeFeasibleSolution
    */
   if( SCIPsolGetHeur(sol) == NULL || SCIPsolGetHeur(sol) == conshdlrdata->trysolheur )
      return SCIP_OKAY;

   SCIPdebugMsg(scip, "caught new sol event %" SCIP_EVENTTYPE_FORMAT " from heur <%s>\n", SCIPeventGetType(event), SCIPheurGetName(SCIPsolGetHeur(sol)));

   SCIP_CALL( addTightEstimatorCuts(scip, conshdlr, SCIPconshdlrGetConss(conshdlr), SCIPconshdlrGetNConss(conshdlr), sol) );

   return SCIP_OKAY;
}

/** tightens the bounds of the auxiliary variable associated with an expression (or original variable if being a variable-expression) according to given bounds
 *
 *  The given bounds may very well be the exprs activity (when called from forwardPropExpr()), but can also be some
 *  tighter bounds (when called from SCIPtightenExprIntervalNonlinear()).
 *
 *  Nothing will happen if SCIP is not in presolve or solve.
 */
static
SCIP_RETCODE tightenAuxVarBounds(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_EXPR*            expr,               /**< expression whose auxvar is to be tightened */
   SCIP_INTERVAL         bounds,             /**< bounds to be used for tightening (must not be empty) */
   SCIP_Bool*            cutoff,             /**< buffer to store whether a cutoff was detected */
   int*                  ntightenings        /**< buffer to add the total number of tightenings, or NULL */
   )
{
   SCIP_VAR* var;
   SCIP_Bool tightenedlb;
   SCIP_Bool tightenedub;
   SCIP_Bool force;

   assert(scip != NULL);
   assert(conshdlr != NULL);
   assert(expr != NULL);
   assert(cutoff != NULL);

   /* the given bounds must not be empty (we could cope, but we shouldn't be called in this situation) */
   assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, bounds));

   *cutoff = FALSE;

   /* do not tighten variable in problem stage (important for unittests)
    * TODO put some kind of #ifdef UNITTEST around this
    */
   if( SCIPgetStage(scip) < SCIP_STAGE_INITPRESOLVE && SCIPgetStage(scip) > SCIP_STAGE_SOLVING )
      return SCIP_OKAY;

   var = SCIPgetExprAuxVarNonlinear(expr);
   if( var == NULL )
      return SCIP_OKAY;

   /* force tightening if conshdlrdata says so or it would mean fixing the variable */
   force = SCIPconshdlrGetData(conshdlr)->forceboundtightening || SCIPisEQ(scip, bounds.inf, bounds.sup);

   /* try to tighten lower bound of (auxiliary) variable */
   SCIP_CALL( SCIPtightenVarLb(scip, var, bounds.inf, force, cutoff, &tightenedlb) );
   if( tightenedlb )
   {
      if( ntightenings != NULL )
         ++*ntightenings;
      SCIPdebugMsg(scip, "tightened lb on auxvar <%s> to %.15g (forced:%u)\n", SCIPvarGetName(var), SCIPvarGetLbLocal(var), force);
   }
   if( *cutoff )
   {
      SCIPdebugMsg(scip, "cutoff when tightening lb on auxvar <%s> to %.15g\n", SCIPvarGetName(var), bounds.inf);
      return SCIP_OKAY;
   }

   /* try to tighten upper bound of (auxiliary) variable */
   SCIP_CALL( SCIPtightenVarUb(scip, var, bounds.sup, force, cutoff, &tightenedub) );
   if( tightenedub )
   {
      if( ntightenings != NULL )
         ++*ntightenings;
      SCIPdebugMsg(scip, "tightened ub on auxvar <%s> to %.15g (forced:%u)\n", SCIPvarGetName(var), SCIPvarGetUbLocal(var), force);
   }
   if( *cutoff )
   {
      SCIPdebugMsg(scip, "cutoff when tightening ub on auxvar <%s> to %.15g\n", SCIPvarGetName(var), bounds.sup);
      return SCIP_OKAY;
   }

   /* TODO expr->activity should have been reevaluated now due to boundchange-events, but it used to relax bounds
    * that seems unnecessary and we could easily undo this here, e.g.,
    * if( tightenedlb ) expr->activity.inf = bounds.inf
    */

   return SCIP_OKAY;
}

/** propagate bounds of the expressions in a given expression tree (that is, updates activity intervals)
 *  and tries to tighten the bounds of the auxiliary variables accordingly
 */
static
SCIP_RETCODE forwardPropExpr(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_EXPR*            rootexpr,           /**< expression */
   SCIP_Bool             tightenauxvars,     /**< should the bounds of auxiliary variables be tightened? */
   SCIP_Bool*            infeasible,         /**< buffer to store whether the problem is infeasible (NULL if not needed) */
   int*                  ntightenings        /**< buffer to store the number of auxiliary variable tightenings (NULL if not needed) */
   )
{
   SCIP_EXPRITER* it;
   SCIP_EXPR* expr;
   SCIP_EXPR_OWNERDATA* ownerdata;
   SCIP_CONSHDLRDATA* conshdlrdata;

   assert(scip != NULL);
   assert(rootexpr != NULL);

   if( infeasible != NULL )
      *infeasible = FALSE;
   if( ntightenings != NULL )
      *ntightenings = 0;

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);

   /* if value is valid and empty, then we cannot improve, so do nothing */
   if( SCIPexprGetActivityTag(rootexpr) >= conshdlrdata->lastboundrelax && SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(rootexpr)) )
   {
      SCIPdebugMsg(scip, "stored activity of root expr is empty and valid (activitytag >= lastboundrelax (%" SCIP_LONGINT_FORMAT ")), skip forwardPropExpr -> cutoff\n", conshdlrdata->lastboundrelax);

      if( infeasible != NULL )
         *infeasible = TRUE;

      /* just update tag to curboundstag */
      SCIPexprSetActivity(rootexpr, SCIPexprGetActivity(rootexpr), conshdlrdata->curboundstag);

      return SCIP_OKAY;
   }

   /* if value is up-to-date, then nothing to do */
   if( SCIPexprGetActivityTag(rootexpr) == conshdlrdata->curboundstag )
   {
      SCIPdebugMsg(scip, "activitytag of root expr equals curboundstag (%" SCIP_LONGINT_FORMAT "), skip forwardPropExpr\n", conshdlrdata->curboundstag);

      assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(rootexpr))); /* handled in previous if() */

      return SCIP_OKAY;
   }

   ownerdata = SCIPexprGetOwnerData(rootexpr);
   assert(ownerdata != NULL);

   /* if activity of rootexpr is not used, but expr participated in detect (nenfos >= 0), then we do nothing
    * it seems wrong to be called for such an expression (unless we are in detect at the moment), so I add a SCIPABORT()
    * during detect, we are in some in-between state where we may want to eval activity
    * on exprs that we did not notify about their activity usage
    */
   if( ownerdata->nenfos >= 0 && ownerdata->nactivityusesprop == 0 && ownerdata->nactivityusessepa == 0 && !conshdlrdata->indetect)
   {
#ifdef DEBUG_PROP
      SCIPdebugMsg(scip, "root expr activity is not used but enfo initialized, skip inteval\n");
#endif
      SCIPABORT();
      return SCIP_OKAY;
   }

   SCIP_CALL( SCIPcreateExpriter(scip, &it) );
   SCIP_CALL( SCIPexpriterInit(it, rootexpr, SCIP_EXPRITER_DFS, TRUE) );
   SCIPexpriterSetStagesDFS(it, SCIP_EXPRITER_VISITINGCHILD | SCIP_EXPRITER_LEAVEEXPR);

   for( expr = SCIPexpriterGetCurrent(it); !SCIPexpriterIsEnd(it);  )
   {
      switch( SCIPexpriterGetStageDFS(it) )
      {
         case SCIP_EXPRITER_VISITINGCHILD :
         {
            /* skip child if it has been evaluated already */
            SCIP_EXPR* child;

            child = SCIPexpriterGetChildExprDFS(it);
            if( conshdlrdata->curboundstag == SCIPexprGetActivityTag(child) )
            {
               if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(child)) && infeasible != NULL )
                  *infeasible = TRUE;

               expr = SCIPexpriterSkipDFS(it);
               continue;
            }

            break;
         }

         case SCIP_EXPRITER_LEAVEEXPR :
         {
            SCIP_INTERVAL activity;

            /* we should not have entered this expression if its activity was already up to date */
            assert(SCIPexprGetActivityTag(expr) < conshdlrdata->curboundstag);

            ownerdata = SCIPexprGetOwnerData(expr);
            assert(ownerdata != NULL);

            /* for var exprs where varevents are catched, activity is updated immediately when the varbound has been changed
             * so we can assume that the activity is up to date for all these variables
             * UNLESS we changed the method used to evaluate activity of variable expressions
             *   or we currently use global bounds (varevents are catched for local bound changes only)
             */
            if( SCIPisExprVar(scip, expr) && ownerdata->filterpos >= 0 &&
                SCIPexprGetActivityTag(expr) >= conshdlrdata->lastvaractivitymethodchange && !conshdlrdata->globalbounds )
            {
#ifndef NDEBUG
               SCIP_INTERVAL exprhdlrinterval;

               SCIP_CALL( SCIPcallExprInteval(scip, expr, &exprhdlrinterval, conshdlrdata->intevalvar, conshdlrdata) );
               assert(SCIPisRelEQ(scip, exprhdlrinterval.inf, SCIPexprGetActivity(expr).inf));
               assert(SCIPisRelEQ(scip, exprhdlrinterval.sup, SCIPexprGetActivity(expr).sup));
#endif
#ifdef DEBUG_PROP
               SCIPdebugMsg(scip, "skip interval evaluation of expr for var <%s> [%g,%g]\n", SCIPvarGetName(SCIPgetVarExprVar(expr)), SCIPexprGetActivity(expr).inf, SCIPexprGetActivity(expr).sup);
#endif
               SCIPexprSetActivity(expr, SCIPexprGetActivity(expr), conshdlrdata->curboundstag);

               break;
            }

            if( SCIPexprGetActivityTag(expr) < conshdlrdata->lastboundrelax )
            {
               /* start with entire activity if current one is invalid */
               SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &activity);
            }
            else if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(expr)) )
            {
               /* If already empty, then don't try to compute even better activity.
                * If cons_nonlinear were alone, then we should have noted that we are infeasible
                * so an assert(infeasible == NULL || *infeasible) should work here.
                * However, after reporting a cutoff due to expr->activity being empty,
                * SCIP may wander to a different node and call propagation again.
                * If no bounds in a nonlinear constraint have been relaxed when switching nodes
                * (so expr->activitytag >= conshdlrdata->lastboundrelax), then
                * we will still have expr->activity being empty, but will have forgotten
                * that we found infeasibility here before (!2221#note_134120).
                * Therefore we just set *infeasibility=TRUE here and stop.
                */
               if( infeasible != NULL )
                  *infeasible = TRUE;
               SCIPdebugMsg(scip, "expr %p already has empty activity -> cutoff\n", (void*)expr);
               break;
            }
            else
            {
               /* start with current activity, since it is valid */
               activity = SCIPexprGetActivity(expr);
            }

            /* if activity of expr is not used, but expr participated in detect (nenfos >= 0), then do nothing */
            if( ownerdata->nenfos >= 0 && ownerdata->nactivityusesprop == 0 && ownerdata->nactivityusessepa == 0 && !conshdlrdata->indetect )
            {
#ifdef DEBUG_PROP
               SCIPdebugMsg(scip, "expr %p activity is not used but enfo initialized, skip inteval\n", (void*)expr);
#endif
               break;
            }

#ifdef DEBUG_PROP
            SCIPdebugMsg(scip, "interval evaluation of expr %p ", (void*)expr);
            SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
            SCIPdebugMsgPrint(scip, ", current activity = [%.20g, %.20g]\n", SCIPexprGetActivity(expr).inf, SCIPexprGetActivity(expr).sup);
#endif

            /* run interval eval of nonlinear handlers or expression handler */
            if( ownerdata->nenfos > 0 )
            {
               SCIP_NLHDLR* nlhdlr;
               SCIP_INTERVAL nlhdlrinterval;
               int e;

               /* for expressions with enforcement, nlhdlrs take care of interval evaluation */
               for( e = 0; e < ownerdata->nenfos && !SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, activity); ++e )
               {
                  /* skip nlhdlr if it does not want to participate in activity computation */
                  if( (ownerdata->enfos[e]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_ACTIVITY) == 0 )
                     continue;

                  nlhdlr = ownerdata->enfos[e]->nlhdlr;
                  assert(nlhdlr != NULL);

                  /* skip nlhdlr if it does not provide interval evaluation (so it may only provide reverse propagation) */
                  if( !SCIPnlhdlrHasIntEval(nlhdlr) )
                     continue;

                  /* let nlhdlr evaluate current expression */
                  nlhdlrinterval = activity;
                  SCIP_CALL( SCIPnlhdlrInteval(scip, nlhdlr, expr, ownerdata->enfos[e]->nlhdlrexprdata,
                     &nlhdlrinterval, conshdlrdata->intevalvar, conshdlrdata) );
#ifdef DEBUG_PROP
                  SCIPdebugMsg(scip, " nlhdlr <%s>::inteval = [%.20g, %.20g]", SCIPnlhdlrGetName(nlhdlr), nlhdlrinterval.inf, nlhdlrinterval.sup);
#endif

                  /* update activity by intersecting with computed activity */
                  SCIPintervalIntersectEps(&activity, SCIPepsilon(scip), activity, nlhdlrinterval);
#ifdef DEBUG_PROP
                  SCIPdebugMsgPrint(scip, " -> new activity: [%.20g, %.20g]\n", activity.inf, activity.sup);
#endif
               }
            }
            else
            {
               /* for node without enforcement (before or during detect), call the callback of the exprhdlr directly */
               SCIP_INTERVAL exprhdlrinterval = activity;
               SCIP_CALL( SCIPcallExprInteval(scip, expr, &exprhdlrinterval, conshdlrdata->intevalvar, conshdlrdata) );
#ifdef DEBUG_PROP
               SCIPdebugMsg(scip, " exprhdlr <%s>::inteval = [%.20g, %.20g]", SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), exprhdlrinterval.inf, exprhdlrinterval.sup);
#endif

               /* update expr->activity by intersecting with computed activity */
               SCIPintervalIntersectEps(&activity, SCIPepsilon(scip), activity, exprhdlrinterval);
#ifdef DEBUG_PROP
               SCIPdebugMsgPrint(scip, " -> new activity: [%.20g, %.20g]\n", activity.inf, activity.sup);
#endif
            }

            /* if expression is integral, then we try to tighten the interval bounds a bit
             * this should undo the addition of some unnecessary safety added by use of nextafter() in interval arithmetics, e.g., when doing pow()
             * it would be ok to use ceil() and floor(), but for safety we use SCIPceil and SCIPfloor for now
             * do this only if using boundtightening-inteval and not in redundancy check (there we really want to relax all variables)
             * boundtightening-inteval does not relax integer variables, so can omit expressions without children
             * (constants should be ok, too)
             */
            if( SCIPexprIsIntegral(expr) && conshdlrdata->intevalvar == intEvalVarBoundTightening && SCIPexprGetNChildren(expr) > 0 )
            {
               if( activity.inf > -SCIP_INTERVAL_INFINITY )
                  activity.inf = SCIPceil(scip, activity.inf);
               if( activity.sup <  SCIP_INTERVAL_INFINITY )
                  activity.sup = SCIPfloor(scip, activity.sup);
#ifdef DEBUG_PROP
               SCIPdebugMsg(scip, " applying integrality: [%.20g, %.20g]\n", activity.inf, activity.sup);
#endif
            }

            /* mark the current node to be infeasible if either the lower/upper bound is above/below +/- SCIPinfinity()
             * TODO this is a problem if dual-presolve fixed a variable to +/- infinity
             */
            if( SCIPisInfinity(scip, activity.inf) || SCIPisInfinity(scip, -activity.sup) )
            {
               SCIPdebugMsg(scip, "cut off due to activity [%g,%g] beyond infinity\n", activity.inf, activity.sup);
               SCIPintervalSetEmpty(&activity);
            }

            /* now finally store activity in expr */
            SCIPexprSetActivity(expr, activity, conshdlrdata->curboundstag);

            if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, activity) )
            {
               if( infeasible != NULL )
                  *infeasible = TRUE;
            }
            else if( tightenauxvars && ownerdata->auxvar != NULL )
            {
               SCIP_Bool tighteninfeasible;

               SCIP_CALL( tightenAuxVarBounds(scip, conshdlr, expr, activity, &tighteninfeasible, ntightenings) );
               if( tighteninfeasible )
               {
                  if( infeasible != NULL )
                     *infeasible = TRUE;
                  SCIPintervalSetEmpty(&activity);
                  SCIPexprSetActivity(expr, activity, conshdlrdata->curboundstag);
               }
            }

            break;
         }

         default:
            /* you should never be here */
            SCIPerrorMessage("unexpected iterator stage\n");
            SCIPABORT();
            break;
      }

      expr = SCIPexpriterGetNext(it);
   }

   SCIPfreeExpriter(&it);

   return SCIP_OKAY;
}

/** returns whether intersecting `oldinterval` with `newinterval` would provide a properly smaller interval
 *
 * If `subsetsufficient` is TRUE, then the intersection being smaller than oldinterval is sufficient.
 *
 * If `subsetsufficient` is FALSE, then we require
 *  - a change from an unbounded interval to a bounded one, or
 *  - or a change from an unfixed (width > epsilon) to a fixed interval, or
 *  - a minimal tightening of one of the interval bounds as defined by SCIPis{Lb,Ub}Better().
 */
static
SCIP_Bool isIntervalBetter(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_Bool             subsetsufficient,   /**< whether the intersection being a proper subset of oldinterval is sufficient */
   SCIP_INTERVAL         newinterval,        /**< new interval */
   SCIP_INTERVAL         oldinterval         /**< old interval */
   )
{
   assert(scip != NULL);
   assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, newinterval));
   assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, oldinterval));

   if( subsetsufficient )
      /* oldinterval \cap newinterval < oldinterval iff not oldinterval is subset of newinterval */
      return !SCIPintervalIsSubsetEQ(SCIP_INTERVAL_INFINITY, oldinterval, newinterval);

   /* check whether lower bound of interval becomes finite */
   if( oldinterval.inf <= -SCIP_INTERVAL_INFINITY && newinterval.inf > -SCIP_INTERVAL_INFINITY )
      return TRUE;

   /* check whether upper bound of interval becomes finite */
   if( oldinterval.sup >=  SCIP_INTERVAL_INFINITY && newinterval.sup >  SCIP_INTERVAL_INFINITY )
      return TRUE;

   /* check whether intersection will have width <= epsilon, if oldinterval doesn't have yet */
   if( !SCIPisEQ(scip, oldinterval.inf, oldinterval.sup) && SCIPisEQ(scip, MAX(oldinterval.inf, newinterval.inf), MIN(oldinterval.sup, newinterval.sup)) )
      return TRUE;

   /* check whether lower bound on interval will be better by SCIP's quality measures for boundchanges */
   if( SCIPisLbBetter(scip, newinterval.inf, oldinterval.inf, oldinterval.sup) )
      return TRUE;

   /* check whether upper bound on interval will be better by SCIP's quality measures for boundchanges */
   if( SCIPisUbBetter(scip, newinterval.sup, oldinterval.inf, oldinterval.sup) )
      return TRUE;

   return FALSE;
}

/** propagates bounds for each sub-expression in the `reversepropqueue` by starting from the root expressions
 *
 *  The expression will be traversed in breadth first search by using this queue.
 *
 *  @note Calling this function requires feasible intervals for each sub-expression; this is guaranteed by calling
 *  forwardPropExpr() before calling this function.
 *
 *  @note Calling this function with `*infeasible` = TRUE will only empty the queue.
 */
static
SCIP_RETCODE reversePropQueue(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_Bool*            infeasible,         /**< buffer to update whether an expression's bounds were propagated to an empty interval */
   int*                  ntightenings        /**< buffer to store the number of (variable) tightenings */
   )
{
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_EXPR* expr;
   SCIP_EXPR_OWNERDATA* ownerdata;

   assert(infeasible != NULL);
   assert(ntightenings != NULL);

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);

   *ntightenings = 0;

   /* main loop that calls reverse propagation for expressions on the queue
    * when reverseprop finds a tightening for an expression, then that expression is added to the queue (within the reverseprop call)
    */
   while( !SCIPqueueIsEmpty(conshdlrdata->reversepropqueue) && !(*infeasible) )
   {
      SCIP_INTERVAL propbounds;
      int e;

      expr = (SCIP_EXPR*) SCIPqueueRemove(conshdlrdata->reversepropqueue);
      assert(expr != NULL);

      ownerdata = SCIPexprGetOwnerData(expr);
      assert(ownerdata != NULL);

      assert(ownerdata->inpropqueue);
      /* mark that the expression is not in the queue anymore */
      ownerdata->inpropqueue = FALSE;

      /* since the expr was in the propagation queue, the propbounds should belong to current propagation and should not be empty
       * (propbounds being entire doesn't make much sense, so assert this for now, too, but that could be removed)
       */
      assert(ownerdata->propboundstag == conshdlrdata->curpropboundstag);
      assert(!SCIPintervalIsEntire(SCIP_INTERVAL_INFINITY, ownerdata->propbounds));
      assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, ownerdata->propbounds));

      /* this intersects propbounds with activity and auxvar bounds
       * I doubt this would be much helpful, since propbounds are already subset of activity and we also propagate
       * auxvar bounds separately, so disabling this for now
       */
#ifdef SCIP_DISABLED_CODE
      propbounds = SCIPgetExprBoundsNonlinear(scip, expr);
      if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, propbounds) )
      {
         *infeasible = TRUE;
         break;
      }
#else
      propbounds = ownerdata->propbounds;
#endif

      if( ownerdata->nenfos > 0 )
      {
         /* for nodes with enforcement, call reverse propagation callbacks of nlhdlrs */
         for( e = 0; e < ownerdata->nenfos && !*infeasible; ++e )
         {
            SCIP_NLHDLR* nlhdlr;
            int nreds;

            /* skip nlhdlr if it does not want to participate in activity computation */
            if( (ownerdata->enfos[e]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_ACTIVITY) == 0 )
               continue;

            nlhdlr = ownerdata->enfos[e]->nlhdlr;
            assert(nlhdlr != NULL);

            /* call the reverseprop of the nlhdlr */
#ifdef SCIP_DEBUG
            SCIPdebugMsg(scip, "call reverse propagation for ");
            SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
            SCIPdebugMsgPrint(scip, " in [%g,%g] using nlhdlr <%s>\n", propbounds.inf, propbounds.sup, SCIPnlhdlrGetName(nlhdlr));
#endif

            nreds = 0;
            SCIP_CALL( SCIPnlhdlrReverseprop(scip, conshdlr, nlhdlr, expr, ownerdata->enfos[e]->nlhdlrexprdata, propbounds, infeasible, &nreds) );
            assert(nreds >= 0);
            *ntightenings += nreds;
         }
      }
      else if( SCIPexprhdlrHasReverseProp(SCIPexprGetHdlr(expr)) )
      {
         /* if expr without enforcement (before detect), call reverse propagation callback of exprhdlr directly */
         SCIP_INTERVAL* childrenbounds;
         int c;

#ifdef SCIP_DEBUG
         SCIPdebugMsg(scip, "call reverse propagation for ");
         SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
         SCIPdebugMsgPrint(scip, " in [%g,%g] using exprhdlr <%s>\n", SCIPexprGetActivity(expr).inf, SCIPexprGetActivity(expr).sup, SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)));
#endif

         /* if someone added an expr without nlhdlr into the reversepropqueue, then this must be because its enfo hasn't
          * been initialized in detectNlhdlr yet (nenfos < 0)
          */
         assert(ownerdata->nenfos < 0);

         SCIP_CALL( SCIPallocBufferArray(scip, &childrenbounds, SCIPexprGetNChildren(expr)) );
         for( c = 0; c < SCIPexprGetNChildren(expr); ++c )
            childrenbounds[c] = SCIPgetExprBoundsNonlinear(scip, SCIPexprGetChildren(expr)[c]);

         /* call the reverseprop of the exprhdlr */
         SCIP_CALL( SCIPcallExprReverseprop(scip, expr, propbounds, childrenbounds, infeasible) );

         if( !*infeasible )
            for( c = 0; c < SCIPexprGetNChildren(expr); ++c )
            {
               SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, SCIPexprGetChildren(expr)[c], childrenbounds[c], infeasible, ntightenings) );
            }

         SCIPfreeBufferArray(scip, &childrenbounds);
      }
   }

   /* reset inpropqueue for all remaining expr's in queue (can happen in case of early stop due to infeasibility) */
   while( !SCIPqueueIsEmpty(conshdlrdata->reversepropqueue) )
   {
      expr = (SCIP_EXPR*) SCIPqueueRemove(conshdlrdata->reversepropqueue);
      assert(expr != NULL);

      ownerdata = SCIPexprGetOwnerData(expr);
      assert(ownerdata != NULL);

      /* mark that the expression is not in the queue anymore */
      ownerdata->inpropqueue = FALSE;
   }

   return SCIP_OKAY;
}

/** calls domain propagation for a given set of constraints
 *
 *  The algorithm alternates calls of forward and reverse propagation.
 *  Forward propagation ensures that activity of expressions is up to date.
 *  Reverse propagation tries to derive tighter variable bounds by reversing the activity computation, using the constraints
 *  [lhs,rhs] interval as starting point.
 *
 *  The propagation algorithm works as follows:
 *   1. apply forward propagation (update activities) for all constraints not marked as propagated
 *   2. if presolve or propauxvars is disabled: collect expressions for which the constraint sides provide tighter bounds
 *      if solve and propauxvars is enabled: collect expressions for which auxvars (including those in root exprs)
 *      provide tighter bounds
 *   3. apply reverse propagation to all collected expressions; don't explore
 *      sub-expressions which have not changed since the beginning of the propagation loop
 *   4. if we have found enough tightenings go to 1, otherwise leave propagation loop
 *
 *  @note After calling forward propagation for a constraint, we mark this constraint as propagated. This flag might be
 *  reset during the reverse propagation when we find a bound tightening of a variable expression contained in the
 *  constraint. Resetting this flag is done in the EVENTEXEC callback of the event handler
 *
 *  TODO should we distinguish between expressions where activity information is used for separation and those where not,
 *    e.g., try less to propagate on convex constraints?
 */
static
SCIP_RETCODE propConss(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS**           conss,              /**< constraints to propagate */
   int                   nconss,             /**< total number of constraints */
   SCIP_Bool             force,              /**< force tightening even if below bound strengthening tolerance */
   SCIP_RESULT*          result,             /**< pointer to store the result */
   int*                  nchgbds             /**< buffer to add the number of changed bounds */
   )
{
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_CONSDATA* consdata;
   SCIP_EXPR_OWNERDATA* ownerdata;
   SCIP_Bool cutoff = FALSE;
   SCIP_INTERVAL conssides;
   int ntightenings;
   int roundnr;
   SCIP_EXPRITER* revpropcollectit = NULL;
   int i;

   assert(scip != NULL);
   assert(conshdlr != NULL);
   assert(conss != NULL);
   assert(nconss >= 0);
   assert(result != NULL);
   assert(nchgbds != NULL);
   assert(*nchgbds >= 0);

   /* no constraints to propagate */
   if( nconss == 0 )
   {
      *result = SCIP_DIDNOTRUN;
      return SCIP_OKAY;
   }

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);
   assert(conshdlrdata->intevalvar == intEvalVarBoundTightening);
   assert(!conshdlrdata->globalbounds);

   *result = SCIP_DIDNOTFIND;
   roundnr = 0;

   /* tightenAuxVarBounds() needs to know whether boundtightenings are to be forced */
   conshdlrdata->forceboundtightening = force;

   /* invalidate all propbounds (probably not needed) */
   ++conshdlrdata->curpropboundstag;

   /* create iterator that we will use if we need to look at all auxvars */
   if( conshdlrdata->propauxvars )
   {
      SCIP_CALL( SCIPcreateExpriter(scip, &revpropcollectit) );
   }

   /* main propagation loop */
   do
   {
      SCIPdebugMsg(scip, "start propagation round %d\n", roundnr);

      assert(SCIPqueueIsEmpty(conshdlrdata->reversepropqueue));

      /* apply forward propagation (update expression activities)
       * and add promising root expressions into queue for reversepropagation
       */
      for( i = 0; i < nconss; ++i )
      {
         consdata = SCIPconsGetData(conss[i]);
         assert(consdata != NULL);

         /* skip deleted, non-active, or propagation-disabled constraints */
         if( SCIPconsIsDeleted(conss[i]) || !SCIPconsIsActive(conss[i]) || !SCIPconsIsPropagationEnabled(conss[i]) )
            continue;

         /* skip already propagated constraints, i.e., constraints where no (original) variable has changed and thus
          * activity didn't change
          */
         if( consdata->ispropagated )
            continue;

         /* update activities in expression */
         SCIPdebugMsg(scip, "call forwardPropExpr() for constraint <%s> (round %d): ", SCIPconsGetName(conss[i]), roundnr);
         SCIPdebugPrintCons(scip, conss[i], NULL);

         ntightenings = 0;
         SCIP_CALL( forwardPropExpr(scip, conshdlr, consdata->expr, TRUE, &cutoff, &ntightenings) );
         assert(cutoff || !SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(consdata->expr)));

         if( cutoff )
         {
            SCIPdebugMsg(scip, " -> cutoff in forwardPropExpr (due to domain error or auxvar tightening) of constraint <%s>\n", SCIPconsGetName(conss[i]));
            *result = SCIP_CUTOFF;
            break;
         }

         ownerdata = SCIPexprGetOwnerData(consdata->expr);

         /* TODO for a constraint that only has an auxvar for consdata->expr (e.g., convex quadratic), we could also just do the if(TRUE)-branch */
         if( !conshdlrdata->propauxvars || ownerdata->auxvar == NULL )
         {
            /* check whether constraint sides (relaxed by epsilon) or auxvar bounds provide a tightening
             *   (if we have auxvar (not in presolve), then bounds of the auxvar are initially set to constraint sides,
             *   so taking auxvar bounds is enough)
             */
            if( ownerdata->auxvar == NULL )
            {
               /* relax sides by SCIPepsilon() and handle infinite sides */
               SCIP_Real lhs = SCIPisInfinity(scip, -consdata->lhs) ? -SCIP_INTERVAL_INFINITY : consdata->lhs - conshdlrdata->conssiderelaxamount;
               SCIP_Real rhs = SCIPisInfinity(scip,  consdata->rhs) ?  SCIP_INTERVAL_INFINITY : consdata->rhs + conshdlrdata->conssiderelaxamount;
               SCIPintervalSetBounds(&conssides, lhs, rhs);
            }
            else
            {
               conssides = intEvalVarBoundTightening(scip, ownerdata->auxvar, (void*)conshdlrdata);
            }
            SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, consdata->expr, conssides, &cutoff, &ntightenings) );
         }
         else
         {
            /* check whether bounds of any auxvar used in constraint provides a tightening
             *   (for the root expression, bounds of auxvar are initially set to constraint sides)
             * but skip exprs that have an auxvar, but do not participate in propagation
             */
            SCIP_EXPR* expr;

            assert(revpropcollectit != NULL);
            SCIP_CALL( SCIPexpriterInit(revpropcollectit, consdata->expr, SCIP_EXPRITER_BFS, FALSE) );
            for( expr = SCIPexpriterGetCurrent(revpropcollectit); !SCIPexpriterIsEnd(revpropcollectit) && !cutoff; expr = SCIPexpriterGetNext(revpropcollectit) )
            {
               ownerdata = SCIPexprGetOwnerData(expr);
               assert(ownerdata != NULL);

               if( ownerdata->auxvar == NULL )
                  continue;

               if( ownerdata->nactivityusesprop == 0 && ownerdata->nactivityusessepa == 0 )
                  continue;

               conssides = intEvalVarBoundTightening(scip, ownerdata->auxvar, (void*)conshdlrdata);
               SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, expr, conssides, &cutoff, &ntightenings) );
            }
         }

         if( cutoff )
         {
            SCIPdebugMsg(scip, " -> cutoff after intersect with conssides of constraint <%s>\n", SCIPconsGetName(conss[i]));
            *result = SCIP_CUTOFF;
            break;
         }

         assert(ntightenings >= 0);
         if( ntightenings > 0 )
         {
            *nchgbds += ntightenings;
            *result = SCIP_REDUCEDDOM;
         }

         /* mark constraint as propagated; this will be reset via the event system when we find a variable tightening */
         consdata->ispropagated = TRUE;
      }

      /* apply backward propagation (if cutoff is TRUE, then this call empties the queue) */
      SCIP_CALL( reversePropQueue(scip, conshdlr, &cutoff, &ntightenings) );
      assert(ntightenings >= 0);
      assert(SCIPqueueIsEmpty(conshdlrdata->reversepropqueue));

      if( cutoff )
      {
         SCIPdebugMsg(scip, " -> cutoff\n");
         *result = SCIP_CUTOFF;
         break;
      }

      if( ntightenings > 0 )
      {
         *nchgbds += ntightenings;
         *result = SCIP_REDUCEDDOM;
      }
   }
   while( ntightenings > 0 && ++roundnr < conshdlrdata->maxproprounds );

   if( conshdlrdata->propauxvars )
   {
      SCIPfreeExpriter(&revpropcollectit);
   }

   conshdlrdata->forceboundtightening = FALSE;

   /* invalidate propbounds in all exprs, so noone accidentally uses them outside propagation */
   ++conshdlrdata->curpropboundstag;

   return SCIP_OKAY;
}

/** calls the reverseprop callbacks of all nlhdlrs in all expressions in all constraints using activity as bounds
 *
 * This is meant to propagate any domain restrictions on functions onto variable bounds, if possible.
 *
 * Assumes that activities are still valid and curpropboundstag does not need to be increased.
 * Therefore, a good place to call this function is immediately after propConss() or after forwardPropExpr() if outside propagation.
 */
static
SCIP_RETCODE propExprDomains(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS**           conss,              /**< constraints to propagate */
   int                   nconss,             /**< total number of constraints */
   SCIP_RESULT*          result,             /**< pointer to store the result */
   int*                  nchgbds             /**< buffer to add the number of changed bounds */
   )
{
   SCIP_CONSDATA* consdata;
   SCIP_EXPRITER* it;
   SCIP_EXPR* expr;
   SCIP_EXPR_OWNERDATA* ownerdata;
   SCIP_Bool cutoff = FALSE;
   int ntightenings;
   int c;
   int e;

   assert(scip != NULL);
   assert(conshdlr != NULL);
   assert(conss != NULL);
   assert(nconss >= 0);
   assert(result != NULL);
   assert(nchgbds != NULL);
   assert(*nchgbds >= 0);

   assert(SCIPconshdlrGetData(conshdlr)->intevalvar == intEvalVarBoundTightening);
   assert(!SCIPconshdlrGetData(conshdlr)->globalbounds);
   assert(SCIPqueueIsEmpty(SCIPconshdlrGetData(conshdlr)->reversepropqueue));

   *result = SCIP_DIDNOTFIND;

   SCIP_CALL( SCIPcreateExpriter(scip, &it) );
   SCIP_CALL( SCIPexpriterInit(it, NULL, SCIP_EXPRITER_DFS, FALSE) );

   for( c = 0; c < nconss && !cutoff; ++c )
   {
      /* skip deleted, non-active, or propagation-disabled constraints */
      if( SCIPconsIsDeleted(conss[c]) || !SCIPconsIsActive(conss[c]) || !SCIPconsIsPropagationEnabled(conss[c]) )
         continue;

      consdata = SCIPconsGetData(conss[c]);
      assert(consdata != NULL);

      for( expr = SCIPexpriterRestartDFS(it, consdata->expr); !SCIPexpriterIsEnd(it) && !cutoff; expr = SCIPexpriterGetNext(it) )
      {
         ownerdata = SCIPexprGetOwnerData(expr);
         assert(ownerdata != NULL);

         /* call reverseprop for those nlhdlr that participate in this expr's activity computation
          * this will propagate the current activity
          */
         for( e = 0; e < ownerdata->nenfos; ++e )
         {
            SCIP_NLHDLR* nlhdlr;
            assert(ownerdata->enfos[e] != NULL);

            nlhdlr = ownerdata->enfos[e]->nlhdlr;
            assert(nlhdlr != NULL);
            if( (ownerdata->enfos[e]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_ACTIVITY) == 0 )
               continue;

            SCIPdebugMsg(scip, "propExprDomains calling reverseprop for expression %p [%g,%g]\n", (void*)expr,
                  SCIPexprGetActivity(expr).inf, SCIPexprGetActivity(expr).sup);
            ntightenings = 0;
            SCIP_CALL( SCIPnlhdlrReverseprop(scip, conshdlr, nlhdlr, expr, ownerdata->enfos[e]->nlhdlrexprdata,
                    SCIPexprGetActivity(expr), &cutoff, &ntightenings) );

            if( cutoff )
            {
               /* stop everything if we detected infeasibility */
               SCIPdebugMsg(scip, "detect infeasibility for constraint <%s> during reverseprop()\n", SCIPconsGetName(conss[c]));
               *result = SCIP_CUTOFF;
               break;
            }

            assert(ntightenings >= 0);
            if( ntightenings > 0 )
            {
               *nchgbds += ntightenings;
               *result = SCIP_REDUCEDDOM;
            }
         }
      }
   }

   /* apply backward propagation (if cutoff is TRUE, then this call empties the queue) */
   SCIP_CALL( reversePropQueue(scip, conshdlr, &cutoff, &ntightenings) );
   assert(ntightenings >= 0);

   if( cutoff )
   {
      SCIPdebugMsg(scip, " -> cutoff\n");
      *result = SCIP_CUTOFF;
   }
   else if( ntightenings > 0 )
   {
      *nchgbds += ntightenings;
      *result = SCIP_REDUCEDDOM;
   }

   SCIPfreeExpriter(&it);

   /* invalidate propbounds in all exprs, so noone accidentally uses them outside propagation */
   ++SCIPconshdlrGetData(conshdlr)->curpropboundstag;

   return SCIP_OKAY;
}

/** propagates variable locks through expression and adds locks to variables */
static
SCIP_RETCODE propagateLocks(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EXPR*            expr,               /**< expression */
   int                   nlockspos,          /**< number of positive locks */
   int                   nlocksneg           /**< number of negative locks */
   )
{
   SCIP_EXPR_OWNERDATA* ownerdata;
   SCIP_EXPRITER* it;
   SCIP_EXPRITER_USERDATA ituserdata;

   assert(expr != NULL);

   /* if no locks, then nothing to propagate */
   if( nlockspos == 0 && nlocksneg == 0 )
      return SCIP_OKAY;

   SCIP_CALL( SCIPcreateExpriter(scip, &it) );
   SCIP_CALL( SCIPexpriterInit(it, expr, SCIP_EXPRITER_DFS, TRUE) );
   SCIPexpriterSetStagesDFS(it, SCIP_EXPRITER_ENTEREXPR | SCIP_EXPRITER_VISITINGCHILD | SCIP_EXPRITER_LEAVEEXPR);
   assert(SCIPexpriterGetCurrent(it) == expr); /* iterator should not have moved */

   /* store locks in root node */
   ituserdata.intvals[0] = nlockspos;
   ituserdata.intvals[1] = nlocksneg;
   SCIPexpriterSetCurrentUserData(it, ituserdata);

   while( !SCIPexpriterIsEnd(it) )
   {
      /* collect locks */
      ituserdata = SCIPexpriterGetCurrentUserData(it);
      nlockspos = ituserdata.intvals[0];
      nlocksneg = ituserdata.intvals[1];

      ownerdata = SCIPexprGetOwnerData(expr);

      switch( SCIPexpriterGetStageDFS(it) )
      {
         case SCIP_EXPRITER_ENTEREXPR:
         {
            if( SCIPisExprVar(scip, expr) )
            {
               /* if a variable, then also add nlocksneg/nlockspos via SCIPaddVarLocks() */
               SCIP_CALL( SCIPaddVarLocks(scip, SCIPgetVarExprVar(expr), nlocksneg, nlockspos) );
            }

            /* add locks to expression */
            ownerdata->nlockspos += nlockspos;
            ownerdata->nlocksneg += nlocksneg;

            /* add monotonicity information if expression has been locked for the first time */
            if( ownerdata->nlockspos == nlockspos && ownerdata->nlocksneg == nlocksneg && SCIPexprGetNChildren(expr) > 0
               && SCIPexprhdlrHasMonotonicity(SCIPexprGetHdlr(expr)) )
            {
               int i;

               assert(ownerdata->monotonicity == NULL);
               assert(ownerdata->monotonicitysize == 0);

               SCIP_CALL( SCIPallocBlockMemoryArray(scip, &ownerdata->monotonicity, SCIPexprGetNChildren(expr)) );
               ownerdata->monotonicitysize = SCIPexprGetNChildren(expr);

               /* store the monotonicity for each child */
               for( i = 0; i < SCIPexprGetNChildren(expr); ++i )
               {
                  SCIP_CALL( SCIPcallExprMonotonicity(scip, expr, i, &ownerdata->monotonicity[i]) );
               }
            }
            break;
         }

         case SCIP_EXPRITER_LEAVEEXPR :
         {
            /* remove monotonicity information if expression has been unlocked */
            if( ownerdata->nlockspos == 0 && ownerdata->nlocksneg == 0 && ownerdata->monotonicity != NULL )
            {
               assert(ownerdata->monotonicitysize > 0);
               /* keep this assert for checking whether someone changed an expression without updating locks properly */
               assert(ownerdata->monotonicitysize == SCIPexprGetNChildren(expr));

               SCIPfreeBlockMemoryArray(scip, &ownerdata->monotonicity, ownerdata->monotonicitysize);
               ownerdata->monotonicitysize = 0;
            }
            break;
         }

         case SCIP_EXPRITER_VISITINGCHILD :
         {
            SCIP_MONOTONE monotonicity;

            /* get monotonicity of child */
            /* NOTE: the monotonicity stored in an expression might be different from the result obtained by
             * SCIPcallExprMonotonicity
             */
            monotonicity = ownerdata->monotonicity != NULL ? ownerdata->monotonicity[SCIPexpriterGetChildIdxDFS(it)] : SCIP_MONOTONE_UNKNOWN;

            /* compute resulting locks of the child expression */
            switch( monotonicity )
            {
               case SCIP_MONOTONE_INC:
                  ituserdata.intvals[0] = nlockspos;
                  ituserdata.intvals[1] = nlocksneg;
                  break;
               case SCIP_MONOTONE_DEC:
                  ituserdata.intvals[0] = nlocksneg;
                  ituserdata.intvals[1] = nlockspos;
                  break;
               case SCIP_MONOTONE_UNKNOWN:
                  ituserdata.intvals[0] = nlockspos + nlocksneg;
                  ituserdata.intvals[1] = nlockspos + nlocksneg;
                  break;
               case SCIP_MONOTONE_CONST:
                  ituserdata.intvals[0] = 0;
                  ituserdata.intvals[1] = 0;
                  break;
            }
            /* set locks in child expression */
            SCIPexpriterSetChildUserData(it, ituserdata);

            break;
         }

         default :
            /* you should never be here */
            SCIPABORT();
            break;
      }

      expr = SCIPexpriterGetNext(it);
   }

   SCIPfreeExpriter(&it);

   return SCIP_OKAY;
}

/** main function for adding locks to expressions and variables
 *
 * Locks for a nonlinear constraint are used to update locks for all sub-expressions and variables.
 * Locks of expressions depend on the monotonicity of expressions w.r.t. their children, e.g.,
 * consider the constraint \f$x^2 \leq 1\f$ with \f$x \in [-2,-1]\f$ implies an up-lock for the root
 * expression (pow) and a down-lock for its child \f$x\f$ because \f$x^2\f$ is decreasing on [-2,-1].
 * Since the monotonicity (and thus the locks) might also depend on variable bounds, the function remembers
 * the computed monotonicity information of each expression until all locks of an expression have been removed,
 * which implies that updating the monotonicity information during the next locking of this expression does not
 * break existing locks.
 *
 * @note When modifying the structure of an expression, e.g., during simplification, it is necessary to remove all
 *       locks from an expression and repropagating them after the structural changes have been applied.
 *       Because of existing common sub-expressions, it might be necessary to remove the locks of all constraints
 *       to ensure that an expression is unlocked (see canonicalizeConstraints() for an example)
 */
static
SCIP_RETCODE addLocks(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons,               /**< nonlinear constraint */
   int                   nlockspos,          /**< number of positive rounding locks */
   int                   nlocksneg           /**< number of negative rounding locks */
   )
{
   SCIP_CONSDATA* consdata;

   assert(cons != NULL);

   if( nlockspos == 0 && nlocksneg == 0 )
      return SCIP_OKAY;

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);

   /* no constraint sides -> nothing to lock */
   if( SCIPisInfinity(scip, consdata->rhs) && SCIPisInfinity(scip, -consdata->lhs) )
      return SCIP_OKAY;

   /* remember locks */
   consdata->nlockspos += nlockspos;
   consdata->nlocksneg += nlocksneg;

   assert(consdata->nlockspos >= 0);
   assert(consdata->nlocksneg >= 0);

   /* compute locks for lock propagation */
   if( !SCIPisInfinity(scip, consdata->rhs) && !SCIPisInfinity(scip, -consdata->lhs) )
   {
      SCIP_CALL( propagateLocks(scip, consdata->expr, nlockspos + nlocksneg, nlockspos + nlocksneg));
   }
   else if( !SCIPisInfinity(scip, consdata->rhs) )
   {
      SCIP_CALL( propagateLocks(scip, consdata->expr, nlockspos, nlocksneg));
   }
   else
   {
      assert(!SCIPisInfinity(scip, -consdata->lhs));
      SCIP_CALL( propagateLocks(scip, consdata->expr, nlocksneg, nlockspos));
   }

   return SCIP_OKAY;
}

/** create a nonlinear row representation of a nonlinear constraint and stores them in consdata */
static
SCIP_RETCODE createNlRow(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons                /**< nonlinear constraint */
   )
{
   SCIP_CONSDATA* consdata;

   assert(scip != NULL);
   assert(cons != NULL);

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);
   assert(consdata->expr != NULL);

   if( consdata->nlrow != NULL )
   {
      SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
   }

   /* better curvature info will be set in initSolve() just before nlrow is added to NLP */
   SCIP_CALL( SCIPcreateNlRow(scip, &consdata->nlrow, SCIPconsGetName(cons), 0.0,
         0, NULL, NULL, NULL, consdata->lhs, consdata->rhs, SCIP_EXPRCURV_UNKNOWN) );

   if( SCIPisExprSum(scip, consdata->expr) )
   {
      /* if root is a sum, then split into linear and nonlinear terms */
      SCIP_EXPR* nonlinpart;
      SCIP_EXPR* child;
      SCIP_Real* coefs;
      int i;

      coefs = SCIPgetCoefsExprSum(consdata->expr);

      /* constant term of sum */
      SCIP_CALL( SCIPchgNlRowConstant(scip, consdata->nlrow, SCIPgetConstantExprSum(consdata->expr)) );

      /* a sum-expression that will hold the nonlinear terms and be passed to the nlrow eventually */
      SCIP_CALL( SCIPcreateExprSum(scip, &nonlinpart, 0, NULL, NULL, 0.0, exprownerCreate, (void*)SCIPconsGetHdlr(cons)) );

      for( i = 0; i < SCIPexprGetNChildren(consdata->expr); ++i )
      {
         child = SCIPexprGetChildren(consdata->expr)[i];
         if( SCIPisExprVar(scip, child) )
         {
            /* linear term */
            SCIP_CALL( SCIPaddLinearCoefToNlRow(scip, consdata->nlrow, SCIPgetVarExprVar(child), coefs[i]) );
         }
         else
         {
            /* nonlinear term */
            SCIP_CALL( SCIPappendExprSumExpr(scip, nonlinpart, child, coefs[i]) );
         }
      }

      if( SCIPexprGetNChildren(nonlinpart) > 0 )
      {
         /* add expression to nlrow (this will make a copy) */
         SCIP_CALL( SCIPsetNlRowExpr(scip, consdata->nlrow, nonlinpart) );
      }
      SCIP_CALL( SCIPreleaseExpr(scip, &nonlinpart) );
   }
   else
   {
      SCIP_CALL( SCIPsetNlRowExpr(scip, consdata->nlrow, consdata->expr) );
   }

   return SCIP_OKAY;
}

/** compares enfodata by enforcement priority of nonlinear handler
 *
 * If handlers have same enforcement priority, then compare by detection priority, then by name.
 */
static
SCIP_DECL_SORTPTRCOMP(enfodataCmp)
{
   SCIP_NLHDLR* h1;
   SCIP_NLHDLR* h2;

   assert(elem1 != NULL);
   assert(elem2 != NULL);

   h1 = ((EXPRENFO*)elem1)->nlhdlr;
   h2 = ((EXPRENFO*)elem2)->nlhdlr;

   assert(h1 != NULL);
   assert(h2 != NULL);

   if( SCIPnlhdlrGetEnfoPriority(h1) != SCIPnlhdlrGetEnfoPriority(h2) )
      return SCIPnlhdlrGetEnfoPriority(h1) - SCIPnlhdlrGetEnfoPriority(h2);

   if( SCIPnlhdlrGetDetectPriority(h1) != SCIPnlhdlrGetDetectPriority(h2) )
      return SCIPnlhdlrGetDetectPriority(h1) - SCIPnlhdlrGetDetectPriority(h2);

   return strcmp(SCIPnlhdlrGetName(h1), SCIPnlhdlrGetName(h2));
}

/** install nlhdlrs in one expression */
static
SCIP_RETCODE detectNlhdlr(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EXPR*            expr,               /**< expression for which to run detection routines */
   SCIP_CONS*            cons                /**< constraint for which expr == consdata->expr, otherwise NULL */
   )
{
   SCIP_EXPR_OWNERDATA* ownerdata;
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_NLHDLR_METHOD enforcemethodsallowed;
   SCIP_NLHDLR_METHOD enforcemethods;
   SCIP_NLHDLR_METHOD enforcemethodsnew;
   SCIP_NLHDLR_METHOD nlhdlrenforcemethods;
   SCIP_NLHDLR_METHOD nlhdlrparticipating;
   SCIP_NLHDLREXPRDATA* nlhdlrexprdata;
   int enfossize;  /* allocated length of expr->enfos array */
   int h;

   assert(expr != NULL);

   ownerdata = SCIPexprGetOwnerData(expr);
   assert(ownerdata != NULL);

   conshdlrdata = SCIPconshdlrGetData(ownerdata->conshdlr);
   assert(conshdlrdata != NULL);
   assert(conshdlrdata->auxvarid >= 0);
   assert(!conshdlrdata->indetect);

   /* there should be no enforcer yet and detection should not even have considered expr yet */
   assert(ownerdata->nenfos < 0);
   assert(ownerdata->enfos == NULL);

   /* check which enforcement methods are required by setting flags in enforcemethods for those that are NOT required
    * - if no auxiliary variable is used, then do not need sepabelow or sepaabove
    * - if auxiliary variable is used, but nobody positively (up) locks expr -> only need to enforce expr >= auxvar -> no need for underestimation
    * - if auxiliary variable is used, but nobody negatively (down) locks expr -> only need to enforce expr <= auxvar -> no need for overestimation
    * - if no one uses activity, then do not need activity methods
    */
   enforcemethods = SCIP_NLHDLR_METHOD_NONE;
   if( ownerdata->nauxvaruses == 0 )
      enforcemethods |= SCIP_NLHDLR_METHOD_SEPABOTH;
   else
   {
      if( ownerdata->nlockspos == 0 )  /* no need for underestimation */
         enforcemethods |= SCIP_NLHDLR_METHOD_SEPABELOW;
      if( ownerdata->nlocksneg == 0 )  /* no need for overestimation */
         enforcemethods |= SCIP_NLHDLR_METHOD_SEPAABOVE;
   }
   if( ownerdata->nactivityusesprop == 0 && ownerdata->nactivityusessepa == 0 )
      enforcemethods |= SCIP_NLHDLR_METHOD_ACTIVITY;

   /* it doesn't make sense to have been called on detectNlhdlr, if the expr isn't used for anything */
   assert(enforcemethods != SCIP_NLHDLR_METHOD_ALL);

   /* all methods that have not been flagged above are the ones that we want to be handled by nlhdlrs */
   enforcemethodsallowed = ~enforcemethods & SCIP_NLHDLR_METHOD_ALL;

   ownerdata->nenfos = 0;
   enfossize = 2;
   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &ownerdata->enfos, enfossize) );
   conshdlrdata->indetect = TRUE;

   SCIPdebugMsg(scip, "detecting nlhdlrs for %s expression %p (%s); requiring%s%s%s\n",
      cons != NULL ? "root" : "non-root", (void*)expr, SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)),
      (enforcemethods & SCIP_NLHDLR_METHOD_SEPABELOW) != 0 ? "" : " sepabelow",
      (enforcemethods & SCIP_NLHDLR_METHOD_SEPAABOVE) != 0 ? "" : " sepaabove",
      (enforcemethods & SCIP_NLHDLR_METHOD_ACTIVITY) != 0 ? "" : " activity");

   for( h = 0; h < conshdlrdata->nnlhdlrs; ++h )
   {
      SCIP_NLHDLR* nlhdlr;

      nlhdlr = conshdlrdata->nlhdlrs[h];
      assert(nlhdlr != NULL);

      /* skip disabled nlhdlrs */
      if( !SCIPnlhdlrIsEnabled(nlhdlr) )
         continue;

      /* call detect routine of nlhdlr */
      nlhdlrexprdata = NULL;
      enforcemethodsnew = enforcemethods;
      nlhdlrparticipating = SCIP_NLHDLR_METHOD_NONE;
      conshdlrdata->registerusesactivitysepabelow = FALSE;  /* SCIPregisterExprUsageNonlinear() as called by detect may set this to TRUE */
      conshdlrdata->registerusesactivitysepaabove = FALSE;  /* SCIPregisterExprUsageNonlinear() as called by detect may set this to TRUE */
      SCIP_CALL( SCIPnlhdlrDetect(scip, ownerdata->conshdlr, nlhdlr, expr, cons, &enforcemethodsnew, &nlhdlrparticipating, &nlhdlrexprdata) );

      /* nlhdlr might have claimed more than needed: clean up sepa flags */
      nlhdlrparticipating &= enforcemethodsallowed;

      /* detection is only allowed to augment to nlhdlrenforcemethods, so previous enforcemethods must still be set */
      assert((enforcemethodsnew & enforcemethods) == enforcemethods);

      /* Because of the previous assert, nlhdlrenforcenew ^ enforcemethods are the methods enforced by this nlhdlr.
       * They are also cleaned up here to ensure that only the needed methods are claimed.
       */
      nlhdlrenforcemethods = (enforcemethodsnew ^ enforcemethods) & enforcemethodsallowed;

      /* nlhdlr needs to participate for the methods it is enforcing */
      assert((nlhdlrparticipating & nlhdlrenforcemethods) == nlhdlrenforcemethods);

      if( nlhdlrparticipating == SCIP_NLHDLR_METHOD_NONE )
      {
         /* nlhdlr might not have detected anything, or all set flags might have been removed by
          * clean up; in the latter case, we may need to free nlhdlrexprdata */

         /* free nlhdlr exprdata, if there is any and there is a method to free this data */
         if( nlhdlrexprdata != NULL )
         {
            SCIP_CALL( SCIPnlhdlrFreeexprdata(scip, nlhdlr, expr, &nlhdlrexprdata) );
         }
         /* nlhdlr cannot have added an enforcement method if it doesn't participate (actually redundant due to previous asserts) */
         assert(nlhdlrenforcemethods == SCIP_NLHDLR_METHOD_NONE);

         SCIPdebugMsg(scip, "nlhdlr <%s> detect unsuccessful\n", SCIPnlhdlrGetName(nlhdlr));

         continue;
      }

      SCIPdebugMsg(scip, "nlhdlr <%s> detect successful; sepabelow: %s, sepaabove: %s, activity: %s\n",
         SCIPnlhdlrGetName(nlhdlr),
         ((nlhdlrenforcemethods & SCIP_NLHDLR_METHOD_SEPABELOW) != 0) ? "enforcing" : ((nlhdlrparticipating & SCIP_NLHDLR_METHOD_SEPABELOW) != 0) ? "participating" : "no",
         ((nlhdlrenforcemethods & SCIP_NLHDLR_METHOD_SEPAABOVE) != 0) ? "enforcing" : ((nlhdlrparticipating & SCIP_NLHDLR_METHOD_SEPAABOVE) != 0) ? "participating" : "no",
         ((nlhdlrenforcemethods & SCIP_NLHDLR_METHOD_ACTIVITY) != 0) ? "enforcing" : ((nlhdlrparticipating & SCIP_NLHDLR_METHOD_ACTIVITY) != 0) ? "participating" : "no");

      /* store nlhdlr and its data */
      SCIP_CALL( SCIPensureBlockMemoryArray(scip, &ownerdata->enfos, &enfossize, ownerdata->nenfos+1) );
      SCIP_CALL( SCIPallocBlockMemory(scip, &ownerdata->enfos[ownerdata->nenfos]) );
      ownerdata->enfos[ownerdata->nenfos]->nlhdlr = nlhdlr;
      ownerdata->enfos[ownerdata->nenfos]->nlhdlrexprdata = nlhdlrexprdata;
      ownerdata->enfos[ownerdata->nenfos]->nlhdlrparticipation = nlhdlrparticipating;
      ownerdata->enfos[ownerdata->nenfos]->issepainit = FALSE;
      ownerdata->enfos[ownerdata->nenfos]->sepabelowusesactivity = conshdlrdata->registerusesactivitysepabelow;
      ownerdata->enfos[ownerdata->nenfos]->sepaaboveusesactivity = conshdlrdata->registerusesactivitysepaabove;
      ownerdata->nenfos++;

      /* update enforcement flags */
      enforcemethods = enforcemethodsnew;
   }

   conshdlrdata->indetect = FALSE;

   /* stop if an enforcement method is missing but we are already in solving stage
    * (as long as the expression provides its callbacks, the default nlhdlr should have provided all enforcement methods)
    */
   if( enforcemethods != SCIP_NLHDLR_METHOD_ALL && SCIPgetStage(scip) == SCIP_STAGE_SOLVING )
   {
      SCIPerrorMessage("no nonlinear handler provided some of the required enforcement methods\n");
      return SCIP_ERROR;
   }

   assert(ownerdata->nenfos > 0);

   /* sort nonlinear handlers by enforcement priority, in decreasing order */
   if( ownerdata->nenfos > 1 )
      SCIPsortDownPtr((void**)ownerdata->enfos, enfodataCmp, ownerdata->nenfos);

   /* resize enfos array to be nenfos long */
   SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &ownerdata->enfos, enfossize, ownerdata->nenfos) );

   return SCIP_OKAY;
}

/** detect nlhdlrs that can handle the expressions */
static
SCIP_RETCODE detectNlhdlrs(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS**           conss,              /**< constraints for which to run nlhdlr detect */
   int                   nconss              /**< total number of constraints */
   )
{
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_CONSDATA* consdata;
   SCIP_EXPR* expr;
   SCIP_EXPR_OWNERDATA* ownerdata;
   SCIP_EXPRITER* it;
   int i;

   assert(conss != NULL || nconss == 0);
   assert(nconss >= 0);
   assert(SCIPgetStage(scip) == SCIP_STAGE_PRESOLVING || SCIPgetStage(scip) == SCIP_STAGE_INITSOLVE || SCIPgetStage(scip) == SCIP_STAGE_SOLVING);  /* should only be called in presolve or initsolve or consactive */

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);

   SCIP_CALL( SCIPcreateExpriter(scip, &it) );
   SCIP_CALL( SCIPexpriterInit(it, NULL, SCIP_EXPRITER_DFS, TRUE) );

   if( SCIPgetStage(scip) == SCIP_STAGE_SOLVING && SCIPgetDepth(scip) != 0 )
   {
      /* ensure that activities are recomputed w.r.t. the global variable bounds if CONSACTIVE is called in a local node;
       * for example, this happens if globally valid nonlinear constraints are added during the tree search
       */
      SCIPincrementCurBoundsTagNonlinear(conshdlr, TRUE);
      conshdlrdata->globalbounds = TRUE;
      conshdlrdata->lastvaractivitymethodchange = conshdlrdata->curboundstag;
   }

   for( i = 0; i < nconss; ++i )
   {
      assert(conss != NULL && conss[i] != NULL);

      consdata = SCIPconsGetData(conss[i]);
      assert(consdata != NULL);
      assert(consdata->expr != NULL);

      /* if a constraint is separated, we currently need it to be initial, too
       * this is because INITLP will create the auxiliary variables that are used for any separation
       * TODO we may relax this with a little more programming effort when required, see also TODO in INITLP
       */
      assert((!SCIPconsIsSeparated(conss[i]) && !SCIPconsIsEnforced(conss[i])) || SCIPconsIsInitial(conss[i]));

      ownerdata = SCIPexprGetOwnerData(consdata->expr);
      assert(ownerdata != NULL);

      /* because of common sub-expressions it might happen that we already detected a nonlinear handler and added it to the expr
       * then we would normally skip to run DETECT again
       * HOWEVER: most likely we have been running DETECT with cons == NULL, which may interest less nlhdlrs
       * thus, if expr is the root expression, we rerun DETECT
       */
      if( ownerdata->nenfos > 0 )
      {
         SCIP_CALL( freeEnfoData(scip, consdata->expr, FALSE) );
         assert(ownerdata->nenfos < 0);
      }

      /* if constraint will be enforced, and we are in solve, then ensure auxiliary variable for root expression
       *   this way we can treat the root expression like any other expression when enforcing via separation
       * if constraint will be propagated, then register activity usage of root expression
       * this can trigger a call to forwardPropExpr, for which we better have the indetect flag set
       */
      conshdlrdata->indetect = TRUE;
      SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, consdata->expr,
         SCIPgetStage(scip) >= SCIP_STAGE_INITSOLVE && (SCIPconsIsSeparated(conss[i]) || SCIPconsIsEnforced(conss[i])),
         SCIPconsIsPropagated(conss[i]),
         FALSE, FALSE) );
      conshdlrdata->indetect = FALSE;

      /* compute integrality information for all subexpressions */
      SCIP_CALL( SCIPcomputeExprIntegrality(scip, consdata->expr) );

      /* run detectNlhdlr on all expr where required */
      for( expr = SCIPexpriterRestartDFS(it, consdata->expr); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
      {
         ownerdata = SCIPexprGetOwnerData(expr);
         assert(ownerdata != NULL);

         /* skip exprs that we already looked at */
         if( ownerdata->nenfos >= 0 )
            continue;

         /* if there is use of the auxvar, then someone requires that
          *   auxvar == expr (or auxvar >= expr or auxvar <= expr) or we are at the root expression (expr==consdata->expr)
          *   thus, we need to find nlhdlrs that separate or estimate
          * if there is use of the activity, then there is someone requiring that
          *   activity of this expression is updated; this someone would also benefit from better bounds on the activity of this expression
          *   thus, we need to find nlhdlrs that do interval-evaluation
          */
         if( ownerdata->nauxvaruses > 0 || ownerdata->nactivityusesprop > 0 || ownerdata->nactivityusessepa > 0 )
         {
            SCIP_CALL( detectNlhdlr(scip, expr, expr == consdata->expr ? conss[i] : NULL) );

            assert(ownerdata->nenfos >= 0);
         }
         else
         {
            /* remember that we looked at this expression during detectNlhdlrs
             * even though we have not actually run detectNlhdlr, because no nlhdlr showed interest in this expr,
             * in some situations (forwardPropExpr, to be specific) we will have to distinguish between exprs for which
             * we have not initialized enforcement yet (nenfos < 0) and expressions which are just not used in enforcement (nenfos == 0)
             */
            ownerdata->nenfos = 0;
         }
      }

      /* include this constraint into the next propagation round because the added nlhdlr may do find tighter bounds now */
      if( SCIPconsIsPropagated(conss[i]) )
         consdata->ispropagated = FALSE;
   }

   if( SCIPgetStage(scip) == SCIP_STAGE_SOLVING && SCIPgetDepth(scip) != 0 )
   {
      /* ensure that the local bounds are used again when reevaluating the expressions later;
       * this is only needed if CONSACTIVE is called in a local node (see begin of this function)
       */
      SCIPincrementCurBoundsTagNonlinear(conshdlr, FALSE);
      conshdlrdata->globalbounds = FALSE;
      conshdlrdata->lastvaractivitymethodchange = conshdlrdata->curboundstag;
   }
   else
   {
      /* ensure that all activities (except for var-exprs) are reevaluated since better methods may be available now */
      SCIPincrementCurBoundsTagNonlinear(conshdlr, FALSE);
   }

   SCIPfreeExpriter(&it);

   return SCIP_OKAY;
}

/** initializes (pre)solving data of constraints
 *
 * This initializes data in a constraint that is used for separation, propagation, etc, and assumes that expressions will
 * not be modified.
 * In particular, this function
 * - runs the detection method of nlhldrs
 * - looks for unlocked linear variables
 * - checks curvature (if not in presolve)
 * - creates and add row to NLP (if not in presolve)
 *
 * This function can be called in presolve and solve and can be called several times with different sets of constraints,
 * e.g., it should be called in INITSOL and for constraints that are added during solve.
 */
static
SCIP_RETCODE initSolve(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS**           conss,              /**< constraints */
   int                   nconss              /**< number of constraints */
   )
{
   int c;

   for( c = 0; c < nconss; ++c )
   {
      /* check for a linear variable that can be increase or decreased without harming feasibility */
      findUnlockedLinearVar(scip, conss[c]);

      if( SCIPgetStage(scip) == SCIP_STAGE_INITSOLVE || SCIPgetStage(scip) == SCIP_STAGE_SOLVING )
      {
         SCIP_CONSDATA* consdata;
         SCIP_Bool success = FALSE;

         consdata = SCIPconsGetData(conss[c]);  /*lint !e613*/
         assert(consdata != NULL);
         assert(consdata->expr != NULL);

         if( !SCIPconshdlrGetData(conshdlr)->assumeconvex )
         {
            /* call the curvature detection algorithm of the convex nonlinear handler
             * Check only for those curvature that may result in a convex inequality, i.e.,
             * whether f(x) is concave when f(x) >= lhs and/or f(x) is convex when f(x) <= rhs.
             * Also we can assume that we are nonlinear, so do not check for convex if already concave.
             */
            if( !SCIPisInfinity(scip, -consdata->lhs) )
            {
               SCIP_CALL( SCIPhasExprCurvature(scip, consdata->expr, SCIP_EXPRCURV_CONCAVE, &success, NULL) );
               if( success )
                  consdata->curv = SCIP_EXPRCURV_CONCAVE;
            }
            if( !success && !SCIPisInfinity(scip, consdata->rhs) )
            {
               SCIP_CALL( SCIPhasExprCurvature(scip, consdata->expr, SCIP_EXPRCURV_CONVEX, &success, NULL) );
               if( success )
                  consdata->curv = SCIP_EXPRCURV_CONVEX;
            }
         }
         else
         {
            if( !SCIPisInfinity(scip, -consdata->lhs) && !SCIPisInfinity(scip, consdata->rhs) )
            {
               SCIPwarningMessage(scip, "Nonlinear constraint <%s> has finite left- and right-hand side, but constraints/nonlinear/assumeconvex is enabled.\n", SCIPconsGetName(conss[c]));
               consdata->curv = SCIP_EXPRCURV_LINEAR;
            }
            else
            {
               consdata->curv = !SCIPisInfinity(scip, consdata->rhs) ? SCIP_EXPRCURV_CONVEX : SCIP_EXPRCURV_CONCAVE;
            }
         }
         SCIPdebugMsg(scip, "root curvature of constraint %s = %d\n", SCIPconsGetName(conss[c]), consdata->curv);

         /* add nlrow representation to NLP, if NLP had been constructed */
         if( SCIPisNLPConstructed(scip) && SCIPconsIsActive(conss[c]) )
         {
            if( consdata->nlrow == NULL )
            {
               SCIP_CALL( createNlRow(scip, conss[c]) );
               assert(consdata->nlrow != NULL);
            }
            SCIPnlrowSetCurvature(consdata->nlrow, consdata->curv);
            SCIP_CALL( SCIPaddNlRow(scip, consdata->nlrow) );
         }
      }
   }

   /* register non linear handlers */
   SCIP_CALL( detectNlhdlrs(scip, conshdlr, conss, nconss) );

   return SCIP_OKAY;
}

/** deinitializes (pre)solving data of constraints
 *
 * This removes the initialization data created in initSolve().
 *
 * This function can be called in presolve and solve.
 *
 * TODO At the moment, it should not be called for a constraint if there are other constraints
 * that use the same expressions but still require their nlhdlr.
 * We should probably only decrement the auxvar and activity usage for the root expr and then
 * proceed as in detectNlhdlrs(), i.e., free enfo data only where none is used.
 */
static
SCIP_RETCODE deinitSolve(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS**           conss,              /**< constraints */
   int                   nconss              /**< number of constraints */
   )
{
   SCIP_EXPRITER* it;
   SCIP_EXPR* expr;
   SCIP_CONSDATA* consdata;
   SCIP_Bool rootactivityvalid;
   int c;

   SCIP_CALL( SCIPcreateExpriter(scip, &it) );
   SCIP_CALL( SCIPexpriterInit(it, NULL, SCIP_EXPRITER_DFS, FALSE) );
   SCIPexpriterSetStagesDFS(it, SCIP_EXPRITER_LEAVEEXPR);

   /* call deinitialization callbacks of expression and nonlinear handlers
    * free nonlinear handlers information from expressions
    * remove auxiliary variables and nactivityuses counts from expressions
    */
   for( c = 0; c < nconss; ++c )
   {
      assert(conss != NULL);
      assert(conss[c] != NULL);

      consdata = SCIPconsGetData(conss[c]);
      assert(consdata != NULL);
      assert(consdata->expr != NULL);

      /* check and remember whether activity in root is valid */
      rootactivityvalid = SCIPexprGetActivityTag(consdata->expr) >= SCIPconshdlrGetData(conshdlr)->lastboundrelax;

      for( expr = SCIPexpriterRestartDFS(it, consdata->expr); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
      {
         SCIPdebugMsg(scip, "exitsepa and free nonlinear handler data for expression %p\n", (void*)expr);

         /* remove nonlinear handlers in expression and their data and auxiliary variables; reset activityusage count */
         SCIP_CALL( freeEnfoData(scip, expr, TRUE) );

         /* remove quadratic info */
         SCIPfreeExprQuadratic(scip, expr);

         if( rootactivityvalid )
         {
            /* ensure activity is valid if consdata->expr activity is valid
             * this is mainly to ensure that we do not leave invalid activities in parts of the expression tree where activity was not used,
             * e.g., an expr's activity was kept up to date by a nlhdlr, but without using some childs activity
             * so this childs activity would be invalid, which can generate confusion
             */
            SCIP_CALL( SCIPevalExprActivity(scip, expr) );
         }
      }

      if( consdata->nlrow != NULL )
      {
         /* remove row from NLP, if still in solving
          * if we are in exitsolve, the whole NLP will be freed anyway
          */
         if( SCIPgetStage(scip) == SCIP_STAGE_SOLVING )
         {
            SCIP_CALL( SCIPdelNlRow(scip, consdata->nlrow) );
         }

         SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
      }

      /* forget about linear variables that can be increased or decreased without harming feasibility */
      consdata->linvardecr = NULL;
      consdata->linvarincr = NULL;

      /* forget about curvature */
      consdata->curv = SCIP_EXPRCURV_UNKNOWN;
   }

   SCIPfreeExpriter(&it);

   return SCIP_OKAY;
}

/** helper method to decide whether a given expression is product of at least two binary variables */
static
SCIP_Bool isBinaryProduct(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EXPR*            expr                /**< expression */
   )
{
   int i;

   assert(expr != NULL);

   /* check whether the expression is a product */
   if( !SCIPisExprProduct(scip, expr) )
      return FALSE;

   /* don't consider products with a coefficient != 1 and products with a single child
    * simplification will take care of this expression later
    */
   if( SCIPexprGetNChildren(expr) <= 1 || SCIPgetCoefExprProduct(expr) != 1.0 )
      return FALSE;

   for( i = 0; i < SCIPexprGetNChildren(expr); ++i )
   {
      SCIP_EXPR* child;
      SCIP_VAR* var;
      SCIP_Real ub;
      SCIP_Real lb;

      child = SCIPexprGetChildren(expr)[i];
      assert(child != NULL);

      if( !SCIPisExprVar(scip, child) )
         return FALSE;

      var = SCIPgetVarExprVar(child);
      lb = SCIPvarGetLbLocal(var);
      ub = SCIPvarGetUbLocal(var);

      /* check whether variable is integer and has [0,1] as variable bounds */
      if( !SCIPvarIsIntegral(var) || !SCIPisEQ(scip, lb, 0.0) || !SCIPisEQ(scip, ub, 1.0) )
         return FALSE;
   }

   return TRUE;
}

/** helper method to collect all bilinear binary product terms */
static
SCIP_RETCODE getBilinearBinaryTerms(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EXPR*            sumexpr,            /**< sum expression */
   SCIP_VAR**            xs,                 /**< array to collect first variable of each bilinear binary product */
   SCIP_VAR**            ys,                 /**< array to collect second variable of each bilinear binary product */
   int*                  childidxs,          /**< array to store the index of the child of each stored bilinear binary product */
   int*                  nterms              /**< pointer to store the total number of bilinear binary terms */
   )
{
   int i;

   assert(sumexpr != NULL);
   assert(SCIPisExprSum(scip, sumexpr));
   assert(xs != NULL);
   assert(ys != NULL);
   assert(childidxs != NULL);
   assert(nterms != NULL);

   *nterms = 0;

   for( i = 0; i < SCIPexprGetNChildren(sumexpr); ++i )
   {
      SCIP_EXPR* child;

      child = SCIPexprGetChildren(sumexpr)[i];
      assert(child != NULL);

      if( SCIPexprGetNChildren(child) == 2 && isBinaryProduct(scip, child) )
      {
         SCIP_VAR* x = SCIPgetVarExprVar(SCIPexprGetChildren(child)[0]);
         SCIP_VAR* y = SCIPgetVarExprVar(SCIPexprGetChildren(child)[1]);

         assert(x != NULL);
         assert(y != NULL);

         if( x != y )
         {
            xs[*nterms] = x;
            ys[*nterms] = y;
            childidxs[*nterms] = i;
            ++(*nterms);
         }
      }
   }

   return SCIP_OKAY;
}

/** helper method to reformulate \f$x_i \sum_j c_{ij} x_j\f$ */
static
SCIP_RETCODE reformulateFactorizedBinaryQuadratic(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS*            cons,               /**< constraint */
   SCIP_VAR*             facvar,             /**< variable that has been factorized */
   SCIP_VAR**            vars,               /**< variables of sum_j c_ij x_j */
   SCIP_Real*            coefs,              /**< coefficients of sum_j c_ij x_j */
   int                   nvars,              /**< total number of variables in sum_j c_ij x_j */
   SCIP_EXPR**           newexpr,            /**< pointer to store the new expression */
   int*                  naddconss           /**< pointer to update the total number of added constraints (might be NULL) */
   )
{
   SCIP_VAR* auxvar;
   SCIP_CONS* newcons;
   SCIP_Real minact = 0.0;
   SCIP_Real maxact = 0.0;
   SCIP_Bool integral = TRUE;
   char name [SCIP_MAXSTRLEN];
   int i;

   assert(facvar != NULL);
   assert(vars != NULL);
   assert(nvars > 1);
   assert(newexpr != NULL);

   /* compute minimum and maximum activity of sum_j c_ij x_j */
   /* TODO could compute minact and maxact for facvar=0 and facvar=1 separately, taking implied bounds into account, allowing for possibly tighter big-M's below */
   for( i = 0; i < nvars; ++i )
   {
      minact += MIN(coefs[i], 0.0);
      maxact += MAX(coefs[i], 0.0);
      integral = integral && SCIPisIntegral(scip, coefs[i]);
   }
   assert(minact <= maxact);

   /* create and add auxiliary variable */
   (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s", SCIPconsGetName(cons), SCIPvarGetName(facvar));
   SCIP_CALL( SCIPcreateVarBasic(scip, &auxvar, name, minact, maxact, 0.0, integral ? SCIP_VARTYPE_IMPLINT : SCIP_VARTYPE_CONTINUOUS) );
   SCIP_CALL( SCIPaddVar(scip, auxvar) );

   /* create and add z - maxact x <= 0 */
   if( !SCIPisZero(scip, maxact) )
   {
      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s_1", SCIPconsGetName(cons), SCIPvarGetName(facvar));
      SCIP_CALL( SCIPcreateConsBasicVarbound(scip, &newcons, name, auxvar, facvar, -maxact, -SCIPinfinity(scip), 0.0) );
      SCIP_CALL( SCIPaddCons(scip, newcons) );
      SCIP_CALL( SCIPreleaseCons(scip, &newcons) );
      if( naddconss != NULL )
         ++(*naddconss);
   }

   /* create and add  0 <= z - minact x */
   if( !SCIPisZero(scip, minact) )
   {
      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s_2", SCIPconsGetName(cons), SCIPvarGetName(facvar));
      SCIP_CALL( SCIPcreateConsBasicVarbound(scip, &newcons, name, auxvar, facvar, -minact, 0.0, SCIPinfinity(scip)) );
      SCIP_CALL( SCIPaddCons(scip, newcons) );
      SCIP_CALL( SCIPreleaseCons(scip, &newcons) );
      if( naddconss != NULL )
         ++(*naddconss);
   }

   /* create and add minact <= sum_j c_j x_j - z + minact x_i */
   (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s_3", SCIPconsGetName(cons), SCIPvarGetName(facvar));
   SCIP_CALL( SCIPcreateConsBasicLinear(scip, &newcons, name, nvars, vars, coefs, minact, SCIPinfinity(scip)) );
   SCIP_CALL( SCIPaddCoefLinear(scip, newcons, auxvar, -1.0) );
   if( !SCIPisZero(scip, minact) )
   {
      SCIP_CALL( SCIPaddCoefLinear(scip, newcons, facvar, minact) );
   }
   SCIP_CALL( SCIPaddCons(scip, newcons) );
   SCIP_CALL( SCIPreleaseCons(scip, &newcons) );
   if( naddconss != NULL )
      ++(*naddconss);

   /* create and add sum_j c_j x_j - z + maxact x_i <= maxact */
   (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s_4", SCIPconsGetName(cons), SCIPvarGetName(facvar));
   SCIP_CALL( SCIPcreateConsBasicLinear(scip, &newcons, name, nvars, vars, coefs, -SCIPinfinity(scip), maxact) );
   SCIP_CALL( SCIPaddCoefLinear(scip, newcons, auxvar, -1.0) );
   if( !SCIPisZero(scip, maxact) )
   {
      SCIP_CALL( SCIPaddCoefLinear(scip, newcons, facvar, maxact) );
   }
   SCIP_CALL( SCIPaddCons(scip, newcons) );
   SCIP_CALL( SCIPreleaseCons(scip, &newcons) );
   if( naddconss != NULL )
      ++(*naddconss);

   /* create variable expression */
   SCIP_CALL( createExprVar(scip, conshdlr, newexpr, auxvar) );

   /* release auxvar */
   SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );

   return SCIP_OKAY;
}

/** helper method to generate an expression for a sum of products of binary variables; note that the method captures the generated expression */
static
SCIP_RETCODE getFactorizedBinaryQuadraticExpr(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS*            cons,               /**< constraint */
   SCIP_EXPR*            sumexpr,            /**< expression */
   int                   minterms,           /**< minimum number of terms in a the sum of x_i sum_j c_j x_j */
   SCIP_EXPR**           newexpr,            /**< pointer to store the expression that represents the binary quadratic */
   int*                  naddconss           /**< pointer to update the total number of added constraints (might be NULL) */
   )
{
   SCIP_EXPR** exprs = NULL;
   SCIP_VAR** tmpvars = NULL;
   SCIP_VAR** vars = NULL;
   SCIP_VAR** xs = NULL;
   SCIP_VAR** ys = NULL;
   SCIP_Real* exprcoefs = NULL;
   SCIP_Real* tmpcoefs = NULL;
   SCIP_Real* sumcoefs;
   SCIP_Bool* isused  = NULL;
   int* childidxs = NULL;
   int* count = NULL;
   int nchildren;
   int nexprs = 0;
   int nterms;
   int nvars;
   int ntotalvars;
   int i;

   assert(sumexpr != NULL);
   assert(minterms > 1);
   assert(newexpr != NULL);

   *newexpr = NULL;

   /* check whether sumexpr is indeed a sum */
   if( !SCIPisExprSum(scip, sumexpr) )
      return SCIP_OKAY;

   nchildren = SCIPexprGetNChildren(sumexpr);
   sumcoefs = SCIPgetCoefsExprSum(sumexpr);
   nvars = SCIPgetNVars(scip);
   ntotalvars = SCIPgetNTotalVars(scip);

   /* check whether there are enough terms available */
   if( nchildren < minterms )
      return SCIP_OKAY;

   /* allocate memory */
   SCIP_CALL( SCIPallocBufferArray(scip, &xs, nchildren) );
   SCIP_CALL( SCIPallocBufferArray(scip, &ys, nchildren) );
   SCIP_CALL( SCIPallocBufferArray(scip, &childidxs, nchildren) );

   /* collect all bilinear binary product terms */
   SCIP_CALL( getBilinearBinaryTerms(scip, sumexpr, xs, ys, childidxs, &nterms) );

   /* check whether there are enough terms available */
   if( nterms < minterms )
      goto TERMINATE;

   /* store how often each variable appears in a bilinear binary product */
   SCIP_CALL( SCIPduplicateBufferArray(scip, &vars, SCIPgetVars(scip), nvars) );
   SCIP_CALL( SCIPallocClearBufferArray(scip, &count, ntotalvars) );
   SCIP_CALL( SCIPallocClearBufferArray(scip, &isused, nchildren) );

   SCIP_CALL( SCIPallocBufferArray(scip, &exprs, nchildren) );
   SCIP_CALL( SCIPallocBufferArray(scip, &exprcoefs, nchildren) );
   SCIP_CALL( SCIPallocBufferArray(scip, &tmpvars, MIN(nterms, nvars)) );
   SCIP_CALL( SCIPallocBufferArray(scip, &tmpcoefs, MIN(nterms, nvars)) );

   for( i = 0; i < nterms; ++i )
   {
      int xidx;
      int yidx;

      assert(xs[i] != NULL);
      assert(ys[i] != NULL);

      xidx = SCIPvarGetIndex(xs[i]);
      assert(xidx < ntotalvars);
      yidx = SCIPvarGetIndex(ys[i]);
      assert(yidx < ntotalvars);

      ++count[xidx];
      ++count[yidx];

      SCIPdebugMsg(scip, "increase counter for %s to %d\n", SCIPvarGetName(xs[i]), count[xidx]);
      SCIPdebugMsg(scip, "increase counter for %s to %d\n", SCIPvarGetName(ys[i]), count[yidx]);
   }

   /* sort variables; don't change order of count array because it depends on problem indices */
   {
      int* tmpcount;

      SCIP_CALL( SCIPduplicateBufferArray(scip, &tmpcount, count, nvars) );
      SCIPsortDownIntPtr(tmpcount, (void**)vars, nvars);
      SCIPfreeBufferArray(scip, &tmpcount);
   }

   for( i = 0; i < nvars; ++i )
   {
      SCIP_VAR* facvar = vars[i];
      int ntmpvars = 0;
      int j;

      /* skip candidate if there are not enough terms left */
      if( count[SCIPvarGetIndex(vars[i])] < minterms )
         continue;

      SCIPdebugMsg(scip, "consider facvar = %s with count = %d\n", SCIPvarGetName(facvar), count[SCIPvarGetIndex(vars[i])]);

      /* collect variables for x_i * sum_j c_ij x_j */
      for( j = 0; j < nterms; ++j )
      {
         int childidx = childidxs[j];
         assert(childidx >= 0 && childidx < nchildren);

         if( !isused[childidx] && (xs[j] == facvar || ys[j] == facvar) )
         {
            SCIP_Real coef;
            int xidx;
            int yidx;

            coef = sumcoefs[childidx];
            assert(coef != 0.0);

            /* collect corresponding variable */
            tmpvars[ntmpvars] = (xs[j] == facvar) ? ys[j] : xs[j];
            tmpcoefs[ntmpvars] = coef;
            ++ntmpvars;

            /* update counters */
            xidx = SCIPvarGetIndex(xs[j]);
            assert(xidx < ntotalvars);
            yidx = SCIPvarGetIndex(ys[j]);
            assert(yidx < ntotalvars);
            --count[xidx];
            --count[yidx];
            assert(count[xidx] >= 0);
            assert(count[yidx] >= 0);

            /* mark term to be used */
            isused[childidx] = TRUE;
         }
      }
      assert(ntmpvars >= minterms);
      assert(SCIPvarGetIndex(facvar) < ntotalvars);
      assert(count[SCIPvarGetIndex(facvar)] == 0); /* facvar should not appear in any other bilinear term */

      /* create required constraints and store the generated expression */
      SCIP_CALL( reformulateFactorizedBinaryQuadratic(scip, conshdlr, cons, facvar, tmpvars, tmpcoefs, ntmpvars, &exprs[nexprs], naddconss) );
      exprcoefs[nexprs] = 1.0;
      ++nexprs;
   }

   /* factorization was only successful if at least one expression has been generated */
   if( nexprs > 0 )
   {
      int nexprsold = nexprs;

      /* add all children of the sum that have not been used */
      for( i = 0; i < nchildren; ++i )
      {
         if( !isused[i] )
         {
            exprs[nexprs] = SCIPexprGetChildren(sumexpr)[i];
            exprcoefs[nexprs] = sumcoefs[i];
            ++nexprs;
         }
      }

      /* create a new sum expression */
      SCIP_CALL( SCIPcreateExprSum(scip, newexpr, nexprs, exprs, exprcoefs, SCIPgetConstantExprSum(sumexpr), exprownerCreate, (void*)conshdlr) );

      /* release all expressions that have been generated by reformulateFactorizedBinaryQuadratic() */
      for( i = 0; i < nexprsold; ++i )
      {
         SCIP_CALL( SCIPreleaseExpr(scip, &exprs[i]) );
      }
   }

TERMINATE:
   /* free memory */
   SCIPfreeBufferArrayNull(scip, &tmpcoefs);
   SCIPfreeBufferArrayNull(scip, &tmpvars);
   SCIPfreeBufferArrayNull(scip, &exprcoefs);
   SCIPfreeBufferArrayNull(scip, &exprs);
   SCIPfreeBufferArrayNull(scip, &vars);
   SCIPfreeBufferArrayNull(scip, &isused);
   SCIPfreeBufferArrayNull(scip, &count);
   SCIPfreeBufferArray(scip, &childidxs);
   SCIPfreeBufferArray(scip, &ys);
   SCIPfreeBufferArray(scip, &xs);

   return SCIP_OKAY;
}

/** helper method to create an AND constraint or varbound constraints for a given binary product expression */
static
SCIP_RETCODE getBinaryProductExprDo(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_EXPR*            prodexpr,           /**< product expression */
   SCIP_EXPR**           newexpr,            /**< pointer to store the expression that represents the product */
   int*                  naddconss,          /**< pointer to update the total number of added constraints (might be NULL) */
   SCIP_Bool             empathy4and         /**< whether to use an AND constraint, if possible */
   )
{
   SCIP_VAR** vars;
   SCIP_CONS* cons;
   SCIP_Real* coefs;
   SCIP_VAR* w;
   char name[SCIP_MAXSTRLEN];
   int nchildren;
   int i;

   assert(conshdlr != NULL);
   assert(prodexpr != NULL);
   assert(SCIPisExprProduct(scip, prodexpr));
   assert(newexpr != NULL);

   nchildren = SCIPexprGetNChildren(prodexpr);
   assert(nchildren >= 2);

   /* memory to store the variables of the variable expressions (+1 for w) */
   SCIP_CALL( SCIPallocBufferArray(scip, &vars, nchildren + 1) );
   SCIP_CALL( SCIPallocBufferArray(scip, &coefs, nchildren + 1) );

   /* prepare the names of the variable and the constraints */
   (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform");
   for( i = 0; i < nchildren; ++i )
   {
      vars[i] = SCIPgetVarExprVar(SCIPexprGetChildren(prodexpr)[i]);
      coefs[i] = 1.0;
      assert(vars[i] != NULL);
      (void) strcat(name, "_");
      (void) strcat(name, SCIPvarGetName(vars[i]));
   }

   /* create and add variable */
   SCIP_CALL( SCIPcreateVarBasic(scip, &w, name, 0.0, 1.0, 0.0, SCIP_VARTYPE_IMPLINT) );
   SCIP_CALL( SCIPaddVar(scip, w) );
   SCIPdebugMsg(scip, "  created auxiliary variable %s\n", name);

   /* use variable bound constraints if it is a bilinear product and there is no empathy for an AND constraint */
   if( nchildren == 2 && !empathy4and )
   {
      SCIP_VAR* x = vars[0];
      SCIP_VAR* y = vars[1];

      assert(x != NULL);
      assert(y != NULL);
      assert(x != y);

      /* create and add x - w >= 0 */
      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s_1", SCIPvarGetName(x), SCIPvarGetName(y));
      SCIP_CALL( SCIPcreateConsBasicVarbound(scip, &cons, name, x, w, -1.0, 0.0, SCIPinfinity(scip)) );
      SCIP_CALL( SCIPaddCons(scip, cons) );
      SCIP_CALL( SCIPreleaseCons(scip, &cons) );

      /* create and add y - w >= 0 */
      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s_2", SCIPvarGetName(x), SCIPvarGetName(y));
      SCIP_CALL( SCIPcreateConsBasicVarbound(scip, &cons, name, y, w, -1.0, 0.0, SCIPinfinity(scip)) );
      SCIP_CALL( SCIPaddCons(scip, cons) );
      SCIP_CALL( SCIPreleaseCons(scip, &cons) );

      /* create and add x + y - w <= 1 */
      vars[2] = w;
      coefs[2] = -1.0;
      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s_3", SCIPvarGetName(x), SCIPvarGetName(y));
      SCIP_CALL( SCIPcreateConsBasicLinear(scip, &cons, name, 3, vars, coefs, -SCIPinfinity(scip), 1.0) );
      SCIP_CALL( SCIPaddCons(scip, cons) );
      SCIP_CALL( SCIPreleaseCons(scip, &cons) );

      /* update number of added constraints */
      if( naddconss != NULL )
         *naddconss += 3;
   }
   else
   {
      /* create, add, and release AND constraint */
      SCIP_CALL( SCIPcreateConsBasicAnd(scip, &cons, name, w, nchildren, vars) );
      SCIP_CALL( SCIPaddCons(scip, cons) );
      SCIP_CALL( SCIPreleaseCons(scip, &cons) );
      SCIPdebugMsg(scip, "  create AND constraint\n");

      /* update number of added constraints */
      if( naddconss != NULL )
         *naddconss += 1;
   }

   /* create variable expression */
   SCIP_CALL( createExprVar(scip, conshdlr, newexpr, w) );

   /* release created variable */
   SCIP_CALL( SCIPreleaseVar(scip, &w) );

   /* free memory */
   SCIPfreeBufferArray(scip, &coefs);
   SCIPfreeBufferArray(scip, &vars);

   return SCIP_OKAY;
}

/** helper method to generate an expression for the product of binary variables; note that the method captures the generated expression */
static
SCIP_RETCODE getBinaryProductExpr(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_HASHMAP*         exprmap,            /**< map to remember generated variables for visited product expressions */
   SCIP_EXPR*            prodexpr,           /**< product expression */
   SCIP_EXPR**           newexpr,            /**< pointer to store the expression that represents the product */
   int*                  naddconss,          /**< pointer to update the total number of added constraints (might be NULL) */
   int*                  nchgcoefs           /**< pointer to update the total number of changed coefficients (might be NULL) */
   )
{
   SCIP_CONSHDLRDATA* conshdlrdata;
   int nchildren;

   assert(prodexpr != NULL);
   assert(newexpr != NULL);

   *newexpr = NULL;

   /* only consider products of binary variables */
   if( !isBinaryProduct(scip, prodexpr) )
      return SCIP_OKAY;

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);
   nchildren = SCIPexprGetNChildren(prodexpr);
   assert(nchildren >= 2);

   /* check whether there is already an expression that represents the product */
   if( SCIPhashmapExists(exprmap, (void*)prodexpr) )
   {
      *newexpr = (SCIP_EXPR*) SCIPhashmapGetImage(exprmap, (void*)prodexpr);
      assert(*newexpr != NULL);

      /* capture expression */
      SCIPcaptureExpr(*newexpr);
   }
   else
   {
      SCIPdebugMsg(scip, "  product expression %p has been considered for the first time\n", (void*)prodexpr);

      if( nchildren == 2 )
      {
         SCIP_CLIQUE** xcliques;
         SCIP_VAR* x;
         SCIP_VAR* y;
         SCIP_Bool found_clique = FALSE;
         int c;

         /* get variables from the product expression */
         x = SCIPgetVarExprVar(SCIPexprGetChildren(prodexpr)[0]);
         assert(x != NULL);
         y = SCIPgetVarExprVar(SCIPexprGetChildren(prodexpr)[1]);
         assert(y != NULL);
         assert(x != y);

         /* first try to find a clique containing both variables */
         xcliques = SCIPvarGetCliques(x, TRUE);

         /* look in cliques containing x */
         for( c = 0; c < SCIPvarGetNCliques(x, TRUE); ++c )
         {
            if( SCIPcliqueHasVar(xcliques[c], y, TRUE) ) /* x + y <= 1 => x*y = 0 */
            {
               /* create zero value expression */
               SCIP_CALL( SCIPcreateExprValue(scip, newexpr, 0.0, exprownerCreate, (void*)conshdlr) );

               if( nchgcoefs != NULL )
                  *nchgcoefs += 1;

               found_clique = TRUE;
               break;
            }

            if( SCIPcliqueHasVar(xcliques[c], y, FALSE) ) /* x + (1-y) <= 1 => x*y = x */
            {
               /* create variable expression for x */
               SCIP_CALL( createExprVar(scip, conshdlr, newexpr, x) );

               if( nchgcoefs != NULL )
                  *nchgcoefs += 2;

               found_clique = TRUE;
               break;
            }
         }

         if( !found_clique )
         {
            xcliques = SCIPvarGetCliques(x, FALSE);

            /* look in cliques containing complement of x */
            for( c = 0; c < SCIPvarGetNCliques(x, FALSE); ++c )
            {
               if( SCIPcliqueHasVar(xcliques[c], y, TRUE) ) /* (1-x) + y <= 1 => x*y = y */
               {
                  /* create variable expression for y */
                  SCIP_CALL( createExprVar(scip, conshdlr, newexpr, y) );

                  if( nchgcoefs != NULL )
                     *nchgcoefs += 1;

                  found_clique = TRUE;
                  break;
               }

               if( SCIPcliqueHasVar(xcliques[c], y, FALSE) ) /* (1-x) + (1-y) <= 1 => x*y = x + y - 1 */
               {
                  /* create sum expression */
                  SCIP_EXPR* sum_children[2];
                  SCIP_Real sum_coefs[2];
                  SCIP_CALL( createExprVar(scip, conshdlr, &sum_children[0], x) );
                  SCIP_CALL( createExprVar(scip, conshdlr, &sum_children[1], y) );
                  sum_coefs[0] = 1.0;
                  sum_coefs[1] = 1.0;
                  SCIP_CALL( SCIPcreateExprSum(scip, newexpr, 2, sum_children, sum_coefs, -1.0, exprownerCreate, (void*)conshdlr) );

                  SCIP_CALL( SCIPreleaseExpr(scip, &sum_children[0]) );
                  SCIP_CALL( SCIPreleaseExpr(scip, &sum_children[1]) );

                  if( nchgcoefs != NULL )
                     *nchgcoefs += 3;

                  found_clique = TRUE;
                  break;
               }
            }
         }

         /* if the variables are not in a clique, do standard linearization */
         if( !found_clique )
         {
            SCIP_CALL( getBinaryProductExprDo(scip, conshdlr, prodexpr, newexpr, naddconss, conshdlrdata->reformbinprodsand) );
         }
      }
      else
      {
         /* linearize binary product using an AND constraint because nchildren > 2 */
         SCIP_CALL( getBinaryProductExprDo(scip, conshdlr, prodexpr, newexpr, naddconss, conshdlrdata->reformbinprodsand) );
      }

      /* hash variable expression */
      SCIP_CALL( SCIPhashmapInsert(exprmap, (void*)prodexpr, *newexpr) );
   }

   return SCIP_OKAY;
}

/** helper function to replace binary products in a given constraint */
static
SCIP_RETCODE replaceBinaryProducts(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS*            cons,               /**< constraint */
   SCIP_HASHMAP*         exprmap,            /**< map to remember generated variables for visited product expressions */
   SCIP_EXPRITER*        it,                 /**< expression iterator */
   int*                  naddconss,          /**< pointer to update the total number of added constraints (might be NULL) */
   int*                  nchgcoefs           /**< pointer to update the total number of changed coefficients (might be NULL) */
   )
{
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_CONSDATA* consdata;
   SCIP_EXPR* expr;

   assert(conshdlr != NULL);
   assert(cons != NULL);
   assert(exprmap != NULL);
   assert(it != NULL);

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);
   assert(consdata->expr != NULL);

   SCIPdebugMsg(scip, "  check constraint %s\n", SCIPconsGetName(cons));

   for( expr = SCIPexpriterRestartDFS(it, consdata->expr); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
   {
      SCIP_EXPR* newexpr = NULL;
      SCIP_EXPR* childexpr;
      int childexpridx;

      childexpridx = SCIPexpriterGetChildIdxDFS(it);
      assert(childexpridx >= 0 && childexpridx < SCIPexprGetNChildren(expr));
      childexpr = SCIPexpriterGetChildExprDFS(it);
      assert(childexpr != NULL);

      /* try to factorize variables in a sum expression that contains several products of binary variables */
      if( conshdlrdata->reformbinprodsfac > 1 )
      {
         SCIP_CALL( getFactorizedBinaryQuadraticExpr(scip, conshdlr, cons, childexpr, conshdlrdata->reformbinprodsfac, &newexpr, naddconss) );
      }

      /* try to create an expression that represents a product of binary variables */
      if( newexpr == NULL )
      {
         SCIP_CALL( getBinaryProductExpr(scip, conshdlr, exprmap, childexpr, &newexpr, naddconss, nchgcoefs) );
      }

      if( newexpr != NULL )
      {
         assert(naddconss == NULL || *naddconss > 0 || nchgcoefs == NULL || *nchgcoefs > 0);

         /* replace product expression */
         SCIP_CALL( SCIPreplaceExprChild(scip, expr, childexpridx, newexpr) );

         /* note that the expression has been captured by getBinaryProductExpr and SCIPreplaceExprChild */
         SCIP_CALL( SCIPreleaseExpr(scip, &newexpr) );

         /* mark the constraint to not be simplified anymore */
         consdata->issimplified = FALSE;
      }
   }

   return SCIP_OKAY;
}

/** reformulates products of binary variables during presolving in the following way:
 *
 * Let \f$\sum_{i,j} Q_{ij} x_i x_j\f$ be a subexpression that only contains binary variables.
 * Each term \f$x_i x_j\f$ is reformulated with the help of an extra (implicit integer) variable \f$z_{ij}\f$ in {0,1}:
 * \f[
 *    z_{ij} \leq x_i, \qquad z_{ij} \leq x_j, \qquad x_i + x_j - z_{ij} \leq 1.
 * \f]
 *
 * Before reformulating \f$x_i x_j\f$ in this way, it is checked whether there is a clique that contains \f$x_i\f$ and \f$x_j\f$.
 * These cliques allow for a better reformulation. There are four cases:
 *
 *    1. \f$x_i + x_j \leq 1\f$ implies that \f$x_i x_j = 0\f$
 *    2. \f$x_i + (1 - x_j) \leq 1\f$ implies \f$x_i x_j = x_i\f$
 *    3. \f$(1 - x_i) + x_j \leq 1\f$ implies \f$x_i x_j = x_j\f$
 *    4. \f$(1 - x_i) + (1 - x_j) \leq 1\f$ implies \f$x_i x_j = x_i + x_j - 1\f$
 *
 * The reformulation using \f$z_{ij}\f$ or the cliques is implemented in getBinaryProductExpr().
 *
 * Introducing too many extra variables and constraints can have a negative impact on the performance (e.g., due to
 * slow probing). For this reason, it is checked in getFactorizedBinaryQuadraticExpr() whether \f$\sum_{i,j} Q_{ij} x_i x_j\f$
 * contains large (&ge; `reformbinprodsfac` parameter) lower sums of the form \f$x_i \sum_j Q_{ij} x_j\f$.
 * Such a lower sum is reformulated with only one extra variable w_i:
 * \f{align}{
 *    \text{maxact} & := \sum_j \max(0, Q_{ij}), \\
 *    \text{minact} & := \sum_j \min(0, Q_{ij}), \\
 *    \text{minact}\, x_i & \leq w_i, \\
 *    w_i &\leq \text{maxact}\, x_i, \\
 *    \text{minact} &\leq \sum_j Q_{ij} x_j - w_i + \text{minact}\, x_i \\
 *    \text{maxact} &\geq \sum_j Q_{ij} x_j - w_i + \text{maxact}\, x_i
 * \f}
 * We mark \f$w_i\f$ to be implicit integer if all \f$Q_{ij}\f$ are integer. After each replacement of a lower sum, it
 * is checked whether there are enough terms left to factorize other binary variables. Lower sums with a larger number
 * of terms are prioritized.
 */
static
SCIP_RETCODE presolveBinaryProducts(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS**           conss,              /**< constraints */
   int                   nconss,             /**< total number of constraints */
   int*                  naddconss,          /**< pointer to store the total number of added constraints (might be NULL) */
   int*                  nchgcoefs           /**< pointer to store the total number of changed coefficients (might be NULL) */
   )
{
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_HASHMAP* exprmap;
   SCIP_EXPRITER* it;
   int c;

   assert(conshdlr != NULL);

   /* no nonlinear constraints or binary variables -> skip */
   if( nconss == 0 || SCIPgetNBinVars(scip) == 0 )
      return SCIP_OKAY;
   assert(conss != NULL);

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);

   /* create expression hash map */
   SCIP_CALL( SCIPhashmapCreate(&exprmap, SCIPblkmem(scip), SCIPgetNVars(scip)) );

   /* create expression iterator */
   SCIP_CALL( SCIPcreateExpriter(scip, &it) );
   SCIP_CALL( SCIPexpriterInit(it, NULL, SCIP_EXPRITER_DFS, FALSE) );
   SCIPexpriterSetStagesDFS(it, SCIP_EXPRITER_VISITINGCHILD);

   SCIPdebugMsg(scip, "call presolveBinaryProducts()\n");

   for( c = 0; c < nconss; ++c )
   {
      SCIP_CONSDATA* consdata;
      SCIP_EXPR* newexpr = NULL;

      assert(conss[c] != NULL);

      consdata = SCIPconsGetData(conss[c]);
      assert(consdata != NULL);

      /* try to reformulate the root expression */
      if( conshdlrdata->reformbinprodsfac > 1 )
      {
         SCIP_CALL( getFactorizedBinaryQuadraticExpr(scip, conshdlr, conss[c], consdata->expr, conshdlrdata->reformbinprodsfac, &newexpr, naddconss) );
      }

      /* release the root node if another expression has been found */
      if( newexpr != NULL )
      {
         SCIP_CALL( SCIPreleaseExpr(scip, &consdata->expr) );
         consdata->expr = newexpr;

         /* mark constraint to be not simplified anymore */
         consdata->issimplified = FALSE;
      }

      /* replace each product of binary variables separately */
      SCIP_CALL( replaceBinaryProducts(scip, conshdlr, conss[c], exprmap, it, naddconss, nchgcoefs) );
   }

   /* free memory */
   SCIPhashmapFree(&exprmap);
   SCIPfreeExpriter(&it);

   return SCIP_OKAY;
}

/** scales the sides of the constraint \f$\ell \leq \sum_i c_i f_i(x) \leq r\f$.
 *
 *  Let \f$n_+\f$ the number of positive coefficients \f$c_i\f$ and \f$n_-\f$ be the number of negative coefficients.
 *  Then scale by -1 if
 *  - \f$n_+ < n_-\f$, or
 *  - \f$n_+ = n_-\f$ and \f$r = \infty\f$.
 */
static
SCIP_RETCODE scaleConsSides(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< nonlinear constraint handler */
   SCIP_CONS*            cons,               /**< nonlinear constraint */
   SCIP_Bool*            changed             /**< buffer to store if the expression of cons changed */
   )
{
   SCIP_CONSDATA* consdata;
   int i;

   assert(cons != NULL);

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);

   if( SCIPisExprSum(scip, consdata->expr) )
   {
      SCIP_Real* coefs;
      SCIP_Real constant;
      int nchildren;
      int counter = 0;

      coefs = SCIPgetCoefsExprSum(consdata->expr);
      constant = SCIPgetConstantExprSum(consdata->expr);
      nchildren = SCIPexprGetNChildren(consdata->expr);

      /* handle special case when constraint is l <= -f(x) <= r and f(x) not a sum: simplfy ensures f is not a sum */
      if( nchildren == 1 && constant == 0.0 && coefs[0] == -1.0 )
      {
         SCIP_EXPR* expr;
         expr = consdata->expr;

         consdata->expr = SCIPexprGetChildren(expr)[0];
         assert(!SCIPisExprSum(scip, consdata->expr));

         SCIPcaptureExpr(consdata->expr);

         SCIPswapReals(&consdata->lhs, &consdata->rhs);
         consdata->lhs = -consdata->lhs;
         consdata->rhs = -consdata->rhs;

         SCIP_CALL( SCIPreleaseExpr(scip, &expr) );
         *changed = TRUE;
         return SCIP_OKAY;
      }

      /* compute n_+ - n_i */
      for( i = 0; i < nchildren; ++i )
         counter += coefs[i] > 0 ? 1 : -1;

      if( counter < 0 || (counter == 0 && SCIPisInfinity(scip, consdata->rhs)) )
      {
         SCIP_EXPR* expr;
         SCIP_Real* newcoefs;

         /* allocate memory */
         SCIP_CALL( SCIPallocBufferArray(scip, &newcoefs, nchildren) );

         for( i = 0; i < nchildren; ++i )
            newcoefs[i] = -coefs[i];

         /* create a new sum expression */
         SCIP_CALL( SCIPcreateExprSum(scip, &expr, nchildren, SCIPexprGetChildren(consdata->expr), newcoefs, -constant, exprownerCreate, (void*)conshdlr) );

         /* replace expression in constraint data and scale sides */
         SCIP_CALL( SCIPreleaseExpr(scip, &consdata->expr) );
         consdata->expr = expr;
         SCIPswapReals(&consdata->lhs, &consdata->rhs);
         consdata->lhs = -consdata->lhs;
         consdata->rhs = -consdata->rhs;

         /* free memory */
         SCIPfreeBufferArray(scip, &newcoefs);

         *changed = TRUE;
      }
   }

   return SCIP_OKAY;
}

/** forbid multiaggrations of variables that appear nonlinear in constraints */
static
SCIP_RETCODE forbidNonlinearVariablesMultiaggration(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS**           conss,              /**< constraints */
   int                   nconss              /**< number of constraints */
   )
{
   SCIP_EXPRITER* it;
   SCIP_CONSDATA* consdata;
   SCIP_EXPR* expr;
   int c;

   assert(scip != NULL);
   assert(conshdlr != NULL);

   if( !SCIPconshdlrGetData(conshdlr)->forbidmultaggrnlvar )
      return SCIP_OKAY;

   SCIP_CALL( SCIPcreateExpriter(scip, &it) );
   SCIP_CALL( SCIPexpriterInit(it, NULL, SCIP_EXPRITER_DFS, FALSE) );

   for( c = 0; c < nconss; ++c )
   {
      consdata = SCIPconsGetData(conss[c]);
      assert(consdata != NULL);

      /* if root expression is sum, then forbid multiaggregation only for variables that are not in linear terms of sum,
       *   i.e., skip children of sum that are variables
       */
      if( SCIPisExprSum(scip, consdata->expr) )
      {
         int i;
         SCIP_EXPR* child;
         for( i = 0; i < SCIPexprGetNChildren(consdata->expr); ++i )
         {
            child = SCIPexprGetChildren(consdata->expr)[i];

            /* skip variable expression, as they correspond to a linear term */
            if( SCIPisExprVar(scip, child) )
               continue;

            for( expr = SCIPexpriterRestartDFS(it, child); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
               if( SCIPisExprVar(scip, expr) )
               {
                  SCIP_CALL( SCIPmarkDoNotMultaggrVar(scip, SCIPgetVarExprVar(expr)) );
               }
         }
      }
      else
      {
         for( expr = SCIPexpriterRestartDFS(it, consdata->expr); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
            if( SCIPisExprVar(scip, expr) )
            {
               SCIP_CALL( SCIPmarkDoNotMultaggrVar(scip, SCIPgetVarExprVar(expr)) );
            }
      }
   }

   SCIPfreeExpriter(&it);

   return SCIP_OKAY;
}

/** simplifies expressions and replaces common subexpressions for a set of constraints
 * @todo put the constant to the constraint sides
 */
static
SCIP_RETCODE canonicalizeConstraints(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS**           conss,              /**< constraints */
   int                   nconss,             /**< total number of constraints */
   SCIP_PRESOLTIMING     presoltiming,       /**< presolve timing (SCIP_PRESOLTIMING_ALWAYS if not in presolving) */
   SCIP_Bool*            infeasible,         /**< buffer to store whether infeasibility has been detected */
   int*                  ndelconss,          /**< counter to add number of deleted constraints, or NULL */
   int*                  naddconss,          /**< counter to add number of added constraints, or NULL */
   int*                  nchgcoefs           /**< counter to add number of changed coefficients, or NULL */
   )
{
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_CONSDATA* consdata;
   int* nlockspos;
   int* nlocksneg;
   SCIP_Bool havechange;
   int i;

   assert(scip != NULL);
   assert(conshdlr != NULL);
   assert(conss != NULL);
   assert(nconss > 0);
   assert(infeasible != NULL);

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);

   /* update number of canonicalize calls */
   ++(conshdlrdata->ncanonicalizecalls);

   SCIP_CALL( SCIPstartClock(scip, conshdlrdata->canonicalizetime) );

   *infeasible = FALSE;

   /* set havechange to TRUE in the first call of canonicalize; otherwise we might not replace common subexpressions */
   havechange = conshdlrdata->ncanonicalizecalls == 1;

   /* free nonlinear handlers information from expressions */  /* TODO can skip this in first presolve round */
   SCIP_CALL( deinitSolve(scip, conshdlr, conss, nconss) );

   /* allocate memory for storing locks of each constraint */
   SCIP_CALL( SCIPallocBufferArray(scip, &nlockspos, nconss) );
   SCIP_CALL( SCIPallocBufferArray(scip, &nlocksneg, nconss) );

   /* unlock all constraints */
   for( i = 0; i < nconss; ++i )
   {
      assert(conss[i] != NULL);

      consdata = SCIPconsGetData(conss[i]);
      assert(consdata != NULL);

      /* remember locks */
      nlockspos[i] = consdata->nlockspos;
      nlocksneg[i] = consdata->nlocksneg;

      /* remove locks */
      SCIP_CALL( addLocks(scip, conss[i], -consdata->nlockspos, -consdata->nlocksneg) );
      assert(consdata->nlockspos == 0);
      assert(consdata->nlocksneg == 0);
   }

#ifndef NDEBUG
   /* check whether all locks of each expression have been removed */
   for( i = 0; i < nconss; ++i )
   {
      SCIP_EXPR* expr;
      SCIP_EXPRITER* it;

      SCIP_CALL( SCIPcreateExpriter(scip, &it) );

      consdata = SCIPconsGetData(conss[i]);
      assert(consdata != NULL);

      SCIP_CALL( SCIPexpriterInit(it, consdata->expr, SCIP_EXPRITER_RTOPOLOGIC, TRUE) );
      for( expr = consdata->expr; !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
      {
         assert(expr != NULL);
         assert(SCIPexprGetOwnerData(expr)->nlocksneg == 0);
         assert(SCIPexprGetOwnerData(expr)->nlockspos == 0);
      }
      SCIPfreeExpriter(&it);
   }
#endif

   /* reformulate products of binary variables */
   if( conshdlrdata->reformbinprods && SCIPgetStage(scip) == SCIP_STAGE_PRESOLVING
      && (presoltiming & SCIP_PRESOLTIMING_EXHAUSTIVE) )
   {
      int tmpnaddconss = 0;
      int tmpnchgcoefs = 0;

      /* call this function before simplification because expressions might not be simplified after reformulating
       * binary products; the detection of some nonlinear handlers might assume that expressions are simplified
       */
      SCIP_CALL( presolveBinaryProducts(scip, conshdlr, conss, nconss, &tmpnaddconss, &tmpnchgcoefs) );

      /* update counters */
      if( naddconss != NULL )
         *naddconss = tmpnaddconss;
      if( nchgcoefs != NULL )
         *nchgcoefs = tmpnchgcoefs;

      /* check whether at least one expression has changed */
      if( tmpnaddconss + tmpnchgcoefs > 0 )
         havechange = TRUE;
   }

   for( i = 0; i < nconss; ++i )
   {
      consdata = SCIPconsGetData(conss[i]);
      assert(consdata != NULL);

      /* call simplify for each expression */
      if( !consdata->issimplified && consdata->expr != NULL )
      {
         SCIP_EXPR* simplified;
         SCIP_Bool changed;

         changed = FALSE;
         SCIP_CALL( SCIPsimplifyExpr(scip, consdata->expr, &simplified, &changed, infeasible, exprownerCreate, (void*)conshdlr) );
         consdata->issimplified = TRUE;

         if( changed )
            havechange = TRUE;

         /* If root expression changed, then we need to take care updating the locks as well (the consdata is the one holding consdata->expr "as a child").
          * If root expression did not change, some subexpression may still have changed, but the locks were taking care of in the corresponding SCIPreplaceExprChild() call.
          */
         if( simplified != consdata->expr )
         {
            assert(changed);

            /* release old expression */
            SCIP_CALL( SCIPreleaseExpr(scip, &consdata->expr) );

            /* store simplified expression */
            consdata->expr = simplified;
         }
         else
         {
            /* The simplify captures simplified in any case, also if nothing has changed.
             * Therefore, we have to release it here.
             */
            SCIP_CALL( SCIPreleaseExpr(scip, &simplified) );
         }

         if( *infeasible )
            break;

         /* scale constraint sides */
         SCIP_CALL( scaleConsSides(scip, conshdlr, conss[i], &changed) );

         if( changed )
            havechange = TRUE;

         /* handle constant root expression; either the problem is infeasible or the constraint is redundant */
         if( SCIPisExprValue(scip, consdata->expr) )
         {
            SCIP_Real value = SCIPgetValueExprValue(consdata->expr);
            if( (!SCIPisInfinity(scip, -consdata->lhs) && SCIPisFeasNegative(scip, value - consdata->lhs)) ||
                (!SCIPisInfinity(scip,  consdata->rhs) && SCIPisFeasPositive(scip, value - consdata->rhs)) )
            {
               SCIPdebugMsg(scip, "<%s> with constant expression found infeasible\n", SCIPconsGetName(conss[i]));
               SCIPdebugPrintCons(scip, conss[i], NULL);
               *infeasible = TRUE;
               break;
            }
            else
            {
               SCIP_CALL( addLocks(scip, conss[i], nlockspos[i], nlocksneg[i]) );
               SCIP_CALL( SCIPdelCons(scip, conss[i]) );
               if( ndelconss != NULL )
                  ++*ndelconss;
               havechange = TRUE;
            }
         }
      }
   }

   /* replace common subexpressions */
   if( havechange && !*infeasible )
   {
      SCIP_CONS** consssorted;
      SCIP_EXPR** rootexprs;
      SCIP_Bool replacedroot;

      SCIP_CALL( SCIPallocBufferArray(scip, &rootexprs, nconss) );
      for( i = 0; i < nconss; ++i )
         rootexprs[i] = SCIPconsGetData(conss[i])->expr;

      SCIP_CALL( SCIPreplaceCommonSubexpressions(scip, rootexprs, nconss, &replacedroot) );

      /* update pointer to root expr in constraints, if any has changed
       * SCIPreplaceCommonSubexpressions will have released the old expr and captures the new one
       */
      if( replacedroot )
         for( i = 0; i < nconss; ++i )
            SCIPconsGetData(conss[i])->expr = rootexprs[i];

      SCIPfreeBufferArray(scip, &rootexprs);

      /* TODO this is a possibly expensive way to update the variable expressions stored inside an expression which might have
       * been changed after simplification; now we completely recollect all variable expression and variable events
       */

      /* Each variable stores the constraints for which it catched varbound events sorted by the constraint index.
       * Thus, for performance reasons, it is better to call dropVarEvents in descending order of constraint index.
       */
      SCIP_CALL( SCIPduplicateBufferArray(scip, &consssorted, conss, nconss) );
      SCIPsortPtr((void**)consssorted, compIndexConsNonlinear, nconss);

      for( i = nconss-1; i >= 0; --i )
      {
         assert(i == 0 || compIndexConsNonlinear((void*)consssorted[i-1], (void*)consssorted[i]) < 0);
         if( SCIPconsIsDeleted(consssorted[i]) )
            continue;

         SCIP_CALL( dropVarEvents(scip, conshdlrdata->eventhdlr, consssorted[i]) );
         SCIP_CALL( freeVarExprs(scip, SCIPconsGetData(consssorted[i])) );
      }
      for( i = 0; i < nconss; ++i )
      {
         if( SCIPconsIsDeleted(consssorted[i]) )
            continue;

         SCIP_CALL( storeVarExprs(scip, conshdlr, SCIPconsGetData(consssorted[i])) );
         SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, consssorted[i]) );
      }

      SCIPfreeBufferArray(scip, &consssorted);

      /* forbid multiaggregation for nonlinear variables again (in case new variables appeared now)
       * a multiaggregation of a nonlinear variable can yield to a large increase in expressions due to
       * expanding terms in simplify, e.g. ,(sum_i x_i)^2, so we just forbid these
       */
      SCIP_CALL( forbidNonlinearVariablesMultiaggration(scip, conshdlr, conss, nconss) );
   }

   /* restore locks */
   for( i = 0; i < nconss; ++i )
   {
      if( SCIPconsIsDeleted(conss[i]) )
         continue;

      SCIP_CALL( addLocks(scip, conss[i], nlockspos[i], nlocksneg[i]) );
   }

   /* run nlhdlr detect if in presolving stage (that is, not in exitpre)
    * TODO can we skip this in presoltiming fast?
    */
   if( SCIPgetStage(scip) == SCIP_STAGE_PRESOLVING && !*infeasible )
   {
      /* reset one of the number of detections counter to count only current presolving round */
      for( i = 0; i < conshdlrdata->nnlhdlrs; ++i )
         SCIPnlhdlrResetNDetectionslast(conshdlrdata->nlhdlrs[i]);

      SCIP_CALL( initSolve(scip, conshdlr, conss, nconss) );
   }

   /* free allocated memory */
   SCIPfreeBufferArray(scip, &nlocksneg);
   SCIPfreeBufferArray(scip, &nlockspos);

   SCIP_CALL( SCIPstopClock(scip, conshdlrdata->canonicalizetime) );

   return SCIP_OKAY;
}

/** merges constraints that have the same root expression */
static
SCIP_RETCODE presolveMergeConss(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS**           conss,              /**< constraints to process */
   int                   nconss,             /**< number of constraints */
   SCIP_Bool*            success             /**< pointer to store whether at least one constraint could be deleted */
   )
{
   SCIP_HASHMAP* expr2cons;
   SCIP_Bool* updatelocks;
   int* nlockspos;
   int* nlocksneg;
   int c;

   assert(success != NULL);

   *success = FALSE;

   /* not enough constraints available */
   if( nconss <= 1 )
      return SCIP_OKAY;

   SCIP_CALL( SCIPhashmapCreate(&expr2cons, SCIPblkmem(scip), nconss) );
   SCIP_CALL( SCIPallocClearBufferArray(scip, &updatelocks, nconss) );
   SCIP_CALL( SCIPallocBufferArray(scip, &nlockspos, nconss) );
   SCIP_CALL( SCIPallocBufferArray(scip, &nlocksneg, nconss) );

   for( c = 0; c < nconss; ++c )
   {
      SCIP_CONSDATA* consdata;

      /* ignore deleted constraints */
      if( SCIPconsIsDeleted(conss[c]) )
         continue;

      consdata = SCIPconsGetData(conss[c]);
      assert(consdata != NULL);

      /* add expression to the hash map if not seen so far */
      if( !SCIPhashmapExists(expr2cons, (void*)consdata->expr) )
      {
         SCIP_CALL( SCIPhashmapInsertInt(expr2cons, (void*)consdata->expr, c) );
      }
      else
      {
         SCIP_CONSDATA* imgconsdata;
         int idx;

         idx = SCIPhashmapGetImageInt(expr2cons, (void*)consdata->expr);
         assert(idx >= 0 && idx < nconss);

         imgconsdata = SCIPconsGetData(conss[idx]);
         assert(imgconsdata != NULL);
         assert(imgconsdata->expr == consdata->expr);

         SCIPdebugMsg(scip, "merge constraint %g <= %s <= %g with %g <= %s <= %g\n", consdata->lhs,
            SCIPconsGetName(conss[c]), consdata->rhs, imgconsdata->lhs, SCIPconsGetName(conss[idx]), imgconsdata->rhs);

         /* check whether locks need to be updated */
         if( !updatelocks[idx] && ((SCIPisInfinity(scip, -imgconsdata->lhs) && !SCIPisInfinity(scip, -consdata->lhs))
            || (SCIPisInfinity(scip, imgconsdata->rhs) && !SCIPisInfinity(scip, consdata->rhs))) )
         {
            nlockspos[idx] = imgconsdata->nlockspos;
            nlocksneg[idx] = imgconsdata->nlocksneg;
            SCIP_CALL( addLocks(scip, conss[idx], -imgconsdata->nlockspos, -imgconsdata->nlocksneg) );
            updatelocks[idx] = TRUE;
         }

         /* update constraint sides */
         imgconsdata->lhs = MAX(imgconsdata->lhs, consdata->lhs);
         imgconsdata->rhs = MIN(imgconsdata->rhs, consdata->rhs);

         /* delete constraint */
         SCIP_CALL( SCIPdelCons(scip, conss[c]) );
         *success = TRUE;
      }
   }

   /* restore locks of updated constraints */
   if( *success )
   {
      for( c = 0; c < nconss; ++c )
      {
         if( updatelocks[c] )
         {
            SCIP_CALL( addLocks(scip, conss[c], nlockspos[c], nlocksneg[c]) );
         }
      }
   }

   /* free memory */
   SCIPfreeBufferArray(scip, &nlocksneg);
   SCIPfreeBufferArray(scip, &nlockspos);
   SCIPfreeBufferArray(scip, &updatelocks);
   SCIPhashmapFree(&expr2cons);

   return SCIP_OKAY;
}

/** interval evaluation of variables as used in redundancy check
 *
 * Returns local variable bounds of a variable, relaxed by feastol, as interval.
 */
static
SCIP_DECL_EXPR_INTEVALVAR(intEvalVarRedundancyCheck)
{  /*lint --e{715}*/
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_INTERVAL interval;
   SCIP_Real lb;
   SCIP_Real ub;

   assert(scip != NULL);
   assert(var != NULL);

   conshdlrdata = (SCIP_CONSHDLRDATA*)intevalvardata;
   assert(conshdlrdata != NULL);

   if( conshdlrdata->globalbounds )
   {
      lb = SCIPvarGetLbGlobal(var);
      ub = SCIPvarGetUbGlobal(var);
   }
   else
   {
      lb = SCIPvarGetLbLocal(var);
      ub = SCIPvarGetUbLocal(var);
   }
   assert(lb <= ub);  /* can SCIP ensure by now that variable bounds are not contradicting? */

   /* relax variable bounds, if there are bounds and variable is not fixed
    * (actually some assert complains if trying SCIPisRelEQ if both bounds are at different infinity)
    */
   if( !(SCIPisInfinity(scip, -lb) && SCIPisInfinity(scip, ub)) && !SCIPisRelEQ(scip, lb, ub) )
   {
      if( !SCIPisInfinity(scip, -lb) )
         lb -= SCIPfeastol(scip);

      if( !SCIPisInfinity(scip, ub) )
         ub += SCIPfeastol(scip);
   }

   /* convert SCIPinfinity() to SCIP_INTERVAL_INFINITY */
   lb = -infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, -lb);
   ub =  infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY,  ub);
   assert(lb <= ub);

   SCIPintervalSetBounds(&interval, lb, ub);

   return interval;
}

/** removes constraints that are always feasible or very simple
 *
 * Checks whether the activity of constraint functions is a subset of the constraint sides (relaxed by feastol).
 * To compute the activity, we use forwardPropExpr(), but relax variable bounds by feastol, because solutions to be checked
 * might violate variable bounds by up to feastol, too.
 * This is the main reason why the redundancy check is not done in propConss(), which relaxes variable bounds by epsilon only.
 *
 * Also removes constraints of the form lhs &le; variable &le; rhs.
 *
 * @todo it would be sufficient to check constraints for which we know that they are not currently violated by a valid solution
 *
 * @note This could should not run during solving, because the forwardProp takes the bounds of auxiliary variables into account.
 * For the root expression, these bounds are already set to the constraint sides, so that the activity of every expression
 * would appear as if the constraint is redundant.
 */
static
SCIP_RETCODE presolveRedundantConss(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   SCIP_CONS**           conss,              /**< constraints to propagate */
   int                   nconss,             /**< total number of constraints */
   SCIP_Bool*            cutoff,             /**< pointer to store whether infeasibility has been identified */
   int*                  ndelconss,          /**< buffer to add the number of deleted constraints */
   int*                  nchgbds             /**< buffer to add the number of variable bound tightenings */
   )
{
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_CONSDATA* consdata;
   SCIP_INTERVAL activity;
   SCIP_INTERVAL sides;
   int i;

   assert(scip != NULL);
   assert(conshdlr != NULL);
   assert(conss != NULL);
   assert(nconss >= 0);
   assert(cutoff != NULL);
   assert(ndelconss != NULL);
   assert(nchgbds != NULL);

   /* no constraints to check */
   if( nconss == 0 )
      return SCIP_OKAY;

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);

   /* increase curboundstag and set lastvaractivitymethodchange
    * we do this here to trigger a reevaluation of all variable bounds, since we will relax variable bounds
    * for the redundancy check differently than for domain propagation
    * we also update lastboundrelax to ensure activites of all expressions are indeed reevaluated
    */
   ++conshdlrdata->curboundstag;
   assert(conshdlrdata->curboundstag > 0);
   conshdlrdata->lastvaractivitymethodchange = conshdlrdata->curboundstag;
   conshdlrdata->lastboundrelax = conshdlrdata->curboundstag;
   conshdlrdata->intevalvar = intEvalVarRedundancyCheck;

   SCIPdebugMsg(scip, "checking %d constraints for redundancy\n", nconss);

   *cutoff = FALSE;
   for( i = 0; i < nconss; ++i )
   {
      if( !SCIPconsIsActive(conss[i]) || SCIPconsIsDeleted(conss[i]) )
         continue;

      consdata = SCIPconsGetData(conss[i]);
      assert(consdata != NULL);

      /* handle constant expressions separately: either the problem is infeasible or the constraint is redundant */
      if( SCIPisExprValue(scip, consdata->expr) )
      {
         SCIP_Real value = SCIPgetValueExprValue(consdata->expr);

         if( (!SCIPisInfinity(scip, -consdata->lhs) && value < consdata->lhs - SCIPfeastol(scip)) ||
             (!SCIPisInfinity(scip,  consdata->rhs) && value > consdata->rhs + SCIPfeastol(scip)) )
         {
            SCIPdebugMsg(scip, "constant constraint <%s> is infeasible: %g in [%g,%g] ", SCIPconsGetName(conss[i]), value, consdata->lhs, consdata->rhs);
            *cutoff = TRUE;

            goto TERMINATE;
         }

         SCIPdebugMsg(scip, "constant constraint <%s> is redundant: %g in [%g,%g] ", SCIPconsGetName(conss[i]), value, consdata->lhs, consdata->rhs);

         SCIP_CALL( SCIPdelConsLocal(scip, conss[i]) );
         ++*ndelconss;

         continue;
      }

      /* handle variable expressions separately: tighten variable bounds to constraint sides, then remove constraint (now redundant) */
      if( SCIPisExprVar(scip, consdata->expr) )
      {
         SCIP_VAR* var;
         SCIP_Bool tightened;

         var = SCIPgetVarExprVar(consdata->expr);
         assert(var != NULL);

         SCIPdebugMsg(scip, "variable constraint <%s> can be made redundant: <%s>[%g,%g] in [%g,%g]\n", SCIPconsGetName(conss[i]), SCIPvarGetName(var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var), consdata->lhs, consdata->rhs);

         /* ensure that variable bounds are within constraint sides */
         if( !SCIPisInfinity(scip, -consdata->lhs) )
         {
            SCIP_CALL( SCIPtightenVarLb(scip, var, consdata->lhs, TRUE, cutoff, &tightened) );

            if( tightened )
               ++*nchgbds;

            if( *cutoff )
               goto TERMINATE;
         }

         if( !SCIPisInfinity(scip, consdata->rhs) )
         {
            SCIP_CALL( SCIPtightenVarUb(scip, var, consdata->rhs, TRUE, cutoff, &tightened) );

            if( tightened )
               ++*nchgbds;

            if( *cutoff )
               goto TERMINATE;
         }

         /* delete the (now) redundant constraint locally */
         SCIP_CALL( SCIPdelConsLocal(scip, conss[i]) );
         ++*ndelconss;

         continue;
      }

      /* reevaluate expression activity, now using intEvalVarRedundancyCheck
       * we relax variable bounds by feastol here, as solutions that are checked later can also violate
       * variable bounds by up to feastol
       * (relaxing fixed variables seems to be too much, but they would be removed by presolve soon anyway)
       */
      SCIPdebugMsg(scip, "call forwardPropExpr() for constraint <%s>: ", SCIPconsGetName(conss[i]));
      SCIPdebugPrintCons(scip, conss[i], NULL);

      SCIP_CALL( forwardPropExpr(scip, conshdlr, consdata->expr, FALSE, cutoff, NULL) );
      assert(*cutoff || !SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(consdata->expr)));

      /* it is unlikely that we detect infeasibility by doing forward propagation */
      if( *cutoff )
      {
         SCIPdebugMsg(scip, " -> cutoff\n");
         goto TERMINATE;
      }

      assert(SCIPexprGetActivityTag(consdata->expr) == conshdlrdata->curboundstag);
      activity = SCIPexprGetActivity(consdata->expr);

      /* relax sides by feastol
       * we could accept every solution that violates constraints up to feastol as redundant, so this is the most permissive we can be
       */
      SCIPintervalSetBounds(&sides,
         SCIPisInfinity(scip, -consdata->lhs) ? -SCIP_INTERVAL_INFINITY : consdata->lhs - SCIPfeastol(scip),
         SCIPisInfinity(scip,  consdata->rhs) ?  SCIP_INTERVAL_INFINITY : consdata->rhs + SCIPfeastol(scip));

      if( SCIPintervalIsSubsetEQ(SCIP_INTERVAL_INFINITY, activity, sides) )
      {
         SCIPdebugMsg(scip, " -> redundant: activity [%g,%g] within sides [%g,%g]\n", activity.inf, activity.sup, consdata->lhs, consdata->rhs);

         SCIP_CALL( SCIPdelConsLocal(scip, conss[i]) );
         ++*ndelconss;

         continue;
      }

      SCIPdebugMsg(scip, " -> not redundant: activity [%g,%g] not within sides [%g,%g]\n", activity.inf, activity.sup, consdata->lhs, consdata->rhs);
   }

TERMINATE:
   /* make sure all activities are reevaluated again, since we relaxed bounds in a different way */
   ++conshdlrdata->curboundstag;
   conshdlrdata->lastvaractivitymethodchange = conshdlrdata->curboundstag;
   conshdlrdata->lastboundrelax = conshdlrdata->curboundstag;
   conshdlrdata->intevalvar = intEvalVarBoundTightening;

   return SCIP_OKAY;
}

/** tries to automatically convert a nonlinear constraint into a more specific and more specialized constraint */
static
SCIP_RETCODE presolveUpgrade(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler data structure */
   SCIP_CONS*            cons,               /**< source constraint to try to convert */
   SCIP_Bool*            upgraded,           /**< buffer to store whether constraint was upgraded */
   int*                  nupgdconss,         /**< buffer to increase if constraint was upgraded */
   int*                  naddconss           /**< buffer to increase with number of additional constraints created during upgrade */
   )
{
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_CONSDATA* consdata;
   SCIP_CONS** upgdconss;
   int upgdconsssize;
   int nupgdconss_;
   int i;

   assert(scip != NULL);
   assert(conshdlr != NULL);
   assert(cons != NULL);
   assert(!SCIPconsIsModifiable(cons));
   assert(upgraded   != NULL);
   assert(nupgdconss != NULL);
   assert(naddconss  != NULL);

   *upgraded = FALSE;

   nupgdconss_ = 0;

   conshdlrdata = SCIPconshdlrGetData(conshdlr);
   assert(conshdlrdata != NULL);

   /* if there are no upgrade methods, we can stop */
   if( conshdlrdata->nconsupgrades == 0 )
      return SCIP_OKAY;

   upgdconsssize = 2;
   SCIP_CALL( SCIPallocBufferArray(scip, &upgdconss, upgdconsssize) );

   /* call the upgrading methods */
   SCIPdebugMsg(scip, "upgrading nonlinear constraint <%s> (up to %d upgrade methods): ", SCIPconsGetName(cons), conshdlrdata->nconsupgrades);
   SCIPdebugPrintCons(scip, cons, NULL);

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);

   /* try all upgrading methods in priority order in case the upgrading step is enable  */
   for( i = 0; i < conshdlrdata->nconsupgrades; ++i )
   {
      if( !conshdlrdata->consupgrades[i]->active )
         continue;

      assert(conshdlrdata->consupgrades[i]->consupgd != NULL);

      SCIP_CALL( conshdlrdata->consupgrades[i]->consupgd(scip, cons, consdata->nvarexprs, &nupgdconss_, upgdconss, upgdconsssize) );

      while( nupgdconss_ < 0 )
      {
         /* upgrade function requires more memory: resize upgdconss and call again */
         assert(-nupgdconss_ > upgdconsssize);
         upgdconsssize = -nupgdconss_;
         SCIP_CALL( SCIPreallocBufferArray(scip, &upgdconss, -nupgdconss_) );

         SCIP_CALL( conshdlrdata->consupgrades[i]->consupgd(scip, cons, consdata->nvarexprs, &nupgdconss_, upgdconss, upgdconsssize) );

         assert(nupgdconss_ != 0);
      }

      if( nupgdconss_ > 0 )
      {
         /* got upgrade */
         int j;

         SCIPdebugMsg(scip, " -> upgraded to %d constraints:\n", nupgdconss_);

         /* add the upgraded constraints to the problem and forget them */
         for( j = 0; j < nupgdconss_; ++j )
         {
            SCIPdebugMsgPrint(scip, "\t");
            SCIPdebugPrintCons(scip, upgdconss[j], NULL);

            SCIP_CALL( SCIPaddCons(scip, upgdconss[j]) );      /*lint !e613*/
            SCIP_CALL( SCIPreleaseCons(scip, &upgdconss[j]) ); /*lint !e613*/
         }

         /* count the first upgrade constraint as constraint upgrade and the remaining ones as added constraints */
         *nupgdconss += 1;
         *naddconss += nupgdconss_ - 1;
         *upgraded = TRUE;

         /* delete upgraded constraint */
         SCIPdebugMsg(scip, "delete constraint <%s> after upgrade\n", SCIPconsGetName(cons));
         SCIP_CALL( SCIPdelCons(scip, cons) );

         break;
      }
   }

   SCIPfreeBufferArray(scip, &upgdconss);

   return SCIP_OKAY;
}

/** returns whether the variable of a given variable expression is a candidate for presolveSingleLockedVars(), i.e.,
 *  the variable is only contained in a single nonlinear constraint, has no objective coefficient, has finite
 *  variable bounds, and is not binary
 */
static
SCIP_Bool isSingleLockedCand(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EXPR*            expr                /**< variable expression */
   )
{
   SCIP_VAR* var;
   SCIP_EXPR_OWNERDATA* ownerdata;

   assert(SCIPisExprVar(scip, expr));

   var = SCIPgetVarExprVar(expr);
   assert(var != NULL);

   ownerdata = SCIPexprGetOwnerData(expr);
   assert(ownerdata != NULL);

   return SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) == ownerdata->nlocksneg
      && SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) == ownerdata->nlockspos
      && ownerdata->nconss == 1 && SCIPisZero(scip, SCIPvarGetObj(var))
      && !SCIPisInfinity(scip, -SCIPvarGetLbGlobal(var)) && !SCIPisInfinity(scip, SCIPvarGetUbGlobal(var))
      && SCIPvarGetType(var) != SCIP_VARTYPE_BINARY
      && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
}

/** removes all variable expressions that are contained in a given expression from a hash map */
static
SCIP_RETCODE removeSingleLockedVars(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EXPR*            expr,               /**< expression */
   SCIP_EXPRITER*        it,                 /**< expression iterator */
   SCIP_HASHMAP*         exprcands           /**< map to hash variable expressions */
   )
{
   SCIP_EXPR* e;

   for( e = SCIPexpriterRestartDFS(it, expr); !SCIPexpriterIsEnd(it); e = SCIPexpriterGetNext(it) )
   {
      if( SCIPisExprVar(scip, e) && SCIPhashmapExists(exprcands, (void*)e) )
      {
         SCIP_CALL( SCIPhashmapRemove(exprcands, (void*)e) );
      }
   }

   return SCIP_OKAY;
}

/** presolving method to fix a variable \f$x_i\f$ to one of its bounds if the variable is only contained in a single
 *  nonlinear constraint g(x) &le; rhs (&ge; lhs) if g() is concave (convex) in \f$x_i\f$
 *
 *  If a continuous variable has bounds [0,1], then the variable type is changed to be binary.
 *  Otherwise, a bound disjunction constraint is added.
 *
 *  @todo the same reduction can be applied if g(x) is not concave, but monotone in \f$x_i\f$ for g(x) &le; rhs
 *  @todo extend this to cases where a variable can appear in a monomial with an exponent, essentially relax