#ifndef TACHYON_GEM_H
#define TACHYON_GEM_H

/*Drawing Layouts:
	Reingold-Tilford
	Random 3D Layout
	Barycentric
	Fruchterman-Reingold
	GEM Force-directed
	Ring
	UserDefined
*/

#include "Vector2.h"

/*typedef Vector2<double> arrow;

struct graph_embedding_data {
	arrow currentPosition;		// E
	arrow lastImpulse;			// p
	double localTemperature;	// t
	double skewGauge;			// d
};

typedef std::vector<graph_embedding_data> gem_data;

double gem_t_global(const gem_data &vertex_data) {
	double t_total(0);
	for (gem_data::const_iterator i = vertex_data.begin(); i != vertex_data.end(); ++i)
		t_total += i->localTemperature;
	return t_total / vertex_data.size();
}

arrow gem_position_sum(const gem_data &vertex_data) {
	arrow c(0);
	for (gem_data::const_iterator i = vertex_data.begin(); i != vertex_data.end(); ++i)
		c += i->currentPosition;
	return c;
}

template <typename Graph>
double gem_functor(Graph &g, uint32_t v) {
	uint32_t degree = boost::degree(v, g);
	return degree * (1 + degree / 2.0);
}

template <typename Graph>
arrow gem_impulse(Graph &g, const gem_data &vertex_data, uint32_t v, arrow c) {
	static const double e_des(128);
	static const double gravity(0.0625);

	const graph_embedding_data &data = vertex_data[v];

	arrow p = (c / vertex_data.size() - data.currentPosition) * gravity * gem_functor(g, v);
	p += arrow(rand() * 64.0 / RAND_MAX - 32, rand() * 64.0 / RAND_MAX - 32);

	for (gem_data::const_iterator u = vertex_data.begin(); u != vertex_data.end(); ++u) {
		arrow delta = data.currentPosition - u->currentPosition;
		if (delta != arrow(0)) {
			double multiple = e_des / delta.Length();
			p += delta * multiple * multiple;
		}
	}

	boost::graph_traits<Graph>::out_edge_iterator oei, oedge_end;
	for (boost::tie(oei, oedge_end) = boost::out_edges(boost::vertex(v, g), g); oei != oedge_end; ++oei) {
		arrow delta = data.currentPosition - vertex_data[boost::target(*oei, g)].currentPosition;
		double multiple = delta.Length() / e_des;
		p -= delta * multiple * multiple / gem_functor(g, v);
	}

	boost::graph_traits<Graph>::in_edge_iterator iei, iedge_end;
	for (boost::tie(iei, iedge_end) = boost::in_edges(boost::vertex(v, g), g); iei != iedge_end; ++iei) {
		arrow delta = data.currentPosition - vertex_data[boost::source(*iei, g)].currentPosition;
		double multiple = delta.Length() / e_des;
		p -= delta * multiple * multiple / gem_functor(g, v);
	}

	return p;
}*/

/*double gem_incidence(arrow v1, arrow v2) {
	return acos(v1.Dot(v2) / (v1.Length() * v2.Length()));

	/* // c^2 = a^2 + b^2 -2*a*b*cos(C);
	double a = v1.Length();
	double b = v2.Length();
	double c = (v1-v2).Length();

	double dude = pow(c, 2) - pow(a, 2) - pow(b, 2) / (-2 * a * b);
	return acos(dude);*/
//}

