view tests/pointset_test/genpoints2.c @ 395:bc7a821004bb api-inversion

Invert audioDB::status / audiodb_status().

To do that without breaking abstractions, we actually need a new field
in the status structure, storing the size of the data region.
Previously, this was computed in the audioDB::status request from the
database header, but I'm assuming that "user" code doesn't have access
to such internals.  While we're at it, name some intermediate values in
audioDB::status() so that I don't get confused.

Here's the thing, though: we need to make sure that the adb_t * that we
have from audiodb_open() or audiodb_create() is propagated all the way
through into the C++ routines that implement library functions -- in
particular those which actually write to the database; otherwise we
won't have a consistent view in memory of the header on-disk (as the adb
header that will have been written to disk won't be the same as the one
in memory).

We can do that, by altering the "API" audioDB constructors to take the
adb_t * argument, and setting the adb field in the audioDB object that
we've already introduced to that.  But now we need to be careful a
couple of times: if we have one, then audioDB::initTables() mustn't
stomp on it; also, if we're only constructing an audioDB instance to
fulfil an API request, we mustn't audiodb_close() the one we have when
we destroy the audioDB object, because the adb_t * is the one we have
passed in and are going to reuse in later calls to the API.

The good news is that we can be careful in just these ways with minimal
code.  The really good news is that once the inversion is complete, all
of this horribleness will automatically go away (as there will be no
code which constructs audioDB objects to fulfil API functions).  Hooray!
It's almost like it was all planned this way.

#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

double randn();
double randbl();

/* genpoints count radius^2 */
int main(int argc, char *argv[]) {
  if (argc < 3) {
    fprintf(stderr, "usage: %s count radius^2 [dim]\n", argv[0]);
    exit(1);
  }
  long int count = strtol(argv[1], NULL, 0);
  double rsquared = strtod(argv[2], NULL);
  long int dim = 3;
  if(argc > 3)
    dim = strtol(argv[3], NULL, 0);
  
  // Generate *count* Gaussian Random vectors in R^*dim*
  // sitting on the *rdashed*-sphere

  srandom(time(NULL));

  int i,j;
  for (i = 0; i < count + 1; i++) {
    // Normed Gaussian random vectors are distributed uniformly on unit sphere
    double* coords = malloc(dim * sizeof(double));
    double nmsq = 0.0;

    for (j = 0; j < dim; j++){
      if(i < count)
	coords[j] = randn();
      else
	coords[j] = 0.0;
      nmsq += coords[j]*coords[j];
    }

    double nm2 = 0.0;
    if(i < count){
      nm2 = sqrt(rsquared/nmsq);
      // Place on rdash-sphere
      for (j = 0; j < dim; j++)
	coords[j] *= nm2;
    }
    // Translate to (0,0,...,1)
    coords[dim-1]+=1.0; 

    // Compute distance to (0,0,...,1)
    nmsq = 0.0;
    for (j = 0; j < dim-1; j++){
      nmsq += coords[j]*coords[j];
    }
    // Save last value to distance calulcation to query(0,0,...,1)
    double nth = coords[dim-1];
    // Output to ASCII terminal
    printf("(");
    for(j = 0; j < dim; j++)
      printf("%8.3f ", coords[j]);
    printf(") d = %8.3f\n", sqrt(nmsq + (nth-1)*(nth-1)));
    

    // Save single feature vector
    char name[40];
    if(i < count)
      snprintf(name, 39, i<10?"testfeature0%d":"testfeature%d", i);
    else
      snprintf(name, 39, "queryfeature");
    /* assumes $PWD is right */
    int fd = open(name, O_CREAT|O_TRUNC|O_WRONLY, S_IWUSR|S_IRUSR|S_IRGRP|S_IROTH);

    write(fd, &dim, sizeof(int));
    for(j = 0; j < dim; j++)
      write(fd, coords + j, sizeof(double));
    close(fd);

    free(coords);
  }
  exit(0);
}

// Genereate U[0,1]
double randbl(){
  return (   (double)rand() / ((double)(RAND_MAX)+(double)(1)) );
}

// Generate z ~ N(0,1)
double randn(){
// Box-Muller
  double x1, x2;
  do{
    x1 = randbl();
  } while (x1 == 0); // cannot take log of 0
  x2 = randbl();
  double z = sqrt(-2.0 * log(x1)) * cos(2.0 * M_PI * x2);
  return z;
}