#include <string>
#include <vector>
#include <math.h>
#include "rating.h"

using namespace std;

double stdnormal_inv(double p) {
	/* Calculates approximation to the inverse of the cumulative
		standard normal distribution function
	
		http://www.pitt.edu/~wpilib/statfaq/gaussfaq.html */
	
	if (p <= 0) return -1e300;
	if (p >= 1) return 1e300;
	
	double t;
	
	if (p <= .5) {
		t = sqrt(-2.0 * log(p));
		return -t + (2.515517 + t * (0.802853 + t * 0.010328)) / (1 + t * (1.432788 + t * (0.189269 + t * 0.001308)));
	} else {
		t = sqrt(-2.0 * log(1 - p));
		return t - (2.515517 + t * (0.802853 + t * 0.010328)) / (1 + t * (1.432788 + t * (0.189269 + t * 0.001308)));
	}
}

#ifdef NOT_DEFINED
/*
* The inverse standard normal distribution.
*
*   Author:      Peter J. Acklam <pjacklam@online.no>
*   URL:         http://home.online.no/~pjacklam
*
* This function is based on the MATLAB code from the address above,
* translated to C, and adapted for our purposes.
*/
double stdnormal_inv(double p)
{
	const double a[6] = {
		-3.969683028665376e+01,  2.209460984245205e+02,
		-2.759285104469687e+02,  1.383577518672690e+02,
		-3.066479806614716e+01,  2.506628277459239e+00
	};
	const double b[5] = {
		-5.447609879822406e+01,  1.615858368580409e+02,
		-1.556989798598866e+02,  6.680131188771972e+01,
		-1.328068155288572e+01
	};
	const double c[6] = {
		-7.784894002430293e-03, -3.223964580411365e-01,
		-2.400758277161838e+00, -2.549732539343734e+00,
		4.374664141464968e+00,  2.938163982698783e+00
	};
	const double d[4] = {
		7.784695709041462e-03,  3.224671290700398e-01,
		2.445134137142996e+00,  3.754408661907416e+00
	};
	
	register double q, t, u;
	
	if (isnan(p) || p > 1.0 || p < 0.0) return NAN;
	if (p == 0.0) return -HUGE_VAL;
	if (p == 1.0) return  HUGE_VAL;
	q = p <? 1-p;
	if (q > 0.02425) {
		/* Rational approximation for central region. */
		u = q-0.5;
		t = u*u;
		u = u*(((((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4])*t+a[5])
		/(((((b[0]*t+b[1])*t+b[2])*t+b[3])*t+b[4])*t+1);
	} else {
		/* Rational approximation for tail region. */
		t = sqrt(-2*log(q));
		u = (((((c[0]*t+c[1])*t+c[2])*t+c[3])*t+c[4])*t+c[5])
		/((((d[0]*t+d[1])*t+d[2])*t+d[3])*t+1);
	}
	/* The relative error of the approximation has absolute value less
	than 1.15e-9.  One iteration of Halley's rational method (third
	order) gives full machine precision... */
	t = stdnormal_cdf(u)-q;    /* error */
	t = t*M_SQRT2PI*exp(u*u/2);   /* f(u)/df(u) */
	u = u-t/(1+u*t/2);     /* Halley's method */
	
	return (p > 0.5 ? -u : u);
}
#endif

double winp(int Rating1, int Vol1, int Rating2, int Vol2) {
	return erf((Rating1-Rating2)/sqrt(2.0*(Vol1*Vol1 + Vol2*Vol2)))*0.5 + 0.5;
}

int calcRatings(vector<coder> &c) {
	int sum = 0;
	long long s2 = 0, v2 = 0;
	
	// compute sums
	for (int i = 0; i < c.size(); i++) {
		sum += c[i].rating;
		s2 += c[i].rating*c[i].rating;
		v2 += c[i].vol*c[i].vol;
	}
	
	// now calculate competition factor, etc
	double avg, cf;
	
	avg = sum/(double)c.size();
	cf = sqrt(v2/(double)c.size() + (s2 - c.size()*avg*avg)/(c.size()-1));
	
	/* calculate expected rank */
	vector<double> erank(c.size(), 0.5);
	for (int i = 0; i < c.size(); i++)
		for (int j = 0; j < c.size(); j++)
			erank[i] += winp(c[i].rating, c[i].vol, c[j].rating, c[j].vol);
	
	for (int i = 0; i < c.size(); i++) {
		double eperf = -stdnormal_inv((erank[i] - 0.5)/c.size());
		
		/* not using ties for this simulation */
		double aperf = -stdnormal_inv((i+0.5)/c.size());
		
		double perfas = c[i].rating + cf*(aperf-eperf);
		
		double wt = 1/(1-(.42/(c[i].tp + 1)+.18))-1;
		if (c[i].rating >= 2500)
			wt *= .8;
		else if (c[i].rating >= 2000)
			wt *= .9;
			
		int cap = 150 + 1500/(1 + c[i].tp);
		
		int rating = (int)round((c[i].rating + wt*perfas)/(1+wt));
		
		rating <?= c[i].rating+cap;
		rating >?= c[i].rating-cap;
		
		c[i].vol = (int)round(sqrt( (rating-c[i].rating)*(rating-c[i].rating)/wt +
			c[i].vol*c[i].vol/(wt+1)));
		c[i].rating = rating;
	}
	return 0;
}