/* nag_ip_bb (h02bbc) Example Program.
 *
 * Copyright 2014 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 <string.h>
#include <nag_stdlib.h>
#include <nag_string.h>
#include <nagh02.h>

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

#define A(I, J) a[(I) *tda + J]

int main(void)
{
  static double ruser[1] = {-1.0};
  Integer     exit_status = 0;
  Integer     i, j, m, n, nbnd, tda;
  char        **crnames = 0, *names = 0;
  double      *a = 0, *bl = 0, *bu = 0, *cvec = 0, objf, red_bnd, *x = 0;
  Nag_Boolean *intvar = 0, *intvar2 = 0;
  char        nag_enum_arg[40];
  Nag_Comm    comm;
  Nag_H02_Opt options;
  NagError    fail;

  INIT_FAIL(fail);

  printf("nag_ip_bb (h02bbc) Example Program Results\n");

  /* For communication with user-supplied functions: */
  comm.user = ruser;
  scanf(" %*[^\n]"); /* Skip heading */

  /* Read the problem dimensions */
  scanf(" %*[^\n]");
  scanf("%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(nbnd*9, char)))
        {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
      tda = n;
    }
  else
    {
      printf("Invalid n or m.\n");
      exit_status = 1;
      return exit_status;
    }
  /* Read names */
  scanf(" %*[^\n]");
  nbnd = n+m;
  for (i = 0; i < nbnd; ++i)
    {
      scanf("%8s", &names[9*i]);
      crnames[i] = &names[9*i];
    }
  /* Read objective coefficients */
  scanf(" %*[^\n]");
  for (i = 0; i < n; ++i)
    scanf("%lf", &cvec[i]);

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

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

  /* Read which variables are integer */
  scanf(" %*[^\n]");
  for (i = 0; i < n; ++i)
    {
      scanf("%39s", nag_enum_arg);
      /* intvar = Nag_TRUE if integer variable, Nag_FALSE if not */
      intvar[i] = (Nag_Boolean) nag_enum_name_to_value(nag_enum_arg);
    }

  /* Read the initial estimate of x */
  scanf(" %*[^\n]");
  for (i = 0; i < n; ++i)
    scanf("%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. */
  fflush(stdout);
  nag_ip_bb(n, m, a, tda, bl, bu, intvar, cvec, (double *) 0, 0,
            NULLFN, x, &objf, &options, &comm, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("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 */
  scanf(" %*[^\n]");
  scanf("%lf", &red_bnd);
  bl[n] -= red_bnd;

  printf("\nSolve modified problem - use different tree search.\n");
  printf("---------------------------------------------------\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;

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

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

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

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

  /* nag_ip_bb (h02bbc), see above. */
  fflush(stdout);
  nag_ip_bb(n, m, a, tda, bl, bu, intvar, cvec, (double *) 0, 0,
            NULLFN, x, &objf, &options, &comm, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_ip_bb (h02bbc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  printf("\n***IP objective value = %16.7e\n", objf);
  /* nag_ip_free (h02xzc).
   * Free NAG allocated memory from option structures
   */
  nag_ip_free(&options, "", &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("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 */

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

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

  options.print_level = Nag_NoPrint;
  printf("***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, &comm, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_ip_bb (h02bbc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  printf("***LP objective value = %16.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. */
  fflush(stdout);
  nag_ip_free(&options, "", &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("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;

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

  /* nag_ip_bb (h02bbc), see above. */
  fflush(stdout);
  nag_ip_bb(n, m, a, tda, bl, bu, intvar, cvec, (double *) 0, 0,
            qphess, x, &objf, &options, &comm, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("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)
    {
      printf("Error from nag_ip_free (h02xzc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

 END:
  NAG_FREE(a);
  NAG_FREE(cvec);
  NAG_FREE(bl);
  NAG_FREE(bu);
  NAG_FREE(x);
  NAG_FREE(intvar);
  NAG_FREE(intvar2);
  NAG_FREE(crnames);
  NAG_FREE(names);

  return exit_status;
}

static void NAG_CALL qphess(Integer n, Integer jthcol, const double h[],
                            Integer tdh, const double x[], double hx[],
                            Nag_Comm *comm)
{
  Integer i;
  if (comm->user[0] == -1.0)
    {
      printf("(User-supplied callback qphess, first invocation.)\n");
      fflush(stdout);
      comm->user[0] = 0.0;
    }
  /* In this qphess function the Hessian is defined implicitly */
  if (jthcol==0) {
    for (i = 0; i < n; ++i)
      hx[i] = 2.0*x[i];
  }
  else {
    for (i = 0; i < n; ++i)
      hx[i] = (i == jthcol - 1 ? 2.0 : 0.0);
  }
} /* qphess */