#pragma warning (disable: 4786)

#include <stdio.h>
#include <sys/timeb.h>
#include <math.h>
#include <string>
#include <map>
#include "strtokenizer.h"
#include "utils.h"
#include "model.h"

using namespace std; 

#define MAX(x, y) (((x) > (y))? (x): (y))
#define MIN(x, y) (((x) < (y))? (x): (y))

// Generates an arbitary seed for the random number generator. 
int utils::get_seed(void) {
	unsigned int seed; 
	struct timeb timeBuf; 

	ftime(&timeBuf); 
	seed = (((unsigned int)timeBuf.time&0xFFFF) + (unsigned int)timeBuf.millitm)^(unsigned int)timeBuf.millitm; 

	return (unsigned int)seed; 
}

// Sets the seed of the random number generator to a specific value.
void utils::set_seed(int seed) {
	srand(seed); 
}

// Returns a (float) pseudo random number in the interval [0, 1]. 
double utils::get_rand(void) {
	return (double)rand() / RAND_MAX; 
}

void utils::print_error (char *fmt, ...) {
	va_list args; 
	
	va_start(args, fmt); 
	fprintf(stderr, "ERROR: ");  
	vfprintf(stderr, fmt, args); 
	va_end(args); 
}

int utils::read_and_parse(string filename, model *pmodel) {
    // open file <model>.others to read:
    // alpha=?
    // beta=?
    // ntopics=?	the number of topics
    // ndocs=?		the number of documents
    // nwords=?		the vocabulary size
	// nauthors=?	the number of unique authors
    // liter=?		the Gibbs sampling iteration at which the model was saved 
    
    FILE *fin; 

    if ((fin = fopen(filename.c_str(), "r")) == NULL) {
		print_error("Cannot open file: %s\n", filename.c_str()); 
		return EXIT_FAILURE; 
    }
    
    char buff[BUFF_SIZE_SHORT]; 
    string line; 
    
    while (fgets(buff, BUFF_SIZE_SHORT - 1, fin)) {
		line = buff; 
		strtokenizer strtok(line, "= \t\r\n"); 
		int count = strtok.count_tokens(); 
	
		if (count != 2) { // invalid, ignore this line
			continue; 
		}
	
		if (strtok.token(0) == "alpha") {
			pmodel->alpha = atof(strtok.token(1).c_str()); 	
		} else if (strtok.token(0) == "beta") {	    
			pmodel->beta = atof(strtok.token(1).c_str()); 
		} else if (strtok.token(0) == "ntopics") {
			pmodel->K = atoi(strtok.token(1).c_str()); 
		} else if (strtok.token(0) == "ndocs") {	   
			pmodel->M = atoi(strtok.token(1).c_str()); 
		} else if (strtok.token(0) == "nwords") {
			pmodel->V = atoi(strtok.token(1).c_str()); 
		} else if (strtok.token(0) == "nauthors") {
			pmodel->A = atoi(strtok.token(1).c_str());
		} else if (strtok.token(0) == "liter") {
			pmodel->liter = atoi(strtok.token(1).c_str()); 
		} else if (strtok.token(0) == "seed") {
			pmodel->seed = atoi(strtok.token(1).c_str()); 
		} else {
			;	// do nothing
		}
    }
    
    fclose(fin); 
    
    return EXIT_SUCCESS; 
}

string utils::generate_model_name(int iter) {
    string model_name = "model-";

    char buff[BUFF_SIZE_SHORT];
    
    if (0 <= iter && iter < 10) {
		sprintf(buff, "0000%d", iter);
    } else if (10 <= iter && iter < 100) {
		sprintf(buff, "000%d", iter);
    } else if (100 <= iter && iter < 1000) {
		sprintf(buff, "00%d", iter);
    } else if (1000 <= iter && iter < 10000) {
		sprintf(buff, "0%d", iter);
    } else {
		sprintf(buff, "%d", iter);
    }
    
	model_name += (iter >= 0)? buff: "final"; 
    
    return model_name;
}

