/* nag_rand_bb_init (g05xac) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 24, 2013.
 */
#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, Integer ntimes, 
                     Integer d, double *b, 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)*pdc + I-1]
  Integer exit_status = 0;
  NagError fail;
  /*  Scalars */
  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;
  INIT_FAIL(fail);

  /* Parameters which determine the bridge */
  ntimes = 10;
  t0 = 0.0;
  npaths = 2;
  a = 0;
  nmove = 0;
  d = 3;
  pdz = d*(ntimes+1-a);
  pdb = d*(ntimes+1);
  pdc = d;
  /* 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(pdc * d, double)) ||
          !( z = NAG_ALLOC(pdz * npaths, double)) || 
          !( b = NAG_ALLOC(pdb * 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);

  /* 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);

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

  /* We want the following covariance matrix*/
  C( 1,1 ) = 6.0;
  C( 2,1 ) = 1.0;
  C( 3,1 ) = -0.2;
  C( 1,2 ) = 1.0;
  C( 2,2 ) = 5.0;
  C( 3,2 ) = 0.3;
  C( 1,3 ) = -0.2;
  C( 2,3 ) = 0.3;
  C( 3,3 ) = 4.0;
  /* nag_rand_bb uses the Cholesky factorization of the covariance matrix C */
  /* f07fdc. Cholesky factorization of real positive definite matrix */ 
  nag_dpotrf(Nag_ColMajor, Nag_Lower, d, c, pdc, &fail);
  CHECK_FAIL("nag_dpotrf",fail);

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

  for(i=0; i<d; i++) start[i] = 0.0;
  /* g05xbc. Generate paths for a free or non-free Wiener process using the */
  /* Brownian bridge algorithm   */ 
  nag_rand_bb(Nag_RowMajor, npaths, d, start, a, term, z, pdz, c, pdc, b, pdb,
              rcomm, &fail);
  CHECK_FAIL("nag_rand_bb",fail);

  /* Display the results*/
  display_results(Nag_RowMajor, 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(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);

  /* g05skc.  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, 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_init (g05xac) Example Program Results\n\n");
  for ( p=1; p<=npaths; p++)
    {
      printf("Wiener Path ");
      printf("%1ld ", p);
      printf(",  ");
      printf("%1ld ",  ntimes + 1);
      printf(" time steps, ");
      printf("%1ld ",  d);
      printf(" dimensions \n");

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

    
      for (i=1; i<= ntimes+1; i++)
        {
              printf("%2ld ", i);
              for (k=1; k<=d; k++)
                  {
                      printf("%10.4f", B(p, k + (i-1)*d));

                  }
              printf("\n");

        }
      printf("\n");
    }
}