/********************************************************/
/*  Xpress-BCL C# Example Problems                      */
/*  ==============================                      */
/*                                                      */
/*  file xbfixbv.cs                                     */
/*  ```````````````                                     */
/*  Example for the use of Xpress-BCL                   */
/*  (Using the complete Coco Problem, as in xbcoco.cxx, */
/*  this program implements a binary fixing heuristic)  */
/*                                                      */
/*  (c) 2008 Fair Isaac Corporation                     */
/*      authors: S.Heipcke, D.Brett.                    */
/********************************************************/

using System;
using System.Text;
using System.IO;
using Optimizer;
using BCL;


namespace Examples
{

    public class TestAdvFixBV
    {
        const double TOL = 0,001;

        const int PHASE = 5;
        /* Phase = 4: Mines may open/closed freely; when closed save 20000 per month
        * Phase = 5: Once closed always closed; larger saving */

        const int NP = 2;            /* Number of products (p) */
        const int NF = 2;            /*           factories (f) */
        const int NR = 2;            /*           raw materials (r) */
        const int NT = 4;            /*           time periods (t) */

        /****DATA****/
        double[,] REV =      /* Unit selling price of product p in period t */
        {{400, 380, 405, 350},
        {410, 397, 412, 397}};
        double[,] CMAK =     /* Unit cost to make product p at factory f */
        {{150, 153},
        { 75,  68}};
        double[,] CBUY =     /* Unit cost to buy raw material r in period t */
        {{100,  98,  97, 100},
        {200, 195, 198, 200}};
        double[] COPEN =        /* Fixed cost of factory f being open for one period */
        {50000, 63000};
        double CPSTOCK = 2;   /* Unit cost to store any product p */
        double CRSTOCK = 1;   /* Unit cost to store any raw material r */
        double[,] REQ =      /* Requirement by unit of prod. p for raw material r */
        {{1.0, 0,5},
        {1.3, 0,4}};
        double[,] MXSELL =   /* Max. amount of p that can be sold in period t */
        {{650, 600, 500, 400},
        {600, 500, 300, 250}};
        double[] MXMAKE =       /* Max. amount factory f can make over all products */
        {400, 500};
        double MXRSTOCK = 300;  /* Max. amount of r that can be stored each f and t */
        double[,] PSTOCK0 =  /* Initial product p stock level at factory f */
        {{50, 100},
        {50,  50}};
        double[,] RSTOCK0 =  /* Initial raw material r stock level at factory f*/
        {{100, 150},
        { 50, 100}};

        XPRBprob pb = new XPRBprob("Coco");    /* Initialize a new problem in BCL */

        /****VARIABLES****/
        XPRBvar[,] openm = new XPRBvar[NF,NT];

        /***********************************************************************/

        void modCoco()
        {
            XPRBvar[,,] make = new XPRBvar[NP,NF,NT];
            XPRBvar[,,] sell = new XPRBvar[NP,NF,NT];
            XPRBvar[,,] pstock = new XPRBvar[NP,NF,NT+1];
            XPRBvar[,,] buy = new XPRBvar[NR,NF,NT];
            XPRBvar[,,] rstock = new XPRBvar[NR,NF,NT+1];
            XPRBexpr lobj, lc;
            int p,f,r,t;

            /****VARIABLES****/
            for(p=0;p<NP;p++)
                for(f=0;f<NF;f++)
                {
                    for(t=0;t<NT;t++)
                    {
                        make[p,f,t] = pb.newVar("make_p" + (p+1) + "_f" + (f+1));
                        /* Amount of prod. p to make at factory f in period t */
                        sell[p,f,t] = pb.newVar("sell_p" + (p+1) + "_f" + (f+1));
                        /* Amount of prod. p sold from factory f in period t */
                    }
                    for(t=0;t<NT+1;t++)
                        pstock[p,f,t] = pb.newVar("pstock_p" + (p+1) + "_f" + (f+1));
                    /* Stock level of prod. p at factory f at start of period t */
                }

            for(r=0;r<NR;r++)
                for(f=0;f<NF;f++)
                {
                    for(t=0;t<NT;t++)
                        buy[r,f,t] = pb.newVar("buy_r" + (r+1) + "_f" + (f+1));
                    /* Amount of raw material r bought for factory f in period t */
                    for(t=0;t<NT+1;t++)
                        rstock[r,f,t] = pb.newVar("rstock_r" + (r+1) + "_f" + (f+1));
                    /* Stock level of raw mat. r at factory f at start of per. t */
                }

            for(f=0;f<NF;f++)
                for(t=0;t<NT;t++)
                    openm[f,t] = pb.newVar("open_f" + (f+1), BCLconstant.XPRB_BV);
                    /* 1 if factory f is open in period t, else 0 */

            /****OBJECTIVE****/
            lobj = new XPRBexpr();
            for(f=0;f<NF;f++)        /* Objective: maximize total profit */
            {
                for(p=0;p<NP;p++)
                {
                    for(t=0;t<NT;t++)
                        lobj += REV[p,t]*sell[p,f,t] - CMAK[p,f]*make[p,f,t];
                    for(t=1;t<NT+1;t++)  lobj -= CPSTOCK*pstock[p,f,t];
                }
                if(PHASE==4)
                    for(t=0;t<NT;t++)  lobj += (20000-COPEN[f])*openm[f,t];
                else if(PHASE==5)
                    for(t=0;t<NT;t++)  lobj -= COPEN[f]*openm[f,t];
                for(r=0;r<NR;r++)
                {
                    for(t=0;t<NT;t++)  lobj -= CBUY[r,t]*buy[r,f,t];
                    for(t=1;t<NT+1;t++)  lobj -= CRSTOCK*rstock[r,f,t];
                }
            }
            pb.setObj(pb.newCtr("OBJ",lobj));   /* Set objective function */

            /****CONSTRAINTS****/
            for(p=0;p<NP;p++)        /* Product stock balance */
                for(f=0;f<NF;f++)
                    for(t=0;t<NT;t++)
                        pb.newCtr("PBal", pstock[p,f,t]+make[p,f,t] == sell[p,f,t]+pstock[p,f,t+1] );

            for(r=0;r<NR;r++)        /* Raw material stock balance */
                for(f=0;f<NF;f++)
                    for(t=0;t<NT;t++)
                    {
                        lc = new XPRBexpr(0);
                        for(p=0;p<NP;p++)  lc += REQ[p,r]*make[p,f,t];
                        pb.newCtr("RBal", rstock[r,f,t]+buy[r,f,t] == lc+rstock[r,f,t+1]);
                    }

            for(p=0;p<NP;p++)
                for(t=0;t<NT;t++)
                {                       /* Limit on the amount of product p to be sold */
                    lc = new XPRBexpr(0);
                    for(f=0;f<NF;f++)  lc += sell[p,f,t];
                    pb.newCtr("MxSell", lc <= MXSELL[p,t]);
                }

            for(f=0;f<NF;f++)
                for(t=0;t<NT;t++)
                {                       /* Capacity limit at factory f */
                    lc = new XPRBexpr(0);
                    for(p=0;p<NP;p++)  lc += make[p,f,t];
                    pb.newCtr("MxMake", lc <= MXMAKE[f]*openm[f,t]);
                }

            for(f=0;f<NF;f++)
                for(t=1;t<NT+1;t++)
                {                       /* Raw material stock limit */
                    lc = new XPRBexpr(0);
                    for(r=0;r<NR;r++) lc += rstock[r,f,t];
                    pb.newCtr("MxRStock", lc <= MXRSTOCK);
                }

            if(PHASE==5)
            for(f=0;f<NF;f++)
               for(t=0;t<NT-1;t++)    /* Once closed, always closed */
                  pb.newCtr("Closed", openm[f,t+1] <= openm[f,t]);

            /****BOUNDS****/
            for(p=0;p<NP;p++)
               for(f=0;f<NF;f++)
                  pstock[p,f,0].fix(PSTOCK0[p,f]);    /* Initial product levels */

            for(r=0;r<NR;r++)
                for(f=0;f<NF;f++)
                   rstock[r,f,0].fix(RSTOCK0[r,f]);    /* Initial raw mat. levels */

            if(PHASE<=3)
            for(f=0;f<NF;f++)
                for(t=0;t<NT;t++)
                    openm[f,t].fix(1);
        }

