comparison sample.cpp @ 280:3be15407e814

Merge sampling branch (-r361:405, though I hope that the branch is now
finished) onto trunk.  API developers take note.

Things still to clear up:
* whether the threshold distance it currently reports bears any relation
  to reality;
* if not, how to bring it a bit more into alignment;
* minor code cleanup issues in sample.cpp;
* incorporating --absolute-threshold handling into sampling;
* writing suitable test cases.

+#include "audioDB.h"
+
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_rng.h>
+
+static
+double yfun(double d) {
+  return gsl_sf_log(d) - gsl_sf_psi(d);
+}
+
+static
+double yinv(double y) {
+  double a = 1.0e-5;
+  double b = 1000.0;
+
+  double ay = yfun(a);
+  double by = yfun(b);
+
+  double c, cy;
+
+  /* FIXME: simple binary search; there's probably some clever solver
+     in gsl somewhere which is less sucky. */
+  while ((b - a) > 1.0e-5) {
+    c = (a + b) / 2;
+    cy = yfun(c);
+    if (cy > y) {
+      a = c;
+      ay = cy;
+    } else {
+      b = c;
+      by = cy;
+    }
+  }
+
+  return c;
+}
+
+unsigned audioDB::random_track(unsigned *propTable, unsigned total, gsl_rng *rng) {
+  /* FIXME: make this O(1) by using the alias-rejection method, or
+     some other sensible method of sampling from a discrete
+     distribution. */
+  double thing = gsl_rng_uniform(rng);
+  unsigned sofar = 0;
+  for (unsigned int i = 0; i < dbH->numFiles; i++) {
+    sofar += propTable[i];
+    if (thing < ((double) sofar / (double) total)) {
+      return i;
+    }
+  }
+  error("fell through in random_track()");
+
+  /* FIXME: decorate error's declaration so that this isn't necessary */
+  return 0;
+}
+
+void audioDB::sample(const char *dbName) {
+  initTables(dbName, 0);
+
+  gsl_rng *rng = gsl_rng_alloc(gsl_rng_mt19937);
+
+  /* FIXME: in Real Life we'll want to initialize the RNG using
+     /dev/random or the current time or something, like this:
+
+     unsigned int seed;
+     int fd = open("/dev/urandom", O_RDONLY);
+     read(fd, &seed, 4);
+     
+     gsl_rng_set(rng, seed);
+  */
+
+  // build track offset table (FIXME: cut'n'pasted from query.cpp)
+  off_t *trackOffsetTable = new off_t[dbH->numFiles];
+  unsigned cumTrack=0;
+  for(unsigned int k = 0; k < dbH->numFiles; k++){
+    trackOffsetTable[k] = cumTrack;
+    cumTrack += trackTable[k] * dbH->dim;
+  }
+
+  unsigned *propTable = new unsigned[dbH->numFiles];
+  unsigned total = 0;
+  unsigned count = 0;
+
+  for (unsigned int i = 0; i < dbH->numFiles; i++) {
+    /* what kind of a stupid language doesn't have binary max(), let
+       alone nary? */
+    unsigned int prop = trackTable[i] - sequenceLength + 1;
+    prop = prop > 0 ? prop : 0;
+    if (prop > 0) 
+      count++;
+    propTable[i] = prop;
+    total += prop;
+  }
+
+  if (total == 0) {
+    error("no sequences of this sequence length in the database", dbName);
+  }
+
+  unsigned int vlen = dbH->dim * sequenceLength;
+  double *v1 = new double[vlen];
+  double *v2 = new double[vlen];
+  double v1norm, v2norm, v1v2;
+
+  double sumdist = 0;
+  double sumlogdist = 0;
+
+  for (unsigned int i = 0; i < nsamples;) {
+    unsigned track1 = random_track(propTable, total, rng);
+    unsigned track2 = random_track(propTable, total, rng);
+
+    if(track1 == track2)
+      continue;
+
+    unsigned i1 = gsl_rng_uniform_int(rng, propTable[track1]);
+    unsigned i2 = gsl_rng_uniform_int(rng, propTable[track2]);
+
+    VERB_LOG(1, "%d %d, %d %d | ", track1, i1, track2, i2);
+
+    /* FIXME: this seeking, reading and distance calculation should
+       share more code with the query loop */
+    lseek(dbfid, dbH->dataOffset + trackOffsetTable[track1] * sizeof(double) + i1 * dbH->dim * sizeof(double), SEEK_SET);
+    read(dbfid, v1, dbH->dim * sequenceLength * sizeof(double));
+
+    lseek(dbfid, dbH->dataOffset + trackOffsetTable[track2] * sizeof(double) + i2 * dbH->dim * sizeof(double), SEEK_SET);
+    read(dbfid, v2, dbH->dim * sequenceLength * sizeof(double));
+
+    v1norm = 0;
+    v2norm = 0;
+    v1v2 = 0;
+
+    for (unsigned int j = 0; j < vlen; j++) {
+      v1norm += v1[j]*v1[j];
+      v2norm += v2[j]*v2[j];
+      v1v2 += v1[j]*v2[j];
+    }
+
+    /* FIXME: we must deal with infinities better than this; there
+       could be all sorts of NaNs from arbitrary features.  Best
+       include power thresholds or something... */
+    if(isfinite(v1norm) && isfinite(v2norm) && isfinite(v1v2)) {
+
+      VERB_LOG(1, "%f %f %f | ", v1norm, v2norm, v1v2);
+      /* assume normalizedDistance == true for now */
+      /* FIXME: not convinced that the statistics we calculated in
+	 TASLP paper are technically valid for normalizedDistance */
+
+      double dist = 2 - 2 * v1v2 / sqrt(v1norm * v2norm);
+      // double dist = v1norm + v2norm - 2*v1v2;
+      
+      VERB_LOG(1, "%f %f\n", dist, log(dist));
+      sumdist += dist;
+      sumlogdist += log(dist);
+      i++;
+    } else {
+      VERB_LOG(1, "infinity/NaN found: %f %f %f\n", v1norm, v2norm, v1v2);
+    }
+  }
+
+  /* FIXME: the mean isn't really what we should be reporting here */
+  unsigned meanN = total / count;
+
+  double sigma2 = sumdist / (sequenceLength * dbH->dim * nsamples);
+  double d = 2 * yinv(log(sumdist/nsamples) - sumlogdist/nsamples);
+
+  std::cout << "Summary statistics" << std::endl;
+  std::cout << "number of samples: " << nsamples << std::endl;
+  std::cout << "sum of distances (S): " << sumdist << std::endl;
+  std::cout << "sum of log distances (L): " << sumlogdist << std::endl;
+
+  /* FIXME: we'll also want some more summary statistics based on
+     propTable, for the minimum-of-X estimate */
+  std::cout << "mean number of applicable sequences (N): " << meanN << std::endl;
+  std::cout << std::endl;
+  std::cout << "Estimated parameters" << std::endl;
+  std::cout << "sigma^2: " << sigma2 << "; ";
+  std::cout << "Msigma^2: " << sumdist / nsamples << std::endl;
+  std::cout << "d: " << d << std::endl;
+
+  double logw = (2 / d) * gsl_sf_log(-gsl_sf_log(0.99));
+  double logxthresh = gsl_sf_log(sumdist / nsamples) + logw
+    - (2 / d) * gsl_sf_log(meanN)
+    - gsl_sf_log(d/2)
+    - (2 / d) * gsl_sf_log(2 / d)
+    + (2 / d) * gsl_sf_lngamma(d / 2);
+
+  std::cout << "track xthresh: " << exp(logxthresh) << std::endl;
+
+  delete[] propTable;
+  delete[] v1;
+  delete[] v2;
+}