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

  file folioheur.c
  ````````````````
  Modeling a small MIP problem
  to perform portfolio optimization.
   -- Heuristic solution --

  (c) 2008-2024 Fair Isaac Corporation
      author: S.Heipcke, Dec. 2003, rev. Mar. 2011
********************************************************/

#include <stdio.h>
#include "xprb.h"
#include "xprs.h"

#define MAXNUM 4                  /* Max. number of different assets */

#define NSHARES 10                /* Number of shares */
#define NRISK 5                   /* Number of high-risk shares */
#define NNA 4                     /* Number of North-American shares */

void solveHeur(void);

double RET[] = {5,17,26,12,8,9,7,6,31,21};  /* Estimated return in investment */
int RISK[] = {1,2,3,8,9};         /* High-risk values among shares */
int NA[] = {0,1,2,3};             /* Shares issued in N.-America */

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

int main(int argc, char **argv)
{
 int s;
 XPRBctr Risk,Na,Return,Cap,Num;

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

/* Create the decision variables (including upper bounds for `frac') */
 for(s=0;s<NSHARES;s++)
 {
  frac[s] = XPRBnewvar(prob, XPRB_PL, "frac", 0, 0.3);
  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, 1.0/3);

/* Minimum amount of North-American values */
 Na = XPRBnewctr(prob, "NA", XPRB_G);
 for(s=0;s<NNA;s++) XPRBaddterm(Na, frac[NA[s]], 1);
 XPRBaddterm(Na, NULL, 0.5);

/* 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++) XPRBnewprec(prob, "Link", frac[s], 0, buy[s]);

/* Solve problem heuristically */
 XPRBsetsense(prob, XPRB_MAXIM);
 solveHeur();

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

/* Solution printing */
 if(XPRBgetmipstat(prob)==XPRB_MIP_SOLUTION ||
    XPRBgetmipstat(prob)==XPRB_MIP_OPTIMAL)
 {
  printf("Exact solution: Total return: %g\n", XPRBgetobjval(prob));
  for(s=0;s<NSHARES;s++)
   printf(" %d : %g%%\n", s, XPRBgetsol(frac[s])*100);
 }
 else
  printf("Heuristic solution is optimal.\n");

 return 0;
}

void solveHeur(void)
{
 XPRBbasis basis;
 int s, ifgsol;
 double solval, fsol[NSHARES],TOL;

 XPRSsetintcontrol(XPRBgetXPRSprob(prob), XPRS_CUTSTRATEGY, 0);
                                   /* Disable automatic cuts */
 XPRSsetintcontrol(XPRBgetXPRSprob(prob), XPRS_PRESOLVE, 0);
                                   /* Switch presolve off */
 XPRSgetdblcontrol(XPRBgetXPRSprob(prob), XPRS_FEASTOL, &TOL);
                                   /* Get feasibility tolerance */

 XPRBmipoptimize(prob, "l");       /* Solve the LP-problem */
 basis=XPRBsavebasis(prob);        /* Save the current basis */

/* Fix all variables `buy' for which `frac' is at 0 or at a relatively
   large value */
 for(s=0;s<NSHARES;s++)
 {
  fsol[s]=XPRBgetsol(frac[s]);     /* Get the solution values of `frac' */
  if(fsol[s] < TOL)  XPRBsetub(buy[s], 0);
  else if(fsol[s] > 0.2-TOL)  XPRBsetlb(buy[s], 1);
 }

 XPRBmipoptimize(prob, "c");       /* Continue solving the MIP-problem */
 ifgsol=0;
 if(XPRBgetmipstat(prob)==XPRB_MIP_SOLUTION ||
    XPRBgetmipstat(prob)==XPRB_MIP_OPTIMAL)
 {                                 /* If an integer feas. solution was found */
  ifgsol=1;
  solval=XPRBgetobjval(prob);      /* Get the value of the best solution */
  printf("Heuristic solution: Total return: %g\n", XPRBgetobjval(prob));
  for(s=0;s<NSHARES;s++)
   printf(" %d : %g%%\n", s, XPRBgetsol(frac[s])*100);
 }

 XPRSpostsolve(XPRBgetXPRSprob(prob));  /* Re-initialize the global search */

/* Reset variables to their original bounds */
 for(s=0;s<NSHARES;s++)
  if((fsol[s] < TOL) || (fsol[s] > 0.2-TOL))
  {
   XPRBsetlb(buy[s], 0);
   XPRBsetub(buy[s], 1);
  }

 XPRBloadbasis(basis);      /* Load the saved basis: bound changes are
                               immediately passed on from BCL to the
                               Optimizer if the problem has not been modified
                               in any other way, so that there is no need to
                               reload the matrix */
 XPRBdelbasis(basis);       /* No need to store the saved basis any longer */
 if(ifgsol==1)
  XPRSsetdblcontrol(XPRBgetXPRSprob(prob), XPRS_MIPABSCUTOFF, solval+TOL);
                            /* Set the cutoff to the best known solution */
}