1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the program and library             */
4    	/*         SCIP --- Solving Constraint Integer Programs                      */
5    	/*                                                                           */
6    	/*  Copyright 2002-2022 Zuse Institute Berlin                                */
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   cutsel_ensemble.c
26   	 * @ingroup DEFPLUGINS_CUTSEL
27   	 * @brief  ensemble cut selector
28   	 * @author Mark Turner
29   	 *
30   	 * @todo separator hard limit on density is inappropriate for MINLP. Need to relax hard limit in case of all cuts dense
31   	 * @todo penalising via parallelism is overly costly if many cuts. Hash cuts before and find appropriate groups?
32   	 */
33   	
34   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
35   	
36   	#include <assert.h>
37   	
38   	#include "scip/scip_cutsel.h"
39   	#include "scip/scip_cut.h"
40   	#include "scip/scip_lp.h"
41   	#include "scip/scip_randnumgen.h"
42   	#include "scip/cutsel_ensemble.h"
43   	
44   	
45   	#define CUTSEL_NAME              "ensemble"
46   	#define CUTSEL_DESC              "weighted sum of many terms with optional filtering and penalties"
47   	#define CUTSEL_PRIORITY           7000
48   	
49   	#define RANDSEED                  0x5EED
50   	
51   	#define DEFAULT_MINSCORE                0.0  /**< minimum score s.t. a cut can be selected */
52   	#define DEFAULT_EFFICACYWEIGHT          0.75 /**< weight of normed-efficacy in score calculation */
53   	#define DEFAULT_DIRCUTOFFDISTWEIGHT     0.0  /**< weight of normed-directed cutoff distance in score calculation */
54   	#define DEFAULT_OBJPARALWEIGHT          0.25 /**< weight of objective parallelism in score calculation */
55   	#define DEFAULT_INTSUPPORTWEIGHT        0.45 /**< weight of integral support in cut score calculation */
56   	#define DEFAULT_EXPIMPROVWEIGHT         0.1  /**< weight of normed-expected improvement in cut score calculation */
57   	#define DEFAULT_PSCOSTWEIGHT            0.75 /**< weight of normalised pseudo-costs in cut score calculation */
58   	#define DEFAULT_NLOCKSWEIGHT            0.25 /**< weight of normalised number of locks in cut score calculation */
59   	#define DEFAULT_MAXSPARSITYBONUS        0.5  /**< score given to a cut with complete sparsity */
60   	#define DEFAULT_SPARSITYENDBONUS        0.2  /**< the density at which a cut no longer receives additional score */
61   	#define DEFAULT_GOODNUMERICBONUS        0.0  /**< bonus provided for good numerics */
62   	#define DEFAULT_MAXCOEFRATIOBONUS       10000 /**< maximum coefficient ratio of cut for which numeric bonus is given */
63   	#define DEFAULT_PENALISELOCKS           TRUE /**< whether having less locks should be rewarded instead of more */
64   	#define DEFAULT_PENALISEOBJPARAL        TRUE /**< whether objective parallelism should be penalised not rewarded */
65   	#define DEFAULT_FILTERPARALCUTS         FALSE /**< should cuts be filtered so no two parallel cuts are added */
66   	#define DEFAULT_MAXPARAL                0.95 /**< threshold for when two cuts are considered parallel to each other */
67   	#define DEFAULT_PENALISEPARALCUTS       TRUE /**< should two parallel cuts be penalised instead of outright filtered */
68   	#define DEFAULT_PARALPENALTY            0.25 /**< penalty for weaker of two parallel cuts if penalising parallel cuts */
69   	#define DEFAULT_FILTERDENSECUTS         TRUE /**< should cuts over a given density threshold be filtered */
70   	#define DEFAULT_MAXCUTDENSITY           0.425 /**< max allowed cut density if filtering dense cuts */
71   	#define DEFAULT_MAXNONZEROROOTROUND     4.5  /**< max nonzeros per round (root). Gets multiplied by num LP cols */
72   	#define DEFAULT_MAXNONZEROTREEROUND     9.5  /**< max nonzeros per round (tree). Gets multiplied by num LP cols */
73   	#define DEFAULT_MAXCUTS                 200  /**< maximum number of cuts that can be considered by this cut selector */
74   	#define DEFAULT_MAXNUMVARS              50000 /**< maximum number of variables that a problem can have while calling this cut selector */
75   	
76   	/*
77   	 * Data structures
78   	 */
79   	
80   	/** cut selector data */
81   	struct SCIP_CutselData
82   	{
83   	   SCIP_RANDNUMGEN*      randnumgen;         /**< random generator for tie-breaking */
84   	   SCIP_Real             minscore;           /**< minimum score for a cut to be added to the LP */
85   	   SCIP_Real             objparalweight;     /**< weight of objective parallelism in cut score calculation */
86   	   SCIP_Real             efficacyweight;     /**< weight of normed-efficacy in cut score calculation */
87   	   SCIP_Real             dircutoffdistweight;/**< weight of normed-directed cutoff distance in cut score calculation */
88   	   SCIP_Real             expimprovweight;    /**< weight of normed-expected improvement in cut score calculation */
89   	   SCIP_Real             intsupportweight;   /**< weight of integral support in cut score calculation */
90   	   SCIP_Real             pscostweight;       /**< weight of normalised pseudo-costs in cut score calculation */
91   	   SCIP_Real             locksweight;        /**< weight of normed-number of active locks in cut score calculation */
92   	   SCIP_Real             maxsparsitybonus;   /**< weight of maximum sparsity reward in cut score calculation */
93   	   SCIP_Real             goodnumericsbonus;  /**< weight of good numeric bonus in cut score calculation */
94   	   SCIP_Real             endsparsitybonus;   /**< max sparsity value for which a bonus is applied */
95   	   SCIP_Real             maxparal;           /**< threshold for when two cuts are considered parallel to each other */
96   	   SCIP_Real             paralpenalty;       /**< penalty for weaker of two parallel cuts if penalising parallel cuts */
97   	   SCIP_Real             maxcutdensity;      /**< max allowed cut density if filtering dense cuts */
98   	   SCIP_Real             maxnonzerorootround;/**< max nonzeros per round (root). Gets multiplied by num LP cols */
99   	   SCIP_Real             maxnonzerotreeround;/**< max nonzeros per round (tree). Gets multiplied by num LP cols */
100  	   SCIP_Bool             filterparalcuts;   /**< should cuts be filtered so no two parallel cuts are added */
101  	   SCIP_Bool             penaliseparalcuts;  /**< should two parallel cuts be penalised instead of outright filtered */
102  	   SCIP_Bool             filterdensecuts;    /**< should cuts over a given density threshold be filtered */
103  	   SCIP_Bool             penaliselocks;      /**< whether the number of locks should be penalised instead of rewarded */
104  	   SCIP_Bool             penaliseobjparal;   /**< whether objective parallelism should be penalised */
105  	   int                   maxcoefratiobonus;  /**< maximum coefficient ratio for which numeric bonus is applied */
106  	   int                   maxcuts;            /**< maximum number of cuts that can be considered by this cut selector */
107  	   int                   maxnumvars;         /**< maximum number of variables that a problem can have while calling this cut selector */
108  	};
109  	
110  	
111  	/*
112  	 * Local methods
113  	 */
114  	
115  	/** returns the maximum score of cuts; if scores is not NULL, then stores the individual score of each cut in scores */
116  	static
117  	SCIP_RETCODE scoring(
118  	   SCIP*                 scip,               /**< SCIP data structure */
119  	   SCIP_ROW**            cuts,               /**< array with cuts to score */
120  	   SCIP_CUTSELDATA*      cutseldata,         /**< cut selector data */
121  	   SCIP_Real*            scores,             /**< array to store the score of cuts or NULL */
122  	   SCIP_Bool             root,               /**< whether we are at the root node or not */
123  	   int                   ncuts               /**< number of cuts in cuts array */
124  	)
125  	{
126  	   SCIP_Real* effs;
127  	   SCIP_Real* dcds;
128  	   SCIP_Real* exps;
129  	   SCIP_Real* cutdensities;
130  	   SCIP_Real* cutlocks;
131  	   SCIP_Real* pscosts;
132  	   SCIP_SOL* sol;
133  	   SCIP_Real maxdcd = 0.0;
134  	   SCIP_Real maxeff = 0.0;
135  	   SCIP_Real maxexp = 0.0;
136  	   SCIP_Real maxpscost = 0.0;
137  	   SCIP_Real maxlocks = 0.0;
138  	   SCIP_Real ncols;
139  	
140  	   /* Get the solution that we use for directed cutoff distance calculations. Get the number of columns too */
141  	   sol = SCIPgetBestSol(scip);
142  	   ncols = SCIPgetNLPCols(scip);
143  	
144  	   /* Initialise all array information that we're going to use for scoring */
145  	   SCIP_CALL(SCIPallocBufferArray(scip, &effs, ncuts));
146  	   SCIP_CALL(SCIPallocBufferArray(scip, &dcds, ncuts));
147  	   SCIP_CALL(SCIPallocBufferArray(scip, &exps, ncuts));
148  	   SCIP_CALL(SCIPallocBufferArray(scip, &cutdensities, ncuts));
149  	   SCIP_CALL(SCIPallocBufferArray(scip, &cutlocks, ncuts));
150  	   SCIP_CALL(SCIPallocBufferArray(scip, &pscosts, ncuts));
151  	
152  	
153  	   /* Populate the number of cut locks, the pseudo-cost scores, and the cut densities */
154  	   for (int i = 0; i < ncuts; ++i )
155  	   {
156  	      SCIP_COL** cols;
157  	      SCIP_Real* cutvals;
158  	      SCIP_Real sqrcutnorm;
159  	      SCIP_Real ncutcols;
160  	      SCIP_Real cutalpha;
161  	
162  	      cols = SCIProwGetCols(cuts[i]);
163  	      cutvals = SCIProwGetVals(cuts[i]);
164  	      sqrcutnorm = MAX(SCIPsumepsilon(scip), SQR(SCIProwGetNorm(cuts[i]))); /*lint !e666*/
165  	      cutalpha = -SCIPgetRowFeasibility(scip, cuts[i]) / sqrcutnorm;
166  	      ncutcols = SCIProwGetNNonz(cuts[i]);
167  	      cutdensities[i] = ncutcols / ncols;
168  	      cutlocks[i] = 0;
169  	      pscosts[i] = 0;
170  	
171  	
172  	      for ( int j = 0; j < (int) ncutcols; ++j )
173  	      {
174  	         SCIP_VAR* colvar;
175  	         SCIP_Real colval;
176  	         SCIP_Real l1dist;
177  	
178  	         colval = SCIPcolGetPrimsol(cols[j]);
179  	         colvar = SCIPcolGetVar(cols[j]);
180  	         /* Get the number of active locks feature in the cut */
181  	         if( ! SCIPisInfinity(scip, SCIProwGetRhs(cuts[i])) && cutvals[j] > 0.0 )
182  	            cutlocks[i] += SCIPvarGetNLocksUp(colvar);
183  	         if( ! SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i])) && cutvals[j] < 0.0 )
184  	            cutlocks[i] += SCIPvarGetNLocksUp(colvar);
185  	         if( ! SCIPisInfinity(scip, SCIProwGetRhs(cuts[i])) && cutvals[j] < 0.0 )
186  	            cutlocks[i] += SCIPvarGetNLocksDown(colvar);
187  	         if( ! SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i])) && cutvals[j] > 0.0 )
188  	            cutlocks[i] += SCIPvarGetNLocksDown(colvar);
189  	
190  	         /* Get the L1 distance from the projection onto the cut and the LP solution in the variable direction */
191  	         l1dist = ABS(colval - (cutalpha * cutvals[j]));
192  	         pscosts[i] += SCIPgetVarPseudocostScore(scip, colvar, colval) * l1dist;
193  	      }
194  	      cutlocks[i] = cutlocks[i] / ncutcols; /*lint !e414*/
195  	
196  	      if( cutlocks[i] > maxlocks )
197  	         maxlocks = cutlocks[i];
198  	
199  	      if( pscosts[i] > maxpscost )
200  	         maxpscost = pscosts[i];
201  	   }
202  	   
203  	   /* account for the case where maxlocks or maxpscost is 0 */
204  	   maxpscost = MAX(maxpscost, SCIPepsilon(scip)); /*lint !e666*/
205  	   maxlocks = MAX(maxlocks, 1);
206  	
207  	   for ( int i = 0; i < ncuts; i++ )
208  	   {
209  	      cutlocks[i] = cutlocks[i] / maxlocks; /*lint !e414*/
210  	      /* if locks are penalized, we complement the corresponding score */
211  	      if( cutseldata->penaliselocks )
212  	         cutlocks[i] = 1 - cutlocks[i];
213  	      pscosts[i] = pscosts[i] / maxpscost; /*lint !e414*/
214  	   }
215  	
216  	
217  	   /* Get the arrays / maximums of directed cutoff distance, efficacy, and expected improvement values. */
218  	   if ( sol != NULL && root )
219  	   {
220  	      for ( int i = 0; i < ncuts; i++ )
221  	      {
222  	         dcds[i] = SCIPgetCutLPSolCutoffDistance(scip, sol, cuts[i]);
223  	         maxdcd = MAX(maxdcd, dcds[i]);
224  	      }
225  	   }
226  	
227  	   for ( int i = 0; i < ncuts; ++i )
228  	   {
229  	      effs[i] = SCIPgetCutEfficacy(scip, NULL, cuts[i]);
230  	      exps[i] = effs[i] * SCIPgetRowObjParallelism(scip, cuts[i]);
231  	      maxeff = MAX(maxeff, effs[i]);
232  	      maxexp = MAX(maxexp, exps[i]);
233  	   }
234  	
235  	   /* Now score the cuts */
236  	   for ( int i = 0; i < ncuts; ++i )
237  	   {
238  	      SCIP_Real score;
239  	      SCIP_Real scaleddcd;
240  	      SCIP_Real scaledeff;
241  	      SCIP_Real scaledexp;
242  	      SCIP_Real objparallelism;
243  	      SCIP_Real intsupport;
244  	      SCIP_Real density;
245  	      SCIP_Real dynamism;
246  	      SCIP_Real mincutval;
247  	      SCIP_Real maxcutval;
248  	      SCIP_Real pscost;
249  	      SCIP_Real cutlock;
250  	
251  	      /* Get the integer support */
252  	      intsupport = SCIPgetRowNumIntCols(scip, cuts[i]) / (SCIP_Real) SCIProwGetNNonz(cuts[i]);
253  	      intsupport *= cutseldata->intsupportweight;
254  	
255  	      /* Get the objective parallelism and orthogonality */
256  	      if( ! cutseldata->penaliseobjparal )
257  	         objparallelism = cutseldata->objparalweight * SCIPgetRowObjParallelism(scip, cuts[i]);
258  	      else
259  	         objparallelism = cutseldata->objparalweight * (1 - SCIPgetRowObjParallelism(scip, cuts[i]));
260  	
261  	      /* Get the density score */
262  	      density = (cutseldata->maxsparsitybonus / cutseldata->endsparsitybonus) * -1 * cutdensities[i];
263  	      density += cutseldata->maxsparsitybonus;
264  	      density = MAX(density, 0.0);
265  	
266  	      /* Get the normalised pseudo-cost and number of locks score */
267  	      if( root )
268  	         pscost = 0.0;
269  	      else
270  	         pscost = cutseldata->pscostweight * pscosts[i];
271  	      cutlock = cutseldata->locksweight * cutlocks[i];
272  	
273  	      /* Get the dynamism (good numerics) score */
274  	      maxcutval = SCIPgetRowMaxCoef(scip, cuts[i]);
275  	      mincutval = SCIPgetRowMinCoef(scip, cuts[i]);
276  	      mincutval = mincutval > 0.0 ? mincutval : 1.0;
277  	      dynamism = cutseldata->maxcoefratiobonus >= maxcutval / mincutval ? cutseldata->goodnumericsbonus : 0.0;
278  	
279  	      /* Get the dcd / eff / exp score */
280  	      if ( sol != NULL && root )
281  	      {
282  	         if ( SCIPisSumLE(scip, dcds[i], 0.0))
283  	            scaleddcd = 0.0;
284  	         else
285  	            scaleddcd = cutseldata->dircutoffdistweight * SQR(LOG1P(dcds[i]) / LOG1P(maxdcd)); /*lint !e666*/
286  	      }
287  	      else
288  	      {
289  	         scaleddcd = 0.0;
290  	      }
291  	
292  	      if ( SCIPisSumLE(scip, exps[i], 0.0))
293  	         scaledexp = 0.0;
294  	      else
295  	         scaledexp = cutseldata->expimprovweight * SQR(LOG1P(exps[i]) / LOG1P(maxexp)); /*lint !e666*/
296  	
297  	      if ( SCIPisSumLE(scip, effs[i], 0.0))
298  	      {
299  	         scaledeff = 0.0;
300  	      }
301  	      else
302  	      {
303  	         if ( sol != NULL && root )
304  	            scaledeff = cutseldata->efficacyweight * SQR(LOG1P(effs[i]) / LOG1P(maxeff)); /*lint !e666*/
305  	         else
306  	            scaledeff = (cutseldata->efficacyweight + cutseldata->dircutoffdistweight) * SQR(LOG1P(effs[i]) / LOG1P(maxeff)); /*lint !e666*/
307  	      }
308  	
309  	      /* Combine all scores and introduce some minor randomness */
310  	      score = scaledeff + scaleddcd + scaledexp + objparallelism + intsupport + density + dynamism + pscost + cutlock;
311  	
312  	      score += SCIPrandomGetReal(cutseldata->randnumgen, 0.0, 1e-6);
313  	
314  	      scores[i] = score;
315  	   }
316  	
317  	   SCIPfreeBufferArray(scip, &effs);
318  	   SCIPfreeBufferArray(scip, &dcds);
319  	   SCIPfreeBufferArray(scip, &exps);
320  	   SCIPfreeBufferArray(scip, &cutdensities);
321  	   SCIPfreeBufferArray(scip, &cutlocks);
322  	   SCIPfreeBufferArray(scip, &pscosts);
323  	
324  	   return SCIP_OKAY;
325  	
326  	}
327  	
328  	
329  	/** move the cut with the highest score to the first position in the array; there must be at least one cut */
330  	static
331  	void selectBestCut(
332  	   SCIP_ROW**            cuts,               /**< array with cuts to perform selection algorithm */
333  	   SCIP_Real*            scores,             /**< array with scores of cuts to perform selection algorithm */
334  	   int                   ncuts               /**< number of cuts in given array */
335  	)
336  	{
337  	   int bestpos;
338  	   SCIP_Real bestscore;
339  	
340  	   assert(ncuts > 0);
341  	   assert(cuts != NULL);
342  	   assert(scores != NULL);
343  	
344  	   bestscore = scores[0];
345  	   bestpos = 0;
346  	
347  	   for( int i = 1; i < ncuts; ++i )
348  	   {
349  	      if( scores[i] > bestscore )
350  	      {
351  	         bestpos = i;
352  	         bestscore = scores[i];
353  	      }
354  	   }
355  	
356  	   SCIPswapPointers((void**) &cuts[bestpos], (void**) &cuts[0]);
357  	   SCIPswapReals(&scores[bestpos], &scores[0]);
358  	}
359  	
360  	/** filters the given array of cuts to enforce a maximum parallelism constraint
361  	 *  w.r.t the given cut; moves filtered cuts to the end of the array and returns number of selected cuts */
362  	static
363  	int filterWithParallelism(
364  	   SCIP_ROW*             cut,                /**< cut to filter orthogonality with */
365  	   SCIP_ROW**            cuts,               /**< array with cuts to perform selection algorithm */
366  	   SCIP_Real*            scores,             /**< array with scores of cuts to perform selection algorithm */
367  	   int                   ncuts,              /**< number of cuts in given array */
368  	   SCIP_Real             maxparallel         /**< maximal parallelism for all cuts that are not good */
369  	)
370  	{
371  	
372  	   assert( cut != NULL );
373  	   assert( ncuts == 0 || cuts != NULL );
374  	   assert( ncuts == 0 || scores != NULL );
375  	
376  	   for( int i = ncuts - 1; i >= 0; --i )
377  	   {
378  	      SCIP_Real thisparallel;
379  	
380  	      thisparallel = SCIProwGetParallelism(cut, cuts[i], 'e');
381  	
382  	      if( thisparallel > maxparallel )
383  	      {
384  	         --ncuts;
385  	         SCIPswapPointers((void**) &cuts[i], (void**) &cuts[ncuts]);
386  	         SCIPswapReals(&scores[i], &scores[ncuts]);
387  	      }
388  	   }
389  	
390  	   return ncuts;
391  	}
392  	
393  	/** penalises any cut too parallel to cut by reducing the parallel cut's score. */
394  	static
395  	int penaliseWithParallelism(
396  	   SCIP*                 scip,               /**< SCIP data structure */
397  	   SCIP_ROW*             cut,                /**< cut to filter orthogonality with */
398  	   SCIP_ROW**            cuts,               /**< array with cuts to perform selection algorithm */
399  	   SCIP_Real*            scores,             /**< array with scores of cuts to perform selection algorithm */
400  	   int                   ncuts,              /**< number of cuts in given array */
401  	   SCIP_Real             maxparallel,        /**< maximal parallelism for all cuts that are not good */
402  	   SCIP_Real             paralpenalty        /**< penalty for weaker of two parallel cuts if penalising parallel cuts */
403  	)
404  	{
405  	
406  	   assert( cut != NULL );
407  	   assert( ncuts == 0 || cuts != NULL );
408  	   assert( ncuts == 0 || scores != NULL );
409  	
410  	   for( int i = ncuts - 1; i >= 0; --i )
411  	   {
412  	      SCIP_Real thisparallel;
413  	
414  	      thisparallel = SCIProwGetParallelism(cut, cuts[i], 'e');
415  	
416  	      /* Filter cuts that are absolutely parallel still. Otherwise penalise if closely parallel */
417  	      if( thisparallel > 1 - SCIPsumepsilon(scip) )
418  	      {
419  	         --ncuts;
420  	         SCIPswapPointers((void**) &cuts[i], (void**) &cuts[ncuts]);
421  	         SCIPswapReals(&scores[i], &scores[ncuts]);
422  	      }
423  	      else if( thisparallel > maxparallel )
424  	      {
425  	         scores[i] -= paralpenalty;
426  	      }
427  	   }
428  	
429  	   return ncuts;
430  	}
431  	
432  	/** filters the given array of cuts to enforce a maximum density constraint,
433  	 *  Moves filtered cuts to the end of the array and returns number of selected cuts */
434  	static
435  	int filterWithDensity(
436  	   SCIP*                 scip,               /**< SCIP data structure */
437  	   SCIP_ROW**            cuts,               /**< array with cuts to perform selection algorithm */
438  	   SCIP_Real             maxdensity,         /**< maximum density s.t. a cut is not filtered */
439  	   int                   ncuts               /**< number of cuts in given array */
440  	)
441  	{
442  	   SCIP_Real ncols;
443  	
444  	   assert( ncuts == 0 || cuts != NULL );
445  	
446  	   ncols = SCIPgetNLPCols(scip);
447  	
448  	   for( int i = ncuts - 1; i >= 0; --i )
449  	   {
450  	      SCIP_Real nvals;
451  	
452  	      nvals = SCIProwGetNNonz(cuts[i]);
453  	
454  	      if( maxdensity < nvals / ncols )
455  	      {
456  	         --ncuts;
457  	         SCIPswapPointers((void**) &cuts[i], (void**) &cuts[ncuts]);
458  	      }
459  	   }
460  	
461  	   return ncuts;
462  	}
463  	
464  	
465  	/*
466  	 * Callback methods of cut selector
467  	 */
468  	
469  	
470  	/** copy method for cut selector plugin (called when SCIP copies plugins) */
471  	static
472  	SCIP_DECL_CUTSELCOPY(cutselCopyEnsemble)
473  	{  /*lint --e{715}*/
474  	   assert(scip != NULL);
475  	   assert(cutsel != NULL);
476  	   assert(strcmp(SCIPcutselGetName(cutsel), CUTSEL_NAME) == 0);
477  	
478  	   /* call inclusion method of cut selector */
479  	   SCIP_CALL( SCIPincludeCutselEnsemble(scip) );
480  	
481  	   return SCIP_OKAY;
482  	}
483  	
484  	/** destructor of cut selector to free user data (called when SCIP is exiting) */
485  	/**! [SnippetCutselFreeEnsemble] */
486  	static
487  	SCIP_DECL_CUTSELFREE(cutselFreeEnsemble)
488  	{  /*lint --e{715}*/
489  	   SCIP_CUTSELDATA* cutseldata;
490  	
491  	   cutseldata = SCIPcutselGetData(cutsel);
492  	
493  	   SCIPfreeBlockMemory(scip, &cutseldata);
494  	
495  	   SCIPcutselSetData(cutsel, NULL);
496  	
497  	   return SCIP_OKAY;
498  	}
499  	/**! [SnippetCutselFreeEnsemble] */
500  	
501  	/** initialization method of cut selector (called after problem was transformed) */
502  	static
503  	SCIP_DECL_CUTSELINIT(cutselInitEnsemble)
504  	{  /*lint --e{715}*/
505  	   SCIP_CUTSELDATA* cutseldata;
506  	
507  	   cutseldata = SCIPcutselGetData(cutsel);
508  	   assert(cutseldata != NULL);
509  	
510  	   SCIP_CALL( SCIPcreateRandom(scip, &(cutseldata)->randnumgen, RANDSEED, TRUE) );
511  	
512  	   return SCIP_OKAY;
513  	}
514  	
515  	/** deinitialization method of cut selector (called before transformed problem is freed) */
516  	static
517  	SCIP_DECL_CUTSELEXIT(cutselExitEnsemble)
518  	{  /*lint --e{715}*/
519  	   SCIP_CUTSELDATA* cutseldata;
520  	
521  	   cutseldata = SCIPcutselGetData(cutsel);
522  	   assert(cutseldata != NULL);
523  	   assert(cutseldata->randnumgen != NULL);
524  	
525  	   SCIPfreeRandom(scip, &cutseldata->randnumgen);
526  	
527  	   return SCIP_OKAY;
528  	}
529  	
530  	/** cut selection method of cut selector */
531  	static
532  	SCIP_DECL_CUTSELSELECT(cutselSelectEnsemble)
533  	{  /*lint --e{715}*/
534  	   SCIP_CUTSELDATA* cutseldata;
535  	
536  	   assert(cutsel != NULL);
537  	   assert(result != NULL);
538  	
539  	   cutseldata = SCIPcutselGetData(cutsel);
540  	   assert(cutseldata != NULL);
541  	   
542  	   if( ncuts > cutseldata->maxcuts || SCIPgetNVars(scip) > cutseldata->maxnumvars )
543  	   {
544  	      *result = SCIP_DIDNOTFIND;
545  	      return SCIP_OKAY;
546  	   }
547  	   
548  	   *result = SCIP_SUCCESS;
549  	
550  	   SCIP_CALL( SCIPselectCutsEnsemble(scip, cuts, forcedcuts, cutseldata, root, ncuts, nforcedcuts,
551  	      maxnselectedcuts, nselectedcuts) );
552  	
553  	   return SCIP_OKAY;
554  	}
555  	
556  	
557  	/*
558  	 * cut selector specific interface methods
559  	 */
560  	
561  	/** creates the ensemble cut selector and includes it in SCIP */
562  	SCIP_RETCODE SCIPincludeCutselEnsemble(
563  	   SCIP*                 scip                /**< SCIP data structure */
564  	)
565  	{
566  	   SCIP_CUTSELDATA* cutseldata;
567  	   SCIP_CUTSEL* cutsel;
568  	
569  	   /* create ensemble cut selector data */
570  	   SCIP_CALL( SCIPallocBlockMemory(scip, &cutseldata) );
571  	   BMSclearMemory(cutseldata);
572  	
573  	   SCIP_CALL( SCIPincludeCutselBasic(scip, &cutsel, CUTSEL_NAME, CUTSEL_DESC, CUTSEL_PRIORITY, cutselSelectEnsemble,
574  	                                     cutseldata) );
575  	
576  	   assert(cutsel != NULL);
577  	
578  	   /* set non fundamental callbacks via setter functions */
579  	   SCIP_CALL( SCIPsetCutselCopy(scip, cutsel, cutselCopyEnsemble) );
580  	
581  	   SCIP_CALL( SCIPsetCutselFree(scip, cutsel, cutselFreeEnsemble) );
582  	   SCIP_CALL( SCIPsetCutselInit(scip, cutsel, cutselInitEnsemble) );
583  	   SCIP_CALL( SCIPsetCutselExit(scip, cutsel, cutselExitEnsemble) );
584  	
585  	   /* add ensemble cut selector parameters */
586  	   SCIP_CALL( SCIPaddRealParam(scip,
587  	                               "cutselection/" CUTSEL_NAME "/efficacyweight",
588  	      "weight of normed-efficacy in cut score calculation",
589  	      &cutseldata->efficacyweight, FALSE, DEFAULT_EFFICACYWEIGHT, 0.0, SCIP_INVALID/10.0, NULL, NULL) );
590  	
591  	   SCIP_CALL( SCIPaddRealParam(scip,
592  	                               "cutselection/" CUTSEL_NAME "/dircutoffdistweight",
593  	      "weight of normed-directed cutoff distance in cut score calculation",
594  	      &cutseldata->dircutoffdistweight, FALSE, DEFAULT_DIRCUTOFFDISTWEIGHT, 0.0, SCIP_INVALID/10.0, NULL, NULL) );
595  	
596  	   SCIP_CALL( SCIPaddRealParam(scip,
597  	                               "cutselection/" CUTSEL_NAME "/objparalweight",
598  	      "weight of objective parallelism in cut score calculation",
599  	      &cutseldata->objparalweight, FALSE, DEFAULT_OBJPARALWEIGHT, 0.0, SCIP_INVALID/10.0, NULL, NULL) );
600  	
601  	   SCIP_CALL( SCIPaddRealParam(scip,
602  	                               "cutselection/" CUTSEL_NAME "/intsupportweight",
603  	      "weight of integral support in cut score calculation",
604  	      &cutseldata->intsupportweight, FALSE, DEFAULT_INTSUPPORTWEIGHT, 0.0, SCIP_INVALID/10.0, NULL, NULL) );
605  	
606  	   SCIP_CALL( SCIPaddRealParam(scip,
607  	                               "cutselection/" CUTSEL_NAME "/expimprovweight",
608  	      "weight of normed-expected obj improvement in cut score calculation",
609  	      &cutseldata->expimprovweight, FALSE, DEFAULT_EXPIMPROVWEIGHT, 0.0, SCIP_INVALID/10.0, NULL, NULL) );
610  	
611  	   SCIP_CALL( SCIPaddRealParam(scip,
612  	                               "cutselection/" CUTSEL_NAME "/minscore",
613  	      "minimum score s.t. a cut can be added",
614  	      &cutseldata->minscore, FALSE, DEFAULT_MINSCORE, -SCIP_INVALID/10.0, SCIP_INVALID/10.0, NULL, NULL) );
615  	
616  	   SCIP_CALL( SCIPaddRealParam(scip,
617  	                               "cutselection/" CUTSEL_NAME "/pscostweight",
618  	      "weight of normed-pseudo-costs in cut score calculation",
619  	      &cutseldata->pscostweight, FALSE, DEFAULT_PSCOSTWEIGHT, 0.0, SCIP_INVALID/10.0, NULL, NULL) );
620  	
621  	   SCIP_CALL( SCIPaddRealParam(scip,
622  	                               "cutselection/" CUTSEL_NAME "/locksweight",
623  	      "weight of normed-num-locks in cut score calculation",
624  	      &cutseldata->locksweight, FALSE, DEFAULT_NLOCKSWEIGHT, 0.0, SCIP_INVALID/10.0, NULL, NULL) );
625  	
626  	   SCIP_CALL( SCIPaddRealParam(scip,
627  	                               "cutselection/" CUTSEL_NAME "/maxsparsitybonus",
628  	      "weight of maximum sparsity reward in cut score calculation",
629  	      &cutseldata->maxsparsitybonus, FALSE, DEFAULT_MAXSPARSITYBONUS, 0.0, SCIP_INVALID/10.0, NULL, NULL) );
630  	
631  	   SCIP_CALL( SCIPaddRealParam(scip,
632  	                               "cutselection/" CUTSEL_NAME "/goodnumericsbonus",
633  	      "weight of good numerics bonus (ratio of coefficients) in cut score calculation",
634  	      &cutseldata->goodnumericsbonus, FALSE, DEFAULT_GOODNUMERICBONUS, 0.0, SCIP_INVALID/10.0, NULL, NULL) );
635  	
636  	   SCIP_CALL( SCIPaddRealParam(scip,
637  	                               "cutselection/" CUTSEL_NAME "/endsparsitybonus",
638  	      "max sparsity value for which a bonus is applied in cut score calculation",
639  	      &cutseldata->endsparsitybonus, FALSE, DEFAULT_SPARSITYENDBONUS, 0.0, SCIP_INVALID/10.0, NULL, NULL) );
640  	
641  	   SCIP_CALL( SCIPaddRealParam(scip,
642  	                               "cutselection/" CUTSEL_NAME "/maxparal",
643  	      "threshold for when two cuts are considered parallel to each other",
644  	      &cutseldata->maxparal, FALSE, DEFAULT_MAXPARAL, 0.0, 1.0, NULL, NULL) );
645  	
646  	   SCIP_CALL( SCIPaddRealParam(scip,
647  	                               "cutselection/" CUTSEL_NAME "/paralpenalty",
648  	      "penalty for weaker of two parallel cuts if penalising parallel cuts",
649  	      &cutseldata->paralpenalty, TRUE, DEFAULT_PARALPENALTY, 0.0, SCIP_INVALID/10.0, NULL, NULL) );
650  	
651  	   SCIP_CALL( SCIPaddRealParam(scip,
652  	                               "cutselection/" CUTSEL_NAME "/maxcutdensity",
653  	      "max allowed cut density if filtering dense cuts",
654  	      &cutseldata->maxcutdensity, TRUE, DEFAULT_MAXCUTDENSITY, 0.0, 1.0, NULL, NULL) );
655  	
656  	   SCIP_CALL( SCIPaddRealParam(scip,
657  	                               "cutselection/" CUTSEL_NAME "/maxnonzerorootround",
658  	      "max non-zeros per round applied cuts (root). multiple num LP cols.",
659  	      &cutseldata->maxnonzerorootround, FALSE, DEFAULT_MAXNONZEROROOTROUND, 0.0, SCIP_INVALID/10.0, NULL, NULL) );
660  	
661  	   SCIP_CALL( SCIPaddRealParam(scip,
662  	                               "cutselection/" CUTSEL_NAME "/maxnonzerotreeround",
663  	      "max non-zeros per round applied cuts (tree). multiple num LP cols.",
664  	      &cutseldata->maxnonzerotreeround, FALSE, DEFAULT_MAXNONZEROTREEROUND, 0.0, SCIP_INVALID/10.0, NULL, NULL) );
665  	
666  	   SCIP_CALL( SCIPaddBoolParam(scip,
667  	                               "cutselection/" CUTSEL_NAME "/filterparalcuts",
668  	      "should cuts be filtered so no two parallel cuts are added",
669  	      &cutseldata->filterparalcuts, FALSE, DEFAULT_FILTERPARALCUTS, NULL, NULL) );
670  	
671  	   SCIP_CALL( SCIPaddBoolParam(scip,
672  	                               "cutselection/" CUTSEL_NAME "/penaliseparalcuts",
673  	      "should two parallel cuts be penalised instead of outright filtered",
674  	      &cutseldata->penaliseparalcuts, TRUE, DEFAULT_PENALISEPARALCUTS, NULL, NULL) );
675  	
676  	   SCIP_CALL( SCIPaddBoolParam(scip,
677  	                               "cutselection/" CUTSEL_NAME "/filterdensecuts",
678  	      "should cuts over a given density threshold be filtered",
679  	      &cutseldata->filterdensecuts, TRUE, DEFAULT_FILTERDENSECUTS, NULL, NULL) );
680  	
681  	   SCIP_CALL( SCIPaddBoolParam(scip,
682  	                               "cutselection/" CUTSEL_NAME "/penaliselocks",
683  	      "should the number of locks be penalised instead of rewarded",
684  	      &cutseldata->penaliselocks, TRUE, DEFAULT_PENALISELOCKS, NULL, NULL) );
685  	
686  	   SCIP_CALL( SCIPaddBoolParam(scip,
687  	                               "cutselection/" CUTSEL_NAME "/penaliseobjparal",
688  	      "should objective parallelism be penalised instead of rewarded",
689  	      &cutseldata->penaliseobjparal, TRUE, DEFAULT_PENALISEOBJPARAL, NULL, NULL) );
690  	
691  	   SCIP_CALL( SCIPaddIntParam(scip,
692  	                               "cutselection/" CUTSEL_NAME "/maxcoefratiobonus",
693  	      "max coefficient ratio for which numeric bonus is applied.",
694  	      &cutseldata->maxcoefratiobonus, TRUE, DEFAULT_MAXCOEFRATIOBONUS, 1, 1000000, NULL, NULL) );
695  	   
696  	   SCIP_CALL( SCIPaddIntParam(scip,
697  	                              "cutselection/" CUTSEL_NAME "/maxcuts",
698  	      "max number of cuts such that cut selector is applied.",
699  	      &cutseldata->maxcuts, TRUE, DEFAULT_MAXCUTS, 1, 1000000, NULL, NULL) );
700  	   
701  	   SCIP_CALL( SCIPaddIntParam(scip,
702  	                              "cutselection/" CUTSEL_NAME "/maxnumvars",
703  	      "max number of variables such that cut selector is applied.",
704  	      &cutseldata->maxnumvars, TRUE, DEFAULT_MAXNUMVARS, 1, 1000000, NULL, NULL) );
705  	
706  	   return SCIP_OKAY;
707  	
708  	}
709  	
710  	/** perform a cut selection algorithm for the given array of cuts
711  	 *
712  	 *  This is the selection method of the ensemble cut selector. It uses a weighted sum of normalised efficacy,
713  	 *  normalised directed cutoff distance, normalised expected improvements, objective parallelism,
714  	 *  integer support, sparsity, dynamism, pseudo-costs, and variable locks.
715  	 *  In addition to the weighted sum score, there are optionally parallelism-based filtering and penalties,
716  	 *  and density filtering.
717  	 *  There are also additional budget constraints on the number of cuts that should be added.
718  	 *  The input cuts array gets re-sorted such that the selected cuts come first and the remaining ones are the end.
719  	 */
720  	SCIP_RETCODE SCIPselectCutsEnsemble(
721  	   SCIP*                 scip,               /**< SCIP data structure */
722  	   SCIP_ROW**            cuts,               /**< array with cuts to perform selection algorithm */
723  	   SCIP_ROW**            forcedcuts,         /**< array with forced cuts */
724  	   SCIP_CUTSELDATA*      cutseldata,         /**< cut selector data */
725  	   SCIP_Bool             root,               /**< whether we are at the root node or not */
726  	   int                   ncuts,              /**< number of cuts in cuts array */
727  	   int                   nforcedcuts,        /**< number of forced cuts */
728  	   int                   maxselectedcuts,    /**< maximal number of cuts from cuts array to select */
729  	   int*                  nselectedcuts       /**< pointer to return number of selected cuts from cuts array */
730  	)
731  	{
732  	   SCIP_Real* scores;
733  	   SCIP_Real* origscoresptr;
734  	   SCIP_Real nonzerobudget;
735  	   SCIP_Real budgettaken = 0.0;
736  	   SCIP_Real ncols;
737  	
738  	   assert(cuts != NULL && ncuts > 0);
739  	   assert(forcedcuts != NULL || nforcedcuts == 0);
740  	   assert(nselectedcuts != NULL);
741  	
742  	   *nselectedcuts = 0;
743  	   ncols = SCIPgetNLPCols(scip);
744  	
745  	   /* filter dense cuts */
746  	   if( cutseldata->filterdensecuts )
747  	   {
748  	      ncuts = filterWithDensity(scip, cuts, cutseldata->maxcutdensity, ncuts);
749  	      if( ncuts == 0 )
750  	         return SCIP_OKAY;
751  	   }
752  	
753  	   /* Initialise the score array */
754  	   SCIP_CALL( SCIPallocBufferArray(scip, &scores, ncuts) );
755  	   origscoresptr = scores;
756  	
757  	   /* compute scores of cuts */
758  	   SCIP_CALL( scoring(scip, cuts, cutseldata, scores, root, ncuts) );
759  	
760  	   /* perform cut selection algorithm for the cuts */
761  	
762  	   /* forced cuts are going to be selected so use them to filter cuts */
763  	   for( int i = 0; i < nforcedcuts && ncuts > 0; ++i )
764  	   {
765  	      if( cutseldata->filterparalcuts )
766  	         ncuts = filterWithParallelism(forcedcuts[i], cuts, scores, ncuts, cutseldata->maxparal);
767  	      else if( cutseldata->penaliseparalcuts )
768  	         ncuts = penaliseWithParallelism(scip, forcedcuts[i], cuts, scores, ncuts, cutseldata->maxparal, cutseldata->paralpenalty);
769  	   }
770  	
771  	   /* Get the budget depending on if we are the root or not */
772  	   nonzerobudget = root ? cutseldata->maxnonzerorootround : cutseldata->maxnonzerotreeround;
773  	
774  	   /* now greedily select the remaining cuts */
775  	   while( ncuts > 0 )
776  	   {
777  	      SCIP_ROW* selectedcut;
778  	
779  	      selectBestCut(cuts, scores, ncuts);
780  	      selectedcut = cuts[0];
781  	
782  	      /* if the best cut of the remaining cuts is considered bad, we discard it and all remaining cuts */
783  	      if( scores[0] < cutseldata->minscore )
784  	         break;
785  	
786  	      ++(*nselectedcuts);
787  	
788  	      /* if the maximal number of cuts was selected, we can stop here */
789  	      if( *nselectedcuts == maxselectedcuts )
790  	         break;
791  	
792  	      /* increase the non-zero budget counter of added cuts */
793  	      budgettaken += SCIProwGetNNonz(cuts[0]) / ncols;
794  	
795  	      /* move the pointers to the next position and filter the remaining cuts to enforce the maximum parallelism constraint */
796  	      ++cuts;
797  	      ++scores;
798  	      --ncuts;
799  	
800  	      if( cutseldata->filterparalcuts && ncuts > 0)
801  	         ncuts = filterWithParallelism(selectedcut, cuts, scores, ncuts, cutseldata->maxparal);
802  	      else if( cutseldata->penaliseparalcuts && ncuts > 0 )
803  	         ncuts = penaliseWithParallelism(scip, selectedcut, cuts, scores, ncuts, cutseldata->maxparal, cutseldata->paralpenalty);
804  	
805  	      /* Filter out all remaining cuts that would go over the non-zero budget threshold */
806  	      if( nonzerobudget - budgettaken < 1 && ncuts > 0 )
807  	         ncuts = filterWithDensity(scip, cuts, nonzerobudget - budgettaken, ncuts);
808  	
809  	   }
810  	
811  	   SCIPfreeBufferArray(scip, &origscoresptr);
812  	
813  	   return SCIP_OKAY;
814  	}
815