Scippy

SCIP

Solving Constraint Integer Programs

cons_orbitope.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2014 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not email to scip@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file cons_orbitope.c
17  * @brief constraint handler for (partitioning/packing) orbitope constraints w.r.t. the full symmetric group
18  * @author Timo Berthold
19  * @author Marc Pfetsch
20  *
21  * The type of constraints of this constraint handler is described in cons_orbitope.h.
22  *
23  * The details of the method implemented here are described in the following papers.
24  *
25  * Packing and Partitioning Orbitopes@n
26  * Volker Kaibel and Marc E. Pfetsch,@n
27  * Math. Program. 114, No. 1, 1-36 (2008)
28  *
29  * Among other things, this paper describes so-called shifted column inequalities of the following
30  * form \f$x(S) \leq x(B)\f$, where \f$S\f$ is a so-called shifted column and \f$B\f$ is a so-called
31  * bar. These inequalities can be used to handle symmetry and they are separated in this constraint
32  * handler. We use the linear time separation algorithm of the paper.@par
33  *
34  * Orbitopal Fixing@n
35  * Volker Kaibel, Matthias Peinhardt, and Marc E. Pfetsch,@n
36  * Discrete Optimization 8, No. 4, 595-610 (2011)
37  * (A preliminary version appears in Proc. IPCO 2007.)
38  *
39  * In this paper a linear time propagation algorithm is described, a variant of which is implemented
40  * here. The implemented variant does not run in linear time, but is very fast in practice.
41  *
42  * <table>
43  * <caption>translation table</caption>
44  * <tr><td>here</td><td>paper</td></tr>
45  * <tr><td></td><td></td></tr>
46  * <tr><td>nspcons </td><td>p </td></tr>
47  * <tr><td>nblocks </td><td>q </td></tr>
48  * <tr><td>vars </td><td>x </td></tr>
49  * <tr><td>vals </td><td>A^\\star</td></tr>
50  * <tr><td>weights </td><td>\\omega </td></tr>
51  * <tr><td>cases </td><td>\\tau </td></tr>
52  * <tr><td>fixtriangle </td><td>-- </td></tr>
53  * <tr><td>resolveprop </td><td>-- </td></tr>
54  * <tr><td>firstnonzeros</td><td>\\mu </td></tr>
55  * <tr><td>lastones </td><td>\\alpha </td></tr>
56  * <tr><td>frontiersteps</td><td>\\Gamma </td></tr>
57  * </table>
58  */
59 
60 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
61 
62 #include <assert.h>
63 #include <string.h>
64 #include <ctype.h>
65 
66 #include "scip/cons_orbitope.h"
67 
68 /* constraint handler properties */
69 #define CONSHDLR_NAME "orbitope"
70 #define CONSHDLR_DESC "symmetry breaking constraint handler relying on (partitioning/packing) orbitopes"
71 #define CONSHDLR_SEPAPRIORITY +40100 /**< priority of the constraint handler for separation */
72 #define CONSHDLR_ENFOPRIORITY -1005200 /**< priority of the constraint handler for constraint enforcing */
73 #define CONSHDLR_CHECKPRIORITY -1005200 /**< priority of the constraint handler for checking feasibility */
74 #define CONSHDLR_SEPAFREQ 5 /**< frequency for separating cuts; zero means to separate only in the root node */
75 #define CONSHDLR_PROPFREQ -1 /**< frequency for propagating domains; zero means only preprocessing propagation */
76 #define CONSHDLR_EAGERFREQ -1 /**< frequency for using all instead of only the useful constraints in separation,
77  * propagation and enforcement, -1 for no eager evaluations, 0 for first only */
78 #define CONSHDLR_MAXPREROUNDS -1 /**< maximal number of presolving rounds the constraint handler participates in (-1: no limit) */
79 #define CONSHDLR_DELAYSEPA FALSE /**< should separation method be delayed, if other separators found cuts? */
80 #define CONSHDLR_DELAYPROP FALSE /**< should propagation method be delayed, if other propagators found reductions? */
81 #define CONSHDLR_DELAYPRESOL FALSE /**< should presolving method be delayed, if other presolvers found reductions? */
82 #define CONSHDLR_NEEDSCONS TRUE /**< should the constraint handler be skipped, if no constraints are available? */
83 
84 #define CONSHDLR_PROP_TIMING SCIP_PROPTIMING_BEFORELP
85 
86 
87 /*
88  * Data structures
89  */
90 
91 /** constraint data for orbitope constraints */
92 struct SCIP_ConsData
93 {
94  SCIP_VAR*** vars; /**< matrix of variables on which the symmetry acts */
95  SCIP_VAR** tmpvars; /**< temporary storage for variables */
96  SCIP_Real** vals; /**< LP-solution for those variables */
97  SCIP_Real* tmpvals; /**< temporary storage for values */
98  SCIP_Real** weights; /**< SC weight table */
99  int** cases; /**< indicator of the SC cases */
100  int nspcons; /**< number of set partitioning/packing constraints <=> p */
101  int nblocks; /**< number of symmetric variable blocks <=> q */
102  SCIP_Bool ispart; /**< whether we deal with the partitioning case (packing otherwise) */
103  SCIP_Bool resolveprop; /**< should propagation be resolved? */
104  SCIP_Bool istrianglefixed; /**< has the upper right triangle already globally been fixed to zero? */
105 };
106 
107 
108 /*
109  * Local methods
110  */
111 
112 /** frees an orbitope constraint data */
113 static
114 SCIP_RETCODE consdataFree(
115  SCIP* scip, /**< SCIP data structure */
116  SCIP_CONSDATA** consdata /**< pointer to orbitope constraint data */
117  )
118 {
119  int i;
120  int p;
121  int q;
122 
123  assert( consdata != NULL );
124  assert( *consdata != NULL );
125 
126  p = (*consdata)->nspcons;
127  q = (*consdata)->nblocks;
128  for (i = 0; i < p; ++i)
129  {
130  SCIPfreeBlockMemoryArrayNull(scip, &((*consdata)->cases[i]), q); /*lint !e866*/
131  SCIPfreeBlockMemoryArrayNull(scip, &((*consdata)->vars[i]), q); /*lint !e866*/
132  SCIPfreeBlockMemoryArrayNull(scip, &((*consdata)->weights[i]), q); /*lint !e866*/
133  SCIPfreeBlockMemoryArrayNull(scip, &((*consdata)->vals[i]), q); /*lint !e866*/
134  }
135 
136  SCIPfreeBlockMemoryArrayNull(scip, &((*consdata)->cases), p);
137  SCIPfreeBlockMemoryArrayNull(scip, &((*consdata)->vars), p);
138  SCIPfreeBlockMemoryArrayNull(scip, &((*consdata)->weights), p);
139  SCIPfreeBlockMemoryArrayNull(scip, &((*consdata)->vals), p);
140 
141  SCIPfreeBlockMemoryArrayNull(scip, &((*consdata)->tmpvals), p + q);
142  SCIPfreeBlockMemoryArrayNull(scip, &((*consdata)->tmpvars), p + q);
143 
144  SCIPfreeBlockMemory(scip, consdata);
145 
146  return SCIP_OKAY;
147 }
148 
149 
150 /** creates orbitope constraint data */
151 static
152 SCIP_RETCODE consdataCreate(
153  SCIP* scip, /**< SCIP data structure */
154  SCIP_CONSDATA** consdata, /**< pointer to store constraint data */
155  SCIP_VAR*** vars, /**< variables array, must have size nspcons x nblocks */
156  int nspcons, /**< number of set partitioning (packing) constraints <=> p */
157  int nblocks, /**< number of symmetric variable blocks <=> q */
158  SCIP_Bool ispart, /**< deal with the partitioning case (packing otherwise) */
159  SCIP_Bool resolveprop /**< should propagation be resolved? */
160  )
161 {
162  int i;
163 
164  assert(consdata != NULL);
165 
166  SCIP_CALL( SCIPallocBlockMemory(scip, consdata) );
167 
168  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*consdata)->vals, nspcons) );
169  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*consdata)->weights, nspcons) );
170  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*consdata)->vars, nspcons) );
171  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*consdata)->cases, nspcons) );
172 
173  for (i = 0; i < nspcons; ++i)
174  {
175  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*consdata)->vals[i], nblocks) ); /*lint !e866*/
176  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*consdata)->weights[i], nblocks) ); /*lint !e866*/
177  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->vars[i], vars[i], nblocks) ); /*lint !e866*/
178  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*consdata)->cases[i], nblocks) ); /*lint !e866*/
179  }
180 
181  (*consdata)->tmpvals = NULL;
182  (*consdata)->tmpvars = NULL;
183  (*consdata)->nspcons = nspcons;
184  (*consdata)->nblocks = nblocks;
185  (*consdata)->ispart = ispart;
186  (*consdata)->resolveprop = resolveprop;
187  (*consdata)->istrianglefixed = FALSE;
188 
189  /* get transformed variables, if we are in the transformed problem */
190  if ( SCIPisTransformed(scip) )
191  {
192  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*consdata)->tmpvals, nspcons + nblocks) );
193  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*consdata)->tmpvars, nspcons + nblocks) );
194 
195  for (i = 0; i < nspcons; ++i)
196  {
197  SCIP_CALL( SCIPgetTransformedVars(scip, (*consdata)->nblocks, (*consdata)->vars[i], (*consdata)->vars[i]) );
198  }
199  }
200 
201  return SCIP_OKAY;
202 }
203 
204 
205 #ifdef PRINT_MATRIX
206 /** debug method, prints variable matrix */
207 static
208 void printMatrix(
209  SCIP* scip, /**< SCIP data structure */
210  SCIP_CONSDATA* consdata /**< the constraint data */
211  )
212 {
213  int i;
214  int j;
215 
216  assert( consdata != NULL );
217  assert( consdata->nspcons > 0 );
218  assert( consdata->nblocks > 0 );
219  assert( consdata->vars != NULL );
220 
221  for (j = 0; j < consdata->nblocks; ++j)
222  SCIPinfoMessage(scip, NULL, "-");
223 
224  SCIPinfoMessage(scip, NULL, "\n");
225  for (i = 0; i < consdata->nspcons; ++i)
226  {
227  for (j = 0; j < consdata->nblocks; ++j)
228  {
229  if ( SCIPvarGetUbLocal(consdata->vars[i][j]) - SCIPvarGetLbLocal(consdata->vars[i][j]) < 0.5 )
230  SCIPinfoMessage(scip, NULL, "%1.0f", REALABS(SCIPvarGetUbLocal(consdata->vars[i][j])));
231  else
232  SCIPinfoMessage(scip, NULL, " ");
233  }
234  SCIPinfoMessage(scip, NULL, "|\n");
235  }
236  for (j = 0; j < consdata->nblocks; ++j)
237  SCIPinfoMessage(scip, NULL, "-");
238  SCIPinfoMessage(scip, NULL, "\n");
239 }
240 #endif
241 
242 
243 #ifdef SHOW_SCI
244 /** Print SCI in nice form for debugging */
245 static
246 SCIP_RETCODE printSCI(
247  SCIP* scip, /**< SCIP pointer */
248  int p, /**< number of rows */
249  int q, /**< number of columns */
250  int** cases, /**< SCI dynamic programming table */
251  int i, /**< row position of bar */
252  int j /**< column position of bar */
253  )
254 {
255  int k;
256  int l;
257  int** M;
258  int p1;
259  int p2;
260 
261  SCIP_CALL( SCIPallocBufferArray(scip, &M, p) );
262  for (k = 0; k < p; ++k)
263  {
264  SCIP_CALL( SCIPallocBufferArray(scip, &M[k], q) ); /*lint !e866*/
265  for (l = 0; l < q; ++l)
266  M[k][l] = 0;
267  }
268 
269  /* first add bar */
270  for (l = j; l < q; ++l)
271  {
272  assert( M[i][l] == 0 );
273  M[i][l] = 1;
274  }
275 
276  /* then add shifted column */
277  p1 = i-1;
278  p2 = j-1;
279  do
280  {
281  assert( cases[p1][p2] != -1 );
282  assert( p1 >= 0 && p1 < i );
283  assert( p2 >= 0 && p2 < j );
284 
285  /* if case 1 */
286  if ( cases[p1][p2] == 1 )
287  --p2; /* decrease column */
288  else
289  {
290  /* case 2 or 3: */
291  assert( cases[p1][p2] == 2 || cases[p1][p2] == 3 );
292  assert( M[p1][p2] == 0 );
293  M[p1][p2] = -1;
294  if ( cases[p1][p2] == 3 )
295  break;
296  }
297  --p1; /* decrease row */
298  }
299  while ( p1 >= 0 ); /* should always be true, i.e. the break should end the loop */
300  assert( cases[p1][p2] == 3 );
301 
302  /* now output matrix M */
303  for (l = 0; l < q; ++l)
304  SCIPinfoMessage(scip, NULL, "-");
305  SCIPinfoMessage(scip, NULL, "\n");
306 
307  for (k = 0; k < p; ++k)
308  {
309  for (l = 0; l < q; ++l)
310  {
311  if ( l > k )
312  SCIPinfoMessage(scip, NULL, "*");
313  else
314  {
315  switch (M[k][l])
316  {
317  case 1:
318  SCIPinfoMessage(scip, NULL, "+");
319  break;
320  case -1:
321  SCIPinfoMessage(scip, NULL, "-");
322  break;
323  case 0:
324  SCIPinfoMessage(scip, NULL, "#");
325  break;
326  default:
327  SCIPerrorMessage("unexpected matrix entry <%d>: should be -1, 0 or +1\n", M[k][l]);
328  SCIPABORT();
329  }
330  }
331  }
332  SCIPinfoMessage(scip, NULL, "\n");
333  }
334 
335  for (l = 0; l < q; ++l)
336  SCIPinfoMessage(scip, NULL, "-");
337  SCIPinfoMessage(scip, NULL, "\n");
338 
339  for (k = 0; k < p; ++k)
340  SCIPfreeBufferArray(scip, &M[k]);
341  SCIPfreeBufferArray(scip, &M);
342 
343  return SCIP_OKAY;
344 }
345 #endif
346 
347 
348 /** copies the variables values from the solution to the constraint data structure */
349 static
350 void copyValues(
351  SCIP* scip, /**< the SCIP data structure */
352  SCIP_CONSDATA* consdata, /**< the constraint data */
353  SCIP_SOL* sol /**< a primal solution or NULL for the current LP optimum */
354  )
355 {
356  int i;
357  int j;
358 
359  assert( scip != NULL );
360  assert( consdata != NULL );
361  assert( consdata->nspcons > 0 );
362  assert( consdata->nblocks > 0 );
363  assert( consdata->vars != NULL );
364  assert( consdata->vals != NULL );
365 
366  for (i = 0; i < consdata->nspcons; ++i)
367  {
368  for (j = 0; j < consdata->nblocks; ++j)
369  consdata->vals[i][j] = SCIPgetSolVal(scip, sol, consdata->vars[i][j]);
370  }
371 }
372 
373 
374 /** compute the dynamic programming table for SC
375  *
376  * Build up dynamic programming table in order to find SCs with minimum weight.
377  *
378  * The values of the minimal SCIs are stored in @a weights.
379  * The array @a cases[i][j] stores which of the cases were applied to get @a weights[i][j].
380  * Here, 3 means that we have reached the upper limit.
381  *
382  * We assume that the upper right triangle is fixed to 0. Hence we can perform the computation a
383  * bit more efficient.
384  */
385 static
386 void computeSCTable(
387  SCIP* scip, /**< SCIP pointer */
388  int nspcons, /**< number of set partitioning (packing) constraints <=> p */
389  int nblocks, /**< number of symmetric variable blocks <=> q */
390  SCIP_Real** weights, /**< SC weight table */
391  int** cases, /**< indicator of the SC cases */
392  SCIP_Real** vals /**< current solution */
393  )
394 {
395  SCIP_Real minvalue;
396  int diagsize;
397  int i;
398  int j;
399 
400  assert( weights != NULL );
401  assert( cases != NULL );
402  assert( vals != NULL );
403 
404 #ifndef NDEBUG
405  /* for debugging */
406  for (i = 0; i < nspcons; ++i)
407  {
408  for (j = 0; j < nblocks; ++j)
409  {
410  if ( i >= j )
411  {
412  weights[i][j] = -1.0;
413  cases[i][j] = -1;
414  }
415  }
416  }
417 #endif
418 
419  /* initialize diagonal */
420  minvalue = vals[0][0];
421  weights[0][0] = minvalue;
422  cases[0][0] = 3;
423 
424  /* get last row of triangle */
425  diagsize = nblocks;
426  if ( nspcons < nblocks )
427  diagsize = nspcons;
428 
429  for (j = 1; j < diagsize; ++j)
430  {
431  /* use LT to move entry as far to the left as possible */
432  if ( SCIPisLT(scip, vals[j][j], minvalue) )
433  {
434  minvalue = vals[j][j];
435  cases[j][j] = 3;
436  }
437  else
438  cases[j][j] = 1;
439  weights[j][j] = minvalue;
440  }
441 
442  /* initialize first column */
443  for (i = 1; i < nspcons; ++i)
444  {
445  weights[i][0] = weights[i-1][0] + vals[i][0];
446  cases[i][0] = 2; /* second case */
447  }
448 
449  /* build the table */
450  for (i = 2; i < nspcons; ++i)
451  {
452  for (j = 1; j < nblocks && j < i; ++j)
453  {
454  SCIP_Real weightleft;
455  SCIP_Real weightright;
456 
457  assert( cases[i-1][j] != -1 );
458  assert( cases[i-1][j-1] != -1 );
459 
460  weightleft = weights[i-1][j-1];
461  weightright = vals[i][j] + weights[i-1][j];
462 
463  /* For first column: cannot take left possibility */
464  if ( SCIPisLT(scip, weightleft, weightright) )
465  {
466  weights[i][j] = weightleft;
467  cases[i][j] = 1;
468  }
469  else
470  {
471  weights[i][j] = weightright;
472  cases[i][j] = 2;
473  }
474  }
475  }
476 }
477 
478 
479 /** fix upper right triangle if necessary */
480 static
481 SCIP_RETCODE fixTriangle(
482  SCIP* scip, /**< SCIP data structure */
483  SCIP_CONS* cons, /**< constraint to be processed */
484  SCIP_Bool* infeasible, /**< pointer to store TRUE, if the node can be cut off */
485  int* nfixedvars /**< pointer to add up the number of found domain reductions */
486  )
487 {
488  SCIP_CONSDATA* consdata;
489  SCIP_VAR*** vars;
490  SCIP_Bool fixedglobal;
491  SCIP_Bool fixed;
492  int diagsize;
493  int nspcons;
494  int nblocks;
495  int i;
496  int j;
497 
498  assert( scip != NULL );
499  assert( cons != NULL );
500  assert( infeasible != NULL );
501  assert( nfixedvars != NULL );
502 
503  consdata = SCIPconsGetData(cons);
504  assert( consdata != NULL );
505  assert( consdata->nspcons > 0 );
506  assert( consdata->nblocks > 0 );
507  assert( consdata->vars != NULL );
508 
509  *infeasible = FALSE;
510  *nfixedvars = 0;
511 
512  if ( consdata->istrianglefixed )
513  return SCIP_OKAY;
514 
515  nspcons = consdata->nspcons;
516  nblocks = consdata->nblocks;
517  vars = consdata->vars;
518  fixedglobal = TRUE;
519 
520  /* get last row of triangle */
521  diagsize = nblocks;
522  if ( nspcons < nblocks )
523  diagsize = nspcons;
524 
525  /* fix variables to 0 */
526  for (i = 0; i < diagsize; ++i)
527  {
528  for (j = i+1; j < nblocks; ++j)
529  {
530  /* fix variable, if not in the root the fixation is local */
531  SCIP_CALL( SCIPfixVar(scip, vars[i][j], 0.0, infeasible, &fixed) );
532 
533  if ( *infeasible )
534  {
535  SCIPdebugMessage("The problem is infeasible: some variable in the upper right triangle is fixed to 1.\n");
536  return SCIP_OKAY;
537  }
538 
539  if ( fixed )
540  ++(*nfixedvars);
541 
542  if ( SCIPvarGetUbGlobal(vars[i][j]) > 0.5 )
543  fixedglobal = FALSE;
544  }
545  }
546  if ( *nfixedvars > 0 )
547  {
548  SCIPdebugMessage("<%s>: %s fixed upper right triangle to 0 (fixed vars: %d).\n", SCIPconsGetName(cons), fixedglobal ? "globally" : "locally", *nfixedvars);
549  }
550  else
551  {
552  SCIPdebugMessage("<%s>: Upper right triangle already fixed to 0.\n", SCIPconsGetName(cons));
553  }
554 
555  if ( fixedglobal )
556  consdata->istrianglefixed = TRUE;
557 
558  return SCIP_OKAY;
559 }
560 
561 
562 /** separates shifted column inequalities according to the solution stored in consdata->vals */
563 static
564 SCIP_RETCODE separateSCIs(
565  SCIP* scip, /**< the SCIP data structure */
566  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
567  SCIP_CONS* cons, /**< constraint */
568  SCIP_CONSDATA* consdata, /**< the constraint data */
569  SCIP_Bool* infeasible, /**< whether we detected infeasibility */
570  int* nfixedvars, /**< number of variables fixed */
571  int* ncuts /**< pointer to store number of separated SCIs */
572  )
573 {
574  SCIP_Real** vals;
575  SCIP_Real** weights;
576  SCIP_Real* tmpvals;
577  SCIP_VAR*** vars;
578  SCIP_VAR** tmpvars;
579  int** cases;
580  int nspcons;
581  int nblocks;
582  int i;
583  int j;
584  int l;
585 
586  assert( scip != NULL );
587  assert( conshdlr != NULL );
588  assert( cons != NULL );
589  assert( infeasible != NULL);
590  assert( nfixedvars != NULL );
591  assert( ncuts != NULL );
592 
593  assert( consdata != NULL );
594  assert( consdata->nspcons > 0 );
595  assert( consdata->nblocks > 0 );
596  assert( consdata->vars != NULL );
597  assert( consdata->vals != NULL );
598  assert( consdata->tmpvars != NULL );
599  assert( consdata->tmpvals != NULL );
600  assert( consdata->weights != NULL );
601  assert( consdata->cases != NULL );
602 
603  *infeasible = FALSE;
604  *nfixedvars = 0;
605 
606  nspcons = consdata->nspcons;
607  nblocks = consdata->nblocks;
608  vars = consdata->vars;
609  vals = consdata->vals;
610  tmpvars = consdata->tmpvars;
611  tmpvals = consdata->tmpvals;
612  weights = consdata->weights;
613  cases = consdata->cases;
614 
615  /* check for upper right triangle */
616  if ( ! consdata->istrianglefixed )
617  {
618  SCIP_CALL( fixTriangle(scip, cons, infeasible, nfixedvars) );
619  if ( *infeasible )
620  return SCIP_OKAY;
621  if ( *nfixedvars > 0 )
622  return SCIP_OKAY;
623  }
624 
625  /* compute table if necessary (i.e., not computed before) */
626  computeSCTable(scip, nspcons, nblocks, weights, cases, vals);
627 
628  /* loop through rows */
629  for (i = 1; i < nspcons && ! (*infeasible); ++i)
630  {
631  SCIP_Real bar; /* value of bar: */
632  int lastcolumn; /* last column considered as part of the bar */
633 
634  bar = 0.0;
635  lastcolumn = nblocks - 1;
636  if ( lastcolumn > i )
637  lastcolumn = i;
638 
639  /* traverse row from right to left: */
640  /* j >= 2, since for j = 1 we look at column 0, which is uninteresting due to the one at position (0,0) */
641  for (j = lastcolumn; j > 1; --j)
642  {
643  bar += vals[i][j];
644 
645  /* check whether weights[i-1][j-1] < bar (<=> bar - weights[i-1][j-1] > 0), i.e. cut is violated) */
646  if ( SCIPisEfficacious(scip, bar - weights[i-1][j-1]) )
647  {
648  SCIP_Real weight;
649  SCIP_ROW* row;
650 #ifdef SCIP_DEBUG
651  char name[SCIP_MAXSTRLEN];
652 #endif
653  int nvars;
654  int p1;
655  int p2;
656 
657  nvars = 0;
658  p1 = i-1;
659  p2 = j-1;
660  weight = 0.0;
661 
662  /* first add bar */
663  for (l = j; l <= lastcolumn; ++l)
664  {
665  tmpvars[nvars] = vars[i][l];
666  tmpvals[nvars] = 1.0;
667  nvars++;
668  }
669 
670  /* then add shifted column */
671  do
672  {
673  assert( cases[p1][p2] != -1 );
674  assert( p1 >= 0 && p1 < i );
675  assert( p2 >= 0 && p2 < j );
676 
677  /* if case 1 */
678  if (cases[p1][p2] == 1)
679  p2--; /* decrease column */
680  else
681  {
682  /* case 2 or 3: */
683  assert( cases[p1][p2] == 2 || cases[p1][p2] == 3 );
684  tmpvars[nvars] = vars[p1][p2];
685  tmpvals[nvars] = -1.0;
686  nvars++;
687  weight += vals[p1][p2];
688  if ( cases[p1][p2] == 3 )
689  break;
690  }
691  p1--; /* decrease row */
692  }
693  while ( p1 >= 0 ); /* should always be true, i.e. the break should end the loop */
694  assert( cases[p1][p2] == 3 );
695 
696  /* generate cut */
697 #ifdef SCIP_DEBUG
698  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "sci_%d_%d", i, j);
699  SCIP_CALL( SCIPcreateEmptyRowCons(scip, &row, conshdlr, name, -SCIPinfinity(scip), 0.0, FALSE, FALSE, TRUE) );
700 #else
701  SCIP_CALL( SCIPcreateEmptyRowCons(scip, &row, conshdlr, "", -SCIPinfinity(scip), 0.0, FALSE, FALSE, TRUE) );
702 #endif
703  SCIP_CALL( SCIPaddVarsToRow(scip, row, nvars, tmpvars, tmpvals) );
704  /*SCIP_CALL( SCIPprintRow(scip, row, NULL) ); */
705  SCIP_CALL( SCIPaddCut(scip, NULL, row, FALSE, infeasible) );
706  SCIP_CALL( SCIPreleaseRow(scip, &row) );
707  ++(*ncuts);
708 
709 #ifdef SHOW_SCI
710  SCIP_CALL( printSCI(scip, nspcons, nblocks, cases, i, j) );
711 #endif
712 
713  assert( SCIPisSumEQ(scip, weights[i-1][j-1], weight) );
714  }
715  }
716  }
717  return SCIP_OKAY;
718 }
719 
720 
721 /** propagation method for a single orbitope constraint */
722 static
723 SCIP_RETCODE propagateCons(
724  SCIP* scip, /**< SCIP data structure */
725  SCIP_CONS* cons, /**< constraint to be processed */
726  SCIP_Bool* infeasible, /**< pointer to store TRUE, if the node can be cut off */
727  int* nfixedvars /**< pointer to add up the number of found domain reductions */
728  )
729 {
730  SCIP_CONSDATA* consdata;
731  SCIP_VAR*** vars;
732  SCIP_Bool ispart;
733  int* firstnonzeros;
734  int* lastones;
735  int* frontiersteps;
736  int lastoneprevrow;
737  int nspcons;
738  int nblocks;
739  int nsteps;
740  int i;
741  int j;
742 
743  assert( scip != NULL );
744  assert( cons != NULL );
745  assert( infeasible != NULL );
746  assert( nfixedvars != NULL );
747 
748  consdata = SCIPconsGetData(cons);
749  assert( consdata != NULL );
750  assert( consdata->nspcons > 0 );
751  assert( consdata->nblocks > 0 );
752  assert( consdata->vars != NULL );
753 
754  *nfixedvars = 0;
755 
756  nspcons = consdata->nspcons;
757  nblocks = consdata->nblocks;
758  vars = consdata->vars;
759  ispart = consdata->ispart;
760 
761  /* fix upper right triangle if still necessary */
762  if ( ! consdata->istrianglefixed )
763  {
764  int nfixed = 0;
765  SCIP_CALL( fixTriangle(scip, cons, infeasible, &nfixed) );
766  *nfixedvars += nfixed;
767  }
768 
769  /* prepare further propagation */
770  SCIP_CALL( SCIPallocBufferArray(scip, &firstnonzeros, nspcons) );
771  SCIP_CALL( SCIPallocBufferArray(scip, &lastones, nspcons) );
772  SCIP_CALL( SCIPallocBufferArray(scip, &frontiersteps, nblocks) );
773 
774 #ifdef PRINT_MATRIX
775  SCIPdebugMessage("Matrix:\n");
776  printMatrix(scip, consdata);
777 #endif
778 
779  /* propagate */
780  lastoneprevrow = 0;
781  lastones[0] = 0;
782 
783  if ( ! ispart )
784  {
785  /* packing case: if entry (0,0) is fixed to 0 */
786  if ( SCIPvarGetUbLocal(vars[0][0]) < 0.5 )
787  {
788  lastoneprevrow = -1;
789  lastones[0] = -1;
790  }
791  }
792  nsteps = 0;
793 
794  for (i = 1; i < nspcons; ++i)
795  {
796  int lastcolumn;
797  int firstnonzeroinrow;
798  int lastoneinrow;
799  SCIP_Bool infrontier;
800 
801  /* last column considered as part of the bar: */
802  lastcolumn = nblocks - 1;
803  if ( lastcolumn > i )
804  lastcolumn = i;
805 
806  /* find first position not fixed to 0 (partitioning) or fixed to 1 (packing) */
807  firstnonzeroinrow = -1;
808  for (j = 0; j <= lastcolumn; ++j)
809  {
810  if ( ispart )
811  {
812  /* partitioning case: if variable is not fixed to 0 */
813  if ( SCIPvarGetUbLocal(vars[i][j]) > 0.5 )
814  {
815  firstnonzeroinrow = j;
816  break;
817  }
818  }
819  else
820  {
821  /* packing case: if variable is fixed to 1 */
822  if ( SCIPvarGetLbLocal(vars[i][j]) > 0.5 )
823  {
824  firstnonzeroinrow = j;
825  break;
826  }
827  }
828  }
829  /* if all variables are fixed to 0 in the partitioning case - should not happen */
830  if ( firstnonzeroinrow == -1 && ispart )
831  {
832  SCIPdebugMessage(" -> Infeasible node: all variables in row %d are fixed to 0.\n", i);
833  *infeasible = TRUE;
834  /* conflict should be analyzed by setppc constraint handler */
835  goto TERMINATE;
836  }
837  firstnonzeros[i] = firstnonzeroinrow;
838  assert( !ispart || firstnonzeroinrow >= 0 );
839  assert( -1 <= firstnonzeroinrow && firstnonzeroinrow <= lastcolumn );
840 
841  /* compute rightmost possible position for a 1 */
842  lastoneinrow = -1;
843  assert( !ispart || 0 <= lastoneprevrow );
844  assert( lastoneprevrow <= lastcolumn );
845 
846  /* if we are at right border or if entry in column lastoneprevrow+1 is fixed to 0 */
847  infrontier = FALSE;
848  if ( lastoneprevrow == nblocks-1 || SCIPvarGetUbLocal(vars[i][lastoneprevrow+1]) < 0.5 )
849  lastoneinrow = lastoneprevrow;
850  else
851  {
852  lastoneinrow = lastoneprevrow + 1;
853  frontiersteps[nsteps++] = i;
854  infrontier = TRUE;
855  }
856 
857  /* store lastoneinrow */
858  assert( !ispart || 0 <= lastoneinrow );
859  assert( lastoneinrow <= lastcolumn );
860  lastones[i] = lastoneinrow;
861 
862  /* check whether we are infeasible */
863  if ( firstnonzeroinrow > lastoneinrow )
864  {
865  int k;
866 
867 #ifdef SCIP_DEBUG
868  if ( ispart )
869  {
870  SCIPdebugMessage(" -> Infeasible node: row %d, leftmost nonzero at %d, rightmost 1 at %d\n",
871  i, firstnonzeroinrow, lastoneinrow);
872  }
873  else
874  {
875  SCIPdebugMessage(" -> Infeasible node: row %d, 1 at %d, rightmost position for 1 at %d\n",
876  i, firstnonzeroinrow, lastoneinrow);
877  }
878 #endif
879  /* check if conflict analysis is applicable */
880  if ( SCIPisConflictAnalysisApplicable(scip) )
881  {
882  /* conflict analysis only applicable in SOLVING stage */
883  assert( SCIPgetStage(scip) == SCIP_STAGE_SOLVING || SCIPinProbing(scip) );
884 
885  /* perform conflict analysis */
886  SCIP_CALL( SCIPinitConflictAnalysis(scip) );
887 
888  if ( ispart )
889  {
890  /* add bounds (variables fixed to 0) that result in the first nonzero entry */
891  for (j = 0; j <= lastcolumn; ++j)
892  {
893  /* add varaibles in row up to the first variable fixed to 0 */
894  if ( SCIPvarGetUbLocal(vars[i][j]) > 0.5 )
895  break;
896 
897  assert( SCIPvarGetUbLocal(vars[i][j]) < 0.5 );
898  SCIP_CALL( SCIPaddConflictBinvar(scip, vars[i][j]) );
899  }
900  }
901  else
902  {
903  /* add bounds that result in the last one - check top left entry for packing case */
904  if ( lastones[0] == -1 )
905  {
906  assert( SCIPvarGetUbLocal(vars[0][0]) < 0.5 );
907  SCIP_CALL( SCIPaddConflictBinvar(scip, vars[0][0]) );
908  }
909 
910  /* mark variable fixed to 1 */
911  assert( SCIPvarGetLbLocal(vars[i][firstnonzeroinrow]) > 0.5 );
912  SCIP_CALL( SCIPaddConflictBinvar(scip, vars[i][firstnonzeroinrow]) );
913  }
914 
915  /* add bounds that result in the last one - pass through rows */
916  for (k = 1; k < i; ++k)
917  {
918  int l;
919  l = lastones[k] + 1;
920 
921  /* if the frontier has not moved and we are not beyond the matrix boundaries */
922  if ( l <= nblocks-1 && l <= k && lastones[k-1] == lastones[k] )
923  {
924  assert( SCIPvarGetUbLocal(vars[k][l]) < 0.5 );
925  SCIP_CALL( SCIPaddConflictBinvar(scip, vars[k][l]) );
926  }
927  }
928  SCIP_CALL( SCIPanalyzeConflictCons(scip, cons, NULL) );
929  }
930 
931  *infeasible = TRUE;
932  goto TERMINATE;
933  }
934 
935  /* fix entries beyond the last possible position for a 1 in the row to 0 (see Lemma 1 in the paper) */
936  for (j = lastoneinrow+1; j <= lastcolumn; ++j)
937  {
938  /* if the entry is not yet fixed to 0 */
939  if ( SCIPvarGetUbLocal(vars[i][j]) > 0.5 )
940  {
941  SCIP_Bool tightened;
942  int inferInfo;
943 
944  SCIPdebugMessage(" -> Fixing entry (%d,%d) to 0.\n", i, j);
945 
946  tightened = FALSE;
947 
948  /* fix variable to 0 and store position of (i,lastoneinrow+1) for conflict resolution */
949  inferInfo = i * nblocks + lastoneinrow + 1;
950  /* correction according to Lemma 1 in the paper (second part): store (i,lastoneinrow+2) */
951  if ( !infrontier )
952  ++inferInfo;
953  SCIP_CALL( SCIPinferBinvarCons(scip, vars[i][j], FALSE, cons, inferInfo, infeasible, &tightened) );
954 
955  /* if entry is fixed to one -> infeasible node */
956  if ( *infeasible )
957  {
958  SCIPdebugMessage(" -> Infeasible node: row %d, 1 in column %d beyond rightmost position %d\n", i, j, lastoneinrow);
959  /* check if conflict analysis is applicable */
960  if( SCIPisConflictAnalysisApplicable(scip) )
961  {
962  int k;
963 
964  /* conflict analysis only applicable in SOLVING stage */
965  assert(SCIPgetStage(scip) == SCIP_STAGE_SOLVING || SCIPinProbing(scip));
966 
967  /* perform conflict analysis */
968  SCIP_CALL( SCIPinitConflictAnalysis(scip) );
969 
970  /* add current bound */
971  SCIP_CALL( SCIPaddConflictBinvar(scip, vars[i][j]) );
972 
973  /* add bounds that result in the last one - check top left entry for packing case */
974  if ( ! ispart && lastones[0] == -1 )
975  {
976  assert( SCIPvarGetUbLocal(vars[0][0]) < 0.5 );
977  SCIP_CALL( SCIPaddConflictBinvar(scip, vars[0][0]) );
978  }
979 
980  /* add bounds that result in the last one - pass through rows */
981  for (k = 1; k < i; ++k)
982  {
983  int l;
984  l = lastones[k] + 1;
985 
986  /* if the frontier has not moved and we are not beyond the matrix boundaries */
987  if ( l <= nblocks-1 && l <= k && lastones[k-1] == lastones[k] )
988  {
989  assert( SCIPvarGetUbLocal(vars[k][l]) < 0.5 );
990  SCIP_CALL( SCIPaddConflictBinvar(scip, vars[k][l]) );
991  }
992  }
993  SCIP_CALL( SCIPanalyzeConflictCons(scip, cons, NULL) );
994  }
995 
996  goto TERMINATE;
997  }
998  if ( tightened )
999  ++(*nfixedvars);
1000  }
1001  }
1002 
1003  lastoneprevrow = lastoneinrow;
1004  }
1005 
1006  /* check whether fixing any entry to 0 results in a contradiction -> loop through rows in frontiersteps (a.k.a. gamma) */
1007  for (j = 0; j < nsteps; ++j)
1008  {
1009  int s;
1010  int lastoneinrow;
1011 
1012  s = frontiersteps[j];
1013  lastoneinrow = lastones[s];
1014  /* note for packing case: if we are in a frontier step then lastoneinrow >= 0 */
1015  assert( 0 <= lastoneinrow && lastoneinrow < nblocks );
1016 
1017  /* if entry is not fixed */
1018  if ( SCIPvarGetLbLocal(vars[s][lastoneinrow]) < 0.5 && SCIPvarGetUbLocal(vars[s][lastoneinrow]) > 0.5 )
1019  {
1020  int betaprev;
1021  betaprev = lastoneinrow - 1;
1022 
1023  /* loop through rows below s */
1024  for (i = s+1; i < nspcons; ++i)
1025  {
1026  int beta;
1027  beta = -2;
1028 
1029  if ( betaprev == nblocks-1 || SCIPvarGetUbLocal(vars[i][betaprev+1]) < 0.5 )
1030  beta = betaprev;
1031  else
1032  beta = betaprev + 1;
1033  assert( -1 <= beta && beta < nblocks );
1034 
1035  if ( firstnonzeros[i] > beta )
1036  {
1037  SCIP_Bool tightened;
1038  int inferInfo;
1039 
1040  /* can fix (s,lastoneinrow) (a.k.a (s,alpha)) to 1
1041  * (do not need to fix other entries to 0, since they will be
1042  * automatically fixed by SCIPtightenVarLb.)
1043  */
1044  assert( SCIPvarGetLbLocal(vars[s][lastoneinrow]) < 0.5 );
1045  SCIPdebugMessage(" -> Fixing entry (%d,%d) to 1.\n", s, lastoneinrow);
1046 
1047  tightened = FALSE;
1048 
1049  /* store position (i,firstnonzeros[i]) */
1050  inferInfo = nblocks * nspcons + i * nblocks + firstnonzeros[i];
1051  SCIP_CALL( SCIPinferBinvarCons(scip, vars[s][lastoneinrow], TRUE, cons, inferInfo, infeasible, &tightened) );
1052 
1053  assert( !(*infeasible) );
1054  if ( tightened )
1055  ++(*nfixedvars);
1056  break;
1057  }
1058  betaprev = beta;
1059  }
1060  }
1061  }
1062 
1063  TERMINATE:
1064  SCIPfreeBufferArray(scip, &frontiersteps);
1065  SCIPfreeBufferArray(scip, &lastones);
1066  SCIPfreeBufferArray(scip, &firstnonzeros);
1067 
1068  return SCIP_OKAY;
1069 }
1070 
1071 
1072 /** Propagation conflict resolving method of propagator
1073  *
1074  * In this function we use that the propagation method above implicitly propagates SCIs, i.e., every
1075  * fixing can also be gotten via an SCI-fixing.
1076  *
1077  * Since the storage of an integer is not enough to store the complete information about the fixing
1078  * nor a complete shifted column, we have to use the linear time algorithm for SCIs.
1079  *
1080  * The inferinfo integer is set as follows:
1081  *
1082  * - If a shifted column is fixed to 0 and the corresponding bar does not necessarily has value 1
1083  * then we fix these entries to 0 and inferinfo is i * nblocks + j, where (i,j) is the leader of the
1084  * bar. The SCI depends on whether i is in Gamma or not (see Lemma 1 in the paper and the comments
1085  * above).
1086  *
1087  * - If a bar has value 1 and the shifted column has one entry that is not fixed, it can be fixed to
1088  * 1 and inferinfo is (nspcons*nblocks) + i * nblocks + j, where (i,j) is the leader of the bar; see
1089  * Proposition 1 (2c).
1090  */
1091 static
1092 SCIP_RETCODE resolvePropagation(
1093  SCIP* scip, /**< SCIP data structure */
1094  SCIP_CONS* cons, /**< constraint that inferred the bound change */
1095  SCIP_VAR* infervar, /**< variable that was deduced */
1096  int inferinfo, /**< inference information */
1097  SCIP_BOUNDTYPE boundtype, /**< the type of the changed bound (lower or upper bound) */
1098  SCIP_BDCHGIDX* bdchgidx, /**< bound change index (time stamp of bound change), or NULL for current time */
1099  SCIP_RESULT* result /**< pointer to store the result of the propagation conflict resolving call */
1100  )
1101 { /*lint --e{715}*/
1102  SCIP_CONSDATA* consdata;
1103  SCIP_Real** vals;
1104  SCIP_Real** weights;
1105  SCIP_VAR*** vars;
1106  SCIP_Bool ispart;
1107  int** cases;
1108 
1109  int i;
1110  int j;
1111  int nspcons;
1112  int nblocks;
1113 
1114  assert( scip != NULL );
1115  assert( cons != NULL );
1116  assert( result != NULL );
1117 
1118  consdata = SCIPconsGetData(cons);
1119  assert( consdata != NULL );
1120  assert( consdata->nspcons > 0 );
1121  assert( consdata->nblocks > 0 );
1122  assert( consdata->vars != NULL );
1123  assert( consdata->vals != NULL );
1124  assert( consdata->weights != NULL );
1125  assert( consdata->cases != NULL );
1126  assert( consdata->istrianglefixed );
1127 
1128  *result = SCIP_DIDNOTFIND;
1129  if ( ! consdata->resolveprop )
1130  return SCIP_OKAY;
1131 
1132  nspcons = consdata->nspcons;
1133  nblocks = consdata->nblocks;
1134  vars = consdata->vars;
1135  vals = consdata->vals;
1136  weights = consdata->weights;
1137  ispart = consdata->ispart;
1138  cases = consdata->cases;
1139 
1140  SCIPdebugMessage("Propagation resolution method of orbitope constraint using orbitopal fixing\n");
1141 
1142  /* fill table */
1143  for (i = 0; i < nspcons; ++i)
1144  {
1145  int lastcolumn;
1146 
1147  /* last column considered as part of the bar: */
1148  lastcolumn = nblocks - 1;
1149  if ( lastcolumn > i )
1150  lastcolumn = i;
1151  for (j = 0; j <= lastcolumn; ++j)
1152  {
1153  /* if the variable was fixed to zero at conflict time */
1154  if ( SCIPvarGetUbAtIndex(vars[i][j], bdchgidx, FALSE) < 0.5 )
1155  vals[i][j] = 0.0;
1156  else
1157  {
1158  /* if the variable was fixed to one at conflict time */
1159  if ( SCIPvarGetLbAtIndex(vars[i][j], bdchgidx, FALSE) > 0.5 )
1160  vals[i][j] = 2.0;
1161  else
1162  vals[i][j] = 1.0;
1163  }
1164  }
1165  }
1166 
1167 #ifdef PRINT_MATRIX
1168  SCIPdebugMessage("Matrix:\n");
1169  printMatrix(scip, consdata);
1170 #endif
1171 
1172  /* computation of table: this now minimizes the value of the shifted column */
1173  assert( consdata->istrianglefixed );
1174  computeSCTable(scip, nspcons, nblocks, weights, cases, vals);
1175 
1176  /* if we fixed variables in the bar to zero */
1177  assert( inferinfo >= 0 && inferinfo < 2 * nspcons * nblocks );
1178  if ( inferinfo < nspcons * nblocks )
1179  {
1180  int p1;
1181  int p2;
1182 #ifdef SCIP_DEBUG
1183  char str[SCIP_MAXSTRLEN];
1184  char tmpstr[SCIP_MAXSTRLEN];
1185 #endif
1186 
1187  i = (int) (inferinfo / nblocks);
1188  j = inferinfo % nblocks;
1189  assert( 0 <= i && i < nspcons );
1190  assert( 0 <= j && j < nblocks );
1191 
1192  /* find SCI with value 0 */
1193  assert( weights[i-1][j-1] < 0.5 );
1194 
1195  SCIPdebugMessage(" -> reason for x[%d][%d] = ... = x[%d][%d] = 0 was the following SC:\n", i, j, i, MIN(i,nblocks));
1196 #ifdef SCIP_DEBUG
1197  str[0] = '\0';
1198 #endif
1199 
1200  p1 = i-1;
1201  p2 = j-1;
1202  do
1203  {
1204  assert( cases[p1][p2] != -1 );
1205  assert( p1 >= 0 && p1 < i );
1206  assert( p2 >= 0 && p2 < j );
1207 
1208  /* if case 1 */
1209  if ( cases[p1][p2] == 1 )
1210  --p2; /* decrease column */
1211  else
1212  {
1213  /* case 2 or 3: */
1214  assert( cases[p1][p2] == 2 || cases[p1][p2] == 3 );
1215  assert( SCIPvarGetUbAtIndex(vars[p1][p2], bdchgidx, FALSE) < 0.5 );
1216  SCIP_CALL( SCIPaddConflictUb(scip, vars[p1][p2], bdchgidx) );
1217  *result = SCIP_SUCCESS;
1218 
1219 #ifdef SCIP_DEBUG
1220  (void) SCIPsnprintf(tmpstr, SCIP_MAXSTRLEN, " (%d,%d)", p1, p2);
1221  (void) strncat(str, tmpstr, SCIP_MAXSTRLEN);
1222 #endif
1223 
1224  if ( cases[p1][p2] == 3 )
1225  break;
1226  }
1227  --p1; /* decrease row */
1228  }
1229  while ( p1 >= 0 ); /* should always be true, i.e. the break should end the loop */
1230  assert( cases[p1][p2] == 3 );
1231 
1232 #ifdef SCIP_DEBUG
1233  SCIPdebugMessage("%s\n", str);
1234 #endif
1235  }
1236  else
1237  {
1238  int k;
1239  int p1;
1240  int p2;
1241 #ifndef NDEBUG
1242  int pos1;
1243  int pos2;
1244 #endif
1245 #ifdef SCIP_DEBUG
1246  char str[SCIP_MAXSTRLEN];
1247  char tmpstr[SCIP_MAXSTRLEN];
1248 #endif
1249 
1250  /* if we fixed a variable in the SC to 1 */
1251  inferinfo -= nspcons * nblocks;
1252  i = (int) inferinfo / nblocks;
1253  j = inferinfo % nblocks;
1254  assert( 0 <= i && i < nspcons );
1255  assert( 0 <= j && j < nblocks );
1256 
1257  /* In rare cases it might happen that we fixed a variable to 1, but the node later becomes infeasible by globally
1258  * fixing variables to 0. In this case, it might happen that we find a SC with value 0 instead of 1. We then
1259  * cannot use this SC to repropagate (and do not know how to reconstruct the original reasoning). */
1260  if ( weights[i-1][j-1] > 0.5 && weights[i-1][j-1] < 1.5 )
1261  {
1262  SCIPdebugMessage(" -> reason for x[%d][%d] = 1 was the following SC:\n", i, j);
1263 #ifdef SCIP_DEBUG
1264  (void) SCIPsnprintf(str, SCIP_MAXSTRLEN, "SC:");
1265 #endif
1266 
1267  p1 = i-1;
1268  p2 = j-1;
1269 #ifndef NDEBUG
1270  pos1 = -1;
1271  pos2 = -1;
1272 #endif
1273  do
1274  {
1275  assert( cases[p1][p2] != -1 );
1276  assert( p1 >= 0 && p1 < i );
1277  assert( p2 >= 0 && p2 < j );
1278 
1279  /* if case 1 */
1280  if ( cases[p1][p2] == 1 )
1281  --p2; /* decrease column */
1282  else
1283  {
1284  /* case 2 or 3: reason are formed by variables in SC fixed to 0 */
1285  assert( cases[p1][p2] == 2 || cases[p1][p2] == 3 );
1286  if ( SCIPvarGetUbAtIndex(vars[p1][p2], bdchgidx, FALSE) < 0.5 )
1287  {
1288  SCIP_CALL( SCIPaddConflictUb(scip, vars[p1][p2], bdchgidx) );
1289  *result = SCIP_SUCCESS;
1290 
1291 #ifdef SCIP_DEBUG
1292  (void) SCIPsnprintf(tmpstr, SCIP_MAXSTRLEN, " (%d,%d)", p1, p2);
1293  (void) strncat(str, tmpstr, SCIP_MAXSTRLEN);
1294 #endif
1295  }
1296 #ifndef NDEBUG
1297  else
1298  {
1299  assert( SCIPvarGetLbAtIndex(vars[p1][p2], bdchgidx, FALSE) < 0.5 );
1300  assert( pos1 == -1 && pos2 == -1 );
1301  pos1 = p1;
1302  pos2 = p2;
1303  }
1304 #endif
1305  if ( cases[p1][p2] == 3 )
1306  break;
1307  }
1308  --p1; /* decrease row */
1309  }
1310  while ( p1 >= 0 ); /* should always be true, i.e., the break should end the loop */
1311  assert( cases[p1][p2] == 3 );
1312  assert( pos1 >= 0 && pos2 >= 0 );
1313 
1314  /* distinguish partitioning/packing */
1315  if ( ispart )
1316  {
1317  /* partitioning case */
1318 #ifdef SCIP_DEBUG
1319  (void) SCIPsnprintf(tmpstr, SCIP_MAXSTRLEN, " before bar: ");
1320  (void) strncat(str, tmpstr, SCIP_MAXSTRLEN);
1321 #endif
1322  /* add variables before the bar in the partitioning case */
1323  for (k = 0; k < j; ++k)
1324  {
1325  assert( SCIPvarGetUbAtIndex(vars[i][k], bdchgidx, FALSE) < 0.5 );
1326  SCIP_CALL( SCIPaddConflictUb(scip, vars[i][k], bdchgidx) );
1327  *result = SCIP_SUCCESS;
1328 #ifdef SCIP_DEBUG
1329  (void) SCIPsnprintf(tmpstr, SCIP_MAXSTRLEN, " (%d,%d)", i, k);
1330  (void) strncat(str, tmpstr, SCIP_MAXSTRLEN);
1331 #endif
1332  }
1333 
1334 #ifdef SCIP_DEBUG
1335  SCIPdebugMessage("%s\n", str);
1336 #endif
1337  }
1338  else
1339  {
1340  /* packing case */
1341  int lastcolumn;
1342 
1343  /* last column considered as part of the bar: */
1344  lastcolumn = nblocks - 1;
1345  if ( lastcolumn > i )
1346  lastcolumn = i;
1347 
1348  /* search for variable in the bar that is fixed to 1 in the packing case */
1349  for (k = j; k <= lastcolumn; ++k)
1350  {
1351  if ( SCIPvarGetLbAtIndex(vars[i][k], bdchgidx, FALSE) > 0.5 )
1352  {
1353  SCIP_CALL( SCIPaddConflictLb(scip, vars[i][k], bdchgidx) );
1354  *result = SCIP_SUCCESS;
1355  SCIPdebugMessage(" and variable x[%d][%d] fixed to 1.\n", i, k);
1356  break;
1357  }
1358  }
1359  }
1360  }
1361  }
1362 
1363  return SCIP_OKAY;
1364 }
1365 
1366 
1367 /*
1368  * Callback methods of constraint handler
1369  */
1370 
1371 /** copy method for constraint handler plugins (called when SCIP copies plugins) */
1372 static
1373 SCIP_DECL_CONSHDLRCOPY(conshdlrCopyOrbitope)
1374 { /*lint --e{715}*/
1375  assert(scip != NULL);
1376  assert(conshdlr != NULL);
1377  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
1378 
1379  /* call inclusion method of constraint handler */
1380  SCIP_CALL( SCIPincludeConshdlrOrbitope(scip) );
1381 
1382  *valid = TRUE;
1383 
1384  return SCIP_OKAY;
1385 }
1386 
1387 /** frees specific constraint data */
1388 static
1389 SCIP_DECL_CONSDELETE(consDeleteOrbitope)
1390 { /*lint --e{715}*/
1391  assert(conshdlr != NULL);
1392  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
1393 
1394  SCIP_CALL( consdataFree(scip, consdata) );
1395 
1396  return SCIP_OKAY;
1397 }
1398 
1399 /** transforms constraint data into data belonging to the transformed problem */
1400 static
1401 SCIP_DECL_CONSTRANS(consTransOrbitope)
1402 { /*lint --e{715}*/
1403  SCIP_CONSDATA* sourcedata;
1404  SCIP_CONSDATA* targetdata;
1405 
1406  assert(conshdlr != NULL);
1407  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
1408  assert(SCIPgetStage(scip) == SCIP_STAGE_TRANSFORMING);
1409  assert(sourcecons != NULL);
1410  assert(targetcons != NULL);
1411 
1412  sourcedata = SCIPconsGetData(sourcecons);
1413  assert(sourcedata != NULL);
1414 
1415  /* create linear constraint data for target constraint */
1416  SCIP_CALL( consdataCreate(scip, &targetdata, sourcedata->vars, sourcedata->nspcons, sourcedata->nblocks,
1417  sourcedata->ispart, sourcedata->resolveprop) );
1418 
1419  /* create target constraint */
1420  SCIP_CALL( SCIPcreateCons(scip, targetcons, SCIPconsGetName(sourcecons), conshdlr, targetdata,
1421  SCIPconsIsInitial(sourcecons), SCIPconsIsSeparated(sourcecons), SCIPconsIsEnforced(sourcecons),
1422  SCIPconsIsChecked(sourcecons), SCIPconsIsPropagated(sourcecons),
1423  SCIPconsIsLocal(sourcecons), SCIPconsIsModifiable(sourcecons),
1424  SCIPconsIsDynamic(sourcecons), SCIPconsIsRemovable(sourcecons), SCIPconsIsStickingAtNode(sourcecons)) );
1425 
1426  return SCIP_OKAY;
1427 }
1428 
1429 /** separation method of constraint handler for LP solutions */
1430 static
1431 SCIP_DECL_CONSSEPALP(consSepalpOrbitope)
1432 { /*lint --e{715}*/
1433  assert( scip != NULL );
1434  assert( conshdlr != NULL );
1435  assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
1436  assert( result != NULL );
1437 
1438  *result = SCIP_DIDNOTRUN;
1439 
1440  /* if solution is not integer */
1441  if ( SCIPgetNLPBranchCands(scip) > 0 )
1442  {
1443  SCIP_Bool infeasible;
1444  int nfixedvars = 0;
1445  int ncuts = 0;
1446  int c;
1447 
1448  *result = SCIP_DIDNOTFIND;
1449  infeasible = FALSE;
1450 
1451  /* loop through constraints */
1452  for (c = 0; c < nusefulconss && ! infeasible; c++)
1453  {
1454  SCIP_CONSDATA* consdata;
1455 
1456  assert( conss[c] != NULL );
1457 
1458  /* get data of constraint */
1459  consdata = SCIPconsGetData(conss[c]);
1460  assert( consdata != NULL );
1461 
1462  /* get solution */
1463  copyValues(scip, consdata, NULL);
1464  SCIPdebugMessage("Separating SCIs for orbitope constraint <%s> ...\n", SCIPconsGetName(conss[c]));
1465 
1466  /* separate */
1467  SCIP_CALL( separateSCIs(scip, conshdlr, conss[c], consdata, &infeasible, &nfixedvars, &ncuts) );
1468  }
1469  if ( infeasible )
1470  {
1471  SCIPdebugMessage("Infeasible node.\n");
1472  *result = SCIP_CUTOFF;
1473  }
1474  else if ( nfixedvars > 0 )
1475  {
1476  SCIPdebugMessage("Fixed %d variables.\n", nfixedvars);
1477  *result = SCIP_REDUCEDDOM;
1478  }
1479  else if ( ncuts > 0 )
1480  {
1481  SCIPdebugMessage("Separated %d SCIs.\n", ncuts);
1482  *result = SCIP_SEPARATED;
1483  }
1484  else
1485  {
1486  SCIPdebugMessage("No violated SCI found.\n");
1487  }
1488  }
1489 
1490  return SCIP_OKAY;
1491 }
1492 
1493 /** separation method of constraint handler for arbitrary primal solutions */
1494 static
1495 SCIP_DECL_CONSSEPASOL(consSepasolOrbitope)
1496 { /*lint --e{715}*/
1497  SCIP_Bool infeasible = FALSE;
1498  int nfixedvars = 0;
1499  int ncuts = 0;
1500  int c;
1501 
1502  assert( scip != NULL );
1503  assert( conshdlr != NULL );
1504  assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
1505  assert( result != NULL );
1506 
1507  *result = SCIP_DIDNOTFIND;
1508 
1509  /* loop through constraints */
1510  for (c = 0; c < nusefulconss && ! infeasible; c++)
1511  {
1512  SCIP_CONSDATA* consdata;
1513 
1514  /* get data of constraint */
1515  assert( conss[c] != NULL );
1516  consdata = SCIPconsGetData(conss[c]);
1517  assert( consdata != NULL );
1518 
1519  /* get solution */
1520  copyValues(scip, consdata, sol);
1521  SCIPdebugMessage("Separating SCIs (solution) for orbitope constraint <%s> \n", SCIPconsGetName(conss[c]));
1522 
1523  /* separate */
1524  SCIP_CALL( separateSCIs(scip, conshdlr, conss[c], consdata, &infeasible, &nfixedvars, &ncuts) );
1525  }
1526 
1527  if ( infeasible )
1528  {
1529  SCIPdebugMessage("Infeasible node.\n");
1530  *result = SCIP_CUTOFF;
1531  }
1532  else if ( nfixedvars > 0 )
1533  {
1534  SCIPdebugMessage("Fixed %d variables.\n", nfixedvars);
1535  *result = SCIP_REDUCEDDOM;
1536  }
1537  else if ( ncuts > 0 )
1538  {
1539  SCIPdebugMessage("Separated %d SCIs.\n", ncuts);
1540  *result = SCIP_SEPARATED;
1541  }
1542  else
1543  {
1544  SCIPdebugMessage("No violated SCI found.\n");
1545  }
1546 
1547  return SCIP_OKAY;
1548 }
1549 
1550 
1551 /** constraint enforcing method of constraint handler for LP solutions */
1552 static
1553 SCIP_DECL_CONSENFOLP(consEnfolpOrbitope)
1554 { /*lint --e{715}*/
1555  SCIP_Bool infeasible = FALSE;
1556  int nfixedvars = 0;
1557  int ncuts = 0;
1558  int c;
1559 
1560  assert( scip != NULL );
1561  assert( conshdlr != NULL );
1562  assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
1563  assert( result != NULL );
1564 
1565  *result = SCIP_FEASIBLE;
1566 
1567  /* we have a negative priority, so we should come after the integrality conshdlr */
1568  assert( SCIPgetNLPBranchCands(scip) == 0 );
1569 
1570  /* loop through constraints */
1571  for (c = 0; c < nusefulconss && ! infeasible; c++)
1572  {
1573  SCIP_CONSDATA* consdata;
1574 
1575  assert( conss[c] != NULL );
1576 
1577  /* get data of constraint */
1578  consdata = SCIPconsGetData(conss[c]);
1579  assert( consdata != NULL );
1580 
1581  /* get solution */
1582  copyValues(scip, consdata, NULL);
1583  SCIPdebugMessage("Enforcing for orbitope constraint <%s>\n", SCIPconsGetName(conss[c]));
1584 
1585  /* separate */
1586  SCIP_CALL( separateSCIs(scip, conshdlr, conss[c], consdata, &infeasible, &nfixedvars, &ncuts) );
1587  }
1588 
1589  if ( infeasible )
1590  {
1591  SCIPdebugMessage("Infeasible node.\n");
1592  *result = SCIP_CUTOFF;
1593  }
1594  else if ( nfixedvars > 0 )
1595  {
1596  SCIPdebugMessage("Fixed %d variables.\n", nfixedvars);
1597  *result = SCIP_REDUCEDDOM;
1598  }
1599  else if ( ncuts > 0 )
1600  {
1601  SCIPdebugMessage("Separated %d SCIs during enforcement.\n", ncuts);
1602  *result = SCIP_SEPARATED;
1603  }
1604  else
1605  {
1606  SCIPdebugMessage("No violated SCI found during enforcement.\n");
1607  }
1608 
1609  return SCIP_OKAY;
1610 }
1611 
1612 
1613 /** constraint enforcing method of constraint handler for pseudo solutions */
1614 static
1615 SCIP_DECL_CONSENFOPS(consEnfopsOrbitope)
1616 { /*lint --e{715}*/
1617  int c;
1618 
1619  assert( scip != NULL );
1620  assert( conshdlr != NULL );
1621  assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
1622  assert( result != NULL );
1623 
1624  *result = SCIP_FEASIBLE;
1625  if ( objinfeasible || solinfeasible )
1626  return SCIP_OKAY;
1627 
1628  /* loop through constraints */
1629  for (c = 0; c < nconss; ++c)
1630  {
1631  SCIP_CONSDATA* consdata;
1632  SCIP_Real** weights;
1633  SCIP_Real** vals;
1634  SCIP_CONS* cons;
1635  int** cases;
1636  int nspcons;
1637  int nblocks;
1638  int i;
1639  int j;
1640 
1641  /* get data of constraint */
1642  cons = conss[c];
1643  assert( cons != 0 );
1644  consdata = SCIPconsGetData(cons);
1645 
1646  assert( consdata != NULL );
1647  assert( consdata->nspcons > 0 );
1648  assert( consdata->nblocks > 0 );
1649  assert( consdata->vals != NULL );
1650  assert( consdata->weights != NULL );
1651  assert( consdata->cases != NULL );
1652 
1653  /* check for upper right triangle */
1654  if ( ! consdata->istrianglefixed )
1655  {
1656  SCIP_Bool infeasible = FALSE;
1657  int nfixedvars = 0;
1658 
1659  SCIP_CALL( fixTriangle(scip, cons, &infeasible, &nfixedvars) );
1660  if ( infeasible )
1661  {
1662  *result = SCIP_CUTOFF;
1663  return SCIP_OKAY;
1664  }
1665  if ( nfixedvars > 0 )
1666  {
1667  *result = SCIP_REDUCEDDOM;
1668  return SCIP_OKAY;
1669  }
1670  }
1671 
1672  nspcons = consdata->nspcons;
1673  nblocks = consdata->nblocks;
1674  vals = consdata->vals;
1675  weights = consdata->weights;
1676  cases = consdata->cases;
1677 
1678  /* get solution */
1679  copyValues(scip, consdata, NULL);
1680  SCIPdebugMessage("Enforcing (pseudo solutions) for orbitope constraint <%s>\n", SCIPconsGetName(conss[c]));
1681 
1682  /* compute table */
1683  assert( consdata->istrianglefixed );
1684  computeSCTable(scip, nspcons, nblocks, weights, cases, vals);
1685 
1686  /* loop through rows */
1687  for (i = 1; i < nspcons; ++i)
1688  {
1689  SCIP_Real bar = 0.0;
1690  int lastcolumn;
1691 
1692  lastcolumn = nblocks - 1;
1693 
1694  /* last column considered as part of the bar: */
1695  if ( lastcolumn > i )
1696  lastcolumn = i;
1697 
1698  /* traverse row from right to left */
1699  for (j = lastcolumn; j > 1; --j)
1700  {
1701  bar += vals[i][j];
1702  assert( SCIPisIntegral(scip, vals[i][j]) );
1703 
1704  /* check whether weights[i-1][j-1] < bar (<=> bar - weights[i-1][j-1] > 0), i.e. cut is violated) */
1705  if ( SCIPisGT(scip, bar - weights[i-1][j-1], 0.0) )
1706  {
1707  SCIPdebugMessage("Solution is infeasible.\n");
1708  *result = SCIP_INFEASIBLE;
1709  return SCIP_OKAY;
1710  }
1711  }
1712  }
1713  }
1714 
1715  return SCIP_OKAY;
1716 }
1717 
1718 
1719 /** feasibility check method of constraint handler for integral solutions */
1720 static
1721 SCIP_DECL_CONSCHECK(consCheckOrbitope)
1722 { /*lint --e{715}*/
1723  int c;
1724 
1725  assert( scip != NULL );
1726  assert( conshdlr != NULL );
1727  assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
1728  assert( result != NULL );
1729 
1730  *result = SCIP_FEASIBLE;
1731 
1732  /* loop through constraints */
1733  for (c = 0; c < nconss; ++c)
1734  {
1735  SCIP_CONSDATA* consdata;
1736  SCIP_VAR*** vars;
1737  SCIP_Real** vals;
1738  SCIP_Real** weights;
1739  int** cases;
1740  int nspcons;
1741  int nblocks;
1742  int i;
1743  int j;
1744 
1745  /* get data of constraint */
1746  assert( conss[c] != 0 );
1747  consdata = SCIPconsGetData(conss[c]);
1748 
1749  assert( consdata != NULL );
1750  assert( consdata->nspcons > 0 );
1751  assert( consdata->nblocks > 0 );
1752  assert( consdata->vars != NULL );
1753  assert( consdata->vals != NULL );
1754  assert( consdata->weights != NULL );
1755  assert( consdata->cases != NULL );
1756 
1757  nspcons = consdata->nspcons;
1758  nblocks = consdata->nblocks;
1759  vars = consdata->vars;
1760  vals = consdata->vals;
1761  weights = consdata->weights;
1762  cases = consdata->cases;
1763 
1764  /* get solution */
1765  copyValues(scip, consdata, sol);
1766  SCIPdebugMessage("Checking orbitope constraint <%s> ...\n", SCIPconsGetName(conss[c]));
1767 
1768  /* check upper right triangle (if not yet fixed to zero or in debug mode */
1769 #ifdef NDEBUG
1770  if ( ! consdata->istrianglefixed )
1771 #endif
1772  {
1773  int diagsize;
1774 
1775  /* get last row of triangle */
1776  diagsize = nblocks;
1777  if ( nspcons < nblocks )
1778  diagsize = nspcons;
1779 
1780  /* check variables */
1781  for (i = 0; i < diagsize; ++i)
1782  {
1783  for (j = i+1; j < nblocks; ++j)
1784  {
1785  if ( ! SCIPisFeasZero(scip, vals[i][j]) )
1786  {
1787  if ( printreason )
1788  SCIPinfoMessage(scip, NULL, "variable x[%d][%d] = %f on upper right nonzero.\n", i, j, vals[i][j]);
1789  *result = SCIP_INFEASIBLE;
1790  return SCIP_OKAY;
1791  }
1792  }
1793  }
1794  }
1795 
1796  /* compute table */
1797  computeSCTable(scip, nspcons, nblocks, weights, cases, vals);
1798 
1799  /* loop through rows */
1800  for (i = 1; i < nspcons; ++i)
1801  {
1802  SCIP_Real bar;
1803  int lastcolumn;
1804 
1805  lastcolumn = nblocks - 1;
1806  bar = 0.0;
1807  /* last column considered as part of the bar: */
1808  if ( lastcolumn > i )
1809  lastcolumn = i;
1810 
1811  /* traverse row from right to left */
1812  for (j = lastcolumn; j > 1; --j)
1813  {
1814  bar += vals[i][j];
1815  assert( SCIPisFeasIntegral(scip, vals[i][j]) );
1816 
1817  /* check whether weights[i-1][j-1] < bar (<=> bar - weights[i-1][j-1] > 0), i.e. cut is violated) */
1818  if ( SCIPisGT(scip, bar - weights[i-1][j-1], 0.0) )
1819  {
1820  SCIPdebugMessage("Solution is infeasible.\n");
1821  *result = SCIP_INFEASIBLE;
1822 
1823  if ( printreason )
1824  {
1825  int l;
1826  int p1;
1827  int p2;
1828 
1829  SCIPinfoMessage(scip, NULL, "violated SCI: bar(");
1830 
1831  /* first output bar */
1832  for (l = j; l < nblocks; ++l)
1833  SCIPinfoMessage(scip, NULL, "<%s> (%f)", SCIPvarGetName(vars[i][l]), consdata->vals[i][l]);
1834 
1835  SCIPinfoMessage(scip, NULL, ") SC(");
1836 
1837  /* output shifted column */
1838  p1 = i-1;
1839  p2 = j-1;
1840  do
1841  {
1842  assert( cases[p1][p2] != -1 );
1843  assert( p1 >= 0 && p1 < i );
1844  assert( p2 >= 0 && p2 < j );
1845 
1846  /* if case 1 */
1847  if (cases[p1][p2] == 1)
1848  --p2; /* decrease column */
1849  else
1850  {
1851  /* case 2 or 3: */
1852  assert( cases[p1][p2] == 2 || cases[p1][p2] == 3 );
1853  SCIPinfoMessage(scip, NULL, "<%s> (%f)", SCIPvarGetName(vars[p1][p2]), consdata->vals[p1][p2]);
1854  if ( cases[p1][p2] == 3 )
1855  break;
1856  }
1857  --p1; /* decrease row */
1858  }
1859  while ( p1 >= 0 ); /* should always be true, i.e. the break should end the loop */
1860  assert( cases[p1][p2] == 3 );
1861 
1862  SCIPinfoMessage(scip, NULL, ")");
1863  }
1864 
1865  return SCIP_OKAY;
1866  }
1867  }
1868  }
1869  }
1870  SCIPdebugMessage("Solution is feasible.\n");
1871 
1872  return SCIP_OKAY;
1873 }
1874 
1875 
1876 /** domain propagation method of constraint handler */
1877 static
1878 SCIP_DECL_CONSPROP(consPropOrbitope)
1879 { /*lint --e{715}*/
1880  SCIP_Bool infeasible = FALSE;
1881  int nfixedvars = 0;
1882  int c;
1883 
1884  assert( scip != NULL );
1885  assert( conshdlr != NULL );
1886  assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
1887  assert( result != NULL );
1888 
1889  *result = SCIP_DIDNOTRUN;
1890 
1891  /* propagate all useful constraints */
1892  for (c = 0; c < nusefulconss && !infeasible; ++c)
1893  {
1894  assert( conss[c] != 0 );
1895 
1896  SCIPdebugMessage("Propagation of orbitope constraint <%s> ...\n", SCIPconsGetName(conss[c]));
1897 
1898  SCIP_CALL( propagateCons(scip, conss[c], &infeasible, &nfixedvars) );
1899  }
1900 
1901  /* return the correct result */
1902  if ( infeasible )
1903  {
1904  *result = SCIP_CUTOFF;
1905  SCIPdebugMessage("Propagation via orbitopal fixing proved node to be infeasible.\n");
1906  }
1907  else if ( nfixedvars > 0 )
1908  {
1909  *result = SCIP_REDUCEDDOM;
1910  SCIPdebugMessage("Propagated %d variables via orbitopal fixing.\n", nfixedvars);
1911  }
1912  else if ( nusefulconss > 0 )
1913  {
1914  *result = SCIP_DIDNOTFIND;
1915  SCIPdebugMessage("Propagation via orbitopal fixing did not find anything.\n");
1916  }
1917 
1918  return SCIP_OKAY;
1919 }
1920 
1921 
1922 /** presolving method of constraint handler */
1923 static
1924 SCIP_DECL_CONSPRESOL(consPresolOrbitope)
1925 { /*lint --e{715}*/
1926  SCIP_Bool infeasible = FALSE;
1927  int noldfixedvars;
1928  int c;
1929 
1930  assert( scip != NULL );
1931  assert( conshdlr != NULL );
1932  assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
1933  assert( result != NULL );
1934 
1935  *result = SCIP_DIDNOTRUN;
1936  noldfixedvars = *nfixedvars;
1937 
1938  /* propagate all useful constraints */
1939  for (c = 0; c < nconss && !infeasible; ++c)
1940  {
1941  int nfixed = 0;
1942 
1943  assert( conss[c] != 0 );
1944 
1945  SCIPdebugMessage("Presolving of orbitope constraint <%s> ...\n", SCIPconsGetName(conss[c]));
1946 
1947  SCIP_CALL( propagateCons(scip, conss[c], &infeasible, &nfixed) );
1948  *nfixedvars += nfixed;
1949  }
1950 
1951  if ( infeasible )
1952  {
1953  *result = SCIP_CUTOFF;
1954  SCIPdebugMessage("Presolving detected infeasibility.\n");
1955  }
1956  else if ( *nfixedvars > noldfixedvars )
1957  {
1958  *result = SCIP_SUCCESS;
1959  }
1960  else if ( nconss > 0 )
1961  {
1962  *result = SCIP_DIDNOTFIND;
1963  SCIPdebugMessage("Presolving via orbitopal fixing did not find anything.\n");
1964  }
1965 
1966  return SCIP_OKAY;
1967 }
1968 
1969 
1970 /** propagation conflict resolving method of constraint handler */
1971 static
1972 SCIP_DECL_CONSRESPROP(consRespropOrbitope)
1973 { /*lint --e{715}*/
1974  assert( scip != NULL );
1975  assert( cons != NULL );
1976  assert( infervar != NULL );
1977  assert( bdchgidx != NULL );
1978  assert( result != NULL );
1979 
1980  SCIP_CALL( resolvePropagation(scip, cons, infervar, inferinfo, boundtype, bdchgidx, result) );
1981 
1982  return SCIP_OKAY;
1983 }
1984 
1985 
1986 /** variable rounding lock method of constraint handler */
1987 static
1988 SCIP_DECL_CONSLOCK(consLockOrbitope)
1989 { /*lint --e{715}*/
1990  SCIP_CONSDATA* consdata;
1991  SCIP_VAR*** vars;
1992  int i;
1993  int j;
1994  int nspcons;
1995  int nblocks;
1996 
1997  assert( scip != NULL );
1998  assert( conshdlr != NULL );
1999  assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
2000  assert( cons != NULL );
2001 
2002  consdata = SCIPconsGetData(cons);
2003  assert( consdata != NULL );
2004  assert( consdata->nspcons > 0 );
2005  assert( consdata->nblocks > 0 );
2006  assert( consdata->vars != NULL );
2007 
2008  SCIPdebugMessage("Locking method for orbitope constraint handler\n");
2009 
2010  nspcons = consdata->nspcons;
2011  nblocks = consdata->nblocks;
2012  vars = consdata->vars;
2013 
2014  /* add up locks and down locks on each variable */
2015  for (i = 0; i < nspcons; ++i)
2016  {
2017  for (j = 0; j < nblocks; ++j)
2018  SCIP_CALL( SCIPaddVarLocks(scip, vars[i][j], nlockspos + nlocksneg, nlockspos + nlocksneg) );
2019  }
2020 
2021  return SCIP_OKAY;
2022 }
2023 
2024 
2025 /** constraint display method of constraint handler */
2026 static
2027 SCIP_DECL_CONSPRINT(consPrintOrbitope)
2028 {
2029  SCIP_CONSDATA* consdata;
2030  SCIP_VAR*** vars;
2031  int i;
2032  int j;
2033  int nspcons;
2034  int nblocks;
2035 
2036  assert( scip != NULL );
2037  assert( conshdlr != NULL );
2038  assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
2039  assert( cons != NULL );
2040 
2041  consdata = SCIPconsGetData(cons);
2042  assert( consdata != NULL );
2043  assert( consdata->nspcons > 0 );
2044  assert( consdata->nblocks > 0 );
2045  assert( consdata->vars != NULL );
2046 
2047  nspcons = consdata->nspcons;
2048  nblocks = consdata->nblocks;
2049  vars = consdata->vars;
2050 
2051  SCIPdebugMessage("Printing method for orbitope constraint handler\n");
2052 
2053  if ( consdata->ispart )
2054  SCIPinfoMessage(scip, file, "partOrbitope(");
2055  else
2056  SCIPinfoMessage(scip, file, "packOrbitope(");
2057 
2058  for (i = 0; i < nspcons; ++i)
2059  {
2060  for (j = 0; j < nblocks; ++j)
2061  {
2062  if ( j > 0 )
2063  SCIPinfoMessage(scip, file, ",");
2064  SCIPinfoMessage(scip, file, "%s", SCIPvarGetName(vars[i][j]));
2065  }
2066  if ( i < nspcons-1 )
2067  SCIPinfoMessage(scip, file, ".");
2068  }
2069  SCIPinfoMessage(scip, file, ")");
2070 
2071  return SCIP_OKAY;
2072 }
2073 
2074 
2075 /** constraint copying method of constraint handler */
2076 static
2077 SCIP_DECL_CONSCOPY(consCopyOrbitope)
2078 {
2079  SCIP_CONSDATA* sourcedata;
2080  SCIP_VAR*** sourcevars;
2081  SCIP_VAR*** vars;
2082  int i;
2083  int j;
2084  int nspcons;
2085  int nblocks;
2086 
2087  assert( scip != NULL );
2088  assert( cons != NULL );
2089  assert( sourcescip != NULL );
2090  assert( sourceconshdlr != NULL );
2091  assert( strcmp(SCIPconshdlrGetName(sourceconshdlr), CONSHDLR_NAME) == 0 );
2092  assert( sourcecons != NULL );
2093  assert( varmap != NULL );
2094 
2095  *valid = TRUE;
2096 
2097  SCIPdebugMessage("Copying method for orbitope constraint handler.\n");
2098 
2099  sourcedata = SCIPconsGetData(sourcecons);
2100  assert(sourcedata != NULL);
2101  assert(sourcedata->nspcons > 0);
2102  assert(sourcedata->nblocks > 0);
2103  assert(sourcedata->vars != NULL);
2104 
2105  nspcons = sourcedata->nspcons;
2106  nblocks = sourcedata->nblocks;
2107  sourcevars = sourcedata->vars;
2108 
2109  SCIP_CALL( SCIPallocBufferArray(scip, &vars, nspcons) );
2110  for (i = 0; i < nspcons && *valid; ++i)
2111  {
2112  SCIP_CALL( SCIPallocBufferArray(scip, &(vars[i]), nblocks) ); /*lint !e866*/
2113 
2114  for (j = 0; j < nblocks && *valid; ++j)
2115  {
2116  SCIP_CALL( SCIPgetVarCopy(sourcescip, scip, sourcevars[i][j], &vars[i][j], varmap, consmap, global, valid) );
2117  assert(!(*valid) || vars[i][j] != NULL);
2118  }
2119  }
2120 
2121  /* only create the target constraint, if all variables could be copied */
2122  if ( *valid )
2123  {
2124  /* create copied constraint */
2125  if ( name == NULL )
2126  name = SCIPconsGetName(sourcecons);
2127 
2128  SCIP_CALL( SCIPcreateConsOrbitope(scip, cons, name,
2129  vars, sourcedata->ispart, nspcons, nblocks, sourcedata->resolveprop,
2130  initial, separate, enforce, check, propagate, local, modifiable, dynamic, removable, stickingatnode) );
2131  }
2132 
2133  for (i = 0; i < nspcons; ++i)
2134  SCIPfreeBufferArray(scip, &vars[i]);
2135  SCIPfreeBufferArray(scip, &vars);
2136 
2137  return SCIP_OKAY;
2138 }
2139 
2140 
2141 /** constraint parsing method of constraint handler */
2142 static
2143 SCIP_DECL_CONSPARSE(consParseOrbitope)
2144 { /*lint --e{715}*/
2145  const char* s;
2146  SCIP_Bool ispart;
2147  char varname[SCIP_MAXSTRLEN];
2148  SCIP_VAR*** vars;
2149  SCIP_VAR* var;
2150  int nspcons;
2151  int maxnspcons;
2152  int nblocks;
2153  int maxnblocks;
2154  int k;
2155  int j;
2156 
2157  assert( success != NULL );
2158 
2159  *success = TRUE;
2160  s = str;
2161 
2162  /* skip white space */
2163  while ( *s != '\0' && isspace((unsigned char)*s) )
2164  ++s;
2165 
2166  ispart = FALSE;
2167  if ( strncmp(s, "partOrbitope(", 13) == 0 )
2168  ispart = TRUE;
2169  else
2170  {
2171  if ( strncmp(s, "packOrbitope(", 13) != 0 )
2172  {
2173  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "Syntax error - expected \"partOrbitope\" or \"packOrbitope\": %s\n", s);
2174  *success = FALSE;
2175  return SCIP_OKAY;
2176  }
2177  }
2178  s += 13;
2179 
2180  /* loop through string */
2181  nspcons = 0;
2182  nblocks = 0;
2183  maxnspcons = 10;
2184  maxnblocks = 10;
2185 
2186  SCIP_CALL( SCIPallocBufferArray(scip, &vars, maxnspcons) );
2187  SCIP_CALL( SCIPallocBufferArray(scip, &(vars[0]), maxnblocks) );
2188 
2189  j = 0;
2190  do
2191  {
2192  /* find variable name */
2193  k = 0;
2194  while ( *s != '\0' && ! isspace((unsigned char)*s) && *s != ',' && *s != '.' && *s != ')' )
2195  varname[k++] = *s++;
2196  varname[k] = '\0';
2197 
2198  /* get variable */
2199  var = SCIPfindVar(scip, varname);
2200  if ( var == NULL )
2201  {
2202  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "unknown variable <%s>\n", varname);
2203  *success = FALSE;
2204  return SCIP_OKAY;
2205  }
2206  vars[nspcons][j++] = var;
2207 
2208  if ( j > nblocks )
2209  {
2210  int newsize;
2211 
2212  if ( nspcons > 0 )
2213  {
2214  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "variables per row do not match.\n");
2215  *success = FALSE;
2216  return SCIP_OKAY;
2217  }
2218 
2219  nblocks = j;
2220  newsize = SCIPcalcMemGrowSize(scip, nblocks);
2221  SCIP_CALL( SCIPreallocBufferArray(scip, &(vars[nspcons]), newsize) ); /*lint !e866*/
2222  maxnblocks = newsize;
2223  assert( nblocks <= maxnblocks );
2224  }
2225 
2226  /* skip white space and ',' */
2227  while ( *s != '\0' && ( isspace((unsigned char)*s) || *s == ',' ) )
2228  ++s;
2229 
2230  /* begin new row if required */
2231  if ( *s == '.' )
2232  {
2233  ++nspcons;
2234  ++s;
2235 
2236  if ( nspcons >= maxnspcons )
2237  {
2238  int newsize;
2239 
2240  newsize = SCIPcalcMemGrowSize(scip, nspcons+1);
2241  SCIP_CALL( SCIPreallocBufferArray(scip, &vars, newsize) );
2242  maxnspcons = newsize;
2243  }
2244  assert(nspcons < maxnspcons);
2245 
2246  SCIP_CALL( SCIPallocBufferArray(scip, &(vars[nspcons]), nblocks) ); /*lint !e866*/
2247  j = 0;
2248  }
2249  }
2250  while ( *s != ')' );
2251  ++nspcons;
2252 
2253  SCIP_CALL( SCIPcreateConsOrbitope(scip, cons, name, vars, ispart, nspcons, nblocks, TRUE,
2254  initial, separate, enforce, check, propagate, local, modifiable, dynamic, removable, stickingatnode) );
2255 
2256  for (k = 0; k < nspcons; ++k)
2257  SCIPfreeBufferArray(scip, &vars[k]);
2258  SCIPfreeBufferArray(scip, &vars);
2259 
2260  return SCIP_OKAY;
2261 }
2262 
2263 
2264 /** constraint method of constraint handler which returns the variables (if possible) */
2265 static
2266 SCIP_DECL_CONSGETVARS(consGetVarsOrbitope)
2267 { /*lint --e{715}*/
2268  SCIP_CONSDATA* consdata;
2269 
2270  assert( cons != NULL );
2271  assert( success != NULL );
2272  assert( vars != NULL );
2273 
2274  consdata = SCIPconsGetData(cons);
2275  assert( consdata != NULL );
2276 
2277  if ( varssize < consdata->nblocks * consdata->nspcons )
2278  (*success) = FALSE;
2279  else
2280  {
2281  int cnt = 0;
2282  int i;
2283  int j;
2284 
2285  for (i = 0; i < consdata->nspcons; ++i)
2286  {
2287  for (j = 0; j < consdata->nblocks; ++j)
2288  vars[cnt++] = consdata->vars[i][j];
2289  }
2290  (*success) = TRUE;
2291  }
2292 
2293  return SCIP_OKAY;
2294 }
2295 
2296 
2297 /** constraint method of constraint handler which returns the number of variables (if possible) */
2298 static
2299 SCIP_DECL_CONSGETNVARS(consGetNVarsOrbitope)
2300 { /*lint --e{715}*/
2301  SCIP_CONSDATA* consdata;
2302 
2303  assert( cons != NULL );
2304 
2305  consdata = SCIPconsGetData(cons);
2306  assert( consdata != NULL );
2307 
2308  (*nvars) = consdata->nblocks * consdata->nspcons;
2309  (*success) = TRUE;
2310 
2311  return SCIP_OKAY;
2312 }
2313 
2314 
2315 /*
2316  * constraint specific interface methods
2317  */
2318 
2319 /** creates the handler for orbitope constraints and includes it in SCIP */
2320 SCIP_RETCODE SCIPincludeConshdlrOrbitope(
2321  SCIP* scip /**< SCIP data structure */
2322  )
2323 {
2324  SCIP_CONSHDLRDATA* conshdlrdata;
2325  SCIP_CONSHDLR* conshdlr;
2326 
2327  /* create orbitope constraint handler data */
2328  conshdlrdata = NULL;
2329 
2330  /* include constraint handler */
2331  SCIP_CALL( SCIPincludeConshdlrBasic(scip, &conshdlr, CONSHDLR_NAME, CONSHDLR_DESC,
2332  CONSHDLR_ENFOPRIORITY, CONSHDLR_CHECKPRIORITY,
2333  CONSHDLR_EAGERFREQ, CONSHDLR_NEEDSCONS,
2334  consEnfolpOrbitope, consEnfopsOrbitope, consCheckOrbitope, consLockOrbitope,
2335  conshdlrdata) );
2336  assert(conshdlr != NULL);
2337 
2338  /* set non-fundamental callbacks via specific setter functions */
2339  SCIP_CALL( SCIPsetConshdlrCopy(scip, conshdlr, conshdlrCopyOrbitope, consCopyOrbitope) );
2340  SCIP_CALL( SCIPsetConshdlrDelete(scip, conshdlr, consDeleteOrbitope) );
2341  SCIP_CALL( SCIPsetConshdlrGetVars(scip, conshdlr, consGetVarsOrbitope) );
2342  SCIP_CALL( SCIPsetConshdlrGetNVars(scip, conshdlr, consGetNVarsOrbitope) );
2343  SCIP_CALL( SCIPsetConshdlrParse(scip, conshdlr, consParseOrbitope) );
2344  SCIP_CALL( SCIPsetConshdlrPresol(scip, conshdlr, consPresolOrbitope, CONSHDLR_MAXPREROUNDS, CONSHDLR_DELAYPRESOL) );
2345  SCIP_CALL( SCIPsetConshdlrPrint(scip, conshdlr, consPrintOrbitope) );
2346  SCIP_CALL( SCIPsetConshdlrProp(scip, conshdlr, consPropOrbitope, CONSHDLR_PROPFREQ, CONSHDLR_DELAYPROP,
2347  CONSHDLR_PROP_TIMING) );
2348  SCIP_CALL( SCIPsetConshdlrResprop(scip, conshdlr, consRespropOrbitope) );
2349  SCIP_CALL( SCIPsetConshdlrSepa(scip, conshdlr, consSepalpOrbitope, consSepasolOrbitope, CONSHDLR_SEPAFREQ,
2350  CONSHDLR_SEPAPRIORITY, CONSHDLR_DELAYSEPA) );
2351  SCIP_CALL( SCIPsetConshdlrTrans(scip, conshdlr, consTransOrbitope) );
2352 
2353  return SCIP_OKAY;
2354 }
2355 
2356 
2357 /** creates and captures a orbitope constraint
2358  *
2359  * @note the constraint gets captured, hence at one point you have to release it using the method SCIPreleaseCons()
2360  */
2361 SCIP_RETCODE SCIPcreateConsOrbitope(
2362  SCIP* scip, /**< SCIP data structure */
2363  SCIP_CONS** cons, /**< pointer to hold the created constraint */
2364  const char* name, /**< name of constraint */
2365  SCIP_VAR*** vars, /**< matrix of variables on which the symmetry acts */
2366  SCIP_Bool ispart, /**< whether we deal with the partitioning case (packing otherwise) */
2367  int nspcons, /**< number of set partitioning/packing constraints <=> p */
2368  int nblocks, /**< number of symmetric variable blocks <=> q */
2369  SCIP_Bool resolveprop, /**< should propagation be resolved? */
2370  SCIP_Bool initial, /**< should the LP relaxation of constraint be in the initial LP?
2371  * Usually set to TRUE. Set to FALSE for 'lazy constraints'. */
2372  SCIP_Bool separate, /**< should the constraint be separated during LP processing?
2373  * Usually set to TRUE. */
2374  SCIP_Bool enforce, /**< should the constraint be enforced during node processing?
2375  * TRUE for model constraints, FALSE for additional, redundant constraints. */
2376  SCIP_Bool check, /**< should the constraint be checked for feasibility?
2377  * TRUE for model constraints, FALSE for additional, redundant constraints. */
2378  SCIP_Bool propagate, /**< should the constraint be propagated during node processing?
2379  * Usually set to TRUE. */
2380  SCIP_Bool local, /**< is constraint only valid locally?
2381  * Usually set to FALSE. Has to be set to TRUE, e.g., for branching constraints. */
2382  SCIP_Bool modifiable, /**< is constraint modifiable (subject to column generation)?
2383  * Usually set to FALSE. In column generation applications, set to TRUE if pricing
2384  * adds coefficients to this constraint. */
2385  SCIP_Bool dynamic, /**< is constraint subject to aging?
2386  * Usually set to FALSE. Set to TRUE for own cuts which
2387  * are separated as constraints. */
2388  SCIP_Bool removable, /**< should the relaxation be removed from the LP due to aging or cleanup?
2389  * Usually set to FALSE. Set to TRUE for 'lazy constraints' and 'user cuts'. */
2390  SCIP_Bool stickingatnode /**< should the constraint always be kept at the node where it was added, even
2391  * if it may be moved to a more global node?
2392  * Usually set to FALSE. Set to TRUE to for constraints that represent node data. */
2393  )
2394 {
2395  SCIP_CONSHDLR* conshdlr;
2396  SCIP_CONSDATA* consdata;
2397 
2398  /* find the orbitope constraint handler */
2399  conshdlr = SCIPfindConshdlr(scip, CONSHDLR_NAME);
2400  if ( conshdlr == NULL )
2401  {
2402  SCIPerrorMessage("orbitope constraint handler not found\n");
2403  return SCIP_PLUGINNOTFOUND;
2404  }
2405 
2406  assert( nspcons > 0 );
2407  assert( nblocks > 0 );
2408 
2409  /* run some checks */
2410 #ifndef NDEBUG
2411  {
2412  SCIP_Real obj;
2413  int i;
2414  int j;
2415  for (i = 0; i < nspcons; ++i)
2416  {
2417  /* initialize obj to infinity */
2418  obj = SCIPinfinity(scip);
2419  for (j = 0; j < nblocks; ++j)
2420  {
2421  SCIP_Bool fixedZero;
2422  SCIP_VAR* var;
2423 
2424  var = vars[i][j];
2425  assert(var != NULL);
2426 
2427  /* all variables need to be binary */
2428  assert( SCIPvarGetType(var) == SCIP_VARTYPE_BINARY );
2429 
2430  /* fixed variables have obj = 0; for variables fixed to 0, we assume that there is no
2431  problem (but we cannot always check it, e.g., when in the original problem
2432  variables were fixed and this problem was copied.) */
2433  fixedZero = ( SCIPisZero(scip, SCIPvarGetLbGlobal(var)) && SCIPisZero(scip, SCIPvarGetUbGlobal(var)) );
2434 
2435  /* check whether all variables in a row have the same objective */
2436  if ( ! fixedZero && SCIPisInfinity(scip, obj) )
2437  obj = SCIPvarGetObj(var);
2438  else
2439  {
2440  assert( fixedZero || SCIPisEQ(scip, obj, SCIPvarGetObj(var)) );
2441  }
2442  }
2443  }
2444  }
2445 #endif
2446 
2447  /* create constraint data */
2448  SCIP_CALL( consdataCreate(scip, &consdata, vars, nspcons, nblocks, ispart, resolveprop) );
2449 
2450  /* create constraint */
2451  SCIP_CALL( SCIPcreateCons(scip, cons, name, conshdlr, consdata, initial, separate, enforce, check, propagate,
2452  local, modifiable, dynamic, removable, stickingatnode) );
2453 
2454  return SCIP_OKAY;
2455 }
2456 
2457 /** creates and captures an orbitope constraint
2458  * in its most basic variant, i. e., with all constraint flags set to their default values
2459  *
2460  * @note the constraint gets captured, hence at one point you have to release it using the method SCIPreleaseCons()
2461  */
2462 SCIP_RETCODE SCIPcreateConsBasicOrbitope(
2463  SCIP* scip, /**< SCIP data structure */
2464  SCIP_CONS** cons, /**< pointer to hold the created constraint */
2465  const char* name, /**< name of constraint */
2466  SCIP_VAR*** vars, /**< matrix of variables on which the symmetry acts */
2467  SCIP_Bool ispart, /**< whether we deal with the partitioning case (packing otherwise) */
2468  int nspcons, /**< number of set partitioning/packing constraints <=> p */
2469  int nblocks, /**< number of symmetric variable blocks <=> q */
2470  SCIP_Bool resolveprop /**< should propagation be resolved? */
2471  )
2472 {
2473  SCIP_CALL( SCIPcreateConsOrbitope(scip, cons, name, vars, ispart, nspcons, nblocks, resolveprop,
2474  TRUE, TRUE, TRUE, TRUE, TRUE,
2475  FALSE, FALSE, FALSE, FALSE, FALSE) );
2476 
2477  return SCIP_OKAY;
2478 }
2479