#include "MpUtil.h"
#include "NashEquilibrium.cpp"

double Rho;
double Lambda;
double Sigma;
double qL;
double qR;

int Func1(const Vector<double> &y, Vector<double> &Util)
{

	double PL;
	if(y(0) < y(1))
	{
		PL = NormalCdf((Lambda + Sqr(y(1)) - Sqr(y(0))) / (Two * (y(1) - y(0)) * Sigma));
	}
	else if(y(1) < y(0))
	{
		PL = One - NormalCdf((Lambda + Sqr(y(1)) - Sqr(y(0))) / (Two * (y(1) - y(0)) * Sigma));
	}
	else
	{
		     if(Lambda > Zero) PL = One;
		else if(Lambda < Zero) PL = Zero;
		else                   PL = Half;
	}

	Util(0) =  PL * (Rho - Sqr(y(0) - qL)) - (One - PL) *  Sqr(y(1) - qL);
	Util(1) = -PL *  Sqr(y(0) - qR)        + (One - PL) * (Rho - Sqr(y(1) - qR));

	return 0;
}

int main()
{
	Print("Begin");

	Rho = 1.0;
	Lambda = 0.0;
	Sigma = 0.5;
	qL = -3.0;
	qR =  3.0;

	NashEquilibrium ne(Func1,2);
	ne.Option = 2;
	Set(ne.Low,-5.0);
	Set(ne.High,5.0);
	ne.K = 101;
	ne.Solve();
	Print("ne.x",ne.x);
	Print("ne.Check()",ne.Check());

	/*
	ne.GraphUtil("util.dat");
	ne.GraphBestResponse("br.dat");

	int I = 101;
	Vector<double> yLStar(I);	yLStar.Name = "yLStar";
	Vector<double> yRStar(I);	yRStar.Name = "yRStar";
	Vector<double> Grid1 = Grid(-2.0,2.0,I);
	Grid1.Name = "Rho";

	for(int i=0; i<I; i++)
	{
		Lambda = Grid1(i);
		ne.Solve();
		yLStar(i) = ne.x(0);
		yRStar(i) = ne.x(1);
	}

	Output("LambdaPlot.dat",Grid1,yLStar,yRStar);

	*/

	Print("End");
	return 0;
}