/********************************************************************/
/**                                                                **/
/**   exapmle_5.c                                                  **/
/**                                                                **/
/**   A universal generator for multivariate T-concave             **/
/**   distributions                                                **/
/**                                                                **/
/**   generate multivariate random points                          **/
/**                                                                **/
/**   domain:  IR^n or rectangle                                   **/
/**   mode:    (m_1,...m_N)    (see below)                         **/
/**                                                                **/
/**   densities:                                                   **/
/**   (1)  normal         exp( -\sum_{i=1}^N a_i x_i^2             **/
/**   (2)  "exponential"  exp( -\sum_{i=1}^N a_i |x_i|             **/
/**   (3)  gugelhupf      max( 0, \prod_{i=1}^N (1-a_i*x_i^2) )    **/
/**   (4)  gugelhupf      max( 0, \prod_{i=1}^N (1-|x_i|^{a_i})    **/
/**   (5)  thimple        max( 0, \sum_{i=1}^N (1-a_i*x_i^2) )     **/
/**                                                                **/
/**   all a_i > 0                                                  **/
/**                                                                **/
/**   authors: Josef Leydold                                       **/
/**   last modified: Mon Mar  9 13:38:00 MET 1998                  **/
/**                                                                **/
/********************************************************************/              
/**                                                                **/
/**   Remark:                                                      **/
/**   get_mode does not work when started in the plain region      **/
/**   of densities (3), (4) and (5) !!                             **/
/**   (Sorry no error message)                                     **/
/**                                                                **/
/**   The algorithm does work if the center for the construction   **/
/**   of the cones is in this plain region !!                      **/
/**   (Sorry no error message)                                     **/
/**                                                                **/
/********************************************************************/              
/* the header file */
#include "tdrmv.h"

/*------------------------------------------------------------------*/
/* Compiling options in config.h                                    */
/*                                                                  */
/*    N          2, ..., 14                                         */
/*    MODE       1                                                  */
/*    RECTANGLE  0 or 1                                             */
/*    GAMMA      1 or 2  (is set to one if RECTANGLE == 1)          */
/*                                                                  */
/*------------------------------------------------------------------*/

#define SAMPLE_SIZE   1000      /* sample size                      */

#define DENSITY  1   /* select a density                            */

#define ROTATION 1   /* 1 ... Rotate coordinate system              */
                     /* 0 ... do not rotate                         */

/* mode of unrestricted density function */
/* static const double mode[] = {0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1}; */
static const double mode[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};

/* bounds of rectangle */
static const double lower[] = {0.1,0.1,0.1,0.1,0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1};
static const double upper[] = {0.9,0.9,0.9,0.9,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.};


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

/* some parameters for setting coefficients */
#define LOWESTCOEFF 2.     /* lowest value for coefficients */
#define RATIOCOEFF  1.5    /* bound for ratio between highest and lowest coefficients */
                           /* (all coefficients are are sampled from uniform distribution */
                           /* in intervall [ LOWESTCOEFF, LOWESTCOEFF * RATIOCOEFF] */

/* initialize uniform random number generator */
#define INIR1       2
#define INIR2       100

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

/* print random point */
static void print_rpoint(double *v, int n);

/* make rotation of coordinate system */
static void rotate(double *coord, double *rotcoord);  /* rotates and moves point */
static void invrotate(double *coord);                 /* rotates point back into original coordinate system */
static void initialize_rotation(void);                /* set angles for rotation randomly */

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

/*******************************************************************

  Main

  ******************************************************************/

