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

   file ComplexUserFunctions.c
   ``````````````````````````
    Implement a problem using complex user functions

    This problem demonstrates how to

    1. add user functions implemented locally to an SLP problem,
    2. track how function evaluations happen,
    3. define user functions that return derivatives,
    4. solve a problem iteration by iteration,
    5. optionally write out the linearization of each iteration for analysis.

 Solve

 Minimize t                              (opt = 27)
   s.t.
   t = y ^ 2 + z - v + w ^ 2
   x ^ 2 + t ^ 2 <= 730                  (= 27*27+1)
   and
   1 <= x <= 2
   2 <= y <= 3
   3 <= z <= 4
   4 <= v <= 5
   5 <= w <= 6

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

#include <stdlib.h>
#include <stdio.h>
#include <stdarg.h>
#include <signal.h>
#include <string.h>
#include <math.h>
#include <float.h>

#include <xprs.h>
#include <xslp.h>

#define WRITE_LINEARIZATIONS 0

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:  /* dialogue */
    case 1:  /* information */
             printf("%s\n", msg);
             break;
    default: /* exiting - buffers need flushing */
             fflush(stdout);
             break;
  }
}

// user function implementing
// y^2 + z - v + w^2
// with derivatives
double XSLP_CC MyObj(double *InputValues, int *FunctionInfo, double *Deltas, double *ReturnArray)
{
  int i;
  // Information contained by FunctionInfo
  int CallFlag, nInput, nOutput, nDelta, nInStr, nOutStr, SLPUF, nInst;
  // User function name
  char UFName[256];
  // the XSLP problem object
  XSLPprob sprob;
  // the input values
  double y = 0;
  double z = 0;
  double v = 0;
  double w = 0;
  // value of the function
  double Value = 0;
  // the return array counter
  int nRes = 0;

  XSLPgetfuncinfo(FunctionInfo, &CallFlag, &nInput, &nOutput, &nDelta, &nInStr, &nOutStr, &SLPUF, &nInst);

  // retrieve the XSLP object pointer
  XSLPgetfuncobject(FunctionInfo, XSLP_XSLPPROBLEM, (void **) &sprob);

  // set up local variables
  if (nInput) {
    y = InputValues[0];
    z = InputValues[1];
    v = InputValues[2];
    w = InputValues[3];
  }

  // calculate the function value
  if (nInput && nOutput) {
    Value = y*y + z - v + w*w;

    // for functions returning multiple values, the first value in the return array is the calculated value at the given point
    ReturnArray[nRes++] = Value;

    // retrieve the user function's name
    XSLPitemname(sprob, XSLP_FUN, SLPUF, UFName);
    // report calculation
    printf("%s%f,%f,%f,%f) = %f\n",UFName,y,z,v,w,Value);
  }

  // return deltas (y^2 + z - v + w^2)
  // derivatives are expected to follow only for those for which the Delta array holds a nonzero value
  if (nInput && nDelta) {
    double dtmp = 1;
    if (Deltas[0]) ReturnArray[nRes++] = 2*y;
    if (Deltas[1]) ReturnArray[nRes++] = 1;
    if (Deltas[2]) ReturnArray[nRes++] = -1;
    if (Deltas[3]) ReturnArray[nRes++] = 2*w;
    // report derivatives
    printf("Diff(%s%f,%f,%f,%f)) = %f,%f,%f,%f\n",UFName,y,z,v,w,2*y,dtmp,-dtmp,2*w);
  }

  // return value indicates if the function call was successful
  return 0;
}

