#include <stdlib.h>
#include <ctype.h>
#include <math.h>
#include <assert.h>
#include <vlcutils/error.h>
#include <vlcutils/heap.h>
#include <vlcutils/set.h>
#include <vlcutils/bplus.h>
#include <vlcutils/bitmap.h>
#include "definitions.h"
#include "freqdist.h"
#include "index.h"
#include "matchtable.h"

struct matchtable {
     const struct index *index;
     unsigned char *data;
     int index_size, query_size, page_size, window_size;
     int marked; /* Number of entries marked */
};

static struct matchtable *make_matchtable(const struct index *index, int qlen)
{
     struct matchtable *mt;
     size_t data_size;
     int page_size;
     MBR last_box;

     assert(index);
     assert(qlen >= 0);
     page_size = index->window_size;
     mt = malloc(sizeof(struct matchtable));
     if (!mt) fatal_perror("malloc");
     mt->index = index;
     last_box = index->elements[index->num_boxes - 1];
     mt->query_size = ceil(qlen / (page_size * 1.0));
     mt->index_size = ceil((last_box.start + last_box.capacity + mt->window_size)
			   / (page_size * 1.0));
     fprintf(stderr, "Match table size: %d × %d (index × query) = %d entries ", 
	    mt->index_size, mt->query_size, mt->index_size * mt->query_size);
     mt->index_size = ceil(mt->index_size / 8.0);
     mt->window_size = index->window_size;
     mt->page_size = page_size;
     data_size = mt->index_size * mt->query_size;
     fprintf(stderr, "(%lu bytes)\n", (long unsigned)data_size);
     mt->data = calloc(data_size, 1);
     if (!mt->data) fatal_perror("make_matchtable(): calloc(%d)", data_size);
     mt->marked = 0;
     return mt;
}

void free_matchtable(struct matchtable *mt)
{
     assert(mt);
     assert(mt->data);
     free(mt->data);
     mt->data = NULL;
     free(mt);
}

static void matchtable_mark(struct matchtable *mt, int box_number, 
			    int query_page)
{
     int i, i_max, j;
     MBR box;
     int value, shift;

     assert(mt);
     box = mt->index->elements[box_number];
     i = box.start / mt->page_size;
     i_max = (box.start + box.capacity + mt->window_size - 2) / mt->page_size;
     for (;;) {
	  j = query_page;
	  assert(0 <= j && j < mt->query_size);
	  assert(0 <= i && i < mt->index_size * 8);
	  value = mt->data[j * mt->index_size + i / 8];
	  shift = 7 - i % 8;
	  if (!(value & 1 << shift))
	       mt->marked++;
	  value |= 1 << shift;
	  mt->data[j * mt->index_size + i / 8] = value;
	  i++;
	  if (i > i_max)
	       break;
     }
}

static struct matchtable *make_matchtable_boxes(const struct index *index1,
                                                const struct index *index2)
{
     struct matchtable *mt;
     size_t data_size;
     int page_size;
     MBR last_box;

     assert(index1 && index2);
     assert(index1->window_size == index2->window_size);
     page_size = index1->window_size;
     mt = malloc(sizeof(struct matchtable));
     if (!mt)
	  fatal_perror("malloc");
     mt->window_size = index1->window_size;
     mt->page_size = page_size;
     mt->index = index1;
     last_box = index1->elements[index1->num_boxes - 1];
     mt->index_size = ceil((last_box.start + last_box.capacity +
                            mt->window_size) / (page_size * 1.0));
     last_box = index2->elements[index2->num_boxes - 1];
     mt->query_size = ceil((last_box.start + last_box.capacity +
                            mt->window_size) / (page_size * 1.0));
     fprintf(stderr, "Match table size: %d × %d (index × query) = %d entries", 
	    mt->index_size, mt->query_size, mt->index_size * mt->query_size);
     mt->index_size = ceil(mt->index_size / 8.0);
     data_size = mt->index_size * mt->query_size;
     fprintf(stderr, " (%lu bytes)\n", (long unsigned)data_size);
     mt->data = calloc(data_size, 1);
     if (!mt->data)
	  fatal_perror("make_matchtable(): calloc(%d)", data_size);
     mt->marked = 0;
     return mt;
}

