/********************************************************
  BCL Example Problems
  ====================

  file xbfixbvls.c
  ````````````````
  Using the complete Coco Problem, as in xbcoco.c,
  this program implements a binary fixing heuristic.
  - Model version saving and loading a start solution -

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

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

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

#define NP 2            /* Number of products (p) */
#define NF 2            /*           factories (f) */
#define NR 2            /*           raw materials (r) */
#define NT 4            /*           time periods (t) */

#define TOL 5.0E-4

/****DATA****/
double REV[][NT] =      /* Unit selling price of product p in period t */
        {{400, 380, 405, 350},
         {410, 397, 412, 397}};
double CMAK[][NF] =     /* Unit cost to make product p at factory f */
        {{150, 153},
         { 75,  68}};
double CBUY[][NT] =     /* 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.0;   /* Unit cost to store any product p */
double CRSTOCK = 1.0;   /* Unit cost to store any raw material r */
double REQ[][NR] =      /* Requirement by unit of prod. p for raw material r */
        {{1.0, 0.5},
         {1.3, 0.4}};
double MXSELL[][NT] =   /* 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[][NF] =  /* Initial product p stock level at factory f */
        {{50, 100},
         {50,  50}};
double RSTOCK0[][NF] =  /* Initial raw material r stock level at factory f*/
        {{100, 150},
         { 50, 100}};

/****VARIABLES****/
XPRBvar openm[NF][NT];

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

void modcoco(XPRBprob prob)
{
 XPRBvar make[NP][NF][NT], sell[NP][NF][NT], pstock[NP][NF][NT+1],
       buy[NR][NF][NT], rstock[NR][NF][NT+1];
 XPRBctr ctr;
 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]=XPRBnewvar(prob,XPRB_PL, XPRBnewname("make_p%d_f%d",p+1,f+1),
        0, XPRB_INFINITY);
        /* Amount of prod. p to make at factory f in period t */
    sell[p][f][t]=XPRBnewvar(prob,XPRB_PL, XPRBnewname("sell_p%d_f%d",p+1,f+1),
        0, XPRB_INFINITY);
        /* Amount of prod. p sold from factory f in period t */
   }
   for(t=0;t<NT+1;t++)
    pstock[p][f][t]=XPRBnewvar(prob,XPRB_PL, XPRBnewname("pstock_p%d_f%d",p+1,f+1),
        0, XPRB_INFINITY);
        /* 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]=XPRBnewvar(prob,XPRB_PL, XPRBnewname("buy_r%d_f%d",r+1,f+1),
        0, XPRB_INFINITY);
        /* Amount of raw material r bought for factory f in period t */
   for(t=0;t<NT+1;t++)
    rstock[r][f][t]=XPRBnewvar(prob,XPRB_PL, XPRBnewname("rstock_r%d_f%d",r+1,f+1),
        0, XPRB_INFINITY);
        /* 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]=XPRBnewvar(prob,XPRB_BV, XPRBnewname("open_f%d",f+1), 0, 1);
        /* 1 if factory f is open in period t, else 0 */

/****OBJECTIVE****/
 ctr = XPRBnewctr(prob,"OBJ",XPRB_N);
 for(f=0;f<NF;f++)                  /* Objective: maximize total profit */
 {
  for(p=0;p<NP;p++)
  {
   for(t=0;t<NT;t++)
   {
    XPRBaddterm(ctr, sell[p][f][t], REV[p][t]);
    XPRBaddterm(ctr, make[p][f][t], -1*CMAK[p][f]);
   }
   for(t=1;t<NT+1;t++)
    XPRBaddterm(ctr, pstock[p][f][t], -1*CPSTOCK);
  }
  if(PHASE==4)
   for(t=0;t<NT;t++)
    XPRBaddterm(ctr, openm[f][t], 20000-COPEN[f]);
  else if(PHASE==5)
   for(t=0;t<NT;t++)
    XPRBaddterm(ctr, openm[f][t], -1*COPEN[f]);
  for(r=0;r<NR;r++)
  {
   for(t=0;t<NT;t++)
    XPRBaddterm(ctr, buy[r][f][t], -1*CBUY[r][t]);
   for(t=1;t<NT+1;t++)
    XPRBaddterm(ctr, rstock[r][f][t], -1*CRSTOCK);
  }
 }
 XPRBsetobj(prob,ctr);                       /* Select objective function */

