/* nag_rand_bb (g05xbc) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 */
#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg05.h>
#include <nagf07.h>

int get_z(Nag_OrderType order, Integer ntimes, Integer d, Integer a,
          Integer npaths, double *z, Integer pdz);
void display_results(Nag_OrderType order, Integer npaths, Integer ntimes,
                     Integer d, double *b, Integer pdb);

#define CHECK_FAIL(name,fail) if(fail.code != NE_NOERROR) { \
    printf("Error from %s.\n%s\n",name,fail.message); \
    exit_status = -1; goto END; }

int main(void)
{
  /*  Scalars */
  Integer exit_status = 0;
  double t0, tend;
  Integer a, d, pdb, pdc, pdz, nmove, npaths, ntimes, i;
  /*  Arrays */
  double *b = 0, *c = 0, *intime = 0, *rcomm = 0, *start = 0, *term = 0,
         *times = 0, *z = 0;
  Integer *move = 0;
  /* Nag Types */
  NagError fail;
  Nag_OrderType order;

  INIT_FAIL(fail);

  /* Parameters which determine the bridge. */
  ntimes = 10;
  t0 = 0.0;
  npaths = 2;
  /* Create a non-free bridge. */
  a = 1;
  nmove = 0;
  d = 3;
#ifdef NAG_COLUMN_MAJOR
  order = Nag_ColMajor;
  pdz = npaths;
  pdb = npaths;
#else
  order = Nag_RowMajor;
  pdz = d * (ntimes + 1 - a);
  pdb = d * (ntimes + 1);
#endif
  pdc = d;
#define C(I,J) c[(J-1)*pdc+I-1]

  /* Allocate memory */
  if (!(intime = NAG_ALLOC((ntimes), double)) ||
      !(times = NAG_ALLOC((ntimes), double)) ||
      !(rcomm = NAG_ALLOC((12 * (ntimes + 1)), double)) ||
      !(start = NAG_ALLOC(d, double)) ||
      !(term = NAG_ALLOC(d, double)) ||
      !(c = NAG_ALLOC(d * pdc, double)) ||
      !(z = NAG_ALLOC(d * (ntimes + 1 - a) * npaths, double)) ||
      !(b = NAG_ALLOC(d * (ntimes + 1) * npaths, double)) ||
      !(move = NAG_ALLOC(nmove, Integer))
         )
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Fix the time points at which the bridge is required */
  for (i = 0; i < ntimes; i++)
    intime[i] = t0 + (double) (i + 1);
  tend = t0 + (double) (ntimes + 1);

  /* Create a Brownian bridge construction order out of a set of input times
   * using nag_rand_bb_make_bridge_order (g05xec).
   */
  nag_rand_bb_make_bridge_order(Nag_RLRoundDown, t0, tend, ntimes, intime,
                                nmove, move, times, &fail);
  CHECK_FAIL("nag_rand_bb_make_bridge_order", fail);

  /* Initialize the Brownian bridge generator using 
   * nag_rand_bb_init (g05xac).
   */
  nag_rand_bb_init(t0, tend, times, ntimes, rcomm, &fail);
  CHECK_FAIL("nag_rand_bb_init (g05xac)", fail);

  /* We want the following covariance matrix ... */
  C(1, 1) = 6.0;
  C(2, 1) = C(1, 2) = 1.0;
  C(3, 1) = C(1, 3) = -0.2;
  C(2, 2) = 5.0;
  C(3, 2) = C(2, 3) = 0.3;
  C(3, 3) = 4.0;
  /* Cholesky factorize the covariance matrix C, as required by
   * nag_rand_bb (g05xbc), using nag_dpotrf (f07fdc).
   */
  nag_dpotrf(Nag_ColMajor, Nag_Lower, d, c, pdc, &fail);
  CHECK_FAIL("nag_dpotrf", fail);

  /* Generate the random numbers z. */
  if (get_z(order, ntimes, d, a, npaths, z, pdz) != 0) {
    printf("Error generating random numbers\n");
    exit_status = -1;
    goto END;
  }
  /* Give start and terminal values of pinned bridge */
  start[0] = start[1] = start[2] = 0.0;
  term[0] = 1.0;
  term[1] = 0.5;
  term[2] = 0.0;

  /* Generate paths for a free or non-free Wiener process using the
   * Brownian bridge algorithm: nag_rand_bb (g05xbc).
   */
  nag_rand_bb(order, npaths, d, start, a, term, z, pdz, c, pdc,
              b, pdb, rcomm, &fail);
  CHECK_FAIL("nag_rand_bb", fail);

  /* Display the results */
  display_results(order, npaths, ntimes, d, b, pdb);
END:
  NAG_FREE(b);
  NAG_FREE(c);
  NAG_FREE(intime);
  NAG_FREE(rcomm);
  NAG_FREE(start);
  NAG_FREE(term);
  NAG_FREE(times);
  NAG_FREE(z);
  NAG_FREE(move);
  return exit_status;
}

