/* ------------------------------------------------------------------------- *
 *
 * Examples for using Algorithms HITRO-box and HITRO-plate
 *
 * Roman Karawatzki, Josef Leydold, and Klaus Poetzelberger (2006).
 * Automatic Markov Chain Monte Carlo Procedures for Sampling from 
 * Multivariate Distributions.
 *
 * (c) 2006, Josef Leydold
 *
 * ------------------------------------------------------------------------- *
 *
 * Sampling from Multinormal distribution with independent components,
 * mean (0,...,0), and variance 'i' for the 'i'-th component.
 *
 * ------------------------------------------------------------------------- *
 *
 * Usage:
 *   cc hitro.c example.c -o example -lm
 *
 * ------------------------------------------------------------------------- */

/* ------------------------------------------------------------------------- *
 *
 * WARNING: These algorithms are Markov chain sampling algorithm, i.e.
 * they produce *correlated* random points.
 * Thus it they might fail a goodness-of-fit test (e.g. the chisquare
 * goodness-of-fit test). This is in particular the case when "Gibbs
 * sampling" is selected. The autocorrelation can be reduced by
 * "thinning" the generated Markov chain, i.e., by using only each
 * 't'-th element. We made good experiences when 't' was choosen as the
 * number of variables.
 *
 * ------------------------------------------------------------------------- */

 
#define N (3)        /* dimension   */
#define S (100)      /* sample size */

/* ------------------------------------------------------------------------- */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "hitro.h"

/* ------------------------------------------------------------------------- *
 * [ uniform generator ]
 *
 * Combined multiple recursive generator.
 * The two components of the generator are
 *
 *   x_{1,k}  = (2^{22} x_{1,k-2} + (2^7 +1)x_{1,k-3}) mod (2^{31}-1)
 *   x_{2,k}  = (2^{15} x_{2,k-1} + (2^{15} +1)x_{2,k-3} mod (2^{31}-21069)
 *
 * x_{1,k} and x_{2,k} are combined together and the result is
 * multiplied by 1/(2^{31}-1) to have a number between 0 and 1.
 *
 * REFERENCE:
 * L'Ecuyer, P. and R. Touzin (2000): Fast Combined Multiple Recursive
 *    Generators with Multipliers of the Form a = ±2^q±2^r.
 *    in: J.A. Jones, R.R. Barton, K. Kang, and P.A. Fishwick (eds.),
 *    Proc. 2000 Winter Simulation Conference, 683-689.
 *
 * Copyright for uniform generator by Renee Touzin.
 *
 */

/* set seed (must not be 0!) */
#define SEED  (1245889L)

/* status variables */
static unsigned long x10 = SEED;
static unsigned long x11 = SEED;
static unsigned long x12 = SEED;
static unsigned long x20 = SEED;
static unsigned long x21 = SEED;
static unsigned long x22 = SEED;

static double uniform_MRG31k3p (void)
     /* Combined multiple recursive generator.                               */
     /* Copyright (c) 2002 Renee Touzin.                                     */
{

# define m1      2147483647
# define m2      2147462579
# define norm    4.656612873077393e-10
# define mask11  511
# define mask12  16777215
# define mask20  65535
 
  register unsigned long y1, y2;  /* For intermediate results */
  
  /* First component */
  y1 = ( (((x11 & mask11) << 22) + (x11 >> 9))
         + (((x12 & mask12) << 7)  + (x12 >> 24)) );
  if (y1 > m1) y1 -= m1;
  y1 += x12;
  if (y1 > m1) y1 -= m1;
  x12 = x11;  x11 = x10;  x10 = y1;
 
  /* Second component */
  y1 = ((x20 & mask20) << 15) + 21069 * (x20 >> 16);
  if (y1 > m2) y1 -= m2;
  y2 = ((x22 & mask20) << 15) + 21069 * (x22 >> 16);
  if (y2 > m2) y2 -= m2;
  y2 += x22;
  if (y2 > m2) y2 -= m2;
  y2 += y1;
  if (y2 > m2) y2 -= m2;
  x22 = x21;  x21 = x20;  x20 = y2;

  /* Combination */
  if (x10 <= x20)
    return ((x10 - x20 + m1) * norm);
  else 
    return ((x10 - x20) * norm);

} /* end of uniform_MRG31k3p() */
 
/*---------------------------------------------------------------------------*
 * [ PDF of a multinormal distribution ]
 *
 * mean(X_i)    = 0  for all i = 1, ..., n
 * Var(X_i)     = i  for all i = 1, ..., n
 * Cov(X_i,X_j) = 0  for all i!=j
 *
 */

static double pdf_normal (int dim, double *x) 
{
  int i;
  double f=0.0;
  
  for (i=0; i<dim; i++)
    f += x[i]*x[i] / (i+1.0);

  return exp(-f/2.0);
} /* end of pdf_normal() */

