/* AUTORIGHTS
Copyright (C) 2007 Princeton University
      
This file is part of Ferret Toolkit.

Ferret Toolkit is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software Foundation,
Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_vector.h>
#include "LSH.h"
#include "perturb.h"

#define MISS_GAMMA	1

#define MAX_M	64
#define MAX_T	600

static inline double p_col_helper (double x)
{
	double result;
	result = 2.0*gsl_cdf_ugaussian_P(x) - 1.0;
	result += sqrt(2.0/M_PI) * (exp(-x*x/2.0)-1.0)/x;
	return result; 
}

static inline double p_col_p_helper (double x, double k)
{
	return gsl_cdf_ugaussian_P((1.0 + k) * x) - gsl_cdf_ugaussian_P(k * x);
}

int LSH_recall_init (LSH_recall_t *recall, int d_step, float d_min, float d_max,
		int M, int L, int T, float W)
{
	ptb_vec_t *ptb_score;
	ptb_vec_t *ptb_set;
	float **table; 
	double dist, r, p;
	int i, j, k;
	int l;

	ptb_score = gen_score(M);
	assert(ptb_score != NULL);
	ptb_set = type_calloc(ptb_vec_t, T);
	assert(ptb_set != NULL);
	gen_perturb_set(ptb_score, ptb_set, M, T);

	recall->d_step = d_step;
	recall->d_min = d_min;
	recall->d_max = d_max;
	recall->T = T;
	recall->table = table = type_matrix_alloc(float, T + 1, d_step);

	for (i = 0; i < d_step; i++)
	{
		dist = d_min + (d_max - d_min) * i / d_step;
		dist = W / sqrt(dist);
		p = p_col_helper(dist);
		table[0][i] = log(1.0 - pow(p, M));
		for (j = 1; j <= T; j++)
		{
			r = 1.0;
			l = 0;
			for (k = 0; k < 2 * M; k++)
			if (ptb_set[j-1].set & (1 << k))
			{
				r *= p_col_p_helper(dist, ptb_score[k].key1);
				l++;
			}
			assert(l > 0);
			l = M - l;
			r *= pow(p, l);
			table[j][i] = table[j-1][i] + log(1.0-r);
		}
		for (j = 0; j <= T; j++)
		{
			table[j][i] = 1.0 - exp(L * table[j][i]);
		}
	}

	free(ptb_score);
	free(ptb_set);
	return 0;
}

int LSH_recall_dump (LSH_recall_t *recall, CASS_FILE *fout)
{
	int ret;
	ret = cass_write_size(&recall->d_step, 1, fout);
	ret += cass_write_size(&recall->T, 1, fout);
	ret += cass_write_float(&recall->d_min, 1, fout);
	ret += cass_write_float(&recall->d_max, 1, fout);

	if (ret != 4) return CASS_ERR_IO;

	if (type_matrix_dump_stream(float, fout, recall->T + 1, recall->d_step, recall->table) != 0) return CASS_ERR_IO;

	return 0;
}

int LSH_recall_load (LSH_recall_t *recall, CASS_FILE *fin)
{
	int ret;
	cass_size_t row, col;
	ret = cass_read_size(&recall->d_step, 1, fin);
	ret += cass_read_size(&recall->T, 1, fin);
	ret += cass_read_float(&recall->d_min, 1, fin);
	ret += cass_read_float(&recall->d_max, 1, fin);

	if (ret != 4) return CASS_ERR_IO;

	type_matrix_load_stream(float, fin, &row, &col, &recall->table);

	assert(row == recall->T +1);
	assert(col == recall->d_step);

	return 0;
}