#include <iostream>
#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 dcut 1

using namespace std;

void _dco() {
	cout << "dco" << endl;

	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++) {
		x[i] = 2.0;
		dco::ga1s<double>::global_tape->register_variable(x[i]);
	}

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

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

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

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

	for (int i = 0; i < N; i++) {
		cout << "x[" << i << "]_(1) = " << dco::derivative(x[i]) << endl;
	}

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

void _dco_fuzzy() {
	cout << "dco_fuzzy" << endl;

	dco::ga1s< fuzzy<double, ncuts> >::global_tape = dco::ga1s< fuzzy<double, ncuts> >::tape_t::create();

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

	for (int i = 0; i < N; i++) {
		dco::value(x[i]).set_cut(0, interval<double>(2.0));
		dco::value(x[i]).set_cut(1, interval<double>(1.9, 2.1));
		dco::ga1s< fuzzy<double, ncuts> >::global_tape->register_variable(x[i]);
	}

	dco::value(y).set_cut(0, interval<double>(1.0));
	dco::value(y).set_cut(1, interval<double>(1.0));
	
	for (int i = 0; i < N; i++) {
		y = y * x[i];
	}

	cout << "y = " << dco::value(y) << endl;

	dco::ga1s< fuzzy<double, ncuts> >::global_tape->register_output_variable(y);
	dco::derivative(y).set_cut(0, interval<double>(1.0));
	dco::derivative(y).set_cut(1, interval<double>(1.0));

	cout << "y_(1)" << " = " << dco::derivative(y) << endl;

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

	for (int i = 0; i < N; i++) {
		cout << "x[" << i << "]_(1) = " << dco::derivative(x[i]) << endl;
	}

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

void _fuzzy() {
	double h = 0.00001;

	cout << "fuzzy" << endl;

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

	for (int i = 0; i < N; i++) {
		x[i].set_cut(0, interval<double>(2.0));
		x[i].set_cut(1, interval<double>(1.9, 2.1));
	}

	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;

	for (int i = 0; i < N; i++) {
		x[i].set_cut(0, interval<double>(2.0));
		x[i].set_cut(1, interval<double>(1.9, 2.1));
	}

	x[0].set_cut(dcut, interval<double>(1.9 + h));

	y2.set_cut(0, interval<double>(1.0));
	y2.set_cut(1, interval<double>(1.0));
	
	for (int i = 0; i < N; i++) {
		y2 = y2 * x[i];
	}

	cout << "y2 = " << y2 << endl;

	cout << "lower dy / dx = " << ((y2.get_cut(dcut).lower() - y.get_cut(dcut).lower()) / h) << endl;
	cout << "upper dy / dx = " << ((y2.get_cut(dcut).upper() - y.get_cut(dcut).upper()) / h) << endl;

	delete [] x;
}

void _fuzzy_dco() {
	cout << "fuzzy_dco" << endl;

	dco::ga1s<double>::global_tape = dco::ga1s<double>::tape_t::create();

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

	for (int i = 0; i < N; i++) {
		dco::ga1s<double>::type l1, l2, u1, u2;
		l1 = u1 = 2.0;
		l2 = 1.9; u2 = 2.1;

		dco::ga1s<double>::global_tape->register_variable(l1);
		dco::ga1s<double>::global_tape->register_variable(u1);
		dco::ga1s<double>::global_tape->register_variable(l2);
		dco::ga1s<double>::global_tape->register_variable(u2);

		x[i].set_cut(0, interval<dco::ga1s<double>::type>(l1, u1));
		x[i].set_cut(1, interval<dco::ga1s<double>::type>(l2, u2));
	}

	y.set_cut(0, interval<dco::ga1s<double>::type>(1.0));
	y.set_cut(1, interval<dco::ga1s<double>::type>(1.0));
	
	for (int i = 0; i < N; i++) {
		y = y * x[i];
	}

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

	dco::ga1s<double>::type l1;

	l1 = y.get_cut(dcut).lower();

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

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

	for (int i = 0; i < N; i++) {
		cout << "lower x[" << i << "]_(1) = " << dco::derivative(x[i].get_cut(dcut).lower()) << endl;
		cout << "upper x[" << i << "]_(1) = " << dco::derivative(x[i].get_cut(dcut).upper()) << endl;
	}

	dco::ga1s<double>::tape_t::remove(dco::ga1s<double>::global_tape);

	delete [] x;
}

int main(void) {
	_dco();
	_fuzzy();
	_fuzzy_dco();
	_dco_fuzzy();

	return 0;
}
