#include <functional>
#include <iostream>
#include <math.h>
#include <vector>
#include <stdlib.h>
#include <time.h>

using namespace std;

const double zero = 0.0;
const double nan = 1/zero;
const double pi = acos(-1);
const double minDouble = 	0.00000000000001;


//
//	ApproxFunc is a template class that represents both a function and an approximation
//  It is a generic solution to performing multidimensional approximations using a recursive 
// algorithm.
//
template <class approx, class multiDFunc>
class ApproxFunc : public approx, public multiDFunc
{
public:
	ApproxFunc(double lowerBound, double upperBound, multiDFunc* function, int dimensions=2)	
		: approx(lowerBound, upperBound,function), multiDFunc(*function),  current_d(dimensions), store_vals(1)
	{
		f=this;
		if (a >=b) cerr << "Lowerbound " << a << " is greater than or equal to upper bound " << b << endl;
	};
	virtual double calculate()
	{
		double ret;
		if (current_d == 1)
		{
			store_vals=0;
			ret = approx::calculate();
			store_vals=1;
		}
		else
		{
			current_d--;
			ret =  approx::calculate();
			current_d++;
		}
		return ret;
	}
	virtual double operator()(double x_)
	{
		if (store_vals)
		{
			x[current_d-1]=x_;
			return calculate();
		}
		else 
		{
			return multiDFunc::operator ()(x_);
		}
	};
protected:
	int current_d;
	int store_vals;
};


class Function : public unary_function<double,double>
{
public:
	Function() {};
	virtual double operator()(double x){ return x;};
	virtual Function& getDerivative() { return *this;};
};

//
//	MultiDFunction is used as a base class for functions that wish to be incorporated into
//  the ApproxFunc template.  It stores a vector of double values that represent chosen
// points in the previous dimensions of the recursion.   These are available during the
// actual calculation of a point.
//
class MultiDFunction :public Function
{
public:
	MultiDFunction(int dimension) :x(d-1), d(dimension) {};
	virtual double operator()(double x, double y, ...) = 0;
	~MultiDFunction() {};
protected:
	int d; 
	vector<double> x;
};

class ElectrostaticPotential: public MultiDFunction
{
public:
	ElectrostaticPotential(double at_x, double at_y, double singularAvoidance = minDouble)
		: xp(at_x), yp(at_y), sinc(singularAvoidance), MultiDFunction(2) 
	{};
	ElectrostaticPotential(const ElectrostaticPotential& ep)
		: xp(ep.get_xp()), yp(ep.get_yp()), sinc(ep.get_singularAvoidance()), MultiDFunction(2)
	{};
	virtual double operator()(double val)
	{
		if ((fabs(x[0] - xp) < sinc) && (fabs(val-yp) < sinc)) 
			return  1/(sqrt(2)*sinc); //top off singularities //nan;  
		return sqrt(1/((x[0]-xp)*(x[0]-xp)+(val-yp)*(val-yp)));
	};
	virtual double operator()(double x, double y, ...)
	{
		if ((fabs(x - xp) < sinc) && (fabs(y-yp) < sinc)) 
		{
			cerr << "singularity found in ElectrostaticFunction at (" << x << "," << y << ")" << endl;
			return  1/(sqrt(2)*sinc); //top off singularities //nan;  
		}
		return sqrt(1/((x-xp)*(x-xp)+(y-yp)*(y-yp)));
	};

	double get_xp() const {return xp;};
	double get_yp() const {return yp;};
	double get_singularAvoidance() const {return sinc;};
protected:
	double xp;
	double yp;
	double sinc;
};


class Approximation
{
public:
	Approximation(double lowerBound, double upperBound, Function* function)
		:a(lowerBound), b(upperBound), f(function)
	{};

	virtual double calculate(double precision=0, int maxTries = 3000)
	{
		cerr << "Default approximation not implemented" << endl;
		return -1;
	}
protected:
	double a;
	double b;
	Function* f;
	//TBD: set precision & maxTries?
};

class MonteCarlo : public Approximation
{
public:
	MonteCarlo(double lowerBound, double upperBound, Function * function)	
		: Approximation(lowerBound, upperBound,function)
	{
		if (a >=b) cerr << "Lowerbound " << a << " is greater than or equal to upper bound " << b << endl;
	};

	virtual double calculate(double precision=0.07, int maxTries = 300000)
	{
		double sum=0, val;
		double r;
		double n;
		double N = 1/(precision*precision);
		double randmult = (b-a)/RAND_MAX;

		for (n=0;n < N;n++)
		{
			r = rand()*randmult;
			//cerr << n << "\t" << r << endl; //test quality of random number generator
			val = (*f)(r+a);
			if (val == nan)
			{
				n--;
				cerr << "Monte Carlo got  NAN (" << (*(ElectrostaticPotential*)f).get_xp() << "," << (*(ElectrostaticPotential*)f).get_yp() << ") " << endl;
			}
			else sum += val;
		}
		return (sum/n)*(b-a);
	}
};

class SimpsonsApproximation : public Approximation
{
public:
	SimpsonsApproximation(double lowerBound, double upperBound, Function* function)	
		: Approximation(lowerBound, upperBound,function)
	{
		if (a >=b) cerr << "Lowerbound " << a << " is greater than or equal to upper bound " << b << endl;
	};

