/* nag_glopt_bnd_mcs_solve (e05jbc) Example Program.
 *
 * Copyright 2006 Numerical Algorithms Group.
 *
 * Mark 9, 2009.
 */
#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage05.h>

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL objfun(Integer n, const double x[], double *f,
                            Integer nstate, Nag_Comm *comm, Integer *inform);
static void NAG_CALL monit(Integer n, Integer ncall, const double xbest[],
                           const Integer icount[], Integer sdlist,
                           const double list[], const Integer numpts[],
                           const Integer initpt[], Integer nbaskt,
                           const double xbaskt[], const double boxl[],
                           const double boxu[], Integer nstate, Nag_Comm *comm,
                           Integer *inform);
static void NAG_CALL output_current_box(Integer n, const double boxl[],
                                        const double boxu[]);
#ifdef __cplusplus
}
#endif

int main(void)
{
  /* Scalars */
  double            obj;
  Integer           exit_status, i, n, plot, sdlist;
  Nag_BoundType     boundenum;
  Nag_MCSInitMethod initmethodenum;
  /* Arrays */
  char              bound[16], initmethod[18];
  double            *bl = 0, *bu = 0, *list = 0, *x = 0;
  Integer           *initpt = 0, *numpts = 0;
  Integer           iuser[1];
  /* Nag Types */
  Nag_E05State      state;
  NagError          fail;
  Nag_Comm          comm;

  comm.iuser = iuser;

  exit_status = 0;
  INIT_FAIL(fail);

  printf("nag_glopt_bnd_mcs_solve (e05jbc) Example Program Results\n");

  comm.p = (Pointer) "stdout";

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  /* Read n and sdlist from data file */
  scanf("%ld%ld%*[^\n] ", &n, &sdlist);
  if (n > 0 && sdlist > 0)
    {
      /* Allocate memory */
      if (!(bl = NAG_ALLOC(n, double)) ||
          !(bu = NAG_ALLOC(n, double)) ||
          !(list = NAG_ALLOC(n*sdlist, double)) ||
          !(x = NAG_ALLOC(n, double)) ||
          !(initpt = NAG_ALLOC(n, Integer)) ||
          !(numpts = NAG_ALLOC(n, Integer)))
        {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }

      /* Read in bound (and bl and bu if necessary) */
      scanf("%s%*[^\n] ", bound);

      /*
       * nag_enum_name_to_value (x04nac).
       * Converts NAG enum member name to value
       */
      boundenum = (Nag_BoundType) nag_enum_name_to_value(bound);

      if (boundenum == Nag_Bounds)
      /* Read in the whole of each bound */
        {
          for (i = 0; i < n; ++i)
            scanf("%lf", &bl[i]);
          scanf("%*[^\n] ");

          for (i = 0; i < n; ++i)
            scanf("%lf", &bu[i]);
          scanf("%*[^\n] ");
        }
      else if (boundenum == Nag_BoundsEqual)
      /* Bounds are uniform: read in only the first entry of each */
        {
          scanf("%lf%*[^\n] ", &bl[0]);
          scanf("%lf%*[^\n] ", &bu[0]);
        }

      /* Read in initmethod */
      scanf("%s%*[^\n] ", initmethod);

      /*
       * nag_enum_name_to_value (x04nac).
       * Converts NAG enum member name to value
       */
      initmethodenum = (Nag_MCSInitMethod) nag_enum_name_to_value(initmethod);

      /* Read in plot. Its value determines whether monit displays
         information on the current search box */
      scanf("%ld%*[^\n] ", &plot);

      /* Communicate plot through to monit */
      iuser[0] = plot;

      /* Call nag_glopt_bnd_mcs_init (e05jac) to initialize
       * nag_glopt_bnd_mcs_solve (e05jbc). */
      /* Its first argument is a legacy argument and has no significance. */
      nag_glopt_bnd_mcs_init(n,&state, &fail);

      if (fail.code != NE_NOERROR)
        {
          printf("Initialization of nag_glopt_bnd_mcs_init (e05jac) "
                  "failed.\n%s\n", fail.message);
          exit_status = 1;
          goto END;
        }

      /* Solve the problem. */
      /* nag_glopt_bnd_mcs_solve (e05jbc).
       * Global optimization by multilevel coordinate search, simple bounds.
       */
      nag_glopt_bnd_mcs_solve(n, objfun, boundenum, initmethodenum, bl, bu,
                              sdlist, list, numpts, initpt, monit, x, &obj,
                              &state, &comm, &fail);

      if (fail.code == NE_NOERROR)
        {
          printf("Final objective value = %11.5f\n", obj);
          printf("Global optimum X = ");

          for (i = 0; i < n; ++i)
            printf("%9.5f", x[i]);
          printf("\n");
        }
      else
        {
          printf(
                  "Error message from nag_glopt_bnd_mcs_solve (e05jbc).\n%s\n",
                  fail.message);
          exit_status = 1;
        }
    }
 END:
  if (bl) NAG_FREE(bl);
  if (bu) NAG_FREE(bu);
  if (list) NAG_FREE(list);
  if (x) NAG_FREE(x);
  if (initpt) NAG_FREE(initpt);
  if (numpts) NAG_FREE(numpts);

  return exit_status;
}

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

  *inform = 0;

  if (*inform >= 0)
  /* Here we're prepared to evaluate objfun at the current x */
    {
      if (nstate == 1)
      /* This is the first call to objfun */
        {
          printf("\n(objfun was just called for the first time)\n");
        }

      *f = (
        3.0*pow((1.0-x[0]), 2)*exp(-pow(x[0], 2)-pow((x[1]+1), 2))
        - (10.0*(x[0]/5.0-pow(x[0], 3)-pow(x[1], 5))*
           exp(-pow(x[0], 2)-pow(x[1], 2)))
        - 1.0/3.0*exp(-pow((x[0]+1.0), 2)-pow(x[1], 2))
        );
    }
}

