Migrating from using the 'glpk-cut-log' repo to using an official GPLK release + a set of patches (#4775)

This PR attempts to migrate from @timothy-king's repo for glpk-cut-log to using a set of patches against the official release that 'glpk-cut-log' was based off of (4.52).

This is in preparation of trying to rework these patches to be against the latest GPLK release (4.65). If we can do this, then it makes it easier for CVC4's dependancy on GPLK to be able to follow upstream (rather than being stuck on a release that is 6.5 years old!).

I have added GPLK as an option for CI, but I do not know if we actually want this. My concern with adding this to CI is that we're then not testing non-GPL builds, which doesn't seem ideal.

However, before starting to rework the patches to be against 4.65, I want to make sure that things are actually working/stable against 4.52; so having at least one CI target that does use GPLK would be great.

Signed-off-by: Andrew V. Jones andrew.jones@vector.com

2 files changed, 1539 insertions, 2 deletions

diff --git a/contrib/get-glpk-cut-log b/contrib/get-glpk-cut-log
index 16acf97ae..951e10620 100755
--- a/contrib/get-glpk-cut-log
+++ b/contrib/get-glpk-cut-log
@@ -2,15 +2,24 @@
 #
 source "$(dirname "$0")/get-script-header.sh"
 
+# get the full path to the contrib dir; needs to be the full path so we can
+# refer the patch from within the 'glpk-cut-log' build directory
+contrib_dir="$( cd "$( dirname "${BASH_SOURCE[0]}" )" >/dev/null 2>&1 && pwd -P)"
+
+# Get the full path the patch we wish to apply
+patch_file=${contrib_dir}/glpk-cut-log.patch
+
 GLPK_DIR="$DEPS_DIR/glpk-cut-log"
-commit=b420454e732f4b3d229c552ef7cd46fec75fe65c
+version="4.52"
 
 check_dep_dir "$GLPK_DIR"
 setup_dep \
-  "https://github.com/timothy-king/glpk-cut-log/archive/$commit.tar.gz" \
+  "https://ftp.gnu.org/gnu/glpk/glpk-${version}.tar.gz" \
   "$GLPK_DIR"
 cd "$GLPK_DIR"
 
+patch -p1 < ${patch_file}
+
 libtoolize
 aclocal
 autoheader
diff --git a/contrib/glpk-cut-log.patch b/contrib/glpk-cut-log.patch
new file mode 100644
index 000000000..9ef01886c
--- /dev/null
+++ b/contrib/glpk-cut-log.patch
@@ -0,0 +1,1528 @@
+@@ This patch is taken from https://github.com/timothy-king/glpk-cut-log and
+@@ has the following license:
+@@
+@@ GLPK is free software: you can redistribute it and/or modify it under the
+@@ terms of the GNU General Public License as published by the Free Software
+@@ Foundation, either version 3 of the License, or (at your option) any later
+@@ version.
+@@
+From cf84e3855d8c5676557daaca434b42b2fc17dc29 Mon Sep 17 00:00:00 2001
+From: Tim King <taking@cs.nyu.edu>
+Date: Sat, 14 Dec 2013 16:03:35 -0500
+Subject: [PATCH] Adding a function for returning the iteration count.
+
+Adding GLP_ICUTADDED callback:
+- Adds a new callback location that is called after every cut is added to the pool.
+- During this callback users can request a copy of this cut and query what kind of cut it is.
+- GMI cuts also support returning what row in the tableau it came from. Users can use this to replay how the cut was generated.
+
+Logging MIR Cuts
+- Changed the IOSAUX to be defined as a linear sum of input rows.
+  This is visible from the interface
+- Logging the multipliers used to create the aggregate row in MIR.
+- Logging the delta selected by final successful cmir_sep call.
+- Logging the cset selected by final successful cmir_sep call.
+- The delta and cset are not yet user visible.
+
+The extra mir logging information is now user visible.
+
+Adding GLP_ICUTSELECT callback:
+- ios_process_cuts marks the cuts as selected after turning the cut into a row.
+- This callback happens after ios_process_cuts marks and before the pool is
+  cleared by ios_clear_pool.
+
+Cleaning up the git directory to not track autogenerated files.
+
+Branch logging
+- Adds to the integer interface a callback function for when branches are made.
+  This is called after branch_on.
+- Added a public function glp_ios_branch_log.
+  This returns the structural variable for the branch as well as the value
+  branched on, and the ids of the nodes for the down and up branches.
+  If both branches are infeasible, the ids are -1 for both dn and up.
+  If one branch is infeasible (and no branch is done), this returns -1 for
+  the infeasible branch and the current node id for the up branch.
+- Each node now has a unique ordinal associated with it.
+- Added a public function to convert the id of an active node to the unique
+  node ordinal.
+- Added a callback for when a node is closed due to being linearly infeasible.
+
+Improved the mir cut logging facilities to now include:
+- the rid id used for virtual lower/upper bound selection,
+- and the subst map,
+- Changed the size of the cset array to be the same as the new arrays (m+n)
+- Fixed a bug in returning the cset selected.
+
+Fixing a bug selecting the correct branch.
+
+Added a callback for tracking row deletion.
+
+Adding notes to node freeze and revive functions.
+
+Commenting out a few debugging printfs for gomory cuts.
+
+Changing an overly aggressive assert in MIR generate when an assignment is outside of the bounds to skipping the application of the mir rule.
+
+Removing some printfs.
+
+Weakening assert to a failure.
+
+Adding a stability failure limit to simplex.
+
+Updating the installation instructions.
+
+---
+diff --git a/configure.ac b/configure.ac
+index 0141a8a..ffb9ae4 100644
+--- a/configure.ac
++++ b/configure.ac
+@@ -1,7 +1,7 @@
+ dnl Process this file with autoconf to produce a configure script
+ 
+ AC_INIT([GLPK], [4.52], [bug-glpk@gnu.org])
+-
++AM_INIT_AUTOMAKE([subdir-objects])
+ AC_CONFIG_SRCDIR([src/glpk.h])
+ 
+ AC_CONFIG_MACRO_DIR([m4])
+diff --git a/src/Makefile.am b/src/Makefile.am
+index b39efa6..2a025bc 100644
+--- a/src/Makefile.am
++++ b/src/Makefile.am
+@@ -20,7 +20,10 @@ libglpk_la_LDFLAGS = \
+ -version-info 36:0:1 \
+ -export-symbols-regex '^glp_*'
+ 
++
++# all of the cut log sources are listed first
+ libglpk_la_SOURCES = \
++cutlog01.c \
+ glpapi01.c \
+ glpapi02.c \
+ glpapi03.c \
+diff --git a/src/cutlog01.c b/src/cutlog01.c
+new file mode 100644
+index 0000000..9a85bb7
+--- /dev/null
++++ b/src/cutlog01.c
+@@ -0,0 +1,452 @@
++/* cutlog01.c (api extension routines) */
++
++
++#include "glpapi.h"
++#include "glpios.h"
++
++int glp_get_it_cnt(glp_prob *P){
++  if(P == NULL){
++    return 0;
++  }else{
++    return P->it_cnt;
++  }
++}
++
++int glp_ios_get_cut(glp_tree *T, int i, int* ind, double* val, int* klass, int* type, double* rhs){
++  xassert(T != NULL);
++
++  IOSCUT* cut;
++  int len;
++  IOSAIJ* aij;
++  glp_prob* prob;
++
++  if (T->reason != GLP_ICUTADDED){
++    xerror("glp_ios_get_cut: not called during cut added.\n");
++  }
++  cut = ios_find_row(T->local, i);
++  if ( cut == NULL ) {
++    xerror("glp_ios_get_cut: called with an invalid index.");
++  }
++  len = 0;
++  for(len = 0, aij = cut->ptr; aij != NULL; aij = aij->next)
++  {
++    len++;
++    if(ind != NULL){ ind[len] = aij->j; }
++    if(val != NULL){ val[len] = aij->val; }
++  }
++  if(klass != NULL){ *klass = cut->klass; }
++  if(type != NULL){ *type = cut->type; }
++  if(rhs != NULL){ *rhs = cut->rhs; }
++  return len;
++}
++
++
++IOSAUX *ios_create_aux(int n, const int rows[], const double coeffs[]){
++  IOSAUX *aux;
++  int i;
++  aux = xmalloc(sizeof(IOSAUX));
++  aux->nrows = n;
++  aux->rows = xcalloc(1+n, sizeof(int));
++  aux->mult = xcalloc(1+n, sizeof(double));
++  aux->selected = -1;
++  aux->mir_cset = NULL;
++
++  for ( i = 1; i <= n; i++)
++  {  aux->rows[i] = rows[i];
++     aux->mult[i] = coeffs[i];
++  }
++
++  return aux;
++}
++
++void ios_delete_aux(IOSAUX *aux){
++  xassert(aux != NULL);
++  xfree(aux->rows);
++  xfree(aux->mult);
++  if( aux->mir_cset != NULL ){
++    xfree(aux->mir_cset);
++     xfree(aux->mir_subst);
++     xfree(aux->mir_vlb_rows);
++     xfree(aux->mir_vub_rows);
++  }
++  xfree(aux);
++  return;
++}
++
++static void cut_set_aux(IOSCUT *cut, int n,
++                        const int rows[], const double coeffs[]){
++  int i;
++  xassert( cut != NULL );
++  if( cut->aux != NULL ) {
++    ios_delete_aux(cut-> aux);
++  }
++
++  cut->aux = ios_create_aux(n, rows, coeffs);
++  xassert( cut->aux->nrows == n );
++}
++
++static void cut_set_aux_mir(IOSAUX *aux, double delta, int m, int n,
++                            const char cset[], const char subst[],
++                            const int vlbrs[], const int vubrs[]){
++  int i;
++  xassert( aux != NULL );
++  if ( aux->mir_cset != NULL )
++  {  xfree(aux->mir_cset);
++     xfree(aux->mir_subst);
++     xfree(aux->mir_vlb_rows);
++     xfree(aux->mir_vub_rows);
++  }
++
++  aux->mir_cset     = xcalloc(1+n+m, sizeof(char));
++  aux->mir_subst    = xcalloc(1+n+m, sizeof(char));
++  aux->mir_vlb_rows = xcalloc(1+n+m, sizeof(int));
++  aux->mir_vub_rows = xcalloc(1+n+m, sizeof(int));
++
++  for ( i = 1; i <= n+m; i++)
++  {  aux->mir_cset[i] = cset[i];
++     aux->mir_subst[i] = subst[i];
++     aux->mir_vlb_rows[i] = vlbrs[i];
++     aux->mir_vub_rows[i] = vubrs[i];
++  }
++
++  aux->mir_delta = delta;
++}
++
++void ios_cut_set_aux_mir(glp_tree *T, int ord, double delta,
++                         const char cset[], const char subst[],
++                         const int vlbrs[], const int vubrs[]){
++  int m, n;
++  IOSCUT *cut;
++  glp_prob *mip;
++  mip = T->mip;
++  m = mip->m;
++  n = mip->n;
++  cut = ios_find_row(T->local, ord);
++  xassert(cut != NULL);
++  cut_set_aux_mir(cut->aux, delta, m, n, cset, subst, vlbrs, vubrs);
++}
++
++void ios_cut_set_single_aux(glp_tree *T, int ord, int j){
++  IOSCUT *cut;
++  cut = ios_find_row(T->local, ord);
++  xassert(cut != NULL);
++
++  /* set up arrays */
++  int ind[1+1];
++  double coeffs[1+1];
++  ind[1] = j;
++  coeffs[1] = +1.0;
++
++  /* call general procedure */
++  cut_set_aux(cut, 1, ind, coeffs);
++}
++
++void ios_cut_set_aux(glp_tree *T, int ord, int n,
++                     const int rows[], const double coeffs[]){
++  IOSCUT *cut;
++  cut = ios_find_row(T->local, ord);
++  xassert(cut != NULL);
++  cut_set_aux(cut, n, rows, coeffs);
++}
++
++int glp_ios_cut_get_aux_nrows(glp_tree *tree, int ord){
++  IOSCUT *cut;
++  IOSAUX *aux;
++  if (tree->reason != GLP_ICUTADDED){
++    xerror("glp_ios_cut_get_aux_nrows: not called during cut added.\n");
++  }
++  cut = ios_find_row(tree->local, ord);
++  if ( cut == NULL ){
++    xerror("glp_ios_cut_get_aux_nrows: not called on a valid cut.\n");
++  }
++  aux = cut->aux;
++  return (aux == NULL) ? 0 : aux->nrows;
++}
++
++void glp_ios_cut_get_aux_rows(glp_tree *tree, int ord,
++                              int rows[], double coeffs[]){
++  IOSCUT *cut;
++  IOSAUX *aux;
++  int j, nrows;
++  if (tree->reason != GLP_ICUTADDED){
++    xerror("glp_ios_cut_get_aux_rows: not called during cut added.\n");
++  }
++  cut = ios_find_row(tree->local, ord);
++  if ( cut == NULL ){
++    xerror("glp_ios_cut_get_aux_rows: not called on a valid cut.\n");
++  }
++  aux = cut->aux;
++  if( aux != NULL ){
++    nrows = aux->nrows;
++    for ( j = 1; j <= nrows; j++ )
++    {  if ( rows != NULL ) { rows[j] = aux->rows[j]; }
++       if ( coeffs != NULL ) { coeffs[j] = aux->mult[j]; }
++    }
++  }
++  return;
++}
++
++
++void glp_ios_cut_get_mir_subst(glp_tree *tree, int ord, char subst[]);
++/* gets mir cut substition information. */
++void glp_ios_cut_get_mir_virtual_rows(glp_tree *tree, int ord,
++                                      int vlb[], int vub[]);
++/* gets mir cut virtual bounds rows. */
++
++void glp_ios_cut_get_mir_cset(glp_tree *tree, int ord, char *cset){
++  glp_prob *mip;
++  IOSCUT *cut;
++  IOSAUX *aux;
++  int j, n, m;
++  if ( tree->reason != GLP_ICUTADDED ){
++    xerror("glp_ios_cut_get_aux_mir: not called during cut added.\n");
++  }
++  cut = ios_find_row(tree->local, ord);
++  if ( cut == NULL ){
++    xerror("glp_ios_cut_get_aux_mir: not called on a cut.\n");
++  }
++  if ( cut->klass != GLP_RF_MIR ){
++    xerror("glp_ios_cut_get_aux_mir: not called on a mir cut.\n");
++  }
++  aux = cut->aux;
++  mip = tree->mip;
++  m = mip->m;
++  n = mip->n;
++
++  if( cset != NULL ){
++    for ( j=1; j <= n+m; j++ ){
++      cset[j] = (aux->mir_cset == NULL) ? 0 : aux->mir_cset[j];
++    }
++  }
++}
++void glp_ios_cut_get_mir_subst(glp_tree *tree, int ord, char *subst){
++  glp_prob *mip;
++  IOSCUT *cut;
++  IOSAUX *aux;
++  int j, n, m;
++  if ( tree->reason != GLP_ICUTADDED ){
++    xerror("glp_ios_cut_get_aux_mir: not called during cut added.\n");
++  }
++  cut = ios_find_row(tree->local, ord);
++  if ( cut == NULL ){
++    xerror("glp_ios_cut_get_aux_mir: not called on a cut.\n");
++  }
++  if ( cut->klass != GLP_RF_MIR ){
++    xerror("glp_ios_cut_get_aux_mir: not called on a mir cut.\n");
++  }
++  aux = cut->aux;
++  mip = tree->mip;
++  m = mip->m;
++  n = mip->n;
++
++  if( subst != NULL ){
++    for ( j=1; j <= n+m; j++ ){
++      subst[j] = (aux->mir_subst == NULL) ? 0 : aux->mir_subst[j];
++    }
++  }
++}
++void glp_ios_cut_get_mir_virtual_rows(glp_tree *tree, int ord, int vlb_rows[], int vub_rows[]){
++  glp_prob *mip;
++  IOSCUT *cut;
++  IOSAUX *aux;
++  int j, n, m;
++  if ( tree->reason != GLP_ICUTADDED ){
++    xerror("glp_ios_cut_get_aux_mir: not called during cut added.\n");
++  }
++  cut = ios_find_row(tree->local, ord);
++  if ( cut == NULL ){
++    xerror("glp_ios_cut_get_aux_mir: not called on a cut.\n");
++  }
++  if ( cut->klass != GLP_RF_MIR ){
++    xerror("glp_ios_cut_get_aux_mir: not called on a mir cut.\n");
++  }
++  aux = cut->aux;
++  mip = tree->mip;
++  m = mip->m;
++  n = mip->n;
++
++  for ( j=1; j <= n+m; j++ ){
++    vlb_rows[j] = (aux->mir_vlb_rows == NULL) ? 0 : aux->mir_vlb_rows[j];
++    vub_rows[j] = (aux->mir_vub_rows == NULL) ? 0 : aux->mir_vub_rows[j];
++  }
++}
++double glp_ios_cut_get_mir_delta(glp_tree *tree, int ord){
++  glp_prob *mip;
++  IOSCUT *cut;
++  IOSAUX *aux;
++  int j, n, m;
++  if ( tree->reason != GLP_ICUTADDED ){
++    xerror("glp_ios_cut_get_aux_mir: not called during cut added.\n");
++  }
++  cut = ios_find_row(tree->local, ord);
++  if ( cut == NULL ){
++    xerror("glp_ios_cut_get_aux_mir: not called on a cut.\n");
++  }
++  if ( cut->klass != GLP_RF_MIR ){
++    xerror("glp_ios_cut_get_aux_mir: not called on a mir cut.\n");
++  }
++  aux = cut->aux;
++  return aux->mir_delta;
++}
++
++
++void ios_cut_set_selected(IOSCUT *cut, int sel){
++#ifdef CUT_DEBUG
++  static int i = 0;
++  ++i;
++  printf("ios_cut_set_selected: %d %d %p\n", i, sel, cut);
++#endif
++
++  IOSAUX *aux;
++  aux = cut->aux;
++  if ( aux != NULL ){
++    aux->selected = sel;
++  }
++}
++
++int glp_ios_selected_cuts(glp_tree *tree, int ords[], int sel[]){
++  int len, j, N, s;
++  IOSPOOL* pool;
++  IOSCUT* cut;
++  IOSAUX* aux;
++  if ( tree == NULL ){
++    xerror("glp_ios_selected_cuts: not called with a valid tree.\n");
++  }
++  if ( tree->reason != GLP_ICUTSELECT ){
++    xerror("glp_ios_selected_cuts: not called during cut selected.\n");
++  }
++
++  pool = tree->local;
++  if ( pool == NULL ){
++    xerror("glp_ios_selected_cuts: called on a malformed tree.\n");
++  }
++
++  for (len = 0, j = 1, cut = pool->head; cut != NULL; cut = cut->next, j++)
++  {  aux = cut->aux;
++#ifdef CUT_DEBUG
++     printf("glp_ios_selected_cuts: %d %p\n", j, cut);
++#endif
++     if ( aux != NULL )
++     { s = aux->selected;
++       if ( s >= 0 )
++       {  len++;
++          if (ords != NULL) { ords[len] = j; }
++          if (sel != NULL)  { sel[len] = s; }
++       }
++     }
++  }
++  return len;
++}
++
++int glp_ios_branch_log(glp_tree *tree, double *br_val, int* parent, int* dn, int* up){
++  IOSNPD *node;
++  int br_result, br_var;
++  int p, d, u;
++  double v;
++  glp_prob *mip;
++  if ( tree == NULL ){
++    xerror("glp_ios_branch_log: not called with a valid tree \n");
++  }
++  if ( tree->reason != GLP_LI_BRANCH ){
++    xerror("glp_ios_branch_log: not called during GLP_LI_BRANCH \n");
++  }
++  mip = tree->mip;
++  if ( mip == NULL ){
++    xerror("glp_ios_branch_log: not called with a valid tree\n");
++  }
++  br_result = tree->br_result;
++  br_var = tree->br_var;
++  switch(br_result){
++  case 0:
++    p = tree->br_node;
++    node = tree->slot[p].node;
++    break;
++  case 1:
++  case 2:
++    node = tree->curr;
++    p = node->p;
++    break;
++  default:
++    xerror("glp_ios_branch_log: br_result is not properly set \n");
++  }
++  if( node == NULL ){
++    xerror("glp_ios_branch_log: not called with a valid tree \n");
++  }
++  switch(br_result){
++  case 0:
++    v = node->br_val;
++    d = tree->dn_child;
++    u = tree->up_child;
++    break;
++  case 1:
++    v = mip->col[br_var]->prim;
++    if(tree->br_to_up){
++      d = -1;
++      u = p;
++    }else{
++      d = p;
++      u = -1;
++    }
++    break;
++  case 2:
++    v = mip->col[br_var]->prim;
++    d = -1;
++    u = -1;
++    break;
++  default:
++    xerror("glp_ios_branch_log: not called with a valid tree \n");
++  }
++
++  if(br_val != NULL){ *br_val = v; }
++  if(parent != NULL){ *parent = p; }
++  if(dn != NULL){ *dn = d; }
++  if(up != NULL){ *up = u; }
++
++  return br_var;
++}
++
++int glp_ios_node_ord(glp_tree *tree, int p){
++  IOSNPD *node;
++  if ( tree == NULL ){
++    xerror("glp_ios_node_ord: not called with a valid tree.\n");
++  }
++  if (!(1 <= p && p <= tree->nslots)){
++    xerror("glp_ios_node_ord: not called with a valid p.\n");
++  }
++  node = tree->slot[p].node;
++  return node->ord;
++}
++
++void ios_cb_rows_deleted(glp_tree *T, int n, const int* rows){
++  if (T->parm->cb_func != NULL)
++  {
++    xassert(T->reason == 0);
++    xassert(T->deleting_rows == NULL);
++    xassert(T->num_deleting_rows == 0);
++    T->num_deleting_rows = n;
++    T->deleting_rows = rows;
++
++    T->reason = GLP_LI_DELROW;
++    T->parm->cb_func(T, T->parm->cb_info);
++    T->reason = 0;
++    T->num_deleting_rows = 0;
++    T->deleting_rows = NULL;
++  }
++}
++
++int glp_ios_rows_deleted(glp_tree *tree, int* rows){
++  if ( tree == NULL ){
++    xerror("glp_ios_rows_deleted: not called with a valid tree.\n");
++  }
++  if ( tree->reason != GLP_LI_DELROW ){
++    xerror("glp_ios_rows_deleted: not called with a valid reason.\n");
++  }
++
++  int j;
++  if(rows != NULL){
++    for(j=1; j <= tree->num_deleting_rows; j++){
++      rows[j] = tree->deleting_rows[j];
++    }
++  }
++  return tree->num_deleting_rows;
++}
+diff --git a/src/glpapi06.c b/src/glpapi06.c
+index 743d6be..05d0498 100644
+--- a/src/glpapi06.c
++++ b/src/glpapi06.c
+@@ -488,6 +488,7 @@ void glp_init_smcp(glp_smcp *parm)
+       parm->out_frq = 500;
+       parm->out_dly = 0;
+       parm->presolve = GLP_OFF;
++      parm->stability_lmt = 200;
+       return;
+ }
+ 
+diff --git a/src/glpapi13.c b/src/glpapi13.c
+index 926dda4..55adf44 100644
+--- a/src/glpapi13.c
++++ b/src/glpapi13.c
+@@ -453,8 +453,14 @@ void glp_ios_row_attr(glp_tree *tree, int i, glp_attr *attr)
+ 
+ int glp_ios_pool_size(glp_tree *tree)
+ {     /* determine current size of the cut pool */
+-      if (tree->reason != GLP_ICUTGEN)
+-         xerror("glp_ios_pool_size: operation not allowed\n");
++      switch(tree->reason)
++      {  case GLP_ICUTGEN:
++         case GLP_ICUTADDED:
++         case GLP_ICUTSELECT:
++           break;
++         default:
++           xerror("glp_ios_pool_size: operation not allowed\n");
++      }
+       xassert(tree->local != NULL);
+       return tree->local->size;
+ }
+diff --git a/src/glpios.h b/src/glpios.h
+index 9e2d6b2..3b31901 100644
+--- a/src/glpios.h
++++ b/src/glpios.h
+@@ -36,6 +36,9 @@ typedef struct IOSAIJ IOSAIJ;
+ typedef struct IOSPOOL IOSPOOL;
+ typedef struct IOSCUT IOSCUT;
+ 
++
++typedef struct IOSAUX IOSAUX;
++
+ struct glp_tree
+ {     /* branch-and-bound tree */
+       int magic;
+@@ -217,6 +220,19 @@ struct glp_tree
+          GLP_NO_BRNCH - use general selection technique */
+       int child;
+       /* subproblem reference number corresponding to br_sel */
++
++      /* start of cut log extras */
++      int br_result;
++      int br_to_up;
++      int br_node;
++      /* subproblem reference number for the just branched node.*/
++      int dn_child;
++      /* subproblem reference number for the just created down node */
++      int up_child;
++      /* subproblem reference number for the just created up node */
++
++      const int* deleting_rows;
++      int num_deleting_rows;
+ };
+ 
+ struct IOSLOT
+@@ -229,6 +245,8 @@ struct IOSLOT
+ 
+ struct IOSNPD
+ {     /* node subproblem descriptor */
++      int ord;
++      /* this is a unique ordinal for each subproblem  */
+       int p;
+       /* subproblem reference number (it is the index to corresponding
+          slot, i.e. slot[p] points to this descriptor) */
+@@ -393,12 +411,51 @@ struct IOSCUT
+          GLP_FX: sum a[j] * x[j]  = b */
+       double rhs;
+       /* cut right-hand side */
++
++      IOSAUX *aux;
++      /* cut auxillary source information */
++
+       IOSCUT *prev;
+       /* pointer to previous cut */
+       IOSCUT *next;
+       /* pointer to next cut */
+ };
+ 
++struct IOSAUX
++{
++      /* aux (auxillary source information for each cut)
++       * Each cut operates on a sum of rows
++       * Each cut operates on a row that can be described using
++       * a current row or a previous cut.
++       * To generalize, we assume each row is a sum of two rows:
++       *   row[r]* r_mult + pool[c] * c_mult
++       */
++      int nrows;
++      int *rows;
++      double *mult;
++
++      int selected;
++      /* when < 0 this has not yet been turned into a row
++         when >=0 this is the id of the row added. */
++
++      double mir_delta;
++      /* the delta selected by a mir cut */
++
++      char *mir_cset;
++      /* complimented set */
++      /* if this is NULL, it is implicitly all 0s */
++
++      char *mir_subst;
++      /* the substition vectors */
++
++      int *mir_vlb_rows;
++      /* the substition vectors */
++
++      int *mir_vub_rows;
++      /* the substition vectors */
++};
++
++
+ #define ios_create_tree _glp_ios_create_tree
+ glp_tree *ios_create_tree(glp_prob *mip, const glp_iocp *parm);
+ /* create branch-and-bound tree */
+@@ -615,6 +672,31 @@ int ios_choose_node(glp_tree *T);
+ int ios_choose_var(glp_tree *T, int *next);
+ /* select variable to branch on */
+ 
++/* functions added to retrieve information */
++IOSAUX *ios_create_aux(int n, const int rows[], const double coeffs[]);
++/* creates an aux with n rows */
++
++void ios_delete_aux(IOSAUX *aux);
++/* deletes an aux */
++
++void ios_cut_set_single_aux(glp_tree *T, int ord, int j);
++/* sets the aux of row ord to be a single row j */
++
++
++void ios_cut_set_aux(glp_tree *T, int ord, int n,
++                     const int rows[], const double coeffs[]);
++/* sets an arbitrary aux sum */
++
++void ios_cut_set_aux_mir(glp_tree *T, int ord, double delta,
++                         const char cset[], const char subst[],
++                         const int vlbrs[], const int vubrs[]);
++/* sets the extra mir information */
++
++void ios_cut_set_selected(IOSCUT *cut, int i);
++/* the cut has been added as row i */
++
++void ios_cb_rows_deleted(glp_tree *T, int n, const int* rows);
++
+ #endif
+ 
+ /* eof */
+diff --git a/src/glpios01.c b/src/glpios01.c
+index 70798fd..d9e5703 100644
+--- a/src/glpios01.c
++++ b/src/glpios01.c
+@@ -159,6 +159,16 @@ glp_tree *ios_create_tree(glp_prob *mip, const glp_iocp *parm)
+       tree->stop = 0;
+       /* create the root subproblem, which initially is identical to
+          the original MIP */
++
++      tree->br_result = 0;
++      tree->br_to_up = 0;
++      tree->br_node = 0;
++      tree->dn_child = 0;
++      tree->up_child = 0;
++
++      tree->deleting_rows = NULL;
++      tree->num_deleting_rows = 0;
++
+       new_node(tree, NULL);
+       return tree;
+ }
+@@ -278,6 +288,7 @@ void ios_revive_node(glp_tree *tree, int p)
+             double *val;
+             ind = xcalloc(1+n, sizeof(int));
+             val = xcalloc(1+n, sizeof(double));
++            /* maintains the row order during revival */
+             for (r = node->r_ptr; r != NULL; r = r->next)
+             {  i = glp_add_rows(mip, 1);
+                glp_set_row_name(mip, i, r->name);
+@@ -457,6 +468,13 @@ void ios_freeze_node(glp_tree *tree)
+             double *val;
+             ind = xcalloc(1+n, sizeof(int));
+             val = xcalloc(1+n, sizeof(double));
++            /* Added rows are stored in the same order!
++             * This is done by 2 reversals.
++             * - Going from i = m down pred_m.
++             * - Storing the list by a stack "push"
++             *   node->r_ptr = r;
++             *   r->next = node->r_ptr;
++             */
+             for (i = m; i > pred_m; i--)
+             {  GLPROW *row = mip->row[i];
+                IOSROW *r;
+@@ -501,6 +519,8 @@ void ios_freeze_node(glp_tree *tree)
+             xassert(nrs > 0);
+             num = xcalloc(1+nrs, sizeof(int));
+             for (i = 1; i <= nrs; i++) num[i] = root_m + i;
++            /* To not call ios_cb_rows_deleted here.
++               These rows have been saved earlier. */
+             glp_del_rows(mip, nrs, num);
+             xfree(num);
+          }
+@@ -636,6 +656,9 @@ static IOSNPD *new_node(glp_tree *tree, IOSNPD *parent)
+       tree->a_cnt++;
+       tree->n_cnt++;
+       tree->t_cnt++;
++
++      node->ord = tree->t_cnt;
++
+       /* increase the number of child subproblems */
+       if (parent == NULL)
+          xassert(p == 1);
+@@ -818,6 +841,8 @@ void ios_delete_tree(glp_tree *tree)
+          xassert(nrs > 0);
+          num = xcalloc(1+nrs, sizeof(int));
+          for (i = 1; i <= nrs; i++) num[i] = tree->orig_m + i;
++         /* Do not call ios_cb_rows_deleted here.
++            This does not help log information. */
+          glp_del_rows(mip, nrs, num);
+          xfree(num);
+       }
+@@ -1417,6 +1442,7 @@ int ios_add_row(glp_tree *tree, IOSPOOL *pool,
+       cut->rhs = rhs;
+       cut->prev = pool->tail;
+       cut->next = NULL;
++      cut->aux = NULL;
+       if (cut->prev == NULL)
+          pool->head = cut;
+       else
+@@ -1517,6 +1543,11 @@ void ios_del_row(glp_tree *tree, IOSPOOL *pool, int i)
+          cut->ptr = aij->next;
+          dmp_free_atom(tree->pool, aij, sizeof(IOSAIJ));
+       }
++
++      if (cut->aux != NULL){
++        ios_delete_aux(cut->aux);
++      }
++
+       dmp_free_atom(tree->pool, cut, sizeof(IOSCUT));
+       pool->size--;
+       return;
+@@ -1535,6 +1566,10 @@ void ios_clear_pool(glp_tree *tree, IOSPOOL *pool)
+             cut->ptr = aij->next;
+             dmp_free_atom(tree->pool, aij, sizeof(IOSAIJ));
+          }
++
++         if (cut->aux != NULL){
++           ios_delete_aux(cut->aux);
++         }
+          dmp_free_atom(tree->pool, cut, sizeof(IOSCUT));
+       }
+       pool->size = 0;
+diff --git a/src/glpios03.c b/src/glpios03.c
+index 80d701b..03e9208 100644
+--- a/src/glpios03.c
++++ b/src/glpios03.c
+@@ -358,12 +358,12 @@ static void fix_by_red_cost(glp_tree *T)
+ *  2 - both branches are hopeless and have been pruned; new subproblem
+ *      selection is needed to continue the search. */
+ 
+-static int branch_on(glp_tree *T, int j, int next)
++static int branch_on(glp_tree *T, int j, int next, int clone[], int* to_up)
+ {     glp_prob *mip = T->mip;
+       IOSNPD *node;
+       int m = mip->m;
+       int n = mip->n;
+-      int type, dn_type, up_type, dn_bad, up_bad, p, ret, clone[1+2];
++      int type, dn_type, up_type, dn_bad, up_bad, p, ret;
+       double lb, ub, beta, new_ub, new_lb, dn_lp, up_lp, dn_bnd, up_bnd;
+       /* determine bounds and value of x[j] in optimal solution to LP
+          relaxation of the current subproblem */
+@@ -431,6 +431,7 @@ static int branch_on(glp_tree *T, int j, int next)
+          else
+             xassert(mip != mip);
+          ret = 1;
++         *to_up = 0; /* up is bad. Do not go to up. */
+          goto done;
+       }
+       else if (dn_bad)
+@@ -449,6 +450,7 @@ static int branch_on(glp_tree *T, int j, int next)
+          else
+             xassert(mip != mip);
+          ret = 1;
++         *to_up = 1; /* down is bad. Go to up. */
+          goto done;
+       }
+       /* both down- and up-branches seem to be hopeful */
+@@ -702,7 +704,8 @@ static void remove_cuts(glp_tree *T)
+          }
+       }
+       if (cnt > 0)
+-      {  glp_del_rows(T->mip, cnt, num);
++      {  ios_cb_rows_deleted(T, cnt, num);
++         glp_del_rows(T->mip, cnt, num);
+ #if 0
+          xprintf("%d inactive cut(s) removed\n", cnt);
+ #endif
+@@ -784,6 +787,7 @@ static void display_cut_info(glp_tree *T)
+ 
+ int ios_driver(glp_tree *T)
+ {     int p, curr_p, p_stat, d_stat, ret;
++      int branch_clones[1 + 2];
+ #if 1 /* carry out to glp_tree */
+       int pred_p = 0;
+       /* if the current subproblem has been just created due to
+@@ -1010,6 +1014,12 @@ more: /* minor loop starts here */
+          if (T->parm->msg_lev >= GLP_MSG_DBG)
+             xprintf("LP relaxation has no feasible solution\n");
+          /* prune the branch */
++         if (T->parm->cb_func != NULL)
++         {  xassert(T->reason == 0);
++            T->reason = GLP_LI_CLOSE;
++            T->parm->cb_func(T, T->parm->cb_info);
++            T->reason = 0;
++         }
+          goto fath;
+       }
+       else
+@@ -1219,6 +1229,14 @@ more: /* minor loop starts here */
+          ios_process_cuts(T);
+          T->reason = 0;
+       }
++      /* if the local cut pool is not empty and the callback func is there,
++         this gives the callback the chance to see what was selected. */
++      if (T->parm->cb_func != NULL && T->local->size > 0)
++      {  xassert(T->reason == 0);
++         T->reason = GLP_ICUTSELECT;
++         T->parm->cb_func(T, T->parm->cb_info);
++         T->reason = 0;
++      }
+       /* clear the local cut pool */
+       ios_clear_pool(T, T->local);
+       /* perform re-optimization, if necessary */
+@@ -1255,7 +1273,34 @@ more: /* minor loop starts here */
+          T->br_var = ios_choose_var(T, &T->br_sel);
+       /* perform actual branching */
+       curr_p = T->curr->p;
+-      ret = branch_on(T, T->br_var, T->br_sel);
++      xassert(T->br_to_up == 0);
++      ret = branch_on(T, T->br_var, T->br_sel, branch_clones, &T->br_to_up);
++      if (T->parm->cb_func != NULL)
++      {  xassert(T->reason == 0);
++         xassert(T->br_node == 0);
++         xassert(T->dn_child == 0);
++         xassert(T->up_child == 0);
++         // record a branch here
++         T->reason = GLP_LI_BRANCH;
++         // at this point T->br_var is the branching variable
++         T->br_node = curr_p;
++         T->br_result = ret;
++         if(ret == 0){
++           T->dn_child = branch_clones[1];
++           T->up_child = branch_clones[2];
++         }
++         T->parm->cb_func(T, T->parm->cb_info);
++         T->reason = 0;
++         T->br_node = 0;
++         T->dn_child = 0;
++         T->up_child = 0;
++         if (T->stop)
++         {  ret = GLP_ESTOP;
++            goto done;
++         }
++      }
++      T->br_to_up = 0;
++
+       T->br_var = T->br_sel = 0;
+       if (ret == 0)
+       {  /* both branches have been created */
+diff --git a/src/glpios05.c b/src/glpios05.c
+index b9322b9..caf69d7 100644
+--- a/src/glpios05.c
++++ b/src/glpios05.c
+@@ -65,6 +65,9 @@ static void gen_cut(glp_tree *tree, struct worka *worka, int j)
+       double *phi = worka->phi;
+       int i, k, len, kind, stat;
+       double lb, ub, alfa, beta, ksi, phi1, rhs;
++      int input_j;
++      input_j = j;
++
+       /* compute row of the simplex tableau, which (row) corresponds
+          to specified basic variable xB[i] = x[m+j]; see (23) */
+       len = glp_eval_tab_row(mip, m+j, ind, val);
+@@ -73,6 +76,7 @@ static void gen_cut(glp_tree *tree, struct worka *worka, int j)
+          if it would be computed with formula (27); it is assumed that
+          beta[i] is fractional enough */
+       beta = mip->col[j]->prim;
++
+       /* compute cut coefficients phi and right-hand side rho, which
+          correspond to formula (30); dense format is used, because rows
+          of the simplex tableau is usually dense */
+@@ -104,6 +108,7 @@ static void gen_cut(glp_tree *tree, struct worka *worka, int j)
+          xassert(stat != GLP_BS);
+          /* determine row coefficient ksi[i,j] at xN[j]; see (23) */
+          ksi = val[j];
++         /* printf("%d %4.15f ", k, ksi); */
+          /* if ksi[i,j] is too large in the magnitude, do not generate
+             the cut */
+          if (fabs(ksi) > 1e+05) goto fini;
+@@ -135,6 +140,7 @@ static void gen_cut(glp_tree *tree, struct worka *worka, int j)
+                /* y[j] is integer */
+                if (fabs(alfa - floor(alfa + 0.5)) < 1e-10)
+                {  /* alfa[i,j] is close to nearest integer; skip it */
++                  /* printf("(skip)"); */
+                   goto skip;
+                }
+                else if (f(alfa) <= f(beta))
+@@ -170,6 +176,13 @@ static void gen_cut(glp_tree *tree, struct worka *worka, int j)
+          }
+ skip:    ;
+       }
++      /* printf("\n"); */
++      /* for (i = 1; i <= m+n; i++) */
++      /* { */
++      /*   printf("%i %f, ", i, phi[i]); */
++      /* } */
++      /* printf("\n"); */
++
+       /* now the cut has the form sum_k phi[k] * x[k] >= rho, where cut
+          coefficients are stored in the array phi in dense format;
+          x[1,...,m] are auxiliary variables, x[m+1,...,m+n] are struc-
+@@ -218,8 +231,20 @@ skip:    ;
+       ios_add_cut_row(tree, pool, GLP_RF_GMI, len, ind, val, GLP_LO,
+          rhs);
+ #else
+-      glp_ios_add_row(tree, NULL, GLP_RF_GMI, 0, len, ind, val, GLP_LO,
+-         rhs);
++      int ord;
++      ord = glp_ios_add_row(tree, NULL, GLP_RF_GMI, 0, len, ind, val,
++                            GLP_LO, rhs);
++      ios_cut_set_single_aux(tree, ord, input_j);
++
++      /* printf("ord: % d beta %f\n", ord, beta); */
++
++      /** callback for a cut being added to the cut pool */
++      if (tree->parm->cb_func != NULL)
++      {  xassert(tree->reason == GLP_ICUTGEN);
++         tree->reason = GLP_ICUTADDED;
++         tree->parm->cb_func(tree, tree->parm->cb_info);
++         tree->reason = GLP_ICUTGEN;
++      }
+ #endif
+ fini: return;
+ }
+diff --git a/src/glpios06.c b/src/glpios06.c
+index 2bd0453..4bd3cf5 100644
+--- a/src/glpios06.c
++++ b/src/glpios06.c
+@@ -100,6 +100,29 @@ struct MIR
+       /* sparse vector of cutting plane coefficients, alpha[k] */
+       double cut_rhs;
+       /* right-hand size of the cutting plane, beta */
++
++      /*-------------------------------------------------------------*/
++      /* Extras I've added to reproduce a cut externally */
++      double cut_delta;
++      /* the delta used for generating the cut */
++      char *cut_cset; /* char cut_vec[1+m+n]; */
++      /* cut_cset[k], 1 <= k <= m+n, is set to true if structural
++         variable x[k] was complemented in the cut:
++         0 - x[k] has been not been complemented
++         non 0 - x[k] has been complemented */
++      int *vlb_rows; /* int vlb_rows[1+m+n]; */
++      /* vlb_rows[k], 1 <= k <= m+n,
++       * vlb_rows[k] <= 0 if virtual lower bound has not been set
++       * vlb_rows[k] = r if virtual lower bound was set using the row r
++       */
++      int *vub_rows; /* int vub_rows[1+m+n]; */
++      /* vub_rows[k], 1 <= k <= m+n,
++       * vub_rows[k] <= 0 if virtual upper bound has not been set
++       * vub_rows[k] = r if virtual upper bound was set using the row r
++       */
++
++      double *agg_coeffs;
++      /* coefficients used to multiply agg_coeffs */
+ };
+ 
+ /***********************************************************************
+@@ -146,6 +169,7 @@ static void set_row_attrib(glp_tree *tree, struct MIR *mir)
+                xassert(row != row);
+          }
+          mir->vlb[k] = mir->vub[k] = 0;
++         mir->vlb_rows[k] = mir->vub_rows[k] = 0;
+       }
+       return;
+ }
+@@ -181,6 +205,7 @@ static void set_col_attrib(glp_tree *tree, struct MIR *mir)
+                xassert(col != col);
+          }
+          mir->vlb[k] = mir->vub[k] = 0;
++         mir->vlb_rows[k] = mir->vub_rows[k] = 0;
+       }
+       return;
+ }
+@@ -230,6 +255,7 @@ static void set_var_bounds(glp_tree *tree, struct MIR *mir)
+             {  /* set variable lower bound for x1 */
+                mir->lb[k1] = - a2 / a1;
+                mir->vlb[k1] = k2;
++               mir->vlb_rows[k1]  = i;
+                /* the row should not be used */
+                mir->skip[i] = 1;
+             }
+@@ -240,6 +266,7 @@ static void set_var_bounds(glp_tree *tree, struct MIR *mir)
+             {  /* set variable upper bound for x1 */
+                mir->ub[k1] = - a2 / a1;
+                mir->vub[k1] = k2;
++               mir->vub_rows[k1]  = i;
+                /* the row should not be used */
+                mir->skip[i] = 1;
+             }
+@@ -313,6 +340,13 @@ void *ios_mir_init(glp_tree *tree)
+       mir->subst = xcalloc(1+m+n, sizeof(char));
+       mir->mod_vec = ios_create_vec(m+n);
+       mir->cut_vec = ios_create_vec(m+n);
++
++      /* added */
++      mir->cut_cset = xcalloc(1+m+n, sizeof(char));
++      mir->vlb_rows = xcalloc(1+m+n, sizeof(int));
++      mir->vub_rows = xcalloc(1+m+n, sizeof(int));
++      mir->agg_coeffs = xcalloc(1+MAXAGGR, sizeof(double));
++
+       /* set global row attributes */
+       set_row_attrib(tree, mir);
+       /* set global column attributes */
+@@ -405,6 +439,9 @@ static void initial_agg_row(glp_tree *tree, struct MIR *mir, int i)
+       mir->skip[i] = 2;
+       mir->agg_cnt = 1;
+       mir->agg_row[1] = i;
++
++      mir->agg_coeffs[1] = +1.0;
++
+       /* use x[i] - sum a[i,j] * x[m+j] = 0, where x[i] is auxiliary
+          variable of row i, x[m+j] are structural variables */
+       ios_clear_vec(mir->agg_vec);
+@@ -784,19 +821,19 @@ static int cmir_cmp(const void *p1, const void *p2)
+ 
+ static double cmir_sep(const int n, const double a[], const double b,
+       const double u[], const double x[], const double s,
+-      double alpha[], double *beta, double *gamma)
++      double alpha[], double *beta, double *gamma,
++      double* delta, char cset[])
+ {     int fail, j, k, nv, v;
+-      double delta, eps, d_try[1+3], r, r_best;
+-      char *cset;
++      double eps, d_try[1+3], r, r_best;
+       struct vset *vset;
+       /* allocate working arrays */
+-      cset = xcalloc(1+n, sizeof(char));
++      //cset = xcalloc(1+n, sizeof(char));
+       vset = xcalloc(1+n, sizeof(struct vset));
+       /* choose initial C */
+       for (j = 1; j <= n; j++)
+          cset[j] = (char)(x[j] >= 0.5 * u[j]);
+       /* choose initial delta */
+-      r_best = delta = 0.0;
++      r_best = (*delta) = 0.0;
+       for (j = 1; j <= n; j++)
+       {  xassert(a[j] != 0.0);
+          /* if x[j] is close to its bounds, skip it */
+@@ -809,16 +846,16 @@ static double cmir_sep(const int n, const double a[], const double b,
+          /* compute violation */
+          r = - (*beta) - (*gamma) * s;
+          for (k = 1; k <= n; k++) r += alpha[k] * x[k];
+-         if (r_best < r) r_best = r, delta = fabs(a[j]);
++         if (r_best < r) r_best = r, (*delta) = fabs(a[j]);
+       }
+       if (r_best < 0.001) r_best = 0.0;
+       if (r_best == 0.0) goto done;
+-      xassert(delta > 0.0);
++      xassert((*delta) > 0.0);
+       /* try to increase violation by dividing delta by 2, 4, and 8,
+          respectively */
+-      d_try[1] = delta / 2.0;
+-      d_try[2] = delta / 4.0;
+-      d_try[3] = delta / 8.0;
++      d_try[1] = (*delta) / 2.0;
++      d_try[2] = (*delta) / 4.0;
++      d_try[3] = (*delta) / 8.0;
+       for (j = 1; j <= 3; j++)
+       {  /* construct c-MIR inequality */
+          fail = cmir_ineq(n, a, b, u, cset, d_try[j], alpha, beta,
+@@ -827,7 +864,7 @@ static double cmir_sep(const int n, const double a[], const double b,
+          /* compute violation */
+          r = - (*beta) - (*gamma) * s;
+          for (k = 1; k <= n; k++) r += alpha[k] * x[k];
+-         if (r_best < r) r_best = r, delta = d_try[j];
++         if (r_best < r) r_best = r, (*delta) = d_try[j];
+       }
+       /* build subset of variables lying strictly between their bounds
+          and order it by nondecreasing values of |x[j] - u[j]/2| */
+@@ -849,7 +886,7 @@ static double cmir_sep(const int n, const double a[], const double b,
+          /* replace x[j] by its complement or vice versa */
+          cset[j] = (char)!cset[j];
+          /* construct c-MIR inequality */
+-         fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
++         fail = cmir_ineq(n, a, b, u, cset, (*delta), alpha, beta, gamma);
+          /* restore the variable */
+          cset[j] = (char)!cset[j];
+          /* do not replace the variable in case of failure */
+@@ -860,10 +897,11 @@ static double cmir_sep(const int n, const double a[], const double b,
+          if (r_best < r) r_best = r, cset[j] = (char)!cset[j];
+       }
+       /* construct the best c-MIR inequality chosen */
+-      fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
++      fail = cmir_ineq(n, a, b, u, cset, (*delta), alpha, beta, gamma);
+       xassert(!fail);
++
+ done: /* free working arrays */
+-      xfree(cset);
++
+       xfree(vset);
+       /* return to the calling routine */
+       return r_best;
+@@ -874,7 +912,8 @@ static double generate(struct MIR *mir)
+       int m = mir->m;
+       int n = mir->n;
+       int j, k, kk, nint;
+-      double s, *u, *x, *alpha, r_best = 0.0, b, beta, gamma;
++      double s, *u, *x, *alpha, r_best = 0.0, b, beta, gamma, delta;
++      char *cset;
+       ios_copy_vec(mir->cut_vec, mir->mod_vec);
+       mir->cut_rhs = mir->mod_rhs;
+       /* remove small terms, which can appear due to substitution of
+@@ -923,6 +962,7 @@ static double generate(struct MIR *mir)
+       u = xcalloc(1+nint, sizeof(double));
+       x = xcalloc(1+nint, sizeof(double));
+       alpha = xcalloc(1+nint, sizeof(double));
++      cset = xcalloc(1+nint, sizeof(char));
+       /* determine u and x */
+       for (j = 1; j <= nint; j++)
+       {  k = mir->cut_vec->ind[j];
+@@ -936,6 +976,7 @@ static double generate(struct MIR *mir)
+             x[j] = mir->ub[k] - mir->x[k];
+          else
+             xassert(k != k);
++         if(!(x[j] >= -0.001)) { goto skip; }
+          xassert(x[j] >= -0.001);
+          if (x[j] < 0.0) x[j] = 0.0;
+       }
+@@ -965,6 +1006,7 @@ static double generate(struct MIR *mir)
+          }
+          else
+             xassert(k != k);
++         if(!(x >= -0.001)) { goto skip; }
+          xassert(x >= -0.001);
+          if (x < 0.0) x = 0.0;
+          s -= mir->cut_vec->val[j] * x;
+@@ -973,7 +1015,7 @@ static double generate(struct MIR *mir)
+       /* apply heuristic to obtain most violated c-MIR inequality */
+       b = mir->cut_rhs;
+       r_best = cmir_sep(nint, mir->cut_vec->val, b, u, x, s, alpha,
+-         &beta, &gamma);
++         &beta, &gamma, &delta, cset);
+       if (r_best == 0.0) goto skip;
+       xassert(r_best > 0.0);
+       /* convert to raw cut */
+@@ -988,10 +1030,22 @@ static double generate(struct MIR *mir)
+ #if _MIR_DEBUG
+       ios_check_vec(mir->cut_vec);
+ #endif
++      /* added */
++      mir->cut_delta = delta;
++      // this is not a great place for resetting the array,
++      // but it should be sufficient
++      for (j = 1; j <= n+m;  j++)
++         mir->cut_cset[j] = 0;
++      for (j = 1; j <= nint; j++)
++      {  k = mir->cut_vec->ind[j];
++         xassert(m <= k  && k <= m+n);
++         mir->cut_cset[k] = cset[j];
++      }
+ skip: /* free working arrays */
+       xfree(u);
+       xfree(x);
+       xfree(alpha);
++      xfree(cset);
+ done: return r_best;
+ }
+ 
+@@ -1199,8 +1253,25 @@ static void add_cut(glp_tree *tree, struct MIR *mir)
+       ios_add_cut_row(tree, pool, GLP_RF_MIR, len, ind, val, GLP_UP,
+          mir->cut_rhs);
+ #else
+-      glp_ios_add_row(tree, NULL, GLP_RF_MIR, 0, len, ind, val, GLP_UP,
++      int ord;
++      ord = glp_ios_add_row(tree, NULL, GLP_RF_MIR, 0, len, ind, val, GLP_UP,
+          mir->cut_rhs);
++      ios_cut_set_aux(tree, ord, mir->agg_cnt, mir->agg_row, mir->agg_coeffs);
++      ios_cut_set_aux_mir(tree, ord, mir->cut_delta,
++                          mir->cut_cset, mir->subst,
++                          mir->vlb_rows, mir->vub_rows);
++
++      /** callback for a cut being added to the cut pool */
++      /* printf("mir tree parm %p %d\n", tree->parm->cb_func, ord); */
++      /* printf("  agg_rhs %f\n", mir->agg_rhs); */
++      /* printf("  mod_rhs %f\n", mir->mod_rhs); */
++      /* printf("  cut_rhs %f\n", mir->cut_rhs); */
++      if (tree->parm->cb_func != NULL)
++      {  xassert(tree->reason == GLP_ICUTGEN);
++         tree->reason = GLP_ICUTADDED;
++         tree->parm->cb_func(tree, tree->parm->cb_info);
++         tree->reason = GLP_ICUTGEN;
++      }
+ #endif
+       xfree(ind);
+       xfree(val);
+@@ -1216,6 +1287,7 @@ static int aggregate_row(glp_tree *tree, struct MIR *mir)
+       IOSVEC *v;
+       int ii, j, jj, k, kk, kappa = 0, ret = 0;
+       double d1, d2, d, d_max = 0.0;
++      double guass_coeff;
+       /* choose appropriate structural variable in the aggregated row
+          to be substituted */
+       for (j = 1; j <= mir->agg_vec->nnz; j++)
+@@ -1300,8 +1372,9 @@ static int aggregate_row(glp_tree *tree, struct MIR *mir)
+       xassert(j != 0);
+       jj = v->pos[kappa];
+       xassert(jj != 0);
+-      ios_linear_comb(mir->agg_vec,
+-         - mir->agg_vec->val[j] / v->val[jj], v);
++      guass_coeff = - mir->agg_vec->val[j] / v->val[jj];
++      mir->agg_coeffs[mir->agg_cnt] = guass_coeff;
++      ios_linear_comb(mir->agg_vec, guass_coeff, v);
+       ios_delete_vec(v);
+       ios_set_vj(mir->agg_vec, kappa, 0.0);
+ #if _MIR_DEBUG
+@@ -1440,6 +1513,13 @@ void ios_mir_term(void *gen)
+       xfree(mir->subst);
+       ios_delete_vec(mir->mod_vec);
+       ios_delete_vec(mir->cut_vec);
++
++      /* added */
++      xfree(mir->cut_cset);
++      xfree(mir->vlb_rows);
++      xfree(mir->vub_rows);
++      xfree(mir->agg_coeffs);
++
+       xfree(mir);
+       return;
+ }
+diff --git a/src/glpios07.c b/src/glpios07.c
+index 0892550..9f742e6 100644
+--- a/src/glpios07.c
++++ b/src/glpios07.c
+@@ -543,6 +543,14 @@ void ios_cov_gen(glp_tree *tree)
+          /* add the cut to the cut pool */
+          glp_ios_add_row(tree, NULL, GLP_RF_COV, 0, len, ind, val,
+             GLP_UP, val[0]);
++
++         /** callback for a cut being added to the cut pool */
++         if (tree->parm->cb_func != NULL)
++         {  xassert(tree->reason == GLP_ICUTGEN);
++            tree->reason = GLP_ICUTADDED;
++            tree->parm->cb_func(tree, tree->parm->cb_info);
++            tree->reason = GLP_ICUTGEN;
++         }
+       }
+       /* free working arrays */
+       xfree(ind);
+diff --git a/src/glpios08.c b/src/glpios08.c
+index 2d2fcd5..3aad808 100644
+--- a/src/glpios08.c
++++ b/src/glpios08.c
+@@ -129,6 +129,15 @@ void ios_clq_gen(glp_tree *T, void *G_)
+       /* add cut inequality to local cut pool */
+       glp_ios_add_row(T, NULL, GLP_RF_CLQ, 0, len, ind, val, GLP_UP,
+          rhs);
++
++      /** callback for a cut being added to the cut pool */
++      if (T->parm->cb_func != NULL)
++      {  xassert(T->reason == GLP_ICUTGEN);
++         T->reason = GLP_ICUTADDED;
++         T->parm->cb_func(T, T->parm->cb_info);
++         T->reason = GLP_ICUTGEN;
++      }
++
+ skip: /* free working arrays */
+       tfree(ind);
+       tfree(val);
+diff --git a/src/glpios11.c b/src/glpios11.c
+index c40e9a5..82d5698 100644
+--- a/src/glpios11.c
++++ b/src/glpios11.c
+@@ -193,6 +193,9 @@ void ios_process_cuts(glp_tree *T)
+          glp_set_mat_row(T->mip, i, len, ind, val);
+          xassert(cut->type == GLP_LO || cut->type == GLP_UP);
+          glp_set_row_bnds(T->mip, i, cut->type, cut->rhs, cut->rhs);
++
++         /* setting this as selected */
++         ios_cut_set_selected(cut, i);
+       }
+       /* free working arrays */
+       xfree(info);
+diff --git a/src/glpk.h b/src/glpk.h
+index 75c292c..e117782 100644
+--- a/src/glpk.h
++++ b/src/glpk.h
+@@ -130,6 +130,7 @@ typedef struct
+       int out_frq;            /* spx.out_frq */
+       int out_dly;            /* spx.out_dly (milliseconds) */
+       int presolve;           /* enable/disable using LP presolver */
++      int stability_lmt;      /* maximum number of check stability failures before stopping */
+       double foo_bar[36];     /* (reserved) */
+ } glp_smcp;
+ 
+@@ -229,6 +230,11 @@ typedef struct
+ #define GLP_IBRANCH     0x05  /* request for branching */
+ #define GLP_ISELECT     0x06  /* request for subproblem selection */
+ #define GLP_IPREPRO     0x07  /* request for preprocessing */
++#define GLP_ICUTADDED   0x08  /* cut was added to the pool */
++#define GLP_ICUTSELECT  0x09  /* cuts were selected as rows */
++#define GLP_LI_BRANCH   0x10  /* a branch was made */
++#define GLP_LI_CLOSE    0x11  /* an active node was closed */
++#define GLP_LI_DELROW   0x12  /* an active node was closed */
+ 
+ /* branch selection indicator: */
+ #define GLP_NO_BRNCH       0  /* select no branch */
+@@ -1057,6 +1063,54 @@ int glp_top_sort(glp_graph *G, int v_num);
+ int glp_wclique_exact(glp_graph *G, int v_wgt, double *sol, int v_set);
+ /* find maximum weight clique with exact algorithm */
+ 
++/*******************************************/
++/***               CUT LOG               ***/
++/*******************************************/
++
++int glp_get_it_cnt(glp_prob *P);
++/* get the iteration count of the current problem */
++
++
++int glp_ios_get_cut(glp_tree *T, int i, int ind[], double val[], int* klass, int* type, double* rhs);
++/* determine reason for calling the callback routine */
++
++int glp_ios_cut_get_aux_nrows(glp_tree *tree, int ord);
++/* gets the number of rows used to generate a cut. */
++
++
++void glp_ios_cut_get_aux_rows(glp_tree *tree, int ord,
++                              int rows[], double coeffs[]);
++/* gets a cut as an input sequence of rows times coefficients. */
++
++
++void glp_ios_cut_get_mir_cset(glp_tree *tree, int ord, char cset[]);
++/* gets mir cut complement set. */
++double glp_ios_cut_get_mir_delta(glp_tree *tree, int ord);
++/* gets mir cut delta. */
++void glp_ios_cut_get_mir_subst(glp_tree *tree, int ord, char subst[]);
++/* gets mir cut substition information. */
++void glp_ios_cut_get_mir_virtual_rows(glp_tree *tree, int ord,
++                                      int vlb_rows[], int vub_rows[]);
++/* gets mir cut virtual bounds rows. */
++
++int glp_ios_selected_cuts(glp_tree *tree, int ords[], int sel[]);
++/* gets the list of selected cuts.
++   Can only be called when GLP_ICUTSELECT */
++
++
++int glp_ios_branch_log(glp_tree *tree, double *val, int* parent, int* dn, int* up);
++/* can only be called when GLP_LI_BRANCH.
++ * If id is non-null, returns the id of the structural variable branched upon.
++ * If val is non-null, it is set to the value branched upon.
++ * If parent is non-null, it is set to node id of the node branched upon.
++ * If dn is non-null, it is set to node id of the newly created down node.
++ * If up is non-null, it is set to node id of the newly created up node.
++ */
++
++int glp_ios_node_ord(glp_tree *tree, int node_p);
++
++int glp_ios_rows_deleted(glp_tree *tree, int* rows);
++
+ #ifdef __cplusplus
+ }
+ #endif
+diff --git a/src/glpspx01.c b/src/glpspx01.c
+index 5c17114..caaab27 100644
+--- a/src/glpspx01.c
++++ b/src/glpspx01.c
+@@ -238,6 +238,9 @@ struct csa
+       double *work2; /* double work2[1+m]; */
+       double *work3; /* double work3[1+m]; */
+       double *work4; /* double work4[1+m]; */
++
++      /** Things Tim has added. */
++      int stability_failures;
+ };
+ 
+ static const double kappa = 0.10;
+@@ -400,6 +403,8 @@ static void init_csa(struct csa *csa, glp_prob *lp)
+       csa->refct = 0;
+       memset(&refsp[1], 0, (m+n) * sizeof(char));
+       for (j = 1; j <= n; j++) gamma[j] = 1.0;
++
++      csa->stability_failures = 0;
+       return;
+ }
+ 
+@@ -2647,6 +2652,11 @@ loop: /* main loop starts here */
+          if (check_stab(csa, parm->tol_bnd))
+          {  /* there are excessive bound violations due to round-off
+                errors */
++            csa->stability_failures++;
++            if (csa->stability_failures >= parm->stability_lmt){
++              ret = GLP_EINSTAB;
++              goto done;
++            }
+             if (parm->msg_lev >= GLP_MSG_ERR)
+                xprintf("Warning: numerical instability (primal simplex,"
+                   " phase %s)\n", csa->phase == 1 ? "I" : "II");
+diff --git a/src/glpspx02.c b/src/glpspx02.c
+index 19152a6..de3d677 100644
+--- a/src/glpspx02.c
++++ b/src/glpspx02.c
+@@ -288,6 +288,9 @@ struct csa
+       double *work2; /* double work2[1+m]; */
+       double *work3; /* double work3[1+m]; */
+       double *work4; /* double work4[1+m]; */
++
++      /* count the number of stability failures */
++      int stability_failures;
+ };
+ 
+ static const double kappa = 0.10;
+@@ -501,6 +504,8 @@ static void init_csa(struct csa *csa, glp_prob *lp)
+       csa->refct = 0;
+       memset(&refsp[1], 0, (m+n) * sizeof(char));
+       for (i = 1; i <= m; i++) gamma[i] = 1.0;
++
++      csa->stability_failures = 0;
+       return;
+ }
+ 
+@@ -2741,7 +2746,13 @@ loop: /* main loop starts here */
+          /* make sure that the current basic solution remains dual
+             feasible */
+          if (check_stab(csa, parm->tol_dj) != 0)
+-         {  if (parm->msg_lev >= GLP_MSG_ERR)
++         {
++            csa->stability_failures++;
++            if (csa->stability_failures >= parm->stability_lmt){
++              ret = GLP_EINSTAB;
++              goto done;
++            }
++            if (parm->msg_lev >= GLP_MSG_ERR)
+                xprintf("Warning: numerical instability (dual simplex, p"
+                   "hase %s)\n", csa->phase == 1 ? "I" : "II");
+ #if 1
+@@ -3023,6 +3034,10 @@ loop: /* main loop starts here */
+       /* accuracy check based on the pivot element */
+       {  double piv1 = csa->tcol_vec[csa->p]; /* more accurate */
+          double piv2 = csa->trow_vec[csa->q]; /* less accurate */
++         if(piv1 == 0.0){
++           ret = GLP_EFAIL;
++           goto done;
++         }
+          xassert(piv1 != 0.0);
+          if (fabs(piv1 - piv2) > 1e-8 * (1.0 + fabs(piv1)) ||
+              !(piv1 > 0.0 && piv2 > 0.0 || piv1 < 0.0 && piv2 < 0.0))
+-- 
+2.26.2