Scippy

SCIP

Solving Constraint Integer Programs

prop_symmetry.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (c) 2002-2024 Zuse Institute Berlin (ZIB) */
7 /* */
8 /* Licensed under the Apache License, Version 2.0 (the "License"); */
9 /* you may not use this file except in compliance with the License. */
10 /* You may obtain a copy of the License at */
11 /* */
12 /* http://www.apache.org/licenses/LICENSE-2.0 */
13 /* */
14 /* Unless required by applicable law or agreed to in writing, software */
15 /* distributed under the License is distributed on an "AS IS" BASIS, */
16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17 /* See the License for the specific language governing permissions and */
18 /* limitations under the License. */
19 /* */
20 /* You should have received a copy of the Apache-2.0 license */
21 /* along with SCIP; see the file LICENSE. If not visit scipopt.org. */
22 /* */
23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24 
25 /**@file prop_symmetry.c
26  * @ingroup DEFPLUGINS_PROP
27  * @brief propagator for handling symmetries
28  * @author Marc Pfetsch
29  * @author Thomas Rehn
30  * @author Christopher Hojny
31  * @author Fabian Wegscheider
32  * @author Jasper van Doornmalen
33  *
34  * This propagator combines the following symmetry handling functionalities:
35  * - It allows to compute symmetries of the problem and to store this information in adequate form. The symmetry
36  * information can be accessed through external functions.
37  * - It implements various methods to handle the symmetries:
38  * - orbital reduction, which generalizes orbital fixing. See symmetry_orbital.c
39  * - (dynamic) orbitopal reduction, which generalizes (dynamic) orbital fixing. See symmetry_orbitopal.c
40  * - static orbitopal fixing (for binary variable domains) for full orbitopes. See cons_orbitope.c
41  * - static orbitopal fixing (for binary variable domains) for packing-partitioning orbitopes. See cons_orbitope.c
42  * - (dynamic) lexicographic reduction. See symmetry_lexred.c
43  * - static lexicographic fixing for binary variable domains (i.e., symresack propagation). See cons_symresack.c
44  * - static lexicographic fixing for binary variable domains on involutions (i.e., orbisacks). See cons_orbisack.c
45  * - Symmetry breaking inequalities based on the Schreier-Sims Table (i.e., SST cuts).
46  * - Strong and weak symmetry breaking inequalities.
47  *
48  *
49  * @section SYMCOMP Symmetry Computation
50  *
51  * The generic functionality of the compute_symmetry.h interface is used.
52  * We do not copy symmetry information, since it is not clear how this information transfers. Moreover, copying
53  * symmetry might inhibit heuristics. But note that solving a sub-SCIP might then happen without symmetry information!
54  *
55  *
56  * @section SYMBREAK Symmetry handling by the (unified) symmetry handling constraints
57  *
58  * Many common methods are captured by a framework that dynamifies symmetry handling constraints. The ideas are
59  * described in@n
60  * J. van Doornmalen, C. Hojny, "A Unified Framework for Symmetry Handling", preprint, 2023,
61  * https://doi.org/10.48550/arXiv.2211.01295.
62  *
63  * This paper shows that various symmetry handling methods are compatible under certain conditions, and provides
64  * generalizations to common symmetry handling constraints from binary variable domains to arbitrary variable domains.
65  * This includes symresack propagation, orbitopal fixing, and orbital fixing, that are generalized to
66  * lexicographic reduction, orbitopal reduction and orbital reduction, respectively. For a description and
67  * implementation, see symmetry_lexred.c, symmetry_orbitopal.c and symmetry_orbital.c, respectively.
68  * The static counterparts on binary variable domains are cons_symresack.c and cons_orbisack.c for lexicographic
69  * reduction (cf. symresack propagation), and cons_orbitope.c and cons_orbisack.c for orbitopal reduction
70  * (cf. orbitopal fixing). We refer to the description of tryAddSymmetryHandlingMethods for the order in which these
71  * methods are applied.
72  *
73  * @section SST Cuts derived from the Schreier Sims table
74  *
75  * SST cuts have been introduced by@n
76  * D. Salvagnin: Symmetry Breaking Inequalities from the Schreier-Sims table. CPAIOR 2018 Proceedings, 521-529, 2018.
77  *
78  * These inequalities are computed as follows. Throughout these procedure a set of so-called leaders is maintained.
79  * Initially the set of leaders is empty. In a first step, select a variable \f$x_i\f$ and compute its orbit w.r.t.
80  * the symmetry group of the mixed-integer program. For each variable \f$x_j\f$ in the orbit of \f$x_i\f$, the
81  * inequality \f$x_i \geq x_j\f$ is a valid symmetry handling inequality, which can be added to the mixed-integer
82  * program. We call \f$x_i\f$ the leader of this inequality. Add the leader \f$x_i\f$ to the set of leaders and
83  * compute the pointwise stabilizer of the leader set. In the next step, select a new variable, compute its orbit
84  * w.r.t. the stabilizer group of the leaders, add the inequalities based on this orbit, and add the new leader
85  * to the set of leaders. This procedure is iterated until the pointwise stabilizer group of the leaders has become
86  * trivial.
87  *
88  * @todo Possibly turn off propagator in subtrees.
89  * @todo Check application of conflict resolution.
90  * @todo Check whether one should switch the role of 0 and 1
91  * @todo Implement stabilizer computation?
92  * @todo Implement isomorphism pruning?
93  * @todo Implement particular preprocessing rules
94  * @todo Separate permuted cuts (first experiments not successful)
95  * @todo Allow the computation of local symmetries
96  * @todo Order rows of orbitopes (in particular packing/partitioning) w.r.t. cliques in conflict graph.
97  * @todo A dynamic variant for packing-partitioning orbitopal structures
98  * @todo A dynamic variant for suborbitopes
99  */
100 /* #define SCIP_OUTPUT */
101 /* #define SCIP_OUTPUT_COMPONENT */
102 /* #define SCIP_DISPLAY_SYM_CHECK */
103 
104 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
105 
106 #include <scip/cons_linear.h>
107 #include <scip/cons_knapsack.h>
108 #include <scip/cons_varbound.h>
109 #include <scip/cons_setppc.h>
110 #include <scip/cons_and.h>
111 #include <scip/cons_logicor.h>
112 #include <scip/cons_or.h>
113 #include <scip/cons_orbitope.h>
114 #include <scip/cons_symresack.h>
115 #include <scip/cons_xor.h>
116 #include <scip/cons_linking.h>
118 #include <scip/cons_indicator.h>
119 #include <scip/cons_nonlinear.h>
120 #include <scip/cons_sos1.h>
121 #include <scip/cons_sos2.h>
122 #include <scip/expr_pow.h>
123 #include <scip/expr_product.h>
124 #include <scip/pub_expr.h>
125 #include <scip/misc.h>
127 
128 #include <scip/prop_symmetry.h>
130 #include <scip/event_shadowtree.h>
131 #include <scip/symmetry.h>
132 #include <scip/symmetry_graph.h>
133 #include <scip/symmetry_orbitopal.h>
134 #include <scip/symmetry_orbital.h>
135 #include <scip/symmetry_lexred.h>
136 
137 #include <math.h>
138 #include <string.h>
139 
140 /* propagator properties */
141 #define PROP_NAME "symmetry"
142 #define PROP_DESC "propagator for handling symmetry"
143 #define PROP_TIMING SCIP_PROPTIMING_BEFORELP /**< propagation timing mask */
144 #define PROP_PRIORITY -1000000 /**< propagator priority */
145 #define PROP_FREQ 1 /**< propagator frequency */
146 #define PROP_DELAY FALSE /**< should propagation method be delayed, if other propagators found reductions? */
147 
148 #define PROP_PRESOL_PRIORITY -10000000 /**< priority of the presolving method (>= 0: before, < 0: after constraint handlers) */
149 #define PROP_PRESOLTIMING SCIP_PRESOLTIMING_EXHAUSTIVE /* timing of the presolving method (fast, medium, or exhaustive) */
150 #define PROP_PRESOL_MAXROUNDS -1 /**< maximal number of presolving rounds the presolver participates in (-1: no limit) */
151 
152 
153 /* default parameter values for symmetry computation */
154 #define DEFAULT_MAXGENERATORS 1500 /**< limit on the number of generators that should be produced within symmetry detection (0 = no limit) */
155 #define DEFAULT_CHECKSYMMETRIES FALSE /**< Should all symmetries be checked after computation? */
156 #define DEFAULT_DISPLAYNORBITVARS FALSE /**< Should the number of variables affected by some symmetry be displayed? */
157 #define DEFAULT_USECOLUMNSPARSITY FALSE /**< Should the number of conss a variable is contained in be exploited in symmetry detection? */
158 #define DEFAULT_DOUBLEEQUATIONS FALSE /**< Double equations to positive/negative version? */
159 #define DEFAULT_COMPRESSSYMMETRIES TRUE /**< Should non-affected variables be removed from permutation to save memory? */
160 #define DEFAULT_COMPRESSTHRESHOLD 0.5 /**< Compression is used if percentage of moved vars is at most the threshold. */
161 #define DEFAULT_SYMFIXNONBINARYVARS FALSE /**< Disabled parameter */
162 #define DEFAULT_ENFORCECOMPUTESYMMETRY FALSE /**< always compute symmetries, even if they cannot be handled */
163 #define DEFAULT_SYMTYPE (int) SYM_SYMTYPE_PERM /**< type of symmetries to be computed */
164 
165 /* default parameters for linear symmetry constraints */
166 #define DEFAULT_CONSSADDLP TRUE /**< Should the symmetry breaking constraints be added to the LP? */
167 #define DEFAULT_ADDSYMRESACKS TRUE /**< Add inequalities for symresacks for each generator? */
168 #define DEFAULT_DETECTDOUBLELEX TRUE /**< Should we check whether the components of the symmetry group can be handled by double lex matrices? */
169 #define DEFAULT_DETECTORBITOPES TRUE /**< Should we check whether the components of the symmetry group can be handled by orbitopes? */
170 #define DEFAULT_DETECTSUBGROUPS TRUE /**< Should we try to detect orbitopes in subgroups of the symmetry group? */
171 #define DEFAULT_ADDWEAKSBCS TRUE /**< Should we add weak SBCs for enclosing orbit of symmetric subgroups? */
172 #define DEFAULT_ADDSTRONGSBCS FALSE /**< Should we add strong SBCs for enclosing orbit of symmetric subgroups if orbitopes are not used? */
173 #define DEFAULT_ADDCONSSTIMING 2 /**< timing of adding constraints (0 = before presolving, 1 = during presolving, 2 = after presolving) */
174 #define DEFAULT_MAXNCONSSSUBGROUP 500000 /**< Maximum number of constraints up to which subgroup structures are detected */
175 #define DEFAULT_USEDYNAMICPROP TRUE /**< whether dynamic propagation should be used for full orbitopes */
176 #define DEFAULT_PREFERLESSROWS TRUE /**< Shall orbitopes with less rows be preferred in detection? */
177 
178 /* default parameters for symmetry computation */
179 #define DEFAULT_SYMCOMPTIMING 2 /**< timing of symmetry computation (0 = before presolving, 1 = during presolving, 2 = at first call) */
180 #define DEFAULT_PERFORMPRESOLVING 0 /**< Run orbital fixing during presolving? (disabled parameter) */
181 #define DEFAULT_RECOMPUTERESTART 0 /**< Recompute symmetries after a restart has occurred? (0 = never) */
182 
183 /* default parameters for Schreier Sims constraints */
184 #define DEFAULT_SSTTIEBREAKRULE 1 /**< index of tie break rule for selecting orbit for Schreier Sims constraints? */
185 #define DEFAULT_SSTLEADERRULE 0 /**< index of rule for selecting leader variables for Schreier Sims constraints? */
186 #define DEFAULT_SSTLEADERVARTYPE 14 /**< bitset encoding which variable types can be leaders (1: binary; 2: integer; 4: impl. int; 8: continuous);
187  * if multiple types are allowed, take the one with most affected vars */
188 #define DEFAULT_ADDCONFLICTCUTS TRUE /**< Should Schreier Sims constraints be added if we use a conflict based rule? */
189 #define DEFAULT_SSTADDCUTS TRUE /**< Should Schreier Sims constraints be added? */
190 #define DEFAULT_SSTMIXEDCOMPONENTS TRUE /**< Should Schreier Sims constraints be added if a symmetry component contains variables of different types? */
192 /* output table properties */
193 #define TABLE_NAME_SYMMETRY "symmetry"
194 #define TABLE_DESC_SYMMETRY "symmetry handling statistics"
195 #define TABLE_POSITION_SYMMETRY 7001 /**< the position of the statistics table */
196 #define TABLE_EARLIEST_SYMMETRY SCIP_STAGE_SOLVING /**< output of the statistics table is only printed from this stage onwards */
198 
199 /* other defines */
200 #define MAXGENNUMERATOR 64000000 /**< determine maximal number of generators by dividing this number by the number of variables */
201 #define COMPRESSNVARSLB 25000 /**< lower bound on the number of variables above which compression could be performed */
203 /* macros for getting activeness of symmetry handling methods */
204 #define ISSYMRETOPESACTIVE(x) (((unsigned) x & SYM_HANDLETYPE_SYMBREAK) != 0)
205 #define ISORBITALREDUCTIONACTIVE(x) (((unsigned) x & SYM_HANDLETYPE_ORBITALREDUCTION) != 0)
206 #define ISSSTACTIVE(x) (((unsigned) x & SYM_HANDLETYPE_SST) != 0)
208 #define ISSSTBINACTIVE(x) (((unsigned) x & SCIP_SSTTYPE_BINARY) != 0)
209 #define ISSSTINTACTIVE(x) (((unsigned) x & SCIP_SSTTYPE_INTEGER) != 0)
210 #define ISSSTIMPLINTACTIVE(x) (((unsigned) x & SCIP_SSTTYPE_IMPLINT) != 0)
211 #define ISSSTCONTACTIVE(x) (((unsigned) x & SCIP_SSTTYPE_CONTINUOUS) != 0)
213 /* enable symmetry statistics */
214 #define SYMMETRY_STATISTICS 1
216 /** propagator data */
217 struct SCIP_PropData
218 {
219  /* symmetry group information */
220  int npermvars; /**< number of variables for permutations */
221  int nbinpermvars; /**< number of binary variables for permutations */
222  SCIP_VAR** permvars; /**< variables on which permutations act */
223  int nperms; /**< number of permutations */
224  int nmaxperms; /**< maximal number of permutations (needed for freeing storage) */
225  int** perms; /**< pointer to store permutation generators as (nperms x npermvars) matrix */
226  int** permstrans; /**< pointer to store transposed permutation generators as (npermvars x nperms) matrix */
227  SCIP_HASHMAP* permvarmap; /**< map of variables to indices in permvars array */
228  int nmovedpermvars; /**< number of variables moved by any permutation */
229  int nmovedbinpermvars; /**< number of binary variables moved by any permutation */
230  int nmovedintpermvars; /**< number of integer variables moved by any permutation */
231  int nmovedimplintpermvars; /**< number of implicitly integer variables moved by any permutation */
232  int nmovedcontpermvars; /**< number of continuous variables moved by any permutation */
233  SCIP_HASHMAP* customsymopnodetypes; /**< types of operator nodes introduced
234  * by a user for symmetry detection */
235  int nopnodetypes; /**< current number of operator node types used for symmetry detection */
236  SCIP_Real* permvardomaincenter; /**< center of variable domains (needed for signed permutations) */
237  int symtype; /**< type of symmetries to be computed */
238 
239  /* components of symmetry group */
240  int ncomponents; /**< number of components of symmetry group */
241  int ncompblocked; /**< number of components that have been blocked */
242  int* components; /**< array containing the indices of permutations sorted by components */
243  int* componentbegins; /**< array containing in i-th position the first position of
244  * component i in components array */
245  int* vartocomponent; /**< array containing for each permvar the index of the component it is
246  * contained in (-1 if not affected) */
247  unsigned* componentblocked; /**< array to store which symmetry methods have been applied to a component using
248  * the same bitset as for misc/usesymmetry */
249  SCIP_Bool* componenthassignedperm; /**< array to indicate whether a component has a signed permutation */
250 
251  /* further symmetry information */
252  int nmovedvars; /**< number of variables moved by some permutation */
253  SCIP_Real log10groupsize; /**< log10 of size of symmetry group */
254  SCIP_Bool binvaraffected; /**< whether binary variables are affected by some symmetry */
255 
256  /* for symmetry computation */
257  int maxgenerators; /**< limit on the number of generators that should be produced within symmetry detection (0 = no limit) */
258  SCIP_Bool checksymmetries; /**< Should all symmetries be checked after computation? */
259  SCIP_Bool displaynorbitvars; /**< Whether the number of variables in non-trivial orbits shall be computed */
260  SCIP_Bool compresssymmetries; /**< Should non-affected variables be removed from permutation to save memory? */
261  SCIP_Real compressthreshold; /**< Compression is used if percentage of moved vars is at most the threshold. */
262  SCIP_Bool compressed; /**< Whether symmetry data has been compressed */
263  SCIP_Bool computedsymmetry; /**< Have we already tried to compute symmetries? */
264  int usesymmetry; /**< encoding of active symmetry handling methods (for debugging) */
265  SCIP_Bool usecolumnsparsity; /**< Should the number of conss a variable is contained in be exploited in symmetry detection? */
266  SCIP_Bool doubleequations; /**< Double equations to positive/negative version? */
267  SCIP_Bool enforcecomputesymmetry; /**< always compute symmetries, even if they cannot be handled */
268 
269  /* for symmetry constraints */
270  SCIP_Bool triedaddsymmethods; /**< whether we already tried to add symmetry handling methods */
271  SCIP_Bool conssaddlp; /**< Should the symmetry breaking constraints be added to the LP? */
272  SCIP_Bool addsymresacks; /**< Add symresack constraints for each generator? */
273  int addconsstiming; /**< timing of adding constraints (0 = before presolving, 1 = during presolving, 2 = after presolving) */
274  SCIP_CONS** genorbconss; /**< list of generated orbitope/orbisack/symresack constraints */
275  SCIP_CONS** genlinconss; /**< list of generated linear constraints */
276  int ngenorbconss; /**< number of generated orbitope/orbisack/symresack constraints */
277  int genorbconsssize; /**< size of generated orbitope/orbisack/symresack constraints array */
278  int ngenlinconss; /**< number of generated linear constraints */
279  int genlinconsssize; /**< size of linear constraints array */
280  int nsymresacks; /**< number of symresack constraints */
281  SCIP_Bool detectdoublelex; /**< Should we check whether the components of the symmetry group can be handled by double lex matrices? */
282  SCIP_Bool detectorbitopes; /**< Should we check whether the components of the symmetry group can be handled by orbitopes? */
283  SCIP_Bool detectsubgroups; /**< Should we try to detect orbitopes in subgroups of the symmetry group? */
284  SCIP_Bool addweaksbcs; /**< Should we add weak SBCs for enclosing orbit of symmetric subgroups? */
285  SCIP_Bool addstrongsbcs; /**< Should we add strong SBCs for enclosing orbit of symmetric subgroups if orbitopes are not used? */
286  int norbitopes; /**< number of orbitope constraints */
287  SCIP_Bool* isnonlinvar; /**< array indicating whether variables appear non-linearly */
288  SCIP_CONSHDLR* conshdlr_nonlinear; /**< nonlinear constraint handler */
289  int maxnconsssubgroup; /**< maximum number of constraints up to which subgroup structures are detected */
290  SCIP_Bool usedynamicprop; /**< whether dynamic propagation should be used for full orbitopes */
291  SCIP_Bool preferlessrows; /**< Shall orbitopes with less rows be preferred in detection? */
292 
293  /* data necessary for symmetry computation order */
294  int recomputerestart; /**< Recompute symmetries after a restart has occurred? (0 = never, 1 = always, 2 = if symmetry reduction found) */
295  int symcomptiming; /**< timing for computation symmetries (0 = before presolving, 1 = during presolving, 2 = at first call) */
296  int lastrestart; /**< last restart for which symmetries have been computed */
297  SCIP_Bool symfoundreduction; /**< whether symmetry handling propagation has found a reduction since the last time computing symmetries */
298 
299  /* data necessary for Schreier Sims constraints */
300  SCIP_CONS** sstconss; /**< list of generated schreier sims conss */
301  int nsstconss; /**< number of generated schreier sims conss */
302  int maxnsstconss; /**< maximum number of conss in sstconss */
303  int sstleaderrule; /**< rule to select leader */
304  int ssttiebreakrule; /**< tie break rule for leader selection */
305  int sstleadervartype; /**< bitset encoding which variable types can be leaders;
306  * if multiple types are allowed, take the one with most affected vars */
307  int* leaders; /**< index of orbit leaders in permvars */
308  int nleaders; /**< number of orbit leaders in leaders array */
309  int maxnleaders; /**< maximum number of leaders in leaders array */
310  SCIP_Bool addconflictcuts; /**< Should Schreier Sims constraints be added if we use a conflict based rule? */
311  SCIP_Bool sstaddcuts; /**< Should Schreier Sims constraints be added? */
312  SCIP_Bool sstmixedcomponents; /**< Should Schreier Sims constraints be added if a symmetry component contains variables of different types? */
313 
314  SCIP_EVENTHDLR* shadowtreeeventhdlr;/**< pointer to event handler for shadow tree */
315  SCIP_ORBITOPALREDDATA* orbitopalreddata; /**< container for the orbitopal reduction data */
316  SCIP_ORBITALREDDATA* orbitalreddata; /**< container for orbital reduction data */
317  SCIP_LEXREDDATA* lexreddata; /**< container for lexicographic reduction propagation */
318 };
319 
320 /** conflict data structure for SST cuts */
321 struct SCIP_ConflictData
322 {
323  SCIP_VAR* var; /**< variable belonging to node */
324  int orbitidx; /**< orbit of variable w.r.t. current stabilizer subgroup
325  * or -1 if not affected by symmetry */
326  int nconflictinorbit; /**< number of variables the node's var is in conflict with */
327  int orbitsize; /**< size of the variable's orbit */
328  int posinorbit; /**< position of variable in its orbit */
329  SCIP_Bool active; /**< whether variable has not been fixed by Schreier Sims code */
330  SCIP_CLIQUE** cliques; /**< List of setppc constraints. */
331  int ncliques; /**< Number of setppc constraints. */
332 };
333 typedef struct SCIP_ConflictData SCIP_CONFLICTDATA;
335 
336 /** compare function for sorting an array by the addresses of its members */
337 static
338 SCIP_DECL_SORTPTRCOMP(sortByPointerValue)
339 {
340  /* @todo move to misc.c? */
341  if ( elem1 < elem2 )
342  return -1;
343  else if ( elem1 > elem2 )
344  return +1;
345  return 0;
346 }
347 
348 
349 /** checks whether two arrays that are sorted with the same comparator have a common element */
350 static
352  void** arr1, /**< first array */
353  int narr1, /**< number of elements in first array */
354  void** arr2, /**< second array */
355  int narr2, /**< number of elements in second array */
356  SCIP_DECL_SORTPTRCOMP((*compfunc)) /**< comparator function that was used to sort arr1 and arr2; must define a total ordering */
357  )
358 {
359  /* @todo move to misc.c? */
360  int it1;
361  int it2;
362  int cmp;
363 
364  assert( arr1 != NULL || narr1 == 0 );
365  assert( narr1 >= 0 );
366  assert( arr2 != NULL || narr2 == 0 );
367  assert( narr2 >= 0 );
368  assert( compfunc != NULL );
369 
370  /* there is no overlap if one of the two arrays is empty */
371  if ( narr1 <= 0 )
372  return FALSE;
373  if ( narr2 <= 0 )
374  return FALSE;
375 
376  it1 = 0;
377  it2 = 0;
378 
379  while ( TRUE ) /*lint !e716*/
380  {
381  cmp = compfunc(arr1[it1], arr2[it2]);
382  if ( cmp < 0 )
383  {
384  /* comparison function determines arr1[it1] < arr2[it2]
385  * increase iterator for arr1
386  */
387  if ( ++it1 >= narr1 )
388  break;
389  continue;
390  }
391  else if ( cmp > 0 )
392  {
393  /* comparison function determines arr1[it1] > arr2[it2]
394  * increase iterator for arr2
395  */
396  if ( ++it2 >= narr2 )
397  break;
398  continue;
399  }
400  else
401  {
402  /* the entries arr1[it1] and arr2[it2] are the same with respect to the comparison function */
403  assert( cmp == 0 );
404  return TRUE;
405  }
406  }
407 
408  /* no overlap detected */
409  assert( it1 >= narr1 || it2 >= narr2 );
410  return FALSE;
411 }
412 
413 
414 /*
415  * Display dialog callback methods
416  */
417 
418 /** displays the cycle of a symmetry */
419 static
421  SCIP* scip, /**< SCIP pointer */
422  int* perm, /**< symmetry */
423  SYM_SYMTYPE symtype, /**< type of symmetry */
424  int baseidx, /**< variable index for which cycle is computed */
425  SCIP_Bool* covered, /**< allocated array to store covered variables */
426  int nvars, /**< number of (non-negated) variables in symmetry */
427  SCIP_VAR** vars /**< variables on which symmetry acts */
428  )
429 {
430  char* string;
431  int varidx;
432  int j;
433 
434  assert( scip != NULL );
435  assert( perm != NULL );
436  assert( 0 <= baseidx );
437  assert( (symtype == SYM_SYMTYPE_PERM && baseidx < nvars) ||
438  (symtype == SYM_SYMTYPE_SIGNPERM && baseidx < 2 * nvars) );
439  assert( covered != NULL );
440 
441  /* skip fixed points or elements already covered in previous cycle */
442  if ( perm[baseidx] == baseidx || covered[baseidx] )
443  return SCIP_OKAY;
444 
445  varidx = baseidx >= nvars ? baseidx - nvars : baseidx;
446  string = (char*) SCIPvarGetName(vars[varidx]);
447  SCIPinfoMessage(scip, NULL, " (%s<%s>", baseidx >= nvars ? "negated " : "", string);
448  j = perm[baseidx];
449  covered[baseidx] = TRUE;
450  while ( j != baseidx )
451  {
452  covered[j] = TRUE;
453  varidx = j >= nvars ? j - nvars : j;
454  string = (char*) SCIPvarGetName(vars[varidx]);
455  SCIPinfoMessage(scip, NULL, ",%s<%s>", j >= nvars ? "negated " : "", string);
456  j = perm[j];
457  }
458  SCIPinfoMessage(scip, NULL, ")\n");
459 
460  return SCIP_OKAY;
461 }
462 
463 /** displays symmetry information without taking components into account */
464 static
466  SCIP* scip, /**< SCIP pointer */
467  SCIP_PROPDATA* propdata /**< propagator data */
468  )
469 {
470  SCIP_Bool* covered;
471  SYM_SYMTYPE symtype;
472  int* perm;
473  int permlen;
474  int npermvars;
475  int i;
476  int p;
477 
478  assert( scip != NULL );
479  assert( propdata != NULL );
480  assert( propdata->nperms > 0 );
481  assert( propdata->permvars != NULL );
482  assert( propdata->npermvars > 0 );
483 
484  symtype = (SYM_SYMTYPE) propdata->symtype;
485  npermvars = propdata->npermvars;
486  permlen = symtype == SYM_SYMTYPE_PERM ? npermvars : 2 * npermvars;
487 
488  if ( symtype == SYM_SYMTYPE_SIGNPERM )
489  SCIPinfoMessage(scip, NULL, "Display permutations as signed permutations (allowing translations)\n");
490 
491  SCIP_CALL( SCIPallocClearBufferArray(scip, &covered, permlen) );
492 
493  for (p = 0; p < propdata->nperms; ++p)
494  {
495  SCIPinfoMessage(scip, NULL, "Permutation %d:\n", p);
496  perm = propdata->perms[p];
497 
498  for (i = 0; i < permlen; ++i)
499  {
500  SCIP_CALL( displayCycleOfSymmetry(scip, perm, symtype, i, covered, npermvars, propdata->permvars) );
501  }
502 
503  for (i = 0; i < permlen; ++i)
504  covered[i] = FALSE;
505  }
506 
507  SCIPfreeBufferArray(scip, &covered);
508 
509  return SCIP_OKAY;
510 }
511 
512 /** displays symmetry information taking components into account */
513 static
515  SCIP* scip, /**< SCIP pointer */
516  SCIP_PROPDATA* propdata /**< propagator data */
517  )
518 {
519  SCIP_Bool* covered;
520  SYM_SYMTYPE symtype;
521  int* perm;
522  int comppermlen;
523  int permlen;
524  int npermvars;
525  int i;
526  int p;
527  int c;
528 
529  assert( scip != NULL );
530  assert( propdata != NULL );
531  assert( propdata->nperms > 0 );
532  assert( propdata->permvars != NULL );
533  assert( propdata->npermvars > 0 );
534  assert( propdata->ncomponents > 0 );
535 
536  symtype = (SYM_SYMTYPE) propdata->symtype;
537  npermvars = propdata->npermvars;
538  permlen = symtype == SYM_SYMTYPE_PERM ? npermvars : 2 * npermvars;
539 
540  SCIP_CALL( SCIPallocClearBufferArray(scip, &covered, permlen) );
541 
542  for (c = 0; c < propdata->ncomponents; ++c)
543  {
544  int cnt;
545 
546  SCIPinfoMessage(scip, NULL, "Display symmetries of component %d.\n", c);
547  if ( propdata->componenthassignedperm[c] )
548  SCIPinfoMessage(scip, NULL, " Symmetries are displayed as signed permutations (allowing translations).\n");
549  else
550  SCIPinfoMessage(scip, NULL, " Symmetries are displayed as permutations.\n");
551 
552  comppermlen = propdata->componenthassignedperm[c] ? 2 * npermvars : npermvars;
553 
554  for (p = propdata->componentbegins[c], cnt = 0; p < propdata->componentbegins[c + 1]; ++p, ++cnt)
555  {
556  SCIPinfoMessage(scip, NULL, "Permutation %d:\n", cnt);
557  perm = propdata->perms[propdata->components[p]];
558 
559  for (i = 0; i < comppermlen; ++i)
560  {
561  SCIP_CALL( displayCycleOfSymmetry(scip, perm, symtype, i, covered, npermvars, propdata->permvars) );
562  }
563 
564  for (i = 0; i < comppermlen; ++i)
565  covered[i] = FALSE;
566  }
567  }
568 
569  SCIPfreeBufferArray(scip, &covered);
570 
571  return SCIP_OKAY;
572 }
573 
574 /** dialog execution method for the display symmetry information command */
575 static
576 SCIP_DECL_DIALOGEXEC(dialogExecDisplaySymmetry)
577 { /*lint --e{715}*/
578  SCIP_PROPDATA* propdata;
579 
580  /* add your dialog to history of dialogs that have been executed */
581  SCIP_CALL( SCIPdialoghdlrAddHistory(dialoghdlr, dialog, NULL, FALSE) );
582 
583  propdata = (SCIP_PROPDATA*)SCIPdialogGetData(dialog);
584  assert( propdata != NULL );
585 
586  if ( propdata->nperms == -1 )
587  {
588  SCIPinfoMessage(scip, NULL, "Cannot display symmetries. Symmetries have not been computed yet.\n");
589  }
590  else if ( propdata->nperms == 0 )
591  {
592  SCIPinfoMessage(scip, NULL, "Cannot display symmetries. No symmetries detected.\n");
593  }
594  else if ( propdata->ncomponents < 0 )
595  {
597  }
598  else
599  {
600  SCIP_CALL( displaySymmetriesWithComponents(scip, propdata) );
601  }
602 
603  /* next dialog will be root dialog again */
604  *nextdialog = SCIPdialoghdlrGetRoot(dialoghdlr);
605 
606  return SCIP_OKAY;
607 }
608 
609 
610 /*
611  * Table callback methods
612  */
613 
614 /** table data */
615 struct SCIP_TableData
616 {
617  SCIP_PROPDATA* propdata; /** pass data of propagator for table output function */
618 };
619 
620 
621 /** output method of symmetry propagator statistics table to output file stream 'file' */
622 static
623 SCIP_DECL_TABLEOUTPUT(tableOutputSymmetry)
624 {
625  SCIP_TABLEDATA* tabledata;
626  int nred;
627  int ncutoff;
628  SCIP_Real time;
629 
630  assert( scip != NULL );
631  assert( table != NULL );
632 
633  tabledata = SCIPtableGetData(table);
634  assert( tabledata != NULL );
635  assert( tabledata->propdata != NULL );
636 
637  if ( tabledata->propdata->orbitopalreddata || tabledata->propdata->orbitalreddata
638  || tabledata->propdata->lexreddata )
639  {
640  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, file, "Symmetry :\n");
641  if ( tabledata->propdata->orbitopalreddata )
642  {
643  SCIP_CALL( SCIPorbitopalReductionGetStatistics(scip, tabledata->propdata->orbitopalreddata, &nred, &ncutoff) );
644  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, file, " orbitopal reduction: %10d reductions applied,"
645  " %10d cutoffs\n", nred, ncutoff);
646  }
647  if ( tabledata->propdata->orbitalreddata )
648  {
649  SCIP_CALL( SCIPorbitalReductionGetStatistics(scip, tabledata->propdata->orbitalreddata, &nred, &ncutoff) );
650  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, file, " orbital reduction: %10d reductions applied,"
651  " %10d cutoffs\n", nred, ncutoff);
652  }
653  if ( tabledata->propdata->lexreddata )
654  {
655  SCIP_CALL( SCIPlexicographicReductionGetStatistics(scip, tabledata->propdata->lexreddata, &nred, &ncutoff) );
656  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, file, " lexicographic red: %10d reductions applied,"
657  " %10d cutoffs\n", nred, ncutoff);
658  }
659  if ( tabledata->propdata->shadowtreeeventhdlr )
660  {
661  time = SCIPgetShadowTreeEventHandlerExecutionTime(scip, tabledata->propdata->shadowtreeeventhdlr);
662  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, file, " shadow tree time : %10.2f s\n", time);
663  }
664  }
665 
666  return SCIP_OKAY;
667 }
668 
669 
670 /** destructor of statistics table to free user data (called when SCIP is exiting) */
671 static
672 SCIP_DECL_TABLEFREE(tableFreeSymmetry)
673 {
674  SCIP_TABLEDATA* tabledata;
675  tabledata = SCIPtableGetData(table);
676  assert( tabledata != NULL );
677 
678  SCIPfreeBlockMemory(scip, &tabledata);
679 
680  return SCIP_OKAY;
681 }
682 
683 
684 
685 /*
686  * local data structures
687  */
688 
689 /** data structure to store arrays used for sorting colored component types */
691 {
692  int* components; /**< array of components */
693  int* colors; /**< array of colors */
694 };
697 /** sorts variable indices according to their corresponding component in the graph
698  *
699  * Variables are sorted first by the color of their component and then by the component index.
700  *
701  * result:
702  * < 0: ind1 comes before (is better than) ind2
703  * = 0: both indices have the same value
704  * > 0: ind2 comes after (is worse than) ind2
705  */
706 static
707 SCIP_DECL_SORTINDCOMP(SYMsortGraphCompVars)
708 {
709  SYM_SORTGRAPHCOMPVARS* data;
710 
711  data = (SYM_SORTGRAPHCOMPVARS*) dataptr;
712 
713  if ( data->colors[ind1] < data->colors[ind2] )
714  return -1;
715  else if ( data->colors[ind1] > data->colors[ind2] )
716  return 1;
717 
718  if ( data->components[ind1] < data->components[ind2] )
719  return -1;
720  if ( data->components[ind1] > data->components[ind2] )
721  return 1;
722 
723  return 0;
724 }
725 
726 /** compares two symmetry detection graphs
727  *
728  * Graphs are compared by their number of consnodes, then their constypes, then by their lhs,
729  * then by their rhs, then by their total number of nodes, then by the number of operator nodes,
730  * then by their number of value nodes, and then by their number of edges.
731  *
732  * result:
733  * < 0: ind1 comes before (is better than) ind2
734  * = 0: both indices have the same value
735  * > 0: ind2 comes after (is worse than) ind2
736  */
737 static
738 int compareSymgraphs(
739  SCIP* scip, /**< SCIP pointer (or NULL) */
740  SYM_GRAPH* G1, /**< first graph in comparison */
741  SYM_GRAPH* G2 /**< second graph in comparison */
742  )
743 {
744  int c;
745  int perm1;
746  int perm2;
747 
748  /* compare the number of constraint nodes */
749  if ( G1->nconsnodes < G2->nconsnodes )
750  return -1;
751  if ( G1->nconsnodes > G2->nconsnodes )
752  return 1;
753 
754  /* compare the constraint nodes of the two graphs */
755  for (c = 0; c < G1->nconsnodes; ++c)
756  {
757  perm1 = G1->consnodeperm[c];
758  perm2 = G2->consnodeperm[c];
759 
760  if ( SCIPconsGetHdlr(G1->conss[perm1]) < SCIPconsGetHdlr(G2->conss[perm2]) )
761  return -1;
762  if ( SCIPconsGetHdlr(G1->conss[perm1]) > SCIPconsGetHdlr(G2->conss[perm2]) )
763  return 1;
764 
765  /* compare using SCIP functions when SCIP is available */
766  if ( scip != NULL )
767  {
768  if ( SCIPisLT(scip, G1->lhs[perm1], G2->lhs[perm2]) )
769  return -1;
770  if ( SCIPisGT(scip, G1->lhs[perm1], G2->lhs[perm2]) )
771  return 1;
772 
773  if ( SCIPisLT(scip, G1->rhs[perm1], G2->rhs[perm2]) )
774  return -1;
775  if ( SCIPisGT(scip, G1->rhs[perm1], G2->rhs[perm2]) )
776  return 1;
777  }
778  else
779  {
780  if ( G1->lhs[perm1] < G2->lhs[perm2] )
781  return -1;
782  if ( G1->lhs[perm1] > G2->lhs[perm2] )
783  return 1;
784 
785  if ( G1->rhs[perm1] < G2->rhs[perm2] )
786  return -1;
787  if ( G1->rhs[perm1] > G2->rhs[perm2] )
788  return 1;
789  }
790  }
791 
792  /* compare number of remaining node types */
793  if ( G1->nnodes < G2->nnodes )
794  return -1;
795  if ( G1->nnodes > G2->nnodes )
796  return 1;
797 
798  if ( G1->nopnodes < G2->nopnodes )
799  return -1;
800  if ( G1->nopnodes > G2->nopnodes )
801  return 1;
802 
803  if ( G1->nvalnodes < G2->nvalnodes )
804  return -1;
805  if ( G1->nvalnodes > G2->nvalnodes )
806  return 1;
807 
808  if ( G1->nedges < G2->nedges )
809  return -1;
810  if ( G1->nedges > G2->nedges )
811  return 1;
812 
813  return 0;
814 }
815 
816 /** sorts symmetry detection graphs
817  *
818  * Graphs are sorted by their number of consnodes, then their constypes, then by their lhs,
819  * then by their rhs, then by their total number of nodes, then by the number of operator nodes,
820  * then by their number of value nodes, and then by their number of edges.
821  *
822  * result:
823  * < 0: ind1 comes before (is better than) ind2
824  * = 0: both indices have the same value
825  * > 0: ind2 comes after (is worse than) ind2
826  */
827 static
828 SCIP_DECL_SORTINDCOMP(SYMsortSymgraphs)
829 {
830  SYM_GRAPH** data;
831  SYM_GRAPH* G1;
832  SYM_GRAPH* G2;
833 
834  data = (SYM_GRAPH**) dataptr;
835  G1 = data[ind1];
836  G2 = data[ind2];
837 
838  return compareSymgraphs(NULL, G1, G2);
839 }
840 
841 /*
842  * Local methods
843  */
844 
845 #ifndef NDEBUG
846 /** checks that symmetry data is all freed */
847 static
849  SCIP_PROPDATA* propdata /**< propagator data */
850  )
851 {
852  assert( propdata->permvarmap == NULL );
853  assert( propdata->genorbconss == NULL );
854  assert( propdata->genlinconss == NULL );
855  assert( propdata->ngenlinconss == 0 );
856  assert( propdata->ngenorbconss == 0 );
857  assert( propdata->genorbconsssize == 0 );
858  assert( propdata->genlinconsssize == 0 );
859  assert( propdata->sstconss == NULL );
860  assert( propdata->leaders == NULL );
861 
862  assert( propdata->permvardomaincenter == NULL );
863  assert( propdata->permvars == NULL );
864  assert( propdata->perms == NULL );
865  assert( propdata->permstrans == NULL );
866  assert( propdata->npermvars == 0 );
867  assert( propdata->nbinpermvars == 0 );
868  assert( propdata->nperms == -1 || propdata->nperms == 0 );
869  assert( propdata->nmaxperms == 0 );
870  assert( propdata->nmovedpermvars == -1 );
871  assert( propdata->nmovedbinpermvars == 0 );
872  assert( propdata->nmovedintpermvars == 0 );
873  assert( propdata->nmovedimplintpermvars == 0 );
874  assert( propdata->nmovedcontpermvars == 0 );
875  assert( propdata->nmovedvars == -1 );
876  assert( propdata->binvaraffected == FALSE );
877  assert( propdata->isnonlinvar == NULL );
878 
879  assert( propdata->componenthassignedperm == NULL );
880  assert( propdata->componentblocked == NULL );
881  assert( propdata->componentbegins == NULL );
882  assert( propdata->components == NULL );
883  assert( propdata->ncomponents == -1 );
884  assert( propdata->ncompblocked == 0 );
885 
886  return TRUE;
887 }
888 #endif
889 
890 
891 /** resets symmetry handling propagators that depend on the branch-and-bound tree structure */
892 static
894  SCIP* scip, /**< SCIP pointer */
895  SCIP_PROPDATA* propdata /**< propagator data */
896  )
897 {
898  assert( scip != NULL );
899  assert( propdata != NULL );
900 
901  /* propagators managed by a different file */
902  if ( propdata->orbitalreddata != NULL )
903  {
904  SCIP_CALL( SCIPorbitalReductionReset(scip, propdata->orbitalreddata) );
905  }
906  if ( propdata->orbitopalreddata != NULL )
907  {
908  SCIP_CALL( SCIPorbitopalReductionReset(scip, propdata->orbitopalreddata) );
909  }
910  if ( propdata->lexreddata != NULL )
911  {
912  SCIP_CALL( SCIPlexicographicReductionReset(scip, propdata->lexreddata) );
913  }
914 
915  return SCIP_OKAY;
916 }
917 
918 
919 /** frees symmetry data */
920 static
922  SCIP* scip, /**< SCIP pointer */
923  SCIP_PROPDATA* propdata /**< propagator data */
924  )
925 {
926  int i;
927 
928  assert( scip != NULL );
929  assert( propdata != NULL );
930 
931  SCIP_CALL( resetDynamicSymmetryHandling(scip, propdata) );
932 
933  if ( propdata->permvarmap != NULL )
934  {
935  SCIPhashmapFree(&propdata->permvarmap);
936  }
937 
938  /* release all variables contained in permvars array */
939  for (i = 0; i < propdata->npermvars; ++i)
940  {
941  assert( propdata->permvars[i] != NULL );
942  SCIP_CALL( SCIPreleaseVar(scip, &propdata->permvars[i]) );
943  }
944 
945  /* free permstrans matrix*/
946  if ( propdata->permstrans != NULL )
947  {
948  assert( propdata->nperms > 0 );
949  assert( propdata->permvars != NULL );
950  assert( propdata->npermvars > 0 );
951  assert( propdata->nmaxperms > 0 );
952 
953  for (i = 0; i < propdata->npermvars; ++i)
954  {
955  SCIPfreeBlockMemoryArray(scip, &propdata->permstrans[i], propdata->nmaxperms);
956  }
957  SCIPfreeBlockMemoryArray(scip, &propdata->permstrans, propdata->npermvars);
958  }
959 
960  /* free data of added orbitope/orbisack/symresack constraints */
961  if ( propdata->genorbconss != NULL )
962  {
963  assert( propdata->ngenorbconss > 0 );
964 
965  /* release constraints */
966  while ( propdata->ngenorbconss > 0 )
967  {
968  assert( propdata->genorbconss[propdata->ngenorbconss - 1] != NULL );
969  SCIP_CALL( SCIPreleaseCons(scip, &propdata->genorbconss[--propdata->ngenorbconss]) );
970  }
971  assert( propdata->ngenorbconss == 0 );
972 
973  /* free pointers to symmetry group and binary variables */
974  SCIPfreeBlockMemoryArray(scip, &propdata->genorbconss, propdata->genorbconsssize);
975  propdata->genorbconsssize = 0;
976  }
977 
978  /* free data of added constraints */
979  if ( propdata->genlinconss != NULL )
980  {
981  /* release constraints */
982  for (i = 0; i < propdata->ngenlinconss; ++i)
983  {
984  assert( propdata->genlinconss[i] != NULL );
985  SCIP_CALL( SCIPreleaseCons(scip, &propdata->genlinconss[i]) );
986  }
987 
988  /* free pointers to symmetry group and binary variables */
989  SCIPfreeBlockMemoryArray(scip, &propdata->genlinconss, propdata->genlinconsssize);
990  propdata->ngenlinconss = 0;
991  propdata->genlinconsssize = 0;
992  }
993 
994  if ( propdata->sstconss != NULL )
995  {
996  assert( propdata->nsstconss > 0 );
997 
998  /* release constraints */
999  for (i = 0; i < propdata->nsstconss; ++i)
1000  {
1001  assert( propdata->sstconss[i] != NULL );
1002  SCIP_CALL( SCIPreleaseCons(scip, &propdata->sstconss[i]) );
1003  }
1004 
1005  /* free pointers to symmetry group and binary variables */
1006  SCIPfreeBlockMemoryArray(scip, &propdata->sstconss, propdata->maxnsstconss);
1007  propdata->sstconss = NULL;
1008  propdata->nsstconss = 0;
1009  propdata->maxnsstconss = 0;
1010  }
1011 
1012  if ( propdata->leaders != NULL )
1013  {
1014  assert( propdata->maxnleaders > 0 );
1015 
1016  SCIPfreeBlockMemoryArray(scip, &propdata->leaders, propdata->maxnleaders);
1017  propdata->maxnleaders = 0;
1018  propdata->leaders = NULL;
1019  propdata->nleaders = 0;
1020  }
1021 
1022  /* free components */
1023  if ( propdata->ncomponents > 0 )
1024  {
1025  assert( propdata->componentblocked != NULL );
1026  assert( propdata->vartocomponent != NULL );
1027  assert( propdata->componentbegins != NULL );
1028  assert( propdata->components != NULL );
1029 
1030  SCIPfreeBlockMemoryArray(scip, &propdata->componenthassignedperm, propdata->ncomponents);
1031  SCIPfreeBlockMemoryArray(scip, &propdata->componentblocked, propdata->ncomponents);
1032  SCIPfreeBlockMemoryArray(scip, &propdata->vartocomponent, propdata->npermvars);
1033  SCIPfreeBlockMemoryArray(scip, &propdata->componentbegins, propdata->ncomponents + 1);
1034  SCIPfreeBlockMemoryArray(scip, &propdata->components, propdata->nperms);
1035 
1036  propdata->ncomponents = -1;
1037  propdata->ncompblocked = 0;
1038  }
1039 
1040  /* free main symmetry data */
1041  if ( propdata->nperms > 0 )
1042  {
1043  int permlen;
1044 
1045  assert( propdata->permvars != NULL );
1046 
1047  if ( (SYM_SYMTYPE) propdata->symtype == SYM_SYMTYPE_SIGNPERM )
1048  permlen = 2 * propdata->npermvars;
1049  else
1050  permlen = propdata->npermvars;
1051 
1052  SCIPfreeBlockMemoryArray(scip, &propdata->permvars, propdata->npermvars);
1053  SCIPfreeBlockMemoryArray(scip, &propdata->permvardomaincenter, propdata->npermvars);
1054 
1055  if ( propdata->perms != NULL )
1056  {
1057  for (i = 0; i < propdata->nperms; ++i)
1058  {
1059  SCIPfreeBlockMemoryArray(scip, &propdata->perms[i], permlen);
1060  }
1061  SCIPfreeBlockMemoryArray(scip, &propdata->perms, propdata->nmaxperms);
1062  }
1063 
1064  SCIPfreeBlockMemoryArrayNull(scip, &propdata->isnonlinvar, propdata->npermvars);
1065 
1066  propdata->npermvars = 0;
1067  propdata->nbinpermvars = 0;
1068  propdata->nperms = -1;
1069  propdata->nmaxperms = 0;
1070  propdata->nmovedpermvars = -1;
1071  propdata->nmovedbinpermvars = 0;
1072  propdata->nmovedintpermvars = 0;
1073  propdata->nmovedimplintpermvars = 0;
1074  propdata->nmovedcontpermvars = 0;
1075  propdata->nmovedvars = -1;
1076  propdata->log10groupsize = -1.0;
1077  propdata->binvaraffected = FALSE;
1078  propdata->isnonlinvar = NULL;
1079  }
1080  propdata->nperms = -1;
1081 
1082  assert( checkSymmetryDataFree(propdata) );
1083 
1084  propdata->computedsymmetry = FALSE;
1085  propdata->compressed = FALSE;
1086 
1087  return SCIP_OKAY;
1088 }
1089 
1090 
1091 /** makes sure that the constraint array (potentially NULL) of given array size is sufficiently large */
1092 static
1094  SCIP* scip, /**< SCIP pointer */
1095  SCIP_CONS*** consarrptr, /**< constraint array pointer */
1096  int* consarrsizeptr, /**< constraint array size pointer */
1097  int consarrsizereq /**< constraint array size required */
1098  )
1099 {
1100  int newsize;
1101 
1102  assert( scip != NULL );
1103  assert( consarrptr != NULL );
1104  assert( consarrsizeptr != NULL );
1105  assert( consarrsizereq > 0 );
1106  assert( *consarrsizeptr >= 0 );
1107  assert( (*consarrsizeptr == 0) == (*consarrptr == NULL) );
1108 
1109  /* array is already sufficiently large */
1110  if ( consarrsizereq <= *consarrsizeptr )
1111  return SCIP_OKAY;
1112 
1113  /* compute new size */
1114  newsize = SCIPcalcMemGrowSize(scip, consarrsizereq);
1115  assert( newsize > *consarrsizeptr );
1116 
1117  /* allocate or reallocate */
1118  if ( *consarrptr == NULL )
1119  {
1120  SCIP_CALL( SCIPallocBlockMemoryArray(scip, consarrptr, newsize) );
1121  }
1122  else
1123  {
1124  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, consarrptr, *consarrsizeptr, newsize) );
1125  }
1126 
1127  *consarrsizeptr = newsize;
1128 
1129  return SCIP_OKAY;
1130 }
1131 
1132 /** set symmetry data */
1133 static
1135  SCIP* scip, /**< SCIP pointer */
1136  SYM_SYMTYPE symtype, /**< type of symmetries in perms */
1137  SCIP_VAR** vars, /**< vars present at time of symmetry computation */
1138  int nvars, /**< number of vars present at time of symmetry computation */
1139  int nbinvars, /**< number of binary vars present at time of symmetry computation */
1140  SCIP_VAR*** permvars, /**< pointer to permvars array */
1141  int* npermvars, /**< pointer to store number of permvars */
1142  int* nbinpermvars, /**< pointer to store number of binary permvars */
1143  SCIP_Real** permvardomaincenter, /**< pointer to store center points of variable domains */
1144  int** perms, /**< permutations matrix (nperms x nvars) */
1145  int nperms, /**< number of permutations */
1146  int* nmovedvars, /**< pointer to store number of vars affected by symmetry (if usecompression) or NULL */
1147  SCIP_Bool* binvaraffected, /**< pointer to store whether a binary variable is affected by symmetry */
1148  SCIP_Bool usecompression, /**< whether symmetry data shall be compressed */
1149  SCIP_Real compressthreshold, /**< if percentage of moved vars is at most threshold, compression is done */
1150  SCIP_Bool* compressed /**< pointer to store whether compression has been performed */
1151  )
1152 {
1153  SCIP_Real ub;
1154  SCIP_Real lb;
1155  int i;
1156  int p;
1157 
1158  assert( scip != NULL );
1159  assert( vars != NULL );
1160  assert( nvars > 0 );
1161  assert( permvars != NULL );
1162  assert( npermvars != NULL );
1163  assert( nbinpermvars != NULL );
1164  assert( perms != NULL );
1165  assert( nperms > 0 );
1166  assert( binvaraffected != NULL );
1167  assert( SCIPisGE(scip, compressthreshold, 0.0) );
1168  assert( SCIPisLE(scip, compressthreshold, 1.0) );
1169  assert( compressed != NULL );
1170 
1171  /* set default return values */
1172  *permvars = vars;
1173  *npermvars = nvars;
1174  *nbinpermvars = nbinvars;
1175  *binvaraffected = FALSE;
1176  *compressed = FALSE;
1177 
1178  /* if we possibly perform compression */
1179  if ( usecompression && SCIPgetNVars(scip) >= COMPRESSNVARSLB )
1180  {
1181  SCIP_Real percentagemovedvars;
1182  int* labelmovedvars;
1183  int* labeltopermvaridx;
1184  int nbinvarsaffected = 0;
1185 
1186  assert( nmovedvars != NULL );
1187 
1188  *nmovedvars = 0;
1189 
1190  /* detect number of moved vars and label moved vars */
1191  SCIP_CALL( SCIPallocBufferArray(scip, &labelmovedvars, nvars) );
1192  SCIP_CALL( SCIPallocBufferArray(scip, &labeltopermvaridx, nvars) );
1193  for (i = 0; i < nvars; ++i)
1194  {
1195  labelmovedvars[i] = -1;
1196 
1197  for (p = 0; p < nperms; ++p)
1198  {
1199  if ( perms[p][i] != i )
1200  {
1201  labeltopermvaridx[*nmovedvars] = i;
1202  labelmovedvars[i] = (*nmovedvars)++;
1203 
1204  if ( SCIPvarIsBinary(vars[i]) )
1205  ++nbinvarsaffected;
1206  break;
1207  }
1208  }
1209  }
1210 
1211  if ( nbinvarsaffected > 0 )
1212  *binvaraffected = TRUE;
1213 
1214  /* check whether compression should be performed */
1215  percentagemovedvars = (SCIP_Real) *nmovedvars / (SCIP_Real) nvars;
1216  if ( *nmovedvars > 0 && SCIPisLE(scip, percentagemovedvars, compressthreshold) )
1217  {
1218  /* remove variables from permutations that are not affected by any permutation */
1219  for (p = 0; p < nperms; ++p)
1220  {
1221  /* iterate over labels and adapt permutation (possibly taking signed permutations into account) */
1222  for (i = 0; i < *nmovedvars; ++i)
1223  {
1224  assert( i <= labeltopermvaridx[i] );
1225  if ( perms[p][labeltopermvaridx[i]] >= nvars )
1226  {
1227  int imgvaridx;
1228 
1229  assert( symtype == SYM_SYMTYPE_SIGNPERM );
1230 
1231  imgvaridx = perms[p][labeltopermvaridx[i]] - nvars;
1232  perms[p][i] = labelmovedvars[imgvaridx] + *nmovedvars;
1233 
1234  assert( 0 <= perms[p][i] && perms[p][i] < 2 * (*nmovedvars) );
1235  }
1236  else
1237  perms[p][i] = labelmovedvars[perms[p][labeltopermvaridx[i]]];
1238  }
1239 
1240  if ( symtype == SYM_SYMTYPE_SIGNPERM )
1241  {
1242  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &perms[p], 2 * nvars, 2 * (*nmovedvars)) );
1243  }
1244  else
1245  {
1246  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &perms[p], nvars, *nmovedvars) );
1247  }
1248  }
1249 
1250  /* remove variables from permvars array that are not affected by any symmetry */
1251  SCIP_CALL( SCIPallocBlockMemoryArray(scip, permvars, *nmovedvars) );
1252  for (i = 0; i < *nmovedvars; ++i)
1253  {
1254  (*permvars)[i] = vars[labeltopermvaridx[i]];
1255  }
1256  *npermvars = *nmovedvars;
1257  *nbinpermvars = nbinvarsaffected;
1258  *compressed = TRUE;
1259 
1260  SCIPfreeBlockMemoryArray(scip, &vars, nvars);
1261  }
1262  SCIPfreeBufferArray(scip, &labeltopermvaridx);
1263  SCIPfreeBufferArray(scip, &labelmovedvars);
1264  }
1265  else
1266  {
1267  /* detect whether binary variable is affected by symmetry and count number of binary permvars */
1268  for (i = 0; i < nbinvars; ++i)
1269  {
1270  for (p = 0; p < nperms && ! *binvaraffected; ++p)
1271  {
1272  if ( perms[p][i] != i )
1273  {
1274  if ( SCIPvarIsBinary(vars[i]) )
1275  *binvaraffected = TRUE;
1276  break;
1277  }
1278  }
1279  }
1280  }
1281 
1282  /* store center points of variable domains */
1283  SCIP_CALL( SCIPallocBlockMemoryArray(scip, permvardomaincenter, *npermvars) );
1284  for (i = 0; i < *npermvars; ++i)
1285  {
1286  ub = SCIPvarGetUbGlobal((*permvars)[i]);
1287  lb = SCIPvarGetLbGlobal((*permvars)[i]);
1288 
1289  (*permvardomaincenter)[i] = 0.5 * (ub + lb);
1290  }
1291 
1292  return SCIP_OKAY;
1293 }
1294 
1295 /** returns whether a constraint handler can provide required symmetry information */
1296 static
1298  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
1299  SYM_SYMTYPE symtype /**< type of symmetries for which information are needed */
1300  )
1301 {
1302  assert( conshdlr != NULL );
1303 
1304  switch ( symtype )
1305  {
1306  case SYM_SYMTYPE_PERM:
1307  return SCIPconshdlrSupportsPermsymDetection(conshdlr);
1308  default:
1309  assert( symtype == SYM_SYMTYPE_SIGNPERM );
1311  } /*lint !e788*/
1312 }
1313 
1314 /** returns whether all constraint handlers with constraints can provide symmetry information */
1315 static
1317  SCIP* scip, /**< SCIP pointer */
1318  SYM_SYMTYPE symtype /**< type of symmetries for which information are needed */
1319  )
1320 {
1321  SCIP_CONSHDLR** conshdlrs;
1322  SCIP_CONSHDLR* conshdlr;
1323  int nconshdlrs;
1324  int c;
1325 
1326  conshdlrs = SCIPgetConshdlrs(scip);
1327  assert( conshdlrs != NULL );
1328 
1329  nconshdlrs = SCIPgetNConshdlrs(scip);
1330  for (c = 0; c < nconshdlrs; ++c)
1331  {
1332  conshdlr = conshdlrs[c];
1333  assert( conshdlr != NULL );
1334 
1335  if ( ! conshdlrCanProvideSymInformation(conshdlr, symtype) && SCIPconshdlrGetNConss(conshdlr) > 0 )
1336  {
1337  char name[SCIP_MAXSTRLEN];
1338 
1339  if ( symtype == SYM_SYMTYPE_PERM )
1340  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "CONSGETPERMSYMGRAPH");
1341  else
1342  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "CONSGETSIGNEDPERMSYMGRAPH");
1343 
1345  " Symmetry detection interrupted: constraints of type %s do not provide symmetry information.\n"
1346  " If symmetries shall be detected, implement the %s callback.\n",
1347  SCIPconshdlrGetName(conshdlr), name);
1348 
1349  return FALSE;
1350  }
1351  }
1352 
1353  /* check whether all expressions provide sufficient symmetry information */
1354  conshdlr = SCIPfindConshdlr(scip, "nonlinear");
1355  if ( conshdlr != NULL && SCIPconshdlrGetNConss(conshdlr) > 0 )
1356  {
1357  SCIP_EXPRHDLR* exprhdlr;
1358 
1359  for (c = 0; c < SCIPgetNExprhdlrs(scip); ++c)
1360  {
1361  SCIP_Bool found = FALSE;
1362  exprhdlr = SCIPgetExprhdlrs(scip)[c];
1363 
1364  if ( SCIPexprhdlrHasGetSymData(exprhdlr) )
1365  continue;
1366 
1367  /* check whether exprhdlr is known by SCIP (and handles symmetries correctly) */
1368  if ( strcmp(SCIPexprhdlrGetName(exprhdlr), "var") == 0
1369  || strcmp(SCIPexprhdlrGetName(exprhdlr), "sum") == 0
1370  || strcmp(SCIPexprhdlrGetName(exprhdlr), "product") == 0
1371  || strcmp(SCIPexprhdlrGetName(exprhdlr), "val") == 0
1372  || strcmp(SCIPexprhdlrGetName(exprhdlr), "pow") == 0
1373  || strcmp(SCIPexprhdlrGetName(exprhdlr), "signpow") == 0
1374  || strcmp(SCIPexprhdlrGetName(exprhdlr), "exp") == 0
1375  || strcmp(SCIPexprhdlrGetName(exprhdlr), "log") == 0
1376  || strcmp(SCIPexprhdlrGetName(exprhdlr), "abs") == 0
1377  || strcmp(SCIPexprhdlrGetName(exprhdlr), "sin") == 0
1378  || strcmp(SCIPexprhdlrGetName(exprhdlr), "cos") == 0
1379  || strcmp(SCIPexprhdlrGetName(exprhdlr), "entropy") == 0
1380  || strcmp(SCIPexprhdlrGetName(exprhdlr), "erf") == 0
1381  || strcmp(SCIPexprhdlrGetName(exprhdlr), "varidx") == 0 )
1382  found = TRUE;
1383 
1384  /* there exists an unknown expression handler that does not provide symmetry information */
1385  if ( ! found )
1386  {
1387  SCIPwarningMessage(scip, "Expression handler %s does not implement the EXPRGETSYMDATA callback.\n"
1388  "Computed symmetries might be incorrect if the expression uses different constants or assigns\n"
1389  "different coefficients to its children.\n", SCIPexprhdlrGetName(SCIPgetExprhdlrs(scip)[c]));
1390  }
1391  }
1392  }
1393 
1394  return TRUE;
1395 }
1396 
1397 /** provides estimates for the number of nodes and edges in a symmetry detection graph */
1398 static
1400  SCIP* scip, /**< SCIP pointer */
1401  int* nopnodes, /**< pointer to store estimate for number of operator nodes */
1402  int* nvalnodes, /**< pointer to store estimate for number of value nodes */
1403  int* nconsnodes, /**< pointer to store estimate for number of constraint nodes */
1404  int* nedges /**< pointer to store estimate for number of edges */
1405  )
1406 {
1407  SCIP_CONS** conss;
1408  SCIP_Bool success;
1409  int nvars;
1410  int nconss;
1411  int num;
1412  int c;
1413 
1414  assert( scip != NULL );
1415  assert( nopnodes != NULL );
1416  assert( nvalnodes != NULL );
1417  assert( nconsnodes != NULL );
1418  assert( nedges != NULL );
1419 
1420  nvars = SCIPgetNVars(scip);
1421  nconss = SCIPgetNConss(scip);
1422  conss = SCIPgetConss(scip);
1423  assert( conss != NULL || nconss == 0 );
1424 
1425  *nconsnodes = nconss;
1426 
1427  /* get estimate from different types of constraints */
1428  *nopnodes = 0;
1429  *nvalnodes = 0;
1430  for (c = 0; c < nconss; ++c)
1431  {
1432  if ( strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(conss[c])), "bounddisjunction") == 0 )
1433  {
1434  SCIP_CALL( SCIPgetConsNVars(scip, conss[c], &num, &success) );
1435 
1436  if ( success )
1437  {
1438  *nopnodes += num;
1439  *nvalnodes += num;
1440  }
1441  }
1442  else if ( strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(conss[c])), "indicator") == 0 )
1443  {
1444  *nopnodes += 3;
1445  *nvalnodes += 1;
1446  }
1447  else if ( strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(conss[c])), "nonlinear") == 0 )
1448  {
1449  SCIP_CALL( SCIPgetConsNVars(scip, conss[c], &num, &success) );
1450 
1451  /* use binary trees as a proxy for an expression tree */
1452  if ( success )
1453  {
1454  int depth;
1455  int numnodes;
1456  int expval;
1457 
1458  depth = (int) log2((double) num);
1459  expval = (int) exp2((double) (depth + 1));
1460  numnodes = MIN(expval, 100);
1461 
1462  *nopnodes += numnodes;
1463  *nvalnodes += MAX((int) 0.1 * numnodes, 1);
1464  }
1465  }
1466  else if ( strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(conss[c])), "SOS1") == 0 )
1467  {
1468  SCIP_CALL( SCIPgetConsNVars(scip, conss[c], &num, &success) );
1469 
1470  if ( success )
1471  *nopnodes += num;
1472  }
1473  else if ( strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(conss[c])), "SOS2") == 0 )
1474  {
1475  SCIP_CALL( SCIPgetConsNVars(scip, conss[c], &num, &success) );
1476 
1477  if ( success )
1478  *nopnodes += num - 1;
1479  }
1480  }
1481 
1482  /* use a staggered scheme for the number of edges since this can become large
1483  *
1484  * In most cases, edges represent variable coefficients from linear constraints.
1485  * For this reason, use number of variables as proxy.
1486  */
1487  if ( nvars <= 100000 )
1488  *nedges = 100 * nvars;
1489  else if ( nvars <= 1000000 )
1490  *nedges = 32 * nvars;
1491  else if ( nvars <= 16700000 )
1492  *nedges = 16 * nvars;
1493  else
1494  *nedges = INT_MAX / 10;
1495 
1496  return SCIP_OKAY;
1497 }
1498 
1499 /** checks whether computed symmetries are indeed symmetries */
1500 static
1502  SCIP* scip, /**< SCIP pointer */
1503  SYM_SYMTYPE symtype, /**< type of symmetries to be checked */
1504  int** perms, /**< array of permutations */
1505  int nperms, /**< number of permutations */
1506  int npermvars, /**< number of variables permutations act on */
1507  SYM_SPEC fixedtype /**< variable types that must be fixed by symmetries */
1508  )
1509 {
1510  SYM_GRAPH** graphs;
1511  SCIP_CONS** conss;
1512  SCIP_VAR** symvars;
1513  SCIP_Bool success;
1514  int* graphperm;
1515  int* groupbegins;
1516  int ngroups = 1;
1517  int nsymvars;
1518  int nconss;
1519  int p;
1520  int c;
1521  int g;
1522 #ifdef SCIP_DISPLAY_SYM_CHECK
1523  int permlen;
1524  SCIP_Bool* covered;
1525 #endif
1526 
1527  assert( scip != NULL );
1528  assert( perms != NULL );
1529  assert( nperms > 0 );
1530  assert( npermvars > 0 );
1531 
1532  /* get symmetry detection graphs for all constraints */
1533  nconss = SCIPgetNConss(scip);
1534  conss = SCIPgetConss(scip);
1535  assert( conss != NULL );
1536 
1537  symvars = SCIPgetVars(scip);
1538  nsymvars = SCIPgetNVars(scip);
1539  assert( nsymvars == npermvars );
1540 
1541  SCIP_CALL( SCIPallocBufferArray(scip, &graphs, nconss) );
1542 
1543  for (c = 0; c < nconss; ++c)
1544  {
1545  SCIP_CALL( SCIPcreateSymgraph(scip, symtype, &graphs[c], symvars, nsymvars, 10, 10, 1, 100) );
1546 
1547  switch ( symtype )
1548  {
1549  case SYM_SYMTYPE_PERM:
1550  SCIP_CALL( SCIPgetConsPermsymGraph(scip, conss[c], graphs[c], &success) );
1551  break;
1552  default:
1553  assert( symtype == SYM_SYMTYPE_SIGNPERM );
1554  SCIP_CALL( SCIPgetConsSignedPermsymGraph(scip, conss[c], graphs[c], &success) );
1555  } /*lint !e788*/
1556 
1557  SCIP_CALL( SCIPcomputeSymgraphColors(scip, graphs[c], fixedtype) );
1558 
1559  assert( success );
1560  }
1561 
1562  /* sort graphs for quicker comparisons */
1563  SCIP_CALL( SCIPallocBufferArray(scip, &graphperm, nconss) );
1564  SCIP_CALL( SCIPallocBufferArray(scip, &groupbegins, nconss + 1) );
1565  for (c = 0; c < nconss; ++c)
1566  {
1567  SCIP_CALL( SCIPcreateSymgraphConsnodeperm(scip, graphs[c]) );
1568  }
1569 
1570  SCIPsort(graphperm, SYMsortSymgraphs, graphs, nconss);
1571 
1572  groupbegins[0] = 0;
1573  for (c = 1; c < nconss; ++c)
1574  {
1575  if ( compareSymgraphs(scip, graphs[graphperm[c]], graphs[graphperm[c-1]]) != 0 )
1576  groupbegins[ngroups++] = c;
1577  }
1578  groupbegins[ngroups] = nconss;
1579 
1580  /* remove information from symmetry detection graph that is not needed anymore */
1581  for (c = 0; c < nconss; ++c)
1582  {
1583  SCIP_CALL( SCIPfreeSymgraphConsnodeperm(scip, graphs[c]) );
1584  }
1585 
1586 #ifdef SCIP_DISPLAY_SYM_CHECK
1587  permlen = symtype == SYM_SYMTYPE_SIGNPERM ? 2 * npermvars : npermvars;
1588  SCIP_CALL( SCIPallocClearBufferArray(scip, &covered, permlen) );
1589 #endif
1590 
1591  /* iterate over all permutations and check whether they define symmetries */
1592  for (p = 0; p < nperms; ++p)
1593  {
1594  SYM_GRAPH* graph;
1595  SCIP_Bool found = TRUE;
1596  int d;
1597 #ifdef SCIP_DISPLAY_SYM_CHECK
1598  int i;
1599 
1600  SCIPinfoMessage(scip, NULL, "Check whether permutation %d is a symmetry:\n", p);
1601  for (i = 0; i < permlen; ++i)
1602  {
1603  SCIP_CALL( displayCycleOfSymmetry(scip, perms[p], symtype, i, covered, npermvars, SCIPgetVars(scip)) );
1604  }
1605 
1606  for (i = 0; i < permlen; ++i)
1607  covered[i] = FALSE;
1608  SCIPinfoMessage(scip, NULL, "Check whether every constraint has a symmetric counterpart.\n");
1609 #endif
1610 
1611  /* for every constraint, create permuted graph by copying nodes and edges */
1612  for (g = 0; g < ngroups; ++g)
1613  {
1614  for (c = groupbegins[g]; c < groupbegins[g+1]; ++c)
1615  {
1616 #ifdef SCIP_DISPLAY_SYM_CHECK
1617  SCIPinfoMessage(scip, NULL, "Check whether constraint %d has a symmetric counterpart:\n",
1618  graphperm[c]);
1619  SCIP_CALL( SCIPprintCons(scip, conss[graphperm[c]], NULL) );
1620  SCIPinfoMessage(scip, NULL, "\n");
1621 #endif
1622  SCIP_CALL( SCIPcopySymgraph(scip, &graph, graphs[graphperm[c]], perms[p], fixedtype) );
1623 
1624  /* if adapted graph is equivalent to original graph, we don't need to check further graphs */
1625  if ( SYMcheckGraphsAreIdentical(scip, symtype, graph, graphs[graphperm[c]]) )
1626  {
1627 #ifdef SCIP_DISPLAY_SYM_CHECK
1628  SCIPinfoMessage(scip, NULL, "\tconstraint is symmetric to itself\n");
1629 #endif
1630  SCIP_CALL( SCIPfreeSymgraph(scip, &graph) );
1631  continue;
1632  }
1633 
1634  /* check whether graph has an isomorphic counterpart */
1635  found = FALSE;
1636  for (d = groupbegins[g]; d < groupbegins[g+1] && ! found; ++d)
1637  {
1638  found = SYMcheckGraphsAreIdentical(scip, symtype, graph, graphs[graphperm[d]]);
1639 
1640 #ifdef SCIP_DISPLAY_SYM_CHECK
1641  SCIPinfoMessage(scip, NULL, "\tconstraint is %ssymmetric to constraint %d\n\t", !found ? "not " : "", d);
1642  SCIP_CALL( SCIPprintCons(scip, conss[graphperm[d]], NULL) );
1643  SCIPinfoMessage(scip, NULL, "\n");
1644 #endif
1645  }
1646 
1647  SCIP_CALL( SCIPfreeSymgraph(scip, &graph) );
1648 
1649  if ( ! found )
1650  {
1651 #ifdef SCIP_DISPLAY_SYM_CHECK
1652  SCIPfreeBufferArray(scip, &covered);
1653 #endif
1654  SCIPerrorMessage("permutation %d is not a symmetry\n", p);
1655  return SCIP_ERROR;
1656  }
1657  }
1658  }
1659  }
1660 
1661 #ifdef SCIP_DISPLAY_SYM_CHECK
1662  SCIPfreeBufferArray(scip, &covered);
1663 #endif
1664 
1665  SCIPfreeBufferArray(scip, &groupbegins);
1666  SCIPfreeBufferArray(scip, &graphperm);
1667 
1668  for (c = nconss - 1; c >= 0; --c)
1669  {
1670  SCIP_CALL( SCIPfreeSymgraph(scip, &graphs[c]) );
1671  }
1672  SCIPfreeBufferArray(scip, &graphs);
1673 
1674  return SCIP_OKAY;
1675 }
1676 
1677 /** computes symmetry group of a CIP */
1678 static
1680  SCIP* scip, /**< SCIP pointer */
1681  SYM_SYMTYPE symtype, /**< type of symmetries to be computed */
1682  SCIP_Bool compresssymmetries, /**< Should non-affected variables be removed from permutation to save memory? */
1683  SCIP_Real compressthreshold, /**< if percentage of moved vars is at most threshold, compression is done */
1684  int maxgenerators, /**< maximal number of generators constructed (= 0 if unlimited) */
1685  SYM_SPEC fixedtype, /**< variable types that must be fixed by symmetries */
1686  SCIP_Bool checksymmetries, /**< Should all symmetries be checked after computation? */
1687  SCIP_VAR*** permvars, /**< pointer to permvars array */
1688  int* npermvars, /**< pointer to store number of permvars */
1689  int* nbinpermvars, /**< pointer to store number of binary permvars */
1690  SCIP_Real** permvardomaincenter, /**< pointer to store center points of variable domains */
1691  int*** perms, /**< pointer to store permutation matrix (nperms x nvars) */
1692  int* nperms, /**< pointer to store number of permutations */
1693  int* nmaxperms, /**< pointer to store maximal number of permutations
1694  * (needed for freeing storage) */
1695  int* nmovedvars, /**< pointer to store number of vars affected
1696  * by symmetry (if usecompression) or NULL */
1697  SCIP_Bool* binvaraffected, /**< pointer to store whether a binary variable is affected by symmetry */
1698  SCIP_Bool* compressed, /**< pointer to store whether compression has been performed */
1699  SCIP_Real* log10groupsize, /**< pointer to store log10 of size of group */
1700  SCIP_Real* symcodetime, /**< pointer to store the time for symmetry code */
1701  SCIP_Bool* success /**< pointer to store whether symmetry computation was successful */
1702  )
1703 {
1704  SCIP_CONS** conss;
1705  SYM_GRAPH* graph;
1706  int nconsnodes = 0;
1707  int nvalnodes = 0;
1708  int nopnodes = 0;
1709  int nedges = 0;
1710  int nconss;
1711  int c;
1712 
1713  assert( scip != NULL );
1714  assert( permvars != NULL );
1715  assert( npermvars != NULL );
1716  assert( nbinpermvars != NULL );
1717  assert( perms != NULL );
1718  assert( nperms != NULL );
1719  assert( nmaxperms != NULL );
1720  assert( nmovedvars != NULL );
1721  assert( binvaraffected != NULL );
1722  assert( compressed != NULL );
1723  assert( log10groupsize != NULL );
1724  assert( symcodetime != NULL );
1725  assert( success != NULL );
1726 
1727  /* init pointers */
1728  *permvars = NULL;
1729  *npermvars = 0;
1730  *nbinpermvars = 0;
1731  *perms = NULL;
1732  *nperms = 0;
1733  *nmaxperms = 0;
1734  *nmovedvars = -1;
1735  *binvaraffected = FALSE;
1736  *compressed = FALSE;
1737  *log10groupsize = 0;
1738  *success = FALSE;
1739  *symcodetime = 0.0;
1740 
1741  /* check whether all constraints can provide symmetry information */
1742  if ( ! conshdlrsCanProvideSymInformation(scip, symtype) )
1743  return SCIP_OKAY;
1744 
1745  /* get symmetry detection graphs from constraints */
1746  conss = SCIPgetConss(scip);
1747  assert( conss != NULL );
1748 
1749  nconss = SCIPgetNConss(scip);
1750 
1751  /* exit if no constraints or no variables are available */
1752  if ( nconss == 0 || SCIPgetNVars(scip) == 0 )
1753  {
1754  *success = TRUE;
1755  return SCIP_OKAY;
1756  }
1757 
1758  /* get an estimate for the number of nodes and edges */
1759  SCIP_CALL( estimateSymgraphSize(scip, &nopnodes, &nvalnodes, &nconsnodes, &nedges) );
1760 
1761  /* create graph */
1762  SCIP_CALL( SCIPcreateSymgraph(scip, symtype, &graph, SCIPgetVars(scip), SCIPgetNVars(scip),
1763  nopnodes, nvalnodes, nconsnodes, nedges) );
1764 
1765  *success = TRUE;
1766  for (c = 0; c < nconss && *success; ++c)
1767  {
1768  if ( symtype == SYM_SYMTYPE_PERM )
1769  {
1770  SCIP_CALL( SCIPgetConsPermsymGraph(scip, conss[c], graph, success) );
1771  }
1772  else
1773  {
1774  assert( symtype == SYM_SYMTYPE_SIGNPERM );
1775  SCIP_CALL( SCIPgetConsSignedPermsymGraph(scip, conss[c], graph, success) );
1776  }
1777 
1778  /* terminate early if graph could not be returned */
1779  if ( ! *success )
1780  {
1781  SCIP_CALL( SCIPfreeSymgraph(scip, &graph) );
1782 
1783  return SCIP_OKAY;
1784  }
1785  }
1786 
1787  SCIP_CALL( SCIPcomputeSymgraphColors(scip, graph, fixedtype) );
1788 
1789  /* terminate early in case all variables are different */
1790  if ( (symtype == SYM_SYMTYPE_PERM && SCIPgetSymgraphNVarcolors(graph) == SCIPgetNVars(scip))
1791  || (symtype == SYM_SYMTYPE_SIGNPERM && SCIPgetSymgraphNVarcolors(graph) == 2 * SCIPgetNVars(scip)) )
1792  {
1793  SCIP_CALL( SCIPfreeSymgraph(scip, &graph) );
1794  return SCIP_OKAY;
1795  }
1796 
1797  /*
1798  * actually compute symmetries
1799  */
1800  SCIP_CALL( SYMcomputeSymmetryGenerators(scip, maxgenerators, graph, nperms, nmaxperms,
1801  perms, log10groupsize, symcodetime) );
1802 
1803  if ( checksymmetries && *nperms > 0 )
1804  {
1805  SCIP_CALL( checkSymmetriesAreSymmetries(scip, symtype, *perms, *nperms, SCIPgetNVars(scip), fixedtype) );
1806  }
1807 
1808  /* potentially store symmetries */
1809  if ( *nperms > 0 )
1810  {
1811  SCIP_VAR** vars;
1812  int nvars;
1813 
1814  nvars = SCIPgetNVars(scip);
1815  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &vars, SCIPgetVars(scip), nvars) ); /*lint !e666*/
1816 
1817  SCIP_CALL( setSymmetryData(scip, symtype, vars, nvars, SCIPgetNBinVars(scip), permvars, npermvars, nbinpermvars,
1818  permvardomaincenter, *perms, *nperms, nmovedvars, binvaraffected,
1819  compresssymmetries, compressthreshold, compressed) );
1820  }
1821 
1822  /* free symmetry graph */
1823  SCIP_CALL( SCIPfreeSymgraph(scip, &graph) );
1824 
1825  return SCIP_OKAY;
1826 }
1827 
1828 /** returns whether a symmetry is a non-standard permutation */
1829 static
1831  SCIP* scip, /**< SCIP instance */
1832  int* symmetry, /**< a symmetry encoded as a signed permutation */
1833  SCIP_VAR** vars, /**< array of variables the symmetry acts on */
1834  int nvars /**< number of variables in vars */
1835  )
1836 {
1837  int v;
1838 
1839  assert( symmetry != NULL );
1840  assert( vars != NULL );
1841  assert( nvars > 0 );
1842 
1843  for (v = 0; v < nvars; ++v)
1844  {
1845  /* the symmetry is signed */
1846  if ( symmetry[v] >= nvars )
1847  return TRUE;
1848 
1849  /* the domain of symmetric variables is different */
1850  if ( !SCIPisEQ(scip, SCIPvarGetLbLocal(vars[v]), SCIPvarGetLbLocal(vars[symmetry[v]]))
1851  || !SCIPisEQ(scip, SCIPvarGetUbLocal(vars[v]), SCIPvarGetUbLocal(vars[symmetry[v]])) )
1852  {
1853  assert( SCIPisEQ(scip, SCIPvarGetUbLocal(vars[v]) - SCIPvarGetLbLocal(vars[v]),
1854  SCIPvarGetUbLocal(vars[symmetry[v]]) - SCIPvarGetLbLocal(vars[symmetry[v]])) );
1855  return TRUE;
1856  }
1857  }
1858 
1859  return FALSE;
1860 }
1861 
1862 /** checks whether component contains non-standard permutations
1863  *
1864  * If all symmetries are standard permutations, stores them as such.
1865  */
1866 static
1868  SCIP* scip, /**< SCIP instance */
1869  SCIP_PROPDATA* propdata /**< propagator data */
1870  )
1871 {
1872  int* components;
1873  int* componentbegins;
1874  int ncomponents;
1875  int i;
1876  int c;
1877 
1878  assert( scip != NULL );
1879  assert( propdata != NULL );
1880  assert( propdata->ncomponents > 0 );
1881  assert( propdata->components != NULL );
1882  assert( propdata->componentbegins != NULL );
1883 
1884  components = propdata->components;
1885  componentbegins = propdata->componentbegins;
1886  ncomponents = propdata->ncomponents;
1887 
1888  SCIP_CALL( SCIPallocClearBlockMemoryArray(scip, &(propdata->componenthassignedperm), ncomponents) );
1889 
1890  /* stop if no non-standard permutations can exist */
1891  if ( (SYM_SYMTYPE) propdata->symtype == SYM_SYMTYPE_PERM )
1892  return SCIP_OKAY;
1893 
1894  /* for each component, check whether it has a non-standard permutation */
1895  for (c = 0; c < ncomponents; ++c)
1896  {
1897  for (i = componentbegins[c]; i < componentbegins[c + 1]; ++i)
1898  {
1899  if ( isNonstandardPerm(scip, propdata->perms[components[i]], propdata->permvars, propdata->npermvars) )
1900  {
1901  propdata->componenthassignedperm[c] = TRUE;
1902  break;
1903  }
1904  }
1905  }
1906 
1907  return SCIP_OKAY;
1908 }
1909 
1910 /** ensures that the symmetry components are already computed */
1911 static
1913  SCIP* scip, /**< SCIP instance */
1914  SCIP_PROPDATA* propdata /**< propagator data */
1915  )
1916 {
1917  assert( scip != NULL );
1918  assert( propdata != NULL );
1919 
1920  /* symmetries must have been determined */
1921  assert( propdata->nperms >= 0 );
1922 
1923  /* stop if already computed */
1924  if ( propdata->ncomponents >= 0 )
1925  return SCIP_OKAY;
1926 
1927  /* compute components */
1928  assert( propdata->ncomponents == -1 );
1929  assert( propdata->components == NULL );
1930  assert( propdata->componentbegins == NULL );
1931  assert( propdata->vartocomponent == NULL );
1932 
1933 #ifdef SCIP_OUTPUT_COMPONENT
1934  SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, " (%.1fs) component computation started\n", SCIPgetSolvingTime(scip));
1935 #endif
1936 
1937  SCIP_CALL( SCIPcomputeComponentsSym(scip, (SYM_SYMTYPE) propdata->symtype, propdata->perms, propdata->nperms,
1938  propdata->permvars, propdata->npermvars, FALSE, &propdata->components, &propdata->componentbegins,
1939  &propdata->vartocomponent, &propdata->componentblocked, &propdata->ncomponents) );
1940 
1941 #ifdef SCIP_OUTPUT_COMPONENT
1942  SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, " (%.1fs) component computation finished\n", SCIPgetSolvingTime(scip));
1943 #endif
1944 
1945  assert( propdata->components != NULL );
1946  assert( propdata->componentbegins != NULL );
1947  assert( propdata->ncomponents > 0 );
1948 
1949  /* structure of symmetries can be simplified if they are standard permutations */
1950  SCIP_CALL( checkComponentsForNonstandardPerms(scip, propdata) );
1951  assert( propdata->componenthassignedperm != NULL );
1952 
1953  return SCIP_OKAY;
1954 }
1955 
1956 
1957 /** ensures that permvarmap is initialized */
1958 static
1960  SCIP* scip, /**< SCIP instance */
1961  SCIP_PROPDATA* propdata /**< propagator data */
1962  )
1963 {
1964  int v;
1965 
1966  assert( scip != NULL );
1967  assert( propdata != NULL );
1968 
1969  /* symmetries must have been determined */
1970  assert( propdata->nperms >= 0 );
1971 
1972  /* stop if already computed */
1973  if ( propdata->permvarmap != NULL )
1974  return SCIP_OKAY;
1975 
1976  /* create hashmap for storing the indices of variables */
1977  SCIP_CALL( SCIPhashmapCreate(&propdata->permvarmap, SCIPblkmem(scip), propdata->npermvars) );
1978 
1979  /* insert variables into hashmap */
1980  for (v = 0; v < propdata->npermvars; ++v)
1981  {
1982  SCIP_CALL( SCIPhashmapInsertInt(propdata->permvarmap, propdata->permvars[v], v) );
1983  }
1984 
1985  return SCIP_OKAY;
1986 }
1987 
1988 
1989 /** ensures that permstrans is initialized */
1990 static
1992  SCIP* scip, /**< SCIP instance */
1993  SCIP_PROPDATA* propdata /**< propagator data */
1994  )
1995 {
1996  int v;
1997  int p;
1998 
1999  assert( scip != NULL );
2000  assert( propdata != NULL );
2001 
2002  /* symmetries must have been determined */
2003  assert( propdata->nperms >= 0 );
2004 
2005  /* stop if already computed */
2006  if ( propdata->permstrans != NULL )
2007  return SCIP_OKAY;
2008 
2009  /* transpose symmetries matrix here */
2010  assert( propdata->permstrans == NULL );
2011  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &propdata->permstrans, propdata->npermvars) );
2012  for (v = 0; v < propdata->npermvars; ++v)
2013  {
2014  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(propdata->permstrans[v]), propdata->nmaxperms) );
2015  for (p = 0; p < propdata->nperms; ++p)
2016  propdata->permstrans[v][p] = propdata->perms[p][v];
2017  }
2018 
2019  return SCIP_OKAY;
2020 }
2021 
2022 
2023 /** ensures that movedpermvarscounts is initialized */
2024 static
2026  SCIP* scip, /**< SCIP instance */
2027  SCIP_PROPDATA* propdata /**< propagator data */
2028  )
2029 {
2030  int v;
2031  int p;
2032 
2033  assert( scip != NULL );
2034  assert( propdata != NULL );
2035 
2036  /* symmetries must have been determined */
2037  assert( propdata->nperms >= 0 );
2038 
2039  /* stop if already computed */
2040  if ( propdata->nmovedpermvars >= 0 )
2041  return SCIP_OKAY;
2042  assert( propdata->nmovedpermvars == -1 );
2043 
2044  propdata->nmovedpermvars = 0;
2045  propdata->nmovedbinpermvars = 0;
2046  propdata->nmovedintpermvars = 0;
2047  propdata->nmovedimplintpermvars = 0;
2048  propdata->nmovedcontpermvars = 0;
2049 
2050  for (p = 0; p < propdata->nperms; ++p)
2051  {
2052  for (v = 0; v < propdata->npermvars; ++v)
2053  {
2054  if ( propdata->perms[p][v] != v )
2055  {
2056  ++propdata->nmovedpermvars;
2057 
2058  switch ( SCIPvarGetType(propdata->permvars[v]) )
2059  {
2060  case SCIP_VARTYPE_BINARY:
2061  ++propdata->nmovedbinpermvars;
2062  break;
2063  case SCIP_VARTYPE_INTEGER:
2064  ++propdata->nmovedintpermvars;
2065  break;
2066  case SCIP_VARTYPE_IMPLINT:
2067  ++propdata->nmovedimplintpermvars;
2068  break;
2070  ++propdata->nmovedcontpermvars;
2071  break;
2072  default:
2073  SCIPerrorMessage("Variable provided with unknown vartype\n");
2074  return SCIP_ERROR;
2075  }
2076  }
2077  }
2078  }
2079 
2080  return SCIP_OKAY;
2081 }
2082 
2083 
2084 /** returns whether any allowed symmetry handling method is effective for the problem instance */
2085 static
2087  SCIP* scip, /**< SCIP instance */
2088  SCIP_PROPDATA* propdata /**< propagator data */
2089  )
2090 {
2091  /* must always compute symmetry if it is enforced */
2092  if ( propdata->enforcecomputesymmetry )
2093  return TRUE;
2094 
2095  /* for dynamic symmetry handling or orbital reduction, branching must be possible */
2096  if ( propdata->usedynamicprop || ISORBITALREDUCTIONACTIVE(propdata->usesymmetry) )
2097  {
2098  /* @todo a proper test whether variables can be branched on or not */
2099  if ( SCIPgetNBinVars(scip) > 0 )
2100  return TRUE;
2101  if ( SCIPgetNIntVars(scip) > 0 )
2102  return TRUE;
2103  /* continuous variables can be branched on if nonlinear constraints exist */
2104  if ( ( SCIPgetNContVars(scip) > 0 || SCIPgetNImplVars(scip) > 0 )
2105  && SCIPconshdlrGetNActiveConss(propdata->conshdlr_nonlinear) > 0 )
2106  return TRUE;
2107  }
2108 
2109  /* for SST, matching leadervartypes */
2110  if ( ISSSTACTIVE(propdata->usesymmetry) )
2111  {
2112  if ( ISSSTBINACTIVE(propdata->sstleadervartype) && SCIPgetNBinVars(scip) > 0 ) /*lint !e641*/
2113  return TRUE;
2114  if ( ISSSTINTACTIVE(propdata->sstleadervartype) && SCIPgetNIntVars(scip) > 0 ) /*lint !e641*/
2115  return TRUE;
2116  if ( ISSSTIMPLINTACTIVE(propdata->sstleadervartype) && SCIPgetNImplVars(scip) > 0 ) /*lint !e641*/
2117  return TRUE;
2118  if ( ISSSTCONTACTIVE(propdata->sstleadervartype) && SCIPgetNContVars(scip) > 0 ) /*lint !e641*/
2119  return TRUE;
2120  }
2121 
2122  /* for static symmetry handling constraints, binary variables must be present */
2123  if ( ISSYMRETOPESACTIVE(propdata->usesymmetry) )
2124  {
2125  if ( SCIPgetNBinVars(scip) > 0 )
2126  return TRUE;
2127  }
2128 
2129  /* if all tests above fail, then the symmetry handling methods cannot achieve anything */
2130  return FALSE;
2131 }
2132 
2133 /** determines symmetry */
2134 static
2136  SCIP* scip, /**< SCIP instance */
2137  SCIP_PROPDATA* propdata, /**< propagator data */
2138  SYM_SPEC symspecrequire, /**< symmetry specification for which we need to compute symmetries */
2139  SYM_SPEC symspecrequirefixed /**< symmetry specification of variables which must be fixed by symmetries */
2140  )
2141 { /*lint --e{641}*/
2142  SCIP_Bool successful;
2143  SCIP_Real symcodetime = 0.0;
2144  int maxgenerators;
2145  unsigned int type = 0;
2146  int nvars;
2147  int i;
2148 
2149  assert( scip != NULL );
2150  assert( propdata != NULL );
2151  assert( propdata->usesymmetry >= 0 );
2152 
2153  /* do not compute symmetry if reoptimization is enabled */
2154  if ( SCIPisReoptEnabled(scip) )
2155  return SCIP_OKAY;
2156 
2157  /* do not compute symmetry if Benders decomposition enabled */
2158  if ( SCIPgetNActiveBenders(scip) > 0 )
2159  return SCIP_OKAY;
2160 
2161  /* skip symmetry computation if no graph automorphism code was linked */
2162  if ( ! SYMcanComputeSymmetry() )
2163  {
2165  " Deactivated symmetry handling methods, since SCIP was built without symmetry detector (SYM=none).\n");
2166 
2167  return SCIP_OKAY;
2168  }
2169 
2170  /* do not compute symmetry if there are active pricers */
2171  if ( SCIPgetNActivePricers(scip) > 0 )
2172  return SCIP_OKAY;
2173 
2174  /* avoid trivial cases */
2175  nvars = SCIPgetNVars(scip);
2176  if ( nvars <= 0 )
2177  return SCIP_OKAY;
2178 
2179  /* do not compute symmetry if we cannot handle it */
2180  if ( !testSymmetryComputationRequired(scip, propdata) )
2181  return SCIP_OKAY;
2182 
2183  /* determine symmetry specification */
2184  if ( SCIPgetNBinVars(scip) > 0 )
2185  type |= (int) SYM_SPEC_BINARY;
2186  if ( SCIPgetNIntVars(scip) > 0 )
2187  type |= (int) SYM_SPEC_INTEGER;
2188  /* count implicit integer variables as real variables, since we cannot currently handle integral variables well */
2189  if ( SCIPgetNContVars(scip) > 0 || SCIPgetNImplVars(scip) > 0 )
2190  type |= (int) SYM_SPEC_REAL;
2191 
2192  /* skip symmetry computation if required variables are not present */
2193  if ( ! (type & symspecrequire) )
2194  {
2196  " (%.1fs) symmetry computation skipped: type (bin %c, int %c, cont %c) does not match requirements (bin %c, int %c, cont %c).\n",
2197  SCIPgetSolvingTime(scip),
2198  SCIPgetNBinVars(scip) > 0 ? '+' : '-',
2199  SCIPgetNIntVars(scip) > 0 ? '+' : '-',
2200  SCIPgetNContVars(scip) + SCIPgetNImplVars(scip) > 0 ? '+' : '-',
2201  (symspecrequire & (int) SYM_SPEC_BINARY) != 0 ? '+' : '-',
2202  (symspecrequire & (int) SYM_SPEC_INTEGER) != 0 ? '+' : '-',
2203  (symspecrequire & (int) SYM_SPEC_REAL) != 0 ? '+' : '-');
2204 
2205  return SCIP_OKAY;
2206  }
2207 
2208  /* skip computation if symmetry has already been computed */
2209  if ( propdata->computedsymmetry )
2210  return SCIP_OKAY;
2211 
2212  assert( propdata->npermvars == 0 );
2213  assert( propdata->permvars == NULL );
2214  assert( propdata->nperms < 0 );
2215  assert( propdata->nmaxperms == 0 );
2216  assert( propdata->perms == NULL );
2217 
2218  /* output message */
2220  " (%.1fs) symmetry computation started: requiring (bin %c, int %c, cont %c), (fixed: bin %c, int %c, cont %c)\n",
2221  SCIPgetSolvingTime(scip),
2222  (symspecrequire & (int) SYM_SPEC_BINARY) != 0 ? '+' : '-',
2223  (symspecrequire & (int) SYM_SPEC_INTEGER) != 0 ? '+' : '-',
2224  (symspecrequire & (int) SYM_SPEC_REAL) != 0 ? '+' : '-',
2225  (symspecrequirefixed & (int) SYM_SPEC_BINARY) != 0 ? '+' : '-',
2226  (symspecrequirefixed & (int) SYM_SPEC_INTEGER) != 0 ? '+' : '-',
2227  (symspecrequirefixed & (int) SYM_SPEC_REAL) != 0 ? '+' : '-');
2228 
2229  /* output warning if we want to fix certain symmetry parts that we also want to compute */
2230  if ( symspecrequire & symspecrequirefixed )
2231  SCIPwarningMessage(scip, "Warning: some required symmetries must be fixed.\n");
2232 
2233  /* determine maximal number of generators depending on the number of variables */
2234  maxgenerators = propdata->maxgenerators;
2235  maxgenerators = MIN(maxgenerators, MAXGENNUMERATOR / nvars);
2236 
2237  /* actually compute (global) symmetry */
2238  SCIP_CALL( computeSymmetryGroup(scip, (SYM_SYMTYPE) propdata->symtype,
2239  propdata->compresssymmetries, propdata->compressthreshold,
2240  maxgenerators, symspecrequirefixed, propdata->checksymmetries, &propdata->permvars,
2241  &propdata->npermvars, &propdata->nbinpermvars, &propdata->permvardomaincenter,
2242  &propdata->perms, &propdata->nperms, &propdata->nmaxperms,
2243  &propdata->nmovedvars, &propdata->binvaraffected, &propdata->compressed,
2244  &propdata->log10groupsize, &symcodetime, &successful) );
2245 
2246  /* mark that we have computed the symmetry group */
2247  propdata->computedsymmetry = TRUE;
2248 
2249  /* store restart level */
2250  propdata->lastrestart = SCIPgetNRuns(scip);
2251 
2252  /* return if not successful */
2253  if ( ! successful )
2254  {
2255  assert( checkSymmetryDataFree(propdata) );
2256  SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, " (%.1fs) could not compute symmetry\n", SCIPgetSolvingTime(scip));
2257 
2258  return SCIP_OKAY;
2259  }
2260 
2261  /* return if no symmetries found */
2262  if ( propdata->nperms == 0 )
2263  {
2264  assert( checkSymmetryDataFree(propdata) );
2265  SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, " (%.1fs) no symmetry present (symcode time: %.2f)\n", SCIPgetSolvingTime(scip), symcodetime);
2266 
2267  return SCIP_OKAY;
2268  }
2269  assert( propdata->nperms > 0 );
2270  assert( propdata->npermvars > 0 );
2271 
2272  /* display statistics */
2273  SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, " (%.1fs) symmetry computation finished: %d generators found (max: ",
2274  SCIPgetSolvingTime(scip), propdata->nperms);
2275 
2276  /* display statistics: maximum number of generators */
2277  if ( maxgenerators == 0 )
2279  else
2280  SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, "%d", maxgenerators);
2281 
2282  /* display statistics: log10 group size, number of affected vars*/
2283  SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, ", log10 of symmetry group size: %.1f", propdata->log10groupsize);
2284 
2285  if ( propdata->displaynorbitvars )
2286  {
2287  if ( propdata->nmovedvars == -1 )
2288  {
2289  SCIP_CALL( SCIPdetermineNVarsAffectedSym(scip, propdata->perms, propdata->nperms, propdata->permvars,
2290  propdata->npermvars, &(propdata->nmovedvars)) );
2291  }
2292  SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, ", number of affected variables: %d)\n", propdata->nmovedvars);
2293  }
2294  SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, ") (symcode time: %.2f)\n", symcodetime);
2295 
2296  /* capture all variables while they are in the permvars array */
2297  for (i = 0; i < propdata->npermvars; ++i)
2298  {
2299  SCIP_CALL( SCIPcaptureVar(scip, propdata->permvars[i]) );
2300  }
2301 
2302  return SCIP_OKAY;
2303 }
2304 
2305 
2306 /*
2307  * Functions for symmetry constraints
2308  */
2309 
2310 
2311 /** Checks whether given set of 2-cycle permutations forms an orbitope and if so, builds the variable index matrix.
2312  *
2313  * If @p activevars == NULL, then the function assumes all permutations of the component are active and therefore all
2314  * moved vars are considered.
2315  *
2316  * We need to keep track of the number of generated columns, because we might not be able to detect all orbitopes.
2317  * For example (1,2), (2,3), (3,4), (3,5) defines the symmetric group on {1,2,3,4,5}, but the generators we expect
2318  * in our construction need shape (1,2), (2,3), (3,4), (4,5).
2319  *
2320  * @pre @p orbitopevaridx has to be an initialized 2D array of size @p ntwocycles x @p nperms
2321  * @pre @p columnorder has to be an initialized array of size nperms
2322  * @pre @p nusedelems has to be an initialized array of size npermvars
2323  */
2324 static
2326  SCIP* scip, /**< SCIP instance */
2327  SCIP_VAR** permvars, /**< array of all permutation variables */
2328  int npermvars, /**< number of permutation variables */
2329  int** perms, /**< array of all permutations of the symmetry group */
2330  int* activeperms, /**< indices of the relevant permutations in perms */
2331  int ntwocycles, /**< number of 2-cycles in the permutations */
2332  int nactiveperms, /**< number of active permutations */
2333  int** orbitopevaridx, /**< pointer to store variable index matrix */
2334  int* columnorder, /**< pointer to store column order */
2335  int* nusedelems, /**< pointer to store how often each element was used */
2336  int* nusedcols, /**< pointer to store number of columns used in orbitope (or NULL) */
2337  SCIP_Shortbool* rowisbinary, /**< pointer to store which rows are binary (or NULL) */
2338  SCIP_Bool* isorbitope, /**< buffer to store result */
2339  SCIP_Shortbool* activevars /**< bitset to store whether a variable is active (or NULL) */
2340  )
2341 { /*lint --e{571}*/
2342  SCIP_Bool* usedperm;
2343  SCIP_Bool foundperm = FALSE;
2344  int nusedperms = 0;
2345  int nfilledcols;
2346  int coltoextend;
2347  int ntestedperms = 0;
2348  int row = 0;
2349  int j;
2350 
2351  assert( scip != NULL );
2352  assert( permvars != NULL );
2353  assert( perms != NULL );
2354  assert( activeperms != NULL );
2355  assert( orbitopevaridx != NULL );
2356  assert( columnorder != NULL );
2357  assert( nusedelems != NULL );
2358  assert( isorbitope != NULL );
2359  assert( nactiveperms > 0 );
2360  assert( ntwocycles > 0 );
2361  assert( npermvars > 0 );
2362  assert( activevars == NULL || (0 <= nactiveperms && nactiveperms < npermvars) );
2363 
2364  *isorbitope = TRUE;
2365  if ( nusedcols != NULL )
2366  *nusedcols = 0;
2367 
2368  /* whether a permutation was considered to contribute to orbitope */
2369  SCIP_CALL( SCIPallocClearBufferArray(scip, &usedperm, nactiveperms) );
2370 
2371  /* fill first two columns of orbitopevaridx matrix */
2372 
2373  /* look for the first active permutation which moves an active variable */
2374  while ( ! foundperm )
2375  {
2376  int permidx;
2377 
2378  assert( ntestedperms < nactiveperms );
2379 
2380  permidx = activeperms[ntestedperms];
2381 
2382  for (j = 0; j < npermvars; ++j)
2383  {
2384  if ( activevars != NULL && ! activevars[j] )
2385  continue;
2386 
2387  assert( activevars == NULL || activevars[perms[permidx][j]] );
2388 
2389  /* avoid adding the same 2-cycle twice */
2390  if ( perms[permidx][j] > j )
2391  {
2392  assert( SCIPvarIsBinary(permvars[j]) == SCIPvarIsBinary(permvars[perms[permidx][j]]) );
2393 
2394  if ( rowisbinary != NULL && SCIPvarIsBinary(permvars[j]) )
2395  rowisbinary[row] = TRUE;
2396 
2397  orbitopevaridx[row][0] = j;
2398  orbitopevaridx[row++][1] = perms[permidx][j];
2399  ++(nusedelems[j]);
2400  ++(nusedelems[perms[permidx][j]]);
2401 
2402  foundperm = TRUE;
2403  }
2404 
2405  if ( row == ntwocycles )
2406  break;
2407  }
2408 
2409  ++ntestedperms;
2410  }
2411 
2412  /* in the subgroup case it might happen that a generator has less than ntwocycles many 2-cyles */
2413  if ( row != ntwocycles )
2414  {
2415  *isorbitope = FALSE;
2416  SCIPfreeBufferArray(scip, &usedperm);
2417  return SCIP_OKAY;
2418  }
2419 
2420  usedperm[ntestedperms - 1] = TRUE;
2421  ++nusedperms;
2422  columnorder[0] = 0;
2423  columnorder[1] = 1;
2424  nfilledcols = 2;
2425 
2426  /* extend orbitopevaridx matrix to the left, i.e., iteratively find new permutations that
2427  * intersect the last added left column in each row in exactly one entry, starting with
2428  * column 0 */
2429  coltoextend = 0;
2430  for (j = ntestedperms; j < nactiveperms; ++j)
2431  { /* lint --e{850} */
2432  SCIP_Bool success = FALSE;
2433  SCIP_Bool infeasible = FALSE;
2434 
2435  if ( nusedperms == nactiveperms )
2436  break;
2437 
2438  if ( usedperm[j] )
2439  continue;
2440 
2441  SCIP_CALL( SCIPextendSubOrbitope(orbitopevaridx, ntwocycles, nfilledcols, coltoextend,
2442  perms[activeperms[j]], TRUE, &nusedelems, permvars, NULL, &success, &infeasible) );
2443 
2444  if ( infeasible )
2445  {
2446  *isorbitope = FALSE;
2447  break;
2448  }
2449  else if ( success )
2450  {
2451  usedperm[j] = TRUE;
2452  ++nusedperms;
2453  coltoextend = nfilledcols;
2454  columnorder[nfilledcols++] = -1; /* mark column to be filled from the left */
2455  j = 0; /*lint !e850*/ /* reset j since previous permutations can now intersect with the latest added column */
2456  }
2457  }
2458 
2459  if ( ! *isorbitope ) /*lint !e850*/
2460  {
2461  SCIPfreeBufferArray(scip, &usedperm);
2462  return SCIP_OKAY;
2463  }
2464 
2465  coltoextend = 1;
2466  for (j = ntestedperms; j < nactiveperms; ++j)
2467  { /*lint --e(850)*/
2468  SCIP_Bool success = FALSE;
2469  SCIP_Bool infeasible = FALSE;
2470 
2471  if ( nusedperms == nactiveperms )
2472  break;
2473 
2474  if ( usedperm[j] )
2475  continue;
2476 
2477  SCIP_CALL( SCIPextendSubOrbitope(orbitopevaridx, ntwocycles, nfilledcols, coltoextend,
2478  perms[activeperms[j]], FALSE, &nusedelems, permvars, NULL, &success, &infeasible) );
2479 
2480  if ( infeasible )
2481  {
2482  *isorbitope = FALSE;
2483  break;
2484  }
2485  else if ( success )
2486  {
2487  usedperm[j] = TRUE;
2488  ++nusedperms;
2489  coltoextend = nfilledcols;
2490  columnorder[nfilledcols] = 1; /* mark column to be filled from the right */
2491  ++nfilledcols;
2492  j = 0; /*lint !e850*/ /* reset j since previous permutations can now intersect with the latest added column */
2493  }
2494  }
2495 
2496  if ( activevars == NULL && nusedperms < nactiveperms ) /*lint !e850*/
2497  *isorbitope = FALSE;
2498 
2499  if ( nusedcols != NULL )
2500  *nusedcols = nfilledcols;
2501  assert( ! *isorbitope || activevars == NULL || nusedperms < nfilledcols );
2502 
2503  SCIPfreeBufferArray(scip, &usedperm);
2504 
2505  return SCIP_OKAY;
2506 }
2507 
2508 /** choose an order in which the generators should be added for subgroup detection */
2509 static
2511  SCIP* scip, /**< SCIP instance */
2512  SCIP_PROPDATA* propdata, /**< pointer to data of symmetry propagator */
2513  int compidx, /**< index of component */
2514  int** genorder, /**< (initialized) buffer to store the resulting order of generator */
2515  int* ntwocycleperms /**< pointer to store the number of 2-cycle permutations in component compidx */
2516  )
2517 {
2518  int** perms;
2519  int* components;
2520  int* componentbegins;
2521  int* ntwocycles;
2522  int npermvars;
2523  int npermsincomp;
2524  int i;
2525 
2526  assert( scip != NULL );
2527  assert( propdata != NULL );
2528  assert( compidx >= 0 );
2529  assert( compidx < propdata->ncomponents );
2530  assert( genorder != NULL );
2531  assert( *genorder != NULL );
2532  assert( ntwocycleperms != NULL );
2533  assert( propdata->computedsymmetry );
2534  assert( propdata->nperms > 0 );
2535  assert( propdata->perms != NULL );
2536  assert( propdata->npermvars > 0 );
2537  assert( propdata->ncomponents > 0 );
2538  assert( propdata->components != NULL );
2539  assert( propdata->componentbegins != NULL );
2540 
2541  perms = propdata->perms;
2542  npermvars = propdata->npermvars;
2543  components = propdata->components;
2544  componentbegins = propdata->componentbegins;
2545  npermsincomp = componentbegins[compidx + 1] - componentbegins[compidx];
2546  *ntwocycleperms = npermsincomp;
2547 
2548  SCIP_CALL( SCIPallocBufferArray(scip, &ntwocycles, npermsincomp) );
2549 
2550  for (i = 0; i < npermsincomp; ++i)
2551  {
2552  int* perm;
2553  int nbincycles;
2554 
2555  perm = perms[components[componentbegins[compidx] + i]];
2556 
2557  SCIP_CALL( SCIPisInvolutionPerm(perm, propdata->permvars, npermvars, &(ntwocycles[i]), &nbincycles, FALSE) );
2558 
2559  /* we skip permutations which do not purely consist of 2-cycles */
2560  if ( ntwocycles[i] == 0 )
2561  {
2562  /* we change the number of two cycles for this perm so that it will be sorted to the end */
2563  if ( propdata->preferlessrows )
2564  ntwocycles[i] = npermvars;
2565  else
2566  ntwocycles[i] = 0;
2567  --(*ntwocycleperms);
2568  }
2569  else if ( ! propdata->preferlessrows )
2570  ntwocycles[i] = - ntwocycles[i];
2571  }
2572 
2573  SCIPsortIntInt(ntwocycles, *genorder, npermsincomp);
2574 
2575  SCIPfreeBufferArray(scip, &ntwocycles);
2576 
2577  return SCIP_OKAY;
2578 }
2579 
2580 
2581 /** builds the graph for symmetric subgroup detection from the given permutation of generators
2582  *
2583  * After execution, @p graphcomponents contains all permvars sorted by their color and component,
2584  * @p graphcompbegins points to the indices where new components in @p graphcomponents start and
2585  * @p compcolorbegins points to the indices where new colors in @p graphcompbegins start.
2586 */
2587 static
2589  SCIP* scip, /**< SCIP instance */
2590  SCIP_PROPDATA* propdata, /**< pointer to data of symmetry propagator */
2591  int* genorder, /**< order in which the generators should be considered */
2592  int ntwocycleperms, /**< number of 2-cycle permutations in this component */
2593  int compidx, /**< index of the component */
2594  int** graphcomponents, /**< buffer to store the components of the graph (ordered var indices) */
2595  int** graphcompbegins, /**< buffer to store the indices of each new graph component */
2596  int** compcolorbegins, /**< buffer to store at which indices a new color begins */
2597  int* ngraphcomponents, /**< pointer to store the number of graph components */
2598  int* ncompcolors, /**< pointer to store the number of different colors */
2599  int** usedperms, /**< buffer to store the indices of permutations that were used */
2600  int* nusedperms, /**< pointer to store the number of used permutations in the graph */
2601  int usedpermssize, /**< initial size of usedperms */
2602  SCIP_Shortbool* permused /**< initialized buffer to store which permutations have been used
2603  * (identified by index in component) */
2604  )
2605 {
2606  SCIP_DISJOINTSET* vartocomponent;
2607  SCIP_DISJOINTSET* comptocolor;
2608  int** perms;
2609  int* components;
2610  int* componentbegins;
2611  int* componentslastperm;
2612  SYM_SORTGRAPHCOMPVARS graphcompvartype;
2613  int npermvars;
2614  int nextcolor;
2615  int nextcomp;
2616  int j;
2617  int k;
2618 
2619  assert( scip != NULL );
2620  assert( propdata != NULL );
2621  assert( graphcomponents != NULL );
2622  assert( graphcompbegins != NULL );
2623  assert( compcolorbegins != NULL );
2624  assert( ngraphcomponents != NULL );
2625  assert( ncompcolors != NULL );
2626  assert( genorder != NULL );
2627  assert( usedperms != NULL );
2628  assert( nusedperms != NULL );
2629  assert( usedpermssize > 0 );
2630  assert( permused != NULL );
2631  assert( ntwocycleperms >= 0 );
2632  assert( compidx >= 0 );
2633  assert( compidx < propdata->ncomponents );
2634  assert( propdata->computedsymmetry );
2635  assert( propdata->nperms > 0 );
2636  assert( propdata->perms != NULL );
2637  assert( propdata->npermvars > 0 );
2638  assert( propdata->ncomponents > 0 );
2639  assert( propdata->components != NULL );
2640  assert( propdata->componentbegins != NULL );
2641  assert( ! propdata->componentblocked[compidx] );
2642 
2643  perms = propdata->perms;
2644  npermvars = propdata->npermvars;
2645  components = propdata->components;
2646  componentbegins = propdata->componentbegins;
2647  *nusedperms = 0;
2648 
2649  assert( ntwocycleperms <= componentbegins[compidx + 1] - componentbegins[compidx] );
2650 
2651  SCIP_CALL( SCIPcreateDisjointset(scip, &vartocomponent, npermvars) );
2652  SCIP_CALL( SCIPcreateDisjointset(scip, &comptocolor, npermvars) );
2653  SCIP_CALL( SCIPallocBufferArray( scip, &componentslastperm, npermvars) );
2654 
2655  for (k = 0; k < npermvars; ++k)
2656  componentslastperm[k] = -1;
2657 
2658  for (j = 0; j < ntwocycleperms; ++j)
2659  {
2660  int* perm;
2661  int firstcolor = -1;
2662 
2663  /* use given order of generators */
2664  perm = perms[components[componentbegins[compidx] + genorder[j]]];
2665  assert( perm != NULL );
2666 
2667  /* iteratively handle each swap of perm until an invalid one is found or all edges have been added */
2668  for (k = 0; k < npermvars; ++k)
2669  {
2670  int comp1;
2671  int comp2;
2672  int color1;
2673  int color2;
2674  int img;
2675 
2676  img = perm[k];
2677  assert( perm[img] == k );
2678 
2679  if ( img <= k )
2680  continue;
2681 
2682  comp1 = SCIPdisjointsetFind(vartocomponent, k);
2683  comp2 = SCIPdisjointsetFind(vartocomponent, img);
2684 
2685  if ( comp1 == comp2 )
2686  {
2687  /* another permutation has already merged these variables into one component; store its color */
2688  if ( firstcolor < 0 )
2689  {
2690  assert( SCIPdisjointsetFind(comptocolor, comp1) == SCIPdisjointsetFind(comptocolor, comp2) );
2691  firstcolor = SCIPdisjointsetFind(comptocolor, comp1);
2692  }
2693  componentslastperm[comp1] = j;
2694  continue;
2695  }
2696 
2697  /* if it is the second time that the component is used for this generator,
2698  * it is not guaranteed that the group acts like the symmetric group, so skip it
2699  */
2700  if ( componentslastperm[comp1] == j || componentslastperm[comp2] == j )
2701  break;
2702 
2703  color1 = SCIPdisjointsetFind(comptocolor, comp1);
2704  color2 = SCIPdisjointsetFind(comptocolor, comp2);
2705 
2706  /* a generator is not allowed to connect two components of the same color, since they depend on each other */
2707  if ( color1 == color2 )
2708  break;
2709 
2710  componentslastperm[comp1] = j;
2711  componentslastperm[comp2] = j;
2712 
2713  if ( firstcolor < 0 )
2714  firstcolor = color1;
2715  }
2716 
2717  /* if the generator is invalid, delete the newly added edges, go to next generator */
2718  if ( k < npermvars )
2719  continue;
2720 
2721  /* if the generator only acts on already existing components, we don't have to store it */
2722  if ( firstcolor == -1 )
2723  continue;
2724 
2725  /* check whether we need to resize */
2726  if ( *nusedperms >= usedpermssize )
2727  {
2728  int newsize = SCIPcalcMemGrowSize(scip, (*nusedperms) + 1);
2729  assert( newsize > usedpermssize );
2730 
2731  SCIP_CALL( SCIPreallocBufferArray(scip, usedperms, newsize) );
2732 
2733  usedpermssize = newsize;
2734  }
2735 
2736  (*usedperms)[*nusedperms] = components[componentbegins[compidx] + genorder[j]];
2737  ++(*nusedperms);
2738  permused[genorder[j]] = TRUE;
2739 
2740  /* if the generator can be added, update the datastructures for graph components and colors */
2741  for (k = 0; k < npermvars; ++k)
2742  {
2743  int comp1;
2744  int comp2;
2745  int color1;
2746  int color2;
2747  int img;
2748 
2749  img = perm[k];
2750  assert( perm[img] == k );
2751 
2752  if ( img <= k )
2753  continue;
2754 
2755  comp1 = SCIPdisjointsetFind(vartocomponent, k);
2756  comp2 = SCIPdisjointsetFind(vartocomponent, img);
2757 
2758  /* components and colors don't have to be updated if the components are the same */
2759  if ( comp1 == comp2 )
2760  continue;
2761 
2762  color1 = SCIPdisjointsetFind(comptocolor, comp1);
2763  color2 = SCIPdisjointsetFind(comptocolor, comp2);
2764 
2765  if ( color1 != color2 )
2766  {
2767  SCIPdisjointsetUnion(comptocolor, firstcolor, color1, TRUE);
2768  SCIPdisjointsetUnion(comptocolor, firstcolor, color2, TRUE);
2769  }
2770 
2771  SCIPdisjointsetUnion(vartocomponent, comp1, comp2, FALSE);
2772 
2773  assert( SCIPdisjointsetFind(vartocomponent, k) == SCIPdisjointsetFind(vartocomponent, img) );
2774  assert( SCIPdisjointsetFind(comptocolor, SCIPdisjointsetFind(vartocomponent, k)) == firstcolor );
2775  assert( SCIPdisjointsetFind(comptocolor, SCIPdisjointsetFind(vartocomponent, img)) == firstcolor );
2776  }
2777  }
2778 
2779  SCIP_CALL( SCIPallocBlockMemoryArray(scip, graphcomponents, npermvars) );
2780  SCIP_CALL( SCIPallocBufferArray(scip, &(graphcompvartype.components), npermvars) );
2781  SCIP_CALL( SCIPallocBufferArray(scip, &(graphcompvartype.colors), npermvars) );
2782 
2783  /*
2784  * At this point, we have built the colored graph. Now we transform the information in the
2785  * disjoint sets to the arrays graphcomponents, graphcompbegins, and compcolorbegins (see above).
2786  */
2787 
2788  /* build the struct graphcompvartype which is used to sort the graphcomponents array */
2789  for (j = 0; j < npermvars; ++j)
2790  {
2791  int comp;
2792 
2793  comp = SCIPdisjointsetFind(vartocomponent, j);
2794 
2795  graphcompvartype.components[j] = comp;
2796  graphcompvartype.colors[j] = SCIPdisjointsetFind(comptocolor, comp);
2797 
2798  (*graphcomponents)[j] = j;
2799  }
2800 
2801  /* sort graphcomponents first by color, then by component */
2802  SCIPsort(*graphcomponents, SYMsortGraphCompVars, (void*) &graphcompvartype, npermvars);
2803 
2804  *ngraphcomponents = SCIPdisjointsetGetComponentCount(vartocomponent);
2805  *ncompcolors = SCIPdisjointsetGetComponentCount(comptocolor);
2806  SCIP_CALL( SCIPallocBlockMemoryArray(scip, graphcompbegins, (*ngraphcomponents) + 1) );
2807  SCIP_CALL( SCIPallocBlockMemoryArray(scip, compcolorbegins, (*ncompcolors) + 1) );
2808 
2809  nextcolor = 1;
2810  nextcomp = 1;
2811  (*graphcompbegins)[0] = 0;
2812  (*compcolorbegins)[0] = 0;
2813 
2814  /* find the starting indices of new components and new colors */
2815  for (j = 1; j < npermvars; ++j)
2816  {
2817  int idx1;
2818  int idx2;
2819 
2820  idx1 = (*graphcomponents)[j];
2821  idx2 = (*graphcomponents)[j-1];
2822 
2823  assert( graphcompvartype.colors[idx1] >= graphcompvartype.colors[idx2] );
2824 
2825  if ( graphcompvartype.components[idx1] != graphcompvartype.components[idx2] )
2826  {
2827  (*graphcompbegins)[nextcomp] = j;
2828 
2829  if ( graphcompvartype.colors[idx1] > graphcompvartype.colors[idx2] )
2830  {
2831  (*compcolorbegins)[nextcolor] = nextcomp;
2832  ++nextcolor;
2833  }
2834 
2835  ++nextcomp;
2836  }
2837  }
2838  assert( nextcomp == *ngraphcomponents );
2839  assert( nextcolor == *ncompcolors );
2840 
2841  (*compcolorbegins)[nextcolor] = *ngraphcomponents;
2842  (*graphcompbegins)[nextcomp] = npermvars;
2843 
2844  SCIPfreeBufferArray(scip, &(graphcompvartype.colors));
2845  SCIPfreeBufferArray(scip, &(graphcompvartype.components));
2846  SCIPfreeBufferArray(scip, &componentslastperm);
2847  SCIPfreeDisjointset(scip, &comptocolor);
2848  SCIPfreeDisjointset(scip, &vartocomponent);
2849 
2850  return SCIP_OKAY;
2851 }
2852 
2853 /** adds an orbitope constraint for a suitable color of the subgroup graph */
2854 static
2856  SCIP* scip, /**< SCIP instance */
2857  SCIP_PROPDATA* propdata, /**< pointer to data of symmetry propagator */
2858  int* usedperms, /**< array of the permutations that build the orbitope */
2859  int nusedperms, /**< number of permutations in usedperms */
2860  int* compcolorbegins, /**< array indicating where a new graphcolor begins */
2861  int* graphcompbegins, /**< array indicating where a new graphcomponent begins */
2862  int* graphcomponents, /**< array of all variable indices sorted by color and comp */
2863  int graphcoloridx, /**< index of the graph color */
2864  int nrows, /**< number of rows in the orbitope */
2865  int ncols, /**< number of columns in the orbitope */
2866  int* firstvaridx, /**< buffer to store the index of the largest variable (or NULL) */
2867  int* compidxfirstrow, /**< buffer to store the comp index for the first row (or NULL) */
2868  int** lexorder, /**< pointer to array storing lexicographic order defined by sub orbitopes */
2869  int* nvarslexorder, /**< number of variables in lexicographic order */
2870  int* maxnvarslexorder, /**< maximum number of variables in lexicographic order */
2871  SCIP_Bool mayinteract, /**< whether orbitope's symmetries might interact with other symmetries */
2872  SCIP_Bool* success /**< whether the orbitope could be added */
2873  )
2874 { /*lint --e{571}*/
2875  char name[SCIP_MAXSTRLEN];
2876  SCIP_VAR*** orbitopevarmatrix;
2877  SCIP_Shortbool* activevars;
2878  int** orbitopevaridx;
2879  int* columnorder;
2880  int* nusedelems;
2881  SCIP_CONS* cons;
2882  SCIP_Bool isorbitope;
2883  SCIP_Bool infeasible = FALSE;
2884 #ifndef NDEBUG
2885  int nactivevars = 0;
2886 #endif
2887  int ngencols = 0;
2888  int k;
2889 
2890  assert( scip != NULL );
2891  assert( propdata != NULL );
2892  assert( usedperms != NULL );
2893  assert( compcolorbegins != NULL );
2894  assert( graphcompbegins != NULL );
2895  assert( graphcomponents != NULL );
2896  assert( nusedperms > 0 );
2897  assert( nrows > 0 );
2898  assert( ncols > 0 );
2899  assert( lexorder != NULL );
2900  assert( nvarslexorder != NULL );
2901  assert( maxnvarslexorder != NULL );
2902 
2903  *success = FALSE;
2904 
2905  /* create hashset to mark variables */
2906  SCIP_CALL( SCIPallocClearBufferArray(scip, &activevars, propdata->npermvars) );
2907 
2908  /* orbitope matrix for indices of variables in permvars array */
2909  SCIP_CALL( SCIPallocBufferArray(scip, &orbitopevaridx, nrows) );
2910  for (k = 0; k < nrows; ++k)
2911  {
2912  SCIP_CALL( SCIPallocBufferArray(scip, &orbitopevaridx[k], ncols) ); /*lint !e866*/
2913  }
2914 
2915  /* order of columns of orbitopevaridx */
2916  SCIP_CALL( SCIPallocBufferArray(scip, &columnorder, ncols) );
2917  for (k = 0; k < ncols; ++k)
2918  columnorder[k] = ncols + 1;
2919 
2920  /* count how often an element was used in the potential orbitope */
2921  SCIP_CALL( SCIPallocClearBufferArray(scip, &nusedelems, propdata->npermvars) );
2922 
2923  /* mark variables in this subgroup orbitope */
2924  for (k = compcolorbegins[graphcoloridx]; k < compcolorbegins[graphcoloridx+1]; ++k)
2925  {
2926  SCIP_VAR* firstvar;
2927  int compstart;
2928  int l;
2929 
2930  compstart = graphcompbegins[k];
2931  firstvar = propdata->permvars[graphcomponents[compstart]];
2932 
2933  if ( ! SCIPvarIsBinary(firstvar) )
2934  continue;
2935 
2936  for (l = 0; l < ncols; ++l)
2937  {
2938  int varidx;
2939 
2940  varidx = graphcomponents[compstart + l];
2941  assert( ! activevars[varidx] );
2942 
2943  activevars[varidx] = TRUE;
2944 #ifndef NDEBUG
2945  ++nactivevars;
2946 #endif
2947  }
2948  }
2949  assert( nactivevars == nrows * ncols );
2950 
2951  /* build the variable index matrix for the orbitope
2952  *
2953  * It is possible that we find an orbitope, but not using all possible columns. For example
2954  * (1,2), (2,3), (3,4), (3,5) defines the symmetric group on {1,2,3,4,5}, but the generators
2955  * we expect in our construction need shape (1,2), (2,3), (3,4), (4,5). For this reason,
2956  * we need to store how many columns have been generated.
2957  *
2958  * @todo ensure compatibility with more general generators
2959  */
2960  SCIP_CALL( checkTwoCyclePermsAreOrbitope(scip, propdata->permvars, propdata->npermvars,
2961  propdata->perms, usedperms, nrows, nusedperms, orbitopevaridx, columnorder,
2962  nusedelems, &ngencols, NULL, &isorbitope, activevars) );
2963 
2964  /* it might happen that we cannot detect the orbitope if it is generated by permutations with different
2965  * number of 2-cycles.
2966  */
2967  if ( ! isorbitope )
2968  {
2969  SCIPfreeBufferArray(scip, &nusedelems);
2970  SCIPfreeBufferArray(scip, &columnorder);
2971  for (k = nrows - 1; k >= 0; --k)
2972  {
2973  SCIPfreeBufferArray(scip, &orbitopevaridx[k]);
2974  }
2975  SCIPfreeBufferArray(scip, &orbitopevaridx);
2976  SCIPfreeBufferArray(scip, &activevars);
2977 
2978  return SCIP_OKAY;
2979  }
2980 
2981  /* There are three possibilities for the structure of columnorder:
2982  * 1) [0, 1, -1, -1, ..., -1]
2983  * 2) [0, 1, 1, 1, ..., 1]
2984  * 3) [0, 1, -1, -1, ...., -1, 1, 1, ..., 1]
2985  *
2986  * The '1'-columns will be added to the matrix first and in the last 2
2987  * cases the method starts from the right. So to store the variable index
2988  * that will be in the upper-left corner, we need either the entryin the
2989  * second column (case 1) or the entry in the last column (cases 2 and 3).
2990  */
2991  if ( firstvaridx != NULL )
2992  {
2993  if ( columnorder[ngencols-1] > -1 )
2994  *firstvaridx = orbitopevaridx[0][ngencols-1];
2995  else
2996  *firstvaridx = orbitopevaridx[0][1];
2997  }
2998 
2999  /* find corresponding graphcomponent of first variable (needed for weak sbcs) */
3000  if ( compidxfirstrow != NULL && firstvaridx != NULL )
3001  {
3002  *compidxfirstrow = -1;
3003 
3004  for (k = compcolorbegins[graphcoloridx]; k < compcolorbegins[graphcoloridx+1] && (*compidxfirstrow) < 0; ++k)
3005  {
3006  SCIP_VAR* firstvar;
3007  int compstart;
3008  int l;
3009 
3010  compstart = graphcompbegins[k];
3011  firstvar = propdata->permvars[graphcomponents[compstart]];
3012 
3013  if ( ! SCIPvarIsBinary(firstvar) )
3014  continue;
3015 
3016  /* iterate over all columns (elements in orbit), because we cannot see from ngencols which columns
3017  * have been left out
3018  */
3019  for (l = 0; l < ncols; ++l)
3020  {
3021  if ( graphcomponents[compstart + l] == *firstvaridx )
3022  {
3023  *compidxfirstrow = k;
3024  break;
3025  }
3026  }
3027  }
3028  assert( *compidxfirstrow > -1 );
3029  }
3030 
3031  /* prepare orbitope variable matrix */
3032  SCIP_CALL( SCIPallocBufferArray(scip, &orbitopevarmatrix, nrows) );
3033  for (k = 0; k < nrows; ++k)
3034  {
3035  SCIP_CALL( SCIPallocBufferArray(scip, &orbitopevarmatrix[k], ngencols) );
3036  }
3037 
3038  /* build the matrix containing the actual variables of the orbitope */
3039  SCIP_CALL( SCIPgenerateOrbitopeVarsMatrix(scip, &orbitopevarmatrix, nrows, ngencols,
3040  propdata->permvars, propdata->npermvars, orbitopevaridx, columnorder,
3041  nusedelems, NULL, &infeasible, TRUE, lexorder, nvarslexorder, maxnvarslexorder) );
3042 
3043  assert( ! infeasible );
3044  assert( firstvaridx == NULL || propdata->permvars[*firstvaridx] == orbitopevarmatrix[0][0] );
3045 
3046  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "suborbitope_%d_%d", graphcoloridx, propdata->norbitopes);
3047 
3048  SCIP_CALL( SCIPcreateConsOrbitope(scip, &cons, name, orbitopevarmatrix,
3049  SCIP_ORBITOPETYPE_FULL, nrows, ngencols, FALSE, mayinteract, FALSE, FALSE, propdata->conssaddlp,
3050  TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
3051 
3052  SCIP_CALL( SCIPaddCons(scip, cons) );
3053  *success = TRUE;
3054 
3055  /* do not release constraint here - will be done later */
3056  SCIP_CALL( ensureDynamicConsArrayAllocatedAndSufficientlyLarge(scip, &propdata->genorbconss,
3057  &propdata->genorbconsssize, propdata->ngenorbconss + 1) );
3058  propdata->genorbconss[propdata->ngenorbconss++] = cons;
3059  ++propdata->norbitopes;
3060 
3061  for (k = nrows - 1; k >= 0; --k)
3062  SCIPfreeBufferArray(scip, &orbitopevarmatrix[k]);
3063  SCIPfreeBufferArray(scip, &orbitopevarmatrix);
3064  SCIPfreeBufferArray(scip, &nusedelems);
3065  SCIPfreeBufferArray(scip, &columnorder);
3066  for (k = nrows - 1; k >= 0; --k)
3067  SCIPfreeBufferArray(scip, &orbitopevaridx[k]);
3068  SCIPfreeBufferArray(scip, &orbitopevaridx);
3069  SCIPfreeBufferArray(scip, &activevars);
3070 
3071  return SCIP_OKAY;
3072 }
3073 
3074 /** adds strong SBCs for a suitable color of the subgroup graph */
3075 static
3077  SCIP* scip, /**< SCIP instance */
3078  SCIP_PROPDATA* propdata, /**< pointer to data of symmetry propagator */
3079  int* graphcompbegins, /**< array indicating where a new graphcomponent begins */
3080  int* graphcomponents, /**< array of all variable indices sorted by color and comp */
3081  int graphcompidx, /**< index of the graph component */
3082  SCIP_Bool storelexorder, /**< whether the lexicographic order induced by the orbitope shall be stored */
3083  int** lexorder, /**< pointer to array storing lexicographic order defined by sub orbitopes */
3084  int* nvarsorder, /**< number of variables in lexicographic order */
3085  int* maxnvarsorder /**< maximum number of variables in lexicographic order */
3086  )
3087 {
3088  int k;
3089 
3090  assert( scip != NULL );
3091  assert( propdata != NULL );
3092  assert( graphcompbegins != NULL );
3093  assert( graphcomponents != NULL );
3094  assert( graphcompidx >= 0 );
3095  assert( ! storelexorder || lexorder != NULL );
3096  assert( ! storelexorder || nvarsorder != NULL );
3097  assert( ! storelexorder || maxnvarsorder != NULL );
3098 
3099  /* possibly store lexicographic order defined by strong SBCs */
3100  if ( storelexorder )
3101  {
3102  if ( *maxnvarsorder == 0 )
3103  {
3104  *maxnvarsorder = graphcompbegins[graphcompidx + 1] - graphcompbegins[graphcompidx + 1];
3105  *nvarsorder = 0;
3106 
3107  SCIP_CALL( SCIPallocBlockMemoryArray(scip, lexorder, *maxnvarsorder) );
3108  }
3109  else
3110  {
3111  assert( *nvarsorder == *maxnvarsorder );
3112 
3113  *maxnvarsorder += graphcompbegins[graphcompidx + 1] - graphcompbegins[graphcompidx + 1];
3114 
3115  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, lexorder, *nvarsorder, *maxnvarsorder) );
3116  }
3117 
3118  (*lexorder)[*nvarsorder++] = graphcomponents[graphcompbegins[graphcompidx]];
3119  }
3120 
3121  /* add strong SBCs (lex-max order) for chosen graph component */
3122  for (k = graphcompbegins[graphcompidx]+1; k < graphcompbegins[graphcompidx+1]; ++k)
3123  {
3124  char name[SCIP_MAXSTRLEN];
3125  SCIP_CONS* cons;
3126  SCIP_VAR* vars[2];
3127  SCIP_Real vals[2] = {1, -1};
3128 
3129  vars[0] = propdata->permvars[graphcomponents[k-1]];
3130  vars[1] = propdata->permvars[graphcomponents[k]];
3131 
3132  if ( storelexorder )
3133  (*lexorder)[*nvarsorder++] = graphcomponents[k];
3134 
3135  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "strong_sbcs_%s_%s", SCIPvarGetName(vars[0]), SCIPvarGetName(vars[1]));
3136 
3137  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, name, 2, vars, vals, 0.0,
3138  SCIPinfinity(scip), propdata->conssaddlp, propdata->conssaddlp, TRUE,
3139  TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
3140 
3141  SCIP_CALL( SCIPaddCons(scip, cons) );
3142 
3143 #ifdef SCIP_MORE_DEBUG
3144  SCIP_CALL( SCIPprintCons(scip, cons, NULL) );
3145  SCIPinfoMessage(scip, NULL, "\n");
3146 #endif
3147 
3148  /* check whether we need to resize */
3149  SCIP_CALL( ensureDynamicConsArrayAllocatedAndSufficientlyLarge(scip, &propdata->genlinconss,
3150  &propdata->genlinconsssize, propdata->ngenlinconss + 1) );
3151  propdata->genlinconss[propdata->ngenlinconss] = cons;
3152  ++propdata->ngenlinconss;
3153  }
3154 
3155  return SCIP_OKAY;
3156 }
3157 
3158 /** adds weak SBCs for a suitable color of the subgroup graph */
3159 static
3161  SCIP* scip, /**< SCIP instance */
3162  SCIP_PROPDATA* propdata, /**< pointer to data of symmetry propagator */
3163  int* compcolorbegins, /**< array indicating where a new graphcolor begins */
3164  int* graphcompbegins, /**< array indicating where a new graphcomponent begins */
3165  int* graphcomponents, /**< array of all variable indices sorted by color and comp */
3166  int ncompcolors, /**< number of colors in the graph */
3167  int* chosencomppercolor, /**< array indicating which comp was handled per color */
3168  int* firstvaridxpercolor,/**< array indicating the largest variable per color */
3169  int symgrpcompidx, /**< index of the component of the symmetry group */
3170  int* naddedconss, /**< buffer to store the number of added constraints */
3171  SCIP_Bool storelexorder, /**< whether the lexicographic order induced by the orbitope shall be stored */
3172  int** lexorder, /**< pointer to array storing lexicographic order defined by sub orbitopes */
3173  int* nvarsorder, /**< number of variables in lexicographic order */
3174  int* maxnvarsorder /**< maximum number of variables in lexicographic order */
3175  )
3176 { /*lint --e{571}*/
3177  SCIP_HASHMAP* varsinlexorder;
3178  SCIP_Shortbool* usedvars;
3179  SCIP_VAR* vars[2];
3180  SCIP_Real vals[2] = {1, -1};
3181  SCIP_Shortbool* varfound;
3182  int* orbit[2];
3183  int orbitsize[2] = {1, 1};
3184  int activeorb = 0;
3185  int chosencolor = -1;
3186  int j;
3187  int k;
3188 
3189  assert( scip != NULL );
3190  assert( propdata != NULL );
3191  assert( compcolorbegins != NULL );
3192  assert( graphcompbegins != NULL );
3193  assert( graphcomponents != NULL );
3194  assert( firstvaridxpercolor != NULL );
3195  assert( chosencomppercolor != NULL );
3196  assert( naddedconss != NULL );
3197  assert( symgrpcompidx >= 0 );
3198  assert( symgrpcompidx < propdata->ncomponents );
3199  assert( ! storelexorder || lexorder != NULL );
3200  assert( ! storelexorder || nvarsorder != NULL );
3201  assert( ! storelexorder || maxnvarsorder != NULL );
3202 
3203  *naddedconss = 0;
3204 
3205  SCIP_CALL( SCIPallocCleanBufferArray(scip, &usedvars, propdata->npermvars) );
3206  SCIP_CALL( SCIPallocClearBufferArray(scip, &varfound, propdata->npermvars) );
3207  SCIP_CALL( SCIPallocBufferArray(scip, &orbit[0], propdata->npermvars) );
3208  SCIP_CALL( SCIPallocBufferArray(scip, &orbit[1], propdata->npermvars) );
3209 
3210  /* Store the entries in lexorder in a hashmap, for fast lookups. */
3211  if ( lexorder == NULL || *lexorder == NULL )
3212  {
3213  /* Lexorder does not exist, so do not create hashmap. */
3214  varsinlexorder = NULL;
3215  }
3216  else
3217  {
3218  assert( *maxnvarsorder >= 0 );
3219  assert( *nvarsorder >= 0 );
3220 
3221  SCIP_CALL( SCIPhashmapCreate(&varsinlexorder, SCIPblkmem(scip), *maxnvarsorder) );
3222 
3223  for (k = 0; k < *nvarsorder; ++k)
3224  {
3225  /* add element from lexorder to hashmap.
3226  * Use insert, as duplicate entries in lexorder is not permitted. */
3227  assert((*lexorder)[k] >= 0);
3228  assert( ! SCIPhashmapExists(varsinlexorder, (void*) (size_t) (*lexorder)[k]) ); /* Use int as pointer */
3229  SCIP_CALL( SCIPhashmapInsertInt(varsinlexorder, (void*) (size_t) (*lexorder)[k], k) );
3230  }
3231  }
3232 
3233  /* We will store the newest and the largest orbit and activeorb will be used to mark at which entry of the array
3234  * orbit the newly computed one will be stored. */
3235  if ( ncompcolors > 0 )
3236  {
3237  SCIP_CALL( ensureSymmetryPermstransComputed(scip, propdata) );
3238  }
3239  for (j = 0; j < ncompcolors; ++j)
3240  {
3241  int graphcomp;
3242  int graphcompsize;
3243  int varidx;
3244 
3245  /* skip color for which we did not add anything */
3246  if ( chosencomppercolor[j] < 0 )
3247  continue;
3248 
3249  assert( firstvaridxpercolor[j] >= 0 );
3250 
3251  graphcomp = chosencomppercolor[j];
3252  graphcompsize = graphcompbegins[graphcomp+1] - graphcompbegins[graphcomp];
3253  varidx = firstvaridxpercolor[j];
3254  assert(varidx >= 0);
3255 
3256  /* if the first variable was already contained in another orbit or if there are no variables left anyway, skip the
3257  * component */
3258  if ( varfound[varidx] || graphcompsize == propdata->npermvars )
3259  continue;
3260 
3261  /* If varidx is in lexorder, then it must be the first entry of lexorder. */
3262  if ( varsinlexorder != NULL
3263  && SCIPhashmapExists(varsinlexorder, (void*) (size_t) varidx)
3264  && lexorder != NULL && *lexorder != NULL && *maxnvarsorder > 0 && *nvarsorder > 0
3265  && (*lexorder)[0] != varidx )
3266  continue;
3267 
3268  /* mark all variables that have been used in strong SBCs */
3269  for (k = graphcompbegins[graphcomp]; k < graphcompbegins[graphcomp+1]; ++k)
3270  {
3271  assert( 0 <= graphcomponents[k] && graphcomponents[k] < propdata->npermvars );
3272 
3273  usedvars[graphcomponents[k]] = TRUE;
3274  }
3275 
3276  SCIP_CALL( SCIPcomputeOrbitVar(scip, propdata->npermvars, propdata->perms,
3277  propdata->permstrans, propdata->components, propdata->componentbegins,
3278  usedvars, varfound, varidx, symgrpcompidx,
3279  orbit[activeorb], &orbitsize[activeorb]) );
3280 
3281  assert( orbit[activeorb][0] == varidx );
3282 
3283  if ( orbitsize[activeorb] > orbitsize[1 - activeorb] ) /*lint !e514*/
3284  {
3285  /* if the new orbit is larger then the old largest one, flip activeorb */
3286  activeorb = 1 - activeorb;
3287  chosencolor = j;
3288  }
3289 
3290  /* reset array */
3291  for (k = graphcompbegins[graphcomp]; k < graphcompbegins[graphcomp+1]; ++k)
3292  usedvars[graphcomponents[k]] = FALSE;
3293  }
3294 
3295  /* check if we have found at least one non-empty orbit */
3296  if ( chosencolor > -1 )
3297  {
3298  /* flip activeorb again to avoid confusion, it is then at the largest orbit */
3299  activeorb = 1 - activeorb;
3300 
3301  assert( orbit[activeorb][0] == firstvaridxpercolor[chosencolor] );
3302  vars[0] = propdata->permvars[orbit[activeorb][0]];
3303 
3304  assert( chosencolor > -1 );
3305  SCIPdebugMsg(scip, " adding %d weak sbcs for enclosing orbit of color %d.\n", orbitsize[activeorb]-1, chosencolor);
3306 
3307  *naddedconss = orbitsize[activeorb] - 1;
3308 
3309  /* add weak SBCs for rest of enclosing orbit */
3310  for (j = 1; j < orbitsize[activeorb]; ++j)
3311  {
3312  SCIP_CONS* cons;
3313  char name[SCIP_MAXSTRLEN];
3314 
3315  vars[1] = propdata->permvars[orbit[activeorb][j]];
3316 
3317  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "weak_sbcs_%d_%s_%s", symgrpcompidx, SCIPvarGetName(vars[0]), SCIPvarGetName(vars[1]));
3318 
3319  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, name, 2, vars, vals, 0.0,
3320  SCIPinfinity(scip), propdata->conssaddlp, propdata->conssaddlp, TRUE,
3321  TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
3322 
3323  SCIP_CALL( SCIPaddCons(scip, cons) );
3324 
3325 #ifdef SCIP_MORE_DEBUG
3326  SCIP_CALL( SCIPprintCons(scip, cons, NULL) );
3327  SCIPinfoMessage(scip, NULL, "\n");
3328 #endif
3329 
3330  /* check whether we need to resize */
3331  SCIP_CALL( ensureDynamicConsArrayAllocatedAndSufficientlyLarge(scip, &propdata->genlinconss,
3332  &propdata->genlinconsssize, propdata->ngenlinconss + 1) );
3333  propdata->genlinconss[propdata->ngenlinconss] = cons;
3334  ++propdata->ngenlinconss;
3335  }
3336 
3337  /* possibly store lexicographic order defined by weak SBCs */
3338  if ( storelexorder )
3339  {
3340  int varidx;
3341 
3342  varidx = orbit[activeorb][0];
3343  assert(varidx >= 0);
3344 
3345  if ( *maxnvarsorder == 0 )
3346  {
3347  *maxnvarsorder = 1;
3348  *nvarsorder = 0;
3349 
3350  SCIP_CALL( SCIPallocBlockMemoryArray(scip, lexorder, *maxnvarsorder) );
3351  (*lexorder)[(*nvarsorder)++] = varidx;
3352  }
3353  else
3354  {
3355  assert( *nvarsorder == *maxnvarsorder );
3356  assert( varsinlexorder != NULL );
3357  assert( lexorder != NULL );
3358  assert( *lexorder != NULL );
3359 
3360  /* the leader of the weak inequalities has to be the first element in the lexicographic order */
3361  if ( varidx == (*lexorder)[0] )
3362  {
3363  /* lexorder is already ok!! */
3364  assert( SCIPhashmapExists(varsinlexorder, (void*) (size_t) varidx) );
3365  }
3366  else
3367  {
3368  /* Then varidx must not be in the lexorder,
3369  * We must add it at the front of the array, and maintain the current order. */
3370  assert( ! SCIPhashmapExists(varsinlexorder, (void*) (size_t) varidx) );
3371 
3372  ++(*maxnvarsorder);
3373  ++(*nvarsorder);
3374 
3375  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, lexorder, *nvarsorder, *maxnvarsorder) );
3376 
3377  /* Shift array by one position to the right */
3378  for (k = *maxnvarsorder - 1; k >= 1; --k)
3379  (*lexorder)[k] = (*lexorder)[k - 1];
3380 
3381  (*lexorder)[0] = varidx;
3382  }
3383  }
3384  }
3385  }
3386  else
3387  SCIPdebugMsg(scip, " no further weak sbcs are valid\n");
3388 
3389  SCIPfreeBufferArray(scip, &orbit[1]);
3390  SCIPfreeBufferArray(scip, &orbit[0]);
3391  if ( varsinlexorder != NULL )
3392  SCIPhashmapFree(&varsinlexorder);
3393  SCIPfreeBufferArray(scip, &varfound);
3394  SCIPfreeCleanBufferArray(scip, &usedvars);
3395 
3396  return SCIP_OKAY;
3397 }
3398 
3399 
3400 /** temporarily adapt symmetry data to new variable order given by Schreier Sims */
3401 static
3403  SCIP* scip, /**< SCIP instance */
3404  int** origperms, /**< permutation matrix w.r.t. original variable ordering */
3405  int** modifiedperms, /**< memory for permutation matrix w.r.t. new variable ordering */
3406  int nperms, /**< number of permutations */
3407  SCIP_VAR** origpermvars, /**< array of permutation vars w.r.t. original variable ordering */
3408  SCIP_VAR** modifiedpermvars, /**< memory for array of permutation vars w.r.t. new variable ordering */
3409  int npermvars, /**< length or modifiedpermvars array */
3410  int* leaders, /**< leaders of Schreier Sims constraints */
3411  int nleaders /**< number of leaders */
3412  )
3413 {
3414  int* permvaridx;
3415  int* posinpermvar;
3416  int leader;
3417  int curposleader;
3418  int varidx;
3419  int lidx;
3420  int i;
3421  int l;
3422  int p;
3423 
3424  assert( scip != NULL );
3425  assert( origperms != NULL );
3426  assert( modifiedperms != NULL );
3427  assert( nperms > 0 );
3428  assert( origpermvars != NULL );
3429  assert( modifiedpermvars != NULL );
3430  assert( npermvars > 0 );
3431  assert( leaders != NULL );
3432  assert( nleaders > 0 );
3433 
3434  /* initialize map from position in lexicographic order to index of original permvar */
3435  SCIP_CALL( SCIPallocBufferArray(scip, &permvaridx, npermvars) );
3436  for (i = 0; i < npermvars; ++i)
3437  permvaridx[i] = i;
3438 
3439  /* initialize map from permvaridx to its current position in the reordered permvars array */
3440  SCIP_CALL( SCIPallocBufferArray(scip, &posinpermvar, npermvars) );
3441  for (i = 0; i < npermvars; ++i)
3442  posinpermvar[i] = i;
3443 
3444  /* Iterate over leaders and put the l-th leader to the l-th position of the lexicographic order.
3445  * We do this by swapping the l-th leader with the element at position l of the current permvars array. */
3446  for (l = 0; l < nleaders; ++l)
3447  {
3448  leader = leaders[l];
3449  curposleader = posinpermvar[leader];
3450  varidx = permvaridx[curposleader];
3451  lidx = permvaridx[l];
3452 
3453  /* swap the permvar at position l with the l-th leader */
3454  permvaridx[curposleader] = lidx;
3455  permvaridx[l] = varidx;
3456 
3457  /* update the position map */
3458  posinpermvar[lidx] = curposleader;
3459  posinpermvar[leader] = l;
3460  }
3461 
3462  /* update the permvars array to new variable order */
3463  for (i = 0; i < npermvars; ++i)
3464  modifiedpermvars[i] = origpermvars[permvaridx[i]];
3465 
3466  /* update the permutation to the new variable order */
3467  for (p = 0; p < nperms; ++p)
3468  {
3469  for (i = 0; i < npermvars; ++i)
3470  modifiedperms[p][i] = posinpermvar[origperms[p][permvaridx[i]]];
3471  }
3472 
3473  SCIPfreeBufferArray(scip, &permvaridx);
3474  SCIPfreeBufferArray(scip, &posinpermvar);
3475 
3476  return SCIP_OKAY;
3477 }
3478 
3479 
3480 /* returns the number of found orbitopes with at least three columns per graph component or 0
3481  * if the found orbitopes do not satisfy certain criteria for being used
3482  */
3483 static
3485  SCIP_VAR** permvars, /**< array of variables affected by symmetry */
3486  int* graphcomponents, /**< array of graph components */
3487  int* graphcompbegins, /**< array indicating starting position of graph components */
3488  int* compcolorbegins, /**< array indicating starting positions of potential orbitopes */
3489  int ncompcolors, /**< number of components encoded in compcolorbegins */
3490  int symcompsize /**< size of symmetry component for that we detect suborbitopes */
3491  )
3492 {
3493  SCIP_Bool oneorbitopecriterion = FALSE;
3494  SCIP_Bool multorbitopecriterion = FALSE;
3495  int norbitopes = 0;
3496  int j;
3497 
3498  assert( graphcompbegins != NULL );
3499  assert( compcolorbegins != NULL );
3500  assert( ncompcolors >= 0 );
3501  assert( symcompsize > 0 );
3502 
3503  for (j = 0; j < ncompcolors; ++j)
3504  {
3505  SCIP_VAR* firstvar;
3506  int largestcompsize = 0;
3507  int nbinrows= 0;
3508  int k;
3509 
3510  /* skip trivial components */
3511  if ( graphcompbegins[compcolorbegins[j+1]] - graphcompbegins[compcolorbegins[j]] < 2 )
3512  continue;
3513 
3514  /* check whether components of this color build an orbitope (with > 2 columns) */
3515  for (k = compcolorbegins[j]; k < compcolorbegins[j+1]; ++k)
3516  {
3517  int compsize;
3518 
3519  compsize = graphcompbegins[k+1] - graphcompbegins[k];
3520 
3521  /* the first component that we are looking at for this color */
3522  if ( largestcompsize < 1 )
3523  {
3524  if ( compsize < 3 )
3525  break;
3526 
3527  largestcompsize = compsize;
3528  }
3529  else if ( compsize != largestcompsize )
3530  break;
3531 
3532  firstvar = permvars[graphcomponents[graphcompbegins[k]]];
3533 
3534  /* count number of binary orbits (comps) */
3535  if ( SCIPvarIsBinary(firstvar) )
3536  ++nbinrows;
3537  }
3538 
3539  /* we have found an orbitope */
3540  if ( k == compcolorbegins[j+1] )
3541  {
3542  SCIP_Real threshold;
3543  int ncols;
3544 
3545  ++norbitopes;
3546  ncols = graphcompbegins[compcolorbegins[j] + 1] - graphcompbegins[compcolorbegins[j]];
3547 
3548  threshold = 0.7 * (SCIP_Real) symcompsize;
3549 
3550  /* check whether criteria for adding orbitopes are satisfied */
3551  if ( nbinrows <= 2 * ncols || (nbinrows <= 8 * ncols && nbinrows < 100) )
3552  multorbitopecriterion = TRUE;
3553  else if ( nbinrows <= 3 * ncols || (SCIP_Real) nbinrows * ncols >= threshold )
3554  oneorbitopecriterion = TRUE;
3555  }
3556  }
3557 
3558  if ( (norbitopes == 1 && oneorbitopecriterion) || (norbitopes >= 2 && multorbitopecriterion) )
3559  return norbitopes;
3560 
3561  return 0;
3562 }
3563 
3564 
3565 /** checks whether subgroups of the components are symmetric groups and adds SBCs for them */
3566 static
3568  SCIP* scip, /**< SCIP instance */
3569  SCIP_PROPDATA* propdata, /**< pointer to data of symmetry propagator */
3570  int cidx /**< index of component which shall be handled */
3571  )
3572 {
3573  int* genorder;
3574  int p;
3575 #ifdef SCIP_DEBUG
3576  int norbitopes = 0;
3577  int nstrongsbcs = 0;
3578  int nweaksbcs = 0;
3579 #endif
3580  int** modifiedperms;
3581  SCIP_VAR** modifiedpermvars;
3582  int* nvarsincomponent;
3583 
3584  int* graphcomponents;
3585  int* graphcompbegins;
3586  int* compcolorbegins;
3587  int* chosencomppercolor = NULL;
3588  int* firstvaridxpercolor = NULL;
3589  int* usedperms;
3590  int usedpermssize;
3591  int ngraphcomponents;
3592  int ncompcolors;
3593  int ntwocycleperms;
3594  int npermsincomp;
3595  int nusedperms;
3596  int ntrivialcolors = 0;
3597  int j;
3598  int* lexorder = NULL;
3599  int nvarslexorder = 0;
3600  int maxnvarslexorder = 0;
3601  SCIP_Shortbool* permused;
3602  SCIP_Bool allpermsused = FALSE;
3603  SCIP_Bool handlednonbinarysymmetry = FALSE;
3604  int norbitopesincomp;
3605 
3606  assert( scip != NULL );
3607  assert( propdata != NULL );
3608  assert( propdata->computedsymmetry );
3609  assert( propdata->nperms >= 0 );
3610  assert( 0 <= cidx && cidx < propdata->ncomponents );
3611  assert( propdata->components != NULL );
3612  assert( propdata->componentbegins != NULL );
3613 
3614  /* exit if no symmetry is present or component is blocked */
3615  if ( propdata->nperms == 0 || propdata->componentblocked[cidx] )
3616  return SCIP_OKAY;
3617 
3618  /* exit if instance is too large */
3619  if ( SCIPgetNConss(scip) > propdata->maxnconsssubgroup )
3620  return SCIP_OKAY;
3621 
3622  assert( propdata->nperms > 0 );
3623  assert( propdata->perms != NULL );
3624  assert( propdata->npermvars > 0 );
3625  assert( propdata->permvars != NULL );
3626 
3627  /* create array for permutation order */
3628  SCIP_CALL( SCIPallocBufferArray(scip, &genorder, propdata->nperms) );
3629 
3630  /* create arrays for modified permutations in case we adapt the lexicographic order because of suborbitopes */
3631  SCIP_CALL( SCIPallocBufferArray(scip, &modifiedperms, propdata->nperms) );
3632  for (p = 0; p < propdata->nperms; ++p)
3633  {
3634  SCIP_CALL( SCIPallocBufferArray(scip, &modifiedperms[p], propdata->npermvars) );
3635  }
3636  SCIP_CALL( SCIPallocBufferArray(scip, &modifiedpermvars, propdata->npermvars) );
3637 
3638  SCIP_CALL( SCIPallocClearBufferArray(scip, &nvarsincomponent, propdata->npermvars) );
3639  for (p = 0; p < propdata->npermvars; ++p)
3640  {
3641  if ( propdata->vartocomponent[p] >= 0 )
3642  ++nvarsincomponent[propdata->vartocomponent[p]];
3643  }
3644 
3645  SCIPdebugMsg(scip, "starting subgroup detection routine for component %d\n", cidx);
3646 
3647  npermsincomp = propdata->componentbegins[cidx + 1] - propdata->componentbegins[cidx];
3648 
3649  /* set the first npermsincomp entries of genorder; the others are not used for this component */
3650  for (j = 0; j < npermsincomp; ++j)
3651  genorder[j] = j;
3652 
3653  SCIP_CALL( chooseOrderOfGenerators(scip, propdata, cidx, &genorder, &ntwocycleperms) );
3654 
3655  assert( ntwocycleperms >= 0 );
3656  assert( ntwocycleperms <= npermsincomp );
3657 
3658  SCIPdebugMsg(scip, "component %d has %d permutations consisting of 2-cycles\n", cidx, ntwocycleperms);
3659 
3660 #ifdef SCIP_MORE_DEBUG
3661  SCIP_Bool* used;
3662  int perm;
3663  int p;
3664  int k;
3665 
3666  SCIP_CALL( SCIPallocBufferArray(scip, &used, propdata->npermvars) );
3667  for (p = propdata->componentbegins[cidx]; p < propdata->componentbegins[cidx+1]; ++p)
3668  {
3669  perm = propdata->components[p];
3670 
3671  for (k = 0; k < propdata->npermvars; ++k)
3672  used[k] = FALSE;
3673 
3674  SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, "permutation %d\n", perm);
3675 
3676  for (k = 0; k < propdata->npermvars; ++k)
3677  {
3678  if ( used[k] )
3679  continue;
3680 
3681  j = propdata->perms[perm][k];
3682 
3683  if ( k == j )
3684  continue;
3685 
3686  SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, "(%s,", SCIPvarGetName(propdata->permvars[k]));
3687  used[k] = TRUE;
3688  while (j != k)
3689  {
3690  SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, "%s,", SCIPvarGetName(propdata->permvars[j]));
3691  used[j] = TRUE;
3692 
3693  j = propdata->perms[perm][j];
3694  }
3696  }
3698  }
3699 
3700  SCIPfreeBufferArray(scip, &used);
3701 #endif
3702 
3703  if ( ntwocycleperms < 2 )
3704  {
3705  SCIPdebugMsg(scip, " --> skip\n");
3706  goto FREEBASICMEM;
3707  }
3708 
3709  usedpermssize = ntwocycleperms / 2;
3710  SCIP_CALL( SCIPallocBufferArray(scip, &usedperms, usedpermssize) );
3711  SCIP_CALL( SCIPallocClearBufferArray(scip, &permused, npermsincomp) );
3712 
3713  SCIP_CALL( buildSubgroupGraph(scip, propdata, genorder, ntwocycleperms, cidx,
3714  &graphcomponents, &graphcompbegins, &compcolorbegins, &ngraphcomponents,
3715  &ncompcolors, &usedperms, &nusedperms, usedpermssize, permused) );
3716 
3717  SCIPdebugMsg(scip, " created subgroup detection graph using %d of the permutations\n", nusedperms);
3718 
3719  if ( nusedperms == npermsincomp )
3720  allpermsused = TRUE;
3721 
3722  assert( graphcomponents != NULL );
3723  assert( graphcompbegins != NULL );
3724  assert( compcolorbegins != NULL );
3725  assert( ngraphcomponents > 0 );
3726  assert( ncompcolors > 0 );
3727  assert( nusedperms <= ntwocycleperms );
3728  assert( ncompcolors < propdata->npermvars );
3729 
3730  if ( nusedperms == 0 )
3731  {
3732  SCIPdebugMsg(scip, " -> skipping component, since less no permutation was used\n");
3733 
3734  SCIPfreeBufferArray(scip, &permused);
3735  SCIPfreeBufferArray(scip, &usedperms);
3736 
3737  goto FREEBASICMEM;
3738  }
3739 
3740  SCIPdebugMsg(scip, " number of different colors in the graph: %d\n", ncompcolors);
3741 
3742  if ( propdata->addstrongsbcs || propdata->addweaksbcs )
3743  {
3744  SCIP_CALL( SCIPallocBufferArray(scip, &chosencomppercolor, ncompcolors) );
3745  SCIP_CALL( SCIPallocBufferArray(scip, &firstvaridxpercolor, ncompcolors) );
3746 
3747  /* Initialize the arrays with -1 to encode that we have not added orbitopes/strong SBCs
3748  * yet. In case we do not modify this entry, no weak inequalities are added based on
3749  * this component.
3750  */
3751  for (j = 0; j < ncompcolors; ++j)
3752  {
3753  chosencomppercolor[j] = -1;
3754  firstvaridxpercolor[j] = -1;
3755  }
3756  }
3757 
3758  norbitopesincomp = getNOrbitopesInComp(propdata->permvars, graphcomponents, graphcompbegins, compcolorbegins,
3759  ncompcolors, nvarsincomponent[cidx]);
3760 
3761  /* if there is just one orbitope satisfying the requirements, handle the full component by symresacks */
3762  if ( norbitopesincomp == 1 )
3763  {
3764  int k;
3765 
3766  for (k = 0; k < npermsincomp; ++k)
3767  {
3768  SCIP_CONS* cons;
3769  char name[SCIP_MAXSTRLEN];
3770 
3771  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "symresack_comp%d_perm%d", cidx, k);
3772 
3773  SCIP_CALL( SCIPcreateSymbreakCons(scip, &cons, name,
3774  propdata->perms[propdata->components[propdata->componentbegins[cidx] + k]],
3775  propdata->permvars, propdata->npermvars, FALSE,
3776  propdata->conssaddlp, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
3777  SCIP_CALL( SCIPaddCons(scip, cons));
3778 
3779  /* do not release constraint here - will be done later */
3780  SCIP_CALL( ensureDynamicConsArrayAllocatedAndSufficientlyLarge(scip, &propdata->genorbconss,
3781  &propdata->genorbconsssize, propdata->ngenorbconss + 1) );
3782  propdata->genorbconss[propdata->ngenorbconss++] = cons;
3783  ++propdata->nsymresacks;
3784 
3785  if ( ! propdata->componentblocked[cidx] )
3786  {
3787  propdata->componentblocked[cidx] |= SYM_HANDLETYPE_SYMBREAK;
3788  ++propdata->ncompblocked;
3789  }
3790 
3791  SCIPdebugMsg(scip, " add symresack for permutation %d of component %d\n", k, cidx);
3792  }
3793 
3794  goto FREEMEMORY;
3795  }
3796 
3797  for (j = 0; j < ncompcolors; ++j)
3798  {
3799  int nbinarycomps = 0;
3800  int largestcolorcomp = -1;
3801  int largestcompsize = 0;
3802  int k;
3803  SCIP_Bool isorbitope = TRUE;
3804  SCIP_Bool orbitopeadded = FALSE;
3805  SCIP_Bool useorbitope;
3806 #ifdef SCIP_DEBUG
3807  SCIP_Bool binaffected = FALSE;
3808  SCIP_Bool intaffected = FALSE;
3809  SCIP_Bool contaffected = FALSE;
3810 #endif
3811 
3812  /* skip trivial components */
3813  if ( graphcompbegins[compcolorbegins[j+1]] - graphcompbegins[compcolorbegins[j]] < 2 )
3814  {
3815  if( chosencomppercolor != NULL )
3816  chosencomppercolor[j] = -1;
3817 
3818  ++ntrivialcolors;
3819  continue;
3820  }
3821 
3822  SCIPdebugMsg(scip, " color %d has %d components with overall %d variables\n", j, compcolorbegins[j+1] - compcolorbegins[j],
3823  graphcompbegins[compcolorbegins[j+1]] - graphcompbegins[compcolorbegins[j]]);
3824 
3825  /* check whether components of this color might build an orbitope (with > 2 columns) */
3826  for (k = compcolorbegins[j]; k < compcolorbegins[j+1]; ++k)
3827  {
3828  SCIP_VAR* firstvar;
3829  int compsize;
3830 
3831  compsize = graphcompbegins[k+1] - graphcompbegins[k];
3832 
3833  /* the first component that we are looking at for this color */
3834  if ( largestcompsize < 1 )
3835  {
3836  if ( compsize < 3 )
3837  {
3838  isorbitope = FALSE;
3839  break;
3840  }
3841 
3842  largestcompsize = compsize;
3843  largestcolorcomp = k;
3844  }
3845  else if ( compsize != largestcompsize )
3846  {
3847  /* variable orbits (compsize) have not the same size, cannot define orbitope */
3848  isorbitope = FALSE;
3849  break;
3850  }
3851 
3852  firstvar = propdata->permvars[graphcomponents[graphcompbegins[k]]];
3853 
3854  /* count number of binary orbits (comps) */
3855  if ( SCIPvarIsBinary(firstvar) )
3856  ++nbinarycomps;
3857  }
3858 
3859 #ifdef SCIP_DEBUG
3860  for (k = compcolorbegins[j]; k < compcolorbegins[j+1]; ++k)
3861  {
3862  SCIP_VAR* firstvar;
3863 
3864  firstvar = propdata->permvars[graphcomponents[graphcompbegins[k]]];
3865 
3866  if ( SCIPvarIsBinary(firstvar) )
3867  binaffected = TRUE;
3868  else if (SCIPvarIsIntegral(firstvar) )
3869  intaffected = TRUE;
3870  else
3871  contaffected = TRUE;
3872  }
3873 
3874  SCIPdebugMsg(scip, " affected types (bin,int,cont): (%d,%d,%d)\n", binaffected, intaffected, contaffected);
3875 #endif
3876 
3877  /* only use the orbitope if there are binary rows */
3878  useorbitope = FALSE;
3879  if ( norbitopesincomp > 0 && nbinarycomps > 0 )
3880  useorbitope = TRUE;
3881 
3882  if ( isorbitope && useorbitope )
3883  {
3884  int firstvaridx;
3885  int chosencomp;
3886 
3887  SCIPdebugMsg(scip, " detected an orbitope with %d rows and %d columns\n", nbinarycomps, largestcompsize);
3888 
3889  assert( nbinarycomps > 0 );
3890  assert( largestcompsize > 2 );
3891 
3892  /* add the orbitope constraint for this color
3893  *
3894  * It might happen that we cannot generate the orbitope matrix if the orbitope is not generated by permutations
3895  * all having the same number of 2-cycles, e.g., the orbitope generated by (1,2)(4,5), (2,3), (5,6).
3896  */
3897  SCIP_CALL( addOrbitopeSubgroup(scip, propdata, usedperms, nusedperms, compcolorbegins,
3898  graphcompbegins, graphcomponents, j, nbinarycomps, largestcompsize, &firstvaridx, &chosencomp,
3899  &lexorder, &nvarslexorder, &maxnvarslexorder, allpermsused, &orbitopeadded) );
3900 
3901  if ( orbitopeadded )
3902  {
3903  if ( propdata->addstrongsbcs || propdata->addweaksbcs )
3904  {
3905  assert( chosencomppercolor != NULL );
3906  assert( firstvaridxpercolor != NULL );
3907 
3908  /* adapt the first variable per color to be compatible with the created orbiope (upper left variable) */
3909  assert( compcolorbegins[j] <= chosencomp && chosencomp < compcolorbegins[j+1] );
3910  assert( 0 <= firstvaridx && firstvaridx < propdata->npermvars );
3911 
3912  chosencomppercolor[j] = chosencomp;
3913  firstvaridxpercolor[j] = firstvaridx;
3914  }
3915 
3916  if ( ! propdata->componentblocked[cidx] )
3917  {
3918  propdata->componentblocked[cidx] |= SYM_HANDLETYPE_SYMBREAK;
3919  ++propdata->ncompblocked;
3920  }
3921 
3922 #ifdef SCIP_DEBUG
3923  ++norbitopes;
3924 #endif
3925  }
3926  }
3927 
3928  /* if no (useable) orbitope was found, possibly add strong SBCs */
3929  if ( propdata->addstrongsbcs && ! orbitopeadded )
3930  {
3931  assert( largestcolorcomp >= 0 );
3932  assert( largestcolorcomp < ngraphcomponents );
3933  assert( largestcompsize > 0 );
3934 
3935  if( propdata->addweaksbcs )
3936  {
3937  assert( chosencomppercolor != NULL );
3938  assert( firstvaridxpercolor != NULL );
3939 
3940  chosencomppercolor[j] = largestcolorcomp;
3941  firstvaridxpercolor[j] = graphcomponents[graphcompbegins[largestcolorcomp]];
3942  }
3943 
3944  SCIPdebugMsg(scip, " choosing component %d with %d variables and adding strong SBCs\n",
3945  largestcolorcomp, graphcompbegins[largestcolorcomp+1] - graphcompbegins[largestcolorcomp]);
3946 
3947  /* add the strong SBCs for the corresponding component */
3948  SCIP_CALL( addStrongSBCsSubgroup(scip, propdata, graphcompbegins, graphcomponents, largestcolorcomp,
3949  propdata->addsymresacks, &lexorder, &nvarslexorder, &maxnvarslexorder) );
3950 
3951  /* store whether symmetries on non-binary symmetries have been handled */
3952  if ( ! SCIPvarIsBinary(propdata->permvars[graphcomponents[graphcompbegins[largestcolorcomp]]]) )
3953  handlednonbinarysymmetry = TRUE;
3954 
3955  if ( ! propdata->componentblocked[cidx] )
3956  {
3957  propdata->componentblocked[cidx] |= SYM_HANDLETYPE_SYMBREAK;
3958  ++propdata->ncompblocked;
3959  }
3960 
3961 #ifdef SCIP_DEBUG
3962  nstrongsbcs += graphcompbegins[largestcolorcomp+1] - graphcompbegins[largestcolorcomp] - 1;
3963 #endif
3964  }
3965  else if ( ! orbitopeadded )
3966  {
3967  SCIPdebugMsg(scip, " no useable orbitope found and no SBCs added\n");
3968 
3969  /* mark the color as not handled */
3970  if ( propdata->addweaksbcs )
3971  {
3972  assert( chosencomppercolor != NULL );
3973  chosencomppercolor[j] = -1; /*lint !e613*/
3974  }
3975  }
3976  }
3977 
3978  SCIPdebugMsg(scip, " skipped %d trivial colors\n", ntrivialcolors);
3979 
3980  /* possibly add weak SBCs for enclosing orbit of first component */
3981  if ( propdata->addweaksbcs && propdata->componentblocked[cidx] && nusedperms < npermsincomp )
3982  {
3983  int naddedconss;
3984 
3985  assert( firstvaridxpercolor != NULL );
3986  assert( chosencomppercolor != NULL );
3987 
3988  SCIP_CALL( addWeakSBCsSubgroup(scip, propdata, compcolorbegins, graphcompbegins,
3989  graphcomponents, ncompcolors, chosencomppercolor, firstvaridxpercolor,
3990  cidx, &naddedconss, propdata->addsymresacks, &lexorder, &nvarslexorder, &maxnvarslexorder) );
3991 
3992  assert( naddedconss < propdata->npermvars );
3993 
3994 #ifdef SCIP_DEBUG
3995  nweaksbcs += naddedconss;
3996 #endif
3997  }
3998  else
3999  SCIPdebugMsg(scip, " don't add weak sbcs because all generators were used or the settings forbid it\n");
4000 
4001  /* if suborbitopes or strong group actions have been found, potentially add symresacks adapted to
4002  * variable order given by lexorder if no symmetries on non-binary variables have been handled
4003  */
4004  if ( nvarslexorder > 0 && propdata->addsymresacks && ! handlednonbinarysymmetry )
4005  {
4006  int k;
4007 
4008  SCIP_CALL( adaptSymmetryDataSST(scip, propdata->perms, modifiedperms, propdata->nperms,
4009  propdata->permvars, modifiedpermvars, propdata->npermvars, lexorder, nvarslexorder) );
4010 
4011  for (k = 0; k < npermsincomp; ++k)
4012  {
4013  SCIP_CONS* cons;
4014  char name[SCIP_MAXSTRLEN];
4015  int* symresackperm;
4016  SCIP_Bool actsonbinary = FALSE;
4017 
4018  /* skip permutations that have been used to build an orbitope */
4019  if ( permused[k] )
4020  continue;
4021 
4022  /* skip permutations that do not act on binary variables */
4023  symresackperm = modifiedperms[propdata->components[propdata->componentbegins[cidx] + k]];
4024  for (j = 0; j < propdata->nperms && !actsonbinary; ++j)
4025  {
4026  if ( symresackperm[j] != j && SCIPvarIsBinary(modifiedpermvars[j]) )
4027  actsonbinary = TRUE;
4028  }
4029 
4030  if ( ! actsonbinary )
4031  continue;
4032 
4033  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "symresack_comp%d_perm%d", cidx, k);
4034 
4035  SCIP_CALL( SCIPcreateSymbreakCons(scip, &cons, name,
4036  symresackperm, modifiedpermvars, propdata->npermvars, FALSE,
4037  propdata->conssaddlp, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
4038  SCIP_CALL( SCIPaddCons(scip, cons));
4039 
4040  /* do not release constraint here - will be done later */
4041  SCIP_CALL( ensureDynamicConsArrayAllocatedAndSufficientlyLarge(scip, &propdata->genorbconss,
4042  &propdata->genorbconsssize, propdata->ngenorbconss + 1) );
4043  propdata->genorbconss[propdata->ngenorbconss++] = cons;
4044  ++propdata->nsymresacks;
4045 
4046  if ( ! propdata->componentblocked[cidx] )
4047  {
4048  propdata->componentblocked[cidx] |= SYM_HANDLETYPE_SYMBREAK;
4049  ++propdata->ncompblocked;
4050  }
4051 
4052  SCIPdebugMsg(scip, " add symresack for permutation %d of component %d adapted to suborbitope lexorder\n", k, cidx);
4053  }
4054  }
4055 
4056  FREEMEMORY:
4057  SCIPfreeBlockMemoryArrayNull(scip, &lexorder, maxnvarslexorder);
4058 
4059  SCIPfreeBufferArrayNull(scip, &firstvaridxpercolor);
4060  SCIPfreeBufferArrayNull(scip, &chosencomppercolor);
4061  SCIPfreeBlockMemoryArrayNull(scip, &compcolorbegins, ncompcolors + 1);
4062  SCIPfreeBlockMemoryArrayNull(scip, &graphcompbegins, ngraphcomponents + 1);
4063  SCIPfreeBlockMemoryArrayNull(scip, &graphcomponents, propdata->npermvars);
4064  SCIPfreeBufferArrayNull(scip, &permused);
4065  SCIPfreeBufferArrayNull(scip, &usedperms);
4066 
4067 #ifdef SCIP_DEBUG
4068  SCIPdebugMsg(scip, "total number of added (sub-)orbitopes: %d\n", norbitopes);
4069  SCIPdebugMsg(scip, "total number of added strong sbcs: %d\n", nstrongsbcs);
4070  SCIPdebugMsg(scip, "total number of added weak sbcs: %d\n", nweaksbcs);
4071 #endif
4072 
4073  FREEBASICMEM:
4074  SCIPfreeBufferArray(scip, &nvarsincomponent);
4075 
4076  SCIPfreeBufferArray(scip, &modifiedpermvars);
4077  for (p = propdata->nperms - 1; p >= 0; --p)
4078  {
4079  SCIPfreeBufferArray(scip, &modifiedperms[p]);
4080  }
4081  SCIPfreeBufferArray(scip, &modifiedperms);
4082  SCIPfreeBufferArray(scip, &genorder);
4083 
4084  return SCIP_OKAY;
4085 }
4086 
4087 
4088 /*
4089  * Functions for symmetry constraints
4090  */
4091 
4092 
4093 /** update symmetry information of conflict graph */
4094 static
4096  SCIP* scip, /**< SCIP instance */
4097  SCIP_CONFLICTDATA* varconflicts, /**< conflict structure */
4098  SCIP_VAR** conflictvars, /**< variables encoded in conflict structure */
4099  int nconflictvars, /**< number of nodes/vars in conflict structure */
4100  int* orbits, /**< array of non-trivial orbits */
4101  int* orbitbegins, /**< array containing begin positions of new orbits in orbits array */
4102  int norbits /**< number of non-trivial orbits */
4103  )
4104 {
4105  int i;
4106  int j;
4107  int ii;
4108  int jj;
4109  int r; /* r from orbit, the orbit index. */
4110 
4111  assert( scip != NULL );
4112  assert( varconflicts != NULL );
4113  assert( conflictvars != NULL );
4114  assert( nconflictvars > 0 );
4115  assert( orbits != NULL );
4116  assert( orbitbegins != NULL );
4117  assert( norbits >= 0 );
4118 
4119  /* initialize/reset variable information of nodes in conflict graph */
4120  for (i = 0; i < nconflictvars; ++i)
4121  {
4122  /* (re-)set node data */
4123  varconflicts[i].orbitidx = -1;
4124  varconflicts[i].nconflictinorbit = 0;
4125  varconflicts[i].orbitsize = -1;
4126  varconflicts[i].posinorbit = -1;
4127  }
4128 
4129  /* add orbit information to nodes of conflict graph */
4130  for (r = 0; r < norbits; ++r)
4131  {
4132  int posinorbit = 0;
4133  int orbitsize;
4134 
4135  orbitsize = orbitbegins[r + 1] - orbitbegins[r];
4136  assert( orbitsize >= 0 );
4137 
4138  for (i = orbitbegins[r]; i < orbitbegins[r + 1]; ++i)
4139  {
4140  int pos;
4141 
4142  /* get variable and position in conflict graph */
4143  pos = orbits[i];
4144  assert( pos < nconflictvars );
4145  assert( varconflicts[pos].var == conflictvars[pos] );
4146 
4147  varconflicts[pos].orbitidx = r;
4148  varconflicts[pos].nconflictinorbit = 0;
4149  varconflicts[pos].orbitsize = orbitsize;
4150  varconflicts[pos].posinorbit = posinorbit++;
4151  }
4152 
4153  /* determine nconflictsinorbit
4154  *
4155  * For each pair of active variables in this orbit, check if it is part of a conflict clique.
4156  * Use that we store the cliques of this type in varconflicts[pos].cliques.
4157  * These lists are sorted (by the address of the constraint), so we only need to check for each i, j in the orbit
4158  * whether they are contained in the same clique.
4159  */
4160  for (i = orbitbegins[r]; i < orbitbegins[r + 1]; ++i)
4161  {
4162  ii = orbits[i];
4163  assert( varconflicts[ii].orbitidx == r );
4164 
4165  /* skip inactive variables */
4166  if ( ! varconflicts[ii].active )
4167  continue;
4168 
4169  for (j = i + 1; j < orbitbegins[r + 1]; ++j)
4170  {
4171  jj = orbits[j];
4172  assert( varconflicts[jj].orbitidx == r );
4173 
4174  /* skip inactive variables */
4175  if ( ! varconflicts[jj].active )
4176  continue;
4177 
4178  /* Check if i and j are overlapping in some clique, where only one of the two could have value 1.
4179  * Use that cliques are sorted by the constraint address.
4180  *
4181  * @todo A better sorted order would be: First constraints with large variables (higher hitting probability)
4182  * and then by a unique constraint identifier (address, or conspos).
4183  */
4184  if ( checkSortedArraysHaveOverlappingEntry((void**)varconflicts[ii].cliques,
4185  varconflicts[ii].ncliques, (void**)varconflicts[jj].cliques, varconflicts[jj].ncliques,
4186  sortByPointerValue) )
4187  {
4188  /* there is overlap! */
4189  ++varconflicts[ii].nconflictinorbit;
4190  ++varconflicts[jj].nconflictinorbit;
4191  }
4192  }
4193  }
4194  }
4195 
4196  return SCIP_OKAY;
4197 }
4198 
4199 
4200 /** create conflict graph either for symmetric or for all variables
4201  *
4202  * This routine just creates the graph, but does not add (symmetry) information to its nodes.
4203  * This has to be done separately by the routine updateSymInfoConflictGraphSST().
4204  *
4205  * The function returns with varconflicts as NULL when we do not create it.
4206  */
4207 static
4209  SCIP* scip, /**< SCIP instance */
4210  SCIP_CONFLICTDATA** varconflicts, /**< pointer to store the variable conflict data */
4211  SCIP_VAR** conflictvars, /**< array of variables to encode in conflict graph */
4212  int nconflictvars, /**< number of vars to encode in conflict graph */
4213  SCIP_HASHMAP* conflictvarmap /**< map of variables to indices in conflictvars array */
4214  )
4215 {
4216  SCIP_CLIQUE** cliques;
4217  SCIP_VAR** cliquevars;
4218  SCIP_CLIQUE* clique;
4219  int* tmpncliques;
4220  int ncliques;
4221  int ncliquevars;
4222  int node;
4223  int c;
4224  int i;
4225 
4226 #ifdef SCIP_DEBUG
4227  int varncliques = 0;
4228 #endif
4229 
4230  assert( scip != NULL );
4231  assert( varconflicts != NULL );
4232  assert( conflictvars != NULL );
4233  assert( nconflictvars > 0 );
4234 
4235  /* we set the pointer of varconflicts to NULL to illustrate that we didn't generate it */
4236  *varconflicts = NULL;
4237 
4238  /* get cliques for creating conflict structure */
4239 
4240  cliques = SCIPgetCliques(scip);
4241  ncliques = SCIPgetNCliques(scip);
4242  if ( ncliques == 0 )
4243  {
4244  SCIPdebugMsg(scip, "No cliques present --> construction of conflict structure aborted.\n");
4245  return SCIP_OKAY;
4246  }
4247 
4248  /* construct variable conflicts */
4249  SCIPdebugMsg(scip, "Construction of conflict structure:\n");
4250  SCIP_CALL( SCIPallocBlockMemoryArray(scip, varconflicts, nconflictvars) );
4251  for (i = 0; i < nconflictvars; ++i)
4252  {
4253  (*varconflicts)[i].ncliques = 0;
4254  (*varconflicts)[i].active = TRUE;
4255  (*varconflicts)[i].var = conflictvars[i];
4256  /* set remaining variable conflictdata at neutral entries */
4257  (*varconflicts)[i].cliques = NULL;
4258  (*varconflicts)[i].orbitidx = -1;
4259  (*varconflicts)[i].nconflictinorbit = 0;
4260  (*varconflicts)[i].orbitsize = -1;
4261  (*varconflicts)[i].posinorbit = -1;
4262  }
4263 
4264  /* Store, for each variable, the conflict cliques it is contained in.
4265  * In three steps:
4266  * (1.) Count the number of cliques it's contained in, per var, then
4267  * (2.) Create the array of this size, and
4268  * (3.) Fill the array with the cliques.
4269  * Starting with (1.):
4270  */
4271  for (c = 0; c < ncliques; ++c)
4272  {
4273  clique = cliques[c];
4274  assert( clique != NULL );
4275 
4276  cliquevars = SCIPcliqueGetVars(clique);
4277  ncliquevars = SCIPcliqueGetNVars(clique);
4278  assert( cliquevars != NULL );
4279  assert( ncliquevars > 0 );
4280 
4281  SCIPdebugMsg(scip, "\tIdentify edges for clique ID: %d; Index: %d).\n", SCIPcliqueGetId(clique),
4282  SCIPcliqueGetIndex(clique));
4283 
4284  /* for all variables, list which cliques it is part of */
4285  for (i = 0; i < ncliquevars; ++i)
4286  {
4287  node = SCIPhashmapGetImageInt(conflictvarmap, cliquevars[i]);
4288 
4289  /* skip variables not in the conflictvars array (so not in hashmap, too) */
4290  if ( node == INT_MAX )
4291  continue;
4292  assert( node >= 0 );
4293  assert( node < nconflictvars );
4294 
4295  assert( (*varconflicts)[node].var == cliquevars[i] );
4296  (*varconflicts)[node].active = TRUE;
4297  (*varconflicts)[node].ncliques++;
4298  }
4299  }
4300 
4301  /* (2.) allocate the arrays */
4302  for (i = 0; i < nconflictvars; ++i)
4303  {
4304  assert( (*varconflicts)[i].ncliques >= 0 );
4305  assert( (*varconflicts)[i].cliques == NULL );
4306  if ( (*varconflicts)[i].ncliques > 0 )
4307  {
4308  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*varconflicts)[i].cliques, (*varconflicts)[i].ncliques) );
4309  }
4310  }
4311 
4312  /* (3.) fill the clique constraints */
4313  SCIP_CALL( SCIPallocClearBufferArray(scip, &tmpncliques, nconflictvars) );
4314  for (c = 0; c < ncliques; ++c)
4315  {
4316  clique = cliques[c];
4317  assert( clique != NULL );
4318 
4319  cliquevars = SCIPcliqueGetVars(clique);
4320  ncliquevars = SCIPcliqueGetNVars(clique);
4321  assert( cliquevars != NULL );
4322  assert( ncliquevars > 0 );
4323 
4324  SCIPdebugMsg(scip, "\tAdd edges for clique ID: %d; Index: %d).\n", SCIPcliqueGetId(clique),
4325  SCIPcliqueGetIndex(clique));
4326 
4327  /* for all variables, list which cliques it is part of */
4328  for (i = 0; i < ncliquevars; ++i)
4329  {
4330  node = SCIPhashmapGetImageInt(conflictvarmap, cliquevars[i]);
4331 
4332  /* skip variables not in the conflictvars array (so not in hashmap, too) */
4333  if ( node == INT_MAX )
4334  continue;
4335 
4336  assert( node >= 0 );
4337  assert( node < nconflictvars );
4338  assert( (*varconflicts)[node].var == cliquevars[i] );
4339 
4340  /* add clique to the cliques */
4341  assert( tmpncliques[node] < (*varconflicts)[node].ncliques );
4342  assert( (*varconflicts)[node].cliques != NULL );
4343  (*varconflicts)[node].cliques[tmpncliques[node]++] = clique;
4344 
4345 #ifdef SCIP_DEBUG
4346  varncliques++;
4347 #endif
4348  }
4349  }
4350 
4351  /* sort the variable cliques by the address, so checkSortedArraysHaveOverlappingEntry can detect intersections */
4352  for (i = 0; i < nconflictvars; ++i)
4353  {
4354  SCIPsortPtr((void**)(*varconflicts)[i].cliques, sortByPointerValue, (*varconflicts)[i].ncliques);
4355  }
4356 
4357 #ifndef NDEBUG
4358  for (i = 0; i < nconflictvars; ++i)
4359  {
4360  assert( tmpncliques[i] == (*varconflicts)[i].ncliques );
4361  }
4362 #endif
4363 
4364  SCIPfreeBufferArray(scip, &tmpncliques);
4365 
4366 #ifdef SCIP_DEBUG
4367  SCIPdebugMsg(scip, "Construction of conflict graph terminated; %d variable-clique combinations detected.\n",
4368  varncliques);
4369 #endif
4370 
4371  return SCIP_OKAY;
4372 }
4373 
4374 /** frees conflict graph */
4375 static
4377  SCIP* scip, /**< SCIP instance */
4378  SCIP_CONFLICTDATA** varconflicts, /**< conflict graph */
4379  int nvars /**< number of nodes in conflict graph */
4380  )
4381 {
4382  int i;
4383  int n;
4384 
4385  assert( scip != NULL );
4386  assert( varconflicts != NULL );
4387  assert( *varconflicts != NULL );
4388  assert( nvars >= 0 );
4389 
4390  for (i = nvars - 1; i >= 0; --i)
4391  {
4392  n = (*varconflicts)[i].ncliques;
4393  SCIPfreeBlockMemoryArray(scip, &(*varconflicts)[i].cliques, n);
4394  }
4395  SCIPfreeBlockMemoryArray(scip, varconflicts, nvars);
4396 
4397  return SCIP_OKAY;
4398 }
4399 
4400 
4401 /** adds symresack constraints */
4402 static
4404  SCIP* scip, /**< SCIP instance */
4405  SCIP_PROPDATA* propdata, /**< data of symmetry propagator */
4406  int cidx /**< index of component to be handled */
4407  )
4408 { /*lint --e{641}*/
4409  int* components;
4410  int* componentbegins;
4411  SCIP_VAR** permvars;
4412  SCIP_Bool conssaddlp;
4413  int** modifiedperms = NULL;
4414  SCIP_VAR** modifiedpermvars = NULL;
4415  int** perms;
4416  int nsymresackcons = 0;
4417  int npermvars;
4418  int nperms;
4419  int i;
4420  int p;
4421 
4422  assert( scip != NULL );
4423  assert( propdata != NULL );
4424  assert( propdata->npermvars >= 0 );
4425  assert( propdata->nbinpermvars >= 0 );
4426 
4427  /* if no symmetries on binary variables are present */
4428  if ( propdata->nbinpermvars == 0 )
4429  {
4430  assert( propdata->binvaraffected == 0 );
4431  return SCIP_OKAY;
4432  }
4433 
4434  perms = propdata->perms;
4435  nperms = propdata->nperms;
4436  permvars = propdata->permvars;
4437  npermvars = propdata->npermvars;
4438  conssaddlp = propdata->conssaddlp;
4439  components = propdata->components;
4440  componentbegins = propdata->componentbegins;
4441 
4442  assert( nperms <= 0 || perms != NULL );
4443  assert( permvars != NULL );
4444  assert( npermvars > 0 );
4445  assert( components != NULL );
4446  assert( componentbegins != NULL );
4447  assert( 0 <= cidx && cidx < propdata->ncomponents );
4448 
4449  /* exit if component is already blocked by incompatible methods */
4450  if ( propdata->componentblocked[cidx] & (~SYM_HANDLETYPE_SST) )
4451  return SCIP_OKAY;
4452  if ( (propdata->componentblocked[cidx] & SYM_HANDLETYPE_SST) )
4453  {
4454  /* the leader must be binary for compatability */
4455  if ( (ISSSTINTACTIVE(propdata->sstleadervartype)
4456  || ISSSTIMPLINTACTIVE(propdata->sstleadervartype)
4457  || ISSSTCONTACTIVE(propdata->sstleadervartype)) )
4458  return SCIP_OKAY;
4459  }
4460 
4461  /* skip component if it has signed permutations */
4462  if ( propdata->componenthassignedperm[cidx] )
4463  return SCIP_OKAY;
4464 
4465  /* adapt natural variable order to a variable order that is compatible with Schreier Sims constraints */
4466  if ( propdata->nleaders > 0 && ISSSTBINACTIVE(propdata->sstleadervartype) )
4467  {
4468  SCIP_CALL( SCIPallocBufferArray(scip, &modifiedperms, nperms) );
4469  for (p = 0; p < nperms; ++p)
4470  {
4471  SCIP_CALL( SCIPallocBufferArray(scip, &modifiedperms[p], npermvars) );
4472  }
4473  SCIP_CALL( SCIPallocBufferArray(scip, &modifiedpermvars, npermvars) );
4474 
4475  for (i = 0; i < npermvars; ++i)
4476  modifiedpermvars[i] = permvars[i];
4477 
4478  SCIP_CALL( adaptSymmetryDataSST(scip, perms, modifiedperms, nperms, permvars, modifiedpermvars, npermvars,
4479  propdata->leaders, propdata->nleaders) );
4480  }
4481 
4482  /* loop through perms in component cidx and add symresack constraints */
4483  for (p = componentbegins[cidx]; p < componentbegins[cidx + 1]; ++p)
4484  {
4485  SCIP_CONS* cons;
4486  int permidx;
4487  char name[SCIP_MAXSTRLEN];
4488 
4489  permidx = components[p];
4490 
4491  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "symbreakcons_component%d_perm%d", cidx, permidx);
4492 
4493  /* adapt permutation to leader */
4494  if ( propdata->nleaders > 0 && ISSSTBINACTIVE(propdata->sstleadervartype) )
4495  {
4496  assert( (propdata->componentblocked[cidx] & SYM_HANDLETYPE_SST) != 0 );
4497  assert( modifiedperms != NULL );
4498  assert( modifiedpermvars != NULL );
4499 
4500  SCIP_CALL( SCIPcreateSymbreakCons(scip, &cons, name, modifiedperms[permidx], modifiedpermvars, npermvars, FALSE,
4501  conssaddlp, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
4502  }
4503  else
4504  {
4505  SCIP_CALL( SCIPcreateSymbreakCons(scip, &cons, name, perms[permidx], permvars, npermvars, FALSE,
4506  conssaddlp, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
4507  }
4508  propdata->componentblocked[cidx] |= SYM_HANDLETYPE_SYMBREAK;
4509  SCIP_CALL( SCIPaddCons(scip, cons) );
4510 
4511  /* do not release constraint here - will be done later */
4512  SCIP_CALL( ensureDynamicConsArrayAllocatedAndSufficientlyLarge(scip, &propdata->genorbconss,
4513  &propdata->genorbconsssize, propdata->ngenorbconss + 1) );
4514  propdata->genorbconss[propdata->ngenorbconss++] = cons;
4515  ++propdata->nsymresacks;
4516  ++nsymresackcons;
4517  }
4518 
4519  if ( propdata->nleaders > 0 && ISSSTBINACTIVE(propdata->sstleadervartype) )
4520  {
4521  assert( modifiedperms != NULL );
4522  assert( modifiedpermvars != NULL );
4523 
4524  SCIPfreeBufferArray(scip, &modifiedpermvars);
4525  for (p = nperms - 1; p >= 0; --p)
4526  {
4527  SCIPfreeBufferArray(scip, &modifiedperms[p]);
4528  }
4529  SCIPfreeBufferArray(scip, &modifiedperms);
4530  }
4531 
4532  SCIPdebugMsg(scip, "Added %d symresack constraints.\n", nsymresackcons);
4533 
4534  return SCIP_OKAY;
4535 }
4536 
4537 
4538 /** add Schreier Sims constraints for a specific orbit and update Schreier Sims table */
4539 static
4541  SCIP* scip, /**< SCIP instance */
4542  SCIP_CONFLICTDATA* varconflicts, /**< conflict graph or NULL if useconflictgraph == FALSE */
4543  SCIP_PROPDATA* propdata, /**< data of symmetry propagator */
4544  SCIP_VAR** permvars, /**< permvars array */
4545  int* orbits, /**< symmetry orbits */
4546  int* orbitbegins, /**< array storing begin position for each orbit */
4547  int orbitidx, /**< index of orbit for Schreier Sims constraints */
4548  int orbitleaderidx, /**< index of leader variable for Schreier Sims constraints */
4549  SCIP_Shortbool* orbitvarinconflict, /**< indicator whether orbitvar is in conflict with orbit leader */
4550  int norbitvarinconflict,/**< number of variables in conflict with orbit leader */
4551  int* nchgbds /**< pointer to store number of bound changes (or NULL) */
4552  )
4553 { /*lint --e{613,641}*/
4554  SCIP_CONS* cons;
4555  char name[SCIP_MAXSTRLEN];
4556  SCIP_VAR* vars[2];
4557  SCIP_Real vals[2];
4558  int orbitsize;
4559  int posleader;
4560  int poscur;
4561  int ncuts = 0;
4562  SCIP_Bool addcuts = FALSE;
4563  int i;
4564 #ifndef NDEBUG
4565  int j;
4566 #endif
4567 
4568  assert( scip != NULL );
4569  assert( propdata != NULL );
4570  assert( permvars != NULL );
4571  assert( orbits != NULL );
4572  assert( orbitbegins != NULL );
4573  assert( orbitidx >= 0 );
4574  assert( orbitleaderidx >= 0 );
4575  assert( orbitvarinconflict != NULL || varconflicts == NULL );
4576  assert( norbitvarinconflict >= 0 );
4577  assert( nchgbds != NULL );
4578 
4579  orbitsize = orbitbegins[orbitidx + 1] - orbitbegins[orbitidx];
4580 
4581  /* variables in conflict with leader are fixed and not treated by a cut; trailing -1 to not count the leader */
4582  if ( propdata->sstaddcuts )
4583  addcuts = TRUE;
4584  else if ( propdata->sstleaderrule == SCIP_LEADERRULE_MAXCONFLICTSINORBIT
4585  || propdata->ssttiebreakrule == SCIP_LEADERTIEBREAKRULE_MAXCONFLICTSINORBIT )
4586  addcuts = propdata->addconflictcuts;
4587 
4588  if ( addcuts )
4589  ncuts = orbitsize - norbitvarinconflict - 1;
4590 
4591  /* (re-)allocate memory for Schreier Sims constraints and leaders */
4592  if ( ncuts > 0 )
4593  {
4594  if ( propdata->nsstconss == 0 )
4595  {
4596  assert( propdata->sstconss == NULL );
4597  assert( propdata->maxnsstconss == 0 );
4598  propdata->maxnsstconss = 2 * ncuts;
4599  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(propdata->sstconss), propdata->maxnsstconss) );
4600  }
4601  else if ( propdata->nsstconss + ncuts > propdata->maxnsstconss )
4602  {
4603  int newsize;
4604 
4605  newsize = SCIPcalcMemGrowSize(scip, propdata->maxnsstconss + 2 * ncuts);
4606  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(propdata->sstconss),
4607  propdata->maxnsstconss, newsize) );
4608  propdata->maxnsstconss = newsize;
4609  }
4610  }
4611 
4612  if ( propdata->nleaders == 0 )
4613  {
4614  propdata->maxnleaders = MIN(propdata->nperms, propdata->npermvars);
4615  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(propdata->leaders), propdata->maxnleaders) );
4616  }
4617  assert( propdata->nleaders < propdata->maxnleaders );
4618 
4619  /* add Schreier Sims constraints vars[0] >= vars[1], where vars[0] is always the leader */
4620  posleader = orbitbegins[orbitidx] + orbitleaderidx;
4621  vars[0] = permvars[orbits[posleader]];
4622  vals[0] = -1.0;
4623  vals[1] = 1.0;
4624  propdata->leaders[propdata->nleaders++] = orbits[posleader];
4625  *nchgbds = 0;
4626  for (i = 0, poscur = orbitbegins[orbitidx]; i < orbitsize; ++i, ++poscur)
4627  {
4628  if ( i == orbitleaderidx )
4629  {
4630  assert( orbitvarinconflict == NULL || ! orbitvarinconflict[i] );
4631  continue;
4632  }
4633 
4634  vars[1] = permvars[orbits[poscur]];
4635 #ifndef NDEBUG
4636  for (j = 0; j < propdata->nleaders - 1; ++j)
4637  {
4638  assert( propdata->leaders[j] != orbits[poscur] );
4639  }
4640 #endif
4641 
4642  /* if the i-th variable in the orbit is in a conflict with the leader, fix it to 0 */
4643  if ( varconflicts != NULL )
4644  {
4645  if ( orbitvarinconflict[i] )
4646  {
4647  assert( SCIPvarIsBinary(vars[1]) );
4648  assert( SCIPvarGetLbLocal(vars[1]) < 0.5 );
4649  assert( varconflicts != NULL );
4650 
4651  /* if variable is fixed */
4652  if ( SCIPvarGetUbLocal(vars[1]) > 0.5 )
4653  {
4654  SCIP_CALL( SCIPchgVarUb(scip, vars[1], 0.0) );
4655  ++(*nchgbds);
4656 
4657  /* deactivate the fixed variable (cannot contribute to a conflict anymore) */
4658  assert( varconflicts[orbits[poscur]].active );
4659  varconflicts[orbits[poscur]].active = FALSE;
4660  }
4661 
4662  /* reset value */
4663  orbitvarinconflict[i] = FALSE;
4664  }
4665  else if ( addcuts )
4666  {
4667  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "SSTcut_%d_%d", orbits[posleader], orbits[poscur]);
4668  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, name, 2, vars, vals, - SCIPinfinity(scip), 0.0,
4670 
4671  SCIP_CALL( SCIPaddCons(scip, cons) );
4672  propdata->sstconss[propdata->nsstconss++] = cons;
4673  }
4674  }
4675  else if ( addcuts )
4676  {
4677  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "SSTcut_%d_%d", orbits[posleader], orbits[poscur]);
4678  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, name, 2, vars, vals, - SCIPinfinity(scip), 0.0,
4680 
4681  SCIP_CALL( SCIPaddCons(scip, cons) );
4682  propdata->sstconss[propdata->nsstconss++] = cons;
4683  }
4684  }
4685 
4686  return SCIP_OKAY;
4687 }
4688 
4689 
4690 /** selection rule of next orbit/leader in orbit for Schreier Sims constraints */
4691 static
4693  SCIP* scip, /**< SCIP instance */
4694  SCIP_CONFLICTDATA* varconflicts, /**< variable conflicts structure, or NULL if we do not use it */
4695  SCIP_VAR** conflictvars, /**< variables encoded in conflict graph */
4696  int nconflictvars, /**< number of variables encoded in conflict graph */
4697  int* orbits, /**< orbits of stabilizer subgroup, expressed in terms of conflictvars */
4698  int* orbitbegins, /**< array storing the begin position of each orbit in orbits */
4699  int norbits, /**< number of orbits */
4700  int leaderrule, /**< rule to select leader */
4701  int tiebreakrule, /**< tie break rule to select leader */
4702  SCIP_VARTYPE leadervartype, /**< variable type of leader */
4703  int* orbitidx, /**< pointer to index of selected orbit */
4704  int* leaderidx, /**< pointer to leader in orbit */
4705  SCIP_Shortbool* orbitvarinconflict, /**< array to store whether a var in the orbit is conflicting with leader */
4706  int* norbitvarinconflict,/**< pointer to store number of vars in the orbit in conflict with leader */
4707  SCIP_Bool* success /**< pointer to store whether orbit cut be selected successfully */
4708  )
4709 {
4710  int varidx;
4711  int orbitcriterion;
4712  int curcriterion = INT_MIN;
4713  int orbitsize;
4714  int i;
4715  int leader = -1;
4716 
4717  assert( scip != NULL );
4718  assert( conflictvars != NULL );
4719  assert( nconflictvars > 0 );
4720  assert( orbits != NULL );
4721  assert( orbitbegins != NULL );
4722  assert( norbits > 0 );
4723  assert( orbitidx != NULL );
4724  assert( leaderidx != NULL );
4725  assert( orbitvarinconflict != NULL || varconflicts == NULL );
4726  assert( norbitvarinconflict != NULL );
4727  assert( success != NULL );
4728 
4729  *orbitidx = 0;
4730  *leaderidx = 0;
4731  *norbitvarinconflict = 0;
4732  *success = FALSE;
4733 
4734  /* terminate if leader or tiebreak rule cannot be checked */
4735  if ( varconflicts == NULL && (leaderrule == (int) SCIP_LEADERRULE_MAXCONFLICTSINORBIT
4736  || tiebreakrule == (int) SCIP_LEADERTIEBREAKRULE_MAXCONFLICTSINORBIT) )
4737  return SCIP_OKAY;
4738 
4739  /* select the leader and its orbit */
4740  if ( leaderrule == (int) SCIP_LEADERRULE_FIRSTINORBIT || leaderrule == (int) SCIP_LEADERRULE_LASTINORBIT )
4741  {
4742  orbitcriterion = INT_MIN;
4743 
4744  /* iterate over orbits and select the first one that meets the tiebreak rule */
4745  for (i = 0; i < norbits; ++i)
4746  {
4747  /* skip orbits containing vars different to the leader's vartype */
4748  /* Conflictvars is permvars! */
4749  if ( SCIPvarGetType(conflictvars[orbits[orbitbegins[i]]]) != leadervartype )
4750  continue;
4751 
4752  if ( tiebreakrule == (int) SCIP_LEADERTIEBREAKRULE_MINORBIT )
4753  curcriterion = orbitbegins[i] - orbitbegins[i + 1];
4754  else if ( tiebreakrule == (int) SCIP_LEADERTIEBREAKRULE_MAXORBIT )
4755  curcriterion = orbitbegins[i + 1] - orbitbegins[i];
4756  else
4757  {
4758  assert( tiebreakrule == (int) SCIP_LEADERTIEBREAKRULE_MAXCONFLICTSINORBIT );
4759 
4760  /* get first or last active variable in orbit */
4761  if ( leaderrule == (int) SCIP_LEADERRULE_FIRSTINORBIT )
4762  {
4763  int cnt;
4764 
4765  cnt = orbitbegins[i];
4766 
4767  do
4768  {
4769  varidx = orbits[cnt++];
4770  }
4771  while ( SCIPvarGetProbindex(conflictvars[varidx]) == -1 && cnt < orbitbegins[i + 1]);
4772  }
4773  else
4774  {
4775  int cnt;
4776 
4777  cnt = orbitbegins[i + 1] - 1;
4778 
4779  do
4780  {
4781  varidx = orbits[cnt--];
4782  }
4783  while ( SCIPvarGetProbindex(conflictvars[varidx]) == -1 && cnt >= orbitbegins[i]);
4784  }
4785 
4786  /* skip inactive variables */
4787  if ( SCIPvarGetProbindex(conflictvars[varidx]) == -1 )
4788  continue;
4789 
4790  assert( varconflicts[varidx].orbitidx == i );
4791  /* coverity[var_deref_op] */
4792  curcriterion = varconflicts[varidx].nconflictinorbit;
4793  }
4794 
4795  /* update selected orbit */
4796  if ( curcriterion > orbitcriterion )
4797  {
4798  orbitcriterion = curcriterion;
4799  *orbitidx = i;
4800  *success = TRUE;
4801 
4802  if ( leaderrule == (int) SCIP_LEADERRULE_FIRSTINORBIT )
4803  *leaderidx = 0;
4804  else
4805  *leaderidx = orbitbegins[i + 1] - orbitbegins[i] - 1;
4806  }
4807  }
4808 
4809  /* store variables in conflict with leader */
4810  if ( *success && varconflicts != NULL )
4811  {
4812  leader = orbits[orbitbegins[*orbitidx] + *leaderidx];
4813  assert( leader < nconflictvars );
4814 
4815  if ( tiebreakrule == (int) SCIP_LEADERTIEBREAKRULE_MAXCONFLICTSINORBIT
4816  && varconflicts[leader].ncliques > 0 )
4817  {
4818  /* count how many active variables in the orbit conflict with "leader"
4819  * This is only needed if there are possible conflicts.
4820  */
4821  int varmapid;
4822 
4823  orbitsize = orbitbegins[*orbitidx + 1] - orbitbegins[*orbitidx];
4824  assert( varconflicts != NULL );
4825  assert( leader >= 0 && leader < nconflictvars );
4826 
4827  assert( orbitvarinconflict != NULL );
4828 
4829  for (i = 0; i < orbitsize; ++i)
4830  {
4831  /* skip the leader */
4832  if ( i == *leaderidx )
4833  continue;
4834 
4835  /* get variable index in conflict graph */
4836  varmapid = orbits[orbitbegins[*orbitidx] + i];
4837 
4838  /* only active variables */
4839  if ( ! varconflicts[varmapid].active )
4840  continue;
4841 
4842  /* check if leader and var have overlap */
4843  if ( checkSortedArraysHaveOverlappingEntry((void**)varconflicts[leader].cliques,
4844  varconflicts[leader].ncliques, (void**)varconflicts[varmapid].cliques,
4845  varconflicts[varmapid].ncliques, sortByPointerValue) )
4846  {
4847  /* there is overlap! */
4848  orbitvarinconflict[i] = TRUE;
4849  ++(*norbitvarinconflict);
4850  }
4851  }
4852  }
4853  }
4854  }
4855  else
4856  {
4857  /* only three possible values for leaderrules, so it must be MAXCONFLICTSINORBIT
4858  * In this case, the code must have computed the conflict graph.
4859  */
4860  assert( leaderrule == (int) SCIP_LEADERRULE_MAXCONFLICTSINORBIT );
4861  assert( varconflicts != NULL );
4862 
4863  orbitcriterion = 0;
4864 
4865  /* iterate over variables and select the first one that meets the tiebreak rule */
4866  for (i = 0; i < nconflictvars; ++i)
4867  {
4868  /* skip vars different to the leader's vartype */
4869  if ( SCIPvarGetType(conflictvars[i]) != leadervartype )
4870  continue;
4871 
4872  /* skip variables not affected by symmetry */
4873  /* coverity[var_deref_op] */
4874  if ( varconflicts[i].orbitidx == -1 )
4875  continue;
4876 
4877  curcriterion = varconflicts[i].nconflictinorbit;
4878 
4879  if ( curcriterion > orbitcriterion )
4880  {
4881  orbitcriterion = curcriterion;
4882  *orbitidx = varconflicts[i].orbitidx;
4883  *leaderidx = varconflicts[i].posinorbit;
4884  *success = TRUE;
4885  }
4886  }
4887 
4888  /* store variables in conflict with leader */
4889  leader = orbits[orbitbegins[*orbitidx] + *leaderidx];
4890  assert( leader < nconflictvars );
4891  assert( norbitvarinconflict != NULL );
4892 
4893  if ( *success && varconflicts[leader].ncliques > 0 )
4894  {
4895  /* count how many active variables in the orbit conflict with leader */
4896  int varmapid;
4897 
4898  orbitsize = orbitbegins[*orbitidx + 1] - orbitbegins[*orbitidx];
4899  assert( varconflicts != NULL );
4900  assert( leader >= 0 && leader < nconflictvars );
4901 
4902  assert( orbitvarinconflict != NULL );
4903 
4904  for (i = 0; i < orbitsize; ++i)
4905  {
4906  /* skip the leader */
4907  if ( i == *leaderidx )
4908  continue;
4909 
4910  /* get variable index in conflict graph */
4911  varmapid = orbits[orbitbegins[*orbitidx] + i];
4912  /* only active variables */
4913  if ( ! varconflicts[varmapid].active )
4914  continue;
4915 
4916  /* check if leader and var have overlap */
4917  if ( checkSortedArraysHaveOverlappingEntry((void**)varconflicts[leader].cliques,
4918  varconflicts[leader].ncliques, (void**)varconflicts[varmapid].cliques,
4919  varconflicts[varmapid].ncliques, sortByPointerValue) )
4920  {
4921  /* there is overlap! */
4922  orbitvarinconflict[i] = TRUE;
4923  ++(*norbitvarinconflict);
4924  }
4925  }
4926  }
4927  }
4928 
4929  return SCIP_OKAY;
4930 }
4931 
4932 
4933 /** add Schreier Sims constraints to the problem */
4934 static
4936  SCIP* scip, /**< SCIP instance */
4937  SCIP_PROPDATA* propdata, /**< data of symmetry propagator */
4938  SCIP_Bool onlywithcontvars, /**< only handle components that contain continuous variables with SST */
4939  int* nchgbds, /**< pointer to store number of bound changes (or NULL) */
4940  int cidx /**< index of component which shall be handled */
4941  )
4942 { /*lint --e{641}*/
4943  SCIP_CONFLICTDATA* varconflicts = NULL;
4944  SCIP_HASHMAP* permvarmap;
4945  SCIP_VAR** permvars;
4946  int** permstrans;
4947  int npermvars;
4948  int nmovedpermvars;
4949  int nmovedbinpermvars;
4950  int nmovedintpermvars;
4951  int nmovedimplintpermvars;
4952  int nmovedcontpermvars;
4953  int nperms;
4954 
4955  int* orbits;
4956  int* orbitbegins;
4957  int norbits;
4958  int* components;
4959  int* componentbegins;
4960  int* vartocomponent;
4961  int ncomponents;
4962  unsigned* componentblocked;
4963 
4964  int orbitidx;
4965  int orbitleaderidx;
4966  SCIP_Shortbool* orbitvarinconflict = NULL;
4967  int norbitvarinconflict;
4968  SCIP_Shortbool* inactiveperms;
4969  int ninactiveperms;
4970  int posleader;
4971  int leaderrule;
4972  int tiebreakrule;
4973  int leadervartype;
4974  SCIP_VARTYPE selectedtype = SCIP_VARTYPE_CONTINUOUS;
4975  int nvarsselectedtype;
4976  SCIP_Bool conflictgraphcreated = FALSE;
4977  SCIP_Bool mixedcomponents;
4978  int norbitleadercomponent;
4979  int* perm;
4980  SCIP_VARTYPE vartype;
4981 
4982  int i;
4983  int c;
4984  int p;
4985  SCIP_Bool success = TRUE;
4986 
4987  assert( scip != NULL );
4988  assert( propdata != NULL );
4989  assert( propdata->computedsymmetry );
4990 
4991  permvars = propdata->permvars;
4992  npermvars = propdata->npermvars;
4993  nperms = propdata->nperms;
4994  assert( permvars != NULL );
4995  assert( npermvars > 0 );
4996  assert( nperms > 0 );
4997 
4998  SCIP_CALL( ensureSymmetryPermvarmapComputed(scip, propdata) );
4999  permvarmap = propdata->permvarmap;
5000  assert( permvarmap != NULL );
5001 
5002  SCIP_CALL( ensureSymmetryPermstransComputed(scip, propdata) );
5003  permstrans = propdata->permstrans;
5004  assert( permstrans != NULL );
5005 
5006  components = propdata->components;
5007  componentbegins = propdata->componentbegins;
5008  componentblocked = propdata->componentblocked;
5009  vartocomponent = propdata->vartocomponent;
5010  ncomponents = propdata->ncomponents;
5011 
5012  assert( components != NULL );
5013  assert( componentbegins != NULL );
5014  assert( vartocomponent != NULL );
5015  assert( componentblocked != NULL );
5016  assert( ncomponents > 0 );
5017  assert( 0 <= cidx && cidx < ncomponents );
5018 
5019  /* exit if component is blocked */
5020  if ( componentblocked[cidx] )
5021  return SCIP_OKAY;
5022 
5023  /* skip component if it has signed permutations */
5024  if ( propdata->componenthassignedperm[cidx] )
5025  return SCIP_OKAY;
5026 
5027  leaderrule = propdata->sstleaderrule;
5028  tiebreakrule = propdata->ssttiebreakrule;
5029  leadervartype = propdata->sstleadervartype;
5030  mixedcomponents = propdata->sstmixedcomponents;
5031 
5032  /* if not already computed, get number of affected vars */
5034  nmovedpermvars = propdata->nmovedpermvars;
5035  nmovedbinpermvars = propdata->nmovedbinpermvars;
5036  nmovedintpermvars = propdata->nmovedintpermvars;
5037  nmovedimplintpermvars = propdata->nmovedimplintpermvars;
5038  nmovedcontpermvars = propdata->nmovedcontpermvars;
5039  assert( nmovedpermvars > 0 ); /* nperms > 0 implies this */
5040 
5041  /* determine the leader's vartype */
5042  nvarsselectedtype = 0;
5043  if ( ISSSTBINACTIVE(leadervartype) && nmovedbinpermvars > nvarsselectedtype )
5044  {
5045  selectedtype = SCIP_VARTYPE_BINARY;
5046  nvarsselectedtype = nmovedbinpermvars;
5047  }
5048 
5049  if ( ISSSTINTACTIVE(leadervartype) && nmovedintpermvars > nvarsselectedtype )
5050  {
5051  selectedtype = SCIP_VARTYPE_INTEGER;
5052  nvarsselectedtype = nmovedintpermvars;
5053  }
5054 
5055  if ( ISSSTIMPLINTACTIVE(leadervartype) && nmovedimplintpermvars > nvarsselectedtype )
5056  {
5057  selectedtype = SCIP_VARTYPE_IMPLINT;
5058  nvarsselectedtype = nmovedimplintpermvars;
5059  }
5060 
5061  if ( ISSSTCONTACTIVE(leadervartype) && nmovedcontpermvars > nvarsselectedtype )
5062  {
5063  selectedtype = SCIP_VARTYPE_CONTINUOUS;
5064  nvarsselectedtype = nmovedcontpermvars;
5065  }
5066 
5067  /* terminate if no variables of a possible leader type is affected */
5068  if ( nvarsselectedtype == 0 )
5069  return SCIP_OKAY;
5070 
5071  /* ignore this component if no continuous variables are contained */
5072  if ( onlywithcontvars )
5073  {
5074  for (p = componentbegins[cidx]; p < componentbegins[cidx + 1]; ++p)
5075  {
5076  perm = propdata->perms[p];
5077  for (i = 0; i < propdata->npermvars; ++i)
5078  {
5079  if ( perm[i] == i )
5080  continue;
5081  vartype = SCIPvarGetType(propdata->permvars[i]);
5082  if ( vartype == SCIP_VARTYPE_CONTINUOUS || vartype == SCIP_VARTYPE_IMPLINT )
5083  goto COMPONENTOK;
5084  }
5085  }
5086  /* loop terminated naturally, so component does not have continuous or implicitly integer variables. */
5087  return SCIP_OKAY;
5088 
5089  COMPONENTOK:
5090  ;
5091  }
5092 
5093  /* @todo online create the conflict graph for the variable in the current component */
5094  /* possibly create conflict graph; graph is not created if no cliques are present */
5095  if ( selectedtype == SCIP_VARTYPE_BINARY && (leaderrule == SCIP_LEADERRULE_MAXCONFLICTSINORBIT
5096  || tiebreakrule == SCIP_LEADERTIEBREAKRULE_MAXCONFLICTSINORBIT) )
5097  {
5098  SCIP_CALL( createConflictGraphSST(scip, &varconflicts, permvars, npermvars, permvarmap) );
5099  conflictgraphcreated = varconflicts != NULL;
5100  }
5101 
5102  /* allocate data structures necessary for orbit computations and conflict graph */
5103  SCIP_CALL( SCIPallocBufferArray(scip, &inactiveperms, nperms) );
5104  SCIP_CALL( SCIPallocBufferArray(scip, &orbits, npermvars) );
5105  SCIP_CALL( SCIPallocBufferArray(scip, &orbitbegins, npermvars) );
5106 
5107  if ( conflictgraphcreated )
5108  {
5109  SCIP_CALL( SCIPallocClearBufferArray(scip, &orbitvarinconflict, npermvars) );
5110  }
5111 
5112  SCIPdebugMsg(scip, "Start selection of orbits and leaders for Schreier Sims constraints.\n");
5113  SCIPdebugMsg(scip, "orbitidx\tleaderidx\torbitsize\n");
5114 
5115  if ( nchgbds != NULL )
5116  *nchgbds = 0;
5117 
5118  /* initialize array indicating whether permutations shall not be considered for orbit permutations */
5119  for (c = 0; c < ncomponents; ++c)
5120  {
5121  for (p = componentbegins[c]; p < componentbegins[c + 1]; ++p)
5122  {
5123  if ( c == cidx )
5124  inactiveperms[components[p]] = FALSE;
5125  else
5126  inactiveperms[components[p]] = TRUE;
5127  }
5128  }
5129  ninactiveperms = nperms - componentbegins[cidx + 1] + componentbegins[cidx];
5130 
5131  /* as long as the stabilizer is non-trivial, add Schreier Sims constraints */
5132  norbitleadercomponent = 0;
5133  while ( ninactiveperms < nperms )
5134  {
5135  int nchanges = 0;
5136 
5137  /* compute orbits w.r.t. active perms */
5138  SCIP_CALL( SCIPcomputeOrbitsFilterSym(scip, npermvars, permstrans, nperms, inactiveperms,
5139  orbits, orbitbegins, &norbits, components, componentbegins, vartocomponent,
5140  componentblocked, ncomponents, nmovedpermvars) );
5141 
5142  /* stop if we require pure components and a component contains variables of different types */
5143  if ( ! mixedcomponents )
5144  {
5145  for (p = 0; p < norbits; ++p)
5146  {
5147  /* stop if the first element of an orbits has the wrong vartype */
5148  if ( SCIPvarGetType(permvars[orbits[orbitbegins[p]]]) != selectedtype )
5149  {
5150  success = FALSE;
5151  break;
5152  }
5153  }
5154  }
5155 
5156  if ( ! success )
5157  break;
5158 
5159  /* update symmetry information of conflict graph */
5160  if ( conflictgraphcreated )
5161  {
5162  SCIP_CALL( updateSymInfoConflictGraphSST(scip, varconflicts, permvars, npermvars, orbits, orbitbegins,
5163  norbits) );
5164  }
5165 
5166  /* possibly adapt the leader and tie-break rule */
5167  if ( leaderrule == SCIP_LEADERRULE_MAXCONFLICTSINORBIT && ! conflictgraphcreated )
5168  leaderrule = SCIP_LEADERRULE_FIRSTINORBIT;
5169  if ( leaderrule == SCIP_LEADERRULE_MAXCONFLICTSINORBIT && selectedtype != SCIP_VARTYPE_BINARY )
5170  leaderrule = SCIP_LEADERRULE_FIRSTINORBIT;
5171  if ( tiebreakrule == SCIP_LEADERTIEBREAKRULE_MAXCONFLICTSINORBIT && ! conflictgraphcreated )
5172  tiebreakrule = SCIP_LEADERTIEBREAKRULE_MAXORBIT;
5173  if ( tiebreakrule == SCIP_LEADERTIEBREAKRULE_MAXCONFLICTSINORBIT && selectedtype != SCIP_VARTYPE_BINARY )
5174  tiebreakrule = SCIP_LEADERTIEBREAKRULE_MAXORBIT;
5175 
5176  /* select orbit and leader */
5177  SCIP_CALL( selectOrbitLeaderSSTConss(scip, varconflicts, permvars, npermvars, orbits, orbitbegins,
5178  norbits, propdata->sstleaderrule, propdata->ssttiebreakrule, selectedtype, &orbitidx, &orbitleaderidx,
5179  orbitvarinconflict, &norbitvarinconflict, &success) );
5180 
5181  if ( ! success )
5182  break;
5183 
5184  assert( 0 <= orbitidx && orbitidx < norbits );
5185  assert( 0 <= orbitleaderidx && orbitleaderidx < orbitbegins[orbitidx + 1] - orbitbegins[orbitidx] );
5186  SCIPdebugMsg(scip, "%d\t\t%d\t\t%d\n", orbitidx, orbitleaderidx, orbitbegins[orbitidx + 1] - orbitbegins[orbitidx]);
5187 
5188  /* add Schreier Sims constraints for the selected orbit and update Schreier Sims table */
5189  SCIP_CALL( addSSTConssOrbitAndUpdateSST(scip, varconflicts, propdata, permvars,
5190  orbits, orbitbegins, orbitidx, orbitleaderidx, orbitvarinconflict, norbitvarinconflict, &nchanges) );
5191 
5192  ++norbitleadercomponent;
5193 
5194  if ( nchgbds != NULL )
5195  *nchgbds += nchanges;
5196 
5197  /* deactivate permutations that move the orbit leader */
5198  posleader = orbits[orbitbegins[orbitidx] + orbitleaderidx];
5199  for (p = 0; p < nperms; ++p)
5200  {
5201  if ( inactiveperms[p] )
5202  continue;
5203 
5204  if ( permstrans[posleader][p] != posleader )
5205  {
5206  inactiveperms[p] = TRUE;
5207  ++ninactiveperms;
5208  }
5209  }
5210  }
5211 
5212  /* if Schreier Sims constraints have been added, store that Schreier Sims has been used for this component */
5213  if ( norbitleadercomponent > 0 )
5214  componentblocked[cidx] |= SYM_HANDLETYPE_SST;
5215 
5216  if ( conflictgraphcreated )
5217  {
5218  SCIPfreeBufferArray(scip, &orbitvarinconflict);
5219  }
5220  SCIPfreeBufferArray(scip, &orbitbegins);
5221  SCIPfreeBufferArray(scip, &orbits);
5222  if ( varconflicts != NULL )
5223  {
5224  /* nconflictvars at construction is npermvars */
5225  SCIP_CALL( freeConflictGraphSST(scip, &varconflicts, npermvars) );
5226  }
5227  SCIPfreeBufferArray(scip, &inactiveperms);
5228 
5229  return SCIP_OKAY;
5230 }
5231 
5232 
5233 /** orbitopal reduction */
5234 static
5236  SCIP* scip, /**< SCIP instance */
5237  SCIP_PROPDATA* propdata, /**< propdata */
5238  int id, /**< ID for orbitope constraint (needed for name) */
5239  int** varidxmatrix, /**< matrix containing variable indices in orbitope matrix */
5240  int nrows, /**< number of rows of orbitope */
5241  int ncols, /**< number of columns of orbitope */
5242  SCIP_Bool* success /**< pointer to store whether orbitope could be added successfully */
5243  )
5244 {
5245  char name[SCIP_MAXSTRLEN];
5246  int i;
5247  int j;
5248 
5249  SCIP_Bool ispporbitope;
5250  SCIP_VAR*** varmatrix;
5251  SCIP_Bool* pprows;
5252  int npprows;
5253  SCIP_ORBITOPETYPE type;
5254 
5255  assert( scip != NULL );
5256  assert( propdata != NULL );
5257  assert( propdata->usedynamicprop );
5258  assert( varidxmatrix != NULL );
5259  assert( nrows > 0 );
5260  assert( ncols > 0 );
5261  assert( success != NULL );
5262 
5263  *success = FALSE;
5264 
5265  /* add linear constraints x_1 >= x_2 >= ... >= x_ncols for single-row orbitopes */
5266  if ( nrows == 1 )
5267  {
5268  /* restrict to the packing and partitioning rows */
5269  SCIP_CONS* cons;
5270  SCIP_VAR* consvars[2];
5271  SCIP_Real conscoefs[2] = { -1.0, 1.0 };
5272 
5273  /* for all adjacent column pairs, add linear constraint */
5274  SCIP_CALL( ensureDynamicConsArrayAllocatedAndSufficientlyLarge(scip, &propdata->genlinconss,
5275  &propdata->genlinconsssize, propdata->ngenlinconss + ncols - 1) );
5276  for (i = 0; i < ncols - 1; ++i)
5277  {
5278  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "orbitope_1row_comp_%d_col%d", id, i);
5279 
5280  consvars[0] = propdata->permvars[varidxmatrix[0][i]];
5281  consvars[1] = propdata->permvars[varidxmatrix[0][i + 1]];
5282 
5283  /* enforce, but do not check */
5284  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, name, 2, consvars, conscoefs, -SCIPinfinity(scip), 0.0,
5285  propdata->conssaddlp, propdata->conssaddlp, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE ) );
5286 
5287  SCIP_CALL( SCIPaddCons(scip, cons) );
5288  propdata->genlinconss[propdata->ngenlinconss++] = cons;
5289  }
5290 
5291  *success = TRUE;
5292  return SCIP_OKAY;
5293  }
5294 
5295  /* for only 2 columns, the the component can be completely handled by lexicographic reduction */
5296  if ( ncols == 2 && propdata->lexreddata != NULL )
5297  {
5298  int* orbisackperm;
5299 
5300  /* If the component is an orbitope with 2 columns, then there is 1 generator of order 2. */
5301  orbisackperm = propdata->perms[propdata->components[propdata->componentbegins[id]]];
5302 
5303  SCIP_CALL( SCIPlexicographicReductionAddPermutation(scip, propdata->lexreddata,
5304  propdata->permvars, propdata->npermvars, orbisackperm, (SYM_SYMTYPE) propdata->symtype,
5305  propdata->permvardomaincenter, TRUE, success) );
5306  if ( *success )
5307  return SCIP_OKAY;
5308  }
5309 
5310  /* create orbitope variable matrix */
5311  SCIP_CALL( SCIPallocBufferArray(scip, &varmatrix, nrows) );
5312  for (i = 0; i < nrows; ++i)
5313  {
5314  SCIP_CALL( SCIPallocBufferArray(scip, &varmatrix[i], ncols) );
5315  for (j = 0; j < ncols; ++j)
5316  varmatrix[i][j] = propdata->permvars[varidxmatrix[i][j]];
5317  }
5318 
5319  pprows = NULL;
5320  SCIP_CALL( SCIPisPackingPartitioningOrbitope(scip, varmatrix, nrows, ncols, &pprows, &npprows, &type) );
5321 
5322  /* does it have at least 3 packing-partitioning rows? */
5323  ispporbitope = npprows >= 3; /* (use same magic number as cons_orbitope.c) */
5324 
5325  if ( ispporbitope ) /* @todo if it's a pporbitope, we do it statically right now. */
5326  {
5327  /* restrict to the packing and partitioning rows */
5328  SCIP_CONS* cons;
5329  SCIP_VAR*** ppvarsarrayonlypprows;
5330  int r;
5331 
5332  assert( pprows != NULL );
5333 
5334  SCIP_CALL( SCIPallocBufferArray(scip, &ppvarsarrayonlypprows, npprows) );
5335 
5336  r = 0;
5337  for (i = 0; i < nrows; ++i)
5338  {
5339  if ( pprows[i] )
5340  {
5341  assert( r < npprows );
5342  ppvarsarrayonlypprows[r++] = varmatrix[i];
5343  }
5344  }
5345  assert( r == npprows );
5346 
5347  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "orbitope_pp_comp_%d", id);
5348  SCIP_CALL( SCIPcreateConsOrbitope(scip, &cons, name, ppvarsarrayonlypprows, SCIP_ORBITOPETYPE_PACKING,
5349  npprows, ncols, FALSE, FALSE, FALSE, FALSE, propdata->conssaddlp,
5350  TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
5351 
5352  SCIP_CALL( SCIPaddCons(scip, cons) );
5353 
5354  /* check whether we need to resize */
5355  SCIP_CALL( ensureDynamicConsArrayAllocatedAndSufficientlyLarge(scip, &propdata->genlinconss,
5356  &propdata->genlinconsssize, propdata->ngenlinconss + 1) );
5357  /* @todo we add orbitopes to the dynamically sized array `genlinconss` instead of `genorbconss` to ensure
5358  * compatability with the static orbitope function, which allocates this array statically
5359  */
5360  propdata->genlinconss[propdata->ngenlinconss++] = cons;
5361  *success = TRUE;
5362 
5363  SCIPfreeBufferArray(scip, &ppvarsarrayonlypprows);
5364  }
5365  else
5366  {
5367  /* use orbitopal reduction for component */
5368  SCIP_COLUMNORDERING columnordering;
5369  SCIP_VAR** orbitopevarmatrix;
5370  int nelem;
5371  int pos = 0;
5372 
5373  /* variable array */
5374  nelem = nrows * ncols;
5375  SCIP_CALL( SCIPallocBufferArray(scip, &orbitopevarmatrix, nelem) );
5376  for (i = 0; i < nrows; ++i)
5377  {
5378  for (j = 0; j < ncols; ++j)
5379  orbitopevarmatrix[pos++] = varmatrix[i][j];
5380  }
5381 
5382  /* get column ordering */
5383  columnordering = SCIPorbitopalReductionGetDefaultColumnOrdering(propdata->orbitopalreddata);
5384 
5385  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "orbitope_full_comp_%d", id);
5386  SCIP_CALL( SCIPorbitopalReductionAddOrbitope(scip, propdata->orbitopalreddata,
5387  SCIP_ROWORDERING_BRANCHING, columnordering,
5388  orbitopevarmatrix, nrows, ncols, success) );
5389  *success = TRUE;
5390 
5391  SCIPfreeBufferArray(scip, &orbitopevarmatrix);
5392  }
5393 
5394  SCIPfreeBlockMemoryArrayNull(scip, &pprows, nrows);
5395 
5396  for (i = nrows - 1; i >= 0; --i)
5397  {
5398  SCIPfreeBufferArray(scip, &varmatrix[i]);
5399  }
5400  SCIPfreeBufferArray(scip, &varmatrix);
5401 
5402  return SCIP_OKAY;
5403 }
5404 
5405 
5406 /** applies pp-orbitope upgrade if at least 50% of the permutations in a component correspond to pp-orbisacks */
5407 static
5409  SCIP* scip, /**< SCIP instance */
5410  SCIP_PROPDATA* propdata, /**< propdata */
5411  int** componentperms, /**< permutations in the component */
5412  int componentsize, /**< number of permutations in the component */
5413  SCIP_Bool hassignedperm, /**< whether the component has a signed permutation */
5414  SCIP_Bool* success /**< whether the packing partitioning upgrade succeeded */
5415  )
5416 {
5417  int c;
5418  int i;
5419  int j;
5420  int p;
5421  int* perm;
5422  SCIP_CONSHDLR* setppcconshdlr;
5423  SCIP_CONS** setppcconss;
5424  SCIP_CONS* cons;
5425  SCIP_CONS** setppconsssort;
5426  int nsetppconss;
5427  int nsetppcvars;
5428  SCIP_VAR** setppcvars;
5429  int nsetppcconss;
5430  int** pporbisackperms;
5431  int npporbisackperms;
5432  SCIP_VAR* var;
5433  int varid;
5434  SCIP_CONS*** permvarssetppcconss;
5435  int* npermvarssetppcconss;
5436  int* maxnpermvarssetppcconss;
5437  int maxntwocycles;
5438  int ntwocycles;
5439 
5440  assert( scip != NULL );
5441  assert( propdata != NULL );
5442  assert( componentperms != NULL );
5443  assert( componentsize > 0 );
5444  assert( success != NULL );
5445 
5446  /* we did not upgrade yet */
5447  *success = FALSE;
5448 
5449  /* currently, we cannot handle signed permutations */
5450  if ( hassignedperm )
5451  return SCIP_OKAY;
5452 
5453  setppcconshdlr = SCIPfindConshdlr(scip, "setppc");
5454  if ( setppcconshdlr == NULL )
5455  return SCIP_OKAY;
5456 
5457  nsetppcconss = SCIPconshdlrGetNConss(setppcconshdlr);
5458  if ( nsetppcconss == 0 )
5459  return SCIP_OKAY;
5460 
5461  setppcconss = SCIPconshdlrGetConss(setppcconshdlr);
5462  assert( setppcconss != NULL );
5463 
5464  SCIP_CALL( ensureSymmetryPermvarmapComputed(scip, propdata) );
5465 
5466  /* collect non-covering constraints and sort by pointer for easy intersection finding */
5467  SCIP_CALL( SCIPallocBufferArray(scip, &setppconsssort, nsetppcconss) );
5468  nsetppconss = 0;
5469  for (c = 0; c < nsetppcconss; ++c)
5470  {
5471  cons = setppcconss[c];
5472 
5473  /* only packing or partitioning constraints, no covering types */
5474  if ( SCIPgetTypeSetppc(scip, cons) == SCIP_SETPPCTYPE_COVERING )
5475  continue;
5476 
5477  setppconsssort[nsetppconss++] = cons;
5478  }
5479  SCIPsortPtr((void**) setppconsssort, sortByPointerValue, nsetppcconss);
5480 
5481  /* For each permvar, introduce an array of setppc constraints (initially NULL) for each variable,
5482  * and populate it with the setppc constraints that it contains. This array follows the ordering by cons ptr address.
5483  */
5484  SCIP_CALL( SCIPallocCleanBufferArray(scip, &permvarssetppcconss, propdata->npermvars) );
5485  SCIP_CALL( SCIPallocCleanBufferArray(scip, &npermvarssetppcconss, propdata->npermvars) );
5486  SCIP_CALL( SCIPallocCleanBufferArray(scip, &maxnpermvarssetppcconss, propdata->npermvars) );
5487  for (c = 0; c < nsetppconss; ++c)
5488  {
5489  assert( c >= 0 );
5490  assert( c < nsetppconss );
5491  cons = setppconsssort[c];
5492  assert( cons != NULL );
5493 
5494  setppcvars = SCIPgetVarsSetppc(scip, cons);
5495  nsetppcvars = SCIPgetNVarsSetppc(scip, cons);
5496 
5497  for (i = 0; i < nsetppcvars; ++i)
5498  {
5499  var = setppcvars[i];
5500  assert( var != NULL );
5501  varid = SCIPhashmapGetImageInt(propdata->permvarmap, (void*) var);
5502  assert( varid == INT_MAX || varid < propdata->npermvars );
5503  assert( varid >= 0 );
5504  if ( varid < propdata->npermvars )
5505  {
5507  &(permvarssetppcconss[varid]), &maxnpermvarssetppcconss[varid], npermvarssetppcconss[varid] + 1) );
5508  assert( npermvarssetppcconss[varid] < maxnpermvarssetppcconss[varid] );
5509  permvarssetppcconss[varid][npermvarssetppcconss[varid]++] = cons;
5510  }
5511  }
5512  }
5513 
5514  /* for all permutations, test involutions on binary variables and test if they are captured by setppc conss */
5515  SCIP_CALL( SCIPallocBufferArray(scip, &pporbisackperms, componentsize) );
5516  maxntwocycles = 0;
5517  npporbisackperms = 0;
5518  for (p = 0; p < componentsize; ++p)
5519  {
5520  perm = componentperms[p];
5521  ntwocycles = 0;
5522 
5523  /* check if the binary orbits are involutions */
5524  for (i = 0; i < propdata->npermvars; ++i)
5525  {
5526  j = perm[i];
5527 
5528  /* ignore fixed points in permutation */
5529  if ( i == j )
5530  continue;
5531  /* only check for situations where i and j are binary variables */
5532  assert( SCIPvarGetType(propdata->permvars[i]) == SCIPvarGetType(propdata->permvars[j]) );
5533  if ( SCIPvarGetType(propdata->permvars[i]) != SCIP_VARTYPE_BINARY )
5534  continue;
5535  /* the permutation must be an involution on binary variables */
5536  if ( perm[j] != i )
5537  goto NEXTPERMITER;
5538  /* i and j are a two-cycle, so we find this once for i and once for j. Only handle this once for i < j. */
5539  if ( i > j )
5540  continue;
5541  /* disqualify permutation if i and j are not in a common set packing constraint */
5542  if ( !checkSortedArraysHaveOverlappingEntry((void**) permvarssetppcconss[i], npermvarssetppcconss[i],
5543  (void**) permvarssetppcconss[j], npermvarssetppcconss[j], sortByPointerValue) )
5544  goto NEXTPERMITER;
5545  ++ntwocycles;
5546  }
5547 
5548  /* The permutation qualifies if all binary variables are either a reflection or in a 2-cycle. There must be at
5549  * least one binary 2-cycle, because otherwise the permutation is the identity, or it permutes
5550  * nonbinary variables.
5551  */
5552  if ( ntwocycles > 0 )
5553  {
5554  pporbisackperms[npporbisackperms++] = perm;
5555  if ( ntwocycles > maxntwocycles )
5556  maxntwocycles = ntwocycles;
5557  }
5558 
5559  NEXTPERMITER:
5560  ;
5561  }
5562 
5563  /* if at least 50% of such permutations are packing-partitioning type, apply packing upgrade */
5564  if ( npporbisackperms * 2 >= componentsize )
5565  {
5566  char name[SCIP_MAXSTRLEN];
5567  SCIP_VAR** ppvarsblock;
5568  SCIP_VAR*** ppvarsmatrix;
5569  SCIP_VAR** row;
5570  int nrows;
5571 
5572  assert( npporbisackperms > 0 );
5573  assert( maxntwocycles > 0 );
5574 
5575  /* instead of allocating and re-allocating multiple times, recycle the ppvars array */
5576  SCIP_CALL( SCIPallocBufferArray(scip, &ppvarsblock, 2 * maxntwocycles) );
5577  SCIP_CALL( SCIPallocBufferArray(scip, &ppvarsmatrix, maxntwocycles) );
5578  for (i = 0; i < maxntwocycles; ++i)
5579  ppvarsmatrix[i] = &(ppvarsblock[2 * i]);
5580 
5581  /* for each of these perms, create the packing orbitope matrix and add constraint*/
5582  for (p = 0; p < npporbisackperms; ++p)
5583  {
5584  perm = pporbisackperms[p];
5585 
5586  /* populate ppvarsmatrix */
5587  nrows = 0;
5588  for (i = 0; i < propdata->npermvars; ++i)
5589  {
5590  j = perm[i];
5591 
5592  /* ignore fixed points in permutation, and only consider rows with i < j */
5593  if ( i >= j )
5594  continue;
5595  /* only for situations where i and j are binary variables */
5596  assert( SCIPvarGetType(propdata->permvars[i]) == SCIPvarGetType(propdata->permvars[j]) );
5597  if ( SCIPvarGetType(propdata->permvars[i]) != SCIP_VARTYPE_BINARY )
5598  continue;
5599  assert( perm[j] == i );
5600  assert( checkSortedArraysHaveOverlappingEntry((void**) permvarssetppcconss[i], npermvarssetppcconss[i],
5601  (void**) permvarssetppcconss[j], npermvarssetppcconss[j], sortByPointerValue) );
5602 
5603  assert( nrows < maxntwocycles );
5604  row = ppvarsmatrix[nrows++];
5605  row[0] = propdata->permvars[i];
5606  row[1] = propdata->permvars[j];
5607  assert( row[0] != row[1] );
5608  }
5609  assert( nrows > 0 );
5610 
5611  /* create constraint, use same parameterization as in orbitope packing partitioning checker */
5612  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "orbitope_pp_upgrade_lexred%d", p);
5613  SCIP_CALL( SCIPcreateConsOrbitope(scip, &cons, name, ppvarsmatrix, SCIP_ORBITOPETYPE_PACKING, nrows, 2,
5614  FALSE, FALSE, FALSE, FALSE,
5615  propdata->conssaddlp, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
5616 
5617  SCIP_CALL( ensureDynamicConsArrayAllocatedAndSufficientlyLarge(scip, &propdata->genlinconss,
5618  &propdata->genlinconsssize, propdata->ngenlinconss + 1) );
5619  /* @todo we add orbitopes to the dynamically sized array `genlinconss` instead of `genorbconss` to ensure
5620  * compatability with the static orbitope function, which allocates this array statically
5621  */
5622  propdata->genlinconss[propdata->ngenlinconss++] = cons;
5623  SCIP_CALL( SCIPaddCons(scip, cons) );
5624  }
5625 
5626  SCIPfreeBufferArray(scip, &ppvarsmatrix);
5627  SCIPfreeBufferArray(scip, &ppvarsblock);
5628 
5629  *success = TRUE;
5630  }
5631 
5632  /* free pp orbisack array */
5633  SCIPfreeBufferArray(scip, &pporbisackperms);
5634 
5635  /* clean the non-clean arrays */
5636  for (varid = 0; varid < propdata->npermvars; ++varid)
5637  {
5638  assert( (permvarssetppcconss[varid] == NULL) == (maxnpermvarssetppcconss[varid] == 0) );
5639  assert( npermvarssetppcconss[varid] >= 0 );
5640  assert( maxnpermvarssetppcconss[varid] >= 0 );
5641  assert( npermvarssetppcconss[varid] <= maxnpermvarssetppcconss[varid] );
5642  if ( npermvarssetppcconss[varid] == 0 )
5643  continue;
5644  SCIPfreeBlockMemoryArray(scip, &permvarssetppcconss[varid], maxnpermvarssetppcconss[varid]);
5645  permvarssetppcconss[varid] = NULL;
5646  npermvarssetppcconss[varid] = 0;
5647  maxnpermvarssetppcconss[varid] = 0;
5648  }
5649  SCIPfreeCleanBufferArray(scip, &maxnpermvarssetppcconss);
5650  SCIPfreeCleanBufferArray(scip, &npermvarssetppcconss);
5651  SCIPfreeCleanBufferArray(scip, &permvarssetppcconss);
5652  SCIPfreeBufferArray(scip, &setppconsssort);
5653 
5654  return SCIP_OKAY;
5655 }
5656 
5657 
5658 /** dynamic permutation lexicographic reduction */
5659 static
5661  SCIP* scip, /**< SCIP instance */
5662  SCIP_PROPDATA* propdata, /**< propdata */
5663  int cidx /**< index of component */
5664  )
5665 {
5666  int componentsize;
5667  int** componentperms;
5668  int p;
5669 
5670  SCIP_Bool checkorbired;
5671  SCIP_Bool checklexred;
5672  SCIP_Bool success;
5673  SCIP_PARAM* checkpporbisack;
5674 
5675  assert( scip != NULL );
5676  assert( propdata != NULL );
5677  assert( ISORBITALREDUCTIONACTIVE(propdata->usesymmetry)
5678  || (
5679  ISSYMRETOPESACTIVE(propdata->usesymmetry)
5680  && propdata->usedynamicprop
5681  && propdata->addsymresacks
5682  ) );
5683  assert( propdata->nperms > 0 );
5684  assert( 0 <= cidx && cidx < propdata->ncomponents );
5685  assert( propdata->componentblocked != NULL );
5686 
5687  /* exit if component is already blocked */
5688  if ( propdata->componentblocked[cidx] )
5689  return SCIP_OKAY;
5690 
5691  /* in this function orbital reduction or dynamic lexicographic reduction propagation must be enabled */
5692  checkorbired = ISORBITALREDUCTIONACTIVE(propdata->usesymmetry);
5693  checklexred = ISSYMRETOPESACTIVE(propdata->usesymmetry) && propdata->usedynamicprop && propdata->addsymresacks;
5694  assert( checkorbired || checklexred );
5695 
5697  assert( propdata->nmovedpermvars );
5698 
5699  /* collect the permutations of this component */
5700  componentsize = propdata->componentbegins[cidx + 1] - propdata->componentbegins[cidx];
5701  SCIP_CALL( SCIPallocBufferArray(scip, &componentperms, componentsize) );
5702  for (p = 0; p < componentsize; ++p)
5703  componentperms[p] = propdata->perms[propdata->components[propdata->componentbegins[cidx] + p]];
5704 
5705  /* check if many component permutations contain many packing partitioning orbisacks
5706  *
5707  * 1. Get the checkpporbisack param from the parameter hashset. This returns NULL if it is not initialized,
5708  * likely because the orbisack constraint handler is not loaded.
5709  * 2. If the param is not NULL, then we only do the packing-partitioning upgrade step if its value is TRUE.
5710  * Packing-partitioning orbitopes are only implemented for binary orbitopes, so binary variables must be moved.
5711  */
5712  checkpporbisack = SCIPgetParam(scip, "constraints/orbisack/checkpporbisack");
5713  if ( ( checkpporbisack == NULL || SCIPparamGetBool(checkpporbisack) == TRUE ) && propdata->nmovedbinpermvars > 0 )
5714  {
5716  componentperms, componentsize, propdata->componenthassignedperm[cidx], &success) );
5717 
5718  if ( success )
5719  {
5720  propdata->componentblocked[cidx] |= SYM_HANDLETYPE_SYMBREAK;
5721  goto FINISHCOMPONENT;
5722  }
5723  }
5724 
5725  /* handle component permutations with orbital reduction */
5726  if ( checkorbired && !propdata->componenthassignedperm[cidx] )
5727  {
5728  SCIP_CALL( SCIPorbitalReductionAddComponent(scip, propdata->orbitalreddata,
5729  propdata->permvars, propdata->npermvars, componentperms, componentsize, &success) );
5730  if ( success )
5731  propdata->componentblocked[cidx] |= SYM_HANDLETYPE_ORBITALREDUCTION;
5732  }
5733 
5734  /* handle component permutations with the dynamic lexicographic reduction propagator */
5735  if ( checklexred )
5736  {
5737  /* handle every permutation in the component with the dynamic lexicographic reduction propagator */
5738  for (p = 0; p < componentsize; ++p)
5739  {
5740  assert( componentperms[p] != NULL );
5741  SCIP_CALL( SCIPlexicographicReductionAddPermutation(scip, propdata->lexreddata,
5742  propdata->permvars, propdata->npermvars, componentperms[p],
5743  (SYM_SYMTYPE) propdata->symtype, propdata->permvardomaincenter, TRUE, &success) );
5744  if ( success )
5745  propdata->componentblocked[cidx] |= SYM_HANDLETYPE_SYMBREAK;
5746  }
5747  }
5748 
5749  FINISHCOMPONENT:
5750  /* if it got blocked here */
5751  if ( propdata->componentblocked[cidx] )
5752  ++propdata->ncompblocked;
5753 
5754  SCIPfreeBufferArray(scip, &componentperms);
5755 
5756  return SCIP_OKAY;
5757 }
5758 
5759 
5760 /** displays statistics on the used symmetry handling methods */
5761 static
5763  SCIP* scip, /**< SCIP instance */
5764  SCIP_PROPDATA* propdata /**< data of symmetry propagator */
5765  )
5766 {
5767  int ncomponentshandled;
5768  int i;
5769 
5770  assert( scip != NULL );
5771  assert( propdata != NULL );
5772 
5773  SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, "dynamic symmetry handling statistics:\n");
5774  if ( propdata->orbitopalreddata )
5775  {
5776  SCIP_CALL( SCIPorbitopalReductionPrintStatistics(scip, propdata->orbitopalreddata) );
5777  }
5778  if ( propdata->orbitalreddata )
5779  {
5780  SCIP_CALL( SCIPorbitalReductionPrintStatistics(scip, propdata->orbitalreddata) );
5781  }
5782  if ( propdata->lexreddata )
5783  {
5784  SCIP_CALL( SCIPlexicographicReductionPrintStatistics(scip, propdata->lexreddata) );
5785  }
5786  if ( propdata->ncomponents >= 0 )
5787  {
5788  /* report the number of handled components
5789  *
5790  * Since SST is compatible with static symresacks, the propdata->ncompblocked counter is not the number of
5791  * handled components. Compute this statistic based on the componentblocked array.
5792  */
5793  ncomponentshandled = 0;
5794  for (i = 0; i < propdata->ncomponents; ++i)
5795  {
5796  if ( propdata->componentblocked[i] )
5797  ++ncomponentshandled;
5798  }
5799  assert( propdata->ncompblocked <= ncomponentshandled );
5800  assert( ncomponentshandled <= propdata->ncomponents );
5801  SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, "handled %d out of %d symmetry components\n",
5802  ncomponentshandled, propdata->ncomponents);
5803  }
5804 
5805  return SCIP_OKAY;
5806 }
5807 
5808 /** handles orbitope action by static or dynamic symmetry handling methods */
5809 static
5811  SCIP* scip, /**< SCIP instance */
5812  SCIP_PROPDATA* propdata, /**< data of symmetry propagator */
5813  int id, /**< ID of orbitope (used for constraint name) */
5814  int** varidxmatrix, /**< matrix containing variable indices of orbitope */
5815  int nrows, /**< number of rows of matrix */
5816  int ncols, /**< number of columns of matrix */
5817  SCIP_Bool* success /**< pointer to store whether orbitope could be added successfully */
5818  )
5819 {
5820  assert( scip != NULL );
5821  assert( propdata != NULL );
5822  assert( varidxmatrix != NULL );
5823  assert( nrows > 0 );
5824  assert( ncols > 0 );
5825  assert( success != NULL );
5826 
5827  *success = FALSE;
5828 
5829  /* dynamic propagation */
5830  if ( propdata->usedynamicprop )
5831  {
5832  SCIP_CALL( addOrbitopesDynamic(scip, propdata, id, varidxmatrix, nrows, ncols, success) );
5833  }
5834  /* static variant only for binary variables */
5835  else if ( propdata->binvaraffected )
5836  {
5837  char name[SCIP_MAXSTRLEN];
5838  SCIP_VAR*** orbitopematrix;
5839  SCIP_CONS* cons;
5840  int i;
5841  int j;
5842  int nbinrows = 0;
5843 
5844  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "orbitope_component_%d", id);
5845 
5846  SCIP_CALL( SCIPallocBufferArray(scip, &orbitopematrix, nrows) );
5847  for (i = 0; i < nrows; ++i)
5848  {
5849  /* skip rows without binary variables */
5850  if ( ! SCIPvarIsBinary(propdata->permvars[varidxmatrix[i][0]]) )
5851  continue;
5852 
5853  SCIP_CALL( SCIPallocBufferArray(scip, &orbitopematrix[nbinrows], ncols) );
5854  for (j = 0; j < ncols; ++j)
5855  {
5856  assert( SCIPvarIsBinary(propdata->permvars[varidxmatrix[i][j]]) );
5857  orbitopematrix[nbinrows][j] = propdata->permvars[varidxmatrix[i][j]];
5858  }
5859  ++nbinrows;
5860  }
5861 
5862  if ( nbinrows > 0 )
5863  {
5864  SCIP_CALL( SCIPcreateConsOrbitope(scip, &cons, name, orbitopematrix, SCIP_ORBITOPETYPE_FULL,
5865  nbinrows, ncols, propdata->usedynamicprop /* @todo disable */, FALSE, FALSE, FALSE,
5866  propdata->conssaddlp, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
5867 
5868  SCIP_CALL( SCIPaddCons(scip, cons) );
5869 
5870  /* do not release constraint here - will be done later */
5871  SCIP_CALL( ensureDynamicConsArrayAllocatedAndSufficientlyLarge(scip, &propdata->genorbconss,
5872  &propdata->genorbconsssize, propdata->ngenorbconss + 1) );
5873  propdata->genorbconss[propdata->ngenorbconss++] = cons;
5874  ++propdata->norbitopes;
5875 
5876  *success = TRUE;
5877  }
5878 
5879  for (i = nbinrows - 1; i >= 0; --i)
5880  {
5881  SCIPfreeBufferArray(scip, &orbitopematrix[i]);
5882  }
5883  SCIPfreeBufferArray(scip, &orbitopematrix);
5884  }
5885 
5886  return SCIP_OKAY;
5887 }
5888 
5889 /** handles binary double lex matrix by adding static orbitope constraints
5890  *
5891  * @todo Extend method to general variable types and dynamic variable orders.
5892  */
5893 static
5895  SCIP* scip, /**< SCIP instance */
5896  SCIP_PROPDATA* propdata, /**< data of symmetry propagator */
5897  int id, /**< ID of double lex matrix (used for constraint names) */
5898  int** varidxmatrix, /**< matrix containing variable indices of double lex matrix */
5899  int nrows, /**< number of rows of matrix */
5900  int ncols, /**< number of columns of matrix */
5901  int* rowsbegin, /**< array indicating where a new row block begins */
5902  int* colsbegin, /**< array indicating where a new column block begins */
5903  int nrowblocks, /**< number of row blocks */
5904  int ncolblocks, /**< number of column blocks */
5905  SCIP_Bool* success /**< pointer to store whether orbitope could be added successfully */
5906  )
5907 {
5908  char name[SCIP_MAXSTRLEN];
5909  SCIP_VAR*** orbitopematrix;
5910  SCIP_CONS* cons;
5911  int maxdim;
5912  int i;
5913  int p;
5914  int j;
5915  int col;
5916  int nbinrows;
5917 
5918  assert( scip != NULL );
5919  assert( propdata != NULL );
5920  assert( varidxmatrix != NULL );
5921  assert( nrows > 0 );
5922  assert( ncols > 0 );
5923  assert( rowsbegin != NULL );
5924  assert( colsbegin != NULL );
5925  assert( nrowblocks > 0 );
5926  assert( ncolblocks > 0 );
5927  assert( success != NULL );
5928 
5929  /* ensure that we can store orbitope constraints in probdata */
5930  SCIP_CALL( ensureDynamicConsArrayAllocatedAndSufficientlyLarge(scip, &propdata->genorbconss,
5931  &propdata->genorbconsssize, propdata->ngenorbconss + nrowblocks + ncolblocks) );
5932 
5933  maxdim = MAX(nrows, ncols);
5934  SCIP_CALL( SCIPallocBufferArray(scip, &orbitopematrix, maxdim) );
5935  for (i = 0; i < maxdim; ++i)
5936  {
5937  SCIP_CALL( SCIPallocBufferArray(scip, &orbitopematrix[i], maxdim) );
5938  }
5939 
5940  /* add orbitopes corresponding to column blocks of doublelexmatrix */
5941  for (p = 0; p < ncolblocks; ++p)
5942  {
5943  nbinrows = 0;
5944  for (i = 0; i < nrows; ++i)
5945  {
5946  /* skip rows that do not contain binary variables */
5947  if ( ! SCIPvarIsBinary(propdata->permvars[varidxmatrix[i][colsbegin[p]]]) )
5948  continue;
5949 
5950  for (col = 0, j = colsbegin[p]; j < colsbegin[p + 1]; ++j, ++col)
5951  {
5952  assert( SCIPvarIsBinary(propdata->permvars[varidxmatrix[i][j]]) );
5953  orbitopematrix[nbinrows][col] = propdata->permvars[varidxmatrix[i][j]];
5954  }
5955  ++nbinrows;
5956  }
5957 
5958  if ( nbinrows > 0 )
5959  {
5960  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "doublelex_cols_%d_%d", id, p);
5961  SCIP_CALL( SCIPcreateConsOrbitope(scip, &cons, name, orbitopematrix, SCIP_ORBITOPETYPE_FULL,
5962  nrows, colsbegin[p + 1] - colsbegin[p], FALSE, FALSE, TRUE, FALSE,
5963  propdata->conssaddlp, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
5964  SCIP_CALL( SCIPaddCons(scip, cons) );
5965  propdata->genorbconss[(propdata->ngenorbconss)++] = cons;
5966  /* do not release constraint here - will be done later */
5967  }
5968  }
5969 
5970  /* add orbitopes corresponding to row blocks of doublelexmatrix */
5971  for (p = 0; p < nrowblocks; ++p)
5972  {
5973  nbinrows = 0;
5974  for (i = 0; i < ncols; ++i)
5975  {
5976  /* skip rows that do not contain binary variables */
5977  if ( ! SCIPvarIsBinary(propdata->permvars[varidxmatrix[rowsbegin[p]][i]]) )
5978  continue;
5979 
5980  for (col = 0, j = rowsbegin[p]; j < rowsbegin[p + 1]; ++j, ++col)
5981  {
5982  assert( SCIPvarIsBinary(propdata->permvars[varidxmatrix[j][i]]) );
5983  orbitopematrix[nbinrows][col] = propdata->permvars[varidxmatrix[j][i]];
5984  }
5985  ++nbinrows;
5986  }
5987 
5988  if ( nbinrows > 0 )
5989  {
5990  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "doublelex_rows_%d_%d", id, p);
5991  SCIP_CALL( SCIPcreateConsOrbitope(scip, &cons, name, orbitopematrix, SCIP_ORBITOPETYPE_FULL,
5992  ncols, rowsbegin[p + 1] - rowsbegin[p], FALSE, FALSE, TRUE, FALSE,
5993  propdata->conssaddlp, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
5994  SCIP_CALL( SCIPaddCons(scip, cons) );
5995  propdata->genorbconss[(propdata->ngenorbconss)++] = cons;
5996  /* do not release constraint here - will be done later */
5997  }
5998  }
5999 
6000  for (i = maxdim - 1; i >= 0; --i)
6001  {
6002  SCIPfreeBufferArray(scip, &orbitopematrix[i]);
6003  }
6004  SCIPfreeBufferArray(scip, &orbitopematrix);
6005 
6006  return SCIP_OKAY;
6007 }
6008 
6009 /** tries to handle symmetries of single lex matrices (orbitopes) or double lex matrices */
6010 static
6012  SCIP* scip, /**< SCIP instance */
6013  SCIP_PROPDATA* propdata, /**< data of symmetry propagator */
6014  SCIP_Bool detectsinglelex, /**< whether single lex matrices shall be detected */
6015  int cidx /**< index of component */
6016  )
6017 {
6018  int** lexmatrix = NULL;
6019  int* lexrowsbegin = NULL;
6020  int* lexcolsbegin = NULL;
6021  int nrows;
6022  int ncols;
6023  int nrowmatrices;
6024  int ncolmatrices;
6025  int** perms;
6026  int compsize;
6027  int i;
6028  int p;
6029  SCIP_Bool isorbitope;
6030  SCIP_Bool success = FALSE;
6031 
6032  assert( scip != NULL );
6033  assert( propdata != NULL );
6034  assert( 0 <= cidx && cidx < propdata->ncomponents );
6035 
6036  /* exit if component is already blocked */
6037  if ( propdata->componentblocked[cidx] )
6038  return SCIP_OKAY;
6039 
6040  /* exit if component has non-standard permutations */
6041  if ( propdata->componenthassignedperm[cidx] )
6042  return SCIP_OKAY;
6043 
6044  /* exit if polyhedral methods are disabled when looking for double lex matrices */
6045  if ( !ISSYMRETOPESACTIVE(propdata->usesymmetry) && !detectsinglelex )
6046  return SCIP_OKAY;
6047 
6048  /* get permutations of component */
6049  compsize = propdata->componentbegins[cidx + 1] - propdata->componentbegins[cidx];
6050  SCIP_CALL( SCIPallocBufferArray(scip, &perms, compsize) );
6051  for (p = 0, i = propdata->componentbegins[cidx]; i < propdata->componentbegins[cidx + 1]; ++i)
6052  perms[p++] = propdata->perms[propdata->components[i]];
6053 
6054  SCIP_CALL( SCIPdetectSingleOrDoubleLexMatrices(scip, detectsinglelex, perms, compsize, propdata->npermvars,
6055  &success, &isorbitope, &lexmatrix, &nrows, &ncols,
6056  &lexrowsbegin, &lexcolsbegin, &nrowmatrices, &ncolmatrices) );
6057 
6058  SCIPfreeBufferArray(scip, &perms);
6059 
6060  /* possibly handle double lex matrix or orbitope */
6061  if ( success )
6062  {
6063  assert( lexmatrix != NULL );
6064  assert( nrows > 0 );
6065  assert( ncols > 0 );
6066 
6067  if ( isorbitope )
6068  {
6069  SCIP_CALL( handleOrbitope(scip, propdata, cidx, lexmatrix, nrows, ncols, &success) );
6070  }
6071  else
6072  {
6073  SCIP_Bool hasbinaryvar = FALSE;
6074 
6075  /* check whether a binary variable is contained in the matrix */
6076  for (i = 0; i < nrows && !hasbinaryvar; ++i)
6077  {
6078  for (p = 0; p < ncols; ++p)
6079  {
6080  if ( SCIPvarIsBinary(propdata->permvars[lexmatrix[i][p]]) )
6081  {
6082  hasbinaryvar = TRUE;
6083  break;
6084  }
6085  }
6086  }
6087 
6088  if ( hasbinaryvar )
6089  {
6090  SCIP_CALL( handleDoublelLexMatrix(scip, propdata, cidx, lexmatrix, nrows, ncols,
6091  lexrowsbegin, lexcolsbegin, nrowmatrices, ncolmatrices, &success) );
6092  }
6093  else
6094  success = FALSE;
6095  }
6096 
6097  /* free memory not needed anymore */
6098