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

  file d1wagon2.c
  ``````````````
   Load balancing of train wagons
   (second version, using heuristic solution as
    start solution for MIP)

  (c) 2014-2024 Fair Isaac Corporation
      author: L. Bertacco, September 2014
********************************************************/

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

#define NBOXES (sizeof(WEIGHT)/sizeof(*WEIGHT)) /* Number of boxes                       */
#define NWAGONS 3                               /* Number of wagons                      */

                                                /* Box weights                           */
int WEIGHT[] = { 34, 6, 8, 17, 16, 5, 13, 21, 25, 31, 14, 13, 33, 9, 25, 25 };
int WMAX = 100;                                 /* Weight limit of the wagons            */

int HeurSol[NBOXES];                            /* Heuristic solution: for each box      */
                                                /* specifies in which wagon it is loaded */

/****VARIABLES****/
XPRBvar load[NBOXES][NWAGONS];                  /* 1 if box loaded on wagon, 0 otherwise */
XPRBvar maxweight;                              /* Weight of the heaviest wagon load     */

XPRBprob prob;

void XPRS_CC solnotify(XPRSprob my_prob, void* my_object, const char* solname, int status);

void d1w2_model(XPRBprob prob)
{
  int b, w;
  double sum_weights = 0;
  XPRBctr ctr;

  /****VARIABLES****/

  /* Create load[box,wagon] (binary) */
  for (b = 0; b < NBOXES; b++) for(w = 0; w < NWAGONS; w++)
    load[b][w] = XPRBnewvar(prob, XPRB_BV, XPRBnewname("load_%d_%d",b+1,w+1), 0, 1);

  /* Create maxweight (continuous with lb=ceil((sum(b in BOXES) WEIGHT(b))/NBOXES) */
  for (b = 0; b < NBOXES; b++) sum_weights += WEIGHT[b];
  maxweight = XPRBnewvar(prob, XPRB_PL, "maxweight", ceil(sum_weights/NBOXES), XPRB_INFINITY);

  /****CONSTRAINTS****/

  /* Every box into one wagon: forall(b in BOXES) sum(w in WAGONS) load(b,w) = 1 */
  for (b = 0; b < NBOXES; b++) {
    ctr = XPRBnewctr(prob, NULL, XPRB_E);
    for (w = 0; w < NWAGONS; w++) XPRBaddterm(ctr, load[b][w], 1);
    XPRBaddterm(ctr, NULL, 1);
  }

  /* Limit the weight loaded into every wagon: forall(w in WAGONS) sum(b in BOXES) WEIGHT(b)*load(b,w) <= maxweight */
  for (w = 0; w < NWAGONS; w++) {
    ctr = XPRBnewctr(prob, NULL, XPRB_L);
    for (b = 0; b < NBOXES; b++) XPRBaddterm(ctr, load[b][w], WEIGHT[b]);
    XPRBaddterm(ctr, maxweight, -1);
  }

  /****OBJECTIVE****/

  ctr = XPRBnewctr(prob, "MaxWeight", XPRB_N);
  XPRBaddterm(ctr, maxweight, 1);
  XPRBsetobj(prob, ctr);
  XPRBsetsense(prob, XPRB_MINIM);
}

void d1w2_solve(XPRBprob prob)
{
  int b, w;
  int statmip;
  XPRBsol sol;

  /* Alternative to lower bound on maxweight: adapt the optimizer cutoff value  */
  /* XPRSsetdblcontrol(XPRBgetXPRSprob(prob), XPRS_MIPADDCUTOFF, -0.99999); */

  /* Comment out the following line to enable the optimizer log */
  XPRSsetintcontrol(XPRBgetXPRSprob(prob), XPRS_OUTPUTLOG, 0);

  /* Create a BCL solution from the heuristic solution we have found */
  sol = XPRBnewsol(prob);
  /* Set the solution values for all discrete variables that are non-zero */
  for (b = 0; b < NBOXES; b++) XPRBsetsolvar(sol, load[b][HeurSol[b]], 1);

  /* It is possible, but not necessary, to set values for ALL discrete vars  */
  /* by uncommenting the following line. In this case, the usersolnotify     */
  /* callback would return status equal to 2 (instead of 3), as the solution */
  /* would be feasible without the need of a local search.                   */
  /* for (b=0; b<NBOXES; b++) for (w=0; w<NWAGONS; w++) XPRBsetsolvar(sol, load[b][w], w==HeurSol[b]); */

  XPRBaddmipsol(prob, sol, "heurSol"); /* Send the solution to the optimizer */
  XPRBdelsol(sol);                     /* Free the solution object */
  /* Request notification of solution status after processing */
  XPRSaddcbusersolnotify(XPRBgetXPRSprob(prob), solnotify, NULL, 0);

  /* Parameter settings to make use of loaded solution */
  XPRSsetdblcontrol(XPRBgetXPRSprob(prob), XPRS_HEURSEARCHEFFORT, 2);
  XPRSsetintcontrol(XPRBgetXPRSprob(prob), XPRS_HEURSEARCHROOTSELECT, 31);
  XPRSsetintcontrol(XPRBgetXPRSprob(prob), XPRS_HEURSEARCHTREESELECT, 19);

  XPRBmipoptimize(prob,"c");       /* Solve the MIP problem */
  statmip = XPRBgetmipstat(prob); /* Get the problem status */
  if (statmip == XPRB_MIP_SOLUTION || statmip == XPRB_MIP_OPTIMAL) { /* An integer solution has been found */
    printf("Optimal solution:\n Max weight: %g\n", XPRBgetobjval(prob));
    for (w = 0; w < NWAGONS; w++) {
      int tot_weight = 0;
      printf(" %d:", w + 1);
      for (b = 0; b < NBOXES; b++) if (XPRBgetsol(load[b][w]) > .5) {
        printf(" %d", b+1);
        tot_weight += WEIGHT[b];
      }
      printf(" (total weight: %d)\n", tot_weight);
    }
  }
}

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

