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   branch_distribution.c
26   	 * @ingroup DEFPLUGINS_BRANCH
27   	 * @ingroup BRANCHINGRULES
28   	 * @brief  probability based branching rule based on an article by J. Pryor and J.W. Chinneck
29   	 * @author Gregor Hendel
30   	 *
31   	 * The distribution branching rule selects a variable based on its impact on row activity distributions. More formally,
32   	 * let \f$ a(x) = a_1 x_1 + \dots + a_n x_n \leq b \f$ be a row of the LP. Let further \f$ l_i, u_i \in R\f$ denote the
33   	 * (finite) lower and upper bound, respectively, of the \f$ i \f$-th variable \f$x_i\f$.
34   	 * Viewing every variable \f$x_i \f$ as (continuously) uniformly distributed within its bounds, we can approximately
35   	 * understand the row activity \f$a(x)\f$ as a Gaussian random variate with mean value \f$ \mu = E[a(x)] = \sum_i a_i\frac{l_i + u_i}{2}\f$
36   	 * and variance \f$ \sigma^2 = \sum_i a_i^2 \sigma_i^2 \f$, with \f$ \sigma_i^2 = \frac{(u_i - l_i)^2}{12}\f$ for
37   	 * continuous and \f$ \sigma_i^2 = \frac{(u_i - l_i + 1)^2 - 1}{12}\f$ for discrete variables.
38   	 * With these two parameters, we can calculate the probability to satisfy the row in terms of the cumulative distribution
39   	 * of the standard normal distribution: \f$ P(a(x) \leq b) = \Phi(\frac{b - \mu}{\sigma})\f$.
40   	 *
41   	 * The impact of a variable on the probability to satisfy a constraint after branching can be estimated by altering
42   	 * the variable contribution to the sums described above. In order to keep the tree size small,
43   	 * variables are preferred which decrease the probability
44   	 * to satisfy a row because it is more likely to create infeasible subproblems.
45   	 *
46   	 * The selection of the branching variable is subject to the parameter @p scoreparam. For both branching directions,
47   	 * an individual score is calculated. Available options for scoring methods are:
48   	 * - @b d: select a branching variable with largest difference in satisfaction probability after branching
49   	 * - @b l: lowest cumulative probability amongst all variables on all rows (after branching).
50   	 * - @b h: highest cumulative probability amongst all variables on all rows (after branching).
51   	 * - @b v: highest number of votes for lowering row probability for all rows a variable appears in.
52   	 * - @b w: highest number of votes for increasing row probability for all rows a variable appears in.
53   	 *
54   	 * If the parameter @p usescipscore is set to @a TRUE, a single branching score is calculated from the respective
55   	 * up and down scores as defined by the SCIP branching score function (see advanced branching parameter @p scorefunc),
56   	 * otherwise, the variable with the single highest score is selected, and the maximizing direction is assigned
57   	 * higher branching priority.
58   	 *
59   	 * The original idea of probability based branching schemes appeared in:
60   	 *
61   	 * J. Pryor and J.W. Chinneck:@n
62   	 * Faster Integer-Feasibility in Mixed-Integer Linear Programs by Branching to Force Change@n
63   	 * Computers and Operations Research, vol. 38, 2011, p. 1143–1152@n
64   	 * (Paper: <a href="http://www.sce.carleton.ca/faculty/chinneck/docs/PryorChinneck.pdf">PDF</a>).
65   	 */
66   	
67   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
68   	
69   	#include "scip/branch_distribution.h"
70   	#include "scip/pub_branch.h"
71   	#include "scip/pub_event.h"
72   	#include "scip/pub_lp.h"
73   	#include "scip/pub_message.h"
74   	#include "scip/pub_misc.h"
75   	#include "scip/pub_var.h"
76   	#include "scip/scip_branch.h"
77   	#include "scip/scip_event.h"
78   	#include "scip/scip_general.h"
79   	#include "scip/scip_lp.h"
80   	#include "scip/scip_message.h"
81   	#include "scip/scip_mem.h"
82   	#include "scip/scip_numerics.h"
83   	#include "scip/scip_param.h"
84   	#include "scip/scip_pricer.h"
85   	#include "scip/scip_prob.h"
86   	#include "scip/scip_probing.h"
87   	#include <string.h>
88   	
89   	
90   	#define BRANCHRULE_NAME            "distribution"
91   	#define BRANCHRULE_DESC            "branching rule based on variable influence on cumulative normal distribution of row activities"
92   	#define BRANCHRULE_PRIORITY        0
93   	#define BRANCHRULE_MAXDEPTH        -1
94   	#define BRANCHRULE_MAXBOUNDDIST    1.0
95   	
96   	#define SCOREPARAM_VALUES          "dhlvw"
97   	#define DEFAULT_SCOREPARAM         'v'
98   	#define DEFAULT_PRIORITY           2.0
99   	#define SQRTOFTWO                  1.4142136
100  	#define SQUARED(x) ((x) * (x))
101  	#define DEFAULT_ONLYACTIVEROWS     FALSE     /**< should only rows which are active at the current node be considered? */
102  	#define DEFAULT_USEWEIGHTEDSCORE   FALSE     /**< should the branching score weigh up- and down-scores of a variable */
103  	
104  	/* branching rule event handler to catch variable bound changes */
105  	#define EVENTHDLR_NAME             "eventhdlr_distribution"    /**< event handler name */
106  	#define EVENT_DISTRIBUTION         SCIP_EVENTTYPE_BOUNDCHANGED /**< the event tyoe to be handled by this event handler */
107  	
108  	/*
109  	 * Data structures
110  	 */
111  	
112  	/** branching rule data */
113  	struct SCIP_BranchruleData
114  	{
115  	   SCIP_EVENTHDLR*       eventhdlr;          /**< event handler pointer */
116  	   SCIP_VAR**            updatedvars;        /**< variables to process bound change events for */
117  	   SCIP_Real*            rowmeans;           /**< row activity mean values for all rows */
118  	   SCIP_Real*            rowvariances;       /**< row activity variances for all rows */
119  	   SCIP_Real*            currentubs;         /**< variable upper bounds as currently saved in the row activities */
120  	   SCIP_Real*            currentlbs;         /**< variable lower bounds as currently saved in the row activities */
121  	   int*                  rowinfinitiesdown;  /**< count the number of variables with infinite bounds which allow for always
122  	                                              *   repairing the constraint right hand side */
123  	   int*                  rowinfinitiesup;    /**< count the number of variables with infinite bounds which allow for always
124  	                                              *   repairing the constraint left hand side */
125  	   int*                  varposs;            /**< array of variable positions in the updated variables array */
126  	   int*                  varfilterposs;      /**< array of event filter positions for variable events */
127  	   int                   nupdatedvars;       /**< the current number of variables with pending bound changes */
128  	   int                   memsize;            /**< memory size of current arrays, needed for dynamic reallocation */
129  	   int                   varpossmemsize;     /**< memory size of updated vars and varposs array */
130  	   char                  scoreparam;         /**< parameter how the branch score is calculated */
131  	   SCIP_Bool             onlyactiverows;     /**< should only rows which are active at the current node be considered? */
132  	   SCIP_Bool             usescipscore;       /**< should the branching use SCIP's branching score function */
133  	};
134  	
135  	struct SCIP_EventhdlrData
136  	{
137  	   SCIP_BRANCHRULEDATA*  branchruledata;     /**< the branching rule data to access distribution arrays */
138  	};
139  	
140  	/*
141  	 * Local methods
142  	 */
143  	
144  	/** ensure that maxindex + 1 rows can be represented in data arrays; memory gets reallocated with 10% extra space
145  	 *  to save some time for future allocations */
146  	static
147  	SCIP_RETCODE branchruledataEnsureArraySize(
148  	   SCIP*                 scip,               /**< SCIP data structure */
149  	   SCIP_BRANCHRULEDATA*  branchruledata,     /**< branchruledata */
150  	   int                   maxindex            /**< row index at hand (size must be at least this large) */
151  	   )
152  	{
153  	   int newsize;
154  	   int r;
155  	
156  	   /* maxindex fits in current array -> nothing to do */
157  	   if( maxindex < branchruledata->memsize )
158  	      return SCIP_OKAY;
159  	
160  	   /* new memory size is the max index + 1 plus 10% additional space */
161  	   newsize = (int)SCIPfeasCeil(scip, (maxindex + 1) * 1.1);
162  	   assert(newsize > branchruledata->memsize);
163  	   assert(branchruledata->memsize >= 0);
164  	
165  	   /* alloc memory arrays for row information */
166  	   if( branchruledata->memsize == 0 )
167  	   {
168  	      SCIP_VAR** vars;
169  	      int v;
170  	      int nvars;
171  	
172  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->rowinfinitiesdown, newsize) );
173  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->rowinfinitiesup, newsize) );
174  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->rowmeans, newsize) );
175  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->rowvariances, newsize) );
176  	
177  	      assert(SCIPgetStage(scip) == SCIP_STAGE_SOLVING);
178  	
179  	      vars = SCIPgetVars(scip);
180  	      nvars = SCIPgetNVars(scip);
181  	
182  	      assert(nvars > 0);
183  	
184  	      /* allocate variable update event processing array storage */
185  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->varfilterposs, nvars) );
186  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->varposs, nvars) );
187  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->updatedvars, nvars) );
188  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->currentubs, nvars) );
189  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->currentlbs, nvars) );
190  	
191  	      branchruledata->varpossmemsize = nvars;
192  	      branchruledata->nupdatedvars = 0;
193  	
194  	      /* init variable event processing data */
195  	      for( v = 0; v < nvars; ++v )
196  	      {
197  	         assert(SCIPvarIsActive(vars[v]));
198  	         assert(SCIPvarGetProbindex(vars[v]) == v);
199  	
200  	         /* set up variable events to catch bound changes */
201  	         SCIP_CALL( SCIPcatchVarEvent(scip, vars[v], EVENT_DISTRIBUTION, branchruledata->eventhdlr, NULL, &(branchruledata->varfilterposs[v])) );
202  	         assert(branchruledata->varfilterposs[v] >= 0);
203  	
204  	         branchruledata->varposs[v] = -1;
205  	         branchruledata->updatedvars[v] = NULL;
206  	         branchruledata->currentlbs[v] = SCIP_INVALID;
207  	         branchruledata->currentubs[v] = SCIP_INVALID;
208  	      }
209  	   }
210  	   else
211  	   {
212  	      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &branchruledata->rowinfinitiesdown, branchruledata->memsize, newsize) );
213  	      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &branchruledata->rowinfinitiesup, branchruledata->memsize, newsize) );
214  	      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &branchruledata->rowmeans, branchruledata->memsize, newsize) );
215  	      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &branchruledata->rowvariances, branchruledata->memsize, newsize) );
216  	   }
217  	
218  	   /* loop over extended arrays and invalidate data to trigger initialization of this row when necessary */
219  	   for( r = branchruledata->memsize; r < newsize; ++r )
220  	   {
221  	      branchruledata->rowmeans[r] = SCIP_INVALID;
222  	      branchruledata->rowvariances[r] = SCIP_INVALID;
223  	      branchruledata->rowinfinitiesdown[r] = 0;
224  	      branchruledata->rowinfinitiesup[r] = 0;
225  	   }
226  	
227  	   /* adjust memsize */
228  	   branchruledata->memsize = newsize;
229  	
230  	   return SCIP_OKAY;
231  	}
232  	
233  	/* update the variables current lower and upper bound */
234  	static
235  	void branchruledataUpdateCurrentBounds(
236  	   SCIP*                 scip,               /**< SCIP data structure */
237  	   SCIP_BRANCHRULEDATA*  branchruledata,     /**< branchrule data */
238  	   SCIP_VAR*             var                 /**< the variable to update current bounds */
239  	   )
240  	{
241  	   int varindex;
242  	   SCIP_Real lblocal;
243  	   SCIP_Real ublocal;
244  	
245  	   assert(var != NULL);
246  	
247  	   varindex = SCIPvarGetProbindex(var);
248  	   assert(0 <= varindex && varindex < branchruledata->varpossmemsize);
249  	   lblocal = SCIPvarGetLbLocal(var);
250  	   ublocal = SCIPvarGetUbLocal(var);
251  	
252  	   assert(SCIPisFeasLE(scip, lblocal, ublocal));
253  	
254  	   branchruledata->currentlbs[varindex] = lblocal;
255  	   branchruledata->currentubs[varindex] = ublocal;
256  	}
257  	
258  	/** calculate the variable's distribution parameters (mean and variance) for the bounds specified in the arguments.
259  	 *  special treatment of infinite bounds necessary */
260  	void SCIPvarCalcDistributionParameters(
261  	   SCIP*                 scip,               /**< SCIP data structure */
262  	   SCIP_Real             varlb,              /**< variable lower bound */
263  	   SCIP_Real             varub,              /**< variable upper bound */
264  	   SCIP_VARTYPE          vartype,            /**< type of the variable */
265  	   SCIP_Real*            mean,               /**< pointer to store mean value */
266  	   SCIP_Real*            variance            /**< pointer to store the variance of the variable uniform distribution */
267  	   )
268  	{
269  	   assert(mean != NULL);
270  	   assert(variance != NULL);
271  	
272  	   /* calculate mean and variance of variable uniform distribution before and after branching */
273  	   if( SCIPisInfinity(scip, varub) || SCIPisInfinity(scip, -varlb) )
274  	   {
275  	      /* variables with infinite bounds are not kept in the row activity variance */
276  	      *variance = 0.0;
277  	
278  	      if( !SCIPisInfinity(scip, varub) )
279  	         *mean = varub;
280  	      else if( !SCIPisInfinity(scip, -varlb) )
281  	         *mean = varlb;
282  	      else
283  	         *mean = 0.0;
284  	   }
285  	   else
286  	   {
287  	      /* if the variable is continuous, we assume a continuous uniform distribution, otherwise a discrete one */
288  	      if( vartype == SCIP_VARTYPE_CONTINUOUS )
289  	         *variance = SQUARED(varub - varlb);
290  	      else
291  	         *variance = SQUARED(varub - varlb + 1) - 1;
292  	      *variance /= 12.0;
293  	      *mean = (varub + varlb) * .5;
294  	   }
295  	
296  	   assert(!SCIPisNegative(scip, *variance));
297  	}
298  	
299  	/** calculates the cumulative distribution P(-infinity <= x <= value) that a normally distributed
300  	 *  random variable x takes a value between -infinity and parameter \p value.
301  	 *
302  	 *  The distribution is given by the respective mean and deviation. This implementation
303  	 *  uses the error function SCIPerf().
304  	 */
305  	SCIP_Real SCIPcalcCumulativeDistribution(
306  	   SCIP*                 scip,               /**< current SCIP */
307  	   SCIP_Real             mean,               /**< the mean value of the distribution */
308  	   SCIP_Real             variance,           /**< the square of the deviation of the distribution */
309  	   SCIP_Real             value               /**< the upper limit of the calculated distribution integral */
310  	   )
311  	{
312  	   SCIP_Real normvalue;
313  	   SCIP_Real std;
314  	
315  	   /* we need to calculate the standard deviation from the variance */
316  	   assert(!SCIPisNegative(scip, variance));
317  	   if( SCIPisFeasZero(scip, variance) )
318  	      std = 0.0;
319  	   else
320  	      std = sqrt(variance);
321  	
322  	   /* special treatment for zero variance */
323  	   if( SCIPisFeasZero(scip, std) )
324  	   {
325  	      if( SCIPisFeasLE(scip, value, mean) )
326  	         return 1.0;
327  	      else
328  	         return 0.0;
329  	   }
330  	   assert( std != 0.0 ); /* for lint */
331  	
332  	   /* scale and translate to standard normal distribution. Factor sqrt(2) is needed for SCIPerf() function */
333  	   normvalue = (value - mean)/(std * SQRTOFTWO);
334  	
335  	   SCIPdebugMsg(scip, " Normalized value %g = ( %g - %g ) / (%g * 1.4142136)\n", normvalue, value, mean, std);
336  	
337  	   /* calculate the cumulative distribution function for normvalue. For negative normvalues, we negate
338  	    * the normvalue and use the oddness of the SCIPerf()-function; special treatment for values close to zero.
339  	    */
340  	   if( SCIPisFeasZero(scip, normvalue) )
341  	      return .5;
342  	   else if( normvalue > 0 )
343  	   {
344  	      SCIP_Real erfresult;
345  	
346  	      erfresult = SCIPerf(normvalue);
347  	      return  erfresult / 2.0 + 0.5;
348  	   }
349  	   else
350  	   {
351  	      SCIP_Real erfresult;
352  	
353  	      erfresult = SCIPerf(-normvalue);
354  	
355  	      return 0.5 - erfresult / 2.0;
356  	   }
357  	}
358  	
359  	/** calculates the probability of satisfying an LP-row under the assumption
360  	 *  of uniformly distributed variable values.
361  	 *
362  	 *  For inequalities, we use the cumulative distribution function of the standard normal
363  	 *  distribution PHI(rhs - mu/sqrt(sigma2)) to calculate the probability
364  	 *  for a right hand side row with mean activity mu and variance sigma2 to be satisfied.
365  	 *  Similarly, 1 - PHI(lhs - mu/sqrt(sigma2)) is the probability to satisfy a left hand side row.
366  	 *  For equations (lhs==rhs), we use the centeredness measure p = min(PHI(lhs'), 1-PHI(lhs'))/max(PHI(lhs'), 1 - PHI(lhs')),
367  	 *  where lhs' = lhs - mu / sqrt(sigma2).
368  	 */
369  	SCIP_Real SCIProwCalcProbability(
370  	   SCIP*                 scip,               /**< current scip */
371  	   SCIP_ROW*             row,                /**< the row */
372  	   SCIP_Real             mu,                 /**< the mean value of the row distribution */
373  	   SCIP_Real             sigma2,             /**< the variance of the row distribution */
374  	   int                   rowinfinitiesdown,  /**< the number of variables with infinite bounds to DECREASE activity */
375  	   int                   rowinfinitiesup     /**< the number of variables with infinite bounds to INCREASE activity */
376  	   )
377  	{
378  	   SCIP_Real rowprobability;
379  	   SCIP_Real lhs;
380  	   SCIP_Real rhs;
381  	   SCIP_Real lhsprob;
382  	   SCIP_Real rhsprob;
383  	
384  	   lhs = SCIProwGetLhs(row);
385  	   rhs = SCIProwGetRhs(row);
386  	
387  	   lhsprob = 1.0;
388  	   rhsprob = 1.0;
389  	
390  	   /* use the cumulative distribution if the row contains no variable to repair every infeasibility */
391  	   if( !SCIPisInfinity(scip, rhs) && rowinfinitiesdown == 0 )
392  	      rhsprob = SCIPcalcCumulativeDistribution(scip, mu, sigma2, rhs);
393  	
394  	   /* use the cumulative distribution if the row contains no variable to repair every infeasibility
395  	    * otherwise the row can always be made feasible by increasing activity far enough
396  	    */
397  	   if( !SCIPisInfinity(scip, -lhs) && rowinfinitiesup == 0 )
398  	      lhsprob = 1.0 - SCIPcalcCumulativeDistribution(scip, mu, sigma2, lhs);
399  	
400  	   assert(SCIPisFeasLE(scip, lhsprob, 1.0) && SCIPisFeasGE(scip, lhsprob, 0.0));
401  	   assert(SCIPisFeasLE(scip, rhsprob, 1.0) && SCIPisFeasGE(scip, rhsprob, 0.0));
402  	
403  	   /* use centeredness measure for equations; for inequalities, the minimum of the two probabilities is the
404  	    * probability to satisfy the row */
405  	   if( SCIPisFeasEQ(scip, lhs, rhs) )
406  	   {
407  	      SCIP_Real minprobability;
408  	      SCIP_Real maxprobability;
409  	
410  	      minprobability = MIN(rhsprob, lhsprob);
411  	      maxprobability = MAX(lhsprob, rhsprob);
412  	      rowprobability = minprobability / maxprobability;
413  	   }
414  	   else
415  	      rowprobability = MIN(rhsprob, lhsprob);
416  	
417  	   SCIPdebug( SCIPprintRow(scip, row, NULL) );
418  	   SCIPdebugMsg(scip, " Row %s, mean %g, sigma2 %g, LHS %g, RHS %g has probability %g to be satisfied\n",
419  	      SCIProwGetName(row), mu, sigma2, lhs, rhs, rowprobability);
420  	
421  	   assert(SCIPisFeasGE(scip, rowprobability, 0.0) && SCIPisFeasLE(scip, rowprobability, 1.0));
422  	
423  	   return rowprobability;
424  	}
425  	
426  	/** calculates the initial mean and variance of the row activity normal distribution.
427  	 *
428  	 *  The mean value \f$ \mu \f$ is given by \f$ \mu = \sum_i=1^n c_i * (lb_i +ub_i) / 2 \f$ where
429  	 *  \f$n \f$ is the number of variables, and \f$ c_i, lb_i, ub_i \f$ are the variable coefficient and
430  	 *  bounds, respectively. With the same notation, the variance \f$ \sigma^2 \f$ is given by
431  	 *  \f$ \sigma^2 = \sum_i=1^n c_i^2 * \sigma^2_i \f$, with the variance being
432  	 *  \f$ \sigma^2_i = ((ub_i - lb_i + 1)^2 - 1) / 12 \f$ for integer variables and
433  	 *  \f$ \sigma^2_i = (ub_i - lb_i)^2 / 12 \f$ for continuous variables.
434  	 */
435  	static
436  	void rowCalculateGauss(
437  	   SCIP*                 scip,               /**< SCIP data structure */
438  	   SCIP_BRANCHRULEDATA*  branchruledata,     /**< the branching rule data */
439  	   SCIP_ROW*             row,                /**< the row for which the gaussian normal distribution has to be calculated */
440  	   SCIP_Real*            mu,                 /**< pointer to store the mean value of the gaussian normal distribution */
441  	   SCIP_Real*            sigma2,             /**< pointer to store the variance value of the gaussian normal distribution */
442  	   int*                  rowinfinitiesdown,  /**< pointer to store the number of variables with infinite bounds to DECREASE activity */
443  	   int*                  rowinfinitiesup     /**< pointer to store the number of variables with infinite bounds to INCREASE activity */
444  	   )
445  	{
446  	   SCIP_COL** rowcols;
447  	   SCIP_Real* rowvals;
448  	   int nrowvals;
449  	   int c;
450  	
451  	   assert(scip != NULL);
452  	   assert(row != NULL);
453  	   assert(mu != NULL);
454  	   assert(sigma2 != NULL);
455  	   assert(rowinfinitiesup != NULL);
456  	   assert(rowinfinitiesdown != NULL);
457  	
458  	   rowcols = SCIProwGetCols(row);
459  	   rowvals = SCIProwGetVals(row);
460  	   nrowvals = SCIProwGetNNonz(row);
461  	
462  	   assert(nrowvals == 0 || rowcols != NULL);
463  	   assert(nrowvals == 0 || rowvals != NULL);
464  	
465  	   *mu = SCIProwGetConstant(row);
466  	   *sigma2 = 0.0;
467  	   *rowinfinitiesdown = 0;
468  	   *rowinfinitiesup = 0;
469  	
470  	   /* loop over nonzero row coefficients and sum up the variable contributions to mu and sigma2 */
471  	   for( c = 0; c < nrowvals; ++c )
472  	   {
473  	      SCIP_VAR* colvar;
474  	      SCIP_Real colval;
475  	      SCIP_Real colvarlb;
476  	      SCIP_Real colvarub;
477  	      SCIP_Real squarecoeff;
478  	      SCIP_Real varvariance;
479  	      SCIP_Real varmean;
480  	      int varindex;
481  	
482  	      assert(rowcols[c] != NULL);
483  	      colvar = SCIPcolGetVar(rowcols[c]);
484  	      assert(colvar != NULL);
485  	
486  	      colval = rowvals[c];
487  	      colvarlb = SCIPvarGetLbLocal(colvar);
488  	      colvarub = SCIPvarGetUbLocal(colvar);
489  	
490  	      varmean = 0.0;
491  	      varvariance = 0.0;
492  	      varindex = SCIPvarGetProbindex(colvar);
493  	      assert((branchruledata->currentlbs[varindex] == SCIP_INVALID) == (branchruledata->currentubs[varindex] == SCIP_INVALID)); /*lint !e777*/
494  	
495  	      /* variable bounds need to be watched from now on */
496  	      if( branchruledata->currentlbs[varindex] == SCIP_INVALID ) /*lint !e777*/
497  	         branchruledataUpdateCurrentBounds(scip, branchruledata, colvar);
498  	
499  	      assert(!SCIPisInfinity(scip, colvarlb));
500  	      assert(!SCIPisInfinity(scip, -colvarub));
501  	      assert(SCIPisFeasLE(scip, colvarlb, colvarub));
502  	
503  	      /* variables with infinite bounds are skipped for the calculation of the variance; they need to
504  	       * be accounted for by the counters for infinite row activity decrease and increase and they
505  	       * are used to shift the row activity mean in case they have one nonzero, but finite bound */
506  	      if( SCIPisInfinity(scip, -colvarlb) || SCIPisInfinity(scip, colvarub) )
507  	      {
508  	         if( SCIPisInfinity(scip, colvarub) )
509  	         {
510  	         /* an infinite upper bound gives the row an infinite maximum activity or minimum activity, if the coefficient is
511  	          * positive or negative, resp.
512  	          */
513  	            if( colval < 0.0 )
514  	               ++(*rowinfinitiesdown);
515  	            else
516  	               ++(*rowinfinitiesup);
517  	         }
518  	
519  	         /* an infinite lower bound gives the row an infinite maximum activity or minimum activity, if the coefficient is
520  	          * negative or positive, resp.
521  	          */
522  	         if( SCIPisInfinity(scip, -colvarlb) )
523  	         {
524  	            if( colval > 0.0 )
525  	               ++(*rowinfinitiesdown);
526  	            else
527  	               ++(*rowinfinitiesup);
528  	         }
529  	      }
530  	      SCIPvarCalcDistributionParameters(scip, colvarlb, colvarub, SCIPvarGetType(colvar), &varmean, &varvariance);
531  	
532  	      /* actual values are updated; the contribution of the variable to mu is the arithmetic mean of its bounds */
533  	      *mu += colval * varmean;
534  	
535  	      /* the variance contribution of a variable is c^2 * (u - l)^2 / 12.0 for continuous and c^2 * ((u - l + 1)^2 - 1) / 12.0 for integer */
536  	      squarecoeff = SQUARED(colval);
537  	      *sigma2 += squarecoeff * varvariance;
538  	
539  	      assert(!SCIPisFeasNegative(scip, *sigma2));
540  	   }
541  	
542  	   SCIPdebug( SCIPprintRow(scip, row, NULL) );
543  	   SCIPdebugMsg(scip, "  Row %s has a mean value of %g at a sigma2 of %g \n", SCIProwGetName(row), *mu, *sigma2);
544  	}
545  	
546  	/** update the up- and downscore of a single variable after calculating the impact of branching on a
547  	 *  particular row, depending on the chosen score parameter
548  	 */
549  	SCIP_RETCODE SCIPupdateDistributionScore(
550  	   SCIP*                 scip,               /**< current SCIP pointer */
551  	   SCIP_Real             currentprob,        /**< the current probability */
552  	   SCIP_Real             newprobup,          /**< the new probability if branched upwards */
553  	   SCIP_Real             newprobdown,        /**< the new probability if branched downwards */
554  	   SCIP_Real*            upscore,            /**< pointer to store the new score for branching up */
555  	   SCIP_Real*            downscore,          /**< pointer to store the new score for branching down */
556  	   char                  scoreparam          /**< parameter to determine the way the score is calculated */
557  	   )
558  	{
559  	   assert(scip != NULL);
560  	   assert(SCIPisFeasGE(scip, currentprob, 0.0) && SCIPisFeasLE(scip, currentprob, 1.0));
561  	   assert(SCIPisFeasGE(scip, newprobup, 0.0) && SCIPisFeasLE(scip, newprobup, 1.0));
562  	   assert(SCIPisFeasGE(scip, newprobdown, 0.0) && SCIPisFeasLE(scip, newprobdown, 1.0));
563  	   assert(upscore != NULL);
564  	   assert(downscore != NULL);
565  	
566  	   /* update up and downscore depending on score parameter */
567  	   switch( scoreparam )
568  	   {
569  	   case 'l' :
570  	      /* 'l'owest cumulative probability */
571  	      if( SCIPisGT(scip, 1.0 - newprobup, *upscore) )
572  	         *upscore = 1.0 - newprobup;
573  	      if( SCIPisGT(scip, 1.0 - newprobdown, *downscore) )
574  	         *downscore = 1.0 - newprobdown;
575  	      break;
576  	
577  	   case 'd' :
578  	      /* biggest 'd'ifference currentprob - newprob */
579  	      if( SCIPisGT(scip, currentprob - newprobup, *upscore) )
580  	         *upscore = currentprob - newprobup;
581  	      if( SCIPisGT(scip, currentprob - newprobdown, *downscore) )
582  	         *downscore = currentprob - newprobdown;
583  	      break;
584  	
585  	   case 'h' :
586  	      /* 'h'ighest cumulative probability */
587  	      if( SCIPisGT(scip, newprobup, *upscore) )
588  	         *upscore = newprobup;
589  	      if( SCIPisGT(scip, newprobdown, *downscore) )
590  	         *downscore = newprobdown;
591  	      break;
592  	
593  	   case 'v' :
594  	      /* 'v'otes lowest cumulative probability */
595  	      if( SCIPisLT(scip, newprobup, newprobdown) )
596  	         *upscore += 1.0;
597  	      else if( SCIPisGT(scip, newprobup, newprobdown) )
598  	         *downscore += 1.0;
599  	      break;
600  	
601  	   case 'w' :
602  	      /* votes highest cumulative probability */
603  	      if( SCIPisGT(scip, newprobup, newprobdown) )
604  	         *upscore += 1.0;
605  	      else if( SCIPisLT(scip, newprobup, newprobdown) )
606  	         *downscore += 1.0;
607  	      break;
608  	
609  	   default :
610  	      SCIPerrorMessage(" ERROR! No branching scheme selected! Exiting  method.\n");
611  	      return SCIP_INVALIDCALL;
612  	   }
613  	
614  	   return SCIP_OKAY;
615  	}
616  	
617  	/** calculate the branching score of a variable, depending on the chosen score parameter */
618  	static
619  	SCIP_RETCODE calcBranchScore(
620  	   SCIP*                 scip,               /**< current SCIP */
621  	   SCIP_BRANCHRULEDATA*  branchruledata,     /**< branch rule data */
622  	   SCIP_VAR*             var,                /**< candidate variable */
623  	   SCIP_Real             lpsolval,           /**< current fractional LP-relaxation solution value  */
624  	   SCIP_Real*            upscore,            /**< pointer to store the variable score when branching on it in upward direction */
625  	   SCIP_Real*            downscore,          /**< pointer to store the variable score when branching on it in downward direction */
626  	   char                  scoreparam          /**< the score parameter of this branching rule */
627  	   )
628  	{
629  	   SCIP_COL* varcol;
630  	   SCIP_ROW** colrows;
631  	   SCIP_Real* rowvals;
632  	   SCIP_Real varlb;
633  	   SCIP_Real varub;
634  	   SCIP_Real squaredbounddiff; /* current squared difference of variable bounds (ub - lb)^2 */
635  	   SCIP_Real newub;            /* new upper bound if branching downwards */
636  	   SCIP_Real newlb;            /* new lower bound if branching upwards */
637  	   SCIP_Real squaredbounddiffup; /* squared difference after branching upwards (ub - lb')^2 */
638  	   SCIP_Real squaredbounddiffdown; /* squared difference after branching downwards (ub' - lb)^2 */
639  	   SCIP_Real currentmean;      /* current mean value of variable uniform distribution */
640  	   SCIP_Real meanup;           /* mean value of variable uniform distribution after branching up */
641  	   SCIP_Real meandown;         /* mean value of variable uniform distribution after branching down*/
642  	   SCIP_VARTYPE vartype;
643  	   int ncolrows;
644  	   int i;
645  	
646  	   SCIP_Bool onlyactiverows; /* should only rows which are active at the current node be considered? */
647  	
648  	   assert(scip != NULL);
649  	   assert(var != NULL);
650  	   assert(upscore != NULL);
651  	   assert(downscore != NULL);
652  	   assert(!SCIPisIntegral(scip, lpsolval));
653  	   assert(SCIPvarGetStatus(var) == SCIP_VARSTATUS_COLUMN);
654  	
655  	   varcol = SCIPvarGetCol(var);
656  	   assert(varcol != NULL);
657  	
658  	   colrows = SCIPcolGetRows(varcol);
659  	   rowvals = SCIPcolGetVals(varcol);
660  	   ncolrows = SCIPcolGetNNonz(varcol);
661  	   varlb = SCIPvarGetLbLocal(var);
662  	   varub = SCIPvarGetUbLocal(var);
663  	   assert(SCIPisFeasLT(scip, varlb, varub));
664  	   vartype = SCIPvarGetType(var);
665  	
666  	   /* calculate mean and variance of variable uniform distribution before and after branching */
667  	   currentmean = 0.0;
668  	   squaredbounddiff = 0.0;
669  	   SCIPvarCalcDistributionParameters(scip, varlb, varub, vartype, &currentmean, &squaredbounddiff);
670  	
671  	   newlb = SCIPfeasCeil(scip, lpsolval);
672  	   newub = SCIPfeasFloor(scip, lpsolval);
673  	
674  	   /* calculate the variable's uniform distribution after branching up and down, respectively. */
675  	   squaredbounddiffup = 0.0;
676  	   meanup = 0.0;
677  	   SCIPvarCalcDistributionParameters(scip, newlb, varub, vartype, &meanup, &squaredbounddiffup);
678  	
679  	   /* calculate the distribution mean and variance for a variable with finite lower bound */
680  	   squaredbounddiffdown = 0.0;
681  	   meandown = 0.0;
682  	   SCIPvarCalcDistributionParameters(scip, varlb, newub, vartype, &meandown, &squaredbounddiffdown);
683  	
684  	   /* initialize the variable's up and down score */
685  	   *upscore = 0.0;
686  	   *downscore = 0.0;
687  	
688  	   onlyactiverows = branchruledata->onlyactiverows;
689  	
690  	   /* loop over the variable rows and calculate the up and down score */
691  	   for( i = 0; i < ncolrows; ++i )
692  	   {
693  	      SCIP_ROW* row;
694  	      SCIP_Real changedrowmean;
695  	      SCIP_Real rowmean;
696  	      SCIP_Real rowvariance;
697  	      SCIP_Real changedrowvariance;
698  	      SCIP_Real currentrowprob;
699  	      SCIP_Real newrowprobup;
700  	      SCIP_Real newrowprobdown;
701  	      SCIP_Real squaredcoeff;
702  	      SCIP_Real rowval;
703  	      int rowinfinitiesdown;
704  	      int rowinfinitiesup;
705  	      int rowpos;
706  	
707  	      row = colrows[i];
708  	      rowval = rowvals[i];
709  	      assert(row != NULL);
710  	
711  	      /* we access the rows by their index */
712  	      rowpos = SCIProwGetIndex(row);
713  	
714  	      /* skip non-active rows if the user parameter was set this way */
715  	      if( onlyactiverows && SCIPisSumPositive(scip, SCIPgetRowLPFeasibility(scip, row)) )
716  	         continue;
717  	
718  	      /* call method to ensure sufficient data capacity */
719  	      SCIP_CALL( branchruledataEnsureArraySize(scip, branchruledata, rowpos) );
720  	
721  	      /* calculate row activity distribution if this is the first candidate to appear in this row */
722  	      if( branchruledata->rowmeans[rowpos] == SCIP_INVALID ) /*lint !e777*/
723  	      {
724  	         rowCalculateGauss(scip, branchruledata, row, &branchruledata->rowmeans[rowpos], &branchruledata->rowvariances[rowpos],
725  	               &branchruledata->rowinfinitiesdown[rowpos], &branchruledata->rowinfinitiesup[rowpos]);
726  	      }
727  	
728  	      /* retrieve the row distribution parameters from the branch rule data */
729  	      rowmean = branchruledata->rowmeans[rowpos];
730  	      rowvariance = branchruledata->rowvariances[rowpos];
731  	      rowinfinitiesdown = branchruledata->rowinfinitiesdown[rowpos];
732  	      rowinfinitiesup = branchruledata->rowinfinitiesdown[rowpos];
733  	      assert(!SCIPisNegative(scip, rowvariance));
734  	
735  	      currentrowprob = SCIProwCalcProbability(scip, row, rowmean, rowvariance,
736  	            rowinfinitiesdown, rowinfinitiesup);
737  	
738  	      /* get variable's current expected contribution to row activity */
739  	      squaredcoeff = SQUARED(rowval);
740  	
741  	      /* first, get the probability change for the row if the variable is branched on upwards. The probability
742  	       * can only be affected if the variable upper bound is finite
743  	       */
744  	      if( !SCIPisInfinity(scip, varub) )
745  	      {
746  	         int rowinftiesdownafterbranch;
747  	         int rowinftiesupafterbranch;
748  	
749  	         /* calculate how branching would affect the row parameters */
750  	         changedrowmean = rowmean + rowval * (meanup - currentmean);
751  	         changedrowvariance = rowvariance + squaredcoeff * (squaredbounddiffup - squaredbounddiff);
752  	         changedrowvariance = MAX(0.0, changedrowvariance);
753  	
754  	         rowinftiesdownafterbranch = rowinfinitiesdown;
755  	         rowinftiesupafterbranch = rowinfinitiesup;
756  	
757  	         /* account for changes of the row's infinite bound contributions */
758  	         if( SCIPisInfinity(scip, -varlb) && rowval < 0.0 )
759  	            rowinftiesupafterbranch--;
760  	         if( SCIPisInfinity(scip, -varlb) && rowval > 0.0 )
761  	            rowinftiesdownafterbranch--;
762  	
763  	         assert(rowinftiesupafterbranch >= 0);
764  	         assert(rowinftiesdownafterbranch >= 0);
765  	         newrowprobup = SCIProwCalcProbability(scip, row, changedrowmean, changedrowvariance, rowinftiesdownafterbranch,
766  	               rowinftiesupafterbranch);
767  	      }
768  	      else
769  	         newrowprobup = currentrowprob;
770  	
771  	      /* do the same for the other branching direction */
772  	      if( !SCIPisInfinity(scip, varlb) )
773  	      {
774  	         int rowinftiesdownafterbranch;
775  	         int rowinftiesupafterbranch;
776  	
777  	         changedrowmean = rowmean + rowval * (meandown - currentmean);
778  	         changedrowvariance = rowvariance + squaredcoeff * (squaredbounddiffdown - squaredbounddiff);
779  	         changedrowvariance = MAX(0.0, changedrowvariance);
780  	
781  	         rowinftiesdownafterbranch = rowinfinitiesdown;
782  	         rowinftiesupafterbranch = rowinfinitiesup;
783  	
784  	         /* account for changes of the row's infinite bound contributions */
785  	         if( SCIPisInfinity(scip, varub) && rowval > 0.0 )
786  	            rowinftiesupafterbranch -= 1;
787  	         if( SCIPisInfinity(scip, varub) && rowval < 0.0 )
788  	            rowinftiesdownafterbranch -= 1;
789  	
790  	         assert(rowinftiesdownafterbranch >= 0);
791  	         assert(rowinftiesupafterbranch >= 0);
792  	         newrowprobdown = SCIProwCalcProbability(scip, row, changedrowmean, changedrowvariance, rowinftiesdownafterbranch,
793  	               rowinftiesupafterbranch);
794  	      }
795  	      else
796  	         newrowprobdown = currentrowprob;
797  	
798  	      /* update the up and down score depending on the chosen scoring parameter */
799  	      SCIP_CALL( SCIPupdateDistributionScore(scip, currentrowprob, newrowprobup, newrowprobdown, upscore, downscore, scoreparam) );
800  	
801  	      SCIPdebugMsg(scip, "  Variable %s changes probability of row %s from %g to %g (branch up) or %g;\n",
802  	         SCIPvarGetName(var), SCIProwGetName(row), currentrowprob, newrowprobup, newrowprobdown);
803  	      SCIPdebugMsg(scip, "  -->  new variable score: %g (for branching up), %g (for branching down)\n",
804  	         *upscore, *downscore);
805  	   }
806  	
807  	   return SCIP_OKAY;
808  	}
809  	
810  	/** free branchrule data */
811  	static
812  	void branchruledataFreeArrays(
813  	   SCIP*                 scip,               /**< SCIP data structure */
814  	   SCIP_BRANCHRULEDATA*  branchruledata      /**< branching rule data */
815  	   )
816  	{
817  	   assert(branchruledata->memsize == 0 || branchruledata->rowmeans != NULL);
818  	   assert(branchruledata->memsize >= 0);
819  	
820  	   if( branchruledata->memsize > 0 )
821  	   {
822  	      SCIPfreeBlockMemoryArray(scip, &branchruledata->rowmeans, branchruledata->memsize);
823  	      SCIPfreeBlockMemoryArray(scip, &branchruledata->rowvariances, branchruledata->memsize);
824  	      SCIPfreeBlockMemoryArray(scip, &branchruledata->rowinfinitiesup, branchruledata->memsize);
825  	      SCIPfreeBlockMemoryArray(scip, &branchruledata->rowinfinitiesdown, branchruledata->memsize);
826  	
827  	      SCIPfreeBlockMemoryArray(scip, &branchruledata->varfilterposs, branchruledata->varpossmemsize);
828  	      SCIPfreeBlockMemoryArray(scip, &branchruledata->varposs, branchruledata->varpossmemsize);
829  	      SCIPfreeBlockMemoryArray(scip, &branchruledata->updatedvars, branchruledata->varpossmemsize);
830  	      SCIPfreeBlockMemoryArray(scip, &branchruledata->currentubs, branchruledata->varpossmemsize);
831  	      SCIPfreeBlockMemoryArray(scip, &branchruledata->currentlbs, branchruledata->varpossmemsize);
832  	
833  	      branchruledata->memsize = 0;
834  	   }
835  	}
836  	
837  	/** add variable to the bound change event queue; skipped if variable is already in there, or if variable has
838  	 *  no row currently watched
839  	 */
840  	static
841  	void branchruledataAddBoundChangeVar(
842  	   SCIP_BRANCHRULEDATA*  branchruledata,     /**< branchrule data */
843  	   SCIP_VAR*             var                 /**< the variable whose bound changes need to be processed */
844  	   )
845  	{
846  	   int varindex;
847  	   int varpos;
848  	
849  	   assert(var != NULL);
850  	
851  	   varindex = SCIPvarGetProbindex(var);
852  	   assert(-1 <= varindex && varindex < branchruledata->varpossmemsize);
853  	
854  	   /* if variable is not active, it should not be watched */
855  	   if( varindex == -1 )
856  	      return;
857  	   varpos = branchruledata->varposs[varindex];
858  	   assert(varpos < branchruledata->nupdatedvars);
859  	
860  	   /* nothing to do if variable is already in the queue */
861  	   if( varpos >= 0 )
862  	   {
863  	      assert(branchruledata->updatedvars[varpos] == var);
864  	
865  	      return;
866  	   }
867  	
868  	   /* if none of the variables rows was calculated yet, variable needs not to be watched */
869  	   assert((branchruledata->currentlbs[varindex] == SCIP_INVALID) == (branchruledata->currentubs[varindex] == SCIP_INVALID)); /*lint !e777*/
870  	   if( branchruledata->currentlbs[varindex] == SCIP_INVALID ) /*lint !e777*/
871  	      return;
872  	
873  	   /* add the variable to the branch rule data of variables to process updates for */
874  	   assert(branchruledata->varpossmemsize > branchruledata->nupdatedvars);
875  	   varpos = branchruledata->nupdatedvars;
876  	   branchruledata->updatedvars[varpos] = var;
877  	   branchruledata->varposs[varindex] = varpos;
878  	   ++branchruledata->nupdatedvars;
879  	}
880  	
881  	/** returns the next unprocessed variable (last in, first out) with pending bound changes, or NULL */
882  	static
883  	SCIP_VAR* branchruledataPopBoundChangeVar(
884  	   SCIP_BRANCHRULEDATA*  branchruledata      /**< branchrule data */
885  	   )
886  	{
887  	   SCIP_VAR* var;
888  	   int varpos;
889  	   int varindex;
890  	
891  	   assert(branchruledata->nupdatedvars >= 0);
892  	
893  	   /* return if no variable is currently pending */
894  	   if( branchruledata->nupdatedvars == 0 )
895  	      return NULL;
896  	
897  	   varpos = branchruledata->nupdatedvars - 1;
898  	   var = branchruledata->updatedvars[varpos];
899  	   assert(var != NULL);
900  	   varindex = SCIPvarGetProbindex(var);
901  	   assert(0 <= varindex && varindex < branchruledata->varpossmemsize);
902  	   assert(varpos == branchruledata->varposs[varindex]);
903  	
904  	   branchruledata->varposs[varindex] = -1;
905  	   branchruledata->nupdatedvars--;
906  	
907  	   return var;
908  	}
909  	
910  	/** process a variable from the queue of changed variables */
911  	static
912  	SCIP_RETCODE varProcessBoundChanges(
913  	   SCIP*                 scip,               /**< SCIP data structure */
914  	   SCIP_BRANCHRULEDATA*  branchruledata,     /**< branchrule data */
915  	   SCIP_VAR*             var                 /**< the variable whose bound changes need to be processed */
916  	   )
917  	{
918  	   SCIP_ROW** colrows;
919  	   SCIP_COL* varcol;
920  	   SCIP_Real* colvals;
921  	   SCIP_Real oldmean;
922  	   SCIP_Real newmean;
923  	   SCIP_Real oldvariance;
924  	   SCIP_Real newvariance;
925  	   SCIP_Real oldlb;
926  	   SCIP_Real newlb;
927  	   SCIP_Real oldub;
928  	   SCIP_Real newub;
929  	   SCIP_VARTYPE vartype;
930  	   int ncolrows;
931  	   int r;
932  	   int varindex;
933  	
934  	   /* skip event execution if SCIP is in Probing mode because these bound changes will be undone anyway before branching
935  	    * rule is called again
936  	    */
937  	   assert(!SCIPinProbing(scip));
938  	
939  	   assert(var != NULL);
940  	   varcol = SCIPvarGetCol(var);
941  	   assert(varcol != NULL);
942  	   colrows = SCIPcolGetRows(varcol);
943  	   colvals = SCIPcolGetVals(varcol);
944  	   ncolrows = SCIPcolGetNNonz(varcol);
945  	
946  	   varindex = SCIPvarGetProbindex(var);
947  	
948  	   oldlb = branchruledata->currentlbs[varindex];
949  	   oldub = branchruledata->currentubs[varindex];
950  	
951  	   /* skip update if the variable has never been subject of previously calculated row activities */
952  	   assert((oldlb == SCIP_INVALID) == (oldub == SCIP_INVALID)); /*lint !e777*/
953  	   if( oldlb == SCIP_INVALID ) /*lint !e777*/
954  	      return SCIP_OKAY;
955  	
956  	   newlb = SCIPvarGetLbLocal(var);
957  	   newub = SCIPvarGetUbLocal(var);
958  	
959  	   /* skip update if the bound change events have cancelled out */
960  	   if( SCIPisFeasEQ(scip, oldlb, newlb) && SCIPisFeasEQ(scip, oldub, newub) )
961  	      return SCIP_OKAY;
962  	
963  	   /* calculate old and new variable distribution mean and variance */
964  	   oldvariance = 0.0;
965  	   newvariance = 0.0;
966  	   oldmean = 0.0;
967  	   newmean = 0.0;
968  	   vartype = SCIPvarGetType(var);
969  	   SCIPvarCalcDistributionParameters(scip, oldlb, oldub, vartype, &oldmean, &oldvariance);
970  	   SCIPvarCalcDistributionParameters(scip, newlb, newub, vartype, &newmean, &newvariance);
971  	
972  	   /* loop over all rows of this variable and update activity distribution */
973  	   for( r = 0; r < ncolrows; ++r )
974  	   {
975  	      int rowpos;
976  	
977  	      assert(colrows[r] != NULL);
978  	      rowpos = SCIProwGetIndex(colrows[r]);
979  	      assert(rowpos >= 0);
980  	
981  	      SCIP_CALL( branchruledataEnsureArraySize(scip, branchruledata, rowpos) );
982  	
983  	      /* only consider rows for which activity distribution was already calculated */
984  	      if( branchruledata->rowmeans[rowpos] != SCIP_INVALID ) /*lint !e777*/
985  	      {
986  	         SCIP_Real coeff;
987  	         SCIP_Real coeffsquared;
988  	         /*lint -e777*/
989  	         assert(branchruledata->rowvariances[rowpos] != SCIP_INVALID
990  	               && SCIPisFeasGE(scip, branchruledata->rowvariances[rowpos], 0.0));
991  	
992  	         coeff = colvals[r];
993  	         coeffsquared = SQUARED(coeff);
994  	
995  	         /* update variable contribution to row activity distribution */
996  	         branchruledata->rowmeans[rowpos] += coeff * (newmean - oldmean);
997  	         branchruledata->rowvariances[rowpos] += coeffsquared * (newvariance - oldvariance);
998  	         branchruledata->rowvariances[rowpos] = MAX(0.0, branchruledata->rowvariances[rowpos]);
999  	
1000 	         /* account for changes of the infinite contributions to row activities */
1001 	         if( coeff > 0.0 )
1002 	         {
1003 	            /* if the coefficient is positive, upper bounds affect activity up */
1004 	            if( SCIPisInfinity(scip, newub) && !SCIPisInfinity(scip, oldub) )
1005 	               ++branchruledata->rowinfinitiesup[rowpos];
1006 	            else if( !SCIPisInfinity(scip, newub) && SCIPisInfinity(scip, oldub) )
1007 	               --branchruledata->rowinfinitiesup[rowpos];
1008 	
1009 	            if( SCIPisInfinity(scip, newlb) && !SCIPisInfinity(scip, oldlb) )
1010 	               ++branchruledata->rowinfinitiesdown[rowpos];
1011 	            else if( !SCIPisInfinity(scip, newlb) && SCIPisInfinity(scip, oldlb) )
1012 	               --branchruledata->rowinfinitiesdown[rowpos];
1013 	         }
1014 	         else if( coeff < 0.0 )
1015 	         {
1016 	            if( SCIPisInfinity(scip, newub) && !SCIPisInfinity(scip, oldub) )
1017 	               ++branchruledata->rowinfinitiesdown[rowpos];
1018 	            else if( !SCIPisInfinity(scip, newub) && SCIPisInfinity(scip, oldub) )
1019 	               --branchruledata->rowinfinitiesdown[rowpos];
1020 	
1021 	            if( SCIPisInfinity(scip, newlb) && !SCIPisInfinity(scip, oldlb) )
1022 	               ++branchruledata->rowinfinitiesup[rowpos];
1023 	            else if( !SCIPisInfinity(scip, newlb) && SCIPisInfinity(scip, oldlb) )
1024 	               --branchruledata->rowinfinitiesup[rowpos];
1025 	         }
1026 	         assert(branchruledata->rowinfinitiesdown[rowpos] >= 0);
1027 	         assert(branchruledata->rowinfinitiesup[rowpos] >= 0);
1028 	      }
1029 	   }
1030 	
1031 	   /* store the new local bounds in the data */
1032 	   branchruledataUpdateCurrentBounds(scip, branchruledata, var);
1033 	
1034 	   return SCIP_OKAY;
1035 	}
1036 	
1037 	/** destructor of event handler to free user data (called when SCIP is exiting) */
1038 	static
1039 	SCIP_DECL_EVENTFREE(eventFreeDistribution)
1040 	{
1041 	   SCIP_EVENTHDLRDATA* eventhdlrdata;
1042 	
1043 	   eventhdlrdata = SCIPeventhdlrGetData(eventhdlr);
1044 	   assert(eventhdlrdata != NULL);
1045 	
1046 	   SCIPfreeBlockMemory(scip, &eventhdlrdata);
1047 	   SCIPeventhdlrSetData(eventhdlr, NULL);
1048 	
1049 	   return SCIP_OKAY;
1050 	}
1051 	
1052 	/*
1053 	 * Callback methods of branching rule
1054 	 */
1055 	
1056 	/** copy method for branchrule plugins (called when SCIP copies plugins) */
1057 	static
1058 	SCIP_DECL_BRANCHCOPY(branchCopyDistribution)
1059 	{  /*lint --e{715}*/
1060 	   assert(scip != NULL);
1061 	
1062 	   SCIP_CALL( SCIPincludeBranchruleDistribution(scip) );
1063 	
1064 	   return SCIP_OKAY;
1065 	}
1066 	
1067 	/** solving process deinitialization method of branching rule (called before branch and bound process data is freed) */
1068 	static
1069 	SCIP_DECL_BRANCHEXITSOL(branchExitsolDistribution)
1070 	{
1071 	   SCIP_BRANCHRULEDATA* branchruledata;
1072 	
1073 	   assert(branchrule != NULL);
1074 	   assert(strcmp(SCIPbranchruleGetName(branchrule), BRANCHRULE_NAME) == 0);
1075 	
1076 	   branchruledata = SCIPbranchruleGetData(branchrule);
1077 	   assert(branchruledata != NULL);
1078 	
1079 	   /* drop variable events at the end of branch and bound process (cannot be used after restarts, anyway) */
1080 	   if( branchruledata->varfilterposs != NULL)
1081 	   {
1082 	      SCIP_VAR** vars;
1083 	      int nvars;
1084 	      int v;
1085 	
1086 	      vars = SCIPgetVars(scip);
1087 	      nvars = SCIPgetNVars(scip);
1088 	
1089 	      assert(nvars > 0);
1090 	      for( v = 0; v < nvars; ++v )
1091 	      {
1092 	         SCIP_CALL( SCIPdropVarEvent(scip, vars[v], EVENT_DISTRIBUTION, branchruledata->eventhdlr, NULL, branchruledata->varfilterposs[v]) );
1093 	      }
1094 	   }
1095 	
1096 	   /* free row arrays when branch-and-bound data is freed */
1097 	   branchruledataFreeArrays(scip, branchruledata);
1098 	
1099 	   return SCIP_OKAY;
1100 	}
1101 	
1102 	/** destructor of branching rule to free user data (called when SCIP is exiting) */
1103 	static
1104 	SCIP_DECL_BRANCHFREE(branchFreeDistribution)
1105 	{
1106 	   SCIP_BRANCHRULEDATA* branchruledata;
1107 	
1108 	   assert(branchrule != NULL);
1109 	   assert(strcmp(SCIPbranchruleGetName(branchrule), BRANCHRULE_NAME) == 0);
1110 	
1111 	   branchruledata = SCIPbranchruleGetData(branchrule);
1112 	   assert(branchruledata != NULL);
1113 	
1114 	   /* free internal arrays first */
1115 	   branchruledataFreeArrays(scip, branchruledata);
1116 	   SCIPfreeBlockMemory(scip, &branchruledata);
1117 	   SCIPbranchruleSetData(branchrule, NULL);
1118 	
1119 	   return SCIP_OKAY;
1120 	}
1121 	
1122 	/** branching execution method for fractional LP solutions */
1123 	static
1124 	SCIP_DECL_BRANCHEXECLP(branchExeclpDistribution)
1125 	{  /*lint --e{715}*/
1126 	   SCIP_BRANCHRULEDATA* branchruledata;
1127 	   SCIP_VAR** lpcands;
1128 	   SCIP_VAR* bestcand;
1129 	   SCIP_NODE* downchild;
1130 	   SCIP_NODE* upchild;
1131 	
1132 	   SCIP_Real* lpcandssol;
1133 	
1134 	   SCIP_Real bestscore;
1135 	   SCIP_BRANCHDIR bestbranchdir;
1136 	   int nlpcands;
1137 	   int c;
1138 	
1139 	   assert(branchrule != NULL);
1140 	   assert(strcmp(SCIPbranchruleGetName(branchrule), BRANCHRULE_NAME) == 0);
1141 	   assert(scip != NULL);
1142 	   assert(result != NULL);
1143 	
1144 	   *result = SCIP_DIDNOTRUN;
1145 	
1146 	   SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, NULL, NULL, &nlpcands, NULL) );
1147 	
1148 	   if( nlpcands == 0 )
1149 	      return SCIP_OKAY;
1150 	
1151 	   if( SCIPgetNActivePricers(scip) > 0 )
1152 	      return SCIP_OKAY;
1153 	
1154 	   branchruledata = SCIPbranchruleGetData(branchrule);
1155 	
1156 	   /* if branching rule data arrays were not initialized before (usually the first call of the branching execution),
1157 	    * allocate arrays with an initial capacity of the number of LP rows */
1158 	   if( branchruledata->memsize == 0 )
1159 	   {
1160 	      int nlprows;
1161 	
1162 	      /* get LP rows data */
1163 	      nlprows = SCIPgetNLPRows(scip);
1164 	
1165 	      /* initialize arrays only if there are actual LP rows to operate on */
1166 	      if( nlprows > 0 )
1167 	      {
1168 	         SCIP_CALL( branchruledataEnsureArraySize(scip, branchruledata, nlprows) );
1169 	      }
1170 	      else /* if there are no LP rows, branching rule cannot be used */
1171 	         return SCIP_OKAY;
1172 	   }
1173 	
1174 	   /* process pending bound change events */
1175 	   while( branchruledata->nupdatedvars > 0 )
1176 	   {
1177 	      SCIP_VAR* nextvar;
1178 	
1179 	      /* pop the next variable from the queue and process its bound changes */
1180 	      nextvar = branchruledataPopBoundChangeVar(branchruledata);
1181 	      assert(nextvar != NULL);
1182 	      SCIP_CALL( varProcessBoundChanges(scip, branchruledata, nextvar) );
1183 	   }
1184 	
1185 	   bestscore = -1;
1186 	   bestbranchdir = SCIP_BRANCHDIR_AUTO;
1187 	   bestcand = NULL;
1188 	
1189 	   /* invalidate the current row distribution data which is initialized on the fly when looping over the candidates */
1190 	
1191 	   /* loop over candidate variables and calculate their score in changing the cumulative
1192 	    * probability of fulfilling each of their constraints */
1193 	   for( c = 0; c < nlpcands; ++c )
1194 	   {
1195 	      SCIP_Real upscore;
1196 	      SCIP_Real downscore;
1197 	      SCIP_VAR* lpcand;
1198 	      int varindex;
1199 	
1200 	      lpcand = lpcands[c];
1201 	      assert(lpcand != NULL);
1202 	
1203 	      varindex = SCIPvarGetProbindex(lpcand);
1204 	
1205 	      /* in debug mode, ensure that all bound process events which occurred in the mean time have been captured
1206 	       * by the branching rule event system
1207 	       */
1208 	      assert(SCIPisFeasLE(scip, SCIPvarGetLbLocal(lpcand), SCIPvarGetUbLocal(lpcand)));
1209 	      assert(0 <= varindex && varindex < branchruledata->varpossmemsize);
1210 	
1211 	      /*lint -e777*/
1212 	      assert((branchruledata->currentlbs[varindex] == SCIP_INVALID) == (branchruledata->currentubs[varindex] == SCIP_INVALID));
1213 	      assert((branchruledata->currentlbs[varindex] == SCIP_INVALID)
1214 	            || SCIPisFeasEQ(scip, SCIPvarGetLbLocal(lpcand), branchruledata->currentlbs[varindex]));
1215 	      assert((branchruledata->currentubs[varindex] == SCIP_INVALID)
1216 	                  || SCIPisFeasEQ(scip, SCIPvarGetUbLocal(lpcand), branchruledata->currentubs[varindex]));
1217 	
1218 	      /* if the branching rule has not captured the variable bounds yet, this can be done now */
1219 	      if( branchruledata->currentlbs[varindex] == SCIP_INVALID ) /*lint !e777*/
1220 	      {
1221 	         branchruledataUpdateCurrentBounds(scip, branchruledata, lpcand);
1222 	      }
1223 	
1224 	      upscore = 0.0;
1225 	      downscore = 0.0;
1226 	
1227 	      /* loop over candidate rows and determine the candidate up- and down- branching score w.r.t. the score parameter */
1228 	      SCIP_CALL( calcBranchScore(scip, branchruledata, lpcand, lpcandssol[c],
1229 	            &upscore, &downscore, branchruledata->scoreparam) );
1230 	
1231 	      /* if weighted scoring is enabled, use the branching score method of SCIP to weigh up and down score */
1232 	      if( branchruledata->usescipscore )
1233 	      {
1234 	         SCIP_Real score;
1235 	
1236 	         score = SCIPgetBranchScore(scip, lpcand, downscore, upscore);
1237 	
1238 	         /* select the candidate with the highest branching score */
1239 	         if( score > bestscore )
1240 	         {
1241 	            bestscore = score;
1242 	            bestcand = lpcand;
1243 	            /* prioritize branching direction with the higher score */
1244 	            if( upscore > downscore )
1245 	               bestbranchdir = SCIP_BRANCHDIR_UPWARDS;
1246 	            else
1247 	               bestbranchdir = SCIP_BRANCHDIR_DOWNWARDS;
1248 	         }
1249 	      }
1250 	      else
1251 	      {
1252 	         /* no weighted score; keep candidate which has the single highest score in one direction */
1253 	         if( upscore > bestscore && upscore > downscore )
1254 	         {
1255 	            bestscore = upscore;
1256 	            bestbranchdir = SCIP_BRANCHDIR_UPWARDS;
1257 	            bestcand = lpcand;
1258 	         }
1259 	         else if( downscore > bestscore )
1260 	         {
1261 	            bestscore = downscore;
1262 	            bestbranchdir = SCIP_BRANCHDIR_DOWNWARDS;
1263 	            bestcand = lpcand;
1264 	         }
1265 	      }
1266 	
1267 	      SCIPdebugMsg(scip, "  Candidate %s has score down %g and up %g \n", SCIPvarGetName(lpcand), downscore, upscore);
1268 	      SCIPdebugMsg(scip, "  Best candidate: %s, score %g, direction %d\n", SCIPvarGetName(bestcand), bestscore, bestbranchdir);
1269 	   }
1270 	   assert(!SCIPisFeasIntegral(scip, SCIPvarGetSol(bestcand, TRUE)));
1271 	   assert(bestbranchdir == SCIP_BRANCHDIR_DOWNWARDS || bestbranchdir == SCIP_BRANCHDIR_UPWARDS);
1272 	   assert(bestcand != NULL);
1273 	
1274 	   SCIPdebugMsg(scip, "  Branching on variable %s with bounds [%g, %g] and solution value <%g>\n", SCIPvarGetName(bestcand),
1275 	      SCIPvarGetLbLocal(bestcand), SCIPvarGetUbLocal(bestcand), SCIPvarGetLPSol(bestcand));
1276 	
1277 	   /* branch on the best candidate variable */
1278 	   SCIP_CALL( SCIPbranchVar(scip, bestcand, &downchild, NULL, &upchild) );
1279 	
1280 	   assert(downchild != NULL);
1281 	   assert(upchild != NULL);
1282 	
1283 	   if( bestbranchdir == SCIP_BRANCHDIR_UPWARDS )
1284 	   {
1285 	      SCIP_CALL( SCIPchgChildPrio(scip, upchild, DEFAULT_PRIORITY) );
1286 	      SCIPdebugMsg(scip, "  Changing node priority of up-child.\n");
1287 	   }
1288 	   else
1289 	   {
1290 	      assert(bestbranchdir == SCIP_BRANCHDIR_DOWNWARDS);
1291 	      SCIP_CALL( SCIPchgChildPrio(scip, downchild, DEFAULT_PRIORITY) );
1292 	      SCIPdebugMsg(scip, "  Changing node priority of down-child.\n");
1293 	   }
1294 	
1295 	   *result = SCIP_BRANCHED;
1296 	
1297 	   return SCIP_OKAY;
1298 	}
1299 	
1300 	/** event execution method of distribution branching which handles bound change events of variables */
1301 	static
1302 	SCIP_DECL_EVENTEXEC(eventExecDistribution)
1303 	{  /*lint --e{715}*/
1304 	   SCIP_BRANCHRULEDATA* branchruledata;
1305 	   SCIP_EVENTHDLRDATA* eventhdlrdata;
1306 	   SCIP_VAR* var;
1307 	
1308 	   assert(eventhdlr != NULL);
1309 	   eventhdlrdata = SCIPeventhdlrGetData(eventhdlr);
1310 	   assert(eventhdlrdata != NULL);
1311 	
1312 	   branchruledata = eventhdlrdata->branchruledata;
1313 	   var = SCIPeventGetVar(event);
1314 	
1315 	   /* add the variable to the queue of unprocessed variables; method itself ensures that every variable is added at most once */
1316 	   branchruledataAddBoundChangeVar(branchruledata, var);
1317 	
1318 	   return SCIP_OKAY;
1319 	}
1320 	
1321 	
1322 	/*
1323 	 * branching rule specific interface methods
1324 	 */
1325 	
1326 	/** creates the distribution branching rule and includes it in SCIP */
1327 	SCIP_RETCODE SCIPincludeBranchruleDistribution(
1328 	   SCIP*                 scip                /**< SCIP data structure */
1329 	   )
1330 	{
1331 	   SCIP_BRANCHRULE* branchrule = NULL;
1332 	   SCIP_BRANCHRULEDATA* branchruledata;
1333 	   SCIP_EVENTHDLRDATA* eventhdlrdata;
1334 	
1335 	   /* create distribution branching rule data */
1336 	   branchruledata = NULL;
1337 	   SCIP_CALL( SCIPallocBlockMemory(scip, &branchruledata) );
1338 	
1339 	   branchruledata->memsize = 0;
1340 	   branchruledata->rowmeans = NULL;
1341 	   branchruledata->rowvariances = NULL;
1342 	   branchruledata->rowinfinitiesdown = NULL;
1343 	   branchruledata->rowinfinitiesup = NULL;
1344 	   branchruledata->varfilterposs = NULL;
1345 	   branchruledata->currentlbs = NULL;
1346 	   branchruledata->currentubs = NULL;
1347 	
1348 	   /* create event handler first to finish branch rule data */
1349 	   eventhdlrdata = NULL;
1350 	   SCIP_CALL( SCIPallocBlockMemory(scip, &eventhdlrdata) );
1351 	   eventhdlrdata->branchruledata = branchruledata;
1352 	
1353 	   branchruledata->eventhdlr = NULL;
1354 	   SCIP_CALL( SCIPincludeEventhdlrBasic(scip, &branchruledata->eventhdlr, EVENTHDLR_NAME,
1355 	         "event handler for dynamic acitivity distribution updating",
1356 	         eventExecDistribution, eventhdlrdata) );
1357 	   assert( branchruledata->eventhdlr != NULL);
1358 	   SCIP_CALL( SCIPsetEventhdlrFree(scip, branchruledata->eventhdlr, eventFreeDistribution) );
1359 	
1360 	   /* include branching rule */
1361 	   SCIP_CALL( SCIPincludeBranchruleBasic(scip, &branchrule, BRANCHRULE_NAME, BRANCHRULE_DESC, BRANCHRULE_PRIORITY, BRANCHRULE_MAXDEPTH,
1362 		 BRANCHRULE_MAXBOUNDDIST, branchruledata) );
1363 	
1364 	   assert(branchrule != NULL);
1365 	   SCIP_CALL( SCIPsetBranchruleCopy(scip, branchrule, branchCopyDistribution) );
1366 	   SCIP_CALL( SCIPsetBranchruleFree(scip, branchrule, branchFreeDistribution) );
1367 	   SCIP_CALL( SCIPsetBranchruleExitsol(scip, branchrule, branchExitsolDistribution) );
1368 	   SCIP_CALL( SCIPsetBranchruleExecLp(scip, branchrule, branchExeclpDistribution) );
1369 	
1370 	   /* add distribution branching rule parameters */
1371 	   SCIP_CALL( SCIPaddCharParam(scip, "branching/" BRANCHRULE_NAME "/scoreparam",
1372 	         "the score;largest 'd'ifference, 'l'owest cumulative probability,'h'ighest c.p., 'v'otes lowest c.p., votes highest c.p.('w') ",
1373 	         &branchruledata->scoreparam, TRUE, DEFAULT_SCOREPARAM, SCOREPARAM_VALUES, NULL, NULL) );
1374 	
1375 	   SCIP_CALL( SCIPaddBoolParam(scip, "branching/" BRANCHRULE_NAME "/onlyactiverows",
1376 	         "should only rows which are active at the current node be considered?",
1377 	         &branchruledata->onlyactiverows, TRUE, DEFAULT_ONLYACTIVEROWS, NULL, NULL) );
1378 	
1379 	   SCIP_CALL( SCIPaddBoolParam(scip, "branching/" BRANCHRULE_NAME "/weightedscore",
1380 	         "should the branching score weigh up- and down-scores of a variable",
1381 	         &branchruledata->usescipscore, TRUE, DEFAULT_USEWEIGHTEDSCORE, NULL, NULL) );
1382 	
1383 	   return SCIP_OKAY;
1384 	}
1385