/****CONSTRAINTS****/
 for(p=0;p<NP;p++)
  for(f=0;f<NF;f++)
   for(t=0;t<NT;t++)
   {
    ctr=XPRBnewctr(prob,"PBal",XPRB_E);      /* Product stock balance */
    XPRBaddterm(ctr, pstock[p][f][t], 1);
    XPRBaddterm(ctr, make[p][f][t], 1);
    XPRBaddterm(ctr, sell[p][f][t], -1);
    XPRBaddterm(ctr, pstock[p][f][t+1], -1);
   }

 for(r=0;r<NR;r++)
  for(f=0;f<NF;f++)
   for(t=0;t<NT;t++)
   {
    ctr=XPRBnewctr(prob,"RBal",XPRB_E);      /* Raw material stock balance */
    XPRBaddterm(ctr, rstock[r][f][t], 1);
    XPRBaddterm(ctr, buy[r][f][t], 1);
    XPRBaddterm(ctr, rstock[r][f][t+1], -1);
    for(p=0;p<NP;p++)
     XPRBaddterm(ctr, make[p][f][t], -1*REQ[p][r]);
   }

 for(p=0;p<NP;p++)
  for(t=0;t<NT;t++)
  {         /* Limit on the amount of product p to be sold */
   ctr=XPRBnewctr(prob,"MxSell",XPRB_L);
   for(f=0;f<NF;f++)
    XPRBaddterm(ctr, sell[p][f][t], 1);
   XPRBaddterm(ctr, NULL, MXSELL[p][t]);
  }

 for(f=0;f<NF;f++)
  for(t=0;t<NT;t++)
  {
   ctr=XPRBnewctr(prob,"MxMake",XPRB_L);     /* Capacity limit at factory f */
   for(p=0;p<NP;p++)
    XPRBaddterm(ctr, make[p][f][t], 1);
   XPRBaddterm(ctr, openm[f][t], -1*MXMAKE[f]);
  }

 for(f=0;f<NF;f++)
  for(t=1;t<NT+1;t++)
  {
   ctr=XPRBnewctr(prob,"MxRStock",XPRB_L);   /* Raw material stock limit */
   for(r=0;r<NR;r++)
    XPRBaddterm(ctr, rstock[r][f][t], 1);
   XPRBaddterm(ctr, NULL, MXRSTOCK);
  }

 if(PHASE==5)
  for(f=0;f<NF;f++)
   for(t=0;t<NT-1;t++)
   {
    ctr=XPRBnewctr(prob,"Closed",XPRB_L);    /* Once closed, always closed */
    XPRBaddterm(ctr, openm[f][t+1], 1);
    XPRBaddterm(ctr, openm[f][t], -1);
   }

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

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

/**************************************************************************/
/*  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(XPRBprob prob)
{
 XPRBbasis basis;
 int f,t, ifgsol, ncol, status;
 double solval, osol[NF][2];
 double *allsol;

/****SOLVING + OUTPUT****/
 XPRSsetintcontrol(XPRBgetXPRSprob(prob), XPRS_CUTSTRATEGY, 0);
                              /* Disable automatic cuts - we use a heuristic */
 XPRSsetintcontrol(XPRBgetXPRSprob(prob), XPRS_PRESOLVE, 0);
 XPRBsetsense(prob,XPRB_MAXIM);  /* Choose the sense of the optimization */
 XPRBmipoptimize(prob,"l");      /* Solve the LP-problem */
 basis=XPRBsavebasis(prob);      /* 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]=XPRBgetsol(openm[f][t]);
      /* Fix all variables which are integer feasible: */
   if(osol[f][t] < TOL)  XPRBsetub(openm[f][t], 0.0);
   else if(1-osol[f][t] < TOL)  XPRBsetlb(openm[f][t], 1.0);
  }

 XPRBmipoptimize(prob,"c");    /* Continue solving the MIP-problem */
 ifgsol=0;
 if(XPRBgetmipstat(prob)==XPRS_MIP_SOLUTION ||
    XPRBgetmipstat(prob)==XPRS_MIP_OPTIMAL)
 {                             /* If a global solution was found */
  ifgsol=1;
  solval=XPRBgetobjval(prob);  /* Get the value of the best solution */

                               /* Save the MIP solution values */
  XPRSgetintattrib(XPRBgetXPRSprob(prob), XPRS_ORIGINALCOLS, &ncol);
  allsol = calloc(ncol,sizeof(double));
  for(f=0;f<NF;f++)
   for(t=0;t<NT;t++)
    allsol[XPRBgetcolnum(openm[f][t])] = XPRBgetsol(openm[f][t]);
 }

/* XPRSpostsolve(XPRBgetXPRSprob(prob));  / * 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))
   {
    XPRBsetlb(openm[f][t], 0.0);
    XPRBsetub(openm[f][t], 1.0);
   }

 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)
 {                            /* Load the previously saved MIP solution */
  status = XPRBloadmipsol(prob, allsol, ncol, 1);
  if(status!=0) printf("Loading MIP solution failed\n");
  free(allsol);               /* Not needed any longer */
 }
 XPRBmipoptimize(prob,"c");   /* Solve the MIP-problem */

 if(XPRBgetmipstat(prob)==XPRS_MIP_SOLUTION ||
    XPRBgetmipstat(prob)==XPRS_MIP_OPTIMAL)
  printf("Best integer solution: %g\n",XPRBgetobjval(prob));
 else
  printf("Best integer solution: %g\n",solval);
}

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

int main(int argc, char **argv)
{
 XPRBprob prob;

 prob=XPRBnewprob("Coco");    /* Initialize a new problem in BCL */
 modcoco(prob);               /* Model the problem */
 solvecoco(prob);             /* Solve the problem */

 return 0;
}