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

  file xbmt2.c
  ````````````
  Implementing multiple threads with BCL and
  Xpress-Optimizer.

  Win32 version
  
  (c) 2008-2024 Fair Isaac Corporation
      author: Y. Colombani, 2001, rev. Mar. 2011
********************************************************/

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

/* Constants for contr1() */
#define District 6                    /* Number of districts */
#define Contract 10                   /* Number of contracts */
const int OUTPUT[] = {50, 40, 10, 20, 70, 50};  /* Max. output per district */
const int COST[]   = {50, 20, 25, 30, 45, 40};  /* Cost per district */
const int VOLUME[]   = {20, 10, 30, 15, 20, 30, 10, 50, 10, 20};  
                                      /* Volume of contracts */

/* Constants for portf() */
#define NVal 30                       /* Total number of values */
#define LIMIT 20                      /* Maximum number to be chosen */

#define QFILE XPRBDATAPATH "/portf/pfqcost.dat" /* Quadratic cost coeff.s */
#define BFILE XPRBDATAPATH "/portf/pfubds.dat"  /* Upper bounds on percentages */
#define CFILE XPRBDATAPATH "/portf/pflcost.dat" /* Linear cost coefficients */

/* The problems to solve */
DWORD WINAPI contr1(LPVOID param);
DWORD WINAPI portf(LPVOID param);

/**********************/
/* Where things start */
/**********************/
int main(int argc,char *argv[])
{
 HANDLE tosolve1, tosolve2;
 XPRBprob modcontr1, modfolio;

 /* Initialize BCL and create two problems (separate initialization
    of BCL only for clarity's sake, it could as well be left to
    be done during the creation of the first problem) */
 XPRBinit();
 modcontr1=XPRBnewprob("Contr1");
 modfolio=XPRBnewprob("Folio");

 /* Start the building+solving in parallel */
 tosolve1=CreateThread(NULL, 0, contr1, modcontr1,0,0);
 tosolve2=CreateThread(NULL, 0, portf,  modfolio ,0,0);

 /* Wait for results */
 WaitForSingleObject(tosolve1, INFINITE);
 WaitForSingleObject(tosolve2, INFINITE);

 /* Clear up everything (not really required here since the program 
    terminates immediately afterwards) */
 XPRBdelprob(modcontr1);
 XPRBdelprob(modfolio);
 XPRBfinish();
 return 0;
}

/*********************/
/* Problem 'contr1'  */
/*********************/
DWORD WINAPI contr1(LPVOID param)
{
 int d,c;
 XPRBprob bprob;
 XPRSprob oprob;
 XPRBctr c1,c2,cobj;
 XPRBvar x[District][Contract];    /* Variables indicating whether a project 
                                      is chosen */
 XPRBvar y[District][Contract];    /* Quantities allocated to contractors */
 int i, ncol, len, stat, offset;
 double *sol, val;
 char *names;
 
 bprob=(XPRBprob)param;

/**** VARIABLES ****/
 for(d=0;d<District;d++)
  for(c=0;c<Contract;c++)
  {
   x[d][c] = XPRBnewvar(bprob,XPRB_BV,XPRBnewname("x_d%d",d+1),0,1);
   y[d][c] = XPRBnewvar(bprob,XPRB_SC,XPRBnewname("q_d%d",d+1),0,OUTPUT[d]);
   XPRBsetlim(y[d][c],5);
  } 

/****OBJECTIVE****/
 cobj = XPRBnewctr(bprob,"OBJ",XPRB_N);  /* Define objective: total cost */
 for(d=0;d<District;d++)
  for(c=0;c<Contract;c++)
   XPRBaddterm(cobj,y[d][c],COST[d]);     
 XPRBsetobj(bprob,cobj);                 /* Set objective function */ 
 
/**** CONSTRAINTS ****/
 for(c=0;c<Contract;c++)
 {
  c1=XPRBnewctr(bprob,"Size",XPRB_G); /* "Size": cover the required contract volume */
  c2=XPRBnewctr(bprob,"Min",XPRB_G);  /* "Min": at least 2 districts per contract */
  for(d=0;d<District;d++)
  {
   XPRBaddterm(c1,y[d][c],1);
   XPRBaddterm(c2,x[d][c],1);
  }
  XPRBaddterm(c1,NULL,VOLUME[c]);
  XPRBaddterm(c2,NULL,2); 
 } 
 
 for(d=0;d<District;d++)        /* Do not exceed max. output of any district */
 {
  c1=XPRBnewctr(bprob,"Output",XPRB_L);
  for(c=0;c<Contract;c++)
   XPRBaddterm(c1,y[d][c],1);
  XPRBaddterm(c1,NULL,OUTPUT[d]);
 } 
 
 for(d=0;d<District;d++)        /* If a contract is allocated to a district,
                                   then at least 1 unit is allocated to it */
  for(c=0;c<Contract;c++)
   XPRBnewprec(bprob,"XY",x[d][c],0,y[d][c]);

/****SOLVING + OUTPUT****/
 XPRBloadmat(bprob);            /* Load the matrix explicitely */
 oprob=XPRBgetXPRSprob(bprob);
 XPRSchgobjsense(oprob, XPRS_OBJ_MINIMIZE);  /* Set sense to minimization */
 XPRSmipoptimize(oprob, "");                 /* Solve the MIP problem */

 XPRSgetintattrib(oprob, XPRS_MIPSTATUS, &stat);  /* Get MIP status */
 if((stat==XPRS_MIP_SOLUTION) || (stat==XPRS_MIP_OPTIMAL))
                                /* Test whether an integer solution was found */
 {
  XPRSgetdblattrib(oprob, XPRS_MIPOBJVAL, &val);   
  printf("Problem %s Objective: %g\n", XPRBgetprobname(bprob), val);
  XPRSgetintattrib(oprob, XPRS_ORIGINALCOLS, &ncol); 
  sol = (double *)malloc(ncol * sizeof(double));
  XPRSgetmipsol(oprob, sol, NULL);   /* Get the primal solution values */
  XPRSgetnamelist(oprob, 2, NULL, 0, &len, 0, ncol-1);    
                        /* Get number of bytes required for retrieving names */
  names = (char *)malloc(len*sizeof(char));
  XPRSgetnamelist(oprob, 2, names, len, NULL, 0, ncol-1); 
                                     /* Get the variable names */
  offset=0;
  for(i=0; i<ncol; i++) {            /* Print out the solution */
   if(sol[i]!=0)
    printf("%s: %g, ", names+offset, sol[i]);
   offset += strlen(names+offset)+1;
  }   
  printf("\n");  
  free(sol);
  free(names);
 }
  
 return 0;
} 