static void matchtable_mark_boxes(struct matchtable *mt, MBR box1, MBR box2)
{
     int i, i_max, j, j_max;
     int value, shift;

     assert(mt);
     i_max = (box1.start + box1.capacity + mt->window_size - 2) 
          / mt->page_size;
     j_max = (box2.start + box2.capacity + mt->window_size - 2) 
          / mt->page_size;
     for (j = box2.start / mt->page_size; j <= j_max; j++) {
          for (i = box1.start / mt->page_size; i <= i_max; i++) {
               assert(0 <= j && j < mt->query_size);
               assert(0 <= i && i < mt->index_size * 8);
               value = mt->data[j * mt->index_size + i / 8];
               shift = 7 - i % 8;
               if (!(value & 1 << shift))
                    mt->marked++;
               value |= 1 << shift;
               mt->data[j * mt->index_size + i / 8] = value;
          }
     }
}

void fprint_matchtable(FILE *stream, struct matchtable *mt)
{
  int i, j;
  int value;

  /* -1 is just junk to fill an unused field. */
  fprintf(stream, "%d %d\n", -1, mt->page_size);
  for (j = 0; j < mt->query_size; j++) {
       for (i = 0; i < mt->index_size; i++) {
	    value = mt->data[j * mt->index_size + i];
	    if (value & (1 << 7)) fprintf(stream, "%d %d\n", i * 8 + 0, j);
	    if (value & (1 << 6)) fprintf(stream, "%d %d\n", i * 8 + 1, j);
	    if (value & (1 << 5)) fprintf(stream, "%d %d\n", i * 8 + 2, j);
	    if (value & (1 << 4)) fprintf(stream, "%d %d\n", i * 8 + 3, j);
	    if (value & (1 << 3)) fprintf(stream, "%d %d\n", i * 8 + 4, j);
	    if (value & (1 << 2)) fprintf(stream, "%d %d\n", i * 8 + 5, j);
	    if (value & (1 << 1)) fprintf(stream, "%d %d\n", i * 8 + 6, j);
	    if (value & (1 << 0)) fprintf(stream, "%d %d\n", i * 8 + 7, j);
       }
  }
}

void fwrite_matchtable(FILE *stream, struct matchtable *mt)
{
     if (fwrite(mt->data, mt->query_size * mt->index_size, 1, stream) != 1)
          fatal_perror("fwrite");
}

/*************************************************************************/

/* Events for the plane sweep. */
typedef enum { e_none, e_start, e_end, e_point } event_type;

struct event {
     event_type type;
     int point[ALPH_LEN];
     int box_number;       /* Also used for page number (for point events) */
};

struct event_queue {
     int size;
     int fill_pointer;
     struct event *events;
};

void init_event_queue(struct event_queue *eq, int size)
{
     eq->size = size;
     eq->fill_pointer = 0;
     eq->events = malloc(size * sizeof(struct event));
     if (!eq->events) fatal_perror("init_event_queue(): malloc");
 }

void event_queue_realloc(struct event_queue *eq)
{
     assert(eq && eq->events);
     eq->size = eq->fill_pointer;
     eq->events = realloc(eq->events, eq->size * sizeof(struct event));
     if (!eq->events && eq->size)
          fatal_perror("event_queue_realloc(): realloc");
}

static void event_queue_add(struct event_queue *eq, event_type type, 
			    int point[ALPH_LEN], int box_number, int grow)
{
     int i;

     assert(eq->size > eq->fill_pointer);
     eq->events[eq->fill_pointer].type = type;
     for (i = 0; i < ALPH_LEN; i++)
	  eq->events[eq->fill_pointer].point[i] = point[i];
     eq->events[eq->fill_pointer].point[0] += grow;
     eq->events[eq->fill_pointer].box_number = box_number;
     eq->fill_pointer++;
}

/* Read query from disk and add points to the event queue.  eq is the
 * event queue, and qlen is the number of A, C, G, T in the file. */
