/* nag_glopt_bnd_pso (e05sac) Example Program.
 *
 * Copyright 2011, Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */
#include <string.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage05.h>
#include <nagx02.h>
#include <nagx04.h>

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL objfun_schwefel(Integer *mode, Integer ndim,
                                     const double x[], double *objf,
                                     double vecout[], Integer nstate,
                                     Nag_Comm *comm);
static void NAG_CALL monmod(Integer ndim, Integer npar, double x[],
                            const double xb[], double fb, const double xbest[],
                            const double fbest[], const Integer itt[],
                            Nag_Comm *comm, Integer *inform);
static void display_result(Integer ndim, const double xb[], const double 
                           x_target[], double fb, double f_target,
                           const Integer itt[], Integer  inform);
static void display_option(const char *optstr, Nag_VariableType optype,
                           Integer ivalue, double rvalue, const char *cvalue);
static void get_known_solution(Integer ndim, double x_target[], 
                               double *f_target);
#ifdef __cplusplus
}
#endif

/* Global constants - set the behaviour of the monitoring function.*/
static const Integer detail_level = 0;
static const Integer report_freq = 100;

int main(void)
{
  /* This example program demonstrates how to use
   * nag_glopt_bnd_pso (e05sac) in standard execution, and with a
   * selection of coupled local minimizers.
   * The non-default option 'Repeatability = On' is used here, giving
   * repeatable results.
   */
  /* Scalars */
  Integer          ndim = 2, npar = 5;
  Integer          exit_status = 0, lcvalue = 17;
  Integer          liopts = 100, lopts = 100;
  double           fb, f_target, rvalue;
  Integer          i, inform, ivalue;
  /* Arrays */
  char             cvalue[17], optstr[81];
  double           *bl=0, *bu=0, opts[100], *xb=0, *x_target=0;
  Integer          iopts[100], itt[6];
  /* NAG types */
  Nag_VariableType optype;
  NagError         fail;
  Nag_Comm         *comm;

  /* Print advisory information.*/
  printf("nag_glopt_bnd_pso (e05sac) Example Program Results\n\n");
  printf("Minimization of the Schwefel function.\n\n");

  /* Allocate memory.*/
  if (!(bl = NAG_ALLOC(ndim, double)) || !(bu = NAG_ALLOC(ndim, double)) ||
      !(xb = NAG_ALLOC(ndim, double)) || !(x_target = NAG_ALLOC(ndim, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Store the known solution of the problem for a comparison.*/
  get_known_solution(ndim, x_target, &f_target);
  
  /* Set box bounds for problem.*/
  for (i = 0; i < ndim; i++)
    {
      bl[i] = -500.0;
      bu[i] = 500.0;
    }

  /* Initialize fail structures.*/
  INIT_FAIL(fail);
  
  /* Initialize comm.*/
  comm = NAGCOMM_NULL;

  /* Initialize the option arrays for nag_glopt_bnd_pso (e05sac)
   * using nag_glopt_opt_set (e05zkc).
   */
  nag_glopt_opt_set("Initialize = e05sac", iopts, liopts, opts, lopts,
                    NAGERR_DEFAULT);

  /* Query some default option values.*/
  printf("  Default Option Queries:\n\n");
  /* nag_glopt_opt_get (e05zlc).
     Option getting routine for nag_glopt_bnd_pso (e05sac).
   */
  strcpy(optstr, "Boundary");
  nag_glopt_opt_get(optstr, &ivalue, &rvalue, cvalue, lcvalue, &optype,
                    iopts, opts, NAGERR_DEFAULT);
  display_option(optstr, optype, ivalue, rvalue, cvalue);

  strcpy(optstr, "Maximum Iterations Completed");
  nag_glopt_opt_get(optstr, &ivalue, &rvalue, cvalue, lcvalue, &optype,
                    iopts, opts, NAGERR_DEFAULT);
  display_option(optstr, optype, ivalue, rvalue, cvalue);

  strcpy(optstr, "Distance Tolerance");
  nag_glopt_opt_get(optstr, &ivalue, &rvalue, cvalue, lcvalue, &optype,
                    iopts, opts, NAGERR_DEFAULT);
  display_option(optstr, optype, ivalue, rvalue, cvalue);

  /* ------------------------------------------------------------------*/

  printf("\n1. Solution without using coupled local minimizer.\n\n");

  /* Set various options to non-default values if required.*/
  nag_glopt_opt_set("Repeatability = On", iopts, liopts, opts, lopts,
                    NAGERR_DEFAULT);
  nag_glopt_opt_set("Verify Gradients = Off", iopts, liopts, opts, lopts,
                    NAGERR_DEFAULT);
  nag_glopt_opt_set("Boundary = Hyperspherical", iopts, liopts, opts, lopts,
                    NAGERR_DEFAULT);
  nag_glopt_opt_set("Maximum iterations static = 150", iopts, liopts, opts,
                    lopts, NAGERR_DEFAULT);;
  nag_glopt_opt_set("Repulsion Initialize = 30", iopts, liopts, opts, lopts,
                    NAGERR_DEFAULT);
  nag_glopt_opt_set("Repulsion Finalize = 30", iopts, liopts, opts, lopts,
                    NAGERR_DEFAULT);

  /* nag_glopt_bnd_pso (e05sac).
   * Global optimization using particle swarm algorithm (PSO),
   * bound constraints only.
   */
  nag_glopt_bnd_pso(ndim, npar, xb, &fb, bl, bu, objfun_schwefel, monmod,
                    iopts, opts, comm, itt, &inform, &fail);
  /* It is essential to test fail.code on exit.*/
  switch (fail.code)
    {
    case NE_NOERROR:
    case NW_FAST_SOLUTION:
    case NW_SOLUTION_NOT_GUARANTEED:
      /* No errors, best found solution at xb returned in fb.*/
      display_result(ndim, xb, x_target, fb, f_target, itt, inform);
      break;
    case NE_USER_STOP:
      /* Exit flag set in objfun or monmod and returned in inform.*/
      display_result(ndim, xb, x_target, fb, f_target, itt, inform);
      break;
    default: /* An error was detected.*/
      exit_status = 1;
      printf("Error from nag_glopt_bnd_pso (e05sac)\n%s\n", fail.message);
      goto END;
    }

  /* ------------------------------------------------------------------*/

  printf("2. Solution using coupled local minimizer "
         "nag_opt_simplex_easy (e04cbc).\n\n");

  /* Initialize fail structures.*/
  INIT_FAIL(fail);
  
  /* Set an objective target.*/
  sprintf(optstr, "Target Objective Value = %32.16e", f_target);
  nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT);
  sprintf(optstr, "Target Objective Tolerance = %32.16e", 1.0e-5);
  nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT);
  sprintf(optstr, "Target Objective Safeguard = %32.16e", 1.0e-8);
  nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT);

  /* Set the local minimizer to be nag_opt_simplex_easy (e04cbc)
   * and set corresponding options.
   */
  nag_glopt_opt_set("Local Minimizer = e04cbc", iopts, liopts, opts, lopts,
                    NAGERR_DEFAULT);
  nag_glopt_opt_set("Local Interior Iterations = 10", iopts, liopts, opts,
                    lopts, NAGERR_DEFAULT);
  nag_glopt_opt_set("Local Exterior Iterations = 20", iopts, liopts, opts,
                    lopts, NAGERR_DEFAULT);
  sprintf(optstr, "Local Interior Tolerance = %32.16e", 1.0e-4);
  nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT);
  sprintf(optstr, "Local Exterior Tolerance = %32.16e", 1.0e-4);
  nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT);

  /* Search for the global optimum.*/
  nag_glopt_bnd_pso(ndim, npar, xb, &fb, bl, bu, objfun_schwefel, monmod,
                    iopts, opts, comm, itt, &inform, &fail);
  /* It is essential to test fail.code on exit.*/
  switch (fail.code)
    {
    case NE_NOERROR:
    case NW_FAST_SOLUTION:
    case NW_SOLUTION_NOT_GUARANTEED:
      /* No errors, best found solution at xb returned in fb.*/
      display_result(ndim, xb, x_target, fb, f_target, itt, inform);
      break;
    case NE_USER_STOP:
      /* Exit flag set in objfun or monmod and returned in inform.*/
      display_result(ndim, xb, x_target, fb, f_target, itt, inform);
      break;
    default: /* An error was detected.*/
      exit_status = 1;
      printf("Error from nag_glopt_bnd_pso (e05sac)\n%s\n", fail.message);
      goto END;
    }

  /* -----------------------------------------------------------------*/

  printf("3. Solution using coupled local minimizer "
         "nag_opt_conj_grad (e04dgc).\n\n");

  /* Initialize fail structures.*/
  INIT_FAIL(fail);
  
  /* set the local minimizer to be nag_opt_conj_grad (e04dgc)
   * and set corresponding options.
   */
  nag_glopt_opt_set("Local Minimizer = e04dgc", iopts, liopts, opts, lopts,
                    NAGERR_DEFAULT);
  nag_glopt_opt_set("Local Interior Iterations = 5", iopts, liopts, opts,
                    lopts, NAGERR_DEFAULT);
  nag_glopt_opt_set("Local Exterior Iterations = 20", iopts, liopts, opts,
                    lopts, NAGERR_DEFAULT);

  /* Search for the global optimum.*/
  nag_glopt_bnd_pso(ndim, npar, xb, &fb, bl, bu, objfun_schwefel, monmod,
                    iopts, opts, comm, itt, &inform, &fail);
  /* It is essential to test fail.code on exit.*/
  switch (fail.code)
    {
    case NE_NOERROR:
    case NW_FAST_SOLUTION:
    case NW_SOLUTION_NOT_GUARANTEED:
      /* No errors, best found solution at xb returned in fb.*/
      display_result(ndim, xb, x_target, fb, f_target, itt, inform);
      break;
    case NE_USER_STOP:
      /* Exit flag set in objfun or monmod and returned in inform.*/
      display_result(ndim, xb, x_target, fb, f_target, itt, inform);
      break;
    default: /* An error was detected.*/
      exit_status = 1;
      printf("Error from nag_glopt_bnd_pso (e05sac)\n%s\n", fail.message);
      goto END;
    }

 END:
  /* Clean up memory.*/
  if (bl) NAG_FREE(bl);
  if (bu) NAG_FREE(bu);
  if (xb) NAG_FREE(xb);
  if (x_target) NAG_FREE(x_target);

  return exit_status;
}

