43 #include "scip/scipdefplugins.h"
44 #include "scip/cons_linear.h"
45 #include "scip/misc.h"
48 #define HEUR_NAME "gcgfeaspump"
49 #define HEUR_DESC "objective feasibility pump 2.0"
50 #define HEUR_DISPCHAR 'F'
51 #define HEUR_PRIORITY -1000000
53 #define HEUR_FREQOFS 0
54 #define HEUR_MAXDEPTH -1
55 #define HEUR_TIMING SCIP_HEURTIMING_AFTERPLUNGE
56 #define HEUR_USESSUBSCIP TRUE
58 #define DEFAULT_MAXLPITERQUOT 0.01
59 #define DEFAULT_MAXLPITEROFS 1000
60 #define DEFAULT_MAXSOLS 10
62 #define DEFAULT_MAXLOOPS 10000
63 #define DEFAULT_MAXSTALLLOOPS 10
64 #define DEFAULT_MINFLIPS 10
65 #define DEFAULT_CYCLELENGTH 3
66 #define DEFAULT_PERTURBFREQ 100
67 #define DEFAULT_OBJFACTOR 1.0
69 #define DEFAULT_ALPHADIFF 1.0
70 #define DEFAULT_USEFP20 FALSE
71 #define DEFAULT_PERTSOLFOUND TRUE
72 #define DEFAULT_STAGE3 FALSE
73 #define DEFAULT_NEIGHBORHOODSIZE 18
74 #define DEFAULT_COPYCUTS TRUE
78 #define MINLPITER 5000
80 #define DEFAULT_RANDSEED 13
119 SCIP_HASHMAP** varmapfw,
125 SCIP_VAR** tmpsubvars;
131 SCIP_CALL( SCIPcreate(divingscip) );
134 SCIP_CALL( SCIPhashmapCreate(varmapfw, SCIPblkmem(*divingscip), SCIPgetNVars(scip)) );
138 SCIP_CALL( SCIPcopy(scip, *divingscip, *varmapfw, NULL,
"gcgfeaspump", FALSE, FALSE, FALSE, TRUE, success) );
143 SCIP_CALL( SCIPcopyCuts(scip, *divingscip, *varmapfw, NULL, FALSE, NULL) );
147 SCIP_CALL( SCIPgetVarsData(*divingscip, &tmpsubvars, &nsubvars, NULL, NULL, NULL, NULL) );
148 SCIP_CALL( SCIPduplicateBufferArray(*divingscip, &subvars, tmpsubvars, nsubvars) );
149 for( i = 0; i < nsubvars; ++i )
151 SCIP_Bool infeasible;
152 SCIP_CALL( SCIPchgVarType(*divingscip, subvars[i], SCIP_VARTYPE_CONTINUOUS, &infeasible) );
155 SCIPfreeBufferArray(*divingscip, &subvars);
158 SCIP_CALL( SCIPsetBoolParam(*divingscip,
"misc/catchctrlc", FALSE) );
161 SCIP_CALL( SCIPsetIntParam(*divingscip,
"display/verblevel", 0) );
164 SCIP_CALL( SCIPsetSeparating(*divingscip, SCIP_PARAMSETTING_OFF, TRUE) );
167 SCIP_CALL( SCIPsetPresolving(*divingscip, SCIP_PARAMSETTING_FAST, TRUE) );
170 SCIP_CALL( SCIPsetHeuristics(*divingscip, SCIP_PARAMSETTING_OFF, TRUE) );
173 if( !SCIPisParamFixed(*divingscip,
"conflict/enable") )
175 SCIP_CALL( SCIPsetBoolParam(*divingscip,
"conflict/enable", FALSE) );
179 SCIP_CALL( SCIPsetLongintParam(*divingscip,
"limits/nodes", 1LL) );
189 SCIP_HASHMAP* varmapfw,
198 assert(scip != NULL);
199 assert(divingscip != NULL);
200 assert(varmapfw != NULL);
201 assert(lpsol != NULL);
204 subsol = SCIPgetBestSol(divingscip);
205 assert(subsol != NULL);
208 SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
211 for( i = 0; i < nvars; ++i )
215 subvar = (SCIP_VAR*) SCIPhashmapGetImage(varmapfw, vars[i]);
216 assert(subvar != NULL);
218 SCIP_CALL( SCIPsetSolVal(scip, lpsol, vars[i], SCIPgetSolVal(divingscip, subsol, subvar)) );
237 assert(scip != NULL);
238 assert(lpsol != NULL);
239 assert(nfracs != NULL);
242 SCIP_CALL( SCIPgetVarsData(scip, &vars, NULL, &nbinvars, &nintvars, NULL, NULL) );
246 for( i = 0; i < nbinvars + nintvars; ++i )
247 if( !SCIPisFeasIntegral(scip, SCIPgetSolVal(scip, lpsol, vars[i])) )
258 SCIP_HASHMAP** varmapfw,
264 SCIP_CALL( SCIPcreate(probingscip) );
267 SCIP_CALL( SCIPhashmapCreate(varmapfw, SCIPblkmem(*probingscip), SCIPgetNVars(scip)) );
271 SCIP_CALL( SCIPcopy(scip, *probingscip, *varmapfw, NULL,
"gcgfeaspump_probing", FALSE, FALSE, FALSE, TRUE, success) );
276 SCIP_CALL( SCIPcopyCuts(scip, *probingscip, *varmapfw, NULL, FALSE, NULL) );
285 SCIP_VAR** mostfracvars,
286 SCIP_Real* mostfracvals,
295 assert(mostfracvars != NULL);
296 assert(mostfracvals != NULL);
297 assert(nflipcands != NULL);
304 if( *nflipcands >= maxnflipcands )
306 if( frac <= mostfracvals[*nflipcands-1] )
313 for( i = *nflipcands; i > 0 && mostfracvals[i-1] < frac; i-- )
315 mostfracvars[i] = mostfracvars[i-1];
316 mostfracvals[i] = mostfracvals[i-1];
318 assert(0 <= i && i <= *nflipcands && *nflipcands < maxnflipcands);
321 mostfracvars[i] = var;
322 mostfracvals[i] = frac;
333 SCIP_HASHMAP* varmapfw,
334 SCIP_HEURDATA* heurdata,
335 SCIP_VAR** mostfracvars,
343 for( i = 0; i < nflipcands; i++ )
349 SCIP_Real newobjcoeff;
350 SCIP_Real orgobjcoeff;
352 var = mostfracvars[i];
353 solval = SCIPvarGetLPSol(var);
354 orgobjcoeff = SCIPvarGetObj(var);
355 frac = SCIPfeasFrac(scip, solval);
359 newobjcoeff = (1.0 - alpha) + alpha * orgobjcoeff;
360 solval = SCIPfeasFloor(scip, solval);
364 newobjcoeff = - (1.0 - alpha) + alpha * orgobjcoeff;
365 solval = SCIPfeasCeil(scip, solval);
368 SCIP_CALL( SCIPsetSolVal(scip, heurdata->roundedsol, var, solval) );
369 divingvar = (SCIP_VAR*) SCIPhashmapGetImage(varmapfw, var);
370 SCIP_CALL( SCIPchgVarObj(divingscip, divingvar, newobjcoeff) );
382 SCIP_HASHMAP* varmapfw,
383 SCIP_HEURDATA* heurdata,
392 for( i = 0; i < nbinandintvars; i++ )
397 SCIP_Real orgobjcoeff;
402 solval = SCIPvarGetLPSol(var);
403 orgobjcoeff = SCIPvarGetObj(var);
404 frac = SCIPfeasFrac(scip, solval);
405 flipprob = -0.3 + SCIPrandomGetReal(heurdata->randnumgen, 0.0, 1.0);
408 if( MIN(frac, 1.0-frac) + MAX(flipprob, 0.0) > 0.5 )
410 SCIP_Real newobjcoeff;
415 newobjcoeff = (1.0 - alpha) + alpha * orgobjcoeff;
416 solval = SCIPfeasFloor(scip, solval);
420 newobjcoeff = - (1.0 - alpha) + alpha * orgobjcoeff;
421 solval = SCIPfeasCeil(scip, solval);
423 SCIP_CALL( SCIPsetSolVal(scip, heurdata->roundedsol, var, solval) );
424 divingvar = (SCIP_VAR*) SCIPhashmapGetImage(varmapfw, var);
425 SCIP_CALL( SCIPchgVarObj(divingscip, divingvar, newobjcoeff) );
437 SCIP_HASHMAP* varmapfw,
439 SCIP_Real neighborhoodsize
451 char consname[SCIP_MAXSTRLEN];
453 (void) SCIPsnprintf(consname, SCIP_MAXSTRLEN,
"%s_localbranchcons", SCIPgetProbName(scip));
456 SCIP_CALL( SCIPgetVarsData(scip, &vars, NULL, &nbinvars, NULL, NULL, NULL) );
458 SCIP_CALL( SCIPallocBufferArray(scip, &consvars, nbinvars) );
459 SCIP_CALL( SCIPallocBufferArray(scip, &consvals, nbinvars) );
463 rhs = neighborhoodsize;
466 for( i = 0; i < nbinvars; i++ )
470 solval = SCIPgetSolVal(scip, bestsol, vars[i]);
471 assert( SCIPisFeasIntegral(scip,solval) );
474 if( SCIPisFeasEQ(scip,solval,1.0) )
482 consvars[i] = (SCIP_VAR*) SCIPhashmapGetImage(varmapfw, vars[i]);
483 SCIP_CALL( SCIPchgVarObj(probingscip, consvars[i], consvals[i]) );
484 assert( SCIPvarGetType(consvars[i]) == SCIP_VARTYPE_BINARY );
488 SCIP_CALL( SCIPcreateConsLinear(probingscip, &cons, consname, nbinvars, consvars, consvals,
489 lhs, rhs, FALSE, FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
490 SCIP_CALL( SCIPaddCons(probingscip, cons) );
491 SCIP_CALL( SCIPreleaseCons(probingscip, &cons) );
494 SCIPfreeBufferArray(scip, &consvals);
495 SCIPfreeBufferArray(scip, &consvars);
505 SCIP_HASHMAP* varmapfw,
515 SCIP_Real* subsolvals;
519 assert(scip != NULL);
520 assert(subscip != NULL);
523 SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
528 assert(nvars <= SCIPgetNOrigVars(subscip));
531 SCIP_CALL( SCIPallocBufferArray(scip, &subvars, nvars) );
532 for( i = 0; i < nvars; i++ )
533 subvars[i] = (SCIP_VAR*) SCIPhashmapGetImage(varmapfw, vars[i]);
535 SCIP_CALL( SCIPallocBufferArray(scip, &subsolvals, nvars) );
537 nsubsols = SCIPgetNSols(subscip);
538 subsols = SCIPgetSols(subscip);
541 for( i = 0; i < nsubsols && !(*success); ++i )
544 SCIP_CALL( SCIPgetSolVals(subscip, subsols[i], nvars, subvars, subsolvals) );
547 SCIP_CALL( SCIPcreateSol(scip, &newsol, heur) );
548 SCIP_CALL( SCIPsetSolVals(scip, newsol, nvars, vars, subsolvals) );
551 SCIP_CALL( SCIPtrySolFree(scip, &newsol, FALSE, FALSE, TRUE, TRUE, TRUE, success) );
554 SCIPfreeBufferArray(scip, &subvars);
555 SCIPfreeBufferArray(scip, &subsolvals);
564 SCIP_HEURDATA* heurdata;
566 assert(heur != NULL);
567 assert(strcmp(SCIPheurGetName(heur),
HEUR_NAME) == 0);
568 assert(scip != NULL);
571 heurdata = SCIPheurGetData(heur);
572 assert(heurdata != NULL);
573 SCIPfreeMemory(scip, &heurdata);
574 SCIPheurSetData(heur, NULL);
584 SCIP_HEURDATA* heurdata;
586 assert(heur != NULL);
587 assert(strcmp(SCIPheurGetName(heur),
HEUR_NAME) == 0);
590 heurdata = SCIPheurGetData(heur);
591 assert(heurdata != NULL);
594 SCIP_CALL( SCIPcreateSol(scip, &heurdata->sol, heur) );
595 SCIP_CALL( SCIPcreateSol(scip, &heurdata->roundedsol, heur) );
598 heurdata->nlpiterations = 0;
599 heurdata->nsuccess = 0;
602 SCIP_CALL( SCIPcreateRandom(scip, &heurdata->randnumgen,
612 SCIP_HEURDATA* heurdata;
614 assert(heur != NULL);
615 assert(strcmp(SCIPheurGetName(heur),
HEUR_NAME) == 0);
618 heurdata = SCIPheurGetData(heur);
619 assert(heurdata != NULL);
622 SCIPfreeRandom(scip, &heurdata->randnumgen);
625 SCIP_CALL( SCIPfreeSol(scip, &heurdata->sol) );
626 SCIP_CALL( SCIPfreeSol(scip, &heurdata->roundedsol) );
635 SCIP_Longint maxnlpiterations,
636 SCIP_Longint nsolsfound,
640 if( nstallloops <= 1 )
642 if( nsolsfound == 0 )
643 return 4*maxnlpiterations;
645 return 2*maxnlpiterations;
648 return maxnlpiterations;
656 SCIP_HEURDATA* heurdata;
658 SCIP_SOL** lastroundedsols;
659 SCIP_SOL* closestsol;
660 SCIP_Real* lastalphas;
664 SCIP_HASHMAP* varmapfwdive;
665 SCIP_HASHMAP* varmapfw;
669 SCIP_VAR** pseudocands;
670 SCIP_VAR** tmppseudocands;
671 SCIP_VAR** mostfracvars;
674 SCIP_Real* mostfracvals;
675 SCIP_Real newobjcoeff;
676 SCIP_Real orgobjcoeff;
683 SCIP_Real scalingfactor;
684 SCIP_Real mindistance;
686 SCIP_Longint nlpiterations;
687 SCIP_Longint maxnlpiterations;
688 SCIP_Longint nsolsfound;
690 SCIP_Longint nbestsolsfound;
710 SCIP_RETCODE retcode;
712 assert(heur != NULL);
713 assert(strcmp(SCIPheurGetName(heur),
HEUR_NAME) == 0);
714 assert(scip != NULL);
715 assert(result != NULL);
719 assert(masterprob != NULL);
721 *result = SCIP_DELAYED;
726 if( !SCIPisRelaxSolValid(scip) )
730 if( SCIPgetStage(masterprob) > SCIP_STAGE_SOLVING || SCIPgetLPSolstat(masterprob) != SCIP_LPSOLSTAT_OPTIMAL )
733 *result = SCIP_DIDNOTRUN;
736 if( SCIPgetDepth(scip) == 0 && SCIPheurGetNCalls(heur) > 0 )
740 heurdata = SCIPheurGetData(heur);
741 assert(heurdata != NULL);
744 if( heurdata->maxsols >= 0 && SCIPgetNSolsFound(scip) > heurdata->maxsols && SCIPgetNPricers(scip) == 0 )
748 SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, &nbinvars, &nintvars, NULL, NULL) );
749 nfracs = SCIPgetNExternBranchCands(scip);
750 assert(0 <= nfracs && nfracs <= nbinvars + nintvars);
755 nlpiterations = SCIPgetNLPIterations(scip);
756 ncalls = SCIPheurGetNCalls(heur);
757 nsolsfound = 10*SCIPheurGetNBestSolsFound(heur) + heurdata->nsuccess;
758 maxnlpiterations = (SCIP_Longint)((1.0 + 10.0*(nsolsfound+1.0)/(ncalls+1.0)) * heurdata->maxlpiterquot * nlpiterations);
759 maxnlpiterations += heurdata->maxlpiterofs;
762 if( heurdata->nlpiterations >= maxnlpiterations )
766 if( SCIPheurGetNCalls(heur) == 0 && SCIPgetNSolsFound(scip) == 0 && SCIPgetDepth(scip) == 0 )
767 maxnlpiterations += nlpiterations;
770 maxnlpiterations = MAX(maxnlpiterations, heurdata->nlpiterations +
MINLPITER);
773 maxflips = 3*heurdata->minflips;
774 maxloops = (heurdata->maxloops == -1 ? INT_MAX : heurdata->maxloops);
775 maxstallloops = (heurdata->maxstallloops == -1 ? INT_MAX : heurdata->maxstallloops);
777 SCIPdebugMessage(
"executing GCG feasibility pump heuristic, nlpiters=%"SCIP_LONGINT_FORMAT
", maxnlpit:%"SCIP_LONGINT_FORMAT
", maxflips:%d \n",
778 nlpiterations, maxnlpiterations, maxflips);
780 *result = SCIP_DIDNOTFIND;
785 if( heurdata->usefp20 )
787 SCIP_CALL(
setupProbingSCIP(scip, &probingscip, &varmapfw, heurdata->copycuts, &success) );
791 if( SCIPisParamFixed(probingscip,
"heuristics/"HEUR_NAME"/freq") )
793 SCIPwarningMessage(scip,
"unfixing parameter heuristics/"HEUR_NAME"/freq in probingscip of "HEUR_NAME" heuristic to avoid recursive calls\n");
794 SCIP_CALL( SCIPunfixParam(probingscip,
"heuristics/"HEUR_NAME"/freq") );
796 SCIP_CALL( SCIPsetIntParam(probingscip,
"heuristics/"HEUR_NAME"/freq", -1) );
799 SCIP_CALL( SCIPsetBoolParam(probingscip,
"misc/catchctrlc", FALSE) );
803 SCIP_CALL( SCIPsetIntParam(probingscip,
"display/verblevel", 0) );
807 SCIP_CALL( SCIPsetLongintParam(probingscip,
"limits/nodes", 1LL) );
808 if( SCIPisParamFixed(probingscip,
"lp/solvefreq") )
810 SCIPwarningMessage(scip,
"unfixing parameter lp/solvefreq in probingscip of "HEUR_NAME" heuristic to avoid recursive calls\n");
811 SCIP_CALL( SCIPunfixParam(probingscip,
"lp/solvefreq") );
813 SCIP_CALL( SCIPsetIntParam(probingscip,
"lp/solvefreq", -1) );
816 SCIP_CALL( SCIPsetPresolving(probingscip, SCIP_PARAMSETTING_FAST, TRUE) );
817 retcode = SCIPsolve(probingscip);
821 if( retcode != SCIP_OKAY )
824 SCIP_CALL( retcode );
826 SCIPwarningMessage(scip,
"Error while solving subproblem in feaspump heuristic; sub-SCIP terminated with code <%d>\n", retcode);
829 SCIPhashmapFree(&varmapfw);
830 SCIP_CALL( SCIPfree(&probingscip) );
834 if( SCIPgetStage(probingscip) != SCIP_STAGE_SOLVING )
836 SCIP_STATUS probingstatus = SCIPgetStatus(probingscip);
838 if( probingstatus == SCIP_STATUS_OPTIMAL )
840 assert( SCIPgetNSols(probingscip) > 0 );
841 SCIP_CALL(
createNewSols(scip, probingscip, varmapfw, heur, &success) );
843 *result = SCIP_FOUNDSOL;
847 SCIPhashmapFree(&varmapfw);
848 SCIP_CALL( SCIPfree(&probingscip) );
851 SCIP_CALL( SCIPsetLongintParam(probingscip,
"limits/nodes", 2LL) );
854 SCIP_CALL( SCIPstartProbing(probingscip) );
855 SCIP_CALL( SCIPnewProbingNode(probingscip) );
857 SCIPdebugMessage(
"successfully copied SCIP instance -> feasibility pump 2.0 can be used.\n");
862 SCIP_CALL( SCIPallocBufferArray(scip, &mostfracvars, maxflips) );
863 SCIP_CALL( SCIPallocBufferArray(scip, &mostfracvals, maxflips) );
864 SCIP_CALL( SCIPallocBufferArray(scip, &lastroundedsols, heurdata->cyclelength) );
865 SCIP_CALL( SCIPallocBufferArray(scip, &lastalphas, heurdata->cyclelength) );
866 SCIP_CALL( SCIPallocBufferArray(scip, &cycles, heurdata->cyclelength) );
868 for( j = 0; j < heurdata->cyclelength; j++ )
870 SCIP_CALL( SCIPcreateSol(scip, &lastroundedsols[j], heur) );
874 if( heurdata->stage3 )
876 SCIP_CALL( SCIPcreateSol(scip, &closestsol, heur) );
880 SCIP_CALL(
setupDivingSCIP(scip, &divingscip, &varmapfwdive, heurdata->copycuts, &success) );
883 nsolsfound = SCIPgetNBestSolsFound(scip);
884 if( heurdata->objfactor == 1.0 )
885 objfactor = MIN(1.0 - 0.1 / (SCIP_Real)(1 + nsolsfound), 0.999);
887 objfactor = heurdata->objfactor;
890 objnorm = SCIPgetObjNorm(scip);
891 objnorm = MAX(objnorm, 1.0);
892 scalingfactor = SQRT((SCIP_Real)(nbinvars + nintvars)) / objnorm;
898 nbestsolsfound = SCIPgetNBestSolsFound(scip);
899 bestnfracs = INT_MAX;
900 mindistance = SCIPinfinity(scip);
902 SCIP_CALL( SCIPlinkRelaxSol(scip, heurdata->sol) );
903 SCIP_CALL( SCIPlinkRelaxSol(scip, heurdata->roundedsol) );
908 && nloops < maxloops && nstallloops < maxstallloops
909 && !SCIPisStopped(scip) )
912 SCIP_Real* pseudocandsfrac;
914 SCIP_Longint nlpiterationsleft;
915 SCIP_Longint iterlimit;
921 SCIPdebugMessage(
"feasibility pump loop %d: %d fractional variables (alpha: %.4f, stall: %d/%d)\n",
922 nloops, nfracs, alpha, nstallloops, maxstallloops);
926 SCIP_CALL( SCIProundSol(scip, heurdata->sol, &success) );
931 SCIP_CALL( SCIPtrySol(scip, heurdata->sol, FALSE, FALSE, FALSE, FALSE, FALSE, &success) );
933 *result = SCIP_FOUNDSOL;
937 maxnflipcands = SCIPrandomGetInt(heurdata->randnumgen, MIN(nfracs/2+1, heurdata->minflips), MIN(nfracs, maxflips));
941 SCIP_CALL( SCIPgetPseudoBranchCands(scip, &tmppseudocands, &npseudocands, NULL) );
942 SCIP_CALL( SCIPduplicateBufferArray(scip, &pseudocands, tmppseudocands, npseudocands) );
945 if( heurdata->usefp20 )
947 SCIP_CALL( SCIPallocBufferArray(scip, &pseudocandsfrac, npseudocands) );
949 for( i = 0; i < npseudocands; i++ )
951 frac = SCIPfeasFrac(scip, SCIPgetSolVal(scip, heurdata->roundedsol, pseudocands[i]));
952 pseudocandsfrac[i] = MIN(frac, 1.0-frac);
953 if( SCIPvarGetType(pseudocands[i]) == SCIP_VARTYPE_BINARY )
954 pseudocandsfrac[i] -= 10.0;
956 SCIPsortRealPtr(pseudocandsfrac, (
void**)pseudocands, npseudocands);
957 SCIPfreeBufferArray(scip, &pseudocandsfrac);
959 SCIPdebugMessage(
"iteratively fix and propagate variables\n");
962 for( i = 0; i < npseudocands; i++ )
965 SCIP_Bool infeasible;
966 SCIP_Longint ndomreds;
968 var = pseudocands[i];
969 orgobjcoeff = SCIPvarGetObj(var);
972 solval = SCIPgetSolVal(scip, heurdata->roundedsol, var);
973 frac = SCIPfeasFrac(scip, solval);
975 solval = SCIPfloor(scip, solval+0.5);
978 if( heurdata->usefp20 )
980 SCIP_VAR* probingvar;
984 probingvar = (SCIP_VAR*) SCIPhashmapGetImage(varmapfw, var);
985 lb = SCIPvarGetLbLocal(probingvar);
986 ub = SCIPvarGetUbLocal(probingvar);
988 solval = MAX(solval, lb);
989 solval = MIN(solval, ub);
992 if( !SCIPisFeasEQ(probingscip, lb, ub) )
994 assert(SCIPisFeasLE(probingscip, lb, ub));
995 SCIP_CALL( SCIPnewProbingNode(probingscip) );
997 SCIP_CALL( SCIPfixVarProbing(probingscip, probingvar, solval) );
998 SCIPdebugMessage(
"try to fix variable <%s> (domain [%f,%f] to %f\n",SCIPvarGetName(probingvar), lb, ub,
1000 SCIP_CALL( SCIPpropagateProbing(probingscip, 3, &infeasible, &ndomreds) );
1001 SCIPdebugMessage(
" -> reduced %"SCIP_LONGINT_FORMAT
" domains\n", ndomreds);
1005 SCIPdebugMessage(
" -> infeasible!\n");
1006 SCIP_CALL( SCIPbacktrackProbing(probingscip, SCIPgetProbingDepth(probingscip)-1) );
1011 SCIPdebugMessage(
"variable <%s> is already fixed to %f\n",SCIPvarGetName(probingvar), solval);
1015 assert(SCIPisIntegral(scip,solval));
1016 SCIP_CALL( SCIPsetSolVal(scip, heurdata->roundedsol, var, solval) );
1019 if( SCIPisFeasZero(scip, frac) )
1025 lb = SCIPvarGetLbLocal(var);
1026 ub = SCIPvarGetUbLocal(var);
1027 if( SCIPisFeasEQ(scip, solval, lb) )
1028 newobjcoeff = (1.0 - alpha)/scalingfactor + alpha * orgobjcoeff;
1029 else if( SCIPisFeasEQ(scip, solval, ub) )
1030 newobjcoeff = - (1.0 - alpha)/scalingfactor + alpha * orgobjcoeff;
1032 newobjcoeff = alpha * orgobjcoeff;
1037 insertFlipCand(mostfracvars, mostfracvals, &nflipcands, maxnflipcands, var, frac);
1040 newobjcoeff = - (1.0 - alpha)/scalingfactor + alpha * orgobjcoeff;
1042 newobjcoeff = (1.0 - alpha)/scalingfactor + alpha * orgobjcoeff;
1046 divingvar = (SCIP_VAR*) SCIPhashmapGetImage(varmapfwdive, var);
1047 SCIP_CALL( SCIPchgVarObj(divingscip, divingvar, newobjcoeff) );
1050 if( heurdata->usefp20 )
1052 SCIP_CALL( SCIPbacktrackProbing(probingscip, 1) );
1056 for( i = nbinvars+nintvars; i < nvars; i++ )
1058 SCIP_VAR* divingvar;
1059 divingvar = (SCIP_VAR*) SCIPhashmapGetImage(varmapfwdive, vars[i]);
1060 SCIP_CALL( SCIPchgVarObj(divingscip, divingvar, alpha * SCIPvarGetObj(vars[i])) );
1063 SCIPfreeBufferArray(scip, &pseudocands);
1066 minimum = MIN(heurdata->cyclelength, nloops-1);
1067 for( j = 0; j < heurdata->cyclelength; j++ )
1068 cycles[j] = (nloops > j+1) && (REALABS(lastalphas[j] - alpha) < heurdata->alphadiff);
1071 for( i = 0; i < nbinvars+nintvars; i++ )
1073 solval = SCIPgetSolVal(scip, heurdata->roundedsol, vars[i]);
1076 for( j = 0; j < minimum; j++ )
1078 oldsolval = SCIPgetSolVal(scip, lastroundedsols[j], vars[i]);
1079 cycles[j] = cycles[j] && SCIPisFeasEQ(scip, solval, oldsolval);
1086 assert(heurdata->perturbfreq > 0);
1087 if( nloops % heurdata->perturbfreq == 0 || (heurdata->pertsolfound && SCIPgetNBestSolsFound(scip) > nbestsolsfound) )
1089 SCIPdebugMessage(
" -> random perturbation\n");
1090 SCIP_CALL(
handleCycle(scip, divingscip, varmapfwdive, heurdata, vars, nintvars+nbinvars, alpha) );
1091 nbestsolsfound = SCIPgetNBestSolsFound(scip);
1095 minimum = MIN(heurdata->cyclelength, nloops-1);
1097 for( j = 0; j < minimum; j++ )
1105 SCIPdebugMessage(
" -> avoiding 1-cycle: flipping %d candidates\n", nflipcands);
1106 SCIP_CALL(
handle1Cycle(scip, divingscip, varmapfwdive, heurdata, mostfracvars, nflipcands, alpha) );
1110 SCIPdebugMessage(
" -> avoiding %d-cycle by random flip\n", j+1);
1111 SCIP_CALL(
handleCycle(scip, divingscip, varmapfwdive, heurdata, vars, nintvars+nbinvars, alpha) );
1120 iterlimit = MAX((
int)nlpiterationsleft,
MINLPITER);
1121 SCIP_CALL( SCIPsetLongintParam(divingscip,
"lp/iterlim", iterlimit) );
1122 SCIPdebugMessage(
" -> solve LP with iteration limit %"SCIP_LONGINT_FORMAT
"\n", iterlimit);
1124 if( heurdata->stage3 )
1126 SCIP_CALL( SCIPunlinkSol(scip, heurdata->roundedsol) );
1130 retcode = SCIPsolve(divingscip);
1135 if( retcode != SCIP_OKAY )
1138 SCIP_CALL( retcode );
1140 SCIPwarningMessage(scip,
"Error while solving subproblem in Feasibility Pump heuristic; sub-SCIP terminated with code <%d>\n",retcode);;
1141 SCIPwarningMessage(scip,
"This does not affect the remaining solution procedure --> continue\n");
1145 heurdata->nlpiterations += SCIPgetNLPIterations(divingscip);
1146 SCIPdebugMessage(
" -> number of iterations: %"SCIP_LONGINT_FORMAT
"/%"SCIP_LONGINT_FORMAT
", status=%d\n",
1147 heurdata->nlpiterations,
adjustedMaxNLPIterations(maxnlpiterations, nsolsfound, nstallloops), SCIPgetStatus(divingscip));
1150 if( SCIPgetStage(divingscip) != SCIP_STAGE_SOLVED || SCIPgetBestSol(divingscip) == NULL )
1152 SCIPdebugMessage(
" -> solstat is %d\n", SCIPgetStatus(divingscip));
1153 SCIPdebugMessage(
" -> diving LP was not solved to optimality --> abort heuristic\n");
1158 SCIP_CALL(
getDivingLPSol(scip, divingscip, varmapfwdive, heurdata->sol) );
1160 if( heurdata->stage3 )
1164 assert(closestsol != NULL);
1168 for( i = 0; i < nbinvars+nintvars; i++ )
1170 SCIP_Real roundedval;
1173 roundedval = SCIPgetSolVal(scip, heurdata->roundedsol, vars[i]);
1174 lpval = SCIPgetSolVal(scip, heurdata->sol, vars[i]);
1175 distance += REALABS(roundedval - lpval);
1179 if( SCIPisLT(scip, distance, mindistance) )
1181 for( i = 0; i < nbinvars+nintvars; i++ )
1183 assert(SCIPisIntegral(scip,SCIPgetSolVal(scip, heurdata->roundedsol, vars[i])));
1184 SCIP_CALL( SCIPsetSolVal(scip, closestsol, vars[i], SCIPgetSolVal(scip, heurdata->roundedsol, vars[i])) );
1186 mindistance = distance;
1191 tmpsol = lastroundedsols[heurdata->cyclelength-1];
1192 for( j = heurdata->cyclelength-1; j > 0; j-- )
1194 lastroundedsols[j] = lastroundedsols[j-1];
1195 lastalphas[j] = lastalphas[j-1];
1197 lastroundedsols[0] = heurdata->roundedsol;
1198 lastalphas[0] = alpha;
1199 heurdata->roundedsol = tmpsol;
1201 SCIP_CALL(
getDivingLPSol(scip, divingscip, varmapfwdive, heurdata->roundedsol) );
1204 SCIP_CALL(
getNFracs(scip, heurdata->sol, &nfracs) );
1205 if( nfracs < bestnfracs )
1207 bestnfracs = nfracs;
1214 SCIP_CALL( SCIPfreeTransform(divingscip) );
1216 SCIPdebugMessage(
" -> loop finished: %d fractional variables (stall: %d/%d, iterations: %"SCIP_LONGINT_FORMAT
"/%"SCIP_LONGINT_FORMAT
")\n",
1217 nfracs, nstallloops, maxstallloops, heurdata->nlpiterations,
adjustedMaxNLPIterations(maxnlpiterations, nsolsfound, nstallloops));
1225 SCIP_CALL( SCIPtrySol(scip, heurdata->sol, FALSE, FALSE, FALSE, FALSE, FALSE, &success) );
1227 *result = SCIP_FOUNDSOL;
1231 SCIPhashmapFree(&varmapfwdive);
1232 SCIP_CALL( SCIPfree(&divingscip) );
1235 if( heurdata->usefp20 )
1237 SCIP_CALL( SCIPendProbing(probingscip) );
1242 if( heurdata->stage3 && (*result != SCIP_FOUNDSOL) && SCIPisLE(scip, mindistance, (SCIP_Real) heurdata->neighborhoodsize) )
1245 SCIP_Real timelimit;
1246 SCIP_Real memorylimit;
1248 assert(closestsol != NULL);
1249 assert(!SCIPisInfinity(scip, mindistance) || nloops == 0);
1252 if( heurdata->usefp20 )
1254 assert(probingscip != NULL);
1255 SCIP_CALL( SCIPfreeTransform(probingscip) );
1259 assert(probingscip == NULL);
1260 SCIP_CALL(
setupProbingSCIP(scip, &probingscip, &varmapfw, heurdata->copycuts, &success) );
1264 SCIP_CALL( SCIPgetRealParam(scip,
"limits/time", &timelimit) );
1265 if( !SCIPisInfinity(scip, timelimit) )
1266 timelimit -= SCIPgetSolvingTime(scip);
1267 SCIP_CALL( SCIPgetRealParam(scip,
"limits/memory", &memorylimit) );
1270 if( !SCIPisInfinity(scip, memorylimit) )
1272 memorylimit -= SCIPgetMemUsed(scip)/1048576.0;
1273 memorylimit -= SCIPgetMemExternEstim(scip)/1048576.0;
1277 if( timelimit > 0.0 && memorylimit > 2.0*SCIPgetMemExternEstim(scip)/1048576.0 )
1280 SCIP_CALL( SCIPsetBoolParam(probingscip,
"misc/catchctrlc", FALSE) );
1284 SCIP_CALL( SCIPsetIntParam(probingscip,
"display/verblevel", 0) );
1287 SCIP_CALL( SCIPsetLongintParam(probingscip,
"limits/nodes", 1000LL) );
1288 SCIP_CALL( SCIPsetLongintParam(probingscip,
"limits/stallnodes", 100LL) );
1289 SCIP_CALL( SCIPsetRealParam(probingscip,
"limits/time", timelimit) );
1290 SCIP_CALL( SCIPsetRealParam(probingscip,
"limits/memory", memorylimit) );
1293 SCIP_CALL( SCIPsetSubscipsOff(probingscip, TRUE) );
1294 if( SCIPisParamFixed(probingscip,
"heuristics/"HEUR_NAME"/freq") )
1296 SCIPwarningMessage(scip,
"unfixing parameter heuristics/"HEUR_NAME"/freq in probingscip of "HEUR_NAME" heuristic to avoid recursive calls\n");
1297 SCIP_CALL( SCIPunfixParam(probingscip,
"heuristics/"HEUR_NAME"/freq") );
1299 SCIP_CALL( SCIPsetIntParam(probingscip,
"heuristics/feaspump/freq", -1) );
1302 if( !SCIPisParamFixed(probingscip,
"heuristics/octane/freq") )
1304 SCIP_CALL( SCIPsetIntParam(probingscip,
"heuristics/octane/freq", -1) );
1306 if( !SCIPisParamFixed(probingscip,
"heuristics/objpscostdiving/freq") )
1308 SCIP_CALL( SCIPsetIntParam(probingscip,
"heuristics/objpscostdiving/freq", -1) );
1310 if( !SCIPisParamFixed(probingscip,
"heuristics/rootsoldiving/freq") )
1312 SCIP_CALL( SCIPsetIntParam(probingscip,
"heuristics/rootsoldiving/freq", -1) );
1316 SCIP_CALL( SCIPsetSeparating(probingscip, SCIP_PARAMSETTING_OFF, TRUE) );
1319 SCIP_CALL( SCIPsetPresolving(probingscip, SCIP_PARAMSETTING_FAST, TRUE) );
1322 if( SCIPfindNodesel(probingscip,
"estimate") != NULL && !SCIPisParamFixed(probingscip,
"nodeselection/estimate/stdpriority") )
1324 SCIP_CALL( SCIPsetIntParam(probingscip,
"nodeselection/estimate/stdpriority", INT_MAX/4) );
1328 if( SCIPfindBranchrule(probingscip,
"inference") != NULL && !SCIPisParamFixed(probingscip,
"branching/inference/priority") )
1330 SCIP_CALL( SCIPsetIntParam(probingscip,
"branching/inference/priority", INT_MAX/4) );
1334 if( !SCIPisParamFixed(probingscip,
"conflict/useprop") )
1336 SCIP_CALL( SCIPsetBoolParam(probingscip,
"conflict/useprop", FALSE) );
1338 if( !SCIPisParamFixed(probingscip,
"conflict/useinflp") )
1340 SCIP_CALL( SCIPsetBoolParam(probingscip,
"conflict/useinflp", FALSE) );
1342 if( !SCIPisParamFixed(probingscip,
"conflict/useboundlp") )
1344 SCIP_CALL( SCIPsetBoolParam(probingscip,
"conflict/useboundlp", FALSE) );
1346 if( !SCIPisParamFixed(probingscip,
"conflict/usesb") )
1348 SCIP_CALL( SCIPsetBoolParam(probingscip,
"conflict/usesb", FALSE) );
1350 if( !SCIPisParamFixed(probingscip,
"conflict/usepseudo") )
1352 SCIP_CALL( SCIPsetBoolParam(probingscip,
"conflict/usepseudo", FALSE) );
1356 mindistance = SCIPceil(scip, 2.2*mindistance);
1364 retcode = SCIPsolve(probingscip);
1365 if( retcode != SCIP_OKAY )
1367 SCIPwarningMessage(scip,
"Error while solving sub-SCIP in stage 3 of feasibility pump heuristic; sub-SCIP terminated with code <%d>\n",
1371 SCIP_CALL( SCIPsolve(probingscip) );
1374 if( SCIPgetNSols(probingscip) > 0 )
1377 SCIP_CALL(
createNewSols(scip, probingscip, varmapfw, heur, &success) );
1379 *result = SCIP_FOUNDSOL;
1384 if( *result == SCIP_FOUNDSOL )
1385 heurdata->nsuccess++;
1388 if( varmapfw != NULL )
1389 SCIPhashmapFree(&varmapfw);
1391 if( probingscip != NULL )
1393 SCIP_CALL( SCIPfree(&probingscip) );
1396 if( heurdata->stage3 )
1398 SCIP_CALL( SCIPfreeSol(scip, &closestsol) );
1402 for( j = 0; j < heurdata->cyclelength; j++ )
1404 SCIP_CALL( SCIPfreeSol(scip, &lastroundedsols[j]) );
1407 SCIPfreeBufferArray(scip, &cycles);
1408 SCIPfreeBufferArray(scip, &lastalphas);
1409 SCIPfreeBufferArray(scip, &lastroundedsols);
1410 SCIPfreeBufferArray(scip, &mostfracvals);
1411 SCIPfreeBufferArray(scip, &mostfracvars);
1413 SCIPdebugMessage(
"feasibility pump finished [%d iterations done].\n", nloops);
1428 SCIP_HEURDATA* heurdata;
1432 SCIP_CALL( SCIPallocMemory(scip, &heurdata) );
1435 SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
1439 assert(heur != NULL);
1442 SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeGcgfeaspump) );
1443 SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitGcgfeaspump) );
1444 SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitGcgfeaspump) );
1447 SCIP_CALL( SCIPaddRealParam(scip,
1449 "maximal fraction of diving LP iterations compared to node LP iterations",
1451 SCIP_CALL( SCIPaddRealParam(scip,
1453 "factor by which the regard of the objective is decreased in each round, 1.0 for dynamic",
1455 SCIP_CALL( SCIPaddRealParam(scip,
1457 "threshold difference for the convex parameter to perform perturbation",
1459 SCIP_CALL( SCIPaddIntParam(scip,
1461 "additional number of allowed LP iterations",
1463 SCIP_CALL( SCIPaddIntParam(scip,
1465 "total number of feasible solutions found up to which heuristic is called (-1: no limit)",
1466 &heurdata->maxsols, TRUE,
DEFAULT_MAXSOLS, -1, INT_MAX, NULL, NULL) );
1467 SCIP_CALL( SCIPaddIntParam(scip,
1469 "maximal number of pumping loops (-1: no limit)",
1471 SCIP_CALL( SCIPaddIntParam(scip,
1473 "maximal number of pumping rounds without fractionality improvement (-1: no limit)",
1475 SCIP_CALL( SCIPaddIntParam(scip,
1477 "minimum number of random variables to flip, if a 1-cycle is encountered",
1479 SCIP_CALL( SCIPaddIntParam(scip,
1481 "maximum length of cycles to be checked explicitly in each round",
1483 SCIP_CALL( SCIPaddIntParam(scip,
1485 "number of iterations until a random perturbation is forced",
1487 SCIP_CALL( SCIPaddIntParam(scip,
"heuristics/"HEUR_NAME"/neighborhoodsize",
1488 "radius (using Manhattan metric) of the neighborhood to be searched in stage 3",
1490 SCIP_CALL( SCIPaddBoolParam(scip,
1492 "should an iterative round-and-propagate scheme be used to find the integral points?",
1494 SCIP_CALL( SCIPaddBoolParam(scip,
1496 "should a random perturbation be performed if a feasible solution was found?",
1498 SCIP_CALL( SCIPaddBoolParam(scip,
1500 "should we solve a local branching sub-MIP if no solution could be found?",
1502 SCIP_CALL( SCIPaddBoolParam(scip,
"heuristics/"HEUR_NAME"/copycuts",
1503 "should all active cuts from cutpool be copied to constraints in subproblem?",