// user function implementing
// x ^ 2 + t ^ 2
// derivatives will be calculated using tangential approximation, note the perturbed calls in the output
double XSLP_CC Ball(double *InputValues, int *FunctionInfo)
{
  // Information contained by FunctionInfo
  int CallFlag, nInput, nOutput, nDelta, nInStr, nOutStr, SLPUF, nInst;
  // User function name
  char UFName[256];
  // the XSLP problem object
  XSLPprob sprob;
  // input values
  double x = 0;
  double t = 0;
  // value of the function
  double Value = 0;
  // the return array counter
  int nRes = 0;

  XSLPgetfuncinfo(FunctionInfo, &CallFlag, &nInput, &nOutput, &nDelta, &nInStr, &nOutStr, &SLPUF, &nInst);

  // set up local variables
  if (nInput) {
    x = InputValues[0];
    t = InputValues[1];
  }

  if (nInput) {

    Value = x*x + t*t;

    // retrieve the XSLP object pointer
    XSLPgetfuncobject(FunctionInfo, XSLP_XSLPPROBLEM, (void **) &sprob);

    // retrieve the user function's name
    XSLPitemname(sprob, XSLP_FUN, SLPUF, UFName);
    // report calculation
    printf("%s(%f,%f) = %f\n",UFName,x,t,Value);
  }

  // for functions returning a single value, the return value is the calculated value of the function
  return Value;
}

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

  XPRSprob xprob = NULL;
  XSLPprob sprob = NULL;
  int ReturnValue = 0;

  char buffer[1000];
  int loop;

  char *sProblem=NULL;        /* Problem name */

  /* Initialise Optimizer */
  ReturnValue = XPRSinit(NULL);
  XPRSgetbanner(buffer); printf("%s\n",buffer);
  if (ReturnValue == 8) {
    printf("Unable to initialize XPRESS\n");
    goto ErrorReturn;
  }

  XSLPinit();
  XSLPgetbanner(buffer); printf("%s\n",buffer);

  /* Initialise problem object */
  ReturnValue = XPRScreateprob(&xprob);
  if (ReturnValue != 0 && ReturnValue != 32) {
    printf("Unable to create XPRESS problem\n");
    XPRSfree();
    goto ErrorReturn;
  }

  // set messaging callback
  if (ReturnValue=XPRSsetcbmessage(xprob,Message,NULL)) goto ErrorReturn;

  // create SLP problem
  if (ReturnValue=XSLPcreateprob(&sprob, &xprob)) goto ErrorReturn;

  // use XSLP
  if (ReturnValue=XSLPsetintcontrol(sprob, XSLP_SOLVER, 0)) goto ErrorReturn;

