diff options
Diffstat (limited to 'contrib/glpk-cut-log.patch')
-rw-r--r-- | contrib/glpk-cut-log.patch | 1528 |
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 |