Scippy

SCIP

Solving Constraint Integer Programs

heur_intshifting.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 heur_intshifting.c
17  * @brief LP rounding heuristic that tries to recover from intermediate infeasibilities, shifts integer variables, and
18  * solves a final LP to calculate feasible values for continuous variables
19  * @author Tobias Achterberg
20  */
21 
22 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
23 
24 #include <assert.h>
25 #include <string.h>
26 
27 #include "scip/heur_intshifting.h"
28 
29 
30 #define HEUR_NAME "intshifting"
31 #define HEUR_DESC "LP rounding heuristic with infeasibility recovering and final LP solving"
32 #define HEUR_DISPCHAR 'i'
33 #define HEUR_PRIORITY -10000
34 #define HEUR_FREQ 10
35 #define HEUR_FREQOFS 0
36 #define HEUR_MAXDEPTH -1
37 #define HEUR_TIMING SCIP_HEURTIMING_AFTERLPPLUNGE
38 #define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */
39 
40 #define MAXSHIFTINGS 50 /**< maximal number of non improving shiftings */
41 #define WEIGHTFACTOR 1.1
42 
43 
44 /* locally defined heuristic data */
45 struct SCIP_HeurData
46 {
47  SCIP_SOL* sol; /**< working solution */
48  SCIP_Longint lastlp; /**< last LP number where the heuristic was applied */
49  unsigned int randseed; /**< seed value for random number generator */
50 };
51 
52 
53 /*
54  * local methods
55  */
56 
57 /** update row violation arrays after a row's activity value changed */
58 static
59 void updateViolations(
60  SCIP* scip, /**< SCIP data structure */
61  SCIP_ROW* row, /**< LP row */
62  SCIP_ROW** violrows, /**< array with currently violated rows */
63  int* violrowpos, /**< position of LP rows in violrows array */
64  int* nviolrows, /**< pointer to the number of currently violated rows */
65  SCIP_Real oldminactivity, /**< old minimal activity value of LP row */
66  SCIP_Real oldmaxactivity, /**< old maximal activity value of LP row */
67  SCIP_Real newminactivity, /**< new minimal activity value of LP row */
68  SCIP_Real newmaxactivity /**< new maximal activity value of LP row */
69  )
70 {
71  SCIP_Real lhs;
72  SCIP_Real rhs;
73  SCIP_Bool oldviol;
74  SCIP_Bool newviol;
75 
76  assert(violrows != NULL);
77  assert(violrowpos != NULL);
78  assert(nviolrows != NULL);
79 
80  lhs = SCIProwGetLhs(row);
81  rhs = SCIProwGetRhs(row);
82  oldviol = (SCIPisFeasLT(scip, oldmaxactivity, lhs) || SCIPisFeasGT(scip, oldminactivity, rhs));
83  newviol = (SCIPisFeasLT(scip, newmaxactivity, lhs) || SCIPisFeasGT(scip, newminactivity, rhs));
84  if( oldviol != newviol )
85  {
86  int rowpos;
87  int violpos;
88 
89  rowpos = SCIProwGetLPPos(row);
90  assert(rowpos >= 0);
91 
92  if( oldviol )
93  {
94  /* the row violation was repaired: remove row from violrows array, decrease violation count */
95  violpos = violrowpos[rowpos];
96  assert(0 <= violpos && violpos < *nviolrows);
97  assert(violrows[violpos] == row);
98  violrowpos[rowpos] = -1;
99  if( violpos != *nviolrows-1 )
100  {
101  violrows[violpos] = violrows[*nviolrows-1];
102  violrowpos[SCIProwGetLPPos(violrows[violpos])] = violpos;
103  }
104  (*nviolrows)--;
105  }
106  else
107  {
108  /* the row is now violated: add row to violrows array, increase violation count */
109  assert(violrowpos[rowpos] == -1);
110  violrows[*nviolrows] = row;
111  violrowpos[rowpos] = *nviolrows;
112  (*nviolrows)++;
113  }
114  }
115 }
116 
117 /** update row activities after a variable's solution value changed */
118 static
119 SCIP_RETCODE updateActivities(
120  SCIP* scip, /**< SCIP data structure */
121  SCIP_Real* minactivities, /**< LP row minimal activities */
122  SCIP_Real* maxactivities, /**< LP row maximal activities */
123  SCIP_ROW** violrows, /**< array with currently violated rows */
124  int* violrowpos, /**< position of LP rows in violrows array */
125  int* nviolrows, /**< pointer to the number of currently violated rows */
126  int nlprows, /**< number of rows in current LP */
127  SCIP_VAR* var, /**< variable that has been changed */
128  SCIP_Real oldsolval, /**< old solution value of variable */
129  SCIP_Real newsolval /**< new solution value of variable */
130  )
131 {
132  SCIP_COL* col;
133  SCIP_ROW** colrows;
134  SCIP_Real* colvals;
135  SCIP_Real delta;
136  int ncolrows;
137  int r;
138 
139  assert(minactivities != NULL);
140  assert(maxactivities != NULL);
141  assert(nviolrows != NULL);
142  assert(0 <= *nviolrows && *nviolrows <= nlprows);
143  assert(SCIPvarGetType(var) == SCIP_VARTYPE_BINARY || SCIPvarGetType(var) == SCIP_VARTYPE_INTEGER);
144 
145  delta = newsolval - oldsolval;
146  col = SCIPvarGetCol(var);
147  colrows = SCIPcolGetRows(col);
148  colvals = SCIPcolGetVals(col);
149  ncolrows = SCIPcolGetNLPNonz(col);
150  assert(ncolrows == 0 || (colrows != NULL && colvals != NULL));
151 
152  for( r = 0; r < ncolrows; ++r )
153  {
154  SCIP_ROW* row;
155  int rowpos;
156 
157  row = colrows[r];
158  rowpos = SCIProwGetLPPos(row);
159  assert(-1 <= rowpos && rowpos < nlprows);
160 
161  if( rowpos >= 0 && !SCIProwIsLocal(row) )
162  {
163  SCIP_Real oldminactivity;
164  SCIP_Real oldmaxactivity;
165  SCIP_Real newminactivity;
166  SCIP_Real newmaxactivity;
167 
168  assert(SCIProwIsInLP(row));
169 
170  /* update row activities */
171  oldminactivity = minactivities[rowpos];
172  oldmaxactivity = maxactivities[rowpos];
173 
174  if( !SCIPisInfinity(scip, -oldminactivity) )
175  {
176  newminactivity = oldminactivity + delta * colvals[r];
177  minactivities[rowpos] = newminactivity;
178  }
179  else
180  newminactivity = -SCIPinfinity(scip);
181  if( !SCIPisInfinity(scip, oldmaxactivity) )
182  {
183  newmaxactivity = oldmaxactivity + delta * colvals[r];
184  maxactivities[rowpos] = newmaxactivity;
185  }
186  else
187  newmaxactivity = SCIPinfinity(scip);
188 
189  /* update row violation arrays */
190  updateViolations(scip, row, violrows, violrowpos, nviolrows, oldminactivity, oldmaxactivity,
191  newminactivity, newmaxactivity);
192  }
193  }
194 
195  return SCIP_OKAY;
196 }
197 
198 /** returns an integer variable, that pushes activity of the row in the given direction with minimal negative impact on
199  * other rows;
200  * if variables have equal impact, chooses the one with best objective value improvement in corresponding direction;
201  * prefer fractional integers over other variables in order to become integral during the process;
202  * shifting in a direction is forbidden, if this forces the objective value over the upper bound, or if the variable
203  * was already shifted in the opposite direction
204  */
205 static
206 SCIP_RETCODE selectShifting(
207  SCIP* scip, /**< SCIP data structure */
208  SCIP_SOL* sol, /**< primal solution */
209  SCIP_ROW* row, /**< LP row */
210  SCIP_Real rowactivity, /**< activity of LP row */
211  int direction, /**< should the activity be increased (+1) or decreased (-1)? */
212  SCIP_Real* nincreases, /**< array with weighted number of increasings per variables */
213  SCIP_Real* ndecreases, /**< array with weighted number of decreasings per variables */
214  SCIP_Real increaseweight, /**< current weight of increase/decrease updates */
215  SCIP_VAR** shiftvar, /**< pointer to store the shifting variable, returns NULL if impossible */
216  SCIP_Real* oldsolval, /**< pointer to store old solution value of shifting variable */
217  SCIP_Real* newsolval /**< pointer to store new (shifted) solution value of shifting variable */
218  )
219 {
220  SCIP_COL** rowcols;
221  SCIP_Real* rowvals;
222  int nrowcols;
223  SCIP_Real activitydelta;
224  SCIP_Real bestshiftscore;
225  SCIP_Real bestdeltaobj;
226  int c;
227 
228  assert(direction == +1 || direction == -1);
229  assert(nincreases != NULL);
230  assert(ndecreases != NULL);
231  assert(shiftvar != NULL);
232  assert(oldsolval != NULL);
233  assert(newsolval != NULL);
234 
235  /* get row entries */
236  rowcols = SCIProwGetCols(row);
237  rowvals = SCIProwGetVals(row);
238  nrowcols = SCIProwGetNLPNonz(row);
239 
240  /* calculate how much the activity must be shifted in order to become feasible */
241  activitydelta = (direction == +1 ? SCIProwGetLhs(row) - rowactivity : SCIProwGetRhs(row) - rowactivity);
242  assert((direction == +1 && SCIPisPositive(scip, activitydelta))
243  || (direction == -1 && SCIPisNegative(scip, activitydelta)));
244 
245  /* select shifting variable */
246  bestshiftscore = SCIP_REAL_MAX;
247  bestdeltaobj = SCIPinfinity(scip);
248  *shiftvar = NULL;
249  *newsolval = 0.0;
250  *oldsolval = 0.0;
251  for( c = 0; c < nrowcols; ++c )
252  {
253  SCIP_COL* col;
254  SCIP_VAR* var;
255  SCIP_Real val;
256  SCIP_Real solval;
257  SCIP_Real shiftval;
258  SCIP_Real shiftscore;
259  SCIP_Bool isfrac;
260  SCIP_Bool increase;
261  int probindex;
262 
263  col = rowcols[c];
264  var = SCIPcolGetVar(col);
265  val = rowvals[c];
266  assert(!SCIPisZero(scip, val));
267  solval = SCIPgetSolVal(scip, sol, var);
268 
269  /* only accept integer variables */
270  if( SCIPvarGetType(var) != SCIP_VARTYPE_BINARY && SCIPvarGetType(var) != SCIP_VARTYPE_INTEGER )
271  continue;
272 
273  isfrac = !SCIPisFeasIntegral(scip, solval);
274  increase = (direction * val > 0.0);
275  probindex = SCIPvarGetProbindex(var);
276 
277  /* calculate the score of the shifting (prefer smaller values) */
278  if( isfrac )
279  shiftscore = increase ? -1.0 / (SCIPvarGetNLocksUp(var) + 1.0) :
280  -1.0 / (SCIPvarGetNLocksDown(var) + 1.0);
281  else
282  {
283  if( increase )
284  shiftscore = ndecreases[probindex]/increaseweight;
285  else
286  shiftscore = nincreases[probindex]/increaseweight;
287  shiftscore += 1.0;
288  }
289 
290  if( shiftscore <= bestshiftscore )
291  {
292  SCIP_Real deltaobj;
293 
294  if( !increase )
295  {
296  /* shifting down */
297  assert(direction * val < 0.0);
298  if( isfrac )
299  shiftval = SCIPfeasFloor(scip, solval);
300  else
301  {
302  SCIP_Real lb;
303 
304  assert(activitydelta/val < 0.0);
305  shiftval = solval + activitydelta/val;
306  assert(shiftval <= solval); /* may be equal due to numerical digit erasement in the subtraction */
307  shiftval = SCIPfeasFloor(scip, shiftval);
308  lb = SCIPvarGetLbGlobal(var);
309  shiftval = MAX(shiftval, lb);
310  }
311  }
312  else
313  {
314  /* shifting up */
315  assert(direction * val > 0.0);
316  if( isfrac )
317  shiftval = SCIPfeasCeil(scip, solval);
318  else
319  {
320  SCIP_Real ub;
321 
322  assert(activitydelta/val > 0.0);
323  shiftval = solval + activitydelta/val;
324  assert(shiftval >= solval); /* may be equal due to numerical digit erasement in the subtraction */
325  shiftval = SCIPfeasCeil(scip, shiftval);
326  ub = SCIPvarGetUbGlobal(var);
327  shiftval = MIN(shiftval, ub);
328  }
329  }
330 
331  if( SCIPisEQ(scip, shiftval, solval) )
332  continue;
333 
334  deltaobj = SCIPvarGetObj(var) * (shiftval - solval);
335  if( shiftscore < bestshiftscore || deltaobj < bestdeltaobj )
336  {
337  bestshiftscore = shiftscore;
338  bestdeltaobj = deltaobj;
339  *shiftvar = var;
340  *oldsolval = solval;
341  *newsolval = shiftval;
342  }
343  }
344  }
345 
346  return SCIP_OKAY;
347 }
348 
349 /** returns a fractional variable, that has most impact on rows in opposite direction, i.e. that is most crucial to
350  * fix in the other direction;
351  * if variables have equal impact, chooses the one with best objective value improvement in corresponding direction;
352  * shifting in a direction is forbidden, if this forces the objective value over the upper bound
353  */
354 static
355 SCIP_RETCODE selectEssentialRounding(
356  SCIP* scip, /**< SCIP data structure */
357  SCIP_SOL* sol, /**< primal solution */
358  SCIP_Real minobj, /**< minimal objective value possible after shifting remaining fractional vars */
359  SCIP_VAR** lpcands, /**< fractional variables in LP */
360  int nlpcands, /**< number of fractional variables in LP */
361  SCIP_VAR** shiftvar, /**< pointer to store the shifting variable, returns NULL if impossible */
362  SCIP_Real* oldsolval, /**< old (fractional) solution value of shifting variable */
363  SCIP_Real* newsolval /**< new (shifted) solution value of shifting variable */
364  )
365 {
366  SCIP_VAR* var;
367  SCIP_Real solval;
368  SCIP_Real shiftval;
369  SCIP_Real obj;
370  SCIP_Real deltaobj;
371  SCIP_Real bestdeltaobj;
372  int maxnlocks;
373  int nlocks;
374  int v;
375 
376  assert(shiftvar != NULL);
377  assert(oldsolval != NULL);
378  assert(newsolval != NULL);
379 
380  /* select shifting variable */
381  maxnlocks = -1;
382  bestdeltaobj = SCIPinfinity(scip);
383  *shiftvar = NULL;
384  for( v = 0; v < nlpcands; ++v )
385  {
386  var = lpcands[v];
387  assert(SCIPvarGetType(var) == SCIP_VARTYPE_BINARY || SCIPvarGetType(var) == SCIP_VARTYPE_INTEGER);
388 
389  solval = SCIPgetSolVal(scip, sol, var);
390  if( !SCIPisFeasIntegral(scip, solval) )
391  {
392  obj = SCIPvarGetObj(var);
393 
394  /* shifting down */
395  nlocks = SCIPvarGetNLocksUp(var);
396  if( nlocks >= maxnlocks )
397  {
398  shiftval = SCIPfeasFloor(scip, solval);
399  deltaobj = obj * (shiftval - solval);
400  if( (nlocks > maxnlocks || deltaobj < bestdeltaobj) && minobj - obj < SCIPgetCutoffbound(scip) )
401  {
402  maxnlocks = nlocks;
403  bestdeltaobj = deltaobj;
404  *shiftvar = var;
405  *oldsolval = solval;
406  *newsolval = shiftval;
407  }
408  }
409 
410  /* shifting up */
411  nlocks = SCIPvarGetNLocksDown(var);
412  if( nlocks >= maxnlocks )
413  {
414  shiftval = SCIPfeasCeil(scip, solval);
415  deltaobj = obj * (shiftval - solval);
416  if( (nlocks > maxnlocks || deltaobj < bestdeltaobj) && minobj + obj < SCIPgetCutoffbound(scip) )
417  {
418  maxnlocks = nlocks;
419  bestdeltaobj = deltaobj;
420  *shiftvar = var;
421  *oldsolval = solval;
422  *newsolval = shiftval;
423  }
424  }
425  }
426  }
427 
428  return SCIP_OKAY;
429 }
430 
431 /** adds a given value to the fractionality counters of the rows in which the given variable appears */
432 static
433 void addFracCounter(
434  int* nfracsinrow, /**< array to store number of fractional variables per row */
435  int nlprows, /**< number of rows in LP */
436  SCIP_VAR* var, /**< variable for which the counting should be updated */
437  int incval /**< value that should be added to the corresponding array entries */
438  )
439 {
440  SCIP_COL* col;
441  SCIP_ROW** rows;
442  int nrows;
443  int r;
444 
445  col = SCIPvarGetCol(var);
446  rows = SCIPcolGetRows(col);
447  nrows = SCIPcolGetNLPNonz(col);
448  for( r = 0; r < nrows; ++r )
449  {
450  int rowidx;
451 
452  rowidx = SCIProwGetLPPos(rows[r]);
453  assert(0 <= rowidx && rowidx < nlprows);
454  nfracsinrow[rowidx] += incval;
455  assert(nfracsinrow[rowidx] >= 0);
456  }
457 }
458 
459 
460 /*
461  * Callback methods
462  */
463 
464 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */
465 static
466 SCIP_DECL_HEURCOPY(heurCopyIntshifting)
467 { /*lint --e{715}*/
468  assert(scip != NULL);
469  assert(heur != NULL);
470  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
471 
472  /* call inclusion method of primal heuristic */
473  SCIP_CALL( SCIPincludeHeurIntshifting(scip) );
474 
475  return SCIP_OKAY;
476 }
477 
478 
479 /** initialization method of primal heuristic (called after problem was transformed) */
480 static
481 SCIP_DECL_HEURINIT(heurInitIntshifting) /*lint --e{715}*/
482 { /*lint --e{715}*/
483  SCIP_HEURDATA* heurdata;
484 
485  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
486  assert(SCIPheurGetData(heur) == NULL);
487 
488  /* create heuristic data */
489  SCIP_CALL( SCIPallocMemory(scip, &heurdata) );
490  SCIP_CALL( SCIPcreateSol(scip, &heurdata->sol, heur) );
491  heurdata->lastlp = -1;
492  heurdata->randseed = 0;
493  SCIPheurSetData(heur, heurdata);
494 
495  return SCIP_OKAY;
496 }
497 
498 /** deinitialization method of primal heuristic (called before transformed problem is freed) */
499 static
500 SCIP_DECL_HEUREXIT(heurExitIntshifting) /*lint --e{715}*/
501 { /*lint --e{715}*/
502  SCIP_HEURDATA* heurdata;
503 
504  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
505 
506  /* free heuristic data */
507  heurdata = SCIPheurGetData(heur);
508  assert(heurdata != NULL);
509  SCIP_CALL( SCIPfreeSol(scip, &heurdata->sol) );
510  SCIPfreeMemory(scip, &heurdata);
511  SCIPheurSetData(heur, NULL);
512 
513  return SCIP_OKAY;
514 }
515 
516 /** solving process initialization method of primal heuristic (called when branch and bound process is about to begin) */
517 static
518 SCIP_DECL_HEURINITSOL(heurInitsolIntshifting)
519 {
520  SCIP_HEURDATA* heurdata;
521 
522  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
523 
524  heurdata = SCIPheurGetData(heur);
525  assert(heurdata != NULL);
526  heurdata->lastlp = -1;
527 
528  return SCIP_OKAY;
529 }
530 
531 
532 /** execution method of primal heuristic */
533 static
534 SCIP_DECL_HEUREXEC(heurExecIntshifting) /*lint --e{715}*/
535 { /*lint --e{715}*/
536  SCIP_HEURDATA* heurdata;
537  SCIP_SOL* sol;
538  SCIP_VAR** lpcands;
539  SCIP_Real* lpcandssol;
540  SCIP_ROW** lprows;
541  SCIP_Real* minactivities;
542  SCIP_Real* maxactivities;
543  SCIP_ROW** violrows;
544  SCIP_Real* nincreases;
545  SCIP_Real* ndecreases;
546  int* violrowpos;
547  int* nfracsinrow;
548  SCIP_Real increaseweight;
549  SCIP_Real obj;
550  SCIP_Real bestshiftval;
551  SCIP_Real minobj;
552  int nlpcands;
553  int nlprows;
554  int nvars;
555  int nfrac;
556  int nviolrows;
557  int nprevviolrows;
558  int minnviolrows;
559  int nnonimprovingshifts;
560  int c;
561  int r;
562  SCIP_Longint nlps;
563  SCIP_Longint ncalls;
564  SCIP_Longint nsolsfound;
565  SCIP_Longint nnodes;
566 
567  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
568  assert(scip != NULL);
569  assert(result != NULL);
570  assert(SCIPhasCurrentNodeLP(scip));
571 
572  *result = SCIP_DIDNOTRUN;
573 
574  /* do not call heuristic of node was already detected to be infeasible */
575  if( nodeinfeasible )
576  return SCIP_OKAY;
577 
578  /* don't call heuristic, if no continuous variables are present
579  * -> in this case, it is equivalent to shifting heuristic
580  */
581  if( SCIPgetNContVars(scip) == 0 )
582  return SCIP_OKAY;
583 
584  /* only call heuristic, if an optimal LP solution is at hand */
585  if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
586  return SCIP_OKAY;
587 
588  /* only call heuristic, if the LP objective value is smaller than the cutoff bound */
589  if( SCIPisGE(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip)) )
590  return SCIP_OKAY;
591 
592  /* get heuristic data */
593  heurdata = SCIPheurGetData(heur);
594  assert(heurdata != NULL);
595 
596  /* don't call heuristic, if we have already processed the current LP solution */
597  nlps = SCIPgetNLPs(scip);
598  if( nlps == heurdata->lastlp )
599  return SCIP_OKAY;
600  heurdata->lastlp = nlps;
601 
602  /* don't call heuristic, if it was not successful enough in the past */
603  ncalls = SCIPheurGetNCalls(heur);
604  nsolsfound = 10*SCIPheurGetNBestSolsFound(heur) + SCIPheurGetNSolsFound(heur);
605  nnodes = SCIPgetNNodes(scip);
606  if( nnodes % (ncalls/(nsolsfound+1)+1) != 0 ) /*?????????? ncalls/100 */
607  return SCIP_OKAY;
608 
609  /* get fractional variables, that should be integral */
610  /* todo check if heuristic should include implicit integer variables for its calculations */
611  SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, NULL, &nlpcands, NULL, NULL) );
612  nfrac = nlpcands;
613 
614  /* only call heuristic, if LP solution is fractional */
615  if( nfrac == 0 )
616  return SCIP_OKAY;
617 
618  *result = SCIP_DIDNOTFIND;
619 
620  /* get LP rows */
621  SCIP_CALL( SCIPgetLPRowsData(scip, &lprows, &nlprows) );
622 
623  SCIPdebugMessage("executing intshifting heuristic: %d LP rows, %d fractionals\n", nlprows, nfrac);
624 
625  /* get memory for activities, violated rows, and row violation positions */
626  nvars = SCIPgetNVars(scip);
627  SCIP_CALL( SCIPallocBufferArray(scip, &minactivities, nlprows) );
628  SCIP_CALL( SCIPallocBufferArray(scip, &maxactivities, nlprows) );
629  SCIP_CALL( SCIPallocBufferArray(scip, &violrows, nlprows) );
630  SCIP_CALL( SCIPallocBufferArray(scip, &violrowpos, nlprows) );
631  SCIP_CALL( SCIPallocBufferArray(scip, &nfracsinrow, nlprows) );
632  SCIP_CALL( SCIPallocBufferArray(scip, &nincreases, nvars) );
633  SCIP_CALL( SCIPallocBufferArray(scip, &ndecreases, nvars) );
634  BMSclearMemoryArray(nfracsinrow, nlprows);
635  BMSclearMemoryArray(nincreases, nvars);
636  BMSclearMemoryArray(ndecreases, nvars);
637 
638  /* get the minimal and maximal activity for all globally valid rows for continuous variables in their full range;
639  * these are the values of a*x' with x' being the LP solution for integer variables and the lower or upper bound
640  * for the continuous variables
641  */
642  nviolrows = 0;
643  for( r = 0; r < nlprows; ++r )
644  {
645  SCIP_ROW* row;
646 
647  row = lprows[r];
648  assert(SCIProwGetLPPos(row) == r);
649 
650  if( !SCIProwIsLocal(row) )
651  {
652  SCIP_COL** cols;
653  SCIP_Real* vals;
654  int nnonz;
655  SCIP_Bool mininf;
656  SCIP_Bool maxinf;
657 
658  mininf = FALSE;
659  maxinf = FALSE;
660  minactivities[r] = 0.0;
661  maxactivities[r] = 0.0;
662  cols = SCIProwGetCols(row);
663  vals = SCIProwGetVals(row);
664  nnonz = SCIProwGetNNonz(row);
665  for( c = 0; c < nnonz && !(mininf && maxinf); ++c )
666  {
667  SCIP_VAR* var;
668 
669  var = SCIPcolGetVar(cols[c]);
670  if( SCIPvarGetType(var) == SCIP_VARTYPE_BINARY || SCIPvarGetType(var) == SCIP_VARTYPE_INTEGER )
671  {
672  SCIP_Real act;
673 
674  act = vals[c] * SCIPcolGetPrimsol(cols[c]);
675  minactivities[r] += act;
676  maxactivities[r] += act;
677  }
678  else if( vals[c] > 0.0 )
679  {
680  SCIP_Real lb;
681  SCIP_Real ub;
682 
683  lb = SCIPvarGetLbGlobal(var);
684  ub = SCIPvarGetUbGlobal(var);
685  if( SCIPisInfinity(scip, -lb) )
686  mininf = TRUE;
687  else
688  minactivities[r] += vals[c] * lb;
689  if( SCIPisInfinity(scip, ub) )
690  maxinf = TRUE;
691  else
692  maxactivities[r] += vals[c] * ub;
693  }
694  else if( vals[c] < 0.0 )
695  {
696  SCIP_Real lb;
697  SCIP_Real ub;
698 
699  lb = SCIPvarGetLbGlobal(var);
700  ub = SCIPvarGetUbGlobal(var);
701  if( SCIPisInfinity(scip, ub) )
702  mininf = TRUE;
703  else
704  minactivities[r] += vals[c] * ub;
705  if( SCIPisInfinity(scip, -lb) )
706  maxinf = TRUE;
707  else
708  maxactivities[r] += vals[c] * lb;
709  }
710 
711  if( mininf )
712  minactivities[r] = -SCIPinfinity(scip);
713  if( maxinf )
714  maxactivities[r] = SCIPinfinity(scip);
715  }
716 
717  if( SCIPisFeasLT(scip, maxactivities[r], SCIProwGetLhs(row))
718  || SCIPisFeasGT(scip, minactivities[r], SCIProwGetRhs(row)) )
719  {
720  violrows[nviolrows] = row;
721  violrowpos[r] = nviolrows;
722  nviolrows++;
723  }
724  else
725  violrowpos[r] = -1;
726  }
727  }
728 
729  /* calc the current number of fractional variables in rows */
730  for( c = 0; c < nlpcands; ++c )
731  addFracCounter(nfracsinrow, nlprows, lpcands[c], +1);
732 
733  /* get the working solution from heuristic's local data */
734  sol = heurdata->sol;
735  assert(sol != NULL);
736 
737  /* copy the current LP solution to the working solution */
738  SCIP_CALL( SCIPlinkLPSol(scip, sol) );
739 
740  /* calculate the minimal objective value possible after rounding fractional variables */
741  minobj = SCIPgetSolTransObj(scip, sol);
742  assert(minobj < SCIPgetCutoffbound(scip));
743  for( c = 0; c < nlpcands; ++c )
744  {
745  obj = SCIPvarGetObj(lpcands[c]);
746  bestshiftval = obj > 0.0 ? SCIPfeasFloor(scip, lpcandssol[c]) : SCIPfeasCeil(scip, lpcandssol[c]);
747  minobj += obj * (bestshiftval - lpcandssol[c]);
748  }
749 
750  /* try to shift remaining variables in order to become/stay feasible */
751  nnonimprovingshifts = 0;
752  minnviolrows = INT_MAX;
753  increaseweight = 1.0;
754  while( (nfrac > 0 || nviolrows > 0) && nnonimprovingshifts < MAXSHIFTINGS && !SCIPisStopped(scip) )
755  {
756  SCIP_VAR* shiftvar;
757  SCIP_Real oldsolval;
758  SCIP_Real newsolval;
759  SCIP_Bool oldsolvalisfrac;
760  int probindex;
761 
762  SCIPdebugMessage("intshifting heuristic: nfrac=%d, nviolrows=%d, obj=%g (best possible obj: %g), cutoff=%g\n",
763  nfrac, nviolrows, SCIPgetSolOrigObj(scip, sol), SCIPretransformObj(scip, minobj),
764  SCIPretransformObj(scip, SCIPgetCutoffbound(scip)));
765 
766  nprevviolrows = nviolrows;
767 
768  /* choose next variable to process:
769  * - if a violated row exists, shift a variable decreasing the violation, that has least impact on other rows
770  * - otherwise, shift a variable, that has strongest devastating impact on rows in opposite direction
771  */
772  shiftvar = NULL;
773  oldsolval = 0.0;
774  newsolval = 0.0;
775  if( nviolrows > 0 && (nfrac == 0 || nnonimprovingshifts < MAXSHIFTINGS-1) )
776  {
777  SCIP_ROW* row;
778  int rowidx;
779  int rowpos;
780  int direction;
781 
782  rowidx = -1;
783  rowpos = -1;
784  row = NULL;
785  if( nfrac > 0 )
786  {
787  for( rowidx = nviolrows-1; rowidx >= 0; --rowidx )
788  {
789  row = violrows[rowidx];
790  rowpos = SCIProwGetLPPos(row);
791  assert(violrowpos[rowpos] == rowidx);
792  if( nfracsinrow[rowpos] > 0 )
793  break;
794  }
795  }
796  if( rowidx == -1 )
797  {
798  rowidx = SCIPgetRandomInt(0, nviolrows-1, &heurdata->randseed);
799  row = violrows[rowidx];
800  rowpos = SCIProwGetLPPos(row);
801  assert(0 <= rowpos && rowpos < nlprows);
802  assert(violrowpos[rowpos] == rowidx);
803  assert(nfracsinrow[rowpos] == 0);
804  }
805  assert(violrowpos[rowpos] == rowidx);
806 
807  SCIPdebugMessage("intshifting heuristic: try to fix violated row <%s>: %g <= [%g,%g] <= %g\n",
808  SCIProwGetName(row), SCIProwGetLhs(row), minactivities[rowpos], maxactivities[rowpos], SCIProwGetRhs(row));
809  SCIPdebug( SCIP_CALL( SCIPprintRow(scip, row, NULL) ) );
810 
811  /* get direction in which activity must be shifted */
812  assert(SCIPisFeasLT(scip, maxactivities[rowpos], SCIProwGetLhs(row))
813  || SCIPisFeasGT(scip, minactivities[rowpos], SCIProwGetRhs(row)));
814  direction = SCIPisFeasLT(scip, maxactivities[rowpos], SCIProwGetLhs(row)) ? +1 : -1;
815 
816  /* search an integer variable that can shift the activity in the necessary direction */
817  SCIP_CALL( selectShifting(scip, sol, row, direction == +1 ? maxactivities[rowpos] : minactivities[rowpos],
818  direction, nincreases, ndecreases, increaseweight, &shiftvar, &oldsolval, &newsolval) );
819  }
820 
821  if( shiftvar == NULL && nfrac > 0 )
822  {
823  SCIPdebugMessage("intshifting heuristic: search rounding variable and try to stay feasible\n");
824  SCIP_CALL( selectEssentialRounding(scip, sol, minobj, lpcands, nlpcands, &shiftvar, &oldsolval, &newsolval) );
825  }
826 
827  /* check, whether shifting was possible */
828  if( shiftvar == NULL || SCIPisEQ(scip, oldsolval, newsolval) )
829  {
830  SCIPdebugMessage("intshifting heuristic: -> didn't find a shifting variable\n");
831  break;
832  }
833 
834  assert(SCIPvarGetType(shiftvar) == SCIP_VARTYPE_BINARY || SCIPvarGetType(shiftvar) == SCIP_VARTYPE_INTEGER);
835 
836  SCIPdebugMessage("intshifting heuristic: -> shift var <%s>[%g,%g], type=%d, oldval=%g, newval=%g, obj=%g\n",
837  SCIPvarGetName(shiftvar), SCIPvarGetLbGlobal(shiftvar), SCIPvarGetUbGlobal(shiftvar), SCIPvarGetType(shiftvar),
838  oldsolval, newsolval, SCIPvarGetObj(shiftvar));
839 
840  /* update row activities of globally valid rows */
841  SCIP_CALL( updateActivities(scip, minactivities, maxactivities, violrows, violrowpos, &nviolrows, nlprows,
842  shiftvar, oldsolval, newsolval) );
843  if( nviolrows >= nprevviolrows )
844  nnonimprovingshifts++;
845  else if( nviolrows < minnviolrows )
846  {
847  minnviolrows = nviolrows;
848  nnonimprovingshifts = 0;
849  }
850 
851  /* store new solution value and decrease fractionality counter */
852  SCIP_CALL( SCIPsetSolVal(scip, sol, shiftvar, newsolval) );
853 
854  /* update fractionality counter and minimal objective value possible after shifting remaining variables */
855  oldsolvalisfrac = !SCIPisFeasIntegral(scip, oldsolval);
856  obj = SCIPvarGetObj(shiftvar);
857  if( oldsolvalisfrac )
858  {
859  assert(SCIPisFeasIntegral(scip, newsolval));
860  nfrac--;
861  nnonimprovingshifts = 0;
862  minnviolrows = INT_MAX;
863  addFracCounter(nfracsinrow, nlprows, shiftvar, -1);
864 
865  /* the rounding was already calculated into the minobj -> update only if rounding in "wrong" direction */
866  if( obj > 0.0 && newsolval > oldsolval )
867  minobj += obj;
868  else if( obj < 0.0 && newsolval < oldsolval )
869  minobj -= obj;
870  }
871  else
872  {
873  /* update minimal possible objective value */
874  minobj += obj * (newsolval - oldsolval);
875  }
876 
877  /* update increase/decrease arrays */
878  if( !oldsolvalisfrac )
879  {
880  probindex = SCIPvarGetProbindex(shiftvar);
881  assert(0 <= probindex && probindex < nvars);
882  increaseweight *= WEIGHTFACTOR;
883  if( newsolval < oldsolval )
884  ndecreases[probindex] += increaseweight;
885  else
886  nincreases[probindex] += increaseweight;
887  if( increaseweight >= 1e+09 )
888  {
889  int i;
890 
891  for( i = 0; i < nvars; ++i )
892  {
893  nincreases[i] /= increaseweight;
894  ndecreases[i] /= increaseweight;
895  }
896  increaseweight = 1.0;
897  }
898  }
899 
900  SCIPdebugMessage("intshifting heuristic: -> nfrac=%d, nviolrows=%d, obj=%g (best possible obj: %g)\n",
901  nfrac, nviolrows, SCIPgetSolOrigObj(scip, sol), SCIPretransformObj(scip, minobj));
902  }
903 
904  /* check, if the new solution is potentially feasible and solve the LP to calculate values for the continuous
905  * variables
906  */
907  if( nfrac == 0 && nviolrows == 0 )
908  {
909  SCIP_VAR** vars;
910  SCIP_Bool lperror;
911  int nintvars;
912  int v;
913 #ifdef NDEBUG
914  SCIP_RETCODE retstat;
915 #endif
916 
917  SCIPdebugMessage("shifted solution is potentially feasible -> solve LP to fix continuous variables\n");
918 
919  /* start diving to calculate the LP relaxation */
920  SCIP_CALL( SCIPstartDive(scip) );
921 
922  /* set the bounds of the variables: fixed for integers, global bounds for continuous */
923  vars = SCIPgetVars(scip);
924  nintvars = SCIPgetNBinVars(scip) + SCIPgetNIntVars(scip);
925  for( v = 0; v < nvars; ++v )
926  {
927  if( SCIPvarGetStatus(vars[v]) == SCIP_VARSTATUS_COLUMN )
928  {
929  SCIP_CALL( SCIPchgVarLbDive(scip, vars[v], SCIPvarGetLbGlobal(vars[v])) );
930  SCIP_CALL( SCIPchgVarUbDive(scip, vars[v], SCIPvarGetUbGlobal(vars[v])) );
931  }
932  }
933  for( v = 0; v < nintvars; ++v ) /* apply this after global bounds to not cause an error with intermediate empty domains */
934  {
935  if( SCIPvarGetStatus(vars[v]) == SCIP_VARSTATUS_COLUMN )
936  {
937  SCIP_Real solval;
938  solval = SCIPgetSolVal(scip, sol, vars[v]);
939  SCIP_CALL( SCIPchgVarLbDive(scip, vars[v], solval) );
940  SCIP_CALL( SCIPchgVarUbDive(scip, vars[v], solval) );
941  }
942  }
943 
944  /* solve LP */
945  SCIPdebugMessage(" -> old LP iterations: %"SCIP_LONGINT_FORMAT"\n", SCIPgetNLPIterations(scip));
946 
947  /* Errors in the LP solver should not kill the overall solving process, if the LP is just needed for a heuristic.
948  * Hence in optimized mode, the return code is caught and a warning is printed, only in debug mode, SCIP will stop.
949  */
950 #ifdef NDEBUG
951  retstat = SCIPsolveDiveLP(scip, -1, &lperror, NULL);
952  if( retstat != SCIP_OKAY )
953  {
954  SCIPwarningMessage(scip, "Error while solving LP in Intshifting heuristic; LP solve terminated with code <%d>\n",retstat);
955  }
956 #else
957  SCIP_CALL( SCIPsolveDiveLP(scip, -1, &lperror, NULL) );
958 #endif
959 
960  SCIPdebugMessage(" -> new LP iterations: %"SCIP_LONGINT_FORMAT"\n", SCIPgetNLPIterations(scip));
961  SCIPdebugMessage(" -> error=%u, status=%d\n", lperror, SCIPgetLPSolstat(scip));
962 
963  /* check if this is a feasible solution */
964  if( !lperror && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL )
965  {
966  SCIP_Bool stored;
967 
968  /* copy the current LP solution to the working solution */
969  SCIP_CALL( SCIPlinkLPSol(scip, sol) );
970 
971  /* check solution for feasibility, and add it to solution store if possible
972  * neither integrality nor feasibility of LP rows has to be checked, because this is already
973  * done in the intshifting heuristic itself and due to the LP resolve
974  */
975  SCIP_CALL( SCIPtrySol(scip, sol, FALSE, FALSE, FALSE, FALSE, &stored) );
976 
977  if( stored )
978  {
979  SCIPdebugMessage("found feasible shifted solution:\n");
980  SCIPdebug( SCIP_CALL( SCIPprintSol(scip, sol, NULL, FALSE) ) );
981  *result = SCIP_FOUNDSOL;
982  }
983  }
984 
985  /* terminate the diving */
986  SCIP_CALL( SCIPendDive(scip) );
987  }
988 
989  /* free memory buffers */
990  SCIPfreeBufferArray(scip, &ndecreases);
991  SCIPfreeBufferArray(scip, &nincreases);
992  SCIPfreeBufferArray(scip, &nfracsinrow);
993  SCIPfreeBufferArray(scip, &violrowpos);
994  SCIPfreeBufferArray(scip, &violrows);
995  SCIPfreeBufferArray(scip, &maxactivities);
996  SCIPfreeBufferArray(scip, &minactivities);
997 
998  return SCIP_OKAY;
999 }
1000 
1001 
1002 /*
1003  * heuristic specific interface methods
1004  */
1005 
1006 /** creates the intshifting heuristic with infeasibility recovering and includes it in SCIP */
1007 SCIP_RETCODE SCIPincludeHeurIntshifting(
1008  SCIP* scip /**< SCIP data structure */
1009  )
1010 {
1011  SCIP_HEUR* heur;
1012 
1013  /* include primal heuristic */
1014  SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
1015  HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS,
1016  HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecIntshifting, NULL) );
1017 
1018  assert(heur != NULL);
1019 
1020  /* set non-NULL pointers to callback methods */
1021  SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyIntshifting) );
1022  SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitIntshifting) );
1023  SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitIntshifting) );
1024  SCIP_CALL( SCIPsetHeurInitsol(scip, heur, heurInitsolIntshifting) );
1025 
1026 
1027  return SCIP_OKAY;
1028 }
1029 
1030