42 #include "scip/misc.h"
45 #define HEUR_NAME "gcgshifting"
46 #define HEUR_DESC "LP rounding heuristic on original variables with infeasibility recovering also using continuous variables"
47 #define HEUR_DISPCHAR 's'
48 #define HEUR_PRIORITY -5000
50 #define HEUR_FREQOFS 0
51 #define HEUR_MAXDEPTH -1
52 #define HEUR_TIMING SCIP_HEURTIMING_AFTERNODE
53 #define HEUR_USESSUBSCIP FALSE
55 #define MAXSHIFTINGS 50
56 #define WEIGHTFACTOR 1.1
57 #define DEFAULT_RANDSEED 31
83 SCIP_Real oldactivity,
92 assert(violrows != NULL);
93 assert(violrowpos != NULL);
94 assert(nviolrows != NULL);
96 lhs = SCIProwGetLhs(row);
97 rhs = SCIProwGetRhs(row);
98 oldviol = (SCIPisFeasLT(scip, oldactivity, lhs) || SCIPisFeasGT(scip, oldactivity, rhs));
99 newviol = (SCIPisFeasLT(scip, newactivity, lhs) || SCIPisFeasGT(scip, newactivity, rhs));
100 if( oldviol != newviol )
104 rowpos = SCIProwGetLPPos(row);
112 violpos = violrowpos[rowpos];
113 assert(0 <= violpos && violpos < *nviolrows);
114 assert(violrows[violpos] == row);
115 violrowpos[rowpos] = -1;
116 if( violpos != *nviolrows-1 )
118 violrows[violpos] = violrows[*nviolrows-1];
119 violrowpos[SCIProwGetLPPos(violrows[violpos])] = violpos;
126 assert(violrowpos[rowpos] == -1);
127 violrows[*nviolrows] = row;
128 violrowpos[rowpos] = *nviolrows;
138 SCIP_Real* activities,
155 assert(activities != NULL);
156 assert(nviolrows != NULL);
157 assert(0 <= *nviolrows && *nviolrows <= nlprows);
159 delta = newsolval - oldsolval;
160 col = SCIPvarGetCol(var);
161 colrows = SCIPcolGetRows(col);
162 colvals = SCIPcolGetVals(col);
163 ncolrows = SCIPcolGetNLPNonz(col);
164 assert(ncolrows == 0 || (colrows != NULL && colvals != NULL));
166 for( r = 0; r < ncolrows; ++r )
172 rowpos = SCIProwGetLPPos(row);
173 assert(-1 <= rowpos && rowpos < nlprows);
175 if( rowpos >= 0 && !SCIProwIsLocal(row) )
177 SCIP_Real oldactivity;
179 assert(SCIProwIsInLP(row));
182 oldactivity = activities[rowpos];
183 if( !SCIPisInfinity(scip, -oldactivity) && !SCIPisInfinity(scip, oldactivity) )
185 SCIP_Real newactivity;
187 newactivity = oldactivity + delta * colvals[r];
188 if( SCIPisInfinity(scip, newactivity) )
189 newactivity = SCIPinfinity(scip);
190 else if( SCIPisInfinity(scip, -newactivity) )
191 newactivity = -SCIPinfinity(scip);
192 activities[rowpos] = newactivity;
195 updateViolations(scip, row, violrows, violrowpos, nviolrows, oldactivity, newactivity);
214 SCIP_Real rowactivity,
216 SCIP_Real* nincreases,
217 SCIP_Real* ndecreases,
218 SCIP_Real increaseweight,
220 SCIP_Real* oldsolval,
227 SCIP_Real activitydelta;
228 SCIP_Real bestshiftscore;
229 SCIP_Real bestdeltaobj;
232 assert(direction == +1 || direction == -1);
233 assert(nincreases != NULL);
234 assert(ndecreases != NULL);
235 assert(shiftvar != NULL);
236 assert(oldsolval != NULL);
237 assert(newsolval != NULL);
240 rowcols = SCIProwGetCols(row);
241 rowvals = SCIProwGetVals(row);
242 nrowcols = SCIProwGetNLPNonz(row);
245 activitydelta = (direction == +1 ? SCIProwGetLhs(row) - rowactivity : SCIProwGetRhs(row) - rowactivity);
246 assert((direction == +1 && SCIPisPositive(scip, activitydelta))
247 || (direction == -1 && SCIPisNegative(scip, activitydelta)));
250 bestshiftscore = SCIP_REAL_MAX;
251 bestdeltaobj = SCIPinfinity(scip);
255 for( c = 0; c < nrowcols; ++c )
261 SCIP_Real shiftscore;
267 var = SCIPcolGetVar(col);
269 assert(!SCIPisZero(scip, val));
270 solval = SCIPgetSolVal(scip, sol, var);
272 isinteger = (SCIPvarGetType(var) == SCIP_VARTYPE_BINARY || SCIPvarGetType(var) == SCIP_VARTYPE_INTEGER);
273 isfrac = isinteger && !SCIPisFeasIntegral(scip, solval);
274 increase = (direction * val > 0.0);
278 shiftscore = increase ? -1.0 / (SCIPvarGetNLocksUp(var) + 1.0) :
279 -1.0 / (SCIPvarGetNLocksDown(var) + 1.0);
283 probindex = SCIPvarGetProbindex(var);
286 shiftscore = ndecreases[probindex]/increaseweight;
288 shiftscore = nincreases[probindex]/increaseweight;
293 if( shiftscore <= bestshiftscore )
301 assert(direction * val < 0.0);
303 shiftval = SCIPfeasFloor(scip, solval);
308 assert(activitydelta/val < 0.0);
309 shiftval = solval + activitydelta/val;
310 assert(shiftval <= solval);
311 if( SCIPvarIsIntegral(var) )
312 shiftval = SCIPfeasFloor(scip, shiftval);
313 lb = SCIPvarGetLbGlobal(var);
314 shiftval = MAX(shiftval, lb);
320 assert(direction * val > 0.0);
322 shiftval = SCIPfeasCeil(scip, solval);
327 assert(activitydelta/val > 0.0);
328 shiftval = solval + activitydelta/val;
329 assert(shiftval >= solval);
330 if( SCIPvarIsIntegral(var) )
331 shiftval = SCIPfeasCeil(scip, shiftval);
332 ub = SCIPvarGetUbGlobal(var);
333 shiftval = MIN(shiftval, ub);
337 if( SCIPisEQ(scip, shiftval, solval) )
340 deltaobj = SCIPvarGetObj(var) * (shiftval - solval);
341 if( shiftscore < bestshiftscore || deltaobj < bestdeltaobj )
343 bestshiftscore = shiftscore;
344 bestdeltaobj = deltaobj;
347 *newsolval = shiftval;
368 SCIP_Real* oldsolval,
372 SCIP_Real bestdeltaobj;
376 assert(shiftvar != NULL);
377 assert(oldsolval != NULL);
378 assert(newsolval != NULL);
382 bestdeltaobj = SCIPinfinity(scip);
384 for( v = 0; v < nlpcands; ++v )
390 assert(SCIPvarGetType(var) == SCIP_VARTYPE_BINARY || SCIPvarGetType(var) == SCIP_VARTYPE_INTEGER);
392 solval = SCIPgetSolVal(scip, sol, var);
393 if( !SCIPisFeasIntegral(scip, solval) )
400 obj = SCIPvarGetObj(var);
403 nlocks = SCIPvarGetNLocksUp(var);
404 if( nlocks >= maxnlocks )
406 shiftval = SCIPfeasFloor(scip, solval);
407 deltaobj = obj * (shiftval - solval);
408 if( (nlocks > maxnlocks || deltaobj < bestdeltaobj) && minobj - obj < SCIPgetCutoffbound(scip) )
411 bestdeltaobj = deltaobj;
414 *newsolval = shiftval;
419 nlocks = SCIPvarGetNLocksDown(var);
420 if( nlocks >= maxnlocks )
422 shiftval = SCIPfeasCeil(scip, solval);
423 deltaobj = obj * (shiftval - solval);
424 if( (nlocks > maxnlocks || deltaobj < bestdeltaobj) && minobj + obj < SCIPgetCutoffbound(scip) )
427 bestdeltaobj = deltaobj;
430 *newsolval = shiftval;
453 col = SCIPvarGetCol(var);
454 rows = SCIPcolGetRows(col);
455 nrows = SCIPcolGetNLPNonz(col);
456 for( r = 0; r < nrows; ++r )
460 rowidx = SCIProwGetLPPos(rows[r]);
461 assert(0 <= rowidx && rowidx < nlprows);
462 nfracsinrow[rowidx] += incval;
463 assert(nfracsinrow[rowidx] >= 0);
474 #define heurCopyGcgshifting NULL
477 #define heurFreeGcgshifting NULL
484 SCIP_HEURDATA* heurdata;
486 assert(strcmp(SCIPheurGetName(heur),
HEUR_NAME) == 0);
487 assert(SCIPheurGetData(heur) == NULL);
490 SCIP_CALL( SCIPallocMemory(scip, &heurdata) );
491 SCIP_CALL( SCIPcreateSol(scip, &heurdata->sol, heur) );
492 heurdata->lastlp = -1;
495 SCIP_CALL( SCIPcreateRandom(scip, &heurdata->randnumgen,
498 SCIPheurSetData(heur, heurdata);
507 SCIP_HEURDATA* heurdata;
509 assert(strcmp(SCIPheurGetName(heur),
HEUR_NAME) == 0);
512 heurdata = SCIPheurGetData(heur);
513 assert(heurdata != NULL);
514 SCIP_CALL( SCIPfreeSol(scip, &heurdata->sol) );
517 SCIPfreeRandom(scip, &heurdata->randnumgen);
519 SCIPfreeMemory(scip, &heurdata);
520 SCIPheurSetData(heur, NULL);
529 SCIP_HEURDATA* heurdata;
531 assert(strcmp(SCIPheurGetName(heur),
HEUR_NAME) == 0);
533 heurdata = SCIPheurGetData(heur);
534 assert(heurdata != NULL);
535 heurdata->lastlp = -1;
542 #define heurExitsolGcgshifting NULL
550 SCIP_HEURDATA* heurdata;
553 SCIP_Real* lpcandssol;
555 SCIP_Real* activities;
557 SCIP_Real* nincreases;
558 SCIP_Real* ndecreases;
561 SCIP_Real increaseweight;
570 int nnonimprovingshifts;
575 SCIP_Longint nsolsfound;
578 assert(strcmp(SCIPheurGetName(heur),
HEUR_NAME) == 0);
579 assert(scip != NULL);
580 assert(result != NULL);
584 assert(masterprob != NULL);
586 *result = SCIP_DIDNOTRUN;
591 if( !SCIPisRelaxSolValid(scip) )
593 SCIPdebugMessage(
"skipping GCG shifting: invalid relaxation solution\n");
598 if( SCIPgetStage(masterprob) > SCIP_STAGE_SOLVING || SCIPgetLPSolstat(masterprob) != SCIP_LPSOLSTAT_OPTIMAL )
602 heurdata = SCIPheurGetData(heur);
603 assert(heurdata != NULL);
606 nlps = SCIPgetNLPs(masterprob);
607 if( nlps == heurdata->lastlp )
609 heurdata->lastlp = nlps;
612 ncalls = SCIPheurGetNCalls(heur);
613 nsolsfound = 10*SCIPheurGetNBestSolsFound(heur) + SCIPheurGetNSolsFound(heur);
614 nnodes = SCIPgetNNodes(scip);
615 if( nnodes % ((ncalls/100)/(nsolsfound+1)+1) != 0 )
619 SCIP_CALL( SCIPgetExternBranchCands(scip, &lpcands, &lpcandssol, NULL, &nlpcands, NULL, NULL, NULL, NULL) );
626 *result = SCIP_DIDNOTFIND;
629 SCIP_CALL( SCIPgetLPRowsData(scip, &lprows, &nlprows) );
631 SCIPdebugMessage(
"executing GCG shifting heuristic: %d LP rows, %d fractionals\n", nlprows, nfrac);;
634 nvars = SCIPgetNVars(scip);
635 SCIP_CALL( SCIPallocBufferArray(scip, &activities, nlprows) );
636 SCIP_CALL( SCIPallocBufferArray(scip, &violrows, nlprows) );
637 SCIP_CALL( SCIPallocBufferArray(scip, &violrowpos, nlprows) );
638 SCIP_CALL( SCIPallocBufferArray(scip, &nfracsinrow, nlprows) );
639 SCIP_CALL( SCIPallocBufferArray(scip, &nincreases, nvars) );
640 SCIP_CALL( SCIPallocBufferArray(scip, &ndecreases, nvars) );
641 BMSclearMemoryArray(nfracsinrow, nlprows);
642 BMSclearMemoryArray(nincreases, nvars);
643 BMSclearMemoryArray(ndecreases, nvars);
649 for( r = 0; r < nlprows; ++r )
654 assert(SCIProwGetLPPos(row) == r);
656 if( !SCIProwIsLocal(row) )
659 if( SCIPisFeasLT(scip, activities[r], SCIProwGetLhs(row) )
660 || SCIPisFeasGT(scip, activities[r], SCIProwGetRhs(row)) )
662 violrows[nviolrows] = row;
663 violrowpos[r] = nviolrows;
672 for( c = 0; c < nlpcands; ++c )
680 SCIP_CALL( SCIPlinkRelaxSol(scip, sol) );
683 minobj = SCIPgetSolTransObj(scip, sol);
687 if( minobj >= SCIPgetCutoffbound(scip) )
689 *result = SCIP_DELAYED;
690 SCIPfreeBufferArray(scip, &ndecreases);
691 SCIPfreeBufferArray(scip, &nincreases);
692 SCIPfreeBufferArray(scip, &nfracsinrow);
693 SCIPfreeBufferArray(scip, &violrowpos);
694 SCIPfreeBufferArray(scip, &violrows);
695 SCIPfreeBufferArray(scip, &activities);
698 for( c = 0; c < nlpcands; ++c )
700 SCIP_Real bestshiftval;
702 obj = SCIPvarGetObj(lpcands[c]);
703 bestshiftval = obj > 0.0 ? SCIPfeasFloor(scip, lpcandssol[c]) : SCIPfeasCeil(scip, lpcandssol[c]);
704 minobj += obj * (bestshiftval - lpcandssol[c]);
708 nnonimprovingshifts = 0;
709 minnviolrows = INT_MAX;
710 increaseweight = 1.0;
711 while( (nfrac > 0 || nviolrows > 0) && nnonimprovingshifts <
MAXSHIFTINGS )
716 SCIP_Bool oldsolvalisfrac;
719 SCIPdebugMessage(
"GCG shifting heuristic: nfrac=%d, nviolrows=%d, obj=%g (best possible obj: %g), cutoff=%g\n",
720 nfrac, nviolrows, SCIPgetSolOrigObj(scip, sol), SCIPretransformObj(scip, minobj),
721 SCIPretransformObj(scip, SCIPgetCutoffbound(scip)));
723 nprevviolrows = nviolrows;
732 if( nviolrows > 0 && (nfrac == 0 || nnonimprovingshifts <
MAXSHIFTINGS-1) )
744 for( rowidx = nviolrows-1; rowidx >= 0; --rowidx )
746 row = violrows[rowidx];
747 rowpos = SCIProwGetLPPos(row);
748 assert(violrowpos[rowpos] == rowidx);
749 if( nfracsinrow[rowpos] > 0 )
755 rowidx = SCIPrandomGetInt(heurdata->randnumgen, 0, nviolrows-1);
756 row = violrows[rowidx];
757 rowpos = SCIProwGetLPPos(row);
758 assert(0 <= rowpos && rowpos < nlprows);
759 assert(violrowpos[rowpos] == rowidx);
760 assert(nfracsinrow[rowpos] == 0);
762 assert(violrowpos[rowpos] == rowidx);
764 SCIPdebugMessage(
"GCG shifting heuristic: try to fix violated row <%s>: %g <= %g <= %g\n",
765 SCIProwGetName(row), SCIProwGetLhs(row), activities[rowpos], SCIProwGetRhs(row));
766 SCIPdebug( SCIP_CALL( SCIPprintRow(scip, row, NULL) ) );
769 assert(SCIPisFeasLT(scip, activities[rowpos], SCIProwGetLhs(row))
770 || SCIPisFeasGT(scip, activities[rowpos], SCIProwGetRhs(row)));
771 direction = SCIPisFeasLT(scip, activities[rowpos], SCIProwGetLhs(row)) ? +1 : -1;
774 SCIP_CALL(
selectShifting(scip, sol, row, activities[rowpos], direction,
775 nincreases, ndecreases, increaseweight, &shiftvar, &oldsolval, &newsolval) );
778 if( shiftvar == NULL && nfrac > 0 )
780 SCIPdebugMessage(
"GCG shifting heuristic: search rounding variable and try to stay feasible\n");
785 if( shiftvar == NULL || SCIPisEQ(scip, oldsolval, newsolval) )
787 SCIPdebugMessage(
"GCG shifting heuristic: -> didn't find a shifting variable\n");
791 SCIPdebugMessage(
"GCG shifting heuristic: -> shift var <%s>[%g,%g], type=%d, oldval=%g, newval=%g, obj=%g\n",
792 SCIPvarGetName(shiftvar), SCIPvarGetLbGlobal(shiftvar), SCIPvarGetUbGlobal(shiftvar), SCIPvarGetType(shiftvar),
793 oldsolval, newsolval, SCIPvarGetObj(shiftvar));
796 SCIP_CALL(
updateActivities(scip, activities, violrows, violrowpos, &nviolrows, nlprows,
797 shiftvar, oldsolval, newsolval) );
798 if( nviolrows >= nprevviolrows )
799 nnonimprovingshifts++;
800 else if( nviolrows < minnviolrows )
802 minnviolrows = nviolrows;
803 nnonimprovingshifts = 0;
807 SCIP_CALL( SCIPsetSolVal(scip, sol, shiftvar, newsolval) );
810 oldsolvalisfrac = !SCIPisFeasIntegral(scip, oldsolval)
811 && (SCIPvarGetType(shiftvar) == SCIP_VARTYPE_BINARY || SCIPvarGetType(shiftvar) == SCIP_VARTYPE_INTEGER);
812 obj = SCIPvarGetObj(shiftvar);
813 if( (SCIPvarGetType(shiftvar) == SCIP_VARTYPE_BINARY || SCIPvarGetType(shiftvar) == SCIP_VARTYPE_INTEGER )
816 assert(SCIPisFeasIntegral(scip, newsolval));
818 nnonimprovingshifts = 0;
819 minnviolrows = INT_MAX;
823 if( obj > 0.0 && newsolval > oldsolval )
825 else if( obj < 0.0 && newsolval < oldsolval )
831 minobj += obj * (newsolval - oldsolval);
835 if( !oldsolvalisfrac )
839 probindex = SCIPvarGetProbindex(shiftvar);
840 assert(0 <= probindex && probindex < nvars);
842 if( newsolval < oldsolval )
843 ndecreases[probindex] += increaseweight;
845 nincreases[probindex] += increaseweight;
846 if( increaseweight >= 1e+09 )
850 for( i = 0; i < nvars; ++i )
852 nincreases[i] /= increaseweight;
853 ndecreases[i] /= increaseweight;
855 increaseweight = 1.0;
859 SCIPdebugMessage(
"gcg shifting heuristic: -> nfrac=%d, nviolrows=%d, obj=%g (best possible obj: %g)\n",
860 nfrac, nviolrows, SCIPgetSolOrigObj(scip, sol), SCIPretransformObj(scip, minobj));
864 if( nfrac == 0 && nviolrows == 0 )
873 SCIP_CALL( SCIPtrySol(scip, sol, FALSE, FALSE, FALSE, FALSE, TRUE, &stored) );
877 SCIPdebugMessage(
"found feasible shifted solution:\n");
878 SCIPdebug(SCIPprintSol(scip, sol, NULL, FALSE));
879 *result = SCIP_FOUNDSOL;
884 SCIPfreeBufferArray(scip, &ndecreases);
885 SCIPfreeBufferArray(scip, &nincreases);
886 SCIPfreeBufferArray(scip, &nfracsinrow);
887 SCIPfreeBufferArray(scip, &violrowpos);
888 SCIPfreeBufferArray(scip, &violrows);
889 SCIPfreeBufferArray(scip, &activities);