Repairing infeasibility

/***********************************************************************
   Xpress Optimizer Examples
   =========================

   file repair.c
   `````````````
   Demonstrates the FICO repairinfeas utility.

   Prints a relaxation summary and creates the relaxed subproblem.

   (c) 2017 Fair Isaac Corporation       
***********************************************************************/

#include 
#include 
#include 

#include "xprs.h"

void errormsg(const char *sSubName,int nLineNo,int nErrCode);

XPRSprob xprob;
void XPRS_CC Message(XPRSprob my_prob, void* my_object, const char *msg, int len, int msgtype);
int  XPRS_CC XPRSgetinfeasibilitybreakers(XPRSprob xprob, const double *x, const double *slacks, double constraints[], double bounds[], const double tolarence);
int  XPRS_CC XPRSgeneratefeasibilityreport(XPRSprob xprob, const double *constraintviolations, const double *boundviolations);
int  XPRS_CC XPRSgetrowbounds(XPRSprob xprob, const int row, double *lb, double *ub);
int  XPRS_CC XPRSaplyrepairinfeasresulttoproblem(XPRSprob xprob, const double *constraintviolations, const double *boundviolations);
int  XPRS_CC XPRSsetrowbounds(XPRSprob xprob, const int i, const double lb, const double ub);

int main(int argc, char **argv) {

  int i,j;
  int nReturn = 0;                    /* Return value of Optimizer subroutines */
  int Status;                         /* Status of Optimizer subroutines */
  int IIScount;                       /* Number of IISs found     */
  int nOptimizerVersion;              /* Optimizer version number */
  char *sProblem = NULL;              /* Problem name */
  int nExpiry;
  int nrows, ncols;
  char banner[256];
  double *x = NULL, *slacks = NULL;
  double *constraintviolations = NULL, *boundviolations = NULL;
  int nglent, nsets;

  double *lrp_array = NULL;   // Array of size ROWS containing the preferences for relaxing the less or equal side of row.  
  double *grp_array = NULL;   // Array of size ROWS containing the preferences for relaxing the greater or equal side of a row.  
  double *lbp_array = NULL;   // Array of size COLS containing the preferences for relaxing lower bounds.  
  double *ubp_array = NULL;   // Array of size COLS containing preferences for relaxing upper bounds.  

  /* Parse the arguments to the program. */
  if (argc != 2) {
    printf("syntax: repair  \n");
    return(1);
  }
  sProblem = argv[1];

  /* Initialize Optimizer */   
  nReturn = XPRSinit(NULL);
  XPRSgetbanner(banner); printf("%s\n",banner);
  if (nReturn == 8) {
    printf("Unable to initialize XPRESS\n");
    return(1);
  }

  /* Initialize problem object */
  nReturn = XPRScreateprob(&xprob);
  if (nReturn != 0 && nReturn != 32) {
    printf("Unable to create XPRESS problem\n");
    XPRSfree();
    return(1);
  }

  /* Get and display the Optimizer version number */
  if (nReturn = XPRSgetintcontrol(xprob, XPRS_VERSION, &nOptimizerVersion))
    errormsg("Unable to retrive optimizer version number",__LINE__,nReturn);

  printf("Xpress Optimizer Subroutine Library Release %.2f\n\n",
         (float)nOptimizerVersion/100);

  // forward messages to screen
  XPRSsetcbmessage(xprob,Message,NULL);

  /* Load problem file */
  if (nReturn = XPRSreadprob(xprob, sProblem, "")) 
    errormsg("Cannot input lp file",__LINE__,nReturn);  

  /* Get problem size */
  XPRSgetintattrib(xprob,XPRS_ORIGINALROWS,&nrows);
  XPRSgetintattrib(xprob,XPRS_ORIGINALCOLS,&ncols);

  // Allocate memory for the preference arrays
  lrp_array = malloc(nrows*sizeof(double)); if (!lrp_array) errormsg("Out of memory",__LINE__,-1);
  grp_array = malloc(nrows*sizeof(double)); if (!grp_array) errormsg("Out of memory",__LINE__,-1);
  lbp_array = malloc(ncols*sizeof(double)); if (!lbp_array) errormsg("Out of memory",__LINE__,-1);
  ubp_array = malloc(ncols*sizeof(double)); if (!ubp_array) errormsg("Out of memory",__LINE__,-1);

  // Allocate memory for the solution and the infeasibility breakers
  x                    = malloc(ncols*sizeof(double)); if (!x) errormsg("Out of memory",__LINE__,-1);
  slacks               = malloc(nrows*sizeof(double)); if (!slacks) errormsg("Out of memory",__LINE__,-1);  
  boundviolations      = malloc(ncols*sizeof(double)); if (!boundviolations) errormsg("Out of memory",__LINE__,-1);
  constraintviolations = malloc(nrows*sizeof(double)); if (!constraintviolations) errormsg("Out of memory",__LINE__,-1);

  // Set relaxation values
  for (i=0; i= XPRS_PLUSINFINITY) {
    rowtype = 'N';
    nReturn += XPRSchgrowtype(xprob, 1, &i, &rowtype);
  } else
  if (lb <= XPRS_MINUSINFINITY) {
    rowtype = 'L';
    nReturn += XPRSchgrowtype(xprob, 1, &i, &rowtype);
    nReturn += XPRSchgrhs(xprob, 1, &i, &ub);
  } else 
  if (ub >= XPRS_PLUSINFINITY) {
    rowtype = 'G';
    nReturn += XPRSchgrowtype(xprob, 1, &i, &rowtype);
    nReturn += XPRSchgrhs(xprob, 1, &i, &lb);
  } else 
  if (lb == ub) {
    rowtype = 'E';
    nReturn += XPRSchgrowtype(xprob, 1, &i, &rowtype);
    nReturn += XPRSchgrhs(xprob, 1, &i, &ub);
  } else {    
    rowtype = 'L';
    nReturn += XPRSchgrowtype(xprob, 1, &i, &rowtype);
    nReturn += XPRSchgrhs(xprob, 1, &i, &ub);
    range = ub-lb;
    nReturn += XPRSchgrhsrange(xprob, 1, &i, &ub);
  }

  return(nReturn);

}

