/*
 * lsh.c - experiments for the LSH article.
 *
 *   ./lsh test                          self-check of the multi-probe generator
 *   ./lsh curve  TRIALS SEED            hyperplane S-curves and p-stable p(u), Monte Carlo
 *   ./lsh knn    DIR K W L T Q SEED     E2LSH + query-directed multi-probe on SIFT1M
 *
 * All outputs are counts or frequencies; nothing depends on the clock.
 * Build: gcc -O2 -Wall -Wextra -o lsh lsh.c -lm
 */
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/* ---------- xoshiro256** seeded by splitmix64 ---------- */

static uint64_t rs[4];

static uint64_t splitmix64(uint64_t *x)
{
    uint64_t z = (*x += 0x9e3779b97f4a7c15ULL);
    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
    z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
    return z ^ (z >> 31);
}

static void rng_seed(uint64_t seed)
{
    for (int i = 0; i < 4; i++) rs[i] = splitmix64(&seed);
}

static inline uint64_t rotl(uint64_t x, int k) { return (x << k) | (x >> (64 - k)); }

static uint64_t rng_next(void)
{
    uint64_t r = rotl(rs[1] * 5, 7) * 9, t = rs[1] << 17;
    rs[2] ^= rs[0]; rs[3] ^= rs[1]; rs[1] ^= rs[2]; rs[0] ^= rs[3];
    rs[2] ^= t; rs[3] = rotl(rs[3], 45);
    return r;
}

/* uniform in [0, 1) */
static double rng_u01(void) { return (double)(rng_next() >> 11) * 0x1.0p-53; }

/* standard normal, Box-Muller */
static double rng_gauss(void)
{
    static int have;
    static double spare;
    if (have) { have = 0; return spare; }
    double u1 = 1.0 - rng_u01(), u2 = rng_u01(); /* u1 in (0, 1] */
    double r = sqrt(-2.0 * log(u1)), t = 2.0 * M_PI * u2;
    spare = r * sin(t); have = 1;
    return r * cos(t);
}

static void *xmalloc(size_t n)
{
    void *p = malloc(n);
    if (!p) { fprintf(stderr, "out of memory (%zu bytes)\n", n); exit(1); }
    return p;
}

static double dot(const double *a, const double *b, int d)
{
    double s = 0;
    for (int i = 0; i < d; i++) s += a[i] * b[i];
    return s;
}

/* random unit vector u, and unit e orthogonal to u */
static void random_frame(double *u, double *e, int d)
{
    for (int i = 0; i < d; i++) { u[i] = rng_gauss(); e[i] = rng_gauss(); }
    double nu = sqrt(dot(u, u, d));
    for (int i = 0; i < d; i++) u[i] /= nu;
    double pr = dot(u, e, d);
    for (int i = 0; i < d; i++) e[i] -= pr * u[i];
    double ne = sqrt(dot(e, e, d));
    for (int i = 0; i < d; i++) e[i] /= ne;
}

/* ---------- mode curve ---------- */

#define CD 128          /* dimension used for the Monte Carlo */
#define NANG 37         /* angles 0, 5, ..., 180 degrees */
#define MAXH 256        /* hyperplanes drawn per trial */
#define NU 41           /* p-stable: distance/w = 0, 0.1, ..., 4.0 */
#define PS_REP 20       /* p-stable samples per trial */

static const int CFG_K[] = {1, 4, 8, 8};
static const int CFG_L[] = {1, 4, 8, 32};
#define NCFG 4

/* Datar et al. 2004 with a Gaussian projection: X = a.(p - q) ~ N(0, u^2),
 * p(u) = int_0^w (2/u) phi(t/u) (1 - t/w) dt, closed form with s = w/u. */
static double pstable_p(double u, double w)
{
    if (u == 0) return 1.0;
    double s = w / u;
    return 1.0 - erfc(s / M_SQRT2) - 2.0 / (sqrt(2.0 * M_PI) * s) * (1.0 - exp(-s * s / 2));
}