void event_queue_read_query(struct event_queue *eq, FILE *input, 
			    int qlen, int page_size, int window_size)
{
     char c, *qs;
     int i, j, k, m, wPerPage, delta, qnum, start1;
     int q[ALPH_LEN];

     assert(input);
     /* Read the complete query into memory. */
     qs = (char*)malloc(sizeof(char) * qlen);
     if (!qs)
	  fatal_perror("malloc");
     j = 0;
     while (fscanf(input, "%c", &qs[j]) > 0) 
	  j++;
     qlen = j;

     delta = window_size;
     qnum = floor((qlen - window_size + delta - 1) / delta);
     wPerPage = floor(page_size / window_size);
     if (wPerPage == 0)
	  wPerPage = 1;

     start1 = k = 0; 
     while (k < qnum) {
	  /* Consume a page size portion of the query. */
	  for (i = 0; i < wPerPage; i++){
	       /* Initialize event */
	       for(j = 0; j < ALPH_LEN; j++)
		    q[j] = 0;
	       /* map k'th subquery into integer space */
	       j = start1;
	       for (m = 0; m < window_size; m++) {
		    c = qs[j];
		    j++;
		    increment(c, q);
	       }
	       event_queue_add(eq, e_point, q, k, 0);
	       /* Rearrange the boundaries for the next subquery. */
	       start1 = start1 + delta;
	       if (k == qnum - 1) {
		    if (start1 > qlen)
			 start1 = qlen - window_size;
	       }
	       k++;
	       if (k == qnum)
		    break;
	  }
     }
     free(qs);
}

static int compare_events_by_x(const void *datum1, const void *datum2)
{
     int x1, x2;

     x1 = ((struct event *)datum1)->point[0];
     x2 = ((struct event *)datum2)->point[0];
     return x1 < x2 ? -1 : (x1 == x2 ? 0 : 1);
}

#if 0
static void print_point(int q[ALPH_LEN], FILE *stream)
{
     int i;
     assert(ALPH_LEN > 0);
     fprintf(stream, "[%d", q[0]);
     for (i = 1; i < ALPH_LEN; i++)
	  fprintf(stream, ", %d", q[i]);
     fprintf(stream, "]");
}
#endif

/* Returns the minimum distance between a point and the intersection
 * of an MBR with the SI-plane with respect to the frequency of each
 * letter.
 */
static __inline__ int infdist4(int q[], MBR box, int maxdist)
{
     int posDist, negDist, sum;
     int q0, q1, q2, q3, bl0, bl1, bl2, bl3, bh0, bh1, bh2, bh3, qw;

     sum = posDist = negDist = 0;

     q1 = q[1]; bl1 = box.lower[1]; bh1 = box.higher[1];
     if (bl1 > q1 + maxdist || bh1 < q1 - maxdist)
       return maxdist + 1;
     q2 = q[2]; bl2 = box.lower[2]; bh2 = box.higher[2];
     if (bl2 > q2 + maxdist || bh2 < q2 - maxdist)
       return maxdist + 1;
     q3 = q[3]; bl3 = box.lower[3]; bh3 = box.higher[3];
     if (bl3 > q3 + maxdist || bh3 < q3 - maxdist)
       return maxdist + 1;

     q0 = q[0]; bl0 = box.lower[0]; bh0 = box.higher[0]; 
     if (q0 < bl0) {
	  posDist += bl0 - q0;
	  sum += bl0;
     }
     else if (q0 > bh0) {
	  negDist -= bh0 - q0;
	  sum += bh0;
     }
     else
	  sum += q0;

     if (q1 < bl1) {
	  posDist += bl1 - q1;
	  sum += bl1;
     }
     else if (q1 > bh1) {
	  negDist -= bh1 - q1;
	  sum += bh1;
     }
     else
	  sum += q1;

     if (q2 < bl2) {
	  posDist += bl2 - q2;
	  sum += bl2;
     }
     else if (q2 > bh2) {
	  negDist -= bh2 - q2;
	  sum += bh2;
     }
     else
	  sum += q2;

     if (q3 < bl3) {
	  posDist += bl3 - q3;
	  sum += bl3;
     }
     else if (q3 > bh3) {
	  negDist -= bh3 - q3;
	  sum += bh3;
     }
     else
	  sum += q3;

     qw = q0 + q1 + q2 + q3;
     if (sum > qw) {
	  posDist = 4 * posDist;
	  negDist = 4 * negDist + 4 + sum - qw;
     }
     else {
	  posDist = 4 * posDist + 4 + qw - sum;
	  negDist = 4 * negDist;
     }

     return (posDist > negDist) ? posDist : negDist;
}

