/********************************************************************/ /** **/ /** tdrmv.h **/ /** **/ /** A universal generator for multivariate T-concave **/ /** distributions **/ /** **/ /** header file **/ /** **/ /** author: Josef Leydold **/ /** last modified: Mon Mar 9 13:38:00 MET 1998 **/ /** **/ /********************************************************************/ /*------------------------------------------------------------------*/ /** config file **/ #include "config.h" /** the libraries **/ #include #include #include #include #include #include /*------------------------------------------------------------------*/ /* check compiling options */ #if RECTANGLE == 1 #undef MODE #define MODE 1 #undef GAMMA #define GAMMA 2 #endif /*------------------------------------------------------------------*/ /** global variables **/ /* file handle for log file */ FILE *LOG; /** debugging constanst **/ #define DB_VERTICES 2 #define DB_CONES 4 #define DB_EDGES 8 #define DB_CPARAMS 16 #define DB_RPOINT 32 #define DB_GUIDE 64 #define DB_GAMMA 128 #define DB_VOLUME 1024 #define DB_G 2048 /*------------------------------------------------------------------*/ /** Macros **/ /* write message, file name and line number and exit */ #define FATAL(str) fatal((str),__FILE__,__LINE__) /* write message, file name and line number */ #define WARNING(str) warning((str),__FILE__,__LINE__) /* maximum of two numbers */ #define MAX(x,y) (((x)>(y))?(x):(y)) /* minimum of two numbers */ #define MIN(x,y) (((x)<(y))?(x):(y)) /** Transformation **/ /* this version supports T(x) = log(x) only */ #define T(x) (log(x)) /* transformation function */ #define T_deriv(x) (1./(x)) /* first derivative of transform funtion */ #define T_inv(x) (exp(x)) /* inverse of transform function */ /*------------------------------------------------------------------*/ /** Structures **/ typedef struct s_constr /* parameters for constructing hat function */ { #if MODE == 1 double mode[N]; /* mode of distribution */ #endif #if RECTANGLE == 1 double rl[N]; /* lower bounds for rectangle */ double ru[N]; /* upper bounds for rectangle */ #endif int max_cones; /* maximum number of cones (at least 2^(N+T_STEPS_MIN) */ int steps_min; /* minimum number of triangulation steps */ int step_tp; /* triangulation step when optimal touching points is calculated */ double reject_tp_abs; /* absolute value of volume Hi for which touching point is rejected */ double reject_tp_rel; /* max/min ration of volume Hi for which touching point is rejected */ double reject_tp_percentile; /* percentile for rejection */ int vertex_alloc_block; /* number of vertices to be allocated at once */ int cone_alloc_block; /* number of cones to be allocated at once */ double find_tp_start; /* starting point for finding optimal touching point */ int find_tp_steps_r; /* max number of steps for finding proper touching point, --> Infinity */ int find_tp_steps_l; /* max number of steps for finding proper touching point, --> 0 */ double find_tp_step_size; /* size of each step for finding proper touching point */ int find_tp_steps_lb; /* max number of steps for finding lower bound of bracket */ int find_tp_steps_ub; /* max number of steps for finding upper bound of bracket */ double find_tp_tol; /* acceptable tolerance for Brent's algorithm for finding minimum */ double find_tp_steps_brent; /* max number of interation steps for Brent's algorithm */ #if MODE == 1 double hooke_rho_begin; /* stepsize geometric shrink */ double hooke_rho_second; /* stepsize geometric shrink (second run) */ double hooke_epsilon; /* ending value of stepsize */ int hooke_imax; /* max # of iterations */ double mode_to_boundary; /* move mode to boundary if |mode - boundary| / length < MODE_TO_BOUNDARY */ #endif int guide_rel_size; /* size of guide table */ double volume_density; /* volume below density (for debugging only) */ } C_CONSTR; typedef struct s_vertex /* data for one vertex */ { struct s_vertex *next; /* pointer to next vertex in list */ int index; /* index of vertex */ double coord[N]; /* coordinates of spanning vector(norm = 1), "vertex" */ double norm; /* norm of vertex */ } VERTEX; typedef struct s_cone /* a cone */ { struct s_cone *next; /* pointer to next cone in list */ #if DEBUG > 0 int index; /* index of cone */ #endif int level; /* level of triangulation */ VERTEX *v[N]; /* list of vertices of the cone */ double center[N]; /* barycenter of cone */ double detf; /* determinant / (N-1)! for cone */ double alpha; /* parameters for hat function */ double beta; double gv[N]; /* for all vertices v */ double ai; /* coefficient for marginal density */ double tp; /* coordinate of touching point */ double Hi; /* volume under hat in cone */ double Hsum; /* accumulated sum of volumes */ #if RECTANGLE == 1 double height; /* height of pyramid */ double tdrg_vol; /* volume below hat of gamma density */ #endif #if DEBUG > 1000 double fp; /* value of density at touching point */ double g[N]; /* direction of sweep-plane */ double Vi; /* volume of triangle of distance 1 */ #endif } CONE; typedef struct s_root /* root of tree of cones */ { double (*density)(); /* pointer to density function */ void (*grad_density)(); /* pointer to density function */ #if MODE == 1 double mode[N]; /* mode of distribution */ #endif #if RECTANGLE == 1 /* rectangle moved to mode == origin */ double rl[N]; /* lower bounds for rectangle */ double ru[N]; /* upper bounds for rectangle */ #endif CONE *c; /* root of list of cones */ VERTEX *v; /* root of list of vertices */ double Htot; /* total volume below hat */ double Hi_median; /* median of the Hi's */ double Hi_percentile; /* percentile of the Hi's */ int max_v; /* maximum number of vertices */ VERTEX *last_v; /* pointer to last vertex in list */ int max_c; /* maximum number of cones */ CONE *last_c; /* pointer to last cone in list */ int steps_min; /* minimum number of triangulation steps */ int steps_max; /* maximum number of triangulation steps */ int step_tp; /* triangulation step, when touching point is calculated */ CONE **guide; /* pointer to guide table */ int size_guide; /* size of guide table */ C_CONSTR *hc; /* pointer to construction parameters for hat function */ #if RECTANGLE == 1 double max_gamma; /* maximum value for gamma variaties */ #endif #if DEBUG double volume_density; /* volume below density (for debugging only) */ int db_n_rpoint_accept; /* number of accepted random points (for debugging only) */ int db_n_rpoint_reject; /* number of rejected random points (for debugging only) */ #if RECTANGLE == 1 int db_n_rpoint_reject_outside; /* number of rejected random points outside domain (for debugging only) */ #endif #if GAMMA == 2 int db_n_tdrg_accept; /* number of accepted random points (gamma dist.; for debugging only) */ int db_n_tdrg_reject; /* number of rejected random points (gamma dist.; for debugging only) */ #endif #endif #if GAMMA == 2 int tdrg_ntp; /* total number of touching points for constructing hat function of gamma density */ #endif } C_ROOT; /*------------------------------------------------------------------*/ /** Functions for user interface **/ /* get default values for construction of hat function (hat.c) */ C_CONSTR *get_constr_default(void); /* compute cones and hat function (hat.c) */ C_ROOT *get_hat(C_CONSTR *hc, double (*density)(), void (*grad_density)()); /* sample random point (rpoint.c) */ void sample_point(C_ROOT *cr, double *rpoint); #if MODE == 1 /* find mode of given density (mode.c) */ int get_mode(C_CONSTR *hc, double (*density)()); #endif #if DEBUG > 0 /* write statistics into log file and close log file (debug.c) */ void db_close(C_ROOT *cr); #endif