diff options
Diffstat (limited to 'test/monniaux/glpk-4.65/src/proxy/proxy.c')
| -rw-r--r-- | test/monniaux/glpk-4.65/src/proxy/proxy.c | 1073 |
1 files changed, 1073 insertions, 0 deletions
diff --git a/test/monniaux/glpk-4.65/src/proxy/proxy.c b/test/monniaux/glpk-4.65/src/proxy/proxy.c new file mode 100644 index 00000000..7d890003 --- /dev/null +++ b/test/monniaux/glpk-4.65/src/proxy/proxy.c @@ -0,0 +1,1073 @@ +/* proxy.c (proximity search heuristic algorithm) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Author: Giorgio Sartor <0gioker0@gmail.com>. +* +* Copyright (C) 2013, 2016 Andrew Makhorin, Department for Applied +* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights +* reserved. E-mail: <mao@gnu.org>. +* +* GLPK is free software: you can redistribute it and/or modify it +* under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* GLPK is distributed in the hope that it will be useful, but WITHOUT +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +* License for more details. +* +* You should have received a copy of the GNU General Public License +* along with GLPK. If not, see <http://www.gnu.org/licenses/>. +* +************************************************************************ +* +* THIS CODE IS AN IMPLEMENTATION OF THE ALGORITHM PROPOSED IN +* +* M. Fischetti, M. Monaci, +* "Proximity Search for 0-1 Mixed-Integer Convex Programming" +* Technical Report DEI, University of Padua, March 2013. +* +* AVAILABLE AT +* http://www.dei.unipd.it/~fisch/papers/proximity_search.pdf +* +* THE CODE HAS BEEN WRITTEN BY GIORGIO SARTOR, " 0gioker0@gmail.com " +* +* BASIC IDEA: +* +* The initial feasible solution x_tilde is defined. This initial +* solution can be found by an ad-hoc heuristic and proxy can be used to +* refine it by exploiting an underlying MIP model whose solution from +* scratch turned out to be problematic. Otherwise, x_tilde can be found +* by running the GLPK mip solver until a first feasible solution is +* found, setting a conservative time limit of 10 minutes (by default). +* Time limit can be modified passing variable tlim [ms]. +* +* Then the cutoff tolerance "delta" is defined. The default tolerance +* is 1% of the last feasible solution obj value--rounded to integer if +* all the variables and obj coefficients are integer. +* +* Next, the objective function c' x is replaced by the Hamming distance +* between x (the actual obj coefficients) and x_tilde (the given +* solution). Distance is only computed wrt the binary variables. +* +* The GLPK solver is then invoked to hopefully find a new incumbent +* x_star with cost c' x_star <= c' x_tilde - delta. A crucial property +* here is that the root-node solution of the LP relaxation is expected +* to be not too different from x_tilde, as this latter solution would +* be optimal without the cutoff constraint, that for a small delta can +* typically be fulfilled with just local adjustments. +* +* If no new solution x_star is found within the time limit the +* algorithm stops. Of course, if the MIP solver proved infeasibility +* for the given delta, we have that c' x_tilde - delta is a valid lower +* bound (in case of minimazation) on the optimal value of the original +* MIP. +* +* The new solution x_star, if any, is possibly improved by solving a +* simplified problem (refinement) where all binary variables have been +* fixed to their value in x_star so as to find the best solution within +* the neighborhood. +* +* Finally, the approach is reapplied on x_star (that replaces x_tilde) +* so as to recenter the distance Hamming function and by modifying the +* cutoff tolerance delta. +* +* In this way, there will be a series of hopefully not-too-difficult +* sub-MIPs to solve, each leading to an improvement of the incumbent. +* More aggressive policies on the definition of tolerance delta can +* lead to a better performance, but would require an ad-hoc tuning. +* +************************************************************************ +* +* int proxy(glp_prob *lp, double *zstar, double *xstar, +* const double[] initsol, double rel_impr, int tlim, +* int verbose) +* +* lp : GLPK problem pointer to a MIP with binary variables +* +* zstar : the value of objective function of the best solution found +* +* xstar : best solution with components xstar[1],...,xstar[ncols] +* +* initsol : pointer to a initial feasible solution, see +* glp_ios_heur_sol +* If initsol = NULL, the procedure finds the first solution +* by itself. +* +* rel_impr : minimum relative obj improvement to be achieved at each +* internal step; if <= 0.0 a default value of 0.01 (1%) is +* used; for some problems (e.g., set covering with small +* integer costs) a more-conservative choice of 0.001 (0.1%) +* can lead to a better final solution; values larger than +* 0.05 (5%) are typically too aggressive and do not work +* well. +* +* tlim : time limit to find a new solution, in ms. +* If tlim = 0, it is set to its default value, 600000 ms +* +* verbose : if 1 the output is activated. If 0 only errors are +* displayed +* +* The procedure returns -1 if an error occurred, 0 otherwise (possibly, +* time limit) +* +***********************************************************************/ + +/**********************************************************************/ +/* 1. INCLUDE */ +/**********************************************************************/ + +#include "glpk.h" +#include "env.h" +#include "proxy.h" + +/**********************************************************************/ +/* 2. PARAMETERS AND CONSTANTS */ +/**********************************************************************/ + +#define TDAY 86400.0 +#define TRUE 1 +#define FALSE 0 +#define EPS 1e-6 +#define RINF 1e38 +#define MAXVAL 1e20 +#define MINVAL -1e20 +#if 0 /* by gioker */ + #define PROXY_DEBUG +#endif + +/**********************************************************************/ +/* 3. GLOBAL VARIABLES */ +/**********************************************************************/ + +struct csa { + +int integer_obj; /* TRUE if each feasible solution has an + integral cost */ +int b_vars_exist; /* TRUE if there is at least one binary + variable in the problem */ +int i_vars_exist; /* TRUE if there is at least one general + integer variable in the problem */ +const double *startsol; /* Pointer to the initial solution */ + +int *ckind; /* Store the kind of the structural variables + of the problem */ +double *clb; /* Store the lower bound on the structural + variables of the problem */ +double *cub; /* Store the upper bound on the structural + variables of the problem */ +double *true_obj; /* Store the obj coefficients of the problem */ + +int dir; /* Minimization or maximization problem */ +int ncols; /* Number of structural variables of the + problem */ + +time_t GLOtstart; /* starting time of the algorithm */ + +glp_prob *lp_ref; /* glp problem for refining only*/ + +}; + +/**********************************************************************/ +/* 4. FUNCTIONS PROTOTYPES */ +/**********************************************************************/ + +static void callback(glp_tree *tree, void *info); +static void get_info(struct csa *csa, glp_prob *lp); +static int is_integer(struct csa *csa); +static void check_integrality(struct csa *csa); +static int check_ref(struct csa *csa, glp_prob *lp, double *xref); +static double second(void); +static int add_cutoff(struct csa *csa, glp_prob *lp); +static void get_sol(struct csa *csa, glp_prob *lp, double *xstar); +static double elapsed_time(struct csa *csa); +static void redefine_obj(glp_prob *lp, double *xtilde, int ncols, + int *ckind, double *clb, double *cub); +static double update_cutoff(struct csa *csa, glp_prob *lp, + double zstar, int index, double rel_impr); +static double compute_delta(struct csa *csa, double z, + double rel_impr); +static double objval(int ncols, double *x, double *true_obj); +static void array_copy(int begin, int end, double *source, + double *destination); +static int do_refine(struct csa *csa, glp_prob *lp_ref, int ncols, + int *ckind, double *xref, int *tlim, int tref_lim, + int verbose); +static void deallocate(struct csa *csa, int refine); + +/**********************************************************************/ +/* 5. FUNCTIONS */ +/**********************************************************************/ + +int proxy(glp_prob *lp, double *zfinal, double *xfinal, + const double initsol[], double rel_impr, int tlim, + int verbose) + +{ struct csa csa_, *csa = &csa_; + glp_iocp parm; + glp_smcp parm_lp; + size_t tpeak; + int refine, tref_lim, err, cutoff_row, niter, status, i, tout; + double *xref, *xstar, zstar, tela, cutoff, zz; + + memset(csa, 0, sizeof(struct csa)); + + + /********** **********/ + /********** RETRIEVING PROBLEM INFO **********/ + /********** **********/ + + /* getting problem direction (min or max) */ + csa->dir = glp_get_obj_dir(lp); + + /* getting number of variables */ + csa->ncols = glp_get_num_cols(lp); + + /* getting kind, bounds and obj coefficient of each variable + information is stored in ckind, cub, clb, true_obj */ + get_info(csa, lp); + + /* checking if the objective function is always integral */ + check_integrality(csa); + + /* Proximity search cannot be used if there are no binary + variables */ + if (csa->b_vars_exist == FALSE) { + if (verbose) { + xprintf("The problem has not binary variables. Proximity se" + "arch cannot be used.\n"); + } + tfree(csa->ckind); + tfree(csa->clb); + tfree(csa->cub); + tfree(csa->true_obj); + return -1; + } + + /* checking if the problem needs refinement, i.e., not all + variables are binary. If so, the routine creates a copy of the + lp problem named lp_ref and initializes the solution xref to + zero. */ + xref = talloc(csa->ncols+1, double); +#if 0 /* by mao */ + memset(xref, 0, sizeof(double)*(csa->ncols+1)); +#endif + refine = check_ref(csa, lp, xref); +#ifdef PROXY_DEBUG + xprintf("REFINE = %d\n",refine); +#endif + + /* Initializing the solution */ + xstar = talloc(csa->ncols+1, double); +#if 0 /* by mao */ + memset(xstar, 0, sizeof(double)*(csa->ncols+1)); +#endif + + /********** **********/ + /********** FINDING FIRST SOLUTION **********/ + /********** **********/ + + if (verbose) { + xprintf("Applying PROXY heuristic...\n"); + } + + /* get the initial time */ + csa->GLOtstart = second(); + + /* setting the optimization parameters */ + glp_init_iocp(&parm); + glp_init_smcp(&parm_lp); +#if 0 /* by gioker */ + /* Preprocessing should be disabled because the mip passed + to proxy is already preprocessed */ + parm.presolve = GLP_ON; +#endif +#if 1 /* by mao */ + /* best projection backtracking seems to be more efficient to find + any integer feasible solution */ + parm.bt_tech = GLP_BT_BPH; +#endif + + /* Setting the default value of the minimum relative improvement + to 1% */ + if ( rel_impr <= 0.0 ) { + rel_impr = 0.01; + } + + /* Setting the default value of time limit to 10 minutes */ + if (tlim <= 0) { + tlim = INT_MAX; + } + if (verbose) { + xprintf("Proxy's time limit set to %d seconds.\n",tlim/1000); + xprintf("Proxy's relative improvement " + "set to %2.2lf %c.\n",rel_impr*100,37); + } + + parm_lp.tm_lim = tlim; + + parm.mip_gap = 9999999.9; /* to stop the optimization at the first + feasible solution found */ + + /* finding the first solution */ + if (verbose) { + xprintf("Searching for a feasible solution...\n"); + } + + /* verifying the existence of an input starting solution */ + if (initsol != NULL) { + csa->startsol = initsol; + parm.cb_func = callback; + parm.cb_info = csa; + if (verbose) { + xprintf("Input solution found.\n"); + } + } + + tout = glp_term_out(GLP_OFF); + err = glp_simplex(lp,&parm_lp); + glp_term_out(tout); + + status = glp_get_status(lp); + + if (status != GLP_OPT) { + if (verbose) { + xprintf("Proxy heuristic terminated.\n"); + } +#ifdef PROXY_DEBUG + /* For debug only */ + xprintf("GLP_SIMPLEX status = %d\n",status); + xprintf("GLP_SIMPLEX error code = %d\n",err); +#endif + tfree(xref); + tfree(xstar); + deallocate(csa, refine); + return -1; + } + + tela = elapsed_time(csa); + if (tlim-tela*1000 <= 0) { + if (verbose) { + xprintf("Time limit exceeded. Proxy could not " + "find optimal solution to LP relaxation.\n"); + xprintf("Proxy heuristic aborted.\n"); + } + tfree(xref); + tfree(xstar); + deallocate(csa, refine); + return -1; + } + + parm.tm_lim = tlim - tela*1000; + tref_lim = (tlim - tela *1000) / 20; + + tout = glp_term_out(GLP_OFF); + err = glp_intopt(lp, &parm); + glp_term_out(tout); + + status = glp_mip_status(lp); + + /***** If no solution was found *****/ + + if (status == GLP_NOFEAS || status == GLP_UNDEF) { + if (err == GLP_ETMLIM) { + if (verbose) { + xprintf("Time limit exceeded. Proxy could not " + "find an initial integer feasible solution.\n"); + xprintf("Proxy heuristic aborted.\n"); + } + } + else { + if (verbose) { + xprintf("Proxy could not " + "find an initial integer feasible solution.\n"); + xprintf("Proxy heuristic aborted.\n"); + } + } + tfree(xref); + tfree(xstar); + deallocate(csa, refine); + return -1; + } + + /* getting the first solution and its value */ + get_sol(csa, lp,xstar); + zstar = glp_mip_obj_val(lp); + + if (verbose) { + xprintf(">>>>> first solution = %e;\n", zstar); + } + + /* If a feasible solution was found but the time limit is + exceeded */ + if (err == GLP_ETMLIM) { + if (verbose) { + xprintf("Time limit exceeded. Proxy heuristic terminated.\n"); + } + goto done; + } + + tela = elapsed_time(csa); + tpeak = 0; + glp_mem_usage(NULL, NULL, NULL, &tpeak); + if (verbose) { + xprintf("Time used: %3.1lf secs. Memory used: %2.1lf Mb\n", + tela,(double)tpeak/1048576); + xprintf("Starting proximity search...\n"); + } + + /********** **********/ + /********** PREPARING THE PROBLEM FOR PROXY **********/ + /********** **********/ + + /* adding a dummy cutoff constraint */ + cutoff_row = add_cutoff(csa, lp); + + /* proximity search needs minimization direction + even if the problem is a maximization one */ + if (csa->dir == GLP_MAX) { + glp_set_obj_dir(lp, GLP_MIN); + } + + /********** **********/ + /********** STARTING PROXIMITY SEARCH **********/ + /********** **********/ + + + niter = 0; + + while (TRUE) { + niter++; + + /********** CHANGING THE OBJ FUNCTION **********/ + + redefine_obj(lp,xstar, csa->ncols, csa->ckind, csa->clb, + csa->cub); + + /********** UPDATING THE CUTOFF CONSTRAINT **********/ + + cutoff = update_cutoff(csa, lp,zstar, cutoff_row, rel_impr); + +#ifdef PROXY_DEBUG + xprintf("TRUE_OBJ[0] = %f\n",csa->true_obj[0]); + xprintf("ZSTAR = %f\n",zstar); + xprintf("CUTOFF = %f\n",cutoff); +#endif + + /********** SEARCHING FOR A BETTER SOLUTION **********/ + + tela = elapsed_time(csa); + if (tlim-tela*1000 <= 0) { + if (verbose) { + xprintf("Time limit exceeded. Proxy heuristic " + "terminated.\n"); + } + goto done; + } +#ifdef PROXY_DEBUG + xprintf("TELA = %3.1lf\n",tela*1000); + xprintf("TLIM = %3.1lf\n",tlim - tela*1000); +#endif + parm_lp.tm_lim = tlim -tela*1000; + + tout = glp_term_out(GLP_OFF); + err = glp_simplex(lp,&parm_lp); + glp_term_out(tout); + + status = glp_get_status(lp); + + if (status != GLP_OPT) { + if (status == GLP_NOFEAS) { + if (verbose) { + xprintf("Bound exceeded = %f. ",cutoff); + } + } + if (verbose) { + xprintf("Proxy heuristic terminated.\n"); + } +#ifdef PROXY_DEBUG + xprintf("GLP_SIMPLEX status = %d\n",status); + xprintf("GLP_SIMPLEX error code = %d\n",err); +#endif + goto done; + } + + tela = elapsed_time(csa); + if (tlim-tela*1000 <= 0) { + if (verbose) { + xprintf("Time limit exceeded. Proxy heuristic " + "terminated.\n"); + } + goto done; + } + parm.tm_lim = tlim - tela*1000; + parm.cb_func = NULL; +#if 0 /* by gioker */ + /* Preprocessing should be disabled because the mip passed + to proxy is already preprocessed */ + parm.presolve = GLP_ON; +#endif + tout = glp_term_out(GLP_OFF); + err = glp_intopt(lp, &parm); + glp_term_out(tout); + + /********** MANAGEMENT OF THE SOLUTION **********/ + + status = glp_mip_status(lp); + + /***** No feasible solutions *****/ + + if (status == GLP_NOFEAS) { + if (verbose) { + xprintf("Bound exceeded = %f. Proxy heuristic " + "terminated.\n",cutoff); + } + goto done; + } + + /***** Undefined solution *****/ + + if (status == GLP_UNDEF) { + if (err == GLP_ETMLIM) { + if (verbose) { + xprintf("Time limit exceeded. Proxy heuristic " + "terminated.\n"); + } + } + else { + if (verbose) { + xprintf("Proxy terminated unexpectedly.\n"); +#ifdef PROXY_DEBUG + xprintf("GLP_INTOPT error code = %d\n",err); +#endif + } + } + goto done; + } + + /***** Feasible solution *****/ + + if ((status == GLP_FEAS) || (status == GLP_OPT)) { + + /* getting the solution and computing its value */ + get_sol(csa, lp,xstar); + zz = objval(csa->ncols, xstar, csa->true_obj); + + /* Comparing the incumbent solution with the current best + one */ +#ifdef PROXY_DEBUG + xprintf("ZZ = %f\n",zz); + xprintf("ZSTAR = %f\n",zstar); + xprintf("REFINE = %d\n",refine); +#endif + if (((zz<zstar) && (csa->dir == GLP_MIN)) || + ((zz>zstar) && (csa->dir == GLP_MAX))) { + + /* refining (possibly) the solution */ + if (refine) { + + /* copying the incumbent solution in the refinement + one */ + array_copy(1, csa->ncols +1, xstar, xref); + err = do_refine(csa, csa->lp_ref, csa->ncols, + csa->ckind, xref, &tlim, tref_lim, verbose); + if (!err) { + double zref = objval(csa->ncols, xref, + csa->true_obj); + if (((zref<zz) && (csa->dir == GLP_MIN)) || + ((zref>zz) && (csa->dir == GLP_MAX))) { + zz = zref; + /* copying the refinement solution in the + incumbent one */ + array_copy(1, csa->ncols +1, xref, xstar); + } + } + } + zstar = zz; + tela = elapsed_time(csa); + if (verbose) { + xprintf(">>>>> it: %3d: mip = %e; elapsed time " + "%3.1lf sec.s\n", niter,zstar,tela); + } + } + } + } + +done: + tela = elapsed_time(csa); + glp_mem_usage(NULL, NULL, NULL, &tpeak); + if (verbose) { + xprintf("Time used: %3.1lf. Memory used: %2.1lf Mb\n", + tela,(double)tpeak/1048576); + } + + + /* Exporting solution and obj val */ + *zfinal = zstar; + + for (i=1; i < (csa->ncols + 1); i++) { + xfinal[i]=xstar[i]; + } + + /* Freeing allocated memory */ + tfree(xref); + tfree(xstar); + deallocate(csa, refine); + + return 0; +} + +/**********************************************************************/ +static void callback(glp_tree *tree, void *info){ +/**********************************************************************/ + struct csa *csa = info; + switch(glp_ios_reason(tree)) { + case GLP_IHEUR: + glp_ios_heur_sol(tree, csa->startsol); + break; + default: break; + } +} + +/**********************************************************************/ +static void get_info(struct csa *csa, glp_prob *lp) +/**********************************************************************/ +{ + int i; + + /* Storing helpful info of the problem */ + + csa->ckind = talloc(csa->ncols+1, int); +#if 0 /* by mao */ + memset(csa->ckind, 0, sizeof(int)*(csa->ncols+1)); +#endif + csa->clb = talloc(csa->ncols+1, double); +#if 0 /* by mao */ + memset(csa->clb, 0, sizeof(double)*(csa->ncols+1)); +#endif + csa->cub = talloc(csa->ncols+1, double); +#if 0 /* by mao */ + memset(csa->cub, 0, sizeof(double)*(csa->ncols+1)); +#endif + csa->true_obj = talloc(csa->ncols+1, double); +#if 0 /* by mao */ + memset(csa->true_obj, 0, sizeof(double)*(csa->ncols+1)); +#endif + for( i = 1 ; i < (csa->ncols + 1); i++ ) { + csa->ckind[i] = glp_get_col_kind(lp, i); + csa->clb[i] = glp_get_col_lb(lp, i); + csa->cub[i] = glp_get_col_ub(lp, i); + csa->true_obj[i] = glp_get_obj_coef(lp, i); + } + csa->true_obj[0] = glp_get_obj_coef(lp, 0); +} + +/**********************************************************************/ +static int is_integer(struct csa *csa) +/**********************************************************************/ +{ + int i; + csa->integer_obj = TRUE; + for ( i = 1; i < (csa->ncols + 1); i++ ) { + if (fabs(csa->true_obj[i]) > INT_MAX ) { + csa->integer_obj = FALSE; + } + if (fabs(csa->true_obj[i]) <= INT_MAX) { + double tmp, rem; + if (fabs(csa->true_obj[i]) - floor(fabs(csa->true_obj[i])) + < 0.5) { + tmp = floor(fabs(csa->true_obj[i])); + } + else { + tmp = ceil(fabs(csa->true_obj[i])); + } + rem = fabs(csa->true_obj[i]) - tmp; + rem = fabs(rem); + if (rem > EPS) { + csa->integer_obj = FALSE; + } + + } + } + return csa->integer_obj; +} + +/**********************************************************************/ +static void check_integrality(struct csa *csa) +/**********************************************************************/ +{ + /* + Checking if the problem has binary, integer or continuos variables. + integer_obj is TRUE if the problem has no continuous variables + and all the obj coefficients are integer (and < INT_MAX). + */ + + int i; + csa->integer_obj = is_integer(csa); + csa->b_vars_exist = FALSE; + csa->i_vars_exist = FALSE; + for ( i = 1; i < (csa->ncols + 1); i++ ) { + if ( csa->ckind[i] == GLP_IV ){ + csa->i_vars_exist = TRUE; + continue; + } + if ( csa->ckind[i] == GLP_BV ){ + csa->b_vars_exist =TRUE; + continue; + } + csa->integer_obj = FALSE; + } +} + +/**********************************************************************/ +static int check_ref(struct csa *csa, glp_prob *lp, double *xref) +/**********************************************************************/ +{ + /* + checking if the problem has continuos or integer variables. If so, + refinement is prepared. + */ + int refine = FALSE; + int i; + for ( i = 1; i < (csa->ncols + 1); i++ ) { + if ( csa->ckind[i] != GLP_BV) { + refine = TRUE; + break; + } + } + + /* possibly creating a mip clone for refinement only */ + if ( refine ) { + csa->lp_ref = glp_create_prob(); + glp_copy_prob(csa->lp_ref, lp, GLP_ON); + } + + return refine; +} + +/**********************************************************************/ +static double second(void) +/**********************************************************************/ +{ +#if 0 /* by mao */ + return ((double)clock()/(double)CLOCKS_PER_SEC); +#else + return xtime() / 1000.0; +#endif +} + +/**********************************************************************/ +static int add_cutoff(struct csa *csa, glp_prob *lp) +/**********************************************************************/ +{ + /* + Adding a cutoff constraint to set an upper bound (in case of + minimaztion) on the obj value of the next solution, i.e., the next + value of the true obj function that we would like to find + */ + + /* store non-zero coefficients in the objective function */ + int *obj_index = talloc(csa->ncols+1, int); +#if 0 /* by mao */ + memset(obj_index, 0, sizeof(int)*(csa->ncols+1)); +#endif + double *obj_value = talloc(csa->ncols+1, double); +#if 0 /* by mao */ + memset(obj_value, 0, sizeof(double)*(csa->ncols+1)); +#endif + int obj_nzcnt = 0; + int i, irow; + const char *rowname; + for ( i = 1; i < (csa->ncols + 1); i++ ) { + if ( fabs(csa->true_obj[i]) > EPS ) { + obj_nzcnt++; + obj_index[obj_nzcnt] = i; + obj_value[obj_nzcnt] = csa->true_obj[i]; + } + } + + irow = glp_add_rows(lp, 1); + rowname = "Cutoff"; + glp_set_row_name(lp, irow, rowname); + if (csa->dir == GLP_MIN) { + /* minimization problem */ + glp_set_row_bnds(lp, irow, GLP_UP, MAXVAL, MAXVAL); + } + else { + /* maximization problem */ + glp_set_row_bnds(lp, irow, GLP_LO, MINVAL, MINVAL); + } + + glp_set_mat_row(lp, irow, obj_nzcnt, obj_index, obj_value); + + tfree(obj_index); + tfree(obj_value); + + return irow; +} + +/**********************************************************************/ +static void get_sol(struct csa *csa, glp_prob *lp, double *xstar) +/**********************************************************************/ +{ + /* Retrieving and storing the coefficients of the solution */ + + int i; + for (i = 1; i < (csa->ncols +1); i++) { + xstar[i] = glp_mip_col_val(lp, i); + } +} + +/**********************************************************************/ +static double elapsed_time(struct csa *csa) +/**********************************************************************/ +{ + double tela = second() - csa->GLOtstart; + if ( tela < 0 ) tela += TDAY; + return(tela); +} + +/**********************************************************************/ +static void redefine_obj(glp_prob *lp, double *xtilde, int ncols, + int *ckind, double *clb, double *cub) +/**********************************************************************/ + +/* + Redefine the lp objective function obj as the distance-to-integrality + (Hamming distance) from xtilde (the incumbent feasible solution), wrt + to binary vars only + */ + +{ + int j; + double *delta = talloc(ncols+1, double); +#if 0 /* by mao */ + memset(delta, 0, sizeof(double)*(ncols+1)); +#endif + + for ( j = 1; j < (ncols +1); j++ ) { + delta[j] = 0.0; + /* skip continuous variables */ + if ( ckind[j] == GLP_CV ) continue; + + /* skip integer variables that have been fixed */ + if ( cub[j]-clb[j] < 0.5 ) continue; + + /* binary variable */ + if ( ckind[j] == GLP_BV ) { + if ( xtilde[j] > 0.5 ) { + delta[j] = -1.0; + } + else { + delta[j] = 1.0; + } + } + } + + /* changing the obj coeff. for all variables, including continuous + ones */ + for ( j = 1; j < (ncols +1); j++ ) { + glp_set_obj_coef(lp, j, delta[j]); + } + glp_set_obj_coef(lp, 0, 0.0); + + tfree(delta); +} + +/**********************************************************************/ +static double update_cutoff(struct csa *csa, glp_prob *lp, + double zstar, int cutoff_row, + double rel_impr) +/**********************************************************************/ +{ + /* + Updating the cutoff constraint with the value we would like to + find during the next optimization + */ + double cutoff; + zstar -= csa->true_obj[0]; + if (csa->dir == GLP_MIN) { + cutoff = zstar - compute_delta(csa, zstar, rel_impr); + glp_set_row_bnds(lp, cutoff_row, GLP_UP, cutoff, cutoff); + } + else { + cutoff = zstar + compute_delta(csa, zstar, rel_impr); + glp_set_row_bnds(lp, cutoff_row, GLP_LO, cutoff, cutoff); + } + + return cutoff; +} + +/**********************************************************************/ +static double compute_delta(struct csa *csa, double z, double rel_impr) +/**********************************************************************/ +{ + /* Computing the offset for the next best solution */ + + double delta = rel_impr * fabs(z); + if ( csa->integer_obj ) delta = ceil(delta); + + return(delta); +} + +/**********************************************************************/ +static double objval(int ncols, double *x, double *true_obj) +/**********************************************************************/ +{ + /* Computing the true cost of x (using the original obj coeff.s) */ + + int j; + double z = 0.0; + for ( j = 1; j < (ncols +1); j++ ) { + z += x[j] * true_obj[j]; + } + return z + true_obj[0]; +} + +/**********************************************************************/ +static void array_copy(int begin, int end, double *source, + double *destination) +/**********************************************************************/ +{ + int i; + for (i = begin; i < end; i++) { + destination[i] = source[i]; + } +} +/**********************************************************************/ +static int do_refine(struct csa *csa, glp_prob *lp_ref, int ncols, + int *ckind, double *xref, int *tlim, int tref_lim, + int verbose) +/**********************************************************************/ +{ + /* + Refinement is applied when the variables of the problem are not + all binary. Binary variables are fixed to their value and + remaining ones are optimized. If there are only continuos + variables (in addition to those binary) the problem becomes just + an LP. Otherwise, it remains a MIP but of smaller size. + */ + + int j, tout; + double refineStart = second(); + double val, tela, tlimit; + + if ( glp_get_num_cols(lp_ref) != ncols ) { + if (verbose) { + xprintf("Error in Proxy refinement: "); + xprintf("wrong number of columns (%d vs %d).\n", + ncols, glp_get_num_cols(lp_ref)); + } + return 1; + } + + val = -1.0; + + /* fixing all binary variables to their current value in xref */ + for ( j = 1; j < (ncols + 1); j++ ) { + if ( ckind[j] == GLP_BV ) { + val = 0.0; + if ( xref[j] > 0.5 ) val = 1.0; + glp_set_col_bnds(lp_ref, j, GLP_FX, val, val); + } + } + + /* re-optimizing (refining) if some bound has been changed */ + if ( val > -1.0 ) { + glp_iocp parm_ref; + glp_smcp parm_ref_lp; + int err, status; + + glp_init_iocp(&parm_ref); + parm_ref.presolve = GLP_ON; + glp_init_smcp(&parm_ref_lp); + /* + If there are no general integer variable the problem becomes + an LP (after fixing the binary variables) and can be solved + quickly. Otherwise the problem is still a MIP problem and a + timelimit has to be set. + */ + parm_ref.tm_lim = tref_lim; + if (parm_ref.tm_lim > *tlim) { + parm_ref.tm_lim = *tlim; + } + parm_ref_lp.tm_lim = parm_ref.tm_lim; +#ifdef PROXY_DEBUG + xprintf("***** REFINING *****\n"); +#endif + tout = glp_term_out(GLP_OFF); + if (csa->i_vars_exist == TRUE) { + err = glp_intopt(lp_ref, &parm_ref); + } + else { + err = glp_simplex(lp_ref, &parm_ref_lp); + } + glp_term_out(tout); + + if (csa->i_vars_exist == TRUE) { + status = glp_mip_status(lp_ref); + } + else { + status = glp_get_status(lp_ref); + } + +#if 1 /* 29/II-2016 by mao as reported by Chris */ + switch (status) + { case GLP_OPT: + case GLP_FEAS: + break; + default: + status = GLP_UNDEF; + break; + } +#endif + +#ifdef PROXY_DEBUG + xprintf("STATUS REFINING = %d\n",status); +#endif + if (status == GLP_UNDEF) { + if (err == GLP_ETMLIM) { +#ifdef PROXY_DEBUG + xprintf("Time limit exceeded on Proxy refining.\n"); +#endif + return 1; + } + } + for( j = 1 ; j < (ncols + 1); j++ ){ + if (ckind[j] != GLP_BV) { + if (csa->i_vars_exist == TRUE) { + xref[j] = glp_mip_col_val(lp_ref, j); + } + else{ + xref[j] = glp_get_col_prim(lp_ref, j); + } + } + } + } + tela = second() - refineStart; +#ifdef PROXY_DEBUG + xprintf("REFINE TELA = %3.1lf\n",tela*1000); +#endif + return 0; +} +/**********************************************************************/ +static void deallocate(struct csa *csa, int refine) +/**********************************************************************/ +{ + /* Deallocating routine */ + + if (refine) { + glp_delete_prob(csa->lp_ref); + } + + tfree(csa->ckind); + tfree(csa->clb); + tfree(csa->cub); + tfree(csa->true_obj); + +} + +/* eof */ |