struct matchtable *
compute_matchtable(FILE *input, const struct index *index, double err)
{
     int qlen, qnum;
     struct event_queue eq;
     int i, j;
     int delta;
     struct set *status;
     int box_number, query_page_number;
     int maxdist, grow, distance;
     struct matchtable *mt;

     assert(index);
     fprintf(stderr, "window size = %d, box_capacity = %d, num_boxes = %d\n", 
	    index->window_size, index->box_capacity, index->num_boxes);
     maxdist = err * index->window_size;
     /* Add 1 so the box includes points exactly maxdist away. */
     grow = ceil(maxdist / 4.0) + 1;

     qlen = file_length(input);
     delta = index->window_size;
     qnum = floor((qlen - index->window_size + delta - 1) / delta);
     init_event_queue(&eq, qnum + index->num_boxes * 2);
     /* Derive points from the query file and add events to the event queue. */
     event_queue_read_query(&eq, input, qlen, index->window_size, 
			    index->window_size);
     /* Add events for the boxes in the index */
     for (i = 0; i < index->num_boxes; i++) {
	  /* FIXME */
	  event_queue_add(&eq, e_start, index->elements[i].lower, i, -grow);
	  event_queue_add(&eq, e_end, index->elements[i].higher, i, +grow);
     }
     mt = make_matchtable(index, qlen);
     status = make_set(index->num_boxes);

     qsort((void *)eq.events, eq.fill_pointer, sizeof(struct event),
	   compare_events_by_x);
     /* For each query point... */
     for (i = 0; i < eq.fill_pointer; i++) {
	  switch (eq.events[i].type) {
	  case e_start:
	       set_add(status, eq.events[i].box_number);
	       break;
	  case e_end:
	       set_remove(status, eq.events[i].box_number);
	       break;
	  case e_point:
	       for (j = 0; j < status->fill_pointer; j++) {
		    box_number = status->items[j];
		    distance = infdist4(eq.events[i].point, 
				       index->elements[box_number], maxdist);
		    if (distance <= maxdist) {
			 query_page_number = eq.events[i].box_number;
			 matchtable_mark(mt, box_number, query_page_number);
		    }
	       }
	       break;
	  default:
	       error("Unknown event type: %d", eq.events[i].type);
	       abort();
	  }
     }
     assert(status->fill_pointer == 0);
     fprintf(stderr, "Done.  Marked entries: %d\n", mt->marked);
     return mt;
}

struct matchtable *
compute_matchtable_bplus(FILE *input, const struct index *index, double err)
{
     int qlen, qnum;
     struct event_queue eq;
     int i;
     int delta;
     struct bplus end_tree;
     int query_page_number;
     int maxdist, grow, distance;
     struct matchtable *mt;
     MBR *box;
     int py;
     long matches = 0;

     assert(index);
     fprintf(stderr, "window size = %d, box_capacity = %d, num_boxes = %d\n", 
	    index->window_size, index->box_capacity, index->num_boxes);
     maxdist = err * index->window_size;
     /* Add 1 so the box includes points exactly maxdist away. */
     grow = ceil(maxdist / 4.0) + 1;

     qlen = file_length(input);
     delta = index->window_size;
     qnum = floor((qlen - index->window_size + delta - 1) / delta);
     init_event_queue(&eq, qnum + index->num_boxes * 2);
     event_queue_read_query(&eq, input, qlen, index->window_size, 
			    index->window_size);
     for (i = 0; i < index->num_boxes; i++) {
	  event_queue_add(&eq, e_start, index->elements[i].lower, i, -grow);
	  event_queue_add(&eq, e_end, index->elements[i].higher, i, +grow);
     }
     mt = make_matchtable(index, qlen);
     init_bplus(&end_tree, index->num_boxes, 
                (char *)&box->higher[1] - (char *)box);
     qsort((void *)eq.events, eq.fill_pointer, sizeof(struct event),
	   compare_events_by_x);
     /* For each query point... */
     for (i = 0; i < eq.fill_pointer; i++) {
	  switch (eq.events[i].type) {
	  case e_start:
               box = &index->elements[eq.events[i].box_number];
               bplus_insert(&end_tree, box);
	       break;
	  case e_end:
               box = &index->elements[eq.events[i].box_number];
               bplus_delete(&end_tree, box);
	       break;
	  case e_point:
               py = eq.events[i].point[1];
               box = bplus_search(&end_tree, py - maxdist);
               while (box) {
		    distance = infdist4(eq.events[i].point, *box, maxdist);
		    if (distance <= maxdist) {
                         matches++;
			 query_page_number = eq.events[i].box_number;
			 matchtable_mark(mt, box->box_number, 
                                         query_page_number);
		    }
                    box = bplus_search(NULL, 0);
               }
	       break;
	  default:
	       error("Unknown event type: %d", eq.events[i].type);
	       abort();
	  }
     }
     assert(bplus_size(&end_tree) == 0);
     fprintf(stderr, "Done.  Marked entries: %d\n", mt->marked);
     fprintf(stderr, "Matches in B+-tree: %ld\n", matches);
     return mt;
}