/* LPT (Longest processing time) heuristic:     */
/* One at a time, place the heaviest unassigned */
/* box onto the wagon with the least load       */
int weight_cmp(const int *arg1, const int *arg2) { return WEIGHT[*arg2] - WEIGHT[*arg1]; }
double solve_heur()
{
  int b, w, i;
  int ORDERW[NBOXES];                   /* Box indices sorted in decreasing weight order                                              */
  int CurNum[NWAGONS] = { 0 };          /* For each wagon w, this is the number of boxes currently loaded                             */
  int CurWeight[NWAGONS] = { 0 };       /* For each wagon w, this is the current weight, i.e. the sum of weights of loaded boxes      */
  int Load[NWAGONS][NBOXES] = { 0 };    /* Load[w][i] (for i=0..CurNum[w]-1) contains the box index of the i-th box loaded on wagon w */
  double heurobj = 0;                   /* heuristic solution objective value (max wagon weight) */

  /* Copy the box indices into array ORDERW and sort them in decreasing     */
  /* order of box weights (the sorted indices are returned in array ORDERW) */
  for (b = 0; b < NBOXES; b++) ORDERW[b] = b;
  qsort(ORDERW, NBOXES, sizeof(*ORDERW), weight_cmp);

  /* Distribute the loads to the wagons using the LPT heuristic  */
  for (b = 0; b < NBOXES; b++) {
    int v = 0;                          /* Find wagon v with the smallest load */
    for (w = 0; w < NWAGONS; w++) if (CurWeight[w] <= CurWeight[v]) v = w;
    Load[v][CurNum[v]] = ORDERW[b];     /* Add current box to wagon v */
    CurNum[v]++;                        /* Increase the counter of boxes on v */
    CurWeight[v] += WEIGHT[ORDERW[b]];  /* Update current weight of the wagon */
  }

  /* Calculate the solution value */
  for (w = 0; w < NWAGONS; w++) if (CurWeight[w]>heurobj) heurobj = CurWeight[w];

  /* Solution printing */
  printf("Heuristic solution:\n Max weight: %g\n", heurobj);
  for (w = 0; w < NWAGONS; w++) {
    printf(" %d:", w + 1);
    for (i = 0; i < CurNum[w]; i++) printf(" %d", Load[w][i]+1);
    printf(" (total weight: %d)\n", CurWeight[w]);
  }

  /* Save the heuristic solution into the HeurSol array */
  for (w = 0; w < NWAGONS; w++) for (i = 0; i < CurNum[w]; i++) HeurSol[Load[w][i]] = w;

  return heurobj;
}

/* Callback function reporting loaded solution status */
void XPRS_CC solnotify(XPRSprob my_prob, void* my_object, const char* solname, int status)
{
  printf("Optimizer loaded solution %s with status=%d\n", solname, status);
}

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

int main(int argc, char **argv)
{
  XPRBprob prob = XPRBnewprob("d1wagon2"); /* Initialize a new problem in BCL */

  if (solve_heur() <= WMAX) {
    printf("Heuristic solution fits capacity limits\n");
    exit(0);
  }

  d1w2_model(prob);             /* Model the problem */
  d1w2_solve(prob);             /* Solve the problem */

  return 0;
}