static void NAG_CALL objfun_schwefel(Integer *mode, Integer ndim,
                                     const double x[], double *objf,
                                     double vecout[], Integer nstate,
                                     Nag_Comm *comm)
{
  /* Objective function routine returning the schwefel function and
   * its gradient.
   */
  Nag_Boolean evalobjf, evalobjg;
  Integer     i;
  /* Test nstate to indicate what stage of computation has been reached.*/
  switch (nstate)
    {
    case 2:
      /* objfun is called for the very first time. */
      break;
    case 1:
      /* objfun is called on entry to a NAG local minimizer. */
      break;
    default: /* This will be the normal value of nstate. */
      ;
    }
  /* Test mode to determine whether to calculate objf and/or objgrd.*/
  evalobjf = Nag_FALSE;
  evalobjg = Nag_FALSE;
  switch (*mode)
    {
    case 0:
    case 5:
      /* Only the value of the objective function is needed.*/
      evalobjf = Nag_TRUE;
      break;
    case 1:
    case 6:
      /* Only the values of the ndim gradients are required.*/
      evalobjg = Nag_TRUE;
      break;
    case 2:
    case 7:
      /* Both the objective function and the ndim gradients are required.*/
      evalobjf = Nag_TRUE;
      evalobjg = Nag_TRUE;
    }
  if (evalobjf)
    {
      /* Evaluate the objective function.*/
      *objf = 0.0;
      for (i = 0; i < ndim; i++) 
        *objf += x[i]*sin(sqrt(fabs(x[i])));
    }
  if (evalobjg)
    {
      /* Calculate the gradient of the objective function,
       * and return the result in vecout.
       */
      for (i = 0; i < ndim; i++)
        {
          vecout[i] = sqrt(fabs(x[i]));
          vecout[i] = sin(vecout[i]) + 0.5*vecout[i]*cos(vecout[i]);
        }
    }
}