int main()
{
  C_ROOT *hat;         /* pointer to hat */
  C_CONSTR *hc;        /* parameters for construting hat functions */
  double rpoint[N];    /* array for random point */
  int i;

  /* density */
#if DENSITY == 1
  double initialize_normal();      /* initialize density function */
  double normal();                 /* density function */
  void grad_normal();              /* gradient of density */
#elif DENSITY == 2
  double initialize_dexponential();/* initialize density function */
  double dexponential();           /* density function */
  void grad_dexponential();        /* gradient of density */
#elif DENSITY == 3
  double initialize_gugel();       /* initialize density function */
  double gugel();                  /* density function */
  void grad_gugel();               /* gradient of density */
#elif DENSITY == 4
  double initialize_gugel2();      /* initialize density function */
  double gugel2();                 /* density function */
  void grad_gugel2();              /* gradient of density */
#elif DENSITY == 5
  double initialize_thimble();     /* initialize density function */
  double thimble();                /* density function */
  void grad_thimble();             /* gradient of density */
#else
  printf("DENSITY must be set to 1,2,3,4 or 5\n");
  exit(1);
#endif  

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

  /* default parameters for constructing hat function */
  hc = get_constr_default();

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

  /** overwrite default parameters                               **/
  hc->steps_min = 4;

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

  /* volume below density (default is 0. == unknown) */
#if DENSITY == 1
  hc->volume_density = initialize_normal();
#elif DENSITY == 2
  hc->volume_density = initialize_dexponential();
#elif DENSITY == 3
  hc->volume_density = initialize_gugel();
#elif DENSITY == 4
  hc->volume_density = initialize_gugel2();
#elif DENSITY == 5
  hc->volume_density = initialize_thimble();
#endif  

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

#if RECTANGLE == 1
  /* set rectangle */
  for( i=0; i<N; i++ ) {
    hc->rl[i] = lower[i];
    hc->ru[i] = upper[i];
  }
#endif

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

#if MODE == 1
  /* we use the origin as starting point for         */
  /* searching the mode.                             */
  /* since this is the default, we have nothing to   */
  /* do here.                                        */
#endif  

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

  /* hat */
#if DENSITY == 1
  get_mode(hc,normal);                    /* find mode */
  hat = get_hat(hc,normal,grad_normal);   /* hat       */
#elif DENSITY == 2
  get_mode(hc,dexponential);                /* find mode */
  hat = get_hat(hc,dexponential,grad_dexponential);
#elif DENSITY == 3
  get_mode(hc,gugel);                     /* find mode */
  hat = get_hat(hc,gugel,grad_gugel);
#elif DENSITY == 4
  get_mode(hc,gugel2);                    /* find mode */
  hat = get_hat(hc,gugel2,grad_gugel2);
#elif DENSITY == 5
  get_mode(hc,thimble);                   /* find mode */
  hat = get_hat(hc,thimble,grad_thimble);
#endif  

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

  /* sample */
  for( i=0; i<SAMPLE_SIZE; i++ ) {
    sample_point(hat,rpoint);
    print_rpoint(rpoint,N);
  }

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

#if DEBUG > 0
  db_close(hat);
#endif
  
  exit(0);
} /* end of main() */

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

static void print_rpoint(double *v, int n)
/*
  print random point
*/
{
  int i;

  printf("( %f",v[0]);
  for( i=1; i<n; i++ )
    printf(", %f",v[i]);
  printf(" )\n");

} /* end of db_print_rpoint() */

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

/*******************************************************************

  rotation

  ******************************************************************/

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

/* quick and dirty: the parameters are global */
static double rotangle[N-1];      /* storing angles for rotation */
static double sinangle[N-1];      /* sin of angles for rotation */
static double cosangle[N-1];      /* cos of angles for rotation */
static double sininvangle[N-1];   /* sin of angles for rotation into original position */
static double cosinvangle[N-1];   /* cos of angles for rotation into original position */

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

static void rotate(double *coord, double *xc)
/*
  rotate and move coordinates
*/
{
#if ROTATION == 1
  int i;
  double x,y;

  xc[0] = coord[0]-mode[0];

  for( i=0; i<N-1; i++ ) {
    x = xc[0];
    y = coord[i+1]-mode[i+1];
    xc[0] = cosangle[i] * x + sinangle[i] * y;
    xc[i+1] = - sinangle[i] * x + cosangle[i] * y;
  }
#else
  /* no rotation */
  /* just move */
  int i;

  for( i=0; i<N; i++ )
    xc[i] = coord[i]-mode[i];
#endif

} /* end of rotate */

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

static void invrotate(double *xc)
/*
  rotate coordinates into original position
*/
{
#if ROTATION == 1
  int i;
  double x,y;

  for( i=N-2; i>=0; i-- ) {
    x = xc[0];
    y = xc[i+1];
    xc[0] = cosinvangle[i] * x + sininvangle[i] * y;
    xc[i+1] = - sininvangle[i] * x + cosinvangle[i] * y;
  }
#else
  /* nothing to do */
#endif

} /* end of invrotate */

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

#define DPI         6.28318530717959   /* 2*Pi */

