mas01mc@285
|
1 #include <sys/types.h>
|
mas01mc@285
|
2 #include <sys/stat.h>
|
mas01mc@285
|
3 #include <fcntl.h>
|
mas01mc@285
|
4 #include <math.h>
|
mas01mc@285
|
5 #include <stdlib.h>
|
mas01mc@285
|
6 #include <stdio.h>
|
mas01mc@285
|
7
|
mas01mc@285
|
8 double randn();
|
mas01mc@285
|
9 double randbl();
|
mas01mc@285
|
10
|
mas01mc@285
|
11 /* genpoints count radius^2 */
|
mas01mc@285
|
12 int main(int argc, char *argv[]) {
|
mas01mc@285
|
13 if (argc < 3) {
|
mas01mc@285
|
14 fprintf(stderr, "usage: %s count radius^2 [dim]\n", argv[0]);
|
mas01mc@285
|
15 exit(1);
|
mas01mc@285
|
16 }
|
mas01mc@285
|
17 long int count = strtol(argv[1], NULL, 0);
|
mas01mc@285
|
18 double rsquared = strtod(argv[2], NULL);
|
mas01mc@285
|
19 long int dim = 3;
|
mas01mc@285
|
20 if(argc > 3)
|
mas01mc@285
|
21 dim = strtol(argv[3], NULL, 0);
|
mas01mc@285
|
22
|
mas01mc@285
|
23 // Generate *count* Gaussian Random vectors in R^*dim*
|
mas01mc@285
|
24 // sitting on the *rdashed*-sphere
|
mas01mc@285
|
25
|
mas01cr@303
|
26 srandom(time(NULL));
|
mas01mc@285
|
27
|
mas01mc@285
|
28 int i,j;
|
mas01mc@285
|
29 for (i = 0; i < count + 1; i++) {
|
mas01mc@285
|
30 // Normed Gaussian random vectors are distributed uniformly on unit sphere
|
mas01mc@285
|
31 double* coords = malloc(dim * sizeof(double));
|
mas01mc@285
|
32 double nmsq = 0.0;
|
mas01mc@285
|
33
|
mas01mc@285
|
34 for (j = 0; j < dim; j++){
|
mas01mc@285
|
35 if(i < count)
|
mas01mc@285
|
36 coords[j] = randn();
|
mas01mc@285
|
37 else
|
mas01mc@285
|
38 coords[j] = 0.0;
|
mas01mc@285
|
39 nmsq += coords[j]*coords[j];
|
mas01mc@285
|
40 }
|
mas01mc@285
|
41
|
mas01mc@285
|
42 double nm2 = 0.0;
|
mas01mc@285
|
43 if(i < count){
|
mas01mc@285
|
44 nm2 = sqrt(rsquared/nmsq);
|
mas01mc@285
|
45 // Place on rdash-sphere
|
mas01mc@285
|
46 for (j = 0; j < dim; j++)
|
mas01mc@285
|
47 coords[j] *= nm2;
|
mas01mc@285
|
48 }
|
mas01mc@285
|
49 // Translate to (0,0,...,1)
|
mas01mc@285
|
50 coords[dim-1]+=1.0;
|
mas01mc@285
|
51
|
mas01mc@286
|
52 // Compute distance to (0,0,...,1)
|
mas01mc@285
|
53 nmsq = 0.0;
|
mas01mc@286
|
54 for (j = 0; j < dim-1; j++){
|
mas01mc@285
|
55 nmsq += coords[j]*coords[j];
|
mas01mc@285
|
56 }
|
mas01mc@285
|
57 // Save last value to distance calulcation to query(0,0,...,1)
|
mas01mc@285
|
58 double nth = coords[dim-1];
|
mas01mc@285
|
59 // Output to ASCII terminal
|
mas01mc@285
|
60 printf("(");
|
mas01mc@285
|
61 for(j = 0; j < dim; j++)
|
mas01mc@285
|
62 printf("%8.3f ", coords[j]);
|
mas01mc@286
|
63 printf(") d = %8.3f\n", sqrt(nmsq + (nth-1)*(nth-1)));
|
mas01mc@285
|
64
|
mas01mc@285
|
65
|
mas01mc@285
|
66 // Save single feature vector
|
mas01mc@285
|
67 char name[40];
|
mas01mc@285
|
68 if(i < count)
|
mas01mc@290
|
69 snprintf(name, 39, i<10?"testfeature0%d":"testfeature%d", i);
|
mas01mc@285
|
70 else
|
mas01mc@285
|
71 snprintf(name, 39, "queryfeature");
|
mas01mc@285
|
72 /* assumes $PWD is right */
|
mas01mc@285
|
73 int fd = open(name, O_CREAT|O_TRUNC|O_WRONLY, S_IWUSR|S_IRUSR|S_IRGRP|S_IROTH);
|
mas01mc@285
|
74
|
mas01mc@285
|
75 write(fd, &dim, sizeof(int));
|
mas01mc@285
|
76 for(j = 0; j < dim; j++)
|
mas01mc@285
|
77 write(fd, coords + j, sizeof(double));
|
mas01mc@285
|
78 close(fd);
|
mas01mc@285
|
79
|
mas01mc@285
|
80 free(coords);
|
mas01mc@285
|
81 }
|
mas01mc@285
|
82 exit(0);
|
mas01mc@285
|
83 }
|
mas01mc@285
|
84
|
mas01mc@285
|
85 // Genereate U[0,1]
|
mas01mc@285
|
86 double randbl(){
|
mas01mc@285
|
87 return ( (double)rand() / ((double)(RAND_MAX)+(double)(1)) );
|
mas01mc@285
|
88 }
|
mas01mc@285
|
89
|
mas01mc@285
|
90 // Generate z ~ N(0,1)
|
mas01mc@285
|
91 double randn(){
|
mas01mc@285
|
92 // Box-Muller
|
mas01mc@285
|
93 double x1, x2;
|
mas01mc@285
|
94 do{
|
mas01mc@285
|
95 x1 = randbl();
|
mas01mc@285
|
96 } while (x1 == 0); // cannot take log of 0
|
mas01mc@285
|
97 x2 = randbl();
|
mas01mc@285
|
98 double z = sqrt(-2.0 * log(x1)) * cos(2.0 * M_PI * x2);
|
mas01mc@285
|
99 return z;
|
mas01mc@285
|
100 }
|
mas01mc@285
|
101
|