/* nag_glopt_nlp_pso (e05sbc) 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 confun(Integer *mode, Integer ncon, Integer ndim, 
                            Integer tdcj, const Integer needc[], 
                            const double x[], double c[], double cjac[], 
                            Integer nstate, Nag_Comm *comm);

static void NAG_CALL monmod(Integer ndim, Integer ncon, Integer npar,
                            double x[], const double xb[], double fb, 
                            const double cb[], const double xbest[], 
                            const double fbest[], const double cbest[], 
                            const Integer itt[], Nag_Comm *comm,
                            Integer *inform);

static void display_option(const char *optstr, Nag_VariableType optype,
                           Integer ivalue, double rvalue, const char *cvalue);

static void display_result(Integer ndim, Integer ncon, const double xb[],
                           double fb, const double cb[], const Integer itt[],
                           Integer inform);
#ifdef __cplusplus
}
#endif

/* Global constants.*/
/* Set the behaviour of the monitoring function.*/
static const Integer detail_level = 0;
static const Integer report_freq = 100;
/* Known solution for a comparison.*/
static const double  f_target = -731.70709230672696;
static const double  c_scale[] = { 2490.0, 750000.0, 0.1 };
static const double  c_target[] = { 0.0, 0.0, 0.0 };
static const double  x_target[] = { -394.1470221120988, -433.48214189947606 };

