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

  file xbpurch.c
  ``````````````
  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 <stdio.h>
#include <stdlib.h>
#include "xprb.h"

#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 prob)
{
 XPRBctr *refrow, cobj,ctr;
 int s,b;
 XPRBvar *x,**lam;
 XPRBsos sos2;

/****VARIABLES****/
 x=(XPRBvar *)malloc(NS*sizeof(XPRBvar));
 for(s=0;s<NS;s++)
  x[s]=XPRBnewvar(prob,XPRB_PL, "x", 0, XPRB_INFINITY);
                         /* Quantity to purchase from supplier s */

 lam=(XPRBvar **)malloc(NS*sizeof(XPRBvar *));
 for(s=0;s<NS;s++)
 {
  lam[s]=(XPRBvar *)malloc(NB2*sizeof(XPRBvar));
  for(b=0;b<NB2;b++)
   lam[s][b]=XPRBnewvar(prob,XPRB_PL, XPRBnewname("lam_%d",s+1), 0, XPRB_INFINITY);
 }                       /* Weights at breakpoint b for supplier s */

/****OBJECTIVE****/
    /* Objective: minimize the sum of costs*weights */
 cobj = XPRBnewctr(prob,"OBJ",XPRB_N);
 for(s=0;s<NS;s++)
  for(b=0;b<NB2;b++)
   XPRBaddterm(cobj, lam[s][b], C[s][b]);
 XPRBsetobj(prob,cobj);       /* Select objective function */

/****CONSTRAINTS****/
    /* Define x and also order the lam variables by breakpoint quantities */
 refrow=(XPRBctr *)malloc(NS*sizeof(XPRBctr));
 for(s=0;s<NS;s++)
 {
  refrow[s]=XPRBnewctr(prob,"DefX", XPRB_E);
  for(b=0;b<NB2;b++)
   XPRBaddterm(refrow[s], lam[s][b], X[s][b]);
  XPRBaddterm(refrow[s], x[s], -1);
 }

    /* The convexity row (lam sum to 1) */
 for(s=0;s<NS;s++)
 {
  ctr=XPRBnewctr(prob,"Conv", XPRB_E);
  for(b=0;b<NB2;b++)
   XPRBaddterm(ctr, lam[s][b], 1);
  XPRBaddterm(ctr, NULL, 1);
 }

    /* The minimum quantity that must be bought */
 ctr=XPRBnewctr(prob,"MustBuy", XPRB_G);
 for(s=0;s<NS;s++)
  XPRBaddterm(ctr, x[s], 1);
 XPRBaddterm(ctr, NULL, REQ);

/****BOUNDS****/
    /* No more than the maximum percentage from each supplier */
 for(s=0;s<NS;s++) XPRBsetub(x[s], 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. */
 for(s=0;s<NS;s++)
 {
  sos2=XPRBnewsos(prob,"SetLam",XPRB_S2);
  for(b=0;b<NB2;b++)
   XPRBaddsosel(sos2,lam[s][b],X[s][b]+1); /* Add 1 to all SOS coefficients
                                              because BCL does not accept 0 */
 }

/****SOLVING + OUTPUT****/
 XPRBexportprob(prob,XPRB_MPS,"purc");  /* Write out an MPS file */

 XPRBsetsense(prob,XPRB_MINIM);         /* Choose the sense of optimization */
 XPRBmipoptimize(prob,"");              /* Solve the MIP-problem */

 printf("Objective: %g\n", XPRBgetobjval(prob));  /* Get objective value */

 for(s=0;s<NS;s++)                      /* Print out the solution values */
  printf("%s:%g ", XPRBgetvarname(x[s]), XPRBgetsol(x[s]));
 for(s=0;s<NS;s++)
  for(b=0;b<NB2;b++)
   printf("%s:%g ", XPRBgetvarname(lam[s][b]), XPRBgetsol(lam[s][b]));
 printf("\n");
}

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

    /**** Read data from files ****/
void readdata(void)
{
 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=(double **)malloc(NS*sizeof(double *));
 for(s=0;s<NS;s++) UC[s]=(double *)malloc(NB*sizeof(double));

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

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

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

 DELTA=(double *)malloc(NB2*sizeof(double));
 MAXPERC=(double *)malloc(NS*sizeof(double));

        /* 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 prob;

 prob=XPRBnewprob("Purchase");     /* Initialize a new problem in BCL */
 readdata();                       /* Data input from file */
 modpurchase(prob);                /* Formulate and solve the problem */

 return 0;
}