|
1 /* glpios09.c (branching heuristics) */ |
|
2 |
|
3 /*********************************************************************** |
|
4 * This code is part of GLPK (GNU Linear Programming Kit). |
|
5 * |
|
6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
|
7 * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, |
|
8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
|
9 * E-mail: <mao@gnu.org>. |
|
10 * |
|
11 * GLPK is free software: you can redistribute it and/or modify it |
|
12 * under the terms of the GNU General Public License as published by |
|
13 * the Free Software Foundation, either version 3 of the License, or |
|
14 * (at your option) any later version. |
|
15 * |
|
16 * GLPK is distributed in the hope that it will be useful, but WITHOUT |
|
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
|
18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
|
19 * License for more details. |
|
20 * |
|
21 * You should have received a copy of the GNU General Public License |
|
22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
|
23 ***********************************************************************/ |
|
24 |
|
25 #include "glpios.h" |
|
26 |
|
27 /*********************************************************************** |
|
28 * NAME |
|
29 * |
|
30 * ios_choose_var - select variable to branch on |
|
31 * |
|
32 * SYNOPSIS |
|
33 * |
|
34 * #include "glpios.h" |
|
35 * int ios_choose_var(glp_tree *T, int *next); |
|
36 * |
|
37 * The routine ios_choose_var chooses a variable from the candidate |
|
38 * list to branch on. Additionally the routine provides a flag stored |
|
39 * in the location next to suggests which of the child subproblems |
|
40 * should be solved next. |
|
41 * |
|
42 * RETURNS |
|
43 * |
|
44 * The routine ios_choose_var returns the ordinal number of the column |
|
45 * choosen. */ |
|
46 |
|
47 static int branch_first(glp_tree *T, int *next); |
|
48 static int branch_last(glp_tree *T, int *next); |
|
49 static int branch_mostf(glp_tree *T, int *next); |
|
50 static int branch_drtom(glp_tree *T, int *next); |
|
51 |
|
52 int ios_choose_var(glp_tree *T, int *next) |
|
53 { int j; |
|
54 if (T->parm->br_tech == GLP_BR_FFV) |
|
55 { /* branch on first fractional variable */ |
|
56 j = branch_first(T, next); |
|
57 } |
|
58 else if (T->parm->br_tech == GLP_BR_LFV) |
|
59 { /* branch on last fractional variable */ |
|
60 j = branch_last(T, next); |
|
61 } |
|
62 else if (T->parm->br_tech == GLP_BR_MFV) |
|
63 { /* branch on most fractional variable */ |
|
64 j = branch_mostf(T, next); |
|
65 } |
|
66 else if (T->parm->br_tech == GLP_BR_DTH) |
|
67 { /* branch using the heuristic by Dreebeck and Tomlin */ |
|
68 j = branch_drtom(T, next); |
|
69 } |
|
70 else if (T->parm->br_tech == GLP_BR_PCH) |
|
71 { /* hybrid pseudocost heuristic */ |
|
72 j = ios_pcost_branch(T, next); |
|
73 } |
|
74 else |
|
75 xassert(T != T); |
|
76 return j; |
|
77 } |
|
78 |
|
79 /*********************************************************************** |
|
80 * branch_first - choose first branching variable |
|
81 * |
|
82 * This routine looks up the list of structural variables and chooses |
|
83 * the first one, which is of integer kind and has fractional value in |
|
84 * optimal solution to the current LP relaxation. |
|
85 * |
|
86 * This routine also selects the branch to be solved next where integer |
|
87 * infeasibility of the chosen variable is less than in other one. */ |
|
88 |
|
89 static int branch_first(glp_tree *T, int *_next) |
|
90 { int j, next; |
|
91 double beta; |
|
92 /* choose the column to branch on */ |
|
93 for (j = 1; j <= T->n; j++) |
|
94 if (T->non_int[j]) break; |
|
95 xassert(1 <= j && j <= T->n); |
|
96 /* select the branch to be solved next */ |
|
97 beta = glp_get_col_prim(T->mip, j); |
|
98 if (beta - floor(beta) < ceil(beta) - beta) |
|
99 next = GLP_DN_BRNCH; |
|
100 else |
|
101 next = GLP_UP_BRNCH; |
|
102 *_next = next; |
|
103 return j; |
|
104 } |
|
105 |
|
106 /*********************************************************************** |
|
107 * branch_last - choose last branching variable |
|
108 * |
|
109 * This routine looks up the list of structural variables and chooses |
|
110 * the last one, which is of integer kind and has fractional value in |
|
111 * optimal solution to the current LP relaxation. |
|
112 * |
|
113 * This routine also selects the branch to be solved next where integer |
|
114 * infeasibility of the chosen variable is less than in other one. */ |
|
115 |
|
116 static int branch_last(glp_tree *T, int *_next) |
|
117 { int j, next; |
|
118 double beta; |
|
119 /* choose the column to branch on */ |
|
120 for (j = T->n; j >= 1; j--) |
|
121 if (T->non_int[j]) break; |
|
122 xassert(1 <= j && j <= T->n); |
|
123 /* select the branch to be solved next */ |
|
124 beta = glp_get_col_prim(T->mip, j); |
|
125 if (beta - floor(beta) < ceil(beta) - beta) |
|
126 next = GLP_DN_BRNCH; |
|
127 else |
|
128 next = GLP_UP_BRNCH; |
|
129 *_next = next; |
|
130 return j; |
|
131 } |
|
132 |
|
133 /*********************************************************************** |
|
134 * branch_mostf - choose most fractional branching variable |
|
135 * |
|
136 * This routine looks up the list of structural variables and chooses |
|
137 * that one, which is of integer kind and has most fractional value in |
|
138 * optimal solution to the current LP relaxation. |
|
139 * |
|
140 * This routine also selects the branch to be solved next where integer |
|
141 * infeasibility of the chosen variable is less than in other one. |
|
142 * |
|
143 * (Alexander Martin notices that "...most infeasible is as good as |
|
144 * random...".) */ |
|
145 |
|
146 static int branch_mostf(glp_tree *T, int *_next) |
|
147 { int j, jj, next; |
|
148 double beta, most, temp; |
|
149 /* choose the column to branch on */ |
|
150 jj = 0, most = DBL_MAX; |
|
151 for (j = 1; j <= T->n; j++) |
|
152 { if (T->non_int[j]) |
|
153 { beta = glp_get_col_prim(T->mip, j); |
|
154 temp = floor(beta) + 0.5; |
|
155 if (most > fabs(beta - temp)) |
|
156 { jj = j, most = fabs(beta - temp); |
|
157 if (beta < temp) |
|
158 next = GLP_DN_BRNCH; |
|
159 else |
|
160 next = GLP_UP_BRNCH; |
|
161 } |
|
162 } |
|
163 } |
|
164 *_next = next; |
|
165 return jj; |
|
166 } |
|
167 |
|
168 /*********************************************************************** |
|
169 * branch_drtom - choose branching var using Driebeck-Tomlin heuristic |
|
170 * |
|
171 * This routine chooses a structural variable, which is required to be |
|
172 * integral and has fractional value in optimal solution of the current |
|
173 * LP relaxation, using a heuristic proposed by Driebeck and Tomlin. |
|
174 * |
|
175 * The routine also selects the branch to be solved next, again due to |
|
176 * Driebeck and Tomlin. |
|
177 * |
|
178 * This routine is based on the heuristic proposed in: |
|
179 * |
|
180 * Driebeck N.J. An algorithm for the solution of mixed-integer |
|
181 * programming problems, Management Science, 12: 576-87 (1966); |
|
182 * |
|
183 * and improved in: |
|
184 * |
|
185 * Tomlin J.A. Branch and bound methods for integer and non-convex |
|
186 * programming, in J.Abadie (ed.), Integer and Nonlinear Programming, |
|
187 * North-Holland, Amsterdam, pp. 437-50 (1970). |
|
188 * |
|
189 * Must note that this heuristic is time-expensive, because computing |
|
190 * one-step degradation (see the routine below) requires one BTRAN for |
|
191 * each fractional-valued structural variable. */ |
|
192 |
|
193 static int branch_drtom(glp_tree *T, int *_next) |
|
194 { glp_prob *mip = T->mip; |
|
195 int m = mip->m; |
|
196 int n = mip->n; |
|
197 char *non_int = T->non_int; |
|
198 int j, jj, k, t, next, kase, len, stat, *ind; |
|
199 double x, dk, alfa, delta_j, delta_k, delta_z, dz_dn, dz_up, |
|
200 dd_dn, dd_up, degrad, *val; |
|
201 /* basic solution of LP relaxation must be optimal */ |
|
202 xassert(glp_get_status(mip) == GLP_OPT); |
|
203 /* allocate working arrays */ |
|
204 ind = xcalloc(1+n, sizeof(int)); |
|
205 val = xcalloc(1+n, sizeof(double)); |
|
206 /* nothing has been chosen so far */ |
|
207 jj = 0, degrad = -1.0; |
|
208 /* walk through the list of columns (structural variables) */ |
|
209 for (j = 1; j <= n; j++) |
|
210 { /* if j-th column is not marked as fractional, skip it */ |
|
211 if (!non_int[j]) continue; |
|
212 /* obtain (fractional) value of j-th column in basic solution |
|
213 of LP relaxation */ |
|
214 x = glp_get_col_prim(mip, j); |
|
215 /* since the value of j-th column is fractional, the column is |
|
216 basic; compute corresponding row of the simplex table */ |
|
217 len = glp_eval_tab_row(mip, m+j, ind, val); |
|
218 /* the following fragment computes a change in the objective |
|
219 function: delta Z = new Z - old Z, where old Z is the |
|
220 objective value in the current optimal basis, and new Z is |
|
221 the objective value in the adjacent basis, for two cases: |
|
222 1) if new upper bound ub' = floor(x[j]) is introduced for |
|
223 j-th column (down branch); |
|
224 2) if new lower bound lb' = ceil(x[j]) is introduced for |
|
225 j-th column (up branch); |
|
226 since in both cases the solution remaining dual feasible |
|
227 becomes primal infeasible, one implicit simplex iteration |
|
228 is performed to determine the change delta Z; |
|
229 it is obvious that new Z, which is never better than old Z, |
|
230 is a lower (minimization) or upper (maximization) bound of |
|
231 the objective function for down- and up-branches. */ |
|
232 for (kase = -1; kase <= +1; kase += 2) |
|
233 { /* if kase < 0, the new upper bound of x[j] is introduced; |
|
234 in this case x[j] should decrease in order to leave the |
|
235 basis and go to its new upper bound */ |
|
236 /* if kase > 0, the new lower bound of x[j] is introduced; |
|
237 in this case x[j] should increase in order to leave the |
|
238 basis and go to its new lower bound */ |
|
239 /* apply the dual ratio test in order to determine which |
|
240 auxiliary or structural variable should enter the basis |
|
241 to keep dual feasibility */ |
|
242 k = glp_dual_rtest(mip, len, ind, val, kase, 1e-9); |
|
243 if (k != 0) k = ind[k]; |
|
244 /* if no non-basic variable has been chosen, LP relaxation |
|
245 of corresponding branch being primal infeasible and dual |
|
246 unbounded has no primal feasible solution; in this case |
|
247 the change delta Z is formally set to infinity */ |
|
248 if (k == 0) |
|
249 { delta_z = |
|
250 (T->mip->dir == GLP_MIN ? +DBL_MAX : -DBL_MAX); |
|
251 goto skip; |
|
252 } |
|
253 /* row of the simplex table that corresponds to non-basic |
|
254 variable x[k] choosen by the dual ratio test is: |
|
255 x[j] = ... + alfa * x[k] + ... |
|
256 where alfa is the influence coefficient (an element of |
|
257 the simplex table row) */ |
|
258 /* determine the coefficient alfa */ |
|
259 for (t = 1; t <= len; t++) if (ind[t] == k) break; |
|
260 xassert(1 <= t && t <= len); |
|
261 alfa = val[t]; |
|
262 /* since in the adjacent basis the variable x[j] becomes |
|
263 non-basic, knowing its value in the current basis we can |
|
264 determine its change delta x[j] = new x[j] - old x[j] */ |
|
265 delta_j = (kase < 0 ? floor(x) : ceil(x)) - x; |
|
266 /* and knowing the coefficient alfa we can determine the |
|
267 corresponding change delta x[k] = new x[k] - old x[k], |
|
268 where old x[k] is a value of x[k] in the current basis, |
|
269 and new x[k] is a value of x[k] in the adjacent basis */ |
|
270 delta_k = delta_j / alfa; |
|
271 /* Tomlin noticed that if the variable x[k] is of integer |
|
272 kind, its change cannot be less (eventually) than one in |
|
273 the magnitude */ |
|
274 if (k > m && glp_get_col_kind(mip, k-m) != GLP_CV) |
|
275 { /* x[k] is structural integer variable */ |
|
276 if (fabs(delta_k - floor(delta_k + 0.5)) > 1e-3) |
|
277 { if (delta_k > 0.0) |
|
278 delta_k = ceil(delta_k); /* +3.14 -> +4 */ |
|
279 else |
|
280 delta_k = floor(delta_k); /* -3.14 -> -4 */ |
|
281 } |
|
282 } |
|
283 /* now determine the status and reduced cost of x[k] in the |
|
284 current basis */ |
|
285 if (k <= m) |
|
286 { stat = glp_get_row_stat(mip, k); |
|
287 dk = glp_get_row_dual(mip, k); |
|
288 } |
|
289 else |
|
290 { stat = glp_get_col_stat(mip, k-m); |
|
291 dk = glp_get_col_dual(mip, k-m); |
|
292 } |
|
293 /* if the current basis is dual degenerate, some reduced |
|
294 costs which are close to zero may have wrong sign due to |
|
295 round-off errors, so correct the sign of d[k] */ |
|
296 switch (T->mip->dir) |
|
297 { case GLP_MIN: |
|
298 if (stat == GLP_NL && dk < 0.0 || |
|
299 stat == GLP_NU && dk > 0.0 || |
|
300 stat == GLP_NF) dk = 0.0; |
|
301 break; |
|
302 case GLP_MAX: |
|
303 if (stat == GLP_NL && dk > 0.0 || |
|
304 stat == GLP_NU && dk < 0.0 || |
|
305 stat == GLP_NF) dk = 0.0; |
|
306 break; |
|
307 default: |
|
308 xassert(T != T); |
|
309 } |
|
310 /* now knowing the change of x[k] and its reduced cost d[k] |
|
311 we can compute the corresponding change in the objective |
|
312 function delta Z = new Z - old Z = d[k] * delta x[k]; |
|
313 note that due to Tomlin's modification new Z can be even |
|
314 worse than in the adjacent basis */ |
|
315 delta_z = dk * delta_k; |
|
316 skip: /* new Z is never better than old Z, therefore the change |
|
317 delta Z is always non-negative (in case of minimization) |
|
318 or non-positive (in case of maximization) */ |
|
319 switch (T->mip->dir) |
|
320 { case GLP_MIN: xassert(delta_z >= 0.0); break; |
|
321 case GLP_MAX: xassert(delta_z <= 0.0); break; |
|
322 default: xassert(T != T); |
|
323 } |
|
324 /* save the change in the objective fnction for down- and |
|
325 up-branches, respectively */ |
|
326 if (kase < 0) dz_dn = delta_z; else dz_up = delta_z; |
|
327 } |
|
328 /* thus, in down-branch no integer feasible solution can be |
|
329 better than Z + dz_dn, and in up-branch no integer feasible |
|
330 solution can be better than Z + dz_up, where Z is value of |
|
331 the objective function in the current basis */ |
|
332 /* following the heuristic by Driebeck and Tomlin we choose a |
|
333 column (i.e. structural variable) which provides largest |
|
334 degradation of the objective function in some of branches; |
|
335 besides, we select the branch with smaller degradation to |
|
336 be solved next and keep other branch with larger degradation |
|
337 in the active list hoping to minimize the number of further |
|
338 backtrackings */ |
|
339 if (degrad < fabs(dz_dn) || degrad < fabs(dz_up)) |
|
340 { jj = j; |
|
341 if (fabs(dz_dn) < fabs(dz_up)) |
|
342 { /* select down branch to be solved next */ |
|
343 next = GLP_DN_BRNCH; |
|
344 degrad = fabs(dz_up); |
|
345 } |
|
346 else |
|
347 { /* select up branch to be solved next */ |
|
348 next = GLP_UP_BRNCH; |
|
349 degrad = fabs(dz_dn); |
|
350 } |
|
351 /* save the objective changes for printing */ |
|
352 dd_dn = dz_dn, dd_up = dz_up; |
|
353 /* if down- or up-branch has no feasible solution, we does |
|
354 not need to consider other candidates (in principle, the |
|
355 corresponding branch could be pruned right now) */ |
|
356 if (degrad == DBL_MAX) break; |
|
357 } |
|
358 } |
|
359 /* free working arrays */ |
|
360 xfree(ind); |
|
361 xfree(val); |
|
362 /* something must be chosen */ |
|
363 xassert(1 <= jj && jj <= n); |
|
364 #if 1 /* 02/XI-2009 */ |
|
365 if (degrad < 1e-6 * (1.0 + 0.001 * fabs(mip->obj_val))) |
|
366 { jj = branch_mostf(T, &next); |
|
367 goto done; |
|
368 } |
|
369 #endif |
|
370 if (T->parm->msg_lev >= GLP_MSG_DBG) |
|
371 { xprintf("branch_drtom: column %d chosen to branch on\n", jj); |
|
372 if (fabs(dd_dn) == DBL_MAX) |
|
373 xprintf("branch_drtom: down-branch is infeasible\n"); |
|
374 else |
|
375 xprintf("branch_drtom: down-branch bound is %.9e\n", |
|
376 lpx_get_obj_val(mip) + dd_dn); |
|
377 if (fabs(dd_up) == DBL_MAX) |
|
378 xprintf("branch_drtom: up-branch is infeasible\n"); |
|
379 else |
|
380 xprintf("branch_drtom: up-branch bound is %.9e\n", |
|
381 lpx_get_obj_val(mip) + dd_up); |
|
382 } |
|
383 done: *_next = next; |
|
384 return jj; |
|
385 } |
|
386 |
|
387 /**********************************************************************/ |
|
388 |
|
389 struct csa |
|
390 { /* common storage area */ |
|
391 int *dn_cnt; /* int dn_cnt[1+n]; */ |
|
392 /* dn_cnt[j] is the number of subproblems, whose LP relaxations |
|
393 have been solved and which are down-branches for variable x[j]; |
|
394 dn_cnt[j] = 0 means the down pseudocost is uninitialized */ |
|
395 double *dn_sum; /* double dn_sum[1+n]; */ |
|
396 /* dn_sum[j] is the sum of per unit degradations of the objective |
|
397 over all dn_cnt[j] subproblems */ |
|
398 int *up_cnt; /* int up_cnt[1+n]; */ |
|
399 /* up_cnt[j] is the number of subproblems, whose LP relaxations |
|
400 have been solved and which are up-branches for variable x[j]; |
|
401 up_cnt[j] = 0 means the up pseudocost is uninitialized */ |
|
402 double *up_sum; /* double up_sum[1+n]; */ |
|
403 /* up_sum[j] is the sum of per unit degradations of the objective |
|
404 over all up_cnt[j] subproblems */ |
|
405 }; |
|
406 |
|
407 void *ios_pcost_init(glp_tree *tree) |
|
408 { /* initialize working data used on pseudocost branching */ |
|
409 struct csa *csa; |
|
410 int n = tree->n, j; |
|
411 csa = xmalloc(sizeof(struct csa)); |
|
412 csa->dn_cnt = xcalloc(1+n, sizeof(int)); |
|
413 csa->dn_sum = xcalloc(1+n, sizeof(double)); |
|
414 csa->up_cnt = xcalloc(1+n, sizeof(int)); |
|
415 csa->up_sum = xcalloc(1+n, sizeof(double)); |
|
416 for (j = 1; j <= n; j++) |
|
417 { csa->dn_cnt[j] = csa->up_cnt[j] = 0; |
|
418 csa->dn_sum[j] = csa->up_sum[j] = 0.0; |
|
419 } |
|
420 return csa; |
|
421 } |
|
422 |
|
423 static double eval_degrad(glp_prob *P, int j, double bnd) |
|
424 { /* compute degradation of the objective on fixing x[j] at given |
|
425 value with a limited number of dual simplex iterations */ |
|
426 /* this routine fixes column x[j] at specified value bnd, |
|
427 solves resulting LP, and returns a lower bound to degradation |
|
428 of the objective, degrad >= 0 */ |
|
429 glp_prob *lp; |
|
430 glp_smcp parm; |
|
431 int ret; |
|
432 double degrad; |
|
433 /* the current basis must be optimal */ |
|
434 xassert(glp_get_status(P) == GLP_OPT); |
|
435 /* create a copy of P */ |
|
436 lp = glp_create_prob(); |
|
437 glp_copy_prob(lp, P, 0); |
|
438 /* fix column x[j] at specified value */ |
|
439 glp_set_col_bnds(lp, j, GLP_FX, bnd, bnd); |
|
440 /* try to solve resulting LP */ |
|
441 glp_init_smcp(&parm); |
|
442 parm.msg_lev = GLP_MSG_OFF; |
|
443 parm.meth = GLP_DUAL; |
|
444 parm.it_lim = 30; |
|
445 parm.out_dly = 1000; |
|
446 parm.meth = GLP_DUAL; |
|
447 ret = glp_simplex(lp, &parm); |
|
448 if (ret == 0 || ret == GLP_EITLIM) |
|
449 { if (glp_get_prim_stat(lp) == GLP_NOFEAS) |
|
450 { /* resulting LP has no primal feasible solution */ |
|
451 degrad = DBL_MAX; |
|
452 } |
|
453 else if (glp_get_dual_stat(lp) == GLP_FEAS) |
|
454 { /* resulting basis is optimal or at least dual feasible, |
|
455 so we have the correct lower bound to degradation */ |
|
456 if (P->dir == GLP_MIN) |
|
457 degrad = lp->obj_val - P->obj_val; |
|
458 else if (P->dir == GLP_MAX) |
|
459 degrad = P->obj_val - lp->obj_val; |
|
460 else |
|
461 xassert(P != P); |
|
462 /* degradation cannot be negative by definition */ |
|
463 /* note that the lower bound to degradation may be close |
|
464 to zero even if its exact value is zero due to round-off |
|
465 errors on computing the objective value */ |
|
466 if (degrad < 1e-6 * (1.0 + 0.001 * fabs(P->obj_val))) |
|
467 degrad = 0.0; |
|
468 } |
|
469 else |
|
470 { /* the final basis reported by the simplex solver is dual |
|
471 infeasible, so we cannot determine a non-trivial lower |
|
472 bound to degradation */ |
|
473 degrad = 0.0; |
|
474 } |
|
475 } |
|
476 else |
|
477 { /* the simplex solver failed */ |
|
478 degrad = 0.0; |
|
479 } |
|
480 /* delete the copy of P */ |
|
481 glp_delete_prob(lp); |
|
482 return degrad; |
|
483 } |
|
484 |
|
485 void ios_pcost_update(glp_tree *tree) |
|
486 { /* update history information for pseudocost branching */ |
|
487 /* this routine is called every time when LP relaxation of the |
|
488 current subproblem has been solved to optimality with all lazy |
|
489 and cutting plane constraints included */ |
|
490 int j; |
|
491 double dx, dz, psi; |
|
492 struct csa *csa = tree->pcost; |
|
493 xassert(csa != NULL); |
|
494 xassert(tree->curr != NULL); |
|
495 /* if the current subproblem is the root, skip updating */ |
|
496 if (tree->curr->up == NULL) goto skip; |
|
497 /* determine branching variable x[j], which was used in the |
|
498 parent subproblem to create the current subproblem */ |
|
499 j = tree->curr->up->br_var; |
|
500 xassert(1 <= j && j <= tree->n); |
|
501 /* determine the change dx[j] = new x[j] - old x[j], |
|
502 where new x[j] is a value of x[j] in optimal solution to LP |
|
503 relaxation of the current subproblem, old x[j] is a value of |
|
504 x[j] in optimal solution to LP relaxation of the parent |
|
505 subproblem */ |
|
506 dx = tree->mip->col[j]->prim - tree->curr->up->br_val; |
|
507 xassert(dx != 0.0); |
|
508 /* determine corresponding change dz = new dz - old dz in the |
|
509 objective function value */ |
|
510 dz = tree->mip->obj_val - tree->curr->up->lp_obj; |
|
511 /* determine per unit degradation of the objective function */ |
|
512 psi = fabs(dz / dx); |
|
513 /* update history information */ |
|
514 if (dx < 0.0) |
|
515 { /* the current subproblem is down-branch */ |
|
516 csa->dn_cnt[j]++; |
|
517 csa->dn_sum[j] += psi; |
|
518 } |
|
519 else /* dx > 0.0 */ |
|
520 { /* the current subproblem is up-branch */ |
|
521 csa->up_cnt[j]++; |
|
522 csa->up_sum[j] += psi; |
|
523 } |
|
524 skip: return; |
|
525 } |
|
526 |
|
527 void ios_pcost_free(glp_tree *tree) |
|
528 { /* free working area used on pseudocost branching */ |
|
529 struct csa *csa = tree->pcost; |
|
530 xassert(csa != NULL); |
|
531 xfree(csa->dn_cnt); |
|
532 xfree(csa->dn_sum); |
|
533 xfree(csa->up_cnt); |
|
534 xfree(csa->up_sum); |
|
535 xfree(csa); |
|
536 tree->pcost = NULL; |
|
537 return; |
|
538 } |
|
539 |
|
540 static double eval_psi(glp_tree *T, int j, int brnch) |
|
541 { /* compute estimation of pseudocost of variable x[j] for down- |
|
542 or up-branch */ |
|
543 struct csa *csa = T->pcost; |
|
544 double beta, degrad, psi; |
|
545 xassert(csa != NULL); |
|
546 xassert(1 <= j && j <= T->n); |
|
547 if (brnch == GLP_DN_BRNCH) |
|
548 { /* down-branch */ |
|
549 if (csa->dn_cnt[j] == 0) |
|
550 { /* initialize down pseudocost */ |
|
551 beta = T->mip->col[j]->prim; |
|
552 degrad = eval_degrad(T->mip, j, floor(beta)); |
|
553 if (degrad == DBL_MAX) |
|
554 { psi = DBL_MAX; |
|
555 goto done; |
|
556 } |
|
557 csa->dn_cnt[j] = 1; |
|
558 csa->dn_sum[j] = degrad / (beta - floor(beta)); |
|
559 } |
|
560 psi = csa->dn_sum[j] / (double)csa->dn_cnt[j]; |
|
561 } |
|
562 else if (brnch == GLP_UP_BRNCH) |
|
563 { /* up-branch */ |
|
564 if (csa->up_cnt[j] == 0) |
|
565 { /* initialize up pseudocost */ |
|
566 beta = T->mip->col[j]->prim; |
|
567 degrad = eval_degrad(T->mip, j, ceil(beta)); |
|
568 if (degrad == DBL_MAX) |
|
569 { psi = DBL_MAX; |
|
570 goto done; |
|
571 } |
|
572 csa->up_cnt[j] = 1; |
|
573 csa->up_sum[j] = degrad / (ceil(beta) - beta); |
|
574 } |
|
575 psi = csa->up_sum[j] / (double)csa->up_cnt[j]; |
|
576 } |
|
577 else |
|
578 xassert(brnch != brnch); |
|
579 done: return psi; |
|
580 } |
|
581 |
|
582 static void progress(glp_tree *T) |
|
583 { /* display progress of pseudocost initialization */ |
|
584 struct csa *csa = T->pcost; |
|
585 int j, nv = 0, ni = 0; |
|
586 for (j = 1; j <= T->n; j++) |
|
587 { if (glp_ios_can_branch(T, j)) |
|
588 { nv++; |
|
589 if (csa->dn_cnt[j] > 0 && csa->up_cnt[j] > 0) ni++; |
|
590 } |
|
591 } |
|
592 xprintf("Pseudocosts initialized for %d of %d variables\n", |
|
593 ni, nv); |
|
594 return; |
|
595 } |
|
596 |
|
597 int ios_pcost_branch(glp_tree *T, int *_next) |
|
598 { /* choose branching variable with pseudocost branching */ |
|
599 glp_long t = xtime(); |
|
600 int j, jjj, sel; |
|
601 double beta, psi, d1, d2, d, dmax; |
|
602 /* initialize the working arrays */ |
|
603 if (T->pcost == NULL) |
|
604 T->pcost = ios_pcost_init(T); |
|
605 /* nothing has been chosen so far */ |
|
606 jjj = 0, dmax = -1.0; |
|
607 /* go through the list of branching candidates */ |
|
608 for (j = 1; j <= T->n; j++) |
|
609 { if (!glp_ios_can_branch(T, j)) continue; |
|
610 /* determine primal value of x[j] in optimal solution to LP |
|
611 relaxation of the current subproblem */ |
|
612 beta = T->mip->col[j]->prim; |
|
613 /* estimate pseudocost of x[j] for down-branch */ |
|
614 psi = eval_psi(T, j, GLP_DN_BRNCH); |
|
615 if (psi == DBL_MAX) |
|
616 { /* down-branch has no primal feasible solution */ |
|
617 jjj = j, sel = GLP_DN_BRNCH; |
|
618 goto done; |
|
619 } |
|
620 /* estimate degradation of the objective for down-branch */ |
|
621 d1 = psi * (beta - floor(beta)); |
|
622 /* estimate pseudocost of x[j] for up-branch */ |
|
623 psi = eval_psi(T, j, GLP_UP_BRNCH); |
|
624 if (psi == DBL_MAX) |
|
625 { /* up-branch has no primal feasible solution */ |
|
626 jjj = j, sel = GLP_UP_BRNCH; |
|
627 goto done; |
|
628 } |
|
629 /* estimate degradation of the objective for up-branch */ |
|
630 d2 = psi * (ceil(beta) - beta); |
|
631 /* determine d = max(d1, d2) */ |
|
632 d = (d1 > d2 ? d1 : d2); |
|
633 /* choose x[j] which provides maximal estimated degradation of |
|
634 the objective either in down- or up-branch */ |
|
635 if (dmax < d) |
|
636 { dmax = d; |
|
637 jjj = j; |
|
638 /* continue the search from a subproblem, where degradation |
|
639 is less than in other one */ |
|
640 sel = (d1 <= d2 ? GLP_DN_BRNCH : GLP_UP_BRNCH); |
|
641 } |
|
642 /* display progress of pseudocost initialization */ |
|
643 if (T->parm->msg_lev >= GLP_ON) |
|
644 { if (xdifftime(xtime(), t) >= 10.0) |
|
645 { progress(T); |
|
646 t = xtime(); |
|
647 } |
|
648 } |
|
649 } |
|
650 if (dmax == 0.0) |
|
651 { /* no degradation is indicated; choose a variable having most |
|
652 fractional value */ |
|
653 jjj = branch_mostf(T, &sel); |
|
654 } |
|
655 done: *_next = sel; |
|
656 return jjj; |
|
657 } |
|
658 |
|
659 /* eof */ |