1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the program and library             */
4    	/*         SCIP --- Solving Constraint Integer Programs                      */
5    	/*                                                                           */
6    	/*  Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB)                      */
7    	/*                                                                           */
8    	/*  Licensed under the Apache License, Version 2.0 (the "License");          */
9    	/*  you may not use this file except in compliance with the License.         */
10   	/*  You may obtain a copy of the License at                                  */
11   	/*                                                                           */
12   	/*      http://www.apache.org/licenses/LICENSE-2.0                           */
13   	/*                                                                           */
14   	/*  Unless required by applicable law or agreed to in writing, software      */
15   	/*  distributed under the License is distributed on an "AS IS" BASIS,        */
16   	/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17   	/*  See the License for the specific language governing permissions and      */
18   	/*  limitations under the License.                                           */
19   	/*                                                                           */
20   	/*  You should have received a copy of the Apache-2.0 license                */
21   	/*  along with SCIP; see the file LICENSE. If not visit scipopt.org.         */
22   	/*                                                                           */
23   	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24   	
25   	/**@file   prop_nlobbt.c
26   	 * @ingroup DEFPLUGINS_PROP
27   	 * @brief  nlobbt propagator
28   	 * @author Benjamin Mueller
29   	 */
30   	
31   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
32   	
33   	#include "blockmemshell/memory.h"
34   	#include "scip/prop_genvbounds.h"
35   	#include "scip/prop_nlobbt.h"
36   	#include "scip/pub_message.h"
37   	#include "scip/pub_misc.h"
38   	#include "scip/pub_misc_sort.h"
39   	#include "scip/pub_nlp.h"
40   	#include "scip/pub_prop.h"
41   	#include "scip/pub_tree.h"
42   	#include "scip/pub_var.h"
43   	#include "scip/scip_general.h"
44   	#include "scip/scip_lp.h"
45   	#include "scip/scip_mem.h"
46   	#include "scip/scip_message.h"
47   	#include "scip/scip_nlp.h"
48   	#include "scip/scip_nlpi.h"
49   	#include "scip/scip_numerics.h"
50   	#include "scip/scip_param.h"
51   	#include "scip/scip_prob.h"
52   	#include "scip/scip_probing.h"
53   	#include "scip/scip_prop.h"
54   	#include "scip/scip_randnumgen.h"
55   	#include "scip/scip_solvingstats.h"
56   	#include "scip/scip_timing.h"
57   	#include "scip/scip_tree.h"
58   	#include "scip/scip_var.h"
59   	#include <string.h>
60   	
61   	#define PROP_NAME              "nlobbt"
62   	#define PROP_DESC              "propagator template"
63   	#define PROP_PRIORITY          -1100000
64   	#define PROP_FREQ                    -1
65   	#define PROP_DELAY                 TRUE
66   	#define PROP_TIMING            SCIP_PROPTIMING_AFTERLPLOOP
67   	
68   	#define DEFAULT_MINNONCONVEXFRAC   0.20      /**< default minimum (# convex nlrows)/(# nonconvex nlrows) threshold to apply propagator */
69   	#define DEFAULT_MINLINEARFRAC      0.02      /**< default minimum (# convex nlrows)/(# linear nlrows) threshold to apply propagator */
70   	#define DEFAULT_FEASTOLFAC         0.01      /**< default factor for NLP feasibility tolerance */
71   	#define DEFAULT_RELOBJTOLFAC       0.01      /**< default factor for NLP relative objective tolerance */
72   	#define DEFAULT_ADDLPROWS          TRUE      /**< should (non-initial) LP rows be used? */
73   	#define DEFAULT_ITLIMITFACTOR       2.0      /**< multiple of root node LP iterations used as total LP iteration
74   	                                              *   limit for nlobbt (<= 0: no limit ) */
75   	#define DEFAULT_NLPITERLIMIT        500      /**< default iteration limit of NLP solver; 0 for no limit */
76   	#define DEFAULT_NLPTIMELIMIT        0.0      /**< default time limit of NLP solver; 0.0 for no limit */
77   	#define DEFAULT_NLPVERLEVEL           0      /**< verbosity level of NLP solver */
78   	#define DEFAULT_RANDSEED             79      /**< initial random seed */
79   	
80   	/*
81   	 * Data structures
82   	 */
83   	
84   	/* status of bound candidates */
85   	#define UNSOLVED                      1      /**< did not solve LB or UB problem */
86   	#define SOLVEDLB                      2      /**< solved LB problem */
87   	#define SOLVEDUB                      4      /**< solved UB problem */
88   	
89   	/** propagator data */
90   	struct SCIP_PropData
91   	{
92   	   SCIP_NLPI*            nlpi;               /**< nlpi used to create the nlpi problem */
93   	   SCIP_NLPIPROBLEM*     nlpiprob;           /**< nlpi problem representing the convex NLP relaxation */
94   	   SCIP_HASHMAP*         var2nlpiidx;        /**< mapping between variables and nlpi indices */
95   	   SCIP_VAR**            nlpivars;           /**< array containing all variables of the nlpi */
96   	   int                   nlpinvars;          /**< total number of nlpi variables */
97   	   SCIP_Real*            nlscore;            /**< score for each nonlinear variable */
98   	   int*                  status;             /**< array containing a bound status for each candidate */
99   	   SCIP_PROP*            genvboundprop;      /**< genvbound propagator */
100  	   SCIP_RANDNUMGEN*      randnumgen;         /**< random number generator */
101  	   SCIP_Bool             skipprop;           /**< should the propagator be skipped? */
102  	   SCIP_Longint          lastnode;           /**< number of last node where obbt was performed */
103  	   int                   currpos;            /**< current position in the nlpivars array */
104  	
105  	   int                   nlpiterlimit;       /**< iteration limit of NLP solver; 0 for no limit */
106  	   SCIP_Real             nlptimelimit;       /**< time limit of NLP solver; 0.0 for no limit */
107  	   int                   nlpverblevel;       /**< verbosity level of NLP solver */
108  	   SCIP_NLPSTATISTICS    nlpstatistics;      /**< statistics from NLP solver */
109  	
110  	   SCIP_Real             feastolfac;         /**< factor for NLP feasibility tolerance */
111  	   SCIP_Real             relobjtolfac;       /**< factor for NLP relative objective tolerance */
112  	   SCIP_Real             minnonconvexfrac;   /**< minimum (#convex nlrows)/(#nonconvex nlrows) threshold to apply propagator */
113  	   SCIP_Real             minlinearfrac;      /**< minimum (#convex nlrows)/(#linear nlrows) threshold to apply propagator */
114  	   SCIP_Bool             addlprows;          /**< should (non-initial) LP rows be used? */
115  	   SCIP_Real             itlimitfactor;      /**< LP iteration limit for nlobbt will be this factor times total LP
116  	                                              *   iterations in root node */
117  	};
118  	
119  	/*
120  	 * Local methods
121  	 */
122  	
123  	/** clears the propagator data */
124  	static
125  	SCIP_RETCODE propdataClear(
126  	   SCIP*                 scip,               /**< SCIP data structure */
127  	   SCIP_PROPDATA*        propdata            /**< propagator data */
128  	   )
129  	{
130  	   assert(propdata != NULL);
131  	
132  	   if( propdata->nlpiprob != NULL )
133  	   {
134  	      assert(propdata->nlpi != NULL);
135  	
136  	      SCIPfreeBlockMemoryArray(scip, &propdata->status, propdata->nlpinvars);
137  	      SCIPfreeBlockMemoryArray(scip, &propdata->nlscore, propdata->nlpinvars);
138  	      SCIPfreeBlockMemoryArray(scip, &propdata->nlpivars, propdata->nlpinvars);
139  	      SCIPhashmapFree(&propdata->var2nlpiidx);
140  	      SCIP_CALL( SCIPfreeNlpiProblem(scip, propdata->nlpi, &propdata->nlpiprob) );
141  	
142  	      propdata->nlpinvars = 0;
143  	   }
144  	   assert(propdata->nlpinvars == 0);
145  	
146  	   propdata->skipprop = FALSE;
147  	   propdata->currpos = 0;
148  	   propdata->lastnode = -1;
149  	
150  	   return SCIP_OKAY;
151  	}
152  	
153  	/** checks whether it is worth to call nonlinear OBBT procedure */
154  	static
155  	SCIP_Bool isNlobbtApplicable(
156  	   SCIP*                 scip,               /**< SCIP data structure */
157  	   SCIP_PROPDATA*        propdata            /**< propagation data */
158  	   )
159  	{
160  	   SCIP_NLROW** nlrows;
161  	   int nnonconvexnlrows;
162  	   int nconvexnlrows;
163  	   int nlinearnlrows;
164  	   int nnlrows;
165  	   int i;
166  	
167  	   nlrows = SCIPgetNLPNlRows(scip);
168  	   nnlrows = SCIPgetNNLPNlRows(scip);
169  	   nnonconvexnlrows = 0;
170  	   nconvexnlrows = 0;
171  	   nlinearnlrows = 0;
172  	
173  	   for( i = 0; i < nnlrows; ++i )
174  	   {
175  	      if( SCIPnlrowGetExpr(nlrows[i]) == NULL )
176  	         ++nlinearnlrows;
177  	      else if( SCIPnlrowGetCurvature(nlrows[i]) == SCIP_EXPRCURV_CONVEX )
178  	      {
179  	         if( !SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrows[i])) )
180  	            ++nconvexnlrows;
181  	         if( !SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrows[i])) )
182  	            ++nnonconvexnlrows;
183  	      }
184  	      else if( SCIPnlrowGetCurvature(nlrows[i]) == SCIP_EXPRCURV_CONCAVE )
185  	      {
186  	         if( !SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrows[i])) )
187  	            ++nnonconvexnlrows;
188  	         if( !SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrows[i])) )
189  	            ++nconvexnlrows;
190  	      }
191  	      else
192  	      {
193  	         if( !SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrows[i])) )
194  	            ++nnonconvexnlrows;
195  	         if( !SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrows[i])) )
196  	            ++nnonconvexnlrows;
197  	      }
198  	   }
199  	
200  	   SCIPdebugMsg(scip, "nconvex=%d nnonconvex=%d nlinear=%d\n", nconvexnlrows, nnonconvexnlrows, nlinearnlrows);
201  	
202  	   return nconvexnlrows > 0
203  	      && nnonconvexnlrows > 0
204  	      && (SCIPisGE(scip, (SCIP_Real)nconvexnlrows, nnonconvexnlrows * propdata->minnonconvexfrac))
205  	      && (SCIPisGE(scip, (SCIP_Real)nconvexnlrows, nlinearnlrows * propdata->minlinearfrac));
206  	}
207  	
208  	/** filters variables which achieve their lower or dual bound in the current NLP solution */
209  	static
210  	SCIP_RETCODE filterCands(
211  	   SCIP*                 scip,               /**< SCIP data structure */
212  	   SCIP_PROPDATA*        propdata            /**< propagator data */
213  	   )
214  	{
215  	   SCIP_Real* primal;
216  	   int i;
217  	
218  	   assert(propdata->currpos >= 0 && propdata->currpos < propdata->nlpinvars);
219  	   assert(SCIPgetNlpiSolstat(scip, propdata->nlpi, propdata->nlpiprob) <= SCIP_NLPSOLSTAT_FEASIBLE);
220  	
221  	   SCIP_CALL( SCIPgetNlpiSolution(scip, propdata->nlpi, propdata->nlpiprob, &primal, NULL, NULL, NULL, NULL) );
222  	   assert(primal != NULL);
223  	
224  	   /* we skip all candidates which have been processed already, i.e., starting at propdata->currpos + 1 */
225  	   for( i = propdata->currpos + 1; i < propdata->nlpinvars; ++i )
226  	   {
227  	      SCIP_VAR* var;
228  	      SCIP_Real val;
229  	      int varidx;
230  	
231  	      /* only uninteresting variables left -> stop filtering */
232  	      if( SCIPisLE(scip, propdata->nlscore[i], 0.0) )
233  	         break;
234  	
235  	      var = propdata->nlpivars[i];
236  	      assert(var != NULL && SCIPhashmapExists(propdata->var2nlpiidx, (void*)var));
237  	
238  	      varidx = SCIPhashmapGetImageInt(propdata->var2nlpiidx, (void*)var);
239  	      assert(SCIPgetVars(scip)[varidx] == var);
240  	      val = primal[varidx];
241  	
242  	      if( (propdata->status[i] & SOLVEDLB) == 0 && !SCIPisInfinity(scip, -val) /*lint !e641*/
243  	         && SCIPisFeasLE(scip, val, SCIPvarGetLbLocal(var)) )
244  	      {
245  	         SCIPdebugMsg(scip, "filter LB of %s in [%g,%g] with %g\n", SCIPvarGetName(var), SCIPvarGetLbLocal(var),
246  	            SCIPvarGetUbLocal(var), val);
247  	         propdata->status[i] |= SOLVEDLB; /*lint !e641*/
248  	         assert((propdata->status[i] & SOLVEDLB) != 0); /*lint !e641*/
249  	      }
250  	
251  	      if( (propdata->status[i] & SOLVEDUB) == 0 && !SCIPisInfinity(scip, val) /*lint !e641*/
252  	         && SCIPisFeasGE(scip, val, SCIPvarGetUbLocal(var)) )
253  	      {
254  	         SCIPdebugMsg(scip, "filter UB of %s in [%g,%g] with %g\n", SCIPvarGetName(var), SCIPvarGetLbLocal(var),
255  	            SCIPvarGetUbLocal(var), val);
256  	         propdata->status[i] |= SOLVEDUB; /*lint !e641*/
257  	         assert((propdata->status[i] & SOLVEDUB) != 0); /*lint !e641*/
258  	      }
259  	   }
260  	
261  	   return SCIP_OKAY;
262  	}
263  	
264  	/** tries to add a generalized variable bound by exploiting the dual solution of the last NLP solve (see @ref
265  	 *  prop_nlobbt.h for more information)
266  	 */
267  	static
268  	SCIP_RETCODE addGenVBound(
269  	   SCIP*                 scip,               /**< SCIP data structure */
270  	   SCIP_PROPDATA*        propdata,           /**< propagator data */
271  	   SCIP_VAR*             var,                /**< variable used in last NLP solve */
272  	   int                   varidx,             /**< variable index in the propdata->nlpivars array */
273  	   SCIP_BOUNDTYPE        boundtype,          /**< type of bound provided by the genvbound */
274  	   SCIP_Real             cutoffbound         /**< cutoff bound */
275  	   )
276  	{
277  	   SCIP_VAR** lvbvars;
278  	   SCIP_Real* lvbcoefs;
279  	   SCIP_Real* primal;
280  	   SCIP_Real* dual;
281  	   SCIP_Real* alpha;
282  	   SCIP_Real* beta;
283  	   SCIP_Real constant;
284  	   SCIP_Real mu;
285  	   int nlvbvars;
286  	   int i;
287  	
288  	   assert(propdata->genvboundprop != NULL);
289  	   assert(var != NULL);
290  	   assert(varidx >= 0 && varidx < propdata->nlpinvars);
291  	   assert(SCIPgetNlpiSolstat(scip, propdata->nlpi, propdata->nlpiprob) <= SCIP_NLPSOLSTAT_LOCOPT);
292  	
293  	   SCIP_CALL( SCIPgetNlpiSolution(scip, propdata->nlpi, propdata->nlpiprob, &primal, &dual, &alpha, &beta, NULL) );
294  	
295  	   /* not possible to generate genvbound if the duals for the propagated variable do not disappear */
296  	   if( !SCIPisFeasZero(scip, alpha[varidx] - beta[varidx]) )
297  	      return SCIP_OKAY;
298  	
299  	   SCIP_CALL( SCIPallocBufferArray(scip, &lvbcoefs, propdata->nlpinvars) );
300  	   SCIP_CALL( SCIPallocBufferArray(scip, &lvbvars, propdata->nlpinvars) );
301  	   constant = boundtype == SCIP_BOUNDTYPE_LOWER ? primal[varidx] : -primal[varidx];
302  	   mu = 0.0;
303  	   nlvbvars = 0;
304  	
305  	   /* collect coefficients of genvbound */
306  	   for( i = 0; i < propdata->nlpinvars; ++i )
307  	   {
308  	      if( !SCIPisZero(scip, beta[i] - alpha[i]) )
309  	      {
310  	         lvbvars[nlvbvars] = propdata->nlpivars[i];
311  	         lvbcoefs[nlvbvars] = beta[i] - alpha[i];
312  	         ++nlvbvars;
313  	
314  	         constant += (alpha[i] - beta[i]) * primal[i];
315  	      }
316  	   }
317  	
318  	   /* first dual multiplier corresponds to the cutoff row if cutoffbound < SCIPinfinity() */
319  	   if( !SCIPisInfinity(scip, cutoffbound) && SCIPisGT(scip, dual[0], 0.0) )
320  	   {
321  	      mu = dual[0];
322  	      constant += mu * cutoffbound;
323  	   }
324  	
325  	   /* add genvbound to genvbounds propagator */
326  	   if( !SCIPisInfinity(scip, REALABS(constant)) && (nlvbvars > 0 || SCIPisFeasGT(scip, mu, 0.0)) )
327  	   {
328  	      SCIP_CALL( SCIPgenVBoundAdd(scip, propdata->genvboundprop, lvbvars, var, lvbcoefs, nlvbvars, -mu, constant,
329  	            boundtype) );
330  	      SCIPdebugMsg(scip, "add genvbound for %s\n", SCIPvarGetName(var));
331  	   }
332  	
333  	   SCIPfreeBufferArray(scip, &lvbvars);
334  	   SCIPfreeBufferArray(scip, &lvbcoefs);
335  	
336  	   return SCIP_OKAY;
337  	}
338  	
339  	/** sets the objective function, solves the NLP, and tightens the given variable; after calling this function, the
340  	 *  objective function is set to zero
341  	 *
342  	 *  @note function assumes that objective function is zero
343  	 */
344  	static
345  	SCIP_RETCODE solveNlp(
346  	   SCIP*                 scip,               /**< SCIP data structure */
347  	   SCIP_PROPDATA*        propdata,           /**< propagator data */
348  	   SCIP_VAR*             var,                /**< variable to propagate */
349  	   int                   varidx,             /**< variable index in the propdata->nlpivars array */
350  	   SCIP_BOUNDTYPE        boundtype,          /**< minimize or maximize var? */
351  	   SCIP_NLPPARAM*        nlpparam,           /**< NLP solve parameters */
352  	   int*                  nlpiter,            /**< buffer to store the total number of nlp iterations */
353  	   SCIP_RESULT*          result              /**< pointer to store result */
354  	   )
355  	{
356  	   SCIP_Real timelimit;
357  	   SCIP_Real* primal;
358  	   SCIP_Real obj;
359  	   int iterlimit;
360  	
361  	#ifdef SCIP_DEBUG
362  	   SCIP_Real oldlb;
363  	   SCIP_Real oldub;
364  	
365  	   oldlb = SCIPvarGetLbLocal(var);
366  	   oldub = SCIPvarGetUbLocal(var);
367  	#endif
368  	
369  	   assert(var != NULL);
370  	   assert(varidx >= 0 && varidx < propdata->nlpinvars);
371  	   assert(result != NULL && *result != SCIP_CUTOFF);
372  	
373  	   *nlpiter = 0;
374  	
375  	   /* set time and iteration limit */
376  	   SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) );
377  	   if( !SCIPisInfinity(scip, timelimit) )
378  	   {
379  	      timelimit -= SCIPgetSolvingTime(scip);
380  	      if( timelimit <= 0.0 )
381  	      {
382  	         SCIPdebugMsg(scip, "skip NLP solve; no time left\n");
383  	         return SCIP_OKAY;
384  	      }
385  	   }
386  	   if( propdata->nlptimelimit > 0.0 )
387  	      timelimit = MIN(propdata->nlptimelimit, timelimit);
388  	   iterlimit = propdata->nlpiterlimit > 0 ? propdata->nlpiterlimit : INT_MAX;
389  	   nlpparam->timelimit = timelimit;
390  	   nlpparam->iterlimit = iterlimit;
391  	
392  	   /* set corresponding objective coefficient and solve NLP */
393  	   obj = boundtype == SCIP_BOUNDTYPE_LOWER ? 1.0 : -1.0;
394  	   SCIP_CALL( SCIPsetNlpiObjective(scip, propdata->nlpi, propdata->nlpiprob, 1, &varidx, &obj, NULL, 0.0) );
395  	
396  	   SCIPdebugMsg(scip, "solve var=%s boundtype=%d nlscore=%g\n", SCIPvarGetName(var), boundtype,
397  	      propdata->nlscore[propdata->currpos]);
398  	   SCIP_CALL( SCIPsolveNlpiParam(scip, propdata->nlpi, propdata->nlpiprob, *nlpparam) );
399  	   SCIPdebugMsg(scip, "NLP solstat = %d\n", SCIPgetNlpiSolstat(scip, propdata->nlpi, propdata->nlpiprob));
400  	
401  	   /* collect NLP statistics */
402  	   SCIP_CALL( SCIPgetNlpiStatistics(scip, propdata->nlpi, propdata->nlpiprob, &propdata->nlpstatistics) );
403  	   *nlpiter = propdata->nlpstatistics.niterations;
404  	   SCIPdebugMsg(scip, "iterations %d time %g\n", *nlpiter, propdata->nlpstatistics.totaltime);
405  	
406  	   /* filter bound candidates first, otherwise we do not have access to the primal solution values */
407  	   if( SCIPgetNlpiSolstat(scip, propdata->nlpi, propdata->nlpiprob) <= SCIP_NLPSOLSTAT_FEASIBLE )
408  	   {
409  	      SCIP_CALL( filterCands(scip, propdata) );
410  	   }
411  	
412  	   /* try to tighten variable bound */
413  	   if( SCIPgetNlpiSolstat(scip, propdata->nlpi, propdata->nlpiprob) <= SCIP_NLPSOLSTAT_LOCOPT )
414  	   {
415  	      SCIP_Bool tightened;
416  	      SCIP_Bool infeasible;
417  	
418  	      /* try to add a genvbound in the root node */
419  	      if( propdata->genvboundprop != NULL && SCIPgetDepth(scip) == 0 )
420  	      {
421  	         SCIP_CALL( addGenVBound(scip, propdata, var, varidx, boundtype, SCIPgetCutoffbound(scip)) );
422  	      }
423  	
424  	      SCIP_CALL( SCIPgetNlpiSolution(scip, propdata->nlpi, propdata->nlpiprob, &primal, NULL, NULL, NULL, NULL) );
425  	
426  	      if( boundtype == SCIP_BOUNDTYPE_LOWER )
427  	      {
428  	         SCIP_CALL( SCIPtightenVarLb(scip, var, primal[varidx], FALSE, &infeasible, &tightened) );
429  	      }
430  	      else
431  	      {
432  	         SCIP_CALL( SCIPtightenVarUb(scip, var, primal[varidx], FALSE, &infeasible, &tightened) );
433  	      }
434  	
435  	      if( infeasible )
436  	      {
437  	         SCIPdebugMsg(scip, "detect infeasibility after propagating %s\n", SCIPvarGetName(var));
438  	         *result = SCIP_CUTOFF;
439  	      }
440  	      else if( tightened )
441  	      {
442  	         SCIP_Real lb;
443  	         SCIP_Real ub;
444  	
445  	         *result = SCIP_REDUCEDDOM;
446  	
447  	         /* update bounds in NLP */
448  	         lb = SCIPvarGetLbLocal(var);
449  	         ub = SCIPvarGetUbLocal(var);
450  	         SCIP_CALL( SCIPchgNlpiVarBounds(scip, propdata->nlpi, propdata->nlpiprob, 1, &varidx, &lb, &ub) );
451  	
452  	#ifdef SCIP_DEBUG
453  	         SCIPdebugMsg(scip, "tightened bounds of %s from [%g,%g] to [%g,%g]\n", SCIPvarGetName(var), oldlb, oldub, lb, ub);
454  	#endif
455  	      }
456  	   }
457  	
458  	   /* reset objective function */
459  	   obj = 0.0;
460  	   SCIP_CALL( SCIPsetNlpiObjective(scip, propdata->nlpi, propdata->nlpiprob, 1, &varidx, &obj, NULL, 0.0) );
461  	
462  	   return SCIP_OKAY;
463  	}
464  	
465  	/** main method of the propagator
466  	 *
467  	 *  creates a convex NLP relaxation and solves the OBBT-NLPs for each possible candidate;
468  	 *  binary and variables with a small domain will be ignored to reduce the computational cost of the propagator; after
469  	 *  solving each NLP we filter out all variable candidates which are on their lower or upper bound; candidates with a
470  	 *  larger number of occurrences are preferred
471  	 */
472  	static
473  	SCIP_RETCODE applyNlobbt(
474  	   SCIP*                 scip,               /**< SCIP data structure */
475  	   SCIP_PROPDATA*        propdata,           /**< propagation data */
476  	   SCIP_RESULT*          result              /**< pointer to store result */
477  	   )
478  	{
479  	   SCIP_NLPPARAM nlpparam = SCIP_NLPPARAM_DEFAULT(scip);  /*lint !e446*/
480  	   int nlpiterleft;
481  	
482  	   assert(result != NULL);
483  	   assert(!propdata->skipprop);
484  	   assert(SCIPgetNNlpis(scip) > 0);
485  	
486  	   *result = SCIP_DIDNOTRUN;
487  	
488  	   if( propdata->nlpiprob == NULL && !isNlobbtApplicable(scip, propdata) )
489  	   {
490  	      /* do not call the propagator anymore (except after a restart) */
491  	      SCIPdebugMsg(scip, "nlobbt propagator is not applicable\n");
492  	      propdata->skipprop = TRUE;
493  	      return SCIP_OKAY;
494  	   }
495  	
496  	   *result = SCIP_DIDNOTFIND;
497  	
498  	   /* compute NLP iteration limit */
499  	   if( propdata->itlimitfactor > 0.0 )
500  	      nlpiterleft = (int)(propdata->itlimitfactor * SCIPgetNRootLPIterations(scip));
501  	   else
502  	      nlpiterleft = INT_MAX;
503  	
504  	   /* recompute NLP relaxation if the variable set changed */
505  	   if( propdata->nlpiprob != NULL && SCIPgetNVars(scip) != propdata->nlpinvars )
506  	   {
507  	      SCIP_CALL( propdataClear(scip, propdata) );
508  	      assert(propdata->nlpiprob == NULL);
509  	   }
510  	
511  	   /* create or update NLP relaxation */
512  	   if( propdata->nlpiprob == NULL )
513  	   {
514  	      int i;
515  	
516  	      propdata->nlpinvars = SCIPgetNVars(scip);
517  	      propdata->nlpi = SCIPgetNlpis(scip)[0];
518  	      assert(propdata->nlpi != NULL);
519  	
520  	      SCIP_CALL( SCIPhashmapCreate(&propdata->var2nlpiidx, SCIPblkmem(scip), propdata->nlpinvars) );
521  	      SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &propdata->nlpivars, SCIPgetVars(scip), propdata->nlpinvars) ); /*lint !e666*/
522  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &propdata->nlscore, propdata->nlpinvars) );
523  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &propdata->status, propdata->nlpinvars) );
524  	
525  	      SCIP_CALL( SCIPcreateNlpiProblemFromNlRows(scip, propdata->nlpi, &propdata->nlpiprob, "nlobbt-nlp", SCIPgetNLPNlRows(scip), SCIPgetNNLPNlRows(scip),
526  	            propdata->var2nlpiidx, NULL, propdata->nlscore, SCIPgetCutoffbound(scip), FALSE, TRUE) );
527  	
528  	      /* initialize bound status; perturb nlscores by a factor which ensures that zero scores remain zero */
529  	      assert(propdata->randnumgen != NULL);
530  	      for( i = 0; i < propdata->nlpinvars; ++i )
531  	      {
532  	         propdata->status[i] = UNSOLVED; /*lint !e641*/
533  	         propdata->nlscore[i] *= 1.0 + SCIPrandomGetReal(propdata->randnumgen, SCIPfeastol(scip), 2.0 * SCIPfeastol(scip));
534  	      }
535  	
536  	      /* add rows of the LP */
537  	      if( SCIPgetDepth(scip) == 0 )
538  	      {
539  	         SCIP_CALL( SCIPaddNlpiProblemRows(scip, propdata->nlpi, propdata->nlpiprob, propdata->var2nlpiidx,
540  	               SCIPgetLPRows(scip), SCIPgetNLPRows(scip)) );
541  	      }
542  	   }
543  	   else
544  	   {
545  	      SCIP_CALL( SCIPupdateNlpiProblem(scip, propdata->nlpi, propdata->nlpiprob, propdata->var2nlpiidx,
546  	            propdata->nlpivars, propdata->nlpinvars, SCIPgetCutoffbound(scip)) );
547  	   }
548  	
549  	   assert(propdata->nlpiprob != NULL);
550  	   assert(propdata->var2nlpiidx != NULL);
551  	   assert(propdata->nlpivars != NULL);
552  	   assert(propdata->nlscore != NULL);
553  	
554  	   /* sort variables w.r.t. their nlscores if we did not solve any NLP for this node */
555  	   if( propdata->currpos == 0 )
556  	   {
557  	      SCIPsortDownRealIntPtr(propdata->nlscore, propdata->status, (void**)propdata->nlpivars, propdata->nlpinvars);
558  	   }
559  	
560  	   /* set parameters of NLP solver */
561  	   nlpparam.feastol *= propdata->feastolfac;
562  	   nlpparam.opttol = SCIPfeastol(scip) * propdata->relobjtolfac;
563  	   nlpparam.verblevel = (unsigned short)propdata->nlpverblevel;
564  	
565  	   /* main propagation loop */
566  	   while( propdata->currpos < propdata->nlpinvars
567  	      && nlpiterleft > 0
568  	      && SCIPisGT(scip, propdata->nlscore[propdata->currpos], 0.0)
569  	      && *result != SCIP_CUTOFF
570  	      && !SCIPisStopped(scip) )
571  	   {
572  	      SCIP_VAR* var;
573  	      int varidx;
574  	      int iters;
575  	
576  	      var = propdata->nlpivars[propdata->currpos];
577  	      assert(var != NULL);
578  	
579  	      /* skip binary or almost fixed variables */
580  	      if( SCIPvarGetType(var) == SCIP_VARTYPE_BINARY
581  	         || SCIPisFeasEQ(scip, SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var)) )
582  	      {
583  	         ++(propdata->currpos);
584  	         continue;
585  	      }
586  	
587  	      SCIPdebugMsg(scip, "iterations left %d\n", nlpiterleft);
588  	
589  	      /* get index of var in the nlpi */
590  	      assert(SCIPhashmapExists(propdata->var2nlpiidx, (void*)var) );
591  	      varidx = SCIPhashmapGetImageInt(propdata->var2nlpiidx, (void*)var);
592  	      assert(var == SCIPgetVars(scip)[varidx]);
593  	
594  	      /* case: minimize var */
595  	      if( (propdata->status[propdata->currpos] & SOLVEDLB) == 0 ) /*lint !e641*/
596  	      {
597  	         SCIP_CALL( solveNlp(scip, propdata, var, varidx, SCIP_BOUNDTYPE_LOWER, &nlpparam, &iters, result) );
598  	         nlpiterleft -= iters;
599  	      }
600  	
601  	      /* case: maximize var */
602  	      if( *result != SCIP_CUTOFF && (propdata->status[propdata->currpos] & SOLVEDUB) == 0 ) /*lint !e641*/
603  	      {
604  	         SCIP_CALL( solveNlp(scip, propdata, var, varidx, SCIP_BOUNDTYPE_UPPER, &nlpparam, &iters, result) );
605  	         nlpiterleft -= iters;
606  	      }
607  	
608  	      /* update the current position */
609  	      ++(propdata->currpos);
610  	   }
611  	
612  	   return SCIP_OKAY;
613  	}
614  	
615  	/*
616  	 * Callback methods of propagator
617  	 */
618  	
619  	/** destructor of propagator to free user data (called when SCIP is exiting) */
620  	static
621  	SCIP_DECL_PROPFREE(propFreeNlobbt)
622  	{  /*lint --e{715}*/
623  	   SCIP_PROPDATA* propdata;
624  	
625  	   propdata = SCIPpropGetData(prop);
626  	   assert(propdata != NULL);
627  	
628  	   SCIP_CALL( propdataClear(scip, propdata) );
629  	   SCIPfreeBlockMemory(scip, &propdata);
630  	   SCIPpropSetData(prop, NULL);
631  	
632  	   return SCIP_OKAY;
633  	}
634  	
635  	/** solving process initialization method of propagator (called when branch and bound process is about to begin) */
636  	static
637  	SCIP_DECL_PROPINITSOL(propInitsolNlobbt)
638  	{  /*lint --e{715}*/
639  	   SCIP_PROPDATA* propdata;
640  	
641  	   assert(scip != NULL);
642  	   assert(prop != NULL);
643  	
644  	   propdata = SCIPpropGetData(prop);
645  	   assert(propdata != NULL);
646  	
647  	   /* if genvbounds propagator is not available, we cannot create genvbounds */
648  	   propdata->genvboundprop = SCIPfindProp(scip, "genvbounds");
649  	
650  	   SCIP_CALL( SCIPcreateRandom(scip, &propdata->randnumgen,
651  	         DEFAULT_RANDSEED, TRUE) );
652  	   propdata->lastnode = -1;
653  	
654  	   return SCIP_OKAY;
655  	}
656  	
657  	/** solving process deinitialization method of propagator (called before branch and bound process data is freed) */
658  	static
659  	SCIP_DECL_PROPEXITSOL(propExitsolNlobbt)
660  	{  /*lint --e{715}*/
661  	   SCIP_PROPDATA* propdata;
662  	
663  	   propdata = SCIPpropGetData(prop);
664  	   assert(propdata != NULL);
665  	
666  	   SCIPfreeRandom(scip, &propdata->randnumgen);
667  	
668  	   SCIP_CALL( propdataClear(scip, propdata) );
669  	
670  	   return SCIP_OKAY;
671  	}
672  	
673  	/** execution method of propagator */
674  	static
675  	SCIP_DECL_PROPEXEC(propExecNlobbt)
676  	{  /*lint --e{715}*/
677  	   SCIP_PROPDATA* propdata;
678  	
679  	   *result = SCIP_DIDNOTRUN;
680  	
681  	   propdata = SCIPpropGetData(prop);
682  	   assert(propdata != NULL);
683  	
684  	   if( propdata->skipprop || SCIPgetStage(scip) != SCIP_STAGE_SOLVING || SCIPinRepropagation(scip)
685  	      || SCIPinProbing(scip) || SCIPinDive(scip) || !SCIPallowWeakDualReds(scip) || SCIPgetNNlpis(scip) == 0 )
686  	   {
687  	      SCIPdebugMsg(scip, "skip nlobbt propagator\n");
688  	      return SCIP_OKAY;
689  	   }
690  	
691  	   /* only run if LP all columns are in the LP, e.g., do not run if pricers are active */
692  	   if( !SCIPallColsInLP(scip) )
693  	   {
694  	      SCIPdebugMsg(scip, "not all columns in LP, skipping obbt\n");
695  	      return SCIP_OKAY;
696  	   }
697  	
698  	   /* do not run if SCIP does not have constructed an NLP */
699  	   if( !SCIPisNLPConstructed(scip) )
700  	   {
701  	      SCIPdebugMsg(scip, "NLP not constructed, skipping nlobbt\n");
702  	      return SCIP_OKAY;
703  	   }
704  	
705  	   /* consider all variables again if we process a new node */
706  	   if( SCIPnodeGetNumber(SCIPgetCurrentNode(scip)) != propdata->lastnode )
707  	   {
708  	      propdata->lastnode = SCIPnodeGetNumber(SCIPgetCurrentNode(scip));
709  	      propdata->currpos = 0;
710  	   }
711  	
712  	   /* call main procedure of nonlinear OBBT propagator */
713  	   SCIP_CALL( applyNlobbt(scip, propdata, result) );
714  	
715  	   return SCIP_OKAY;
716  	}
717  	
718  	/*
719  	 * propagator specific interface methods
720  	 */
721  	
722  	/** creates the nlobbt propagator and includes it in SCIP */
723  	SCIP_RETCODE SCIPincludePropNlobbt(
724  	   SCIP*                 scip                /**< SCIP data structure */
725  	   )
726  	{
727  	   SCIP_PROPDATA* propdata;
728  	   SCIP_PROP* prop;
729  	
730  	   propdata = NULL;
731  	   prop = NULL;
732  	
733  	   SCIP_CALL( SCIPallocBlockMemory(scip, &propdata) );
734  	   assert(propdata != NULL);
735  	   BMSclearMemory(propdata);
736  	
737  	   SCIP_CALL( SCIPincludePropBasic(scip, &prop, PROP_NAME, PROP_DESC, PROP_PRIORITY, PROP_FREQ, PROP_DELAY, PROP_TIMING,
738  	         propExecNlobbt, propdata) );
739  	   assert(prop != NULL);
740  	
741  	   SCIP_CALL( SCIPsetPropFree(scip, prop, propFreeNlobbt) );
742  	   SCIP_CALL( SCIPsetPropInitsol(scip, prop, propInitsolNlobbt) );
743  	   SCIP_CALL( SCIPsetPropExitsol(scip, prop, propExitsolNlobbt) );
744  	
745  	   SCIP_CALL( SCIPaddRealParam(scip, "propagating/" PROP_NAME "/feastolfac",
746  	         "factor for NLP feasibility tolerance",
747  	         &propdata->feastolfac, TRUE, DEFAULT_FEASTOLFAC, 0.0, 1.0, NULL, NULL) );
748  	
749  	   SCIP_CALL( SCIPaddRealParam(scip, "propagating/" PROP_NAME "/relobjtolfac",
750  	         "factor for NLP relative objective tolerance",
751  	         &propdata->relobjtolfac, TRUE, DEFAULT_RELOBJTOLFAC, 0.0, 1.0, NULL, NULL) );
752  	
753  	   SCIP_CALL( SCIPaddRealParam(scip, "propagating/" PROP_NAME "/minnonconvexfrac",
754  	         "(#convex nlrows)/(#nonconvex nlrows) threshold to apply propagator",
755  	         &propdata->minnonconvexfrac, TRUE, DEFAULT_MINNONCONVEXFRAC, 0.0, SCIPinfinity(scip), NULL, NULL) );
756  	
757  	   SCIP_CALL( SCIPaddRealParam(scip, "propagating/" PROP_NAME "/minlinearfrac",
758  	         "minimum (#convex nlrows)/(#linear nlrows) threshold to apply propagator",
759  	         &propdata->minlinearfrac, TRUE, DEFAULT_MINLINEARFRAC, 0.0, SCIPinfinity(scip), NULL, NULL) );
760  	
761  	   SCIP_CALL( SCIPaddBoolParam(scip, "propagating/" PROP_NAME "/addlprows",
762  	         "should non-initial LP rows be used?",
763  	         &propdata->addlprows, FALSE, DEFAULT_ADDLPROWS, NULL, NULL) );
764  	
765  	   SCIP_CALL( SCIPaddIntParam(scip, "propagating/" PROP_NAME "/nlpiterlimit",
766  	         "iteration limit of NLP solver; 0 for no limit",
767  	         &propdata->nlpiterlimit, TRUE, DEFAULT_NLPITERLIMIT, 0, INT_MAX, NULL, NULL) );
768  	
769  	   SCIP_CALL( SCIPaddRealParam(scip, "propagating/" PROP_NAME "/nlptimelimit",
770  	         "time limit of NLP solver; 0.0 for no limit",
771  	         &propdata->nlptimelimit, TRUE, DEFAULT_NLPTIMELIMIT, 0.0, SCIP_REAL_MAX, NULL, NULL) );
772  	
773  	   SCIP_CALL( SCIPaddIntParam(scip, "propagating/" PROP_NAME "/nlpverblevel",
774  	         "verbosity level of NLP solver",
775  	         &propdata->nlpverblevel, TRUE, DEFAULT_NLPVERLEVEL, 0, 5, NULL, NULL) );
776  	
777  	   SCIP_CALL( SCIPaddRealParam(scip, "propagating/" PROP_NAME "/itlimitfactor",
778  	         "LP iteration limit for nlobbt will be this factor times total LP iterations in root node",
779  	         &propdata->itlimitfactor, TRUE, DEFAULT_ITLIMITFACTOR, 0.0, SCIP_REAL_MAX, NULL, NULL) );
780  	
781  	   return SCIP_OKAY;
782  	}
783