/********************************************************************/ /** **/ /** hat.c **/ /** **/ /** A universal generator for multivariate T-concave **/ /** distributions **/ /** **/ /** construct hat function for given density **/ /** **/ /** author: by Josef Leydold **/ /** last modified: Mon Mar 9 13:38:00 MET 1998 **/ /** **/ /********************************************************************/ /* the header file */ #include "tdrmv.h" /*------------------------------------------------------------------*/ /** structures **/ typedef struct s_edge_table /* hash table for edges */ { int index[2]; /* index of incident vertices */ VERTEX *vertex; /* index of corresponding vertex (=barycenter) */ struct s_edge_table *next; } E_TABLE; typedef struct s_tp_arg /* argument for tp function */ { double x; /* argument */ double f; /* value f(t) */ CONE *c; /* parameters */ double valid; /* bolean to store if x in domain */ } TP_ARG; /*-----------------------------------------------------------------*/ /** error codes **/ #define TP_VALID 0 /* density at touching point not zero */ #define TP_FZERO 1 /* denstiy 0 */ #define TP_GZERO 2 /* grandient of denstiy 0 */ #define TP_UGRAD 3 /* grandient of transformed denstiy does not fit into cone */ #define TP_NAN 9 /* numerical errors (NaN) */ #define TP_LEFT 1 /* minimum in left point */ #define TP_MIDDLE 2 /* minimum in middle point */ #define TP_RIGHT 3 /* minimum in right point */ #define TP_BRACKET 4 /* bracket found */ #define TP_EMPTY 0 /* no proper touching point found */ /*-----------------------------------------------------------------*/ /** global variables **/ static C_ROOT *cr; /* pointer to structure for hat function */ static double fak[N+1]; /* store factorials for 1, ..., N */ /*-----------------------------------------------------------------*/ /** functions **/ /** cones **/ static void get_cones(void); /* get cones and their vertices */ /** initial cones **/ static void get_initial_vertices(void); /* get vertices of initial cones */ static void get_initial_cones(void); /* get initial cones */ static void calc_constants(void); /* calculate some constants */ /** triangulation **/ static int triangulate(int step,int all); /* make one triangulation step */ static void triangulate_cone(CONE *c, int step); /* triangulate a cone */ static VERTEX *new_vertex(VERTEX **vl); /* compute new vertex on edge */ static int number_vertices(int level); /* number of vertices in given triangulation level */ /** data for cones **/ static void cone_center(CONE *c); /* computer center of cone */ static int cone_params(CONE *c); /* compute parameters for hat for a cone */ #if RECTANGLE == 1 static void get_height(CONE *c); /* calculate height of pyramid */ #endif #if DEBUG > 1000 static void cone_volume(CONE *c); /* volume of triangle of distance 1 */ #endif /** hat function **/ static void get_hat_params(void); /* get parameters for hat function for each cone */ static void hat_cone_params(CONE *c, double Hi_bound); /* calculate and check cone parameters */ static void get_median(void); /* compute median and percentile of Hi's */ static void make_guide_table(void); /* create guide table */ /** distance of touching points **/ static void get_tp(void); /* compute optimal distance of touching points for each cone */ static void find_tp(CONE *c); /* find optimal touching point for cone */ static int find_bracket(TP_ARG *a); /* search for proper bracket of minimum of tp_f2min() */ static int find_proper_tp(TP_ARG *a); /* search for proper touching point */ static void tp_min(TP_ARG *a); /* funtion that must be minimized for optimal touching point */ /** hash table for edges **/ #if N > 2 static void etable_new(int size); /* make new hash table */ static void etable_clear(void); /* clear hash table */ static void etable_clear_next(E_TABLE *et); /* clear entries in edge table by recursion */ static VERTEX *etable_find_or_insert(VERTEX **vidx); /* find or insert entry in hash table, return pointer to vertex */ static E_TABLE **p_etable_alloc(int n); /* allocate pointer to edge table */ static E_TABLE *etable_alloc(int n); /* allocate edge table entries */ #endif /** alloc memory **/ static double *double_alloc(int n); /* allocate array of doubles */ static C_CONSTR *hconstr_alloc(int n); /* allocate parameters for constructing hat function */ static C_ROOT *croot_alloc(int n); /* allocate root of cone tree */ static CONE *cone_alloc(int n); /* allocate cone entries */ static CONE *next_cone(void); /* next allocated cone in list */ static VERTEX *vertex_alloc(int n); /* allocate vertex entries */ static VERTEX *next_vertex(void); /* next allocated vertex in list */ static CONE **guide_alloc(int n); /* allocate pointer to array of pointers to cones and initialize */ /** misc **/ static double vec_norm(double *x, int n); /* norm of vector */ static void sort_varray(VERTEX **va); /* sort vertex array */ static int ipower(int n, int k); /* power n^k of integers */ #if RECTANGLE == 1 static double igammaf(int a, double x); /* incomplete gamma function * gamma(a), for integer a, a>=2 */ #endif static double vec_scalarprod(double *, double *, int); /* scalar product of two vectors */ #if N > 2 static void copy_varray(VERTEX **a, VERTEX **b, int n); /* copy first VERTEX array of n elements into second */ #endif /** routines from C math library **/ static double fminbr(TP_ARG *ap, double tol); /* find minimum in an intervall */ /** routines from numerical recipes in C **/ static void hpsort(unsigned long n, double *ra); /* heap sort */ /*-----------------------------------------------------------------*/ /* external functions */ extern void db_initialize(); /* initialize debuging variables */ extern void db_hat(C_ROOT *cr); /* write statistics about hat generation into log file */ #if GAMMA != 1 extern void tdrg_setup(C_ROOT *cr); /* setup for generator tdrg of gamma distributed variates */ #endif extern void warning(char *msg, char *filename, int line); /* write erro message */ extern void fatal(char *msg, char *filename, int line); /* fatal error, exit */ /*-----------------------------------------------------------------*/ /******************************************************************* main ******************************************************************/ /*-----------------------------------------------------------------*/ C_CONSTR *get_constr_default(void) /* get default values for construction of hat function */ { C_CONSTR *hc; /* allocate memory */ hc = hconstr_alloc(1); /** set default values **/ /* the mode */ #if MODE == 1 { int i; for( i=0; imode[i] = 0.; /* mode is set to origin by default */ } #endif /* control generation of cones */ hc->max_cones = MAX_N_CONES; /* maximum number of cones (at least 2^(N+T_STEPS_MIN) */ hc->steps_min = T_STEPS_MIN; /* minimum number of triangulation steps */ hc->step_tp = OPTIMAL_TP_STEP; /* triangulation step when optimal touching points is calculated */ /* reject bad touching points */ hc->reject_tp_abs = REJECT_TP_ABS; /* reject if Hi is greaten than this value */ hc->reject_tp_rel = REJECT_TP_REL; /* reject if Hi > median + reject_tp_rel * (percentile - median) */ hc->reject_tp_percentile = REJECT_TP_PERCENTILE;/* percentile for rejection */ /* control allocation of memory */ hc->cone_alloc_block = CONE_ALLOC_BLOCK; /* number of cones to be allocated at once */ hc->vertex_alloc_block = VERTEX_ALLOC_BLOCK;/* number of vertices to be allocated at once */ /* parameters for finding touching points */ hc->find_tp_start = FIND_TP_START; /* starting point for finding optimal touching point */ hc->find_tp_steps_r = FIND_TP_STEPS_R; /* max number of steps for finding proper touching point, --> Infinity */ hc->find_tp_steps_l = FIND_TP_STEPS_L; /* max number of steps for finding proper touching point, --> 0 */ hc->find_tp_step_size = FIND_TP_STEP_SIZE; /* size of each step for finding proper touching point */ hc->find_tp_steps_lb = FIND_TP_STEPS_LB; /* max number of steps for finding lower bound of bracket */ hc->find_tp_steps_ub = FIND_TP_STEPS_UB; /* max number of steps for finding upper bound of bracket */ hc->find_tp_tol = FIND_TP_TOL; /* acceptable tolerance for Brent's algorithm for finding minimum */ hc->find_tp_steps_brent = FIND_TP_STEPS_BRENT;/* max number of interation steps for Brent's algorithm */ #if MODE == 1 /* parameters for finding mode */ /* see hooke.c for explanation */ hc->hooke_rho_begin = HOOKE_RHO_BEGIN; /* stepsize geometric shrink */ hc->hooke_rho_second = HOOKE_RHO_SECOND; /* stepsize geometric shrink (second run, only if > 0) */ hc->hooke_epsilon = HOOKE_EPSILON; /* ending value of stepsize */ hc->hooke_imax = HOOKE_IMAX; /* max # of iterations */ hc->mode_to_boundary = MODE_TO_BOUNDARY; /* move mode to boundary if |mode - boundary| / length < MODE_TO_BOUNDARY */ #endif #if RECTANGLE == 1 /* set lower bound > upper bound --> thus we can check if the user has set rectangle boundaries */ { int i; for( i=0; irl[i] = 1.; hc->ru[i] = -1.; } } #endif /* size of search table */ hc->guide_rel_size = GUIDE_TABLE_SIZE; /* relative size of guide table compared to number of cones */ hc->volume_density = 0.; /* this is needed for debugging */ return hc; } /* end of get_constr_default() */ /*-----------------------------------------------------------------*/ C_ROOT *get_hat(C_CONSTR *hc, double (*density)(), void (*grad_density)()) /* get hat function */ { #if DEBUG > 0 /* open log file */ db_initialize(); #endif /* alloc memory for root */ cr = croot_alloc(1); /* paramters for density */ cr->hc = hc; /* parameters for constructing hat function */ cr->density = density; /* density function */ cr->grad_density = grad_density; /* gradient of density function */ #if MODE == 1 { int i; for( i=0; imode[i] = hc->mode[i]; /* coordinates of mode of density */ } #endif #if RECTANGLE == 1 /* check bounds of rectangle */ { int i; for( i=0; irl[i] = hc->rl[i] - hc->mode[i]) >= (cr->ru[i] = hc->ru[i] - hc->mode[i]) ) FATAL("rectangle: lower bound >= upper bound"); /* check if given mode is domain */ for( i=0; imode[i] < hc->rl[i] || hc->mode[i] > hc->ru[i] ) FATAL( "mode not in given domain (rectangle)" ); } #endif #if DEBUG > 0 /* initialize counters (for debugging only) */ cr->volume_density = hc->volume_density; /* this is needed for debugging */ cr->db_n_rpoint_accept = 0; cr->db_n_rpoint_reject = 0; #if RECTANGLE == 1 cr->db_n_rpoint_reject_outside = 0; #endif #if GAMMA == 2 cr->db_n_tdrg_accept = 0; cr->db_n_tdrg_reject = 0; #endif #endif /* get cones and its vertices */ get_cones(); /* calculate parameters for hat function */ get_hat_params(); /* we do not need the hash table generated in triangulate() cone any more */ #if N > 2 etable_clear(); #endif #if GAMMA == 2 /* at last we have to call the setup routine for the gamma random number generator */ tdrg_setup(cr); #endif #if DEBUG > 0 /* write data about hat funtion into log file */ db_hat(cr); #endif return cr; } /* end of get_hat() */ /*-----------------------------------------------------------------*/ /******************************************************************* cones ******************************************************************/ /*-----------------------------------------------------------------*/ static void get_cones(void) /* get cones and its vertices */ { int step; /* triangulation steps */ /* constants and parameters */ calc_constants(); /* no triangulation steps yet */ cr->steps_max = -1; /* vertices of initial cones */ get_initial_vertices(); /* initial cones */ get_initial_cones(); /* triangulate until step for computing touching point */ for( step = 1; step <= cr->step_tp; step++ ) triangulate(step,1); /* compute optimal distance of touching points now */ get_tp(); /* continue with triangulation */ for( ; step <= cr->steps_min; step++ ) triangulate(step,1); /* some of cones must be split */ while( triangulate(step,0) ) step++; /* maximum number of triangultions */ cr->steps_max = step-1; } /* end of get_cones() */ /*-----------------------------------------------------------------*/ /******************************************************************* initial cones ******************************************************************/ /*-----------------------------------------------------------------*/ static void get_initial_vertices(void) /* get vertices of initial cones and return number of vertices */ { VERTEX *vt; int i,k; /* number of vertices */ cr->max_v = 2*N; /* alloc memory for vertices */ #if RECTANGLE == 1 cr->v = vertex_alloc(2*N); /* we probably do not need all vertices */ #else cr->v = vertex_alloc(number_vertices(cr->steps_min)); #endif vt = cr->v; /* unit vectors e_k and -e_k */ for( k=0; kcoord)[i] = (i==k) ? 1. : 0.; ((vt+k+N)->coord)[i] = (i==k) ? -1. : 0.; } /* all vectors have norm 1. */ (vt+k)->norm = 1.; (vt+k+N)->norm = 1.; /* index of vertex */ (vt+k)->index = k; (vt+k+N)->index = k+N; #if RECTANGLE == 1 /* all initial vertices which we do not need with a negative integer */ if( (cr->hc)->mode[k] == (cr->hc)->rl[k] ) (vt+k+N)->index = -(vt+k+N)->index -1; /* no cones below lower bound: index -> -index-1 */ if( (cr->hc)->mode[k] == (cr->hc)->ru[k] ) (vt+k)->index = -(vt+k)->index-1; /* no cones above upper bound: index -> -index-1 */ #endif } /* last element in list */ cr->last_v = vt+2*N-1; #if RECTANGLE == 1 /* some of these vertices (at most N) might not be necessary */ /* but it is easier to ignore this case */ #endif } /* end of get_initial_vertices() */ /*-----------------------------------------------------------------*/ static void get_initial_cones(void) /* get initial cones */ { int i,k; int alloc_c; /* number of cones to be allocated */ int max_c; /* maximum number of cones */ CONE *c; double det; #if RECTANGLE == 1 int ok; #endif /* number of cones */ /* we have at least 2^(N+T_STEPS_MIN) cones */ alloc_c = ipower(2,N+cr->steps_min); #if RECTANGLE == 1 /* we might have far less cones, if the mode is at the boundary */ for( i=0; ihc)->mode[i] == (cr->hc)->rl[i] ) || ( (cr->hc)->mode[i] == (cr->hc)->ru[i] ) ) alloc_c /= 2; #endif /* alloc memory for cones */ cr->max_c = 0; /* no cones yet */ cr->c = cone_alloc(alloc_c); c = cr->c; /* factor for a_i (= volume of triangles with vertices in axes) */ det = 1./fak[N-1]; /* we have (at most) 2^N initial cones */ max_c = ipower(2,N); /* we need vertices, index, volume for each cone */ for( k=0; kv)+i+((k>>i)&1)*N)->index < 0 ) { ok = 0; /* cone contains initial vertex outside rectangle */ break; } if( !ok ) /* try next one */ continue; #endif c = (cr->c)+(cr->max_c); ++cr->max_c; /* this is level 0 of triangulation */ c->level = 0; /* each has N edges. */ /* The i-th edge of the cone is either cr->v[i] or cr->v[N+i], */ /* which is equal to (-1)*(cr->v[i]). */ /* WARNING! cr->v is a linked list! */ /* (cr->v)[i] only works for the first 2*N elements !! */ for( i=0; iv)[i] = cr->v + i + ((k>>i)&1)*N; /* we need a sorted list */ sort_varray(c->v); /* determinant and volume of triangle */ c->detf = det; /* touching point not known yet */ c->tp = -1.; /* > 0. if and only if tp is computed !! */ } /* last element in list */ cr->last_c = c; } /* end of get_initial_cones() */ /*-----------------------------------------------------------------*/ static void calc_constants(void) /* setup some constants */ { int i; double tmp; /* first we check parameters */ /* N must be an integer >= 2 */ if( (int) N != N || N < 2 || N > 15 ) FATAL( "N in config.h not an integer >= 2 and <= 15" ); /* some constants */ /* factorial */ fak[0] = 1.; for( i=1; i<=N; i++ ) fak[i] = fak[i-1] * ((double) i); /* we have to check the given constants */ cr->steps_min = MAX( 0, (cr->hc)->steps_min ); /* minimum number of triangulation steps (at least 0 !!) */ cr->steps_max = 0; /* maximum number of triangulation steps (none yet) */ tmp = ipower(2,N+cr->steps_min)+5; if( tmp > (cr->hc)->max_cones ) { WARNING( "number of cones raised to 2^(N+T_STEPS_MIN)" ); (cr->hc)->max_cones = tmp; } /* when should the touching point be computed ? */ cr->step_tp = MAX( 0, (cr->hc)->step_tp ); /* at least 0 !! */ cr->step_tp = MIN( cr->step_tp, cr->steps_min ); /* at most cr->steps_min !! */ } /* end of hat_constants() */ /*-----------------------------------------------------------------*/ /******************************************************************* triangulation ******************************************************************/ /*-----------------------------------------------------------------*/ static int triangulate(int step, int all) /* make one triangulation step */ { int k,nc; CONE *c; #if N > 2 /* We need a hash table for storing the edges */ /* size of table = maximum number of vertices in current triangulation cycle */ /* only at begining of new cycle */ /* length of cycle N-1 */ /* only necessary if N > 2 */ if( step % (N-1) == 1 ) etable_new(number_vertices( (step/(N-1)+1)*(N-1) )); #endif /* number of cones before triangulation */ nc = cr->max_c; /* triangulate every cone */ for( k=0,c=cr->c; ktp < 0. ) { triangulate_cone(c,step); /* triangulate cone */ find_tp(c); /* find new touching point for first new cone */ find_tp(cr->last_c); /* find new touching point for second new cone */ } /* next cone */ c=c->next; } /* return number of new cones */ return (cr->max_c - nc); } /* end of triangulate() */ /*-----------------------------------------------------------------*/ static void triangulate_cone(CONE *c, int step) /* triangulatw a cone */ { CONE *newc; /* new cone */ VERTEX *newv; /* new vertex */ /* find "oldest" edge */ #if N == 2 newv = new_vertex(c->v); #else newv = etable_find_or_insert(c->v); #endif /* construct two new cones */ /* first cone */ newc = next_cone(); /* pointer to new cone */ newc->level = step; /* level */ #if N == 2 (newc->v)[0] = (c->v)[1]; (newc->v)[1] = newv; #else copy_varray( c->v+1, newc->v, N-1 ); /* copy list of vertices to new cone */ (newc->v)[N-1] = newv; /* the new vertex */ #endif newc->detf = c->detf/ (2.*(newv->norm)); /* the determinant of spanning vectors */ newc->tp = c->tp; /* distance of touching point remains unchanged */ /* second cone */ c->level = step; /* level */ #if N > 2 /* moving not necessary if N==2 */ copy_varray( c->v+2, c->v+1, N-2 ); /* shift list of vertices */ #endif (c->v)[N-1] = newv; /* the new vertex */ c->detf = newc->detf; /* the determinant */ } /* end of triangulate_cone() */ /*-----------------------------------------------------------------*/ static VERTEX *new_vertex(VERTEX **vl) /* compute new vertex on edge and return its index */ { int i; VERTEX *newv; /* pointer to new vertex */ /* get an empty vertex */ newv = next_vertex(); /* barycentre of edge */ for( i=0; icoord[i] = 0.5 * ( ((vl[0])->coord)[i] + ((vl[1])->coord)[i] ); /* norm */ newv->norm = vec_norm(newv->coord,N); /* norm --> 1 */ for( i=0; icoord[i] /= newv->norm; return newv; } /* end of new_vertex() */ /*-----------------------------------------------------------------*/ static int number_vertices(int level) /* calculate (approximate) number of vertices in given triangulation level (WARNING! The numbers are found by computer experiments) */ { #if N == 2 /* nothing to declare */ #elif N == 3 static int nv[]={ 6, 13, 18, 40, 66,142,258,538,1026,2098,4098,8290,16386,32962,65538,131458,262146}; static int last = 16; #elif N == 4 static int nv[]={ 8, 19, 25, 32, 80,128,192,456, 824,1408,3120,5968,11008,23264,45600, 87552}; static int last = 15; #elif N == 5 static int nv[]={10, 26, 33, 41, 50,140,220,321, 450,1186,2158,3636, 5890,13970,27130}; static int last = 14; #elif N == 6 static int nv[]={12, 34, 42, 51, 61, 72,224,348, 501, 681, 912,2660, 4896, 8254}; static int last = 13; #elif N == 7 static int nv[]={14, 43, 52, 62, 73, 85, 98,336, 518, 743, 985,1289, 1666}; static int last = 12; #elif N == 8 static int nv[]={16, 53, 63, 74, 86, 99,113,128, 480, 736,1059}; static int last = 10; #elif N == 9 static int nv[]={18, 64, 75, 87,100,114,129,145, 162, 660}; static int last = 9; #elif N ==10 static int nv[]={20, 76, 88,101,115,130,146,163, 181, 200}; static int last = 9; #elif N ==11 static int nv[]={22, 89,102,116,131,147,164,182, 201, 221, 242}; static int last = 10; #else /* N >=12 */ static int nv[]={24,103,117,132,148,165,183,202, 222, 243, 265, 288}; static int last = 11; #endif if( level < 0 ) FATAL( "level < 0" ); #if N == 2 return ipower(N,level+2); #else level = MIN(level,last); return nv[level]; #endif } /* end of number_vertices() */ /*-----------------------------------------------------------------*/ /******************************************************************* data for cones ******************************************************************/ /*-----------------------------------------------------------------*/ static void cone_center(CONE *c) /* computer center of cone and normalize the corresponding vector it */ { int i,k; double norm; norm = 0.; for( i=0; icenter)[i] = 0.; for( k=0; kcenter)[i] += (((c->v)[k])->coord)[i]; /* N * barycenter */ norm += (c->center)[i] * (c->center)[i]; /* norm ^2 */ } /* norm --> 1 */ norm = sqrt(norm); for( i=0; icenter)[i] /= norm; } /* end of cone_center() */ /*-----------------------------------------------------------------*/ static int cone_params(CONE *c) /* compute parameters for hat for a cone */ { double coord[N]; /* coordinates of touching point */ double Tf,f; /* (transformed) density */ double Tgrad[N]; /* gradient of transformed density */ double Tderf; /* T'(f(x)) */ int i; /* aux variable */ #if DEBUG > 1000 double *g = c->g; /* pointer to vector g (g is stored in structure CONE) */ #else double g[N]; /* Vector g */ #endif /* coordinates of touching point */ for( i=0; itp * (c->center)[i]; #if MODE == 1 #if RECTANGLE == 1 /* check if point is in domain */ for( i=0; irl[i] || coord[i] > cr->ru[i] ) return TP_FZERO; #endif { /* move origin (0,0,...,0) into mode of density */ double mcoord[N]; for( i=0; imode[i]; /* density and its gradient */ f = (*cr->density)(mcoord); (*cr->grad_density)(mcoord,Tgrad); } #else /* mode = origin */ /* density and its gradient */ f = (*cr->density)(coord); (*cr->grad_density)(coord,Tgrad); #endif #if DEBUG > 1000 /* we store density in structure */ c->fp = f; #endif /* check density */ if( f < TOLERANCE ) /* f = 0. */ return TP_FZERO; /* transformed density and its gradient */ Tf = T(f); Tderf = T_deriv(f); for( i=0; ialpha = Tf - vec_scalarprod(Tgrad,coord,N); c->beta = vec_norm(Tgrad,N); /* |Tgrad| must not be to small */ if( c->beta < TOLERANCE ) return TP_GZERO; /* compute g and for all vertices of cone */ /* vector g = - grad(T(f)) / |grad(T(f))| */ for( i=0; ibeta; /* for each vertex v of cone and */ /* parameter a1 for volume of cone */ c->ai = c->detf; for( i=0; igv)[i] = vec_scalarprod(g,((c->v)[i])->coord,N ); /* */ if( (c->gv)[i] < TOLERANCE ) return TP_UGRAD; else c->ai /= (c->gv)[i]; } #if RECTANGLE == 1 /* at last calculate height of pyramid */ #if FIND_TP_FUNCTION == 1 get_height(c); /* this is expensive for calculation for every touching point !!! */ /* TODO: approximate gat_height with max_{vertices of rectangle} || vertex - mode || */ #endif #endif /* return error code */ return TP_VALID; } /* end of cone_params() */ /*-----------------------------------------------------------------*/ #if RECTANGLE == 1 /*-----------------------------------------------------------------*/ static void get_height(CONE *c) /* calculate height of pyramid (using simplex algorithm) */ { int i,j,row,ipc,ipr; double pc,pr,ratio; int bu,bl; double A[N+1][N+1]; /* set matrix for simplex algorithm */ for( i=0,row=0; iv[j])->coord[i] > 0 ) bu = 1; if( (c->v[j])->coord[i] < 0 ) bl = 1; } /* coefficients of inequalities */ if( bu ) { /* upper bound */ for( j=0; jv[j])->coord[i]; A[row][N] = cr->ru[i]; } row++; } if( bl ) { /* lower bound */ for( j=0; jv[j])->coord[i]; A[row][N] = - cr->rl[i]; } row++; } } /* objective function */ for( j=0; jgv[j]; A[row][N] = 0.; /* find maximum */ while( 1 ) { /* pivot column */ for( j=0,pc=0.,ipc=-1; jheight = A[row][N]; return; } /* find pivot row */ for( i=0,pr=-1.,ipr=-1; i ratio ) { ipr = i; pr = ratio; } } /* unbounded ? */ if( ipr == -1 ) { WARNING( "unbounded pyramid. This should not happen!" ); c->height = 1.e100; return; } /* make pivot step */ /* the rest */ for( i=0; i<=row; i++ ) if( i!= ipr ) for( j=0; j 1000 /*-----------------------------------------------------------------*/ static void cone_volume(CONE *c) /* volume of triangle of distance 1 */ { int i; c->Vi = c->detf; /* all vertices */ for( i=0; iVi /= vec_scalarprod( ((c->v)[i])->coord, c->center, N ); } /* end of cone_volume() */ /*-----------------------------------------------------------------*/ #endif /*-----------------------------------------------------------------*/ /******************************************************************* hat function ******************************************************************/ /*-----------------------------------------------------------------*/ static void get_hat_params(void) /* get parameters for hat function for each cone */ { CONE *c; /* aux pointer to cones */ double Hi_max; /* max of the Hi's */ double Hi_bound; /* reject Hi above this bound */ /* absolut bound given ? */ Hi_bound = ((cr->hc)->reject_tp_abs > 0.) ? (cr->hc)->reject_tp_abs : -1.; /* initialize */ cr->Htot = 0.; /* accumulated sum of volumes */ Hi_max = 0.; /* maximal Hi is zero at beginn */ /* all cones */ for( c=cr->c; c!=NULL; c=c->next ) { /* first check if calculation of cone parameters is necessary */ /* we only calculate cone parameters if c->level > cr->step_tp */ /* (i.e. when cone has been split after calculation of touching point). */ if( c->level > cr->step_tp ) { cone_center(c); /* center of cone */ if( cone_params(c) != TP_VALID ) c->tp = -1.; /* touching point not valid */ } /* calculate and check parameters for hat for cone (first run) */ hat_cone_params(c,Hi_bound); /* volume below hat */ cr->Htot += c->Hi; /* volume below hat */ c->Hsum = cr->Htot; /* accumulated sum of volumes */ Hi_max = MAX(Hi_max,c->Hi); /* maximal Hi */ if( c == cr->last_c ) break; } /* we need median and percentile for another check of the Hi's */ if( (cr->hc)->reject_tp_rel > 0. ) { /* check required */ /* compute median and percentile of Hi's */ get_median(); /* new bound */ Hi_bound = cr->Hi_median + (cr->hc)->reject_tp_rel * (cr->Hi_percentile - cr->Hi_median); /* checking all cones necessary ? */ if( Hi_max > Hi_bound ) { /* and now check all the cones again */ cr->Htot = 0.; Hi_max = 0.; for( c=cr->c; c!=NULL; c=c->next ) { /* all cones */ hat_cone_params(c,Hi_bound); cr->Htot += c->Hi; /* volume below hat */ c->Hsum = cr->Htot; /* accumulated sum of volumes */ if( c == cr->last_c ) break; } } } /* create guide tbale for finding cones */ make_guide_table(); } /* end of get_hat_params() */ /*-----------------------------------------------------------------*/ static void hat_cone_params(CONE *c, double Hi_bound) /* calculate and check parameters for the hat function for each cone */ { int i; while( 1 ) { /* calculate volume under hat and check its size */ while( c->tp < 0. ) { /* not a proper touching point */ /* we (must) split the cone again */ triangulate_cone(c,c->level+1); /* triangulate cone */ find_tp(c); /* find new touching point for first new cone */ find_tp(cr->last_c); /* find new touching point for second new cone */ #if DEBUG > 0 cr->steps_max = MAX( cr->steps_max, c->level ); #endif } /* calculate the volume */ #if RECTANGLE == 1 get_height(c); /* calculate height of pyramid */ c->Hi = exp(c->alpha) * c->ai * igammaf(N,c->beta*c->height); for( i=0; iHi /= c->beta; #else c->Hi = exp(c->alpha) * c->ai * fak[N-1]; for( i=0; iHi /= c->beta; #endif /* check if volume below hat is not too large */ if( Hi_bound > 0. /* whether to check Hi or not */ && Hi_bound < c->Hi ) /* Hi too large */ c->tp = -1.; else break; } #if DEBUG > 1000 if( DEBUG & DB_VOLUME ) cone_volume(c); /* volume of trinangle of distance 1 */ #endif } /* end of hat_cone_params() */ /*-----------------------------------------------------------------*/ static void get_median(void) /* compute median and percentile of Hi's */ { double *list_Hi; int i; CONE *c; /* alloc memory */ list_Hi = double_alloc(cr->max_c+1); /* we start at 1 because we use hpsort from numerical recipes */ /* copy Hi for all cones */ for( c=cr->c, i=1; c!=NULL; c=c->next, i++ ) { if( i > cr->max_c ) { /* this should not happen */ WARNING("number of cones exceeds cr->max_c (this should not happen)"); break; } list_Hi[i] = c->Hi; if( c == cr->last_c ) break; } if( i < cr->max_c ) /* this should not happen */ WARNING("number of cones less than cr->max_c (this should not happen)"); /* now sort this list */ hpsort((unsigned long)cr->max_c,list_Hi); /* median */ i = (int) (cr->max_c/2); cr->Hi_median = list_Hi[i]; /* percentile */ i = (int) ((cr->hc)->reject_tp_percentile * cr->max_c); i = MIN(i,cr->max_c-3); /* for the case that i >= cr->max_c */ cr->Hi_percentile = list_Hi[i]; /* we do not need this array any more */ free( list_Hi ); } /* end of get_median() */ /*-----------------------------------------------------------------*/ static void make_guide_table(void) /* create guide table */ { int j; CONE *c; /* memory for the guide table */ cr->size_guide = cr->max_c * (cr->hc)->guide_rel_size; cr->guide = guide_alloc(cr->size_guide); /* make table */ for( c=cr->c, j=0; c!=NULL && jsize_guide; j++ ) { while( c->Hsum / cr->Htot < (double) j / cr->size_guide ) { if( c == cr->last_c ) break; c=c->next; } (cr->guide)[j] = c; if( c == cr->last_c ) break; } /* is there an error ? */ if( jsize_guide ) /* this should not happen */ for( ; jsize_guide; j++ ) (cr->guide)[j] = cr->last_c; } /* end of make_guide_table() */ /*-----------------------------------------------------------------*/ /******************************************************************* optimal distance of touching points ******************************************************************/ /*-----------------------------------------------------------------*/ static void get_tp(void) /* compute optimal distance from origin of touching points for each cone */ { int k; CONE *c; for( k=0; kmax_c; k++ ) { /* set pointer to this cone */ c = cr->c+k; if( c->tp <= 0. ) /* find optimal distance for touching point */ find_tp(c); } } /* end of get_tp() */ /*-----------------------------------------------------------------*/ static void find_tp(CONE *c) /* find optimal touching point for cone Brent's method */ { TP_ARG a[3]; /* arguments for tp_min */ /* compute center */ cone_center(c); /* pointer to cone is needed for the function */ a->c = c; (a+1)->c = c; (a+2)->c = c; switch( find_bracket(a) ) { /* searching for intervall that contains minimum */ case TP_BRACKET: /* bracket found */ fminbr(a, (cr->hc)->find_tp_tol);/* use Brent's algorithm */ break; /* c->tp already set by tp_min() */ case TP_LEFT: /* minimum in left point */ c->tp = a->x; break; case TP_MIDDLE: /* minimum in middle point */ tp_min(a+1); break; case TP_RIGHT: /* minimum in right point */ c->tp = (a+2)->x; break; case TP_EMPTY: /* no proper touching point found */ default: c->tp = -1.; return; } } /* end of find_tp() */ /*-----------------------------------------------------------------*/ static int find_bracket(TP_ARG *a) /* search for proper bracket of minimum of function tp_min() returns 0 if no proper bracket could be found 1 otherwise a[0], a[1], a[2] ... left, middle and right point of intervall */ { int i; /* aux variable */ double xleft, xright; /* left boundary of searching region */ /* find proper touching point */ if( !find_proper_tp(a) ) /* no proper point found */ return TP_EMPTY; /** left point of intervall **/ /* initialize boundary of searching region and set starting point */ xleft = a->x; a->x = ((a+1)->x)/2.; /* search */ for( i=1; i <= MAX(1,(cr->hc)->find_tp_steps_lb); i++ ) { tp_min(a); if( !(a->valid) ) { /* a[0] not a proper touching point */ xleft = a->x; /* change boundary of searching region */ a->x += ( (a+1)->x - a->x )*0.5; /* try another one */ } else if( a->f <= (a+1)->f ) { /* a[0] is proper touching point, but ... */ (a+2)->x = (a+1)->x; (a+2)->f = (a+1)->f; (a+2)->valid = 1; (a+1)->x = a->x; (a+1)->f = a->f; (a+1)->valid = 1; a->x = xleft + (a->x - xleft)*0.5; } else /* all right: f(a[0]) > f(a[1]) */ break; } /* search successful ? */ if( !(a->valid) ) /* no proper touching point on left side */ return TP_MIDDLE; /* ????? */ if( a->f <= (a+1)->f ) /* f(left) <= f(middle) */ return TP_LEFT; /** right point of intervall **/ /* initialize a[2] if necessary */ if( (a+2)->x < 0. ) (a+2)->x = 1.1 * (a+1)->x; xright = -1.; /* no right boundary known yet */ xleft = (a+1)->x; /* search */ for( i=1; i <= MAX(1,(cr->hc)->find_tp_steps_ub); i++ ) { tp_min(a+2); if( !((a+2)->valid) ) { /* a[2] not a proper touching point */ xright = (a+2)->x; (a+2)->x = (xleft + (a+2)->x) * 0.5; /* try another one */ } else if( (a+2)->f <= (a+1)->f ) { /* move right */ xleft = (a+2)->x; (a+2)->x = (xright < 0.) ? (a+2)->x * 2. : (xright + (a+2)->x) * 0.5; } else /* all right: f(right) > f(middle) */ break; } /* search successful ? */ if( !((a+2)->valid) ) /* no proper touching point on right side */ return TP_MIDDLE; if( (a+2)->f <= (a+1)->f ) /* f(right) <= f(middle) */ return TP_RIGHT; /* we have found a bracket */ return TP_BRACKET; } /* end of find_bracket() */ /*-----------------------------------------------------------------*/ static int find_proper_tp(TP_ARG *a) /* search for proper touching point returns 0 if no proper point could be found 1 otherwise a[0], a[1], a[2] ... left, middle and right point of intervall the proper touching point is stored on a[1]. (see find_bracket() ) */ { int i; /* aux counter */ /* we store the proper touching point in a[1]. */ /* Thus a[0] and a[2] are no valid points */ a->valid = 0; (a+2)->valid = 0; /** --> infinity **/ #if RECTANGLE == 1 /* TODO: (a+1)->x has to be set to min of (cr->hc)->find_tp_start and c->height if known */ #endif /* initialize boundary of intervall */ a ->x = 0.; /* x[0] must >= 0. */ (a+1)->x = (cr->hc)->find_tp_start; /* starting point for searching proper touching point */ (a+2)->x = -1.; /* not known. marked by setting to -1. */ for( i=1; i <= MAX(1,(cr->hc)->find_tp_steps_r); i++ ) { /* TODO: if f==0 for a point we need not continue */ tp_min(a+1); /* calculate volume function */ if( !((a+1)->valid) ) { /* check validity of touching point */ /* not a proper touching point */ a->x = (a+1)->x; (a+1)->x *= (cr->hc)->find_tp_step_size; } else return 1; } /** --> 0 **/ /* initialize boundary of intervall */ a ->x = 0.; /* x[0] must >= 0. */ (a+1)->x = (cr->hc)->find_tp_start / (cr->hc)->find_tp_step_size;/* starting point for searching proper touching point */ (a+2)->x = (cr->hc)->find_tp_start; /* x[2] must >= x[1]. */ for( i=0; i <= MAX(0,(cr->hc)->find_tp_steps_l); i++ ) { /* TODO: if f==0 for a point we need not continue */ tp_min(a+1); /* calculate volume function */ if( !((a+1)->valid) ) { /* check validity of touching point */ /* not a proper touching point */ (a+2)->x = (a+1)->x; (a+1)->x /= (cr->hc)->find_tp_step_size; } else return 1; } /* no proper touching point found */ return 0; } /* end of find_proper_tp() */ /*-----------------------------------------------------------------*/ static void tp_min(TP_ARG *a) /* funtion that must be minimized for optimal touching point */ { (a->c)->tp = a->x; if( cone_params(a->c) != TP_VALID ) /* something is wrong: beta = 0 and/or <= 0 */ a->valid = 0; else { /*** this line depends on transformation T !! ***/ #if RECTANGLE == 1 #if FIND_TP_FUNCTION == 1 a->f = (a->c)->alpha - N * log((a->c)->beta) + log((a->c)->ai) + log(igammaf(N,(a->c)->beta*(a->c)->height)); #else a->f = (a->c)->alpha - N * log((a->c)->beta) + log((a->c)->ai); #endif #else a->f = (a->c)->alpha - N * log((a->c)->beta) + log((a->c)->ai); #endif /* check for numerical errors (alpha or beta too small) */ a->valid = (isnan(a->f)) ? 0 : 1; } if( !(a->valid) ) /* we mark this case by setting tp = -1 */ (a->c)->tp = -1.; } /* end of tp_min() */ /*-----------------------------------------------------------------*/ /******************************************************************* hash table for edges ******************************************************************/ /*-----------------------------------------------------------------*/ #if N > 2 /*-----------------------------------------------------------------*/ /* common variables */ static E_TABLE **etable; /* pointer to edge table */ static int etable_size; /* size of edge table */ #define hash(x,y) (3*((x)+(y))/2)%etable_size /* hash function */ /*-----------------------------------------------------------------*/ static void etable_new(int size) /* make new hash table */ { /* first clear hash table (if necessary) */ etable_clear(); /* set size of edge table */ etable_size = size; /* make root */ etable = p_etable_alloc(etable_size); } /* end if etable_new() */ /*-----------------------------------------------------------------*/ static void etable_clear(void) /* clear hash table */ { int i; if( etable == NULL ) /* nothing to do */ return; /* clear all branches */ for( i=0; inext != NULL ) etable_clear_next(etable[i]->next); free( etable[i] ); } /* clear root */ free( etable ); etable = NULL; } /* end if etable_clear() */ /*-----------------------------------------------------------------*/ static void etable_clear_next(E_TABLE *et) /* clear branches in edge table by recursion */ { if( et->next == NULL ) free( et ); else etable_clear_next(et->next); } /* end of etable_clear_next() */ /*-----------------------------------------------------------------*/ static VERTEX *etable_find_or_insert(VERTEX **vidx) /* find or insert entry in hash table it returns pointer to edge */ { E_TABLE **pet; /* pointer to table entry */ int idx[2]; /* store indices */ int hidx; /* hash number */ /* hash number */ idx[0] = vidx[0]->index; idx[1] = vidx[1]->index; hidx = hash(idx[0],idx[1]); /* check pointer */ if( etable == NULL ) FATAL( "pointer to edge table is NULL" ); /* get branch of hash table */ pet = etable + hidx; /* now find entry in branch */ while( *pet != NULL ) { if( (*pet)->index[0] == idx[0] && (*pet)->index[1] == idx[1] ) break; /* found ! */ pet = &((*pet)->next); } if( *pet == NULL ) { /* we have not found the index */ /* new entry in hash table */ *pet = etable_alloc(1); /* insert data of new edge */ /* indices of incident vertices */ (*pet)->index[0] = idx[0]; (*pet)->index[1] = idx[1]; /* compute new vertex */ (*pet)->vertex = new_vertex(vidx); } /* return pointer to (new) edge */ return (*pet)->vertex; } /* end of etable_find_or_insert() */ /*-----------------------------------------------------------------*/ static E_TABLE **p_etable_alloc(int n) /* allocate pointer to edge table */ { E_TABLE **pet; pet = (E_TABLE **)malloc((size_t) (n*sizeof(E_TABLE*))); if (!pet) FATAL("allocation failure in p_etable_alloc()"); /* initialize table */ for( n-- ; n>=0; n-- ) pet[n] = NULL; return pet; } /* end of p_etable_alloc() */ /*-----------------------------------------------------------------*/ static E_TABLE *etable_alloc(int n) /* allocate edge table entries */ { E_TABLE *et; et = (E_TABLE *)malloc((size_t) (n*sizeof(E_TABLE))); if (!et) FATAL("allocation failure in etable_alloc()"); /* initialize */ et->next = NULL; return et; } /* end of etable_alloc() */ /*-----------------------------------------------------------------*/ #endif /*-----------------------------------------------------------------*/ /******************************************************************* alloc and free memory for cones, vertices, guide tables ******************************************************************/ /*-----------------------------------------------------------------*/ static double *double_alloc(int n) /* allocate array of doubles */ { double *ad; ad = (double *)malloc((size_t) (n*sizeof(double))); if (!ad) FATAL("allocation failure in double_alloc()"); return ad; } /* end of double_alloc() */ /*-----------------------------------------------------------------*/ static C_CONSTR *hconstr_alloc(int n) /* allocate entries for parameters for constructing hat function */ { C_CONSTR *hc; hc = (C_CONSTR *)malloc((size_t) (n*sizeof(C_CONSTR))); if (!hc) FATAL("allocation failure in hconstr_alloc()"); return hc; } /* end of hconstr_alloc() */ /*-----------------------------------------------------------------*/ static C_ROOT *croot_alloc(int n) /* allocate entries for root of tree of cones */ { C_ROOT *cr; cr = (C_ROOT *)malloc((size_t) (n*sizeof(C_ROOT))); if (!cr) FATAL("allocation failure in croot_alloc()"); return cr; } /* end of croot_alloc() */ /*-----------------------------------------------------------------*/ static CONE *cone_alloc(int n) /* allocate cone entries */ { int i; CONE *c; if( cr->max_c + n > (cr->hc)->max_cones ) { db_hat(cr); db_close(cr); FATAL("maximum number of cones exceeded"); } c = (CONE *)malloc((size_t) (n*sizeof(CONE))); if (!c) FATAL("allocation failure in cone_alloc()"); /* make linked list */ for( i=0; inext = (c+i+1); (c+n-1)->next = NULL; return c; } /* end of cone_alloc() */ /*-----------------------------------------------------------------*/ static CONE *next_cone(void) /* next allocated cone in list */ { /* alloc memory if necessary */ if( (cr->last_c)->next == NULL ) (cr->last_c)->next = cone_alloc( (cr->hc)->cone_alloc_block ); /* move pointer to last entry */ cr->last_c = (cr->last_c)->next; /* and update counter */ ++(cr->max_c); /* return pointer to next vertex */ return cr->last_c; } /* end of next_cone() */ /*-----------------------------------------------------------------*/ static VERTEX *vertex_alloc(int n) /* allocate vertex entries */ { int i; VERTEX *v; v = (VERTEX *)malloc((size_t) (n*sizeof(VERTEX))); if (!v) FATAL("allocation failure in vertex_alloc()"); /* make linked list */ for( i=0; inext = (v+i+1); (v+n-1)->next = NULL; return v; } /* end of vertex_alloc() */ /*-----------------------------------------------------------------*/ static VERTEX *next_vertex(void) /* next allocated vertex in list */ { /* alloc memory if necessary */ if( (cr->last_v)->next == NULL ) (cr->last_v)->next = vertex_alloc( (cr->hc)->vertex_alloc_block ); /* move pointer to last entry */ cr->last_v = (cr->last_v)->next; /* index of vertex */ (cr->last_v)->index = cr->max_v; /* and update counter */ ++(cr->max_v); /* return pointer to next vertex */ return cr->last_v; } /* end of next_vertex() */ /*-----------------------------------------------------------------*/ static CONE **guide_alloc(int n) /* allocate pointer to array of pointers to cones and initialize */ { CONE **ca; ca = (CONE **)malloc((size_t) (n*sizeof(CONE*))); if (!ca) FATAL("allocation failure in guide_alloc()"); /* initialize table */ for( n-- ; n>=0; n-- ) ca[n] = NULL; return ca; } /* end of guide_alloc() */ /*-----------------------------------------------------------------*/ /******************************************************************* misc ******************************************************************/ /*-----------------------------------------------------------------*/ static double vec_norm(double *x, int n) /* norm of vector */ { double norm = 0.; for( --n; n>=0; n-- ) norm += x[n] * x[n]; return sqrt(norm); } /* end of vec_norm() */ /*-----------------------------------------------------------------*/ static void sort_varray(VERTEX **va) /* sort vertex array (insertion sort) */ { int i,j; VERTEX *v; for( i=1; i0 && va[j-1]->index > v->index ) { va[j] = va[j-1]; j--; } va[j] = v; } } /* end of sort_iarray() */ /*-----------------------------------------------------------------*/ static int ipower(int n, int k) /* power n^k of integers */ { int power = 1; for( ; k>0; k-- ) power *= n; return power; } /* end of int_power() */ /*-----------------------------------------------------------------*/ #if RECTANGLE == 1 /*-----------------------------------------------------------------*/ static double igammaf(int a, double x) /* incomplete gamma function * gamma(a) for integer a, a>=2 */ { int k; double sum,potenz; sum = 1.; potenz = 1.; for( k=1; k=0; n-- ) s += x[n]*y[n]; return s; } /* end of vec_scalarprod */ /*-----------------------------------------------------------------*/ #if N > 2 /*-----------------------------------------------------------------*/ static void copy_varray(VERTEX **a, VERTEX **b, int n) /* copy first integer array of n elements into second */ { int i; for( i=0; ic; /* parameters for volume function */ a = ap->x; b = (ap+2)->x; /* set border of range */ if( tol<=0. || b<=a ) /* check parameters */ FATAL( "tol > 0 && b > a" ); v = a + r*(b-a); /* First step - always gold section*/ ar.x = v; tp_min(&ar); fv = ar.f; x = v; w = v; fx=fv; fw=fv; for(i=1; i<=(cr->hc)->find_tp_steps_brent; i++) /* Main iteration loop */ { double range = b-a; /* Range over which the minimum */ /* is seeked for */ double middle_range = (a+b)/2; double tol_act = /* Actual tolerance */ SQRT_EPSILON*fabs(x) + tol/3; double new_step; /* Step at this iteration */ if( fabs(x-middle_range) + range/2 <= 2*tol_act ) return x; /* Acceptable approx. is found */ /* Obtain the gold section step */ new_step = r * ( x= tol_act ) /* If x and w are distinct */ { /* interpolatiom may be tried */ register double p; /* Interpolation step is calcula-*/ register double q; /* ted as p/q; division operation*/ /* is delayed until last moment */ register double t; t = (x-w) * (fx-fv); q = (x-v) * (fx-fw); p = (x-v)*q - (x-w)*t; q = 2*(q-t); if( q>(double)0 ) /* q was calculated with the op-*/ p = -p; /* posite sign; make q positive */ else /* and assign possible minus to */ q = -q; /* p */ if( fabs(p) < fabs(new_step*q) && /* If x+p/q falls in [a,b]*/ p > q*(a-x+2*tol_act) && /* not too close to a and */ p < q*(b-x-2*tol_act) ) /* b, and isn't too large */ new_step = p/q; /* it is accepted */ /* If p/q is too large then the */ /* gold section procedure can */ /* reduce [a,b] range to more */ /* extent */ } if( fabs(new_step) < tol_act ) /* Adjust the step to be not less*/ if( new_step > (double)0 ) /* than tolerance */ new_step = tol_act; else new_step = -tol_act; /* Obtain the next approximation to min */ { /* and reduce the enveloping range */ register double t = x + new_step; /* Tentative point for the min */ register double ft; ar.x = t; tp_min(&ar); ft = ar.f; if( ft <= fx ) { /* t is a better approximation */ if( t < x ) /* Reduce the range so that */ b = x; /* t would fall within it */ else a = x; v = w; w = x; x = t; /* Assign the best approx to x */ fv=fw; fw=fx; fx=ft; } else /* x remains the better approx */ { if( t < x ) /* Reduce the range enclosing x */ a = t; else b = t; if( ft <= fw || w==x ) { v = w; w = t; fv=fw; fw=ft; } else if( ft<=fv || v==x || v==w ) { v = t; fv=ft; } } } /* ----- end-of-block ----- */ } /* ===== End of loop ===== */ /* number of allowed iteration steps exceeded */ return x; } /* end of fminbr() */ #undef SQRT_EPSILON /*-----------------------------------------------------------------*/ /******************************************************************* routines from numerical recipes in C changes by Josef Leydold ******************************************************************/ /*-----------------------------------------------------------------*/ static void hpsort(unsigned long n, double *ra) /* heap sort array ra[1...n] (!!!!) */ { unsigned long i,ir,j,l; double rra; if (n < 2) return; l=(n >> 1)+1; ir=n; for (;;) { if (l > 1) { rra=ra[--l]; } else { rra=ra[ir]; ra[ir]=ra[1]; if (--ir == 1) { ra[1]=rra; break; } } i=l; j=l+l; while (j <= ir) { if (j < ir && ra[j] < ra[j+1]) j++; if (rra < ra[j]) { ra[i]=ra[j]; i=j; j <<= 1; } else j=ir+1; } ra[i]=rra; } } /* (C) Copr. 1986-92 Numerical Recipes Software 75s@Ww+. */ /*-----------------------------------------------------------------*/