#include <assert.h>
#include <stdio.h>
#include <vlcutils/error.h>
#include <vlcutils/heap.h>
#include "freqdist.h"
#include "index.h"

static int cmp_start(const void *d1, const void *d2)
{
     const MBR *r1, *r2;

     r1 = (MBR *)d1;
     r2 = (MBR *)d2;
     return r1->start < r2->start ? -1 : (r1->start == r2->start ? 0 : 1);
}

static void fill_two_boxes(const char *s, int start, int length, int numchar,
			   MBR *box1, MBR *box2, int window_size)
{
     int i, j;
     int p[ALPH_LEN], q[ALPH_LEN];
     MBR r1, r2;
     double volume1, volume2, density1, density2;

     assert(s && box1 && box2);
     assert(start >= 0);
     assert(length >= 2);
     assert(0 < window_size);

     for (i = 0; i < ALPH_LEN; i++) {
	  p[i] = 0;
	  q[i] = 0;
     }

     /* Find the points p and q corresponding to the first and last window. */
     for (i = 0; i < window_size; i++) {
	  increment(s[start + i], p);
	  increment(s[start + length - 1 + i], q);
     }
     /* Initially, r1 contains only p, and r2 contains only q. */
     for (i = 0; i < ALPH_LEN; i++) {
	  r1.lower[i] = r1.higher[i] = p[i];
	  r2.lower[i] = r2.higher[i] = q[i];
     }
     r1.start = start;
     r2.start = start + length - 1;
     r1.capacity = r2.capacity = 1;
     density1 = density2 = 0.0;

     i = r1.start;
     for (;;) {
	  volume1 = 1;
	  volume2 = 1;
	  for (j = 0; j < ALPH_LEN; j++) {
	       volume1 *= (r1.higher[j] - r1.lower[j] + 1) 
		 / (double)(window_size + 1);
	       volume2 *= (r2.higher[j] - r2.lower[j] + 1) 
		 / (double)(window_size + 1);
	  }
	  if (i + 1 == r2.start)
	       break;
	  density1 = r1.capacity / volume1;
	  density2 = r2.capacity / volume2;
	  if (volume1 < volume2) {
	       /* Compute the new point. */
	       decrement(s[i], p);
	       increment(s[i + window_size], p);
	       i++;
	       /* Make sure r1 contains the new point. */
	       for (j = 0; j < ALPH_LEN; j++) {
		    if (p[j] < r1.lower[j]) r1.lower[j] = p[j];
		    if (p[j] > r1.higher[j]) r1.higher[j] = p[j];
	       }
	       r1.capacity++;
	  } else {
	       /* Compute the new point. */
	       r2.start--;
	       decrement(s[r2.start + window_size], q);
	       increment(s[r2.start], q);
	       /* Make sure r2 contains the new point. */
	       for (j = 0; j < ALPH_LEN; j++) {
		    if (q[j] < r2.lower[j]) r2.lower[j] = q[j];
		    if (q[j] > r2.higher[j]) r2.higher[j] = q[j];
	       }
	       r2.capacity++;
	  }
     }
     r1.key = density1;
     r2.key = density2;
     memcpy(box1, &r1, sizeof(MBR));
     memcpy(box2, &r2, sizeof(MBR));
}


static int key_lt_p(const MBR *r1, const MBR *r2)
{
     return r1->key < r2->key;
}

struct index *
make_index_mhistd(FILE *input, int window_size, int num_boxes)
{
     struct index *index;
     int numchar;
     char *s;
     size_t bytes_read;
     struct heap *heap;
     MBR *r;
     int i;

     /* Allocate index */
     index = malloc(sizeof(struct index));
     if (!index)
	  fatal_perror("malloc");
     index->window_size = window_size;
     index->dimensions = ALPH_LEN;
     index->num_boxes = num_boxes;
     index->box_capacity = -1; /* Adaptive */
     index->memory_usage = 
	  sizeof(struct index) + index->num_boxes * sizeof(MBR);
     index->elements = malloc(index->num_boxes * sizeof(MBR));
     if (!index->elements)
	  fatal_perror("malloc");

     printf("Determining input size.\n");
     numchar = file_length(input);

     printf("Reading input.\n");
     s = malloc(numchar);
     if (!s)
       fatal_perror("malloc");
     bytes_read = fread(s, 1, numchar, input);
     if (bytes_read != numchar)
       fatal_error("Read only %d bytes, not the expected %d.\n", bytes_read,
		   numchar);

     printf("Computing boxes...\n");
     fill_two_boxes(s, 0, numchar - window_size + 1, numchar, 
		    &index->elements[0], &index->elements[1], window_size);
     heap = make_heap(num_boxes, key_lt_p);
     heap_insert(&index->elements[0], heap);
     heap_insert(&index->elements[1], heap);

     for (i = 2; i < num_boxes; i++) {
	  r = heap_pop(heap);
	  while (r->capacity <= 2)
	       r = heap_pop(heap);
	  /* printf("...%d (capacity %d)\n", i, r->capacity); */
	  fill_two_boxes(s, r->start, r->capacity, numchar,
			 r, &index->elements[i], window_size);
	  heap_insert(r, heap);
	  heap_insert(&index->elements[i], heap);
     }

     free(s);
     s = NULL;

     /* Sort the index by start position */
     qsort(index->elements, index->num_boxes, sizeof(MBR), cmp_start);
     return index;
}