/********************************************************************/
/**                                                                **/
/**   rpoint.c                                                     **/
/**                                                                **/
/**   A universal generator for multivariate T-concave             **/
/**   distributions                                                **/
/**                                                                **/
/**   get random point                                             **/
/**                                                                **/
/**   author: by Josef Leydold                                     **/
/**   last modified: Mon Mar  9 13:38:00 MET 1998                  **/
/**                                                                **/
/********************************************************************/              

/* the header file */
#include "tdrmv.h"

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

/** local definitions                                              **/

/* internal functions */
/* get point below hat function                                     */
__inline__ static double rp_get_point(C_ROOT *cr, CONE *c, double *rpoint);

/* sample from uniform distribution on standard simplex             */
__inline__ static void sample_uniform(double *u);                   


/* external functions                                               */
/* uniform random number generator */
extern double uniform();       

/* generator for gamma variates                                     */                 
#if GAMMA == 1
extern double gds(double);        /* shape parameters a=N and b = 1 */
#else
extern double tdrg(C_ROOT *cr, CONE *c);      /* a = N and b = beta */                      
#endif

/* for debugging                                                    */
/* print a vector into log file                                     */
extern void db_print_vector(char *begin, double *x, int n, char *end);
/* write error message into log file                                */
extern void warning(char *msg, char *filename, int line);
/* fatal error, exit                                                */
extern void fatal(char *msg, char *filename, int line);

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

void sample_point(C_ROOT *cr, double *rpoint)
/*
  sample random point
*/
{
  CONE *c;       /* cone for generating point */
  double gx;     /* distance of random point */
  double u;      /* uniformly distributed random number */
  double f, h;   /* value of density and hat at random point */
#if MODE == 1
  int i;
#endif
#if RECTANGLE == 1
  int ok;
#endif
#if DEBUG >0 
  double uh;
#endif

  while( 1 ) { /* loop until random point is accepted */

    /*.................................................................*/
    /** find a cone **/

    u = uniform();    /* sample from uniform distribution */

    /* look up in guide table and search for cone */
    for( c =  (cr->guide)[(int) (u*cr->size_guide)]; 
       c!=NULL && c->Hsum < cr->Htot * u;
       c=c->next );
#if DEBUG > 0
    if( c == NULL ) {
      /* something is wrong. This should not happen */
      WARNING( "c==NULL" );
      c=cr->last_c;      /* we simple use the last cone */
    }
    if( DEBUG & DB_RPOINT )
      fprintf(LOG,"H = %g (%g) --> cone = %d",cr->Htot * u,cr->Htot,c->index);
#endif

    /*.................................................................*/
    /** get random point and distance of hyper plane **/
    gx = rp_get_point(cr,c,rpoint);

#if RECTANGLE == 1
    /* is point inside domain (rectangle) ? */
    ok = 1;
    for( i=0; i<N; i++ )
      if( (rpoint[i] < cr->rl[i]) || (rpoint[i] > cr->ru[i]) ) {
	/* point not in rectangle --> reject */
	ok = 0;
	break;
      }
    if( !ok ) {
#if DEBUG >0
      if( DEBUG & DB_RPOINT )
	fprintf(LOG,"\toutside domain\t\t\t\t\t\t-- > reject\n");
      cr->db_n_rpoint_reject++;
      cr->db_n_rpoint_reject_outside++;
#endif
      continue;   /* try again */
    }
#endif

#if MODE == 1
    /** move point into original location */
    for( i=0; i<N; i++ )
      rpoint[i] = rpoint[i] + cr->mode[i];
#endif

    /** value of ... **/
    f = (*cr->density)(rpoint);             /* density */
    h = T_inv( c->alpha - c->beta * gx );   /* hat */

    /*.................................................................*/
    /** accept or reject **/
#if DEBUG >0
    uh = uniform() * h;
    if( DEBUG & DB_RPOINT ) {
      fprintf(LOG,"\that = %g\tf = %g\tu*hat = %g",h,f,uh);
      if( uh <= f )
	fprintf(LOG,"\t-- > accept\n");
      else
	fprintf(LOG,"\t-- > reject\n");
    }
    
    /* check concavity */
    if( h < f * (1.-TOLERANCE) )
      WARNING( "density not T-concave" );
    
    /* accept or reject */
    if( uh <= f ) {
      cr->db_n_rpoint_accept++;
      return;
    }
    else
      cr->db_n_rpoint_reject++;
#else
    if( uniform() * h <= f )
      return;
#endif
  }

} /* end of sample_point() */

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

__inline__ static double rp_get_point(C_ROOT *cr, CONE *c, double *rpoint)
/*
  find hyperplane
*/
{
  int i,j;
  double gx, x;
  double u[N];

  /* get x value for marginal distribution of hat */
#if RECTANGLE == 1
  gx = tdrg(cr,c);
#else
#if GAMMA == 1
  gx = gds((double)N) / (c->beta);
#else
  gx = tdrg(cr,c);
#endif
#endif

  /* uniform random numbers with sum u_i = 1 */
  sample_uniform(u);

  /* clear array */
  for( i=0; i<N; i++ )
    rpoint[i] = 0.;

  /* calculate random point */
  for( j=0; j<N; j++ ) {
    x = gx * u[j] / (c->gv)[j];
    for( i=0; i<N; i++ )
      rpoint[i] += x * ((c->v)[j]->coord)[i];
  }

#if DEBUG > 0
  if( DEBUG & DB_RPOINT ) {
    fprintf(LOG,"   \tgx = %g\n",gx);
    db_print_vector("   rpoint = ",rpoint,N,"\n");
  }
#endif

  return gx;
} /* end of rp_get_point() */

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

__inline__ static void sample_uniform(double *u)
/* sample from uniform distribution on standard simplex */
#if N == 2
{
  u[0] = uniform();
  u[1] = 1. - u[0];
}
/*................................................................*/
#elif N == 3
{
  u[0] = uniform();
  u[1] = uniform();
  if( u[0] > u[1] ) {
    u[2] = u[0];
    u[0] = u[1];
    u[1] = u[2];
  }
  u[2] = 1. - u[1];
  u[1] = u[1] - u[0];
}
/*................................................................*/
#else
{
  int i,j;
  double u_aux;

  /* point in standard simplex 0 <= u[0] <= ... <= u[N-2] <= 1 */

  /* generate N-1 numbers */
  for( i=0; i<N-1; i++ )
    u[i] = uniform();

  /* sort numbers (insertion sort) */
  for( i=1; i<N-1; i++ ) {
    u_aux = u[i];
    j = i;
    while( u[j-1] > u_aux && j>0 ) {
      u[j] = u[j-1];
      j--;
    }
    u[j] = u_aux;
  }

  /* transform to affine coordinates */
  u[N-1] = 1.;
  for( i=N-1; i>0; i-- )
    u[i] = u[i] - u[i-1];

} /* end of sample_uniform() */
#endif

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