-
Notifications
You must be signed in to change notification settings - Fork 4
/
d3sphere.c
85 lines (67 loc) · 2.49 KB
/
d3sphere.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
#include "header.h"
#include "macro.h"
void sphere(char *filename)
{
const counter no_obs=20, no_pts=4000, edge=1000;
const real ratio=edge/UNIMAX;
int fcmpr(const void *u1, const void *u2);
register counter i, j, k;
real d, dmin;
real *x, *y, *z, *p, r3, pvalue;
puts("\n\t|-------------------------------------------------------------|");
puts("\t| THE 3DSPHERES TEST |");
puts("\t|Choose 4000 random points in a cube of edge 1000. At each |");
puts("\t|point, center a sphere large enough to reach the next closest|");
puts("\t|point. Then the volume of the smallest such sphere is (very |");
puts("\t|close to) exponentially distributed with mean 120pi/3. Thus |");
puts("\t|the radius cubed is exponential with mean 30. (The mean is |");
puts("\t|obtained by extensive simulation). The 3DSPHERES test gener-|");
puts("\t|ates 4000 such spheres 20 times. Each min radius cubed leads|");
puts("\t|to a uniform variable by means of 1-exp(-r^3/30.), then a |");
puts("\t| KSTEST is done on the 20 p-values. |");
puts("\t|-------------------------------------------------------------|\n");
printf("\t\t The 3DSPHERES test for file %s\n\n", filename);
printf("\t\tsample no\tr^3\t\tequiv. uni.\n");
x=(real*)malloc(no_pts*sizeof(real));
y=(real*)malloc(no_pts*sizeof(real));
z=(real*)malloc(no_pts*sizeof(real));
p=(real*)malloc(no_obs*sizeof(real));
for(i=1; i<=no_obs; ++i){
dmin=10000000.;
for(j=0; j<no_pts; ++j){
x[j]=ratio*uni(filename);
y[j]=ratio*uni(filename);
z[j]=ratio*uni(filename);
}
qsort(x, no_pts, sizeof(real), fcmpr);
/* to compare with fortran counter part
for(j=0; j<no_pts; ++j){
y[j]=ratio*uni(filename);
z[j]=ratio*uni(filename);
}*/
for(j=0; j<no_pts-1; ++j){
for(k=j+1; k<no_pts; ++k){
d=(x[k]-x[j])*(x[k]-x[j]);
if(d>=dmin) break;
d+=(y[k]-y[j])*(y[k]-y[j])+(z[k]-z[j])*(z[k]-z[j]);
dmin=MIN(dmin, d);
}
}
r3=dmin*sqrt(dmin);
p[i-1]=1-exp(-MIN(r3/30., 20));
printf("\t\t %lu\t\t%.3f\t\t%f\n", i, r3, p[i-1]);
}
uni("close");
free(x); free(y); free(z);
puts("\t--------------------------------------------------------------");
pvalue=KStest(p,no_obs);
printf("\t\tp-value for KS test on those %lu p-values: %f", no_obs, pvalue);
puts("\n");
free(p);
return;
}
/*main()
{
sphere("binc");
return;
}*/