#if 0
struct matchtable *
compute_matchtable_nlj(FILE *input, const struct index *index, double err)
{
     int qlen, qnum;
     struct event_queue eq;
     int i, j;
     int delta;
     int box_number, query_page_number;
     int maxdist, distance;
     struct matchtable *mt;

     assert(index);

     fprintf(stderr, "window size = %d, box_capacity = %d, num_boxes = %d\n", 
	    index->window_size, index->box_capacity, index->num_boxes);
     maxdist = err * index->window_size;

     qlen = file_length(input);
     delta = index->window_size;
     qnum = floor((qlen - index->window_size + delta - 1) / delta);
     init_event_queue(&eq, qnum);
     /* Derive points from the query file and add events to the event queue. */
     event_queue_read_query(&eq, input, qlen, index->window_size, 
			    index->window_size);
     mt = make_matchtable(index, qlen);
     for (i = 0; i < qnum; i++) {
          for (j = 0; j < index->num_boxes; j++) {
               box_number = j;
               distance = infdist4(eq.events[i].point, 
                                   index->elements[box_number], maxdist);
               if (distance <= maxdist) {
                    query_page_number = eq.events[i].box_number;
                    matchtable_mark(mt, box_number, query_page_number);
               }
          }
     }
     fprintf(stderr, "Done.  Marked entries: %d\n", mt->marked);
     return mt;
}
#endif

static int event_lt_x(struct event *e1, struct event *e2)
{
     return e1->point[0] < e2->point[0];
}