/*void gem_update(gem_data &vertex_data, uint32_t v, arrow p, double t_max) {
	static const double pi(3.14159);

	static const double alpha_o(pi);
	static const double alpha_r(pi / 3);

	static const double sigma_o(1.0 / 3.0);
	const double sigma_r(1.0 / (2 * vertex_data.size()));

	graph_embedding_data &data = vertex_data[v];

	if (p != arrow(0)) {
		p = data.localTemperature * p / p.Length();
		data.currentPosition += p;
		// Update c here and save a second.
	}

	if (data.lastImpulse != arrow(0)) {
		double beta = gem_incidence(p, data.lastImpulse);

		double sin_beta = sin(beta);
		double cos_beta = cos(beta);

		if (sin_beta >= sin((pi + alpha_r) / 2))
			data.skewGauge += sigma_r * (sin_beta / fabs(sin_beta));
		if (fabs(cos_beta) >= cos(alpha_o / 2))
			data.localTemperature *= sigma_o * cos_beta;
		data.localTemperature *= 1 - fabs(data.skewGauge);
		data.localTemperature = min(data.localTemperature, t_max);
	}

	data.lastImpulse = p;
}

template <typename Graph>
void graph_embedding(Graph &g, uint32_t r_max, double t_max, double t_min) {
	static const double t_init(100);

	srand(clock());

	gem_data vertex_data(boost::num_vertices(g));

	int b(0);
	boost::graph_traits<Graph>::vertex_iterator vi, vertex_end;
	for (boost::tie(vi, vertex_end) = boost::vertices(g); vi != vertex_end; ++vi) {
		graph_embedding_data &data = vertex_data[*vi];
		data.lastImpulse = 0;
		data.skewGauge = 0;
		data.localTemperature = t_init;
		data.currentPosition = arrow(b += 10, b);
	}

	typedef std::vector<uint32_t> ordering;
	ordering vertex_order(boost::num_vertices(g));

	for (uint32_t i(0); i < vertex_order.size(); ++i)
		vertex_order[i] = i;

	double t_global(t_max);
	for (uint32_t r(0); t_global > t_min && r < r_max; ++r) {
		std::random_shuffle(vertex_order.begin(), vertex_order.end());
		for (ordering::const_iterator v = vertex_order.begin(); v != vertex_order.end(); ++v) {
			arrow p = gem_impulse(g, vertex_data, *v, gem_position_sum(vertex_data));
			gem_update(vertex_data, *v, p, t_max);
			t_global = gem_t_global(vertex_data);
		}
	}

	for (uint32_t iff(0); iff < vertex_data.size(); ++iff) {
		std::wcout << L"Pts(" << (iff*2) << L") = " << vertex_data[iff].currentPosition.x << L" + 100" << std::endl;
		std::wcout << L"Pts(" << (iff*2 + 1) << L") = " << vertex_data[iff].currentPosition.y << L" + 100" << std::endl;
	}
}

template <typename Graph>
void graph_embedding2(Graph &g) {
	graph_embedding(g, 4 * boost::num_vertices(g), 256, 3);
}*/

template <typename T>
inline T min(T v1, T v2) {
	return (v1 > v2) ? v2 : v1;
}

template <typename T>
inline T max(T v1, T v2) {
	return (v1 > v2) ? v1 : v2;
}

typedef int scalar;
typedef Vector2<scalar> arrow;

struct gem_settings {
   float maxtemp;		// 0.01	- 10.0
   float starttemp;		// 0.01	- 10.0
   float finaltemp;		// 0.01	- 0.5
   int maxiter;			// 1	- 100
   float gravity;		// 0.0	- 1.0
   float oscillation;	// 0.0	- 2.0
   float rotation;		// 0.0	- 2.0
   float shake;			// 0.0	- 5.0
};

static const bool randomize(true);
//static const gem_settings ins =	{ 1.00f, 0.3f, 0.05f, 10, 0.05f, 0.4f, 0.5f, 0.2f };
static const gem_settings arr =	{ 1.50f, 1.0f, 0.02f, 100, 0.15f, 0.4f, 0.9f, 0.3f };
static const gem_settings opt =	{ 0.25f, 1.0f, 1.00f,  3, 0.10f, 0.4f, 0.9f, 0.3f };

#define ELEN 128L
#define	ELENSQR		(ELEN*ELEN)
#define	MAXATTRACT	1048576

template <typename Graph>
class graph_embedding {
  private:
	struct vertex_data {
		scalar heat;
		arrow imp;
		int dir;
		int mass;
		arrow pos;
	};

	Graph &g;
	typedef std::vector<vertex_data> VertexData;
	std::vector<vertex_data> data;

	scalar temperature;
	arrow center;
	scalar maxtemp;
	float  oscillation, rotation;

	void init(const scalar starttemp) {
		temperature = 0;
		center = 0;

		boost::graph_traits<Graph>::vertex_iterator vi, vertex_end;
		for (boost::tie(vi, vertex_end) = boost::vertices(g); vi != vertex_end; ++vi) {
			vertex_data &v = data[*vi];

			v.heat = starttemp * ELEN;
			temperature += v.heat * v.heat;

			v.imp = 0;
			v.dir = 0;
			v.mass = 1 + boost::degree(*vi, g) / 3;

			center += v.pos;
		}

		srand(time(NULL));
	}