static int mode_curve(long trials, uint64_t seed)
{
    static double u[CD], e[CD], r[CD], pu[MAXH], pe[MAXH];
    static long single[NANG], hit[NCFG][NANG], ps_hit[NU];
    rng_seed(seed);
    for (long tr = 0; tr < trials; tr++) {
        random_frame(u, e, CD);
        for (int j = 0; j < MAXH; j++) {
            for (int i = 0; i < CD; i++) r[i] = rng_gauss();
            pu[j] = dot(r, u, CD);
            pe[j] = dot(r, e, CD);
        }
        for (int a = 0; a < NANG; a++) {
            /* v = cos(t) u + sin(t) e is a unit vector at angle t from u */
            double t = a * 5.0 * M_PI / 180.0, c = cos(t), s = sin(t);
            unsigned char eq[MAXH];
            for (int j = 0; j < MAXH; j++) {
                eq[j] = (pu[j] >= 0) == (c * pu[j] + s * pe[j] >= 0);
                single[a] += eq[j];
            }
            for (int g = 0; g < NCFG; g++) {
                int k = CFG_K[g], any = 0;
                for (int l = 0; l < CFG_L[g] && !any; l++) {
                    int all = 1;
                    for (int j = l * k; j < (l + 1) * k; j++) all &= eq[j];
                    any |= all;
                }
                hit[g][a] += any;
            }
        }
        /* p-stable, w = 1: h(x) = floor(a.x + b), q = p + dist * e */
        for (int rep = 0; rep < PS_REP; rep++) {
            for (int i = 0; i < CD; i++) r[i] = rng_gauss();
            double b = rng_u01(), sp = 3.0 * dot(r, u, CD), se = dot(r, e, CD);
            for (int m = 0; m < NU; m++) {
                double dist = 0.1 * m;
                ps_hit[m] += floor(sp + b) == floor(sp + dist * se + b);
            }
        }
    }
    printf("# hyperplane LSH, d=%d, trials=%ld, seed=%llu\n", CD, trials,
           (unsigned long long)seed);
    printf("angle_deg,p_mc,p_theory");
    for (int g = 0; g < NCFG; g++) printf(",k%d_L%d_mc,k%d_L%d_theory", CFG_K[g], CFG_L[g], CFG_K[g], CFG_L[g]);
    printf("\n");
    for (int a = 0; a < NANG; a++) {
        double p = 1.0 - a * 5.0 / 180.0;
        printf("%d,%.5f,%.5f", a * 5, (double)single[a] / (trials * MAXH), p);
        for (int g = 0; g < NCFG; g++)
            printf(",%.5f,%.5f", (double)hit[g][a] / trials,
                   1.0 - pow(1.0 - pow(p, CFG_K[g]), CFG_L[g]));
        printf("\n");
    }
    printf("# p-stable (Gaussian) LSH, w=1, trials=%ld\n", trials * PS_REP);
    printf("dist_over_w,p_mc,p_theory\n");
    for (int m = 0; m < NU; m++)
        printf("%.1f,%.5f,%.5f\n", 0.1 * m, (double)ps_hit[m] / (trials * PS_REP),
               pstable_p(0.1 * m, 1.0));
    return 0;
}

/* ---------- E2LSH index: g_t(v) = (h_1(v), ..., h_k(v)), h(v) = floor((a.v + b) / w) ---------- */

#define KMAX 32

typedef struct { uint64_t key; uint32_t id; } Entry;

typedef struct {
    int d, k, L;
    double w;
    float *at;      /* d x (k*L), transposed projection vectors */
    double *b;      /* k*L offsets in [0, w) */
    uint64_t **key; /* per table, sorted bucket keys */
    uint32_t **id;  /* per table, point ids in key order */
    long n;
} Index;

static uint64_t key_of(const int *h, int k)
{
    uint64_t x = 0x243f6a8885a308d3ULL;
    for (int i = 0; i < k; i++) {
        x ^= (uint64_t)(uint32_t)h[i] + 0x9e3779b97f4a7c15ULL + (x << 6) + (x >> 2);
        x = (x ^ (x >> 31)) * 0xbf58476d1ce4e5b9ULL;
    }
    return x ^ (x >> 29);
}

/* f[m] = (a_m . v + b_m) / w for all k*L hash functions */
static void project(const Index *ix, const float *v, double *f, float *acc)
{
    int M = ix->k * ix->L;
    memset(acc, 0, sizeof(float) * M);
    for (int t = 0; t < ix->d; t++) {
        const float *row = ix->at + (size_t)t * M;
        float vt = v[t];
        for (int m = 0; m < M; m++) acc[m] += vt * row[m];
    }
    for (int m = 0; m < M; m++) f[m] = (acc[m] + ix->b[m]) / ix->w;
}