void utils::quicksort(vector<pair<int, double> > &vect, int left, int right) {
    int l_hold, r_hold;
    pair<int, double> pivot;
    
    l_hold = left;
    r_hold = right;    
    int pivotidx = left;
    pivot = vect[pivotidx];

    while (left < right) {
		for (; vect[right].second <= pivot.second && left < right; right--) {
			; 
		}

		if (left != right) {
			vect[left] = vect[right];
			left++;
		}

		for (; vect[left].second >= pivot.second && left < right; left++) {
			; 
		}

		if (left != right) {
			vect[right] = vect[left];
			right--;
		}
    }

    vect[left] = pivot;
    pivotidx = left;
    left = l_hold;
    right = r_hold;
    
    if (left < pivotidx) {
		quicksort(vect, left, pivotidx - 1);
    }
    if (right > pivotidx) {
		quicksort(vect, pivotidx + 1, right);
    }    
}

/*******************************************************************************
* Purpose: 
*	beta_pdf evaluates the Beta PDF. 
* 
* Discussion: 
*	beta_pdf(x, a, b) = x^(a - 1) * (1 - x)^(b - 1) / beta(a, b)
*	a = b = 1 yields the Uniform distribution on [0, 1].
*	a = b = 1/2 yields Arcsin distribution.
*		b = 1 yields the power functin distribution. 
*	a = b -> Infinity tends to the Normal distribtuion. 
* 
* Parameters: 
*	x, double, 0.0 <= x <= 1.0. 
*	a, double, 0.0 <= a. 
*	b, double, 0.0 <= b. 
* 
* Returns: 
*	double, the value of the PDF for Beta distribution. 
*******************************************************************************/
double utils::beta_pdf(double x, double a, double b) {
	double pdf = 0.0; 

	if ((0.0 <= x) && (x <= 1.0)) {
		pdf = pow(x, a - 1.0) * pow(1.0 - x, b - 1.0) / beta(a, b); 
	}

	return pdf; 
}

/*******************************************************************************
* Purpose: 
*	beta returns the value of the Beta function. 
* 
* Discussion: 
*	beta(x, y) = (Gamma(x) * Gamma(y)) / Gamma(x + y)
*	beta(x, y) = beta(y, x)
*	beta(x, y) = Integral (0 <= T <= 1) T^(x - 1) (1 - T)^(y - 1) dT
*
* Parameters: 
*	x, double, 0.0 < x.
*	y, double, 0.0 < y. 
* 
* Returns: 
*	double, the value of the Beta function.
*******************************************************************************/
double utils::beta(double x, double y) {
	if ((x <= 0.0) || (y <= 0.0)) {
		print_error("\nBeta - Fatal error!\n"); 
		print_error("  Both x and y must be greater than 0.\n"); 
		exit(EXIT_FAILURE); 
	}
	
	return (exp(gamma_log(x) + gamma_log(y) - gamma_log(x + y)));
}