static void initialize_rotation(void)
/*
  set angles for rotation randomly
*/
{
  int i;
  extern double uniform();

#if ROTATION == 1
    printf("\nROTATION\nangles:");
#endif

  /* set angles for rotation of points */
  for( i=0; i<N-1; i++ ) {
    rotangle[i] = uniform() * DPI;
    /* rotangle[i] = 0.; */      /* no rotation */
    sinangle[i] = sin( rotangle[i] );
    cosangle[i] = cos( rotangle[i] );
    sininvangle[i] = sin( -rotangle[i] );
    cosinvangle[i] = cos( -rotangle[i] );
#if ROTATION == 1
    printf("\t%g\n",rotangle[i]);
#endif
  }
} /* end of initialize_rotation() */

#undef DPI

/*-----------------------------------------------------------------*/
#if DENSITY == 1
/*-----------------------------------------------------------------*/

/*******************************************************************

  normal distribution

  ******************************************************************/

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

static double normalcoeff[N];  /* coefficients a_i for normal distribution */

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

double initialize_normal(void)
/*
  create random normal density
  returns volume below density function
*/
{
  int i;
  double volume_density;   /* volume below density */
  extern double uniform();

  for( i=0; i<INIR1; i++ )
    uniform();

  printf("NORMAL distribution\ncoefficients:\n");

  /* coefficients of normal distributions */
  normalcoeff[0] = LOWESTCOEFF;
  printf("\t%g\n",normalcoeff[0]);
  for( i=1; i<N; i++ ) {
    normalcoeff[i] = ( (RATIOCOEFF - 1.) * uniform() + 1. ) * LOWESTCOEFF;
    printf("\t%g\n",normalcoeff[i]);
  }

  for( i=0; i<INIR2; i++ )
    uniform();

  /* set angles for rotation of points */
  initialize_rotation();

  /* compute volume below density */
#if RECTANGLE == 1
  volume_density = 0.;
#else
  volume_density = 1.;
  for( i=0; i<N; i++ )
    volume_density *= PI / normalcoeff[i];
#endif
  return( sqrt(volume_density) );

} /* end of initialize_normal() */

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

