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

  file xbpurch.cxx
  ````````````````
  Purchasing problem using SOS-2.

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

/* WARNING. This model is not for the novice, but it does contain many  *
 * useful ideas.                                                        *
 * The formulation is tricky because the discounts are all-quantity,    *
 * so the graph of cost against quantity purchased is discontinuous.    *
 * To maintain sanity in the special ordered set formulation, we must   *
 * coarsen the discontinuity by stopping purchases just at the break    *
 * point. */
 
#include <iostream>
#include <cstdio>
#include "xprb_cpp.h"

using namespace std;
using namespace ::dashoptimization;

#define PARAMSFILE  XPRBDATAPATH "/purchase/params.dat"
#define MAXPERCFILE XPRBDATAPATH "/purchase/maxperc.dat"
#define REQFILE     XPRBDATAPATH "/purchase/required.dat"
#define PRICEFILE   XPRBDATAPATH "/purchase/pricebk.dat"

/****TABLES****/
int NS;             /* Number of suppliers */
int NB;             /* Number of price breaks */
int NB2;            /* Useful parameter */

double **UC;        /* Unit cost */
double **BR;        /* Breakpoints (quantities at which unit cost changes) */
double **X;         /* Coarsened break points */
double **C;         /* Total cost at break points */
double *DELTA;      /* Coarsening factors */
double *MAXPERC;    /* Maximum percentage from each supplier */
double REQ;         /* Total quantity required */

double delta = 0.10;        /* Base coarsening factor */
                    
/***********************************************************************/

void modPurchase(XPRBprob &p)
{
 XPRBexpr lobj, *lr, lc;
 int s,b;
 XPRBvar *x,**lam;
 
/****VARIABLES****/
 x = new XPRBvar[NS];       /* Quantity to purchase from supplier s */
 for(s=0;s<NS;s++)  x[s] = p.newVar("x");

 lam = new XPRBvar*[NS];
 for(s=0;s<NS;s++)
 {
  lam[s] = new XPRBvar[NB2];
  for(b=0;b<NB2;b++)
   lam[s][b] = p.newVar(XPRBnewname("lam_%d",s+1));
 }                          /* Weights at breakpoint b for supplier s */

/****OBJECTIVE****/
 for(s=0;s<NS;s++)          /* Minimize the sum of costs*weights */
  for(b=0;b<NB2;b++)  lobj += C[s][b]*lam[s][b];
 p.setObj(p.newCtr("OBJ",lobj));   /* Set objective function */ 

/****CONSTRAINTS****/  
    /* Define x and also order the lam variables by breakpoint quantities */
 lr = new XPRBexpr[NS];
 for(s=0;s<NS;s++)
 {
  for(b=0;b<NB2;b++)  lr[s] += X[s][b]*lam[s][b];
  p.newCtr("DefX", lr[s] == x[s]);
 }

    /* The convexity constraint (lam sum to 1) */
 for(s=0;s<NS;s++)
 {
  lc=0;
  for(b=0;b<NB2;b++)  lc += lam[s][b];
  p.newCtr("Conv", lc == 1);
 }

    /* The minimum quantity that must be bought */
 lc=0;
 for(s=0;s<NS;s++) lc += x[s];
 p.newCtr("MustBuy", lc >= REQ);

/****BOUNDS****/
    /* No more than the maximum percentage from each supplier */
 for(s=0;s<NS;s++)  x[s].setUB(MAXPERC[s]*REQ/100.0);

/****SETS****/
   /* Define the lam as SOS2 as we can linearly interpolate between the
    * breakpoints. The weights are the DefX coefficients augmented by 1
      because BCL does not accept the weight 0. */
 for(s=0;s<NS;s++) 
 { 
  for(b=0;b<NB2;b++)  lr[s] += lam[s][b];
  p.newSos("SetLam", XPRB_S2, lr[s]);
 } 
   
/****SOLVING + OUTPUT****/
 p.exportProb(XPRB_MPS,"purc");    /* Output to MPS file */
 p.setSense(XPRB_MINIM);           /* Choose the sense of the optimization */
 p.mipOptimize("");                /* Solve the MIP-problem */

 cout << "Objective: " << p.getObjVal() << endl;  /* Get objective value */
    
 for(s=0;s<NS;s++)                 /* Print out the solution values */
  cout << x[s].getName() << ":" << x[s].getSol() << " ";  
 cout << endl;
 for(s=0;s<NS;s++)      
  for(b=0;b<NB2;b++)
   cout << lam[s][b].getName() << ":" << lam[s][b].getSol() << " ";
 cout << endl;
 
 delete [] BR;
 delete [] UC;
 delete [] lr;
 for(s=0;s<NS;s++) delete [] lam[s];
 delete [] lam;
 delete [] x;
}

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

    /**** Read data from files ****/
