/* nag_rand_bb_inc_init (g05xcc) 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(Integer nelements, double *z);
void display_results(Nag_OrderType order, Integer npaths, double t0,
                     double tend, Integer ntimes, double intime[],
                     Integer d, double start[], double bb[],
                     double bd[], Integer pdb);

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

int main(void)
{
#define C(I,J) c[(J-1)*d + I-1]
  Integer exit_status = 0;
  NagError fail;
  /*  Scalars */
  double t0, tend;
  Integer a, d, pdb, pdc, pdz, nmove, npaths, ntimes, i;
  /*  Arrays */
  double *bb = 0, *c = 0, *intime = 0, *rcommb = 0, *start = 0,
         *term = 0, *times = 0, *zb = 0, *rcommd = 0, *zd = 0,
         *diff = 0, *bd = 0;
  Integer *move = 0;
  INIT_FAIL(fail);

  /* Parameters which determine the bridge */
  ntimes = 10;
  t0 = 0.0;
  npaths = 2;
  a = 0;
  nmove = 0;
  d = 2;
  pdz = npaths;
  pdb = npaths;
  pdc = d;
  /* Allocate memory */
  if (!(intime = NAG_ALLOC((ntimes), double)) ||
      !(times = NAG_ALLOC((ntimes), double)) ||
      !(rcommb = NAG_ALLOC((12 * (ntimes + 1)), double)) ||
      !(rcommd = NAG_ALLOC((12 * (ntimes + 1)), double)) ||
      !(start = NAG_ALLOC(d, double)) ||
      !(term = NAG_ALLOC(d, double)) ||
      !(diff = NAG_ALLOC(d, double)) ||
      !(c = NAG_ALLOC(pdc * d, double)) ||
      !(zb = NAG_ALLOC(pdz * d * (ntimes + 1 - a), double)) ||
      !(zd = NAG_ALLOC(pdz * d * (ntimes + 1 - a), double)) ||
      !(bb = NAG_ALLOC(pdb * d * (ntimes + 1), double)) ||
      !(bd = NAG_ALLOC(pdb * d * (ntimes + 1), 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);

  /* g05xec. Creates a Brownian bridge construction order out of a set of */
  /* input times */
  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);

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

  start[0] = 0.0, start[1] = 2.0;
  /* We want the following covariance matrix */
  C(1, 1) = 6.0;
  C(2, 1) = -1.0;
  C(1, 2) = -1.0;
  C(2, 2) = 5.0;
  /* nag_rand_bb uses Cholesky factorization of the covariance matrix C */

  /* f07fdc. Cholesky factorization of real positive definite matrix */
  nag_dpotrf(Nag_ColMajor, Nag_Lower, d, c, d, &fail);
  CHECK_FAIL("nag_dpotrf", fail);

  /* Generate the random numbers */
  if (get_z(npaths * d * (ntimes + 1 - a), zb) != 0) {
    printf("Error generating random numbers\n");
    exit_status = -1;
    goto END;
  }

  /* Copy the random numbers for call to g05xdc */
  for (i = 0; i < npaths * d * (ntimes + 1 - a); i++)
    zd[i] = zb[i];

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

  /* nag_rand_bb_inc_init (g05xcc). Initializes the generator which backs out 
   * the increments of sample paths generated by a Brownian bridge algorithm */
  nag_rand_bb_inc_init(t0, tend, times, ntimes, rcommd, &fail);
  CHECK_FAIL("nag_rand_bb_inc_init", fail);

  /* g05xdc. Backs out the increments from sample paths generated by a */
  /* Brownian bridge algorithm */
  nag_rand_bb_inc(Nag_ColMajor, npaths, d, a, diff, zd, pdz, c, pdc, bd, pdb,
                  rcommd, &fail);
  CHECK_FAIL("nag_rand_bb_inc", fail);

  /* Display the results */
  display_results(Nag_ColMajor, npaths, t0, tend,
                  ntimes, intime, d, start, bb, bd, pdb);
END:
  ;
  NAG_FREE(bb);
  NAG_FREE(bd);
  NAG_FREE(c);
  NAG_FREE(intime);
  NAG_FREE(rcommb);
  NAG_FREE(rcommd);
  NAG_FREE(start);
  NAG_FREE(term);
  NAG_FREE(diff);
  NAG_FREE(times);
  NAG_FREE(zb);
  NAG_FREE(zd);
  NAG_FREE(move);
  return exit_status;
#undef C
}