static int cmp_entry(const void *a, const void *b)
{
    uint64_t x = ((const Entry *)a)->key, y = ((const Entry *)b)->key;
    if (x != y) return x < y ? -1 : 1;
    uint32_t i = ((const Entry *)a)->id, j = ((const Entry *)b)->id;
    return (i > j) - (i < j);
}

static void index_build(Index *ix, const float *data, long n, int d, int k, int L, double w)
{
    int M = k * L;
    ix->d = d; ix->k = k; ix->L = L; ix->w = w; ix->n = n;
    ix->at = xmalloc(sizeof(float) * (size_t)d * M);
    ix->b = xmalloc(sizeof(double) * M);
    for (int m = 0; m < M; m++) {
        for (int t = 0; t < d; t++) ix->at[(size_t)t * M + m] = (float)rng_gauss();
        ix->b[m] = w * rng_u01();
    }
    Entry **tab = xmalloc(sizeof(Entry *) * L);
    for (int l = 0; l < L; l++) tab[l] = xmalloc(sizeof(Entry) * n);
    double *f = xmalloc(sizeof(double) * M);
    float *acc = xmalloc(sizeof(float) * M);
    int h[KMAX];
    for (long i = 0; i < n; i++) {
        project(ix, data + (size_t)i * d, f, acc);
        for (int l = 0; l < L; l++) {
            for (int j = 0; j < k; j++) h[j] = (int)floor(f[l * k + j]);
            tab[l][i].key = key_of(h, k);
            tab[l][i].id = (uint32_t)i;
        }
    }
    ix->key = xmalloc(sizeof(uint64_t *) * L);
    ix->id = xmalloc(sizeof(uint32_t *) * L);
    for (int l = 0; l < L; l++) {
        qsort(tab[l], n, sizeof(Entry), cmp_entry);
        ix->key[l] = xmalloc(sizeof(uint64_t) * n);
        ix->id[l] = xmalloc(sizeof(uint32_t) * n);
        for (long i = 0; i < n; i++) { ix->key[l][i] = tab[l][i].key; ix->id[l][i] = tab[l][i].id; }
        free(tab[l]);
    }
    free(tab); free(f); free(acc);
}

/* bucket [lo, hi) of table l holding key */
static void bucket_range(const Index *ix, int l, uint64_t key, long *lo, long *hi)
{
    const uint64_t *a = ix->key[l];
    long a0 = 0, a1 = ix->n;
    while (a0 < a1) { long mid = (a0 + a1) / 2; if (a[mid] < key) a0 = mid + 1; else a1 = mid; }
    long b1 = a0;
    while (b1 < ix->n && a[b1] == key) b1++;
    *lo = a0; *hi = b1;
}

/* ---------- query-directed probing, Lv et al. VLDB 2007, Sec. 4.4 ----------
 * For hash i of table l, x_i(-1) = f_i - floor(f_i) is the distance (in units of w)
 * to the left boundary and x_i(+1) = 1 - x_i(-1). The 2k single-coordinate moves of a
 * table are sorted by z = x^2; a perturbation set A (positions in that order) has score
 * sum_{j in A} z_j. shift/expand from {0} enumerate every set exactly once in
 * non-decreasing score; one heap serves all L tables. */

typedef struct { double s; int tab, n; unsigned char j[KMAX]; } PSet;

typedef struct {
    int k, L;
    int *coord, *delta; /* [L][2k] */
    double *z;          /* [L][2k], ascending within a table */
    PSet *heap;
    int hn, hcap;
} Prober;

static int pset_less(const PSet *a, const PSet *b) { return a->s < b->s; }

static void heap_push(Prober *p, const PSet *x)
{
    if (p->hn == p->hcap) {
        p->hcap = p->hcap ? 2 * p->hcap : 64;
        PSet *nh = realloc(p->heap, sizeof(PSet) * p->hcap);
        if (!nh) { fprintf(stderr, "out of memory\n"); exit(1); }
        p->heap = nh;
    }
    int i = p->hn++;
    while (i > 0 && pset_less(x, &p->heap[(i - 1) / 2])) { p->heap[i] = p->heap[(i - 1) / 2]; i = (i - 1) / 2; }
    p->heap[i] = *x;
}

