#include "MpUtil.h"

const int N = 100;
Vector<int> Y(N);

const int R1 = 10;
const int R2 = 100;
const int R3 = 1000;
Vector<double> u(R3);
Vector<double> Eps(Max(R1,R2,R3));

const int I = 11;

int Func1(double Alpha, double &F)
{
	F = Zero;
	for(int n=0; n<N; n++)
	{

		double Temp = Zero;
		for(int r=0; r<R1; r++) if(Eps(r) < -Alpha) Temp += One;
		Temp /= (double)R1;

		if(Y(n) == 1)
		{
			if(Temp <= Zero) return 1;
			F -= log(Temp);
		}
		else
		{
			if(Temp >= One) return 1;
			F -= log(One - Temp);
		}
	}
	return 0;
}

int Func2(double Alpha, double &F)
{
	F = Zero;
	for(int n=0; n<N; n++)
	{

		double Temp = Zero;
		for(int r=0; r<R2; r++) if(Eps(r) < -Alpha) Temp += One;
		Temp /= (double)R2;

		if(Y(n) == 1)
		{
			if(Temp <= Zero) return 1;
			F -= log(Temp);
		}
		else
		{
			if(Temp >= One) return 1;
			F -= log(One - Temp);
		}
	}
	return 0;
}

int Func3(double Alpha, double &F)
{
	F = Zero;
	for(int n=0; n<N; n++)
	{

		double Temp = Zero;
		for(int r=0; r<R3; r++) if(Eps(r) < -Alpha) Temp += One;
		Temp /= (double)R3;

		if(Y(n) == 1)
		{
			if(Temp <= Zero) return 1;
			F -= log(Temp);
		}
		else
		{
			if(Temp >= One) return 1;
			F -= log(One - Temp);
		}
	}
	return 0;
}

int Func4(double Alpha, double &F)
{
	F = Zero;
	for(int n=0; n<N; n++)
	{
		double Temp = TrapRule(NormalPdf,-5.0,-Alpha,I);
		if(Y(n) == 1)
		{
			if(Temp <= Zero) return 1;
			F -= log(Temp);
		}
		else
		{
			if(Temp >= One) return 1;
			F -= log(One - Temp);
		}
	}
	return 0;
}

int Func5(double Alpha, double &F)
{
	F = Zero;
	for(int n=0; n<N; n++)
	{
		double Temp = Zero;
		for(int r=0; r<R3; r++) Temp += One - NormalCdf(-Alpha - (Eps(r) - u(r)),0.0,1.0);
		Temp /= (double)R3;

		if(Y(n) == 1)
		{
			if(Temp <= Zero) return 1;
			F -= log(Temp);
		}
		else
		{
			if(Temp >= One) return 1;
			F -= log(One - Temp);
		}
	}
	return 0;
}


int Func6(double Alpha, double &F)
{
	F = Zero;
	for(int n=0; n<N; n++)
	{
		double Temp = NormalCdf(-Alpha);
		if(Y(n) == 1)
		{
			if(Temp <= Zero) return 1;
			F -= log(Temp);
		}
		else
		{
			if(Temp >= One) return 1;
			F -= log(One - Temp);
		}
	}
	return 0;
}

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

	long RandomSeed = -23987324;
	double Alpha0 = 0.5;

	/* Simulate Data */
	for(int n=0; n<N; n++)
	{
		double Temp = Alpha0 + NormalSample(RandomSeed);
		if(Temp >= Zero) Y(n) = 1;
		else             Y(n) = 0;
	}

	/* Random Draws */
	for(int r=0; r<Max(R1,R2,R3); r++)
	{
		Eps(r) = NormalSample(RandomSeed);
	}
	for(int r=0; r<R3; r++)
	{
		u(r) = NormalSample(0.0,1.0,RandomSeed);
	}

	Plot1D(Func1,-5.0,5.0,1001,"plot1.dat");
	Plot1D(Func2,-5.0,5.0,1001,"plot2.dat");
	Plot1D(Func3,-5.0,5.0,1001,"plot3.dat");
	Plot1D(Func4,-5.0,5.0,1001,"plot4.dat");
	Plot1D(Func5,-5.0,5.0,1001,"plot5.dat");
	Plot1D(Func6,-5.0,5.0,1001,"plot6.dat");

	/*
	Min1D m1(Func1);
	m1.PrintLevel = 1;
	int Info1 = m1.Solve(-1.0,0.0,1.0);
	Print("Info1",Info1,"m1.x",m1.x);

	Min1D m2(Func2);
	m2.PrintLevel = 1;
	int Info2 = m2.Solve(-1.0,0.0,1.0);
	Print("Info2",Info2,"m2.x",m2.x);

	Min1D m3(Func3);
	m3.PrintLevel = 1;
	int Info3 = m3.Solve(-1.0,0.0,1.0);
	Print("Info3",Info3,"m3.x",m3.x);

	Min1D m4(Func4);
	m4.PrintLevel = 1;
	int Info4 = m4.Solve(-1.0,0.0,1.0);
	Print("Info4",Info4,"m4.x",m4.x);
	*/


	Print("End");
	return 0;
}