static void NAG_CALL monit(Integer n, Integer ncall, const double xbest[],
                           const Integer icount[], Integer sdlist,
                           const double list[], const Integer numpts[],
                           const Integer initpt[], Integer nbaskt,
                           const double xbaskt[], const double boxl[],
                           const double boxu[], Integer nstate, Nag_Comm *comm,
                           Integer *inform)
{
  /* Scalars */
  Integer i, j;
  Integer plot;

#define XBASKT(I, J) xbaskt[(I-1)*nbaskt + (J-1)]

  *inform = 0;

  if (*inform >= 0)
  /* We are going to allow the iterations to continue */
    {
      /* Extract plot from the communication structure */
      plot = comm->iuser[0];

      if (nstate == 0 || nstate == 1)
      /* When nstate == 1, monit is called for the first time.
       * When nstate == 0, monit is called for the first AND last time.
       * Display a welcome message */
        {
          printf("\n*** Begin monitoring information ***\n\n");

          if (plot != 0 && n == 2)
            printf("<Begin displaying search boxes>\n\n");
        }

      if (plot != 0 && n == 2)
        {
          /* Display the coordinates of the edges of the current search box */
          output_current_box(n, boxl, boxu);
        }

      if (nstate <= 0)
      /* monit is called for the last time */
        {
          if (plot != 0 && n == 2)
            printf("<End displaying search boxes>\n\n");
          printf("Total sub-boxes = %5ld\n", icount[0]);
          printf("Total function evaluations = %5ld\n", ncall);
          printf("Total function evaluations used in local search = %5ld\n",
                  icount[1]);
          printf("Total points used in local search = %5ld\n",
                  icount[2]);
          printf("Total sweeps through levels = %5ld\n", icount[3]);
          printf("Total splits by init. list = %5ld\n", icount[4]);
          printf("Lowest level with nonsplit boxes = %5ld\n",
                  icount[5]);
          printf("Number of candidate minima in the 'shopping basket'"
                  " = %5ld\n", nbaskt);
          printf("Shopping basket:\n");

          for (i = 1; i <= n; ++i)
            {
              printf("xbaskt(%3ld,:) =", i);
              for (j = 1; j <= nbaskt; ++j)
                printf("%9.5f", XBASKT(i, j));
              printf("\n");
            }

          printf("\n");
          printf("*** End monitoring information ***\n");
          printf("\n");
        }
    }
}

static void NAG_CALL output_current_box(Integer n, const double boxl[],
                                        const double boxu[])
{
  printf("%20.15f %20.15f\n", boxl[0], boxl[1]);
  printf("%20.15f %20.15f\n\n", boxl[0], boxu[1]);
  printf("%20.15f %20.15f\n", boxl[0], boxl[1]);
  printf("%20.15f %20.15f\n\n", boxu[0], boxl[1]);
  printf("%20.15f %20.15f\n", boxl[0], boxu[1]);
  printf("%20.15f %20.15f\n\n", boxu[0], boxu[1]);
  printf("%20.15f %20.15f\n", boxu[0], boxl[1]);
  printf("%20.15f %20.15f\n\n", boxu[0], boxu[1]);
}