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

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

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

#include <iostream>
#include "xprb_cpp.h"
#include "xprs.h"

using namespace std;
using namespace ::dashoptimization;

#define MAXNUM 4                   // Max. number of shares to be selected

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

void solveHeur(XPRBprob &p);

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

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;
 XPRBexpr Risk,Na,Return,Cap,Num;

 XPRBprob p("FolioMIPHeur");       // Initialize a new problem in BCL

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

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

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

// Minimum amount of North-American values
 for(s=0;s<NNA;s++) Na += frac[NA[s]];
 p.newCtr(Na >= 0.5);

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

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

// Linking the variables
 for(s=0;s<NSHARES;s++) p.newCtr(frac[s] <= buy[s]);

// Solve problem heuristically
 p.setSense(XPRB_MAXIM);
 solveHeur(p);

// Solve the problem
 p.mipOptimize("");

// Solution printing
 if(p.getMIPStat()==XPRB_MIP_SOLUTION || p.getMIPStat()==XPRB_MIP_OPTIMAL)
 {
  cout << "Exact solution: Total return: " << p.getObjVal() << endl;
  for(s=0;s<NSHARES;s++)
   cout << s << ": " << frac[s].getSol()*100 << "%" << endl;
 }
 else
  cout << "Heuristic solution is optimal." << endl;

 return 0;
}

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

 XPRSsetintcontrol(p.getXPRSprob(), XPRS_CUTSTRATEGY, 0);
                                   // Disable automatic cuts
 XPRSsetintcontrol(p.getXPRSprob(), XPRS_PRESOLVE, 0);
                                   // Switch presolve off
 XPRSgetdblcontrol(p.getXPRSprob(), XPRS_FEASTOL, &TOL);
                                   // Get feasibility tolerance

 p.mipOptimize("l");               // Solve the LP-relaxation
 basis=p.saveBasis();              // 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]=frac[s].getSol();        // Get the solution values of `frac'
  if(fsol[s] < TOL)  buy[s].setUB(0);
  else if(fsol[s] > 0.2-TOL)  buy[s].setLB(1);
 }

 p.mipOptimize("c");               // Solve the MIP-problem
 ifgsol=0;
 if(p.getMIPStat()==XPRB_MIP_SOLUTION || p.getMIPStat()==XPRB_MIP_OPTIMAL)
 {                                 // If an integer feas. solution was found
  ifgsol=1;
  solval=p.getObjVal();            // Get the value of the best solution
  cout << "Heuristic solution: Total return: " << p.getObjVal() << endl;
  for(s=0;s<NSHARES;s++)
   cout << s << ": " << frac[s].getSol()*100 << "%" << endl;
 }

// XPRSpostsolve(p.getXPRSprob());   // 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))
  {
   buy[s].setLB(0);
   buy[s].setUB(1);
  }

 p.loadBasis(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 */
 basis.reset();             // No need to store the saved basis any longer
 if(ifgsol==1)
  XPRSsetdblcontrol(p.getXPRSprob(), XPRS_MIPABSCUTOFF, solval+TOL);
                            // Set the cutoff to the best known solution

}