FICO
FICO Xpress Optimization Examples Repository
FICO Optimization Community FICO Xpress Optimization Home
Back to examples browserPrevious exampleNext example

Repairing infeasibility

Description
Demonstrates the repairinfeas utility, prints a relaxation summary and creates the relaxed subproblem.


Source Files





repair.c

/***********************************************************************
   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 <stdlib.h>
#include <stdio.h>
#include <math.h>

#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 Getinfeasibilitybreakers(XPRSprob xprob, const double *x, const double *slacks, double constraints[], double bounds[], const double tolarence);
int  XPRS_CC Generatefeasibilityreport(XPRSprob xprob, const double *constraintviolations, const double *boundviolations);
int  XPRS_CC Getrowbounds(XPRSprob xprob, const int row, double *lb, double *ub);
int  XPRS_CC Applyrepairinfeasresulttoproblem(XPRSprob xprob, const double *constraintviolations, const double *boundviolations);
int  XPRS_CC Setrowbounds(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 <filename> \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 read 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<nrows; i++) {
    lrp_array[i] = 1;
    grp_array[i] = 1;
  }
  for (j=0; j<ncols; j++) {
    lbp_array[j] = 1;
    ubp_array[j] = 1;
  }

  // Call repairinfeas
  Status = -1; // set a dummy value
  nReturn += XPRSrepairweightedinfeas(xprob, &Status, lrp_array, grp_array, lbp_array, ubp_array, 'n', 0.001, ""); if (nReturn != 0) { errormsg("Out of memory",__LINE__,1); }

  printf("Repairinfeas return code (%i): ", Status);
  switch(Status) {
    case 0: printf("relaxed optimum found\n"); break;
    case 1: printf("relaxed problem is infeasible\n"); break;
    case 2: printf("relaxed problem is unbounded\n"); break;
    case 3: printf("solution of the relaxed problem regarding the original objective is nonoptimal\n"); break;
    case 4: printf("error\n"); break;
    case 5: printf("numerical instability\n"); break;
  }

  XPRSgetglobal(xprob, &nglent, &nsets, NULL, NULL, NULL, NULL, NULL, NULL, NULL);

  // Get solution
  if ((nglent + nsets) == 0) {
    XPRSgetlpsol(xprob, x, slacks, NULL, NULL);
  } else {
    // NOTE: for MIPs, use the getmipsolution function instead
    XPRSgetmipsol(xprob, x, slacks);
  }

  // Get the values of the breaker variables
  Getinfeasibilitybreakers(xprob, x, slacks, constraintviolations, boundviolations, 1.0e-6);

  // Print report
  nReturn += Generatefeasibilityreport(xprob, constraintviolations, boundviolations); if (nReturn != 0) { errormsg("Out of memory",__LINE__,1); }

  // Repair and save problem
  nReturn += Applyrepairinfeasresulttoproblem(xprob, constraintviolations, boundviolations); if (nReturn != 0) { errormsg("Out of memory",__LINE__,1); }
  nReturn += XPRSwriteprob(xprob, "repaired.lp","lp"); if (nReturn != 0) { errormsg("Out of memory",__LINE__,1); }
  printf("Repaired problem is written as repaired.lp.\n");
  printf("(Please note that this matrix might show slightly different behaviour then the repairinfeas problem)\n\n");

  /* Free buffers */
  free(lrp_array);
  free(grp_array);
  free(lbp_array);
  free(ubp_array);

  free(x);
  free(slacks);
  free(boundviolations);
  free(constraintviolations);

  /* Destroy problem object and free environment */
  XPRSdestroyprob(xprob);
  XPRSfree();

  return 0;
}

/*
  This function returns the lower and upper bounds for a row.
*/
int XPRS_CC Getrowbounds(XPRSprob xprob, const int i, double *lb, double *ub) {

  char rowtype;
  double rhs, range;
  int nReturn = 0;

  // the bounds are calculated using the rhs, the range and the type of the row
  nReturn += XPRSgetrhs(xprob, &rhs, i, i);
  nReturn += XPRSgetrhsrange(xprob, &range, i, i);
  nReturn += XPRSgetrowtype(xprob, &rowtype, i, i);

  switch (rowtype) {
    case 'L': {
      *lb = XPRS_MINUSINFINITY;
      *ub = rhs;
    } break;
    case 'E': {
      *lb = rhs;
      *ub = rhs;
    } break;
    case 'G': {
      *lb = rhs;
      *ub = XPRS_PLUSINFINITY;
    } break;
    case 'R': {
      *lb = rhs-range;
      *ub = rhs;
    } break;
    default: {
      *lb = XPRS_MINUSINFINITY;
      *ub = XPRS_PLUSINFINITY;
    }
  }

  return(nReturn);

}