void readData()
{
 FILE *datafile;
 int s,b;  
 double val1,val2;
 
        /* Read the parameter file */
 datafile=fopen(PARAMSFILE,"r");
 XPRBreadlinecb(XPRB_FGETS, datafile, 200, "i", &NS);
 XPRBreadlinecb(XPRB_FGETS, datafile, 200, "i", &NB);
 fclose(datafile);
 NB2 = 2*NB;

        /* Now define the tables that we will use, using the sizes we 
           have just read. */
 UC= new double*[NS];
 for(s=0;s<NS;s++) UC[s]=new double[NB];

 BR = new double*[NS];
 for(s=0;s<NS;s++) BR[s]=new double[NB];

 X=new double*[NS];
 for(s=0;s<NS;s++) X[s]=new double[NB2];

 C=new double*[NS];
 for(s=0;s<NS;s++) C[s]=new double[NB2];

 DELTA=new double[NB2];
 MAXPERC=new double[NS];
 
        /* Define coarsening factors. */
 DELTA[0]=0;
 DELTA[1] = -1*delta;
 for(b=2;b<NB2;b++) DELTA[b] = -1*DELTA[b-1];

        /* Read the price and breakpoint data file */
 for(s=0;s<NS;s++)
  for(b=0;b<NB;b++)
  {
   UC[s][b]=0;
   BR[s][b]=0;
  }
 datafile=fopen(PRICEFILE,"r");
 while(XPRBreadlinecb(XPRB_FGETS, datafile, 200, "i,i,g,g", &s, &b, 
       &val1, &val2)==4)
 {
  UC[s-1][b-1]=val1;
  BR[s-1][b-1]=val2;
 }
 fclose(datafile);

        /* Read the percentages data file */
 datafile=fopen(MAXPERCFILE,"r");
 XPRBreadarrlinecb(XPRB_FGETS, datafile, 200, "g,", MAXPERC, NS);
 fclose(datafile); 
 
        /* Read the required amount data file */
 datafile=fopen(REQFILE,"r");
 XPRBreadlinecb(XPRB_FGETS, datafile, 200, "g", &REQ);
 fclose(datafile); 

        /* We now pick out the data from the tables into which they have 
           been read. */
 for(s=0;s<NS;s++) 
 {
  C[s][0] = 0;      /* First total cost and (new) breakpoint are (0,0) */
  X[s][0] = 0;              
  for(b=1;b<NB2;b++)
  {                          /* Rest calculated... */
   C[s][b] = UC[s][b/2] * BR[s][(b-1)/2];   /* ...unit price*quantity */
   X[s][b] = BR[s][(b-1)/2] + DELTA[b];     /* ...coarsened grids */
  }
 }
}

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

int main(int argc, char **argv)
{
 XPRBprob p("Purchase"); /* Initialize a new problem in BCL */
 readData();             /* Data input from file */
 modPurchase(p);         /* Formulate and solve the problem */
 
 return 0;
}