/***********************************************************************
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);
}
|