/*
  This function sets new lower and upper bounds for a row.
*/
int  XPRS_CC Setrowbounds(XPRSprob xprob, const int i, const double dlb, const double dub) {

  char rowtype;
  double lb, ub, range;
  int nReturn = 0;

  if (dlb < dub) {
    lb = dlb;
    ub = dub;
  } else {
    lb = dub;
    ub = dlb;
  }

  // determine row type and set bounds
  if (lb <= XPRS_MINUSINFINITY && ub >= 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 Getinfeasibilitybreakers(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<nrows; i++) {

    XPRSgetrhs(xprob, &rhs, i, i);
    activity = rhs - slacks[i];
    Getrowbounds(xprob, i, &lb, &ub);

    if (activity < lb-tolarence) {
      constraints[i] = activity-lb;  // a negative value indicates that the lower bound is relaxed
    } else
    if (activity > 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<ncols; j++) {
    XPRSgetlb(xprob, &lb, j, j);
    XPRSgetub(xprob, &ub, j, j);

    if (x[j] < lb-tolarence) {
      bounds[j] = x[j]-lb;  // a negative value indicates that the lower bound is relaxed
    } else
    if (x[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 Generatefeasibilityreport(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;i<nrows;i++) {
   XPRSgetnamelist(xprob, 1, NULL, 0, &len, i, i);
   if (len>NameLength) NameLength=len;
  }
  for(i=0;i<ncols;i++) {
   XPRSgetnamelist(xprob, 2, NULL, 0, &len, i, i);
   if (len>NameLength) 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<nrows; i++) {
    if (constraintviolations[i] != 0) {

      // get constraint name
      XPRSgetnamelist(xprob, 1, name, NameLength, NULL, i, i);
      Getrowbounds(xprob, i, &lb, &ub);
      Bound_to_string(lb, slb); Bound_to_string(ub, sub);

      suminf += fabs(constraintviolations[i]);
      ninf++;

      if (constraintviolations[i] < 0) {
        Bound_to_string(lb+constraintviolations[i], slbshift);
        Bound_to_string(ub, subshift);
      } else {
        Bound_to_string(lb, slbshift);
        Bound_to_string(ub+constraintviolations[i], subshift);
      }

      printf("%.5i  %*s    %*s  %*s   %*s  %*s\n",i,NameLength,name, 7, slb, 7, slbshift, 7, sub, 7, subshift);
    }
  }

  printf("\n");
  printf("The following bound constraints were repaired.\n");
  printf("\n");

  printf("Index   Name%*s   Lower bound      Repaired     Upper bound      Repaired \n",NameLength-4,"");

  for (j=0; j<ncols; j++) {
    if (boundviolations[j] != 0) {

      // get constraint name
      XPRSgetnamelist(xprob, 2, name, NameLength, NULL, j, j);
      XPRSgetlb(xprob, &lb, j, j);
      XPRSgetub(xprob, &ub, j, j);
      Bound_to_string(lb, slb); Bound_to_string(ub, sub);

      suminf += fabs(boundviolations[j]);
      ninf++;

      if (constraintviolations[j] < 0) {
        Bound_to_string(boundviolations[j], slbshift);
        Bound_to_string(ub, subshift);
      } else {
        Bound_to_string(lb, slbshift);
        Bound_to_string(boundviolations[j], subshift);
      }

      printf("%.5i  %*s    %*s  %*s   %*s  %*s\n",j,NameLength,name, 7, slb, 7, slbshift, 7, sub, 7, subshift);
    }
  }

  printf("\n");
  printf("Sum of Violations %f\n",suminf);
  printf("Number of corrections: %i\n", ninf);
  printf("\n");

  free(name);

  return(0);

}

/*
  Relaxes the problem given the values of the feasibility breakers
*/
int  XPRS_CC Applyrepairinfeasresulttoproblem(XPRSprob xprob, const double *constraintviolations, const double *boundviolations) {

  int i,j;
  int nrows, ncols;
  char boundtype;
  double lb,ub;

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

  // apply changes to rows
  for (i=0; i<nrows; i++) {
    if (constraintviolations[i] != 0) {

      // get constraint name
      Getrowbounds(xprob, i, &lb, &ub);

      // apply relaxation
      if (constraintviolations[i] > 0) {
        ub += constraintviolations[i];
      } else {
        lb += constraintviolations[i];
      }

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

  // apply changes to columns
  for (j=0; j<ncols; j++) {
    if (boundviolations[j] != 0) {

      // get constraint name
      XPRSgetlb(xprob, &lb, j, j);
      XPRSgetub(xprob, &ub, j, j);

      if (boundviolations[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);
}





Back to examples browserPrevious exampleNext example