	/*vertex select (void) {
		register vertex	u, v;
		register int 	n;
		static vertex	map[VMAX+1];

		if (iteration == 0)
			for_all_vertices (v)
				map[v] = v;
		n = number_vertices - iteration % number_vertices;
		v = 1 + rand () % n;
		u = map[v]; map[v] = map[n]; map[n] = u;
		return u;
	}*/

	void displace (Graph::vertex_descriptor vi, arrow i) {
		int number_vertices = boost::num_vertices(g);
		register scalar		t, n;
		register arrow		*imp;

		vertex_data &v = data[vi];

		if (i.x || i.y) {
			n = max(abs(i.x), abs(i.y)) / 16384L;
			if (n > 1) {
				i.x /= n;
				i.y /= n;
			}
			t = v.heat;
			n = i.Length();
			i = i * t / n;
			v.pos += i;
			center += i;
			imp = &v.imp;
			n = t * imp->Length();
			if (n) {
				temperature -= t * t;
				t += t * oscillation * i.Dot(*imp) / n;
				t = min(t, maxtemp);
				v.dir += rotation * (i.x * imp->y - i.y * imp->x) / n;
				t -= t * abs(v.dir) / number_vertices;
				t = max(t, 2);
				temperature += t * t;
				v.heat = t;
			}
			*imp = i;
		}
	}

	arrow a_impulse (Graph::vertex_descriptor vi) {
		int number_vertices = boost::num_vertices(g);
		vertex_data &v = data[vi];

		register arrow	i, d;
		register scalar	n;
		arrow p = v.pos;

		n = arr.shake * ELEN;
		i.x = rand () % (2 * n + 1) - n;
		i.y = rand () % (2 * n + 1) - n;

		i += (center / number_vertices - p) * v.mass * arr.gravity;

		boost::graph_traits<Graph>::vertex_iterator ui, vertex_end;
		for (boost::tie(ui, vertex_end) = boost::vertices(g); ui != vertex_end; ++ui) {
			d = p - data[*ui].pos;
			n = d.Dot();
			if (n != 0)
				i += d * ELENSQR / n;
		}

		boost::graph_traits<Graph>::in_edge_iterator iei, iedge_end;
		for (boost::tie(iei, iedge_end) = boost::in_edges(boost::vertex(vi, g), g); iei != iedge_end; ++iei) {
			d = p - data[boost::source(*iei, g)].pos;
			n = d.Dot() / v.mass;
			n = min(n, MAXATTRACT);
			i -= d * n / ELENSQR;
		}

		boost::graph_traits<Graph>::out_edge_iterator oei, oedge_end;
		for (boost::tie(oei, oedge_end) = boost::out_edges(boost::vertex(vi, g), g); oei != oedge_end; ++oei) {
			d = p - data[boost::target(*oei, g)].pos;
			n = d.Dot() / v.mass;
			n = min(n, MAXATTRACT);
			i -= d * n / ELENSQR;
		}

		return i;
	}

	arrow EVdistance (const Graph::edge_descriptor ei, const Graph::vertex_descriptor vi) {
		arrow a = data[boost::source(ei, g)].pos;
		arrow b = data[boost::target(ei, g)].pos;
		arrow c = data[vi].pos;
		scalar	m, n;

		b.x -= a.x; b.y -= a.y; /* b' = b - a */
		m = b.x * (c.x - a.x) + b.y * (c.y - a.y); /* m = <b'|c-a> = <b-a|c-a> */
		n = b.x * b.x + b.y * b.y; /* n = |b'|^2 = |b-a|^2 */
		if (m < 0) m = 0;
		if (m > n) m = n = 1;
		if (m >> 17)	{  	/* prevent integer overflow */
		n /= m >> 16;
		m /= m >> 16;
		}
		a.x += b.x * m / n;	/* a' = m/n b' = a + m/n (b-a) */
		a.y += b.y * m / n;
		return a;
	}