/*
  This function calculates the infeasibilities of a solution (i.e. the values of the feasibility breakers after a repairinfeas call) 
*/
int XPRS_CC XPRSgetinfeasibilitybreakers(XPRSprob xprob, const double *x, const double *slacks, double constraints[], double bounds[], const double tolarence) {

  int i,j,ncols, nrows;  
  double rhs, activity;
  double lb, ub;

  // get problem dimensions
  XPRSgetintattrib(xprob,XPRS_ORIGINALROWS,&nrows);
  XPRSgetintattrib(xprob,XPRS_ORIGINALCOLS,&ncols);

  // check constraints
  for (i=0; i ub+tolarence) {
      constraints[i] = activity-ub;  // a positive value indicates that the upper bound is relaxed
    } else {
      constraints[i] = 0.0;
    }
  }

  // check bounds
  for (j=0; j ub+tolarence) {
      bounds[j] = x[j]-ub;  // a positive value indicates that the upper bound is relaxed
    } else {
      bounds[j] = 0.0;
    }
  }

  return(0);
}

// This function just converts a double to a string for reporting, converting infinity values as "NONE"
void Bound_to_string(double bound, char *string) {
  if (bound > XPRS_MINUSINFINITY && bound < XPRS_PLUSINFINITY) {
    if (bound >= 0)
      sprintf(string, " %.5e", bound);
    else
      sprintf(string, "%.5e", bound);
  } else {
    memcpy(string,"    NONE     \0",14);
  }
}

/*
  Prints a summary of the infeasibily breakers using the already calculated infeasibities
*/
int  XPRS_CC XPRSgeneratefeasibilityreport(XPRSprob xprob, const double *constraintviolations, const double *boundviolations) {
  
  int i,j;
  int nrows, ncols;
  double suminf = 0;                          /* Sum if infeasibility */
  int    ninf = 0;                            /* Number of rows and columns reapired */
  int    len, NameLength;                     /* Name lengths in the LP file */  
  char  *name = NULL;
  double lb,ub;
  char   slb[25], sub[25];                    /* Fixed length name buffers */
  char   slbshift[25], subshift[25];

  /* Get name lengths */
  XPRSgetintattrib(xprob,XPRS_ORIGINALROWS,&nrows);
  XPRSgetintattrib(xprob,XPRS_ORIGINALCOLS,&ncols);

  /* Get the maximum name length */
  NameLength = 0;
  for(i=0;iNameLength) NameLength=len;
  } 
  for(i=0;iNameLength) NameLength=len;
  } 

  /* Allocate name buffer */
  name = malloc(NameLength); 
  if (!name) errormsg("Out of memory",__LINE__,-1);

  printf("\n");
  printf("\n");
  printf("XPRESS-MP FEASIBILITY REPAIR REPORT.\n");
  printf("\n");
  printf("The following constraints were repaired.\n");
  printf("\n");
  
  printf("Index     Name%*s  Lower bound    Repaired        Upper bound    Repaired \n",NameLength-4,"");

  for (i=0; i 0) {
        ub += constraintviolations[i]; 
      } else {
        lb += constraintviolations[i];
      }

      XPRSsetrowbounds(xprob, i, lb, ub);
    }
  }

  // apply changes to columns
  for (j=0; j 0) {
        boundtype = 'U';
        ub += boundviolations[j];
        XPRSchgbounds(xprob, 1, &j, &boundtype, &ub);
      } else {
        boundtype = 'L';
        lb += boundviolations[j];
        XPRSchgbounds(xprob, 1, &j, &boundtype, &lb);
      }
    }
  }

  printf("\n");
  printf("Problem repaired.");
  printf("\n");

  return(0);

}

/* XPRS optimizer message callback */
void XPRS_CC Message(XPRSprob my_prob, void* my_object, 
                     const char *msg, int len, int msgtype)
{
  switch(msgtype)
  {
    case 4:  /* error */
    case 3:  /* warning */
    case 2:  /* not used */
    case 1:  /* information */
             printf("%s\n", msg);
             break;
    default: /* exiting - buffers need flushing */
             fflush(stdout);
             break;
  }
}

void errormsg(const char *sSubName,int nLineNo,int nErrCode)
{
   int nErrNo;   /* Optimizer error number */

   /* Print error message */
   printf("Error in line %d: %s\n",nLineNo,sSubName);

   /* Append the error code, if it exists */
   if (nErrCode!=-1) printf("with error code %d\n",nErrCode);

   /* Append Optimizer error number, if available */
   if (nErrCode==32) {
     XPRSgetintattrib(xprob,XPRS_ERRORCODE,&nErrNo);
     printf("The Optimizer error number is: %d\n",nErrNo);
   }

   /* Free memory, close files and exit */
   XPRSdestroyprob(xprob);
   XPRSfree();
   exit(nErrCode);
}