/* nag_mip_sqp (h02dac) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 */
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagh02.h>

#ifdef __cplusplus
extern "C"
{
#endif
  static void NAG_CALL confun(Integer *mode, Integer ncnln, Integer n,
                              const Integer varcon[], const double x[],
                              double c[], double cjac[], Integer nstate,
                              Nag_Comm *comm);
  static void NAG_CALL objfun(Integer *mode, Integer n,
                              const Integer varcon[], const double x[],
                              double *objmip, double objgrd[], Integer nstate,
                              Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

#define CJAC(I, J) cjac[(J-1)*ncnln+I-1]
#define A(I, J)    a[(J-1)*pda+I-1]

int main(void)
{
  /* Integer scalar and array declarations */
  const Integer liopts = 200, lopts = 100, lcvalue = 40;
  Integer i, j, pda, maxit, n, nclin, ncnln, exit_status = 0;
  Integer iopts[200], p, *varcon = 0, ivalue;

  /* NAG structures and types */
  Nag_VariableType optype;
  NagError fail;
  Nag_Comm comm;

  /* Double scalar and array declarations */
  double acc, accqp, objmip;
  double *a = 0, *ax = 0, *bl = 0, *bu = 0, *c = 0, *cjac = 0;
  double *d = 0, *objgrd = 0, *x = 0, opts[200], rho;
  static double ruser[2] = { -1.0, -1.0 };

  /* Character declarations */
  char cvalue[40];

  /* Initialize the error structure */
  INIT_FAIL(fail);

  printf("nag_mip_sqp (h02dac) Example Program Results\n\n");

  n = 8;
  nclin = 5;
  ncnln = 2;

  pda = nclin;

  if (!(a = NAG_ALLOC(n * pda, double)) ||
      !(d = NAG_ALLOC(nclin, double)) ||
      !(ax = NAG_ALLOC(nclin, double)) ||
      !(bl = NAG_ALLOC(n, double)) ||
      !(bu = NAG_ALLOC(n, double)) ||
      !(varcon = NAG_ALLOC(n + nclin + ncnln, Integer)) ||
      !(x = NAG_ALLOC(n, double)) ||
      !(c = NAG_ALLOC(ncnln, double)) ||
      !(cjac = NAG_ALLOC(ncnln * n, double)) ||
      !(objgrd = NAG_ALLOC(n, double)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  for (i = 0; i < 4; i++) {
    /* Set variable types: continuous then binary */
    varcon[i] = 0;
    varcon[4 + i] = 1;

    /* Set continuous variable bounds */
    bl[i] = 0.0;
    bu[i] = 1.0e3;
  }

  /* Bounds of binary variables need not be provided */
  for (i = 4; i < 8; i++) {
    bl[i] = 0.0;
    bu[i] = 1.0;
  }

  /* Set linear constraint, equality first */
  varcon[n] = 3;
  varcon[n + 1] = varcon[n + 2] = varcon[n + 3] = varcon[n + 4] = 4;

  /* Set Ax=d then Ax>=d */
  for (i = 1; i <= nclin; i++) {
    for (j = 1; j <= n; j++) {
      A(i, j) = 0.0;
    }
  }
  A(1, 1) = A(1, 2) = A(1, 3) = A(1, 4) = 1.0;
  A(2, 1) = -1.0;
  A(2, 5) = 1.0;
  A(3, 2) = -1.0;
  A(3, 6) = 1.0;
  A(4, 3) = -1.0;
  A(4, 7) = 1.0;
  A(5, 4) = -1.0;
  A(5, 8) = 1.0;
  d[0] = 1.0;
  d[1] = d[2] = d[3] = d[4] = 0.0;

  /* Set constraints supplied by CONFUN, equality first */
  varcon[n + nclin] = 3;
  varcon[n + nclin + 1] = 4;

  /* Initialize communication arrays */
  nag_mip_opt_set("Initialize = nag_mip_sqp", iopts, liopts, opts, lopts,
                  &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_mip_opt_set (h02zkc).\n%s\n", fail.message);
    exit_status = -1;
    goto END;
  }

  /* Optimisiation parameters */
  maxit = 500;
  acc = 1.0e-6;

  /* Initial estimate (binary variables need not be given) */
  x[0] = x[1] = x[2] = x[3] = 1.0;
  x[4] = x[5] = x[6] = x[7] = 0.0;

  /* Portfolio parameters */
  p = 3;
  rho = 10.0;
  comm.iuser = &p;
  ruser[0] = rho;
  comm.user = ruser;

  /* Call MINLP solver h02dac (nag_mip_sqp) */
  nag_mip_sqp(n, nclin, ncnln, a, pda, d, ax, bl, bu, varcon, x, confun, c,
              cjac, objfun, objgrd, maxit, acc, &objmip, iopts, opts, &comm,
              &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_mip_sqp (h02dac).\n%s\n", fail.message);
    exit_status = -1;
    goto END;
  }

  /* Query the accuracy of the mixed integer QP solver */
  nag_mip_opt_get("QP Accuracy", &ivalue, &accqp, cvalue, lcvalue, &optype,
                  iopts, opts, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_mip_opt_get (h02zlc).\n%s\n", fail.message);
    exit_status = -1;
    goto END;
  }

  /* Results */
  printf("\nFinal estimate:");
  for (i = 0; i < n; i++) {
    printf("\nx[%4" NAG_IFMT "] = %12.4f", i + 1, x[i]);
  }
  printf("\n\nRequested accuracy of QP subproblems = %12.4g\n", accqp);
  printf("\nOptimised value = %12.4g\n", objmip);

END:
  NAG_FREE(a);
  NAG_FREE(d);
  NAG_FREE(ax);
  NAG_FREE(bl);
  NAG_FREE(bu);
  NAG_FREE(varcon);
  NAG_FREE(x);
  NAG_FREE(c);
  NAG_FREE(cjac);
  NAG_FREE(objgrd);

  return (exit_status);
}

static void NAG_CALL confun(Integer *mode, Integer ncnln, Integer n,
                            const Integer varcon[], const double x[],
                            double c[], double cjac[], Integer nstate,
                            Nag_Comm *comm)
{
  Integer p;
  double rho;

  /* This problem has two nonlinear constraints.
   * As an example of using the mode mechanism,
   * terminate if any other size is supplied.
   */
  if (ncnln != 2) {
    *mode = -1;
    return;
  }

  if (nstate == 1)
    printf("\n(confun was just called for the first time)\n");

  if (*mode == 0) {
    /* Constraints */
    /* Mean return at least rho: */
    rho = comm->user[0];
    c[0] = 8.0 * x[0] + 9.0 * x[1] + 12.0 * x[2] + 7.0 * x[3] - rho;
    /* Maximum of p assets in portfolio: */
    p = *(comm->iuser);
    c[1] = (double) p - x[4] - x[5] - x[6] - x[7];
  }
  else {
    /* Jacobian */
    /* c[0] */
    CJAC(1, 1) = 8.0;
    CJAC(1, 2) = 9.0;
    CJAC(1, 3) = 12.0;
    CJAC(1, 4) = 7.0;
    /* c[1] does not include continuous variables which requires
       that their derivatives are zero */
    CJAC(2, 1) = CJAC(2, 2) = CJAC(2, 3) = CJAC(2, 4) = 0.0;
  }
}

static void NAG_CALL objfun(Integer *mode, Integer n, const Integer varcon[],
                            const double x[], double *objmip, double objgrd[],
                            Integer nstate, Nag_Comm *comm)
{
  /* This is an 8-dimensional problem.
   * As an example of using the mode mechanism,
   * terminate if any other size is supplied.
   */
  if (n != 8) {
    *mode = -1;
    return;
  }

  if (nstate == 1 || comm->user[1] == -1.0) {
    printf("\n(objfun was just called for the first time)\n");
    comm->user[1] = 0.0;
  }

  if (*mode == 0) {
    /* Objective value */
    *objmip =
           x[0] * (4.0 * x[0] + 3.0 * x[1] - x[2]) + x[1] * (3.0 * x[0] +
                                                             6.0 * x[1] +
                                                             x[2]) +
           x[2] * (x[1] - x[0] + 10.0 * x[2]);
  }
  else {
    /* Objective gradients for continuous variables */
    objgrd[0] = 8.0 * x[0] + 6.0 * x[1] - 2.0 * x[2];
    objgrd[1] = 6.0 * x[0] + 12.0 * x[1] + 2.0 * x[2];
    objgrd[2] = 2.0 * (x[1] - x[0]) + 20.0 * x[2];
    objgrd[3] = 0.0;
  }
}