/* nag_opt_nlin_lsq (e04unc) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 *
 * Mark 6 revised, 2000.
 * Mark 7 revised, 2001.
 * Mark 8 revised, 2004.
 *
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <math.h>
#include <nage04.h>

#ifdef __cplusplus
extern "C" {
#endif
  static void objfun(Integer m, Integer n, double x[], double f[],
                     double fjac[], Integer tdfjac, Nag_Comm *comm);
  static void confun(Integer n, Integer ncnlin, Integer needc[], double x[],
                     double conf[], double cjac[], Nag_Comm *comm);
  static void user_print(const Nag_Search_State *st, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
static int ex1(void);
static int ex2(void);
static void objfun(Integer m, Integer n, double x[], double f[],
                   double fjac[], Integer tdfjac, Nag_Comm *comm)
{
#define FJAC(I,J) fjac[(I)*tdfjac + (J)]

  /* Initialized data */
  static double a[44] = {
    8.0, 8.0, 10.0, 10.0, 10.0, 10.0, 12.0, 12.0, 12.0, 12.0, 14.0, 14.0,
    14.0, 16.0, 16.0, 16.0, 18.0, 18.0, 20.0, 20.0, 20.0, 22.0, 22.0, 22.0,
    24.0, 24.0, 24.0, 26.0, 26.0, 26.0, 28.0, 28.0, 30.0, 30.0, 30.0, 32.0,
    32.0, 34.0, 36.0, 36.0, 38.0, 38.0, 40.0, 42.0 };


  /* Local variables */
  double temp;
  Integer i;
  double x0, x1, ai;


  /* Function to evaluate the objective subfunctions
   * and their 1st derivatives.
   */
  x0 = x[0];
  x1 = x[1];
  for (i = 0; i < m; ++i)
    {
      ai = a[i];
      temp = exp(-x1 * (ai - 8.0));
      /* Evaluate objective subfunction f(i+1) only if required */
      if (comm->needf == i+1 || comm->needf == 0)
        f[i] = x0 + (.49 - x0) * temp;
      if (comm->flag == 2)
        {
          FJAC(i,0) = 1.0 - temp;
          FJAC(i,1) = -(.49 - x0) * (ai - 8.0) * temp;
        }
    }
} /* objfun */
static void confun(Integer n, Integer ncnlin, Integer needc[], double x[],
                   double conf[], double cjac[], Nag_Comm *comm)
{
#define CJAC(I,J) cjac[(I)*n + (J)]

  /* Function to evaluate the nonlinear constraints and its 1st derivatives. */

  if (comm->first == Nag_TRUE)
    {
      /* First call to confun.  Set all Jacobian elements to zero.
       *  Note that this will only work when  options.obj_deriv = TRUE
       *  (the default).
       */
      CJAC(0,0) = CJAC(0,1) = 0.0;
    }

  if (needc[0] > 0)
    {
      conf[0] = -0.09 - x[0]*x[1] + 0.49*x[1];

      if (comm->flag == 2)
        {
          CJAC(0,0) = -x[1];
          CJAC(0,1) = -x[0] + 0.49;
        }
    }
} /* confun */


int main(void)
{

  Integer exit_status_ex1=0;
  Integer exit_status_ex2=0;

  Vprintf("nag_opt_nlin_lsq (e04unc) Example Program Results\n");

  exit_status_ex1 = ex1();
  exit_status_ex2 = ex2();

  return exit_status_ex1 == 0 && exit_status_ex2 == 0 ? 0 : 1;

}

#undef FJAC
#define A(I,J) a[(I)*tda + J]
#define FJAC(I,J) fjac[(I)*tdfjac + J]