struct matchtable * 
compute_matchtable_boxes(const struct index *index1, 
                         const struct index *index2, double err)
{
     int window_size;
     struct event *events, *event;
     struct heap *eq;
     int i, j, k;
     struct set *status1, *status2;
     int box_number, event_number;
     int maxdist, grow, distance;
     struct matchtable *mt;
     MBR box, box1, box2;

     assert(index1 && index2);

     fprintf(stderr, "window size = %d, box_capacity = %d, num_boxes = %d\n", 
	    index1->window_size, index1->box_capacity, index1->num_boxes);
     fprintf(stderr, "window size = %d, box_capacity = %d, num_boxes = %d\n", 
	    index2->window_size, index2->box_capacity, index2->num_boxes);
     if (index1->window_size != index2->window_size)
          fatal_error("Indexes do not have the same window size: %d, %d",
                      index1->window_size, index2->window_size);
     window_size = index1->window_size;
     maxdist = err * window_size;
     /* Add 1 so the box includes points exactly maxdist away. */
     grow = ceil(maxdist / 4.0) + 1;

     events = malloc((index1->num_boxes + index2->num_boxes) *
                     sizeof(struct event));
     if (!events)
          fatal_perror("malloc");
     j = 0;
     eq = make_heap(index1->num_boxes + index2->num_boxes, event_lt_x);
     /* Add start events event for the boxes in the indexes */
     for (i = 0; i < index1->num_boxes; i++) {
          events[j].type = e_start;
          for (k = 0; k < ALPH_LEN; k++)
               events[j].point[k] = index1->elements[i].lower[k] - grow;
          events[j].box_number = j;
          heap_insert(&events[j], eq);
          j++;
     }
     for (i = 0; i < index2->num_boxes; i++) {
          events[j].type = e_start;
          for (k = 0; k < ALPH_LEN; k++)
               events[j].point[k] = index2->elements[i].lower[k] - grow;
          events[j].box_number = j;
          heap_insert(&events[j], eq);
          j++;
     }
     mt = make_matchtable_boxes(index1, index2);
     status1 = make_set(index1->num_boxes);
     status2 = make_set(index2->num_boxes);
     /* For each query point... */
     while (heap_size(eq) > 0) {
          event = heap_pop(eq);
	  switch (event->type) {
	  case e_start:
               box_number = event_number = event->box_number;
               if (box_number < index1->num_boxes) {
                    /* Adding a new box from the first index. */
                    set_add(status1, box_number);
                    box = box1 = index1->elements[box_number];
                    for (j = 0; j < status2->fill_pointer; j++) {
                         box2 = index2->elements[status2->items[j]];
                         distance = box_box_distance(box1, box2, window_size);
                         /* if (distance <= maxdist) */
                              matchtable_mark_boxes(mt, box1, box2);
                    }
               } else {
                    /* Adding a new box from the second index. */
                    box_number -= index1->num_boxes;
                    set_add(status2, box_number);
                    box = box2 = index2->elements[box_number];
                    for (j = 0; j < status1->fill_pointer; j++) {
                         box1 = index1->elements[status1->items[j]];
                         distance = box_box_distance(box1, box2, window_size);
                         /* if (distance <= maxdist) */
                              matchtable_mark_boxes(mt, box1, box2);
                    }
               }
               events[event_number].type = e_end;
               for (k = 0; k < ALPH_LEN; k++)
                    events[event_number].point[k] = box.higher[k] + grow;
               heap_insert(&events[event_number], eq);
	       break;
	  case e_end:
               if (event->box_number < index1->num_boxes)
                    set_remove(status1, event->box_number);
               else
                    set_remove(status2, event->box_number - index1->num_boxes);
	       break;
	  default:
	       error("Unknown event type: %d", event->type);
	       abort();
	  }
     }
     assert(status1->fill_pointer == 0);
     assert(status2->fill_pointer == 0);
     free_heap(eq);
     free(events);
     fprintf(stderr, "Done.  Marked entries: %d\n", mt->marked);
     return mt;
}

struct matchtable * 
compute_matchtable_nlj_boxes(const struct index *index1, 
                             const struct index *index2, double err)
{
     int window_size;
     int i, j, k;
     int maxdist, distance;
     long matches = 0;
     struct matchtable *mt;
     int q[ALPH_LEN];
     MBR box1, box2;

     assert(index1 && index2);

     fprintf(stderr, "window size = %d, box_capacity = %d, num_boxes = %d\n", 
	    index1->window_size, index1->box_capacity, index1->num_boxes);
     fprintf(stderr, "window size = %d, box_capacity = %d, num_boxes = %d\n", 
	    index2->window_size, index2->box_capacity, index2->num_boxes);
     if (index1->window_size != index2->window_size)
          fatal_error("Indexes do not have the same window size: %d, %d",
                      index1->window_size, index2->window_size);
     window_size = index1->window_size;
     maxdist = err * window_size;
     mt = make_matchtable_boxes(index1, index2);

     for (i = 0; i < index1->num_boxes; i++)
          for (j = 0; j < index2->num_boxes; j++) {
               box1 = index1->elements[i];
	       maxdist = 0;
	       for (k = 0; k < ALPH_LEN; k++) {
		    q[k] = (box1.lower[k] + box1.higher[k])/2;
		    maxdist += (box1.higher[k] - box1.lower[k])/2;
	       }
	       maxdist = pow(maxdist, 1/ALPH_LEN);
	       maxdist += err * window_size;
               box2 = index2->elements[j];
               distance = infdist(q, box2);
               if (distance <= maxdist) {
		    matches++;
                    matchtable_mark_boxes(mt, box1, box2);
	       }
          }
     fprintf(stderr, "Box--box matches: %ld\n", matches);
     fprintf(stderr, "Done.  Marked entries: %d\n", mt->marked);
     return mt;
}