/**************************/
/* The portfolio example  */
/**************************/
DWORD WINAPI portf(LPVOID param)
{
 /**** DATA ****/
 double Cost[NVal];                   /* Coeff. of lin. part of the obj. */
 double QCost[NVal][NVal];            /* Coeff. of quad. part of the obj. */
 double UBnd[NVal];                   /* Upper bound values */
 double value;
 FILE *datafile;

 int i,j;
 XPRBctr cobj, c;
 XPRBvar x[NVal];                     /* Amount of a value taken into
                                         the portfolio */
 XPRBvar y[NVal];                     /* 1 if value i is chosen, else 0 */
 XPRBprob bprob;
 
 bprob=(XPRBprob)param;

/**** Read data from files ****/
        /* Read the quadratic cost data file */
 memset(QCost,0,NVal*NVal*sizeof(double));      /* Initialize Q to 0 */
 datafile=fopen(QFILE,"r");
 while (XPRBreadlinecb(XPRB_FGETS, datafile, 200, "i,i,g", &i, &j, 
        &value) == 3)
  QCost[i-1][j-1]=value; 
 fclose(datafile);

        /* Read the linear cost data file */
 datafile=fopen(CFILE,"r");
 while (XPRBreadlinecb(XPRB_FGETS, datafile, 200, "i,g", &i, &value) == 2)
  Cost[i-1]=value; 
 fclose(datafile);

        /* Read the bounds data file */
 datafile=fopen(BFILE,"r");
 while (XPRBreadlinecb(XPRB_FGETS, datafile, 200, "i,g", &i, &value) == 2)
  UBnd[i-1]=value;  
 fclose(datafile);

/**** VARIABLES ****/
 for(i=0;i<NVal;i++) 
 {
  x[i] = XPRBnewvar(bprob,XPRB_PL, XPRBnewname("x_%d",i+1), 0, UBnd[i]);
  y[i] = XPRBnewvar(bprob,XPRB_BV, XPRBnewname("y_%d",i+1), 0, 1);
 }
 
/****OBJECTIVE****/
 cobj = XPRBnewctr(bprob,"OBJ",XPRB_N); /* Define objective: total cost */
 for(i=0;i<NVal;i++)
  XPRBaddterm(cobj,x[i],Cost[i]);

 for(i=0;i<NVal;i++)                  /* Define quadratic part of the obj. */
 {
  XPRBaddqterm(cobj,x[i],x[i],QCost[i][i]);
  for(j=i+1;j<NVal;j++)
   XPRBaddqterm(cobj,x[i],x[j], QCost[i][j]);
 }
 XPRBsetobj(bprob,cobj);              /* Set objective function */ 

/**** CONSTRAINTS ****/
 c = XPRBnewctr(bprob,"C1",XPRB_E);   /* Amounts of values chosen must
                                         add up to 100% */
 for(i=0;i<NVal;i++) XPRBaddterm(c, x[i], 1);
 XPRBaddterm(c, NULL, 100);
 
 for(i=0;i<NVal;i++)                  /* Upper limits */
 {
  c = XPRBnewctr(bprob,"UL",XPRB_L);
  XPRBaddterm(c, x[i], 1);
  XPRBaddterm(c, y[i], -UBnd[i]);
 }

 c = XPRBnewctr(bprob,"Card", XPRB_L);   /* Limit on total number of values */
 for(i=0;i<NVal;i++) XPRBaddterm(c, y[i], 1);
 XPRBaddterm(c, NULL, LIMIT); 
  
/****SOLVING + OUTPUT****/
/* XPRBprintprob(bprob); */           /* Print out the problem definition */
 XPRBexportprob(bprob,XPRB_MPS,"Portf"); /* Output matrix in MPS format */
 XPRBexportprob(bprob,XPRB_LP,"Portf");  /* Output matrix in LP format */
  
 XPRBsetsense(bprob,XPRB_MINIM);      /* Choose the sense of the optimization */
 XPRSsetintcontrol(XPRBgetXPRSprob(bprob),XPRS_CUTSTRATEGY,0);
 XPRBlpoptimize(bprob,"");            /* Solve the QP-problem, use 'mipoptimize'
                                         to solve MIQP */

 printf("Problem %s:\n", XPRBgetprobname(bprob));
 for(i=0;i<NVal;i++)
  printf("%s:%g (%g), ", XPRBgetvarname(x[i]), XPRBgetsol(x[i]), 
    XPRBgetsol(y[i]));
 printf("\n");
 return 0;
}