#ifndef _Games_NashEquilibrium
#define _Games_NashEquilibrium

/*

Solve for the pure-strategy Nash equilibrium of a J-player game
where the action space is a single real variable.

The game is described by a function,

Info = Func(x,y)

x is the vector of strategies, y is the returned vector on utilities,
and the function returns a value not equal to zero in case of an error.

A Nash equilibrium is a point x* such that for each j,
y(j) is maximized by x*(j) holding
all the other components of x* fixed.

*/

namespace NashEquilibriumSpace
{
	int J;					// Number of Players
	int j;					// Current Player
	double DerivStep;		// Step Size for Numerical Derivatives
	Vector<double> x;		// Current Vector of Strategies

	/* Function Describing Game */
	int (*Func)(const Vector<double> &, Vector<double> &);

	/* One Dimensional Function */
	double Func_j(double xj);
	double Func_j(double xj)
	{
		x[j] = xj;
		Vector<double> f(J);
		int Info = Func(x,f);
		if(Info) return Inf; // Invalid Point
		return -f[j];
	}
	
	/* One Dimensional Function */
	int Func_j_2(double xj, double &fj);
	int Func_j_2(double xj, double &fj)
	{
		x[j] = xj;
		Vector<double> f(J);
		int Info = Func(x,f);
		if(Info) return 1; // Invalid Point
		fj = -f[j];
		return 0;
	}

	/* Derivative of Function */
	int DFunc(const Vector<double> &x, Vector<double> &df);
	int DFunc(const Vector<double> &x, Vector<double> &df)
	{
		/* For Each Index, Compute First-Order Conditions */
		Vector<double> xCopy = x.Copy();
		Vector<double> fm(J),fp(J);
		for(int j=0; j<J; j++)
		{
			xCopy(j) -= DerivStep;
			Func(xCopy,fm);
			xCopy(j) += Two * DerivStep;
			Func(xCopy,fp);
			xCopy(j) -= DerivStep;
			df(j) = (fp(j) - fm(j)) / (Two * DerivStep);
		}
		return 0;
	}
}



class NashEquilibrium
{

	public:

	int J;				// Number of Players
	Vector<double> x;	// Vector of Strategies
	int (*Func)(const Vector<double> &, Vector<double> &);
	int (*DFunc)(const Vector<double> &, Vector<double> &); 

	double DerivStep;

	int Option; 			// Option (1 = Best Response Iterations, 2 = Guass-Siedel, 3 = Guass Jacobi)
	int Method;				// Optimization Method (1 = Brent's Method, 2 = Golden Section Search)
	double Lambda;			// Dampening Factor
	int MaxIter;			// Maximum Number of Iterations
	double TolX;			// Tolerance for x
	int Min1DPrintLevel;
	int PrintLevel;
	
	int K;					// Size of Grid Used for Each Player
	Vector<double> Low;		// Lower Bound of Grid for Each Player
	Vector<double> High;	// Upper Bound of Grid for Each Player
	
	bool AnalyticalFOC;		// Analytical First-order Conditions Provided

	/* Constructor 1: No Analytical First-order Conditions Provided */
	NashEquilibrium(int (*Func)(const Vector<double> &, Vector<double> &), int J)
	{

		this->J = J;
		this->Func = Func;
		x = Vector<double>(J,Zero);
		DerivStep = MachEps3;

		Option = 1;
		Method = 2;
		Lambda = One;
		MaxIter = 100;
		TolX = 0.00000001;
		Min1DPrintLevel = 1;
		PrintLevel = 2;
		AnalyticalFOC = False;
		
		K = 11;
		Low = Vector<double>(J,-One);	Low.Name  = "NashEquilibrium::Low";
		High = Vector<double>(J,One);	High.Name = "NashEquilibrium::High";
	}

	/* Constructor 2: Analytical First-order Conditions Provided */
	NashEquilibrium(int (*Func)(const Vector<double> &, Vector<double> &), int (*DFunc)(const Vector<double> &, Vector<double> &), int J)
	{
		this->J = J;
		this->Func = Func;
		this->DFunc = DFunc;
		x = Vector<double>(J,Zero);
		DerivStep = MachEps3;

		Option = 1;
		Method = 2;
		Lambda = One;
		MaxIter = 100;
		TolX = 0.00000001;
		Min1DPrintLevel = 1;
		PrintLevel = 2;
		AnalyticalFOC = True;
		
		K = 11;
		Low = Vector<double>(J,-One);	Low.Name  = "NashEquilibrium::Low";
		High = Vector<double>(J,One);	High.Name = "NashEquilibrium::High";
	}