static PSet heap_pop(Prober *p)
{
    PSet top = p->heap[0], last = p->heap[--p->hn];
    int i = 0;
    for (;;) {
        int c = 2 * i + 1;
        if (c >= p->hn) break;
        if (c + 1 < p->hn && pset_less(&p->heap[c + 1], &p->heap[c])) c++;
        if (!pset_less(&p->heap[c], &last)) break;
        p->heap[i] = p->heap[c]; i = c;
    }
    if (p->hn > 0) p->heap[i] = last;
    return top;
}

typedef struct { double z; int coord, delta; } Move;

static int cmp_move(const void *a, const void *b)
{
    const Move *x = a, *y = b;
    if (x->z != y->z) return x->z < y->z ? -1 : 1;
    if (x->coord != y->coord) return x->coord - y->coord;
    return x->delta - y->delta;
}

/* frac[l*k + i] = f - floor(f) for the query */
static void prober_reset(Prober *p, const double *frac)
{
    Move mv[2 * KMAX];
    int k = p->k;
    p->hn = 0;
    for (int l = 0; l < p->L; l++) {
        for (int i = 0; i < k; i++) {
            double x = frac[l * k + i];
            mv[2 * i] = (Move){x * x, i, -1};
            mv[2 * i + 1] = (Move){(1 - x) * (1 - x), i, +1};
        }
        qsort(mv, 2 * k, sizeof(Move), cmp_move);
        for (int j = 0; j < 2 * k; j++) {
            p->coord[l * 2 * k + j] = mv[j].coord;
            p->delta[l * 2 * k + j] = mv[j].delta;
            p->z[l * 2 * k + j] = mv[j].z;
        }
        PSet s = {p->z[l * 2 * k], l, 1, {0}};
        heap_push(p, &s);
    }
}

static int pset_valid(const Prober *p, const PSet *a)
{
    unsigned seen = 0;
    for (int t = 0; t < a->n; t++) {
        unsigned bit = 1u << p->coord[a->tab * 2 * p->k + a->j[t]];
        if (seen & bit) return 0;
        seen |= bit;
    }
    return 1;
}

/* next valid perturbation; writes delta[k] and returns its table, or -1 when exhausted */
static int prober_next(Prober *p, int *delta, double *score)
{
    int k2 = 2 * p->k;
    while (p->hn > 0) {
        PSet a = heap_pop(p);
        const double *z = p->z + a.tab * k2;
        int mx = a.j[a.n - 1];
        if (mx + 1 < k2) {
            PSet s = a;                          /* shift: max(A) -> max(A) + 1 */
            s.j[s.n - 1] = (unsigned char)(mx + 1);
            s.s = a.s - z[mx] + z[mx + 1];
            heap_push(p, &s);
            if (a.n < p->k) {                    /* expand: add max(A) + 1 */
                PSet e = a;
                e.j[e.n++] = (unsigned char)(mx + 1);
                e.s = a.s + z[mx + 1];
                heap_push(p, &e);
            }
        }
        if (!pset_valid(p, &a)) continue;
        memset(delta, 0, sizeof(int) * p->k);
        for (int t = 0; t < a.n; t++)
            delta[p->coord[a.tab * k2 + a.j[t]]] = p->delta[a.tab * k2 + a.j[t]];
        *score = a.s;
        return a.tab;
    }
    return -1;
}

static void prober_init(Prober *p, int k, int L)
{
    memset(p, 0, sizeof(*p));
    p->k = k; p->L = L;
    p->coord = xmalloc(sizeof(int) * 2 * k * L);
    p->delta = xmalloc(sizeof(int) * 2 * k * L);
    p->z = xmalloc(sizeof(double) * 2 * k * L);
}

static void prober_free(Prober *p)
{
    free(p->coord); free(p->delta); free(p->z); free(p->heap);
}

static void index_free(Index *ix)
{
    for (int l = 0; l < ix->L; l++) { free(ix->key[l]); free(ix->id[l]); }
    free(ix->key); free(ix->id); free(ix->at); free(ix->b);
}

/* ---------- mode test: generator order == brute-force enumeration ---------- */

static int cmp_double(const void *a, const void *b)
{
    double x = *(const double *)a, y = *(const double *)b;
    return (x > y) - (x < y);
}

