#include "MpUtil.h"
#include "UncoveredSet.cpp"

/*

Find the Uncovered Set and mixed strategy equilibrium to a two-dimensional
two-candidate model of Downsian competition using a finite grid to approximate
the strategy space.

*/

int main()
{

	long RandomSeed = -197347;
	const int N = 101; 		// Number of Voters
	const int J = Sqr(81);	// Number of Locations
	Timer tm;


	/* Generate Voter Ideal Points from Triangle Distribution */
	Print("1: Generating the Voter Ideal Points");
	tm.Start();
	Vector<double> q1(N);	// Ideal points of voters, coordiante 1
	Vector<double> q2(N);	// Ideal points of voters, coordinate 2
	for(int n=0; n<N; n++)
	{
		q1(n) = UniformSample(RandomSeed);
		q2(n) = UniformSample(RandomSeed);
		if(q2[n] > q1[n]) Swap(q1(n),q2(n));
	}
	tm.Stop();
	Print("");


	/* Create a Finite Grid */
	Print("2: Creating the Grid");
	tm.Start();
	double x1low  = 0.0;
	double x1high = 1.0;
	double x2low  = 0.0;
	double x2high = 1.0;
	Vector<double> x1(J);	// Location, Coordinate 1
	Vector<double> x2(J);	// Location, Coordinate 2
	int Curr = 0;
	for(int j=0; j<sqrt(J); j++) for(int k=0; k<sqrt(J); k++)
	{
		x1[Curr] = x1low + (x1high - x1low) * j / (sqrt(J) - 1.0);
		x2[Curr] = x2low + (x2high - x2low) * k / (sqrt(J) - 1.0);
		Curr++;
	}
	tm.Stop();
	Print("");


	/* Compute the Preference Matrix */
	tm.Start();
	Print("3: Loading the Preference Matrix");
	Matrix<int> A(J,J);
	cout << "Progess (out of " << J << "):\n";
	for(int j=0; j<J; j++)
	{
		if(j % 100 == 0) cout << j << "\n";
		for(int k=0; k<J; k++)
		{
			int vcurr = 0;
	        for(int n=0; n<N; n++)
	        {
	        	if(Sqr(q1[n] - x1[j]) + Sqr(q2[n] - x2[j]) < Sqr(q1[n] - x1[k]) + Sqr(q2[n] - x2[k])) vcurr += 2;
	           	if(Sqr(q1[n] - x1[j]) + Sqr(q2[n] - x2[j]) == Sqr(q1[n] - x1[k]) + Sqr(q2[n] - x2[k])) vcurr += 1;
			}
			if(vcurr >  N) A(j,k) =  1;
			if(vcurr <  N) A(j,k) = -1;
			if(vcurr == N) A(j,k) =  0;
		}
	}
	tm.Stop();
	Print("");


	/* Compute the Uncovered Set */
	Print("4: Computing the Uncovered Set");
	tm.Start();
	Vector<int> UC(J);
	UncoveredSet(A,UC);
	Print("UC",UC);
	tm.Stop();
	Print("");


	/* Compute the Mixed Strategy Equilibrium */
	Print("5: Finding the Mixed Strategy Equilibrium");
	Stack<int> UCStack(J);
	for(int j=0; j<J; j++) if(UC(j)) UCStack.Push(j);
	int NumberUncovered = UCStack.Size();
	Print("NumberUncovered",NumberUncovered);
	Matrix<int> AUC(NumberUncovered,NumberUncovered);
	for(int i=0; i<NumberUncovered; i++) for(int j=0; j<NumberUncovered; j++) AUC(i,j) = A(UCStack(i),UCStack(j));
	Print("AUC",AUC);

	Vector<double> c(NumberUncovered+1);
	for(int j=0; j<NumberUncovered; j++) c[j] = 0.0;
	c[NumberUncovered] = -1.0;
	Matrix<double> Ae(1,NumberUncovered+1);
	for(int j=0; j<NumberUncovered; j++) Ae(0,j) = 1.0;
	Ae(0,NumberUncovered) = 0.0;
	Vector<double> be(1);
	be[0] = 1.0;
	Matrix<double> An(NumberUncovered,NumberUncovered+1);
	for(int i=0; i<NumberUncovered; i++) for(int j=0; j<NumberUncovered; j++) An(i,j) = -AUC(i,j);
	for(int j=0; j<NumberUncovered; j++) An(j,NumberUncovered) = -1.0;
	Vector<double> bn(NumberUncovered);
	for(int j=0; j<NumberUncovered; j++) bn[j] = 0.0;
	Vector<double> lb(NumberUncovered+1,0.0);
	Vector<double> ub(NumberUncovered+1,Inf);
	SymmetricMatrix<double> Q(NumberUncovered+1,Zero);
	MinQuadProg qp(c,Q,Ae,be,An,bn,lb,ub);
	qp.Solve();
	Stack<int> MSEStack(J);
	for(int j=0; j<NumberUncovered; j++) if(abs(qp.Solution(j)) > 0.0000001) MSEStack.Push(UCStack(j));
	Vector<int> MSE(J,0);
	for(int i=0; i<MSEStack.Size(); i++) MSE(MSEStack(i)) = 1;
	tm.Stop();

	ofstream out1("grid.txt");
	ofstream out2("uncovered.txt");
	ofstream out3("mse.txt");
	ofstream out4("voters.txt");

	for(int j=0; j<J; j++) out1 << x1(j) << "\t" << x2(j) << "\n";
	for(int i=0; i<UCStack.Size(); i++) out2 << x1(UCStack(i)) << "\t" << x2(UCStack(i)) << "\n";
	for(int i=0; i<MSEStack.Size(); i++) out3 << x1(MSEStack(i)) << "\t" << x2(MSEStack(i)) << "\n";
	for(int n=0; n<N; n++) out4 << q1(n) << "\t" << q2(n) << "\n";

	Print("End");
	return 0;
}