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_pscost.c
26   	 * @ingroup DEFPLUGINS_BRANCH
27   	 * @brief  pseudo costs branching rule
28   	 * @author Tobias Achterberg
29   	 * @author Stefan Vigerske
30   	 */
31   	
32   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
33   	
34   	#include "blockmemshell/memory.h"
35   	#include "scip/branch_pscost.h"
36   	#include "scip/pub_branch.h"
37   	#include "scip/pub_message.h"
38   	#include "scip/pub_misc.h"
39   	#include "scip/pub_misc_sort.h"
40   	#include "scip/pub_tree.h"
41   	#include "scip/pub_var.h"
42   	#include "scip/scip_branch.h"
43   	#include "scip/scip_mem.h"
44   	#include "scip/scip_message.h"
45   	#include "scip/scip_numerics.h"
46   	#include "scip/scip_param.h"
47   	#include "scip/scip_randnumgen.h"
48   	#include "scip/scip_sol.h"
49   	#include "scip/scip_tree.h"
50   	#include "scip/scip_var.h"
51   	#include <string.h>
52   	
53   	#define BRANCHRULE_NAME          "pscost"
54   	#define BRANCHRULE_DESC          "branching on pseudo cost values"
55   	#define BRANCHRULE_PRIORITY      2000
56   	#define BRANCHRULE_MAXDEPTH      -1
57   	#define BRANCHRULE_MAXBOUNDDIST  1.0
58   	
59   	#define BRANCHRULE_STRATEGIES          "dsuv" /**< possible pseudo cost multiplication strategies for branching on external candidates */
60   	#define BRANCHRULE_STRATEGY_DEFAULT       'u' /**< default pseudo cost multiplication strategy */
61   	#define BRANCHRULE_SCOREMINWEIGHT_DEFAULT 0.8 /**< default weight for minimum of scores of a branching candidate */
62   	#define BRANCHRULE_SCOREMAXWEIGHT_DEFAULT 1.3 /**< default weight for maximum of scores of a branching candidate */
63   	#define BRANCHRULE_SCORESUMWEIGHT_DEFAULT 0.1 /**< default weight for sum of scores of a branching candidate */
64   	#define BRANCHRULE_NCHILDREN_DEFAULT        2 /**< default number of children in n-ary branching */
65   	#define BRANCHRULE_NARYMAXDEPTH_DEFAULT    -1 /**< default maximal depth where to do n-ary branching */
66   	#define BRANCHRULE_NARYMINWIDTH_DEFAULT 0.001 /**< default minimal domain width in children when doing n-ary branching */
67   	#define BRANCHRULE_NARYWIDTHFAC_DEFAULT   2.0 /**< default factor of domain width in n-ary branching */
68   	#define BRANCHRULE_RANDSEED_DEFAULT        47 /**< initial random seed */
69   	
70   	
71   	#define WEIGHTEDSCORING(data, min, max, sum) \
72   	   ((data)->scoreminweight * (min) + (data)->scoremaxweight * (max) + (data)->scoresumweight * (sum))
73   	
74   	/** branching rule data */
75   	struct SCIP_BranchruleData
76   	{
77   	   SCIP_RANDNUMGEN*      randnumgen;         /**< random number generator */
78   	
79   	   char                  strategy;           /**< strategy for computing score of external candidates */
80   	   SCIP_Real             scoreminweight;     /**< weight for minimum of scores of a branching candidate */
81   	   SCIP_Real             scoremaxweight;     /**< weight for maximum of scores of a branching candidate */
82   	   SCIP_Real             scoresumweight;     /**< weight for sum of scores of a branching candidate */
83   	
84   	   char                  updatestrategy;     /**< strategy used to update pseudo costs of continuous variables */
85   	
86   	   int                   nchildren;          /**< targeted number of children in n-ary branching */
87   	   int                   narymaxdepth;       /**< maximal depth where to do n-ary branching, -1 to turn off */
88   	   SCIP_Real             naryminwidth;       /**< minimal domain width in children when doing n-ary branching, relative to global bounds */
89   	   SCIP_Real             narywidthfactor;    /**< factor of domain width in n-ary branching */
90   	};
91   	
92   	/*
93   	 * Local methods
94   	 */
95   	
96   	/** checks if a given branching candidate is better than a previous one and updates the best branching candidate accordingly */
97   	static
98   	SCIP_RETCODE updateBestCandidate(
99   	   SCIP*                 scip,               /**< SCIP data structure */
100  	   SCIP_BRANCHRULEDATA*  branchruledata,     /**< branching rule data */
101  	   SCIP_VAR**            bestvar,            /**< best branching candidate */
102  	   SCIP_Real*            bestbrpoint,        /**< branching point for best branching candidate */
103  	   SCIP_Real*            bestscore,          /**< score of best branching candidate */
104  	   SCIP_Real*            bestrndscore,       /**< random score of the best branching candidate */
105  	   SCIP_VAR*             cand,               /**< branching candidate to consider */
106  	   SCIP_Real             candscoremin,       /**< minimal score of branching candidate */
107  	   SCIP_Real             candscoremax,       /**< maximal score of branching candidate */
108  	   SCIP_Real             candscoresum,       /**< sum of scores of branching candidate */
109  	   SCIP_Real             candrndscore,       /**< random score of branching candidate */
110  	   SCIP_Real             candsol             /**< proposed branching point of branching candidate */
111  	   )
112  	{
113  	   SCIP_Real candbrpoint;
114  	   SCIP_Real branchscore;
115  	
116  	   SCIP_Real deltaminus;
117  	   SCIP_Real deltaplus;
118  	
119  	   SCIP_Real pscostdown;
120  	   SCIP_Real pscostup;
121  	
122  	   char strategy;
123  	
124  	   assert(scip != NULL);
125  	   assert(branchruledata != NULL);
126  	   assert(bestvar != NULL);
127  	   assert(bestbrpoint != NULL);
128  	   assert(bestscore != NULL);
129  	   assert(cand != NULL);
130  	
131  	   /* a branching variable candidate should either be an active problem variable or a multi-aggregated variable */
132  	   assert(SCIPvarIsActive(SCIPvarGetProbvar(cand)) ||
133  	      SCIPvarGetStatus(SCIPvarGetProbvar(cand)) == SCIP_VARSTATUS_MULTAGGR);
134  	
135  	   if( SCIPvarGetStatus(SCIPvarGetProbvar(cand)) == SCIP_VARSTATUS_MULTAGGR )
136  	   {
137  	      /* for a multi-aggregated variable, we call updateBestCandidate function recursively with all variables in the multi-aggregation */
138  	      SCIP_VAR** multvars;
139  	      int nmultvars;
140  	      int i;
141  	      SCIP_Bool success;
142  	      SCIP_Real multvarlb;
143  	      SCIP_Real multvarub;
144  	
145  	      cand = SCIPvarGetProbvar(cand);
146  	      multvars = SCIPvarGetMultaggrVars(cand);
147  	      nmultvars = SCIPvarGetMultaggrNVars(cand);
148  	
149  	      /* if we have a candidate branching point, then first register only aggregation variables
150  	       * for which we can compute a corresponding branching point too (see also comments below)
151  	       * if this fails, then register all (unfixed) aggregation variables, thereby forgetting about candsol
152  	       */
153  	      success = FALSE;
154  	      if( candsol != SCIP_INVALID ) /*lint !e777*/
155  	      {
156  	         SCIP_Real* multscalars;
157  	         SCIP_Real minact;
158  	         SCIP_Real maxact;
159  	         SCIP_Real aggrvarsol;
160  	         SCIP_Real aggrvarsol1;
161  	         SCIP_Real aggrvarsol2;
162  	
163  	         multscalars = SCIPvarGetMultaggrScalars(cand);
164  	
165  	         /* for computing the branching point, we need the current bounds of the multi-aggregated variable */
166  	         minact = SCIPcomputeVarLbLocal(scip, cand);
167  	         maxact = SCIPcomputeVarUbLocal(scip, cand);
168  	
169  	         for( i = 0; i < nmultvars; ++i )
170  	         {
171  	            /* skip fixed variables */
172  	            multvarlb = SCIPcomputeVarLbLocal(scip, multvars[i]);
173  	            multvarub = SCIPcomputeVarUbLocal(scip, multvars[i]);
174  	            if( SCIPisEQ(scip, multvarlb, multvarub) )
175  	               continue;
176  	
177  	            assert(multscalars != NULL);
178  	            assert(multscalars[i] != 0.0);
179  	
180  	            /* we cannot ensure that both the upper bound in the left node and the lower bound in the right node
181  	             * will be candsol by a clever choice for the branching point of multvars[i],
182  	             * but we can try to ensure that at least one of them will be at candsol
183  	             */
184  	            if( multscalars[i] > 0.0 )
185  	            {
186  	               /*    cand >= candsol
187  	                * if multvars[i] >= (candsol - (maxact - multscalars[i] * ub(multvars[i]))) / multscalars[i]
188  	                *                 = (candsol - maxact) / multscalars[i] + ub(multvars[i])
189  	                */
190  	               aggrvarsol1 = (candsol - maxact) / multscalars[i] + multvarub;
191  	
192  	               /*     cand <= candsol
193  	                * if multvars[i] <= (candsol - (minact - multscalar[i] * lb(multvars[i]))) / multscalars[i]
194  	                *                 = (candsol - minact) / multscalars[i] + lb(multvars[i])
195  	                */
196  	               aggrvarsol2 = (candsol - minact) / multscalars[i] + multvarlb;
197  	            }
198  	            else
199  	            {
200  	               /*    cand >= candsol
201  	                * if multvars[i] <= (candsol - (maxact - multscalars[i] * lb(multvars[i]))) / multscalars[i]
202  	                *                 = (candsol - maxact) / multscalars[i] + lb(multvars[i])
203  	                */
204  	               aggrvarsol2 = (candsol - maxact) / multscalars[i] + multvarlb;
205  	
206  	               /*    cand <= candsol
207  	                * if multvars[i] >= (candsol - (minact - multscalar[i] * ub(multvars[i]))) / multscalars[i]
208  	                *                 = (candsol - minact) / multscalars[i] + ub(multvars[i])
209  	                */
210  	               aggrvarsol1 = (candsol - minact) / multscalars[i] + multvarub;
211  	            }
212  	
213  	            /* by the above choice, aggrvarsol1 <= ub(multvars[i]) and aggrvarsol2 >= lb(multvars[i])
214  	             * if aggrvarsol1 <= lb(multvars[i]) or aggrvarsol2 >= ub(multvars[i]), then choose the other one
215  	             * if both are out of bounds, then give up
216  	             * if both are inside bounds, then choose the one closer to 0.0 (someone has better idea???)
217  	             */
218  	            if( SCIPisFeasLE(scip, aggrvarsol1, multvarlb) )
219  	            {
220  	               if( SCIPisFeasGE(scip, aggrvarsol2, multvarub) )
221  	                  continue;
222  	               else
223  	                  aggrvarsol = aggrvarsol2;
224  	            }
225  	            else
226  	            {
227  	               if( SCIPisFeasGE(scip, aggrvarsol2, multvarub) )
228  	                  aggrvarsol = aggrvarsol1;
229  	               else
230  	                  aggrvarsol = REALABS(aggrvarsol1) < REALABS(aggrvarsol2) ? aggrvarsol1 : aggrvarsol2;
231  	            }
232  	            success = TRUE;
233  	
234  	            SCIP_CALL( updateBestCandidate(scip, branchruledata, bestvar, bestbrpoint, bestscore, bestrndscore,
235  	                  multvars[i], candscoremin, candscoremax, candscoresum, candrndscore, aggrvarsol) );
236  	         }
237  	      }
238  	
239  	      if( !success )
240  	         for( i = 0; i < nmultvars; ++i )
241  	         {
242  	            /* skip fixed variables */
243  	            multvarlb = SCIPcomputeVarLbLocal(scip, multvars[i]);
244  	            multvarub = SCIPcomputeVarUbLocal(scip, multvars[i]);
245  	            if( SCIPisEQ(scip, multvarlb, multvarub) )
246  	               continue;
247  	
248  	            SCIP_CALL( updateBestCandidate(scip, branchruledata, bestvar, bestbrpoint, bestscore, bestrndscore,
249  	               multvars[i], candscoremin, candscoremax, candscoresum, candrndscore, SCIP_INVALID) );
250  	         }
251  	
252  	      assert(*bestvar != NULL); /* if all variables were fixed, something is strange */
253  	
254  	      return SCIP_OKAY;
255  	   }
256  	
257  	   /* select branching point for this variable */
258  	   candbrpoint = SCIPgetBranchingPoint(scip, cand, candsol);
259  	   assert(candbrpoint >= SCIPvarGetLbLocal(cand));
260  	   assert(candbrpoint <= SCIPvarGetUbLocal(cand));
261  	
262  	   /* we cannot branch on a huge value for a discrete variable, because we simply cannot enumerate such huge integer values in floating point
263  	    * arithmetics
264  	    */
265  	   if( SCIPvarGetType(cand) != SCIP_VARTYPE_CONTINUOUS && (SCIPisHugeValue(scip, candbrpoint) || SCIPisHugeValue(scip, -candbrpoint)) )
266  	      return SCIP_OKAY;
267  	
268  	   assert(SCIPvarGetType(cand) == SCIP_VARTYPE_CONTINUOUS || !SCIPisIntegral(scip, candbrpoint));
269  	
270  	   if( SCIPvarGetType(cand) == SCIP_VARTYPE_CONTINUOUS )
271  	      strategy = (branchruledata->strategy == 'u' ? branchruledata->updatestrategy : branchruledata->strategy);
272  	   else
273  	      strategy = (branchruledata->strategy == 'u' ? 'l' : branchruledata->strategy);
274  	
275  	   switch( strategy )
276  	   {
277  	   case 'l':
278  	      if( SCIPisInfinity(scip,  SCIPgetSolVal(scip, NULL, cand)) || SCIPgetSolVal(scip, NULL, cand) <= SCIPadjustedVarUb(scip, cand, candbrpoint) )
279  	         deltaminus = 0.0;
280  	      else
281  	         deltaminus = SCIPgetSolVal(scip, NULL, cand) - SCIPadjustedVarUb(scip, cand, candbrpoint);
282  	      if( SCIPisInfinity(scip, -SCIPgetSolVal(scip, NULL, cand)) || SCIPgetSolVal(scip, NULL, cand) >= SCIPadjustedVarLb(scip, cand, candbrpoint) )
283  	         deltaplus = 0.0;
284  	      else
285  	         deltaplus = SCIPadjustedVarLb(scip, cand, candbrpoint) - SCIPgetSolVal(scip, NULL, cand);
286  	      break;
287  	
288  	   case 'd':
289  	      if( SCIPisInfinity(scip, -SCIPvarGetLbLocal(cand)) )
290  	         deltaminus = SCIPisInfinity(scip, candscoremax) ? SCIPinfinity(scip) : WEIGHTEDSCORING(branchruledata, candscoremin, candscoremax, candscoresum);
291  	      else
292  	         deltaminus = SCIPadjustedVarUb(scip, cand, candbrpoint) - SCIPvarGetLbLocal(cand);
293  	
294  	      if( SCIPisInfinity(scip,  SCIPvarGetUbLocal(cand)) )
295  	         deltaplus = SCIPisInfinity(scip, candscoremax) ? SCIPinfinity(scip) : WEIGHTEDSCORING(branchruledata, candscoremin, candscoremax, candscoresum);
296  	      else
297  	         deltaplus = SCIPvarGetUbLocal(cand) - SCIPadjustedVarLb(scip, cand, candbrpoint);
298  	      break;
299  	
300  	   case 's':
301  	      if( SCIPisInfinity(scip, -SCIPvarGetLbLocal(cand)) )
302  	         deltaplus = SCIPisInfinity(scip, candscoremax) ? SCIPinfinity(scip) : WEIGHTEDSCORING(branchruledata, candscoremin, candscoremax, candscoresum);
303  	      else
304  	         deltaplus = SCIPadjustedVarUb(scip, cand, candbrpoint) - SCIPvarGetLbLocal(cand);
305  	
306  	      if( SCIPisInfinity(scip,  SCIPvarGetUbLocal(cand)) )
307  	         deltaminus = SCIPisInfinity(scip, candscoremax) ? SCIPinfinity(scip) : WEIGHTEDSCORING(branchruledata, candscoremin, candscoremax, candscoresum);
308  	      else
309  	         deltaminus = SCIPvarGetUbLocal(cand) - SCIPadjustedVarLb(scip, cand, candbrpoint);
310  	      break;
311  	
312  	   case 'v':
313  	      deltaplus = SCIPisInfinity(scip, candscoremax) ? SCIPinfinity(scip) : WEIGHTEDSCORING(branchruledata, candscoremin, candscoremax, candscoresum);
314  	      deltaminus = deltaplus;
315  	      break;
316  	
317  	   default :
318  	      SCIPerrorMessage("branching strategy %c unknown\n", strategy);
319  	      SCIPABORT();
320  	      return SCIP_INVALIDDATA;  /*lint !e527*/
321  	   }
322  	
323  	   if( SCIPisInfinity(scip, deltaminus) || SCIPisInfinity(scip, deltaplus) )
324  	   {
325  	      branchscore = SCIPinfinity(scip);
326  	   }
327  	   else
328  	   {
329  	      pscostdown  = SCIPgetVarPseudocostVal(scip, cand, -deltaminus);
330  	      pscostup    = SCIPgetVarPseudocostVal(scip, cand,  deltaplus);
331  	      branchscore = SCIPgetBranchScore(scip, cand, pscostdown, pscostup);
332  	      assert(!SCIPisNegative(scip, branchscore));
333  	   }
334  	   SCIPdebugMsg(scip, "branching score variable <%s>[%g,%g] = %g; wscore = %g; type=%d bestbrscore=%g\n",
335  	      SCIPvarGetName(cand), SCIPvarGetLbLocal(cand), SCIPvarGetUbLocal(cand), branchscore, WEIGHTEDSCORING(branchruledata, candscoremin, candscoremax, candscoresum),
336  	      SCIPvarGetType(cand), *bestscore);
337  	
338  	   if( SCIPisInfinity(scip, branchscore) )
339  	      branchscore = 0.9*SCIPinfinity(scip);
340  	
341  	   if( SCIPisSumGT(scip, branchscore, *bestscore) )
342  	   {
343  	      (*bestscore)    = branchscore;
344  	      (*bestrndscore) = candrndscore;
345  	      (*bestvar)      = cand;
346  	      (*bestbrpoint)  = candbrpoint;
347  	      return SCIP_OKAY;
348  	   }
349  	
350  	   /* if score of candidate is worse than bestscore, stay with best candidate */
351  	   if( !SCIPisSumEQ(scip, branchscore, *bestscore) )
352  	      return SCIP_OKAY;
353  	
354  	   if( SCIPisInfinity(scip, -SCIPvarGetLbLocal(*bestvar)) && SCIPisInfinity(scip, SCIPvarGetUbLocal(*bestvar)) )
355  	   {
356  	      /* best candidate is unbounded -> we prefer to branch on it */
357  	      if( SCIPisInfinity(scip, -SCIPvarGetLbLocal(cand)) && SCIPisInfinity(scip, SCIPvarGetUbLocal(cand)) &&
358  	          SCIPrandomGetReal(branchruledata->randnumgen, 0.0, 1.0) <= 0.5
359  	        )
360  	      {
361  	         /* but if the new candidate is also unbounded (thus as good as best candidate),
362  	          * then switch to the candidate with 50% probability to reduce performance variability
363  	          */
364  	         (*bestscore)    = branchscore;
365  	         (*bestrndscore) = candrndscore;
366  	         (*bestvar)      = cand;
367  	         (*bestbrpoint)  = candbrpoint;
368  	      }
369  	
370  	      return SCIP_OKAY;
371  	   }
372  	
373  	   /* best candidate has a finite lower or upper bound -> consider taking the other candidate */
374  	
375  	   if( (SCIPisInfinity(scip, -SCIPvarGetLbLocal(cand))     || SCIPisInfinity(scip, SCIPvarGetUbLocal(cand))) &&
376  	       (SCIPisInfinity(scip, -SCIPvarGetLbLocal(*bestvar)) || SCIPisInfinity(scip, SCIPvarGetUbLocal(*bestvar))) )
377  	   {
378  	      /* both candidates are unbounded, but one side may be finite (for bestcand we know there is one)
379  	       * take the candidate with the larger bound on the bounded side (hope that this avoids branching on always the same variable)
380  	       * this will also prefer unbounded variables over bounded ones
381  	       */
382  	      if( SCIPvarGetUbLocal(cand) > SCIPvarGetUbLocal(*bestvar) || SCIPvarGetLbLocal(cand) < SCIPvarGetLbLocal(*bestvar) )
383  	      {
384  	         /* cand is better than bestvar */
385  	         (*bestscore)    = branchscore;
386  	         (*bestrndscore) = candrndscore;
387  	         (*bestvar)      = cand;
388  	         (*bestbrpoint)  = candbrpoint;
389  	         return SCIP_OKAY;
390  	      }
391  	
392  	      if( SCIPvarGetUbLocal(*bestvar) > SCIPvarGetUbLocal(cand) || SCIPvarGetLbLocal(*bestvar) < SCIPvarGetLbLocal(cand) )
393  	      {
394  	         /* bestvar is better than cand */
395  	         return SCIP_OKAY;
396  	      }
397  	
398  	      /* both are equally good */
399  	   }
400  	
401  	   if( SCIPvarGetType(*bestvar) == SCIPvarGetType(cand) )
402  	   {
403  	      /* if both have the same type, take the one with larger relative diameter */
404  	      if( SCIPrelDiff(SCIPvarGetUbLocal(*bestvar), SCIPvarGetLbLocal(*bestvar)) < SCIPrelDiff(SCIPvarGetUbLocal(cand), SCIPvarGetLbLocal(cand)) )
405  	      {
406  	         /* cand has larger diameter than bestvar*/
407  	         (*bestscore)    = branchscore;
408  	         (*bestrndscore) = candrndscore;
409  	         (*bestvar)      = cand;
410  	         (*bestbrpoint)  = candbrpoint;
411  	         return SCIP_OKAY;
412  	      }
413  	
414  	      if( SCIPrelDiff(SCIPvarGetUbLocal(*bestvar), SCIPvarGetLbLocal(*bestvar)) > SCIPrelDiff(SCIPvarGetUbLocal(cand), SCIPvarGetLbLocal(cand)) )
415  	      {
416  	         /* bestvar has larger diameter than cand */
417  	         return SCIP_OKAY;
418  	      }
419  	   }
420  	
421  	   /* take the one with better type ("more discrete") */
422  	   if( SCIPvarGetType(*bestvar) > SCIPvarGetType(cand) )
423  	   {
424  	      /* cand is more discrete than bestvar */
425  	      (*bestscore)    = branchscore;
426  	      (*bestrndscore) = candrndscore;
427  	      (*bestvar)      = cand;
428  	      (*bestbrpoint)  = candbrpoint;
429  	      return SCIP_OKAY;
430  	   }
431  	   if( SCIPvarGetType(*bestvar) < SCIPvarGetType(cand) )
432  	   {
433  	      /* bestvar is more discrete than cand */
434  	      return SCIP_OKAY;
435  	   }
436  	
437  	   /* cand seems to be as good as the currently best one (bestvar); use the random score as a final tie-breaker */
438  	   if( candrndscore >= (*bestrndscore) )
439  	   {
440  	      (*bestscore)    = branchscore;
441  	      (*bestrndscore) = candrndscore;
442  	      (*bestvar)      = cand;
443  	      (*bestbrpoint)  = candbrpoint;
444  	   }
445  	
446  	   return SCIP_OKAY;
447  	}
448  	
449  	/** selects the branching variable from given candidate array */
450  	static
451  	SCIP_RETCODE selectBranchVar(
452  	   SCIP*                 scip,               /**< SCIP data structure */
453  	   SCIP_BRANCHRULE*      branchrule,         /**< branching rule */
454  	   SCIP_VAR**            cands,              /**< array of branching candidates */
455  	   SCIP_Real*            candssol,           /**< array of candidate solution values */
456  	   SCIP_Real*            candsscore,         /**< array of candidate scores */
457  	   int                   ncands,             /**< the number of candidates */
458  	   SCIP_VAR**            brvar,              /**< pointer to store the selected branching candidate or NULL if none */
459  	   SCIP_Real*            brpoint             /**< pointer to store branching point of selected branching variable */
460  	   )
461  	{ /*lint --e{850}*/ 
462  	   SCIP_BRANCHRULEDATA* branchruledata;
463  	
464  	   SCIP_VAR* cand;
465  	   SCIP_Real candsol;
466  	
467  	   SCIP_Real bestbranchscore;
468  	   SCIP_Real bestrndscore;
469  	
470  	   SCIP_Real scoremin;
471  	   SCIP_Real scoresum;
472  	   SCIP_Real scoremax;
473  	
474  	   SCIP_VAR** candssorted;
475  	   int* candsorigidx;
476  	
477  	   int i;
478  	   int j;
479  	
480  	   assert(brvar   != NULL);
481  	   assert(brpoint != NULL);
482  	
483  	   (*brvar)   = NULL;
484  	   (*brpoint) = SCIP_INVALID;
485  	
486  	   if( ncands == 0 )
487  	      return SCIP_OKAY;
488  	
489  	   branchruledata = SCIPbranchruleGetData(branchrule);
490  	   assert(branchruledata != NULL);
491  	
492  	   /* sort branching candidates (in a copy), such that same variables are on consecutive positions */
493  	   SCIP_CALL( SCIPduplicateBufferArray(scip, &candssorted, cands, ncands) );
494  	   SCIP_CALL( SCIPallocBufferArray(scip, &candsorigidx, ncands) );
495  	   for( i = 0; i < ncands; ++i )
496  	      candsorigidx[i] = i;
497  	
498  	   SCIPsortPtrInt((void**)candssorted, candsorigidx, SCIPvarComp, ncands);
499  	
500  	   bestbranchscore = -1.0;
501  	   bestrndscore = -1.0;
502  	
503  	   for( i = 0; i < ncands; ++i )
504  	   {
505  	      cand = candssorted[i];
506  	
507  	      /* there should be no fixed branching candidates */
508  	      assert(!SCIPisEQ(scip, SCIPvarGetLbLocal(cand), SCIPvarGetUbLocal(cand)));
509  	
510  	      /* compute min, sum, and max of all registered scores for this variables
511  	       * set candsol to a valid value, if someone registered one */
512  	      scoremin = candsscore[candsorigidx[i]];
513  	      scoresum = scoremin;
514  	      scoremax = scoremin;
515  	      candsol  = candssol[candsorigidx[i]];
516  	      for( j = i+1 ; j < ncands && SCIPvarCompare(candssorted[j], cand) == 0; ++j )
517  	      {
518  	         assert(candsscore[candsorigidx[j]] >= 0.0);
519  	         scoresum += candsscore[candsorigidx[j]];
520  	         if( candsscore[candsorigidx[j]] < scoremin )
521  	            scoremin = candsscore[candsorigidx[j]];
522  	         else if( candsscore[candsorigidx[j]] > scoremax )
523  	            scoremax = candsscore[candsorigidx[j]];
524  	
525  	         /* @todo if there are two valid externcandssol available for the same variable, should we take the one closer to the middle of the domain? */
526  	         if( SCIPisInfinity(scip, REALABS(candsol)) )
527  	            candsol = candssol[candsorigidx[j]];
528  	      }
529  	      /* set i to last occurrence of cand in candssorted (instead of first one as before), so in next round we look at another variable */
530  	      i = j-1;
531  	      assert(candssorted[i] == cand);
532  	
533  	      /* check if new candidate is better than previous candidate (if any) */
534  	      SCIP_CALL( updateBestCandidate(scip, branchruledata, brvar, brpoint, &bestbranchscore, &bestrndscore, cand,
535  	            scoremin, scoremax, scoresum, SCIPrandomGetReal(branchruledata->randnumgen, 0.0, 1.0), candsol) );
536  	   }
537  	
538  	   /* there were candidates, but no variable was selected; this can only happen if the branching points are huge values
539  	    * for all non-continuous variables on which we cannot branch
540  	    * @todo delay the node?
541  	    */
542  	   if( (*brvar) == NULL )
543  	   {
544  	      SCIPerrorMessage("no branching could be created: all external candidates have huge bounds\n");
545  	      return SCIP_BRANCHERROR; /*lint !e527*/
546  	   }
547  	
548  	   /* free buffer arrays */
549  	   SCIPfreeBufferArray(scip, &candsorigidx);
550  	   SCIPfreeBufferArray(scip, &candssorted);
551  	
552  	   return SCIP_OKAY;
553  	}
554  	
555  	/*
556  	 * Callback methods
557  	 */
558  	
559  	/** copy method for branchrule plugins (called when SCIP copies plugins) */
560  	static
561  	SCIP_DECL_BRANCHCOPY(branchCopyPscost)
562  	{  /*lint --e{715}*/
563  	   assert(scip != NULL);
564  	   assert(branchrule != NULL);
565  	   assert(strcmp(SCIPbranchruleGetName(branchrule), BRANCHRULE_NAME) == 0);
566  	
567  	   /* call inclusion method of branchrule */
568  	   SCIP_CALL( SCIPincludeBranchrulePscost(scip) );
569  	
570  	   return SCIP_OKAY;
571  	}
572  	
573  	/** destructor of branching rule to free user data (called when SCIP is exiting) */
574  	static
575  	SCIP_DECL_BRANCHFREE(branchFreePscost)
576  	{  /*lint --e{715}*/
577  	   SCIP_BRANCHRULEDATA* branchruledata;
578  	
579  	   /* get branching rule data */
580  	   branchruledata = SCIPbranchruleGetData(branchrule);
581  	   assert(branchruledata != NULL);
582  	
583  	   /* free random number generator */
584  	   SCIPfreeRandom(scip, &branchruledata->randnumgen);
585  	
586  	   /* free branching rule data */
587  	   SCIPfreeBlockMemory(scip, &branchruledata);
588  	   SCIPbranchruleSetData(branchrule, NULL);
589  	
590  	   return SCIP_OKAY;
591  	}
592  	
593  	/** initialization method of branching rule (called after problem was transformed) */
594  	static
595  	SCIP_DECL_BRANCHINIT(branchInitPscost)
596  	{  /*lint --e{715}*/
597  	   SCIP_BRANCHRULEDATA* branchruledata;
598  	
599  	   /* initialize branching rule data */
600  	   branchruledata = SCIPbranchruleGetData(branchrule);
601  	   assert(branchruledata != NULL);
602  	
603  	   SCIPsetRandomSeed(scip, branchruledata->randnumgen, BRANCHRULE_RANDSEED_DEFAULT);
604  	
605  	   return SCIP_OKAY;
606  	}
607  	
608  	/** branching execution method for fractional LP solutions */
609  	static
610  	SCIP_DECL_BRANCHEXECLP(branchExeclpPscost)
611  	{  /*lint --e{715}*/
612  	   SCIP_VAR** lpcands;
613  	   SCIP_Real* lpcandssol;
614  	   SCIP_Real bestscore;
615  	   SCIP_Real bestrootdiff;
616  	   int nlpcands;
617  	   int bestcand;
618  	   int c;
619  	
620  	   assert(branchrule != NULL);
621  	   assert(strcmp(SCIPbranchruleGetName(branchrule), BRANCHRULE_NAME) == 0);
622  	   assert(scip != NULL);
623  	   assert(result != NULL);
624  	
625  	   SCIPdebugMsg(scip, "Execlp method of pscost branching\n");
626  	
627  	   /* get branching candidates */
628  	   SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, NULL, NULL, &nlpcands, NULL) );
629  	   assert(nlpcands > 0);
630  	
631  	   bestcand = -1;
632  	   bestscore = -SCIPinfinity(scip);
633  	   bestrootdiff = 0.0;
634  	   for( c = 0; c < nlpcands; ++c )
635  	   {
636  	      SCIP_Real score;
637  	      SCIP_Real rootsolval;
638  	      SCIP_Real rootdiff;
639  	
640  	      score = SCIPgetVarPseudocostScore(scip, lpcands[c], lpcandssol[c]);
641  	      rootsolval = SCIPvarGetRootSol(lpcands[c]);
642  	      rootdiff = REALABS(lpcandssol[c] - rootsolval);
643  	      if( SCIPisSumGT(scip, score, bestscore) || (SCIPisSumEQ(scip, score, bestscore) && rootdiff > bestrootdiff) )
644  	      {
645  	         bestcand = c;
646  	         bestscore = score;
647  	         bestrootdiff = rootdiff;
648  	      }
649  	   }
650  	   assert(0 <= bestcand && bestcand < nlpcands);
651  	   assert(!SCIPisFeasIntegral(scip, lpcandssol[bestcand]));
652  	   assert(!SCIPisFeasIntegral(scip, SCIPvarGetSol(lpcands[bestcand], TRUE)));
653  	
654  	   /* perform the branching */
655  	   SCIPdebugMsg(scip, " -> %d cands, selected cand %d: variable <%s> (solval=%g, score=%g)\n",
656  	      nlpcands, bestcand, SCIPvarGetName(lpcands[bestcand]), lpcandssol[bestcand], bestscore);
657  	
658  	   /* perform the branching */
659  	   SCIP_CALL( SCIPbranchVar(scip, lpcands[bestcand], NULL, NULL, NULL) );
660  	   *result = SCIP_BRANCHED;
661  	
662  	   return SCIP_OKAY;
663  	}
664  	
665  	
666  	/** branching execution method for external candidates */
667  	static
668  	SCIP_DECL_BRANCHEXECEXT(branchExecextPscost)
669  	{  /*lint --e{715}*/
670  	   SCIP_BRANCHRULEDATA* branchruledata;
671  	   SCIP_VAR** externcands;
672  	   SCIP_Real* externcandssol;
673  	   SCIP_Real* externcandsscore;
674  	   int nprioexterncands;
675  	   SCIP_VAR* brvar;
676  	   SCIP_Real brpoint;
677  	   int nchildren;
678  	
679  	   assert(branchrule != NULL);
680  	   assert(strcmp(SCIPbranchruleGetName(branchrule), BRANCHRULE_NAME) == 0);
681  	   assert(scip != NULL);
682  	   assert(result != NULL);
683  	
684  	   branchruledata = SCIPbranchruleGetData(branchrule);
685  	   assert(branchruledata != NULL);
686  	
687  	   SCIPdebugMsg(scip, "Execext method of pscost branching\n");
688  	
689  	   /* get branching candidates */
690  	   SCIP_CALL( SCIPgetExternBranchCands(scip, &externcands, &externcandssol, &externcandsscore, NULL, &nprioexterncands, NULL, NULL, NULL) );
691  	   assert(nprioexterncands > 0);
692  	
693  	   /* get current update strategy for pseudo costs, if our multiplier rule is 'u' */
694  	   if( branchruledata->strategy == 'u' )
695  	   {
696  	      SCIP_CALL( SCIPgetCharParam(scip, "branching/lpgainnormalize", &branchruledata->updatestrategy) );
697  	   }
698  	
699  	   /* select branching variable */
700  	   SCIP_CALL( selectBranchVar(scip, branchrule, externcands, externcandssol, externcandsscore, nprioexterncands, &brvar, &brpoint) );
701  	
702  	   if( brvar == NULL )
703  	   {
704  	      /* can happen if all candidates were non-continous variables with huge bounds */
705  	      *result = SCIP_DIDNOTFIND;
706  	      return SCIP_OKAY;
707  	   }
708  	
709  	   assert(SCIPvarIsActive(SCIPvarGetProbvar(brvar)));
710  	
711  	   SCIPdebugMsg(scip, "branching on variable <%s>: new intervals: [%g, %g] and [%g, %g]\n",
712  	      SCIPvarGetName(brvar), SCIPvarGetLbLocal(brvar), SCIPadjustedVarUb(scip, brvar, brpoint), SCIPadjustedVarLb(scip, brvar, brpoint), SCIPvarGetUbLocal(brvar));
713  	
714  	   if( branchruledata->nchildren > 2 && SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) <= branchruledata->narymaxdepth )
715  	   {
716  	      /* do n-ary branching */
717  	      SCIP_Real minwidth;
718  	
719  	      minwidth = 0.0;
720  	      if( !SCIPisInfinity(scip, -SCIPvarGetLbGlobal(brvar)) && !SCIPisInfinity(scip, SCIPvarGetUbGlobal(brvar)) )
721  	         minwidth = branchruledata->naryminwidth * (SCIPvarGetUbGlobal(brvar) - SCIPvarGetLbGlobal(brvar));
722  	
723  	      SCIP_CALL( SCIPbranchVarValNary(scip, brvar, brpoint, branchruledata->nchildren, minwidth, branchruledata->narywidthfactor, &nchildren) );
724  	   }
725  	   else
726  	   {
727  	      /* do binary branching */
728  	      SCIP_CALL( SCIPbranchVarValNary(scip, brvar, brpoint, 2, 0.0, 1.0, &nchildren) );
729  	   }
730  	
731  	   if( nchildren > 1 )
732  	   {
733  	      *result = SCIP_BRANCHED;
734  	   }
735  	   else
736  	   {
737  	      /* if there are no children, then variable should have been fixed by SCIPbranchVarVal */
738  	      assert(SCIPisEQ(scip, SCIPvarGetLbLocal(brvar), SCIPvarGetUbLocal(brvar)));
739  	      *result = SCIP_REDUCEDDOM;
740  	   }
741  	
742  	   return SCIP_OKAY;
743  	}
744  	
745  	
746  	/*
747  	 * branching specific interface methods
748  	 */
749  	
750  	/** creates the pseudo cost branching rule and includes it in SCIP */
751  	SCIP_RETCODE SCIPincludeBranchrulePscost(
752  	   SCIP*                 scip                /**< SCIP data structure */
753  	   )
754  	{
755  	   SCIP_BRANCHRULEDATA* branchruledata;
756  	   SCIP_BRANCHRULE* branchrule;
757  	
758  	   /* create pscost branching rule data */
759  	   SCIP_CALL( SCIPallocBlockMemory(scip, &branchruledata) );
760  	
761  	   /* include allfullstrong branching rule */
762  	   SCIP_CALL( SCIPincludeBranchruleBasic(scip, &branchrule, BRANCHRULE_NAME, BRANCHRULE_DESC, BRANCHRULE_PRIORITY,
763  	         BRANCHRULE_MAXDEPTH, BRANCHRULE_MAXBOUNDDIST, branchruledata) );
764  	
765  	   assert(branchrule != NULL);
766  	   /* create a random number generator */
767  	   SCIP_CALL( SCIPcreateRandom(scip, &branchruledata->randnumgen,
768  	         BRANCHRULE_RANDSEED_DEFAULT, TRUE) );
769  	
770  	   /* set non-fundamental callbacks via specific setter functions*/
771  	   SCIP_CALL( SCIPsetBranchruleCopy(scip, branchrule, branchCopyPscost) );
772  	   SCIP_CALL( SCIPsetBranchruleFree(scip, branchrule, branchFreePscost) );
773  	   SCIP_CALL( SCIPsetBranchruleInit(scip, branchrule, branchInitPscost) );
774  	   SCIP_CALL( SCIPsetBranchruleExecLp(scip, branchrule, branchExeclpPscost) );
775  	   SCIP_CALL( SCIPsetBranchruleExecExt(scip, branchrule, branchExecextPscost) );
776  	
777  	   SCIP_CALL( SCIPaddCharParam(scip, "branching/" BRANCHRULE_NAME "/strategy",
778  	         "strategy for utilizing pseudo-costs of external branching candidates (multiply as in pseudo costs 'u'pdate rule, or by 'd'omain reduction, or by domain reduction of 's'ibling, or by 'v'ariable score)",
779  	         &branchruledata->strategy, FALSE, BRANCHRULE_STRATEGY_DEFAULT, BRANCHRULE_STRATEGIES, NULL, NULL) );
780  	
781  	   SCIP_CALL( SCIPaddRealParam(scip, "branching/" BRANCHRULE_NAME "/minscoreweight",
782  	         "weight for minimum of scores of a branching candidate when building weighted sum of min/max/sum of scores",
783  	         &branchruledata->scoreminweight, TRUE, BRANCHRULE_SCOREMINWEIGHT_DEFAULT, -SCIPinfinity(scip), SCIPinfinity(scip), NULL, NULL) );
784  	
785  	   SCIP_CALL( SCIPaddRealParam(scip, "branching/" BRANCHRULE_NAME "/maxscoreweight",
786  	         "weight for maximum of scores of a branching candidate when building weighted sum of min/max/sum of scores",
787  	         &branchruledata->scoremaxweight, TRUE, BRANCHRULE_SCOREMAXWEIGHT_DEFAULT, -SCIPinfinity(scip), SCIPinfinity(scip), NULL, NULL) );
788  	
789  	   SCIP_CALL( SCIPaddRealParam(scip, "branching/" BRANCHRULE_NAME "/sumscoreweight",
790  	         "weight for sum of scores of a branching candidate when building weighted sum of min/max/sum of scores",
791  	         &branchruledata->scoresumweight, TRUE, BRANCHRULE_SCORESUMWEIGHT_DEFAULT, -SCIPinfinity(scip), SCIPinfinity(scip), NULL, NULL) );
792  	
793  	   SCIP_CALL( SCIPaddIntParam(scip, "branching/" BRANCHRULE_NAME "/nchildren",
794  	         "number of children to create in n-ary branching",
795  	         &branchruledata->nchildren, FALSE, BRANCHRULE_NCHILDREN_DEFAULT, 2, INT_MAX, NULL, NULL) );
796  	
797  	   SCIP_CALL( SCIPaddIntParam(scip, "branching/" BRANCHRULE_NAME "/narymaxdepth",
798  	         "maximal depth where to do n-ary branching, -1 to turn off",
799  	         &branchruledata->narymaxdepth, FALSE, BRANCHRULE_NARYMAXDEPTH_DEFAULT, -1, INT_MAX, NULL, NULL) );
800  	
801  	   SCIP_CALL( SCIPaddRealParam(scip, "branching/" BRANCHRULE_NAME "/naryminwidth",
802  	         "minimal domain width in children when doing n-ary branching, relative to global bounds",
803  	         &branchruledata->naryminwidth, FALSE, BRANCHRULE_NARYMINWIDTH_DEFAULT, 0.0, 1.0, NULL, NULL) );
804  	
805  	   SCIP_CALL( SCIPaddRealParam(scip, "branching/" BRANCHRULE_NAME "/narywidthfactor",
806  	         "factor of domain width in n-ary branching when creating nodes with increasing distance from branching value",
807  	         &branchruledata->narywidthfactor, FALSE, BRANCHRULE_NARYWIDTHFAC_DEFAULT, 1.0, SCIP_REAL_MAX, NULL, NULL) );
808  	
809  	   return SCIP_OKAY;
810  	}
811  	
812  	/** selects a branching variable, due to pseudo cost, from the given candidate array and returns this variable together
813  	 *  with a branching point */
814  	SCIP_RETCODE SCIPselectBranchVarPscost(
815  	   SCIP*                 scip,               /**< SCIP data structure */
816  	   SCIP_VAR**            branchcands,        /**< branching candidates */
817  	   SCIP_Real*            branchcandssol,     /**< solution value for the branching candidates */
818  	   SCIP_Real*            branchcandsscore,   /**< array of candidate scores */
819  	   int                   nbranchcands,       /**< number of branching candidates */
820  	   SCIP_VAR**            var,                /**< pointer to store the variable to branch on, or NULL if none */
821  	   SCIP_Real*            brpoint             /**< pointer to store the branching point for the branching variable, will be fractional for a discrete variable */
822  	   )
823  	{
824  	   SCIP_BRANCHRULE* branchrule;
825  	
826  	   assert(scip != NULL);
827  	
828  	   /* find branching rule */
829  	   branchrule = SCIPfindBranchrule(scip, BRANCHRULE_NAME);
830  	   assert(branchrule != NULL);
831  	
832  	   /* select branching variable */
833  	   SCIP_CALL( selectBranchVar(scip, branchrule, branchcands, branchcandssol, branchcandsscore, nbranchcands, var, brpoint) );
834  	
835  	   return SCIP_OKAY;
836  	}
837