1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB) */
7 /* */
8 /* Licensed under the Apache License, Version 2.0 (the "License"); */
9 /* you may not use this file except in compliance with the License. */
10 /* You may obtain a copy of the License at */
11 /* */
12 /* http://www.apache.org/licenses/LICENSE-2.0 */
13 /* */
14 /* Unless required by applicable law or agreed to in writing, software */
15 /* distributed under the License is distributed on an "AS IS" BASIS, */
16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17 /* See the License for the specific language governing permissions and */
18 /* limitations under the License. */
19 /* */
20 /* You should have received a copy of the Apache-2.0 license */
21 /* along with SCIP; see the file LICENSE. If not visit scipopt.org. */
22 /* */
23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24
25 /**@file symmetry.c
26 * @ingroup OTHER_CFILES
27 * @brief methods for handling symmetries
28 * @author Jasper van Doornmalen
29 * @author Christopher Hojny
30 * @author Marc Pfetsch
31 */
32
33 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
34
35 #include "scip/symmetry.h"
36 #include "scip/scip.h"
37 #include "scip/cons_setppc.h"
38 #include "scip/cons_orbitope.h"
39 #include "scip/misc.h"
40 #include "symmetry/struct_symmetry.h"
41 #include "symmetry/type_symmetry.h"
42
43
44 /** compute non-trivial orbits of symmetry group
45 *
46 * The non-trivial orbits of the group action are stored in the array orbits of length npermvars. This array contains
47 * the indices of variables from the permvars array such that variables that are contained in the same orbit appear
48 * consecutively in the orbits array. The variables of the i-th orbit have indices
49 * orbits[orbitbegins[i]], ... , orbits[orbitbegins[i + 1] - 1].
50 * Note that the description of the orbits ends at orbitbegins[norbits] - 1.
51 */
52 SCIP_RETCODE SCIPcomputeOrbitsSym(
53 SCIP* scip, /**< SCIP instance */
54 SCIP_Bool issigend, /**< whether orbits for signed permutations shall be computed */
55 SCIP_VAR** permvars, /**< variables considered in a permutation */
56 int npermvars, /**< length of a permutation array */
57 int** perms, /**< matrix containing in each row a permutation of the symmetry group */
58 int nperms, /**< number of permutations encoded in perms */
59 int* orbits, /**< array of non-trivial orbits */
60 int* orbitbegins, /**< array containing begin positions of new orbits in orbits array */
61 int* norbits /**< pointer to number of orbits currently stored in orbits */
62 )
63 {
64 SCIP_Shortbool* varadded;
65 int permlen;
66 int orbitidx = 0;
67 int i;
68
69 assert( scip != NULL );
70 assert( permvars != NULL );
71 assert( perms != NULL );
72 assert( nperms > 0 );
73 assert( npermvars > 0 );
74 assert( orbits != NULL );
75 assert( orbitbegins != NULL );
76 assert( norbits != NULL );
77
78 permlen = npermvars;
79 if ( issigend )
80 permlen *= 2;
81
82 /* init data structures*/
83 SCIP_CALL( SCIPallocBufferArray(scip, &varadded, permlen) );
84
85 /* initially, every variable is contained in no orbit */
86 for (i = 0; i < permlen; ++i)
87 varadded[i] = FALSE;
88
89 /* find variable orbits */
90 *norbits = 0;
91 for (i = 0; i < permlen; ++i)
92 {
93 int beginorbitidx;
94 int j;
95
96 /* skip variable already contained in an orbit of a previous variable */
97 if ( varadded[i] )
98 continue;
99
100 /* store first variable */
101 beginorbitidx = orbitidx;
102 orbits[orbitidx++] = i;
103 varadded[i] = TRUE;
104
105 /* iterate over variables in curorbit and compute their images */
106 j = beginorbitidx;
107 while ( j < orbitidx )
108 {
109 int curelem;
110 int image;
111 int p;
112
113 curelem = orbits[j];
114
115 for (p = 0; p < nperms; ++p)
116 {
117 image = perms[p][curelem];
118
119 /* found new element of the orbit of i */
120 if ( ! varadded[image] )
121 {
122 orbits[orbitidx++] = image;
123 assert( orbitidx <= permlen );
124 varadded[image] = TRUE;
125 }
126 }
127 ++j;
128 }
129
130 /* if the orbit is trivial, reset storage, otherwise store orbit */
131 if ( orbitidx <= beginorbitidx + 1 )
132 orbitidx = beginorbitidx;
133 else
134 orbitbegins[(*norbits)++] = beginorbitidx;
135 }
136
137 /* store end in "last" orbitbegins entry */
138 assert( *norbits < permlen );
139 orbitbegins[*norbits] = orbitidx;
140
141 #ifdef SCIP_OUTPUT
142 printf("Orbits (total number: %d):\n", *norbits);
143 for (i = 0; i < *norbits; ++i)
144 {
145 int j;
146
147 printf("%d: ", i);
148 for (j = orbitbegins[i]; j < orbitbegins[i+1]; ++j)
149 printf("%d ", orbits[j]);
150 printf("\n");
151 }
152 #endif
153
154 /* free memory */
155 SCIPfreeBufferArray(scip, &varadded);
156
157 return SCIP_OKAY;
158 }
159
160
161 /** compute non-trivial orbits of symmetry group using filtered generators
162 *
163 * The non-trivial orbits of the group action are stored in the array orbits of length npermvars. This array contains
164 * the indices of variables from the permvars array such that variables that are contained in the same orbit appear
165 * consecutively in the orbits array. The variables of the i-th orbit have indices
166 * orbits[orbitbegins[i]], ... , orbits[orbitbegins[i + 1] - 1].
167 * Note that the description of the orbits ends at orbitbegins[norbits] - 1.
168 *
169 * Only permutations that are not inactive (as marked by @p inactiveperms) are used. Thus, one can use this array to
170 * filter out permutations.
171 */
172 SCIP_RETCODE SCIPcomputeOrbitsFilterSym(
173 SCIP* scip, /**< SCIP instance */
174 int npermvars, /**< length of a permutation array */
175 int** permstrans, /**< transposed matrix containing in each column a
176 * permutation of the symmetry group */
177 int nperms, /**< number of permutations encoded in perms */
178 SCIP_Shortbool* inactiveperms, /**< array to store whether permutations are inactive */
179 int* orbits, /**< array of non-trivial orbits */
180 int* orbitbegins, /**< array containing begin positions of new orbits in orbits array */
181 int* norbits, /**< pointer to number of orbits currently stored in orbits */
182 int* components, /**< array containing the indices of permutations sorted by components */
183 int* componentbegins, /**< array containing in i-th position the first position of
184 * component i in components array */
185 int* vartocomponent, /**< array containing for each permvar the index of the component it is
186 * contained in (-1 if not affected) */
187 unsigned* componentblocked, /**< array to store which symmetry methods have been used on a component
188 * using the same bitset information as for misc/usesymmetry */
189 int ncomponents, /**< number of components of symmetry group */
190 int nmovedpermvars /**< number of variables moved by any permutation in a symmetry component
191 * that is handled by orbital fixing */
192 )
193 {
194 SCIP_Shortbool* varadded;
195 int nvaradded = 0;
196 int orbitidx = 0;
197 int i;
198
199 assert( scip != NULL );
200 assert( permstrans != NULL );
201 assert( nperms > 0 );
202 assert( npermvars > 0 );
203 assert( inactiveperms != NULL );
204 assert( orbits != NULL );
205 assert( orbitbegins != NULL );
206 assert( norbits != NULL );
207 assert( components != NULL );
208 assert( componentbegins != NULL );
209 assert( vartocomponent != NULL );
210 assert( ncomponents > 0 );
211 assert( nmovedpermvars > 0 );
212
213 /* init data structures */
214 SCIP_CALL( SCIPallocBufferArray(scip, &varadded, npermvars) );
215
216 /* initially, every variable is contained in no orbit */
217 for (i = 0; i < npermvars; ++i)
218 varadded[i] = FALSE;
219
220 /* find variable orbits */
221 *norbits = 0;
222 for (i = 0; i < npermvars; ++i)
223 {
224 int beginorbitidx;
225 int componentidx;
226 int j;
227
228 /* skip unaffected variables and blocked components */
229 componentidx = vartocomponent[i];
230 if ( componentidx < 0 || componentblocked[componentidx] )
231 continue;
232
233 /* skip variable already contained in an orbit of a previous variable */
234 if ( varadded[i] )
235 continue;
236
237 /* store first variable */
238 beginorbitidx = orbitidx;
239 orbits[orbitidx++] = i;
240 varadded[i] = TRUE;
241 ++nvaradded;
242
243 /* iterate over variables in curorbit and compute their images */
244 j = beginorbitidx;
245 while ( j < orbitidx )
246 {
247 int* pt;
248 int curelem;
249 int image;
250 int p;
251
252 curelem = orbits[j];
253
254 pt = permstrans[curelem];
255 for (p = componentbegins[componentidx]; p < componentbegins[componentidx + 1]; ++p)
256 {
257 int perm;
258
259 perm = components[p];
260
261 if ( ! inactiveperms[perm] )
262 {
263 image = pt[perm];
264 assert( vartocomponent[image] == componentidx );
265
266 /* found new element of the orbit of i */
267 if ( ! varadded[image] )
268 {
269 orbits[orbitidx++] = image;
270 assert( orbitidx <= npermvars );
271 varadded[image] = TRUE;
272 ++nvaradded;
273 }
274 }
275 }
276 ++j;
277 }
278
279 /* if the orbit is trivial, reset storage, otherwise store orbit */
280 if ( orbitidx <= beginorbitidx + 1 )
281 orbitidx = beginorbitidx;
282 else
283 orbitbegins[(*norbits)++] = beginorbitidx;
284
285 /* stop if all variables are covered */
286 if ( nvaradded >= nmovedpermvars )
287 break;
288 }
289
290 /* store end in "last" orbitbegins entry */
291 assert( *norbits < npermvars );
292 orbitbegins[*norbits] = orbitidx;
293
294 #ifdef SCIP_OUTPUT
295 printf("Orbits (total number: %d):\n", *norbits);
296 for (i = 0; i < *norbits; ++i)
297 {
298 int j;
299
300 printf("%d: ", i);
301 for (j = orbitbegins[i]; j < orbitbegins[i+1]; ++j)
302 printf("%d ", orbits[j]);
303 printf("\n");
304 }
305 #endif
306
307 /* free memory */
308 SCIPfreeBufferArray(scip, &varadded);
309
310 return SCIP_OKAY;
311 }
312
313 /** Compute orbit of a given variable and store it in @p orbit. The first entry of the orbit will
314 * be the given variable index and the rest is filled with the remaining variables excluding
315 * the ones specified in @p ignoredvars.
316 *
317 * @pre orbit is an initialized array of size propdata->npermvars
318 * @pre at least one of @p perms and @p permstrans should not be NULL
319 */
320 SCIP_RETCODE SCIPcomputeOrbitVar(
321 SCIP* scip, /**< SCIP instance */
322 int npermvars, /**< number of variables in permvars */
323 int** perms, /**< the generators of the permutation group (or NULL) */
324 int** permstrans, /**< the transposed matrix of generators (or NULL) */
325 int* components, /**< the components of the permutation group */
326 int* componentbegins, /**< array containing the starting index of each component */
327 SCIP_Shortbool* ignoredvars, /**< array indicating which variables should be ignored */
328 SCIP_Shortbool* varfound, /**< bitmap to mark which variables have been added (or NULL) */
329 int varidx, /**< index of variable for which the orbit is requested */
330 int component, /**< component that var is in */
331 int* orbit, /**< array in which the orbit should be stored */
332 int* orbitsize /**< buffer to store the size of the orbit */
333 )
334 { /*lint --e{571}*/
335 SCIP_Shortbool* varadded;
336 int* varstotest;
337 int nvarstotest;
338 int j;
339 int p;
340
341 assert( scip != NULL );
342 assert( perms != NULL || permstrans != NULL );
343 assert( components != NULL );
344 assert( componentbegins != NULL );
345 assert( ignoredvars != NULL );
346 assert( orbit != NULL );
347 assert( orbitsize != NULL );
348 assert( 0 <= varidx && varidx < npermvars );
349 assert( component >= 0 );
350 assert( npermvars > 0 );
351
352 /* init data structures*/
353 SCIP_CALL( SCIPallocClearBufferArray(scip, &varadded, npermvars) );
354 SCIP_CALL( SCIPallocClearBufferArray(scip, &varstotest, npermvars) );
355
356 /* compute and store orbit if it is non-trivial */
357 orbit[0] = varidx;
358 varstotest[0] = varidx;
359 *orbitsize = 1;
360 nvarstotest = 1;
361 varadded[varidx] = TRUE;
362
363 if ( varfound != NULL )
364 varfound[varidx] = TRUE;
365
366 /* iterate over variables in orbit and compute their images */
367 j = 0;
368 while ( j < nvarstotest )
369 {
370 int currvar;
371
372 currvar = varstotest[j++];
373
374 for (p = componentbegins[component]; p < componentbegins[component+1]; ++p)
375 {
376 int image;
377 int comp;
378
379 comp = components[p];
380
381 if ( perms != NULL )
382 image = perms[comp][currvar]; /*lint !e613*/
383 else
384 image = permstrans[currvar][comp];
385
386 /* found new element of the orbit of varidx */
387 if ( ! varadded[image] )
388 {
389 varstotest[nvarstotest++] = image;
390 varadded[image] = TRUE;
391
392 if ( ! ignoredvars[image] )
393 {
394 orbit[(*orbitsize)++] = image;
395
396 if ( varfound != NULL )
397 varfound[image] = TRUE;
398 }
399 }
400 }
401 }
402
403 /* free memory */
404 SCIPfreeBufferArray(scip, &varstotest);
405 SCIPfreeBufferArray(scip, &varadded);
406
407 return SCIP_OKAY;
408 }
409
410 /** compute non-trivial orbits of symmetry group
411 *
412 * The non-trivial orbits of the group action are stored in the array orbits of length npermvars. This array contains
413 * the indices of variables from the permvars array such that variables that are contained in the same orbit appear
414 * consecutively in the orbits array. The variables of the i-th orbit have indices
415 * orbits[orbitbegins[i]], ... , orbits[orbitbegins[i + 1] - 1].
416 * Note that the description of the orbits ends at orbitbegins[norbits] - 1.
417 *
418 * This function is adapted from computeGroupOrbitsFilter().
419 */
420 SCIP_RETCODE SCIPcomputeOrbitsComponentsSym(
421 SCIP* scip, /**< SCIP instance */
422 int npermvars, /**< length of a permutation array */
423 int** permstrans, /**< transposed matrix containing in each column a permutation of the symmetry group */
424 int nperms, /**< number of permutations encoded in perms */
425 int* components, /**< array containing the indices of permutations sorted by components */
426 int* componentbegins, /**< array containing in i-th position the first position of component i in components array */
427 int* vartocomponent, /**< array containing for each permvar the index of the component it is
428 * contained in (-1 if not affected) */
429 int ncomponents, /**< number of components of symmetry group */
430 int* orbits, /**< array of non-trivial orbits */
431 int* orbitbegins, /**< array containing begin positions of new orbits in orbits array */
432 int* norbits, /**< pointer to number of orbits currently stored in orbits */
433 int* varorbitmap /**< array for storing the orbits for each variable */
434 )
435 {
436 SCIP_Shortbool* varadded;
437 int orbitidx = 0;
438 int i;
439
440 assert( scip != NULL );
441 assert( permstrans != NULL );
442 assert( nperms > 0 );
443 assert( npermvars > 0 );
444 assert( components != NULL );
445 assert( componentbegins != NULL );
446 assert( vartocomponent != NULL );
447 assert( ncomponents > 0 );
448 assert( orbits != NULL );
449 assert( orbitbegins != NULL );
450 assert( norbits != NULL );
451 assert( varorbitmap != NULL );
452
453 /* init data structures */
454 SCIP_CALL( SCIPallocBufferArray(scip, &varadded, npermvars) );
455
456 /* initially, every variable is contained in no orbit */
457 for (i = 0; i < npermvars; ++i)
458 {
459 varadded[i] = FALSE;
460 varorbitmap[i] = -1;
461 }
462
463 /* find variable orbits */
464 *norbits = 0;
465 for (i = 0; i < npermvars; ++i)
466 {
467 int beginorbitidx;
468 int componentidx;
469 int j;
470
471 /* skip unaffected variables - note that we also include blocked components */
472 componentidx = vartocomponent[i];
473 if ( componentidx < 0 )
474 continue;
475
476 /* skip variable already contained in an orbit of a previous variable */
477 if ( varadded[i] )
478 continue;
479
480 /* store first variable */
481 beginorbitidx = orbitidx;
482 orbits[orbitidx++] = i;
483 varadded[i] = TRUE;
484 varorbitmap[i] = *norbits;
485
486 /* iterate over variables in curorbit and compute their images */
487 j = beginorbitidx;
488 while ( j < orbitidx )
489 {
490 int* pt;
491 int curelem;
492 int image;
493 int p;
494
495 curelem = orbits[j];
496
497 pt = permstrans[curelem];
498 for (p = componentbegins[componentidx]; p < componentbegins[componentidx + 1]; ++p)
499 {
500 int perm;
501
502 perm = components[p];
503 image = pt[perm];
504 assert( vartocomponent[image] == componentidx );
505
506 /* found new element of the orbit of i */
507 if ( ! varadded[image] )
508 {
509 orbits[orbitidx++] = image;
510 assert( orbitidx <= npermvars );
511 varadded[image] = TRUE;
512 varorbitmap[image] = *norbits;
513 }
514 }
515 ++j;
516 }
517
518 /* if the orbit is trivial, reset storage, otherwise store orbit */
519 if ( orbitidx <= beginorbitidx + 1 )
520 {
521 orbitidx = beginorbitidx;
522 varorbitmap[i] = -1;
523 }
524 else
525 orbitbegins[(*norbits)++] = beginorbitidx;
526 }
527
528 /* store end in "last" orbitbegins entry */
529 assert( *norbits < npermvars );
530 orbitbegins[*norbits] = orbitidx;
531
532 /* free memory */
533 SCIPfreeBufferArray(scip, &varadded);
534
535 return SCIP_OKAY;
536 }
537
538
539 /** Checks whether a permutation is a composition of 2-cycles and in this case determines the number of overall
540 * 2-cycles and binary 2-cycles. It is a composition of 2-cycles iff @p ntwocyclesperm > 0 upon termination.
541 */
542 SCIP_RETCODE SCIPisInvolutionPerm(
543 int* perm, /**< permutation */
544 SCIP_VAR** vars, /**< array of variables perm is acting on */
545 int nvars, /**< number of variables */
546 int* ntwocyclesperm, /**< pointer to store number of 2-cycles or 0 if perm is not an involution */
547 int* nbincyclesperm, /**< pointer to store number of binary cycles */
548 SCIP_Bool earlytermination /**< whether we terminate early if not all affected variables are binary */
549 )
550 {
551 int ntwocycles = 0;
552 int i;
553
554 assert( perm != NULL );
555 assert( vars != NULL );
556 assert( ntwocyclesperm != NULL );
557 assert( nbincyclesperm != NULL );
558
559 *ntwocyclesperm = 0;
560 *nbincyclesperm = 0;
561 for (i = 0; i < nvars; ++i)
562 {
563 assert( 0 <= perm[i] && perm[i] < nvars );
564
565 /* skip fixed points and avoid treating the same 2-cycle twice */
566 if ( perm[i] <= i )
567 continue;
568
569 if ( perm[perm[i]] == i )
570 {
571 if ( SCIPvarIsBinary(vars[i]) && SCIPvarIsBinary(vars[perm[i]]) )
572 ++(*nbincyclesperm);
573 else if ( earlytermination )
574 return SCIP_OKAY;
575
576 ++ntwocycles;
577 }
578 else
579 {
580 /* we do not have only 2-cycles */
581 return SCIP_OKAY;
582 }
583 }
584
585 /* at this point the permutation is a composition of 2-cycles */
586 *ntwocyclesperm = ntwocycles;
587
588 return SCIP_OKAY;
589 }
590
591
592 /** determine number of variables affected by symmetry group */
593 SCIP_RETCODE SCIPdetermineNVarsAffectedSym(
594 SCIP* scip, /**< SCIP instance */
595 int** perms, /**< permutations */
596 int nperms, /**< number of permutations in perms */
597 SCIP_VAR** permvars, /**< variables corresponding to permutations */
598 int npermvars, /**< number of permvars in perms */
599 int* nvarsaffected /**< pointer to store number of all affected variables */
600 )
601 {
602 SCIP_Shortbool* affected;
603 int i;
604 int p;
605
606 assert( scip != NULL );
607 assert( perms != NULL );
608 assert( nperms > 0 );
609 assert( permvars != NULL );
610 assert( npermvars > 0 );
611 assert( nvarsaffected != NULL );
612
613 *nvarsaffected = 0;
614
615 SCIP_CALL( SCIPallocClearBufferArray(scip, &affected, npermvars) );
616
617 /* iterate over permutations and check which variables are affected by some symmetry */
618 for (p = 0; p < nperms; ++p)
619 {
620 for (i = 0; i < npermvars; ++i)
621 {
622 if ( affected[i] )
623 continue;
624
625 if ( perms[p][i] != i )
626 {
627 affected[i] = TRUE;
628 ++(*nvarsaffected);
629 }
630 }
631 }
632 SCIPfreeBufferArray(scip, &affected);
633
634 return SCIP_OKAY;
635 }
636
637
638 /** Given a matrix with nrows and \#perms + 1 columns whose first nfilledcols columns contain entries of variables, this routine
639 * checks whether the 2-cycles of perm intersect each row of column coltoextend in exactly one position. In this case,
640 * we add one column to the suborbitope of the first nfilledcols columns.
641 *
642 * @pre Every non-trivial cycle of perm is a 2-cycle.
643 * @pre perm has nrows many 2-cycles
644 */
645 SCIP_RETCODE SCIPextendSubOrbitope(
646 int** suborbitope, /**< matrix containing suborbitope entries */
647 int nrows, /**< number of rows of suborbitope */
648 int nfilledcols, /**< number of columns of suborbitope which are filled with entries */
649 int coltoextend, /**< index of column that should be extended by perm */
650 int* perm, /**< permutation */
651 SCIP_Bool leftextension, /**< whether we extend the suborbitope to the left */
652 int** nusedelems, /**< pointer to array storing how often an element was used in the orbitope */
653 SCIP_VAR** permvars, /**< permutation vars array */
654 SCIP_Shortbool* rowisbinary, /**< array encoding whether variables in an orbitope row are binary (or NULL) */
655 SCIP_Bool* success, /**< pointer to store whether extension was successful */
656 SCIP_Bool* infeasible /**< pointer to store if the number of intersecting cycles is too small */
657 )
658 {
659 int nintersections = 0;
660 int row;
661 int idx1;
662 int idx2;
663
664 assert( suborbitope != NULL );
665 assert( nrows > 0 );
666 assert( nfilledcols > 0 );
667 assert( coltoextend >= 0 );
668 assert( perm != NULL );
669 assert( nusedelems != NULL );
670 assert( permvars != NULL );
671 assert( success != NULL );
672 assert( infeasible != NULL );
673
674 *success = FALSE;
675 *infeasible = FALSE;
676
677 /* if we try to extend the first orbitope generator by another one,
678 * we can change the order of entries in each row of suborbitope */
679 if ( nfilledcols == 2 )
680 {
681 /* check whether each cycle of perm intersects with a row of suborbitope */
682 for (row = 0; row < nrows; ++row)
683 {
684 idx1 = suborbitope[row][0];
685 idx2 = suborbitope[row][1];
686
687 /* if idx1 or idx2 is affected by perm, we can extend the row of the orbitope */
688 if ( idx1 != perm[idx1] )
689 {
690 /* change order of idx1 and idx2 for right extensions */
691 if ( ! leftextension )
692 {
693 suborbitope[row][0] = idx2;
694 suborbitope[row][1] = idx1;
695 }
696 assert( rowisbinary == NULL || rowisbinary[row] == SCIPvarIsBinary(permvars[perm[idx1]]) );
697
698 suborbitope[row][2] = perm[idx1];
699 ++nintersections;
700
701 ++(*nusedelems)[idx1];
702 ++(*nusedelems)[perm[idx1]];
703
704 /* if an element appears too often in the orbitope matrix */
705 if ( (*nusedelems)[idx1] + (*nusedelems)[perm[idx1]] > 3 )
706 {
707 *infeasible = TRUE;
708 break;
709 }
710 }
711 else if ( idx2 != perm[idx2] )
712 {
713 /* change order of idx1 and idx2 for left extensions */
714 if ( leftextension )
715 {
716 suborbitope[row][0] = idx2;
717 suborbitope[row][1] = idx1;
718 }
719 assert( rowisbinary == NULL || rowisbinary[row] == SCIPvarIsBinary(permvars[perm[idx1]]) );
720
721 suborbitope[row][2] = perm[idx2];
722 ++nintersections;
723
724 ++(*nusedelems)[idx2];
725 ++(*nusedelems)[perm[idx2]];
726
727 /* if an element appears too often in the orbitope matrix */
728 if ( (*nusedelems)[idx2] + (*nusedelems)[perm[idx2]] > 3 )
729 {
730 *infeasible = TRUE;
731 break;
732 }
733 }
734 }
735 }
736 else
737 {
738 /* check whether each cycle of perm intersects with a row of suborbitope */
739 for (row = 0; row < nrows; ++row)
740 {
741 idx1 = suborbitope[row][coltoextend];
742
743 /* if idx1 is affected by perm, we can extend the row of the orbitope */
744 if ( idx1 != perm[idx1] )
745 {
746 assert( rowisbinary == NULL || rowisbinary[row] == SCIPvarIsBinary(permvars[perm[idx1]]) );
747
748 suborbitope[row][nfilledcols] = perm[idx1];
749 ++nintersections;
750
751 ++(*nusedelems)[idx1];
752 ++(*nusedelems)[perm[idx1]];
753
754 /* if an element appears to often in the orbitope matrix */
755 if ( (*nusedelems)[idx1] + (*nusedelems)[perm[idx1]] > 3 )
756 {
757 *infeasible = TRUE;
758 break;
759 }
760 }
761 }
762 }
763
764 /* if there are too few intersection, this is not an orbitope */
765 if ( nintersections > 0 && nintersections < nrows )
766 *infeasible = TRUE;
767 else if ( nintersections == nrows )
768 *success = TRUE;
769
770 return SCIP_OKAY;
771 }
772
773
774 /** compute components of symmetry group */
775 SCIP_RETCODE SCIPcomputeComponentsSym(
776 SCIP* scip, /**< SCIP instance */
777 SYM_SYMTYPE symtype, /**< type of symmetries in perms */
778 int** perms, /**< permutation generators as
779 * (either nperms x npermvars or npermvars x nperms) matrix */
780 int nperms, /**< number of permutations */
781 SCIP_VAR** permvars, /**< variables on which permutations act */
782 int npermvars, /**< number of variables for permutations */
783 SCIP_Bool transposed, /**< transposed permutation generators as (npermvars x nperms) matrix */
784 int** components, /**< array containing the indices of permutations sorted by components */
785 int** componentbegins, /**< array containing in i-th position the first position of
786 * component i in components array */
787 int** vartocomponent, /**< array containing for each permvar the index of the component it is
788 * contained in (-1 if not affected) */
789 unsigned** componentblocked, /**< array to store which symmetry methods have been used on a component
790 * using the same bitset information as for misc/usesymmetry */
791 int* ncomponents /**< pointer to store number of components of symmetry group */
792 )
793 {
794 SCIP_DISJOINTSET* componentstovar = NULL;
795 int* permtovarcomp;
796 int* permtocomponent;
797 int p;
798 int i;
799 int idx;
800
801 assert( scip != NULL );
802 assert( permvars != NULL );
803 assert( npermvars > 0 );
804 assert( perms != NULL );
805 assert( components != NULL );
806 assert( componentbegins != NULL );
807 assert( vartocomponent != NULL );
808 assert( componentblocked != NULL );
809 assert( ncomponents != NULL );
810
811 if ( nperms <= 0 )
812 return SCIP_OKAY;
813
814 SCIP_CALL( SCIPdisjointsetCreate(&componentstovar, SCIPblkmem(scip), npermvars) );
815 *ncomponents = npermvars;
816
817 /* init array that stores for each permutation the representative of its affected variables */
818 SCIP_CALL( SCIPallocBufferArray(scip, &permtovarcomp, nperms) );
819 for (p = 0; p < nperms; ++p)
820 permtovarcomp[p] = -1;
821
822 /* find permutation components and store for each variable an affecting permutation (or -1) */
823 SCIP_CALL( SCIPallocBlockMemoryArray(scip, vartocomponent, npermvars) );
824 for (i = 0; i < npermvars; ++i)
825 {
826 (*vartocomponent)[i] = -1;
827
828 for (p = 0; p < nperms; ++p)
829 {
830 int img;
831
832 img = transposed ? perms[i][p] : perms[p][i];
833
834 /* perm p affects i -> possibly merge var components */
835 if ( img != i )
836 {
837 int component1;
838 int component2;
839 int representative;
840
841 if ( img >= npermvars )
842 {
843 assert( symtype == SYM_SYMTYPE_SIGNPERM );
844 img -= npermvars;
845 assert( 0 <= img && img < npermvars );
846 }
847
848 component1 = SCIPdisjointsetFind(componentstovar, i);
849 component2 = SCIPdisjointsetFind(componentstovar, img);
850 (*vartocomponent)[i] = p;
851 (*vartocomponent)[img] = p;
852
853 /* ensure component1 <= component2 */
854 if ( component2 < component1 )
855 {
856 int swap;
857
858 swap = component1;
859 component1 = component2;
860 component2 = swap;
861 }
862
863 /* init permtovarcomp[p] to component of first moved variable or update the value */
864 if ( permtovarcomp[p] == -1 )
865 {
866 permtovarcomp[p] = component1;
867 representative = component1;
868 }
869 else
870 {
871 permtovarcomp[p] = SCIPdisjointsetFind(componentstovar, permtovarcomp[p]);
872 representative = permtovarcomp[p];
873 }
874
875 /* merge both components if they differ */
876 if ( component1 != component2 )
877 {
878 SCIPdisjointsetUnion(componentstovar, component1, component2, TRUE);
879 --(*ncomponents);
880 }
881
882 /* possibly merge new component and permvartocom[p] and ensure the latter
883 * to have the smallest value */
884 if ( representative != component1 && representative != component2 )
885 {
886 if ( representative > component1 )
887 {
888 SCIPdisjointsetUnion(componentstovar, component1, representative, TRUE);
889 permtovarcomp[p] = component1;
890 }
891 else
892 SCIPdisjointsetUnion(componentstovar, representative, component1, TRUE);
893 --(*ncomponents);
894 }
895 else if ( representative > component1 )
896 {
897 assert( representative == component2 );
898 permtovarcomp[p] = component1;
899 }
900 }
901 }
902
903 /* reduce number of components by singletons */
904 if ( (*vartocomponent)[i] == -1 )
905 --(*ncomponents);
906 }
907 assert( *ncomponents > 0 );
908
909 /* update permvartocomp array to final variable representatives */
910 for (p = 0; p < nperms; ++p)
911 permtovarcomp[p] = SCIPdisjointsetFind(componentstovar, permtovarcomp[p]);
912
913 /* init components array by trivial natural order of permutations */
914 SCIP_CALL( SCIPallocBlockMemoryArray(scip, components, nperms) );
915 for (p = 0; p < nperms; ++p)
916 (*components)[p] = p;
917
918 /* get correct order of components array */
919 SCIPsortIntInt(permtovarcomp, *components, nperms);
920
921 /* determine componentbegins and store components for each permutation */
922 SCIP_CALL( SCIPallocBlockMemoryArray(scip, componentbegins, *ncomponents + 1) );
923 SCIP_CALL( SCIPallocBufferArray(scip, &permtocomponent, nperms) );
924
925 (*componentbegins)[0] = 0;
926 permtocomponent[(*components)[0]] = 0;
927 idx = 0;
928
929 for (p = 1; p < nperms; ++p)
930 {
931 if ( permtovarcomp[p] > permtovarcomp[p - 1] )
932 (*componentbegins)[++idx] = p;
933
934 assert( (*components)[p] >= 0 );
935 assert( (*components)[p] < nperms );
936 permtocomponent[(*components)[p]] = idx;
937 }
938 assert( *ncomponents == idx + 1 );
939 (*componentbegins)[++idx] = nperms;
940
941 /* determine vartocomponent */
942 for (i = 0; i < npermvars; ++i)
943 {
944 int permidx;
945 permidx = (*vartocomponent)[i];
946 assert( -1 <= permidx && permidx < nperms );
947
948 if ( permidx != -1 )
949 {
950 assert( 0 <= permtocomponent[permidx] );
951 assert( permtocomponent[permidx] < *ncomponents );
952
953 (*vartocomponent)[i] = permtocomponent[permidx];
954 }
955 }
956
957 /* init componentblocked */
958 SCIP_CALL( SCIPallocBlockMemoryArray(scip, componentblocked, *ncomponents) );
959 for (i = 0; i < *ncomponents; ++i)
960 (*componentblocked)[i] = 0;
961
962 SCIPfreeBufferArray(scip, &permtocomponent);
963 SCIPfreeBufferArray(scip, &permtovarcomp);
964 SCIPdisjointsetFree(&componentstovar, SCIPblkmem(scip));
965
966 #ifdef SCIP_OUTPUT
967 printf("number of components: %d\n", *ncomponents);
968 for (i = 0; i < *ncomponents; ++i)
969 {
970 printf("Component %d contains the following permutations:\n\t", i);
971 for (p = (*componentbegins)[i]; p < (*componentbegins)[i + 1]; ++p)
972 {
973 printf("%d, ", (*components)[p]);
974 }
975 printf("\n");
976 }
977 #endif
978
979 return SCIP_OKAY;
980 }
981
982
983 /** generate variable matrix for orbitope constraint handler
984 *
985 * @pre if storelexorder is TRUE, then the permutations define an orbitope
986 */
987 SCIP_RETCODE SCIPgenerateOrbitopeVarsMatrix(
988 SCIP* scip, /**< SCIP instance */
989 SCIP_VAR**** vars, /**< pointer to matrix of orbitope variables */
990 int nrows, /**< number of rows of orbitope */
991 int ncols, /**< number of columns of orbitope */
992 SCIP_VAR** permvars, /**< superset of variables that are contained in orbitope */
993 int npermvars, /**< number of variables in permvars array */
994 int** orbitopevaridx, /**< permuted index table of variables in permvars that are contained in orbitope */
995 int* columnorder, /**< permutation to reorder column of orbitopevaridx */
996 int* nusedelems, /**< array storing how often an element was used in the orbitope */
997 SCIP_Shortbool* rowisbinary, /**< array encoding whether a row contains only binary variables (or NULL) */
998 SCIP_Bool* infeasible, /**< pointer to store whether the potential orbitope is not an orbitope */
999 SCIP_Bool storelexorder, /**< whether the lexicographic order induced by the orbitope shall be stored */
1000 int** lexorder, /**< pointer to array storing the lexorder (or NULL) */
1001 int* nvarsorder, /**< pointer to store number of variables in lexorder (or NULL) */
1002 int* maxnvarsorder /**< pointer to store maximum number of variables in lexorder (or NULL) */
1003 )
1004 {
1005 int nfilledcols = 0;
1006 int curcolumn;
1007 int i;
1008 int cnt;
1009 int nvarsorderold = 0;
1010
1011 assert( vars != NULL );
1012 assert( nrows > 0 );
1013 assert( ncols > 0 );
1014 assert( permvars != NULL );
1015 assert( npermvars > 0 );
1016 assert( orbitopevaridx != NULL );
1017 assert( columnorder != NULL );
1018 assert( nusedelems != NULL );
1019 assert( infeasible != NULL );
1020 assert( ! storelexorder || lexorder != NULL );
1021 assert( ! storelexorder || nvarsorder != NULL );
1022 assert( ! storelexorder || maxnvarsorder != NULL );
1023
1024 /* possibly store lexicographic order defined by orbitope
1025 *
1026 * position (i,j) of orbitope has position nrows * j + i in lexicographic order
1027 */
1028 if ( storelexorder )
1029 {
1030 assert( *nvarsorder == *maxnvarsorder );
1031 assert( lexorder != NULL );
1032
1033 *maxnvarsorder += nrows * ncols;
1034 nvarsorderold = *nvarsorder;
1035
1036 if ( *lexorder == NULL )
1037 {
1038 SCIP_CALL( SCIPallocBlockMemoryArray(scip, lexorder, *maxnvarsorder) );
1039 }
1040 else
1041 {
1042 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, lexorder, *nvarsorder, *maxnvarsorder) );
1043 }
1044 }
1045
1046 curcolumn = ncols - 1;
1047
1048 /* start filling vars matrix with the right-most column w.r.t. columnorder */
1049 while ( curcolumn >= 0 && columnorder[curcolumn] >= 0 && ! *infeasible )
1050 {
1051 cnt = 0;
1052 for (i = 0; i < nrows; ++i)
1053 {
1054 /* skip rows containing non-binary variables */
1055 if ( rowisbinary != NULL && ! rowisbinary[i] )
1056 continue;
1057
1058 assert( 0 <= orbitopevaridx[i][curcolumn] && orbitopevaridx[i][curcolumn] < npermvars );
1059 assert( SCIPvarIsBinary(permvars[orbitopevaridx[i][curcolumn]]) );
1060
1061 /* elements in first column of orbitope have to appear exactly once in the orbitope */
1062 if ( nfilledcols == 0 && nusedelems[orbitopevaridx[i][curcolumn]] > 1 )
1063 {
1064 *infeasible = TRUE;
1065 assert( ! storelexorder );
1066 break;
1067 }
1068
1069 if ( storelexorder )
1070 {
1071 (*lexorder)[nvarsorderold + nrows * nfilledcols + cnt] = orbitopevaridx[i][curcolumn];
1072 ++(*nvarsorder);
1073 }
1074 (*vars)[cnt++][nfilledcols] = permvars[orbitopevaridx[i][curcolumn]];
1075 }
1076 --curcolumn;
1077 ++nfilledcols;
1078 }
1079
1080 /* There are three possibilities for the structure of columnorder:
1081 * 1) [0, 1, -1, -1, ..., -1]
1082 * 2) [0, 1, 1, 1, ..., 1]
1083 * 3) [0, 1, -1, -1, ...., -1, 1, 1, ..., 1]
1084 */
1085 /* Either we are in case 1) or case 3), or all columns should have been added to vars in case 2) */
1086 assert( curcolumn > 1 || (curcolumn < 0 && nfilledcols == ncols) );
1087
1088 if ( curcolumn > 1 && ! *infeasible )
1089 {
1090 /* add column with columnorder 1 to vars */
1091 cnt = 0;
1092 for (i = 0; i < nrows; ++i)
1093 {
1094 /* skip rows containing non-binary variables*/
1095 if ( rowisbinary != NULL && ! rowisbinary[i] )
1096 continue;
1097
1098 assert( orbitopevaridx[i][1] < npermvars );
1099 assert( SCIPvarIsBinary(permvars[orbitopevaridx[i][1]]) );
1100
1101 if ( storelexorder )
1102 {
1103 (*lexorder)[nvarsorderold + nrows * nfilledcols + cnt] = orbitopevaridx[i][1];
1104 ++(*nvarsorder);
1105 }
1106 (*vars)[cnt++][nfilledcols] = permvars[orbitopevaridx[i][1]];
1107 }
1108 ++nfilledcols;
1109
1110 /* add column with columnorder 0 to vars */
1111 cnt = 0;
1112 for (i = 0; i < nrows; ++i)
1113 {
1114 /* skip rows containing non-binary variables*/
1115 if ( rowisbinary != NULL && ! rowisbinary[i] )
1116 continue;
1117
1118 assert( orbitopevaridx[i][0] < npermvars );
1119 assert( SCIPvarIsBinary(permvars[orbitopevaridx[i][0]]) );
1120
1121 if ( storelexorder )
1122 {
1123 (*lexorder)[nvarsorderold + nrows * nfilledcols + cnt] = orbitopevaridx[i][0];
1124 ++(*nvarsorder);
1125 }
1126 (*vars)[cnt++][nfilledcols] = permvars[orbitopevaridx[i][0]];
1127 }
1128 ++nfilledcols;
1129
1130 /* add columns with a negative column order to vars */
1131 if ( nfilledcols < ncols )
1132 {
1133 assert( ncols > 2 );
1134
1135 curcolumn = 2;
1136 while ( nfilledcols < ncols && ! *infeasible )
1137 {
1138 assert( columnorder[curcolumn] < 0 );
1139
1140 cnt = 0;
1141 for (i = 0; i < nrows; ++i)
1142 {
1143 /* skip rows containing non-binary variables*/
1144 if ( rowisbinary != NULL && ! rowisbinary[i] )
1145 continue;
1146
1147 assert( orbitopevaridx[i][curcolumn] < npermvars );
1148 assert( SCIPvarIsBinary(permvars[orbitopevaridx[i][curcolumn]]) );
1149
1150 /* elements in last column of orbitope have to appear exactly once in the orbitope */
1151 if ( nfilledcols == ncols - 1 && nusedelems[orbitopevaridx[i][curcolumn]] > 1 )
1152 {
1153 *infeasible = TRUE;
1154 assert( ! storelexorder );
1155 break;
1156 }
1157
1158 if ( storelexorder )
1159 {
1160 (*lexorder)[nvarsorderold + nrows * nfilledcols + cnt] = orbitopevaridx[i][curcolumn];
1161 ++(*nvarsorder);
1162 }
1163 (*vars)[cnt++][nfilledcols] = permvars[orbitopevaridx[i][curcolumn]];
1164 }
1165 ++curcolumn;
1166 ++nfilledcols;
1167 }
1168 }
1169 }
1170
1171 return SCIP_OKAY;
1172 }
1173
1174
1175 /** checks whether an orbitope is a packing or partitioning orbitope; if npprows != NULL,
1176 * count how many rows are contained in packing/partitioning constraints
1177 */
1178 SCIP_RETCODE SCIPisPackingPartitioningOrbitope(
1179 SCIP* scip, /**< SCIP data structure */
1180 SCIP_VAR*** vars, /**< variable matrix of orbitope constraint */
1181 int nrows, /**< pointer to number of rows of variable matrix */
1182 int ncols, /**< number of columns of variable matrix */
1183 SCIP_Bool** pprows, /**< pointer to store which rows are are contained in
1184 * packing/partitioning constraints or NULL if not needed */
1185 int* npprows, /**< pointer to store how many rows are contained
1186 * in packing/partitioning constraints or NULL if not needed */
1187 SCIP_ORBITOPETYPE* type /**< pointer to store type of orbitope constraint after strengthening */
1188 )
1189 {
1190 SCIP_CONSHDLR* setppcconshdlr;
1191 SCIP_CONS** setppcconss;
1192 int nsetppcconss;
1193 int* covered;
1194 int nprobvars;
1195 int* rowidxvar;
1196 int* rowcoveragesetppc;
1197 int* rowsinsetppc;
1198 int ncovered;
1199 int ncoveredpart;
1200 int i;
1201 int j;
1202 int c;
1203
1204 assert( scip != NULL );
1205 assert( vars != NULL );
1206 assert( vars != NULL );
1207 assert( nrows > 0 );
1208 assert( ncols > 0 );
1209 assert( type != NULL );
1210
1211 *type = SCIP_ORBITOPETYPE_FULL;
1212 if ( npprows != NULL )
1213 *npprows = 0;
1214
1215 setppcconshdlr = SCIPfindConshdlr(scip, "setppc");
1216 if ( setppcconshdlr == NULL )
1217 return SCIP_OKAY;
1218
1219 setppcconss = SCIPconshdlrGetConss(setppcconshdlr);
1220 nsetppcconss = SCIPconshdlrGetNConss(setppcconshdlr);
1221
1222 /* we can terminate early if there are not sufficiently many setppc conss
1223 * (for orbitopes treating a full component, we might allow to remove rows
1224 * not contained in setppc cons; for this reason we need the second check)
1225 */
1226 if ( nsetppcconss == 0 || (nsetppcconss < nrows && npprows == NULL ))
1227 return SCIP_OKAY;
1228 assert( setppcconss != NULL );
1229
1230 /* whether a row is contained in packing/partitioning constraint */
1231 SCIP_CALL( SCIPallocClearBufferArray(scip, &covered, nrows) );
1232 ncovered = 0;
1233 ncoveredpart = 0;
1234
1235 nprobvars = SCIPgetNTotalVars(scip);
1236
1237 /* array storing index of orbitope row a variable is contained in */
1238 SCIP_CALL( SCIPallocBufferArray(scip, &rowidxvar, nprobvars) );
1239
1240 for (i = 0; i < nprobvars; ++i)
1241 rowidxvar[i] = -1;
1242
1243 for (i = 0; i < nrows; ++i)
1244 {
1245 for (j = 0; j < ncols; ++j)
1246 {
1247 assert( 0 <= SCIPvarGetIndex(vars[i][j]) && SCIPvarGetIndex(vars[i][j]) < nprobvars );
1248 rowidxvar[SCIPvarGetIndex(vars[i][j])] = i;
1249 }
1250 }
1251
1252 /* storage for number of vars per row that are contained in current setppc cons and
1253 * labels of rows intersecting with current setppc cons
1254 */
1255 SCIP_CALL( SCIPallocCleanBufferArray(scip, &rowcoveragesetppc, nrows) );
1256 SCIP_CALL( SCIPallocClearBufferArray(scip, &rowsinsetppc, nrows) );
1257
1258 /* iterate over set packing and partitioning constraints and check whether the constraint's
1259 * support is a row r of the orbitope (covered[r] = 2) or contains row r (covered[r] = 1)
1260 *
1261 * @todo Check whether we can improve the following loop by using a hash value to check
1262 * whether the setppccons intersects the orbitope matrix
1263 */
1264 for (c = 0; c < nsetppcconss && ncoveredpart < ncols; ++c)
1265 {
1266 int nsetppcvars;
1267 SCIP_VAR** setppcvars;
1268 int nrowintersect = 0;
1269 int nvarsinorbitope;
1270
1271 /* skip covering constraints */
1272 if ( SCIPgetTypeSetppc(scip, setppcconss[c]) == SCIP_SETPPCTYPE_COVERING )
1273 continue;
1274
1275 /* get number of set packing/partitioning variables */
1276 nsetppcvars = SCIPgetNVarsSetppc(scip, setppcconss[c]);
1277
1278 /* constraint does not contain enough variables */
1279 if ( nsetppcvars < ncols )
1280 continue;
1281
1282 setppcvars = SCIPgetVarsSetppc(scip, setppcconss[c]);
1283 assert( setppcvars != NULL );
1284
1285 /* upper bound on variables potentially contained in orbitope */
1286 nvarsinorbitope = nsetppcvars;
1287
1288 /* for each setppc var, check whether it appears in a row of the orbitope and store
1289 * for each row the number of such variables; can be terminated early, if less than
1290 * ncols variables are contained in the orbitope
1291 */
1292 for (i = 0; i < nsetppcvars && nvarsinorbitope >= ncols; ++i)
1293 {
1294 SCIP_VAR* var;
1295 int idx;
1296 int rowidx;
1297
1298 var = setppcvars[i];
1299 idx = SCIPvarGetIndex(var);
1300
1301 assert( 0 <= idx && idx < nprobvars );
1302
1303 rowidx = rowidxvar[idx];
1304
1305 /* skip variables not contained in the orbitope */
1306 if ( rowidx < 0 )
1307 {
1308 --nvarsinorbitope;
1309 continue;
1310 }
1311
1312 /* skip variables corresponding to already treated rows */
1313 if ( covered[rowidx] == 2 || (covered[rowidx] == 1 && (nsetppcvars > ncols || nrowintersect > 1)) )
1314 {
1315 --nvarsinorbitope;
1316 continue;
1317 }
1318
1319 /* store information which rows intersect the setppc cons's support */
1320 if ( rowcoveragesetppc[rowidx] == 0 )
1321 rowsinsetppc[nrowintersect++] = rowidx;
1322 ++(rowcoveragesetppc[rowidx]);
1323
1324 /* we can stop early if not enough variables are left to completely cover one of the rows that
1325 * intersect the setppc cons
1326 */
1327 if ( nsetppcvars - nrowintersect < ncols - 1 )
1328 break;
1329 }
1330
1331 /* store whether rows coincide with set partitioning cons's support or whether
1332 * row is covered by a set packing/partitioning cons's support
1333 */
1334 if ( SCIPgetTypeSetppc(scip, setppcconss[c]) == SCIP_SETPPCTYPE_PARTITIONING
1335 && nrowintersect == 1 && rowcoveragesetppc[rowsinsetppc[0]] == ncols && nsetppcvars == ncols )
1336 {
1337 if ( covered[rowsinsetppc[0]] == 1 )
1338 --ncovered;
1339 covered[rowsinsetppc[0]] = 2;
1340 ++ncoveredpart;
1341 ++ncovered;
1342 }
1343 else
1344 {
1345 for (i = 0; i < nrowintersect; ++i)
1346 {
1347 if ( covered[rowsinsetppc[i]] == 0 && rowcoveragesetppc[rowsinsetppc[i]] >= ncols )
1348 {
1349 covered[rowsinsetppc[i]] = 1;
1350 ++ncovered;
1351 }
1352 }
1353 }
1354
1355 /* reset data */
1356 for (i = 0; i < nrowintersect; ++i)
1357 rowcoveragesetppc[rowsinsetppc[i]] = 0;
1358 }
1359
1360 /* check type of orbitope */
1361 if ( ncovered == nrows )
1362 {
1363 if ( ncoveredpart == nrows )
1364 *type = SCIP_ORBITOPETYPE_PARTITIONING;
1365 else
1366 *type = SCIP_ORBITOPETYPE_PACKING;
1367 }
1368
1369 if ( npprows != NULL )
1370 *npprows = ncovered;
1371
1372 if ( pprows != NULL )
1373 {
1374 SCIP_CALL( SCIPallocBlockMemoryArray(scip, pprows, nrows) );
1375 for (i = 0; i < nrows; ++i)
1376 (*pprows)[i] = covered[i] > 0 ? 1 : 0;
1377 }
1378
1379 SCIPfreeBufferArray(scip, &rowsinsetppc);
1380 SCIPfreeCleanBufferArray(scip, &rowcoveragesetppc);
1381 SCIPfreeBufferArray(scip, &rowidxvar);
1382 SCIPfreeBufferArray(scip, &covered);
1383
1384 return SCIP_OKAY;
1385 }
1386
1387 /** checks whether a (signed) permutation is an involution */
1388 static
1389 SCIP_RETCODE isPermInvolution(
1390 int* perm, /**< permutation */
1391 int permlen, /**< number of original (non-negated) variables in a permutation */
1392 SCIP_Bool* isinvolution, /**< pointer to store whether perm is an involution */
1393 int* ntwocycles /**< pointer to store number of 2-cycles in an involution */
1394 )
1395 {
1396 int v;
1397
1398 assert( perm != NULL );
1399 assert( permlen > 0 );
1400 assert( isinvolution != NULL );
1401 assert( ntwocycles != NULL );
1402
1403 *ntwocycles = 0;
1404 *isinvolution = TRUE;
1405 for (v = 0; v < permlen && *isinvolution; ++v)
1406 {
1407 /* do not handle variables twice */
1408 if ( perm[v] <= v )
1409 continue;
1410
1411 /* detect two cycles */
1412 if ( perm[perm[v]] == v )
1413 ++(*ntwocycles);
1414 else
1415 *isinvolution = FALSE;
1416 }
1417
1418 return SCIP_OKAY;
1419 }
1420
1421 /** checks whether selected permutations define orbitopal symmetries */
1422 static
1423 SCIP_RETCODE detectOrbitopalSymmetries(
1424 SCIP* scip, /**< SCIP pointer */
1425 int** perms, /**< array of permutations */
1426 int* selectedperms, /**< indices of permutations in perm that shall be considered */
1427 int nselectedperms, /**< number of permutations in selectedperms */
1428 int permlen, /**< number of variables in a permutation */
1429 int nrows, /**< number of rows of potential orbitopal symmetries */
1430 SCIP_Bool* success, /**< pointer to store if orbitopal symmetries could be found */
1431 int**** matrices, /**< pointer to store matrices of orbitopal symmetries */
1432 int** ncols, /**< pointer to store number of columns of matrices in matrices */
1433 int* nmatrices /**< pointer to store the number of matrices in matrices */
1434 )
1435 { /*lint --e{771}*/
1436 SCIP_DISJOINTSET* conncomps;
1437 SCIP_DISJOINTSET* compcolors;
1438 int* complastperm;
1439 int* permstoconsider;
1440 int* colorbegins;
1441 int* compidx;
1442 int* colidx;
1443 int* varidx;
1444 int* degrees;
1445 int* perm;
1446 int nposdegree = 0;
1447 int npermstoconsider;
1448 int colorrepresentative1 = -1;
1449 int colorrepresentative2 = -1;
1450 int elemtomove;
1451 int ncurcols;
1452 int curcomp1;
1453 int curcomp2;
1454 int curdeg1;
1455 int curdeg2;
1456 int curcolor1;
1457 int curcolor2;
1458 int ncolors;
1459 int cnt;
1460 int c;
1461 int p;
1462 int v;
1463 int w;
1464
1465 assert( scip != NULL );
1466 assert( perms != NULL );
1467 assert( selectedperms != NULL );
1468 assert( nselectedperms >= 0 );
1469 assert( permlen > 0 );
1470 assert( nrows > 0 || nselectedperms == 0 );
1471 assert( success != NULL );
1472 assert( matrices != NULL );
1473 assert( ncols != NULL );
1474 assert( nmatrices != NULL );
1475
1476 /* initialize data */
1477 *success = TRUE;
1478 *matrices = NULL;
1479 *ncols = NULL;
1480 *nmatrices = 0;
1481
1482 /* we have found the empty set of orbitopal symmetries */
1483 if ( nselectedperms == 0 )
1484 return SCIP_OKAY;
1485
1486 /* build a graph to encode potential orbitopal symmetries
1487 *
1488 * The 2-cycles of a permutation define a set of edges that need to be added simultaneously. We encode
1489 * this as a disjoint set data structure to only encode the connected components of the graph. To ensure
1490 * correctness, we keep track of the degrees of each node and extend a component by a permutation only if
1491 * all nodes to be extended have the same degree. We also make sure that each connected component is a
1492 * path. Although this might not detect all orbitopal symmetries, it seems to cover most of the cases when
1493 * computing symmetries using bliss. We also keep track of which variables are affected by the same
1494 * permutations by coloring the connected components.
1495 */
1496 SCIP_CALL( SCIPcreateDisjointset(scip, &conncomps, permlen) );
1497 SCIP_CALL( SCIPcreateDisjointset(scip, &compcolors, permlen) );
1498 SCIP_CALL( SCIPallocClearBufferArray(scip, °rees, permlen) );
1499 SCIP_CALL( SCIPallocBufferArray(scip, &complastperm, permlen) );
1500 for (v = 0; v < permlen; ++v)
1501 complastperm[v] = -1;
1502
1503 /* try to add edges of permutations to graph */
1504 for (p = 0; p < nselectedperms; ++p)
1505 {
1506 perm = perms[selectedperms[p]];
1507 curdeg1 = -1;
1508 curdeg2 = -1;
1509
1510 /* check whether all potential edges can be added */
1511 for (v = 0; v < permlen; ++v)
1512 {
1513 /* treat each cycle exactly once */
1514 if ( perm[v] <= v )
1515 continue;
1516 w = perm[v];
1517
(1) Event value_overwrite: |
Overwriting previous write to "curcomp1" with value from "SCIPdisjointsetFind(conncomps, v)". |
Also see events: |
[returned_value] |
1518 curcomp1 = SCIPdisjointsetFind(conncomps, v);
1519 curcomp2 = SCIPdisjointsetFind(conncomps, w);
1520
1521 /* an edge is not allowed to connect nodes from the same connected component */
1522 if ( curcomp1 == curcomp2 )
1523 break;
1524
1525 /* permutation p is not allowed to add two edges to the same connected component */
1526 if ( complastperm[curcomp1] == p || complastperm[curcomp2] == p )
1527 break;
1528
1529 /* get colors of nodes */
1530 curcolor1 = SCIPdisjointsetFind(compcolors, v);
1531 curcolor2 = SCIPdisjointsetFind(compcolors, w);
1532
1533 /* an edge is not allowed to connect two nodes with the same color */
1534 if ( curcolor1 == curcolor2 )
1535 break;
1536
1537 if ( curdeg1 == -1 )
1538 {
1539 assert( curdeg2 == -1 );
1540
1541 curdeg1 = degrees[v];
1542 curdeg2 = degrees[w];
1543 colorrepresentative1 = curcolor1;
1544 colorrepresentative2 = curcolor2;
1545
1546 /* stop, we will generate a vertex with degree 3 */
1547 if ( curdeg1 == 2 || curdeg2 == 2 )
1548 break;
1549 }
1550 else
1551 {
1552 /* check whether nodes have compatible degrees */
1553 if ( ! ((curdeg1 == degrees[v] && curdeg2 == degrees[w])
1554 || (curdeg1 == degrees[w] && curdeg2 == degrees[v])) )
1555 break;
1556 assert( colorrepresentative1 >= 0 );
1557 assert( colorrepresentative2 >= 0 || curdeg2 == -1 );
1558
1559 /* check whether all components have compatible colors */
1560 if ( curdeg1 > 0 && curcolor1 != colorrepresentative1 && curcolor2 != colorrepresentative1 )
1561 break;
1562 if ( curdeg2 > 0 && curcolor1 != colorrepresentative2 && curcolor2 != colorrepresentative2 )
1563 break;
1564 }
1565
1566 /* store that permutation p extends the connected components */
1567 complastperm[curcomp1] = p;
1568 complastperm[curcomp2] = p;
1569 }
1570
1571 /* terminate if not all edges can be added */
1572 if ( v < permlen )
1573 {
1574 *success = FALSE;
1575 goto FREEMEMORY;
1576 }
1577 assert( curdeg1 >= 0 && curdeg2 >= 0 );
1578
1579 /* add edges to graph */
1580 for (v = 0; v < permlen; ++v)
1581 {
1582 /* treat each cycle exactly once */
1583 if ( perm[v] <= v )
1584 continue;
1585 w = perm[v];
1586
(2) Event returned_value: |
Assigning value from "SCIPdisjointsetFind(conncomps, v)" to "curcomp1" here, but that stored value is overwritten before it can be used. |
Also see events: |
[value_overwrite] |
1587 curcomp1 = SCIPdisjointsetFind(conncomps, v);
1588 curcomp2 = SCIPdisjointsetFind(conncomps, w);
1589 assert( curcomp1 != curcomp2 );
1590
1591 /* add edge */
1592 SCIPdisjointsetUnion(conncomps, v, w, FALSE);
1593 ++degrees[v];
1594 ++degrees[w];
1595
1596 /* possibly update colors */
1597 curcolor1 = SCIPdisjointsetFind(compcolors, v);
1598 curcolor2 = SCIPdisjointsetFind(compcolors, w);
1599
1600 if ( curcolor1 != curcolor2 )
1601 {
1602 SCIPdisjointsetUnion(compcolors, colorrepresentative1, v, TRUE);
1603 SCIPdisjointsetUnion(compcolors, colorrepresentative1, w, TRUE);
1604 }
1605 }
1606 }
1607
1608 /* find non-trivial components */
1609 for (v = 0; v < permlen; ++v)
1610 {
1611 if ( degrees[v] > 0 )
1612 ++nposdegree;
1613 }
1614
1615 SCIP_CALL( SCIPallocBufferArray(scip, &compidx, nposdegree) );
1616 SCIP_CALL( SCIPallocBufferArray(scip, &colidx, nposdegree) );
1617 SCIP_CALL( SCIPallocBufferArray(scip, &varidx, nposdegree) );
1618
1619 for (v = 0, w = 0; v < permlen; ++v)
1620 {
1621 if ( degrees[v] > 0 )
1622 {
1623 compidx[w] = SCIPdisjointsetFind(conncomps, v);
1624 colidx[w] = SCIPdisjointsetFind(compcolors, v);
1625 #ifndef NDEBUG
1626 if ( w > 0 && compidx[w] == compidx[w-1] )
1627 assert( colidx[w] == colidx[w-1]);
1628 #endif
1629 varidx[w++] = v;
1630 }
1631 }
1632 assert( w == nposdegree );
1633
1634 /* sort variable indices: first by colors, then by components */
1635 SCIP_CALL( SCIPallocBufferArray(scip, &colorbegins, permlen + 1) );
1636
1637 SCIPsortIntIntInt(colidx, compidx, varidx, nposdegree);
1638 w = 0;
1639 ncolors = 0;
1640 colorbegins[0] = 0;
1641 for (v = 1; v < nposdegree; ++v)
1642 {
1643 if ( colidx[v] != colidx[w] )
1644 {
1645 SCIPsortIntInt(&compidx[w], &varidx[w], v - w);
1646 colorbegins[++ncolors] = v;
1647 w = v;
1648 }
1649 }
1650 SCIPsortIntInt(&compidx[w], &varidx[w], nposdegree - w);
1651 colorbegins[++ncolors] = nposdegree;
1652
1653 /* try to find the correct order of variable indices per color class */
1654 SCIP_CALL( SCIPallocBufferArray(scip, &permstoconsider, nselectedperms) );
1655
1656 SCIP_CALL( SCIPallocBlockMemoryArray(scip, matrices, ncolors) );
1657 SCIP_CALL( SCIPallocBlockMemoryArray(scip, ncols, ncolors) );
1658 *nmatrices = ncolors;
1659
1660 for (c = 0; c < ncolors; ++c)
1661 {
1662 /* find an element in the first connected component with degree 1 */
1663 for (v = colorbegins[c]; compidx[v] == compidx[colorbegins[c]]; ++v)
1664 {
1665 assert( v < nposdegree ); /* there should be a node of degree 1 */
1666
1667 if ( degrees[varidx[v]] == 1 )
1668 break;
1669 }
1670 assert( compidx[v] == compidx[colorbegins[c]] );
1671 elemtomove = varidx[v];
1672
1673 /* find the permutations affecting the variables in the first connected component */
1674 npermstoconsider = 0;
1675 for (p = 0; p < nselectedperms; ++p)
1676 {
1677 perm = perms[selectedperms[p]];
1678 for (v = colorbegins[c]; compidx[v] == compidx[colorbegins[c]] && v < nposdegree; ++v)
1679 {
1680 if ( perm[varidx[v]] != varidx[v] )
1681 {
1682 permstoconsider[npermstoconsider++] = selectedperms[p];
1683 break;
1684 }
1685 }
1686 }
1687
1688 /* allocate memory for matrix */
1689 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*matrices)[c], nrows) );
1690 (*ncols)[c] = npermstoconsider + 1;
1691 for (p = 0; p < nrows; ++p)
1692 {
1693 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*matrices)[c][p], (*ncols)[c]) );
1694 }
1695
1696 /* find a permutation that moves the degree-1 element and iteratively extend this to a matrix */
1697 assert( degrees[elemtomove] == 1 );
1698
1699 /* find the first and second column */
1700 for (p = 0; p < npermstoconsider; ++p)
1701 {
1702 perm = perms[permstoconsider[p]];
1703 if ( perm[elemtomove] != elemtomove )
1704 break;
1705 }
1706 assert( p < npermstoconsider );
1707
1708 /* elements moved by perm that have degree 1 are in the first column */
1709 for (v = 0, cnt = 0; v < permlen; ++v)
1710 {
1711 if ( perm[v] > v ) /*lint !e771*/
1712 {
1713 if ( degrees[v] == 1 )
1714 {
1715 (*matrices)[c][cnt][0] = v;
1716 (*matrices)[c][cnt++][1] = perm[v];
1717 }
1718 else
1719 {
1720 (*matrices)[c][cnt][0] = perm[v];
1721 (*matrices)[c][cnt++][1] = v;
1722 }
1723 }
1724 }
1725
1726 /* if the selected permutation do not form orbitopal symmetries */
1727 if ( cnt < nrows )
1728 {
1729 *success = FALSE;
1730 for (p = nrows - 1; p >= 0; --p)
1731 {
1732 SCIPfreeBufferArray(scip, &(*matrices)[c][p]);
1733 }
1734 SCIPfreeBlockMemoryArray(scip, &(*matrices)[c], nrows);
1735 SCIPfreeBlockMemoryArray(scip, matrices, ncolors);
1736 SCIPfreeBlockMemoryArray(scip, ncols, ncolors);
1737 *matrices = NULL;
1738 *ncols = NULL;
1739 goto FREEMOREMEMORY;
1740 }
1741 assert( cnt == nrows );
1742
1743 /* remove p from the list of permutations to be considered */
1744 permstoconsider[p] = permstoconsider[--npermstoconsider];
1745
1746 ncurcols = 1;
1747 while ( npermstoconsider > 0 )
1748 {
1749 elemtomove = (*matrices)[c][0][ncurcols];
1750
1751 /* find permutation moving the elemtomove */
1752 for (p = 0; p < npermstoconsider; ++p)
1753 {
1754 perm = perms[permstoconsider[p]];
1755 if ( perm[elemtomove] != elemtomove )
1756 break;
1757 }
1758 assert( p < npermstoconsider );
1759
1760 /* extend matrix */
1761 for (v = 0; v < nrows; ++v)
1762 {
1763 assert( perm[(*matrices)[c][v][ncurcols]] != (*matrices)[c][v][ncurcols] );
1764 (*matrices)[c][v][ncurcols + 1] = perm[(*matrices)[c][v][ncurcols]];
1765 }
1766 ++ncurcols;
1767 permstoconsider[p] = permstoconsider[--npermstoconsider];
1768 }
1769 }
1770
1771 FREEMOREMEMORY:
1772 SCIPfreeBufferArray(scip, &permstoconsider);
1773 SCIPfreeBufferArray(scip, &colorbegins);
1774 SCIPfreeBufferArray(scip, &varidx);
1775 SCIPfreeBufferArray(scip, &colidx);
1776 SCIPfreeBufferArray(scip, &compidx);
1777
1778 FREEMEMORY:
1779 SCIPfreeBufferArray(scip, &complastperm);
1780 SCIPfreeBufferArray(scip, °rees);
1781 SCIPfreeDisjointset(scip, &compcolors);
1782 SCIPfreeDisjointset(scip, &conncomps);
1783
1784 return SCIP_OKAY;
1785 }
1786
1787 /** checks whether two families of orbitopal symmetries define a double lex matrix, and in case of success, generates matrix
1788 *
1789 * The columns of matrix1 will serve as the columns of the matrix to be generated, the columns of matrix2 will
1790 * serve as rows.
1791 */
1792 static
1793 SCIP_RETCODE isDoublelLexSym(
1794 SCIP* scip, /**< SCIP pointer */
1795 int*** matrices1, /**< first list of matrices associated with orbitopal symmetries */
1796 int nrows1, /**< number of rows of first family of matrices */
1797 int* ncols1, /**< for each matrix in the first family, its number of columns */
1798 int nmatrices1, /**< number of matrices in the first family */
1799 int*** matrices2, /**< second list of matrices associated with orbitopal symmetries */
1800 int nrows2, /**< number of rows of second family of matrices */
1801 int* ncols2, /**< for each matrix in the second family, its number of columns */
1802 int nmatrices2, /**< number of matrices in the second family */
1803 int*** doublelexmatrix, /**< pointer to store combined matrix */
1804 int* nrows, /**< pointer to store number of rows in combined matrix */
1805 int* ncols, /**< pointer to store number of columns in combined matrix */
1806 int** rowsbegin, /**< pointer to store the begin positions of a new lex subset of rows */
1807 int** colsbegin, /**< pointer to store the begin positions of a new lex subset of columns */
1808 SCIP_Bool* success /**< pointer to store whether combined matrix could be generated */
1809 )
1810 {
1811 int* idxtomatrix1;
1812 int* idxtomatrix2;
1813 int* idxtorow1;
1814 int* idxtorow2;
1815 int* idxtocol1;
1816 int* idxtocol2;
1817 int* sortvals;
1818 int elem;
1819 int mat;
1820 int col;
1821 int col2;
1822 int mat2;
1823 int nidx;
1824 int cnt;
1825 int c;
1826 int d;
1827 int i;
1828 int j;
1829
1830 assert( scip != NULL );
1831 assert( matrices1 != NULL );
1832 assert( nrows1 > 0 );
1833 assert( ncols1 != NULL );
1834 assert( nmatrices1 > 0 );
1835 assert( matrices2 != NULL );
1836 assert( nrows2 > 0 || nmatrices2 == 0 );
1837 assert( ncols2 != NULL );
1838 assert( nmatrices2 >= 0 );
1839 assert( doublelexmatrix != NULL );
1840 assert( nrows != NULL );
1841 assert( ncols != NULL );
1842 assert( rowsbegin != NULL );
1843 assert( colsbegin != NULL );
1844 assert( success != NULL );
1845
1846 /* initialize data */
1847 *nrows = nrows1;
1848 *ncols = nrows2;
1849 *success = TRUE;
1850
1851 /* check whether expecteded sizes of matrix match */
1852 for (j = 0, cnt = 0; j < nmatrices1; ++j)
1853 cnt += ncols1[j];
1854 if ( cnt != *ncols )
1855 {
1856 *success = FALSE;
1857 return SCIP_OKAY;
1858 }
1859
1860 for (i = 0, cnt = 0; i < nmatrices2; ++i)
1861 cnt += ncols2[i];
1862 if ( cnt != *nrows )
1863 {
1864 *success = FALSE;
1865 return SCIP_OKAY;
1866 }
1867
1868 /* collect information about entries in matrices */
1869 nidx = nrows1 * nrows2;
1870 SCIP_CALL( SCIPallocBufferArray(scip, &idxtomatrix1, nidx) );
1871 SCIP_CALL( SCIPallocBufferArray(scip, &idxtomatrix2, nidx) );
1872 SCIP_CALL( SCIPallocBufferArray(scip, &idxtorow1, nidx) );
1873 SCIP_CALL( SCIPallocBufferArray(scip, &idxtorow2, nidx) );
1874 SCIP_CALL( SCIPallocBufferArray(scip, &idxtocol1, nidx) );
1875 SCIP_CALL( SCIPallocBufferArray(scip, &idxtocol2, nidx) );
1876
1877 for (c = 0; c < nmatrices1; ++c)
1878 {
1879 for (i = 0; i < nrows1; ++i)
1880 {
1881 for (j = 0; j < ncols1[c]; ++j)
1882 {
1883 idxtomatrix1[matrices1[c][i][j]] = c;
1884 idxtorow1[matrices1[c][i][j]] = i;
1885 idxtocol1[matrices1[c][i][j]] = j;
1886 }
1887 }
1888 }
1889 for (c = 0; c < nmatrices2; ++c)
1890 {
1891 for (i = 0; i < nrows2; ++i)
1892 {
1893 for (j = 0; j < ncols2[c]; ++j)
1894 {
1895 idxtomatrix2[matrices2[c][i][j]] = c;
1896 idxtorow2[matrices2[c][i][j]] = i;
1897 idxtocol2[matrices2[c][i][j]] = j;
1898 }
1899 }
1900 }
1901
1902 /* Find a big matrix such that the columns of this matrix correspond to the columns of matrices in matrices1
1903 * and the rows of this matrix correspond to the columns of matrices in matrices2. In total, this leads to
1904 * a matrix of block shape
1905 *
1906 * A | B | ... | C
1907 * D | E | ... | F
1908 * . | . | | .
1909 * G | H | ... | I
1910 *
1911 * We start by filling the first column of the big matrix by the first column in matrices1. Sort the column
1912 * according to the matrices in matrices2.
1913 */
1914 SCIP_CALL( SCIPallocBufferArray(scip, &sortvals, MAX(*nrows, *ncols)) );
1915 SCIP_CALL( SCIPallocBlockMemoryArray(scip, doublelexmatrix, *nrows) );
1916 for (i = 0; i < *nrows; ++i)
1917 {
1918 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*doublelexmatrix)[i], *ncols) );
1919 (*doublelexmatrix)[i][0] = matrices1[0][i][0];
1920 sortvals[i] = idxtomatrix2[matrices1[0][i][0]];
1921 }
1922 SCIPsortIntPtr(sortvals, (void*) (*doublelexmatrix), *nrows);
1923
1924 /* fill first row of big matrix */
1925 mat = idxtomatrix2[(*doublelexmatrix)[0][0]];
1926 col = idxtocol2[(*doublelexmatrix)[0][0]];
1927 cnt = 0;
1928 for (j = 0; j < *ncols; ++j)
1929 {
1930 /* skip the entry that is already contained in the first column */
1931 if ( matrices2[mat][j][col] == (*doublelexmatrix)[0][0] )
1932 continue;
1933
1934 sortvals[cnt++] = idxtomatrix1[matrices2[mat][j][col]];
1935 (*doublelexmatrix)[0][cnt] = matrices2[mat][j][col];
1936 }
1937 assert( cnt == nrows2 - 1);
1938 SCIPsortIntInt(sortvals, &((*doublelexmatrix)[0][1]), cnt);
1939
1940 /* fill the remaining entries of the big matrix */
1941 for (i = 1; i < *nrows; ++i)
1942 {
1943 for (j = 1; j < *ncols; ++j)
1944 {
1945 /* get the matrices and column/row of the entry */
1946 mat = idxtomatrix1[(*doublelexmatrix)[0][j]];
1947 mat2 = idxtomatrix2[(*doublelexmatrix)[i][0]];
1948 col = idxtocol1[(*doublelexmatrix)[0][j]];
1949 col2 = idxtocol2[(*doublelexmatrix)[i][0]];
1950
1951 /* find the unique element in the col column of matrix mat and the row column of matrix mat2 */
1952 /* @todo improve this by first sorting the columns */
1953 cnt = 0;
1954 elem = -1;
1955 for (c = 0; c < *nrows; ++c)
1956 {
1957 for (d = 0; d < *ncols; ++d)
1958 {
1959 if ( matrices1[mat][c][col] == matrices2[mat2][d][col2] )
1960 {
1961 ++cnt;
1962 elem = matrices1[mat][c][col];
1963 break;
1964 }
1965 }
1966 }
1967
1968 /* stop: the columns do not overlap properly */
1969 if ( cnt != 1 )
1970 {
1971 *success = FALSE;
1972 goto FREEMEMORY;
1973 }
1974 (*doublelexmatrix)[i][j] = elem;
1975 }
1976 }
1977
1978 /* store begin positions of row and column blocks */
1979 SCIP_CALL( SCIPallocBlockMemoryArray(scip, rowsbegin, nmatrices2 + 1) );
1980 SCIP_CALL( SCIPallocBlockMemoryArray(scip, colsbegin, nmatrices1 + 1) );
1981 (*rowsbegin)[0] = 0;
1982 (*colsbegin)[0] = 0;
1983 for (j = 0; j < nmatrices2; ++j)
1984 (*rowsbegin)[j + 1] = (*rowsbegin)[j] + ncols2[j];
1985 for (j = 0; j < nmatrices1; ++j)
1986 (*colsbegin)[j + 1] = (*colsbegin)[j] + ncols1[j];
1987
1988 FREEMEMORY:
1989 SCIPfreeBufferArray(scip, &sortvals);
1990
1991 SCIPfreeBufferArray(scip, &idxtocol2);
1992 SCIPfreeBufferArray(scip, &idxtocol1);
1993 SCIPfreeBufferArray(scip, &idxtorow2);
1994 SCIPfreeBufferArray(scip, &idxtorow1);
1995 SCIPfreeBufferArray(scip, &idxtomatrix2);
1996 SCIPfreeBufferArray(scip, &idxtomatrix1);
1997
1998 if ( !(*success) )
1999 {
2000 for (i = *nrows - 1; i >= 0; --i)
2001 {
2002 SCIPfreeBlockMemoryArray(scip, &(*doublelexmatrix)[i], *ncols);
2003 }
2004 SCIPfreeBlockMemoryArray(scip, doublelexmatrix, *nrows);
2005 SCIPfreeBlockMemoryArray(scip, rowsbegin, nmatrices2 + 1);
2006 SCIPfreeBlockMemoryArray(scip, colsbegin, nmatrices1 + 1);
2007 *doublelexmatrix = NULL;
2008 *rowsbegin = NULL;
2009 *colsbegin = NULL;
2010 }
2011
2012 return SCIP_OKAY;
2013 }
2014
2015 /** detects whether permutations define single or double lex matrices
2016 *
2017 * A single lex matrix is a matrix whose columns can be partitioned into blocks such that the
2018 * columns within each block can be permuted arbitrarily. A double lex matrix is a single lex
2019 * matrix such that also blocks of rows have the aforementioned property.
2020 */
2021 SCIP_RETCODE SCIPdetectSingleOrDoubleLexMatrices(
2022 SCIP* scip, /**< SCIP pointer */
2023 SCIP_Bool detectsinglelex, /**< whether single lex matrices shall be detected */
2024 int** perms, /**< array of permutations */
2025 int nperms, /**< number of permutations in perms */
2026 int permlen, /**< number of variables in a permutation */
2027 SCIP_Bool* success, /**< pointer to store whether structure could be detected */
2028 SCIP_Bool* isorbitope, /**< pointer to store whether detected matrix is orbitopal */
2029 int*** lexmatrix, /**< pointer to store single or double lex matrix */
2030 int* nrows, /**< pointer to store number of rows of lexmatrix */
2031 int* ncols, /**< pointer to store number of columns of lexmatrix */
2032 int** lexrowsbegin, /**< pointer to store array indicating begin of new row-lexmatrix */
2033 int** lexcolsbegin, /**< pointer to store array indicating begin of new col-lexmatrix */
2034 int* nrowmatrices, /**< pointer to store number of single lex row matrices in rows */
2035 int* ncolmatrices /**< pointer to store number of single lex column matrices in rows */
2036 )
2037 {
2038 int*** matricestype1 = NULL;
2039 int*** matricestype2 = NULL;
2040 int* ncolstype1 = NULL;
2041 int* ncolstype2 = NULL;
2042 int nmatricestype1 = 0;
2043 int nmatricestype2 = 0;
2044 int* permstype1;
2045 int* permstype2;
2046 int npermstype1 = 0;
2047 int npermstype2 = 0;
2048 int ncycs1 = -1;
2049 int ncycs2 = -1;
2050 int tmpncycs;
2051 int p;
2052 int i;
2053 SCIP_Bool isinvolution;
2054
2055 assert( scip != NULL );
2056 assert( perms != NULL );
2057 assert( nperms > 0 );
2058 assert( permlen > 0 );
2059 assert( success != NULL );
2060 assert( lexmatrix != NULL );
2061 assert( nrows != NULL );
2062 assert( ncols != NULL );
2063 assert( lexrowsbegin != NULL );
2064 assert( lexcolsbegin != NULL );
2065 assert( nrowmatrices != NULL );
2066 assert( ncolmatrices != NULL );
2067
2068 *success = TRUE;
2069 *isorbitope = FALSE;
2070 *nrowmatrices = 0;
2071 *ncolmatrices = 0;
2072
2073 /* arrays to store the different types of involutions */
2074 SCIP_CALL( SCIPallocBufferArray(scip, &permstype1, nperms) );
2075 SCIP_CALL( SCIPallocBufferArray(scip, &permstype2, nperms) );
2076
2077 /* check whether we can expect lexicographically sorted rows and columns */
2078 for (p = 0; p < nperms; ++p)
2079 {
2080 SCIP_CALL( isPermInvolution(perms[p], permlen, &isinvolution, &tmpncycs) );
2081
2082 /* terminate if not all permutations are involutions */
2083 if ( ! isinvolution )
2084 {
2085 *success = FALSE;
2086 goto FREEMEMORY;
2087 }
2088
2089 /* store number of cycles or terminate if too many different types of involutions */
2090 if ( ncycs1 == -1 || ncycs1 == tmpncycs )
2091 {
2092 ncycs1 = tmpncycs;
2093 permstype1[npermstype1++] = p;
2094 }
2095 else if ( ncycs2 == -1 || ncycs2 == tmpncycs )
2096 {
2097 ncycs2 = tmpncycs;
2098 permstype2[npermstype2++] = p;
2099 }
2100 else
2101 {
2102 *success = FALSE;
2103 goto FREEMEMORY;
2104 }
2105 }
2106
2107 /* for each type, check whether permutations define (disjoint) orbitopal symmetries */
2108 SCIP_CALL( detectOrbitopalSymmetries(scip, perms, permstype1, npermstype1, permlen, ncycs1, success,
2109 &matricestype1, &ncolstype1, &nmatricestype1) );
2110 if ( ! *success )
2111 goto FREEMEMORY;
2112
2113 SCIP_CALL( detectOrbitopalSymmetries(scip, perms, permstype2, npermstype2, permlen, ncycs2, success,
2114 &matricestype2, &ncolstype2, &nmatricestype2) );
2115 if ( ! *success )
2116 goto FREEMEMORY;
2117
2118 /* check whether a double lex matrix is defined */
2119 *success = FALSE;
2120 if ( !detectsinglelex && ncycs2 != -1 )
2121 {
2122 assert( ncycs1 > 0 );
2123
2124 SCIP_CALL( isDoublelLexSym(scip, matricestype1, ncycs1, ncolstype1, nmatricestype1,
2125 matricestype2, ncycs2, ncolstype2, nmatricestype2,
2126 lexmatrix, nrows, ncols, lexrowsbegin, lexcolsbegin, success) );
2127
2128 if ( *success )
2129 {
2130 *nrowmatrices = nmatricestype2;
2131 *ncolmatrices = nmatricestype1;
2132 }
2133 }
2134
2135 /* if no double lex matrix is detected, possibly return orbitope */
2136 if ( !(*success) && ncycs2 == -1 && nmatricestype1 == 1 )
2137 {
2138 SCIP_CALL( SCIPallocBlockMemoryArray(scip, lexmatrix, ncycs1) );
2139 for (i = 0; i < ncycs1; ++i)
2140 {
2141 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*lexmatrix)[i], ncolstype1[0]) );
2142 for (p = 0; p < ncolstype1[0]; ++p)
2143 (*lexmatrix)[i][p] = matricestype1[0][i][p];
2144 }
2145 *nrows = ncycs1;
2146 *ncols = ncolstype1[0];
2147 *success = TRUE;
2148 *isorbitope = TRUE;
2149 }
2150 #ifndef NDEBUG
2151 else if ( !(*success) )
2152 {
2153 assert( *lexmatrix == NULL );
2154 assert( *lexrowsbegin == NULL );
2155 assert( *lexcolsbegin == NULL );
2156 }
2157 #endif
2158
2159 FREEMEMORY:
2160 for (p = nmatricestype2 - 1; p >= 0; --p)
2161 {
2162 for (i = ncycs2 - 1; i >= 0; --i)
2163 {
2164 SCIPfreeBlockMemoryArray(scip, &matricestype2[p][i], ncolstype2[p]);
2165 }
2166 SCIPfreeBlockMemoryArray(scip, &matricestype2[p], ncycs2);
2167 }
2168 SCIPfreeBlockMemoryArrayNull(scip, &matricestype2, nmatricestype2);
2169 SCIPfreeBlockMemoryArrayNull(scip, &ncolstype2, nmatricestype2);
2170 for (p = nmatricestype1 - 1; p >= 0; --p)
2171 {
2172 for (i = ncycs1 - 1; i >= 0; --i)
2173 {
2174 SCIPfreeBlockMemoryArray(scip, &matricestype1[p][i], ncolstype1[p]);
2175 }
2176 SCIPfreeBlockMemoryArray(scip, &matricestype1[p], ncycs1);
2177 }
2178 SCIPfreeBlockMemoryArrayNull(scip, &matricestype1, nmatricestype1);
2179 SCIPfreeBlockMemoryArrayNull(scip, &ncolstype1, nmatricestype1);
2180
2181 SCIPfreeBufferArray(scip, &permstype2);
2182 SCIPfreeBufferArray(scip, &permstype1);
2183
2184 return SCIP_OKAY;
2185 }
2186
2187
2188 /** helper function to test if val1 = val2 while permitting infinity-values */
2189 SCIP_Bool SCIPEQ(
2190 SCIP* scip, /**< SCIP data structure */
2191 SCIP_Real val1, /**< left-hand side value */
2192 SCIP_Real val2 /**< right-hand side value */
2193 )
2194 {
2195 SCIP_Bool inf1;
2196 SCIP_Bool inf2;
2197 SCIP_Bool minf1;
2198 SCIP_Bool minf2;
2199
2200 inf1 = SCIPisInfinity(scip, val1);
2201 inf2 = SCIPisInfinity(scip, val2);
2202 if ( inf1 && inf2 )
2203 return TRUE;
2204 if ( inf1 != inf2 )
2205 return FALSE;
2206 assert( !inf1 );
2207 assert( !inf2 );
2208
2209 minf1 = SCIPisInfinity(scip, -val1);
2210 minf2 = SCIPisInfinity(scip, -val2);
2211 if ( minf1 && minf2 )
2212 return TRUE;
2213 if ( minf1 != minf2 )
2214 return FALSE;
2215 assert( !minf1 );
2216 assert( !minf2 );
2217
2218 return SCIPisEQ(scip, val1, val2);
2219 }
2220
2221
2222 /** helper function to test if val1 <= val2 while permitting infinity-values */
2223 SCIP_Bool SCIPLE(
2224 SCIP* scip, /**< SCIP data structure */
2225 SCIP_Real val1, /**< left-hand side value */
2226 SCIP_Real val2 /**< right-hand side value */
2227 )
2228 {
2229 SCIP_Bool inf1;
2230 SCIP_Bool inf2;
2231 SCIP_Bool minf1;
2232 SCIP_Bool minf2;
2233
2234 inf1 = SCIPisInfinity(scip, val1);
2235 inf2 = SCIPisInfinity(scip, val2);
2236 if ( inf1 && inf2 )
2237 return TRUE;
2238 if ( !inf1 && inf2 )
2239 return TRUE;
2240 if ( inf1 && !inf2 )
2241 return FALSE;
2242 assert( !inf1 );
2243 assert( !inf2 );
2244
2245 minf1 = SCIPisInfinity(scip, -val1);
2246 minf2 = SCIPisInfinity(scip, -val2);
2247 if ( minf1 && minf2 )
2248 return TRUE;
2249 if ( !minf1 && minf2 )
2250 return FALSE;
2251 if ( minf1 && !minf2 )
2252 return TRUE;
2253 assert( !minf1 );
2254 assert( !minf2 );
2255
2256 return SCIPisLE(scip, val1, val2);
2257 }
2258
2259
2260 /** helper function to test if val1 >= val2 while permitting infinity-values */
2261 SCIP_Bool SCIPGE(
2262 SCIP* scip, /**< SCIP data structure */
2263 SCIP_Real val1, /**< left-hand side value */
2264 SCIP_Real val2 /**< right-hand side value */
2265 )
2266 {
2267 SCIP_Bool inf1;
2268 SCIP_Bool inf2;
2269 SCIP_Bool minf1;
2270 SCIP_Bool minf2;
2271
2272 inf1 = SCIPisInfinity(scip, val1);
2273 inf2 = SCIPisInfinity(scip, val2);
2274 if ( inf1 && inf2 )
2275 return TRUE;
2276 if ( !inf1 && inf2 )
2277 return FALSE;
2278 if ( inf1 && !inf2 )
2279 return TRUE;
2280 assert( !inf1 );
2281 assert( !inf2 );
2282
2283 minf1 = SCIPisInfinity(scip, -val1);
2284 minf2 = SCIPisInfinity(scip, -val2);
2285 if ( minf1 && minf2 )
2286 return TRUE;
2287 if ( !minf1 && minf2 )
2288 return TRUE;
2289 if ( minf1 && !minf2 )
2290 return FALSE;
2291 assert( !minf1 );
2292 assert( !minf2 );
2293
2294 return SCIPisGE(scip, val1, val2);
2295 }
2296
2297
2298 /** helper function to test if val1 < val2 while permitting infinity-values */
2299 SCIP_Bool SCIPLT(
2300 SCIP* scip, /**< SCIP data structure */
2301 SCIP_Real val1, /**< left-hand side value */
2302 SCIP_Real val2 /**< right-hand side value */
2303 )
2304 {
2305 SCIP_Bool inf1;
2306 SCIP_Bool inf2;
2307 SCIP_Bool minf1;
2308 SCIP_Bool minf2;
2309
2310 inf1 = SCIPisInfinity(scip, val1);
2311 inf2 = SCIPisInfinity(scip, val2);
2312 if ( inf1 && inf2 )
2313 return FALSE;
2314 if ( !inf1 && inf2 )
2315 return TRUE;
2316 if ( inf1 && !inf2 )
2317 return FALSE;
2318 assert( !inf1 );
2319 assert( !inf2 );
2320
2321 minf1 = SCIPisInfinity(scip, -val1);
2322 minf2 = SCIPisInfinity(scip, -val2);
2323 if ( minf1 && minf2 )
2324 return FALSE;
2325 if ( !minf1 && minf2 )
2326 return FALSE;
2327 if ( minf1 && !minf2 )
2328 return TRUE;
2329 assert( !minf1 );
2330 assert( !minf2 );
2331
2332 return SCIPisLT(scip, val1, val2);
2333 }
2334
2335
2336 /** helper function to test if val1 > val2 while permitting infinity-values */
2337 SCIP_Bool SCIPGT(
2338 SCIP* scip, /**< SCIP data structure */
2339 SCIP_Real val1, /**< left-hand side value */
2340 SCIP_Real val2 /**< right-hand side value */
2341 )
2342 {
2343 SCIP_Bool inf1;
2344 SCIP_Bool inf2;
2345 SCIP_Bool minf1;
2346 SCIP_Bool minf2;
2347
2348 inf1 = SCIPisInfinity(scip, val1);
2349 inf2 = SCIPisInfinity(scip, val2);
2350 if ( inf1 && inf2 )
2351 return FALSE;
2352 if ( !inf1 && inf2 )
2353 return FALSE;
2354 if ( inf1 && !inf2 )
2355 return TRUE;
2356 assert( !inf1 );
2357 assert( !inf2 );
2358
2359 minf1 = SCIPisInfinity(scip, -val1);
2360 minf2 = SCIPisInfinity(scip, -val2);
2361 if ( minf1 && minf2 )
2362 return FALSE;
2363 if ( !minf1 && minf2 )
2364 return TRUE;
2365 if ( minf1 && !minf2 )
2366 return FALSE;
2367 assert( !minf1 );
2368 assert( !minf2 );
2369
2370 return SCIPisGT(scip, val1, val2);
2371 }
2372