/*******************************************************************************
* Purpose: 
*	gamma_log calculates the natural logarithm of Gamma(x) for positive x. 
* 
* Discussion: 
*	The program uses rational functions that theoretically approximate 
*	log(gamma(x)) to at least 18 significant decimal digits. The 
*	approximation for 12 < x is from Hart et al, while approximations 
*	for x < 12.0 are similar to those in Cody and Hillstrom, but are 
*	unpublished. 
* 
*	The accuracy achieved depends on the arithmetic system, the compiler, 
*	intrinsic functions, and proper selection of the machine-dependent 
*	constants. 
* Reference:
*	[1] William Cody, Kenneth Hillstrom, 1967. Chebyshev Approximations for 
*		the Natural Logarithm of the Gamma Function. Mathematics of Computation, 
*		Vol. 21, pp. 198-203. 
*	[2] Kenneth Hillstrom, 1969. ANL/AMD Program ANLC366S, DGAMMA/DLGAMA.
*	[3] Hart, Ward Cheney, Charles Lawson, Maehly, Charles Mesztenyi, John Rice, 
*		Thacher, Witzgall, 1968. Computer Approximations. Wiley and sons, New 
*		York.
*
* Parameters: 
*	x, double, 0.0 < x. 
* 
* Returns: 
*	double, the logarithm of the Gamma function of x. 
*	If x <= 0.0 or if overflow would occur, the program returns the value XINF, 
*	the largest representable floating point number. 
* 
* Machine-dependent constants:
*	BETA   - radix for the floating-point representation.
*	MAXEXP - the smallest positive power of BETA that overflows.
*	XBIG   - largest argument for which LN(GAMMA(X)) is representable in the 
*		machine, i.e., the solution to the equation 
*		LN(GAMMA(XBIG)) = BETA**MAXEXP.
*
*	FRTBIG - Rough estimate of the fourth root of XBIG
*
*	Approximate values for some important machines are:
*                             BETA      MAXEXP         XBIG
*	CRAY-1        (S.P.)        2        8191       9.62E+2461
*	Cyber 180/855
*	  under NOS   (S.P.)        2        1070       1.72E+319
*	IEEE (IBM/XT,
*	  SUN, etc.)  (S.P.)        2         128       4.08E+36
*	IEEE (IBM/XT,
*	  SUN, etc.)  (D.P.)        2        1024       2.55D+305
*	IBM 3033      (D.P.)       16          63       4.29D+73
*	VAX D-Format  (D.P.)        2         127       2.05D+36
*	VAX G-Format  (D.P.)        2        1023       1.28D+305
*
*                           FRTBIG
*	CRAY-1        (S.P.)   3.13E+615
*	Cyber 180/855
*	  under NOS   (S.P.)   6.44E+79
*	IEEE (IBM/XT,
*	  SUN, etc.)  (S.P.)   1.42E+9
*	IEEE (IBM/XT,
*	  SUN, etc.)  (D.P.)   2.25D+76
*	IBM 3033      (D.P.)   2.56D+18
*	VAX D-Format  (D.P.)   1.20D+9
*	VAX G-Format  (D.P.)   1.89D+76
*******************************************************************************/
double utils::gamma_log(double x) {
	double c[7] = {
		-1.910444077728E-03,
		8.4171387781295E-04,
		-5.952379913043012E-04,
		7.93650793500350248E-04,
		-2.777777777777681622553E-03,
		8.333333333333333331554247E-02,
		5.7083835261E-03
	};
	double corr;
	double d1 = - 5.772156649015328605195174E-01;
	double d2 =   4.227843350984671393993777E-01;
	double d4 =   1.791759469228055000094023;
	double frtbig = 1.42E+09; 
	double p1[8] = {
		4.945235359296727046734888,
		2.018112620856775083915565E+02,
		2.290838373831346393026739E+03,
		1.131967205903380828685045E+04,
		2.855724635671635335736389E+04,
		3.848496228443793359990269E+04,
		2.637748787624195437963534E+04,
		7.225813979700288197698961E+03
	};
	double p2[8] = {
		4.974607845568932035012064,
		5.424138599891070494101986E+02,
		1.550693864978364947665077E+04,
		1.847932904445632425417223E+05,
		1.088204769468828767498470E+06,
		3.338152967987029735917223E+06,
		5.106661678927352456275255E+06,
		3.074109054850539556250927E+06
	};
	double p4[8] = {
		1.474502166059939948905062E+04,
		2.426813369486704502836312E+06,
		1.214755574045093227939592E+08,
		2.663432449630976949898078E+09,
		2.940378956634553899906876E+010,
		1.702665737765398868392998E+011,
		4.926125793377430887588120E+011,
		5.606251856223951465078242E+011
	};
	double pnt68 = 0.6796875;
	double q1[8] = {
		6.748212550303777196073036E+01,
		1.113332393857199323513008E+03,
		7.738757056935398733233834E+03,
		2.763987074403340708898585E+04,
		5.499310206226157329794414E+04,
		6.161122180066002127833352E+04,
		3.635127591501940507276287E+04,
		8.785536302431013170870835E+03
	};
	double q2[8] = {
		1.830328399370592604055942E+02,
		7.765049321445005871323047E+03,
		1.331903827966074194402448E+05,
		1.136705821321969608938755E+06,
		5.267964117437946917577538E+06,
		1.346701454311101692290052E+07,
		1.782736530353274213975932E+07,
		9.533095591844353613395747E+06
	};
	double q4[8] = {
		2.690530175870899333379843E+03,
		6.393885654300092398984238E+05,
		4.135599930241388052042842E+07,
		1.120872109616147941376570E+09,
		1.488613728678813811542398E+010,
		1.016803586272438228077304E+011,
		3.417476345507377132798597E+011,
		4.463158187419713286462081E+011
	};
	double res;
	double sqrtpi = 0.9189385332046727417803297;
	double xbig = 4.08E+36;
	double xden;
	double xm1;
	double xm2;
	double xm4;
	double xnum;
	double xsq;

	//  Return immediately if the argument is out of range.
	if ((xbig < x) || (x <= 0.0) ){
		return r8_huge();
	}

	if (x <= r8_epsilon()) {
		res = - log ( x ); 
	} else if (x <= 1.5) {
		if (x < pnt68) {
			corr = -log(x); 
			xm1 = x;
		} else {
			corr = 0.0; 
			xm1 = (x - 0.5) - 0.5; 
		}

		if ((pnt68 <= x) || (x <= 0.5)) { 
			xden = 1.0; 
			xnum = 0.0; 

			for (int i = 0; i < 8; i++) {
				xnum = xnum * xm1 + p1[i]; 
				xden = xden * xm1 + q1[i]; 
			}

			res = corr + (xm1 * (d1 + xm1 * (xnum / xden))); 
		} else { 
			xm2 = (x - 0.5) - 0.5; 
			xden = 1.0; 
			xnum = 0.0; 
			
			for (int i = 0; i < 8; i++) { 
				xnum = xnum * xm2 + p2[i];
				xden = xden * xm2 + q2[i]; 
			}

			res = corr + xm2 * (d2 + xm2 * (xnum / xden)); 
		}
	} else if (x <= 4.0) {
		xm2 = x - 2.0; 
		xden = 1.0; 
		xnum = 0.0; 
		
		for (int i = 0; i < 8; i++) { 
			xnum = xnum * xm2 + p2[i]; 
			xden = xden * xm2 + q2[i]; 
		}

		res = xm2 * (d2 + xm2 * (xnum / xden)); 
	} else if (x <= 12.0) {
		xm4 = x - 4.0; 
		xden = - 1.0; 
		xnum = 0.0; 
		
		for (int i = 0; i < 8; i++) {
			xnum = xnum * xm4 + p4[i]; 
			xden = xden * xm4 + q4[i]; 
		}

		res = d4 + xm4 * ( xnum / xden ); 
	} else {
		res = 0.0; 

		if (x <= frtbig) {
			res = c[6]; 
			xsq = x * x; 

			for (int i = 0; i < 6; i++) { 
				res = res / xsq + c[i]; 
			}
		}

		res = res / x; 
		corr = log (x); 
		res = res + sqrtpi - 0.5 * corr; 
		res = res + x * ( corr - 1.0 ); 
	}

	return res;
}

