Merge branch 'master' into fixCMSBuildPRfixCMSBuildPR

1 files changed, 1528 insertions, 0 deletions

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