static void NAG_CALL monmod(Integer ndim, Integer npar, double x[],
                            const double xb[], double fb, const double xbest[],
                            const double fbest[], const Integer itt[],
                            Nag_Comm *comm, Integer *inform)
{
  Integer i;
  if (detail_level)
    {
      /* Report on the first iteration, and every report_freq iterations.*/
      if (itt[0] == 1 || itt[0]%report_freq == 0)
        {
          printf("* Current global optimum candidate\n");
	  for (i=0; i < ndim; i++)
	    printf("  %2ld  %12.5f\n", i+1, xb[i]);
          printf("* Current global optimum value\n");
          printf("     %13.5f\n\n", fb);
        }
    }
  /* If required set *inform<0 to force exit.*/
  *inform = 0;
}

static void display_option(const char *optstr, Nag_VariableType optype,
                           Integer ivalue, double rvalue, const char *cvalue)
{
  /* Subroutine to query optype and print the appropriate option values.*/
  switch (optype)
    {
    case Nag_Integer:
      printf("%-38s: %13ld\n", optstr, ivalue);
      break;
    case Nag_Real:
      printf("%-38s: %13.4f\n", optstr, rvalue);
      break;
    case Nag_Character:
      printf("%-38s: %13s\n", optstr, cvalue);
      break;
    case Nag_Integer_Additional:
      printf("%-38s: %13ld %16s\n", optstr, ivalue, cvalue);
      break;
    case Nag_Real_Additional:
      printf("%-38s: %13.4f %16s\n", optstr, rvalue, cvalue);
      break;
    default:;
    }
}

