#! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # C/ # C/Dp/ # C/Dp/Drivers/ # C/Dp/Drivers/driver1.c # C/Dp/Drivers/driver2.c # C/Dp/Drivers/driver3.c # C/Dp/Drivers/res1 # C/Dp/Drivers/res2 # C/Dp/Drivers/res3 # C/Dp/Src/ # C/Dp/Src/alc2d.h # C/Dp/Src/src.c # Doc/ # Doc/readme # This archive created: Tue Feb 26 12:26:08 2002 export PATH; PATH=/bin:$PATH if test ! -d 'C' then mkdir 'C' fi cd 'C' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'driver1.c' then echo shar: will not over-write existing file "'driver1.c'" else cat << "SHAR_EOF" > 'driver1.c' /* ALC2DLIB-Version 1.0 implemented by WH 22.3.99, all rights reserved */ /*please report problems or bugs to whoer@statistik.wu-wien.ac.at */ /* or hormannw@boun.edu.tr */ /* This file: main1.c Demonstrates, how to call functions of "alc2d.c" to generate random pairs of two-dimensional log-concave distributions To compile this example on a unix-computer just type cc alc2d.c main1.c -lm */ #include #include #include "alc2d.h" /****************************************************************************/ /* Include your random number generator; the generator here is only suitable for small simulations!!!! */ static unsigned long int urn=0;/*seed of the uniform generator*/ double unif()/*uniform r n Marsaglia m=2^32 a=69069 c=1*/ /*uniform random number generator*/ { urn=(69069*(urn)+1); return(urn/4.294967296e9); } /*******************************************************************/ /*code the logarithm of the density function you want to sapmle from: as this is the first example we take the standard normal distribution; allthough we do not use parameters, we must icnlude them in the prototype as the generation algorithm is expecting them as well*/ /*log of bivariate standard normal density*/ double hstandardnormal(double x[2],double par[]) {return -0.5*(x[0]*x[0]+x[1]*x[1]);} /*partial derivative with respect to x[0]*/ double hnx(double x[2],double par[]) {return(-x[0]);} /*partial derivative with respect to x[1]*/ double hny(double x[2],double par[]) {return(-x[1]);} /*******************************************************************/ int main() { long int i; double x[2],sp[1][2],eq[4][3]; void *t1=NULL; /*first we have to choose a starting-value close to the mode, e.g.:*/ sp[0][0]=0.1; sp[0][1]=0.2; /* As the normal distribution is defined on the whole R^2, we need not define the domain of the distribution, but as we have chosen only one starting-value we have to use a bounded auxiliary domain: In this example we use the square with the corners (-1,-1), (-1,1), (1,-1), and (1,1), as auxiliary domain. It is bordered by the four lines 0=1+x; 0=1+y; 0=-1+x; 0=-1+y lines are stored in the form: 0 = eq[0][0] + eq[0][1]*x + eq[0][2]*y ;*/ eq[0][0]=1.; eq[0][1]=1.;eq[0][2]=0.; eq[1][0]=1.; eq[1][1]=0.;eq[1][2]=1.; eq[2][0]=-1.; eq[2][1]=1.;eq[2][2]=0.; eq[3][0]=-1.; eq[3][1]=0.;eq[3][2]=1.; t1=setup(1,20,sp,0,eq,4,hstandardnormal,hnx,hny,NULL,0); /* 1 ... we have one starting point 20... the algorithm can use up to 20 design points sp... array containing the starting points 0 ... no of equalities for the domain (as it is the whole R^2) eq... array containing the equalities for domain and auxiliary domain 4 ... no of equalities for the auxiliary domain hstandardnormal ... logarithm of the density function hnx... partial derivative with respect to x hny... partial derivative with respect to y NULL.. no array of parameters is needed 0 ... no of parameters of the distribution*/ /* now we can sample from the distribution, the result is stored in x[]*/ for(i=0;i<1000;i++) sample2d(x,t1); printf("The last generated pair was the pair (%f|%f)\n",x[0],x[1]); /* after generation the allocated memory must be freed*/ freesetup(t1); exit(0); } SHAR_EOF fi # end of overwriting check if test -f 'driver2.c' then echo shar: will not over-write existing file "'driver2.c'" else cat << "SHAR_EOF" > 'driver2.c' /* ALC2DLIB-Version 1.0 implemented by WH 22.3.99, all rights reserved */ /*please report problems or bugs to whoer@statistik.wu-wien.ac.at */ /* or hormannw@boun.edu.tr */ /* This file: main2.c Demonstrates, how to call functions of "alc2d.c" to generate random pairs of several two-dimensional log-concave distributions To compile this example on a unix-computer just type cc alc2d.c main2.c -lm */ #include #include #include "alc2d.h" #define MIN(a,b) ((a)<(b)?a:b) /****************************************************************************/ static unsigned long int urn=0;/*seed of the uniform generator*/ double unif()/*uniform r n Marsaglia m=2^32 a=69069 c=1*/ /*uniform random number generator*/ { urn=(69069*(urn)+1); return(urn/4.294967296e9); } /*******************************************************************/ /*Two-dimensional normal distribution*/ /*par[0]..R, par[1]..S1, par[2]..S2*/ double hn(double x[2],double par[]) /*log of bivariate normal density*/ {return (-0.5/(1.-par[0]*par[0])*(x[0]*x[0]*(1./(par[1]*par[1]))- 2.*par[0]*x[0]*x[1]*(1./(par[1]*par[2]))+x[1]*x[1]*(1./(par[2]*par[2]))));} double hnx(double x[2],double par[]) /*partial derivative w.r.t. x*/ {return(-0.5/(1.-par[0]*par[0])* (2*x[0]*(1./(par[1]*par[1]))-2.*par[0]*x[1]*(1./(par[1]*par[2]))));} double hny(double x[2],double par[]) /*partial derivative w.r.t. x*/ { return(-0.5/(1.-par[0]*par[0])* (-2.*par[0]*x[0]*(1./(par[1]*par[2]))+2*x[1]*(1./(par[2]*par[2]))));} /*******************************************************************/ /*"cut normal distribution" of the paper*/ /*log of normal distribution "cut through" by a plane by taking the minimum of the paraboloid and the plane*/ /*par[0]..R, par[1]..S1, par[2]..S2, par[3]...NCA, par[4]...NCB, par[5]...NCC*/ double hnc(double x[2],double par[]) /*log of bivariate cut normal density*/ { double fn; fn=(-0.5/(1.-par[0]*par[0])*(x[0]*x[0]*(1./(par[1]*par[1]))- 2.*par[0]*x[0]*x[1]*(1./(par[1]*par[2]))+ x[1]*x[1]*(1./(par[2]*par[2])))); return MIN(par[5]+par[3]*x[0]+par[4]*x[1],fn); } double hncx(double x[2],double par[]) /*partial derivative w.r.t. x*/ { double fn,fe; fn=(-0.5/(1.-par[0]*par[0])*(x[0]*x[0]*(1./(par[1]*par[1]))- 2.*par[0]*x[0]*x[1]*(1./(par[1]*par[2]))+x[1]*x[1]*(1./(par[2]*par[2])))); fe=par[5]+par[3]*x[0]+par[4]*x[1]; if(fn(par[0])) return(-par[0]*x[0]/r); else return 0.; } double hrefy(double x[],double par[]) /*partial derivative w.r.t. y*/ { double r; r=sqrt(x[0]*x[0]+x[1]*x[1]); if(r>(par[0])) return(-par[0]*x[1]/r); else return 0.; } /*******************************************************************/ /*distribution NS3 of the paper par[0]..a, par[1]..b, par[2]..c, par[3]..d, par[4]..e, if a,b,d,e>0 and 4bd>=c^2 the distribution is log concave on R^2*/ double hq(double a[2],double par[]) /*log of bivariate quartic density*/ { double x2,y2,x,y; x=a[0]; y=a[1]; x2=x*x; y2=y*y; return(-(par[0]*x2*x2+par[1]*x2+par[2]*x*y+par[3]*y2+par[4]*y2*y2)); } double hqx(double a[2],double par[]) /*partial derivative w.r.t. x*/ { double x,y; x=a[0]; y=a[1]; return(-(4*par[0]*x*x*x+2*par[1]*x+par[2]*y)); } double hqy(double a[2],double par[]) /*partial derivative w.r.t. y*/ { double x,y; x=a[0]; y=a[1]; return(-(par[2]*x+2*par[3]*y+4*par[4]*y*y*y)); } /*******************************************************************/ /*bivariate beta-distribution*/ /*par[0]..a1, par[1]..a2, par[2]..a3*/ double hb(double x[2],double par[]) /*log of density of bivariate beta distribution domain of the distribution is the triangle (0,0), (0,1) (1,0)*/ {return (par[0]-1.)*log(x[0])+(par[1]-1.)*log(x[1])+(par[2]-1.)*log(1.-x[0]-x[1]);} double hbx(double x[2],double par[]) /*partial derivative w.r.t. x*/ {return (par[0]-1.)/x[0]-(par[2]-1.)/(1.-x[0]-x[1]);} double hby(double x[2],double par[]) /*Partial derivative w.r.t. y*/ {return (par[1]-1.)/x[1]-(par[2]-1.)/(1.-x[0]-x[1]);} /*function can be used to initialize the domain of the beta-distribution*/ void betainit(double eq[][3]) { /*line x=0*/ eq[0][0]=0.; eq[0][1]=1.; eq[0][2]=0.; /* line y=0 */ eq[1][0]=0.; eq[1][1]=0.; eq[1][2]=1.; /* line x+y=1 */ eq[2][0]=-1.; eq[2][1]=1.; eq[2][2]=1.; } /*******************************************************************/ /*bivariate gamma-distribution Becker and Roux 1981: South African Statistical Journal 15, 1-12*/ /*Caution!!! This distribution is not log-concave for all choices*/ /* of the parameters, allthough it is allways logconcave for a1, a2>1 for x< y and y=0*/ { if(x[0]<=x[1]) return (par[0]-1.)*log(x[0])+(par[1]-1.)*log(par[3]*(x[1]-x[0])+x[0])- (1./par[4]+1./par[5]-par[3]/par[5])*x[0]-par[3]*x[1]/par[5]; else return (par[1]-1.)*log(x[1])+(par[0]-1.)*log(par[2]*(x[0]-x[1])+x[1])- (1./par[4]+1./par[5]-par[2]/par[4])*x[1]-par[2]*x[0]/par[4]; } double hgamx(double x[2],double par[]) /*partial derivative w.r.t. x*/ { if(x[0]<=x[1]) return (par[0]-1.)/x[0]+(par[1]-1.)*(1.-par[3])/(par[3]*(x[1]-x[0])+x[0]) -(1./par[4]+1./par[5]-par[3]/par[5]); else return (par[0]-1.)*par[2]/(par[2]*(x[0]-x[1])+x[1])-par[2]/par[4]; } double hgamy(double x[2],double par[]) /*partial derivative w.r.t. y*/ { if(x[0]<=x[1]) return (par[1]-1.)*par[3]/(par[3]*(x[1]-x[0])+x[0])-par[3]/par[5]; else return (par[1]-1.)/x[1]+(par[0]-1.)*(1.-par[2])/(par[2]*(x[0]-x[1])+x[1]) -(1./par[4]+1./par[5]-par[2]/par[4]); } /*function can be used to initialize the domain of the gamma-distribution*/ void gammainit(double eq[][3]) { /*line x=0*/ eq[0][0]=0.; eq[0][1]=1.; eq[0][2]=0.; /* line y=0 */ eq[1][0]=0.; eq[1][1]=0.; eq[1][2]=1.; } /*******************************************************************/ /* distribution NS1 of the paper */ double hg(double x[2],double par[]) /*log of density of a simple 2-dimensional gamma-distribution domain is the half-plain with x>0 */ {return log(x[0])-x[0]*x[0]-x[0]*x[1]-x[1]*x[1];} double hgx(double x[2],double par[]) /*Partielle Ableitung hb nach x*/ {return 1./x[0]-2.*x[0]-x[1];} double hgy(double x[2],double par[]) /*Partielle Ableitung hb nach y*/ {return -x[0]-2.*x[1];} /*******************************************************************/ #define HD 4 /*Half length of the side of the square used as auxiliary domain*/ int main() { long int i,wid; double x[8][2],sp[1][2],eq[10][3],par[10],pair[2]; void *t1=NULL,*t2=NULL,*t3=NULL,*t4=NULL,*t5=NULL,*t6=NULL,*t7=NULL,*t8=NULL; /*auxiliary domain is the square (-HD,HD)x(-HD,HD) */ eq[0][0]=HD; eq[0][1]=1.; eq[0][2]=0.; eq[1][0]=HD; eq[1][1]=0.; eq[1][2]=1.; eq[2][0]=-HD; eq[2][1]=1.; eq[2][2]=0.; eq[3][0]=-HD; eq[3][1]=0.; eq[3][2]=1.; /*We take a point close the mode as starting point*/ sp[0][0]=0.1; sp[0][1]=0.2122; /*setting the parameters of the normal distribution*/ par[0]=0.5;/*R*/ par[1]=12.;/*S1*/ par[2]=5.;/*S2*/ printf("Setup for the normal distribution with r=%f, s1=%f and s2=%f\n",par[0],par[1],par[2]); t1=setup(1,50,sp,0,eq,4,hn,hnx,hny,par,3); /*t1 is now pointer to the structure that holds all information for sampling from the normal distribution with R=0.5, S1=12 and S2=5 */ /*As second example we generate from a bivariate normal distribution restricted to the domain x>=0 an y>=-x For the definition of the domain we use*/ eq[0][0]=0.; eq[0][1]=1.; eq[0][2]=0.; eq[1][0]=0.; eq[1][1]=1.; eq[1][2]=1.; /*for the auxiliary domain we add the two lines x=3 y=3 */ eq[2][0]=-3.; eq[2][1]=1.; eq[2][2]=0.; eq[3][0]=-3.; eq[3][1]=0.; eq[3][2]=1.; /*We take as starting point*/ sp[0][0]=1.; sp[0][1]=0.; par[0]=0.95;/*R*/ par[1]=10.;/*S1*/ par[2]=10.;/*S2*/ printf("Setup for the normal distribution with r=%f, s1=%f and s2=%f\n",par[0],par[1],par[2]); printf("over the domain : x>=0 and y>-x\n"); t2=setup(1,20,sp,2,eq,4,hn,hnx,hny,par,3); /* (1,20,sp... shows, that we use one starting-point and that we allow 20 design points; ...2,eq,4.. means, that there are 2 equalities for the domain of the distribution and 2 additional equalities for an auxiliary domain t2 is now pointer to the structure that holds all information for sampling from the normal distribution with R=0.95, S1=10 and S2=10 over the restricted domain*/ /*the NS2 distribution is again generated over the whole R2, so we use the same auxiliary domain, as for the unrestricted normal distribution */ /*auxiliary domain is the square (-HD,HD)x(-HD,HD) */ eq[0][0]=HD; eq[0][1]=1.; eq[0][2]=0.; eq[1][0]=HD; eq[1][1]=0.; eq[1][2]=1.; eq[2][0]=-HD; eq[2][1]=1.; eq[2][2]=0.; eq[3][0]=-HD; eq[3][1]=0.; eq[3][2]=1.; /*We take a point close to the mode as starting point*/ sp[0][0]=0.01; sp[0][1]=0.1; par[0]=3.;/*n*/ printf("set-up for distribution NS2 with N=%f\n",par[0]); t3=setup(1,50,sp,0,eq,4,href,hrefx,hrefy,par,1); par[0]=0.05;/*R*/ par[1]=1.;/*S1*/ par[2]=1.;/*S2*/ par[3]=0.5;/*A*/ par[4]=1.;/*B*/ par[5]= -1.;/*C*/ printf("Setup for cut normal-distr.: R=%f S1=%f S2=%f\n",par[0],par[1],par[2]); printf("A %f B %f C %f \n",par[3],par[4],par[5]); t4=setup(1,50,sp,0,eq,4,hnc,hncx,hncy,par,6); par[0]=0.03;/*A*/ par[1]=0.1;/*B*/ par[2]=1.;/*C*/ par[3]=0.5;/*D*/ par[4]=0.01;/*E*/ /*checks concavity*/ if(4*par[1]*par[3]=1 in alc2d.h */ wid=50000; /* First we have to set the repetition counter to 0 */ setwidcount(0); printf("We are generating %ld pairs from the normal distribution\n",wid); for(i=0;i 'driver3.c' /* ALC2DLIB-Version 1.0 implemented by WH 22.3.99, all rights reserved */ /*please report problems or bugs to whoer@statistik.wu-wien.ac.at */ /* or hormannw@boun.edu.tr */ /*This file: mainnort.c */ #define PORTABLE 1 /*0..erf() is included in math-library, 1..erf() is not included*/ #define WID 100000 #define KLAS 30 #define START 17899234 /* starting value for uniform generator */ #define R 0.8 #define S1 1. #define S2 1. #include #include #include "alc2d.h" #ifndef PI #define PI 3.14159265358979323846 #endif #define MIN(a,b) ((a)<(b)?a:b) /****************************************************************************/ static unsigned long int urn=0;/*seed of the uniform generator*/ double unif()/*uniform r n Marsaglia m=2^32 a=69069 c=1*/ /*uniform random number generator*/ { urn=(69069*(urn)+1); return(urn/4.294967296e9); } /*******************************************************************/ /*par[0]..R, par[1]..S1, par[2]..S2*/ double hn(double x[2],double par[]) /*log of bivariate normal density*/ {return (-0.5/(1.-par[0]*par[0])*(x[0]*x[0]*(1./(par[1]*par[1]))- 2.*par[0]*x[0]*x[1]*(1./(par[1]*par[2]))+x[1]*x[1]*(1./(par[2]*par[2]))));} double hnx(double x[2],double par[]) /*Partielle Ableitung hn nach x*/ {return(-0.5/(1.-par[0]*par[0])* (2*x[0]*(1./(par[1]*par[1]))-2.*par[0]*x[1]*(1./(par[1]*par[2]))));} double hny(double x[2],double par[]) /*Partielle Ableitung hn nach y*/ { return(-0.5/(1.-par[0]*par[0])* (-2.*par[0]*x[0]*(1./(par[1]*par[2]))+2*x[1]*(1./(par[2]*par[2]))));} /*******************************************************************/ /*******************************************************************/ #if PORTABLE==0 double dfnormal(x) /*distribution function normal */ double x; { if (x>=0) return((1+erf(x/M_SQRT2))/2); else return(erfc(-x/M_SQRT2)/2); } #else /****************************************************************/ /* Verteilungsfunktion von N(0,1) programmiert nach */ /* Kennedy/Gentle "Statistical Computing" S 90 - 92. */ /****************************************************************/ double dfnormal(x) double x; { int help; double xx,x2,x3,x4,x5,x6,x7,x8,zahler,nenner; static double ep0=242.66795523053175,ep1=21.979261618294152, ep2=6.9963834886191355,ep3= -3.5609843701815385e-2, eq0=215.0588758698612,eq1=91.164905404514901, eq2=15.082797630407787,/*eq3=1,*/ zp0=300.4592610201616005,zq0=300.4592609569832933, zp1=451.9189537118729422,zq1=790.9509253278980272, zp2=339.3208167343436870,zq2=931.3540948506096211, zp3=152.9892850469404039,zq3=638.9802644656311665, zp4=43.16222722205673530,zq4=277.5854447439876434, zp5=7.211758250883093659,zq5=77.00015293522947295, zp6=0.5641955174789739711,zq6=12.78272731962942351, zp7= -1.368648573827167067e-7,/*zq7=1,*/ dp0= -2.99610707703542174e-3,dq0=1.06209230528467918e-2, dp1= -4.94730910623250734e-2,dq1=1.91308926107829841e-1, dp2= -2.26956593539686930e-1,dq2=1.05167510706793207, dp3= -2.78661308609647788e-1,dq3=1.98733201817135256, dp4= -2.23192459734184686e-2/*,dq4=1*/; if(x<0) help=0; else help=1; x=fabs(x)/sqrt(2.); if(x<0.5) { x2=x*x; x4=x2*x2; x6=x4*x2; zahler=ep0+ep1*x2+ep2*x4+ep3*x6; nenner=eq0+eq1*x2+eq2*x4+x6; if(help) return(0.5*(1+x*zahler/nenner)); else return(0.5*(1-x*zahler/nenner)); } else if(x<4) { x2=x*x; x3=x2*x; x4=x3*x; x5=x4*x; x6=x5*x; x7=x6*x; zahler=zp0+zp1*x+zp2*x2+zp3*x3+zp4*x4+zp5*x5+zp6*x6+zp7*x7; nenner=zq0+zq1*x+zq2*x2+zq3*x3+zq4*x4+zq5*x5+zq6*x6+x7; if(help) return(0.5*(2-exp(-x2)*zahler/nenner)); else return(0.5*exp(-x2)*zahler/nenner); } else if(x<50) { xx=x*x; x2=1/xx; x4=x2*x2; x6=x4*x2; x8=x6*x2; zahler=dp0+dp1*x2+dp2*x4+dp3*x6+dp4*x8; nenner=dq0+dq1*x2+dq2*x4+dq3*x6+x8; if(help) return(1-0.5*(exp(-xx)/x)* (1/sqrt(PI)+zahler/(xx*nenner))); else return(0.5*(exp(-xx)/x)*(1/sqrt(PI)+zahler/(xx*nenner))); } else if(help) return(1.); else return(0.); } #endif /**************************************************************** * cumulated distribution function - CHI^2 distribution * * for n degrees of freedom, * * Approximation from Wilson & Hilferty * * Kendall Stuart I, p 371 * ****************************************************************/ double dfchia(x,n) /* approximation */ double x; long n; { double chix,y; y=2./(9.*n); chix=(exp(log(x/n)/3)-1+y)/sqrt(y); return(dfnormal(chix)); } /****************************************/ /*Chi-2 test for uniform distribution*/ double chi2test(long b[]/*obsereved frequencies*/, long l,/*number of classes*/ long wid,/*sample size*/ int printflag/*0..no output, 1.. little, 2..more*/) { long int i; double chi2=0.,erw,pval; erw=(double)wid/l; if (erw<5.) printf("Error chi2test: expected frequency smaller than 5\n"); for (i=0;i=1) { printf(" Chi2-test: samplesize=%ld number of classes= %ld \n",wid,l); printf("Chi2-value %f Approximate P-Value: %f\n",chi2,pval); } return(chi2); } /*******************************************************************/ long int bx[KLAS]; long int by[KLAS]; long int bw[KLAS]; long int bg[KLAS][KLAS]; #define HD 3. int main() { long int j,ix,iy; double x[2],z[2],sp[5][2],ineq[10][3],par[10],w; /*array of the inequalities defining the domain */ /*inequalities are stored as 0 <= ineq[0][0] + ineq[0][1]*x + ineq[0][2]*y */ void *t1=NULL; ineq[0][0]=HD; ineq[0][1]=1.; ineq[0][2]=0.; ineq[1][0]=HD; ineq[1][1]=0.; ineq[1][2]=1.; ineq[2][0]=-HD; ineq[2][1]=1.; ineq[2][2]=0.; ineq[3][0]=-HD; ineq[3][1]=0.; ineq[3][2]=1.; /*sets the seed for the uniform generator*/ urn=START; sp[0][0]=0.1; sp[0][1]=0.2; par[0]=R;par[1]=S1;par[2]=S2; printf("\nsetup and sample of size %d for normal-distr.: R=%f S1=%f S2=%f\n",WID,par[0],par[1],par[2]); t1=setup(1,20,sp,0,ineq,4,hn,hnx,hny,par,3); for(j=0;j 'res1' as volume unbounded setup is restarted using the auxiliary domain points 1, volume below hat 4.135541e+00,(*t).gennr 4 points 2, volume below hat 3.684925e+00,(*t).gennr 8 points 3, volume below hat 3.550357e+00,(*t).gennr 12 points 3, volume below hat 2.732331e+02,(*t).gennr 3 Starting values found, auxiliary domain points dismissed there are 8 allocated pointers with 5400 bytes memory, maximal heapsize up to now 11536 bytes (adr.of first p. 2ca18) points 4, volume below hat 7.927591e+01,(*t).gennr 6 points 5, volume below hat 7.342767e+01,(*t).gennr 16 points 6, volume below hat 3.712564e+01,(*t).gennr 20 points 7, volume below hat 2.049380e+01,(*t).gennr 28 points 8, volume below hat 1.643235e+01,(*t).gennr 36 points 9, volume below hat 1.243682e+01,(*t).gennr 40 points 10, volume below hat 1.162194e+01,(*t).gennr 45 points 12, volume below hat 1.009872e+01,(*t).gennr 57 points 14, volume below hat 8.220588e+00,(*t).gennr 73 points 16, volume below hat 7.744880e+00,(*t).gennr 89 points 18, volume below hat 7.599256e+00,(*t).gennr 105 points 20, volume below hat 7.458044e+00,(*t).gennr 121 Polygon structure was freed The last generated pair was the pair (0.336184|-1.127603) there are 0 allocated pointers with 0 bytes memory, maximal heapsize up to now 19704 bytes (adr.of first p. 0) SHAR_EOF fi # end of overwriting check if test -f 'res2' then echo shar: will not over-write existing file "'res2'" else cat << "SHAR_EOF" > 'res2' Setup for the normal distribution with r=0.500000, s1=12.000000 and s2=5.000000 as volume unbounded setup is restarted using the auxiliary domain points 1, volume below hat 6.408290e+01,(*t).gennr 4 points 2, volume below hat 6.105792e+01,(*t).gennr 8 points 3, volume below hat 5.709573e+01,(*t).gennr 12 points 4, volume below hat 5.690028e+01,(*t).gennr 16 points 5, volume below hat 5.661521e+01,(*t).gennr 22 points 6, volume below hat 5.654689e+01,(*t).gennr 28 points 7, volume below hat 5.617734e+01,(*t).gennr 32 points 7, volume below hat 2.039911e+03,(*t).gennr 24 Starting values found, auxiliary domain points dismissed there are 13 allocated pointers with 14752 bytes memory, maximal heapsize up to now 29472 bytes (adr.of first p. 32260) Setup for the normal distribution with r=0.950000, s1=10.000000 and s2=10.000000 over the domain : x>=0 and y>-x as volume unbounded setup is restarted using the auxiliary domain points 1, volume below hat 1.299823e+01,(*t).gennr 4 points 2, volume below hat 1.212154e+01,(*t).gennr 8 points 2, volume below hat 1.007015e+04,(*t).gennr 3 Starting values found, auxiliary domain points dismissed there are 21 allocated pointers with 20016 bytes memory, maximal heapsize up to now 29472 bytes (adr.of first p. 32260) set-up for distribution NS2 with N=3.000000 as volume unbounded setup is restarted using the auxiliary domain points 1, volume below hat 6.400000e+01,(*t).gennr 2 points 2, volume below hat 5.872184e+01,(*t).gennr 5 points 3, volume below hat 5.355078e+01,(*t).gennr 8 points 4, volume below hat 4.813837e+01,(*t).gennr 11 points 5, volume below hat 4.520246e+01,(*t).gennr 19 points 5, volume below hat 5.039485e+01,(*t).gennr 8 Starting values found, auxiliary domain points dismissed there are 32 allocated pointers with 33088 bytes memory, maximal heapsize up to now 46608 bytes (adr.of first p. 32260) Setup for cut normal-distr.: R=0.050000 S1=1.000000 S2=1.000000 A 0.500000 B 1.000000 C -1.000000 as volume unbounded setup is restarted using the auxiliary domain points 1, volume below hat 2.912920e+02,(*t).gennr 4 points 2, volume below hat 6.386800e+01,(*t).gennr 8 points 3, volume below hat 9.758423e+00,(*t).gennr 14 points 3, volume below hat 1.257848e+01,(*t).gennr 3 Starting values found, auxiliary domain points dismissed there are 41 allocated pointers with 44536 bytes memory, maximal heapsize up to now 56912 bytes (adr.of first p. 32260) setup for NS3-distr.: A=0.030000 B=0.100000 C=0.200000 D 0.500000 E 0.010000 as volume unbounded setup is restarted using the auxiliary domain points 1, volume below hat 6.622113e+01,(*t).gennr 4 points 2, volume below hat 5.686126e+01,(*t).gennr 8 points 3, volume below hat 3.588562e+01,(*t).gennr 14 points 4, volume below hat 2.964845e+01,(*t).gennr 20 points 5, volume below hat 1.952663e+01,(*t).gennr 28 points 6, volume below hat 1.570774e+01,(*t).gennr 34 points 6, volume below hat 1.646918e+01,(*t).gennr 20 Starting values found, auxiliary domain points dismissed there are 53 allocated pointers with 58568 bytes memory, maximal heapsize up to now 73336 bytes (adr.of first p. 32260) setup for beta-distr.: a1=20.000000 a2=3.000000 a3=3.000000 points 1, volume below hat 1.349300e+37,(*t).gennr 2 setup for NS1 (gamma)-distr.: as volume unbounded setup is restarted using the auxiliary domain points 1, volume below hat 2.321723e+02,(*t).gennr 4 points 2, volume below hat 2.992355e+01,(*t).gennr 8 points 3, volume below hat 2.935729e+01,(*t).gennr 14 points 4, volume below hat 8.119297e+00,(*t).gennr 20 points 5, volume below hat 2.847042e+00,(*t).gennr 26 points 6, volume below hat 2.371332e+00,(*t).gennr 32 points 6, volume below hat 2.372066e+00,(*t).gennr 20 Starting values found, auxiliary domain points dismissed there are 71 allocated pointers with 83280 bytes memory, maximal heapsize up to now 97816 bytes (adr.of first p. 32260) setup for gamma-distr.: a1=4.000000 a2=3.000000 l1=3.000000 l2=2.000000 b1=3.000000 b2=1.000000 as volume unbounded setup is restarted using the auxiliary domain points 1, volume below hat 1.697383e+13,(*t).gennr 4 points 2, volume below hat 4.430912e+05,(*t).gennr 8 points 3, volume below hat 5.481583e+03,(*t).gennr 14 points 4, volume below hat 5.023574e+03,(*t).gennr 20 points 5, volume below hat 5.505724e+02,(*t).gennr 26 points 5, volume below hat 5.825382e+02,(*t).gennr 18 Starting values found, auxiliary domain points dismissed there are 82 allocated pointers with 96776 bytes memory, maximal heapsize up to now 110624 bytes (adr.of first p. 32260) We generate 5 matrices , which contain a random pair of each distribution The parameters are those that we have specified above points 8, volume below hat 1.158865e+03,(*t).gennr 28 points 3, volume below hat 7.484789e+02,(*t).gennr 8 points 4, volume below hat 2.136415e+02,(*t).gennr 14 points 2, volume below hat 3.482571e+29,(*t).gennr 6 points 3, volume below hat 1.103226e+27,(*t).gennr 12 points 4, volume below hat 1.504031e+05,(*t).gennr 18 points 5, volume below hat 2.031520e-01,(*t).gennr 24 points 6, volume below hat 8.873361e-07,(*t).gennr 30 points 7, volume below hat 7.990677e-08,(*t).gennr 38 points 7, volume below hat 1.986958e+00,(*t).gennr 28 points 6, volume below hat 4.697298e+02,(*t).gennr 26 matrix 0 normal: (-8.39|2.30) restricted-dnormal:(12.42|10.23) NS2: (1.73|2.07) cut normal: (0.47|-0.21) NS3 (1.53|-0.43) beta: (0.86|0.05) NS1: (1.95|-2.64) Gamma: (6.17|2.11) points 5, volume below hat 1.387881e+02,(*t).gennr 18 points 6, volume below hat 4.807526e+01,(*t).gennr 11 points 7, volume below hat 3.960026e+01,(*t).gennr 14 points 4, volume below hat 7.724189e+00,(*t).gennr 6 points 7, volume below hat 1.475588e+01,(*t).gennr 28 points 8, volume below hat 5.415649e-08,(*t).gennr 46 points 8, volume below hat 1.826199e+00,(*t).gennr 36 matrix 1 normal: (-30.19|-4.96) restricted-dnormal:(2.85|-0.52) NS2: (-0.50|1.56) cut normal: (0.62|2.38) NS3 (0.67|-0.99) beta: (0.74|0.06) NS1: (0.85|-0.87) Gamma: (3.94|2.00) points 9, volume below hat 1.128710e+03,(*t).gennr 33 points 10, volume below hat 1.007647e+03,(*t).gennr 41 points 6, volume below hat 1.300213e+02,(*t).gennr 23 matrix 2 normal: (-14.86|-3.26) restricted-dnormal:(5.25|-3.54) NS2: (-0.08|-2.43) cut normal: (-0.93|0.46) NS3 (1.21|0.41) beta: (0.72|0.11) NS1: (0.51|-0.46) Gamma: (3.35|2.47) points 7, volume below hat 1.214385e+02,(*t).gennr 28 points 8, volume below hat 1.176747e+02,(*t).gennr 36 points 5, volume below hat 6.465530e+00,(*t).gennr 10 points 6, volume below hat 5.498724e+00,(*t).gennr 16 points 9, volume below hat 5.240482e-08,(*t).gennr 52 points 9, volume below hat 1.787739e+00,(*t).gennr 44 matrix 3 normal: (-8.00|-10.60) restricted-dnormal:(7.44|6.99) NS2: (1.24|2.09) cut normal: (-0.74|0.60) NS3 (-1.47|-0.46) beta: (0.65|0.26) NS1: (1.86|-1.77) Gamma: (5.06|3.72) points 8, volume below hat 1.400285e+01,(*t).gennr 36 points 10, volume below hat 1.689096e+00,(*t).gennr 52 matrix 4 normal: (17.64|9.67) restricted-dnormal:(2.65|-0.10) NS2: (0.59|0.92) cut normal: (-0.17|1.40) NS3 (-1.71|-0.34) beta: (0.79|0.09) NS1: (1.27|-1.14) Gamma: (6.69|0.56) If you do not like the messages of the sample2d-function, that gives information about the updated hat, you have to change the define OUTPUT to -1 in alc2d.h We are generating 50000 pairs from the normal distribution points 15, volume below hat 4.940656e+02,(*t).gennr 88 points 20, volume below hat 4.067971e+02,(*t).gennr 120 points 25, volume below hat 3.798256e+02,(*t).gennr 160 points 30, volume below hat 3.698284e+02,(*t).gennr 200 points 35, volume below hat 3.592098e+02,(*t).gennr 240 points 40, volume below hat 3.561572e+02,(*t).gennr 280 points 45, volume below hat 3.514514e+02,(*t).gennr 320 points 50, volume below hat 3.485788e+02,(*t).gennr 360 Polygon structure was freed Observed average acceptance probability d3.1: 0.933132 there are 95 allocated pointers with 101256 bytes memory, maximal heapsize up to now 154728 bytes (adr.of first p. 35398) We are generating 50000 pairs from the restricted normal dist. As we allowed only only 20 design points (calling the set-up) we expect a lower acceptance prob. points 9, volume below hat 1.115039e+02,(*t).gennr 42 points 10, volume below hat 1.108589e+02,(*t).gennr 46 points 12, volume below hat 1.079130e+02,(*t).gennr 64 points 14, volume below hat 1.056893e+02,(*t).gennr 80 points 16, volume below hat 1.033868e+02,(*t).gennr 92 points 18, volume below hat 1.014615e+02,(*t).gennr 108 points 20, volume below hat 9.993151e+01,(*t).gennr 124 Polygon structure was freed Observed average acceptance probability d3.1: 0.930111 The maximal number of corners of one polygon up to now was 9 . there are 81 allocated pointers with 91576 bytes memory, maximal heapsize up to now 154728 bytes (adr.of first p. 38048) We are generating 50000 pairs from the NS2 distribution This function has a big constant part and thus the acceptance-probability is high. points 8, volume below hat 3.899052e+01,(*t).gennr 17 points 9, volume below hat 3.804110e+01,(*t).gennr 20 points 10, volume below hat 3.720609e+01,(*t).gennr 23 points 15, volume below hat 3.613462e+01,(*t).gennr 38 points 20, volume below hat 3.591101e+01,(*t).gennr 53 points 25, volume below hat 3.562934e+01,(*t).gennr 68 points 30, volume below hat 3.550653e+01,(*t).gennr 83 points 35, volume below hat 3.541690e+01,(*t).gennr 98 points 40, volume below hat 3.537068e+01,(*t).gennr 113 points 45, volume below hat 3.535063e+01,(*t).gennr 128 points 50, volume below hat 3.533414e+01,(*t).gennr 143 Polygon structure was freed Observed average acceptance probability d3.1: 0.997685 The maximal number of corners of one polygon up to now was 49 . This high number comes from the big, circle-shaped plateau of the NS2-distr. there are 68 allocated pointers with 77896 bytes memory, maximal heapsize up to now 154728 bytes (adr.of first p. 3adb0) We are generating 50000 pairs from the cut-normal distribution points 7, volume below hat 5.224973e+00,(*t).gennr 24 points 8, volume below hat 5.114996e+00,(*t).gennr 32 points 9, volume below hat 4.837106e+00,(*t).gennr 40 points 10, volume below hat 4.751297e+00,(*t).gennr 48 points 15, volume below hat 4.621735e+00,(*t).gennr 81 points 20, volume below hat 4.442856e+00,(*t).gennr 124 points 25, volume below hat 4.387586e+00,(*t).gennr 164 points 30, volume below hat 4.324068e+00,(*t).gennr 204 points 35, volume below hat 4.294795e+00,(*t).gennr 244 points 40, volume below hat 4.252324e+00,(*t).gennr 284 points 45, volume below hat 4.233574e+00,(*t).gennr 324 points 50, volume below hat 4.223386e+00,(*t).gennr 364 Polygon structure was freed Observed average acceptance probability d3.1: 0.974298 there are 56 allocated pointers with 64048 bytes memory, maximal heapsize up to now 154728 bytes (adr.of first p. 3adf8) We are generating 50000 pairs from the NS3 distribution points 9, volume below hat 1.384203e+01,(*t).gennr 44 points 10, volume below hat 1.356717e+01,(*t).gennr 52 points 15, volume below hat 1.138380e+01,(*t).gennr 84 points 20, volume below hat 1.069460e+01,(*t).gennr 124 points 25, volume below hat 1.040236e+01,(*t).gennr 164 points 30, volume below hat 1.026900e+01,(*t).gennr 204 points 35, volume below hat 1.005843e+01,(*t).gennr 244 points 40, volume below hat 9.901928e+00,(*t).gennr 280 points 45, volume below hat 9.792065e+00,(*t).gennr 316 points 50, volume below hat 9.725706e+00,(*t).gennr 352 Polygon structure was freed Observed average acceptance probability d3.1: 0.942472 there are 42 allocated pointers with 48544 bytes memory, maximal heapsize up to now 154728 bytes (adr.of first p. 3ab10) We are generating 50000 pairs from the beta distribution points 10, volume below hat 5.035690e-08,(*t).gennr 58 points 15, volume below hat 3.813877e-08,(*t).gennr 98 points 20, volume below hat 3.642572e-08,(*t).gennr 138 points 25, volume below hat 3.546012e-08,(*t).gennr 176 points 30, volume below hat 3.443367e-08,(*t).gennr 214 points 35, volume below hat 3.384111e-08,(*t).gennr 254 points 40, volume below hat 3.352087e-08,(*t).gennr 294 points 45, volume below hat 3.319988e-08,(*t).gennr 334 points 50, volume below hat 3.298586e-08,(*t).gennr 374 Polygon structure was freed Observed average acceptance probability d3.1: 0.950914 there are 27 allocated pointers with 31744 bytes memory, maximal heapsize up to now 154728 bytes (adr.of first p. 3a128) We are generating 50000 pairs from the NS1 distribution points 15, volume below hat 1.498010e+00,(*t).gennr 86 points 20, volume below hat 1.410608e+00,(*t).gennr 126 points 25, volume below hat 1.349955e+00,(*t).gennr 162 points 30, volume below hat 1.313204e+00,(*t).gennr 208 points 35, volume below hat 1.296836e+00,(*t).gennr 244 points 40, volume below hat 1.279377e+00,(*t).gennr 284 points 45, volume below hat 1.271058e+00,(*t).gennr 322 points 50, volume below hat 1.262859e+00,(*t).gennr 366 Polygon structure was freed Observed average acceptance probability d3.1: 0.936417 there are 12 allocated pointers with 14424 bytes memory, maximal heapsize up to now 154728 bytes (adr.of first p. 3a170) We are generating 50000 pairs from the Gamma distribution points 7, volume below hat 4.110988e+02,(*t).gennr 34 points 8, volume below hat 4.061815e+02,(*t).gennr 42 points 9, volume below hat 4.030950e+02,(*t).gennr 48 points 10, volume below hat 3.926918e+02,(*t).gennr 52 points 15, volume below hat 3.649642e+02,(*t).gennr 92 points 20, volume below hat 3.563935e+02,(*t).gennr 132 points 25, volume below hat 3.483791e+02,(*t).gennr 177 points 30, volume below hat 3.441924e+02,(*t).gennr 215 points 35, volume below hat 3.420829e+02,(*t).gennr 252 points 40, volume below hat 3.368152e+02,(*t).gennr 292 points 45, volume below hat 3.352364e+02,(*t).gennr 332 points 50, volume below hat 3.327285e+02,(*t).gennr 370 Polygon structure was freed Observed average acceptance probability d3.1: 0.962853 there are 0 allocated pointers with 0 bytes memory, maximal heapsize up to now 154728 bytes (adr.of first p. 0) SHAR_EOF fi # end of overwriting check if test -f 'res3' then echo shar: will not over-write existing file "'res3'" else cat << "SHAR_EOF" > 'res3' setup and sample of size 100000 for normal-distr.: R=0.800000 S1=1.000000 S2=1.000000 as volume unbounded setup is restarted using the auxiliary domain points 1, volume below hat 4.520841e+01,(*t).gennr 4 points 2, volume below hat 2.345216e+01,(*t).gennr 8 points 3, volume below hat 1.768986e+01,(*t).gennr 14 points 4, volume below hat 1.314301e+01,(*t).gennr 20 points 5, volume below hat 8.553789e+00,(*t).gennr 28 points 5, volume below hat 1.119169e+01,(*t).gennr 12 Starting values found, auxiliary domain points dismissed there are 11 allocated pointers with 7088 bytes memory, maximal heapsize up to now 15104 bytes (adr.of first p. 30f50) points 6, volume below hat 9.901565e+00,(*t).gennr 17 points 7, volume below hat 7.526778e+00,(*t).gennr 21 points 8, volume below hat 7.250667e+00,(*t).gennr 29 points 9, volume below hat 5.590759e+00,(*t).gennr 40 points 10, volume below hat 5.028724e+00,(*t).gennr 48 points 12, volume below hat 4.928383e+00,(*t).gennr 60 points 14, volume below hat 4.723411e+00,(*t).gennr 72 points 16, volume below hat 4.576252e+00,(*t).gennr 84 points 18, volume below hat 4.483766e+00,(*t).gennr 100 points 20, volume below hat 4.411470e+00,(*t).gennr 120 Polygon structure was freed Observed average acceptance probability: 0.854942 volume below the hat approximately 4.409552 there are 0 allocated pointers with 0 bytes memory, maximal heapsize up to now 20016 bytes (adr.of first p. 0) Marginal x: Chi2-test: samplesize=100000 number of classes= 30 Chi2-value 29.703800 Approximate P-Value: 0.428885 Marginal y: Chi2-test: samplesize=100000 number of classes= 30 Chi2-value 29.343800 Approximate P-Value: 0.447293 : Chi2-test: samplesize=100000 number of classes= 30 Chi2-value 31.644800 Approximate P-Value: 0.335504 chi2 2 dim: Chi2-test: samplesize=100000 number of classes= 900 Chi2-value 846.746000 Approximate P-Value: 0.892749 SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'alc2d.h' then echo shar: will not over-write existing file "'alc2d.h'" else cat << "SHAR_EOF" > 'alc2d.h' /* ALC2DLIB-Version 1.0 implemented by WH 22.3.99, all rights reserved */ /*please report problems or bugs to whoer@statistik.wu-wien.ac.at */ /* or hormannw@boun.edu.tr */ /*this file: alc2d.h */ /*as the follwing constants are changed from time to time we put them right at the front of the code */ #define OUTPUT 0 /*default: 0 */ /*-1... nothing, 0...little, 1...more, 2..even more output about setup and generation*/ /*the following two defined constants can be used to check the dynamic memory allocation*/ #define MEMCONTROL 1 /*default 0*/ /*0... no controll,1..controll allmost without output,2... more output MEMCONTROL >=1 is allso necessary that the two diagnostic procedures for counting the number of trials and the number of corners work*/ #define MEMOUTPUT 0/*0..none, 2.. everything, 1.. only longer lasting memory*/ /* The algorithm implemented here is a universal algorithm for two-dimensional log-concave distributions. To use the algorithm it is first necessary to call the function "setup" which allocates the necessary memory and computes all constants necessary for generation. Of course set-up needs as input several informations about the distribution. It returns a pointer to the structure of the constants. After calling set-up once, pairs of the two-dimensional distribution are generated by calling the function "sample2d", which needs as input only the pointer to the structure of constants returned by setup. During these calls to "sample2d" the setup procedure is still improving the fit of the hat-function until the maximal number of design-points is reached. If no further samples of the distribution are needed the allocated memory must be freed using "freesetup". The use of these functions is demonstrated in the samplefiles main.c and main1.c Of course a uniform generator is needed. Please include the call of your favorite generator or of the generator included in you C-compiler. We only have included a very small generator in our sample files. As it has a period of only 2^32=4.1e9 numbers it should not be used for big simulations!!!*/ #define UNIFORM unif/*name of your uniform generator; the generator is called without an argument: UNIFORM()*/ /* The input of the below function "setup" is very important, as it must include all information necessary for sampling from the distribution. Before you can use the algorithm to generate variates form the desired 2dimensional distribution the following steps are necessary: 1) Check that the density f(x) is log concave, i.e. log(f(x)) must be a concave function over the entire domain. 2) Code the function h(x)=log(f(x)) in the below style: (All parameters must be included in the double array param[];) double h(double x[2],double param[]) For example for the standard normal distribution: double hstnor(double x[2],double param[]) { return -0.5*(x[0]*x[0]+x[1]*x[1]); } For further examples see the sample files. 3) Code the two partial derivatives of h(x) in the same style as h. 4) Define the domain of the distribution (a possibly open convex polygon) by storing all border lines of this polygon in an array eq[][3]. lines must be stored in the form 0 = eq[0][0] + eq[0][1]*x + eq[0][2]*y; Remark: The domain of the distribution is only fully defined together with the starting point, described in step 5. 5) Provide one (ore more) starting-points, if possible close to the mode of the distribution. Of course all starting points must be inside the domain of the distribution!! 6) If the domain of the distribution is unbounded it can be necessary to add an "auxiliary domain" in order that the set-up algorithm can construct a bounded hat function. This auxiliary domain must be a bounded convex polygon, that is inside the domain of the distribution and must allways contain the mode of the distribution. To define this auxiliary domain add the border lines of the auxiliary domain in the array eq[][3] after the equalities defining the domain of the distribution. In most cases the auxiliary domain can simply be a rectangle containing the mode. (see the sample files for examples) 7) The last information required is the maximal number of design-points, that can be used by the algorithm. Sensible values are between 10 and 30, if the desired samplesize is small to moderate and should go up to about 100 if good generation speed is desired for large samples. */ extern void *setup( int nrsp,/*Number of starting points*/ int nn, /* maximal number of design points*/ double sp[][2],/* array of starting points*/ int nreq, /*number of equalities*/ double eq[][3],/* array of the equalities defining the domain equalities are stored in the form 0 = eq[0][0] + eq[0][1]*x + eq[0][2]*y ;*/ int nrauxun,/*total number of equalities for domain and auxiliary domain we must have nrauxun >=nreq nrauxun== nreq means that no auxiliary domain is provided*/ double (*h)(double x[2],double param[]),/*logarithm of density function*/ double (*hx)(double x[2],double param[]), /*partial derivative (with respect to x) of h*/ double (*hy)(double x[2],double param[]), /*partial derivative (with respect to x) of h*/ double param[],/*current setting of the parameters, NULL if nparam=0*/ int nparam/*number of parameters of the density function*/); extern void sample2d(double res[2],void *t1); /*samples from the 2-dimensional distribution, generated random pair is stored in res[2] t1 must be the pointer that was returned by the "setup" function*/ extern void freesetup(void *t); /*frees all memory that was assigned during the set-up*/ /*The follwing are some internal constants for the algorithm; unless you try to optimize the speed or the memory requirements of the sampling algorithm or of the set-up it is not necessary to change them*/ #define CFGT 20 /* default 20; recommended intervall: 10 to 30 */ /*small values leed to a faster set-up but to a slower sampling-time*/ /*length of guide-table q is determined as CFGT*nn */ #define CFGTSU (0.2) /*default 0.2; must be <=1. */ /*can be used to optimize the speed of the set-up*/ /*length of guide table during setup = CFGT*nn*CFGTSU */ #define NC 10 /*default 10; must be >=4 */ /*A small value of NC can reduce the memory size, but slows down the set-up*/ /*starting number of corners of polygon*/ #define SUNC 10 /*default 10; should be >=5*/ /*for for small values of SUNC the set-up is much faster. High values of SUNC eg. 100 imply a slower set-up but a slightly better fit of the hat. In sampled2d() allways nn/SUNC rejected points together are used as new construction points and a new hat is computed*/ /*the follwing function can be used to find the maximal number of corners (ncmax) of a single polygon that occured in this program.*/ extern int returnncmax(void); /*returns the current value of ncmax The function olny works correctly if MEMCONTOL>=1, otherwise -1 is returned */ /* The following functions can be used to count the number of trials to compute the acceptance probability in the rejection algorithm they only work if MEMCONTROL>=1 is true. In this case for every trial of the rejection-algorithm in sample2d the counter widcount is automatically advanced. PROBLEM: sample2d is also called from the setup function, if the auxiliary domain is utilized to find good starting points. In this case widcount is advanced allthough the user has never called sample2d(). Therefore it is best to reset widcount between the call to setup() and the first call to sample2d(). It must also be considered, that the acceptance probability changes during the first calls to sample2d() as the hat is updated.*/ extern void setwidcount(unsigned long int newval); /*sets the value of widcount to the value passed to the function; is normally used to reset the widcount value to 0 */ extern unsigned long int returnwidcount(void); /*returns the current value of widcount*/ extern void pr_2dlib_def(); /*prints the current settings of the defined constants*/ #ifndef NULL #define NULL (void *) 0 #endif SHAR_EOF fi # end of overwriting check if test -f 'src.c' then echo shar: will not over-write existing file "'src.c'" else cat << "SHAR_EOF" > 'src.c' /* ALC2DLIB-Version 1.0 implemented by WH 22.3.99, all rights reserved */ /*please report your experience to whoer@statistik.wu-wien.ac.at */ /* or hormannw@boun.edu.tr */ /* This file is the main program file alc2d.c for explanations see: alc2d.h */ #include #include #include #include "alc2d.h" extern double UNIFORM(void); /* The following constants should only be changed, if the double type is shorter than 64 bit*/ #define EPS 1.e-13 /*tolerance for linecut*/ #define EPS1 1.e-13 /*tolerance for checksameside*/ #define EPS2 1.e-13 /*tolerance for checkvector point*/ #define EPSTW 1.e-6 /*tolerance for a=0 (numerically tested)*/ #define VOLHATMAX 1.e99 /*large positive value at most HUGE_VAL*/ void pr_2dlib_def(void) /*prints the values of the defined constants*/ { printf("OUTPUT %d SUNC %d MEMCONTROL %d MEMOUTPUT %d VOLHATMAX %e\n", OUTPUT,SUNC,MEMCONTROL,MEMOUTPUT,VOLHATMAX); printf("NC %d CFGT %d EPS %e EPS1 %e EPS2 %e\n",NC,CFGT,EPS,EPS1,EPS2); } /*some output-macros for debugging*/ #if (MEMOUTPUT==2) #define MPR2(a) printf a #else #define MPR2(a) #endif #if (MEMOUTPUT>=1) #define MPR1(a) printf a #else #define MPR1(a) #endif #if (OUTPUT==2) #define OPR2(a) printf a #else #define OPR2(a) #endif #if (OUTPUT>=1) #define OPR1(a) printf a #else #define OPR1(a) #endif #if (OUTPUT>=0) #define OPR0(a) printf a #else #define OPR0(a) #endif /*some macro-definitions*/ #define OK 1239999 #ifndef NULL #define NULL (void *) 0 #endif #define Arrlet2(a,b) do{a[0]=b[0];a[1]=b[1];} while(0) #define Arrlet2zero(a) do{a[0]=0.;a[1]=0.;} while(0) #define Arrlet3(a,b) do{a[0]=b[0];a[1]=b[1];a[2]=b[2];} while(0) #define Arradd2(a,b,c) do{a[0]=b[0]+c[0];a[1]=b[1]+c[1];} while(0) #define Arrsub2(a,b,c) do{a[0]=b[0]-c[0];a[1]=b[1]-c[1];} while(0) #define Arrminuslet2(a,c) do{a[0]= -c[0];a[1]= -c[1];} while(0) #define Arrchangesign2(a) do{a[0]= -a[0];a[1]= -a[1];} while(0) #define Letstru(a,b) memcpy(&(a),&(b),sizeof(b)) #define Sqr(c) ((c)*(c)) #define Sqrdistance(a,b) (Sqr(b[0]-a[0])+Sqr(b[1]-a[1])) #define MIN(a,b) ((a)<(b)?a:b) #define MAX(a,b) ((a)>(b)?a:b) #define Hatvalue(t,c) (t[1]==0.||t[2]==0.?t[0]+t[1]*c[0]+t[2]*c[1]:t[1]*t[2]*(t[0]/(t[1]*t[2])+c[0]/t[2]+c[1]/t[1])) #define Directionvector(a,b) do{a[0]= -b[2];a[1]= b[1];} while(0) /*a[2]..result, b[3]...line, stores the direction vector of b[3] in a[2]*/ typedef struct { int nc/*number of corners of polygon+unbounded edges*/, nca/*number of corners for which memory was allocated up to now*/, update/* 1..update necessary, 0.. no update necessary*/, open/*0..closed pol,1..open,2..twoparallel lines,3..one line,4..empty*/, index/*index number of that polygon in structure TOTAL*/ ; double design[2]/*tangential point*/, tang[3]/*tangential space*/, (*corner)[2]/*pointer to the array of the vertices in consecutive order*/, sc[2]/*corner with the highest z-value*/, s /*value of the tangential plane over sc, i.e. highest value of the hat over the polygon*/, rot[2]/*rotation that rotates the negative gradient of tang in x-direction*/, a/*slope of the hat in the direction of the gradient*/; } POLYGON; /*the convex polygon associated with one point of contact */ typedef struct { int open/*kind of region,0(=normal a negative),1(=normal a positive),2 open*/, ip/*index of polygon containing this generator region*/; double point[2],/*coordinates of the point translated into origin*/ rot[2],/*rotation to rotate negative gradient into (1,0)*/ k[2],/*slopes of the two lines from origin*/ s,/** z-coordinate of point translated into origin*/ a,/* h(x1,x2)=exp(a*x+const) is the hat over rotated generator region*/ b,/*for open<=1 the generator region is in the interval (0,b) for open=2 b is broadness of the parallel strip*/ vol,/**volume of this part*/ vols;/**cumulated volume*/ } GEN; /*All informations about one generator region*/ typedef struct { int n,/* current number of touching points*/ nn,/* maximum number of touching points*/ gennr,/*current number of generator regions*/ ngennr,/*current maximum number of generator regions*/ m,/*current length of guidetable*/ mm,/*maximum length of guide table*/ *q/* pointer to the integer array for the guide table*/; long ok;/*used to check if the correct pointer was passed to sample ok==OK means that the setup was successfull*/ GEN *gen/*[nngen] pointer to the array of generator regions*/; POLYGON *p;/*[nn+1] pointer to the array of POLYGON; only necessary as long as new points of contact are added*/ double (*h)(double x[2],double param[]), (*hx)(double x[2],double param[]),(*hy)(double x[2],double param[]), /*pointers to the functions evaluating the log of the density and the partial derivatives of the density*/ *param,/*pointer to the double array with the current parameter settings*/ oldhatvol;/*total volume below old hat used to check against rounding errors*/ } TOTAL; /*contains all information necessary for set-up and sampling*/ /*******************************************************************/ #if(MEMCONTROL) /*definition of some auxiliary functions used for memory controll*/ typedef struct{void *adress;unsigned long int size;}HEAP; static int nrofap=0/*number of allocated pointers*/, maxsize=0/*maximal size of allocated memory observed so far*/, actsize=0/*current size of allocated memory*/; static HEAP heap[1000]; void checkheap() /*computes the current values of the variables maxsize and actsize*/ { int i; actsize=0; if(nrofap) for(i=0;i=0) Letstru(heap[i],heap[i+1]); } if(found>=0) { nrofap--; free(poi); checkheap(); #if (MEMCONTROL==2) printheap(); #endif } else if(poi!=NULL) printf("ERROR!! memtest: a pointer freed which is not in the list\n"); } void *reallocwh(void *poiold,long int size) /*used instead of realloc() to perform the allocation checks*/ { void *poinew; int i,found= -1; if(size<=0) printf("ERROR!!reallocwh was called with size <=0\n"); for(i=0;i=0) Letstru(heap[i],heap[i+1]); } if(found>=0) nrofap--; else if(poiold!=NULL) printf("ERROR!! memtest: a pointer reallocated which is not in the list\n"); poinew=realloc(poiold,size); if(poinew!=NULL) { heap[nrofap].adress=poinew; heap[nrofap].size=size; nrofap++; } else { printf("Error!!reallocwh memory was not allocated successfully!!\n"); if(size>0) exit(1); } checkheap(); #if (MEMCONTROL==2) printheap(); #endif return(poinew); } #else void *mallocwh(int size) { return(malloc(size));} void freewh(void *poi) { free(poi);} void *reallocwh(void *poiold,int size) { return(realloc(poiold,size));} #endif /**************************************************************************/ /*other utilities for debugging*/ void printgen(GEN *g) /*prints all information of one generator region*/ { printf("open %d ip %d point(%4.3e|%4.3e) rot(%4.3e|%4.3e) k(%4.3e|%4.3e)\n", g->open,g->ip,g->point[0],g->point[1],g->rot[0],g->rot[1],g->k[0],g->k[1]); printf("s %4.3e a %4.3e b %4.3e vol %4.3e vols %4.3e\n", g->s,g->a,g->b,g->vol,g->vols); } void printpoly(POLYGON *poly) /*prints all information of one polygone*/ { int i; printf("n:%d update:%d open:%d design:(%4.3e|%4.3e) tang:(%4.3e|%4.3e|%4.3e)\n", (*poly).nc,(*poly).update,(*poly).open,(*poly).design[0], (*poly).design[1],(*poly).tang[0],(*poly).tang[1],(*poly).tang[2]); for(i=0;i<=((*poly).nc);i++) printf("%d: (%4.3e|%4.3e)\n",i,(*poly).corner[i][0],(*poly).corner[i][1]); } /****************************************************************************/ #if (MEMCONTROL>=1) static int ncmax=0;/*maximal number of corners observed for a polygon*/ int returnncmax(void) /*returns the current value of ncmax*/ { return(ncmax);} #else int returnncmax(void) /*returns -1 as MEMCONTOL==0*/ { printf("Caution! function returnncmax() does not work correctly as MEMCONTROL not >= 1\n"); return(-1);} #endif /***********************************************************************/ double gg2neu(double ab) { double u,h1,h2; long int fak=1,i; u=UNIFORM()*(exp(ab)*(ab-1.)+1.); h1=ab*ab; h2=h1*0.5; /*printf("h2 %f\n",h2);*/ if(u0) {/*Fall a>0*/ ab=a*b; if(ab<1.) return gg2neu(ab)/a; tlx=ab*0.65; al=1.+1./tlx; col=exp(al*(0.-ab))-1.; while(1) { x=ab+log(UNIFORM()*col+1.)/al; h=x-tlx-al*(x-tlx); v=UNIFORM()*tlx; if (v<=x*(1.+h)) return(x/a); if (v<=x*(1.+h*(1.+h*(0.5+h*(1./6.))))) return(x/a); if (v<=x*exp(h)) return(x/a); } } else if(a<0)/*Fall a<0*/ { ab= -a*b; if(abk[1]-ge->k[0]); a=fabs(ge->a); b=fabs(ge->b); u=UNIFORM()*(ak10/a+b); if(u0.); dhelp2=g[0]+p2[0]*g[1]+p2[1]*g[2]; help2=(dhelp2>0.); if (fabs(dhelp1)nc>=p->nca-1)/*da erst unten p->nc erhoeht wird*/ { MPR1(("corner at adress %p with %d points is\n",p->corner,p->nca)); p->nca+=NC; p->corner=reallocwh(p->corner,p->nca*2*sizeof(double)); if(p->corner==NULL) printf("error realloc no space in heap\n"); MPR1(("(realloc)is enlargened to %d points (%d bytes) at adress %p\n",p->nca,p->nca*2*sizeof(double),p->corner)); } for(i=p->nca-1;i>nr;i--) Arrlet2(p->corner[i],p->corner[i-1]); Arrlet2(p->corner[nr],point); p->nc++; #if (MEMCONTROL>=1) if(ncmaxnc) ncmax=p->nc; #endif } /*****************************************************************/ #define POL (*poly) int checkcorner(POLYGON (*poly),double line[3],double helpp[2][2],int a) /*calls the two versions of check same side depending on the value of i, the number of the corner helpp is used to store the intersection points between vector and line, which are computed in the function checksamesinfty()*/ { int n; n=POL.nc; if (POL.open&&a==0) return( checksamesinfty(POL.corner[0],POL.corner[1],POL.design,line,helpp[0])); else if(POL.open&&a==POL.nc-1) return( checksamesinfty(POL.corner[n-1],POL.corner[n-2],POL.design,line,helpp[1])); else return(checksameside(line,POL.corner[a],POL.design)); } /*****************************************************************/ void intersectpolline(POLYGON (*poly),double line[3],double helpp[2][2], int i1,int i2,double intp[2]) /*calls the different versions to compute the intersection with the polygon between the neighbouring corners with numbers i1 and i2 result is stored in intp the result of the intersection of vector with line was allready computed in the function checksamesinfty() and stored in helpp[2][2]*/ { int k; if (i1>i2) {k=i1;i1=i2;i2=k;} if (POL.open) { if(i1==0&&i2==POL.nc-1) { intp[0]= -line[2]; intp[1]= line[1]; return; } else if(i1==0) { Arrlet2(intp,helpp[0]);return;} else if(i1==POL.nc-2) { Arrlet2(intp,helpp[1]);return;} } intersect(POL.corner[i1],POL.corner[i2],line,intp); } /*****************************************************************/ void deletecorners(POLYGON (*poly),int corem[]) /*deletes all corners of the polygon with number i and corem[i]<=0*/ { int i,j=0; while(corem[j]>=1&&j=1) { Arrlet2(POL.corner[j],POL.corner[i]); j++; } POL.nc=j; } /*****************************************************************/ void checkvectorsigns(POLYGON (*poly)) /* checks the direction of the vectors of open polygons and changes them if the touching point is not inside the polygon*/ { double hl[3],hp[2]; int i; /*change vector 0*/ if (POL.nc==3) changelineeq(POL.corner[2],POL.corner[1],1,hl); /* else changelineeq(POL.corner[2],POL.corner[1],0,hl);*/ /*stattdessen fuer Spezialproblem falls viele Punkte sehr eng beisammen*/ else { for(i=2;(i0)&& (Sqrdistance(POL.corner[POL.nc-2],POL.corner[i])<1.e-10);i--); if(i==0) changelineeq(POL.corner[0],POL.corner[POL.nc-2],1,hl); else changelineeq(POL.corner[POL.nc-2],POL.corner[i],0,hl); } Arradd2(hp,POL.corner[POL.nc-2],POL.corner[POL.nc-1]); if(!checksameside(hl,hp,POL.design)) Arrchangesign2(POL.corner[POL.nc-1]); } /*****************************************************************/ void nocornerupdate(POLYGON (*poly),double line[3]) /*a polygon (which has up to now no corners) is updated by the line[]*/ { double hline1[3],hline2[3],help1,help2,help3; POL.update=1;/*polygon is allways considered as updated*/ switch(POL.open) { case 4:/*polygon up to now empty*/ { POL.open=3; Arrlet3(POL.corner[0],line);/*corner is used to store the line*/ break; } case 3:/*up to now only one line stored in corner*/ { Arrlet3(hline1,POL.corner[0]); if (linecut(hline1,line,POL.corner[1])) { POL.nc=3;/*case that the two lines have a cutting point*/ POL.open=1; Directionvector(POL.corner[0],line); Directionvector(POL.corner[2],hline1); checkvectorsigns(poly); } else/*hline1 and line are allmost parallel thus it is computed if both lines or only one line are stored*/ { if (fabs(hline1[1])>EPS) help1=hline1[0]*line[1]/hline1[1]; else help1=hline1[0]*line[2]/hline1[2]; help2= -line[1]*POL.design[0]-line[2]*POL.design[1]; if ((help1help2&&help2>line[0])) { POL.open=2;/*case: both lines are stored*/ /*corner[2] is used to store the 2nd parallel line*/ Arrlet3(POL.corner[2],line); } else if((help2line[0]&&line[0]>help1)) Arrlet3(POL.corner[0],line);/*case: line is stored, hline1 discarde*/ /*else case: as hline1 is kept nothing is changed*/ } break; } case 2:/*up to now two parallel lines stored in corner*/ { Arrlet3(hline1,POL.corner[0]); Arrlet3(hline2,POL.corner[2]); if (linecut(hline1,line,POL.corner[1])) /* case: the new line is not parallel to the two old ones*/ { POL.nc=4; POL.open=1; linecut(hline2,line,POL.corner[2]); /*the direction vector (-line[2],line[1]) is the vector twoards infinity*/ Directionvector(POL.corner[0],hline1); Arrlet2(POL.corner[3],POL.corner[0]); checkvectorsigns(poly); } else/*case: all 3 lines hline1 hline2 and line are allmost parallel thus it is computed if one line has to be replaced*/ { if (fabs(hline1[1])line[0]&&line[0]>help1)||(help3line[0]&&line[0]>help2)||(help31) nocornerupdate(poly,line);/*special case of no corner*/ else { /*cases: POL.open=1 (=open) or 0 (=closed)*/ corem=mallocwh(sizeof(int)*POL.nca); MPR2(("for corem %d byte allocated at adress %p\n",sizeof(int)*POL.nca,corem)); for(i=0;i0) { for(i=0;iip2nr) POL.open=0; if (POL.open) checkvectorsigns(poly); } MPR1(("corem adress %p is freed\n",corem)); freewh(corem); corem=NULL; } } /*****************************************************************/ void rotation(double tang[3],double rot[2],double *a) /* computes a rotation such that the negative gradient of tang[3] points in direction of the positiv x-axis; the result is stored in rot[2] - the length of the gradient is stored in a*/ { double t1,t2,st; t1= tang[1]; t2= tang[2]; st= sqrt(t1*t1+t2*t2); (*a)=-st; if (st<=EPSTW) { (*a)=0.; rot[0]=1.; rot[1]=0.;} else { rot[0]= -t1/st; rot[1]= t2/st;} } /*??perhaps faster when implemented as a macro??*/ void rotate(double rot[2],double vektor[2],double res[2]) /* rotates vektor[2] with rotation rot[2]; result is stored in res[2]*/ { res[0]=rot[0]*vektor[0]-rot[1]*vektor[1]; res[1]=rot[1]*vektor[0]+rot[0]*vektor[1]; } void invrotate(double rot[2],double vektor[2],double res[2]) /*rotates vektor[2] with inverse rotation of rot[2];result is stored in res[2]*/ { res[0]=rot[0]*vektor[0]+rot[1]*vektor[1]; res[1]= -rot[1]*vektor[0]+rot[0]*vektor[1]; } /*****************************************************************/ void volumenc(GEN *ge) /*determines volume of closed generator region (i.e. open<=1)*/ #define GE (*ge) { double eda; if (fabs(GE.a)>EPSTW) { eda=1./GE.a; /* GE.vol= fabs(GE.k[1]-GE.k[0])*exp(GE.s)*eda*( (GE.b-eda)*exp(GE.a*GE.b)+eda); */ GE.vol= (fabs(GE.k[1]-GE.k[0])*(exp(GE.s)*eda*eda+ eda*(GE.b-eda)*exp(GE.a*GE.b+GE.s))); } else GE.vol=exp(GE.s)*fabs(GE.k[1]-GE.k[0])*GE.b*GE.b*0.5; /*case of constant hat*/ } void volumenu(GEN *ge) /*determines volume of unbounded generator region (i.e. open==2)*/ #define GE (*ge) { if (GE.a<-1.e-60) GE.vol=(exp(GE.s)/GE.a)*(fabs(GE.k[1]-GE.k[0])/GE.a-fabs(GE.b)); else { GE.vol= -99999999.; OPR1(("Probelm: Vol infinit, generator region of open==2 with a>=0\n")); } /*Remark vol of "eineck"=exp(s)*fabs(k1-k0)/(a*a) vol of parallel strip= -exp(s)*b/a where b is the broadness of the strip*/ } #undef GE /*****************************************************************/ void ngenplusplus(int *ngen,int ngennr,double vol) /*advances ngen, the current number of generator regions only if vol>0. otherwise ngen is not advanced and thus the generator region is overwritten*/ { if(vol>0.) (*ngen)++; /*due to the if generator regions with vol<=0. are overwritten*/ #if (MEMCONTROL>=1) if(*ngen>ngennr) printf("Error!!!! ngen %d ngennr %d too regions!!\n",*ngen,ngennr); #endif } /*****************************************************************/ void comp1stgen(POLYGON *poly,GEN gen[],int *ngen,int ngennr, double hpoi0[2],double hpoi1[2],double hpoi2[2],int output) /*compute everything of the first generator region of an open or closed triangular. hpoi[0] are the original coordinates of the top, hpoi[1] the translated roated coordinates of the second point, hpoi[2] the coordinates of the rotated vector pointing from the top output ... 0 absolute no output*/ { gen[*ngen].ip=POL.index; gen[*ngen].open=0;/*closed generator region, a<0*/ gen[*ngen].b=hpoi1[0]; if (hpoi1[0]>0.) gen[*ngen].k[0]=hpoi1[1]/hpoi1[0]; else gen[*ngen].k[0]= -99999999.; if (hpoi2[0]>0.) gen[*ngen].k[1]=hpoi2[1]/hpoi2[0]; else {gen[*ngen].k[0]= -99999999.; if(output) printf("!!!!Error problem!! hpoi2[1]/hpoi2[0] %e/%e\n",hpoi2[1],hpoi2[0]); } /*k[0] and k[1] can be interchanged as well*/ Arrlet2(gen[*ngen].point,hpoi0); Arrlet2(gen[*ngen].rot,POL.rot); gen[*ngen].s=Hatvalue(POL.tang,hpoi0); gen[*ngen].a=POL.a; /*due to round-off errors it is possible that b<0. therefore*/ if (gen[*ngen].b0*/ gen[*ngen].s=POL.s+POL.a*point1[0]; Arrminuslet2(gen[*ngen].rot,POL.rot); gen[*ngen].b=point1[0]-point0[0]; if(gen[*ngen].b>0.)gen[*ngen].k[0]= (point1[1]-point0[1])/gen[*ngen].b; else gen[*ngen].k[0]= 99999999.; /*k[0] and k[1] can be interchanged as well*/ if(point1[0]>0.) gen[*ngen].k[1]= point1[1]/point1[0]; else {gen[*ngen].k[1]= 99999999.; if(output) printf("!!!!Error point1[1]/point1[0] %e/%e\n",point1[1],point1[0]); } Arrlet2(gen[*ngen].point,opoint1); gen[*ngen].a= -POL.a; /*due to round-off errors it is possible that b<0. therefore*/ if (gen[*ngen].b<=0.) gen[*ngen].vol=0.; else volumenc(&(gen[*ngen])); ngenplusplus(ngen,ngennr,gen[*ngen].vol); } /*****************************************************************/ int compgeno(POLYGON (*poly),GEN gen[],int *ngen,int ngennr,int output) /*For an infinite biangle the generator regions are computed and stored in gen *ngen ...number of generator regions stored in gen up to now function returns 1 if everything is ok 0 if volume below the hat unbounded occurs output ...0 absolute no output*/ { double hpoi[5][2],helpp[2],hb,hs; /*hpoi[0] are the original coordinates of the top of the biangle, hpoi[1] the translated roated coordinates of the second point, hpoi[2] the coordinates of the rotated vector pointing from the top hpoi[3] the vector pointing from the second point, hpoi[4] the original coordinates of the second point in the case of a one angle the points hpoi[0] and hpoi[1] are identical!*/ if(Hatvalue(POL.tang,POL.corner[1])>=Hatvalue(POL.tang,POL.corner[POL.nc-2])) { Arrlet2(hpoi[0],POL.corner[1]); Arrsub2(helpp,POL.corner[POL.nc-2],hpoi[0]); rotate(POL.rot,helpp,hpoi[1]); rotate(POL.rot,POL.corner[0],hpoi[2]); rotate(POL.rot,POL.corner[POL.nc-1],hpoi[3]); Arrlet2(hpoi[4],POL.corner[POL.nc-2]); } else { Arrlet2(hpoi[0],POL.corner[POL.nc-2]); Arrsub2(helpp,POL.corner[1],hpoi[0]); rotate(POL.rot,helpp,hpoi[1]); rotate(POL.rot,POL.corner[POL.nc-1],hpoi[2]); rotate(POL.rot,POL.corner[0],hpoi[3]); Arrlet2(hpoi[4],POL.corner[1]); } /*the first (bounded) generator region of the biangel is calculated*/ comp1stgen(poly,gen,ngen,ngennr,hpoi[0],hpoi[1],hpoi[2],output); /*the below values of the second generator region are computed using values from the first region; therefore *ngen is increased afterwards this important as the first region is deleted if vol==0 */ if(POL.open==1 && POL.nc==3) hb=0.;/*case of oneangle*/ else if(hpoi[1][0]>EPS) hb= gen[*ngen].b*(gen[*ngen].k[1]-gen[*ngen].k[0]); else hb= -hpoi[1][1]; hs=gen[*ngen].s+POL.a*hpoi[1][0]; ngenplusplus(ngen,ngennr,gen[*ngen].vol); /*the second (unbounded) generator region of the biangel is calculated*/ gen[*ngen].ip=POL.index; gen[*ngen].open=2;/*open generator region*/ gen[*ngen].a= POL.a; gen[*ngen].s=hs; gen[*ngen].b=hb; Arrlet2(gen[*ngen].rot,POL.rot); /*hier darf k[1] und k[0] nicht vertauscht werden!! Warum??*/ if(hpoi[2][0]>1.e-15) gen[*ngen].k[1]=hpoi[2][1]/hpoi[2][0]; else { if(output) { OPR0(("Error!!compgeno2 Volume unbounded probably bad starting values\n")); OPR0(("hpoi[2][1] %e hpoi[2][0] %e \n",hpoi[2][1],hpoi[2][0]));} return(0);} if(hpoi[3][0]>1.e-15) gen[*ngen].k[0]=hpoi[3][1]/hpoi[3][0]; else { if(output) { OPR0(("Error!!compgeno3 Volume unbounded probably bad starting values\n")); OPR0(("hpoi[3][1] %e hpoi[3][0] %e \n",hpoi[3][1],hpoi[3][0]));} return(0); } /*k[0] and k[1] can be interchangend as well*/ Arrlet2(gen[*ngen].point,hpoi[4]); volumenu(&(gen[*ngen])); ngenplusplus(ngen,ngennr,gen[*ngen].vol); return(1); } /*****************************************************************/ #define CP1(i) (i+1)%POL.nc /*cyclic +1 */ #define CM1(i) (POL.nc+i-1)%POL.nc /*cyclic -1 */ int triangulate(POLYGON (*poly),GEN gen[],int *ngen,int ngennr,int output) /*the polygon (open=1 or 0) is decomposed into triangles and the generator regions are computed and stored in gen *ngen...number of generator regions stored in gen up to now function returns 1 if everything is ok 0 if volume below the hat unbounded occurs output...0 absolute no output*/ { int itop,i,poln; double z,zn,xmin,(*ctrp)[2]/*points to corners of transl. rotat. polygon*/, helpp[2],(*corn0p)[2]; if(POL.open>1) return(0);/*as volume is certainly infinity*/ /* this is a trick that allows to make the part of the bounded triangels with only one code for bounded and unbounded polygons. Only the number of corners POL.nc must be changed and the array containing the corners has to start with corner[1] instead of corner[0]*/ if(POL.open==1) { poln=POL.nc-2; corn0p=&(POL.corner[1]);} else { poln=POL.nc; corn0p=&(POL.corner[0]);} /*first the corner itop with the maximal z-value of the hat is searched */ /* if some z-values are equal, the corner with the minimal x-value is taken*/ z= -1.e100; xmin= 1.e99; rotation(POL.tang,POL.rot,&(POL.a)); if(fabs(POL.a)z) { itop=i; z=zn; xmin=corn0p[i][0];} } Arrlet2(POL.sc,corn0p[itop]); POL.s=z; ctrp=mallocwh(sizeof(double)*2*POL.nca); MPR1(("for ctrp %d byte allocated at adress %p\n",sizeof(double)*2*POL.nca, ctrp)); /*array ctrp Corners of Translated Rotated Polygon are computed*/ Arrlet2zero(ctrp[0]); for(i=1;ip void initpol(TOTAL (*t),int nr) /*initialises the polygon if(nrnn) the values of the startpolygon are used as defaults*/ { double hh,hhx,hhy; PO[nr].nca=NC; if(PO[nr].corner!=NULL) freewh(PO[nr].corner); PO[nr].corner=mallocwh(sizeof(double)*2*PO[nr].nca); MPR1(("initpol for PO[%d].corner %d points (%d bytes) at adress %p\n",nr,PO[nr].nca,sizeof(double)*2*PO[nr].nca,PO[nr].corner)); PO[nr].update=1; PO[nr].index=nr; if(nr==t->nn) { PO[nr].open=4;/*empty polygon*/ PO[nr].nc=0; } else { PO[nr].open=PO[t->nn].open; PO[nr].nc=PO[t->nn].nc; memcpy(PO[nr].corner[0],PO[t->nn].corner[0],sizeof(double)*2*MAX(PO[t->nn].nc,2)); /*this maximum for domain a half-plane*/ /* computation of the tangential plane in the point touch*/ hh=(t->h)(PO[nr].design,t->param); hhx=(t->hx)(PO[nr].design,t->param); hhy=(t->hy)(PO[nr].design,t->param); PO[nr].tang[0]=hh-PO[nr].design[0]*hhx-PO[nr].design[1]*hhy; PO[nr].tang[1]=hhx; PO[nr].tang[2]=hhy; } } /****************/ void removetp(TOTAL (*t),int nr) /*removes polygon with number nr out of polygon array and closes the gap; (*t).n is updated aswell! nothing else is changed!*/ { int i; ((*t).n)--; MPR1(("Polygon %d deleted therefore corners at adress %p are freed\n",nr,PO[nr].corner)); /*As the polygon nr is deleted the memory allocated for corner must be freed*/ freewh(PO[nr].corner); PO[nr].corner=NULL; for(i=nr;in;i++) Letstru(PO[i],PO[i+1]); /*PO[*t.n] is marked as empty*/ PO[(*t).n].corner=NULL; PO[(*t).n].nc=0; } /****************/ void initclearpolarr(TOTAL *t,int n,int clear) /*works on the polygon array n.. polygons 0 to n-1 are initialized and possibly cleared clear=0... initialization, clear=1.. clearing and initialization*/ { int i; for(i=0;igennr; t->gennr=0; if(genoldnr>0) for(i=0;igen[i].ip].update==0)/*that means this polygon was not updated*/ { Letstru(t->gen[t->gennr],t->gen[i]); (t->gennr)++; } if(genoldnr==0) /*last addpoint was interrupted with volume infinity, thus all generator regions must be recalculated*/ for(i=0;in;i++) PO[i].update=1; /*computes the (approximate) total number of generatorregions using the formular #gen(1 polygon)=2*(n-2) */ t->ngennr=t->gennr; for(i=0;in;i++) if(PO[i].update) t->ngennr+=2*(PO[i].nc-2); if(t->ngennr<=0) t->ngennr=0; else /*allocate the necessary memory for gen[]*/ { if(t->gen==NULL) { t->gen=mallocwh(sizeof(GEN)*t->ngennr); MPR1(("for t->gen %d byte allocated at adress %p\n",sizeof(GEN)*t->ngennr , t->gen)); } else { t->gen=reallocwh(t->gen,sizeof(GEN)*t->ngennr); MPR1(("for t->gen %d byte REallocated at adress %p\n",sizeof(GEN)*t->ngennr , t->gen)); } } /*triangulate changed polygons(automatically the generator regions of these polygons are stored in gen)*/ for(i=0;in;i++) if(PO[i].update) if(!triangulate(&(PO[i]),t->gen,&(t->gennr),t->ngennr,output)) { t->gennr=0;/*shows that the generator regions were not fully computed*/ return(0); } return(1); } /************************************/ void compute_cumvol_guidetable(TOTAL *t) { double volsh=0.; int i,j; /*cumulated volumes are computed*/ volsh=0.; for(i=0;i<(*t).gennr;i++) { volsh+=(*t).gen[i].vol; (*t).gen[i].vols=volsh; } /*we use the cummulated volume to check if there are apparently numerical ptoblems when computing the volumes*/ if(volsh>t->oldhatvol*(1.+1.e-14))/*USE eps instead of constant*/ { printf("Caution!! new volume below the hat was larger than before1\n"); printf("Perhaps it is better to use fewer touching points!"); } t->oldhatvol=volsh; /*computation of guide-table*/ (*t).q[0]=0; i=0; /*the size of the guide table is proportional to (*t).n but during the setup a smaller guide table is used */ if((*t).n==(*t).nn) (*t).m=(*t).mm; else (*t).m=CFGTSU*(*t).mm; for(j=1;j<(*t).m;j++) { while((*t).gen[i].vols/(*t).gen[(*t).gennr-1].vols<(double)j/(*t).m) i++; (*t).q[j]=i; } } /************************************/ int adddesignpoints(TOTAL *t,int output) /*computes for one or more newly added points of contact the polygones and generator regions, and the necessary updates for the old polygones function returns 1 if everything is ok 0 if volume below the hat unbounded occurs output 0.. absolut kein output*/ { int i,j,nalt=0,initflag=0; long ok; double hger[3]; for(i=0;i < t->n;i++) if(PO[i].corner==NULL) initpol(t,i); else { nalt++; PO[i].update=0;/*up to now no update necessary*/ } /** for the new touching points (it is assumed that the old points are ordered at the beginning of the array) all possible lines (projections of the cutting lines of the tangentialebene with other tangential planes) are computed and possible updates of polygons are made*/ do { initflag=0; for(i=t->n-1;i>=nalt;i--) for(j=0;jn;i++) initpol(t,i); } }while(initflag); #if (OUTPUT>=1) printf("Polygon-update finished!\n"); for(i=0;i<(*t).n;i++) printpoly(&((*t).p[i])); #endif ok=compute_tgen(t,output);/*trinagulates all polygons and computes the generator regions, that are stored in t-gen*/ if(!ok)return(0); #if (OUTPUT>=1) for(i=0;i<(*t).gennr;i++) {printf("%d:",i);printgen(&((*t).gen[i]));} #endif compute_cumvol_guidetable(t); #if (OUTPUT>=1) for(i=0;i<(*t).gennr;i++) {printf("%d:",i);printgen(&((*t).gen[i]));} #endif OPR0(("\n points %d, volume below hat %e,(*t).gennr %d \n",(*t).n,(*t).gen[(*t).gennr-1].vols,(*t).gennr)); return(1); }/*end of adddesignpoints*/ /********************************************************/ void freepolarr(TOTAL *t) { int i; for(i=0;i<=(*t).nn;i++) if(t->p[i].corner!=NULL) { MPR1(("corners von po[%d] at adress %p are freed\n",i,t->p[i].corner)); freewh(t->p[i].corner); t->p[i].corner=NULL; } MPR1(("the polygon array t->p adress %p is freed\n",t->p)); freewh(t->p); t->p=NULL; } /****************/ void freesetup(void *t) /*frees all memory that was assigned during the set-up*/ { TOTAL *t1; if(t!=NULL) { t1=t; if(t1->p!=NULL) freepolarr(t1); MPR1(("structure of generator regions t->gen adress %p is freed\n",t1->gen)); freewh(t1->gen); t1->gen=NULL; MPR1(("structure of guide table t->q adress %p is freed\n",t1->q)); freewh(t1->q); t1->q=NULL; MPR1(("array param at adress %p is freed\n",t1->param)); freewh(t1->param); t1->param=NULL; MPR1(("structure TOTAL at adress %p is freed\n",t1)); freewh(t1); t1=NULL; #if (MEMCONTROL>=1&&OUTPUT>=0) printheap(); #endif } } /************************************/ void sample2d(double res[2],void *t1); /****************/ void *setup(int nrsp,int nn,double sp[][2],int nruneq,double uneq[][3], int nrauxun, double (*h)(double x[2],double param[]), double (*hx)(double x[2],double param[]), double (*hy)(double x[2],double param[]),double param[],int nparam) /* double sp[][2] .. array of the starting points, uneq[][3]..array of the equalities defining the domain of the distribution equalities are stored in the form 0 <= uneq[0][0] + uneq[0][1]*x + uneq[0][2]*y domain is defined by these equalities toegether with the starting points, which define, which of the different regions defined by the inequalities is the domain (that, where the starting points lie inside.) !it is assumed that all starting points in sp[] are inside the same region! nrsp number of starting points nn maximal number of points of contact nruneq number of unequalities for domain of distribution behind these inequalities it can be necessary to store inequalities for an auxiliary domain (which must be a subset of the original domain) nrauxun .. total number of inequalities for auxiliary domain we must have nrauxun >=nruneq nrauxun== nruneq means that no auxiliary domain is provided (*h) log of density function of 2-dim. distribution (can have an arbitrary number of double parameters) (*hx) partial derivative of h with respect to the first variable (*hy) partial derivative of h with respect to the second variable param[] contains the current setting of the parameters of the distribution nparam...... number of parameters of the distribution */ { TOTAL *t,*taux; int i,ok,hic; double help[2]; if(nnnn=nn; t->n=nrsp; t->mm=CFGT*nn; t->ngennr=0; t->gennr=0; t->h=h; t->hx=hx; t->hy=hy; t->param=mallocwh(sizeof(double)*nparam); memcpy(t->param,param,sizeof(double)*nparam); MPR1(("for t->param %d doubles (%d bytes) at adress %p\n",nparam,sizeof(double)*nparam,t->param)); t->oldhatvol=VOLHATMAX; t->gen=NULL; t->q=mallocwh(sizeof(int)*(t->mm));/*for guide table*/ MPR1(("for t->q %d ints (%d bytes) at adress %p\n",t->mm,sizeof(int)*(t->mm),t->q)); t->p=mallocwh(sizeof(POLYGON)*(nn+1))/*in p[nn] the start polygon is stored*/; MPR1(("for t->p %d polygons (%d bytes) at adress %p\n",nn+1,sizeof(POLYGON)*(nn+1),t->p)); initclearpolarr(t,nn+1,0); for(i=0;inn); Arrlet2(PO[nn].design,sp[0]); for(i=0;inruneq) { OPR0(("as volume unbounded setup is restarted using the auxiliary domain\n")); taux=setup(nrsp,nn,sp,nrauxun,uneq,nrauxun,h,hx,hy,t->param,nparam); if(taux!=NULL) { /*now we sample from the distribution over the auxiliary domain to get starting points for the original domain*/ sample2d(help,taux); OPR1(("first auxiliary sampling finished taux-n %d t-n %d\n",taux->n,t->n)); do{ if(taux->n==nn)/*in this case the polygon structure was dismissed*/ { printf("Error!!setup\n"); printf("too few points allowed for auxiliary domain method\n"); freesetup(taux); freesetup(t); return(NULL); } for(i=0;in;i++) Arrlet2(t->p[i].design,taux->p[i].design); OPR1(("touching points transferred\n")); /*generator regions and polygons are cleared to make shure that they are recalculated with the new touching points*/ t->gennr=0; initclearpolarr(t,t->nn,1); OPR1(("clearing finished po[nn]\n")); /*printpoly(&(t->p[t->nn]));*/ t->n=taux->n; ok=adddesignpoints(t,0); OPR1(("adddesignpoints(t) finished\n")); if(!ok) for(hic=0;(taux->nn+1)&&(hic<1000);hic++) { sample2d(help,taux); OPR1(("auxiliary sampling finished taux-n %d t-n %d\n",taux->n,t->n)); } if(hic==1000) { printf("Error!!setup, no additional design points are foun in auxiliary domain\n"); printf("try to use larger auxiliary domain\n"); freesetup(taux); freesetup(t); return(NULL); } }while(!ok); OPR0(("Starting values found, auxiliary domain points dismissed\n")); freesetup(taux); } } /*Controll if computation of has worked correctly*/ if (ok==0||t->gen[t->gennr-1].vols>=VOLHATMAX || t->gen[t->gennr-1].vols<0.) { printf("ErrorSETUP: volume below the hat unbounded!!\n"); printf("you can try:\n"); printf("a) to provide an auxilliary domain before starting set-up\n"); printf("b) use an auxiliary domain and increas nn \n"); printf("c) use more (or better) starting values\n"); freesetup(t); return(0); } else { t->ok = OK; return(t); } } #undef PO #if(MEMCONTROL>=1) static unsigned long int widcount=0; void setwidcount(unsigned long int newval) /*sets the value of widcount*/ { widcount=newval;} unsigned long int returnwidcount(void) /*returns the current value of widcount*/ { return(widcount);} #else /*the functions are defined but return -1 if MEMCONTROLL==1*/ void setwidcount(unsigned long int newval) { printf("Caution! function returnwidcount does not work correctly as MEMCONTROL not >= 1\n"); } unsigned long int returnwidcount(void) { printf("Caution! function returnwidcount does not work correctly as MEMCONTROL not >= 1\n"); return(-1);} #endif /****************/ void sample2d(double res[2],void *t1) /*samples from the 2-dimensional distribution, result is stored in res[2] t1 is the pointer that was returned by setup*/ { TOTAL *t; int j; double zwert,u,v,helpp[2][2]; t=t1; if(t==NULL) { printf("sample2d: es wurde Nullpointer uebergeben!!\n"); return;} else if(t->ok!=OK) { printf("sample2d: es wurde falscher pointer uebergeben!!\n"); return;} /*sampling routine*/ while(1) { /*first the random index j of the region is generated using the guide table*/ u=UNIFORM(); #if (MEMCONTROL>=1) widcount++; #endif j=t->q[(int)(t->m*u)]; u*=t->gen[t->gennr-1].vols; while(t->gen[j].volsgen[j].open<=1) { res[0]=gg2inter(t->gen[j].a,t->gen[j].b); res[1]= UNIFORM()*res[0]*(t->gen[j].k[1]-t->gen[j].k[0])+ res[0]*t->gen[j].k[0]; } else { res[0]=marg2(&(t->gen[j])); res[1]=UNIFORM()*(t->gen[j].b+res[0]*(t->gen[j].k[1]-t->gen[j].k[0]))+ res[0]*t->gen[j].k[0]; } zwert=t->gen[j].s+res[0]*t->gen[j].a; invrotate(t->gen[j].rot,res,helpp[0]); Arradd2(res,helpp[0],t->gen[j].point); /* the acceptance condition is tested*/ v=log(UNIFORM()); if (v<=(t->h)(res,t->param)-zwert) return; /*in the case of rejection if the number of touching points is not exausted the rejected pair is taken as new touching point*/ else if(t->nnn) { t->n++; t->p[t->n-1].design[0]=res[0]; t->p[t->n-1].design[1]=res[1]; /*adddesignpoints is computed only about ??? times*/ if(t->n<10||(t->n%(t->nn/SUNC))==0||t->n==t->nn) adddesignpoints(t,1); /*array of polygons is deleted as the number of points of contact was reached*/ if((t->n==t->nn) &&(t->p!=NULL)) { freepolarr(t); OPR0(("Polygon structure was freed\n")); } } } } SHAR_EOF fi # end of overwriting check cd .. cd .. cd .. if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'readme' then echo shar: will not over-write existing file "'readme'" else cat << "SHAR_EOF" > 'readme' /* ALC2DLIB-Version 1.0 implemented by WH 22.3.99, all rights reserved */ /*please report problems or bugs to whoer@statistik.wu-wien.ac.at */ /* or hormannw@boun.edu.tr */ /*This file readme*/ Random variate generation software ALC2DLIB can generate random pairs from arbitrary log-concave 2-dimensional distributions. It is the implementation of the Algorithm ALC2D described in the paper: "Wolfgang H"ormann: An automatic generator for bivariate log-concave distributions" The software contains the follwing files: readme ... this file alc2d.c ... C-source of the algorithm alc2d.h ... header file also containing explanations for the use of the functions main1.c ... simple demonstration of the use of the algorithm for the standard normal distribution main2.c ... demonstration of the use of the algorithm for several distributions mainnort.c ... Chi2-test for the bivariate normal distribution generated with Algorithm ALC2D The software was developped with the GNU-C-compiler under LINUX. It was tested on a SUN-Workstation using the GNU-C-compiler and under ULTRIX using the ULTRIX C-compiler as-well. SHAR_EOF fi # end of overwriting check cd .. # End of shell archive exit 0