static int ex1(void)
{

  /* Local variables */
  Integer exit_status=0, i, j, m, n, nbnd, nclin, ncnlin, tda, tdfjac;
  NagError fail;
  double *a=0, *bl=0, *bu=0, *f=0, *fjac=0, objf, *x=0, *y=0;

  INIT_FAIL(fail);

  Vprintf("\nExample 1: default options\n");
  Vscanf(" %*[^\n]"); /* Skip heading in data file */
  Vscanf(" %*[^\n]");

  /* Read problem dimensions */
  Vscanf(" %*[^\n]");
  Vscanf("%ld%ld%*[^\n]", &m, &n);
  Vscanf(" %*[^\n]");
  Vscanf("%ld%ld%*[^\n]", &nclin, &ncnlin);

  if (m>0 && n>0 && nclin>=0 && ncnlin>=0)
    {
      nbnd =n + nclin + ncnlin;
      if ( !( x = NAG_ALLOC(n, double)) ||
           !( a = NAG_ALLOC(nclin*n, double)) ||
           !( f = NAG_ALLOC(m, double)) ||
           !( y = NAG_ALLOC(m, double)) ||
           !( fjac = NAG_ALLOC(m*n, double)) ||
           !( bl = NAG_ALLOC(nbnd, double)) ||
           !( bu = NAG_ALLOC(nbnd, double)) )
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
      tda = n;
      tdfjac = n;
    }
  else
    {
      Vprintf("Invalid m or n or nclin or ncnlin.\n");
      exit_status = 1;
      return exit_status;
    }
  /* Read a, y, bl, bu and x from data file */

  if (nclin > 0)
    {
      Vscanf(" %*[^\n]");
      for (i = 0; i < nclin; ++i)
        for (j = 0; j < n; ++j)
          Vscanf("%lf",&A(i,j));
    }

  /* Read the y vector of the objective */
  Vscanf(" %*[^\n]");
  for (i = 0; i < m; ++i)
    Vscanf("%lf",&y[i]);

  /* Read lower bounds */
  Vscanf(" %*[^\n]");
  for (i = 0; i < n + nclin + ncnlin; ++i)
    Vscanf("%lf",&bl[i]);

  /* Read upper bounds */
  Vscanf(" %*[^\n]");
  for (i = 0; i < n + nclin + ncnlin; ++i)
    Vscanf("%lf",&bu[i]);

  /* Read the initial point x */
  Vscanf(" %*[^\n]");
  for (i = 0; i < n; ++i)
    Vscanf("%lf",&x[i]);

  /* Solve the problem */
  /* nag_opt_nlin_lsq (e04unc).
   * Solves nonlinear least-squares problems using the
   * sequential QP method
   */
  nag_opt_nlin_lsq(m, n, nclin, ncnlin, a, tda, bl, bu, y, objfun,
                   confun, x, &objf, f, fjac, tdfjac, E04_DEFAULT,
                   NAGCOMM_NULL, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_opt_nlin_lsq (e04unc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
 END:
  if (x) NAG_FREE(x);
  if (a) NAG_FREE(a);
  if (f) NAG_FREE(f);
  if (y) NAG_FREE(y);
  if (fjac) NAG_FREE(fjac);
  if (bl) NAG_FREE(bl);
  if (bu) NAG_FREE(bu);
  return exit_status;
} /* ex1 */
static void user_print(const Nag_Search_State *st, Nag_Comm *comm)
{
  static const char *states[] = {"IL", "IU", "FR", "LL", "UL", "EQ"};
  Integer i, ixl, ixn;
  Nag_SimSt c_sim, f_sim;


  /* Derivative checking (only handle simple checks here ) */
  if (comm->g_prt)
    {
      /* First ensure that a check has been performed */
      if (st->gprint->g_chk.type != Nag_NoCheck)
        {
          Vprintf("\n\nDerivative Checks\n-----------------\n\n");

          /* Objective check */
          f_sim = st->gprint->f_sim;
          Vprintf("Simple check of objective gradients: ");
          if (f_sim.correct)
            Vprintf("OK\n");
          else
            Vprintf("ERROR!\n");

          Vprintf("Max error was %9.2e in subfunction %1ld\n\n",
                  f_sim.max_error, f_sim.max_subfunction);

          /* Similarly for constraints */
          c_sim = st->gprint->c_sim;
          Vprintf("Simple check of constraint gradients: ");
          if (c_sim.correct)
            Vprintf("OK\n");
          else
            Vprintf("ERROR!\n");

          Vprintf("Max error was %9.2e in constraint %1ld\n",
                  c_sim.max_error, c_sim.max_constraint);
        }
    }

  /* Major iteration output */
  if (comm->it_maj_prt)
    {
      if (st->first)
        {
          Vprintf("\nIterations\n----------\n");
          Vprintf("\n   Maj  Mnr    Step   Merit function\n");
        }
      Vprintf(" %4ld %4ld %8.1e  %14.6e\n", st->iter,
              st->minor_iter, st->step, st->merit);
    }

  /* Solution output */
  if (comm->sol_sqp_prt)
    {
      Vprintf("\nSolution\n--------\n");

      /* Variable results */
      Vprintf("\nVarbl State     Value      Lagr Mult\n");
      for (i = 0; i < st->n; ++i)
        Vprintf("V%4ld  %2s %14.6e %12.4e\n", i+1,
                states[st->state[i] + 2], st->x[i], st->lambda[i]);

      /* Linear constraint results */
      Vprintf("\nL Con State     Value      Lagr Mult\n");
      for (i = st->n; i < st->n + st->nclin; ++i)
        {
          ixl = i - st->n;
          Vprintf("V%4ld  %2s %14.6e %12.4e\n", ixl+1,
                  states[st->state[i] + 2], st->ax[ixl], st->lambda[i]);
        }

      /* Nonlinear constraint results */
      Vprintf("\nN Con State     Value      Lagr Mult\n");
      for (i = st->n + st->nclin; i < st->n + st->nclin + st->ncnlin; ++i)
        {
          ixn = i - st->n - st->nclin;
          Vprintf("V%4ld  %2s %14.6e %12.4e\n", ixn+1,
                  states[st->state[i] + 2], st->cx[ixn], st->lambda[i]);
        }
    }
} /* user_print */


static int ex2(void)
{
  /* Local variables */
  Integer exit_status=0, i, j, m, n, nbnd, nclin, ncnlin, tda, tdfjac;
  NagError fail;
  Nag_E04_Opt options;
  double *a=0, *bl=0, *bu=0, *f=0, *fjac=0, objf, *x=0, *y=0;

  INIT_FAIL(fail);


  Vprintf("\nExample 2: user defined printing option\n");
  Vscanf(" %*[^\n]"); /* Skip heading in data file */

  /* Read problem dimensions */
  Vscanf(" %*[^\n]");
  Vscanf("%ld%ld%*[^\n]", &m, &n);
  Vscanf(" %*[^\n]");
  Vscanf("%ld%ld%*[^\n]", &nclin, &ncnlin);

  if (m>0 && n>0 && nclin>=0 && ncnlin>=0)
    {
      nbnd =  n + nclin + ncnlin;
      if ( !( x = NAG_ALLOC(n, double)) ||
           !( a = NAG_ALLOC(nclin*n, double)) ||
           !( f = NAG_ALLOC(m, double)) ||
           !( y = NAG_ALLOC(m, double)) ||
           !( fjac = NAG_ALLOC(m*n, double)) ||
           !( bl = NAG_ALLOC(nbnd, double)) ||
           !( bu = NAG_ALLOC(nbnd, double)) )
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
      tda = n;
      tdfjac = n;
    }
  else
    {
      Vprintf("Invalid m or n nclin or ncnlin.\n");
      exit_status = 1;
      return exit_status;
    }
  /* Read a, y, bl, bu and x from data file */

  if (nclin > 0)
    {
      Vscanf(" %*[^\n]");
      for (i = 0; i < nclin; ++i)
        for (j = 0; j < n; ++j)
          Vscanf("%lf",&A(i,j));
    }

  /* Read the y vector of the objective */
  Vscanf(" %*[^\n]");
  for (i = 0; i < m; ++i)
    Vscanf("%lf",&y[i]);

  /* Read lower bounds */
  Vscanf(" %*[^\n]");
  for (i = 0; i < n + nclin + ncnlin; ++i)
    Vscanf("%lf",&bl[i]);

  /* Read upper bounds */
  Vscanf(" %*[^\n]");
  for (i = 0; i < n + nclin + ncnlin; ++i)
    Vscanf("%lf",&bu[i]);

  /* Read the initial point x */
  Vscanf(" %*[^\n]");
  for (i = 0; i < n; ++i)
    Vscanf("%lf",&x[i]);

  /* Set an option */
  /* nag_opt_init (e04xxc).
   * Initialization function for option setting
   */
  nag_opt_init(&options);
  options.print_fun = user_print;

  /* Solve the problem */
  /* nag_opt_nlin_lsq (e04unc), see above. */
  nag_opt_nlin_lsq(m, n, nclin, ncnlin, a, tda, bl, bu, y, objfun,
                   confun, x, &objf, f, fjac, tdfjac, &options,
                   NAGCOMM_NULL, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_opt_nlin_lsq (e04unc).\n%s\n", fail.message);
      exit_status = 1;
    }
  /* nag_opt_free (e04xzc).
   * Memory freeing function for use with option setting
   */
  nag_opt_free(&options, "all", &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_opt_free (e04xzc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
 END:
  if (x) NAG_FREE(x);
  if (a) NAG_FREE(a);
  if (f) NAG_FREE(f);
  if (y) NAG_FREE(y);
  if (fjac) NAG_FREE(fjac);
  if (bl) NAG_FREE(bl);
  if (bu) NAG_FREE(bu);
  return exit_status;
} /* ex2 */