/* nag_opt_bounds_bobyqa_func (e04jcc) Example Program.
 *
 * Copyright 2019 Numerical Algorithms Group.
 *
 * Mark 27.0, 2019.
 *
 */

#include <nag.h>
#include <math.h>
#include <stdio.h>
#include <string.h>

#ifdef __cplusplus
extern "C"
{
#endif
  static void NAG_CALL objfun(Integer n, const double x[], double *f,
                              Nag_Comm *comm, Integer *inform);
  static void NAG_CALL monfun(Integer n, Integer nf, const double x[],
                              double f, double rho, Nag_Comm *comm,
                              Integer *inform);
#ifdef __cplusplus
}
#endif

int main(void)
{
  static double ruser[2] = { -1.0, -1.0 };
  Integer exit_status = 0;
  double rhobeg, rhoend, f;
  Integer i, n, nf, npt, maxcal;
  double *bl = 0, *bu = 0, *x = 0;
  NagError fail;
  Nag_Comm comm;

  INIT_FAIL(fail);

  printf("nag_opt_bounds_bobyqa_func (e04jcc) Example Program Results\n");

  /* For communication with user-supplied functions: */
  comm.user = ruser;

  maxcal = 500;
  rhobeg = 1.0e-1;
  rhoend = 1.0e-6;
  n = 4;
  npt = 2 * n + 1;

  if (!(x = NAG_ALLOC(n, double)) ||
      !(bl = NAG_ALLOC(n, double)) || !(bu = NAG_ALLOC(n, double)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Set bounds on variables */
  /* x[2] is not bounded, so we set bl[2] to a large negative
   * number and bu[2] to a large positive number
   */
  bl[0] = 1.0;
  bl[1] = -2.0;
  bl[2] = -1.0e10;
  bl[3] = 1.0;
  bu[0] = 3.0;
  bu[1] = 0.0;
  bu[2] = 1.0e10;
  bu[3] = 3.0;
  x[0] = 3.0;
  x[1] = -1.0;
  x[2] = 0.0;
  x[3] = 1.0;

  /* Call optimization routine */
  /* nag_opt_bounds_bobyqa_func (e04jcc).
     Bound-constrained optimization by quadratic approximations. */
  nag_opt_bounds_bobyqa_func(objfun, n, npt, x, bl, bu, rhobeg, rhoend,
                             monfun, maxcal, &f, &nf, &comm, &fail);

  if (fail.code == NE_NOERROR ||
      fail.code == NE_TOO_MANY_FEVALS ||
      fail.code == NE_TR_STEP_FAILED ||
      fail.code == NE_RESCUE_FAILED || fail.code == NE_USER_STOP) {

    if (fail.code == NE_NOERROR) {
      printf("Successful exit.\n");
    }

    printf("Function value at lowest point found is %11.3f\n", f);
    printf("The corresponding x is:");
    for (i = 0; i <= n - 1; ++i) {
      printf(" %11.3f", x[i]);
    }
    printf("\n");
  }
  else {
    exit_status = 1;
  }

  if (fail.code != NE_NOERROR) {
    printf("%s\n", fail.message);
  }

END:
  NAG_FREE(x);
  NAG_FREE(bl);
  NAG_FREE(bu);

  return exit_status;
}

static void NAG_CALL objfun(Integer n, const double x[], double *f,
                            Nag_Comm *comm, Integer *inform)
{
  /* Routine to evaluate objective function. */

  double a, b, c, d, x1, x2, x3, x4;

  if (comm->user[0] == -1.0) {
    printf("(User-supplied callback objfun, first invocation.)\n");
    comm->user[0] = 0.0;
  }
  *inform = 0;
  x1 = x[0];
  x2 = x[1];
  x3 = x[2];
  x4 = x[3];

  /* Supply a single function value */
  a = x1 + 10.0 * x2;
  b = x3 - x4;
  c = x2 - 2.0 * x3, c *= c;
  d = x1 - x4, d *= d;
  *f = a * a + 5.0 * b * b + c * c + 10.0 * d * d;
}

static void NAG_CALL monfun(Integer n, Integer nf, const double x[], double f,
                            double rho, Nag_Comm *comm, Integer *inform)
{
  /* Monitoring routine */
  Integer j;
  Nag_Boolean verbose;

  if (comm->user[1] == -1.0) {
    printf("(User-supplied callback monfun, first invocation.)\n");
    comm->user[1] = 0.0;
  }
  *inform = 0;

  printf("\nMonitoring: new trust region radius = %13.3e\n", rho);
  verbose = Nag_FALSE; /* Set this to Nag_TRUE to get more detailed output */
  if (verbose)
    {
      printf("Number of function evaluations = %16" NAG_IFMT "\n", nf);
      printf("Current function value = %13.5f\n", f);
      printf("The corresponding x is:\n");
      for (j = 0; j <= n - 1; ++j) {
        printf(" %13.5e", x[j]);
      }
      printf("\n");
    }
}