	virtual double calculate(double precision=0.005, int maxTries=20)
	{
		if (b<=a+precision) return 0;
		double i, sum0, val, val0;
		double h = (b-a)/3;
		double se= val0 = (*f)(a+2*h),	so = (*f)(a+h) + (*f)(b-h),	sd = (*f)(a) + (*f)(b);
		
		double sum;
		
		double singularAvoidance = (b-a)/(3 * pow(2,maxTries));
		while ((so == nan) || (se==nan) || (sd == nan))
		{ 
			//cerr << "Modifying bounds" << endl;
			a += singularAvoidance;
			b -= singularAvoidance;
			h = (b-a)/3;
			se= val0 = (*f)(a+2*h);
			so = (*f)(a+h) + (*f)(b-h);
			sd = (*f)(a) + (*f)(b);
		} 
		
		sum = h*(sd + 2*se + 4*so)/3;

		do 
		{
			h/=2;
			sum0 = sum;
			se+=so;
			so =0;
			for (i = a+h; i < b; i+=2*h)
			{
				val = (*f)(i);
				if (val == nan)
				{
					so+= (val = val0);
					//cerr << "Function returned NAN for i =" << i << " at (" << (*(ElectrostaticPotential*)f).get_xp() << "," << (*(ElectrostaticPotential*)f).get_yp() << ")" << endl;
				}
				else
					so += (val0=val);
			}
			sum =  h*(sd  + 2*se +4*so)/3;
			//cout << sd << "\t" << se << "\t" << so << "\t" << sum << endl;
		}  while ((fabs((sum - sum0)/sum0) > precision) && (--maxTries));
		if (!maxTries) 
			cerr << "maxTries reached at (" << (*(ElectrostaticPotential*)f).get_xp() << "," << (*(ElectrostaticPotential*)f).get_yp() << ") with precision = "  <<
				fabs((sum - sum0)/sum0)  << endl << " and sum = " << sum <<endl ;
		return sum;
	}
};

class MonteCarlo2dApproximation : public Approximation
{
public:
		MonteCarlo2dApproximation(double lowerBound, double upperBound, MultiDFunction* function)	
		: Approximation(lowerBound, upperBound,function)
	{
		if (a >=b) cerr << "Lowerbound " << a << " is greater than or equal to upper bound " << b << endl;
	};


	virtual double calculate(double precision =0.005, int maxTries=-1)
	{
		double total = 0;
		double N = 1/(precision*precision);
		double randmult = (b-a)/RAND_MAX;
		for (int i=0;i<N;i++)
				total += (*(MultiDFunction*)f)((a+rand()*randmult), (a+rand()*randmult));
			
	return total*(b-a)*(b-a)/(N);
	};
};

int main(int argc, char* argv[])
{
	double inc = 0.1, montesinc=0.0000005,simpsonssinc=0.0000005, x, y;
	unsigned long t0,t1,t2,t3;

	//cout.precision(10);

	cout << "# Welcome to the Electrostatic Potential Integration using the Monte Carlo Technique" << endl;
	cout << "# Current Time : " << (t0  =  time(0)) << endl;
	cout << endl;
	
	cout << "# Monte Carlo" << endl;
	for (x=-5.0;x<=5.0;x+=inc)
	{
		for (y=-5.0;y<=5.0;y+=inc)
		{		
			ElectrostaticPotential *f = new ElectrostaticPotential(x,y);
			ApproxFunc<MonteCarlo, ElectrostaticPotential> full_integral(-1.0,1.0, f);
			cout << x << "\t" << y <<"\t" << full_integral.calculate() <<endl;
			delete f;
		}
	};
	cout << "# Monte Carlo ApproxFunc ended at : " << (t1 = time(0)) << endl;

	
	cout << endl << "# 2-D Simpson Approximation" << endl;
	for (x=-5.0;x<=5.0;x+=inc)
	{
		for (y=-5.0;y<=5.0;y+=inc)
		{		
			ElectrostaticPotential *f = new ElectrostaticPotential(x,y,simpsonssinc);	
			ApproxFunc<SimpsonsApproximation, ElectrostaticPotential> full_integral(-1.0,1.0, f);
			cout << x << "\t" << y <<"\t" << full_integral.calculate() <<endl;
			delete f;
		};
	};
	cout << "# Simpsons ApproxFunc ended at : " << (t2 = time(0)) << endl;

	cout <<  endl << "# Monte Carlo Approximation" << endl;
	for (x=-5.0;x<=5.0;x+=inc)
	{
		for (y=-5.0;y<=5.0;y+=inc)
		{		
			ElectrostaticPotential *f = new ElectrostaticPotential(x,y,montesinc);
			MonteCarlo2dApproximation full_integral(-1.0,1.0, f);
			cout << x << "\t" << y <<"\t" << full_integral.calculate(0.005) <<endl;
			delete f;
		}
	};
	cout << endl << "# Execution completed at time : " << (t3 = time(0)) << endl << endl;
	
	cout <<endl << "# MonteCarlo ApproxFunc: " << (t1 -t0) << endl;
	cout << "# Simpsons ApproxFunc: " << (t2 - t1) << endl;
	cout << "# MonteCarlApproximation : " << (t3 - t2) << endl;
	
	return 0;
};