int get_z(Integer nelements, double *z)
{
  NagError fail;
  Integer lseed, lstate, exit_status = 0;
  /*  Arrays */
  Integer seed[1];
  Integer state[80];
  lstate = 80;
  lseed = 1;
  INIT_FAIL(fail);

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

  /* nag_rand_normal.  Generates a vector of pseudorandom numbers from */
  /* a Normal distribution */
  nag_rand_normal(nelements, 0.0, 1.0, state, z, &fail);
  CHECK_FAIL("nag_rand_normal", fail);

END:;
  return exit_status;
}

void display_results(Nag_OrderType order, Integer npaths, double t0,
                     double tend, Integer ntimes, double intime[],
                     Integer d, double start[], double bb[],
                     double bd[], Integer pdb)
{
// Order consistend with Nag_RowMajor
#define BB(I,J) bb[(I-1)*pdb + J-1]
#define BD(I,J) bd[(I-1)*pdb + J-1]
#define CUM(I)  cum[I-1]

  Integer i, p, k;
  double *cum = 0;
  if (!(cum = NAG_ALLOC(d, double)))
  {
    printf("Error allocating memory in display_results\n");
    return;
  }

  printf("nag_rand_bb_inc_init (g05xcc) Example Program Results\n\n");
  for (p = 1; p <= npaths; p++) {
    printf("Weiner Path  ");
    printf("%3" NAG_IFMT " ", p);
    printf(",  ");
    printf("%3" NAG_IFMT " ", ntimes + 1);
    printf(" time steps,  ");
    printf("%3" NAG_IFMT " ", d);
    printf(" dimensions \n");
    printf("      Output of  g05xbc    Output of g05xdc    Sum of g05xdc \n");

    for (k = 1; k <= d; k++)
      CUM(k) = start[k - 1];
    /* Print first point */
    i = 0;
    printf("%2" NAG_IFMT " ", i + 1);
    if (order == Nag_RowMajor) {
      for (k = 1; k <= d; k++)
        CUM(k) += BD(p, k) * (intime[i] - t0);
      for (k = 1; k <= d; k++)
        printf("  %8.4f", BB(p, k));
      for (k = 1; k <= d; k++)
        printf("  %8.4f", BD(p, k));
      for (k = 1; k <= d; k++)
        printf("  %8.4f", CUM(k));
    }
    else {
      for (k = 1; k <= d; k++)
        CUM(k) += BD(k, p) * (intime[i] - t0);
      for (k = 1; k <= d; k++)
        printf("  %8.4f", BB(k, p));
      for (k = 1; k <= d; k++)
        printf("  %8.4f", BD(k, p));
      for (k = 1; k <= d; k++)
        printf("  %8.4f", CUM(k));
    }
    printf("\n");

    /* Print intermediate points */
    for (i = 1; i < ntimes; i++) {
      printf("%2" NAG_IFMT " ", i + 1);
      if (order == Nag_RowMajor) {
        for (k = 1; k <= d; k++)
          CUM(k) += BD(p, i * d + k) * (intime[i] - intime[i - 1]);
        for (k = 1; k <= d; k++)
          printf("  %8.4f", BB(p, i * d + k));
        for (k = 1; k <= d; k++)
          printf("  %8.4f", BD(p, i * d + k));
        for (k = 1; k <= d; k++)
          printf("  %8.4f", CUM(k));
      }
      else {
        for (k = 1; k <= d; k++)
          CUM(k) += BD(i * d + k, p) * (intime[i] - intime[i - 1]);
        for (k = 1; k <= d; k++)
          printf("  %8.4f", BB(i * d + k, p));
        for (k = 1; k <= d; k++)
          printf("  %8.4f", BD(i * d + k, p));
        for (k = 1; k <= d; k++)
          printf("  %8.4f", CUM(k));
      }
      printf("\n");
    }

    /* Print final point */
    i = ntimes;
    printf("%2" NAG_IFMT " ", i + 1);
    if (order == Nag_RowMajor) {
      for (k = 1; k <= d; k++)
        CUM(k) += BD(p, i * d + k) * (tend - intime[i - 1]);
      for (k = 1; k <= d; k++)
        printf("  %8.4f", BB(p, i * d + k));
      for (k = 1; k <= d; k++)
        printf("  %8.4f", BD(p, i * d + k));
      for (k = 1; k <= d; k++)
        printf("  %8.4f", CUM(k));
    }
    else {
      for (k = 1; k <= d; k++)
        CUM(k) += BD(i * d + k, p) * (tend - intime[i - 1]);
      for (k = 1; k <= d; k++)
        printf("  %8.4f", BB(i * d + k, p));
      for (k = 1; k <= d; k++)
        printf("  %8.4f", BD(i * d + k, p));
      for (k = 1; k <= d; k++)
        printf("  %8.4f", CUM(k));
    }
    printf("\n\n");
  }

  NAG_FREE(cum);
}