/********************************************************
  Xpress-BCL C Example Problems
  =============================

  file foliorep.c
  ```````````````
  Modeling a MIP problem
  to perform portfolio optimization.

  Same model as in foliomip3.c.
  -- Infeasible model parameter values --
  -- Repairing infeasibilities --

  (c) 2009-2024 Fair Isaac Corporation
      author: S.Heipcke, June 2009, rev. Oct. 2014
********************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include "xprb.h"
#include "xprs.h"


#define MAXNUM 4                   /* Max. number of different assets */
#define MAXRISK 1.0/3              /* Max. investment into high-risk values */
#define MINREG 0.3                 /* Min. investment per geogr. region */
#define MAXREG 0.5                 /* Max. investment per geogr. region */
#define MAXSEC 0.15                /* Max. investment per ind. sector */
#define MAXVAL 0.2                 /* Max. investment per share */
#define MINVAL 0.1                 /* Min. investment per share */

#define DATAFILE "folio10.cdat"    /* File with problem data */
#define MAXENTRIES 10000


int NSHARES;                       /* Number of shares */
int NRISK;                         /* Number of high-risk shares */
int NREGIONS;                      /* Number of geographical regions */
int NTYPES;                        /* Number of share types */

double *RET;                       /* Estimated return in investment */
int *RISK;                         /* High-risk values among shares */
char **LOC;                        /* Geogr. region of shares */
char **SEC;                        /* Industry sector of shares */

char **SHARES_n;
char **REGIONS_n;
char **TYPES_n;

XPRBprob prob;
XPRBvar *frac;                     /* Fraction of capital used per share */
XPRBvar *buy;                      /* 1 if asset is in portfolio, 0 otherwise */
XPRBctr *allCtr;                   /* array containing all defined contraints */
int allCtrCount = 0;               /* number of constraints in allCtr array */

#include "readfoliodata.c_"

void print_sol_opt(void);
void print_violated(void);

