/* nag_rand_field_2d_generate (g05zsc) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 24, 2013.
 */

#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg05.h>
#include <nagx04.h>

#define NPMAX 5
#define LENST 17
#define LSEED 1
static void read_input_data(Nag_Variogram *cov, Integer *np, double *params,
                            Nag_NormType *norm, double *var, double *xmin,
                            double *xmax, double *ymin, double *ymax,
                            Integer *ns, Integer *maxm, Nag_EmbedScale *corr,
                            Nag_EmbedPad *pad, Integer *s);
static void display_embedding_results(Integer approx, Integer *m, double rho,
                                      double *eig, Integer icount);
static void initialize_state(Integer *state);
static void display_realizations(Integer *ns, Integer s, double *xx,
                                 double *yy,
                                 double *z,
                                 Integer *exit_status);

int main(void)
{
  /*  Scalars */
  Integer        exit_status = 0;
  double         rho, var, xmax, xmin, ymax, ymin;
  Integer        approx, icount, np, s;
  /*  Arrays */
  double         eig[3], params[NPMAX];
  double         *lam = 0, *xx = 0, *yy = 0, *z = 0;
  Integer        m[2], maxm[2], ns[2], state[LENST];
  /* Nag types */
  Nag_NormType   norm;
  Nag_Variogram  cov;
  Nag_EmbedPad   pad;
  Nag_EmbedScale corr;
  NagError       fail;

  INIT_FAIL(fail);

  printf("nag_rand_field_2d_generate (g05zsc) Example Program Results\n\n");
  /* Get problem specifications from data file*/
  read_input_data(&cov, &np, params, &norm, &var, &xmin, &xmax, &ymin, &ymax,
                  ns, maxm, &corr, &pad, &s);
  if (!(lam = NAG_ALLOC(maxm[0]*maxm[1], double)) ||
      !(xx = NAG_ALLOC(ns[0], double)) ||
      !(yy = NAG_ALLOC(ns[1], double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  /* Get square roots of the eigenvalues of the embedding matrix.
   * nag_rand_field_2d_predef_setup (g05zrc).
   * Setup for simulating two-dimensional random fields, preset variogram,
   * circulant embedding method
   */
  nag_rand_field_2d_predef_setup(ns, xmin, xmax, ymin, ymax, maxm, var, cov,
                                 norm, np, params, pad, corr, lam, xx, yy, m,
                                 &approx, &rho, &icount, eig, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_rand_field_2d_predef_setup (g05zrc).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

  display_embedding_results(approx, m, rho, eig, icount);
  /* Initialize state array*/
  initialize_state(state);
  if (!(z = NAG_ALLOC(ns[0]*ns[1]*s, double)))
    {
      printf("Allocation failure\n");
      exit_status = -2;
      goto END;
    }
  /* Compute s random field realisations.
   * nag_rand_field_2d_generate (g05zsc).
   * Generates s realisations of a rwo-dimensional random field by the
   * circulant embedding method.
   */
  nag_rand_field_2d_generate(ns, s, m, lam, rho, state, z, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_rand_field_2d_generate (g05zsc).\n%s\n",
             fail.message);
      exit_status = 2;
      goto END;
    }
  display_realizations(ns, s, xx, yy, z, &exit_status);
 END:
  NAG_FREE(lam);
  NAG_FREE(xx);
  NAG_FREE(yy);
  NAG_FREE(z);
  return exit_status;
}

void read_input_data(Nag_Variogram *cov, Integer *np, double *params,
                     Nag_NormType *norm, double *var, double *xmin,
                     double *xmax, double *ymin, double *ymax, Integer *ns,
                     Integer *maxm, Nag_EmbedScale *corr,
                     Nag_EmbedPad *pad, Integer *s)
{
  Integer j;
  char    nag_enum_arg[40];

  /* Read in covariance function name and convert to value using
   * nag_enum_name_to_value (x04nac).
   */
  scanf("%*[^\n] %39s%*[^\n]", nag_enum_arg);
  *cov = (Nag_Variogram) nag_enum_name_to_value(nag_enum_arg);
  /* Read in parameters */
  scanf("%ld%*[^\n]", np);
  for (j = 0; j < *np; j++)
    scanf("%lf", &params[j]);
  scanf("%*[^\n]");
  /* Read in norm type by name and convert to value */
  scanf(" %39s%*[^\n]", nag_enum_arg);
  *norm = (Nag_NormType) nag_enum_name_to_value(nag_enum_arg);
  /* Read in variance of random field. */
  scanf("%lf%*[^\n]", var);
  /* Read in domain endpoints*/
  scanf("%lf %lf%*[^\n]", xmin, xmax);
  scanf("%lf %lf%*[^\n]", ymin, ymax);
  /* Read in number of sample points in each direction*/
  scanf("%ld %ld%*[^\n]", &ns[0], &ns[1]);
  /* Read in maximum size of embedding matrix*/
  scanf("%ld %ld%*[^\n]", &maxm[0], &maxm[1]);
  /* Read name of scaling in case of approximation and convert to value. */
  scanf(" %39s%*[^\n]", nag_enum_arg);
  *corr = (Nag_EmbedScale) nag_enum_name_to_value(nag_enum_arg);
  /* Read in choice of padding and convert name to value. */
  scanf(" %39s%*[^\n]", nag_enum_arg);
  *pad = (Nag_EmbedPad) nag_enum_name_to_value(nag_enum_arg);
  /* Read in number of realization samples to be generated*/
  scanf("%ld%*[^\n]", s);
}

void display_embedding_results(Integer approx, Integer *m, double rho,
                               double *eig, Integer icount)
{
  Integer j;
  /* Display size of embedding matrix*/
  printf("\nSize of embedding matrix = %ld\n\n", m[0]*m[1]);
  /* Display approximation information if approximation used*/
  if (approx == 1)
    {
      printf("Approximation required\n\n");
      printf("rho = %10.5f\n", rho);
      printf("eig = ");
      for (j = 0; j < 3; j++)
        printf("%10.5f ", eig[j]);
      printf("\nicount = %ld\n", icount);
    }
  else
    {
      printf("Approximation not required\n");
    }
}
void initialize_state(Integer *state)
{
  /*  Scalars */
  Integer  inseed = 14965, lseed = LSEED, subid = 1;
  Integer  lstate;
  /*  Arrays */
  Integer  seed[LSEED];
  /* Nag types */
  NagError fail;

  INIT_FAIL(fail);
  lstate = LENST;
  seed[0] = inseed;
  /* nag_rand_init_repeatable (g05kfc).
   * Initializes a pseudorandom number generator to give a repeatable sequence.
   */
  nag_rand_init_repeatable(Nag_Basic, subid, seed, lseed, state, &lstate,
                           &fail);
}

void display_realizations(Integer *ns, Integer s, double *xx, double *yy,
                          double *z, Integer *exit_status)
{
  /*  Scalars */
  Integer  indent = 0, ncols = 80;
  Integer  i, j, nn;
  /*  Arrays */
  char     **rlabs = 0;
  /* Nag types */
  NagError fail;

  INIT_FAIL(fail);

  nn = ns[0]*ns[1];
  if (!(rlabs = NAG_ALLOC(nn, char *)))
    {
      printf("Allocation failure\n");
      *exit_status = -3;
      goto END;
    }
  /* Set row labels to mesh points (column label is realization number).*/
  for (j = 0; j < ns[1]; j++)
    {
      for (i = 0; i < ns[0]; i++)
        {
          if (!(rlabs[j*ns[0]+i] = NAG_ALLOC(13, char)))
            {
              printf("Allocation failure\n");
              *exit_status = -4;
              goto END;
            }
          if (i == 0)
            {
              sprintf(rlabs[j*ns[0]+i], "%6.1f%6.1f", xx[i], yy[j]);
            }
          else
            {
              sprintf(rlabs[j*ns[0]+i], "%6.1f%6s", xx[i], ".");
            }
        }
    }
  printf("\n");
  fflush(stdout);
  /* Display random field results, z, using the comprehensive real general
   * matrix print routine nag_gen_real_mat_print_comp (x04cbc).
   */
  nag_gen_real_mat_print_comp(Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag,
                              nn, s, z, nn, "%10.5f", "Random field "
                              "realisations (x,y coordinates first):",
                              Nag_CharacterLabels, (const char **) rlabs,
                              Nag_IntegerLabels, NULL, ncols, indent, 0, &fail);
 END:
  for (i = 0; i < nn; i++)
    {
      NAG_FREE(rlabs[i]);
    }
  NAG_FREE(rlabs);
}