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

  file xbcutstk.c
  ```````````````
  Cutting stock problem, solved by column (= cutting
  pattern) generation heuristic looping over the root
  node.

  (c) 2008-2024 Fair Isaac Corporation
      author: S.Heipcke, 2001, rev. Mar. 2014
********************************************************/

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

#define NWIDTHS 5
#define MAXWIDTH 94

#define EPS 1e-6
#define MAXCOL 10

/****DATA****/
double WIDTH[] = {17, 21, 22.5,  24, 29,5};  /* Possible widths */
int DEMAND[] =  {150, 96,   48, 108,  227};  /* Demand per width */
int PATTERNS[NWIDTHS][NWIDTHS];              /* (Basic) cutting patterns */

XPRBvar use[NWIDTHS+MAXCOL];               /* Rolls per pattern */
XPRBctr dem[NWIDTHS];                      /* Demand constraints */
XPRBctr cobj;                              /* Objective function */

double knapsack(int N, double *c, double *a, double R, int *d, int *xbest);

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

void modcutstock(XPRBprob prob)
{
 int i,j;

 for(j=0;j<NWIDTHS;j++) PATTERNS[j][j]=(int)floor(MAXWIDTH/WIDTH[j]); for(j=0;j<NWIDTHS;j++) use[j]=XPRBnewvar(prob,XPRB_UI, cobj=XPRBnewctr(prob,"OBJ",XPRB_N); Minimize total number of rolls for(i=0;i<NWIDTHS;i++) dem[i]=XPRBnewctr(prob,"Demand",XPRB_G); Satisfy the demand per width Column generation loop at top node: solve LP and save basis get solution values generate new (=cutting load modified problem saved void double Objective value int Counters for columns passes Solution variables pat Dual constraints XPRBbasis starttime=XPRBgettime(); npatt=NWIDTHS; for(npass=0;npass<MAXCOL;npass++) Solve Save current objval=XPRBgetobjval(prob); Get objective values: soluse[j]=XPRBgetsol(use[j]); dualdem[i]=XPRBgetdual(dem[i]); printf(%g,=",soluse[j]); printf(" integer knapsack z=min{cx : ax<=r, x in with r=MAXWIDTH, n=NWIDTHS Pass %d:=",(XPRBgettime()-starttime)/1000.0, npass+1);

  if(z < 1+EPS)
  {
   printf(" no profitable column No need to keep any longer else Print pattern: pattern found marginal cost %g=", z-1);
   printf(" Widths distribution:=");
   dw=0;
   for(i=0;i<NWIDTHS;i++)
   {
    printf (" %g:%d=", WIDTH[i], x[i]);
    dw += WIDTH[i]*x[i];
   }
   printf(" Total width: Create a variable this use[npatt]=XPRBnewvar(prob,XPRB_UI, Add var. dw=0;> EPS)
    {
     XPRBaddterm(dem[i], use[npatt], x[i]);
     if((int)ceil((double)DEMAND[i]/x[i]) > dw)
      dw = (int)ceil((double)DEMAND[i]/x[i]);
    }
   XPRBsetub(use[npatt],dw);         /* Change the upper bound on the new var.*/

   npatt++;

   XPRBloadmat(prob);                /* Reload the problem */
   XPRBloadbasis(basis);             /* Load the saved basis */
   XPRBdelbasis(basis);              /* No need to keep the basis any longer */
  }
 }

 XPRBmipoptimize(prob,"c");          /* Solve the MIP */

 printf("(%g sec) Optimal solution: %g rolls, %d patterns\n   ",
   (XPRBgettime()-starttime)/1000.0, XPRBgetobjval(prob), npatt);
 printf("Rolls per pattern: ");
 for(i=0;i<npatt;i++) printf("%g, ", XPRBgetsol(use[i]));
 printf("\n");
}

/**************************************************************************/
/* Integer Knapsack Algorithm for solving the integer knapsack problem    */
/*    z = max{cx : ax <= R, x <= d, x in Z^N}                             */
/*                                                                        */
/* Input data:                                                            */
/*   N:        Number of item types                                       */
/*   c[i]:     Unit profit of item type i, i=1..n                         */
/*   a[i]:     Unit resource use of item type i, i=1..n                   */
/*   R:        Total resource available                                   */
/*   d[i]:     Demand for item type i, i=1..n                             */
/* Return values:                                                         */
/*   xbest[i]: Number of items of type i used in optimal solution         */
/*   zbest:    Value of optimal solution                                  */
/**************************************************************************/
double knapsack(int N, double *c, double *a, double R, int *d, int *xbest)
{
 int j;
 double zbest = 0.0;
 XPRBarrvar x;
 XPRBctr cobj;
 XPRBprob kprob;
/*
 printf ("Solving z = max{cx : ax <= b; x in Z^n}\n   c   =");
 for (j = 0; j < N; j++) printf (" %g,", c[j]);
 printf("\n   a   =");
 for (j = 0; j < N; j++) printf (" %g,", a[j]);
 printf("\n   c/a =");
 for (j = 0; j < N; j++) printf (" %g,", c[j]/a[j]);
 printf("\n   b   = %f\n", R);
*/

 kprob=XPRBnewprob("Knapsack");
 x=XPRBnewarrvar(kprob,N, XPRB_UI, "x", 0, XPRB_INFINITY);
 for(j=0;j<N;j++) XPRBsetub(x[j], d[j]);

 cobj = XPRBnewarrsum(kprob,"OBJ", x, c, XPRB_N, 0);
 XPRBsetobj(kprob,cobj);
 XPRBnewarrsum(kprob,"KnapCtr", x, a, XPRB_L, R);
 XPRBsetsense(kprob,XPRB_MAXIM);
 XPRBmipoptimize(kprob,"");

 zbest = XPRBgetobjval(kprob);
 for(j=0;j<N;j++) xbest[j]=(int)floor(XPRBgetsol(x[j])+0,5);
/*
 printf ("   z = %f\n   x =", zbest);
 for(j=0; j<N; j++) printf (" %d,", xbest[j]);
 printf("\n");
*/

 XPRBdelprob(kprob);

 return (zbest);
}

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

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

 prob=XPRBnewprob("CutStock");         /* Initialize a new problem in BCL */
 modcutstock(prob);                    /* Model the problem */
 solvecutstock(prob);                  /* Solve the problem */

 XPRBdelprob(prob);                    /* Delete the main problem */
 XPRBfree();                           /* Terminate BCL */

 return 0;
}