#include <iostream>
#include <fstream>
#include <cmath>
#include <ctime>
#include "fuzzy_lib/fuzzy_lib.h"
#include "/home/jonathan/dco_cpp/distributions/050315_v3.1/Linux/dco_unscrambled.hpp"

#define N 10
#define ncuts 2
#define mcsamp 100000

using namespace std;

void _monte_carlo(const double *m_x, const double *v_x, double &m_y, double &v_y) {
	double a, b, c, d;
	double mc;
	double *x = new double[N];
	double *y = new double[mcsamp];

	srand(time(0));

	ofstream fout;
	fout.open("mc1.csv");

	double t_v_x_0 = 0.0;

	m_y = 0.0;
	for (int j = 0; j < mcsamp; j++) {
		for (int i = 0; i < N; i++) {
			c = m_x[i];
			d = sqrt(6 * v_x[i]);
			a = c - d;
			b = c + d;

			mc = (double)rand() / RAND_MAX;

			if (mc < (c - a) / (b - a)) {
				x[i] = sqrt(mc * (b - a) * (c - a)) + a;
			} else {
				x[i] = -sqrt(-(mc - 1) * (b - a) * (b - c)) + b;
			}

			fout << x[i] << ", ";
		}

		t_v_x_0 += (x[0] - m_x[0]) * (x[0] - m_x[0]);

		y[j] = 1.0;
		
		for (int i = 0; i < N; i++) {
			y[j] = y[j] * x[i];
		}

		m_y += y[j];
	}

	fout.close();

	m_y /= mcsamp;
	v_y = 0.0;

	for (int j = 0; j < mcsamp; j++) {
		v_y += (y[j] - m_y) * (y[j] - m_y);
	}

	v_y /= mcsamp - 1;

	std::cout << "v_x[0] " << v_x[0] << " t_v_x_0 " << t_v_x_0  / (mcsamp - 1) << std::endl;

	fout.open("mc2.csv");

	for (int j = 0; j < mcsamp; j++) {
		fout << y[j] << ", ";
	}

	fout.close();

	delete [] x;
	delete [] y;
}

void _fuzzy(const double *m_x, const double *v_x, double &m_y, double &v_y) {
	double a, b, c, d;

	fuzzy<double, ncuts> *x = new fuzzy<double, ncuts>[N];
	fuzzy<double, ncuts> y;

	ofstream fout;
	fout.open("abc.csv");

	for (int i = 0; i < N; i++) {
		c = m_x[i];
		d = sqrt(6 * v_x[i]);
		a = c - d;
		b = c + d;

		x[i].set_cut(0, interval<double>(c));
		x[i].set_cut(1, interval<double>(a, b));
	}

	fout << a << ", " << c << ", " << b << endl;

	y.set_cut(0, interval<double>(1.0));
	y.set_cut(1, interval<double>(1.0));

	for (int i = 0; i < N; i++) {
		y = y * x[i];
	}

	cout << "y = " << y << endl;

	a = y.get_cut(1).lower();
	b = y.get_cut(1).upper();
	c = 0.5 * (y.get_cut(0).upper() + y.get_cut(0).lower());

	fout << a << ", " << c << ", " << b << endl;
	fout.close();

	m_y = (a + b + c) / 3;
	v_y = (a * a + b * b + c * c - a * b - a * c - b * c) / 18;

	delete [] x;
}

void _fd_uncertainty(const double *m_x, const double *v_x, double &m_y, double &v_y) {
	double *x = new double[N];
	double y, y2;
	double h = 0.0001;

	for (int i = 0; i < N; i++) {
		x[i] = m_x[i];
	}

	y = 1.0;
	
	for (int i = 0; i < N; i++) {
		y = y * x[i];
	}

	m_y = y;
	v_y = 0.0;

	for (int j = 0; j < N; j++) {
		if (j > 0)
			x[j - 1] -= h;

		x[j] += h;

		y2 = 1.0;
		
		for (int i = 0; i < N; i++) {
			y2 = y2 * x[i];
		}

		double dydxj = (y2 - y) / h;
		v_y += dydxj * dydxj * v_x[j];
	}

	delete [] x;
}

void _dco_uncertainty(const double *m_x, const double *v_x, double &m_y, double &v_y) {
	dco::ga1s<double>::global_tape = dco::ga1s<double>::tape_t::create();

	dco::ga1s<double>::type *x = new dco::ga1s<double>::type[N];
	dco::ga1s<double>::type y;

	for (int i = 0; i < N; i++) {
		dco::value(x[i]) = m_x[i];
		dco::ga1s<double>::global_tape->register_variable(x[i]);
	}

	dco::value(y) = 1.0;
	
	for (int i = 0; i < N; i++) {
		y = y * x[i];
	}

	m_y = dco::value(y);

	dco::ga1s<double>::global_tape->register_output_variable(y);
	dco::derivative(y) = 1.0;

	dco::ga1s<double>::global_tape->interpret_adjoint();

	v_y = 0.0;
	for (int i = 0; i < N; i++) {
		v_y += dco::derivative(x[i]) * dco::derivative(x[i]) * v_x[i];
	}

	ofstream fout;
	fout.open("dco.csv");
	fout << m_x[0] << ", " << v_x[0] << endl;
	fout << m_y << ", " << v_y << endl;
	fout.close();

	dco::ga1s<double>::tape_t::remove(dco::ga1s<double>::global_tape);
	delete [] x;
}

int main(void) {
	double *m_x = new double[N];
	double *v_x = new double[N];
	double m_y, v_y;

	srand(time(0));

	for (int i = 0; i < N; i++) {
		//m_x[i] = 0.5 + 1.0 * ((double)rand() / RAND_MAX);
		//v_x[i] = 0.1 + 0.2 * ((double)rand() / RAND_MAX);
		m_x[i] = 1.0;
		v_x[i] = 0.01;
	}

	cout << "dco_uncertainty" << endl;
	_dco_uncertainty(m_x, v_x, m_y, v_y);
	cout << "m_y = " << m_y << ", v_y = " << v_y << endl;

	cout << "fd_uncertainty" << endl;
	_fd_uncertainty(m_x, v_x, m_y, v_y);
	cout << "m_y = " << m_y << ", v_y = " << v_y << endl;

	cout << "fuzzy" << endl;
	_fuzzy(m_x, v_x, m_y, v_y);
	cout << "m_y = " << m_y << ", v_y = " << v_y << endl;

	cout << "monte_carlo" << endl;
	_monte_carlo(m_x, v_x, m_y, v_y);
	cout << "m_y = " << m_y << ", v_y = " << v_y << endl;

	delete [] m_x;
	delete [] v_x;

	return 0;
}
