/* nag_mip_ilp_dense (h02bbc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.3, 2024.
 *
 */

#include <nag.h>
#include <stdio.h>
#include <string.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_mip_ilp_dense (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("%" NAG_IFMT "%" NAG_IFMT "", &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); /* Initialize options structure */
  options.crnames = crnames;
  options.list = Nag_FALSE;
  options.print_level = Nag_NoPrint;
  /* nag_mip_ilp_dense (h02bbc), see above. */
  fflush(stdout);
  nag_mip_ilp_dense(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_mip_ilp_dense (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.list = Nag_FALSE;
  options.print_level = Nag_NoPrint;

  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_NoPrint\n");

  /* nag_mip_ilp_dense (h02bbc), see above. */
  fflush(stdout);
  nag_mip_ilp_dense(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_mip_ilp_dense (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_mip_ilp_dense (h02bbc), see above. */
  fflush(stdout);
  nag_mip_ilp_dense(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_mip_ilp_dense (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.list = Nag_FALSE;
  options.print_level = Nag_NoPrint;

  /* nag_mip_ilp_dense (h02bbc), see above. */
  nag_mip_ilp_dense(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_mip_ilp_dense (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-initialize 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_FALSE;
  options.print_level = Nag_NoPrint;
  options.prob = Nag_MIQP2;

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

  /* nag_mip_ilp_dense (h02bbc), see above. */
  fflush(stdout);
  nag_mip_ilp_dense(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_mip_ilp_dense (h02bbc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  printf("***MIQP2 objective value = %16.7e\n", objf);

  /* 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 */