static int mode_test(void)
{
    enum { K = 5, L = 3, NV = 242 }; /* 3^5 - 1 non-zero perturbations per table */
    Prober p;
    prober_init(&p, K, L);
    rng_seed(1);
    int fails = 0;
    for (int trial = 0; trial < 2000; trial++) {
        double frac[K * L], brute[L * NV], got[L * NV];
        for (int i = 0; i < K * L; i++) frac[i] = rng_u01();
        int nb = 0;
        for (int l = 0; l < L; l++)
            for (int code = 1; code < NV + 1; code++) {
                double s = 0;
                for (int i = 0, c = code; i < K; i++, c /= 3) {
                    double x = frac[l * K + i];
                    if (c % 3 == 1) s += x * x;
                    if (c % 3 == 2) s += (1 - x) * (1 - x);
                }
                brute[nb++] = s;
            }
        qsort(brute, nb, sizeof(double), cmp_double);
        prober_reset(&p, frac);
        int delta[K], ng = 0, tab;
        unsigned char used[L][NV + 1];
        memset(used, 0, sizeof(used));
        double s;
        while ((tab = prober_next(&p, delta, &s)) >= 0) {
            int code = 0;
            for (int i = K - 1; i >= 0; i--) code = code * 3 + (delta[i] == 0 ? 0 : delta[i] < 0 ? 1 : 2);
            if (code == 0 || used[tab][code]++) { fails++; break; }
            if (ng < L * NV) got[ng] = s;
            ng++;
        }
        if (ng != nb) { fails++; continue; }
        for (int i = 0; i < nb; i++)
            if (fabs(got[i] - brute[i]) > 1e-12) { fails++; break; }
    }
    prober_free(&p);
    printf("probe-order test: 2000 random queries, k=%d, L=%d, %d perturbations each: %s\n",
           K, L, L * NV, fails ? "FAIL" : "ok");
    return fails != 0;
}

/* ---------- mode knn ---------- */

static void *load(const char *dir, const char *name, size_t bytes)
{
    char path[4096];
    snprintf(path, sizeof(path), "%s/%s", dir, name);
    FILE *f = fopen(path, "rb");
    if (!f) { perror(path); exit(1); }
    void *p = xmalloc(bytes);
    if (fread(p, 1, bytes, f) != bytes) { fprintf(stderr, "%s: short read\n", path); exit(1); }
    fclose(f);
    return p;
}

#define NCHK 12
static const int CHK[NCHK] = {0, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024};

/* recall of the true top-10 (and top-1) inside the candidate set; candidates are
 * re-ranked by exact distance, so this equals recall@10 of the final answer */
static void tally(const uint32_t *stamp, uint32_t st, const int32_t *g, double *r1, double *r10)
{
    int c = 0;
    for (int i = 0; i < 10; i++) c += stamp[g[i]] == st;
    *r1 += stamp[g[0]] == st;
    *r10 += c / 10.0;
}