	arrow o_impulse (Graph::vertex_descriptor vi) {
		int number_vertices = boost::num_vertices(g);
		vertex_data &v = data[vi];

		register arrow	i, d;
		register scalar	n;
		arrow p = v.pos;

		n = opt.shake * ELEN;
		i.x = rand () % (2 * n + 1) - n;
		i.y = rand () % (2 * n + 1) - n;

		i += (center / number_vertices - p) * v.mass * opt.gravity;

		boost::graph_traits<Graph>::edge_iterator ei, edge_end;
		for (boost::tie(ei, edge_end) = boost::edges(g); ei != edge_end; ++ei) {
			Graph::vertex_descriptor ui = boost::source(*ei, g);
			Graph::vertex_descriptor wi = boost::target(*ei, g);
			if (ui != vi && wi != vi) {
				d = (data[ui].pos + data[wi].pos) / 2 - p;
				n = d.Dot();
				if (n < 8 * ELENSQR) {
					d = EVdistance(*ei, vi) - p;
					n = d.Dot();
				}
				if (n != 0)
					i -= d * ELENSQR / n;
			} else {
				if (ui == vi)
					ui = wi;
				d = p - data[ui].pos;
				i -= d * min(d.Dot() / v.mass, MAXATTRACT) / ELENSQR;
			}
		}

		return i;
	}

	void arrange() {
		scalar stop_temperature;
		int stop_iteration;

		int number_vertices = boost::num_vertices(g);

		init(arr.starttemp);
		oscillation = arr.oscillation;
		rotation = arr.rotation;
		maxtemp = arr.maxtemp * ELEN;
		stop_temperature = scalar (arr.finaltemp * arr.finaltemp * ELENSQR * number_vertices);
		stop_iteration = arr.maxiter * number_vertices * number_vertices;

		std::vector<Graph::vertex_descriptor> mapping(number_vertices);
		boost::graph_traits<Graph>::vertex_iterator vi, vertex_end;
		for (boost::tie(vi, vertex_end) = boost::vertices(g); vi != vertex_end; ++vi)
			mapping[*vi] = *vi;

		for (int iteration(0); temperature > stop_temperature && iteration < stop_iteration; ) {
			std::random_shuffle(mapping.begin(), mapping.end());
			for (std::vector<Graph::vertex_descriptor>::iterator vi = mapping.begin(); vi != mapping.end(); ++vi, ++iteration)
				displace(*vi, a_impulse(*vi));
		}
	}

	void optimize() {
		scalar stop_temperature;
		int stop_iteration;

		int number_vertices = boost::num_vertices(g);

		init(opt.starttemp);
		oscillation = opt.oscillation;
		rotation = opt.rotation;
		maxtemp = opt.maxtemp * ELEN;
		stop_temperature = scalar (opt.finaltemp * opt.finaltemp * ELENSQR * number_vertices);
		stop_iteration = opt.maxiter * number_vertices * number_vertices;

		std::vector<Graph::vertex_descriptor> mapping(number_vertices);
		boost::graph_traits<Graph>::vertex_iterator vi, vertex_end;
		for (boost::tie(vi, vertex_end) = boost::vertices(g); vi != vertex_end; ++vi)
			mapping[*vi] = *vi;

		for (int iteration(0); temperature > stop_temperature && iteration < stop_iteration; ) {
			std::random_shuffle(mapping.begin(), mapping.end());
			for (std::vector<Graph::vertex_descriptor>::iterator vi = mapping.begin(); vi != mapping.end(); ++vi, ++iteration)
				displace(*vi, o_impulse(*vi));
		}
	}

  public:
	graph_embedding(Graph &g) :
		g(g), data(boost::num_vertices(g))
	{}

	void embed() {
		int number_vertices = boost::num_vertices(g);

		srand(clock());
		if (randomize) {
			long scale = ELEN * sqrt(number_vertices) / 2;
			boost::graph_traits<Graph>::vertex_iterator vi, vertex_end;
			for (boost::tie(vi, vertex_end) = boost::vertices(g); vi != vertex_end; ++vi) {
				data[*vi].pos.x = (scalar) rand() % (2 * scale + 1) - scale / 2;
				data[*vi].pos.y = (scalar) rand() % (2 * scale + 1) - scale / 2;
			}
		} /*else if (ins.finaltemp < ins.starttemp)
			insert();*/

		if (arr.finaltemp < arr.starttemp)
			arrange();
		if (opt.finaltemp < opt.starttemp)
			optimize();

		for (uint32_t iff(0); iff < data.size(); ++iff) {
			std::wcout << L"Pts(" << (iff*2) << L") = " << data[iff].pos.x << std::endl;
			std::wcout << L"Pts(" << (iff*2 + 1) << L") = " << data[iff].pos.y << std::endl;
		}
	}

};

#endif//TACHYON_GEM_H