/*---------------------------------------------------------------------------*/
/* Auxiliary routines                                                        */
/*---------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*
 * [ Print sample ]
 *
 */

static void print_sample( HITRO chain )
{
  int i,k;
  double x[N];

  for (k=0; k<S; k++) {
    hitro_sample( chain, x );
    printf("%13.8f", x[0]);
    for (i=1; i<N; i++)
      printf(",%13.8f", x[i]);
    printf("\n");
  }
  
} /* end of print_sample() */

/*---------------------------------------------------------------------------*
 * [ Check chain ]
 *
 * There are a few tests that check the consistency of the given data.
 * If an error is detected (e.g. dimension < 1) then NULL is returns by
 * the constructor calls.
 */

static void check_chain( HITRO chain )
{
  if (chain == NULL) {
    fprintf(stderr, "ERROR: cannot create HITRO chain\n");
    exit (EXIT_FAILURE);
  }
} /* end of check_chain() */

/*---------------------------------------------------------------------------*/
/* Main                                                                      */
/*---------------------------------------------------------------------------*/

int main() 
{
  HITRO chain;
  int i;
  double x[N];
  double center[N], start[N];
  double llbr[N+1], rubr[N+1];

  /* --- HITRO-plate: simple example --------------------------------------- */
  chain = hitro_plate_create 
    (N,                /* dimension */
     NULL,             /* we use (0,...,0) as center */
     NULL,             /* we use the center as starting point of chain, i.e. (0,...,0) */
     0.,               /* we do not compute the maximum value of the PDF */
     pdf_normal,       /* this is the pointer to the PDF of the multinormal distribution */
     uniform_MRG31k3p, /* this is the pointer to the uniform RNG */
     1.0               /* since our PDF is log concave we can use 'r'=1 */
     );

  /* we should test the result of the constructor. if it fails we get a
     NULL pointer which must not be used. */
  if (chain == NULL) {
    fprintf(stderr, "ERROR: cannot create HITRO chain\n");
    exit (EXIT_FAILURE);
  }

  /* now print a sample */
  printf("[0] HITRO-plate: simple example\n");
  hitro_sample( chain, x );
  for (i=0; i<N; i++)
    printf("%13.8f\t", x[i]);
  printf("\n");
  
  /* at last we delete the HITRO chain */
  hitro_delete(chain);

  /* REMARK: in the following example we use the auxiliary routines 'check_chain' 
     and 'print_sample' to save code lines. */
  chain = hitro_plate_create(N,NULL,NULL,0.,pdf_normal,uniform_MRG31k3p,1.0);
  check_chain( chain );
  printf("[1] HITRO-plate: simple example\n");
  print_sample(chain);
  hitro_delete(chain);

  /* --- HITRO-plate: use upper bound for PDF ------------------------------ */
  /* it is a good idea to provide an upper bound of the PDF (e.g. its value
   * at the mode). This increases the performance of the algorithm.
   */
  chain = hitro_plate_create 
    (N,                /* dimension */
     NULL,             /* we use (0,...,0) as center */
     NULL,             /* we use the center as starting point of chain, i.e. (0,...,0) */
     1.,               /* use PDF(mode) as upper bound for PDF */
     pdf_normal,       /* this is the pointer to the PDF of the multinormal distribution */
     uniform_MRG31k3p, /* this is the pointer to the uniform RNG */
     1.0               /* since our PDF is log concave we can use 'r'=1 */
     );
  check_chain( chain );
  printf("[2] HITRO-plate: use upper bound for PDF\n");
  print_sample(chain);
  hitro_delete(chain);

  /* --- HITRO-plate: 'center' and starting point -------------------------- */
  for (i=0; i<N; i++) center[i] = 1.0/(i+1.0);
  for (i=0; i<N; i++) start[i] = -1.0/(i+1.0);
  chain = hitro_plate_create 
    (N,                /* dimension */
     center,           /* provide 'center' of distribution */
     start,            /* provide starting point of chain */
     1.,               /* use PDF(mode) as upper bound for PDF */
     pdf_normal,       /* this is the pointer to the PDF of the multinormal distribution */
     uniform_MRG31k3p, /* this is the pointer to the uniform RNG */
     1.0               /* since our PDF is log concave we can use 'r'=1 */
     );
  check_chain( chain );
  printf("[3] HITRO-plate: use center and starting point\n");
  print_sample(chain);
  hitro_delete(chain);
  
  /* --- HITRO-plate: use larger value for 'r' ----------------------------- */
  /* sometimes it can be necessary to use values larger than 1 for 'r' to 
   * ensure that the transformed region is bounded.
   * (This is not the case in our example distribution!)
   */
  chain = hitro_plate_create 
    (N,                /* dimension */
     NULL,             /* we use (0,...,0) as center */
     NULL,             /* we use the center as starting point of chain, i.e. (0,...,0) */
     1.,               /* use PDF(mode) as upper bound for PDF */
     pdf_normal,       /* this is the pointer to the PDF of the multinormal distribution */
     uniform_MRG31k3p, /* this is the pointer to the uniform RNG */
     2.0               /* use 'r'=2 */
     );
  check_chain( chain );
  printf("[4] HITRO-plate: use r = 2\n");
  print_sample(chain);
  hitro_delete(chain);

  /* --- HITRO-box: simple example ----------------------------------------- */
  chain = hitro_box_create 
    (N,                /* dimension */
     NULL,             /* we use (0,...,0) as center */
     NULL,             /* we use the center as starting point of chain, i.e. (0,...,0) */
     NULL,             /* we do not know the left lower vertex of the bounding rectangle */
     NULL,             /* we do not know the right upper vertex of the bounding rectangle */
     pdf_normal,       /* this is the pointer to the PDF of the multinormal distribution */
     uniform_MRG31k3p, /* this is the pointer to the uniform RNG */
     1.0,              /* since our PDF is log concave we can use 'r'=1 */
     1,                /* use adaptive sampling from line segments */
     1                 /* use coordinate direction ("Gibbs") sampling */
     );
  check_chain( chain );
  printf("[5] HITRO-box: simple example\n");
  print_sample(chain);
  hitro_delete(chain);

  /* --- HITRO-box: use bounding rectangle --------------------------------- */
  /* it is a good idea to provide the vertices of the bounding rectangle.
   * This increases the performance of the algorithm.
   */
  llbr[0] = 0; rubr[0] = 1;
  for (i=1; i<=N; i++) llbr[i] = -pow(exp(-0.5-N/2.0), 1.0/(N+1.0)) * sqrt(i*(N+1.0));
  for (i=1; i<=N; i++) rubr[i] =  pow(exp(-0.5-N/2.0), 1.0/(N+1.0)) * sqrt(i*(N+1.0));
  chain = hitro_box_create 
    (N,                /* dimension */
     NULL,             /* we use (0,...,0) as center */
     NULL,             /* we use the center as starting point of chain, i.e. (0,...,0) */
     llbr,             /* left lower vertex of bounding rectangle */
     rubr,             /* right upper vertex of bounding rectangle */
     pdf_normal,       /* this is the pointer to the PDF of the multinormal distribution */
     uniform_MRG31k3p, /* this is the pointer to the uniform RNG */
     1.0,              /* since our PDF is log concave we can use 'r'=1 */
     1,                /* use adaptive sampling from line segments */
     1                 /* use coordinate direction ("Gibbs") sampling */
     );
  check_chain( chain );
  printf("[6] HITRO-box: use bounding rectangle\n");
  print_sample(chain);
  hitro_delete(chain);

  /* --- HITRO-box: simple rejection --------------------------------------- */
  /* sometimes it might be necessary to use simple rejection from the line
   * segments (e.g. because the PDF is multimodal and adaptive rejection
   * converges too slowly).
   */
  chain = hitro_box_create 
    (N,                /* dimension */
     NULL,             /* we use (0,...,0) as center */
     NULL,             /* we use the center as starting point of chain, i.e. (0,...,0) */
     llbr,             /* left lower vertex of bounding rectangle */
     rubr,             /* right upper vertex of bounding rectangle */
     pdf_normal,       /* this is the pointer to the PDF of the multinormal distribution */
     uniform_MRG31k3p, /* this is the pointer to the uniform RNG */
     1.0,              /* since our PDF is log concave we can use 'r'=1 */
     0,                /* use simple rejection from line segments */
     1                 /* use coordinate direction ("Gibbs") sampling */
     );
  check_chain( chain );
  printf("[7] HITRO-box: use simple rejection\n");
  print_sample(chain);
  hitro_delete(chain);

  /* --- HITRO-box: use random direction sampling -------------------------- */
  /* One can use random direction sampling instead of Gibbs sampling.
   */
  llbr[0] = 0; rubr[0] = 1;
  chain = hitro_box_create 
    (N,                /* dimension */
     NULL,             /* we use (0,...,0) as center */
     NULL,             /* we use the center as starting point of chain, i.e. (0,...,0) */
     llbr,             /* left lower vertex of bounding rectangle */
     rubr,             /* right upper vertex of bounding rectangle */
     pdf_normal,       /* this is the pointer to the PDF of the multinormal distribution */
     uniform_MRG31k3p, /* this is the pointer to the uniform RNG */
     1.0,              /* since our PDF is log concave we can use 'r'=1 */
     1,                /* use adaptive sampling from line segments */
     0                 /* use random direction sampling */
     );
  check_chain( chain );
  printf("[8] HITRO-box: use random direction sampling\n");
  print_sample(chain);
  hitro_delete(chain);

  /* --- end --------------------------------------------------------------- */
  
  exit (EXIT_SUCCESS);
}

/*---------------------------------------------------------------------------*/