/******************************************************************************
* Purpose: 
*	r8_epsilon returns the R8 round off unit. 
* 
* Discussion: 
*	R8_EPSILON is a number R which is a power of 2 with the property that, 
*	to the precision of the computer's arithmetic, 
*		1 < 1 + R
*	but 
*		1 = (1 + R / 2)
* 
* Parameters: 
*	void
* 
* Returns: 
*	double, the floating point round-off unit. 
******************************************************************************/
double utils::r8_epsilon(void) {
	double r = 1.0; 

	while (1.0 < (double) (1.0 + r)) {
		r /= 2.0; 
	}

	return (2.0 * r);
}

/*******************************************************************************
* Purpose: 
*	r8_huge returns a "huge" R8. 
* 
* Discussion: 
*	HUGE_VAL is the largest representable legal real number, and is usually 
*	defined in math.h, or sometimes in stdlib.h. 
* 
* Parameters: 
*	void
* 
* Returns: 
*	double, a "huge" real value. 
*******************************************************************************/
double utils::r8_huge(void) {
	return HUGE_VAL; 
}

// return the maximal value
double utils::get_max(vector<double> v) {
	double ret = -r8_huge(); 

	for (int i = 0; i < v.size(); i++) {
		ret = MAX(ret, v[i]); 
	}

	return ret; 
}

// return the minimal value
double utils::get_min(vector<double> v) {
	double ret = r8_huge(); 

	for (int i = 0; i < v.size(); i++) {
		ret = MIN(ret, v[i]); 
	}

	return ret; 
}