	void Solve()
	{

		Vector<double> xPrev(J);
		Vector<double> xNew(J);

		if(Option == 1)
		{
		
			/* Set Joint First-Order Conditions Equal to Zero */
			
			if(AnalyticalFOC)
			{
				if(Method == 1)
				{
					/* Newton's Methods */
					SolveNewton< Matrix<double> > sn(DFunc,J);
					sn.NormType = 2;
					sn.Method = 5;
					sn.PrintLevel = PrintLevel;
					Copy(sn.x,x);
					sn.Solve();
					Copy(x,sn.x);
				
				}
				else
				{
					/* Nelder Mead */
					SolveNelderMead sn(DFunc,J);
					sn.TolX = 0.000000000001;
					sn.ScaleX = 0.1;
					sn.NormType = 2;
					sn.PrintLevel = PrintLevel;
					Copy(sn.x,x);
					sn.Solve();
					Copy(x,sn.x);
				}
			}
			else
			{
				if(Method == 1)
				{
				
					/* Newton's Methods */
				
					NashEquilibriumSpace::J         = J;
					NashEquilibriumSpace::Func      = Func;
					NashEquilibriumSpace::DerivStep = DerivStep;

					SolveNewton< Matrix<double> > sn(NashEquilibriumSpace::DFunc,J);
					sn.NormType = 2;
					sn.Method = 5;
					sn.PrintLevel = PrintLevel;
					Copy(sn.x,x);
					sn.Solve();
					Copy(x,sn.x);
				
				}
				else
				{
				
					/* Nelder Mead */
					
					NashEquilibriumSpace::J         = J;
					NashEquilibriumSpace::Func      = Func;
					NashEquilibriumSpace::DerivStep = DerivStep;

					SolveNelderMead sn(NashEquilibriumSpace::DFunc,J);
					sn.TolX = 0.000000000001;
					sn.ScaleX = 0.1;
					sn.NormType = 2;
					sn.PrintLevel = PrintLevel;
					Copy(sn.x,x);
					sn.Solve();
					Copy(x,sn.x);
				}
			}
		}
		else
		{
		
			/* Solve Using Guass-Seidel or Guass Jacobi */

			NashEquilibriumSpace::J     = J;
			NashEquilibriumSpace::x     = Vector<double>(J);
			NashEquilibriumSpace::Func  = Func;

			Min1D m1(NashEquilibriumSpace::Func_j);
			m1.Option = 2;
	
			Vector< Vector<double> > Grid1(J);
			Grid1.Name = "NashEquilibrium::Solve::Grid1";
			for(int j=0; j<J; j++) Grid1(j) = Grid(Low(j),High(j),K);
			
			for(int Iter=1; Iter<=MaxIter; Iter++)
			{
				Copy(xPrev,x);

				if(PrintLevel > 0) Print("Iter",Iter);
				if(PrintLevel > 0) Print("xPrev",x);

				if(Option == 2)
				{
					/* Guass-Siedel */
				
					/* For Each Index.. */
					for(int j=0; j<J; j++)
					{
						/* Apply Max Transformation */
						Copy(NashEquilibriumSpace::x,x);
						NashEquilibriumSpace::j = j;
						m1.Solve(Grid1(j));
						x(j) = m1.x;
					}
				}
				else
				{
					/* Guass-Jacobi */

					/* For Each Index.. */
					for(int j=0; j<J; j++)
					{
						/* Apply Max Transformation */
						Copy(NashEquilibriumSpace::x,xPrev);
						NashEquilibriumSpace::j = j;
						m1.Solve(Grid1(j));
						x(j) = m1.x;
					}
				}

				double xMax = NormL1(x - xPrev);
				xNew = Lambda * x + (One - Lambda) * xPrev;
				if(PrintLevel > 0) Print("x",x);
				if(PrintLevel > 0) Print("xNew",xNew);
				if(PrintLevel > 0) Print("xMax",xMax);
				Copy(x,xNew);

				if(xMax < TolX)
				{
					if(PrintLevel > 0) Print("Succesfully Exiting NashEquilibrium: Within Tolerance");
					if(PrintLevel > 0) Print("x",x);
					if(PrintLevel > 0) Print("xMax",xMax);
					return;
				}
			}

			if(PrintLevel > 0) Print("Unsuccesfully Exiting NashEquilibrium: Too Many Iterations");

		}
	}