/*  Minimize t
    s.t.
    t = y ^ 2 + z - v + w ^ 2
    x ^ 2 + t ^ 2 <= 730                  (27*27+1)
    and
    1 <= x <= 2
    2 <= y <= 3
    3 <= z <= 4
    4 <= v <= 5
    5 <= w <= 6

    Optimum should be 27 with y=2, z=3, v=5, w=5 (forced by objective) and x being 1 (forced by the second constrained) */

  {
   /* Linear data */
   int nRow          = 2;                  /* Number of rows */
   char sRowType[]   = {'E','L'};          /* Row types */
   double dRHS[]     = {0,730};            /* RHS values */
   double *dRange    = NULL;               /* Range values - none in this example */
   char sRowName[]   = "obj\0ball";        /* Row names */

   int nCol          = 7;                  /* Number of columns */
   double dObj[]     = {1,0,0,0,0,0,0};      /* Objective function coefficients */
   double dLowerBd[] = {0,1,2,3,4,5,1};      /* Lower bounds on the columns */
   double dUpperBd[] = {XPRS_PLUSINFINITY,2,3,4,5,6,1};    /* Upper bounds - note macro for infinity */
   char sColName[]   = "t\0x\0y\0z\0v\0w\0="; /* Column names */

   /* Note the existence of the "=" column fixed to 1 and used to add nonlinear terms that are not coefficients to the problem */
   /* Future releases of XSLP (XSLP with optimizer version 21.00 on onwards) will not require this to be done manually */

   /* Matrix data */
   int nColStart[]   = {0,1,1,1,1,1,1,1,1};  /* Start positions of the columns in dMatElem
                                              - note there are nCol+1 entries, the last
                                              indicating where the next column would
                                              start if there was one */
   int *nColElem     = NULL;               /* Number of elements in each column - not
                                              needed thanks to the last (optional) entry
                                              in nColStart */
   int nRowInd[]     = {0};   /* Row indices for the matrix elements */
   double dMatElem[] = {-1};  /* Matrix elements - arranged by column */

   /* Load the matrix into Optimizer */
   if (ReturnValue=XPRSloadlp(xprob,"UserFunctionExample_tracking",nCol,nRow,sRowType,dRHS,dRange,dObj,nColStart,nColElem,nRowInd,dMatElem,dLowerBd,dUpperBd)) goto ErrorReturn;
   /* add names */
   if (ReturnValue=XPRSaddnames(xprob, 1, sRowName, 0, 0)) goto ErrorReturn;
   if (ReturnValue=XPRSaddnames(xprob, 2, sColName, 0, 6)) goto ErrorReturn;
  }

  // Add the user functions to the SLP problem
  {
    int i,nToken;
    int Type[10];
    double Value[10];
    void *Func;

    // Add MyObj as function 1
    nToken = 0;
    if (ReturnValue=XSLPsetstring(sprob,&i,"MyObj")) goto ErrorReturn;
    Type[nToken] = XSLP_STRING;
    Value[nToken++] = (double) i;
    Type[nToken] = XSLP_UFEXETYPE;
    /*
      Octal number describing linkage, evaluation and derivatives.
      Each octal digit corresponds to the following
      digit 1: linkage type (1: DLL, 2:SLX, 3:XLF, 5:Mosel, 6:VB, 7:COM)
      digit 2: evaluation type (0: default, 1: every iteration, 2: if significant change in input)
      digit 3: derivative type (0: default, 1: tangential, 2: forward)
      +
      bit 9: calling style (0: standard, 1: CDECL)
    */
    Value[nToken++] = XSLP_DLL_EXETYPE; /* = (double) 001; linkage is DLL. We will redirect the function later, see below */
    Type[nToken] = XSLP_UFARGTYPE;
    /*
     Octal number describing the arguments of the user function.
     Each octal digit correponds to an argument, indicating its presence and its type
     The possible arguments in order: InputValues, FunctionInfo, InputNames, ReturnNames, Deltas, ReturnArray
     With the corresponding octal digit being:
     0: argument not present, 1: NULL, 2: int *, 3: double *, 4: variant, 6: char *
    */
    Value[nToken++] = (double) 0330023; /* so arguments are: InputValues as (double *), FunctionInfo as (int *), Deltas as (double *) and ReturnArray as (double *) */
    if (ReturnValue=XSLPsetstring(sprob,&i,"DirectC")) goto ErrorReturn; /* not important for internal function */
    Type[nToken] = XSLP_STRING;
    Value[nToken++] = (double) i;
    Type[nToken] = XSLP_EOF;

    if (ReturnValue=XSLPloaduserfuncs(sprob,1,Type,Value)) goto ErrorReturn;
    if (ReturnValue=XSLPaddnames(sprob,XSLP_USERFUNCNAMES,"MyObj",1,1)) goto ErrorReturn;
    Func = MyObj;
    if (ReturnValue=XSLPchguserfuncaddress(sprob,1,&Func)) goto ErrorReturn; /* redirect to internal function */

    // Add Ball as function 2
    nToken = 0;
    if (ReturnValue=XSLPsetstring(sprob,&i,"Ball")) goto ErrorReturn;

    Type[nToken] = XSLP_STRING;
    Value[nToken++] = (double) i;
    Type[nToken] = XSLP_UFEXETYPE;
    Value[nToken++] = XSLP_DLL_EXETYPE + XSLP_2DERIVATIVE; // = (double) 0101; linkage is dll + tangential derivaties
    Type[nToken] = XSLP_UFARGTYPE;
    Value[nToken++] = (double) 0000023; // only InputValues as (double *) and FunctionInfo as (int *) are present as arguments
    XSLPsetstring(sprob,&i,"DirectC"); /* not important for internal function */
    Type[nToken] = XSLP_STRING;
    Value[nToken++] = (double) i;
    Type[nToken] = XSLP_EOF;

    if (ReturnValue=XSLPadduserfuncs(sprob,1,Type,Value)) goto ErrorReturn;
    if (ReturnValue=XSLPaddnames(sprob,XSLP_USERFUNCNAMES,"Ball",2,2)) goto ErrorReturn;
    Func = Ball;
    if (ReturnValue=XSLPchguserfuncaddress(sprob,2,&Func)) goto ErrorReturn; /* redirect to internal function */
  }

  // add the nonlinear functions to the slp problem
  {
    int Col;
    int nToken = 0;
    int nCoef = 0;
    int RowIndex[10];
    int ColIndex[10];
    double Factor[10];
    int Type[100];
    double Value[100];
    int Equal;
    int FormulaStart[10];

    if (ReturnValue=XSLPgetintattrib(sprob, XSLP_EQUALSCOLUMN, &Equal)) goto ErrorReturn;
    // future releases of SLP will not require this...
    if (Equal == 0) {
      if (ReturnValue=XPRSgetindex(xprob, 2, "=", &Equal)) goto ErrorReturn;
    }

    // MyObj(y,z,v,w) = y^2 + z - v + w^2 into row 0
    RowIndex[nCoef] = 0;
    ColIndex[nCoef] = Equal;
    Factor[nCoef] = 1.0;
    FormulaStart[nCoef++] = nToken;

    Type[nToken] = XSLP_FUN; Value[nToken++] = 1;
    Type[nToken] = XSLP_LB;  Value[nToken++] = 0;
    XPRSgetindex(xprob, 2, "y", &Col); Type[nToken] = XSLP_COL; Value[nToken++] = Col;
    Type[nToken] = XSLP_DEL;  Value[nToken++] = XSLP_COMMA;
    XPRSgetindex(xprob, 2, "z", &Col); Type[nToken] = XSLP_COL; Value[nToken++] = Col;
    Type[nToken] = XSLP_DEL;  Value[nToken++] = XSLP_COMMA;
    XPRSgetindex(xprob, 2, "v", &Col); Type[nToken] = XSLP_COL; Value[nToken++] = Col;
    Type[nToken] = XSLP_DEL;  Value[nToken++] = XSLP_COMMA;
    if (ReturnValue=XPRSgetindex(xprob, 2, "w", &Col)) goto ErrorReturn;
    Type[nToken] = XSLP_COL; Value[nToken++] = Col;
    Type[nToken] = XSLP_RB;  Value[nToken++] = 0;
    Type[nToken] = XSLP_EOF; Value[nToken++] = 0;

    // BALL(x,t) = ( x ^ 2 + t ^ 2 ) into row 1
    RowIndex[nCoef] = 1;
    ColIndex[nCoef] = Equal;
    Factor[nCoef] = 1.0;
    FormulaStart[nCoef++] = nToken;

    Type[nToken] = XSLP_FUN; Value[nToken++] = 2;
    Type[nToken] = XSLP_LB;  Value[nToken++] = 0;
    if (ReturnValue=XPRSgetindex(xprob, 2, "x", &Col)) goto ErrorReturn;
    Type[nToken] = XSLP_COL; Value[nToken++] = Col;
    Type[nToken] = XSLP_DEL;  Value[nToken++] = XSLP_COMMA;
    if (ReturnValue=XPRSgetindex(xprob, 2, "t", &Col)) goto ErrorReturn;
    Type[nToken] = XSLP_COL; Value[nToken++] = Col;
    Type[nToken] = XSLP_RB;  Value[nToken++] = 0;
    Type[nToken] = XSLP_EOF; Value[nToken++] = 0;
    if (ReturnValue=XSLPaddcoefs(sprob,nCoef,RowIndex,ColIndex,Factor,FormulaStart,0,Type,Value)) goto ErrorReturn;
  }

  if (WRITE_LINEARIZATIONS) if (ReturnValue=XSLPwriteprob(sprob,"ProblemLoaded","l")) goto ErrorReturn;

  // Solve problem
  {
    int Status;
    int nIterations = 1;

    if (ReturnValue=XSLPsetintcontrol(sprob, XSLP_DELTAZLIMIT , 4)) goto ErrorReturn;

    if (ReturnValue=XSLPsetintcontrol(sprob, XSLP_ITERLIMIT, nIterations)) goto ErrorReturn;
    if (ReturnValue=XSLPminim(sprob,"")) goto ErrorReturn;
    if (ReturnValue=XSLPgetintattrib(sprob, XSLP_STATUS, &Status)) goto ErrorReturn;
    if (ReturnValue=XPRSsetdblcontrol(xprob, XPRS_OUTPUTTOL, 1e-17)) goto ErrorReturn;
    if (WRITE_LINEARIZATIONS) if (ReturnValue=XSLPwriteprob(sprob, "Linearization1.lp","la")) goto ErrorReturn;
    // if stopped because of the iteration limit
    while (Status & XSLP_STATUS_MAXSLPITERATIONS) {
      char str[256];
      nIterations++;
      sprintf(str,"Linearization%i.lp",nIterations);
      if (ReturnValue=XSLPsetintcontrol(sprob, XSLP_ITERLIMIT, nIterations)) goto ErrorReturn;
      if (ReturnValue=XSLPreminim(sprob,"")) goto ErrorReturn;
      if (WRITE_LINEARIZATIONS) if (ReturnValue=XSLPwriteprob(sprob, str,"ap")) goto ErrorReturn;
      if (ReturnValue=XSLPgetintattrib(sprob, XSLP_STATUS, &Status)) goto ErrorReturn;
    }
  }

  // retrieve and report the optimal solution
  {
    int nCol, Col;
    double *dsol;
    if (ReturnValue=XSLPgetintattrib(sprob,XPRS_COLS,&nCol)) goto ErrorReturn;
    dsol = (double *) malloc(nCol*sizeof(double));
    if (!dsol) {
      printf("Out of memory.\n");
      goto ErrorReturn;
    }
    if (ReturnValue=XSLPgetslpsol(sprob,dsol,NULL,NULL,NULL)) goto ErrorReturn;
    printf("Solution:\n");
    if (ReturnValue=XPRSgetindex(xprob, 2, "x", &Col)) goto ErrorReturn; printf("  x = %f\n",dsol[Col]);
    if (ReturnValue=XPRSgetindex(xprob, 2, "y", &Col)) goto ErrorReturn; printf("  y = %f\n",dsol[Col]);
    if (ReturnValue=XPRSgetindex(xprob, 2, "z", &Col)) goto ErrorReturn; printf("  z = %f\n",dsol[Col]);
    if (ReturnValue=XPRSgetindex(xprob, 2, "v", &Col)) goto ErrorReturn; printf("  v = %f\n",dsol[Col]);
    if (ReturnValue=XPRSgetindex(xprob, 2, "w", &Col)) goto ErrorReturn; printf("  w = %f\n",dsol[Col]);
  }


goto NormalReturn;
ErrorReturn:
  printf("\nError %d",ReturnValue);
NormalReturn:
  /* Destroy problem object and free environment */
  XSLPdestroyprob(sprob);
  XPRSdestroyprob(xprob);
  XSLPfree();
  XPRSfree();

  return ReturnValue;
}