double normal(double *coord)
/*
  normal density function
*/
{
  int i;             /* aux counter */
  double f;       /* density function */
  double x[N];

  /* rotate point */
  rotate(coord,x);

  /* gaussian density */
  f=0.;
  for( i=0; i<N; i++ )
    f += normalcoeff[i] * x[i] * x[i];

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

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

void grad_normal(double *coord, double *grad)
/*
  gradient of normal density function
*/
{
  int i;              /* aux counter */
  double f;             /* density at point */
  double x[N];

  /* rotate point */
  rotate(coord,x);

  /* gaussian density */
  f=0.;
  for( i=0; i<N; i++ )
    f += normalcoeff[i] * x[i] * x[i];
  f = exp(-f);

  for( i=0; i<N; i++ )
    grad[i] = (-2.) * normalcoeff[i] * x[i] * f;

  /* inverse rotate gradient */ 
  invrotate(grad);

} /* end of grad_normal() */

/*-----------------------------------------------------------------*/
#elif DENSITY == 2
/*-----------------------------------------------------------------*/

/*******************************************************************

  "exponential" distribution

  ******************************************************************/

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

static double coeff[N];  /* coefficients a_i for distribution */

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

double initialize_dexponential(void)
/*
  create random "exponetial" density
  returns volume below density function
*/
{
  int i;
  double volume_density;   /* volume below density */
  extern double uniform();

  for( i=0; i<INIR1; i++ )
    uniform();

  printf("EXPONENTIAL distribution\ncoefficients:\n");

  /* coefficients of exponential distributions */
  coeff[0] = LOWESTCOEFF;
  printf("\t%g\n",coeff[0]);
  for( i=1; i<N; i++ ) {
    coeff[i] = ( (RATIOCOEFF - 1.) * uniform() + 1. ) * LOWESTCOEFF;
    printf("\t%g\n",coeff[i]);
  }

  for( i=0; i<INIR2; i++ )
    uniform();

  /* set angles for rotation of points */
  initialize_rotation();

  /* compute volume below density */
#if RECTANGLE == 1
  volume_density = 0.;
#else
  volume_density = 1.;
  for( i=0; i<N; i++ )
    volume_density *= 2 / coeff[i];
#endif
  return( volume_density );

} /* end of initialize_dexponential() */

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

double dexponential(double *coord)
/* 
   density function 
*/
{
  int i;             /* aux counter */
  double f;       /* density function */
  double x[N];

  /* rotate point */
  rotate(coord,x);

  /* density */
  f=0.;
  for( i=0; i<N; i++ )
    f += coeff[i] * fabs(x[i]);

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

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

void grad_dexponential(double *coord, double *grad)
/*
  gradient of density function 
*/
{
  int i;              /* aux counter */
  double f;             /* density at point */
  double x[N];

  /* rotate point */
  rotate(coord,x);

  /* density */
  f=0.;
  for( i=0; i<N; i++ )
    f += coeff[i] * fabs(x[i]);
  f = exp(-f);

  /* gradient */
  for( i=0; i<N; i++ ) {
    grad[i] = coeff[i] * f;
    grad[i] *= (x[i] > 0.) ? -1. : 1. ;
  }

  /* inverse rotate gradient */ 
  invrotate(grad);

} /* end of grad_dexponential() */

/*-----------------------------------------------------------------*/
#elif DENSITY == 3
/*-----------------------------------------------------------------*/

/*******************************************************************

  gugelhupf distribution

  ******************************************************************/

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

static double coeff[N];  /* coefficients a_i for distribution */

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

double initialize_gugel(void)
/*
  create random gugelhupf density
  returns volume below density function
*/
{
  int i;
  double volume_density;   /* volume below density */
  extern double uniform();

  for( i=0; i<INIR1; i++ )
    uniform();

  printf("GUGELHUPF distribution\ncoefficients:\n");

  /* coefficients of gugelhupf distributions */
  coeff[0] = LOWESTCOEFF;
  printf("\t%g\n",coeff[0]);
  for( i=1; i<N; i++ ) {
    coeff[i] = ( (RATIOCOEFF - 1.) * uniform() + 1. ) * LOWESTCOEFF;
    printf("\t%g\n",coeff[i]);
  }

  for( i=0; i<INIR2; i++ )
    uniform();

  /* set angles for rotation of points */
  initialize_rotation();

  /* compute volume below density */
#if RECTANGLE == 1
  volume_density = 0.;
#else
  volume_density = 1.;
  for( i=0; i<N; i++ )
    volume_density *= 4 * sqrt( 1. / coeff[i] ) / 3.;
#endif
  return( volume_density );

} /* end of initialize_gugel() */

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

double gugel(double *coord)
/* 
   density function 
*/
{
  int i;             /* aux counter */
  double f,tmp;       /* density function */
  double x[N];

  /* rotate point */
  rotate(coord,x);

  /* gugelhupf density */
  f=1.;
  for( i=0; i<N; i++ ) {
    tmp = 1.-x[i]*x[i]*coeff[i];
    if( tmp <= 0. ) return 0.;
    f *= tmp;
  }

  return( f );
} /* end of gugel() */

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

void grad_gugel(double *coord, double *grad)
/*
  gradient of density function 
*/
{
  int i,j;              /* aux counter */
  double f,tmp;             /* density at point */
  double x[N];

  /* rotate point */
  rotate(coord,x);

  f=1.;
  for( i=0; i<N; i++ ) {
    tmp = 1.-x[i]*x[i]*coeff[i];
    if( tmp <= 0. ) {
      f= 0.;
      break;
    }
    f *= tmp;
  }

  if( f == 0. )
    for( i=0; i<N; i++ )
      grad[i] = 0.;
  else
    for( j=0; j<N; j++ )
      grad[j] = f * (-2*x[j]*coeff[j] ) / (1. - x[j]*x[j]*coeff[j]);
    
  /* inverse rotate gradient */ 
  invrotate(grad);

} /* end of grad_gugel() */

/*-----------------------------------------------------------------*/
#elif DENSITY == 4
/*-----------------------------------------------------------------*/

/*******************************************************************

  gugelhupf distribution

  ******************************************************************/

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

static double coeff[N];  /* coefficients a_i for distribution */

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

double initialize_gugel2(void)
/*
  create random gugelhupf "exponetial" density
  returns volume below density function
*/
{
  int i;
  double volume_density;   /* volume below density */
  extern double uniform();

  for( i=0; i<INIR1; i++ )
    uniform();

  printf("GUGELHUPF 2 distribution\ncoefficients:\n");

  /* coefficients of gugelhupf distributions */
  coeff[0] = LOWESTCOEFF;
  printf("\t%g\n",coeff[0]);
  for( i=1; i<N; i++ ) {
    coeff[i] = ( (RATIOCOEFF - 1.) * uniform() + 1. ) * LOWESTCOEFF;
    printf("\t%g\n",coeff[i]);
  }

  for( i=0; i<INIR2; i++ )
    uniform();

  /* set angles for rotation of points */
  initialize_rotation();

  /* compute volume below density */
#if RECTANGLE == 1
  volume_density = 0.;
#else
  volume_density = 1.;
  for( i=0; i<N; i++ )
    volume_density *= 2*coeff[i] / (1.+coeff[i]);
#endif
  return( volume_density );

} /* end of initialize_gugel2() */

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

double gugel2(double *coord)
/* 
   density function 
*/
{
  int i;             /* aux counter */
  double f,tmp;       /* density function */
  double x[N];

  /* rotate point */
  rotate(coord,x);

  /* gugelhupf density */
  f=1.;
  for( i=0; i<N; i++ ) {
    tmp = 1.- pow( fabs(x[i]), coeff[i]);
    if( tmp <= 0. ) return 0.;
    f *= tmp;
  }

  return( f );
} /* end of gugel2() */

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

void grad_gugel2(double *coord, double *grad)
/*
  gradient of density function 
*/
{
  int i;                /* aux counter */
  double f;             /* density at point */
  double factor[N];     /* factors of density function */
  double x[N];

  /* rotate point */
  rotate(coord,x);

  f=1.;
  for( i=0; i<N; i++ ) {
    factor[i] = 1.- pow( fabs(x[i]), coeff[i]);
    if( factor[i] <= 0. ) {
      f = 0.;
      break;
    }
    f *= factor[i];
  }
  
  /* f == 0 ==> gradient = 0 */
  if( f == 0. ) {
    for( i=0; i<N; i++ )
      grad[i] = 0.;
    return;
  }

  /* compute gradient */
  for( i=0; i<N; i++ )
    grad[i] =  ( (x[i]<0.) ? 1. : -1. ) * f * coeff[i] * pow(fabs(x[i]), coeff[i]-1.) / factor[i];
    
  /* inverse rotate gradient */ 
  invrotate(grad);

} /* end of grad_gugel2() */

/*-----------------------------------------------------------------*/
#elif DENSITY == 5
/*-----------------------------------------------------------------*/

/*******************************************************************

  thimble distribution

  ******************************************************************/

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

static double coeff[N];  /* coefficients a_i for distribution */

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

double initialize_thimble(void)
/*
  create random thimble density
  returns volume below density function
*/
{
  int i;
  double volume_density;   /* volume below density */
  extern double uniform();

  for( i=0; i<INIR1; i++ )
    uniform();

  printf("THIMBLE distribution\ncoefficients:\n");

  /* coefficients of thimble distributions */
  coeff[0] = LOWESTCOEFF;
  printf("\t%g\n",coeff[0]);
  for( i=1; i<N; i++ ) {
    coeff[i] = ( (RATIOCOEFF - 1.) * uniform() + 1. ) * LOWESTCOEFF;
    printf("\t%g\n",coeff[i]);
  }

  for( i=0; i<INIR2; i++ )
    uniform();

  /* set angles for rotation of points */
  initialize_rotation();

  /* compute volume below density */
  volume_density = 0.;   /* not known */
  return( volume_density );

} /* end of initialize_platy() */

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

double thimble(double *coord)
/* 
   density function 
*/
{
  int i;             /* aux counter */
  double f;       /* density function */
  double x[N];

  /* rotate point */
  rotate(coord,x);

  /* thimble density */
  f=0.;
  for( i=0; i<N; i++ ) {
    f += x[i]*x[i]*coeff[i];
  }

  f = (f<=1.) ? 1.-f : 0.;
  
  return( f );
} /* end of thimble() */

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

void grad_thimble(double *coord, double *grad)
/*
  gradient of density function 
*/
{
  int i;              /* aux counter */
  double f;             /* density at point */
  double x[N];

  /* rotate point */
  rotate(coord,x);

  f=0.;
  for( i=0; i<N; i++ ) {
    f += x[i]*x[i]*coeff[i];
  }

  f = (f<=1.) ? 1.-f : 0.;

  if( f == 0. )
    for( i=0; i<N; i++ )
      grad[i] = 0.;
  else
    for( i=0; i<N; i++ )
      grad[i] = -2.*x[i]*coeff[i];
    
  /* inverse rotate gradient */ 
  invrotate(grad);

} /* end of grad_thimble() */

/*-----------------------------------------------------------------*/
#endif
/*-----------------------------------------------------------------*/