contrib/glpk-cut-log.patch

   1 @@ This patch is taken from https://github.com/timothy-king/glpk-cut-log and
   2 @@ has the following license:
   3 @@
   4 @@ GLPK is free software: you can redistribute it and/or modify it under the
   5 @@ terms of the GNU General Public License as published by the Free Software
   6 @@ Foundation, either version 3 of the License, or (at your option) any later
   7 @@ version.
   8 @@
   9 From cf84e3855d8c5676557daaca434b42b2fc17dc29 Mon Sep 17 00:00:00 2001
  10 From: Tim King <taking@cs.nyu.edu>
  11 Date: Sat, 14 Dec 2013 16:03:35 -0500
  12 Subject: [PATCH] Adding a function for returning the iteration count.
  13
  14 Adding GLP_ICUTADDED callback:
  15 - Adds a new callback location that is called after every cut is added to the pool.
  16 - During this callback users can request a copy of this cut and query what kind of cut it is.
  17 - GMI cuts also support returning what row in the tableau it came from. Users can use this to replay how the cut was generated.
  18
  19 Logging MIR Cuts
  20 - Changed the IOSAUX to be defined as a linear sum of input rows.
  21   This is visible from the interface
  22 - Logging the multipliers used to create the aggregate row in MIR.
  23 - Logging the delta selected by final successful cmir_sep call.
  24 - Logging the cset selected by final successful cmir_sep call.
  25 - The delta and cset are not yet user visible.
  26
  27 The extra mir logging information is now user visible.
  28
  29 Adding GLP_ICUTSELECT callback:
  30 - ios_process_cuts marks the cuts as selected after turning the cut into a row.
  31 - This callback happens after ios_process_cuts marks and before the pool is
  32   cleared by ios_clear_pool.
  33
  34 Cleaning up the git directory to not track autogenerated files.
  35
  36 Branch logging
  37 - Adds to the integer interface a callback function for when branches are made.
  38   This is called after branch_on.
  39 - Added a public function glp_ios_branch_log.
  40   This returns the structural variable for the branch as well as the value
  41   branched on, and the ids of the nodes for the down and up branches.
  42   If both branches are infeasible, the ids are -1 for both dn and up.
  43   If one branch is infeasible (and no branch is done), this returns -1 for
  44   the infeasible branch and the current node id for the up branch.
  45 - Each node now has a unique ordinal associated with it.
  46 - Added a public function to convert the id of an active node to the unique
  47   node ordinal.
  48 - Added a callback for when a node is closed due to being linearly infeasible.
  49
  50 Improved the mir cut logging facilities to now include:
  51 - the rid id used for virtual lower/upper bound selection,
  52 - and the subst map,
  53 - Changed the size of the cset array to be the same as the new arrays (m+n)
  54 - Fixed a bug in returning the cset selected.
  55
  56 Fixing a bug selecting the correct branch.
  57
  58 Added a callback for tracking row deletion.
  59
  60 Adding notes to node freeze and revive functions.
  61
  62 Commenting out a few debugging printfs for gomory cuts.
  63
  64 Changing an overly aggressive assert in MIR generate when an assignment is outside of the bounds to skipping the application of the mir rule.
  65
  66 Removing some printfs.
  67
  68 Weakening assert to a failure.
  69
  70 Adding a stability failure limit to simplex.
  71
  72 Updating the installation instructions.
  73
  74 ---
  75 diff --git a/configure.ac b/configure.ac
  76 index 0141a8a..ffb9ae4 100644
  77 --- a/configure.ac
  78 +++ b/configure.ac
  79 @@ -1,7 +1,7 @@
  80  dnl Process this file with autoconf to produce a configure script
  81
  82  AC_INIT([GLPK], [4.52], [bug-glpk@gnu.org])
  83 -
  84 +AM_INIT_AUTOMAKE([subdir-objects])
  85  AC_CONFIG_SRCDIR([src/glpk.h])
  86
  87  AC_CONFIG_MACRO_DIR([m4])
  88 diff --git a/src/Makefile.am b/src/Makefile.am
  89 index b39efa6..2a025bc 100644
  90 --- a/src/Makefile.am
  91 +++ b/src/Makefile.am
  92 @@ -20,7 +20,10 @@ libglpk_la_LDFLAGS = \
  93  -version-info 36:0:1 \
  94  -export-symbols-regex '^glp_*'
  95
  96 +
  97 +# all of the cut log sources are listed first
  98  libglpk_la_SOURCES = \
  99 +cutlog01.c \
 100  glpapi01.c \
 101  glpapi02.c \
 102  glpapi03.c \
 103 diff --git a/src/cutlog01.c b/src/cutlog01.c
 104 new file mode 100644
 105 index 0000000..9a85bb7
 106 --- /dev/null
 107 +++ b/src/cutlog01.c
 108 @@ -0,0 +1,452 @@
 109 +/* cutlog01.c (api extension routines) */
 110 +
 111 +
 112 +#include "glpapi.h"
 113 +#include "glpios.h"
 114 +
 115 +int glp_get_it_cnt(glp_prob *P){
 116 +  if(P == NULL){
 117 +    return 0;
 118 +  }else{
 119 +    return P->it_cnt;
 120 +  }
 121 +}
 122 +
 123 +int glp_ios_get_cut(glp_tree *T, int i, int* ind, double* val, int* klass, int* type, double* rhs){
 124 +  xassert(T != NULL);
 125 +
 126 +  IOSCUT* cut;
 127 +  int len;
 128 +  IOSAIJ* aij;
 129 +  glp_prob* prob;
 130 +
 131 +  if (T->reason != GLP_ICUTADDED){
 132 +    xerror("glp_ios_get_cut: not called during cut added.\n");
 133 +  }
 134 +  cut = ios_find_row(T->local, i);
 135 +  if ( cut == NULL ) {
 136 +    xerror("glp_ios_get_cut: called with an invalid index.");
 137 +  }
 138 +  len = 0;
 139 +  for(len = 0, aij = cut->ptr; aij != NULL; aij = aij->next)
 140 +  {
 141 +    len++;
 142 +    if(ind != NULL){ ind[len] = aij->j; }
 143 +    if(val != NULL){ val[len] = aij->val; }
 144 +  }
 145 +  if(klass != NULL){ *klass = cut->klass; }
 146 +  if(type != NULL){ *type = cut->type; }
 147 +  if(rhs != NULL){ *rhs = cut->rhs; }
 148 +  return len;
 149 +}
 150 +
 151 +
 152 +IOSAUX *ios_create_aux(int n, const int rows[], const double coeffs[]){
 153 +  IOSAUX *aux;
 154 +  int i;
 155 +  aux = xmalloc(sizeof(IOSAUX));
 156 +  aux->nrows = n;
 157 +  aux->rows = xcalloc(1+n, sizeof(int));
 158 +  aux->mult = xcalloc(1+n, sizeof(double));
 159 +  aux->selected = -1;
 160 +  aux->mir_cset = NULL;
 161 +
 162 +  for ( i = 1; i <= n; i++)
 163 +  {  aux->rows[i] = rows[i];
 164 +     aux->mult[i] = coeffs[i];
 165 +  }
 166 +
 167 +  return aux;
 168 +}
 169 +
 170 +void ios_delete_aux(IOSAUX *aux){
 171 +  xassert(aux != NULL);
 172 +  xfree(aux->rows);
 173 +  xfree(aux->mult);
 174 +  if( aux->mir_cset != NULL ){
 175 +    xfree(aux->mir_cset);
 176 +     xfree(aux->mir_subst);
 177 +     xfree(aux->mir_vlb_rows);
 178 +     xfree(aux->mir_vub_rows);
 179 +  }
 180 +  xfree(aux);
 181 +  return;
 182 +}
 183 +
 184 +static void cut_set_aux(IOSCUT *cut, int n,
 185 +                        const int rows[], const double coeffs[]){
 186 +  int i;
 187 +  xassert( cut != NULL );
 188 +  if( cut->aux != NULL ) {
 189 +    ios_delete_aux(cut-> aux);
 190 +  }
 191 +
 192 +  cut->aux = ios_create_aux(n, rows, coeffs);
 193 +  xassert( cut->aux->nrows == n );
 194 +}
 195 +
 196 +static void cut_set_aux_mir(IOSAUX *aux, double delta, int m, int n,
 197 +                            const char cset[], const char subst[],
 198 +                            const int vlbrs[], const int vubrs[]){
 199 +  int i;
 200 +  xassert( aux != NULL );
 201 +  if ( aux->mir_cset != NULL )
 202 +  {  xfree(aux->mir_cset);
 203 +     xfree(aux->mir_subst);
 204 +     xfree(aux->mir_vlb_rows);
 205 +     xfree(aux->mir_vub_rows);
 206 +  }
 207 +
 208 +  aux->mir_cset     = xcalloc(1+n+m, sizeof(char));
 209 +  aux->mir_subst    = xcalloc(1+n+m, sizeof(char));
 210 +  aux->mir_vlb_rows = xcalloc(1+n+m, sizeof(int));
 211 +  aux->mir_vub_rows = xcalloc(1+n+m, sizeof(int));
 212 +
 213 +  for ( i = 1; i <= n+m; i++)
 214 +  {  aux->mir_cset[i] = cset[i];
 215 +     aux->mir_subst[i] = subst[i];
 216 +     aux->mir_vlb_rows[i] = vlbrs[i];
 217 +     aux->mir_vub_rows[i] = vubrs[i];
 218 +  }
 219 +
 220 +  aux->mir_delta = delta;
 221 +}
 222 +
 223 +void ios_cut_set_aux_mir(glp_tree *T, int ord, double delta,
 224 +                         const char cset[], const char subst[],
 225 +                         const int vlbrs[], const int vubrs[]){
 226 +  int m, n;
 227 +  IOSCUT *cut;
 228 +  glp_prob *mip;
 229 +  mip = T->mip;
 230 +  m = mip->m;
 231 +  n = mip->n;
 232 +  cut = ios_find_row(T->local, ord);
 233 +  xassert(cut != NULL);
 234 +  cut_set_aux_mir(cut->aux, delta, m, n, cset, subst, vlbrs, vubrs);
 235 +}
 236 +
 237 +void ios_cut_set_single_aux(glp_tree *T, int ord, int j){
 238 +  IOSCUT *cut;
 239 +  cut = ios_find_row(T->local, ord);
 240 +  xassert(cut != NULL);
 241 +
 242 +  /* set up arrays */
 243 +  int ind[1+1];
 244 +  double coeffs[1+1];
 245 +  ind[1] = j;
 246 +  coeffs[1] = +1.0;
 247 +
 248 +  /* call general procedure */
 249 +  cut_set_aux(cut, 1, ind, coeffs);
 250 +}
 251 +
 252 +void ios_cut_set_aux(glp_tree *T, int ord, int n,
 253 +                     const int rows[], const double coeffs[]){
 254 +  IOSCUT *cut;
 255 +  cut = ios_find_row(T->local, ord);
 256 +  xassert(cut != NULL);
 257 +  cut_set_aux(cut, n, rows, coeffs);
 258 +}
 259 +
 260 +int glp_ios_cut_get_aux_nrows(glp_tree *tree, int ord){
 261 +  IOSCUT *cut;
 262 +  IOSAUX *aux;
 263 +  if (tree->reason != GLP_ICUTADDED){
 264 +    xerror("glp_ios_cut_get_aux_nrows: not called during cut added.\n");
 265 +  }
 266 +  cut = ios_find_row(tree->local, ord);
 267 +  if ( cut == NULL ){
 268 +    xerror("glp_ios_cut_get_aux_nrows: not called on a valid cut.\n");
 269 +  }
 270 +  aux = cut->aux;
 271 +  return (aux == NULL) ? 0 : aux->nrows;
 272 +}
 273 +
 274 +void glp_ios_cut_get_aux_rows(glp_tree *tree, int ord,
 275 +                              int rows[], double coeffs[]){
 276 +  IOSCUT *cut;
 277 +  IOSAUX *aux;
 278 +  int j, nrows;
 279 +  if (tree->reason != GLP_ICUTADDED){
 280 +    xerror("glp_ios_cut_get_aux_rows: not called during cut added.\n");
 281 +  }
 282 +  cut = ios_find_row(tree->local, ord);
 283 +  if ( cut == NULL ){
 284 +    xerror("glp_ios_cut_get_aux_rows: not called on a valid cut.\n");
 285 +  }
 286 +  aux = cut->aux;
 287 +  if( aux != NULL ){
 288 +    nrows = aux->nrows;
 289 +    for ( j = 1; j <= nrows; j++ )
 290 +    {  if ( rows != NULL ) { rows[j] = aux->rows[j]; }
 291 +       if ( coeffs != NULL ) { coeffs[j] = aux->mult[j]; }
 292 +    }
 293 +  }
 294 +  return;
 295 +}
 296 +
 297 +
 298 +void glp_ios_cut_get_mir_subst(glp_tree *tree, int ord, char subst[]);
 299 +/* gets mir cut substition information. */
 300 +void glp_ios_cut_get_mir_virtual_rows(glp_tree *tree, int ord,
 301 +                                      int vlb[], int vub[]);
 302 +/* gets mir cut virtual bounds rows. */
 303 +
 304 +void glp_ios_cut_get_mir_cset(glp_tree *tree, int ord, char *cset){
 305 +  glp_prob *mip;
 306 +  IOSCUT *cut;
 307 +  IOSAUX *aux;
 308 +  int j, n, m;
 309 +  if ( tree->reason != GLP_ICUTADDED ){
 310 +    xerror("glp_ios_cut_get_aux_mir: not called during cut added.\n");
 311 +  }
 312 +  cut = ios_find_row(tree->local, ord);
 313 +  if ( cut == NULL ){
 314 +    xerror("glp_ios_cut_get_aux_mir: not called on a cut.\n");
 315 +  }
 316 +  if ( cut->klass != GLP_RF_MIR ){
 317 +    xerror("glp_ios_cut_get_aux_mir: not called on a mir cut.\n");
 318 +  }
 319 +  aux = cut->aux;
 320 +  mip = tree->mip;
 321 +  m = mip->m;
 322 +  n = mip->n;
 323 +
 324 +  if( cset != NULL ){
 325 +    for ( j=1; j <= n+m; j++ ){
 326 +      cset[j] = (aux->mir_cset == NULL) ? 0 : aux->mir_cset[j];
 327 +    }
 328 +  }
 329 +}
 330 +void glp_ios_cut_get_mir_subst(glp_tree *tree, int ord, char *subst){
 331 +  glp_prob *mip;
 332 +  IOSCUT *cut;
 333 +  IOSAUX *aux;
 334 +  int j, n, m;
 335 +  if ( tree->reason != GLP_ICUTADDED ){
 336 +    xerror("glp_ios_cut_get_aux_mir: not called during cut added.\n");
 337 +  }
 338 +  cut = ios_find_row(tree->local, ord);
 339 +  if ( cut == NULL ){
 340 +    xerror("glp_ios_cut_get_aux_mir: not called on a cut.\n");
 341 +  }
 342 +  if ( cut->klass != GLP_RF_MIR ){
 343 +    xerror("glp_ios_cut_get_aux_mir: not called on a mir cut.\n");
 344 +  }
 345 +  aux = cut->aux;
 346 +  mip = tree->mip;
 347 +  m = mip->m;
 348 +  n = mip->n;
 349 +
 350 +  if( subst != NULL ){
 351 +    for ( j=1; j <= n+m; j++ ){
 352 +      subst[j] = (aux->mir_subst == NULL) ? 0 : aux->mir_subst[j];
 353 +    }
 354 +  }
 355 +}
 356 +void glp_ios_cut_get_mir_virtual_rows(glp_tree *tree, int ord, int vlb_rows[], int vub_rows[]){
 357 +  glp_prob *mip;
 358 +  IOSCUT *cut;
 359 +  IOSAUX *aux;
 360 +  int j, n, m;
 361 +  if ( tree->reason != GLP_ICUTADDED ){
 362 +    xerror("glp_ios_cut_get_aux_mir: not called during cut added.\n");
 363 +  }
 364 +  cut = ios_find_row(tree->local, ord);
 365 +  if ( cut == NULL ){
 366 +    xerror("glp_ios_cut_get_aux_mir: not called on a cut.\n");
 367 +  }
 368 +  if ( cut->klass != GLP_RF_MIR ){
 369 +    xerror("glp_ios_cut_get_aux_mir: not called on a mir cut.\n");
 370 +  }
 371 +  aux = cut->aux;
 372 +  mip = tree->mip;
 373 +  m = mip->m;
 374 +  n = mip->n;
 375 +
 376 +  for ( j=1; j <= n+m; j++ ){
 377 +    vlb_rows[j] = (aux->mir_vlb_rows == NULL) ? 0 : aux->mir_vlb_rows[j];
 378 +    vub_rows[j] = (aux->mir_vub_rows == NULL) ? 0 : aux->mir_vub_rows[j];
 379 +  }
 380 +}
 381 +double glp_ios_cut_get_mir_delta(glp_tree *tree, int ord){
 382 +  glp_prob *mip;
 383 +  IOSCUT *cut;
 384 +  IOSAUX *aux;
 385 +  int j, n, m;
 386 +  if ( tree->reason != GLP_ICUTADDED ){
 387 +    xerror("glp_ios_cut_get_aux_mir: not called during cut added.\n");
 388 +  }
 389 +  cut = ios_find_row(tree->local, ord);
 390 +  if ( cut == NULL ){
 391 +    xerror("glp_ios_cut_get_aux_mir: not called on a cut.\n");
 392 +  }
 393 +  if ( cut->klass != GLP_RF_MIR ){
 394 +    xerror("glp_ios_cut_get_aux_mir: not called on a mir cut.\n");
 395 +  }
 396 +  aux = cut->aux;
 397 +  return aux->mir_delta;
 398 +}
 399 +
 400 +
 401 +void ios_cut_set_selected(IOSCUT *cut, int sel){
 402 +#ifdef CUT_DEBUG
 403 +  static int i = 0;
 404 +  ++i;
 405 +  printf("ios_cut_set_selected: %d %d %p\n", i, sel, cut);
 406 +#endif
 407 +
 408 +  IOSAUX *aux;
 409 +  aux = cut->aux;
 410 +  if ( aux != NULL ){
 411 +    aux->selected = sel;
 412 +  }
 413 +}
 414 +
 415 +int glp_ios_selected_cuts(glp_tree *tree, int ords[], int sel[]){
 416 +  int len, j, N, s;
 417 +  IOSPOOL* pool;
 418 +  IOSCUT* cut;
 419 +  IOSAUX* aux;
 420 +  if ( tree == NULL ){
 421 +    xerror("glp_ios_selected_cuts: not called with a valid tree.\n");
 422 +  }
 423 +  if ( tree->reason != GLP_ICUTSELECT ){
 424 +    xerror("glp_ios_selected_cuts: not called during cut selected.\n");
 425 +  }
 426 +
 427 +  pool = tree->local;
 428 +  if ( pool == NULL ){
 429 +    xerror("glp_ios_selected_cuts: called on a malformed tree.\n");
 430 +  }
 431 +
 432 +  for (len = 0, j = 1, cut = pool->head; cut != NULL; cut = cut->next, j++)
 433 +  {  aux = cut->aux;
 434 +#ifdef CUT_DEBUG
 435 +     printf("glp_ios_selected_cuts: %d %p\n", j, cut);
 436 +#endif
 437 +     if ( aux != NULL )
 438 +     { s = aux->selected;
 439 +       if ( s >= 0 )
 440 +       {  len++;
 441 +          if (ords != NULL) { ords[len] = j; }
 442 +          if (sel != NULL)  { sel[len] = s; }
 443 +       }
 444 +     }
 445 +  }
 446 +  return len;
 447 +}
 448 +
 449 +int glp_ios_branch_log(glp_tree *tree, double *br_val, int* parent, int* dn, int* up){
 450 +  IOSNPD *node;
 451 +  int br_result, br_var;
 452 +  int p, d, u;
 453 +  double v;
 454 +  glp_prob *mip;
 455 +  if ( tree == NULL ){
 456 +    xerror("glp_ios_branch_log: not called with a valid tree \n");
 457 +  }
 458 +  if ( tree->reason != GLP_LI_BRANCH ){
 459 +    xerror("glp_ios_branch_log: not called during GLP_LI_BRANCH \n");
 460 +  }
 461 +  mip = tree->mip;
 462 +  if ( mip == NULL ){
 463 +    xerror("glp_ios_branch_log: not called with a valid tree\n");
 464 +  }
 465 +  br_result = tree->br_result;
 466 +  br_var = tree->br_var;
 467 +  switch(br_result){
 468 +  case 0:
 469 +    p = tree->br_node;
 470 +    node = tree->slot[p].node;
 471 +    break;
 472 +  case 1:
 473 +  case 2:
 474 +    node = tree->curr;
 475 +    p = node->p;
 476 +    break;
 477 +  default:
 478 +    xerror("glp_ios_branch_log: br_result is not properly set \n");
 479 +  }
 480 +  if( node == NULL ){
 481 +    xerror("glp_ios_branch_log: not called with a valid tree \n");
 482 +  }
 483 +  switch(br_result){
 484 +  case 0:
 485 +    v = node->br_val;
 486 +    d = tree->dn_child;
 487 +    u = tree->up_child;
 488 +    break;
 489 +  case 1:
 490 +    v = mip->col[br_var]->prim;
 491 +    if(tree->br_to_up){
 492 +      d = -1;
 493 +      u = p;
 494 +    }else{
 495 +      d = p;
 496 +      u = -1;
 497 +    }
 498 +    break;
 499 +  case 2:
 500 +    v = mip->col[br_var]->prim;
 501 +    d = -1;
 502 +    u = -1;
 503 +    break;
 504 +  default:
 505 +    xerror("glp_ios_branch_log: not called with a valid tree \n");
 506 +  }
 507 +
 508 +  if(br_val != NULL){ *br_val = v; }
 509 +  if(parent != NULL){ *parent = p; }
 510 +  if(dn != NULL){ *dn = d; }
 511 +  if(up != NULL){ *up = u; }
 512 +
 513 +  return br_var;
 514 +}
 515 +
 516 +int glp_ios_node_ord(glp_tree *tree, int p){
 517 +  IOSNPD *node;
 518 +  if ( tree == NULL ){
 519 +    xerror("glp_ios_node_ord: not called with a valid tree.\n");
 520 +  }
 521 +  if (!(1 <= p && p <= tree->nslots)){
 522 +    xerror("glp_ios_node_ord: not called with a valid p.\n");
 523 +  }
 524 +  node = tree->slot[p].node;
 525 +  return node->ord;
 526 +}
 527 +
 528 +void ios_cb_rows_deleted(glp_tree *T, int n, const int* rows){
 529 +  if (T->parm->cb_func != NULL)
 530 +  {
 531 +    xassert(T->reason == 0);
 532 +    xassert(T->deleting_rows == NULL);
 533 +    xassert(T->num_deleting_rows == 0);
 534 +    T->num_deleting_rows = n;
 535 +    T->deleting_rows = rows;
 536 +
 537 +    T->reason = GLP_LI_DELROW;
 538 +    T->parm->cb_func(T, T->parm->cb_info);
 539 +    T->reason = 0;
 540 +    T->num_deleting_rows = 0;
 541 +    T->deleting_rows = NULL;
 542 +  }
 543 +}
 544 +
 545 +int glp_ios_rows_deleted(glp_tree *tree, int* rows){
 546 +  if ( tree == NULL ){
 547 +    xerror("glp_ios_rows_deleted: not called with a valid tree.\n");
 548 +  }
 549 +  if ( tree->reason != GLP_LI_DELROW ){
 550 +    xerror("glp_ios_rows_deleted: not called with a valid reason.\n");
 551 +  }
 552 +
 553 +  int j;
 554 +  if(rows != NULL){
 555 +    for(j=1; j <= tree->num_deleting_rows; j++){
 556 +      rows[j] = tree->deleting_rows[j];
 557 +    }
 558 +  }
 559 +  return tree->num_deleting_rows;
 560 +}
 561 diff --git a/src/glpapi06.c b/src/glpapi06.c
 562 index 743d6be..05d0498 100644
 563 --- a/src/glpapi06.c
 564 +++ b/src/glpapi06.c
 565 @@ -488,6 +488,7 @@ void glp_init_smcp(glp_smcp *parm)
 566        parm->out_frq = 500;
 567        parm->out_dly = 0;
 568        parm->presolve = GLP_OFF;
 569 +      parm->stability_lmt = 200;
 570        return;
 571  }
 572
 573 diff --git a/src/glpapi13.c b/src/glpapi13.c
 574 index 926dda4..55adf44 100644
 575 --- a/src/glpapi13.c
 576 +++ b/src/glpapi13.c
 577 @@ -453,8 +453,14 @@ void glp_ios_row_attr(glp_tree *tree, int i, glp_attr *attr)
 578
 579  int glp_ios_pool_size(glp_tree *tree)
 580  {     /* determine current size of the cut pool */
 581 -      if (tree->reason != GLP_ICUTGEN)
 582 -         xerror("glp_ios_pool_size: operation not allowed\n");
 583 +      switch(tree->reason)
 584 +      {  case GLP_ICUTGEN:
 585 +         case GLP_ICUTADDED:
 586 +         case GLP_ICUTSELECT:
 587 +           break;
 588 +         default:
 589 +           xerror("glp_ios_pool_size: operation not allowed\n");
 590 +      }
 591        xassert(tree->local != NULL);
 592        return tree->local->size;
 593  }
 594 diff --git a/src/glpios.h b/src/glpios.h
 595 index 9e2d6b2..3b31901 100644
 596 --- a/src/glpios.h
 597 +++ b/src/glpios.h
 598 @@ -36,6 +36,9 @@ typedef struct IOSAIJ IOSAIJ;
 599  typedef struct IOSPOOL IOSPOOL;
 600  typedef struct IOSCUT IOSCUT;
 601
 602 +
 603 +typedef struct IOSAUX IOSAUX;
 604 +
 605  struct glp_tree
 606  {     /* branch-and-bound tree */
 607        int magic;
 608 @@ -217,6 +220,19 @@ struct glp_tree
 609           GLP_NO_BRNCH - use general selection technique */
 610        int child;
 611        /* subproblem reference number corresponding to br_sel */
 612 +
 613 +      /* start of cut log extras */
 614 +      int br_result;
 615 +      int br_to_up;
 616 +      int br_node;
 617 +      /* subproblem reference number for the just branched node.*/
 618 +      int dn_child;
 619 +      /* subproblem reference number for the just created down node */
 620 +      int up_child;
 621 +      /* subproblem reference number for the just created up node */
 622 +
 623 +      const int* deleting_rows;
 624 +      int num_deleting_rows;
 625  };
 626
 627  struct IOSLOT
 628 @@ -229,6 +245,8 @@ struct IOSLOT
 629
 630  struct IOSNPD
 631  {     /* node subproblem descriptor */
 632 +      int ord;
 633 +      /* this is a unique ordinal for each subproblem  */
 634        int p;
 635        /* subproblem reference number (it is the index to corresponding
 636           slot, i.e. slot[p] points to this descriptor) */
 637 @@ -393,12 +411,51 @@ struct IOSCUT
 638           GLP_FX: sum a[j] * x[j]  = b */
 639        double rhs;
 640        /* cut right-hand side */
 641 +
 642 +      IOSAUX *aux;
 643 +      /* cut auxillary source information */
 644 +
 645        IOSCUT *prev;
 646        /* pointer to previous cut */
 647        IOSCUT *next;
 648        /* pointer to next cut */
 649  };
 650
 651 +struct IOSAUX
 652 +{
 653 +      /* aux (auxillary source information for each cut)
 654 +       * Each cut operates on a sum of rows
 655 +       * Each cut operates on a row that can be described using
 656 +       * a current row or a previous cut.
 657 +       * To generalize, we assume each row is a sum of two rows:
 658 +       *   row[r]* r_mult + pool[c] * c_mult
 659 +       */
 660 +      int nrows;
 661 +      int *rows;
 662 +      double *mult;
 663 +
 664 +      int selected;
 665 +      /* when < 0 this has not yet been turned into a row
 666 +         when >=0 this is the id of the row added. */
 667 +
 668 +      double mir_delta;
 669 +      /* the delta selected by a mir cut */
 670 +
 671 +      char *mir_cset;
 672 +      /* complimented set */
 673 +      /* if this is NULL, it is implicitly all 0s */
 674 +
 675 +      char *mir_subst;
 676 +      /* the substition vectors */
 677 +
 678 +      int *mir_vlb_rows;
 679 +      /* the substition vectors */
 680 +
 681 +      int *mir_vub_rows;
 682 +      /* the substition vectors */
 683 +};
 684 +
 685 +
 686  #define ios_create_tree _glp_ios_create_tree
 687  glp_tree *ios_create_tree(glp_prob *mip, const glp_iocp *parm);
 688  /* create branch-and-bound tree */
 689 @@ -615,6 +672,31 @@ int ios_choose_node(glp_tree *T);
 690  int ios_choose_var(glp_tree *T, int *next);
 691  /* select variable to branch on */
 692
 693 +/* functions added to retrieve information */
 694 +IOSAUX *ios_create_aux(int n, const int rows[], const double coeffs[]);
 695 +/* creates an aux with n rows */
 696 +
 697 +void ios_delete_aux(IOSAUX *aux);
 698 +/* deletes an aux */
 699 +
 700 +void ios_cut_set_single_aux(glp_tree *T, int ord, int j);
 701 +/* sets the aux of row ord to be a single row j */
 702 +
 703 +
 704 +void ios_cut_set_aux(glp_tree *T, int ord, int n,
 705 +                     const int rows[], const double coeffs[]);
 706 +/* sets an arbitrary aux sum */
 707 +
 708 +void ios_cut_set_aux_mir(glp_tree *T, int ord, double delta,
 709 +                         const char cset[], const char subst[],
 710 +                         const int vlbrs[], const int vubrs[]);
 711 +/* sets the extra mir information */
 712 +
 713 +void ios_cut_set_selected(IOSCUT *cut, int i);
 714 +/* the cut has been added as row i */
 715 +
 716 +void ios_cb_rows_deleted(glp_tree *T, int n, const int* rows);
 717 +
 718  #endif
 719
 720  /* eof */
 721 diff --git a/src/glpios01.c b/src/glpios01.c
 722 index 70798fd..d9e5703 100644
 723 --- a/src/glpios01.c
 724 +++ b/src/glpios01.c
 725 @@ -159,6 +159,16 @@ glp_tree *ios_create_tree(glp_prob *mip, const glp_iocp *parm)
 726        tree->stop = 0;
 727        /* create the root subproblem, which initially is identical to
 728           the original MIP */
 729 +
 730 +      tree->br_result = 0;
 731 +      tree->br_to_up = 0;
 732 +      tree->br_node = 0;
 733 +      tree->dn_child = 0;
 734 +      tree->up_child = 0;
 735 +
 736 +      tree->deleting_rows = NULL;
 737 +      tree->num_deleting_rows = 0;
 738 +
 739        new_node(tree, NULL);
 740        return tree;
 741  }
 742 @@ -278,6 +288,7 @@ void ios_revive_node(glp_tree *tree, int p)
 743              double *val;
 744              ind = xcalloc(1+n, sizeof(int));
 745              val = xcalloc(1+n, sizeof(double));
 746 +            /* maintains the row order during revival */
 747              for (r = node->r_ptr; r != NULL; r = r->next)
 748              {  i = glp_add_rows(mip, 1);
 749                 glp_set_row_name(mip, i, r->name);
 750 @@ -457,6 +468,13 @@ void ios_freeze_node(glp_tree *tree)
 751              double *val;
 752              ind = xcalloc(1+n, sizeof(int));
 753              val = xcalloc(1+n, sizeof(double));
 754 +            /* Added rows are stored in the same order!
 755 +             * This is done by 2 reversals.
 756 +             * - Going from i = m down pred_m.
 757 +             * - Storing the list by a stack "push"
 758 +             *   node->r_ptr = r;
 759 +             *   r->next = node->r_ptr;
 760 +             */
 761              for (i = m; i > pred_m; i--)
 762              {  GLPROW *row = mip->row[i];
 763                 IOSROW *r;
 764 @@ -501,6 +519,8 @@ void ios_freeze_node(glp_tree *tree)
 765              xassert(nrs > 0);
 766              num = xcalloc(1+nrs, sizeof(int));
 767              for (i = 1; i <= nrs; i++) num[i] = root_m + i;
 768 +            /* To not call ios_cb_rows_deleted here.
 769 +               These rows have been saved earlier. */
 770              glp_del_rows(mip, nrs, num);
 771              xfree(num);
 772           }
 773 @@ -636,6 +656,9 @@ static IOSNPD *new_node(glp_tree *tree, IOSNPD *parent)
 774        tree->a_cnt++;
 775        tree->n_cnt++;
 776        tree->t_cnt++;
 777 +
 778 +      node->ord = tree->t_cnt;
 779 +
 780        /* increase the number of child subproblems */
 781        if (parent == NULL)
 782           xassert(p == 1);
 783 @@ -818,6 +841,8 @@ void ios_delete_tree(glp_tree *tree)
 784           xassert(nrs > 0);
 785           num = xcalloc(1+nrs, sizeof(int));
 786           for (i = 1; i <= nrs; i++) num[i] = tree->orig_m + i;
 787 +         /* Do not call ios_cb_rows_deleted here.
 788 +            This does not help log information. */
 789           glp_del_rows(mip, nrs, num);
 790           xfree(num);
 791        }
 792 @@ -1417,6 +1442,7 @@ int ios_add_row(glp_tree *tree, IOSPOOL *pool,
 793        cut->rhs = rhs;
 794        cut->prev = pool->tail;
 795        cut->next = NULL;
 796 +      cut->aux = NULL;
 797        if (cut->prev == NULL)
 798           pool->head = cut;
 799        else
 800 @@ -1517,6 +1543,11 @@ void ios_del_row(glp_tree *tree, IOSPOOL *pool, int i)
 801           cut->ptr = aij->next;
 802           dmp_free_atom(tree->pool, aij, sizeof(IOSAIJ));
 803        }
 804 +
 805 +      if (cut->aux != NULL){
 806 +        ios_delete_aux(cut->aux);
 807 +      }
 808 +
 809        dmp_free_atom(tree->pool, cut, sizeof(IOSCUT));
 810        pool->size--;
 811        return;
 812 @@ -1535,6 +1566,10 @@ void ios_clear_pool(glp_tree *tree, IOSPOOL *pool)
 813              cut->ptr = aij->next;
 814              dmp_free_atom(tree->pool, aij, sizeof(IOSAIJ));
 815           }
 816 +
 817 +         if (cut->aux != NULL){
 818 +           ios_delete_aux(cut->aux);
 819 +         }
 820           dmp_free_atom(tree->pool, cut, sizeof(IOSCUT));
 821        }
 822        pool->size = 0;
 823 diff --git a/src/glpios03.c b/src/glpios03.c
 824 index 80d701b..03e9208 100644
 825 --- a/src/glpios03.c
 826 +++ b/src/glpios03.c
 827 @@ -358,12 +358,12 @@ static void fix_by_red_cost(glp_tree *T)
 828  *  2 - both branches are hopeless and have been pruned; new subproblem
 829  *      selection is needed to continue the search. */
 830
 831 -static int branch_on(glp_tree *T, int j, int next)
 832 +static int branch_on(glp_tree *T, int j, int next, int clone[], int* to_up)
 833  {     glp_prob *mip = T->mip;
 834        IOSNPD *node;
 835        int m = mip->m;
 836        int n = mip->n;
 837 -      int type, dn_type, up_type, dn_bad, up_bad, p, ret, clone[1+2];
 838 +      int type, dn_type, up_type, dn_bad, up_bad, p, ret;
 839        double lb, ub, beta, new_ub, new_lb, dn_lp, up_lp, dn_bnd, up_bnd;
 840        /* determine bounds and value of x[j] in optimal solution to LP
 841           relaxation of the current subproblem */
 842 @@ -431,6 +431,7 @@ static int branch_on(glp_tree *T, int j, int next)
 843           else
 844              xassert(mip != mip);
 845           ret = 1;
 846 +         *to_up = 0; /* up is bad. Do not go to up. */
 847           goto done;
 848        }
 849        else if (dn_bad)
 850 @@ -449,6 +450,7 @@ static int branch_on(glp_tree *T, int j, int next)
 851           else
 852              xassert(mip != mip);
 853           ret = 1;
 854 +         *to_up = 1; /* down is bad. Go to up. */
 855           goto done;
 856        }
 857        /* both down- and up-branches seem to be hopeful */
 858 @@ -702,7 +704,8 @@ static void remove_cuts(glp_tree *T)
 859           }
 860        }
 861        if (cnt > 0)
 862 -      {  glp_del_rows(T->mip, cnt, num);
 863 +      {  ios_cb_rows_deleted(T, cnt, num);
 864 +         glp_del_rows(T->mip, cnt, num);
 865  #if 0
 866           xprintf("%d inactive cut(s) removed\n", cnt);
 867  #endif
 868 @@ -784,6 +787,7 @@ static void display_cut_info(glp_tree *T)
 869
 870  int ios_driver(glp_tree *T)
 871  {     int p, curr_p, p_stat, d_stat, ret;
 872 +      int branch_clones[1 + 2];
 873  #if 1 /* carry out to glp_tree */
 874        int pred_p = 0;
 875        /* if the current subproblem has been just created due to
 876 @@ -1010,6 +1014,12 @@ more: /* minor loop starts here */
 877           if (T->parm->msg_lev >= GLP_MSG_DBG)
 878              xprintf("LP relaxation has no feasible solution\n");
 879           /* prune the branch */
 880 +         if (T->parm->cb_func != NULL)
 881 +         {  xassert(T->reason == 0);
 882 +            T->reason = GLP_LI_CLOSE;
 883 +            T->parm->cb_func(T, T->parm->cb_info);
 884 +            T->reason = 0;
 885 +         }
 886           goto fath;
 887        }
 888        else
 889 @@ -1219,6 +1229,14 @@ more: /* minor loop starts here */
 890           ios_process_cuts(T);
 891           T->reason = 0;
 892        }
 893 +      /* if the local cut pool is not empty and the callback func is there,
 894 +         this gives the callback the chance to see what was selected. */
 895 +      if (T->parm->cb_func != NULL && T->local->size > 0)
 896 +      {  xassert(T->reason == 0);
 897 +         T->reason = GLP_ICUTSELECT;
 898 +         T->parm->cb_func(T, T->parm->cb_info);
 899 +         T->reason = 0;
 900 +      }
 901        /* clear the local cut pool */
 902        ios_clear_pool(T, T->local);
 903        /* perform re-optimization, if necessary */
 904 @@ -1255,7 +1273,34 @@ more: /* minor loop starts here */
 905           T->br_var = ios_choose_var(T, &T->br_sel);
 906        /* perform actual branching */
 907        curr_p = T->curr->p;
 908 -      ret = branch_on(T, T->br_var, T->br_sel);
 909 +      xassert(T->br_to_up == 0);
 910 +      ret = branch_on(T, T->br_var, T->br_sel, branch_clones, &T->br_to_up);
 911 +      if (T->parm->cb_func != NULL)
 912 +      {  xassert(T->reason == 0);
 913 +         xassert(T->br_node == 0);
 914 +         xassert(T->dn_child == 0);
 915 +         xassert(T->up_child == 0);
 916 +         // record a branch here
 917 +         T->reason = GLP_LI_BRANCH;
 918 +         // at this point T->br_var is the branching variable
 919 +         T->br_node = curr_p;
 920 +         T->br_result = ret;
 921 +         if(ret == 0){
 922 +           T->dn_child = branch_clones[1];
 923 +           T->up_child = branch_clones[2];
 924 +         }
 925 +         T->parm->cb_func(T, T->parm->cb_info);
 926 +         T->reason = 0;
 927 +         T->br_node = 0;
 928 +         T->dn_child = 0;
 929 +         T->up_child = 0;
 930 +         if (T->stop)
 931 +         {  ret = GLP_ESTOP;
 932 +            goto done;
 933 +         }
 934 +      }
 935 +      T->br_to_up = 0;
 936 +
 937        T->br_var = T->br_sel = 0;
 938        if (ret == 0)
 939        {  /* both branches have been created */
 940 diff --git a/src/glpios05.c b/src/glpios05.c
 941 index b9322b9..caf69d7 100644
 942 --- a/src/glpios05.c
 943 +++ b/src/glpios05.c
 944 @@ -65,6 +65,9 @@ static void gen_cut(glp_tree *tree, struct worka *worka, int j)
 945        double *phi = worka->phi;
 946        int i, k, len, kind, stat;
 947        double lb, ub, alfa, beta, ksi, phi1, rhs;
 948 +      int input_j;
 949 +      input_j = j;
 950 +
 951        /* compute row of the simplex tableau, which (row) corresponds
 952           to specified basic variable xB[i] = x[m+j]; see (23) */
 953        len = glp_eval_tab_row(mip, m+j, ind, val);
 954 @@ -73,6 +76,7 @@ static void gen_cut(glp_tree *tree, struct worka *worka, int j)
 955           if it would be computed with formula (27); it is assumed that
 956           beta[i] is fractional enough */
 957        beta = mip->col[j]->prim;
 958 +
 959        /* compute cut coefficients phi and right-hand side rho, which
 960           correspond to formula (30); dense format is used, because rows
 961           of the simplex tableau is usually dense */
 962 @@ -104,6 +108,7 @@ static void gen_cut(glp_tree *tree, struct worka *worka, int j)
 963           xassert(stat != GLP_BS);
 964           /* determine row coefficient ksi[i,j] at xN[j]; see (23) */
 965           ksi = val[j];
 966 +         /* printf("%d %4.15f ", k, ksi); */
 967           /* if ksi[i,j] is too large in the magnitude, do not generate
 968              the cut */
 969           if (fabs(ksi) > 1e+05) goto fini;
 970 @@ -135,6 +140,7 @@ static void gen_cut(glp_tree *tree, struct worka *worka, int j)
 971                 /* y[j] is integer */
 972                 if (fabs(alfa - floor(alfa + 0.5)) < 1e-10)
 973                 {  /* alfa[i,j] is close to nearest integer; skip it */
 974 +                  /* printf("(skip)"); */
 975                    goto skip;
 976                 }
 977                 else if (f(alfa) <= f(beta))
 978 @@ -170,6 +176,13 @@ static void gen_cut(glp_tree *tree, struct worka *worka, int j)
 979           }
 980  skip:    ;
 981        }
 982 +      /* printf("\n"); */
 983 +      /* for (i = 1; i <= m+n; i++) */
 984 +      /* { */
 985 +      /*   printf("%i %f, ", i, phi[i]); */
 986 +      /* } */
 987 +      /* printf("\n"); */
 988 +
 989        /* now the cut has the form sum_k phi[k] * x[k] >= rho, where cut
 990           coefficients are stored in the array phi in dense format;
 991           x[1,...,m] are auxiliary variables, x[m+1,...,m+n] are struc-
 992 @@ -218,8 +231,20 @@ skip:    ;
 993        ios_add_cut_row(tree, pool, GLP_RF_GMI, len, ind, val, GLP_LO,
 994           rhs);
 995  #else
 996 -      glp_ios_add_row(tree, NULL, GLP_RF_GMI, 0, len, ind, val, GLP_LO,
 997 -         rhs);
 998 +      int ord;
 999 +      ord = glp_ios_add_row(tree, NULL, GLP_RF_GMI, 0, len, ind, val,
1000 +                            GLP_LO, rhs);
1001 +      ios_cut_set_single_aux(tree, ord, input_j);
1002 +
1003 +      /* printf("ord: % d beta %f\n", ord, beta); */
1004 +
1005 +      /** callback for a cut being added to the cut pool */
1006 +      if (tree->parm->cb_func != NULL)
1007 +      {  xassert(tree->reason == GLP_ICUTGEN);
1008 +         tree->reason = GLP_ICUTADDED;
1009 +         tree->parm->cb_func(tree, tree->parm->cb_info);
1010 +         tree->reason = GLP_ICUTGEN;
1011 +      }
1012  #endif
1013  fini: return;
1014  }
1015 diff --git a/src/glpios06.c b/src/glpios06.c
1016 index 2bd0453..4bd3cf5 100644
1017 --- a/src/glpios06.c
1018 +++ b/src/glpios06.c
1019 @@ -100,6 +100,29 @@ struct MIR
1020        /* sparse vector of cutting plane coefficients, alpha[k] */
1021        double cut_rhs;
1022        /* right-hand size of the cutting plane, beta */
1023 +
1024 +      /*-------------------------------------------------------------*/
1025 +      /* Extras I've added to reproduce a cut externally */
1026 +      double cut_delta;
1027 +      /* the delta used for generating the cut */
1028 +      char *cut_cset; /* char cut_vec[1+m+n]; */
1029 +      /* cut_cset[k], 1 <= k <= m+n, is set to true if structural
1030 +         variable x[k] was complemented in the cut:
1031 +         0 - x[k] has been not been complemented
1032 +         non 0 - x[k] has been complemented */
1033 +      int *vlb_rows; /* int vlb_rows[1+m+n]; */
1034 +      /* vlb_rows[k], 1 <= k <= m+n,
1035 +       * vlb_rows[k] <= 0 if virtual lower bound has not been set
1036 +       * vlb_rows[k] = r if virtual lower bound was set using the row r
1037 +       */
1038 +      int *vub_rows; /* int vub_rows[1+m+n]; */
1039 +      /* vub_rows[k], 1 <= k <= m+n,
1040 +       * vub_rows[k] <= 0 if virtual upper bound has not been set
1041 +       * vub_rows[k] = r if virtual upper bound was set using the row r
1042 +       */
1043 +
1044 +      double *agg_coeffs;
1045 +      /* coefficients used to multiply agg_coeffs */
1046  };
1047
1048  /***********************************************************************
1049 @@ -146,6 +169,7 @@ static void set_row_attrib(glp_tree *tree, struct MIR *mir)
1050                 xassert(row != row);
1051           }
1052           mir->vlb[k] = mir->vub[k] = 0;
1053 +         mir->vlb_rows[k] = mir->vub_rows[k] = 0;
1054        }
1055        return;
1056  }
1057 @@ -181,6 +205,7 @@ static void set_col_attrib(glp_tree *tree, struct MIR *mir)
1058                 xassert(col != col);
1059           }
1060           mir->vlb[k] = mir->vub[k] = 0;
1061 +         mir->vlb_rows[k] = mir->vub_rows[k] = 0;
1062        }
1063        return;
1064  }
1065 @@ -230,6 +255,7 @@ static void set_var_bounds(glp_tree *tree, struct MIR *mir)
1066              {  /* set variable lower bound for x1 */
1067                 mir->lb[k1] = - a2 / a1;
1068                 mir->vlb[k1] = k2;
1069 +               mir->vlb_rows[k1]  = i;
1070                 /* the row should not be used */
1071                 mir->skip[i] = 1;
1072              }
1073 @@ -240,6 +266,7 @@ static void set_var_bounds(glp_tree *tree, struct MIR *mir)
1074              {  /* set variable upper bound for x1 */
1075                 mir->ub[k1] = - a2 / a1;
1076                 mir->vub[k1] = k2;
1077 +               mir->vub_rows[k1]  = i;
1078                 /* the row should not be used */
1079                 mir->skip[i] = 1;
1080              }
1081 @@ -313,6 +340,13 @@ void *ios_mir_init(glp_tree *tree)
1082        mir->subst = xcalloc(1+m+n, sizeof(char));
1083        mir->mod_vec = ios_create_vec(m+n);
1084        mir->cut_vec = ios_create_vec(m+n);
1085 +
1086 +      /* added */
1087 +      mir->cut_cset = xcalloc(1+m+n, sizeof(char));
1088 +      mir->vlb_rows = xcalloc(1+m+n, sizeof(int));
1089 +      mir->vub_rows = xcalloc(1+m+n, sizeof(int));
1090 +      mir->agg_coeffs = xcalloc(1+MAXAGGR, sizeof(double));
1091 +
1092        /* set global row attributes */
1093        set_row_attrib(tree, mir);
1094        /* set global column attributes */
1095 @@ -405,6 +439,9 @@ static void initial_agg_row(glp_tree *tree, struct MIR *mir, int i)
1096        mir->skip[i] = 2;
1097        mir->agg_cnt = 1;
1098        mir->agg_row[1] = i;
1099 +
1100 +      mir->agg_coeffs[1] = +1.0;
1101 +
1102        /* use x[i] - sum a[i,j] * x[m+j] = 0, where x[i] is auxiliary
1103           variable of row i, x[m+j] are structural variables */
1104        ios_clear_vec(mir->agg_vec);
1105 @@ -784,19 +821,19 @@ static int cmir_cmp(const void *p1, const void *p2)
1106
1107  static double cmir_sep(const int n, const double a[], const double b,
1108        const double u[], const double x[], const double s,
1109 -      double alpha[], double *beta, double *gamma)
1110 +      double alpha[], double *beta, double *gamma,
1111 +      double* delta, char cset[])
1112  {     int fail, j, k, nv, v;
1113 -      double delta, eps, d_try[1+3], r, r_best;
1114 -      char *cset;
1115 +      double eps, d_try[1+3], r, r_best;
1116        struct vset *vset;
1117        /* allocate working arrays */
1118 -      cset = xcalloc(1+n, sizeof(char));
1119 +      //cset = xcalloc(1+n, sizeof(char));
1120        vset = xcalloc(1+n, sizeof(struct vset));
1121        /* choose initial C */
1122        for (j = 1; j <= n; j++)
1123           cset[j] = (char)(x[j] >= 0.5 * u[j]);
1124        /* choose initial delta */
1125 -      r_best = delta = 0.0;
1126 +      r_best = (*delta) = 0.0;
1127        for (j = 1; j <= n; j++)
1128        {  xassert(a[j] != 0.0);
1129           /* if x[j] is close to its bounds, skip it */
1130 @@ -809,16 +846,16 @@ static double cmir_sep(const int n, const double a[], const double b,
1131           /* compute violation */
1132           r = - (*beta) - (*gamma) * s;
1133           for (k = 1; k <= n; k++) r += alpha[k] * x[k];
1134 -         if (r_best < r) r_best = r, delta = fabs(a[j]);
1135 +         if (r_best < r) r_best = r, (*delta) = fabs(a[j]);
1136        }
1137        if (r_best < 0.001) r_best = 0.0;
1138        if (r_best == 0.0) goto done;
1139 -      xassert(delta > 0.0);
1140 +      xassert((*delta) > 0.0);
1141        /* try to increase violation by dividing delta by 2, 4, and 8,
1142           respectively */
1143 -      d_try[1] = delta / 2.0;
1144 -      d_try[2] = delta / 4.0;
1145 -      d_try[3] = delta / 8.0;
1146 +      d_try[1] = (*delta) / 2.0;
1147 +      d_try[2] = (*delta) / 4.0;
1148 +      d_try[3] = (*delta) / 8.0;
1149        for (j = 1; j <= 3; j++)
1150        {  /* construct c-MIR inequality */
1151           fail = cmir_ineq(n, a, b, u, cset, d_try[j], alpha, beta,
1152 @@ -827,7 +864,7 @@ static double cmir_sep(const int n, const double a[], const double b,
1153           /* compute violation */
1154           r = - (*beta) - (*gamma) * s;
1155           for (k = 1; k <= n; k++) r += alpha[k] * x[k];
1156 -         if (r_best < r) r_best = r, delta = d_try[j];
1157 +         if (r_best < r) r_best = r, (*delta) = d_try[j];
1158        }
1159        /* build subset of variables lying strictly between their bounds
1160           and order it by nondecreasing values of |x[j] - u[j]/2| */
1161 @@ -849,7 +886,7 @@ static double cmir_sep(const int n, const double a[], const double b,
1162           /* replace x[j] by its complement or vice versa */
1163           cset[j] = (char)!cset[j];
1164           /* construct c-MIR inequality */
1165 -         fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
1166 +         fail = cmir_ineq(n, a, b, u, cset, (*delta), alpha, beta, gamma);
1167           /* restore the variable */
1168           cset[j] = (char)!cset[j];
1169           /* do not replace the variable in case of failure */
1170 @@ -860,10 +897,11 @@ static double cmir_sep(const int n, const double a[], const double b,
1171           if (r_best < r) r_best = r, cset[j] = (char)!cset[j];
1172        }
1173        /* construct the best c-MIR inequality chosen */
1174 -      fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
1175 +      fail = cmir_ineq(n, a, b, u, cset, (*delta), alpha, beta, gamma);
1176        xassert(!fail);
1177 +
1178  done: /* free working arrays */
1179 -      xfree(cset);
1180 +
1181        xfree(vset);
1182        /* return to the calling routine */
1183        return r_best;
1184 @@ -874,7 +912,8 @@ static double generate(struct MIR *mir)
1185        int m = mir->m;
1186        int n = mir->n;
1187        int j, k, kk, nint;
1188 -      double s, *u, *x, *alpha, r_best = 0.0, b, beta, gamma;
1189 +      double s, *u, *x, *alpha, r_best = 0.0, b, beta, gamma, delta;
1190 +      char *cset;
1191        ios_copy_vec(mir->cut_vec, mir->mod_vec);
1192        mir->cut_rhs = mir->mod_rhs;
1193        /* remove small terms, which can appear due to substitution of
1194 @@ -923,6 +962,7 @@ static double generate(struct MIR *mir)
1195        u = xcalloc(1+nint, sizeof(double));
1196        x = xcalloc(1+nint, sizeof(double));
1197        alpha = xcalloc(1+nint, sizeof(double));
1198 +      cset = xcalloc(1+nint, sizeof(char));
1199        /* determine u and x */
1200        for (j = 1; j <= nint; j++)
1201        {  k = mir->cut_vec->ind[j];
1202 @@ -936,6 +976,7 @@ static double generate(struct MIR *mir)
1203              x[j] = mir->ub[k] - mir->x[k];
1204           else
1205              xassert(k != k);
1206 +         if(!(x[j] >= -0.001)) { goto skip; }
1207           xassert(x[j] >= -0.001);
1208           if (x[j] < 0.0) x[j] = 0.0;
1209        }
1210 @@ -965,6 +1006,7 @@ static double generate(struct MIR *mir)
1211           }
1212           else
1213              xassert(k != k);
1214 +         if(!(x >= -0.001)) { goto skip; }
1215           xassert(x >= -0.001);
1216           if (x < 0.0) x = 0.0;
1217           s -= mir->cut_vec->val[j] * x;
1218 @@ -973,7 +1015,7 @@ static double generate(struct MIR *mir)
1219        /* apply heuristic to obtain most violated c-MIR inequality */
1220        b = mir->cut_rhs;
1221        r_best = cmir_sep(nint, mir->cut_vec->val, b, u, x, s, alpha,
1222 -         &beta, &gamma);
1223 +         &beta, &gamma, &delta, cset);
1224        if (r_best == 0.0) goto skip;
1225        xassert(r_best > 0.0);
1226        /* convert to raw cut */
1227 @@ -988,10 +1030,22 @@ static double generate(struct MIR *mir)
1228  #if _MIR_DEBUG
1229        ios_check_vec(mir->cut_vec);
1230  #endif
1231 +      /* added */
1232 +      mir->cut_delta = delta;
1233 +      // this is not a great place for resetting the array,
1234 +      // but it should be sufficient
1235 +      for (j = 1; j <= n+m;  j++)
1236 +         mir->cut_cset[j] = 0;
1237 +      for (j = 1; j <= nint; j++)
1238 +      {  k = mir->cut_vec->ind[j];
1239 +         xassert(m <= k  && k <= m+n);
1240 +         mir->cut_cset[k] = cset[j];
1241 +      }
1242  skip: /* free working arrays */
1243        xfree(u);
1244        xfree(x);
1245        xfree(alpha);
1246 +      xfree(cset);
1247  done: return r_best;
1248  }
1249
1250 @@ -1199,8 +1253,25 @@ static void add_cut(glp_tree *tree, struct MIR *mir)
1251        ios_add_cut_row(tree, pool, GLP_RF_MIR, len, ind, val, GLP_UP,
1252           mir->cut_rhs);
1253  #else
1254 -      glp_ios_add_row(tree, NULL, GLP_RF_MIR, 0, len, ind, val, GLP_UP,
1255 +      int ord;
1256 +      ord = glp_ios_add_row(tree, NULL, GLP_RF_MIR, 0, len, ind, val, GLP_UP,
1257           mir->cut_rhs);
1258 +      ios_cut_set_aux(tree, ord, mir->agg_cnt, mir->agg_row, mir->agg_coeffs);
1259 +      ios_cut_set_aux_mir(tree, ord, mir->cut_delta,
1260 +                          mir->cut_cset, mir->subst,
1261 +                          mir->vlb_rows, mir->vub_rows);
1262 +
1263 +      /** callback for a cut being added to the cut pool */
1264 +      /* printf("mir tree parm %p %d\n", tree->parm->cb_func, ord); */
1265 +      /* printf("  agg_rhs %f\n", mir->agg_rhs); */
1266 +      /* printf("  mod_rhs %f\n", mir->mod_rhs); */
1267 +      /* printf("  cut_rhs %f\n", mir->cut_rhs); */
1268 +      if (tree->parm->cb_func != NULL)
1269 +      {  xassert(tree->reason == GLP_ICUTGEN);
1270 +         tree->reason = GLP_ICUTADDED;
1271 +         tree->parm->cb_func(tree, tree->parm->cb_info);
1272 +         tree->reason = GLP_ICUTGEN;
1273 +      }
1274  #endif
1275        xfree(ind);
1276        xfree(val);
1277 @@ -1216,6 +1287,7 @@ static int aggregate_row(glp_tree *tree, struct MIR *mir)
1278        IOSVEC *v;
1279        int ii, j, jj, k, kk, kappa = 0, ret = 0;
1280        double d1, d2, d, d_max = 0.0;
1281 +      double guass_coeff;
1282        /* choose appropriate structural variable in the aggregated row
1283           to be substituted */
1284        for (j = 1; j <= mir->agg_vec->nnz; j++)
1285 @@ -1300,8 +1372,9 @@ static int aggregate_row(glp_tree *tree, struct MIR *mir)
1286        xassert(j != 0);
1287        jj = v->pos[kappa];
1288        xassert(jj != 0);
1289 -      ios_linear_comb(mir->agg_vec,
1290 -         - mir->agg_vec->val[j] / v->val[jj], v);
1291 +      guass_coeff = - mir->agg_vec->val[j] / v->val[jj];
1292 +      mir->agg_coeffs[mir->agg_cnt] = guass_coeff;
1293 +      ios_linear_comb(mir->agg_vec, guass_coeff, v);
1294        ios_delete_vec(v);
1295        ios_set_vj(mir->agg_vec, kappa, 0.0);
1296  #if _MIR_DEBUG
1297 @@ -1440,6 +1513,13 @@ void ios_mir_term(void *gen)
1298        xfree(mir->subst);
1299        ios_delete_vec(mir->mod_vec);
1300        ios_delete_vec(mir->cut_vec);
1301 +
1302 +      /* added */
1303 +      xfree(mir->cut_cset);
1304 +      xfree(mir->vlb_rows);
1305 +      xfree(mir->vub_rows);
1306 +      xfree(mir->agg_coeffs);
1307 +
1308        xfree(mir);
1309        return;
1310  }
1311 diff --git a/src/glpios07.c b/src/glpios07.c
1312 index 0892550..9f742e6 100644
1313 --- a/src/glpios07.c
1314 +++ b/src/glpios07.c
1315 @@ -543,6 +543,14 @@ void ios_cov_gen(glp_tree *tree)
1316           /* add the cut to the cut pool */
1317           glp_ios_add_row(tree, NULL, GLP_RF_COV, 0, len, ind, val,
1318              GLP_UP, val[0]);
1319 +
1320 +         /** callback for a cut being added to the cut pool */
1321 +         if (tree->parm->cb_func != NULL)
1322 +         {  xassert(tree->reason == GLP_ICUTGEN);
1323 +            tree->reason = GLP_ICUTADDED;
1324 +            tree->parm->cb_func(tree, tree->parm->cb_info);
1325 +            tree->reason = GLP_ICUTGEN;
1326 +         }
1327        }
1328        /* free working arrays */
1329        xfree(ind);
1330 diff --git a/src/glpios08.c b/src/glpios08.c
1331 index 2d2fcd5..3aad808 100644
1332 --- a/src/glpios08.c
1333 +++ b/src/glpios08.c
1334 @@ -129,6 +129,15 @@ void ios_clq_gen(glp_tree *T, void *G_)
1335        /* add cut inequality to local cut pool */
1336        glp_ios_add_row(T, NULL, GLP_RF_CLQ, 0, len, ind, val, GLP_UP,
1337           rhs);
1338 +
1339 +      /** callback for a cut being added to the cut pool */
1340 +      if (T->parm->cb_func != NULL)
1341 +      {  xassert(T->reason == GLP_ICUTGEN);
1342 +         T->reason = GLP_ICUTADDED;
1343 +         T->parm->cb_func(T, T->parm->cb_info);
1344 +         T->reason = GLP_ICUTGEN;
1345 +      }
1346 +
1347  skip: /* free working arrays */
1348        tfree(ind);
1349        tfree(val);
1350 diff --git a/src/glpios11.c b/src/glpios11.c
1351 index c40e9a5..82d5698 100644
1352 --- a/src/glpios11.c
1353 +++ b/src/glpios11.c
1354 @@ -193,6 +193,9 @@ void ios_process_cuts(glp_tree *T)
1355           glp_set_mat_row(T->mip, i, len, ind, val);
1356           xassert(cut->type == GLP_LO || cut->type == GLP_UP);
1357           glp_set_row_bnds(T->mip, i, cut->type, cut->rhs, cut->rhs);
1358 +
1359 +         /* setting this as selected */
1360 +         ios_cut_set_selected(cut, i);
1361        }
1362        /* free working arrays */
1363        xfree(info);
1364 diff --git a/src/glpk.h b/src/glpk.h
1365 index 75c292c..e117782 100644
1366 --- a/src/glpk.h
1367 +++ b/src/glpk.h
1368 @@ -130,6 +130,7 @@ typedef struct
1369        int out_frq;            /* spx.out_frq */
1370        int out_dly;            /* spx.out_dly (milliseconds) */
1371        int presolve;           /* enable/disable using LP presolver */
1372 +      int stability_lmt;      /* maximum number of check stability failures before stopping */
1373        double foo_bar[36];     /* (reserved) */
1374  } glp_smcp;
1375
1376 @@ -229,6 +230,11 @@ typedef struct
1377  #define GLP_IBRANCH     0x05  /* request for branching */
1378  #define GLP_ISELECT     0x06  /* request for subproblem selection */
1379  #define GLP_IPREPRO     0x07  /* request for preprocessing */
1380 +#define GLP_ICUTADDED   0x08  /* cut was added to the pool */
1381 +#define GLP_ICUTSELECT  0x09  /* cuts were selected as rows */
1382 +#define GLP_LI_BRANCH   0x10  /* a branch was made */
1383 +#define GLP_LI_CLOSE    0x11  /* an active node was closed */
1384 +#define GLP_LI_DELROW   0x12  /* an active node was closed */
1385
1386  /* branch selection indicator: */
1387  #define GLP_NO_BRNCH       0  /* select no branch */
1388 @@ -1057,6 +1063,54 @@ int glp_top_sort(glp_graph *G, int v_num);
1389  int glp_wclique_exact(glp_graph *G, int v_wgt, double *sol, int v_set);
1390  /* find maximum weight clique with exact algorithm */
1391
1392 +/*******************************************/
1393 +/***               CUT LOG               ***/
1394 +/*******************************************/
1395 +
1396 +int glp_get_it_cnt(glp_prob *P);
1397 +/* get the iteration count of the current problem */
1398 +
1399 +
1400 +int glp_ios_get_cut(glp_tree *T, int i, int ind[], double val[], int* klass, int* type, double* rhs);
1401 +/* determine reason for calling the callback routine */
1402 +
1403 +int glp_ios_cut_get_aux_nrows(glp_tree *tree, int ord);
1404 +/* gets the number of rows used to generate a cut. */
1405 +
1406 +
1407 +void glp_ios_cut_get_aux_rows(glp_tree *tree, int ord,
1408 +                              int rows[], double coeffs[]);
1409 +/* gets a cut as an input sequence of rows times coefficients. */
1410 +
1411 +
1412 +void glp_ios_cut_get_mir_cset(glp_tree *tree, int ord, char cset[]);
1413 +/* gets mir cut complement set. */
1414 +double glp_ios_cut_get_mir_delta(glp_tree *tree, int ord);
1415 +/* gets mir cut delta. */
1416 +void glp_ios_cut_get_mir_subst(glp_tree *tree, int ord, char subst[]);
1417 +/* gets mir cut substition information. */
1418 +void glp_ios_cut_get_mir_virtual_rows(glp_tree *tree, int ord,
1419 +                                      int vlb_rows[], int vub_rows[]);
1420 +/* gets mir cut virtual bounds rows. */
1421 +
1422 +int glp_ios_selected_cuts(glp_tree *tree, int ords[], int sel[]);
1423 +/* gets the list of selected cuts.
1424 +   Can only be called when GLP_ICUTSELECT */
1425 +
1426 +
1427 +int glp_ios_branch_log(glp_tree *tree, double *val, int* parent, int* dn, int* up);
1428 +/* can only be called when GLP_LI_BRANCH.
1429 + * If id is non-null, returns the id of the structural variable branched upon.
1430 + * If val is non-null, it is set to the value branched upon.
1431 + * If parent is non-null, it is set to node id of the node branched upon.
1432 + * If dn is non-null, it is set to node id of the newly created down node.
1433 + * If up is non-null, it is set to node id of the newly created up node.
1434 + */
1435 +
1436 +int glp_ios_node_ord(glp_tree *tree, int node_p);
1437 +
1438 +int glp_ios_rows_deleted(glp_tree *tree, int* rows);
1439 +
1440  #ifdef __cplusplus
1441  }
1442  #endif
1443 diff --git a/src/glpspx01.c b/src/glpspx01.c
1444 index 5c17114..caaab27 100644
1445 --- a/src/glpspx01.c
1446 +++ b/src/glpspx01.c
1447 @@ -238,6 +238,9 @@ struct csa
1448        double *work2; /* double work2[1+m]; */
1449        double *work3; /* double work3[1+m]; */
1450        double *work4; /* double work4[1+m]; */
1451 +
1452 +      /** Things Tim has added. */
1453 +      int stability_failures;
1454  };
1455
1456  static const double kappa = 0.10;
1457 @@ -400,6 +403,8 @@ static void init_csa(struct csa *csa, glp_prob *lp)
1458        csa->refct = 0;
1459        memset(&refsp[1], 0, (m+n) * sizeof(char));
1460        for (j = 1; j <= n; j++) gamma[j] = 1.0;
1461 +
1462 +      csa->stability_failures = 0;
1463        return;
1464  }
1465
1466 @@ -2647,6 +2652,11 @@ loop: /* main loop starts here */
1467           if (check_stab(csa, parm->tol_bnd))
1468           {  /* there are excessive bound violations due to round-off
1469                 errors */
1470 +            csa->stability_failures++;
1471 +            if (csa->stability_failures >= parm->stability_lmt){
1472 +              ret = GLP_EINSTAB;
1473 +              goto done;
1474 +            }
1475              if (parm->msg_lev >= GLP_MSG_ERR)
1476                 xprintf("Warning: numerical instability (primal simplex,"
1477                    " phase %s)\n", csa->phase == 1 ? "I" : "II");
1478 diff --git a/src/glpspx02.c b/src/glpspx02.c
1479 index 19152a6..de3d677 100644
1480 --- a/src/glpspx02.c
1481 +++ b/src/glpspx02.c
1482 @@ -288,6 +288,9 @@ struct csa
1483        double *work2; /* double work2[1+m]; */
1484        double *work3; /* double work3[1+m]; */
1485        double *work4; /* double work4[1+m]; */
1486 +
1487 +      /* count the number of stability failures */
1488 +      int stability_failures;
1489  };
1490
1491  static const double kappa = 0.10;
1492 @@ -501,6 +504,8 @@ static void init_csa(struct csa *csa, glp_prob *lp)
1493        csa->refct = 0;
1494        memset(&refsp[1], 0, (m+n) * sizeof(char));
1495        for (i = 1; i <= m; i++) gamma[i] = 1.0;
1496 +
1497 +      csa->stability_failures = 0;
1498        return;
1499  }
1500
1501 @@ -2741,7 +2746,13 @@ loop: /* main loop starts here */
1502           /* make sure that the current basic solution remains dual
1503              feasible */
1504           if (check_stab(csa, parm->tol_dj) != 0)
1505 -         {  if (parm->msg_lev >= GLP_MSG_ERR)
1506 +         {
1507 +            csa->stability_failures++;
1508 +            if (csa->stability_failures >= parm->stability_lmt){
1509 +              ret = GLP_EINSTAB;
1510 +              goto done;
1511 +            }
1512 +            if (parm->msg_lev >= GLP_MSG_ERR)
1513                 xprintf("Warning: numerical instability (dual simplex, p"
1514                    "hase %s)\n", csa->phase == 1 ? "I" : "II");
1515  #if 1
1516 @@ -3023,6 +3034,10 @@ loop: /* main loop starts here */
1517        /* accuracy check based on the pivot element */
1518        {  double piv1 = csa->tcol_vec[csa->p]; /* more accurate */
1519           double piv2 = csa->trow_vec[csa->q]; /* less accurate */
1520 +         if(piv1 == 0.0){
1521 +           ret = GLP_EFAIL;
1522 +           goto done;
1523 +         }
1524           xassert(piv1 != 0.0);
1525           if (fabs(piv1 - piv2) > 1e-8 * (1.0 + fabs(piv1)) ||
1526               !(piv1 > 0.0 && piv2 > 0.0 || piv1 < 0.0 && piv2 < 0.0))
1527 --
1528 2.26.2