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 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 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