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

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

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

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include "xprb.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;

#include "readfoliodata.c_"

int main(int argc, char **argv)
{
 int s,r,t;
 double sumf, penalty;
 XPRBprob prob;
 XPRBctr Risk,Return,Cap,Num;
 XPRBctr *MinReg, *MaxReg, *LimSec, LinkL, LinkU;
 XPRBvar *frac;                   /* Fraction of capital used per share */
 XPRBvar *buy;                    /* 1 if asset is in portfolio, 0 otherwise */
 XPRBvar devRisk, devNum, *devMinReg, *devMaxReg, *devSec;
 char *MIPSTATUS[] = {"not loaded", "not optimized", "LP optimized",
                      "unfinished (no solution)",
                      "unfinished (solution found)", "infeasible", "optimal",
                      "unbounded"};

 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 */

/* 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 */
 MinReg = (XPRBctr*)malloc(NREGIONS*sizeof(XPRBctr));
 MaxReg = (XPRBctr*)malloc(NREGIONS*sizeof(XPRBctr));
 for(r=0;r<NREGIONS;r++)
 {
  MinReg[r] = XPRBnewctr(prob, "MinReg", XPRB_G);
  MaxReg[r] = XPRBnewctr(prob, "MaxReg", 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 */
 LimSec = (XPRBctr*)malloc(NTYPES*sizeof(XPRBctr));
 for(t=0;t<NTYPES;t++)
 {
  LimSec[t] = XPRBnewctr(prob, "LimSec", 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 = XPRBnewctr(prob, "LinkU", XPRB_L);
  XPRBaddterm(LinkU, frac[s], 1);
  XPRBaddterm(LinkU, buy[s], -MAXVAL);
  LinkL = XPRBnewctr(prob, "LinkL", XPRB_G);
  XPRBaddterm(LinkL, frac[s], 1);
  XPRBaddterm(LinkL, buy[s], -MINVAL);
 }

/* Solve the problem */
 XPRBsetsense(prob, XPRB_MAXIM);
 XPRBmipoptimize(prob, "");

 printf("Problem status: %s\n", MIPSTATUS[XPRBgetmipstat(prob)]);

 if (XPRBgetmipstat(prob)==XPRB_MIP_INFEAS)
 {
  printf("Original problem infeasible. Adding deviation variables\n");

  /* Define deviation variables and add them to the constraints
     to make problem solvable */
  devRisk = XPRBnewvar(prob, XPRB_PL, "devRisk", 0, XPRB_INFINITY);
  XPRBaddterm(Risk, devRisk, -1);

  devMinReg = (XPRBvar*)malloc(NREGIONS*sizeof(XPRBvar));
  devMaxReg = (XPRBvar*)malloc(NREGIONS*sizeof(XPRBvar));
  for(r=0;r<NREGIONS;r++)
  {                                  /* Only allow small deviations */
   devMinReg[r] = XPRBnewvar(prob, XPRB_PL, "devMinReg", 0, MAXREG/2);
   XPRBaddterm(MinReg[r], devMinReg[r], 1);
   devMaxReg[r] = XPRBnewvar(prob, XPRB_PL, "devMaxReg", 0, MAXREG/2);
   XPRBaddterm(MaxReg[r], devMaxReg[r], -1);
  }

  devSec = (XPRBvar*)malloc(NTYPES*sizeof(XPRBvar));
  for(t=0;t<NTYPES;t++)
  {
   devSec[t] = XPRBnewvar(prob, XPRB_PL, "devSec", 0, MAXSEC/2);
   XPRBaddterm(LimSec[t], devSec[t], -1);
  }

  devNum = XPRBnewvar(prob, XPRB_PL, "devNum", 0, XPRB_INFINITY);
  XPRBaddterm(Num, devNum, -1);


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

  XPRBmipoptimize(prob, "");

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

/* Solution printing */
 printf("Total return: %g\n", XPRBgetobjval(prob));
 for(s=0;s<NSHARES;s++)
  if(XPRBgetsol(buy[s])>0.5)
   printf("  %d : %g%% (%g)\n", s, XPRBgetsol(frac[s])*100, XPRBgetsol(buy[s]));

 return 0;
}