struct matchtable *
compute_matchtable_nlj(FILE *input, const struct index *index, double err)
{
     int qlen, qnum;
     struct event_queue eq;
     MBR box;
     int i, j, delta, query_page_number, maxdist, distance;
     struct matchtable *mt;
     long comparisons = 0;

     assert(index);

     fprintf(stderr, "window size = %d, box_capacity = %d, num_boxes = %d\n", 
	    index->window_size, index->box_capacity, index->num_boxes);
     maxdist = err * index->window_size;

     qlen = file_length(input);
     delta = index->window_size;
     qnum = floor((qlen - index->window_size + delta - 1) / delta);
     init_event_queue(&eq, qnum);
     event_queue_read_query(&eq, input, qlen, index->window_size, 
			    index->window_size);
     mt = make_matchtable(index, qlen);
     for (i = 0; i < index->num_boxes; i++) {
          box = index->elements[i];
          for (j = 0; j < qnum; j++) {
               comparisons++;
               distance = infdist4(eq.events[j].point, box, maxdist);
               if (distance <= maxdist) {
                    query_page_number = eq.events[j].box_number;
                    matchtable_mark(mt, box.box_number, query_page_number);
               }
          }
     }
     fprintf(stderr, "Done.  Marked entries: %d\n", mt->marked);
     fprintf(stderr, "Comparisons: %ld\n", comparisons);
     return mt;
}

/***********************************************************************
 *** PARTITIONING                                                    ***
 ***********************************************************************/

/* Partition a set of boxes along one dimension (d).*/
static void partition1(struct event_queue *eq, const struct index *index, 
                       bitmap **boxmaps, int partitions, int s, int d)
{
     int *histogram, *map, i, j, p, pp, total_depth = 0, depth, target_depth;
     int qnum, factor;
     MBR box;
     bitmap *tmpboxmap;

     qnum = eq[s].fill_pointer;
     eq[s].fill_pointer = 0;
     tmpboxmap = boxmaps[(int)pow(partitions, 3)];
     bitmap_copy(tmpboxmap, boxmaps[s], index->num_boxes);
     assert(bitmap_count(tmpboxmap, index->num_boxes) == 
            bitmap_count(boxmaps[s], index->num_boxes));
     bitmap_clear(boxmaps[s], index->num_boxes);
     factor = pow(partitions, d);
     /* Build histogram. */
     histogram = calloc(index->window_size, sizeof(int));
     if (!histogram) fatal_perror("calloc");
     for (i = 0; i < qnum; i++) {
          assert(eq[s].events);
          histogram[eq[s].events[i].point[d]]++, total_depth++;
     }
     /* Build map */
     map = malloc(index->window_size * sizeof(int));
     if (!map) fatal_perror("malloc");
     target_depth = ceil(total_depth / (partitions * 1.0));
     for (p = 0, depth = 0, j = 0; j < index->window_size; j++) {
          map[j] = p;
          depth += histogram[j];
          if (depth >= target_depth && p + 1 < partitions)
               p++, depth = 0;
     }
     free(histogram);
     /* Initialize lists for points */
     for (p = 0; p < partitions; p++) {
          pp = s + p * factor;
          if (!eq[pp].events)
               init_event_queue(&eq[pp], qnum);
     }
     /* Assign points to partitions */
     for (i = 0; i < qnum; i++) {
          assert(eq[s].events);
          assert(eq[s].events[i].point[d] >= 0 && 
                 eq[s].events[i].point[d] < index->window_size);
          pp = s + map[eq[s].events[i].point[d]] * factor;
          eq[pp].events[eq[pp].fill_pointer++] = eq[s].events[i];
     }
     /* Free up unused memory */
     for (p = 0; p < partitions; p++) {
          pp = s + p * factor;
          event_queue_realloc(&eq[pp]);
     }
     /* Assign boxes to partitions */
     for (i = 0; i < index->num_boxes; i++) {
          if (bitmap_get(tmpboxmap, i)) {
               box = index->elements[i];
               assert(box.lower[d] >= 0 && box.lower[d] < index->window_size &&
                      box.higher[d] >= 0 && box.higher[d] < index->window_size);
               pp = s + map[box.lower[d]] * factor;
               bitmap_set(boxmaps[pp], i);
               if (map[box.lower[d]] != map[box.higher[d]]) {
                    for (p = map[box.lower[d]] + 1; p < map[box.higher[d]]; p++)
                    {
                         pp = s + p * factor;
                         bitmap_set(boxmaps[pp], i);
                    }
                    pp = s + map[box.higher[d]] * factor;
                    bitmap_set(boxmaps[pp], i);
               }
          }
     }
     free(map);
}

