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

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

  Same model as in foliomip3.cpp.
  -- Infeasible model parameter values --
  -- Handling infeasibility through auxiliary variables --

  (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 "xprb_cpp.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;

#include "readfoliodata.c_"


int main(int argc, char **argv)
{
 int s,r,t;
 XPRBprob p("FolioMIP3");          // Initialize a new problem in BCL
 XPRBctr Risk,Return,Num,*MinReg, *MaxReg, *LimSec;
 XPRBexpr le, le2, Cap, LinkL, LinkU;
 XPRBvar *frac;                    // Fraction of capital used per share
 XPRBvar *buy;                     // 1 if asset is in portfolio, 0 otherwise

 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

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

// Limits on geographical distribution
 MinReg = new XPRBctr[NREGIONS];
 MaxReg = new XPRBctr[NREGIONS];
 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(le >= MINREG);
  MaxReg[r] = p.newCtr(le2 <= MAXREG);
 }

// Diversification across industry sectors
 LimSec = new XPRBctr[NTYPES];
 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(le <= MAXSEC);
 }

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

// Limit the total number of assets
 le=0;
 for(s=0;s<NSHARES;s++) le += buy[s];
 Num = p.newCtr(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
 p.setSense(XPRB_MAXIM);
 p.mipOptimize("");

 char *MIPSTATUS[] = {"not loaded", "not optimized", "LP optimized",
                      "unfinished (no solution)",
                      "unfinished (solution found)", "infeasible", "optimal",
                      "unbounded"};

 cout << "Problem status: " << MIPSTATUS[p.getMIPStat()] << endl;

 if (p.getMIPStat()==XPRB_MIP_INFEAS)
 {
  cout << "Original problem infeasible. Adding deviation variables" << endl;

  // Define deviation variables and add them to the constraints
  // to make problem solvable */
  XPRBvar devRisk, devNum, *devMinReg, *devMaxReg, *devSec;
  devRisk = p.newVar("devRisk");
  Risk -= devRisk;

  devMinReg = new XPRBvar[NREGIONS];
  devMaxReg = new XPRBvar[NREGIONS];
  for(r=0;r<NREGIONS;r++)
  {                                  /* Only allow small deviations */
   devMinReg[r] = p.newVar("devMinReg", XPRB_PL, 0, MAXREG/2);
   MinReg[r] += devMinReg[r];
   devMaxReg[r] = p.newVar("devMaxReg", XPRB_PL, 0, MAXREG/2);
   MaxReg[r] -= devMaxReg[r];
  }

  devSec = new XPRBvar[NTYPES];
  for(t=0;t<NTYPES;t++)
  {
   devSec[t] = p.newVar("devSec", XPRB_PL, 0, MAXSEC/2);
   LimSec[t] -= devSec[t];
  }

  devNum = p.newVar("devNum", XPRB_PL, 0, XPRB_INFINITY);
  Num -= devNum;


 /* Resolve the problem with penalty terms added to the objective */
  double penalty = -10;
  Return += devRisk * penalty;
  for(r=0;r<NREGIONS;r++)  Return += (devMinReg[r]+devMaxReg[r]) * penalty;
  for(t=0;t<NTYPES;t++)  Return += devSec[t] * penalty;
  Return += devNum * penalty;
  p.setObj(Return);                    /* Set the new objective function */

  p.mipOptimize("");

  if (p.getMIPStat()== XPRB_MIP_INFEAS)
  {
   cout << "No solution after relaxation" << endl;
   return 0;
  }
  else
  {
   cout << "Constraint violations:" << endl;
   cout << "  Constraint            Activity  Deviation  Bound(s)" << endl;
   printf("  Risk                    %1.2f\t  %1.2f\t  (%1.2f)\n",
          Risk.getAct()+devRisk.getSol(), devRisk.getSol(), MAXRISK);
   for(r=0;r<NREGIONS;r++)
   {
    double sumf=0;
    for(s=0;s<NSHARES;s++) if(LOC[r][s]>0) sumf+=frac[s].getSol();
    printf("  Region %-11s      %1.2f\t  %1.2f\t  (%g-%g)\n", REGIONS_n[r], sumf,
           devMaxReg[r].getSol()-devMinReg[r].getSol(),
           MINREG, MAXREG);
   }
   for(t=0;t<NTYPES;t++)
   {
    printf("  Sector %-14s   %1.2f\t  %1.2f\t  (%g)\n", TYPES_n[t],
           LimSec[t].getAct()+devSec[t].getSol(), devSec[t].getSol(),
           MAXSEC);
   }
   printf("  Number of assets        %1.2f\t  %1.2f\t  (%d)\n",
          Num.getAct()+devNum.getSol(), devNum.getSol(), MAXNUM);
  }
 }

// Solution printing
 cout << "Total return: " << p.getObjVal() << endl;
 for(s=0;s<NSHARES;s++)
  if(buy[s].getSol()>0.5)
   cout << SHARES_n[s] << ": " << frac[s].getSol()*100 << "% (" << buy[s].getSol()
        << ")" << endl;

 return 0;
}