int main(int argc, char **argv)
{
  int s, r, t;
  XPRBctr Return, *Risk, *Cap, *Num;
  XPRBctr *MinReg, *MaxReg, *LimSec, *LinkL, *LinkU;

  readdata(DATAFILE);               /* Data input from file */

  prob = XPRBnewprob("FolioMIP3inf");  /* Initialize a new problem in BCL */

  /* Create the decision variables (including upper bounds for `frac') */
  frac = (XPRBvar*)malloc(NSHARES*sizeof(XPRBvar));
  buy = (XPRBvar*)malloc(NSHARES*sizeof(XPRBvar));
  for (s = 0; s<NSHARES; s++) {
    frac[s] = XPRBnewvar(prob, XPRB_PL, "frac", 0, MAXVAL);
    buy[s] = XPRBnewvar(prob, XPRB_BV, "buy", 0, 1);
  }

  /* Objective: total return */
  Return = XPRBnewctr(prob, "Return", XPRB_N);
  for (s = 0; s<NSHARES; s++) XPRBaddterm(Return, frac[s], RET[s]);
  XPRBsetobj(prob, Return);          /* Set the objective function */

  /* allocate memory for all constraints: Risk(1), Cap(1), Num(1), MinReg(NREGIONS), MaxReg(NREGIONS), LimSec(NTYPES), LinkU(NSHARES), LinkL(NSHARES) */
  allCtr = (XPRBctr*)malloc(sizeof(XPRBctr)* (3 + 2 * NREGIONS + NTYPES + 2 * NSHARES));
  /* init constraint pointers */
  Risk   = &allCtr[allCtrCount++];
  Cap    = &allCtr[allCtrCount++];
  Num    = &allCtr[allCtrCount++];
  MinReg = &allCtr[allCtrCount]; allCtrCount += NREGIONS;
  MaxReg = &allCtr[allCtrCount]; allCtrCount += NREGIONS;
  LimSec = &allCtr[allCtrCount]; allCtrCount += NTYPES;
  LinkU  = &allCtr[allCtrCount]; allCtrCount += NSHARES;
  LinkL  = &allCtr[allCtrCount]; allCtrCount += NSHARES;

  /* Limit the percentage of high-risk values */
  *Risk = XPRBnewctr(prob, "Risk", XPRB_L);
  for (s = 0; s<NRISK; s++) XPRBaddterm(*Risk, frac[RISK[s]], 1);
  XPRBaddterm(*Risk, NULL, MAXRISK);

  /* Limits on geographical distribution */
  for (r = 0; r<NREGIONS; r++) {
    MinReg[r] = XPRBnewctr(prob, XPRBnewname("MinReg(%s)", REGIONS_n[r]), XPRB_G);
    MaxReg[r] = XPRBnewctr(prob, XPRBnewname("MaxReg(%s)", REGIONS_n[r]), XPRB_L);
    for (s = 0; s<NSHARES; s++)
      if (LOC[r][s]>0) {
        XPRBaddterm(MinReg[r], frac[s], 1);
        XPRBaddterm(MaxReg[r], frac[s], 1);
      }
    XPRBaddterm(MinReg[r], NULL, MINREG);
    XPRBaddterm(MaxReg[r], NULL, MAXREG);
  }

  /* Diversification across industry sectors */
  for (t = 0; t<NTYPES; t++) {
    LimSec[t] = XPRBnewctr(prob, XPRBnewname("LimSec(%s)", TYPES_n[t]), XPRB_L);
    for (s = 0; s<NSHARES; s++) if (SEC[t][s]>0) XPRBaddterm(LimSec[t], frac[s], 1);
    XPRBaddterm(LimSec[t], NULL, MAXSEC);
  }

  /* Spend all the capital */
  *Cap = XPRBnewctr(prob, "Cap", XPRB_E);
  for (s = 0; s<NSHARES; s++) XPRBaddterm(*Cap, frac[s], 1);
  XPRBaddterm(*Cap, NULL, 1);

  /* Limit the total number of assets */
  *Num = XPRBnewctr(prob, "Num", XPRB_L);
  for (s = 0; s<NSHARES; s++) XPRBaddterm(*Num, buy[s], 1);
  XPRBaddterm(*Num, NULL, MAXNUM);

  /* Linking the variables */
  for (s = 0; s<NSHARES; s++) {
    LinkU[s] = XPRBnewctr(prob, XPRBnewname("LinkU(%s)", SHARES_n[s]), XPRB_L);
    XPRBaddterm(LinkU[s], frac[s], 1);
    XPRBaddterm(LinkU[s], buy[s], -MAXVAL);
    LinkL[s] = XPRBnewctr(prob, XPRBnewname("LinkL(%s)", SHARES_n[s]), XPRB_G);
    XPRBaddterm(LinkL[s], frac[s], 1);
    XPRBaddterm(LinkL[s], buy[s], -MINVAL);
  }

  /* Solve the problem (LP) */
  XPRBsetmsglevel(prob, 1);
  XPRBsetsense(prob, XPRB_MAXIM);
  XPRBlpoptimize(prob, "");

  if (XPRBgetlpstat(prob) == XPRB_LP_INFEAS) {
    char *rstat[] = { "relaxed optimum found", "relaxed problem infeasible",
      "relaxed problem unbounded", "solution nonoptimal for original objective",
      "error", "numerical instability" };
    double *lrp, *grp; /* Selectors for LEG / GEQ constraints */
    double *lbp, *ubp; /* Selector for lower / upper bounds on vars */
    double delta;
    int ncol, nrow, repstatus;
    XPRSprob op;

    printf("LP infeasible. Start infeasibility repair.\n");

    op = XPRBgetXPRSprob(prob);     /* Retrieve the Optimizer problem */

    /* Must use the weighted infeasibility repair method since
    only some constraints of each type may be relaxed */

    /*
    lrp: (affects = and <= rows)
    ax - aux_var  = b
    ax - aux_var <= b
    grp: (affects = and >= rows)
    ax + aux_var  = b
    ax + aux_var >= b
    lbp:
    x_i + aux_var >= l
    ubp:
    x_i - aux_var <= u
    */

    XPRSgetintattrib(op, XPRS_ORIGINALCOLS, &ncol);
    XPRSgetintattrib(op, XPRS_ORIGINALROWS, &nrow);
    lrp = (double *)calloc(nrow, sizeof(double));
    grp = (double *)calloc(nrow, sizeof(double));
    lbp = (double *)calloc(ncol, sizeof(double));
    ubp = (double *)calloc(ncol, sizeof(double));

    lrp[XPRBgetrownum(*Risk)] = 1;
    for (r = 0; r<NREGIONS; r++) lrp[XPRBgetrownum(MaxReg[r])] = 1;
    for (t = 0; t<NTYPES; t++) lrp[XPRBgetrownum(LimSec[t])] = 1;
    lrp[XPRBgetrownum(*Num)] = 1;
    for (r = 0; r<NREGIONS; r++) grp[XPRBgetrownum(MinReg[r])] = 1;

    for (delta = 0.001; delta < 10; delta *= 10) {
      XPRSrepairweightedinfeas(op, &repstatus, lrp, grp, lbp, ubp, 'r', delta, "");
      printf("delta = %g: Status: %s\n", delta, rstat[repstatus]);
      if (repstatus == 0) {
        XPRBsync(prob, XPRB_XPRS_SOLMIP);
        print_sol_opt();
        print_violated();
      }
    }

    free(lrp);
    free(grp);
    free(lbp);
    free(ubp);
  }

  return 0;
}

/**************************Solution printing****************************/

void print_sol_opt(void)
{
  int s;
  printf("  Total return: %g\n", XPRBgetobjval(prob));
  for (s = 0; s<NSHARES; s++)
    if (XPRBgetsol(buy[s]) > 0.5)
      printf("  %s : %g%% (%g)\n", SHARES_n[s], XPRBgetsol(frac[s])*100, XPRBgetsol(buy[s]));
}

void print_violated(void)
{
  int i;
  char *type=NULL;
  printf(" Violated (relaxed) constraints:\n");
  for (i = 0; i < allCtrCount; i++) {
    double viol, slack = XPRBgetslack(allCtr[i]);
    switch (XPRBgetctrtype(allCtr[i])) {
      case XPRB_E: viol = fabs(slack); type = " ="; break;
      case XPRB_G: viol = slack;       type = ">="; break;
      case XPRB_L: viol = -slack;      type = "<="; break;
      default: printf("  unexpected constraint type\n"); continue;
    }
    if (viol>1e-6) printf("  %s constraint %s by %g\n", type, XPRBgetctrname(allCtr[i]), -slack);
  }
}