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

#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 void ex1(void);
static void 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 == 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)
{
  Vprintf("e04unc Example Program Results\n");
  ex1();
  ex2();
  return EXIT_SUCCESS;
}


static void ex1(void)
{
#define NMAX 10
#define MMAX 50
#define NCLIN 10
#define NCNLIN 10
#define MAXBND NMAX+NCLIN+NCNLIN

  /* Local variables */
  double x[NMAX], a[NCLIN][NMAX];
  double f[MMAX], y[MMAX], fjac[MMAX][NMAX];
  double bl[MAXBND], bu[MAXBND];
  double objf;
  Integer tda, tdfjac;
  Integer i, j, m, n, nclin, ncnlin;
  static NagError fail;

  fail.print = TRUE;

  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 <= MMAX && n <= NMAX && nclin <= NCLIN && ncnlin <= NCNLIN)
    {
      tda = NMAX;
      tdfjac = NMAX;
      /* 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 */
      e04unc(m, n, nclin, ncnlin, &a[0][0], tda, bl, bu, y, objfun,
             confun, x, &objf, f, &fjac[0][0], tdfjac, E04_DEFAULT,
             NAGCOMM_NULL, &fail);
    }

} /* 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 void ex2(void)
{
  /* Local variables */
  double x[NMAX], a[NCLIN][NMAX];
  double f[MMAX], y[MMAX], fjac[MMAX][NMAX];
  double bl[MAXBND], bu[MAXBND];
  double objf;
  Integer tda, tdfjac;
  Integer i, j, m, n, nclin, ncnlin;
  Nag_E04_Opt options;
  static NagError fail, fail1;

  fail.print = TRUE;
  fail1.print = TRUE;

  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 <= MMAX && n <= NMAX && nclin <= NCLIN && ncnlin <= NCNLIN)
    {
      tda = NMAX;
      tdfjac = NMAX;
      /* 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 */
      e04xxc(&options);
      options.print_fun = user_print;

      /* Solve the problem */
      e04unc(m, n, nclin, ncnlin, &a[0][0], tda, bl, bu, y, objfun,
             confun, x, &objf, f, &fjac[0][0], tdfjac, &options,
             NAGCOMM_NULL, &fail);

      e04xzc(&options, "all", &fail1);
    }

} /* ex2 */