int main(void)
{
  /* This example program demonstrates how to use
   * nag_glopt_nlp_pso (e05sbc) in standard execution, and with
   * nag_opt_nlp (e04ucc) as a coupled local minimizer.
   * The non-default option 'Repeatability = On' is used here, giving
   * repeatable results.
   */
  /* Scalars */
  Integer          ncon = 3, ndim = 2, npar = 20;
  Integer          exit_status = 0, lcvalue = 17;
  Integer          liopts = 100, lopts = 100;
  double           fb, rvalue;
  Integer          i, inform, ivalue;
  /* Arrays */
  char             cvalue[17], optstr[81];
  double           opts[100], *bl=0, *bu=0, *cb=0, *xb=0;
  double           *cbest=0, *fbest=0, *xbest=0;
  Integer          iopts[100], itt[7];
  /* NAG types */
  Nag_VariableType optype;
  Nag_Comm         *comm;
  NagError         fail;

  /* Print advisory information.*/
  printf("nag_glopt_nlp_pso (e05sbc) Example Program Results\n\n");
  printf("Minimization of the Schwefel function.\n");
  printf("Subject to one linear and two nonlinear constraints.\n\n");

  /* Allocate memory for arrays.*/
  if (!(bl = NAG_ALLOC(ndim+ncon, double)) || 
      !(bu = NAG_ALLOC(ndim+ncon, double)) ||
      !(cb = NAG_ALLOC(ncon, double)) ||
      !(cbest = NAG_ALLOC(ncon*npar, double)) ||
      !(fbest = NAG_ALLOC(npar, double)) ||
      !(xb = NAG_ALLOC(ndim, double)) || 
      !(xbest = NAG_ALLOC(ndim*npar, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  for (i=0; i < npar; i++) 
    fbest[i] = 0.0;
  for (i=0; i < npar*ndim; i++) 
    xbest[i] = 0.0;
  for (i=0; i < npar*ncon; i++) 
    cbest[i] = 0.0;

  /* Set problem specific values.*/
  /* Set box bounds.*/
  for (i = 0; i < ndim; i++)
    {
      bl[i] = -500.0;
      bu[i] = 500.0;
    }

  /* Set constraint bounds.*/
  bl[ndim]   = -1.0e6;
  bl[ndim+1] = -1.0;
  bl[ndim+2] = -0.9;
  bu[ndim]   = 10.0;
  bu[ndim+1] = 5.0e5;
  bu[ndim+2] = 0.9;

  /* Initialize NagError structure.*/
  INIT_FAIL(fail);

  /* Initialize comm.*/
  comm = NAGCOMM_NULL;

  /* Initialize the option arrays for nag_glopt_nlp_pso (e05sbc)
   * using nag_glopt_opt_set (e05zkc).
   */
  nag_glopt_opt_set("Initialize = e05sbc", 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_nlp_pso (e05sbc).
   */
  ivalue = 0;
  rvalue = 0.0;
  optype = 0;
  strcpy(optstr, "Constraint Norm");
  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.*/
  sprintf(optstr, "Distance Tolerance = %32.16e", rvalue * 0.1);
  nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT);
  nag_glopt_opt_set("Constraint Tolerance = 1.0e-4", iopts, liopts, opts, lopts,
		    NAGERR_DEFAULT);
  nag_glopt_opt_set("Constraint Norm = Euclidean", iopts, liopts, opts, lopts,
		    NAGERR_DEFAULT);
  nag_glopt_opt_set("Repeatability = On", iopts, liopts, opts, lopts,
		    NAGERR_DEFAULT);
  sprintf(optstr, "Target Objective Value = %32.16e", f_target);
  nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT);
  nag_glopt_opt_set("Target Objective Tolerance = 1.0e-4", iopts, liopts, opts,
		    lopts, NAGERR_DEFAULT);

  /* nag_glopt_nlp_pso (e05sbc).
     Global optimization using particle swarm algorithm (PSO), comprehensive.
   */
  nag_glopt_nlp_pso(ndim, ncon, npar, xb, &fb, cb, bl, bu, xbest,
                    fbest, cbest, objfun_schwefel, confun,
                    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, ncon, xb, fb, cb, itt, inform);
      break;
    case NE_USER_STOP:
      /* Exit flag set in objfun, confun or monmod and returned in inform.*/
      display_result(ndim, ncon, xb, fb, cb, itt, inform);
      break;
    default: /* An error was detected.*/
      exit_status = 1;
      printf("Error from nag_glopt_nlp_pso (e05sbc)\n%s\n", fail.message);
      goto END;
    }

  /* ------------------------------------------------------------------*/
  printf("2. Solution using coupled local minimizer nag_opt_nlp (e04ucc).\n\n");

  /*  Set the local minimizer to be nag_opt_nlp (e04ucc) and set corresponding 
   *  options.
   */
  nag_glopt_opt_set("Local Minimizer = e04ucc", iopts, liopts, opts, lopts,
		    NAGERR_DEFAULT);
  nag_glopt_opt_set("Local Interior Major Iterations = 15", iopts, liopts, opts,
		    lopts, NAGERR_DEFAULT);
  nag_glopt_opt_set("Local Interior Minor Iterations = 5", iopts, liopts, opts,
		    lopts, NAGERR_DEFAULT);
  nag_glopt_opt_set("Local Exterior Major Iterations = 50", iopts, liopts, opts,
		    lopts, NAGERR_DEFAULT);
  nag_glopt_opt_set("Local Exterior Minor Iterations = 15", iopts, liopts, opts,
		    lopts, NAGERR_DEFAULT);

  /* Query the option Distance Tolerance*/
  nag_glopt_opt_get("Distance Tolerance", &ivalue, &rvalue, cvalue, lcvalue,
                    &optype, iopts, opts, NAGERR_DEFAULT);

  /* Adjust Distance Tolerance dependent upon its current value*/
  sprintf(optstr, "Distance Tolerance = %32.16e", rvalue*10.0);
  nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT);
  sprintf(optstr, "Local Interior Tolerance = %32.16e", rvalue);
  nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT);
  sprintf(optstr, "Local Exterior Tolerance = %32.16e", rvalue * 1.0e-4);
  nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT);

  /* nag_glopt_nlp_pso (e05sbc).
     Global optimization using particle swarm algorithm (PSO), comprehensive.
  */
  nag_glopt_nlp_pso(ndim, ncon, npar, xb, &fb, cb, bl, bu, xbest,
                    fbest, cbest, objfun_schwefel, confun,
                    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:
    case NW_NOT_FEASIBLE:
      /* nag_glopt_nlp_pso (e05sbc) encountered no errors during
       * operation, and will have returned the best solution found.
       */
      display_result(ndim, ncon, xb, fb, cb, itt, inform);
      break;
    case NE_USER_STOP:
      /* Exit flag set in objfun, confun or monmod and returned in inform.*/
      display_result(ndim, ncon, xb, fb, cb, itt, inform);
      break;
    default: /* An error was detected.*/
      exit_status = 1;
      printf("Error from nag_glopt_nlp_pso (e05sbc)\n%s\n", fail.message);
      goto END;
    }

 END:
  /* Clean up memory.*/
  if (bl) NAG_FREE(bl);
  if (bu) NAG_FREE(bu);
  if (cb) NAG_FREE(cb);
  if (cbest) NAG_FREE(cbest);
  if (fbest) NAG_FREE(fbest);
  if (xb) NAG_FREE(xb);
  if (xbest) NAG_FREE(xbest);

  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 confun(Integer *mode, Integer ncon, Integer ndim, 
                            Integer tdcj, const Integer needc[], 
                            const double x[], double c[], double cjac[], 
                            Integer nstate, Nag_Comm *comm)
{
  /* Supplies constraints.
   * cjac[(k-1)*tdcj + (i-1)] corresponds to dc[k]/dx[i]
   * for k=1,...,ncon and i=1,...,ndim.
   */
  Integer k;
  Nag_Boolean evalc, evalcjac; 

  /* Test nstate to determine whether the local minimizer is being called
   * for the first time from a new start point.
   */
  if (nstate == 1)
    {
      /* Set any constant elements of the Jacobian matrix.*/
      cjac[0] = 3.0;
      cjac[1] = -2.0;
    }
  /* mode: are constraints, derivatives, or both are required? */
  evalc = (*mode == 0 || *mode == 2) ? Nag_TRUE : Nag_FALSE;
  evalcjac = (*mode == 1 || *mode == 2) ? Nag_TRUE : Nag_FALSE;

  for (k = 0; k < ncon; k++)
    {
      if(needc[k] <= 0)
        continue;

      if (evalc == Nag_TRUE)
        {
          /* Constraint values are required.
           * Only those for which needc is non-zero need be set.
           */
          switch (k)
            {
            case 0:
              c[k] = 3.0*x[0] - 2.0*x[1];
              break;
            case 1:
              c[k] = x[0]*x[0] - x[1]*x[1] + 3.0*x[0]*x[1];
              break;
            case 2:
              c[k] = cos(pow((x[0]/200.0), 2) + (x[1]/100.0));
              break;
            default: 
              c[k] = 0.0;
            }
        }
      if (evalcjac == Nag_TRUE)
        {
          /* Constraint derivatives (cjac) are required.*/
          switch (k)
            {
            case 0:
              /* Constant derivatives set when nstate=1 remain throughout
               * the local minimization. 
               */
              break;
            case 1:
              /* If the constraint derivatives are known and are readily
               * calculated, populate cjac when required.
               */
              cjac[k*tdcj] = 2.0*x[0] + 3.0*x[1];
              cjac[k*tdcj+1] = -2.0*x[1] + 3.0*x[0];
              break;
            default:
              /* Any elements of cjac left unaltered will be approximated
               * using finite differences when required.
               */
              ;
            }
        }
    }
  /* If an immediate exit is required, return *mode<0 */
  *mode = 0;
}

static void NAG_CALL monmod(Integer ndim, Integer ncon, Integer npar,
                            double x[], const double xb[], double fb, 
                            const double cb[], const double xbest[], 
                            const double fbest[], const double cbest[], 
                            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", fb);
          printf("* Current constraint violations\n");
          for (i = 0; i < ncon; i++) 
	    printf("  %2ld  %12.5f\n", i+1, cb[i]);
          printf("\n");
        }
    }
  /* 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, Integer ncon, const double xb[],
                           double fb, const double cb[], 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", "Total particles re-initialized",
         itt[5]);
  printf("%-38s: %13ld\n\n", "Total constraints violated",
         itt[6]);
  /* 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;
    case 7:
      printf("Solution Status : Feasible point located\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 constrained objective optimum : %13.3f\n", f_target);
  printf("  Achieved objective value            : %13.3f\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");

  if (ncon > 0)
    {
      printf("  Comparison between scaled constraint violations.\n");
      printf("          c_target            cb\n");
      for (i = 0; i < ncon; i++)
        printf("  %2ld  %12.5f  %12.5f\n", i+1, c_target[i]/c_scale[i],
               cb[i]/c_scale[i]);

      printf("\n");
    }
}