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

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

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

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

import java.io.*;
import java.lang.*;
import java.util.*;
import java.text.DecimalFormat;
import com.dashoptimization.*;

public class folioinfeas
{
 static final String DATAFILE = "folio10.cdat";

 static final int MAXNUM = 4;          /* Max. number of different assets */
 static final double MAXRISK = 1.0/3;  /* Max. investment into high-risk values */
 static final double MINREG = 0.3;     /* Min. investment per geogr. region */
 static final double MAXREG = 0.5;     /* Max. investment per geogr. region */
 static final double MAXSEC = 0.15;    /* Max. investment per ind. sector */
 static final double MAXVAL = 0.2;     /* Max. investment per share */
 static final double MINVAL = 0.1;     /* Min. investment per share */


 static int NSHARES;               /* Number of shares */
 static int NRISK;                 /* Number of high-risk shares */
 static int NREGIONS;              /* Number of geographical regions */
 static int NTYPES;                /* Number of share types */

 static double[] RET;              /* Estimated return in investment  */
 static int[] RISK;                /* High-risk values among shares */
 static boolean LOC[][];           /* Geogr. region of shares */
 static boolean SEC[][];           /* Industry sector of shares */

 static String SHARES_n[];
 static String REGIONS_n[];
 static String TYPES_n[];

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

 private static DecimalFormat form = new DecimalFormat ("0.00");
 
 public static void main(String[] args) throws IOException
 { 
  int s,t,r;
  XPRB bcl;
  XPRBprob p;
  XPRBexpr LinkL, LinkU, le, le2;
  XPRBctr Risk,Return,Cap,Num;
  XPRBctr[] MinReg, MaxReg, LimSec;
  XPRBvar[] frac;                  /* Fraction of capital used per share */
  XPRBvar[] buy;                   /* 1 if asset is in portfolio, 0 otherwise */

  try
  {
   readData();                     /* Read data from file */
  }
  catch(IOException e)
  {
   System.err.println(e.getMessage());
   System.exit(1);
  }

  bcl = new XPRB();                /* Initialize BCL */
  p = bcl.newProb("FolioMIP3inf"); /* Create a new problem in BCL */

/* Create the decision variables */
  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 */
  le = new XPRBexpr();
  for(s=0;s<NSHARES;s++) le.add(frac[s].mul(RET[s])); 
  Return = p.newCtr("Return", le);
  p.setObj(le);                    /* Set the objective function */

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

/* Limits on geographical distribution */
  MinReg = new XPRBctr[NREGIONS];
  MaxReg = new XPRBctr[NREGIONS];
  for(r=0;r<NREGIONS;r++)
  {
   le = new XPRBexpr();
   le2 = new XPRBexpr();
   for(s=0;s<NSHARES;s++)
    if(LOC[r][s])
    {
     le.add(frac[s]); 
     le2.add(frac[s]);
    } 
    MinReg[r] = p.newCtr("MinReg", le.gEql(MINREG));
    MaxReg[r] = p.newCtr("MaxReg", le2.lEql(MAXREG));
  } 

/* Diversification across industry sectors */
  LimSec = new XPRBctr[NTYPES];
  for(t=0;t<NTYPES;t++)
  {
   le = new XPRBexpr();
   for(s=0;s<NSHARES;s++)
    if(SEC[t][s]) le.add(frac[s]); 
   LimSec[t] = p.newCtr("LimSec", le.lEql(MAXSEC));
  } 
 
/* Spend all the capital */
  le = new XPRBexpr();
  for(s=0;s<NSHARES;s++) le.add(frac[s]); 
  Cap = p.newCtr("Cap", le.eql(1));
 
/* Limit the total number of assets */
  le = new XPRBexpr();
  for(s=0;s<NSHARES;s++) le.add(buy[s]);
  Num = p.newCtr("Num", le.lEql(MAXNUM));

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


/* Solve the problem */
  p.setSense(XPRB.MAXIM);
  p.mipOptimize("");

  System.out.println("Problem status: " + MIPSTATUS[p.getMIPStat()]);

 if (p.getMIPStat()==XPRB.MIP_INFEAS)
 {
  System.out.println("Original problem infeasible. Adding deviation variables");

  XPRBvar devRisk, devNum; 
  XPRBvar [] devMinReg, devMaxReg, devSec;

  /* Define deviation variables and add them to the constraints 
     to make problem solvable */
  devRisk = p.newVar("devRisk");
  Risk.add(devRisk.mul(-1));

  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].add(devMinReg[r]);
   devMaxReg[r] = p.newVar("devMaxReg", XPRB.PL, 0, MAXREG/2);
   MaxReg[r].add(devMaxReg[r].mul(-1));
  }

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

  devNum = p.newVar("devNum");
  Num.add(devNum.mul(-1));


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

  p.mipOptimize("");