static int mode_knn(const char *dir, int k, double w, int L, int T, int Q, uint64_t seed)
{
    const long n = 1000000;
    const int d = 128, NQ = 10000, G = 100;
    if (k < 1 || k > KMAX || L < 1 || Q < 1 || Q > NQ || T < 0 || T > CHK[NCHK - 1]) {
        fprintf(stderr, "bad parameters\n");
        return 2;
    }
    float *train = load(dir, "train.f32", sizeof(float) * n * d);
    float *test = load(dir, "test.f32", sizeof(float) * (size_t)NQ * d);
    int32_t *gt = load(dir, "gt.i32", sizeof(int32_t) * (size_t)NQ * G);
    rng_seed(seed);
    Index ix;
    index_build(&ix, train, n, d, k, L, w);

    long nb = 0, mx = 0;
    double sq = 0;
    for (long i = 0; i < n;) {
        long j = i;
        while (j < n && ix.key[0][j] == ix.key[0][i]) j++;
        nb++; if (j - i > mx) mx = j - i; sq += (double)(j - i) * (j - i);
        i = j;
    }
    printf("# knn n=%ld d=%d k=%d w=%g L=%d T=%d Q=%d seed=%llu\n", n, d, k, w, L, T, Q,
           (unsigned long long)seed);
    printf("# table0 buckets=%ld mean_size=%.2f max_size=%ld size_biased_mean=%.1f\n",
           nb, (double)n / nb, mx, sq / n);

    int M = k * L;
    uint32_t *stamp = calloc(n, sizeof(uint32_t));
    double *f = xmalloc(sizeof(double) * M), *frac = xmalloc(sizeof(double) * M);
    float *acc = xmalloc(sizeof(float) * M);
    int *base = xmalloc(sizeof(int) * M), hp[KMAX], delta[KMAX];
    double *sc = calloc(L, sizeof(double)), *s1 = calloc(L, sizeof(double)), *s10 = calloc(L, sizeof(double));
    double mc[NCHK] = {0}, m1[NCHK] = {0}, m10[NCHK] = {0};
    if (!stamp || !sc || !s1 || !s10) { fprintf(stderr, "out of memory\n"); return 1; }
    Prober p;
    prober_init(&p, k, L);

    for (int q = 0; q < Q; q++) {
        uint32_t st = (uint32_t)q + 1;
        const int32_t *g = gt + (size_t)q * G;
        long cand = 0, lo, hi;
        project(&ix, test + (size_t)q * d, f, acc);
        for (int m = 0; m < M; m++) { base[m] = (int)floor(f[m]); frac[m] = f[m] - base[m]; }
        for (int l = 0; l < L; l++) {            /* standard LSH: home bucket per table */
            bucket_range(&ix, l, key_of(base + l * k, k), &lo, &hi);
            for (long i = lo; i < hi; i++) {
                uint32_t id = ix.id[l][i];
                if (stamp[id] != st) { stamp[id] = st; cand++; }
            }
            sc[l] += cand;
            tally(stamp, st, g, &s1[l], &s10[l]);
        }
        int ci = 0, extra = 0;
        if (T > 0) prober_reset(&p, frac);
        for (;;) {                               /* multi-probe on the same L tables */
            while (ci < NCHK && CHK[ci] <= T && CHK[ci] == extra) {
                mc[ci] += cand;
                tally(stamp, st, g, &m1[ci], &m10[ci]);
                ci++;
            }
            if (extra >= T) break;
            double s;
            int tab = prober_next(&p, delta, &s);
            if (tab < 0) { extra = CHK[ci]; continue; } /* exhausted: repeat last state */
            for (int j = 0; j < k; j++) hp[j] = base[tab * k + j] + delta[j];
            bucket_range(&ix, tab, key_of(hp, k), &lo, &hi);
            for (long i = lo; i < hi; i++) {
                uint32_t id = ix.id[tab][i];
                if (stamp[id] != st) { stamp[id] = st; cand++; }
            }
            extra++;
        }
    }
    printf("kind,tables,extra_probes,buckets_probed,candidates,cand_frac,recall1,recall10\n");
    for (int l = 0; l < L; l++)
        printf("std,%d,0,%d,%.1f,%.5f,%.4f,%.4f\n", l + 1, l + 1, sc[l] / Q, sc[l] / Q / n,
               s1[l] / Q, s10[l] / Q);
    for (int c = 0; c < NCHK && CHK[c] <= T; c++)
        printf("mp,%d,%d,%d,%.1f,%.5f,%.4f,%.4f\n", L, CHK[c], L + CHK[c], mc[c] / Q,
               mc[c] / Q / n, m1[c] / Q, m10[c] / Q);
    prober_free(&p);
    index_free(&ix);
    free(stamp); free(f); free(frac); free(acc); free(base); free(sc); free(s1); free(s10);
    free(train); free(test); free(gt);
    return 0;
}

int main(int argc, char **argv)
{
    if (argc >= 2 && !strcmp(argv[1], "test")) return mode_test();
    if (argc == 4 && !strcmp(argv[1], "curve"))
        return mode_curve(atol(argv[2]), strtoull(argv[3], NULL, 10));
    if (argc == 9 && !strcmp(argv[1], "knn"))
        return mode_knn(argv[2], atoi(argv[3]), atof(argv[4]), atoi(argv[5]), atoi(argv[6]),
                        atoi(argv[7]), strtoull(argv[8], NULL, 10));
    fprintf(stderr, "usage: lsh test | lsh curve TRIALS SEED | lsh knn DIR K W L T Q SEED\n");
    return 2;
}