int get_z(Nag_OrderType order, Integer ntimes, Integer d, Integer a,
          Integer npaths, double *z, Integer pdz)
{
  /* Scalars */
  Integer exit_status = 0;
  Integer lseed, lstate, idim, liref, i;
  /* Arrays */
  Integer seed[1], *iref = 0, state[80];
  double *xmean = 0, *stdev = 0;
  /* Nag Types */
  NagError fail;

  INIT_FAIL(fail);

  lstate = 80;
  lseed = 1;
  idim = d * (ntimes + 1 - a);
  liref = 32 * idim + 7;
  if (!(iref = NAG_ALLOC((liref), Integer)) ||
      !(xmean = NAG_ALLOC((idim), double)) ||
      !(stdev = NAG_ALLOC((idim), double)))
  {
    printf("Allocation failure in get_z\n");
    exit_status = -1;
    goto END;
  }

  /* We now need to generate the input pseudorandom numbers. */
  seed[0] = 1023401;
  /* Initialize a pseudorandom number generator to give a repeatable sequence
   * using nag_rand_init_repeatable (g05kfc).
   */
  nag_rand_init_repeatable(Nag_MRG32k3a, 0, seed, lseed, state, &lstate,
                           &fail);
  CHECK_FAIL("nag_rand_init_repeatable (g05kfc)", fail);

  /* Initialize a scrambled quasi-random number generator using
   * nag_quasi_init_scrambled (g05ync).
   */
  nag_quasi_init_scrambled(Nag_QuasiRandom_Sobol, Nag_FaureTezuka, idim,
                           iref, liref, 0, 32, state, &fail);
  CHECK_FAIL("nag_quasi_init_scrambled (g05ync)", fail);

  for (i = 0; i < idim; i++) {
    xmean[i] = 0.0;
    stdev[i] = 1.0;
  }
  /* Generate a (repeatable) Normal quasi-random number sequence using
   * nag_quasi_rand_normal (g05yjc). 
   */
  nag_quasi_rand_normal(order, xmean, stdev, npaths, z, pdz, iref, &fail);
  CHECK_FAIL("nag_quasi_rand_normal (g05yjc)", fail);

END:
  NAG_FREE(iref);
  NAG_FREE(xmean);
  NAG_FREE(stdev);
  return exit_status;
}

void display_results(Nag_OrderType order, Integer npaths, Integer ntimes,
                     Integer d, double *b, Integer pdb)
{
#define B(I,J) (order==Nag_RowMajor ? b[(I-1)*pdb + J-1]:b[(J-1)*pdb + I-1])
  Integer i, p, k;

  printf("nag_rand_bb (g05xbc) Example Program Results\n\n");
  for (p = 1; p <= npaths; p++) {
    printf("Wiener Path  %1" NAG_IFMT ",  %1" NAG_IFMT "", p, ntimes + 1);
    printf(" time steps,  %1" NAG_IFMT " dimensions\n", d);

    for (k = 1; k <= d; k++)
      printf("%10" NAG_IFMT " ", k);
    printf("\n");

    for (i = 0; i < ntimes + 1; i++) {
      printf("%2" NAG_IFMT " ", i + 1);
      for (k = 1; k <= d; k++)
        printf("%10.4f", B(p, k + i * d));
      printf("\n");
    }
    printf("\n");
  }
}