/* nag_rand_field_fracbm_generate (g05ztc) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 */

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

#define NPMAX 4
#define LENST 17
#define LSEED 1
static void read_input_data(double *params, double *xmax, 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 *yy, double *z,
                                 Integer *exit_status);

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

  INIT_FAIL(fail);

  printf("nag_rand_field_fracbm_generate (g05ztc) Example Program Results\n\n");
  /* Get problem specifications from data file */
  read_input_data(params, &xmax, &ns, &maxm, &corr, &pad, &s);
  if (!(lam = NAG_ALLOC(maxm, double)) || !(xx = NAG_ALLOC(ns, double)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  /* Get square roots of the eigenvalues of the embedding matrix.
   * nag_rand_field_1d_predef_setup (g05znc).
   * Setup for simulating one-dimensional random fields, preset variogram,
   * circulant embedding method
   */
  nag_rand_field_1d_predef_setup(ns, xmin, xmax, maxm, var, cov, np,
                                 params, pad, corr, lam, xx, &m, &approx,
                                 &rho, &icount, eig, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_field_1d_predef_setup (g05znc).\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 (!(yy = NAG_ALLOC(ns + 1, double)) ||
      !(z = NAG_ALLOC((ns + 1) * s, double)))
  {
    printf("Allocation failure\n");
    exit_status = -2;
    goto END;
  }
  /* Compute s random field realizations.
   * nag_rand_field_fracbm_generate (g05ztc).
   * Generates s realizations of a one-dimensional random field by the
   * circulant embedding method.
   */
  h = params[0];
  nag_rand_field_fracbm_generate(ns, s, m, xmax, h, lam, rho, state, z, yy,
                                 &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_field_fracbm_generate (g05ztc).\n%s\n",
           fail.message);
    exit_status = 2;
    goto END;
  }
  display_realizations(ns, s, yy, z, &exit_status);
END:
  NAG_FREE(lam);
  NAG_FREE(xx);
  NAG_FREE(yy);
  NAG_FREE(z);
  return exit_status;
}

void read_input_data(double *params, double *xmax, Integer *ns,
                     Integer *maxm, Nag_EmbedScale *corr,
                     Nag_EmbedPad *pad, Integer *s)
{
  char nag_enum_arg[40];
  scanf("%*[^\n]");
  /* Read in Hurst parameter H. */
  scanf("%lf%*[^\n]", &params[0]);
  /* Read in domain endpoint. */
  scanf("%lf%*[^\n]", xmax);
  /* Read in number of sample points. */
  scanf("%" NAG_IFMT "%*[^\n]", ns);
  params[1] = *xmax / ((double) (*ns));
  /* Read in maximum size of embedding matrix. */
  scanf("%" NAG_IFMT "%*[^\n]", maxm);
  /* 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("%" NAG_IFMT "%*[^\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 = %" NAG_IFMT "\n\n", m);
  /* 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 = %" NAG_IFMT "\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 *yy, double *z,
                          Integer *exit_status)
{
  /*  Scalars */
  Integer indent = 0, ncols = 80;
  Integer i, nsp1;
  /*  Arrays */
  char **rlabs = 0;
  /* Nag types */
  NagError fail;

  INIT_FAIL(fail);

  nsp1 = ns + 1;
  if (!(rlabs = NAG_ALLOC(nsp1, char *)))
  {
    printf("Allocation failure\n");
    *exit_status = -3;
    goto END;
  }
  /* Set row labels to mesh points (column label is realization number). */
  for (i = 0; i < nsp1; i++) {
    if (!(rlabs[i] = NAG_ALLOC(7, char)))
    {
      printf("Allocation failure\n");
      *exit_status = -4;
      goto END;
    }
    sprintf(rlabs[i], "%6.1f", yy[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, nsp1, s, z, nsp1, "%10.5f",
                              "Fractional Brownian"
                              " motion realizations (x coordinate first):",
                              Nag_CharacterLabels, (const char **) rlabs,
                              Nag_IntegerLabels, NULL, ncols, indent, 0,
                              &fail);
END:
  for (i = 0; i < nsp1; i++) {
    NAG_FREE(rlabs[i]);
  }
  NAG_FREE(rlabs);
}