1 /*===========================================================================*/
2 /* */
3 /* This file is part of the SYMPHONY MILP Solver Framework. */
4 /* */
5 /* SYMPHONY was jointly developed by Ted Ralphs (ted@lehigh.edu) and */
6 /* Laci Ladanyi (ladanyi@us.ibm.com). */
7 /* */
8 /* (c) Copyright 2000-2019 Ted Ralphs. All Rights Reserved. */
9 /* */
10 /* This software is licensed under the Eclipse Public License. Please see */
11 /* accompanying file for terms. */
12 /* */
13 /*===========================================================================*/
14
15 #include <math.h>
16 #include <memory.h>
17 #include <stdlib.h>
18
19 #include "sym_lp.h"
20 #include "sym_qsort.h"
21 #include "sym_proccomm.h"
22 #include "sym_messages.h"
23 #include "sym_constants.h"
24 #include "sym_macros.h"
25 #include "sym_types.h"
26 #include "sym_lp_solver.h"
27
28 /*===========================================================================*/
29
30 /*===========================================================================*\
31 * This file contains LP functions related to branching.
32 \*===========================================================================*/
33
add_slacks_to_matrix(lp_prob * p,int cand_num,branch_obj ** candidates)34 void add_slacks_to_matrix(lp_prob *p, int cand_num, branch_obj **candidates)
35 {
36 LPdata *lp_data = p->lp_data;
37 int *index;
38 int m = p->lp_data->m;
39 int j, k;
40 branch_obj *can;
41 row_data *newrows;
42 waiting_row **wrows;
43
44 for (j=cand_num-1; j >= 0; j--)
45 if (candidates[j]->type == CANDIDATE_CUT_NOT_IN_MATRIX)
46 break;
47
48 if (j < 0) /* there is nothing to add */
49 return;
50
51 /* We'll create a waiting row for each cut, add them to the matrix,
52 then set their status to be free */
53 wrows = (waiting_row **) malloc(cand_num * sizeof(waiting_row *));
54 /* can't use tmp.p, because that might get resized in add_row_set */
55 for (k=0; j >= 0; j--){
56 can = candidates[j];
57 if (can->type == CANDIDATE_CUT_NOT_IN_MATRIX){
58 wrows[k] = can->row;
59 can->row = NULL;
60 can->position = m + k;
61 can->type = CANDIDATE_CUT_IN_MATRIX;
62 k++;
63 }
64 }
65 add_row_set(p, wrows, k);
66 /* To satisfy the size requirements in free_row_set, the following sizes
67 * are needed: tmp.c:2*m tmp.i1:3*m tmp.d:m */
68 FREE(wrows);
69 index = lp_data->tmp.i1;
70 for (j = 0; j < k; j++)
71 index[j] = m + j;
72 free_row_set(lp_data, k, index);
73 newrows = lp_data->rows + m; /* m is still the old one! */
74 for (j = 0; j < k; j++){
75 newrows[j].ineff_cnt = (MAXINT) >> 1; /* it is slack... */
76 newrows[j].free = TRUE;
77 }
78 }
79
80 /*===========================================================================*/
81
82 /*===========================================================================*\
83 * Ok. So there were violated cuts, either in waiting_rows or among the
84 * slacks (or both). We just have to add those among the slacks to the
85 * waiting_rows; add the appropriate number of cuts to the problem and
86 * continue the node (this second part is just like the end of receive_cuts().
87 \*===========================================================================*/
88
add_violated_slacks(lp_prob * p,int cand_num,branch_obj ** candidates)89 int add_violated_slacks(lp_prob *p, int cand_num, branch_obj **candidates)
90 {
91 LPdata *lp_data = p->lp_data;
92 waiting_row **new_rows;
93 int i, new_row_num = 0;
94
95 /* If there are any violated (former) slack, unpack them and add them
96 * to the set of waiting rows. */
97 if (cand_num > 0){
98 new_rows = (waiting_row **) lp_data->tmp.p1; /* m (actually, candnum<m */
99 for (i=0; i<cand_num; i++){
100 if (candidates[i]->type == VIOLATED_SLACK){
101 new_rows[new_row_num++] = candidates[i]->row;
102 candidates[i]->row = NULL;
103 }
104 }
105 if (new_row_num > 0)
106 add_new_rows_to_waiting_rows(p, new_rows, new_row_num);
107 }
108
109 return( p->waiting_row_num == 0 ? 0 : add_best_waiting_rows(p) );
110 }
111
112 /*===========================================================================*/
113
select_branching_object(lp_prob * p,int * cuts,branch_obj ** candidate)114 int select_branching_object(lp_prob *p, int *cuts, branch_obj **candidate)
115 {
116
117 LPdata *lp_data = p->lp_data;
118 var_desc **vars;
119 row_data *rows;
120 int m;
121 #ifndef MAX_CHILDREN_NUM
122 int maxnum;
123 double *objval, *pobj;
124 int *termcode, *pterm, *feasible, *pfeas, *iterd, *piter;
125 #ifdef COMPILE_FRAC_BRANCHING
126 int *frnum, *pfrnum, **frind, **pfrind;
127 double **frval, **pfrval;
128 #endif
129 #endif
130 int i, j, k, l, branch_var, branch_row, min_ind;
131 double lb, ub, oldobjval, min_obj;
132 cut_data *cut;
133 branch_obj *can, *best_can = NULL;
134 #ifdef COMPILE_FRAC_BRANCHING
135 int *xind;
136 double *xval;
137 #endif
138
139 /* These are the return values from select_candidates_u() */
140 int cand_num = 0, new_vars = 0, *indices;
141 double *values;
142 branch_obj **candidates = NULL;
143 #ifdef STATISTICS
144 int itlim = 0, cnum = 0;
145 #endif
146 int should_use_hot_starts = FALSE, unmark_hs = TRUE;
147 double total_time = 0;
148 double st_time = 0;
149 int total_iters, should_continue;//, max_iter_num;
150 int should_use_rel_br = p->par.should_use_rel_br;
151 double high, low, down_obj, up_obj, best_var_score;
152 int *br_rel_down = p->br_rel_down;
153 int *br_rel_up = p->br_rel_up;
154 double *pcost_down = p->pcost_down;
155 double *pcost_up = p->pcost_up;
156 int rel_threshold = p->par.rel_br_threshold;
157 int best_var;
158 // double lp_time, scaled_by;
159 int max_presolve_iter = 5;
160 const int bc_level = p->bc_level;
161 int strong_br_min_level = p->par.strong_br_min_level;
162
163 used_time(&total_time);
164
165 /*---------------------------------------------------------------------* \
166 * First we call select_candidates_u() to select candidates. It can
167 * -- return with DO_BRANCH and a bunch of candidates, or
168 * -- return with DO_NOT_BRANCH along with a bunch of violated cuts
169 * in the matrix and/or among the slack_cuts, or
170 * -- return with DO_NOT_BRANCH__FATHOMED, i.e., the node can be fathomed.
171 \*------------------------------------------------------------------------*/
172
173 j = select_candidates_u(p, cuts, &new_vars, &cand_num, &candidates);
174
175 switch (j){
176 case DO_NOT_BRANCH__FATHOMED:
177 *candidate = NULL;
178 return(DO_NOT_BRANCH__FATHOMED);
179 case DO_NOT_BRANCH__FEAS_SOL:
180 *candidate = NULL;
181 return(DO_NOT_BRANCH__FEAS_SOL);
182 case DO_NOT_BRANCH:
183 if (cand_num)
184 *cuts += add_violated_slacks(p, cand_num, candidates);
185 /* Free the candidates */
186 if (candidates){
187 for (i=0; i<cand_num; i++){
188 free_candidate(candidates + i);
189 }
190 FREE(candidates);
191 }
192 *candidate = NULL;
193 return(DO_NOT_BRANCH);
194
195 case DO_BRANCH:
196 break;
197
198 case ERROR__NO_BRANCHING_CANDIDATE:
199 *candidate = NULL;
200 return(ERROR__NO_BRANCHING_CANDIDATE);
201 }
202
203 /* OK, now we have to branch. */
204
205 /* First of all, send everything to the cutpool that hasn't been sent
206 before and send the current node description to the TM. */
207 p->comp_times.strong_branching += used_time(&p->tt);
208 #pragma omp critical(cut_pool)
209 send_cuts_to_pool(p, -1);
210 p->comp_times.communication += used_time(&p->tt);
211
212 /* Add all the branching cuts */
213 if (p->par.branch_on_cuts)
214 add_slacks_to_matrix(p, cand_num, candidates);
215 m = lp_data->m;
216 rows = lp_data->rows;
217
218
219 #ifndef MAX_CHILDREN_NUM
220 /* The part below is not needed when we have MAX_CHILDREN_NUM specified */
221 /* Count how many objval/termcode/feasible entry we might need
222 and allocate space for it */
223 for (maxnum = candidates[0]->child_num, j=0, i=1; i<cand_num; i++){
224 if (maxnum < candidates[i]->child_num)
225 maxnum = candidates[i]->child_num;
226 }
227
228 objval = (double *) malloc(maxnum * DSIZE);
229 termcode = (int *) malloc(maxnum * ISIZE);
230 feasible = (int *) malloc(maxnum * ISIZE);
231 iterd = (int *) malloc(maxnum * ISIZE);
232 #ifdef COMPILE_FRAC_BRANCHING
233 frval = (double **) malloc(maxnum * sizeof(double *));
234 pfrval = (double **) malloc(maxnum * sizeof(double *));
235 frind = (int **) malloc(maxnum * sizeof(int *));
236 pfrind = (int **) malloc(maxnum * sizeof(int *));
237 frnum = (int *) malloc(maxnum * ISIZE);
238 pfrnum = (int *) malloc(maxnum * ISIZE);
239 #endif
240 pobj = (double *) malloc(maxnum * DSIZE);
241 pterm = (int *) malloc(maxnum * ISIZE);
242 pfeas = (int *) malloc(maxnum * ISIZE);
243 piter = (int *) malloc(maxnum * ISIZE);
244 #endif
245
246 /* Look at the candidates one-by-one and presolve them. */
247 vars = lp_data->vars;
248 oldobjval = lp_data->objval;
249 st_time = used_time(&total_time);
250 total_iters = 0;
251
252 int *cstat = lp_data->tmp.i1;
253 int *rstat = lp_data->tmp.i2;
254
255 get_basis(lp_data, cstat, rstat);
256
257 if (should_use_rel_br==TRUE) {
258
259 const double lpetol100 = lp_data->lpetol*100;
260 double lpetol = lp_data->lpetol;
261
262 if(!(lp_data->tmp2_size) || lp_data->tmp2_size < 2*lp_data->n){
263 FREE(lp_data->tmp2.i1);
264 FREE(lp_data->tmp2.d);
265 FREE(lp_data->tmp2.c);
266 int tmp_size = 2*lp_data->n;
267 lp_data->tmp2.i1 = (int *)malloc (tmp_size*ISIZE);
268 lp_data->tmp2.d = (double *)malloc ((tmp_size + lp_data->n)*DSIZE);
269 lp_data->tmp2.c = (char *)malloc (tmp_size*CSIZE);
270 }
271 double *bnd_val = lp_data->tmp2.d; //(double *)malloc (2*lp_data->n*DSIZE);
272 int *bnd_ind = lp_data->tmp2.i1; //(int *)malloc (2*lp_data->n*ISIZE);
273 char *bnd_sense = lp_data->tmp2.c; //(char *)malloc (2*lp_data->n*CSIZE);
274 int *up_violation_cnt = NULL, *down_violation_cnt = NULL;
275 int *violation_col_size = NULL;
276 int num_bnd_changes = 0;
277 double xval, floorx, ceilx, var_score;
278 int full_solves = 0, down_is_est, up_is_est, best_down_is_est,
279 best_up_is_est,
280 max_solves_since_impr = p->par.rel_br_cand_threshold,
281 stop_solving = FALSE, both_children_inf = FALSE, rel_up,
282 rel_down, solves_since_impr = 0, best_one_child_inf = FALSE;
283 int max_solves = p->par.rel_br_max_solves;
284 double alpha = p->par.strong_branching_high_low_weight;
285 double one_m_alpha = 1.0 - alpha;
286 int check_first = FALSE;
287 int check_level = 0;
288 int num_up_iters = 0, num_down_iters = 0;
289 int up_status = -1, down_status = -1;
290
291 int check_off = TRUE;
292
293 double *row_lb = lp_data->tmp.d;
294 double *row_ub = lp_data->tmp.d + lp_data->m;
295 char cand_fixed = FALSE;
296
297
298 // experimental - node-presolve
299 if (p->par.use_branching_prep && cand_num > 1){
300 //prep_tighten_bounds(lp_data, &num_bnd_changes, bnd_val, bnd_ind, bnd_sense,
301 // row_ub, row_lb, &cand_fixed);
302
303 if(prep_tighten_bounds(lp_data, &num_bnd_changes, bnd_val, bnd_ind, bnd_sense,
304 row_ub, row_lb, &cand_fixed) == PREP_INFEAS){
305 cand_num = 1;
306 FREE(bnd_val);
307 FREE(bnd_ind);
308 FREE(bnd_sense);
309 FREE(best_can);
310 FREE(candidates);
311 *candidate = NULL;
312 p->lp_stat.prep_nodes_pruned++;
313 set_itlim(lp_data, -1); //both limits should be set for hotstarts
314 return (DO_NOT_BRANCH__FATHOMED);
315 }else if(num_bnd_changes > 0){
316 p->lp_stat.prep_bnd_changes += num_bnd_changes;
317 if(cand_fixed){
318 int c_ind, new_cand_num = 0;
319 int *new_cand_list = lp_data->tmp.i1;
320 for(i = 0; i< cand_num; i++){
321 c_ind = p->br_rel_cand_list[i];
322 xval = lp_data->x[c_ind];
323 if(vars[c_ind]->lb < xval &&
324 vars[c_ind]->ub > xval){
325 new_cand_list[new_cand_num] = c_ind;
326 new_cand_num++;
327 }
328 }
329 cand_num = new_cand_num;
330 memcpy(p->br_rel_cand_list, new_cand_list, ISIZE*cand_num);
331 }
332 }
333
334
335 if(p->par.use_branching_prep){//use_violation){ //}
336 up_violation_cnt = (int *)calloc (lp_data->n,ISIZE);
337 down_violation_cnt = (int *)calloc (lp_data->n,ISIZE);
338 violation_col_size = (int *)calloc(lp_data->n, ISIZE);
339
340
341 const double *si_ub = lp_data->si->getColUpper();
342 const double *si_lb = lp_data->si->getColLower();
343
344 const CoinPackedMatrix * matrix = lp_data->si->getMatrixByCol();
345 const double *matval = matrix->getElements();
346 const int *matind = matrix->getIndices();
347 const int *matbeg = matrix->getVectorStarts();
348 const int *len = matrix->getVectorLengths();
349
350
351 int c_ind, r_ind, col_start, col_end;
352 double coeff;
353
354 const double *r_ub = lp_data->si->getRowUpper();
355 const double *r_lb = lp_data->si->getRowLower();
356 const double inf = lp_data->si->getInfinity();
357 //double new_objval = 0;
358 //double *violation_max_cnt = lp_data->tmp.d + 2*lp_data->m;
359
360 const double *r_act = lp_data->si->getRowActivity();
361 double up_max, down_max, new_act, new_row_lb, new_row_ub;
362 double violation, si_row_ub, si_row_lb;
363 for(i = 0; i < cand_num; i++){
364 up_max = -DBL_MAX;
365 down_max = -DBL_MAX;
366 c_ind = p->br_rel_cand_list[i];
367 col_start = matbeg[c_ind];
368 col_end = col_start + len[c_ind];
369 xval = lp_data->x[c_ind];
370 floorx = floor(xval);
371 ceilx = ceil(xval);
372
373 for(j = col_start; j < col_end; j++){
374 char get_cols_dir = 'R';
375 r_ind = matind[j];
376 coeff = matval[j];
377 if(row_ub[r_ind] < row_lb[r_ind] + 100*lp_data->lpetol)
378 printf("error in row bounds...%i %i %f %f\n", p->bc_index, r_ind, row_lb[r_ind], row_ub[r_ind]);
379 if(row_ub[r_ind] < row_lb[r_ind] + lp_data->lpetol) continue;
380 si_row_ub = r_ub[r_ind];
381 si_row_lb = r_lb[r_ind];
382 violation_col_size[c_ind]++;
383 if(si_row_ub < si_row_lb + lp_data->lpetol){ // 'E'
384 get_cols_dir = 'E';
385 //up_violation_cnt[c_ind]++;
386 //down_violation_cnt[c_ind]++;
387 if(coeff >= 0.0){
388 /* fixing to upper */
389 new_row_lb = row_lb[r_ind] + coeff*(ceilx - si_lb[c_ind]);
390 /* fixing to lower */
391 new_row_ub = row_ub[r_ind] + coeff*(floorx - si_ub[c_ind]);
392 if(new_row_lb > si_row_ub - lp_data->lpetol){
393 get_cols_dir = 'U';
394 up_violation_cnt[c_ind]++;
395 }else if(new_row_ub < si_row_lb + lp_data->lpetol){
396 get_cols_dir = 'D';
397 down_violation_cnt[c_ind]++;
398 }else{
399 up_violation_cnt[c_ind]++;
400 down_violation_cnt[c_ind]++;
401 }
402 }else{
403 /* fixing to upper */
404 new_row_ub = row_ub[r_ind] + coeff*(ceilx - si_lb[c_ind]);
405 /* fixing to lower */
406 new_row_lb = row_lb[r_ind] + coeff*(floorx - si_ub[c_ind]);
407
408 if(new_row_ub < si_row_lb + lp_data->lpetol){
409 get_cols_dir = 'U';
410 up_violation_cnt[c_ind]++;
411 }else if (new_row_lb > si_row_ub - lp_data->lpetol){
412 get_cols_dir = 'D';
413 down_violation_cnt[c_ind]++;
414 }
415 }
416 }else{
417 si_row_ub = MIN(r_ub[r_ind], inf/2);
418 si_row_lb = MAX(r_lb[r_ind], -inf/2);
419
420 new_act = r_act[r_ind] + coeff*(ceilx - xval);
421 violation = MAX(new_act - si_row_ub, si_row_lb - new_act);
422 //if(violation > up_max) up_max = violation;
423 if(violation > lp_data->lpetol) {
424 get_cols_dir = 'U';
425 up_violation_cnt[c_ind]++;
426 }
427
428 new_act = r_act[r_ind] + coeff*(floorx - xval);
429 violation = MAX(new_act - si_row_ub, si_row_lb - new_act);
430 //if(violation > down_max) down_max = violation;
431 if(violation > lp_data->lpetol) {
432 get_cols_dir = 'D';
433 down_violation_cnt[c_ind]++;
434 }
435 }
436 }
437 }
438 }
439 }
440
441 double *x = lp_data->tmp2.d + 2*(lp_data->n); //(double *)malloc (lp_data->n*DSIZE);
442
443 best_var = -1;
444 best_var_score = -SYM_INFINITY;
445 memcpy(x, lp_data->x, lp_data->n*DSIZE);
446
447
448 #ifdef COMPILE_IN_LP
449
450 if(p->par.rel_br_override_default && p->mip->mip_inf && cand_num > 1){
451
452
453 int weighted_iter =
454 p->lp_stat.lp_total_iter_num/(p->lp_stat.lp_calls -
455 p->lp_stat.str_br_lp_calls -
456 p->lp_stat.fp_lp_calls + 1);
457 if(p->mip->nz > 5e3){
458 weighted_iter = (int)
459 ((1.0*weighted_iter * p->mip->nz) / 5e3);
460 }
461
462 if(p->mip->nz > 5e4){
463 rel_threshold = MAX(2, (int)(1.0 * rel_threshold * 5e4/p->mip->nz));
464 }
465
466 if(p->bc_level < 1){
467 if(p->iter_num > 2 && weighted_iter <= 1000){
468 if(p->mip->mip_inf){
469 if(p->mip->mip_inf->prob_type == BINARY_TYPE){
470 strong_br_min_level =
471 MIN(p->par.strong_br_min_level,
472 (int)((p->mip->mip_inf->binary_var_num)/10.0) + 1);
473 }
474 if(p->mip->mip_inf->prob_type == BIN_CONT_TYPE){
475 if(p->mip->mip_inf->bin_var_ratio < 0.1){
476 strong_br_min_level =
477 MIN(MAX(p->par.strong_br_min_level,
478 (int)((p->mip->mip_inf->binary_var_num)/10.0)
479 + 1), 10);
480 if(p->mip->nz < 5e4){
481 strong_br_min_level = MAX(p->par.strong_br_min_level, 8);
482 }
483 }
484 }
485 }
486 }
487 }
488
489 if(weighted_iter * p->bc_index < 5e7){
490 //check_off = FALSE;
491 }
492
493
494 if(p->mip->mip_inf && p->mip->mip_inf->bin_cont_row_num > 0 &&
495 (p->mip->mip_inf->bin_cont_row_num >= p->mip->m ||
496 (p->mip->mip_inf->bin_var_ratio < 0.2) ||
497 p->mip->n - p->mip->mip_inf->cont_var_num <= 100 ||
498 (p->mip->mip_inf->sos_bin_row_ratio < 0.00 &&
499 p->mip->mip_inf->bin_var_ratio < 0.6) ||
500 (p->mip->mip_inf->max_row_ratio < 0.01 &&
501 p->mip->mip_inf->prob_type != BIN_CONT_TYPE))){
502 /* -either we have all continuos rows
503 -less number of bin vars
504 -small number of int vars
505 -large bin but less sos rows
506 -small max_row_size - if rest are all binary, skip to latter
507 */
508
509 if(p->mip->mip_inf->bin_cont_row_num >= p->mip->m){
510 max_solves = MIN(max_solves, 2*cand_num);
511
512 }else if(p->mip->mip_inf->bin_var_ratio < 0.2){
513 max_solves = MIN(max_solves, 2*cand_num);
514 if(p->mip->mip_inf->bin_var_ratio > 0.05){
515 strong_br_min_level = (int)strong_br_min_level/2;
516 }
517 }else{// if(p->mip->mip_inf->max_row_ratio < 0.01){ //}
518 max_solves = MIN(2*max_solves, 2*cand_num);
519 if(p->mip->mip_inf->sos_bin_row_ratio > 0.05){
520 // max_solves = MIN(2*max_solves, 2*cand_num);
521 strong_br_min_level = (int)(2.0*strong_br_min_level);
522 rel_threshold = 2*rel_threshold;
523 }
524 }
525 }else{
526 double imp_avg = 0.0;
527 int backtrack = 0;
528
529 bc_node *node = p->tm->active_nodes[p->proc_index];
530 if(p->bc_level >= 1){
531 while(node->parent){
532 if(node->start_objval > node->parent->end_objval){
533 imp_avg +=
534 fabs(node->start_objval/(node->parent->end_objval +0.0001) - 1.0);
535 }
536 node = node->parent;
537 if(backtrack++ > p->par.rel_br_chain_backtrack) break;
538 }
539 }
540
541 if(backtrack > 0){
542 imp_avg /= backtrack;
543 }
544
545 if(imp_avg > p->par.rel_br_min_imp &&
546 imp_avg < p->par.rel_br_max_imp){
547 if(bc_level <= strong_br_min_level ){
548 max_solves = MIN(3*max_solves, 2*cand_num);
549 }else{
550 max_solves = MIN(2*max_solves, 2*cand_num);
551 }
552 }else{
553
554 int c_cnt = 0;
555 double d_avg = 0.0;
556
557 for (i=0; i<cand_num; i++) {
558 branch_var = p->br_rel_cand_list[i];
559 xval = x[branch_var];
560 floorx = floor(xval);
561 ceilx = ceil(xval);
562 rel_down = br_rel_down[branch_var];
563 rel_up = br_rel_up[branch_var];
564 if(xval - floorx > 0.5){
565 d_avg += ceilx - xval;
566 }else{
567 d_avg += xval - floorx;
568 }
569 }
570
571 d_avg /= cand_num;
572
573 for (i=0; i<cand_num; i++) {
574 branch_var = p->br_rel_cand_list[i];
575 xval = x[branch_var];
576 if(xval - floor(xval) > d_avg && ceil(xval) - xval > d_avg) c_cnt++;
577 else break;
578 }
579
580 if(bc_level < 1){
581 max_solves = (cand_num < p->par.rel_br_override_max_solves ?
582 MIN(p->par.rel_br_override_max_solves/2, cand_num) :
583 p->par.rel_br_override_max_solves);
584 }else if(bc_level < 4){
585 max_solves = MIN((int)(0.75*c_cnt), (int)(0.3 * cand_num) + 1);
586 }else if(bc_level < 20){
587 max_solves = MIN(c_cnt/2, (int)(0.25 * cand_num) + 1);
588 }else if(bc_level < 40){
589 max_solves = MIN(c_cnt/3, (int)(0.20 * cand_num) + 1);
590 }else{
591 max_solves = MAX(c_cnt/4, (int)(0.15 * cand_num) + 1);
592 }
593 }
594
595 max_solves_since_impr = 5;
596
597 //printf("level - set to : %i %i\n", p->bc_level, max_solves);
598 //printf("c_cnt - cand num - max_solves : %i %i %i\n\n",
599 //c_cnt,cand_num, max_solves);
600
601 int int_num = p->mip->n - p->mip->mip_inf->cont_var_num;
602 int max_level = ((p->mip->mip_inf == 0) ? 500 :
603 (int_num)/2);
604 max_level = MIN(500, MAX(100, max_level));
605
606 if(cand_num > 100 && int_num > 500){
607 max_level = MIN(100, max_level);
608 if((p->mip->mip_inf->prob_type == BINARY_TYPE ||
609 p->mip->mip_inf->prob_type == BIN_CONT_TYPE) &&
610 cand_num > 0.05*int_num){
611 max_level /= 2;
612 }
613 }
614
615 if(bc_level > max_level){
616 rel_threshold = max_solves = 0;
617 strong_br_min_level = 1;
618 //cand_num = 1;
619 }
620 //printf("max_level: %i\n", max_level);
621 }
622
623 max_solves = MIN(p->par.rel_br_override_max_solves, max_solves);
624
625 double rel_limit = 0.05;
626 if((p->mip->mip_inf && ((p->mip->mip_inf->mat_density < rel_limit &&
627 p->mip->mip_inf->int_var_ratio > rel_limit &&
628 (p->mip->mip_inf->max_col_ratio > rel_limit ||
629 p->mip->mip_inf->max_row_ratio > rel_limit))||
630 (p->mip->nz > 1e5 && p->mip->mip_inf->mat_density > rel_limit/50) ||
631 (p->mip->mip_inf->max_row_ratio < rel_limit/5 &&
632 p->mip->mip_inf->prob_type != BIN_CONT_TYPE)))){
633 #ifdef __OSI_CLP__
634 lp_data->si->setupForRepeatedUse(2,0);
635 #endif
636 }
637
638 if(p->mip->mip_inf && !check_off &&
639 (p->mip->mip_inf->prob_type == BINARY_TYPE ||
640 p->mip->mip_inf->prob_type == BIN_CONT_TYPE) &&
641 (p->mip->n - p->mip->mip_inf->cont_var_num < 100 ||
642 (p->mip->mip_inf->int_var_ratio > rel_limit &&
643 p->mip->mip_inf->row_density/(p->mip->n + 1) > rel_limit/5))){//
644 }
645
646
647 if(p->mip->mip_inf->binary_sos_row_num > 0){
648 double bin_den = (1.0*p->mip->mip_inf->binary_sos_row_num)/
649 (p->mip->m + 1);
650 if( bin_den > rel_limit && ((bin_den < 10*rel_limit &&
651 p->mip->mip_inf->prob_type != BINARY_TYPE &&
652 p->mip->mip_inf->bin_var_ratio > 10*rel_limit) ||
653 (bin_den < 2*rel_limit &&
654 p->mip->mip_inf->prob_type == BINARY_TYPE))){
655 /* give priority to vars appear in sos rows */
656 int *sos_ind = lp_data->tmp.i1;//(int *)(malloc)(ISIZE*cand_num);
657 int *sos_tot_var = lp_data->tmp.i1+cand_num;//(int *)(malloc)(ISIZE*cand_num);
658 int sos_cnt = 0;
659 for (i=0; i<cand_num; i++) {
660 branch_var = p->br_rel_cand_list[i];
661 //printf("%i %i\n", branch_var, p->mip->mip_inf->cols[branch_var].sos_num);
662 //if(p->mip->mip_inf->cols[branch_var].sos_num > 0.1*p->mip->n){ //}
663 if(p->mip->mip_inf->cols[branch_var].sos_num >= (1.0*p->mip->nz)/(p->mip->m + 1)){
664 sos_tot_var[sos_cnt] = -p->mip->mip_inf->cols[branch_var].sos_num;
665 sos_ind[sos_cnt] = i;
666 sos_cnt++;
667 }
668 }
669 //printf("sos_cnt %i\n", sos_cnt);
670 if(sos_cnt > 0){
671 qsort_ii(sos_tot_var, sos_ind, sos_cnt);
672 int *sos_chosen = lp_data->tmp.i1+cand_num;//(int *)(calloc)(ISIZE,cand_num);
673 int *new_ord = lp_data->tmp.i1+2*cand_num;//(int *)(malloc)(ISIZE*cand_num);
674 memset(sos_chosen, 0, ISIZE*cand_num);
675 for (i=0; i<MIN(max_solves/2 + 1, sos_cnt); i++) {
676 new_ord[i] = p->br_rel_cand_list[sos_ind[i]];
677 sos_chosen[sos_ind[i]] = TRUE;
678 }
679
680 if(i < cand_num){
681 int rest_cnt = 0;
682 for(j = 0; j < cand_num; j++){
683 if(sos_chosen[j]) continue;
684 else new_ord[i+rest_cnt++] = p->br_rel_cand_list[j];
685 }
686 }
687
688 memcpy(p->br_rel_cand_list, new_ord, ISIZE*cand_num);
689 }
690 }
691 }
692 }
693
694 /* order by inf status */
695
696 update_solve_parameters(p);
697
698 if (p->mip->mip_inf){
699 if(1.0*p->mip->mip_inf->cont_var_num/(p->mip->n + 1) < 0.2 ||
700 1.0*p->mip->mip_inf->cont_var_num/(p->mip->n + 1) > 0.8){
701 if(p->bc_level <= 10){
702 max_solves *= 3;
703 max_solves_since_impr *= 2;
704 rel_threshold *=2;
705 }
706 }
707 }
708
709 #endif
710
711 if(cand_num > 1 && !p->par.disable_obj && !p->par.rs_mode_enabled){
712 if(p->par.use_hot_starts && !p->par.branch_on_cuts){
713 should_use_hot_starts = TRUE;
714 }else{
715 should_use_hot_starts = FALSE;
716 }
717
718 if (should_use_hot_starts) {
719 mark_hotstart(lp_data);
720 }
721 }
722
723 if (p->par.max_presolve_iter > 0) {
724 max_presolve_iter = p->par.max_presolve_iter - bc_level;
725
726 #ifdef COMPILE_IN_LP
727 if(p->mip->nz > 5e4){
728 max_presolve_iter = (int)(1.0 * max_presolve_iter * 5e4/p->mip->nz);
729 }
730 #endif
731 max_presolve_iter = MAX(max_presolve_iter, 40);
732 //max_presolve_iter = 40;
733 if(p->par.rs_mode_enabled) max_presolve_iter = 5;
734
735 if (max_presolve_iter < 5) {
736 max_presolve_iter = 5;
737 }
738
739 if(should_use_hot_starts){
740 set_itlim_hotstart(lp_data, max_presolve_iter);
741 }
742 set_itlim(lp_data, max_presolve_iter);
743 }
744
745 char best_is_est = FALSE;
746 char better_cand_found = FALSE;
747 double prog_ratio = fabs(oldobjval)*0.0001;
748
749 //printf("%i %i %i %i %i\n", max_solves, max_solves_since_impr,
750 //rel_threshold,
751 //strong_br_min_level, max_presolve_iter);
752 //printf("first cand: %i \n", p->br_rel_cand_list[0]);
753
754 int str_br_iter_limit = FALSE;
755 double str_br_factor = MAX(10.0, 3.2e6/(1.0*lp_data->m));
756 int str_br_cnt_limit = (int)(lp_data->n*str_br_factor);
757
758 if (1.0*p->lp_stat.str_br_total_iter_num > str_br_cnt_limit) str_br_iter_limit = TRUE;
759
760 #ifdef COMPILE_IN_LP
761 int node_factor = (int)(p->tm->stat.analyzed/50.0);
762 #else
763 int node_factor = 0;
764 #endif
765 double int_factor = 0.5;
766 if (p->mip->mip_inf){
767 int int_var_num = p->mip->n - p->mip->mip_inf->binary_var_num - p->mip->mip_inf->cont_var_num;
768 if (int_var_num < 1 && p->mip->mip_inf->binary_var_num < 500) {
769 int_factor = 0.1;
770 }
771 }
772
773 double init_ratio = MIN(int_factor*((int)((1.0*lp_data->nz)/1e4) + 1), 2.0);
774
775 if (p->bc_index > 0 && p->lp_stat.str_br_lp_calls > 0) {
776 double str_lp_factor = MAX(0.1, init_ratio - node_factor*0.1);
777 int str_iter = p->lp_stat.str_br_total_iter_num;
778 int lp_iter = p->lp_stat.lp_total_iter_num + str_iter;
779 //printf("str_ratio: %.2f\n", (1.0*str_iter)/lp_iter);
780 //printf("str_iter: %i - lp_iter: %i - node_factor: %i - str_ratio: %.2f\n",
781 // str_iter, lp_iter, node_factor, ((1.0*str_iter)/lp_iter));
782 if (((1.0*str_iter)/lp_iter) > str_lp_factor) {
783 str_br_iter_limit = TRUE;
784 }
785 }
786
787 if (
788 #ifdef COMPILE_IN_LP
789 p->tm->stat.analyzed > 5e5 ||
790 #endif
791 p->lp_stat.str_br_total_iter_num > 5e5) {
792 str_br_iter_limit = TRUE;
793 }
794
795 double frac_avg = 0.0;
796 double frac_tol = 1e-5;
797 for (i=0; i<cand_num; i++) {
798 xval = x[p->br_rel_cand_list[i]];
799 frac_avg += MIN(xval - floor(xval), ceil(xval) - xval);
800 }
801 frac_avg = frac_avg/cand_num;
802 if (frac_avg < 1e-2) {
803 frac_tol = 1e-2;
804 }
805 //printf("frac_avg - %f \n", frac_avg);
806
807 for (i=0; i<cand_num; i++) {
808 //printf("cand - %i \n", i);
809 branch_var = p->br_rel_cand_list[i];
810 lb = vars[branch_var]->new_lb;
811 ub = vars[branch_var]->new_ub;
812 xval = x[branch_var];
813 floorx = floor(xval);
814 ceilx = ceil(xval);
815 rel_down = br_rel_down[branch_var];
816 rel_up = br_rel_up[branch_var];
817
818 // ignore the small violations
819 if (best_can != NULL){
820 if (xval - floorx < frac_tol ||
821 ceilx - xval < frac_tol){
822 //printf("xval: %f\n", xval);
823 continue;
824 }
825 }
826
827 if (cand_num < 2 || str_br_iter_limit ||
828 ((rel_down > rel_threshold &&
829 bc_level > strong_br_min_level) &&
830 (i > check_level || (i < check_level + 1 && !check_first ))) ||
831 (p->par.disable_obj) || (stop_solving == TRUE && rel_down > 1)){
832 down_obj = oldobjval + pcost_down[branch_var] * (xval - floorx);
833 down_is_est = TRUE;
834 p->lp_stat.rel_br_pc_down_num++;
835 } else {
836 if (stop_solving == TRUE){
837 continue;
838 }
839 //down_obj = oldobjval;
840 if (strong_branch(p, branch_var, lb, ub, lb, floorx, &down_obj,
841 should_use_hot_starts, &down_status, &num_down_iters, 0, NULL)) {
842 // lp was abandoned
843 continue;
844 }
845 down_is_est = FALSE;
846 if(p->bc_level < p->br_rel_down_min_level[branch_var]){
847 p->br_rel_down_min_level[branch_var] = p->bc_level;
848 }
849 if (down_status == LP_D_INFEASIBLE || down_status == LP_D_OBJLIM ||
850 down_status == LP_D_UNBOUNDED ||
851 (p->has_ub && down_obj > p->ub - p->par.granularity + p->lp_data->lpetol)) {
852 // update bounds
853 bnd_val[num_bnd_changes] = ceilx;
854 bnd_sense[num_bnd_changes] = 'G';
855 bnd_ind[num_bnd_changes] = branch_var;
856 if (p->mip->colname){
857 PRINT(p->par.verbosity, 5,
858 ("Fixing variable index %i (%s) to 1 \n", branch_var,
859 p->mip->colname[p->lp_data->vars[branch_var]->userind]));
860 }
861 num_bnd_changes++;
862 change_lbub(lp_data, branch_var, ceilx, ub);
863 vars[branch_var]->new_lb = ceilx;
864 vars[branch_var]->lb = ceilx;
865 lb = ceilx;
866 p->br_inf_down[branch_var]++;
867 } else {
868 // update pcost
869 pcost_down[branch_var] = (pcost_down[branch_var]*
870 rel_down + (down_obj - oldobjval)/(xval-floorx))/
871 (rel_down + 1);
872 br_rel_down[branch_var]++;
873 p->lp_stat.rel_br_down_update++;
874 }
875 full_solves++;
876 solves_since_impr++;
877 p->lp_stat.rel_br_full_solve_num++;
878 }
879
880 if (cand_num < 2 || str_br_iter_limit ||
881 ((rel_up > rel_threshold &&
882 bc_level > strong_br_min_level) &&
883 (i > check_level || (i < check_level + 1 && !check_first ))) ||
884 (p->par.disable_obj) || (stop_solving == TRUE && rel_up > 1)){
885 up_obj = oldobjval + pcost_up[branch_var] * (ceilx - xval);
886 up_is_est = TRUE;
887 p->lp_stat.rel_br_pc_up_num++;
888 } else {
889 if (stop_solving == TRUE){
890 continue;
891 }
892 //up_obj = oldobjval;
893 if (strong_branch(p, branch_var, lb, ub, ceilx, ub, &up_obj,
894 should_use_hot_starts, &up_status, &num_up_iters, 0, NULL)) {
895 // lp was abandoned
896 continue;
897 }
898 if(p->bc_level < p->br_rel_up_min_level[branch_var]){
899 p->br_rel_up_min_level[branch_var] = p->bc_level;
900 }
901 up_is_est = FALSE;
902 if (up_status == LP_D_INFEASIBLE || up_status == LP_D_OBJLIM ||
903 up_status == LP_D_UNBOUNDED ||
904 (p->has_ub && up_obj > p->ub - p->par.granularity + p->lp_data->lpetol)) {
905 // update bounds
906 bnd_val[num_bnd_changes] = floorx;
907 bnd_sense[num_bnd_changes] = 'L';
908 bnd_ind[num_bnd_changes] = branch_var;
909 if (p->mip->colname){
910 PRINT(p->par.verbosity, 5,
911 ("Fixing variable index %i (%s) to 1 \n", branch_var,
912 p->mip->colname[p->lp_data->vars[branch_var]->userind]));
913 }
914 num_bnd_changes++;
915 change_lbub(lp_data, branch_var, lb, floorx);
916 vars[branch_var]->new_ub = floorx;
917 vars[branch_var]->ub = floorx;
918 ub = floorx;
919 p->br_inf_up[branch_var]++;
920 } else {
921 // update pcost
922 pcost_up[branch_var] = (pcost_up[branch_var]*
923 rel_up + (up_obj - oldobjval)/(ceilx-xval))/
924 (rel_up + 1);
925 br_rel_up[branch_var]++;
926 p->lp_stat.rel_br_up_update++;
927 }
928 full_solves++;
929 solves_since_impr++;
930 p->lp_stat.rel_br_full_solve_num++;
931 }
932
933 if (down_obj > SYM_INFINITY/10 && up_obj > SYM_INFINITY/10) {
934 //printf("d u %f %f\n", down_obj, up_obj);
935 both_children_inf = TRUE;
936 best_can = candidates[0];
937 break;
938 }
939
940 if ((down_obj > SYM_INFINITY/10 || up_obj > SYM_INFINITY/10)) {
941 var_score = MIN(down_obj, up_obj);
942 if(best_can != NULL) {
943 continue;
944 }else{
945 best_one_child_inf = TRUE;
946 }
947 } else {
948 if (down_obj < up_obj) {
949 low = down_obj;
950 high = up_obj;
951 } else {
952 low = up_obj;
953 high = down_obj;
954 }
955 var_score = alpha * low + one_m_alpha * high;
956 }
957
958 double violation_cnt_diff = 0;
959 int inf_cnt_diff = 0;
960 int sos_diff = 0;
961 int frac_cnt_diff = 0, nz_diff = 0;
962
963 if(best_can){
964 inf_cnt_diff = MAX(p->br_inf_up[branch_var], p->br_inf_down[branch_var]) -
965 MAX(p->br_inf_up[best_var], p->br_inf_down[best_var]);
966 }
967
968 if(best_can){
969 if(down_violation_cnt){
970 double cand_v = 0.0, best_v = 0.0;
971 if(violation_col_size[branch_var]){
972 cand_v = 1.0*MAX(down_violation_cnt[branch_var],
973 up_violation_cnt[branch_var])/violation_col_size[branch_var];
974 }
975 if(violation_col_size[best_var]){
976 best_v = 1.0*MAX(down_violation_cnt[best_var],
977 up_violation_cnt[best_var])/violation_col_size[best_var];
978 }
979 violation_cnt_diff = cand_v - best_v;
980 }
981
982 if(p->mip->mip_inf){
983 sos_diff = p->mip->mip_inf->cols[branch_var].sos_num -
984 p->mip->mip_inf->cols[best_var].sos_num;
985
986 //frac_cnt_diff = frac_cnt[branch_var] - frac_cnt[best_var];
987 frac_cnt_diff = (int)(p->var_rank[branch_var] -
988 p->var_rank[best_var]);
989 nz_diff = p->mip->mip_inf->cols[branch_var].nz -
990 p->mip->mip_inf->cols[best_var].nz;
991 }
992 }
993
994 int tot_var_score = 0;
995 if(best_can){
996 better_cand_found = FALSE;
997
998 int s_score = 1, v_score = 32, i_score = 16, f_score = 8, b_score = 4, z_score = 2;
999
1000 double branch_var_frac = fabs(0.5 -(x[best_var] - floorx)) - fabs(0.5 -(x[branch_var] - floorx));
1001 tot_var_score = (sos_diff > 0 ? s_score: (sos_diff < 0 ? -s_score:0)) +
1002 (violation_cnt_diff > 0.0 ? v_score: (violation_cnt_diff < 0.0 ? -v_score:0)) +
1003 (inf_cnt_diff > 0 ? i_score: (inf_cnt_diff < 0 ? - i_score:0)) +
1004 (frac_cnt_diff > 0 ? f_score: (frac_cnt_diff < 0 ? -f_score:0)) +
1005 (branch_var_frac > 0.0 ? b_score: (branch_var_frac < 0.0 ? -b_score:0)) +
1006 (nz_diff > 0 ? z_score : (nz_diff < 0 ? -z_score:0));
1007
1008 //printf("s : v : i : f : b : n : %i %i %i %i %f %i", sos_diff,
1009 // violation_cnt_diff, inf_cnt_diff, frac_cnt_diff, branch_var_frac, nz_diff);
1010
1011 int c_score = 0;
1012 if(!p->par.disable_obj){
1013 char cand_is_est = ((down_is_est && up_is_est) ? TRUE : FALSE);
1014 double score_diff = var_score - best_var_score;
1015 if(score_diff > lpetol100) c_score = 100;
1016 else if(score_diff < -lpetol100) c_score = -100;
1017
1018 if(cand_is_est || best_is_est){
1019 c_score = 0;
1020 if(cand_is_est && best_is_est) {
1021 if(score_diff > lpetol100) c_score = 100;
1022 else if(score_diff < -lpetol100) c_score = -100;
1023 }else{
1024 if(best_is_est){
1025 if(score_diff > lpetol100) c_score = 100;
1026 else if(score_diff < -lpetol100) c_score = -32;
1027 }else{
1028 if(score_diff > lpetol100) c_score = 32;
1029 else if(score_diff < -lpetol100) c_score = -100;
1030 }
1031 }
1032 }
1033
1034 tot_var_score += c_score;
1035 }
1036
1037 if(tot_var_score > 0){
1038 better_cand_found = TRUE;
1039 }
1040 }
1041
1042 if(best_can == NULL || better_cand_found || best_one_child_inf){
1043 //printf("here - %i\n", p->bc_index);
1044 if ( var_score > best_var_score + prog_ratio &&(down_is_est != TRUE ||
1045 up_is_est != TRUE)) {
1046 solves_since_impr = 0;
1047 if(best_can!= NULL){
1048 p->lp_stat.rel_br_impr_num++;
1049 }
1050 }
1051
1052 if(best_can != NULL && best_one_child_inf) {
1053 best_one_child_inf = FALSE;
1054 }
1055
1056 if(down_is_est && up_is_est) best_is_est = TRUE;
1057 else best_is_est = FALSE;
1058
1059 //printf("%f %f\n", var_score, best_var_score);
1060 best_var_score = var_score;
1061 best_var = branch_var;
1062 best_can = candidates[0];
1063 best_can->position = branch_var;
1064 best_can->solutions = NULL;
1065 best_can->sol_inds = NULL;
1066 best_can->sol_sizes = NULL;
1067 best_can->sense[1] = 'L';
1068 best_can->sense[0] = 'G';
1069 if (down_is_est==TRUE) {
1070 best_can->objval[1] = oldobjval;
1071 best_can->iterd[1] = 0;
1072 best_can->termcode[1] = LP_D_ITLIM;
1073 } else {
1074 best_can->objval[1] = down_obj;
1075 best_can->iterd[1] = num_down_iters;
1076 best_can->termcode[1] = down_status;
1077 // added by asm4 because hot starts dont generate a reliable
1078 // bound.
1079 //if (should_use_hot_starts && down_status==LP_D_ITLIM) {
1080 // down_is_est = TRUE;
1081 // best_can->objval[0] = oldobjval;
1082 //}
1083 }
1084 if (up_is_est==TRUE) {
1085 best_can->objval[0] = oldobjval;
1086 best_can->iterd[0] = 0;
1087 best_can->termcode[0] = LP_D_ITLIM;
1088 } else {
1089 best_can->objval[0] = up_obj;
1090 best_can->iterd[0] = num_up_iters;
1091 best_can->termcode[0] = up_status;
1092 // added by asm4 because hot starts dont generate a reliable
1093 // bound.
1094 //if (should_use_hot_starts && up_status==LP_D_ITLIM) {
1095 // up_is_est = TRUE;
1096 // best_can->objval[1] = oldobjval;
1097 //}
1098 }
1099 best_can->is_est[1] = down_is_est;
1100 best_can->is_est[0] = up_is_est;
1101 best_can->rhs[1] = floorx;
1102 best_can->rhs[0] = ceilx;
1103 best_can->value = xval;
1104 best_down_is_est = down_is_est;
1105 best_up_is_est = up_is_est;
1106
1107 if(best_can->objval[0] < best_can->objval[1] + lpetol100 &&
1108 best_can->objval[0] > best_can->objval[1] - lpetol100){
1109 char swap = TRUE;
1110 double objcoef;
1111 get_objcoef(lp_data, branch_var, &objcoef);
1112
1113 if(objcoef > -lpetol) swap = FALSE;
1114 else if(objcoef > -lpetol){
1115 double var_frac_diff = fabs(0.5 -(x[best_var] - floorx)) - fabs(0.5 -(x[branch_var] - floorx));
1116 int var_inf_cnt_diff = p->br_inf_up[branch_var] - p->br_inf_down[branch_var];
1117 int var_violation_cnt_diff = 0;
1118 if(up_violation_cnt){
1119 var_violation_cnt_diff = up_violation_cnt[branch_var] - down_violation_cnt[branch_var];
1120 }
1121 int v_score = 4, i_score = 2, b_score = 1;
1122
1123 int tot_var_score =
1124 (var_violation_cnt_diff > 0.0 ? v_score: (violation_cnt_diff < 0.0 ? -v_score:0)) +
1125 (var_inf_cnt_diff > 0 ? i_score: (inf_cnt_diff < 0 ? - i_score:0)) +
1126 (var_frac_diff > 0.0 ? b_score: (var_frac_diff < 0.0 ? -b_score:0));
1127 if(tot_var_score > 0) swap = FALSE;
1128 }
1129
1130 if(swap){
1131 best_can->sense[0] = 'L';
1132 best_can->sense[1] = 'G';
1133 if (down_is_est==TRUE) {
1134 best_can->objval[0] = oldobjval;
1135 best_can->iterd[0] = 0;
1136 best_can->termcode[0] = LP_D_ITLIM;
1137 } else {
1138 best_can->objval[0] = down_obj;
1139 best_can->iterd[0] = num_down_iters;
1140 best_can->termcode[0] = down_status;
1141 }
1142 if (up_is_est==TRUE) {
1143 best_can->objval[1] = oldobjval;
1144 best_can->iterd[1] = 0;
1145 best_can->termcode[1] = LP_D_ITLIM;
1146 } else {
1147 best_can->objval[1] = up_obj;
1148 best_can->iterd[1] = num_up_iters;
1149 best_can->termcode[1] = up_status;
1150 }
1151 best_can->is_est[0] = down_is_est;
1152 best_can->is_est[1] = up_is_est;
1153 best_can->rhs[0] = floorx;
1154 best_can->rhs[1] = ceilx;
1155 }
1156 }
1157 }
1158
1159 //printf("solves_no_imp %i\n", solves_since_impr);
1160 if ((solves_since_impr > max_solves_since_impr ||
1161 full_solves >= max_solves) || p->par.rs_mode_enabled) {
1162 //printf("breaking because of no gain at iter %d\n", i);
1163 //printf("%i %i %i %i\n", p->bc_level, cand_num, solves_since_impr, full_solves);
1164 stop_solving = TRUE;
1165 }
1166 }
1167 //printf("reliability branching: selected var %d with score %f\n", best_var, best_var_score);
1168
1169 #ifdef COMPILE_IN_LP
1170 if (num_bnd_changes > 0) {
1171 str_br_bound_changes(p, num_bnd_changes, bnd_val, bnd_ind, bnd_sense);
1172 }
1173 #endif
1174
1175
1176 // experimental - sos branching - not tested
1177 if(p->par.use_sos_branching && !both_children_inf && p->mip->mip_inf &&
1178 1.0*p->mip->mip_inf->binary_var_num/(p->mip->n + 1) > 0.5 &&
1179 p->bc_level <= p->par.sos_branching_max_level && p->mip->mip_inf->binary_sos_row_num){
1180
1181 //printf("\nsos row cnt %i", p->mip->mip_inf->binary_sos_row_num);
1182 if (should_use_hot_starts) {
1183 unmark_hs = FALSE;
1184 unmark_hotstart(lp_data);
1185 set_itlim_hotstart(lp_data, -1);
1186 }
1187 double sos_best_var_score = -SYM_INFINITY;
1188 int sos_best_f_cnt = 0;
1189
1190 if (max_presolve_iter < 5) {
1191 max_presolve_iter = 5;
1192 }
1193 set_itlim(lp_data, max_presolve_iter);
1194
1195 //p->mip->mip_inf->cols[best_var].sos_num > 0){ //}
1196 int *l_ind = NULL, *r_ind = NULL;
1197
1198 int col_num = lp_data->n;
1199 int row_num = p->mip->m; // p->base.cutnum ?
1200 int col_ind, row_ind;
1201 int maxmn = MAX(row_num, col_num);
1202 int row_size,row_frac_cnt;
1203
1204 //int *max_frac_ind = lp_data->tmp.i1;
1205 //int *frac_ind = lp_data->tmp.i1 + col_num;
1206 char *col_stat = lp_data->tmp.c + 2*maxmn;
1207
1208 //double *max_frac_val = lp_data->tmp.d + col_num;
1209 //double *frac_val = lp_data->tmp.d + 2*col_num;
1210 int *row_z_cnt = lp_data->tmp.i1;
1211 int *sos_row_size = lp_data->tmp.i1+maxmn;
1212 int *sos_row = NULL;
1213 int sos_row_cnt = 0;
1214 int *row_frac_freq = lp_data->tmp.i1+2*maxmn;
1215
1216 ROWinfo *rows = p->mip->mip_inf->rows;
1217 //COLinfo *cols = p->mip->mip_inf->cols;
1218
1219 int *row_matbeg = p->mip->row_matbeg;
1220 int *row_matind = p->mip->row_matind;
1221 //double *row_matval = p->mip->row_matval;
1222
1223 int *matbeg = p->mip->matbeg;
1224 int *matind = p->mip->matind;
1225 //double *matval = p->mip->matval;
1226 double ub, lb;
1227 memset(row_frac_freq, 0, ISIZE*row_num);
1228 memset(sos_row_size, 0, ISIZE*row_num);
1229 memset(row_z_cnt, 0, ISIZE*row_num);
1230 for(i = 0; i < col_num; i++){
1231 //col_stat[i] = 'N'; // not required
1232 if(vars[i]->is_int){
1233 col_stat[i] = 'I'; //integer
1234 get_ub(lp_data, i, &ub);
1235 get_lb(lp_data, i, &lb);
1236 int col_size = matbeg[i+1] - matbeg[i];
1237 if(ub > lb + lpetol){
1238 col_stat[i] = 'U'; // integer but unfixed yet
1239 if(x[i] - floor(x[i]) > lpetol && ceil(x[i]) - x[i] > lpetol){
1240 col_stat[i] = 'F'; // fractional
1241 for(j = matbeg[i]; j < matbeg[i + 1]; j++){
1242 row_frac_freq[matind[j]]++;
1243 row_z_cnt[matind[j]] -= col_size;
1244 }
1245 }else{
1246 for(j = matbeg[i]; j < matbeg[i + 1]; j++){
1247 row_z_cnt[matind[j]] -= col_size;
1248 }
1249 }
1250 }else{
1251 for(j = matbeg[i]; j < matbeg[i + 1]; j++){
1252 sos_row_size[matind[j]]++;
1253 }
1254 }
1255 }else{
1256 col_stat[i] = 'C'; //continuous
1257 }
1258 }
1259
1260 for(row_ind = 0; row_ind < row_num; row_ind++){
1261 sos_row_size[row_ind] = rows[row_ind].size - sos_row_size[row_ind];
1262 if(rows[row_ind].is_sos_row && row_frac_freq[row_ind] > 1 && sos_row_size[row_ind] > 4){
1263 if(!sos_row){
1264 sos_row = (int*)malloc(ISIZE*row_num);
1265 }
1266 sos_row_size[sos_row_cnt] = -sos_row_size[row_ind] - row_frac_freq[row_ind];
1267 sos_row[sos_row_cnt++] = row_ind;
1268 }
1269 }
1270 //printf("...cnt %i\n", sos_row_cnt);
1271 if(sos_row_cnt > 0){
1272 //qsort_ii(sos_row_size, sos_row, sos_row_cnt);
1273 qsort_ii(row_z_cnt, sos_row, sos_row_cnt);
1274 }
1275 int final_cnt = MIN(5, sos_row_cnt);
1276 double *frac_val = lp_data->tmp.d + maxmn;
1277 int *frac_ind = lp_data->tmp.i1;
1278 char *l_assigned = lp_data->tmp.c;
1279 char *r_assigned = lp_data->tmp.c + maxmn;
1280
1281 for(int k = 0; k < final_cnt; k++){
1282 row_ind = sos_row[k];
1283 int total_f_cnt = 0;
1284 if(rows[row_ind].is_sos_row){
1285 row_frac_cnt = 0;
1286 row_size = 0;
1287 for(i = row_matbeg[row_ind]; i < row_matbeg[row_ind+1]; i++){
1288 col_ind = row_matind[i];
1289 if(col_stat[col_ind] == 'C') {
1290 printf("ERROR in sos branching... - row %i col %i\n",
1291 row_ind, col_ind);
1292 continue;
1293 //exit(0);
1294 }
1295 if(col_stat[col_ind] != 'I') row_size++;
1296 if(col_stat[col_ind] == 'F'){
1297 total_f_cnt += matbeg[col_ind + 1] - matbeg[col_ind];
1298 frac_ind[row_frac_cnt] = col_ind;
1299 frac_val[row_frac_cnt] = floor(x[col_ind]) - x[col_ind];
1300 row_frac_cnt++;
1301 }
1302 //if(col_stat[i] == 'U') bin_cnt[sos_num]--;
1303 //else frac_cnt[sos_num]--;
1304 }
1305
1306 int l_cnt = 0, r_cnt = 0;
1307 if(!l_ind)
1308 l_ind = (int*)malloc(ISIZE*col_num);
1309 if(!r_ind)
1310 r_ind = (int*)malloc(ISIZE*col_num);
1311 qsort_di(frac_val, frac_ind, row_frac_cnt);
1312 l_cnt = row_frac_cnt/2;
1313 r_cnt = row_frac_cnt - l_cnt;
1314 l_ind[0] = frac_ind[0];
1315 r_ind[0] = frac_ind[1];
1316 if(r_cnt > 1){
1317 memcpy(r_ind + 1, frac_ind + 2, ISIZE*(r_cnt - 1));
1318 }
1319 if(l_cnt > 1){
1320 memcpy(l_ind + 1, frac_ind + r_cnt + 1, ISIZE*(l_cnt - 1));
1321 }
1322
1323 int l_assigned_cnt = 0, r_assigned_cnt = 0;
1324
1325 if(row_size > row_frac_cnt){
1326 memset(l_assigned, 0, CSIZE*row_num);
1327 memset(r_assigned, 0, CSIZE*row_num);
1328
1329 for(i = 0; i < l_cnt; i++){
1330 col_ind = l_ind[i];
1331 for(j = matbeg[col_ind]; j < matbeg[col_ind + 1]; j++){
1332 if(!l_assigned[matind[j]]){
1333 l_assigned[matind[j]] = TRUE;
1334 l_assigned_cnt++;
1335 }
1336 }
1337 }
1338 for(i = 0; i < r_cnt; i++){
1339 col_ind = r_ind[i];
1340 for(j = matbeg[col_ind]; j < matbeg[col_ind + 1]; j++){
1341 if(!r_assigned[matind[j]]){
1342 r_assigned[matind[j]] = TRUE;
1343 r_assigned_cnt++;
1344 }
1345 }
1346 }
1347
1348 int bin_l_cnt = 0, bin_r_cnt = 0;
1349 for(i = row_matbeg[row_ind];
1350 i < row_matbeg[row_ind + 1]; i++){
1351 col_ind = row_matind[i];
1352 if(col_stat[col_ind] == 'U'){
1353 bin_l_cnt = bin_r_cnt = 0;
1354 for(j = matbeg[col_ind]; j < matbeg[col_ind + 1]; j++){
1355 if(l_assigned[matind[j]]) bin_l_cnt++;
1356 if(r_assigned[matind[j]]) bin_r_cnt++;
1357 }
1358
1359 if(bin_l_cnt > bin_r_cnt) {
1360 l_ind[l_cnt++] = col_ind;
1361 }else if(bin_l_cnt < bin_r_cnt){
1362 r_ind[r_cnt++] = col_ind;
1363 }else{
1364 if(l_cnt < r_cnt) l_ind[l_cnt++] = col_ind;
1365 else r_ind[r_cnt++] = col_ind;
1366 }
1367 }
1368 }
1369
1370 }
1371
1372 strong_branch(p, 0, 0.0, 0.0, 0.0, 0.0, &down_obj,
1373 FALSE,
1374 &down_status, &num_down_iters, l_cnt,
1375 l_ind);
1376 strong_branch(p, 0, 0.0, 0.0, 0.0, 0.0, &up_obj,
1377 FALSE,
1378 &up_status, &num_up_iters, r_cnt,
1379 r_ind);
1380
1381 if (down_obj > SYM_INFINITY/10 && up_obj > SYM_INFINITY/10) {
1382 both_children_inf = TRUE;
1383 FREE(best_can->sos_ind[0]);
1384 FREE(best_can->sos_ind[1]);
1385 best_can = candidates[0];
1386 break;
1387 }
1388
1389 if (down_obj < up_obj) {
1390 low = down_obj;
1391 high = up_obj;
1392 } else {
1393 low = up_obj;
1394 high = down_obj;
1395 }
1396
1397 double sos_score =
1398 alpha * low + one_m_alpha * high;
1399 int can_iterate = FALSE;
1400 if((best_down_is_est && best_up_is_est) ||
1401 sos_score > SYM_INFINITY/10 || best_var_score < SYM_INFINITY/10) {
1402 can_iterate = TRUE;
1403 //printf("...can iterate \n");
1404 }
1405 if(can_iterate && sos_score > best_var_score - lpetol100 &&
1406 (sos_score > sos_best_var_score + lpetol100 || !(best_can->sos_ind[0]) ||
1407 (sos_score > sos_best_var_score - lpetol100 &&
1408 total_f_cnt > sos_best_f_cnt))){
1409
1410 sos_best_var_score = sos_score;
1411 sos_best_f_cnt = total_f_cnt;
1412
1413 int li = 0;
1414 int ri = 1;
1415
1416 if(l_assigned_cnt < r_assigned_cnt) {
1417 li = 1;
1418 ri = 0;
1419 }
1420
1421 if(!best_can->sos_ind[li]){
1422 best_can->sos_ind[li] = (int*)malloc(ISIZE*l_cnt);
1423 }else if(l_cnt > best_can->sos_cnt[li]){
1424 best_can->sos_ind[li] = (int*)realloc((char *)best_can->sos_ind[li],
1425 ISIZE*l_cnt);
1426 }
1427 if(!best_can->sos_ind[ri]){
1428 //printf("...accepted\n");
1429 best_can->sos_ind[ri] = (int*)malloc(ISIZE*r_cnt);
1430 }else if (r_cnt > best_can->sos_cnt[ri]){
1431 best_can->sos_ind[ri] = (int*)realloc((char*)best_can->sos_ind[ri],
1432 ISIZE*r_cnt);
1433 }
1434
1435 memcpy(best_can->sos_ind[li], l_ind, ISIZE*l_cnt);
1436 memcpy(best_can->sos_ind[ri], r_ind, ISIZE*r_cnt);
1437
1438 best_can->type = SOS1_IMPLICIT;
1439 best_can->sense[li] = 'L';
1440 best_can->sense[ri] = 'G';
1441
1442 best_can->objval[li] = down_obj;
1443 best_can->iterd[li] = num_down_iters;
1444 best_can->termcode[li] = down_status;
1445
1446 best_can->objval[ri] = up_obj;
1447 best_can->iterd[ri] = num_up_iters;
1448 best_can->termcode[ri] = up_status;
1449
1450 best_can->rhs[li] = 0.0;
1451 best_can->rhs[ri] = 0.0;
1452
1453 best_can->sos_cnt[li] = l_cnt;
1454 best_can->sos_cnt[ri] = r_cnt;
1455 }
1456 }
1457 }
1458 FREE(l_ind);
1459 FREE(r_ind);
1460 FREE(sos_row);
1461 }
1462
1463 cand_num = 1;
1464 FREE(up_violation_cnt);
1465 FREE(down_violation_cnt);
1466 FREE(violation_col_size);
1467
1468 if (both_children_inf || num_bnd_changes > 0) {
1469 FREE(best_can);
1470 FREE(candidates);
1471 *candidate = NULL;
1472 if (should_use_hot_starts && unmark_hs){
1473 unmark_hotstart(lp_data);
1474 set_itlim_hotstart(lp_data, -1);
1475 }else{
1476 load_basis(lp_data, rstat, cstat);
1477 }
1478 set_itlim(lp_data, -1); //both limits should be set for hotstarts
1479 if (both_children_inf){
1480 p->lp_stat.str_br_nodes_pruned++;
1481 return (DO_NOT_BRANCH__FATHOMED);
1482 }else{
1483 return (DO_NOT_BRANCH);
1484 }
1485 }
1486
1487 //printf("Branching on %i %c\n", best_can->position, best_can->sense[0]);
1488
1489 } else {
1490 /* do the default symphony branching */
1491
1492 /*
1493 * see if hot-starts should be used. in theory if strong branching is used
1494 * and only variable bounds are changed then, hot-starts should be faster
1495 */
1496
1497 #ifdef __OSI_CLP__
1498 lp_data->si->setupForRepeatedUse(2,0);
1499 #endif
1500 if (p->par.use_hot_starts && !p->par.branch_on_cuts) {
1501 should_use_hot_starts = TRUE;
1502 } else {
1503 should_use_hot_starts = FALSE;
1504 }
1505
1506 /* Set the iteration limit */
1507 if (should_use_hot_starts) {
1508 mark_hotstart(lp_data);
1509 }
1510
1511 if (p->par.max_presolve_iter > 0) {
1512 max_presolve_iter = p->par.max_presolve_iter - bc_level;
1513
1514 if (max_presolve_iter < 5) {
1515 max_presolve_iter = 5;
1516 }
1517 if(should_use_hot_starts){
1518 set_itlim_hotstart(lp_data, max_presolve_iter);
1519 }
1520 set_itlim(lp_data, max_presolve_iter);
1521 }
1522
1523 for (i=0; i<cand_num; i++){
1524 can = candidates[i];
1525
1526 #ifndef MAX_CHILDREN_NUM
1527 can->objval = pobj;
1528 can->termcode = pterm;
1529 can->feasible = pfeas;
1530 can->iterd = piter;
1531 if (p->tm->par.keep_description_of_pruned == KEEP_IN_MEMORY){
1532 can->solutions = (double **) calloc(maxnum, sizeof(double *));
1533 can->sol_inds = (int **) calloc(maxnum, sizeof(int *));
1534 can->sol_size = (int *) calloc(maxnum, ISIZE);
1535 }
1536
1537 #ifdef SENSITIVITY_ANALYSIS
1538 if (p->tm->par.sensitivity_analysis){
1539 can->duals = (double **) calloc(maxnum, sizeof(double *));
1540 }else{
1541 can->duals = NULL;
1542 }
1543 #endif
1544 #ifdef COMPILE_FRAC_BRANCHING
1545 can->frac_num = pfrnum;
1546 can->frac_ind = pfrind;
1547 can->frac_val = pfrval;
1548 #endif
1549
1550 #else
1551 if (p->par.keep_description_of_pruned == KEEP_IN_MEMORY){
1552 can->solutions = (double **) calloc (MAX_CHILDREN_NUM,
1553 sizeof(double *));
1554 can->sol_inds = (int **) calloc(MAX_CHILDREN_NUM,
1555 sizeof(int *));
1556 can->sol_sizes = (int *) calloc(MAX_CHILDREN_NUM, ISIZE);
1557 }
1558 #ifdef SENSITIVITY_ANALYSIS
1559 if (p->par.sensitivity_analysis){
1560 can->duals = (double **) calloc (MAX_CHILDREN_NUM, sizeof(double *));
1561 }else{
1562 can->duals = NULL;
1563 }
1564 #endif
1565 #endif
1566
1567 #ifdef STATISTICS
1568 cnum += can->child_num;
1569 #endif
1570
1571 /* Now depending on the type, adjust ub/lb or rhs/range/sense */
1572 switch (can->type){
1573 case CANDIDATE_VARIABLE:
1574 branch_var = can->position;
1575 #if 0
1576 if (lp_data->status[branch_var] & PERM_FIXED_TO_LB ||
1577 lp_data->status[branch_var] & PERM_FIXED_TO_UB){
1578 if (vars[branch_var]->lb == vars[branch_var]->ub){
1579 printf("Warning -- branching candidate is already fixed. \n");
1580 printf("SYMPHONY has encountered numerical difficulties \n");
1581 printf("With the LP solver. Exiting...\n\n");
1582 }
1583 /* } to unconfuse vi*/
1584 #endif
1585 lb = vars[branch_var]->new_lb;
1586 ub = vars[branch_var]->new_ub;
1587 for (j = 0; j < can->child_num; j++){
1588 switch (can->sense[j]){
1589 case 'E':
1590 change_lbub(lp_data, branch_var, can->rhs[j], can->rhs[j]);
1591 break;
1592 case 'R':
1593 change_lbub(lp_data, branch_var, can->rhs[j],
1594 can->rhs[j] + can->range[j]);
1595 break;
1596 case 'L':
1597 change_lbub(lp_data, branch_var, lb, can->rhs[j]);
1598 break;
1599 case 'G':
1600 change_lbub(lp_data, branch_var, can->rhs[j], ub);
1601 break;
1602 }
1603 check_ub(p);
1604 /* The original basis is in lp_data->lpbas */
1605 if (should_use_hot_starts) {
1606 can->termcode[j] = solve_hotstart(lp_data, can->iterd+j);
1607 total_iters+=*(can->iterd+j);
1608 } else {
1609 load_basis(lp_data, cstat, rstat);
1610 can->termcode[j] = dual_simplex(lp_data, can->iterd+j);
1611 total_iters+=*(can->iterd+j);
1612 }
1613 p->lp_stat.lp_calls++;
1614 p->lp_stat.str_br_lp_calls++;
1615 p->lp_stat.str_br_total_iter_num += *(can->iterd+j);
1616 can->objval[j] = lp_data->objval;
1617 //get_x(lp_data);
1618
1619 #ifdef SENSITIVITY_ANALYSIS
1620 if (p->par.sensitivity_analysis){
1621 get_dj_pi(lp_data);
1622 can->duals[j] = (double *) malloc (DSIZE*p->base.cutnum);
1623 memcpy(can->duals[j], lp_data->dualsol, DSIZE*p->base.cutnum);
1624 }
1625 #endif
1626
1627 if (can->termcode[j] == LP_OPTIMAL){
1628 /* is_feasible_u() fills up lp_data->x, too!! */
1629 switch (is_feasible_u(p, TRUE, FALSE)){
1630
1631 /*NOTE: This is confusing but not all that citical...*/
1632 /*The "feasible" field is only filled out for the
1633 purposes of display (in vbctool) to keep track of
1634 where in the tree the feasible solutions were
1635 found. Since this may not be the actual candidate
1636 branched on, we need to pass this info on to whatever
1637 candidate does get branched on so the that the fact that
1638 a feasible solution was found in presolve can be recorded*/
1639
1640 case IP_FEASIBLE:
1641 can->termcode[j] = LP_OPT_FEASIBLE;
1642 can->feasible[j] = TRUE;
1643 if (p->par.keep_description_of_pruned == KEEP_IN_MEMORY){
1644 can->solutions[j] = (double *) malloc (DSIZE*lp_data->n);
1645 memcpy(can->solutions[j], lp_data->x, DSIZE*lp_data->n);
1646 }
1647 break;
1648
1649 case IP_FEASIBLE_BUT_CONTINUE:
1650 can->termcode[j] = LP_OPT_FEASIBLE_BUT_CONTINUE;
1651 can->feasible[j] = TRUE;
1652 if (p->par.keep_description_of_pruned == KEEP_IN_MEMORY){
1653 can->solutions[j] = (double *) malloc (DSIZE*lp_data->n);
1654 memcpy(can->solutions[j], lp_data->x, DSIZE*lp_data->n);
1655 }
1656 break;
1657
1658 default:
1659 break;
1660 }
1661 } else if (can->termcode[j] == LP_D_OBJLIM ||
1662 can->termcode[j] == LP_D_UNBOUNDED ||
1663 can->termcode[j] == LP_D_INFEASIBLE){
1664 //p->bound_changes_in_iter++;
1665 switch (can->sense[j]){
1666 case 'L':
1667 /* decreasing the ub made the problem inf, so change lb */
1668 //lb = can->rhs[j] + 1;
1669 //vars[can->position]->new_lb = lb;
1670 break;
1671 case 'G':
1672 //ub = can->rhs[j] - 1;
1673 //vars[can->position]->new_ub = ub;
1674 break;
1675 case 'E':
1676 /* problem becomes infeasible */
1677 /* dont know what to do */
1678 break;
1679 }
1680 }
1681 #ifdef COMPILE_FRAC_BRANCHING
1682 else{
1683 if (can->termcode[j] != LP_ABANDONED){
1684 //get_x(lp_data);
1685 }
1686 }
1687 if (can->termcode[j] != LP_ABANDONED){
1688 xind = lp_data->tmp.i1; /* n */
1689 xval = lp_data->tmp.d; /* n */
1690 can->frac_num[j] = collect_fractions(p, lp_data->x, xind, xval);
1691 if (can->frac_num[j] > 0){
1692 can->frac_ind[j] = (int *) malloc(can->frac_num[j] * ISIZE);
1693 can->frac_val[j] = (double *) malloc(can->frac_num[j]*DSIZE);
1694 memcpy(can->frac_ind[j], xind, can->frac_num[j] * ISIZE);
1695 memcpy(can->frac_val[j], xval, can->frac_num[j] * DSIZE);
1696 }
1697 }else{
1698 can->frac_num[j] = 0;
1699 }
1700 #endif
1701 #ifdef STATISTICS
1702 if (can->termcode[j] == LP_D_ITLIM)
1703 itlim++;
1704 #endif
1705 }
1706 change_lbub(lp_data, branch_var, lb, ub);
1707 break;
1708
1709 case CANDIDATE_CUT_IN_MATRIX:
1710 branch_row = can->position;
1711 for (j = 0; j < can->child_num; j++){
1712 change_row(lp_data, branch_row,
1713 can->sense[j], can->rhs[j], can->range[j]);
1714 check_ub(p);
1715 /* The original basis is in lp_data->lpbas */
1716 can->termcode[j] = dual_simplex(lp_data, can->iterd+j);
1717 p->lp_stat.lp_calls++;
1718 p->lp_stat.str_br_lp_calls++;
1719 p->lp_stat.str_br_total_iter_num += *(can->iterd+j);
1720 can->objval[j] = lp_data->objval;
1721
1722
1723 //get_x(lp_data);
1724
1725 #ifdef SENSITIVITY_ANALYSIS
1726 if (p->par.sensitivity_analysis){
1727 get_dj_pi(lp_data);
1728 can->duals[j] = (double *) malloc (DSIZE*p->base.cutnum);
1729 memcpy(can->duals[j], lp_data->dualsol, DSIZE*p->base.cutnum);
1730 }
1731 #endif
1732
1733 if (can->termcode[j] == LP_OPTIMAL){
1734 /* is_feasible_u() fills up lp_data->x, too!! */
1735 switch (is_feasible_u(p, TRUE, FALSE)){
1736
1737 /*NOTE: This is confusing but not all that citical...*/
1738 /*The "feasible" field is only filled out for the
1739 purposes of display (in vbctool) to keep track of
1740 where in the tree the feasible solutions were
1741 found. Since this may not be the actual candidate
1742 branched on, we need to pass this info on to whatever
1743 candidate does get branched on so the that the fact that
1744 a feasible solution was found in presolve can be recorded*/
1745
1746 case IP_FEASIBLE:
1747 can->termcode[j] = LP_OPT_FEASIBLE;
1748 can->feasible[j] = TRUE;
1749 if (p->par.keep_description_of_pruned == KEEP_IN_MEMORY){
1750 can->solutions[j] = (double *) malloc (DSIZE*
1751 lp_data->n);
1752 memcpy(can->solutions[j], lp_data->x, DSIZE*lp_data->n);
1753 }
1754 break;
1755
1756 case IP_FEASIBLE_BUT_CONTINUE:
1757 can->termcode[j] = LP_OPT_FEASIBLE_BUT_CONTINUE;
1758 can->feasible[j] = TRUE;
1759 if (p->par.keep_description_of_pruned == KEEP_IN_MEMORY){
1760 can->solutions[j] = (double *) malloc (DSIZE*
1761 lp_data->n);
1762 memcpy(can->solutions[j], lp_data->x, DSIZE*lp_data->n);
1763 }
1764 break;
1765
1766 default:
1767 break;
1768 }
1769 }
1770 #ifdef COMPILE_FRAC_BRANCHING
1771 else{
1772 if (can->termcode[j] != LP_ABANDONED)
1773 //get_x(lp_data);
1774 }
1775 if (can->termcode[j] != LP_ABANDONED){
1776 xind = lp_data->tmp.i1; /* n */
1777 xval = lp_data->tmp.d; /* n */
1778 can->frac_num[j] = collect_fractions(p, lp_data->x, xind, xval);
1779 if (can->frac_num[j] > 0){
1780 can->frac_ind[j] = (int *) malloc(can->frac_num[j] * ISIZE);
1781 can->frac_val[j] = (double *) malloc(can->frac_num[j]*DSIZE);
1782 memcpy(can->frac_ind[j], xind, can->frac_num[j] * ISIZE);
1783 memcpy(can->frac_val[j], xval, can->frac_num[j] * DSIZE);
1784 }
1785 }else{
1786 can->frac_num[j] = 0;
1787 }
1788 #endif
1789 #ifdef STATISTICS
1790 if (can->termcode[j] == LP_D_ITLIM)
1791 itlim++;
1792 #endif
1793 }
1794 cut = rows[branch_row].cut;
1795 change_row(lp_data, branch_row, cut->sense, cut->rhs, cut->range);
1796 free_row_set(lp_data, 1, &branch_row);
1797 break;
1798 }
1799
1800 switch ((j = compare_candidates_u(p, oldobjval, best_can, can))){
1801 case FIRST_CANDIDATE_BETTER:
1802 case FIRST_CANDIDATE_BETTER_AND_BRANCH_ON_IT:
1803 if (p->par.keep_description_of_pruned == KEEP_IN_MEMORY){
1804 min_ind = -1;
1805 for (k = can->child_num - 1; k >= 0; k--){
1806 if (can->feasible[k]){
1807 if (min_ind < 0){
1808 min_obj = SYM_INFINITY;
1809 for (l = best_can->child_num - 1; l >= 0; l--){
1810 if (best_can->feasible[l] && best_can->objval[k] <
1811 min_obj){
1812 min_obj = best_can->objval[l];
1813 min_ind = l;
1814 }
1815 }
1816 }
1817 if (min_ind > -1){
1818 if(can->objval[k] > best_can->objval[min_ind]){
1819 best_can->feasible[k] = TRUE;
1820 best_can->solutions[k] = can->solutions[k];
1821 can->solutions[k] = 0;
1822 min_ind = -1;
1823 }
1824 }
1825 }
1826 }
1827 } else{
1828 for (k = best_can->child_num - 1; k >= 0; k--){
1829 /* Again, this is only for tracking that there was a feasible
1830 solution discovered in presolve for display purposes */
1831 if (can->feasible[k]){
1832 best_can->feasible[k] = TRUE;
1833 }
1834 }
1835 }
1836 free_candidate(candidates + i);
1837 break;
1838 case SECOND_CANDIDATE_BETTER:
1839 case SECOND_CANDIDATE_BETTER_AND_BRANCH_ON_IT:
1840 #ifndef MAX_CHILDREN_NUM
1841 if (best_can == NULL){
1842 pobj = objval;
1843 pterm = termcode;
1844 pfeas = feasible;
1845 piter = iterd;
1846 #ifdef COMPILE_FRAC_BRANCHING
1847 pfrnum = frnum;
1848 pfrind = frind;
1849 pfrval = frval;
1850 #endif
1851 }else{
1852 pobj = best_can->objval;
1853 pterm = best_can->termcode;
1854 pfeas = best_can->feasible;
1855 piter = best_can->iterd;
1856 #ifdef COMPILE_FRAC_BRANCHING
1857 pfrnum = best_can->frac_num;
1858 pfrind = best_can->frac_ind;
1859 pfrval = best_can->frac_val;
1860 #endif
1861 }
1862 #endif
1863 if (best_can){
1864 if (p->par.keep_description_of_pruned == KEEP_IN_MEMORY){
1865 min_ind = -1;
1866 for (k = best_can->child_num - 1; k >= 0; k--){
1867 if (best_can->feasible[k]){
1868 if (min_ind < 0){
1869 min_obj = SYM_INFINITY;
1870 for (l = can->child_num - 1; l >= 0; l--){
1871 if (can->feasible[l] && can->objval[k] <
1872 min_obj){
1873 min_obj = can->objval[l];
1874 min_ind = l;
1875 }
1876 }
1877 }
1878 if (min_ind > -1){
1879 if(best_can->objval[k] > can->objval[min_ind]){
1880 can->feasible[k] = TRUE;
1881 can->solutions[k] = best_can->solutions[k];
1882 best_can->solutions[k] = 0;
1883 min_ind = -1;
1884 }
1885 }
1886 }
1887 }
1888 }
1889 else{
1890 for (k = can->child_num - 1; k >= 0; k--){
1891 /* Again, this is only for tracking that there was a feasible
1892 solution discovered in presolve for display purposes */
1893 if (best_can->feasible[k]){
1894 can->feasible[k] = TRUE;
1895 }
1896 }
1897 }
1898 free_candidate(&best_can);
1899 }
1900 best_can = can;
1901 candidates[i] = NULL;
1902 break;
1903 }
1904 if ((j & BRANCH_ON_IT)){
1905 break;
1906 }
1907 st_time += used_time(&total_time);
1908
1909 if (p->par.limit_strong_branching_time){
1910 should_continue_strong_branching(p,i,cand_num,st_time,total_iters,
1911 &should_continue);
1912 if (should_continue==FALSE) {
1913 PRINT(p->par.verbosity, 2,
1914 ("too much time in strong branching, breaking\n"));
1915 break;
1916 }
1917 }
1918 }
1919 }
1920 //printf ("total_iters = %d \n",total_iters);
1921 //printf ("candidates evaluated = %d \n",i);
1922
1923 if (should_use_hot_starts && unmark_hs) {
1924 unmark_hotstart(lp_data);
1925 set_itlim_hotstart(lp_data, -1);
1926 }else{
1927 load_basis(lp_data, cstat, rstat);
1928 }
1929 set_itlim(lp_data, -1);
1930
1931 #if 0
1932 if (best_can->type == CANDIDATE_VARIABLE &&
1933 vars[best_can->position]->lb == vars[best_can->position]->ub){
1934 printf("Error -- branching variable is already fixed. \n");
1935 printf("SYMPHONY has encountered numerical difficulties \n");
1936 printf("with the LP solver. Exiting...\n\n");
1937 }
1938
1939 if (best_can->type == CANDIDATE_VARIABLE &&
1940 best_can->rhs[0] < vars[best_can->position]->lb ||
1941 best_can->rhs[1] > vars[best_can->position]->ub){
1942 printf("Warning -- possible illegal branching. \n");
1943 printf("SYMPHONY has encountered possible numerical difficulties \n");
1944 printf("with the LP solver. Exiting...\n\n");
1945 }
1946 #endif
1947
1948 #ifndef MAX_CHILDREN_NUM
1949 FREE(pobj); FREE(pterm); FREE(pfeas); FREE(piter);
1950 # ifdef COMPILE_FRAC_BRANCHING
1951 FREE(pfrnum); FREE(pfrind); FREE(pfrval);
1952 # endif
1953 #endif
1954
1955 if (p->par.max_presolve_iter > 0) {
1956 if (should_use_hot_starts == TRUE) {
1957 set_itlim_hotstart(lp_data, -1);
1958 }
1959 set_itlim(lp_data, -1); // both limits should be set if using hotstarts
1960 }
1961
1962 #ifdef STATISTICS
1963 PRINT(p->par.verbosity, 5,
1964 ("Itlim reached %i times out of %i .\n\n", itlim, cnum));
1965 #endif
1966
1967 if (should_use_rel_br != TRUE) {
1968 for (i++; i<cand_num; i++){
1969 /* Free the remaining candidates */
1970 free_candidate(candidates + i);
1971 }
1972 }
1973 FREE(candidates);
1974
1975 if (p->par.keep_description_of_pruned == KEEP_IN_MEMORY){
1976 indices = lp_data->tmp.i1;
1977 values = lp_data->tmp.d;
1978 for (k = best_can->child_num - 1; k >= 0; k--){
1979 if (best_can->feasible[k]){
1980 best_can->sol_sizes[k] = collect_nonzeros(p,
1981 best_can->solutions[k],
1982 indices, values);
1983 FREE(best_can->solutions[k]);
1984 best_can->sol_inds[k] = (int *) malloc(best_can->sol_sizes[k]*
1985 ISIZE);
1986 best_can->solutions[k] = (double *) malloc(best_can->sol_sizes[k]*
1987 DSIZE);
1988 memcpy(best_can->sol_inds[k], indices, best_can->sol_sizes[k] *
1989 ISIZE);
1990 memcpy(best_can->solutions[k], values, best_can->sol_sizes[k]*
1991 DSIZE);
1992 break;
1993 }
1994 }
1995 }
1996
1997 *candidate = best_can;
1998
1999 p->comp_times.strong_branching += used_time(&p->tt);
2000 send_node_desc(p, NODE_BRANCHED_ON);
2001 p->comp_times.communication += used_time(&p->tt);
2002 return(DO_BRANCH);
2003 }
2004
2005 /*===========================================================================*/
2006
2007 int branch(lp_prob *p, int cuts)
2008 {
2009 LPdata *lp_data = p->lp_data;
2010 branch_obj *can;
2011 char *action;
2012 int branch_var, branch_row, keep;
2013 var_desc *var;
2014 cut_data *cut;
2015 node_desc *desc;
2016 int termcode;
2017 #if defined(DO_TESTS) && defined(COMPILE_IN_LP)
2018 branch_obj *bobj = &(p->tm->active_nodes[p->proc_index]->bobj);
2019 #endif
2020
2021 termcode = select_branching_object(p, &cuts, &can);
2022
2023 if (termcode == ERROR__NO_BRANCHING_CANDIDATE){
2024 return(termcode);
2025 }
2026
2027 if (can == NULL){
2028 if (termcode == DO_NOT_BRANCH__FEAS_SOL) {
2029 /* a better feasible solution was found. return to do reduced cost
2030 * fixing etc.
2031 */
2032 return(FEAS_SOL_FOUND);
2033 }
2034 /* We were either able to fathom the node or found violated cuts
2035 * In any case, send the qualifying cuts to the cutpool */
2036 p->comp_times.strong_branching += used_time(&p->tt);
2037 #pragma omp critical(cut_pool)
2038 send_cuts_to_pool(p, p->par.eff_cnt_before_cutpool);
2039 p->comp_times.communication += used_time(&p->tt);
2040 return (termcode == DO_NOT_BRANCH__FATHOMED ? BRANCHING_INF_NODE : cuts);
2041 }
2042
2043 /*------------------------------------------------------------------------*\
2044 * Now we evaluate can, the best of the candidates.
2045 \*------------------------------------------------------------------------*/
2046 action = lp_data->tmp.c; /* n (good estimate... can->child_num) */
2047 if ((termcode = select_child_u(p, can, action)) < 0)
2048 return(termcode);
2049 if (p->par.verbosity > 4)
2050 print_branch_stat_u(p, can, action);
2051
2052 for (keep = can->child_num-1; keep >= 0; keep--)
2053 if (action[keep] == KEEP_THIS_CHILD) break;
2054
2055 /* Send the branching information to the TM and inquire whether we
2056 should dive */
2057 p->comp_times.strong_branching += used_time(&p->tt);
2058 /* 'keep' may be modified if children are pruned, but we need the original
2059 value */
2060 int old_keep = keep;
2061 send_branching_info(p, can, action, &keep);
2062 p->comp_times.communication += used_time(&p->tt);
2063
2064 /* If we don't dive then return quickly */
2065 if (keep < 0 || p->dive == DO_NOT_DIVE){
2066 free_candidate_completely(&can);
2067 return(FATHOMED_NODE);
2068 }
2069
2070 #if defined(DO_TESTS) && defined(COMPILE_IN_LP)
2071 assert(can->rhs[old_keep] == bobj->rhs[keep]);
2072 #endif
2073
2074 desc = p->desc;
2075 switch (can->type){
2076 case CANDIDATE_VARIABLE:
2077 p->branch_var = can->position;
2078 p->branch_dir = can->sense[old_keep];
2079 var = lp_data->vars[branch_var = can->position];
2080 switch (can->sense[old_keep]){
2081 case 'E':
2082 var->new_lb = var->new_ub = can->rhs[old_keep];
2083 var->lb = var->ub = can->rhs[old_keep]; break;
2084 case 'R':
2085 var->new_lb = can->rhs[old_keep];
2086 var->new_ub = var->lb + can->range[old_keep];
2087 var->lb = can->rhs[old_keep]; var->ub = var->lb + can->range[old_keep]; break;
2088 case 'L':
2089 var->new_ub = can->rhs[old_keep];
2090 var->ub = can->rhs[old_keep]; break;
2091 case 'G':
2092 var->new_lb = can->rhs[old_keep];
2093 var->lb = can->rhs[old_keep]; break;
2094 }
2095 //printf("branching on %i %c %f %f\n", branch_var, can->sense[old_keep], var->lb, var->ub);
2096 change_col(lp_data, branch_var, can->sense[old_keep], var->lb, var->ub);
2097 lp_data->status[branch_var] |= VARIABLE_BRANCHED_ON;
2098 break;
2099 case SOS1_IMPLICIT:
2100 for(int j = 0; j < can->sos_cnt[old_keep]; j++){
2101 branch_var = can->sos_ind[old_keep][j];
2102 change_ub(lp_data, branch_var, 0.0);
2103 lp_data->vars[branch_var]->new_ub = 0.0;
2104 lp_data->vars[branch_var]->ub = 0.0;
2105 //printf("%i ", branch_var);
2106 lp_data->status[branch_var] |= VARIABLE_BRANCHED_ON;
2107 }
2108 //printf("\n");
2109 break;
2110 case CANDIDATE_CUT_IN_MATRIX:
2111 branch_row = can->position;
2112 cut = lp_data->rows[branch_row].cut;
2113 /* To maintain consistency with TM we have to fix a few more things if
2114 we had a non-base, new branching cut */
2115 if (branch_row >= p->base.cutnum && !(cut->branch & CUT_BRANCHED_ON)){
2116 /* insert cut->name into p->desc.cutind.list, and insert SLACK_BASIC
2117 to the same position in p->desc.basis.extrarows.stat */
2118 #ifdef DO_TESTS
2119 if (desc->cutind.size != desc->basis.extrarows.size){
2120 printf("Oops! desc.cutind.size != desc.basis.extrarows.size! \n");
2121 exit(-123);
2122 }
2123 #endif
2124 #ifdef COMPILE_IN_LP
2125 /* Because these cuts are shared with the treemanager, we have to
2126 make a copy before changing them if the LP is compiled in */
2127 cut = (cut_data *) malloc(sizeof(cut_data));
2128 memcpy((char *)cut, (char *)lp_data->rows[branch_row].cut,
2129 sizeof(cut_data));
2130 if (cut->size){
2131 cut->coef = (char *) malloc(cut->size);
2132 memcpy((char *)cut->coef,
2133 (char *)lp_data->rows[branch_row].cut->coef, cut->size);
2134 }
2135 lp_data->rows[branch_row].cut = cut;
2136 #endif
2137 if (desc->cutind.size == 0){
2138 desc->cutind.size = 1;
2139 desc->cutind.list = (int *) malloc(ISIZE);
2140 desc->cutind.list[0] = cut->name;
2141 desc->basis.extrarows.size = 1; /* this must have been 0, too */
2142 desc->basis.extrarows.stat = (int *) malloc(ISIZE);
2143 desc->basis.extrarows.stat[0] = SLACK_BASIC;
2144 }else{
2145 int i, name = cut->name;
2146 int *list;
2147 int *stat;
2148 /* most of the time the one extra element will fit into the
2149 already allocated memory chunk, so it's worth to realloc */
2150 desc->cutind.size++;
2151 list = desc->cutind.list =
2152 (int *) realloc(desc->cutind.list, desc->cutind.size * ISIZE);
2153 desc->basis.extrarows.size++;
2154 stat = desc->basis.extrarows.stat =
2155 (int *) realloc(desc->basis.extrarows.stat,
2156 desc->cutind.size * ISIZE);
2157 for (i = desc->cutind.size - 1; i > 0; i--){
2158 #ifdef DO_TESTS
2159 if (name == list[i-1]){
2160 printf("Oops! name == desc.cutind.list[i] !\n");
2161 exit(-124);
2162 }
2163 #endif
2164 if (name < list[i-1]){
2165 list[i] = list[i-1];
2166 stat[i] = stat[i-1];
2167 }else{
2168 break;
2169 }
2170 }
2171 list[i] = name;
2172 stat[i] = SLACK_BASIC;
2173 }
2174 }
2175 cut->rhs = can->rhs[old_keep];
2176 if ((cut->sense = can->sense[old_keep]) == 'R')
2177 cut->range = can->range[old_keep];
2178 cut->branch = CUT_BRANCHED_ON | can->branch[old_keep];
2179 constrain_row_set(lp_data, 1, &branch_row);
2180 lp_data->rows[branch_row].free = FALSE;
2181 break;
2182 }
2183
2184 /* Since this is a child we dived into, we know that TM stores the stati of
2185 extra vars/rows wrt the parent */
2186 p->desc->basis.extravars.type = WRT_PARENT;
2187 p->desc->basis.extrarows.type = WRT_PARENT;
2188
2189 free_candidate_completely(&can);
2190
2191 /* the new p->bc_index is received in send_branching_info() */
2192 p->bc_level++;
2193 /*p->iter_num = 0;*/
2194
2195 return(NEW_NODE);
2196 }
2197
2198 /*===========================================================================*/
2199
2200 int col_gen_before_branch(lp_prob *p, int *new_vars)
2201 {
2202 our_col_set *new_cols;
2203 int dual_feas;
2204
2205 check_ub(p);
2206 if (! p->has_ub ||
2207 (p->colgen_strategy & BEFORE_BRANCH__DO_NOT_GENERATE_COLS) ||
2208 (p->lp_data->nf_status & NF_CHECK_NOTHING))
2209 return(DO_BRANCH);
2210
2211 PRINT(p->par.verbosity, 2, ("Generating cols before branching.\n"));
2212 p->comp_times.strong_branching += used_time(&p->tt);
2213 new_cols = price_all_vars(p);
2214 p->comp_times.pricing += used_time(&p->tt);
2215 /*price_all_vars sorts by user_ind. We need things sorted by colind */
2216 colind_sort_extra(p);
2217 *new_vars = new_cols->num_vars + new_cols->rel_ub + new_cols->rel_lb;
2218 dual_feas = new_cols->dual_feas;
2219 free_col_set(&new_cols);
2220 check_ub(p);
2221 if (dual_feas == NOT_TDF){
2222 return(DO_NOT_BRANCH);
2223 }else{
2224 if (p->ub - p->par.granularity < p->lp_data->objval ||
2225 p->lp_data->termcode == LP_D_OBJLIM ||
2226 p->lp_data->termcode == LP_OPT_FEASIBLE){
2227 /* If total dual feas and high cost or feasibility ==> fathomable */
2228 PRINT(p->par.verbosity, 1, ("Managed to fathom the node.\n"));
2229 send_node_desc(p, p->lp_data->termcode == LP_OPT_FEASIBLE ?
2230 FEASIBLE_PRUNED : OVER_UB_PRUNED);
2231 p->comp_times.communication += used_time(&p->tt);
2232 return(DO_NOT_BRANCH__FATHOMED);
2233 }else{
2234 return(DO_BRANCH); /* if we got here, then DO_BRANCH */
2235 }
2236 }
2237 return(DO_BRANCH); /* fake return */
2238 }
2239
2240 /*===========================================================================*/
2241
2242 /*****************************************************************************/
2243 /* This is a generic function */
2244 /*****************************************************************************/
2245
2246 void branch_close_to_half(lp_prob *p, int max_cand_num, int *cand_num,
2247 branch_obj ***candidates)
2248 {
2249 LPdata *lp_data = p->lp_data;
2250 double *x = lp_data->x;
2251 //double lpetol100 = lp_data->lpetol*100, lpetol1 = 1 - lpetol100;
2252 double lpetol100 = lp_data->lpetol, lpetol1 = 1 - lpetol100;
2253 int *xind = lp_data->tmp.i1; /* n */
2254 double fracx, *xval = lp_data->tmp.d; /* n */
2255 branch_obj *cand;
2256 int i, j, cnt = 0;
2257 double lim[7] = {.1, .15, .20, .233333, .266667, .3, 1};
2258 var_desc **vars = lp_data->vars;
2259 int should_use_rel_br = p->par.should_use_rel_br;
2260
2261 /* first get the fractional values */
2262 if (should_use_rel_br == TRUE) {
2263 xind = p->br_rel_cand_list;
2264 }
2265 double frac_avg = 0.0;
2266
2267 const CoinPackedMatrix *matrixByCol = lp_data->si->getMatrixByCol();
2268
2269 for (i = lp_data->n-1; i >= 0; i--){
2270 if (vars[i]->is_int){
2271 if (x[i] > vars[i]->new_lb && x[i] < vars[i]->new_ub){
2272 fracx = x[i] - floor(x[i]);
2273 if (fracx > lpetol100 && fracx < lpetol1){
2274 xind[cnt] = i;
2275 int collen = matrixByCol->getVectorSize(i);
2276 //xval[cnt++] = fabs(fracx - .5);
2277 xval[cnt++] = -collen*(0.5 - fabs(fracx - .5));
2278 frac_avg += 0.5 - fabs(fracx - .5);
2279 }
2280 }
2281 }
2282 *cand_num = cnt;
2283 }
2284
2285 #ifdef COMPILE_IN_LP
2286 p->tm->active_nodes[p->proc_index]->frac_cnt = cnt;
2287 p->tm->active_nodes[p->proc_index]->frac_avg = frac_avg;
2288 #endif
2289
2290 if (should_use_rel_br == TRUE) {
2291 *candidates = (branch_obj **) malloc(1 * sizeof(branch_obj *));
2292 cand = (*candidates)[0] = (branch_obj *) calloc(1, sizeof(branch_obj) );
2293 cand->type = CANDIDATE_VARIABLE;
2294 cand->child_num = 2;
2295 cand->sense[0] = 'L';
2296 cand->sense[1] = 'G';
2297 cand->range[0] = cand->range[1] = 0;
2298 qsort_di(xval, xind, cnt);
2299 } else {
2300 qsort_di(xval, xind, cnt);
2301 if (p->bc_level>p->par.strong_br_all_candidates_level ||
2302 p->par.user_set_strong_branching_cand_num) {
2303 for (j = 0, i = 0; i < cnt;){
2304 if (xval[i] > lim[j]){
2305 if (i == 0){
2306 j++; continue;
2307 }else{
2308 break;
2309 }
2310 }else{
2311 i++;
2312 }
2313 }
2314 cnt = i;
2315 *cand_num = MIN(max_cand_num, cnt);
2316 } else {
2317 *cand_num = cnt;
2318 }
2319
2320 if (!*candidates)
2321 *candidates = (branch_obj **) malloc(*cand_num * sizeof(branch_obj *));
2322 for (i=*cand_num-1; i>=0; i--){
2323 cand = (*candidates)[i] = (branch_obj *) calloc(1, sizeof(branch_obj) );
2324 cand->type = CANDIDATE_VARIABLE;
2325 cand->child_num = 2;
2326 cand->position = xind[i];
2327 cand->sense[0] = 'L';
2328 cand->sense[1] = 'G';
2329 cand->rhs[0] = floor(x[xind[i]]);
2330 cand->rhs[1] = cand->rhs[0] + 1;
2331 cand->range[0] = cand->range[1] = 0;
2332 }
2333 }
2334 }
2335
2336 /*===========================================================================*/
2337
2338 /*****************************************************************************/
2339 /* This is a generic function */
2340 /*****************************************************************************/
2341
2342 void branch_close_to_half_and_expensive(lp_prob *p, int max_cand_num,
2343 int *cand_num, branch_obj ***candidates)
2344 {
2345 LPdata *lp_data = p->lp_data;
2346 double *x = lp_data->x;
2347 double lpetol = lp_data->lpetol, lpetol1 = 1 - lpetol;
2348 int *xind = lp_data->tmp.i1; /* n */
2349 double fracx, *xval = lp_data->tmp.d; /* n */
2350 branch_obj *cand;
2351 int i, j, cnt = 0;
2352 double lim[7] = {.1, .15, .20, .233333, .266667, .3, 1};
2353
2354 /* first get the fractional values */
2355 for (i = lp_data->n-1; i >= 0; i--){
2356 fracx = x[i] - floor(x[i]);
2357 if (fracx > lpetol && fracx < lpetol1){
2358 xind[cnt] = i;
2359 xval[cnt++] = fabs(fracx - .5);
2360 }
2361 }
2362
2363 qsort_di(xval, xind, cnt);
2364
2365 for (j=0, i=0; i<cnt; i++){
2366 if (xval[i] > lim[j]){
2367 if (i == 0){
2368 j++; continue;
2369 }else{
2370 break;
2371 }
2372 }
2373 }
2374 cnt = i;
2375
2376 if (max_cand_num >= cnt){
2377 *cand_num = cnt;
2378 }else{
2379 for (i=cnt-1; i>=0; i--){
2380 get_objcoef(p->lp_data, xind[i], xval+i);
2381 xval[i] *= -1;
2382 }
2383 qsort_di(xval, xind, cnt);
2384 *cand_num = max_cand_num;
2385 }
2386
2387 if (!*candidates)
2388 *candidates = (branch_obj **) malloc(*cand_num * sizeof(branch_obj *));
2389 for (i=*cand_num-1; i>=0; i--){
2390 cand = (*candidates)[i] = (branch_obj *) calloc(1, sizeof(branch_obj) );
2391 cand->type = CANDIDATE_VARIABLE;
2392 cand->child_num = 2;
2393 cand->position = xind[i];
2394 cand->sense[0] = 'L';
2395 cand->sense[1] = 'G';
2396 cand->rhs[0] = floor(x[xind[i]]);
2397 cand->rhs[1] = cand->rhs[0] + 1;
2398 cand->range[0] = cand->range[1] = 0;
2399 }
2400 }
2401
2402 /*===========================================================================*/
2403
2404 /*****************************************************************************/
2405 /* This works only for 0/1 problems!!! */
2406 /*****************************************************************************/
2407
2408 void branch_close_to_one_and_cheap(lp_prob *p, int max_cand_num, int *cand_num,
2409 branch_obj ***candidates)
2410 {
2411 LPdata *lp_data = p->lp_data;
2412 double *x = lp_data->x;
2413 double lpetol = lp_data->lpetol, lpetol1 = 1 - lpetol;
2414 int *xind = lp_data->tmp.i1; /* n */
2415 double *xval = lp_data->tmp.d; /* n */
2416 branch_obj *cand;
2417 int i, j, cnt = 0;
2418 double lim[8] = {.1, .2, .25, .3, .333333, .366667, .4, 1};
2419
2420 /* first get the fractional values */
2421 for (i = lp_data->n-1; i >= 0; i--){
2422 if (x[i] > lpetol && x[i] < lpetol1){
2423 xind[cnt] = i;
2424 xval[cnt++] = 1 - x[i];
2425 }
2426 }
2427 qsort_di(xval, xind, cnt);
2428
2429 for (j=0, i=0; i<cnt; i++){
2430 if (xval[i] > lim[j]){
2431 if (i == 0){
2432 j++; continue;
2433 }else{
2434 break;
2435 }
2436 }
2437 }
2438 cnt = i;
2439
2440 if (max_cand_num >= cnt){
2441 *cand_num = cnt;
2442 }else{
2443 for (i=cnt-1; i>=0; i--){
2444 get_objcoef(p->lp_data, xind[i], xval+i);
2445 }
2446 qsort_di(xval, xind, cnt);
2447 *cand_num = max_cand_num;
2448 }
2449
2450 if (!*candidates)
2451 *candidates = (branch_obj **) malloc(*cand_num * sizeof(branch_obj *));
2452 for (i=*cand_num-1; i>=0; i--){
2453 cand = (*candidates)[i] = (branch_obj *) calloc(1, sizeof(branch_obj) );
2454 cand->type = CANDIDATE_VARIABLE;
2455 cand->child_num = 2;
2456 cand->position = xind[i];
2457 cand->sense[0] = 'L';
2458 cand->sense[1] = 'G';
2459 cand->rhs[0] = floor(x[xind[i]]);
2460 cand->rhs[1] = cand->rhs[0] + 1;
2461 cand->range[0] = cand->range[1] = 0;
2462 }
2463 }
2464
2465 /*===========================================================================*/
2466 int should_continue_strong_branching(lp_prob *p, int i, int cand_num,
2467 double st_time, int total_iters,
2468 int *should_continue)
2469 {
2470 double allowed_time = 0;
2471 *should_continue = TRUE;
2472 int min_cands;
2473 int verbosity = p->par.verbosity;
2474 if (p->bc_level<1) {
2475 allowed_time = 20*p->comp_times.lp/p->iter_num;
2476 //allowed_iter = 20*p->lp_stat.lp_total_iter_num/(p->iter_num + 1);
2477 if (allowed_time < 2) {
2478 allowed_time = 2;
2479 }
2480 //allowed_iter = MAX(allowed_iter, 1000);
2481 min_cands = MIN(cand_num,p->par.strong_branching_cand_num_max);
2482 } else {
2483 allowed_time = p->comp_times.lp/2 - p->comp_times.strong_branching;
2484 //allowed_iter = (int)((p->lp_stat.lp_total_iter_num -
2485 // p->lp_stat.str_br_total_iter_num)/2.0);
2486 min_cands = MIN(cand_num,p->par.strong_branching_cand_num_min);
2487 }
2488 PRINT(verbosity,10,("allowed_time = %f\n",allowed_time));
2489 if (st_time/(i+1)*cand_num < allowed_time) {
2490 /* all cands can be evaluated in given time */
2491 *should_continue = TRUE;
2492 } else if (i >= min_cands-1 && st_time>allowed_time) {
2493 /* time is up and min required candidates have been evaluated */
2494 *should_continue = FALSE;
2495 } else if (p->par.user_set_max_presolve_iter == TRUE) {
2496 /* user specified a limit and we wont change it */
2497 *should_continue = TRUE;
2498 } else {
2499 /* we will not be able to evaluate all candidates in given time. we
2500 * reduce the number of iterations */
2501 double min_iters =
2502 (allowed_time-st_time)*total_iters/st_time/(cand_num-i+1);
2503 if (min_iters<10) {
2504 /*
2505 * cant evaluate all candidates in given time with just 10 iters
2506 * as well. we have a choice: increase iters and do min_cands
2507 * or do ten iters and try to evaluate max possible num of cands.
2508 * we like the second option more.
2509 */
2510 min_iters = 10;
2511 }
2512 if (p->par.use_hot_starts && !p->par.branch_on_cuts) {
2513 set_itlim_hotstart(p->lp_data, (int) min_iters);
2514 set_itlim(p->lp_data, (int) min_iters);
2515 } else {
2516 set_itlim(p->lp_data, (int) min_iters);
2517 }
2518 PRINT(verbosity,6, ("iteration limit set to %d\n", (int )min_iters));
2519 *should_continue = TRUE;
2520 }
2521 PRINT(verbosity,29, ("strong branching i = %d\n",i));
2522 return 0;
2523 }
2524
2525 /*===========================================================================*/
2526 int strong_branch(lp_prob *p, int branch_var, double lb, double ub,
2527 double new_lb, double new_ub, double *obj, int should_use_hot_starts,
2528 int *termstatus, int *iterd, int sos_cnt, int *sos_ind)
2529 {
2530 int status = 0;
2531 LPdata *lp_data = p->lp_data;
2532 int *cstat = lp_data->tmp.i1;
2533 int *rstat = lp_data->tmp.i2;
2534
2535 // TODO: LP_ABANDONED
2536 /* change the lb and ub */
2537 if(sos_cnt < 1){
2538 change_lbub(lp_data, branch_var, new_lb, new_ub);
2539 }else{
2540 for(int i = 0; i < sos_cnt; i++){
2541 change_lbub(lp_data, sos_ind[i], 0.0, 0.0);
2542 }
2543 }
2544
2545 // if (p->par.use_hot_starts && !p->par.branch_on_cuts) {
2546 if (should_use_hot_starts) {
2547 *termstatus = solve_hotstart(lp_data, iterd);
2548 } else {
2549 load_basis(lp_data, cstat, rstat);
2550 *termstatus = dual_simplex(lp_data, iterd);
2551 }
2552
2553 if (*termstatus == LP_D_INFEASIBLE || *termstatus == LP_D_OBJLIM ||
2554 *termstatus == LP_D_UNBOUNDED) {
2555 *obj = SYM_INFINITY;
2556 if(sos_cnt < 1){
2557 p->lp_stat.str_br_bnd_changes++;
2558 }
2559 } else {
2560 *obj = lp_data->objval;
2561 // if(lp_data->objval < *obj - lp_data->lpetol){
2562 // printf("dual_simplex error: %i %i\n", p->bc_index, branch_var);
2563 // }else{
2564 // *obj = lp_data->objval;
2565 // }
2566
2567 if (*termstatus == LP_OPTIMAL) {
2568 if (!p->has_ub || *obj < p->ub - p->par.granularity + lp_data->lpetol) {
2569 is_feasible_u(p, TRUE, TRUE);
2570 } else {
2571 *obj = SYM_INFINITY;
2572 *termstatus = LP_D_OBJLIM;
2573 if(sos_cnt < 1){
2574 p->lp_stat.str_br_bnd_changes++;
2575 }
2576 }
2577 } else if (*termstatus == LP_ABANDONED) {
2578 status = LP_ABANDONED;
2579 }
2580 }
2581 p->lp_stat.lp_calls++;
2582 p->lp_stat.str_br_lp_calls++;
2583 p->lp_stat.str_br_total_iter_num += *iterd;
2584 p->lp_stat.num_str_br_cands_in_path++;
2585
2586 if(sos_cnt < 1){
2587 change_lbub(lp_data, branch_var, lb, ub);
2588 }else{
2589 for(int i = 0; i < sos_cnt; i++){
2590 change_lbub(lp_data, sos_ind[i], 0.0, 1.0);
2591 }
2592 }
2593 return status;
2594 }
2595
2596 /*===========================================================================*/
2597 /*===========================================================================*/
2598
2599 int prep_col_fixable(double xval, double aval, double c_lb, double c_ub,
2600 double row_lb, double row_ub, double si_row_lb, double si_row_ub,
2601 double *col_fixed_lb, double *col_fixed_ub, double etol,
2602 double inf)
2603 {
2604
2605 if(xval < c_lb + etol){
2606 if(prep_row_violated(row_lb, row_ub, si_row_lb, si_row_ub,
2607 aval, c_lb, c_ub,
2608 c_lb + 1.0, c_ub, etol, inf)){
2609 *col_fixed_lb = *col_fixed_ub = c_lb;
2610 return TRUE;
2611 }
2612 }else if(xval > c_ub - etol){
2613 if(prep_row_violated(row_lb, row_ub, si_row_lb, si_row_ub,
2614 aval, c_lb, c_ub,
2615 c_lb, c_ub - 1.0, etol, inf)){
2616 *col_fixed_lb = *col_fixed_ub = c_ub;
2617 return TRUE;
2618 }
2619 }else{
2620 /* fractional var
2621 try fixing upper and lower
2622 -might catch infeasibility here
2623 */
2624 double floorx = floor(xval);
2625 double ceilx = ceil(xval);
2626 int status = FALSE;
2627 if(prep_row_violated(row_lb, row_ub, si_row_lb, si_row_ub,
2628 aval, c_lb, c_ub,
2629 c_lb, floorx, etol, inf)){
2630 *col_fixed_lb = ceilx;
2631 *col_fixed_ub = c_ub;
2632 status = TRUE;
2633 }
2634 if(prep_row_violated(row_lb, row_ub, si_row_lb, si_row_ub,
2635 aval, c_lb, c_ub,
2636 ceilx, c_ub, etol, inf)){
2637 if(status){
2638 *col_fixed_lb = c_ub + 1.0;
2639 *col_fixed_ub = c_ub;
2640 }else{
2641 *col_fixed_lb = c_lb;
2642 *col_fixed_ub = floorx;
2643 status = TRUE;
2644 }
2645 }
2646 return status;
2647 }
2648
2649 return FALSE;
2650 }
2651
2652 /*===========================================================================*/
2653 /*===========================================================================*/
2654
2655 int prep_row_violated(double row_lb, double row_ub, double si_row_lb, double si_row_ub,
2656 double aval, double old_col_lb, double old_col_ub,
2657 double new_col_lb, double new_col_ub, double etol,
2658 double inf)
2659 {
2660
2661 if(aval >= 0.0){
2662 if((row_lb > -inf && si_row_ub < inf &&
2663 row_lb + aval*(new_col_lb - old_col_lb) > si_row_ub + etol) ||
2664 (row_ub < inf && si_row_lb > -inf &&
2665 row_ub + aval*(new_col_ub - old_col_ub) < si_row_lb - etol))
2666 return TRUE;
2667 }else{
2668 if((row_lb > -inf && si_row_ub < inf &&
2669 row_lb + aval*(new_col_ub - old_col_ub) > si_row_ub + etol) ||
2670 (row_ub < inf && si_row_lb > -inf &&
2671 row_ub + aval*(new_col_lb - old_col_lb) < si_row_lb - etol))
2672 return TRUE;
2673 }
2674
2675 return FALSE;
2676 }
2677
2678 /*===========================================================================*/
2679 /*===========================================================================*/
2680
2681 int prep_tighten_bounds(LPdata *lp_data, //int cand_num, int *cand_ind,
2682 int *num_changes, double *bnd_val, int *bnd_ind, char *bnd_sense,
2683 double *row_ub, double *row_lb, char *cand_fixed)
2684 {
2685 int j, k, r_ind, c_ind, col_start, col_end;
2686 int n = lp_data->n, m = lp_data->m;
2687 var_desc **vars = lp_data->vars;
2688 double etol = lp_data->lpetol, coeff;
2689 *num_changes = 0;
2690 *cand_fixed = FALSE;
2691
2692 const double *ub = lp_data->si->getColUpper();
2693 const double *lb = lp_data->si->getColLower();
2694 //const double *obj = lp_data->si->getObjCoefficients();
2695
2696 const CoinPackedMatrix * matrix = lp_data->si->getMatrixByCol();
2697 const double *matval = matrix->getElements();
2698 const int *matind = matrix->getIndices();
2699 const int *matbeg = matrix->getVectorStarts();
2700 const int *len = matrix->getVectorLengths();
2701
2702 //double * r_ub = (double *)malloc (lp_data->m*DSIZE);
2703 //double * r_lb = (double *)malloc (lp_data->m*DSIZE);
2704 //memcpy(r_ub, const_cast<double*>(lp_data->si->getRowUpper()), DSIZE*lp_data->m);
2705 //memcpy(r_lb, const_cast<double*>(lp_data->si->getRowLower()), DSIZE*lp_data->m);
2706 const double *si_row_ub = lp_data->si->getRowUpper();
2707 const double *si_row_lb = lp_data->si->getRowLower();
2708 const double inf = lp_data->si->getInfinity();
2709
2710 int iter_cnt = 0, iter_limit = 2;
2711 int bounds_updated;
2712 int col_fixed;
2713 /* get row bounds */
2714 char *row_fixed = lp_data->tmp.c;
2715 char *row_changed = lp_data->tmp.c + m;
2716
2717 for(r_ind = 0; r_ind < m; r_ind++){
2718 row_ub[r_ind] = row_lb[r_ind] = 0.0;
2719 row_fixed[r_ind] = row_changed[r_ind] = TRUE;
2720 }
2721
2722 for(c_ind = 0; c_ind < n; c_ind++){
2723 col_fixed = FALSE;
2724 if(ub[c_ind] < lb[c_ind] + etol) col_fixed = TRUE;
2725 col_start = matbeg[c_ind];
2726 col_end = col_start + len[c_ind];
2727 for(j = col_start; j < col_end; j++){
2728 r_ind = matind[j];
2729 coeff = matval[j];
2730 if(row_ub[r_ind] < inf){
2731 if(coeff >= 0.0 && ub[c_ind] < inf){
2732 row_ub[r_ind] += ub[c_ind]*coeff;
2733 }else if(coeff < 0.0 && lb[c_ind] > -inf){
2734 row_ub[r_ind] += lb[c_ind]*coeff;
2735 }else{
2736 row_ub[r_ind] = inf;
2737 }
2738 }
2739 if(row_lb[r_ind] > -inf){
2740 if(coeff >= 0.0 && lb[c_ind] > -inf){
2741 row_lb[r_ind] += lb[c_ind]*coeff;
2742 }else if(coeff < 0.0 && ub[c_ind] < inf){
2743 row_lb[r_ind] += ub[c_ind]*coeff;
2744 }else{
2745 row_lb[r_ind] = -inf;
2746 }
2747 }
2748 if(row_fixed[r_ind] && !col_fixed) row_fixed[r_ind] = FALSE;
2749 }
2750 }
2751
2752 //return 0;
2753
2754 matrix = lp_data->si->getMatrixByRow();
2755 const double *r_matval = matrix->getElements();
2756 const int *r_matind = matrix->getIndices();
2757 const int *r_matbeg = matrix->getVectorStarts();
2758 const int *r_len = matrix->getVectorLengths();
2759 char row_fix_dir;
2760 int row_start, row_end, rec_row_ind;
2761 double col_fix_lb, col_fix_ub, old_lb, old_ub, rec_coeff;
2762
2763 while(iter_cnt < iter_limit){
2764 bounds_updated = FALSE;
2765 for(r_ind = 0; r_ind < m; r_ind++){
2766 if(!row_fixed[r_ind] && row_changed[r_ind]){
2767 row_changed[r_ind] = FALSE;
2768 row_fix_dir = 'F';
2769 if(row_ub[r_ind] < si_row_lb[r_ind] + etol){
2770 row_fix_dir = 'U';
2771 row_fixed[r_ind] = TRUE;
2772 }else if(row_lb[r_ind] > si_row_ub[r_ind] - etol){
2773 row_fix_dir = 'L';
2774 row_fixed[r_ind] = TRUE;
2775 }
2776 row_start = r_matbeg[r_ind];
2777 row_end = row_start + r_len[r_ind];
2778 for(j = row_start; j < row_end; j++){
2779 c_ind = r_matind[j];
2780 coeff = r_matval[j];
2781 if(vars[c_ind]->ub > vars[c_ind]->lb + etol){
2782 col_fixed = FALSE;
2783 col_fix_lb = col_fix_ub = 0.0;
2784 if(row_fix_dir == 'U'){
2785 continue;
2786 col_fixed = TRUE;
2787 if(coeff >= 0.0) {
2788 col_fix_ub = col_fix_lb = vars[c_ind]->ub;
2789 }else{
2790 col_fix_ub = col_fix_lb = vars[c_ind]->lb;
2791 }
2792 }else if(row_fix_dir == 'L'){
2793 continue;
2794 col_fixed = TRUE;
2795 if(coeff >= 0.0)
2796 col_fix_ub = col_fix_lb = vars[c_ind]->lb;
2797 else
2798 col_fix_ub = col_fix_lb = vars[c_ind]->ub;
2799 }else if(vars[c_ind]->is_int){
2800 /* try fixing if integer*/
2801 //xval = lp_data->x[c_ind];
2802 if(prep_col_fixable(lp_data->x[c_ind], coeff, vars[c_ind]->lb, vars[c_ind]->ub,
2803 row_lb[r_ind], row_ub[r_ind], si_row_lb[r_ind], si_row_ub[r_ind],
2804 &col_fix_lb, &col_fix_ub, etol, inf)){
2805 if(col_fix_lb > col_fix_ub + etol){
2806 return PREP_INFEAS;
2807 }else{
2808 col_fixed = TRUE;
2809 if(*cand_fixed == FALSE && lp_data->x[c_ind] > vars[c_ind]->lb + etol &&
2810 lp_data->x[c_ind] < vars[c_ind]->ub - etol){
2811 *cand_fixed = TRUE;
2812 }
2813 }
2814 }
2815 }
2816 if(col_fixed){
2817 if(col_fix_ub < vars[c_ind]->ub - etol){
2818 bnd_val[*num_changes] = col_fix_ub;
2819 bnd_sense[*num_changes] = 'L';
2820 bnd_ind[*num_changes] = c_ind;
2821 (*num_changes)++;
2822 }else if(col_fix_lb > vars[c_ind]->lb + etol){
2823 bnd_val[*num_changes] = col_fix_lb;
2824 bnd_sense[*num_changes] = 'U';
2825 bnd_ind[*num_changes] = c_ind;
2826 (*num_changes)++;
2827 }else{
2828 printf("error -- prep_tighten_bounds while branching...\n");
2829 return 0;
2830 }
2831
2832 old_lb = vars[c_ind]->lb;
2833 old_ub = vars[c_ind]->ub;
2834 change_lbub(lp_data, c_ind, col_fix_lb, col_fix_ub);
2835 vars[c_ind]->new_ub = col_fix_ub;
2836 vars[c_ind]->ub = col_fix_ub;
2837 vars[c_ind]->new_lb = col_fix_lb;
2838 vars[c_ind]->lb = col_fix_lb;
2839
2840 /*now update row bounds */
2841 col_start = matbeg[c_ind];
2842 col_end = matbeg[c_ind] + len[c_ind];
2843
2844 for(k = col_start; k < col_end; k++){
2845 rec_row_ind = matind[k];
2846 rec_coeff = matval[k];
2847 if(rec_row_ind == r_ind){
2848 if(c_ind > r_matind[r_matbeg[r_ind]])
2849 row_changed[r_ind] = TRUE;
2850 }else{
2851 row_changed[rec_row_ind] = TRUE;
2852 }
2853 if(!bounds_updated && rec_row_ind < r_ind) bounds_updated = TRUE;
2854 if(rec_coeff >= 0.0){
2855 if(row_ub[rec_row_ind] < inf){
2856 row_ub[rec_row_ind] += rec_coeff*(col_fix_ub - old_ub);
2857 }
2858 if(row_lb[rec_row_ind] > -inf){
2859 row_lb[rec_row_ind] += rec_coeff*(col_fix_lb - old_lb);
2860 }
2861 }else{
2862 if(row_ub[rec_row_ind] < inf){
2863 row_ub[rec_row_ind] += rec_coeff*(col_fix_lb - old_lb);
2864 }
2865 if(row_lb[rec_row_ind] > -inf){
2866 row_lb[rec_row_ind] += rec_coeff*(col_fix_ub - old_ub);
2867 }
2868 }
2869 }
2870 }
2871 }
2872 }
2873 }
2874 }
2875 if(!bounds_updated) break;
2876 iter_cnt++;
2877 }
2878
2879 return PREP_MODIFIED;
2880 }
2881
2882
2883
2884
2885
2886