/* nag_ip_bb (h02bbc) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 * Mark 6 revised, 2000.
 * Mark 7 revised, 2001.
 * Mark 8 revised, 2004.
 *
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nag_string.h>
#include <nagh02.h>

static int ex1(void);
static int ex2(void);
#ifdef __cplusplus
extern "C" {
#endif
  static void qphess(Integer n, Integer jthcol, double h[], Integer tdh,
                     double x[], double hx[], Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{
  Integer exit_status_ex1=0;
  Integer exit_status_ex2=0;

  /* Two examples are called, ex1() uses the
   * default settings to solve a problem while
   * ex2() solves another problem with some
   * of the optional parameters set by the user.
   */
  Vprintf("nag_ip_bb (h02bbc) Example Program Results.\n");
  exit_status_ex1 = ex1();
  exit_status_ex2 = ex2();

  return exit_status_ex1 == 0 && exit_status_ex2 == 0 ? 0 : 1;
}

#define A(I,J) a[(I)*tda + J]
static int ex1(void)
{
  /* Local variables */
  Nag_Boolean *intvar=0;
  Integer exit_status=0, i, is_int, j, m, n, nbnd, tda;
  NagError fail;
  double *a=0, *bl=0, *bu=0, *cvec=0, objf, *x=0;

  INIT_FAIL(fail);


  Vprintf("\nExample 1: default options used.\n");
  Vscanf(" %*[^\n]"); /* Skip headings in data file */
  Vscanf(" %*[^\n]");

  /* Read the problem dimensions */
  Vscanf(" %*[^\n]");
  Vscanf("%ld%ld", &m, &n);
  nbnd = n+m;
  if (n>=1 && m>=0)
    {
      if ( !( a = NAG_ALLOC(m*n, double)) ||
           !( cvec = NAG_ALLOC(n, double)) ||
           !( bl = NAG_ALLOC(nbnd, double)) ||
           !( bu = NAG_ALLOC(nbnd, double)) ||
           !( x = NAG_ALLOC(n, double))  ||
           !( intvar = NAG_ALLOC(n, Nag_Boolean)))
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
      tda = n;
    }
  else
    {
      Vprintf("Invalid n or m.\n");
      exit_status = 1;
      return exit_status;
    }
  /* Read objective coefficients */
  Vscanf(" %*[^\n]");
  for (i = 0; i < n; ++i)
    Vscanf("%lf", &cvec[i]);

  /* Read the matrix coefficients */
  Vscanf(" %*[^\n]");
  for (i = 0; i < m; ++i)
    for (j = 0; j < n; ++j)
      Vscanf("%lf",&A(i,j));

  /* Read the bounds */
  Vscanf(" %*[^\n]");
  for (i = 0; i < nbnd; ++i)
    Vscanf("%lf", &bl[i]);
  Vscanf(" %*[^\n]");
  for (i = 0; i < nbnd; ++i)
    Vscanf("%lf", &bu[i]);

  /* Read which variables are integer */
  Vscanf(" %*[^\n]");
  for (i = 0; i < n; ++i)
    {
      Vscanf("%ld", &is_int);
      /* is_int = 1 if integer variable, 0 if not */
      intvar[i] = is_int ? Nag_TRUE : Nag_FALSE;
    }

  /* Read the initial estimate of x */
  Vscanf(" %*[^\n]");
  for (i = 0; i < n; ++i)
    Vscanf("%lf", &x[i]);

  /* nag_ip_bb (h02bbc).
   * Solves integer programming problems using a branch and
   * bound method
   */
  nag_ip_bb(n, m, a, tda, bl, bu, intvar, cvec, (double *)0, 0,
            NULLFN, x, &objf, H02_DEFAULT, NAGCOMM_NULL, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_ip_bb (h02bbc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
 END:
  if (a) NAG_FREE(a);
  if (cvec) NAG_FREE(cvec);
  if (bl) NAG_FREE(bl);
  if (bu) NAG_FREE(bu);
  if (x) NAG_FREE(x);
  if (intvar) NAG_FREE(intvar);
  return exit_status;

} /* ex1 */
static void qphess(Integer n, Integer jthcol, double h[], Integer tdh,
                   double x[], double hx[], Nag_Comm *comm)
{
  Integer i;

  /* In this qphess function the Hessian is defined implicitly */
  for (i = 0; i < n; ++i)
    hx[i] = 2.0*x[i];
} /* qphess */


static int ex2(void)
{
  /* Local variables */
  Nag_Boolean *intvar=0, *intvar2=0;
  Integer exit_status=0, i, is_int, j, m, n, nbnd, tda;
  NagError fail;
  Nag_H02_Opt options;
  char **crnames=0, *names=0;
  double *a=0, *bl=0, *bu=0, *cvec=0, objf, red_bnd, *x=0;

  INIT_FAIL(fail);

  Vprintf("\nExample 2: some options set.\n");
  Vscanf(" %*[^\n]"); /* Skip heading */

  /* Read the problem dimensions */
  Vscanf(" %*[^\n]");
  Vscanf("%ld%ld", &m, &n);
  nbnd = n+m;
  if (n>=1 && m >=0)
    {
      if ( !( a = NAG_ALLOC(m*n, double)) ||
           !( cvec = NAG_ALLOC(n, double)) ||
           !( bl = NAG_ALLOC(nbnd, double)) ||
           !( bu = NAG_ALLOC(nbnd, double)) ||
           !( x = NAG_ALLOC(n, double))  ||
           !( intvar = NAG_ALLOC(n, Nag_Boolean)) ||
           !( intvar2 = NAG_ALLOC(n, Nag_Boolean)) ||
           !( crnames = NAG_ALLOC(nbnd, char *)) ||
           !( names = NAG_ALLOC(9*nbnd, char)) )
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
      tda = n;
    }
  else
    {
      Vprintf("Invalid n or m.\n");
      exit_status = 1;
      return exit_status;
    }
  /* Read names */
  Vscanf(" %*[^\n]");
  nbnd = n+m;
  for (i = 0; i < nbnd; ++i)
    {
      Vscanf("%s", &names[9*i]);
      crnames[i] = &names[9*i];
    }
  /* Read objective coefficients */
  Vscanf(" %*[^\n]");
  for (i = 0; i < n; ++i)
    Vscanf("%lf", &cvec[i]);

  /* Read the matrix coefficients */
  Vscanf(" %*[^\n]");
  for (i = 0; i < m; ++i)
    for (j = 0; j < n; ++j)
      Vscanf("%lf",&A(i,j));

  /* Read the bounds */
  Vscanf(" %*[^\n]");
  for (i = 0; i < nbnd; ++i)
    Vscanf("%lf", &bl[i]);
  Vscanf(" %*[^\n]");
  for (i = 0; i < nbnd; ++i)
    Vscanf("%lf", &bu[i]);

  /* Read which variables are integer */
  Vscanf(" %*[^\n]");
  for (i = 0; i < n; ++i)
    {
      Vscanf("%ld", &is_int);
      /* is_int = 1 if integer variable, 0 if not */
      intvar[i] = is_int ? Nag_TRUE : Nag_FALSE;
    }

  /* Read the initial estimate of x */
  Vscanf(" %*[^\n]");
  for (i = 0; i < n; ++i)
    Vscanf("%lf", &x[i]);

  /* nag_ip_init (h02xxc).
   * Initialize option structure to null values
   */
  nag_ip_init(&options); /* Initialise options structure */
  options.crnames = crnames;
  options.print_level = Nag_Soln;
  /* nag_ip_bb (h02bbc), see above. */
  nag_ip_bb(n, m, a, tda, bl, bu, intvar, cvec, (double *)0, 0,
            NULLFN, x, &objf, &options, NAGCOMM_NULL, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_ip_bb (h02bbc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  /* Now solve a related problem obtained by reducing lower
     bound on a constraint */

  /* Read amount to reduce lower bound on constraint 1 by */
  Vscanf(" %*[^\n]");
  Vscanf("%lf", &red_bnd);
  bl[n] -= red_bnd;

  Vprintf("\nSolve modified problem - use different tree search.\n");
  Vprintf("---------------------------------------------------\n");

  options.list = Nag_FALSE;
  if (red_bnd > 0.0)
    {
      /* We have a valid bound for the objective since this problem
         is less constrained than first one */
      options.int_obj_bound = objf + 1.0e-3;
    }
  options.nodsel = Nag_Deep_Search;
  options.print_level = Nag_Iter;

  Vprintf("***Set options.list = Nag_FALSE\n");
  Vprintf("***Set options.int_obj_bound = %15.7e\n", options.int_obj_bound);
  Vprintf("***Set options.nodsel = Nag_Deep_Search\n");
  Vprintf("***Set options.print_level = Nag_Iter\n");

  /* nag_ip_bb (h02bbc), see above. */
  nag_ip_bb(n, m, a, tda, bl, bu, intvar, cvec, (double *)0, 0,
            NULLFN, x, &objf, &options, NAGCOMM_NULL, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_ip_bb (h02bbc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  Vprintf("\n***IP objective value = %15.7e\n", objf);

  Vprintf("\n\nIllustrate effect of supplying branching directions.\n");
  Vprintf("----------------------------------------------------\n\n");

  options.branch_dir = Nag_Branch_InitX;
  Vprintf("***Set options.branch_dir = Nag_Branch_InitX\n");

  /* nag_ip_bb (h02bbc), see above. */
  nag_ip_bb(n, m, a, tda, bl, bu, intvar, cvec, (double *)0, 0,
            NULLFN, x, &objf, &options, NAGCOMM_NULL, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_ip_bb (h02bbc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  Vprintf("\n***IP objective value = %15.7e\n", objf);
  /* nag_ip_free (h02xzc).
   * Free NAG allocated memory from option structures
   */
  nag_ip_free(&options, "", &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_ip_free (h02xzc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  /* Finally, illustrate solution of an MIQP problem
     - we find the IP solution which is closest in
     least-squares sense to the root node LP solution
     of BB tree */

  Vprintf("\n\nObtain solution of root LP problem.\n");
  Vprintf("-----------------------------------\n\n");

  /* Set all variables non-integer to obtain LP solution */
  for (i = 0; i < n; ++i)
    intvar2[i] = 0;

  options.print_level = Nag_NoPrint;
  Vprintf("***Printout suppressed: options.print_level = Nag_NoPrint\n");

  /* nag_ip_bb (h02bbc), see above. */
  nag_ip_bb(n, m, a, tda, bl, bu, intvar2, cvec, (double *)0, 0,
            NULLFN, x, &objf, &options, NAGCOMM_NULL, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_ip_bb (h02bbc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  Vprintf("***LP objective value = %15.7e\n", objf);

  /* Set linear part of solution */
  for (i = 0; i < n; ++i)
    cvec[i] = -2.0*x[i];

  /* Re-initialise options structure */
  /* nag_ip_free (h02xzc), see above. */
  nag_ip_free(&options, "", &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_ip_free (h02xzc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  /* nag_ip_init (h02xxc), see above. */
  nag_ip_init(&options);
  options.crnames = crnames;

  options.list = Nag_TRUE;
  options.prob = Nag_MIQP2;

  Vprintf("\n\nFinally, solve a related MIQP problem.\n");
  Vprintf("--------------------------------------\n");

  /* nag_ip_bb (h02bbc), see above. */
  nag_ip_bb(n, m, a, tda, bl, bu, intvar, cvec, (double *)0, 0,
            qphess, x, &objf, &options, NAGCOMM_NULL, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_ip_bb (h02bbc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  /* nag_ip_free (h02xzc), see above. */
  nag_ip_free(&options, "", &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_ip_free (h02xzc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
 END:
  if (a) NAG_FREE(a);
  if (cvec) NAG_FREE(cvec);
  if (bl) NAG_FREE(bl);
  if (bu) NAG_FREE(bu);
  if (x) NAG_FREE(x);
  if (intvar) NAG_FREE(intvar);
  if (intvar2) NAG_FREE(intvar2);
  if (crnames) NAG_FREE(crnames);
  if (names) NAG_FREE(names);
  return exit_status;
} /* ex2 */