        /**************************************************************************/
        /*  Solution heuristic:                                                   */
        /*    solve the LP and save the basis                                     */
        /*    fix all open variables which are integer feasible at the relaxation */
        /*    solve the MIP and save the best solution value                      */
        /*    reset all variables to their original bounds                        */
        /*    load the saved basis                                                */
        /*    solve the MIP using the solution value found previously as cutoff   */
        /**************************************************************************/
        void solveCoco()
        {
            XPRBbasis basis;
            int f,t, ifgsol;
            double solval = 0;
            double[,] osol = new double[NF,2];
            XPRSprob xprsp = pb.getXPRSprob();

            /****SOLVING + OUTPUT****/
            xprsp.CutStrategy = 0;
            /* Disable automatic cuts - we use a heuristic */
            xprsp.Presolve = 0;
            /* Switch presolve off */
            pb.setSense(BCLconstant.XPRB_MAXIM);  /* Choose the sense of the optimization */
            pb.lpOptimize();              /* Solve the LP-problem */
            basis=pb.saveBasis();      /* Save the current basis */

            /* For t=0,1 get the solution values of the open variables: */
            for(f=0;f<NF;f++)
                for(t=0;t<2;t++)
                {
                    osol[f,t]=openm[f,t].getSol();
                    /* Fix all variables which are integer feasible: */
                    if(osol[f,t] < TOL)  openm[f,t].setUB(0.0);
                    else if(1-osol[f,t] < TOL)  openm[f,t].setLB(1);
                }

            pb.mipOptimize();             /* Solve the MIP-problem */
            ifgsol=0;
            if(pb.getMIPStat()==4 || pb.getMIPStat()==6)
            {                          /* If a global solution was found */
                ifgsol=1;
                solval=pb.getObjVal();    /* Get the value of the best solution */
            }

            xprsp.InitGlobal(); /* Re-initialize the global search */

            /* Reset variables to their original bounds: */
            for(f=0;f<NF;f++)
                for(t=0;t<2;t++)
                    if((osol[f,t] < TOL) || (1-osol[f,t] < TOL))
                    {
                        openm[f,t].setLB(0.0);
                        openm[f,t].setUB(1);
                    }

            pb.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)
                xprsp.MIPAbsCutoff = solval;
            /* Set the cutoff to the best known solution */
            pb.mipOptimize();             /* Solve the MIP-problem */

            if(pb.getMIPStat()==4 || pb.getMIPStat()==6)
                System.Console.WriteLine("Best integer solution: " + pb.getObjVal());
            else
                System.Console.WriteLine("Best integer solution: " + solval);
        }

        /***********************************************************************/

        public static void Main()
        {
           XPRB.init();
           TestAdvFixBV TestInstance = new TestAdvFixBV();

            TestInstance.modCoco();                 /* Model the problem */
            TestInstance.solveCoco();               /* Solve the problem */

            return;
        }
    }
}