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

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

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

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

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <cmath>
#include <vector>
#include "xprb_cpp.h"
#include "xprs.h"

using namespace std;
using namespace ::dashoptimization;

#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 **SECT;                        // Industry sector of shares

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

XPRBvar *frac;                     // Fraction of capital used per share
XPRBvar *buy;                      // 1 if asset is in portfolio, 0 otherwise

#include "readfoliodata.c_"

void print_sol_opt(XPRBprob &p);
void print_violated(vector<XPRBctr> &ctrs);

int main(int argc, char **argv)
{
  int s, r, t;
  XPRBprob p("FolioMIP3inf");       // Initialize a new problem in BCL
  XPRBctr Return, *Risk, *Num, *MinReg, *MaxReg, *LimSec;
  XPRBexpr le, le2, Cap, LinkL, LinkU;

  readdata(DATAFILE);               // Data input from file

  // Create the decision variables (including upper bounds for `frac')
  frac = new XPRBvar[NSHARES];
  buy = new XPRBvar[NSHARES];
  for (s = 0; s<NSHARES; s++) {
    frac[s] = p.newVar("frac", XPRB_PL, 0, MAXVAL);
    buy[s] = p.newVar("buy", XPRB_BV);
  }

  // Objective: total return
  for (s = 0; s<NSHARES; s++) le += RET[s] * frac[s];
  Return = p.newCtr(le);
  p.setObj(Return);                 // Set the objective function

  // allocate memory for all constraints: Risk(1), Num(1), MinReg(NREGIONS), MaxReg(NREGIONS), LimSec(NTYPES)
  vector<XPRBctr> allCtr(2 + 2 * NREGIONS + NTYPES);
  int allCtrCount = 0;
  /* init constraint pointers */
  Risk = &allCtr[allCtrCount++];
  Num = &allCtr[allCtrCount++];
  MinReg = &allCtr[allCtrCount]; allCtrCount += NREGIONS;
  MaxReg = &allCtr[allCtrCount]; allCtrCount += NREGIONS;
  LimSec = &allCtr[allCtrCount]; allCtrCount += NTYPES;


  // Limit the percentage of high-risk values
  le = 0;
  for (s = 0; s<NRISK; s++) le += frac[RISK[s]];
  *Risk = p.newCtr("Risk", le <= MAXRISK);

  // Limits on geographical distribution
  for (r = 0; r<NREGIONS; r++) {
    le = 0; le2 = 0;
    for (s = 0; s<NSHARES; s++)
    if (LOC[r][s]>0) {
      le += frac[s];
      le2 += frac[s];
    }
    MinReg[r] = p.newCtr(XPRBnewname("MinReg(%s)", REGIONS_n[r]), le >= MINREG);
    MaxReg[r] = p.newCtr(XPRBnewname("MaxReg(%s)", REGIONS_n[r]), le2 <= MAXREG);
  }

  // Diversification across industry sectors
  for (t = 0; t<NTYPES; t++) {
    le = 0;
    for (s = 0; s<NSHARES; s++)
    if (SECT[t][s]>0) le += frac[s];
    LimSec[t] = p.newCtr(XPRBnewname("LimSec(%s)", TYPES_n[t]), le <= MAXSEC);
  }

  // Spend all the capital
  for (s = 0; s<NSHARES; s++) Cap += frac[s];
  p.newCtr("Cap", Cap == 1);

  // Limit the total number of assets
  le = 0;
  for (s = 0; s<NSHARES; s++) le += buy[s];
  *Num = p.newCtr("Num", le <= MAXNUM);

  // Linking the variables
  for (s = 0; s<NSHARES; s++) {
    p.newCtr(frac[s] <= MAXVAL*buy[s]);
    p.newCtr(frac[s] >= MINVAL*buy[s]);
  }

  // Solve the problem (LP)
  p.setMsgLevel(1);
  p.setSense(XPRB_MAXIM);
  p.lpOptimize("");

  if (p.getLPStat() == XPRB_LP_INFEAS) {
    cout << "LP infeasible. Start infeasibility repair." << endl;

    XPRSprob op = p.getXPRSprob();     // 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
    */

    int ncol, nrow, repstatus;
    double *lrp, *grp, *lbp, *ubp;
    XPRSgetintattrib(op, XPRS_ORIGINALCOLS, &ncol);
    XPRSgetintattrib(op, XPRS_ORIGINALROWS, &nrow);
    lrp = new double[nrow];
    grp = new double[nrow];
    lbp = new double[ncol];
    ubp = new double[ncol];
    memset(lrp, 0, nrow*sizeof(double));
    memset(grp, 0, nrow*sizeof(double));
    memset(lbp, 0, ncol*sizeof(double));
    memset(ubp, 0, ncol*sizeof(double));

    lrp[Risk->getRowNum()] = 1;
    for (r = 0; r<NREGIONS; r++) lrp[MaxReg[r].getRowNum()] = 1;
    for (t = 0; t<NTYPES; t++) lrp[LimSec[t].getRowNum()] = 1;
    lrp[Num->getRowNum()] = 1;
    for (r = 0; r<NREGIONS; r++) grp[MinReg[r].getRowNum()] = 1;

    char *rstat[] = { "relaxed optimum found", "relaxed problem infeasible",
      "relaxed problem unbounded", "solution nonoptimal for original objective",
      "error", "numerical instability" };

    for (double delta = 0.001; delta < 10; delta *= 10) {
      XPRSrepairweightedinfeas(op, &repstatus, lrp, grp, lbp, ubp, 'd', delta, "");
      cout << "delta = " << delta << ": Status: " << rstat[repstatus] << endl;
      if (repstatus == 0) {
        p.sync(XPRB_XPRS_SOLMIP);
        print_sol_opt(p);
        print_violated(allCtr);
      }
    }

    delete[] lrp;
    delete[] grp;
    delete[] lbp;
    delete[] ubp;
  }

  return 0;
}

void print_sol_opt(XPRBprob &p)
{
  cout << "  Total return: " << p.getObjVal() << endl;
  for (int s = 0; s<NSHARES; s++)
    if (buy[s].getSol() > 0.5)
      cout << "  " << SHARES_n[s] << " : " << frac[s].getSol() * 100 << "% (" << buy[s].getSol() << ")" << endl;
}

void print_violated(vector<XPRBctr> &ctrs)
{
  char *type = NULL;
  cout << " Violated (relaxed) constraints:" << endl;
  for (vector<XPRBctr>::iterator c = ctrs.begin(); c != ctrs.end(); c++) {
    double viol, slack = c->getSlack();
    switch (c->getType()) {
      case XPRB_E: viol = abs(slack); type = " ="; break;
      case XPRB_G: viol = slack;      type = ">="; break;
      case XPRB_L: viol = -slack;     type = "<="; break;
      default: cout << "  unexpected constraint type" << endl; continue;
    }
    if (viol > 1e-6) cout << "  " << type << " constraint " << c->getName() << " by " << -slack << endl;
  }
}