	void GraphUtil(string FileName)
	{
		/*
		
		Graph the utility of each player,
		varying the player's strategy for Low(j) to High(j),
		and holding the other playres strategies constant.
		Save the results in a specified file.

		*/
		
		Vector<double> fCurr(J);	fCurr.Name = "NashEquilibrium::GraphUtil::fCurr";
		Vector<double> xCurr(J);	xCurr.Name = "NashEquilibrium::GraphUtil::xCurr";
		Matrix<double> XGrid(J,K);	XGrid.Name = "NashEquilibrium::GraphUtil::XGrid";
		Matrix<string> Util(J,K);	Util.Name  = "NashEquilibrium::GraphUtil::Util";
		NashEquilibriumSpace::J = J;
		NashEquilibriumSpace::x = x;
		NashEquilibriumSpace::Func = Func;
		for(int j=0; j<J; j++)
		{
			NashEquilibriumSpace::j = j;
			Copy(xCurr,x);
			for(int k=0; k<K; k++)
			{
				XGrid(j,k) = Low(j) + (High(j) - Low(j)) * k / (K - One);
				xCurr(j) = XGrid(j,k);
				int Info = Func(xCurr,fCurr);
				if(Info) Util(j,k) = "";
				else     Util(j,k) = ToString(fCurr(j));
			}
		}

		/* Write the File */
		ofstream out1(FileName.c_str());
		for(int j=0; j<J-1; j++) out1 << "X(" + ToString(j+1) + ")\tUtil(" + ToString(j+1) + ")\t";
		out1 << "X(" + ToString(J) + ")\tUtil(" + ToString(J) + ")\n";
		for(int k=0; k<K; k++)
		{
			for(int j=0; j<J-1; j++) out1 << XGrid(j,k) << "\t" << Util(j,k) << "\t";
			out1 << XGrid(J-1,k) <<"\t" << Util(J-1,k) << "\n";
		}
	}
	
	
	void GraphBestResponse(string FileName)
	{
	
		/* Graph the Best Response Functions */
	
		if(J != 2)
		{
			Print("Error in NashEquilibrium");
			Print("Graph Best Response Only Works for Two Player Games");
			HaltExecution();
		}

		Vector< Vector<double> > xGrid(J);
		xGrid.Name = "NashEquilibrium::xGrid";
		for(int j=0; j<J; j++)
		{
			xGrid(j) = Grid(Low(j),High(j),K);
			xGrid(j).Name = "xGrid(" + ToString(j+1) + ")";
		}
		
		NashEquilibriumSpace::J = J;
		NashEquilibriumSpace::x = Vector<double>(J);

		Vector< Vector<double> > Br(J);
		Br.Name = "NashEquilibrium::Br";
		for(int j=0; j<J; j++)
		{
			Br(j) = Vector<double>(K);
			Br(j).Name = "Br(" + ToString(j+1) + ")";
		}

		for(int j=0; j<J; j++)
		{
			Copy(NashEquilibriumSpace::x,x);
			for(int k=0; k<K; k++)
			{
				NashEquilibriumSpace::j      = j;
				NashEquilibriumSpace::x(1-j) = xGrid(j)(k);
				NashEquilibriumSpace::Func   = Func;

				Min1D m1(NashEquilibriumSpace::Func_j);
				m1.Solve(xGrid(j));
				Br(j)(k) = m1.x;
			}
		}

		Output(FileName,xGrid(0),Br(0),xGrid(1),Br(1));
	}

	bool Check()
	{
		/* Compute Equilbrium Values */

		Vector<double> fEq(J);
		fEq.Name = "NashEquilibrium1::Check::fEq";
		Func(x,fEq);

		if(PrintLevel >= 1) Print("Checking Equilbrium");
		Print("xEq",x);
		Print("fEq",fEq);

		Vector<double> fCurr(J);					fCurr.Name = "NashEquilibrium::Check::fCurr";
		Vector<double> xCurr(J);					xCurr.Name = "NashEquilibrium::Check::xCurr";
		Vector< Vector<double> > xGrid(J);
		for(int j=0; j<J; j++) xGrid = Grid(Low(j),High(j),K);
	
		bool Info = True;

		for(int j=0; j<J; j++)
		{
			for(int k=0; k<K; k++)
			{
				Copy(xCurr,x);
				xCurr(j) = xGrid(j)(k);
				Func(xCurr,fCurr);

				if(fCurr[j] > fEq[j])
				{
					if(PrintLevel >= 1) Print("This is not an equilbrium: fCurr[j] > fEq[j]");
					if(PrintLevel >= 1) Print("Player",j+1);
					if(PrintLevel >= 1) Print("k",k);
					if(PrintLevel >= 1) Print("xGrid(j)(k)",xGrid(j)(k));
					if(PrintLevel >= 1) Print("fCurr(j)",fCurr(j));
					if(PrintLevel >= 1) Print("fEq(j)",fEq(j));
					Info = False;
				}
			}
		}
		if(PrintLevel >= 1) Print("Checked Equilbrium!");
		
		return Info;
	}


};

#endif