static void partition(struct event_queue *eq, int partitions, 
                      const struct index *index, bitmap **boxmaps)
{
     int p0, p1, pp, sumboxes, sumpoints;

     partition1(eq, index, boxmaps, partitions, 0, 0);
     for (p0 = 0; p0 < partitions; p0++) {
          partition1(eq, index, boxmaps, partitions, p0, 1);
          for (p1 = 0; p1 < partitions; p1++) {
               partition1(eq, index, boxmaps, partitions,
                          p0 + p1 * pow(partitions, 1), 2);
          }
     }
     sumboxes = sumpoints = 0;
     for (pp = 0; pp < pow(partitions, 3); pp++) {
          sumpoints += eq[pp].fill_pointer;
          sumboxes += bitmap_count(boxmaps[pp], index->num_boxes);
     }
     fprintf(stderr, "%d points, %d boxes\n", sumpoints, sumboxes);
}

struct matchtable *
compute_matchtable_partitions(FILE *input, const struct index *index, 
                              double err, int partitions)
{
     int qlen, qnum;
     struct event_queue *eq;
     bitmap **boxmaps;
     struct index *grownindex;
     MBR box;
     int i, j, pp, delta, query_page_number, maxdist, distance;
     struct matchtable *mt;
     long comparisons = 0, markings = 0;

     assert(index);

     fprintf(stderr, "window size = %d, box_capacity = %d, num_boxes = %d\n", 
	    index->window_size, index->box_capacity, index->num_boxes);
     maxdist = err * index->window_size;
     grownindex = index_dup(index);
     /* Add 1 so the box includes points exactly maxdist away. */
     index_grow(grownindex, ceil(maxdist / 4.0) + 1);

     qlen = file_length(input);
     delta = index->window_size;
     qnum = floor((qlen - index->window_size + delta - 1) / delta);
     eq = calloc(pow(partitions, 3), sizeof(struct event_queue));
     if (!eq) fatal_perror("calloc");
     init_event_queue(&eq[0], qnum);
     event_queue_read_query(&eq[0], input, qlen, index->window_size, 
			    index->window_size);
     boxmaps = calloc(pow(partitions, 3) + 1, sizeof(bitmap *));
     if (!boxmaps) fatal_perror("calloc");
     for (i = 0; i < pow(partitions, 3) + 1; i++)
          boxmaps[i] = make_bitmap(index->num_boxes);
     bitmap_setall(boxmaps[0], index->num_boxes);
     fprintf(stderr, "bitmap_count(boxmaps[0], num_boxes) = %d\n",
             bitmap_count(boxmaps[0], index->num_boxes));
     mt = make_matchtable(index, qlen);
     partition(eq, partitions, grownindex, boxmaps);
     for (pp = 0; pp < pow(partitions, 3); pp++) {
          for (i = 0; i < index->num_boxes; i++) {
               if (bitmap_get(boxmaps[pp], i)) {
                    box = index->elements[i];
                    for (j = 0; j < eq[pp].fill_pointer; j++) {
                         comparisons++;
                         distance = infdist4(eq[pp].events[j].point, box,
                                             maxdist);
                         if (distance <= maxdist) {
                              query_page_number = eq[pp].events[j].box_number;
                              markings++;
                              matchtable_mark(mt, box.box_number, 
                                              query_page_number);
                         }
                    }
               }
          }
     }
     for (i = 0; i < pow(partitions, 3); i++)
          free(boxmaps[i]);
     free(boxmaps);
     fprintf(stderr, "Done.  Marked entries: %d\n", mt->marked);
     fprintf(stderr, "%ld comparisons, %ld markings\n", comparisons, markings);
     return mt;
}