static void display_result(Integer ndim, const double xb[], const double 
                           x_target[], double fb, double f_target,
                           const Integer itt[], Integer  inform)
{
  /* Display final results in comparison to known global optimum.*/
  Integer  i;
  
  /* Display final counters.*/
  printf(" Algorithm Statistics\n");
  printf(" --------------------\n");
  printf("%-38s: %13ld\n", "Total complete iterations", itt[0]);
  printf("%-38s: %13ld\n", "Complete iterations since improvement",
         itt[1]);
  printf("%-38s: %13ld\n", "Total particles converged to xb", itt[2]);
  printf("%-38s: %13ld\n", "Total improvements to global optimum",
         itt[3]);
  printf("%-38s: %13ld\n", "Total function evaluations", itt[4]);
  printf("%-38s: %13ld\n\n", "Total particles re-initialized",
         itt[5]);
  /* Display why finalization occurred.*/
  switch (inform)
    {
    case 1:
      printf("Solution Status : Target value achieved\n");
      break;
    case 2:
      printf("Solution Status : Minimum swarm standard deviation obtained\n");
      break;
    case 3:
      printf("Solution Status : Sufficient number of particles converged\n");
      break;
    case 4:
      printf("Solution Status : Maximum static iterations attained\n");
      break;
    case 5:
      printf("Solution Status : Maximum complete iterations attained\n");
      break;
    case 6:
      printf("Solution Status : Maximum function evaluations exceeded\n");
      break;
    default:
      if (inform<0)
        printf("Solution Status : User termination, inform = %16ld\n", 
               inform);
      else
        printf("Solution Status : Termination, an error has been detected\n");
      break;
    }
  /* Display final objective value and location.*/
  printf("  Known objective optimum             : %13.5f\n", f_target);
  printf("  Achieved objective value            : %13.5f\n\n", fb);

  printf("  Comparison between the known optimum and the achieved solution.\n");
  printf("          x_target            xb\n");
  for (i=0; i < ndim; i++)
    printf("  %2ld  %12.2f  %12.2f\n", i+1, x_target[i], xb[i]);
 
  printf("\n");
}

static void get_known_solution(Integer ndim, double x_target[], 
                               double *f_target)
{
  /* Fill in the known solution of multidimensional Schwefel's function. */
  Integer i;

  if (f_target && x_target && ndim>0) 
  {
    *f_target=-418.9828872724337*ndim;
    for (i=0; i < ndim; i++)
      x_target[i]=-420.9687463599820;
  }
}