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

   file tableau.c
   ```````````````
   Read in a user-specified LP and solve it, displaying both the
   initial and optimal tableau.

   Inputs an MPS matrix file and required optimisation sense, and
   proceeds to solve the problem with lpoptimize. The simplex
   algorithm is interrupted to get its intial basis, and a tableau is
   requested with a call to function showtab. Once the solution is
   found, a second call produces the optimal tableau.
   Function showtab retrieves the pivot order of the basic variables,
   along with other problem information, and then constructs (and
   displays) the tableau row-by-row using the backwards transformation,
   btran.

   Arguments:   int argc          number of arguments (2 or 3 expected)
                char argv[0]      program name
                char argv[1]      problem name
                char argv[2]      objective sense ("min" or "max")
                                    if this argument is not provided the
                                    problem will be minimised, by default

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "xprs.h"

/* Calls an Xpress optimizer function and checks the return code.
 * If the call fails then the function
 * - prints a short error message to stderr,
 * - sets variable 'returnCode' to the error,
 * - and branches to label 'cleanup'.
 */
#define CHECK_RETURN(call) do {                         \
    int result_ = call;                                 \
    if ( result_ != 0 ) {                               \
      fprintf(stderr, "Line %d: %s failed with %d\n",   \
              __LINE__, #call, result_);                \
      returnCode = result_;                             \
      goto cleanup;                                     \
    }                                                   \
  } while (0)

#define MIN 1                        /* Sense values */
#define MAX 0

static int showtab(XPRSprob prob, const char *sState, const char *sProblem);
static void XPRS_CC messagecb(XPRSprob cbprob, void* cbdata,
                              const char *msg, int len, int msgtype);

int main(int argc, char *argv[])
{
  int returnCode = 0;
  XPRSprob prob = NULL;             /* The problem instance */
  char   *sProblem;                 /* Problem name */
  char   sLogFile[]="tableau.log";  /* Optimizer log file */
  int    nSense = MIN;              /* Optimisation sense: 1 if min, otherwise
                                       max */
  int    nIteration;                /* Number of simplex iterations */
  double dObjValue;                 /* Objective function value */

  /* Check the number of arguments is correct */
  if (argc < 2 || argc > 3) {
    printf("Usage: tableau <matrix> [min|max]");
    return 1;
  }

  /* Set the problem name */
  sProblem = argv[1];

  /* Set the optimisation sense (default is minim) */
  if (argc == 3) {
    if (strcmp(argv[2],"max") == 0)
      nSense=MAX;
    else if (strcmp(argv[2],"min") != 0) {
      printf("Usage: 'min' or 'max' expected as second argument\n\n");
      return 1;
    }
  }

  /* Initialize the optimizer. */
  if ( XPRSinit("") != 0 ) {
    char message[512];
    XPRSgetlicerrmsg(message, sizeof(message));
    fprintf(stderr, "Licensing error: %s\n", message);
    return 1;
  }

  /* Create a new problem and immediately register a message handler.
   * Once we have a message handler installed, errors will produce verbose
   * error messages on the console and we can limit ourselves to minimal
   * error handling in the code here.
   */
  CHECK_RETURN( XPRScreateprob(&prob) );
  CHECK_RETURN( XPRSaddcbmessage(prob, messagecb, NULL, 0) );

  /* Delete and define the Optimizer log file */
  remove(sLogFile);
  CHECK_RETURN( XPRSsetlogfile(prob, sLogFile) );

  /* Input the matrix */
  CHECK_RETURN( XPRSreadprob(prob, sProblem,"") );

  /* Turn presolve off as we want to display the initial tableau */
  CHECK_RETURN( XPRSsetintcontrol(prob, XPRS_PRESOLVE, 0) );

  /* Set the number of simplex iterations to zero */
  CHECK_RETURN( XPRSsetintcontrol(prob, XPRS_LPITERLIMIT, 0) );

  /* Perform the first step in the simplex algorithm */
  if (nSense == MAX) {
    CHECK_RETURN( XPRSchgobjsense(prob, XPRS_OBJ_MAXIMIZE) );
  }
  CHECK_RETURN( XPRSlpoptimize(prob,"") );

  /* Display the initial tableau */
  CHECK_RETURN( showtab(prob, "Initial", sProblem) );

  /* Continue with the simplex algorithm */
  CHECK_RETURN( XPRSsetintcontrol(prob, XPRS_LPITERLIMIT, 1000000) );

  CHECK_RETURN( XPRSlpoptimize(prob,"") );

  /* Get and display the objective function value and the iteration count */
  CHECK_RETURN( XPRSgetdblattrib(prob, XPRS_LPOBJVAL, &dObjValue) );

  CHECK_RETURN( XPRSgetintattrib(prob, XPRS_SIMPLEXITER, &nIteration) );

  printf("M%simal solution with value %g found in %d iterations.\n",
         (nSense) ? "in" : "ax", dObjValue, nIteration);

  /* Display the optimal tableau */
  CHECK_RETURN( showtab(prob, "Optimal", sProblem) );

 cleanup:
  if (returnCode > 0) {
    /* There was an error with the solver. Get the error code and error message.
     * If prob is still NULL then the error was in XPRScreateprob() and
     * we cannot find more detailed error information.
     */
    if (prob != NULL) {
      int errorCode = -1;
      char errorMessage[512] = {0};
      XPRSgetintattrib(prob, XPRS_ERRORCODE, &errorCode);
      XPRSgetlasterror(prob, errorMessage);
      fprintf(stderr, "Error %d: %s\n", errorCode, errorMessage);
    }
  }

  /* Free the resources (variables are initialized so that this is valid
   * even in case of error).
   */
  XPRSdestroyprob(prob);
  XPRSfree();

  return returnCode;
}


/*****************************************************************************\
 * Name:         showtab
 *
 * Purpose:      Display tableau on screen
 *
 * Arguments:    XPRSprob prob            The problem
 *
 *               const char *sState       Problem state
 *
 *               const char *sProblem     Problem name
 *
 * Return Value: None.
 *
\*****************************************************************************/

int showtab(XPRSprob prob, const char *sState, const char *sProblem)
{
  int returnCode = 0;
  int    nRow;                  /* Number of rows */
  int    nCol;                  /* Number of columns */
  int    nSprR;
  int    nVector;               /* Number of vectors: rows + cols + 1 (for RHS)
                                 */
  int    namlen;                /* Maximum length of row or column name */
  char   *sVectorName = NULL;   /* Row or column name */
  int    nElem;                 /* Number of non-zero matrix elements */
  int    *pRowInd = NULL;       /* Row indices of the matrix elements */
  int    iRow;                  /* Row index holder */
  double *pMatElem = NULL;      /* Values of the matrix elements */
  int    pColStart[2];          /* Start positions of the columns in pMatElem */
  int    *pPivot = NULL;        /* Pivot order of the basic variables */
  int    i, j;                  /* Loop counters */
  double *y = NULL,*z = NULL, d;/* Variables to form tableau */
  double *x = NULL;             /* Primal solution values */
  double *slack = NULL;         /* Slack values */


  /*** Get problem information ***/

  /* Get the number of rows */
  CHECK_RETURN( XPRSgetintattrib(prob, XPRS_ROWS, &nRow) );

  /* Get the number of columns */
  CHECK_RETURN( XPRSgetintattrib(prob, XPRS_COLS, &nCol) );

  /* Get the number of columns */
  CHECK_RETURN( XPRSgetintattrib(prob, XPRS_SPAREROWS, &nSprR) );

  nVector = nRow + nCol + 1;

  /* Get the maximum name length */
  CHECK_RETURN( XPRSgetintattrib(prob, XPRS_NAMELENGTH, &namlen) );
  namlen *= 8;  /* Because name length is returned in 8-byte units */

  /* Allocate memory to the arrays */
  y           = malloc(nRow    * sizeof(*y));
  z           = malloc(nVector * sizeof(*z));
  pRowInd     = malloc(nRow    * sizeof(*pRowInd));
  pMatElem    = malloc(nRow    * sizeof(*pMatElem));
  pPivot      = malloc(nRow    * sizeof(*pPivot));
  x           = malloc(nCol    * sizeof(*x));
  slack       = malloc(nRow    * sizeof(*slack));
  sVectorName = malloc((namlen+1)*sizeof(*sVectorName));

  /* Check that the memory has been successfully allocated */
  if (!y || !z || !pRowInd || !pMatElem || !sVectorName || !pPivot) {
    perror("malloc");
    returnCode = -2;
    goto cleanup;
  }

  /* Retrieve the pivot order of the basic variables */
  CHECK_RETURN( XPRSgetpivotorder(prob, pPivot) );

  /*** Construct and display the tableau ***/

  /* Printer header of matrix names */
  printf("\n%s tableau of Problem %s.\n", sState, sProblem);
  printf("Note: slacks on G type rows range from -infinity to 0\n\n     ");
  for (i=0; i<nRow; i++) {
    /* Get and print the individual row names */
    CHECK_RETURN( XPRSgetnamelist(prob, 1, sVectorName, namlen+1, NULL, i, i) );
    printf(" %-3.3s   ", sVectorName);
  }
  for (i=0; i<nCol; i++) {
    /* Get and print the individual column names */
    CHECK_RETURN( XPRSgetnamelist(prob, 2, sVectorName, namlen+1, NULL, i, i) );
    printf(" %-3.3s   ", sVectorName);
  }
  printf(" RHS\n");

  /* Get the tableau row-by-row using the backwards transformation, btran */

  /* For each row iRow in turn, calculate z = e_irow * B^-1 * A */
  for (iRow=0; iRow<nRow; iRow++) {

    /* Set up the unit vector e_irow */
    for (i=0; i<nRow; i++) y[i]=0.0;
    y[iRow]=1.0;

    /* y = e_irow * B^-1 */
    CHECK_RETURN( XPRSbtran(prob, y) );

    /* Form z = y * A for each slack column of A */
    for (i=0; i<nRow; i++) z[i]=y[i];

    /* Form z = y * A, for each structural column of A */
    for (j=0; j<nCol; j++) {
      CHECK_RETURN( XPRSgetcols(prob, pColStart, pRowInd, pMatElem, nRow,
                                &nElem, j, j) );
      for (d=0.0, i=0; i<nElem; i++)
        d += y[pRowInd[i]] * pMatElem[i];
      z[nRow + j]=d;
    }

    /* Form z for RHS */
    CHECK_RETURN( XPRSgetlpsol(prob, x, slack, NULL, NULL) );

    if (pPivot[iRow]>=nRow)
      z[nVector-1]=x[pPivot[iRow]-nRow-nSprR];
    else
      z[nVector-1]=slack[pPivot[iRow]];

    /* Display single tableau row */
    if (pPivot[iRow]<nRow) {
      /* Pivot is a row */
      CHECK_RETURN( XPRSgetnamelist(prob, 1, sVectorName, namlen+1, NULL,
                                    pPivot[iRow], pPivot[iRow]) );
    }
    else {
      /* Pivot is a column */
      CHECK_RETURN( XPRSgetnamelist(prob, 2, sVectorName, namlen+1, NULL,
                                    pPivot[iRow]-nRow-nSprR,
                                    pPivot[iRow]-nRow-nSprR) );
    }
    printf("%-3.3s:", sVectorName);
    for (j=0; j<nVector; j++) {
      if (z[j]>=0.1 || z[j]<=-0.1)
        printf("%5.1f  ", z[j]);
      else
        printf("       ");
    }
    printf("\n");
  }
  printf("\n");

 cleanup:
  free(y);
  free(z);
  free(pRowInd);
  free(pMatElem);
  free(pPivot);
  free(sVectorName);
  return returnCode;
}

/* XPRS optimizer message callback */
void XPRS_CC messagecb(XPRSprob cbprob, void* cbdata,
                       const char *msg, int len, int msgtype)
{
  (void)cbprob;   /* unused (the problem for which the message is issued) */
  (void)cbdata;   /* unused (the data passed to XPRSaddcbmessage()) */
  switch(msgtype)
  {
  case 4:  /* error */
  case 3:  /* warning */
  case 2:  /* not used */
  case 1:  /* information */
    printf("%*s\n", len, msg);
    break;
  default: /* exiting - buffers need flushing */
    fflush(stdout);
    break;
  }
}