  if (p.getMIPStat()==XPRB.MIP_INFEAS)
  {
   System.out.println("No solution after relaxation");
   return;
  }
  else
  {
   System.out.println("Constraint violations:");
   System.out.println("  Constraint            Activity  Deviation  Bound(s)");
   System.out.println("  Risk                    " +
    form.format(Risk.getAct()+devRisk.getSol()) + "\t  " + 
    form.format(devRisk.getSol()) +"\t  (" + form.format(MAXRISK) + ")");
   for(r=0;r<NREGIONS;r++) 
   {
    double sumf=0;
    for(s=0;s<NSHARES;s++) if(LOC[r][s]) sumf+=frac[s].getSol();
    StringBuffer sb = new StringBuffer(REGIONS_n[r]+"          ");
    sb.setLength(11);
    System.out.println("  Region " + sb.toString() + "      " +
     form.format(sumf) + "\t  " +
     form.format(devMaxReg[r].getSol()-devMinReg[r].getSol()) + "\t  (" + 
     MINREG + "-" + MAXREG + ")");
   }
   for(t=0;t<NTYPES;t++)
   {
    StringBuffer sb = new StringBuffer(TYPES_n[t]+"          ");
    sb.setLength(14);
    System.out.println("  Sector " + sb.toString() + "   " +
     form.format(LimSec[t].getAct()+devSec[t].getSol()) + "\t  " + 
     form.format(devSec[t].getSol()) + "\t  (" + MAXSEC + ")");
   }
   System.out.println("  Number of assets        " +
    form.format(Num.getAct()+devNum.getSol()) + "\t  " +
    form.format(devNum.getSol()) + "\t  (" +MAXNUM + ")");
  }
 }

 
/* Solution printing */
  System.out.println("Total return: " + form.format(p.getObjVal()));
  for(s=0;s<NSHARES;s++) 
   if(buy[s].getSol()>0.5)
    System.out.println(s + ": " + frac[s].getSol()*100 + "% (" +
     buy[s].getSol() + ")");  

/* Delete the problem */
  p.finalize();              
  p=null;

 }

/***********************Data input routines***************************/

/***************************/
/* Input a list of strings */
/***************************/
 private static String [] read_str_list(StreamTokenizer st) throws IOException
 {
  LinkedList<String> l=new LinkedList<String>();

  st.nextToken();                  /* Skip ':' */
  while(st.nextToken()==st.TT_WORD)
  {
   l.addLast(st.sval);
  }

  String a[]=new String[l.size()];
  l.toArray(a);
  return a;
 }

/************************/
/* Input a list of ints */
/************************/
 private static int [] read_int_list(StreamTokenizer st) throws IOException
 {
  LinkedList<Integer> l=new LinkedList<Integer>();

  st.nextToken();                  /* Skip ':' */
  while(st.nextToken()==st.TT_NUMBER)
  {
   l.addLast((int)st.nval);
  }

  int a[]=new int[l.size()];
  for(int i=0;i<l.size();i++)
   a[i]=((Integer)l.get(i)).intValue();
  return a;
 }

/****************************/
/* Input a table of doubles */
/****************************/
 private static void read_dbl_table(StreamTokenizer st,double tbl[]) throws IOException
 {
  int n=0;

  st.nextToken();                  /* Skip ':' */
  while(st.nextToken()==st.TT_NUMBER)
  {
   tbl[n++]=st.nval;
  }
 }

/************************************/
/* Input a sparse table of booleans */
/************************************/
 private static boolean [][] read_bool_table(StreamTokenizer st,int nrow,int ncol) throws IOException
 {
  int i;
  boolean tbl[][]=new boolean[nrow][ncol];

  st.nextToken();                  /* Skip ':' */
  for(int r=0;r<nrow;r++)
  {
   while(st.nextToken()==st.TT_NUMBER)
    tbl[r][(int)st.nval]=true;
  }
  return tbl;
 }

 private static void readData() throws IOException
 {
  int s;
  FileReader datafile=null;
  StreamTokenizer st=null;

  datafile=new FileReader(DATAFILE);/* Open the data file */
  st=new StreamTokenizer(datafile); /* Initialize the stream tokenizer */
  st.commentChar('!');              /* Use the character '!' for comments */
  st.eolIsSignificant(false);       /* Return end-of-line character */
  st.parseNumbers();                /* Read numbers as numbers (not strings) */

  while ( st.nextToken() == st.TT_WORD )
  {
   if ( st.sval.equals("SHARES") && NSHARES==0)
    { SHARES_n=read_str_list(st); NSHARES=SHARES_n.length; }
   else
   if ( st.sval.equals("REGIONS") && NREGIONS==0)
    { REGIONS_n=read_str_list(st); NREGIONS=REGIONS_n.length; }
   else
   if ( st.sval.equals("TYPES") && NTYPES==0)
    { TYPES_n=read_str_list(st); NTYPES=TYPES_n.length; }
   else
   if ( st.sval.equals("RISK") && NRISK==0)
    { RISK=read_int_list(st); NRISK=RISK.length; }
   else
   if ( st.sval.equals("RET") && NSHARES>0)
    {
     RET=new double[NSHARES];
     read_dbl_table(st,RET);
    }
   else
   if ( st.sval.equals("LOC") && NSHARES>0 && NREGIONS>0)
     LOC=read_bool_table(st,NREGIONS,NSHARES);
   else
   if ( st.sval.equals("SEC") && NSHARES>0 && NTYPES>0)
    SEC=read_bool_table(st,NTYPES,NSHARES);
   else
    break;
  }

/*
  for(int i=0;i<NREGIONS;i++) {
   for(int j=0;j<NSHARES;j++)
    System.out.print(" "+LOC[i][j]);
   System.out.println();
  }
*/
  datafile.close();
 }
}