/** @file evaluation.cpp
 *
 *	Definition (implementation) of evaluation module interface.
 *
 *	See evaluation.h for more info.
 *
 *	(C) Ceslav Przywara 2008, MFF UK Prague
 *
 */

// Evaluation requires some advanced math...
#include <cmath>
// ...and floating point operations.
#include <cfloat>

#include "evaluation.h"
// ith_bit() needed
#include "utils.h"

namespace ace {

/** "Override" of std::log10 (return DBL_MIN, if given x is zero).
 *	@param x Value to be logged :)
 *	@return log10 of given x.
 */
double _log10(double x) {
	return (std::abs(x - 0.0) < DBL_EPSILON) ? DBL_MIN : std::log10(x);
}

/* Manning.
double _L(double k, double n, double x) {
	return k*_log10(x) + (n - k)*_log10(1 - x);
}
*/

/**
 */
std::vector<std::vector<size_t> > EvaluationTables::_bits_for_cont_table;
/**
 */
std::vector<std::vector<size_t> > EvaluationTables::_bits_for_inv_table;

/*	Interface item implementation - for info see description in header file.
 */
void EvaluationTables::init() {
	// Requirement: NGram::init() already called.
	//
	_bits_for_cont_table.resize(NGram::types_count());
	_bits_for_inv_table.resize(NGram::types_count());
	//
	for ( ngram_type_t outer = 0; outer < NGram::types_count(); ++outer ) {
		// Do it in the most naive way, it won't have any significant impact
		// on program speed anyway...
		for ( ngram_type_t inner = 0; inner < NGram::types_count(); ++inner ) {
			// Superset exclusive: Inner non-zero bits must be exclusive superset of outer non-zero bits.
			if ( (inner ^ outer) && (bits_in<ngram_type_t>(inner) == bits_in<ngram_type_t>(inner | outer)) ) {
				// if ( ^ ) -> ther're different and if also () == () -> inner bits are superset of outer bits
				_bits_for_cont_table[outer].push_back(inner);
			}
			// Outer non-zero bits must meet inner zero bits (but 0 is excluded)
			if ( inner && !(outer & inner) ) {
				_bits_for_inv_table[outer].push_back(inner);
			}
		}
	}
}

/*	Interface item implementation - for info see description in header file.
 */
// Init contigency table with frequency values!
EvaluationTables::EvaluationTables(const freq_table_t frequencies): _freq(frequencies), _cont(_freq), _inv(), _exp() {
	
	// Now adjust values in contigency table.
	ngram_type_t type = NGram::full_ngram_type();
	do {
		--type; // Values of contigency and *-frequency tables for full N-gram type are always equal.
		// 
		for ( size_t i = 0; i < _bits_for_cont_table[type].size(); ++i ) {
			_cont[type] -= _cont[_bits_for_cont_table[type][i]];
		}
	} while (type != 0);	
	
	// Init inverted frequencies table with contigency table 0-index value.
	_inv.resize(_cont.size(), _cont[0]);
	// Now adjust values in inverted frequencies table.
	for ( ngram_type_t type = 0; type < _freq.size(); ++type ) {
		for ( ngram_type_t i = 0; i < _bits_for_inv_table[type].size(); ++i ) {
			_inv[type] += _cont[_bits_for_inv_table[type][i]];
		}
	}
	// Count expected values (for likelihood ratio and chi square test).
	_exp.resize(_freq.size(), 1.0);
	for ( ngram_type_t type = 0; type < _freq.size(); ++type ) {
		_exp[type] *= _freq[0]; // #N
		for ( ngram_type_t i = 0; i < NGram::n(); i++ ) {
			if ( ith_bit(type, i) ) {
				_exp[type] *= static_cast<double>(_freq[1 << i]) / _freq[0]; // p(x) = #x/#N
			} else {
				_exp[type] *= static_cast<double>(_inv[1 << i]) / _freq[0];
			}
		}
	}
}

namespace eval {

/*	Interface item implementation - for info see description in header file.
 */
double chi_square_test(const EvaluationTables& tables) {

	double chi_square = 0.0;

	for ( ngram_type_t type = 0; type < tables._freq.size(); ++type ) {
		chi_square += std::pow(tables._cont[type] - tables._exp[type], 2) / tables._exp[type];
	}

	return chi_square;
}

/*	Interface item implementation - for info see description in header file.
 */
double log_likelihood_ratio(const EvaluationTables& tables) {

	double likelihood = 0.0;

	for ( ngram_type_t type = 0; type < tables._freq.size(); ++type ) {
		likelihood += tables._cont[type] * _log10(tables._cont[type]/ tables._exp[type]);
	}
	
	return 2 * likelihood; //

	// Based on Manning.
	/*
	double p = static_cast<double>(_tables._freq[1]) / static_cast<double>(_tables._freq[0]),
		p1 = static_cast<double>(_tables._freq[3]) / static_cast<double>(_tables._freq[2]),
		p2 = static_cast<double>(_tables._freq[1] - _tables._freq[3]) / static_cast<double>(_tables._freq[0] - _tables._freq[2]);
	
	double
		L11 = // L( f(X,Y), f(X,*), p )
			_L(_tables._freq[3], _tables._freq[2], p), 
		L12 = // L( f(*,Y)-f(X,Y), f(*,*)-f(X,*), p )
			_L(_tables._freq[1] - _tables._freq[3], _tables._freq[0] - _tables._freq[2], p),
		L21 = // L( f(X,Y), f(X,*), p1 )
			_L(_tables._freq[3], _tables._freq[2], p1),
		L22 = // L( f(*,Y)-f(X,Y), f(*,*)-f(X,*), p2 )
			_L(_tables._freq[1] - _tables._freq[3], _tables._freq[0] - _tables._freq[2], p2);

	//return -2 * ((std::log10(L11) - std::log10(L21)) + (std::log10(L12) - std::log10(L22)));
	return -2 * (L11 - L21 + L12 - L22);
	*/
}

/*	Interface item implementation - for info see description in header file.
 */
double mutual_information(const EvaluationTables& tables) {
	// For given X and Y, the function returns log2( p(XY)/ p(Y)*p(Y) ), where:
	// p(XY) = f(X,Y) / f(*,*)
	// p(X)  = f(X,*) / f(*,*)
	// p(Y)  = f(*,Y) / f(*,*)

	double pxy = static_cast<double>(tables._freq[3]) / static_cast<double>(tables._freq[0]),
		px = static_cast<double>(tables._freq[2]) / static_cast<double>(tables._freq[0]),
		py = static_cast<double>(tables._freq[1]) / static_cast<double>(tables._freq[0]);
	
	// We want log2, but there is no such function in cmath library.
	return _log10(pxy/(px*py)) / std::log10(2.0);
}

/*	Interface item implementation - for info see description in header file.
 */
double pearsons_coefficient(const EvaluationTables& tables) {
	return (tables._cont[3]*tables._cont[0] - tables._cont[2]*tables._cont[1])
		/ std::sqrt(static_cast<double>(
			(tables._cont[3] + tables._cont[2])
			* (tables._cont[3] + tables._cont[1])
			* (tables._cont[0] + tables._cont[2])
			* (tables._cont[0] + tables._cont[1])
		));
}

/*	Interface item implementation - for info see description in header file.
 */
double t_test(const EvaluationTables& tables) {
	//  (f(X,Y) - e(X,Y)) / (f(X,Y)*(1-f(X,Y)/f(*,*)))
	return (tables._freq[3] - tables._exp[3]) / std::sqrt(tables._freq[3]*(1 - static_cast<double>(tables._freq[3])/static_cast<double>(tables._freq[0])));
}

/*	Interface item implementation - for info see description in header file.
 */
double z_score(const EvaluationTables& tables) {
	//  (f(X,Y) - e(X,Y)) / (e(X,Y)*(1-e(X,Y)/f(*,*)))
	return (tables._freq[3] - tables._exp[3]) / std::sqrt(tables._exp[3]*(1 - static_cast<double>(tables._exp[3])/static_cast<double>(tables._freq[0])));
}

} // namespace eval

} // namespace ace
