mas01mc@285: #include <sys/types.h>
mas01mc@285: #include <sys/stat.h>
mas01mc@285: #include <fcntl.h>
mas01mc@285: #include <math.h>
mas01mc@285: #include <stdlib.h>
mas01mc@285: #include <stdio.h>
mas01mc@285: 
mas01mc@285: double randn();
mas01mc@285: double randbl();
mas01mc@285: 
mas01mc@285: /* genpoints count radius^2 */
mas01mc@285: int main(int argc, char *argv[]) {
mas01mc@285:   if (argc < 3) {
mas01mc@285:     fprintf(stderr, "usage: %s count radius^2 [dim]\n", argv[0]);
mas01mc@285:     exit(1);
mas01mc@285:   }
mas01mc@285:   long int count = strtol(argv[1], NULL, 0);
mas01mc@285:   double rsquared = strtod(argv[2], NULL);
mas01mc@285:   long int dim = 3;
mas01mc@285:   if(argc > 3)
mas01mc@285:     dim = strtol(argv[3], NULL, 0);
mas01mc@285:   
mas01mc@285:   // Generate *count* Gaussian Random vectors in R^*dim*
mas01mc@285:   // sitting on the *rdashed*-sphere
mas01mc@285: 
mas01cr@303:   srandom(time(NULL));
mas01mc@285: 
mas01mc@285:   int i,j;
mas01mc@285:   for (i = 0; i < count + 1; i++) {
mas01mc@285:     // Normed Gaussian random vectors are distributed uniformly on unit sphere
mas01mc@285:     double* coords = malloc(dim * sizeof(double));
mas01mc@285:     double nmsq = 0.0;
mas01mc@285: 
mas01mc@285:     for (j = 0; j < dim; j++){
mas01mc@285:       if(i < count)
mas01mc@285: 	coords[j] = randn();
mas01mc@285:       else
mas01mc@285: 	coords[j] = 0.0;
mas01mc@285:       nmsq += coords[j]*coords[j];
mas01mc@285:     }
mas01mc@285: 
mas01mc@285:     double nm2 = 0.0;
mas01mc@285:     if(i < count){
mas01mc@285:       nm2 = sqrt(rsquared/nmsq);
mas01mc@285:       // Place on rdash-sphere
mas01mc@285:       for (j = 0; j < dim; j++)
mas01mc@285: 	coords[j] *= nm2;
mas01mc@285:     }
mas01mc@285:     // Translate to (0,0,...,1)
mas01mc@285:     coords[dim-1]+=1.0; 
mas01mc@285: 
mas01mc@286:     // Compute distance to (0,0,...,1)
mas01mc@285:     nmsq = 0.0;
mas01mc@286:     for (j = 0; j < dim-1; j++){
mas01mc@285:       nmsq += coords[j]*coords[j];
mas01mc@285:     }
mas01mc@285:     // Save last value to distance calulcation to query(0,0,...,1)
mas01mc@285:     double nth = coords[dim-1];
mas01mc@285:     // Output to ASCII terminal
mas01mc@285:     printf("(");
mas01mc@285:     for(j = 0; j < dim; j++)
mas01mc@285:       printf("%8.3f ", coords[j]);
mas01mc@286:     printf(") d = %8.3f\n", sqrt(nmsq + (nth-1)*(nth-1)));
mas01mc@285:     
mas01mc@285: 
mas01mc@285:     // Save single feature vector
mas01mc@285:     char name[40];
mas01mc@285:     if(i < count)
mas01mc@290:       snprintf(name, 39, i<10?"testfeature0%d":"testfeature%d", i);
mas01mc@285:     else
mas01mc@285:       snprintf(name, 39, "queryfeature");
mas01mc@285:     /* assumes $PWD is right */
mas01mc@285:     int fd = open(name, O_CREAT|O_TRUNC|O_WRONLY, S_IWUSR|S_IRUSR|S_IRGRP|S_IROTH);
mas01mc@285: 
mas01mc@285:     write(fd, &dim, sizeof(int));
mas01mc@285:     for(j = 0; j < dim; j++)
mas01mc@285:       write(fd, coords + j, sizeof(double));
mas01mc@285:     close(fd);
mas01mc@285: 
mas01mc@285:     free(coords);
mas01mc@285:   }
mas01mc@285:   exit(0);
mas01mc@285: }
mas01mc@285: 
mas01mc@285: // Genereate U[0,1]
mas01mc@285: double randbl(){
mas01mc@285:   return (   (double)rand() / ((double)(RAND_MAX)+(double)(1)) );
mas01mc@285: }
mas01mc@285: 
mas01mc@285: // Generate z ~ N(0,1)
mas01mc@285: double randn(){
mas01mc@285: // Box-Muller
mas01mc@285:   double x1, x2;
mas01mc@285:   do{
mas01mc@285:     x1 = randbl();
mas01mc@285:   } while (x1 == 0); // cannot take log of 0
mas01mc@285:   x2 = randbl();
mas01mc@285:   double z = sqrt(-2.0 * log(x1)) * cos(2.0 * M_PI * x2);
mas01mc@285:   return z;
mas01mc@285: }
mas01mc@285: